Bug fixes, implicit solvent support, fixes to support a conda build with modeller

Author: Rob-Welch

prepmd/__init__.py

@@ -5,4 +5,5 @@
 from . import util
 from . import prep
 from . import model
-__version__ = "0.2"
+from . import add_modeller_license
+__version__ = "0.3"

prepmd/add_modeller_license.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Change the modeller license key
+"""
+
+import importlib
+from pathlib import Path
+import sys
+
+# note: don't do this
+# i have to do it like this to make the conda build work
+
+def entry_point():
+    key = sys.argv[1]
+    modeller_init_path = Path(importlib.util.find_spec("modeller").origin)
+    modeller_lib = modeller_init_path.parent.parent.absolute() / "modeller" / "config.py"
+    with open(modeller_lib) as file:
+        contents = file.readlines()
+    splitline = contents[1].split("'")
+    contents_1_new = splitline[0] + "'" + key + "'" + splitline[2]
+    contents[1] = contents_1_new
+    with open(modeller_lib, "w") as file:
+        file.writelines(contents)
+    print("Updated modeller license info.")
+
+if __name__ == "__main__":
+    if len(sys.argv) != 2:
+        print("Usage: prep-license LICENSE-KEY")
+        print("Note: this will replace whatever your current license key is.")
+        sys.exit(0)
+    entry_point()

prepmd/get_residues.py

@@ -3,9 +3,12 @@
 """
 Read residue information from structure files
 """
-
+NO_MODELLER = False
+try:
+    from modeller import *
+except:
+    NO_MODELLER = True
 from prepmd import util
-from modeller import *
 
 
 def get_residues_pdb(pdb, code):
@@ -18,6 +21,8 @@ def get_residues_pdb(pdb, code):
     Returns:
         the fasta sequence as a string
     """
+    if NO_MODELLER:
+        raise ImportError("Can't run without MODELLER and a valid license key")
     log.none()
     e = Environ()
     m = Model(e, file=pdb)
@@ -80,5 +85,6 @@ def get_fullseq_pdb(pdb, code):
     fasta_joined = "/".join(fastas)
     chains = ":::::::::"
     if fastas == []:
-        raise IOError("Couldn't get full sequence from contents of "+pdb)
+        raise IOError("Couldn't get full sequence from contents of "+pdb+". "
+                      "Does it contain a sequence?")
     return ">P1;"+code+"_fill\n"+chains+"\n"+fasta_joined+"*"

prepmd/model.py

@@ -3,9 +3,14 @@
 """
 Calls to the MODELLER API
 """
-import modeller
+NO_MODELLER = False
+try:
+    import modeller
+    from modeller.automodel import *
+except:
+    NO_MODELLER = True
+
 from Bio import Align
-from modeller.automodel import *
 import os
 import pathlib
 import shutil
@@ -102,10 +107,12 @@ def get_best_pdb(directory, exts=["pdb", "cif", "mmcif", "mmCif"]):
     Returns:
         path to the file with the highest objective function, a string
     """
+    
     pdbs = []
     for ext in exts:
         pdbs += list(pathlib.Path(directory).glob('*.'+ext))
     scores = {}
+    similarities = {}
     for pdb in pdbs:
         with open(pdb, "r") as file:
             for line in file:
@@ -117,6 +124,26 @@ def get_best_pdb(directory, exts=["pdb", "cif", "mmcif", "mmCif"]):
                 if "_modeller.objective_function" in line:
                     score = line.split()[-1]
                     scores[pdb] = float(score)
+                # pdb
+                if "BEST TEMPLATE" in line:
+                    similarities[pdb] = float(line.split()[-1])
+                # mmcif
+                if "_modeller.best_template_pct_seq_id" in line:
+                    similarities[pdb] = float(line.split()[-1])
+    max_sim = max(similarities.values())
+    if max_sim >= 97:
+        print("Similiarity: "+str(max_sim)+"%.")
+    if max_sim < 97 and max_sim >= 90:
+        print("WARNING: best sequence similarity is "+str(max_sim)+"%. "
+              "Modelled loops could be wrong. Please check them!")
+    if max_sim < 90 and max_sim >= 80:
+        print("SERIOUS WARNING: best sequence similarity is only "+
+              str(max_sim)+"%. Modelled loops may have the wrong sequence." )
+    if max_sim < 80:
+        raise ValueError("Sequences do not match. Cannot reconstruct missing "
+                         "loops. Please double-check your sequence data (by "
+                         "default these are the SEQRES records from the input "
+                         "structure.")
     return str(max(scores, key=scores.get))
 
 
@@ -136,6 +163,9 @@ def fix_missing_residues(code, fastafile, alignmentout, inmodel, outmodel,
     Returns:
         nothing, but writes out outmodel and wdir.
     """
