finished (i think) quasiclassical business

cctk/helper_functions.py

@@ -299,7 +299,7 @@ def compute_RMSD(geometry1, geometry2, checks=True):
         raise ValueError("can't compare two geometries with different lengths!")
 
 #    return np.trace(cdist(geometry1, geometry2)) / len(geometry1)
-    return np.sqrt( np.sum( ( geometry1 - geometry2 ) ** 2) / len(geometry1) )
+    return np.sqrt( np.sum( ( geometry1.view(np.ndarray) - geometry2.view(np.ndarray) ) ** 2) / len(geometry1) )
 
 def get_isotopic_distribution(z):
     """

cctk/molecule.py

@@ -82,6 +82,8 @@ def __init__(self, atomic_numbers, geometry, name=None, bonds=None, charge=0, mu
         self.multiplicity = multiplicity
         self.charge = charge
 
+        self.vibrational_modes = list()
+
         if isinstance(bonds, nx.Graph):
             self.bonds = bonds
         elif isinstance(bonds, (list,np.ndarray,nx.classes.reportviews.EdgeView)):

cctk/parse_gaussian.py

@@ -711,7 +711,7 @@ def read_file_fast(file_text, filename, link1idx, max_len=10000, extended_opt_in
             raise ValueError(f"unexpected # gibbs free energies found!\ngibbs free energies = {gibbs_vals}")
 
 
-        vibrational_modes = parse_modes(block_matches[11])
+        vibrational_modes = parse_modes(block_matches[11], num_atoms=molecules[-1].num_atoms(), hpmodes=re.search("hpmodes", route_card))
         molecules[-1].vibrational_modes = vibrational_modes
 
         frequencies = []
@@ -923,39 +923,67 @@ def parse_hirshfeld(hirshfeld_block):
 
     return cctk.OneIndexedArray(charges), cctk.OneIndexedArray(spins)
 
-def parse_modes(freq_block):
+def parse_modes(freq_block, num_atoms, hpmodes=False):
     freqs = list()
     masses = list()
     force_ks = list()
     displacements = list()
 
-    chunks = freq_block[0].split("\n Freq")
+    chunks = freq_block[0].split("Freq")
+
+    if hpmodes:
+        chunks = chunks[1:]
+        chunks = chunks[:int(len(chunks)/2)]
+
     for chunk in chunks:
         lines = chunk.split("\n")
-        freqs += re.split(" +", lines[0])[2:]
-        masses += re.split(" +", lines[1])[4:]
-        force_ks += re.split(" +", lines[2])[4:]
 
-        current_displacements = [list() for x in re.split(" +", lines[0])[2:]]
+        if hpmodes:
+            num_cols = len(re.split(" +", lines[0])) - 2
+            current_displacements = [np.zeros(shape=(num_atoms, 3)) for x in range(num_cols)]
 
-        for line in lines[5:]:
-            fields = re.split(" +", line)
-            fields = list(filter(None, fields))
+            freqs += list(filter(None, re.split(" +", lines[0])))[2:]
+            masses += list(filter(None, re.split(" +", lines[1])))[3:]
+            force_ks += list(filter(None, re.split(" +", lines[2])))[3:]
 
-            if len(fields) < 5:
-                break
+            for line in lines[6:]:
+                fields = re.split(" +", line)
+                fields = list(filter(None, fields))
+
+                if len(fields) < 5 or fields[0] == "Harmonic":
+                    break
+
+                for col_idx, val in enumerate(fields[3:]):
+                    current_displacements[col_idx][int(fields[1])-1][int(fields[0])-1] = val
+
+            for d in current_displacements:
+                displacements.append(d.view(cctk.OneIndexedArray))
+
+        else:
+            current_displacements = [list() for x in re.split(" +", lines[0])[2:]]
+
+            freqs += re.split(" +", lines[0])[2:]
+            masses += re.split(" +", lines[1])[4:]
+            force_ks += re.split(" +", lines[2])[4:]
+
+            for line in lines[5:]:
+                fields = re.split(" +", line)
+                fields = list(filter(None, fields))
+
+                if len(fields) < 4:
+                    break
 
-            current_displacements[0].append([float(x) for x in fields[2:5]])
+                current_displacements[0].append([float(x) for x in fields[2:5]])
 
-            if len(current_displacements) > 1:
-                current_displacements[1].append([float(x) for x in fields[5:8]])
+                if len(current_displacements) > 1:
+                    current_displacements[1].append([float(x) for x in fields[5:8]])
 
-            if len(current_displacements) > 2:
-                current_displacements[2].append([float(x) for x in fields[8:11]])
+                if len(current_displacements) > 2:
+                    current_displacements[2].append([float(x) for x in fields[8:11]])
 
 
-        for d in current_displacements:
-            displacements.append(cctk.OneIndexedArray(d))
+            for d in current_displacements:
+                displacements.append(cctk.OneIndexedArray(d))
 
     freqs = [float(x) for x in freqs]
     masses = [float(x) for x in masses]

cctk/quasiclassical.py

@@ -7,8 +7,24 @@
 
 import cctk
 
-def get_quasiclassical_vibration(molecule, temperature=298):
+def get_quasiclassical_perturbation(molecule, temperature=298):
+    """
+    Perturbs molecule by treating each mode as a quantum harmonic oscillator and sampling from the distribution appropriate to the temperature.
+
+    This is probably the only useful function in this file.
+
+    Args:
+        molecule (cctk.Molecule): molecule with vibrational modes
+        temperature (float): temperature
+
+    Returns:
+        new ``cctk.Molecule`` object
+        energy above ground state (kcal/mol)
+    """
     assert isinstance(molecule, cctk.Molecule), "need a valid molecule"
+    assert len(molecule.vibrational_modes) > 0, "molecule needs to have vibrational modes (try running a ``freq`` job)"
+
+    assert isinstance(temperature, (int, float)), "temperature must be numeric"
 
     mol = copy.deepcopy(molecule)
     total_PE = 0
@@ -17,12 +33,12 @@ def get_quasiclassical_vibration(molecule, temperature=298):
         PE, KE = apply_vibration(mol, mode, temperature=temperature)
         total_PE += PE
 
-    print(total_PE)
-
-    return mol
+    return mol, total_PE
 
-def apply_vibration(molecule, mode, min_freq=50, temperature=298):
+def apply_vibration(molecule, mode, min_freq=50, temperature=298, verbose=False):
     """
