Skip to content

Commit

Permalink
i have changes here?
Browse files Browse the repository at this point in the history
  • Loading branch information
corinwagen committed Oct 5, 2020
1 parent e452730 commit e8811b1
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 29 deletions.
7 changes: 4 additions & 3 deletions cctk/molecule.py
Original file line number Diff line number Diff line change
Expand Up @@ -1648,19 +1648,20 @@ def optimize(self, inplace=True, nprocs=1):
else:
return optimized

def csearch(self, nprocs=1, constraints=[], logfile=None, noncovalent=False):
def csearch(self, nprocs=1, constraints=[], logfile=None, noncovalent=False, use_tempdir=True):
"""
Optimize molecule at the GFN2-xtb level of theory.
Args:
nprocs (int): number of processors to use
constraints (list): atoms numbers to freeze
noncovalent (Bool): whether or not to use non-covalent settings
noncovalent (bool): whether or not to use non-covalent settings
logfile (str): file to write ongoing ``crest`` output to
use_tempdir (bool): write intermediate files to hidden directory (as opposed to current directory)
Returns
ConformationalEnsemble
"""
import cctk.optimize as opt
assert isinstance(nprocs, int), "nprocs must be int!"
return opt.csearch(self, nprocs=nprocs, constraints=constraints, noncovalent=noncovalent, logfile=logfile)
return opt.csearch(molecule=self, nprocs=nprocs, constraints=constraints, noncovalent=noncovalent, logfile=logfile, use_tempdir=use_tempdir)
61 changes: 35 additions & 26 deletions cctk/optimize.py
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ def run_xtb(molecule, nprocs=1, return_energy=False):
else:
return output_mol

def csearch(molecule, constraints=None, rotamers=False, nprocs=1, noncovalent=False, logfile=None):
def csearch(use_tempdir=True, **kwargs):
"""
Run a conformational search on a molecule using ``crest``.
Expand All @@ -93,45 +93,54 @@ def csearch(molecule, constraints=None, rotamers=False, nprocs=1, noncovalent=Fa
Returns:
cctk.ConformationalEnsemble
"""
assert installed("crest"), "crest must be installed!"

ensemble = None
try:
if use_tempdir:
with tempfile.TemporaryDirectory() as tmpdir:
ensemble = _do_csearch(directory=tmpdir, **kwargs)
else:
ensemble = _do_csearch(directory=os.getcwd(), **kwargs)
except Exception as e:
raise ValueError(f"Error running xtb:\n{e}")

return ensemble

def _do_csearch(molecule, nprocs, logfile, noncovalent, directory, constraints):
assert isinstance(molecule, cctk.Molecule), "need a valid molecule!"
assert isinstance(nprocs, int)
assert isinstance(logfile, str)

assert installed("crest"), "crest must be installed!"
cctk.XYZFile.write_molecule_to_file(f"{directory}/xtb-in.xyz", molecule)

nci = ""
if noncovalent:
nci = "-nci"

ensemble = None
try:
with tempfile.TemporaryDirectory() as tmpdir:
cctk.XYZFile.write_molecule_to_file(f"{tmpdir}/xtb-in.xyz", molecule)
if constraints is not None:
assert isinstance(constraints, list)
assert all(isinstance(n, int) for n in constraints)
command = f"crest xtb-in.xyz --constrain {','.join([str(c) for c in constraints])}"
result = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, cwd=directory, shell=True)
result.check_returncode()

if constraints is not None:
assert isinstance(constraints, list)
assert all(isinstance(n, int) for n in constraints)
command = f"crest xtb-in.xyz --constrain {','.join([str(c) for c in constraints])}"
result = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, cwd=tmpdir, shell=True)
result.check_returncode()
command = f"crest xtb-in.xyz --chrg {molecule.charge} --uhf {molecule.multiplicity - 1} -T {nprocs} {nci}"

command = f"crest xtb-in.xyz --chrg {molecule.charge} --uhf {molecule.multiplicity - 1} -T {nprocs} {nci}"
if logfile:
with open(logfile, "w") as f:
result = sp.run(command, stdout=f, stderr=f, cwd=directory, shell=True)
else:
result = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, cwd=directory, shell=True)

if logfile:
with open(logfile, "w") as f:
result = sp.run(command, stdout=f, stderr=f, cwd=tmpdir, shell=True)
else:
result = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, cwd=tmpdir, shell=True)
result.check_returncode()

result.check_returncode()
if rotamers:
ensemble = cctk.XYZFile.read_ensemble(f"{directory}/crest_rotamers.xyz")
else:
ensemble = cctk.XYZFile.read_ensemble(f"{directory}/crest_conformers.xyz")

if rotamers:
ensemble = cctk.XYZFile.read_ensemble(f"{tmpdir}/crest_rotamers.xyz")
else:
ensemble = cctk.XYZFile.read_ensemble(f"{tmpdir}/crest_conformers.xyz")
return ensemble

except Exception as e:
raise ValueError(f"Error running xtb:\n{e}")

return ensemble

0 comments on commit e8811b1

Please sign in to comment.