+    if NO_MODELLER:
+        raise ImportError("Can't fix without MODELLER and a valid license key")
+
     print("Fixing missing residues...")
 
     if not os.path.isdir(wdir):

prepmd/prep.py

@@ -100,14 +100,23 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None,
     locals_json = json.dumps(locals_copy)
 
     # infer download format from output format
-    if not download_format:
+    if not download_format and not inmodel:
         if (".pdb") in outmodel:
             download_format = "pdb"
             print("No download format specified, downloading PDB.")
         elif (".cif") in outmodel or ".mmcif" in outmodel or (
                 ".pdbx") in outmodel:
             download_format = "mmCif"
             print("No download format specified, downloading mmCif.")
+    
+    #infer download format from input format
+    if not download_format and inmodel:
+        if (".pdb") in inmodel:
+            download_format = "pdb"
+        if ".cif" in inmodel or ".pdbx" in inmodel or ".mmcif" in inmodel:
+            download_format = "mmCif"
+
+        
 
     # check that input/output are specified in the same format
     # i'm not against converting the files but it shouldn't happen implicitly
@@ -126,7 +135,10 @@ def in_string(substr, text): return text == None or substr in text.lower()
         pathlib.Path(workingdir).mkdir(parents=True, exist_ok=True)
 
     if inmodel:
-        shutil.copyfile(inmodel, workingdir+os.path.sep+inmodel)
+        shutil.copyfile(inmodel, workingdir+os.path.sep+code+"."+download_format)
+        # note: modeller requires the filename to be the same as the pdb code
+        # so here we change the filename
+        inmodel = code+"."+download_format.replace("mmCif", "cif")
 
     run_dir = os.getcwd()
     os.chdir(workingdir)

prepmd/run.py

@@ -31,17 +31,22 @@
     "amber14,None": ['amber14-all.xml'],
     "charmm36,None": ['charmm36.xml'],
     "amoeba,None": ['amoeba2018.xml'],
-    "amber19,None": ['amber19-all.xml']
+    "amber19,None": ['amber19-all.xml'],
+    "amber19,implicit": ['amber19-all.xml', 'implicit/gbn2.xml'],
+    "amber14,implicit": ['amber14-all.xml', 'implicit/gbn2.xml'],
+    "charmm36,implicit": ['charmm36.xml', 'implicit/gbn2.xml'],
+    "amoeba,implicit": ['amoeba2018.xml', 'amoeba2009_gk.xml'],  
 }
 
 
 def make_forcefield(ff, solvent=None):
     if ff+","+str(solvent) not in ff_lookup.keys():
+        err = "Solvent "+str(solvent)+" not in list of available solvents."
         validff = set([item[0] for item in list(
             map(methodcaller("split", ","), ff_lookup.keys()))])
         validsolvent = set([item[1] for item in list(
             map(methodcaller("split", ","), ff_lookup.keys()))])
-        err = "Please pick a valid force field and solvent.\n"
+        err += "Please pick a valid force field and solvent.\n"
         err += "Valid force fields are "+str(validff)+"\n"
         err += "Valid solvents are "+str(validsolvent)+"\n"
         raise ValueError(err)
