working on quasicalssical excitations

cctk/__init__.py

@@ -4,6 +4,7 @@
 from .molecule import Molecule
 from .ensemble import Ensemble, ConformationalEnsemble
 from .group import Group
+from .vibrational_mode import VibrationalMode
 
 from .gaussian_file import GaussianJobType, GaussianFile
 from .orca_file import OrcaFile, OrcaJobType

cctk/molecule.py

@@ -36,6 +36,7 @@ class Molecule:
         charge (int): the charge of the molecule
         multiplicity (int): the spin state of the molecule (1 corresponds to singlet, 2 to doublet, 3 to triplet, etc. -- so a multiplicity of 1 is equivalent to S=0)
         checks (bool): whether to check that the constructor parameters are valid
+        vibrational_modes (list of cctk.VibrationalMode): vibrational modes
     """
 
     def __init__(self, atomic_numbers, geometry, name=None, bonds=None, charge=0, multiplicity=1, checks=True):

cctk/parse_gaussian.py

@@ -509,7 +509,7 @@ def extract_success_and_time(lines):
         elapsed_time += fields[0] * 86400 + fields[1] * 3600 + fields[2] * 60 + fields[3]
     return len(success), elapsed_time
 
-def read_file_fast(file_text, filename, link1idx, max_len=1000, extended_opt_info=False):
+def read_file_fast(file_text, filename, link1idx, max_len=10000, extended_opt_info=False):
 
     #### "Make your bottleneck routines fast, everything else clear" - M. Scott Shell, UCSB
     #### Welcome to the fast part!
@@ -554,6 +554,7 @@ def read_file_fast(file_text, filename, link1idx, max_len=1000, extended_opt_inf
         ["Hirshfeld charges, spin densities, dipoles, and CM5 charges", " Hirshfeld charges", 1],
         ["Mulliken charges", "Sum of Mulliken charges", 1],
         ["Electronic spatial extent", "Quadrupole moment", 1], #10
+        ["normal coordinates", "Thermochemistry", 1],
     ]
 
     word_matches = [[] for _ in words]
@@ -709,6 +710,10 @@ def read_file_fast(file_text, filename, link1idx, max_len=1000, extended_opt_inf
         elif len(gibbs_vals) > 1:
             raise ValueError(f"unexpected # gibbs free energies found!\ngibbs free energies = {gibbs_vals}")
 
+
+        vibrational_modes = parse_modes(block_matches[11])
+        molecules[-1].vibrational_modes = vibrational_modes
+
         frequencies = []
         try:
             frequencies += extract_parameter(word_matches[14], 2)
@@ -918,3 +923,51 @@ def parse_hirshfeld(hirshfeld_block):
 
     return cctk.OneIndexedArray(charges), cctk.OneIndexedArray(spins)
 
+def parse_modes(freq_block):
+    freqs = list()
+    masses = list()
+    force_ks = list()
+    displacements = list()
+
+    chunks = freq_block[0].split("\n Freq")
+    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:]]
+
+        for line in lines[5:]:
+            fields = re.split(" +", line)
+            fields = list(filter(None, fields))
+
+            if len(fields) < 5:
+                break
+
+            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) > 2:
+                current_displacements[2].append([float(x) for x in fields[8:11]])
+
+
+        for d in current_displacements:
+            displacements.append(cctk.OneIndexedArray(d))
+
+    freqs = [float(x) for x in freqs]
+    masses = [float(x) for x in masses]
+    force_ks = [float(x) for x in force_ks]
+
+    assert len(freqs) == len(masses)
+    assert len(freqs) == len(force_ks)
+    assert len(freqs) == len(displacements)
+
+    modes = list()
+    for f, m, k, d in zip(freqs, masses, force_ks, displacements):
+        k *= 143.836 # mdyne Å**-1 to kcal/mol Å**-2
+        modes.append(cctk.VibrationalMode(f, m, k, d))
+
+    return modes

cctk/quasiclassical.py

@@ -0,0 +1,78 @@
+"""
+Functions to assist in sampling thermally excited states through quasiclassical approximations.
+"""
+
+import numpy as np
+import math, copy
+
+import cctk
+
+def get_quasiclassical_vibration(molecule, temperature=298):
+    assert isinstance(molecule, cctk.Molecule), "need a valid molecule"
+
+    mol = copy.deepcopy(molecule)
+    total_PE = 0
+
+    for mode in mol.vibrational_modes:
+        PE, KE = apply_vibration(mol, mode, temperature=temperature)
+        total_PE += PE
+
+    print(total_PE)
+
+    return mol
+
+def apply_vibration(molecule, mode, min_freq=50, temperature=298):
+    """
+    Args:
+        molecule (cctk.Molecule)
+        mode (cctk.VibrationalMode)
+
+    Returns:
+        potential energy
+        kinetic energy
+    """
+
+    level = mode.choose_level(temperature)
+    energy = mode.energy(level)
+
+    shift = mode.random_displacement(level)
+    max_shift = mode.classical_turning_point(level)
+    rel_shift = shift/max_shift
+
+    if mode.frequency < min_freq:
+        rel_shift = 0
+
+    mode_coords = mode.displacements
+    molecule.geometry += mode.displacements * rel_shift * max_shift
+
+    potential_energy = 0.5 * mode.force_constant * shift ** 2
+
+    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}\t QC Level {level}\t Shift {rel_shift:.2%} of a potential {max_shift:.2f} Å\tPE = {potential_energy:.2f} kcal/mol")
+
+    return potential_energy, kinetic_energy
+
+
+def get_hermite_polynomial(n):
+    """
+    Returns a ``np.poly1d`` object representing the degree-n Hermite polynomial.
+
+    Adapted from https://scipython.com/blog/the-harmonic-oscillator-wavefunctions/.
+    """
+    assert isinstance(n, int) and n >= 0, "need positive integer"
+
+    Hr = [None] * (n + 1)
+    Hr[0] = np.poly1d([1.,])
+
+    if n > 0:
+        Hr[1] = np.poly1d([2., 0.])
+
+    if n > 1:
+        for v in range(2, n+1):
+            Hr[v] = Hr[1]*Hr[v-1] - 2*(v-1)*Hr[v-2]
+    return Hr[n]
+
+

cctk/vibrational_mode.py

@@ -0,0 +1,171 @@
+import math, copy, re
+import numpy as np
+
+import cctk
+from cctk.quasiclassical import get_hermite_polynomial
+
+# constants
+MAX_QHO_LEVEL = 10000
+MIN_FREQUENCY = 2
+MIN_TEMPERATURE = 10
+MAX_ZPE_RATIO = 0.999999
+
+BOLTZMANN_CONSTANT = 0.001985875 # kcal/mol•K
+
+class VibrationalMode:
+    """
+    Most code adapted from ``jprogdyn``. Displacements will be very low accuracy unless ``freq=hpmodes`` is enabled.
+
+    Values from Gaussian, for now: see https://gaussian.com/vib/.
+
+    Attributes:
+        frequency (float): frequency, in cm-1
+        force_constant (float): force constant, in kcal/mol per Å
+        reduced_mass (float): mass, in amus
+        displacements (cctk.OneIndexedArray): displacements
+
+    """
+    def __init__(self, frequency, force_constant, reduced_mass, displacements):
+        assert isinstance(frequency, float)
+        self.frequency = frequency
+
+        assert isinstance(force_constant, float)
+        self.force_constant = force_constant
+
+        assert isinstance(reduced_mass, float)
+        self.reduced_mass = reduced_mass
+
+        assert isinstance(displacements, cctk.OneIndexedArray)
+        self.displacements = displacements
+
+    def __str__(self):
+        return f"Vibrational mode ({self.frequency:.2f} cm-1, {self.reduced_mass:.2f} amus, {self.force_constant:.2f} kcal/mol Å**-2)"
+
+    def __repr__(self):
+        return f"Vibrational mode ({self.frequency:.2f} cm-1, {self.reduced_mass:.2f} amus, {self.force_constant:.2f} kcal/mol Å**-2)"
+
+    def choose_level(self, temperature=298):
+        if temperature < MIN_TEMPERATURE:
+            return 0
+
+        # zpe_ratio is probability of being in level i vs level i+1, by quantum harmonic oscillator
+        zpe_ratio = math.exp( -2 * self.energy() / BOLTZMANN_CONSTANT * temperature)
+        if zpe_ratio > MAX_ZPE_RATIO:
+            zpe_ratio = MAX_ZPE_RATIO
+
+        # 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
+        P = 1.0 - zpe_ratio
+
+        random = np.random.uniform()
+        level = 0
+        while level < MAX_QHO_LEVEL:
+            if random < P:
+                return level
+            else:
+                P += P * zpe_ratio
+                level += 1
+
+        return level
+
+    def energy(self, level=0):
+        """
+        Calculate energy as a function of level. By default returns zero-point energy (level = 0).
+
+        Args:
+            level (int): which vibrational level the mode is in
+
+        Returns:
+            energy (kcal/mol)
+        """
+        assert isinstance(level, int) and level >= 0, "need positive integer for vibrational level"
+
+        freq = self.frequency
+        if freq < MIN_FREQUENCY:
+            freq = MIN_FREQUENCY
+
+        # 0.5 * h * c * frequency (c in cm/s bc wavenumbers)
+        # 0.5 * (6.626 * 10**-34) * (3 * 10**10) * (6.026 * 10**23) / 4184) = 0.0014305 
+        zpe = 0.0014305 * freq
+        return zpe * (2 * level + 1)
+
+    def random_displacement(self, level=0, method="quasiclassical", max_attempts=1e5):
+        """
+        Args:
+            method (str): "quasiclassical" for now
+            level (int): which vibrational level
+            max_attempts (int): how many tries you get
+
+        Returns:
+            shift
+        """
+        energy = self.energy(level)
+
+        if method == "quasiclassical":
+            min_val = 0
+            max_val = self.quantum_distribution_max(level)
+            max_x = self.classical_turning_point()
+
+            attempts = 0
+            while attempts < max_attempts:
+                x = np.random.uniform(-1 * max_x, max_x)
+                p = self.quantum_distribution_value(x, level)
+
+                y = np.random.uniform(min_val, max_val)
+                if y < p:
+                    return x
+                else:
+                    attempts += 1
+
+            raise ValueError("max_attempts exceeded - can't get a proper initialization for this mode!")
+        else:
+            raise ValueError(f"invalid method {method} - only ``quasiclassical`` implemented currently!")
+
+
+    def quantum_distribution_value(self, x, level=0):
+        """
+        Calculate psi**2 for quantum harmonic oscillator for a given shift in Å.
+
+        Args:
+            x (float): shift in Å
+            level (int): vibrational level
+        """
+        assert isinstance(level, int) and level >= 0, "need positive integer for vibrational level"
+
+        freq = self.frequency
+        if freq < MIN_FREQUENCY:
+            freq = MIN_FREQUENCY
+
+        n = level # brevity is the soul of wit
+        H = get_hermite_polynomial(n)
+
+        # following https://github.com/ekwan/Jprogdyn/blob/master/src/main/java/edu/harvard/chemistry/ekwan/Jprogdyn/HarmonicOscillatorDistribution.java, line 109
+
+        # 4 * pi * 3 * 10**10 / (1000 * 10**23 * 6.022 * 10**23 * 6.626 * 10^-34) = 9.448 * 10 ** -6, take it or leave it
+        omega_term = 9.448 * 10 ** -6 * self.reduced_mass * freq
+        val = math.sqrt(omega_term) * math.exp(-1 * omega_term * math.pi * x ** 2 ) * (H(math.sqrt(omega_term * math.pi) * x) ** 2) / (2 ** n * math.factorial(n))
+        return val
+
+    def quantum_distribution_max(self, level=0, num_pts=1e5):
+        assert isinstance(level, int) and level >= 0, "need positive integer for vibrational level"
+
+        freq = self.frequency
+        if freq < MIN_FREQUENCY:
+            freq = MIN_FREQUENCY
+
+        if level == 0:
+            return self.quantum_distribution_value(0)
+
+        max_x = self.classical_turning_point()
+
+        max_p = 0
+        for x in np.linspace(0, max_x, num_pts):
+            p = self.quantum_distribution_value(x, level)
+            if p > max_p:
+                max_p = p
+
+        return max_p
+
+    def classical_turning_point(self, level=0):
+        return math.sqrt(2 * self.energy(level) / self.force_constant)