update multicomponent tutorial

autosolvate/__main__.py

@@ -27,8 +27,9 @@ def main(args=None):
         print('Usage: autosolvate [OPTIONS]')
         print('  [NO OPTION]                launches GUI')
         print('  boxgen [OPTIONS]           generate initial structure')
-        print('  boxgen_metal[OPTIONS]      generate initial structure for organometallic compounds')
-        print('  mdrun [OPTIONS]            automated QM/MM trajectory generatio')
+        print('  boxgen_metal[OPTIONS]      generate initial structure and forcefield for organometallic compounds')
+        print('  boxgen_multicomponent [OPTIONS] generate initial structure for multicomponent systems')
+        print('  mdrun [OPTIONS]            automated QM/MM trajectory generation')
         print('  clustergen [OPTIONS]       extract microsolvated clusters ')
         print('  -h, --help                 short usage description')
         print()

autosolvate/autosolvate.py

@@ -43,6 +43,7 @@ class AmberParamsBuilder(object):
     1. Generate standard pdb
     2. AnteChamber or Gaussian charge fitting
     3. Tleap create Lib
+    Others work with similar function: acpype
     """
     def __init__(self, xyzfile:str, name = "", resname = "", charge = 0, spinmult = 1, 
                  charge_method="resp", folder = WORKING_DIR, **kwargs):
@@ -136,9 +137,9 @@ def __init__(self,
         ]
 
     def get_solvent(self, solvent:str, slv_xyz:str = "", solvent_frcmod:str = "", solvent_off:str = "", slv_generate:bool = False, slv_count:int = 210*8, solvent_box_name:str = "SLVBOX"):
-        if solvent in AMBER_SOLVENT_DICT:
+        if solvent in AMBER_SOLVENTBOX_DICT:
             # amber solvents
-            self_solvent = AMBER_SOLVENT_DICT[solvent]
+            self_solvent = AMBER_SOLVENTBOX_DICT[solvent]
         elif solvent in custom_solv_dict:
             # solvent data prepared by autosolvate
             solvPrefix = custom_solv_dict[solvent]

autosolvate/data/ch3cn/ch3cn.xyz

@@ -1,5 +1,5 @@
 6
-/home/fren5/AutoSolvae-update/AutoSolvate/autosolvate/data/ch3cn/ch3cn.pdb
+/home/fren5/AutoSolvate-update/AutoSolvate/autosolvate/data/ch3cn/ch3cn.pdb
 N         14.18200       11.28800       13.72000
 C         13.67300       12.18800       13.29900
 C         13.15300       13.38300       12.69800

autosolvate/dockers/_antechamber_docker.py

@@ -24,7 +24,7 @@ def __init__(self,
                  charge_method:         str = 'bcc', 
                  out_format:            str = 'mol2',
                  workfolder:            str = WORKING_DIR,
-                 exeoutfile:            str = None,
+                 exeoutfile:            str = "antechamber.log",
                  eq:                    int = 2,
                  pl:                    int = -1,
 

autosolvate/dockers/_general_docker.py

@@ -52,7 +52,7 @@ def execute(self, cmd):
             self.logger.info("CMD: {}".format(cmd))
             subprocess.run(cmd, shell = True, stdout=exeout, stderr=sys.stdout)
         else:
-            exeout = open(self.exeoutfile, "w")
+            exeout = open(self.exeoutfile, "a")
             self.logger.info("CMD: {}".format(cmd))
             subprocess.run(cmd, shell = True, stdout=exeout, stderr=sys.stdout)
             exeout.close()

autosolvate/dockers/_packmol_docker.py

@@ -112,6 +112,7 @@ def load_solute(self, doc: TextIO, box: SolvatedSystem) -> None:
             # doc.write('{:<15} {}                                 \n'.format('structure', solute.pdb))
             if not os.path.exists(os.path.join(self.workfolder, os.path.basename(solute.pdb))):
                 shutil.copy(solute.pdb, os.path.join(self.workfolder, os.path.basename(solute.pdb)))
+                # this is actually a bug in packmol. The path to the pdb file cannot be too long. One have to copy the file to the working directory
             doc.write('{:<15} {}                                 \n'.format('structure', os.path.basename(solute.pdb)))
             doc.write('{:<15} {:<5}                              \n'.format('number',    1))
             doc.write('{:<10} {posx} {posy} {posz} {com} {com} {com}\n'.format('fixed', posx=solute_pos[0], posy=solute_pos[1], posz=solute_pos[2], com=0.0))

autosolvate/dockers/_tleap_docker.py

@@ -40,6 +40,15 @@ def __init__(self,
     ##### check_system method for different systems
     def check_system_molecule(self, mol:Molecule):
         self.logger.info("Checking system {:s}...".format(mol.name))
+        if mol.amber_solvent:
+            self.logger.info("This is a predefined AMBER solvent.")
+            if mol.name == "water":
+                self.logger.info(f"Using TIP3P water model with predefined residue name {mol.residue_name}")
+            else:
+                self.logger.info("Solvent prep file: {:s}".format(mol.prep))
+                self.logger.info("Solvent frcmod file: {:s}".format(mol.frcmod))
+            return 
+
         suffixs = ["mol2", "lib", "prep", "off"]
         suffixf = [0, 0, 0, 0]
         for i, suffix in enumerate(suffixs):
@@ -199,6 +208,10 @@ def load_solvent_box(self,
                         mol:           SolventBox,
                         check:         bool = False
     ) -> None:
+        if mol.amber_solvent and not mol.name == "water":
+            doc.write('{:<5} {}  \n'.format("loadamberparams", mol.frcmod))
+            return 
+
         if mol.check_exist("off"):
             doc.write('{:<5} {}  \n'.format("loadoff", mol.off))
         if mol.check_exist("lib"):

autosolvate/molecule/molecule.py

@@ -92,7 +92,8 @@ def __setattr__(self, __name: str, __value: Any) -> None:
             else:
                 if __name in self._FILEATTR_FILEPATH_DICT:
                     self._FILEATTR_FILEPATH_DICT[__name] = __value
-                    self.logger.warning("The '{:s}' file of system '{:s}' is set to a non-existent path '{:s}'".format(__name, self.name, __value))
+                    if not (hasattr(self, "amber_solvent") and self.amber_solvent):
+                        self.logger.warning("The '{:s}' file of system '{:s}' is set to a non-existent path '{:s}'".format(__name, self.name, __value))
                 else:
                     self.logger.debug("set the '{:s}' attribute of system '{:s}' as {:s}".format(__name, self.name, __value))
         else:
@@ -120,7 +121,7 @@ def check_exist(self, extensions:List[str]) -> bool:
                 self.logger.debug("The '{:s}' file of system '{:s}' is not defined".format(ext, self.name))
             elif not isinstance(self.__getattribute__(ext), str):
                 flag *= False
-                self.logger.debug("The '{:s}' attribute of system '{:s}' does noe exist".format(ext, self.name))
+                self.logger.debug("The '{:s}' attribute of system '{:s}' does not exist".format(ext, self.name))
             elif not os.path.isfile(self.__getattribute__(ext)):
                 flag *= False
                 self.logger.debug("The '{:s}' file of system '{:s}' does not exist".format(ext, self.name))
@@ -215,60 +216,195 @@ def reference_name(self) -> str:
 # @dataclass 
 class Molecule(System): 
     #constants 
-    _SUPPORT_INPUT_FORMATS = ['pdb', 'xyz'] 
+    _SUPPORT_INPUT_FORMATS = ['pdb', 'xyz', 'mol2', 'prep']
     # other
     def __init__(
             self, 
             xyzfile:        str, 
             charge:         int, 
             multiplicity:   int,
             name            = "",
-            residue_name    = "MOL",
+            residue_name    = "",
             folder          = WORKING_DIR,
+            amber_solvent   = False,
             ) -> None:
         """
         The data class for holding all files of a molecule.
 
         Parameters
         ----------
         xyzfile : str
