From 0526d9368c08e69025d6cd4df98d2918411c8aee Mon Sep 17 00:00:00 2001 From: fangning-ren Date: Tue, 2 Jul 2024 00:12:09 -0400 Subject: [PATCH] enabled TIP3P mixed solvent --- autosolvate/autosolvate.py | 4 +- autosolvate/data/ch3cn/ch3cn.xyz | 2 +- autosolvate/dockers/_tleap_docker.py | 13 +++ autosolvate/molecule/molecule.py | 83 ++++++++++++++-- autosolvate/molecule/solventbox.py | 83 +++++++++------- autosolvate/multicomponent.py | 82 +++++++++------- autosolvate/utils/tools.py | 137 ++++++++++++++++++++------- test_water_acn.ipynb | 85 +++++++++++++++++ 8 files changed, 375 insertions(+), 114 deletions(-) create mode 100644 test_water_acn.ipynb diff --git a/autosolvate/autosolvate.py b/autosolvate/autosolvate.py index 6becb26..51d4acf 100644 --- a/autosolvate/autosolvate.py +++ b/autosolvate/autosolvate.py @@ -136,9 +136,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] diff --git a/autosolvate/data/ch3cn/ch3cn.xyz b/autosolvate/data/ch3cn/ch3cn.xyz index 5caea1b..301b256 100644 --- a/autosolvate/data/ch3cn/ch3cn.xyz +++ b/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 diff --git a/autosolvate/dockers/_tleap_docker.py b/autosolvate/dockers/_tleap_docker.py index 018fa46..b6ac2bf 100644 --- a/autosolvate/dockers/_tleap_docker.py +++ b/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: + 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"): diff --git a/autosolvate/molecule/molecule.py b/autosolvate/molecule/molecule.py index c53ae5c..b1af2b7 100644 --- a/autosolvate/molecule/molecule.py +++ b/autosolvate/molecule/molecule.py @@ -225,6 +225,7 @@ def __init__( name = "", residue_name = "MOL", folder = WORKING_DIR, + amber_solvent = False, ) -> None: """ The data class for holding all files of a molecule. @@ -243,19 +244,24 @@ def __init__( The residue name of the molecule, by default "MOL". 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__() + + if not self.amber_solvent: + self.name = process_system_name(name, xyzfile, support_input_format=Molecule._SUPPORT_INPUT_FORMATS) + self.read_coordinate(xyzfile) def read_coordinate(self, fname:str): ext = os.path.splitext(fname)[-1][1:] @@ -263,12 +269,77 @@ def read_coordinate(self, fname:str): 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) + newpath = self.reference_name + "." + e + subprocess.run(f"obabel -i {ext} {fname} -o {e} -O {newpath} ---errorlevel 0", shell = True) + setattr(self, e, newpath) + + def generate_pdb(self): + if not self.check_exist("pdb") and self.name != "water": + prep2pdb4amber_solvent(self) + self.logger.info(f"AMBER predefined solvent {self.name} is used.") + self.logger.info(f"Converted the predefined prep file {self.prep} to pdb file {self.pdb}") + elif self.name == "water": + 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") + + 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 \ No newline at end of file diff --git a/autosolvate/molecule/solventbox.py b/autosolvate/molecule/solventbox.py index f7f3b69..3690ffe 100644 --- a/autosolvate/molecule/solventbox.py +++ b/autosolvate/molecule/solventbox.py @@ -48,12 +48,13 @@ def __init__(self, """ self.name = process_system_name(name, solventbox, support_input_format=SolventBox._SUPPORT_INPUT_FORMATS, check_exist = False if amber_solvent else True) self.folder = os.path.abspath(folder) - self.frcmod = os.path.abspath(frcmod) self.box_name = box_name self.amber_solvent = amber_solvent - self.solventbox = os.path.abspath(solventbox) if not amber_solvent: + self.frcmod = os.path.abspath(frcmod) + self.solventbox = os.path.abspath(solventbox) + if not os.path.exists(solventbox): logger.error("The pre-built solvent box file {} does not exist".format(solventbox)) raise FileNotFoundError("The solvent box file does not exist") @@ -65,7 +66,9 @@ def __init__(self, self.box_name = self.get_box_name(self.solventbox) self.solventbox = getattr(self, ext) else: - pass + self.solventbox = solventbox + self.frcmod = frcmod + super(SolventBox, self).__init__(name = self.name) self.logger.name = self.__class__.__name__ @@ -86,38 +89,46 @@ def get_box_name(self, fname:str): break return previousname -AMBER_WATER = SolventBox( name='water', - solventbox='solvents.lib', - box_name='TIP3PBOX', - amber_solvent = True - ) - -AMBER_METHANOL = SolventBox( name='methanol', - solventbox='solvents.lib', - frcmod="frcmod.meoh", - box_name='MEOHBOX', - amber_solvent = True - ) - -AMBER_CHLOROFORM = SolventBox(name='chloroform', - solventbox='solvents.lib', - frcmod="frcmod.chcl3", - box_name='CHCL3BOX', - amber_solvent = True - ) - -AMBER_NMA = SolventBox(name='nma', - solventbox='solvents.lib', - frcmod="frcmod.nma", - box_name='NMABOX', - amber_solvent = True - ) + def update(self): + # override the original update method for default system as some parameter for solvent box are predefined then cannot be updated + if self.amber_solvent: + return + super().update() -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, +# available name for pre-equilibrated water box includes: +# POL3BOX, QSPCFWBOX, SPCBOX, SPCFWBOX, TIP3PBOX, TIP3PFBOX, TIP4PBOX, TIP4PEWBOX, OPCBOX, OPC3BOX, TIP5PBOX +AMBER_WATER_BOX = SolventBox(name='water', + solventbox='solvents.lib', + box_name='TIP3PBOX', + amber_solvent=True + ) + +AMBER_METHANOL_BOX = SolventBox(name='methanol', + solventbox='solvents.lib', + frcmod="frcmod.meoh", + box_name='MEOHBOX', + amber_solvent=True + ) + +AMBER_CHLOROFORM_BOX = SolventBox(name='chloroform', + solventbox='solvents.lib', + frcmod="frcmod.chcl3", + box_name='CHCL3BOX', + amber_solvent=True + ) + +AMBER_NMA_BOX = SolventBox(name='nma', + solventbox='solvents.lib', + frcmod="frcmod.nma", + box_name='NMABOX', + amber_solvent=True + ) + +AMBER_SOLVENTBOX_LIST = [AMBER_WATER_BOX, AMBER_METHANOL_BOX, AMBER_CHLOROFORM_BOX, AMBER_NMA_BOX] +AMBER_SOLVENTBOX_DICT = { + AMBER_WATER_BOX.name : AMBER_WATER_BOX, + AMBER_METHANOL_BOX.name : AMBER_METHANOL_BOX, + AMBER_CHLOROFORM_BOX.name : AMBER_CHLOROFORM_BOX, + AMBER_NMA_BOX.name : AMBER_NMA_BOX, } -logging.basicConfig(level = INFO, force = True, handlers=[]) \ No newline at end of file +logging.basicConfig(level=INFO, force=True, handlers=[]) \ No newline at end of file diff --git a/autosolvate/multicomponent.py b/autosolvate/multicomponent.py index d027f52..6909023 100644 --- a/autosolvate/multicomponent.py +++ b/autosolvate/multicomponent.py @@ -121,8 +121,7 @@ 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: - # amber solvents + if solvent in AMBER_SOLVENT_DICT and not slv_generate: self_solvent = AMBER_SOLVENT_DICT[solvent] elif solvent in custom_solv_dict: # solvent data prepared by autosolvate @@ -137,10 +136,7 @@ def get_solvent(self, solvent:str, slv_xyz:str = "", solvent_frcmod:str = "", so self_solvent.frcmod = solvent_frcmod_path self_solvent.prep = solvent_prep_path self_solvent.number = slv_count - elif os.path.exists(solvent_frcmod) and os.path.exists(solvent_off): - # using prebuilt solvent box - self_solvent = SolventBox(solvent_off, solvent_frcmod, name = solvent, folder = self.folder, box_name=solvent_box_name) - elif os.path.exists(slv_xyz) and slv_generate == True: + elif os.path.exists(slv_xyz) and slv_generate: # generate solvent box self_solvent = Molecule(slv_xyz, folder = self.folder) self_solvent.number = slv_count @@ -188,10 +184,9 @@ def __init__(self, folder = WORKING_DIR, cube_size = 54, closeness = 2.0, charge def add_solute(self, xyzfile:str, name="", residue_name="", charge=0, spinmult=1, number = 1, **kwargs): - #use the first three letters of the xyzfile name as the residue name - if not residue_name and len(xyzfile.split(".")[0].split("/")[-1]) >= 3: - residue_name = xyzfile.split(".")[0].split("/")[-1][:3].upper() - + # use the first three letters of the xyzfile name as the residue name if not provided + if not residue_name: + residue_name = try_ones_best_to_get_residue_name(xyzfile, name) molecule = Molecule(xyzfile, charge=charge, multiplicity=spinmult, folder = self.folder, name = name, residue_name=residue_name) if "mol2" in kwargs and os.path.isfile(kwargs["mol2"]): @@ -202,36 +197,54 @@ def add_solute(self, xyzfile:str, name="", residue_name="", charge=0, spinmult=1 molecule.lib = kwargs["lib"] molecule.update() - if not molecule.check_exist("frcmod") or not (molecule.check_exist("mol2") and not molecule.check_exist("lib")): + if not molecule.check_exist("frcmod") or not (molecule.check_exist("mol2") or molecule.check_exist("lib")): for docker in self.single_molecule_pipeline: docker.run(molecule) molecule.number = number self.solutes.append(molecule) - - def add_solvent(self, xyzfile:str, name="", residue_name="", charge=0, spinmult=1, number = 200, **kwargs): - - #use the first three letters of the xyzfile name as the residue name - if not residue_name and len(xyzfile.split(".")[0].split("/")[-1]) >= 3: - residue_name = xyzfile.split(".")[0].split("/")[-1][:3].upper() #get the file name from the absolute path - - molecule = Molecule(xyzfile, charge=charge, multiplicity=spinmult, folder = self.folder, name = name, residue_name=residue_name) - if "mol2" in kwargs and os.path.isfile(kwargs["mol2"]): - molecule.mol2 = kwargs["mol2"] - if "frcmod" in kwargs and os.path.isfile(kwargs["frcmod"]): - molecule.frcmod = kwargs["frcmod"] - if "lib" in kwargs and os.path.isfile(kwargs["lib"]): - molecule.lib = kwargs["lib"] - molecule.update() + def add_solvent(self, xyzfile:str = "", name="", residue_name="", charge=0, spinmult=1, number = 210*8, slv_generate = False, **kwargs): + if name in AMBER_SOLVENT_DICT and not slv_generate: # predefined amber solvent + self_solvent = AMBER_SOLVENT_DICT[name] + self_solvent.folder = self.folder + self_solvent.generate_pdb() + elif name in custom_solv_dict: # custom solvent in autosolvate + # solvent data prepared by autosolvate + solvPrefix = custom_solv_dict[name] + solvent_frcmod_path = pkg_resources.resource_filename('autosolvate', + os.path.join('data',solvPrefix,solvPrefix+".frcmod")) + solvent_prep_path = pkg_resources.resource_filename('autosolvate', + os.path.join('data',solvPrefix,solvPrefix+".prep")) + solvent_pdb_path = pkg_resources.resource_filename('autosolvate', + os.path.join('data',solvPrefix,solvPrefix+".pdb")) + self_solvent = Molecule(solvent_pdb_path, 0, 1, name, residue_name = custom_solv_residue_name[name], folder = self.folder) + self_solvent.frcmod = solvent_frcmod_path + self_solvent.prep = solvent_prep_path + elif os.path.exists(xyzfile) and slv_generate: # user defined solvent + # generate solvent box + residue_name = try_ones_best_to_get_residue_name(xyzfile, name) + self_solvent = Molecule(xyzfile, charge=charge, multiplicity=spinmult, folder = self.folder, name = name, residue_name=residue_name) + if "mol2" in kwargs and os.path.isfile(kwargs["mol2"]): + self_solvent.mol2 = kwargs["mol2"] + if "frcmod" in kwargs and os.path.isfile(kwargs["frcmod"]): + self_solvent.frcmod = kwargs["frcmod"] + if "lib" in kwargs and os.path.isfile(kwargs["lib"]): + self_solvent.lib = kwargs["lib"] + self_solvent.update() - if not molecule.check_exist("frcmod") or not (molecule.check_exist("mol2") and not molecule.check_exist("lib")): - for docker in self.single_molecule_pipeline: - docker.run(molecule) - molecule.number = number - self.solvents.append(molecule) + if not self_solvent.check_exist("frcmod") or not (self_solvent.check_exist("mol2") or self_solvent.check_exist("lib")): + for docker in self.single_molecule_pipeline: + docker.run(self_solvent) + else: + raise ValueError("Solvent not found") + self_solvent.number = number + for solvent in self.solvents: + if solvent.name == self_solvent.name: + raise ValueError(f"Solvent {solvent.name} already exists") + self.solvents.append(self_solvent) def build(self): - system_name = "-".join([m.name for m in self.solutes + self.solvents]) + system_name = "-".join([m.name for m in self.solutes] + [m.name for m in self.solvents]) solute_numbers = [m.number for m in self.solutes] solvent_numbers = [m.number for m in self.solvents] ''' @@ -239,7 +252,10 @@ def build(self): comment left by Patrick Jun 12 2024. the old way to define 'system_name' will cause bug, please have the person who wrote this to fix it. ''' - system_name = 'MYBOX' + ''' + @DEBUG + bug fixed by Fangning Ren on July 1 2024 + ''' system = SolvatedSystem(system_name, solute = self.solutes, solvent = self.solvents, cubesize=self.boxsize, closeness=self.closeness, solute_number = solute_numbers, solvent_number = solvent_numbers, diff --git a/autosolvate/utils/tools.py b/autosolvate/utils/tools.py index 90ed69f..1a88e80 100644 --- a/autosolvate/utils/tools.py +++ b/autosolvate/utils/tools.py @@ -11,6 +11,7 @@ from ..Common import * +# multicomponent utilities def check_multicomponent(filename:str): """ Check if the given file contains multiple molecules @@ -259,6 +260,8 @@ def write_xyz(fname, element, data, energy = ""): f.write(f"{element[i]} {data[i][0]:>12.8f} {data[i][1]:>12.8f} {data[i][2]:>12.8f}\n") f.close() + +# processing PDB def getHeadTail(mol2): r""" Detect start and end of coordinates @@ -306,6 +309,92 @@ def updatePDB(obmol:ob.OBMol, pdbpath:str): pyobmol = pybel.Molecule(ob.OBMol(obmol)) pyobmol.write("pdb", pdbpath, overwrite=True) +def xyz_to_pdb(mol: object) -> None: + r''' + Convert xyz file to pdb file using openbabel + ''' + # os.system('ob -ixyz {} -opdb -O {}'.format(mol.xyz, mol.name+'.pdb')) + obConversion = ob.OBConversion() + obConversion.SetInAndOutFormats("xyz", "pdb") + OBMOL = ob.OBMol() + obConversion.ReadFile(OBMOL, mol.xyz) + obConversion.WriteFile(OBMOL, mol.name+'/'+mol.name+'.pdb') + mol.update() + edit_pdb(mol) + +def edit_pdb(mol: object) -> None: + with open(mol.name+'/'+mol.pdb, 'r') as f: + lines = f.readlines() + with open(mol.name+'/'+mol.pdb, 'w') as f: + for line in lines: + if line.startswith('CONECT'): + continue + else: + f.write(line) + +def edit_system_pdb(box: object) -> None: + with open(box.name+'/'+box.system_pdb, 'r') as f: + lines = f.readlines() + with open(box.name+'/'+box.system_pdb, 'w') as f: + this_resid = 1 + last_resid = 1 + for line in lines: + if 'ATOM' in line: + last_resid = int(this_resid) + this_resid = int(line[22:26]) + if last_resid != this_resid: + f.write("TER\n") + f.write(line) + f.close() + +def prep2pdb4amber_solvent(mol: object) -> None: + r''' + Output a pdb file for a amber predefined solvent (methanol, chloroform, nma) + ''' + outputfolder = mol.folder + resname = mol.residue_name + pdbname = mol.name + '.pdb' + pdbpath = mol.reference_name + '.pdb' + if not os.path.exists(outputfolder): + os.makedirs(outputfolder) + if os.path.exists(pdbpath) and os.path.isfile(pdbpath): + os.remove(pdbpath) + with open(os.path.join(mol.folder, "leap_convert.cmd"), "w") as f: + f.write(f"loadAmberPrep {mol.prep}\n") + f.write(f"singlemol = combine { {mol.residue_name} }\n") + f.write(f"savepdb singlemol {pdbpath}\n") + os.system(f"tleap -f {os.path.join(mol.folder, 'leap_convert.cmd')}") + mol.pdb = pdbpath + +def assign_water_pdb(mol: object) -> None: + '''assign a reference pdb file for water''' + pdbstr = ''' +ATOM 1 O WAT 0 2.537 -0.155 0.000 0.00 0.00 O +ATOM 2 H1 WAT 0 3.074 0.155 0.000 0.00 0.00 H +ATOM 3 H2 WAT 0 2.000 0.155 0.000 0.00 0.00 H +CONECT 1 2 3 +CONECT 2 1 +CONECT 3 1 +END +'''.strip("\n") + outputfolder = mol.folder + if not os.path.exists(outputfolder): + os.makedirs(outputfolder) + pdbpath = mol.reference_name + '.pdb' + with open(pdbpath, 'w') as f: + f.write(pdbstr) + mol.pdb = pdbpath + +def get_residue_name_from_pdb(file_path:str) -> str: + residue_name = None + with open(file_path, 'r') as file: + for line in file: + if line.startswith('ATOM') or line.startswith('HETATM'): + residue_name = line[17:20] + break + return residue_name + + # other functions def process_system_name(name:str, xyzfile:str, support_input_format:Iterable[str] = ("xyz", "pdb", "mol2", "prep", "off", "lib"), check_exist = True): if not os.path.isfile(xyzfile) and check_exist: @@ -374,44 +463,20 @@ def get_list_mol_type(*args: object, mol_type: str) -> list: raise Warning('mol_type is not set') return solvent_list -def xyz_to_pdb(mol: object) -> None: - r''' - Convert xyz file to pdb file using openbabel - ''' - # os.system('ob -ixyz {} -opdb -O {}'.format(mol.xyz, mol.name+'.pdb')) - obConversion = ob.OBConversion() - obConversion.SetInAndOutFormats("xyz", "pdb") - OBMOL = ob.OBMol() - obConversion.ReadFile(OBMOL, mol.xyz) - obConversion.WriteFile(OBMOL, mol.name+'/'+mol.name+'.pdb') - mol.update() - edit_pdb(mol) +def try_ones_best_to_get_residue_name(self, xyzfile:str, name:str): + try: + assert xyzfile.endswith(".pdb"), "Only pdb file is supported for now" + residue_name = get_residue_name_from_pdb(xyzfile) + assert len(residue_name) == 3, "Residue name must be three letters long" + except: + xyzbasename = os.path.splitext(os.path.basename(xyzfile))[0] if not name else name + if len(xyzbasename) >= 3: + residue_name = xyzbasename[:3].upper() + elif len(xyzbasename) < 3: + residue_name = "SLU" + return residue_name -def edit_pdb(mol: object) -> None: - with open(mol.name+'/'+mol.pdb, 'r') as f: - lines = f.readlines() - with open(mol.name+'/'+mol.pdb, 'w') as f: - for line in lines: - if line.startswith('CONECT'): - continue - else: - f.write(line) -def edit_system_pdb(box: object) -> None: - with open(box.name+'/'+box.system_pdb, 'r') as f: - lines = f.readlines() - with open(box.name+'/'+box.system_pdb, 'w') as f: - this_resid = 1 - last_resid = 1 - for line in lines: - if 'ATOM' in line: - last_resid = int(this_resid) - this_resid = int(line[22:26]) - if last_resid != this_resid: - f.write("TER\n") - f.write(line) - f.close() - #handle job submission in class def srun() -> callable: r''' diff --git a/test_water_acn.ipynb b/test_water_acn.ipynb new file mode 100644 index 0000000..c09de38 --- /dev/null +++ b/test_water_acn.ipynb @@ -0,0 +1,85 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "from autosolvate.autosolvate import *\n", + "from autosolvate.multicomponent import *\n", + "import os \n", + "ppath = \"/home/fren5/AutoSolvate-update/AutoSolvate\" \n", + "os.chdir(ppath)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "1 molecule converted\n", + "1 molecule converted\n" + ] + } + ], + "source": [ + "neapath = \"/home/fren5/AutoSolvate-update/AutoSolvate/mixture-test\" \n", + "if not os.path.exists(neapath):\n", + " os.makedirs(neapath)\n", + "os.chdir(neapath)\n", + "\n", + "builder = MixtureBuilder(folder = neapath, cube_size=30, closeness=1.3)\n", + "builder.add_solvent(\"/home/fren5/AutoSolvate-update/AutoSolvate/docs/_data/multicomponent_tutorial/example1/tutorial_step1/acetonitrile.pdb\", \n", + " name = \"acetonitrile\",\n", + " residue_name = \"C3N\",\n", + " charge = 0, \n", + " spinmult = 1, \n", + " number = 200, \n", + ")\n", + "builder.add_solvent(\"/home/fren5/AutoSolvate-update/AutoSolvate/docs/_data/multicomponent_tutorial/example1/tutorial_step1/water.pdb\", \n", + " name = \"water\",\n", + " residue_name = \"WAT\",\n", + " charge = 0,\n", + " spinmult = 1,\n", + " number = 600,\n", + ")\n", + "builder.add_solute(\"/home/fren5/AutoSolvate-update/AutoSolvate/docs/_data/multicomponent_tutorial/example1/tutorial_step1/naphthalene_neutral.xyz\", \n", + " name = \"napthalene\",\n", + " residue_name = \"NAP\",\n", + " charge = 0,\n", + " spinmult = 1,\n", + " number = 1,\n", + ")\n", + "builder.build()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.18" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}