add CC, CI, MP3 energies

cctk/molecule.py

@@ -1638,23 +1638,42 @@ def atomic_symbols(self):
         symbols = {z: get_symbol(z) for z in set(self.atomic_numbers)}
         return [symbols[z] for z in self.atomic_numbers]
 
-    def optimize(self, inplace=True, nprocs=1):
+    def optimize(self, inplace=True, nprocs=1, return_energy=False):
         """
         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
+            return_energy (Bool): whether to return energy or not
         """
         import cctk.optimize as opt
         assert isinstance(nprocs, int), "nprocs must be int!"
-        optimized = opt.optimize_molecule(self, nprocs=nprocs)
+        optimized, energy = opt.optimize_molecule(self, nprocs=nprocs, return_energy=True)
 
         if inplace:
             self.geometry = optimized.geometry
-            return self
+            if return_energy:
+                return self, energy
+            else:
+                return self
         else:
-            return optimized
+            if return_energy:
+                return optimized, energy
+            else:
+                return optimized
+
+    def compute_energy(self, nprocs=1):
+        """
+        Compute energy of molecule at the GFN2-xtb level of theory.
+
+        Args:
+            nprocs (int): number of processors to use
+        """
+        import cctk.optimize as opt
+        assert isinstance(nprocs, int), "nprocs must be int!"
+        energy = opt.get_energy(self, nprocs=nprocs)
+        return energy
 
     def csearch(self, nprocs=1, constraints=[], logfile=None, noncovalent=False, use_tempdir=True):
         """

cctk/optimize.py

@@ -41,9 +41,28 @@ def optimize_molecule(molecule, method=Methods.GFN2_XTB, nprocs=1, return_energy
     assert isinstance(method, Methods), "need a valid molecule!"
 
     if method is Methods.GFN2_XTB:
-        return run_xtb(molecule, nprocs=nprocs, return_energy=return_energy)
+        return run_xtb(molecule, nprocs=nprocs, return_energy=return_energy, opt=True)
 
-def run_xtb(molecule, nprocs=1, return_energy=False):
+def get_energy(molecule, method=Methods.GFN2_XTB, nprocs=1):
+    """
+    Dispatcher method to connect method to molecule.
+
+    Args:
+        molecule (cctk.Molecule):
+        method (Methods):
+        nprocs (int): number of cores to employ
+
+    Returns:
+        energy
+    """
+    assert isinstance(molecule, cctk.Molecule), "need a valid molecule!"
+    assert isinstance(method, Methods), "need a valid molecule!"
+
+    if method is Methods.GFN2_XTB:
+        _, energy = run_xtb(molecule, nprocs=nprocs, return_energy=True, opt=False)
+        return energy
+
+def run_xtb(molecule, nprocs=1, return_energy=False, opt=False):
     """
     Run ``xtb`` in a temporary directory and return the output molecule.
     """
@@ -55,7 +74,11 @@ def run_xtb(molecule, nprocs=1, return_energy=False):
     command = f"xtb --gfn 2 --chrg {molecule.charge} --uhf {molecule.multiplicity - 1}"
     if nprocs > 1:
         command += f" --parallel {nprocs}"
-    command += " xtb-in.xyz --opt tight &> xtb-out.out"
+
+    if opt:
+        command += " xtb-in.xyz --opt tight &> xtb-out.out"
+    else:
+        command += " xtb-in.xyz &> xtb-out.out"
 
     try:
         os.environ["OMP_NUM_THREADS"] = str(nprocs)
@@ -64,8 +87,13 @@ def run_xtb(molecule, nprocs=1, return_energy=False):
             cctk.XYZFile.write_molecule_to_file(f"{tmpdir}/xtb-in.xyz", molecule)
             result = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, cwd=tmpdir, shell=True)
 
-            output_mol = cctk.XYZFile.read_file(f"{tmpdir}/xtbopt.xyz").molecule
-            energy_file = cctk.File.read_file(f"{tmpdir}/xtbopt.log")
+            output_mol, energy_file = None, None
+            if opt:
+                output_mol = cctk.XYZFile.read_file(f"{tmpdir}/xtbopt.xyz").molecule
+                energy_file = cctk.File.read_file(f"{tmpdir}/xtbopt.log")
+            else:
+                print("Need to read energy from xtb-out.out!")
+                raise NotImplementedError
 
             fields = energy_file[1].split()
             energy, gradient = float(fields[1]), float(fields[3])

cctk/parse_gaussian.py

@@ -105,7 +105,9 @@ def read_file_fast(file_text, filename, link1idx, max_len=10000, extended_opt_in
         "Temperature", #15
         "Isotropic",
         "EUMP2",
+        "EUMP3",
         "UMP4(SDTQ)",
+        "Wavefunction amplitudes converged", #20
     ]
 
     #### And here are the blocks of text
@@ -220,8 +222,14 @@ def read_file_fast(file_text, filename, link1idx, max_len=10000, extended_opt_in
     # post-HF methods give weird energies
     if re.search("mp2", route_card, re.IGNORECASE):
         energies = parse_mp2_energies(word_matches[17])
+    elif re.search("mp3", route_card, re.IGNORECASE):
+        energies = parse_mp3_energies(word_matches[18])
     elif re.search("mp4", route_card, re.IGNORECASE):
-        energies = parse_mp4_energies(word_matches[18])
+        energies = parse_mp4_energies(word_matches[19])
+    elif re.search("ccsd", route_card, re.IGNORECASE):
+        energies = parse_cc_energies(word_matches[20])
+    elif re.search("cisd", route_card, re.IGNORECASE):
+        energies = parse_ci_energies(word_matches[20])
 
     f = cctk.GaussianFile(job_types=job_types, route_card=route_card, link0=link0, footer=footer, success=success, elapsed_time=elapsed_time, title=title)
 
@@ -692,6 +700,15 @@ def parse_mp2_energies(lines):
         energies.append(float(energy_str))
     return energies
 
+def parse_mp3_energies(lines):
+    energies = []
+    for line in lines:
+        pieces = list(filter(None, line.split(" ")))
+        energy_str = pieces[3]
+        energy_str = re.sub("D", "E", energy_str)
+        energies.append(float(energy_str))
+    return energies
+
 def parse_mp4_energies(lines):
     energies = []
     for line in lines:
@@ -700,3 +717,15 @@ def parse_mp4_energies(lines):
         energy_str = re.sub("D", "E", energy_str)
         energies.append(float(energy_str))
     return energies
+
+def parse_cc_energies(lines):
+    energies = []
+    for line in lines:
+        pieces = list(filter(None, line.split(" ")))
+        energy_str = pieces[4]
+        energy_str = re.sub("D", "E", energy_str)
+        energies.append(float(energy_str))
+    return energies
+
+def parse_ci_energies(lines):
+    return parse_cc_energies(lines)

cctk/xyz_file.py

@@ -51,6 +51,10 @@ def file_from_lines(cls, lines):
             try:
                 if re.match("[0-9]", pieces[0]):
                     atomic_numbers[index] = int(pieces[0])
+                elif re.match("([A-Za-z])+([0-9])+", pieces[0]):
+                    # mdtraj writes in this format, for some reason
+                    m = re.match("([A-Za-z])+([0-9])+", pieces[0])
+                    atomic_numbers[index] = int(get_number(m.group(1)))
                 else:
                     atomic_numbers[index] = int(get_number(pieces[0]))
                 geometry[index][0] = float(pieces[1])

test/test_opt.py

@@ -10,8 +10,12 @@ def load_molecule(self, path="test/static/test_peptide.xyz"):
 
     def test_basic(self):
         mol = self.load_molecule()
-        mol2 = opt.optimize_molecule(mol, nprocs=4)
+#        e1 = mol.compute_energy()
+#        print(e1)
+
+        mol2, e2 = opt.optimize_molecule(mol, nprocs=4, return_energy=True)
         self.assertTrue(isinstance(mol2, cctk.Molecule))
+        self.assertTrue(e2 + 66.394 < 0.1)
 
         mol.optimize(nprocs=4) #in-place