From Python

The following examples consider a system composed a protein solvated by a mixture of water and glycerol, built with Packmol. The simulations were performed with NAMD with periodic boundary conditions and a NPT ensemble at room temperature and pressure. Molecular pictures were produced with VMD.

Image of the system of the example: a protein solvated by a mixture of glycreol (green) and water, at a concentration of 50%vv. The complete example is available at this repository.

Loading the ComplexMixtures.py file

The Python interface of ComplexMixtures is implemented in the ComplexMixtures.py file. Just download it from the link and save it in a known path.

Installing juliacall

juliacall is a package that allows calling Julia programs from Python. Install it with

pip install juliacall

Installing Julia and underlying packages

Once juliacall is installed, from within Python, execute:

import ComplexMixtures

here we assume that the ComplexMixtures.py file is in the same directory where you launched Python.

How to run this example

The Data directory contains the a pdb file of the system (system.pdb) and a sample from the trajectory (glyc50.dcd), with a few frames. It also contains the result of running the mddf calculation on the complete trajectory, results_glyc50.json. This last file was produced by ComplexMixtures, as indicated in the following examples.

The sample trajectory is provided so that the first example can be run, yet do not expect that the results are the same, as the sampling is much lower in this case. The complete trajectory can be retrieved from this link (3GB file).

We assume that you navigated to the directory of the example, and copied the Python module file to it:

git clone https://github.com/m3g/ComplexMixturesExamples
cd ComplexMixturesExamples/Protein_in_Glycerol/MDDF
cp /path/to/ComplexMixtures.py ./
export JULIA_NUM_THREADS=8

The last line will allow Julia to execute multi-threaded, which will improve a lot the performance on most machines. Set the number of threads to the number of cores of your computer.

Minimum-Distance Distribuion function

Note that the example here follows an identical syntax to the Julia example, except that we qualify the name of the loaded module and implicitly load the PDBTools package.

The script to compute the MDDFs as associated data from within python is, then:

import ComplexMixtures as cm

# Load the pdb file of the system using `PDBTools`:
atoms = cm.readPDB("../Data/system.pdb")

# Create arrays of atoms with the protein and Glycerol atoms, 
# using the `select` function of the `PDBTools` package:
protein = cm.select(atoms,"protein")
glyc = cm.select(atoms,"resname GLYC")

# Setup solute and solvent structures, required for computing the MDDF, 
# with `Selection` function of the `ComplexMixtures` package:
solute = cm.Selection(protein,nmols=1)
solvent = cm.Selection(glyc,natomspermol=14)

# Read and setup the Trajectory structure required for the computations:
trajectory = cm.Trajectory("../Data/glyc50_complete.dcd",solute,solvent)

# Run the calculation and get results:
results = cm.mddf(trajectory)

# Save the reults to recover them later if required
cm.save(results,"./glyc50.json")

options = cm.Options(cutoff=10.)
results = cm.mddf(trajectory,options)

The trajectory that was loaded was for a toy-example. The complete trajectory is available here, but it is a 3GB file. The same procedure above was performed with that file and produced the results_Glyc50.json file, which is available in the Data directory here. We will continue with this file instead.

Produce plots

MDDF and Kirkwood-Buff integrals

import ComplexMixtures as cm
import matplotlib.pyplot as plt

# Load the actual results obtained with the complete simulation:
results = cm.load("../Data/results_glyc50.json")

# Plot MDDF and KB
fig, axs = plt.subplots(2)
axs[0].plot(results.d, results.mddf)
axs[0].set(ylabel="MDDF")

# Plot KB integral
axs[1].plot(results.d, results.kb)
axs[1].set(xlabel="distance / Angs", ylabel="MDDF")

plt.savefig("mddf_kb.png")

Atomic contributions to the MDDF

Selecting the atoms corresponding to the hydroxyl groups, and of the aliphatic carbons of Glycerol. Here we list the types of the atoms as specified by the force-field.

import ComplexMixtures as cm
import matplotlib.pyplot as plt

atoms = cm.readPDB("../Data/system.pdb")
protein = cm.select(atoms,"protein")
glyc = cm.select(atoms,"resname GLYC")
solute = cm.Selection(protein,nmols=1)
solvent = cm.Selection(glyc,natomspermol=14)

# load results
results = cm.load("../Data/results_glyc50.json")

# Select atoms by name
hydroxyls = cm.list(["O1","O2","O3","H1","H2","H3"])
aliphatic = cm.list(["C1","C2","HA","HB","HC","HD"])

# Extract the contributions of the groups above
hydr_contributions = cm.contributions(solvent,results.solvent_atom,hydroxyls)
aliph_contributions = cm.contributions(solvent,results.solvent_atom,aliphatic)

# Plot
plt.plot(results.d, results.mddf)
plt.plot(results.d, hydr_contributions)
plt.plot(results.d, aliph_contributions)
plt.xlabel("distance / Angs")
plt.ylabel("MDDF")
plt.savefig("group_contributions.png")