added crest optimize

cctk/molecule.py

@@ -1647,3 +1647,15 @@ def optimize(self, inplace=True, nprocs=1):
             return self
         else:
             return optimized
+
+    def csearch(self, nprocs=1):
+        """
+        Optimize molecule at the GFN2-xtb level of theory.
+
+        Args:
+            inplace (Bool): whether or not to return a new molecule or simply modify ``self.geometry``
+            nprocs (int): number of processors to use
+        """
+        import cctk.optimize as opt
+        assert isinstance(nprocs, int), "nprocs must be int!"
+        return opt.csearch(self, nprocs=nprocs)

cctk/optimize.py

@@ -70,3 +70,48 @@ def run_xtb(molecule, nprocs=1, return_energy=False):
     else:
         return output_mol
 
+def csearch(molecule, constrain_atoms=None, rotamers=True, nprocs=1):
+    """
+    Run a conformational search on a molecule using ``crest``.
+
+    Args:
+        molecule (cctk.Molecule): molecule of interest
+        constraints (list): list of atom numbers to constrain
+        rotamers (bool): return all rotamers or group into distinct conformers
+        nprocs (int): number of processors to use
+
+    Returns:
+        cctk.ConformationalEnsemble
+    """
+
+    assert isinstance(molecule, cctk.Molecule), "need a valid molecule!"
+    assert isinstance(nprocs, int)
+
+    assert shutil.which("xtb") is not None, "xtb must be installed for csearch to work!"
+    assert shutil.which("crest") is not None, "crest must be installed for csearch to work!"
+
+    ensemble = None
+    try:
+        with tempfile.TemporaryDirectory() as tmpdir:
+            if constrain_atoms is not None:
+                assert isinstance(constrain_atoms, list)
+                assert all(isinstance(n, int) for n in constrain_atoms)
+                command = f"crest --constrain {','.join(constrain_atoms)}"
+                result = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, cwd=tmpdir, shell=True)
+                result.check_returncode()
+
+            cctk.XYZFile.write_molecule_to_file(f"{tmpdir}/xtb-in.xyz", molecule)
+            command = f"crest xtb-in.xyz --chrg {molecule.charge} --uhf {molecule.multiplicity - 1} -T {nprocs}"
+            result = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, cwd=tmpdir, shell=True)
+            result.check_returncode()
+
+            if rotamers:
+                ensemble = cctk.XYZFile.read_ensemble(f"{tmpdir}/crest_rotamers.xyz")
+            else:
+                ensemble = cctk.XYZFile.read_ensemble(f"{tmpdir}/crest_conformers.xyz")
+
+    except Exception as e:
+        raise ValueError(f"Error running xtb:\n{e}")
+
+    return ensemble
+

cctk/xyz_file.py

@@ -114,3 +114,24 @@ def read_trajectory(cls, filename):
             current_lines.append(line)
 
         return files
+
+    @classmethod
+    def read_ensemble(cls, filename, conformational=False):
+        """
+        Alias for read_trajectory.
+        """
+        files = cls.read_trajectory(filename)
+
+        ensemble = None
+        if conformational:
+            ensemble = cctk.ConformationalEnsemble()
+        else:
+            ensemble = cctk.Ensemble()
+
+        for file in files:
+            ensemble.add_molecule(file.molecule)
+
+        return ensemble
+
+
+

test/static/adamantane.mol2

@@ -6,7 +6,7 @@
 #
 
 @<TRIPOS>MOLECULE
-Molecule Name
+title
 26 28
 SMALL
 NO_CHARGES

test/test_core.py

@@ -37,6 +37,11 @@ def test_traj(self):
         self.assertTrue(isinstance(files[0], cctk.XYZFile))
         self.assertEqual(files[0].molecule.num_atoms(), 5)
 
+    def test_ense(self):
+        path = "test/static/methane_traj.xyz"
+        ense = cctk.XYZFile.read_trajectory(path)
+        self.assertEqual(len(ense), 250)
+
 class TestMAE(unittest.TestCase):
     def test_read(self):
         path = "test/static/dodecane_csearch-out.mae"

test/test_opt.py

@@ -1,4 +1,4 @@
-import unittest, sys, os, io, copy, math
+import unittest, sys, os, io, copy, math, shutil
 import numpy as np
 
 import cctk
@@ -18,5 +18,14 @@ def test_basic(self):
         for x1, x2 in zip(np.ravel(mol2.geometry), np.ravel(mol.geometry)):
             self.assertTrue(abs(float(x1)-float(x2)) < 0.1)
 
+    def test_csearch(self):
+        mol = cctk.GaussianFile.read_file("test/static/L-Ala.gjf").get_molecule()
+
+        if shutil.which("crest") is not None:
+            conformers = mol.csearch()
+        else:
+            pass
+
+
 if __name__ == '__main__':
     unittest.main()