:: _api:
cosolvkit API
The CosolvKit API allows for a more advanced, flexible use of CosolvKit where the user can create their own scripts.
Config class and config.json
CosolvKit implements a Config class to handle the list of setup options. In the data folder a template of the config.json file used to setup the system building is provided. It is also possible to overwrite some of the options from the API:
from cosolvkit.config import Config
config = Config.from_config('config.json')
# Modify the receptor option and the radius value
config.receptor = FALSE
config.radius = 8.5
Creating CosolventMolecules without json files
It is possible to avoid to create a cosolvents.json recipe file (although highly recommended).
from cosolvkit.cosolvent_system import CosolventMolecule
cosolvent_molecules = list()
cosolvent_molecules.append(CosolventMolecule(name="benzene",
smiles="C1=CC=CC=C1",
resname="BEN",
concentration=0.25))
cosolvent_molecules.append(CosolventMolecule(name="methanol",
smiles="CO",
# MET is reserved for residues
resname="MNL",
copies=58))
Building a CosolventSystem
Thanks to the flexible API CosolvKit allows the user to instantiate custom CosolventSystem classes with user prepared proteins:
import json
from cosolvkit.cosolvent_system import CosolventSystem
from openmm.app import Modeller
modeller = Modeller(protein_topology, protein_positions)
with open('cosovlents.json') as fi:
cosolvents = json.load(fi)
with open('forcefields.json') as fi:
forcefields = json.load(fi)
cosolvent_system = CosolventSystem(cosolvents=cosolvents,
forcefields=forcefields,
modeller=modeller)
cosolvent_system.build()
Building a CosolventMembraneSystem
A cosolvent system with membrane can be easily built:
import json
from cosolvkit.cosolvent_system import CosolventSystem
from openmm.app import Modeller
modeller = Modeller(protein_topology, protein_positions)
with open('cosovlents.json') as fi:
cosolvents = json.load(fi)
with open('forcefields.json') as fi:
forcefields = json.load(fi)
cosolvent_system = CosolventMembraneSystem(cosolvents=cosolvents,
forcefields=forcefields,
modeller=modeller,
lipid_type="POPC")
# Or if want to pass a different type of lipids (pre-equilibrated patch needed)
cosolvent_system = CosolventMembraneSystem(cosolvents=cosolvents,
forcefields=forcefields,
modeller=modeller,
lipid_patch_path="path/to/the/patch")
cosolvent_system.build()
Adding repulsive forces in case of aggregation events
Aggregation events can be common for some types of cosolvents, if in doubt, we suggest to run a simulation without custom repulsive forces and inspect the RDF profiles (please refer to the original paper for more details). If aggregation is observed, CosolvKit offers the possibility to add a custom repulsive force between specified residues.
import json
from cosolvkit.cosolvent_system import CosolventSystem, CosolventMolecule
from openmm.app import Modeller
modeller = Modeller(protein_topology, protein_positions)
cosolvent_molecules = list()
cosolvent_molecules.append(CosolventMolecule(name="benzene",
smiles="C1=CC=CC=C1",
resname="BEN",
concentration=0.25))
with open('forcefields.json') as fi:
forcefields = json.load(fi)
cosolvent_system = CosolventSystem(cosolvents=cosolvents,
forcefields=forcefields,
modeller=modeller)
cosolvent_system.build()
cosolvent_system.add_repulsive_forces(["BEN"])
# or you can specify epsilon and sigma parameters of the LJ potential
e = 0.05
s = 9
cosolvent_system.add_repulsive_forces(["BEN"], epsilon=e, sigma=s)
Use custom solvent
CosolvKit offers the possibility of using solvents different from water. In case of water the solvation is done by OpenMM, while for custom cosolvents CosolvKit exploits the same method used to place cosolvent molecules to place solvent molecules (if filling the box with solvent can be pretty slow). This feature of CosolvKit is meant to offer flexibility for different advanced tasks.
The solvent can be specified as SMILES string and the number of molecules requested can be specified optionally.
#... Previous code to create cosolvent system
cosolvent_system.build(solvent_smiles="CO", n_solvent_molecules=350)
Saving topologies and the system
Once the cosolvent system is created and parametrized, it has to be saved for the next steps (likely MD simulation). Depending on what MD engine was selected the format of the topology files can change.
#... Previous code to create and parametrize the cosolvent system
cosolvent_system.save_topology(topology=cosolvent_system.modeller.topology,
positions=cosolvent_system.modeller.positions,
system=cosolvent_system.system,
# Gather the md_format from the config file
simulation_format=config.md_format,
forcefield=cosolvent_system.forcefield)
Run MD simulations with CosolvKit
CosolvKit offers a general and standard protocol to run MD simulations that can be used for the majority of the use cases. The flags run_cosolvent_system and run_md in the Config class take care of building the cosolvent system and using the standard MD protocol to run a simulation.
from cosolvkit.simulation import run_simulation
if config.md_format.upper() != "OPENMM":
# Change the next two lines depending on the simulation_format you chose
topo = os.path.join(config.output, "system.prmtop")
pos = os.path.join(config.output, "system.rst7")
# This is for openmm
pdb = None
system = None
else:
topo = None
pos = None
# This is for openmm
pdb = os.path.join(config.output, "system.pdb")
system = os.path.join(config.output, "system.xml")
if config.md_format.upper() == "OPENMM":
print(f"Starting MD simulation from the files: {pdb}, {system}")
else:
print(f"Starting MD simulation from the files: {topo}, {pos}")
run_simulation(
simulation_format = config.md_format,
topology = topo,
positions = pos,
pdb = pdb,
system = system,
warming_steps = 100000,
simulation_steps = 6250000, # 25ns
results_path = config.output, # This should be the name of system being simulated
seed=None
)
Post processing analysis
CosolvKit offers a very basic package to analyze the results of the MD simulations. In particualr, Radial Distribution Functions (RDFs) of the cosolvent atoms and waters are generated with the respective autocorrelation functions. Furthermore, densities of the specified cosolvent molecules are depicted during the simulation and saved as a PyMol session for further analysis (check the pre-print for more examples of the use of cosolvent densities).
# The whole analysis module relies on the Report class
from cosolvkit.analysis import Report
report = Report(log_file="statistics.csv",
traj_file="trajectory.dcd",
top_file="system.prmtop",
cosolvents_path="cosolvents.json")
# Generate RDF and autocorrelation plots
report.generate_report(out_path="results")
# Generate density files
# analysis_selection_string is a string in MDAnalysis format
# to select specific cosolvents for the densities
report.generate_density_maps(out_path="densities",
analysis_selection_string="")
report.generate_pymol_reports(topology="system.prmtop",
trajectory="trajectory.dcd",
density_files=["map_density_BEN.dx"],
# It's possible to specify PyMol selection string to highlight
# specific residues for that particular density
selection_string="",
out_path="results")