-            The file path of the xyz file of the molecule, can also be a pdb file.
+            The file path of the structure file of the molecule, can be in xyz, pdb, mol2 or prep format.
         charge : int
             The charge of the molecule.
         multiplicity : int
             The multiplicity of the molecule.
         name : str, optional
             The name of the molecule, by default "" will be the basename of the xyzfile.
         residue_name : str, optional
-            The residue name of the molecule, by default "MOL".
+            The residue name of the molecule, ignored if the xyzfile is provided as mol2 or prep format, by default will be assigned according to the filename.
         folder : str, optional
             The folder to store all files of the molecule, by default is the current working directory.
+        amber_solvent : bool, optional
+            A flag to indicate whether the molecule will use amber predefined parameters, by default False.
         """
 
-        self.name           = process_system_name(name, xyzfile, support_input_format=Molecule._SUPPORT_INPUT_FORMATS)
         self.folder         = os.path.abspath(folder)
         self.charge         = charge
         self.multiplicity   = multiplicity
         self.spinmult       = multiplicity
         self.residue_name   = residue_name
         self.number         = 0
-        self.read_coordinate(xyzfile)
+        self.amber_solvent = amber_solvent 
+        self.name           = name
         super(Molecule, self).__init__(name = self.name)
         self.logger.name = self.__class__.__name__
-        # super(Molecule, self).__post_init__()
 
-    def read_coordinate(self, fname:str):
-        ext = os.path.splitext(fname)[-1][1:]
-        setattr(self, ext, fname)
+        if not self.amber_solvent:
+            if not self.name:
+                self.name = process_system_name(name, xyzfile, support_input_format=Molecule._SUPPORT_INPUT_FORMATS, check_exist=False)  
+            self.process_input_xyz(xyzfile)
+            self.generate_pdb()
+
+    def process_input_xyz(self, xyzfile:str) -> None:
         for e in Molecule._SUPPORT_INPUT_FORMATS:
-            if e == ext:
-                continue
-            nname = os.path.splitext(fname)[0] + "." + e
-            subprocess.run(f"obabel -i {ext} {fname} -o {e} -O {nname} ---errorlevel 0", shell = True)
-            setattr(self, e, nname)
+            ext = os.path.splitext(xyzfile)[-1]
+            if e in ext:
+                self.__setattr__(e, xyzfile)
+                break
+        else:
+            raise ValueError(f"Unsupported input format for the molecule: {xyzfile}")
+
+    def generate_pdb(self):
+        """
+        Convert the input structure file to pdb file if the provided coordinate file is not in pdb format.
+        If the pdb file for the same molecule already exists, it will be ignored.
+        This behavior is to keep the residue name of the molecule consistent.
+        """
+        if self.name == "water" and self.amber_solvent:
+            assign_water_pdb(self)
+            self.logger.info(f"AMBER predefined water is used.")
+            self.logger.info(f"Write the reference pdb file for {self.name} to {self.pdb}")
+            self.logger.info(f"Water model used: TIP3P")
+        elif self.amber_solvent:
+            prep2pdb(self)
+            self.logger.info(f"Predefined solvent {self.name} is used.")
+            self.logger.info(f"Converted the predefined prep file {self.prep} to pdb file {self.pdb}")
+        elif self.check_exist("prep"):
+            newpath = self.reference_name + ".pdb"
+            if self.check_exist("pdb"):
+                self.logger.warning(f"The existing pdb file {self.pdb} will be ignored as amber prep file is provided.")
+                shutil.move(self.pdb, self.reference_name + "-bak.pdb")
+            self.logger.info(f"Converted the prep file {self.prep} to pdb file {self.pdb}")
+            prep2pdb(self)
+        elif self.check_exist("lib") or self.check_exist("off"):
+            newpath = self.reference_name + ".pdb"
+            if self.check_exist("pdb"):
+                self.logger.warning(f"The existing pdb file {self.pdb} will be ignored as amber lib/off file is provided.")
+                shutil.move(self.pdb, self.reference_name + "-bak.pdb")
+            self.logger.info(f"Converted the lib file {self.lib} to pdb file {self.pdb}")
+            lib2pdb(self)
+        elif self.check_exist("mol2"):
+            newpath = self.reference_name + ".pdb"
+            if self.check_exist("pdb"):
+                self.logger.warning(f"The existing pdb file {self.pdb} will be ignored as mol2 file is provided.")
+                shutil.move(self.pdb, self.reference_name + "-bak.pdb")
+            self.logger.info(f"Converted the mol2 file {self.mol2} to pdb file {self.pdb}")
+            subprocess.run(f"obabel -i mol2 {self.mol2} -o pdb -O {newpath} ---errorlevel 0", shell = True)
+            self.pdb = newpath
+        elif self.check_exist("pdb"):
+            pass
+        elif self.check_exist("xyz"):
+            newpath = self.reference_name + ".pdb"
+            subprocess.run(f"obabel -i xyz {self.xyz} -o pdb -O {newpath} ---errorlevel 0", shell = True)
+            self.pdb = newpath
+
+    def get_residue_name(self):
+        """
+        Get the residue name from the provided pdb, mol2, or prep file.
+        Note that the residue name of the molecule will be automatically assigned if the residue name is not found.
+        Note if the mol2, or prep file is provided, the 'residue_name' attribute will be updated according to these files instead of the argument 'residue_name'.
+        """
+        new_residue_name = ""
+        if self.check_exist("prep"):
+            new_residue_name = extract_residue_name_from_prep(self.prep)
+        elif self.check_exist("lib"):
+            new_residue_name = extract_residue_name_from_lib(self.lib)
+        elif self.check_exist("off"):
+            new_residue_name = extract_residue_name_from_lib(self.off)
+        elif self.check_exist("mol2"):
+            newpath = self.reference_name + "-frommol2.pdb"
+            subprocess.run(f"obabel -i mol2 {self.mol2} -o pdb -O {newpath} ---errorlevel 0", shell = True)
+            new_residue_name = get_residue_name_from_pdb(newpath)
+        elif self.residue_name:
+            new_residue_name = self.residue_name
+        else:
+            self.logger.warning(f"Cannot find the residue name of molecule {self.name}, trying to automatically assign one.")
+            new_residue_name = try_ones_best_to_get_residue_name(self, self.pdb, self.name)
+        if not new_residue_name:
+            new_residue_name = "MOL"
+        if self.residue_name != new_residue_name:
+            self.logger.warning(f"Residue name of {self.name} changed: {self.residue_name} -> {new_residue_name}")
+            self.logger.warning(f"This behavior is to keep the residue name of the molecule consistent.")
+            self.residue_name = new_residue_name
+        self.logger.info(f"Residue name of {self.name} is set to {self.residue_name}")
+
+    def update(self):
+        if not self.amber_solvent:
+            super(Molecule, self).update()
+
+
+
+# if one want to use other water forcefield, one may need to change the residue name of water to other, such as PL3 for POL3 model. 
+AMBER_WATER             = Molecule(name = "water",
+                                   xyzfile = None,
+                                   charge = 0,
+                                   multiplicity = 1,
+                                   residue_name = "WAT",
+                                   folder = WORKING_DIR,
+                                   amber_solvent = True, 
+                                  )
+
+AMBER_METHANOL          = Molecule(name = "methanol",
+                                   xyzfile = None, 
+                                   charge = 0,
+                                   multiplicity = 1,
+                                   residue_name = "MOH",
+                                   folder = WORKING_DIR,
+                                   amber_solvent = True, 
+                                  )
+AMBER_METHANOL.frcmod = "frcmod.meoh"
+AMBER_METHANOL.prep = "meoh.in"
 
+AMBER_CHLOROFORM        = Molecule(name = "chloroform",
+                                   xyzfile = None,
+                                   charge = 0,
+                                   multiplicity = 1,
+                                   residue_name = "CL3",
+                                   folder = WORKING_DIR,
+                                   amber_solvent = True, 
+                                  )
+AMBER_CHLOROFORM.frcmod = "frcmod.chcl3"
+AMBER_CHLOROFORM.prep = "chcl3.in"
 
+AMBER_NMA              = Molecule(name = "nma",
+                                   xyzfile = None,
+                                   charge = 0,
+                                   multiplicity = 1,
+                                   residue_name = "NMA",
+                                   folder = WORKING_DIR,
+                                   amber_solvent = True, 
+                                  )
+AMBER_NMA.frcmod = "frcmod.nma"
+AMBER_NMA.prep = "nma.in"
 
+AMBER_SOLVENT_LIST        = [AMBER_WATER, AMBER_METHANOL, AMBER_CHLOROFORM, AMBER_NMA]  
+AMBER_SOLVENT_DICT        = { 
+                            AMBER_WATER.name:       AMBER_WATER,
+                            AMBER_METHANOL.name:    AMBER_METHANOL,
+                            AMBER_CHLOROFORM.name:  AMBER_CHLOROFORM,
+                            AMBER_NMA.name:         AMBER_NMA,
+                            }
 
 if __name__ == "__main__":
     pass

autosolvate/molecule/molecule_complex.py

@@ -65,8 +65,6 @@ def __init__(self,
         self.residue_name   = "SYS" if not residue_name else residue_name
         self.number         = 0
         self.read_coordinate(xyzfile)
-        super(MoleculeComplex, self).__init__(name = self.name)
-        self.logger.name = self.__class__.__name__
 
         self.mol_obmol          = pybel.readfile("pdb", self.pdb).__next__().OBMol
         self.fragresiduenames   = []    # name of each fragment
@@ -82,19 +80,21 @@ def __init__(self,
 
         self.aminoacidresidues  = AMINO_ACID_RESIDUES
 
+        super(MoleculeComplex, self).__init__(name = self.name)
+        self.logger.name = self.__class__.__name__
+
         if reorder_pdb:
             reorderPDB(self.pdb, self.pdb)
         self.build_molecules()
 
     def read_coordinate(self, fname:str):
         ext = os.path.splitext(fname)[-1][1:]
         setattr(self, ext, fname)
-        for e in Molecule._SUPPORT_INPUT_FORMATS:
-            if e == ext:
-                continue
-            nname = os.path.splitext(fname)[0] + "." + e
-            subprocess.run(f"obabel -i {ext} {fname} -o {e} > {nname}", shell = True)
-            setattr(self, e, nname)
+        if ext != "pdb":
+            subprocess.run(f"obabel -i {ext} {fname} -o pdb -O {self.reference_name}.pdb ---errorlevel 0", shell = True)
+            # subprocess.run(f"obabel -i {ext} {fname} -o pdb -O {self.reference_name}.pdb", shell = True)
+            # os.system(f"obabel -i {ext} {fname} -o pdb -O {self.reference_name}.pdb")
+            self.pdb = f"{self.reference_name}.pdb"
 
     def check_protein_fragment(self, fragment_obmol:ob.OBMol, fragment_index = 0):
         for j in range(fragment_obmol.NumResidues()):
@@ -145,7 +145,7 @@ def getFragments(self):
             fragres = fragment_obmol.GetResidue(0)
             fragresname = fragres.GetName()
             # This is a water molecule or an amino acid moleucle.
-            if fragresname == "HOH" or fragresname in AMINO_ACID_RESIDUES:
+            if fragresname == "HOH" or fragresname == "WAT" or fragresname in AMINO_ACID_RESIDUES:
                 self.fragmols.append(ob.OBMol(fragment_obmol))
                 self.fragresiduenames.append(fragresname)
                 logger.info(f"Fragment {i} is a known molecule with res name {fragresname}. ")