@@ -64,8 +69,8 @@ def check_platforms():
 
 def test_sim(pdb):
     run(pdb, minimised_structure_out=pdb, md_steps=None,
-        integrator="LangevinMiddleIntegrator", solvent=None,
-        pressure=None, test_sim_steps=50)
+        integrator="LangevinMiddleIntegrator",
+        pressure=None, test_sim_steps=250)
 
 
 def run(pdb,
@@ -84,7 +89,7 @@ def run(pdb,
         test_run=True,
         fix_backbone=False,
         constraints="Default",  # Default, None, HBonds
-        solvent="tip4pew",
+        solvent="implicit",
         strip_solvent=False,
         ionic_strength=0.1*molar,
         pressure=1*bar,
@@ -133,7 +138,8 @@ def run(pdb,
         None, HBonds. By default, HBonds will be constrained if the backbone
         is not fixed. A string.
         solvent - solvent to use. Possible values: None (no solvent), tip3p, 
-        tip4pew, spce. A string.
+        tip4pew, spce, implicit (gbn model, equivalent to AMBER igb=8). A
+        string.
         write_solvent: whether to write solvent atoms to the minimised
         structure file. Solvent will always be written to the trajectory. A
         bool.
@@ -195,7 +201,7 @@ def run(pdb,
                              "running a solvated simulation instead.")
 
     if not solvent:
-        print("WARNING: No solvent")
+        print("WARNING: No solvent. Minimised structures may be wrong!!!")
 
     # read in mmcif or pdb
     if ".cif" in pdb or ".mmcif" in pdb:
@@ -209,7 +215,7 @@ def run(pdb,
 
     forcefield = make_forcefield(forcefield, solvent)
 
-    if solvent:
+    if solvent and solvent != "implicit":
         print("Solvating system...")
         modeller.addSolvent(forcefield, model=solvent,
                             ionicStrength=ionic_strength,
@@ -231,7 +237,7 @@ def run(pdb,
 
     # setup non bonded method - if none is chosen, select based on box
     if non_bonded_method == "Default":
-        if solvent is None:
+        if solvent == None or solvent == "implicit":
             non_bonded_method = CutoffNonPeriodic
         else:
             non_bonded_method = PME
@@ -282,6 +288,9 @@ def run(pdb,
         raise ValueError("Integrator must be one of: "
                          "VariableLangevinIntegrator,"
                          " LangevinMiddleIntegrator")
+    
+    # TODO: create metadynamics object here
+    
     simulation = Simulation(modeller.topology, system, integrator)
     simulation.context.setPositions(modeller.positions)
 
@@ -336,6 +345,8 @@ def run(pdb,
                              " minimisation and test run were both skipped.")
 
     if md_steps and traj_out:
+        # todo: if metadynamics, add rmsd reporter here
+        
         if ".dcd" in traj_out.lower():
             simulation.reporters.append(DCDReporter(traj_out, step))
         elif ".xtc" in traj_out.lower():
@@ -350,6 +361,9 @@ def run(pdb,
                                                           step, step=True,
                                                           potentialEnergy=True,
                                                           temperature=True))
+
+    # TODO: if metadynamics, do meta.step here
+
         print("Running production MD...")
         simulation.step(md_steps)
         simulation.saveCheckpoint(checkpoint_output)
@@ -394,7 +408,7 @@ def entry_point():
                         help="Fix the backbone", action="store_true")
     parser.add_argument("-c", "--constraints", help="Constraints to apply. Possible values: None, HBonds. Constraining HBonds can make simulations run faster, at the cost of some accuracy. You can't constrain HBonds if you're also fixing the backbone, for obvious reasons.", default="Default")
     parser.add_argument(
-        "-solv", "--solvent", help="Solvent to use. Possible values: tip3p, tip4pew, spce. If no solvent is specified, none will be used.", default=None)
+        "-solv", "--solvent", help="Solvent to use. Possible values: tip3p, tip4pew, spce, implicit (GBn model, equivalent to AMBER igb=8). If no solvent is specified, implicit will be used.", default="implicit")
     parser.add_argument("-ns", "--strip_solvent",
                         help="Don't write solvent molecules to the minimised output file", action="store_true")
     parser.add_argument("-ion", "--ionic_strength",

pyproject.toml

@@ -5,16 +5,17 @@ requires-python = ">= 3.13.3"
 authors = [{name = "Rob Welch", email = "robert.welch@stfc.ac.uk"}]
 description = "Automatically prepare pdb/mmCif structures for molecular dynamics simulations"
 dependencies = [
-  "requests",
-  "openmm==8.2.0",
-  "pdbfixer==1.11",
-  "gemmi==0.7.1",
-  #"modeller==10.7",
-  "biopython==1.85",
+#  "requests",
+#  "openmm==8.2.0",
+#  "pdbfixer==1.11",
+  #"gemmi==0.7.1",
+#  "modeller==10.7",
+#  "biopython==1.85",
 ]
 readme = "README.md"
 license = "AGPL-3.0-or-later"
 
 [project.scripts]
 prepmd = "prepmd.prep:entry_point"
 runmd = "prepmd.run:entry_point"
+prep-license = "prepmd.add_modeller_license:entry_point"