+    Apply a vibration to molecule ``molecule`` (modified in-place).
+
     Args:
         molecule (cctk.Molecule)
         mode (cctk.VibrationalMode)
@@ -50,12 +66,11 @@ def apply_vibration(molecule, mode, min_freq=50, temperature=298):
     kinetic_energy = energy - potential_energy
 
     # here we could compute atom velocities if we wanted to! initializer lines 440-480
-
-    print(f"Mode {mode.frequency:.2f} ({mode.energy():.2f} kcal/mol)\t QC Level {level}\t Shift {rel_shift:.2%} of a potential {max_shift:.2f} Å\tPE = {potential_energy:.2f} kcal/mol\tk = {mode.force_constant:.2f} kcal/mol Å^-2")
+    if verbose:
+        print(f"Mode {mode.frequency:.2f} ({mode.energy():.2f} kcal/mol)\t QC Level {level}\t Shift {rel_shift:.2%} of a potential {max_shift:.2f} Å\tPE = {potential_energy:.2f} kcal/mol\tk = {mode.force_constant:.2f} kcal/mol Å^-2")
 
     return potential_energy, kinetic_energy
 
-
 def get_hermite_polynomial(n):
     """
     Returns a ``np.poly1d`` object representing the degree-n Hermite polynomial.

cctk/vibrational_mode.py

@@ -55,7 +55,7 @@ def choose_level(self, temperature=298):
 
         # probability of being in state 0 is equal to 1 - zpe_ratio
         # 1 = P(0) + P(1) + P(2) + ... = P + P * zpe_ratio + P * zpe_ratio ** 2 + ...
-        # 1 = P(0) / (1 - zpe_ratio) # geometric series
+        # 1 = P(0) / (1 - zpe_ratio) bc geometric series
         P = 1.0 - zpe_ratio
 
         random = np.random.uniform()
@@ -147,6 +147,9 @@ def quantum_distribution_value(self, x, level=0):
         return val
 
     def quantum_distribution_max(self, level=0, num_pts=1e5):
+        """
+        Returns the maximum value of psi**2 for the quantum harmonic oscillator at a given level.
+        """
         assert isinstance(level, int) and level >= 0, "need positive integer for vibrational level"
 
         freq = self.frequency
@@ -158,6 +161,7 @@ def quantum_distribution_max(self, level=0, num_pts=1e5):
 
         max_x = self.classical_turning_point()
 
+        # there is certainly a better way to do this
         max_p = 0
         for x in np.linspace(0, max_x, num_pts):
             p = self.quantum_distribution_value(x, level)
@@ -167,4 +171,7 @@ def quantum_distribution_max(self, level=0, num_pts=1e5):
         return max_p
 
     def classical_turning_point(self, level=0):
+        """
+        Returns the maximum allowed shift based on modelling the mode as a classical harmonic oscillator (e.g. the point where potential energy is maximum).
+        """
         return math.sqrt(2 * self.energy(level) / self.force_constant)

test/qc_by_mode.py

@@ -0,0 +1,38 @@
+import numpy as np
+import sys
+
+sys.path.insert(0,'/Users/cwagen/code/cctk')
+
+import cctk
+import cctk.quasiclassical as qc
+
+#### This is a script to let you specify vibrational excited states manually, to compare with values from other programs.
+#### This was used to test against Jprogdyn.
+#### Corin Wagen, July 2020
+
+path = "test/static/methane_hpmodes.out"
+file = cctk.GaussianFile.read_file(path)
+
+mol = file.get_molecule()
+
+total_PE = 0
+
+for mode in mol.vibrational_modes:
+    print(mode)
+
+    level = 0
+    energy = mode.energy(level)
+
+    rel_shift = float(input("rel_shift (%):")) / 100
+    max_shift = mode.classical_turning_point(level)
+    shift = rel_shift * max_shift
+
+    mode_coords = mode.displacements
+    mol.geometry += mode.displacements * rel_shift * max_shift
+
+    potential_energy = 0.5 * mode.force_constant * shift ** 2
+    total_PE += potential_energy
+
+    print(f"Mode {mode.frequency:.2f} ({mode.energy():.2f} kcal/mol)\t QC Level {level}\t Shift {rel_shift:.2%} of a potential {max_shift:.2f} Å\tPE = {potential_energy:.2f} kcal/mol\tk = {mode.force_constant:.2f} kcal/mol Å^-2")
+
+cctk.GaussianFile.write_molecule_to_file("test/static/methane_perturbed.gjf", mol, route_card="#p sp hf/sto-3g")

test/static/methane_perturbed.gjf

@@ -0,0 +1,13 @@
+%mem=32GB
+%nprocshared=16
+#p sp hf/sto-3g
+
+title
+
+0 1
+ 6          0.00000000    0.00243584   -0.00822321
+ 1          0.57417619    0.57646316    0.67574662
+ 1         -0.63788247   -0.65552473    0.76005989
+ 1         -0.61897093    0.64186412   -0.69135213
+ 1          0.59476799   -0.59180391   -0.64654410
+

test/static/methane_perturbed_key.gjf

@@ -0,0 +1,10 @@
+#p hf/sto-3g
+
+title
+
+0 1
+   C            0.0073293454       0.0025407821      -0.0081940867
+   H            0.5744416866       0.5757139319       0.6754784436
+   H           -0.6376757688      -0.6555328380       0.7597691638
+   H           -0.6190199611       0.6413481532      -0.6913898112
+   H            0.5949791917      -0.5917801823      -0.6462942546

test/test_freqs.py

@@ -16,10 +16,12 @@ def test_read(self):
         self.assertEqual(len(mol.vibrational_modes), 9)
         self.assertEqual(mol.vibrational_modes[-1].frequency, 3119.6807)
 
-        mol2 = qc.get_quasiclassical_vibration(mol)
+        path = "test/static/methane_hpmodes.out"
+        file = cctk.GaussianFile.read_file(path)
 
-    def test_qho(self):
+    def draw_qho(self):
         from asciichartpy import plot
+
         path = "test/static/methane_normal.out"
         file = cctk.GaussianFile.read_file(path)
         mol = file.get_molecule()
@@ -30,5 +32,37 @@ def test_qho(self):
             tp = mode.classical_turning_point(level=level)
             print(f"LEVEL {level} ({mode.energy(level):.2f} kcal.mol, turning point = {tp:.2f} Å)")
             tp = 5
-            print(plot([mode.quantum_distribution_value(x, level=level) for x in np.linspace(-1 * tp, tp, 100)], {"height": 8}))
+            print(plot([mode.quantum_distribution_value(x, level=level) for x in np.linspace(-1 * tp, tp, 100)], {"height": 10}))
+
+    def test_perturb(self):
+        path = "test/static/methane_hpmodes.out"
+        file = cctk.GaussianFile.read_file(path)
+
+        mol = file.get_molecule()
+
+        energies = list()
+        for i in range(1000):
+            _, e = qc.get_quasiclassical_perturbation(mol)
+            energies.append(e)
+
+        energies = np.array(energies)
+#        from asciichartpy import plot
+#        hist, edges = np.histogram(energies, bins=120)
+#        print(plot(hist, {"height": 10}))
+
+        self.assertTrue(abs(np.mean(energies) - 9.4) < 0.5)
+        self.assertTrue(abs(np.std(energies) - 3) < 0.5)
+
+    def test_final_structure(self):
+        path1 = "test/static/methane_perturbed.gjf"
+        path2 = "test/static/methane_perturbed_key.gjf"
+
+        e = cctk.ConformationalEnsemble()
+
+        e.add_molecule(cctk.GaussianFile.read_file(path1).get_molecule().assign_connectivity())
+        e.add_molecule(cctk.GaussianFile.read_file(path2).get_molecule().assign_connectivity().renumber_to_match(e.molecules[0]))
+
+        e2, before, after = e.align(comparison_atoms="all", compute_RMSD=True)
+
+        self.assertTrue(after[1] < 0.005)