Merge pull request CrawfordGroup#8 from bgpeyton/component

Kick

pycc/rt/rtcc.py

@@ -29,7 +29,14 @@ class rtcc(object):
     mu_tot: NumPy arrays
         1/sqrt(3) * sum of dipole integrals (for isotropic field)
     m: list of NumPy arrays
-        the magnetic dipole integrals for each Cartesian direction (only if dipole = True)
+        the magnetic dipole integrals for each Cartesian direction (only if magnetic = True)
+
+    Parameters
+    ----------
+    magnetic: bool
+        optionally store magnetic dipole integrals (default = False)
+    kick: bool or str
+        optionally isolate 'x', 'y', or 'z' electric field kick (default = False)
 
     Methods
     -------
@@ -46,7 +53,7 @@ class rtcc(object):
     lagrangian()
         Compute the CC Lagrangian energy for a given time t
     """
-    def __init__(self, ccwfn, cclambda, ccdensity, V, magnetic = False):
+    def __init__(self, ccwfn, cclambda, ccdensity, V, magnetic = False, kick = None):
         self.ccwfn = ccwfn
         self.cclambda = cclambda
         self.ccdensity = ccdensity
@@ -59,13 +66,18 @@ def __init__(self, ccwfn, cclambda, ccdensity, V, magnetic = False):
         self.mu = []
         for axis in range(3):
             self.mu.append(C.T @ np.asarray(dipole_ints[axis]) @ C)
-        self.mu_tot = sum(self.mu)/np.sqrt(3.0)  # isotropic field
+        if kick:
+            s_to_i = {"x":0, "y":1, "z":2}
+            self.mu_tot = self.mu[s_to_i[kick.lower()]]
+        else:
+            self.mu_tot = sum(self.mu)/np.sqrt(3.0)  # isotropic field
 
         if magnetic:
             m_ints = mints.ao_angular_momentum()
             self.m = []
             for axis in range(3):
-                self.m.append(C.T @ (np.asarray(m_ints[axis])*-0.5) @ C)
+                m = (C.T @ (np.asarray(m_ints[axis])*-0.5) @ C)
+                self.m.append(m*1.0j)
 
     def f(self, t, y):
         """

pycc/tests/test_007_dipole.py

@@ -52,6 +52,6 @@ def test_dipole_h2_2_cc_pvdz():
     ref = [0, 0, -2.3037968376087573E-5]
     m_x, m_y, m_z = rtcc.dipole(t1, t2, l1, l2, magnetic = True)
 
-    assert (abs(ref[0] - m_x) < 1E-10)
-    assert (abs(ref[1] - m_y) < 1E-10)
-    assert (abs(ref[2] - m_z) < 1E-10)
+    assert (abs(ref[0]*1.0j - m_x) < 1E-10)
+    assert (abs(ref[1]*1.0j - m_y) < 1E-10)
+    assert (abs(ref[2]*1.0j - m_z) < 1E-10)

pycc/tests/test_014_field.py

@@ -0,0 +1,77 @@
+"""
+Test electric field generation 
+"""
+
+# Import package, test suite, and other packages as needed
+import psi4
+import pycc
+import pytest
+import numpy as np
+from ..data.molecules import *
+from pycc.rt.lasers import gaussian_laser
+
+def test_dipole_h2_2_field():
+    """H2 dimer"""
+    psi4.set_memory('2 GiB')
+    psi4.core.set_output_file('output.dat', False)
+    psi4.set_options({'basis': '6-31G',
+                      'scf_type': 'pk',
+                      'mp2_type': 'conv',
+                      'freeze_core': 'false',
+                      'e_convergence': 1e-13,
+                      'd_convergence': 1e-13,
+                      'r_convergence': 1e-13,
+                      'diis': 1})
+    mol = psi4.geometry(moldict["(H2)_2"])
+    rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+    e_conv = 1e-13
+    r_conv = 1e-13
+
+    cc = pycc.ccenergy(rhf_wfn)
+    ecc = cc.solve_ccsd(e_conv, r_conv)
+
+    hbar = pycc.cchbar(cc)
+
+    cclambda = pycc.cclambda(cc, hbar)
+    lecc = cclambda.solve_lambda(e_conv, r_conv)
+
+    ccdensity = pycc.ccdensity(cc, cclambda)
+
+    # narrow Gaussian pulse
+    F_str = 0.01
+    sigma = 0.01
+    center = 0.05
+    V = gaussian_laser(F_str, 0, sigma, center=center)
+    rtcc = pycc.rtcc(cc, cclambda, ccdensity, V, magnetic = True)
+
+    mints = psi4.core.MintsHelper(cc.ref.basisset())
+    dipole_ints = mints.ao_dipole()
+    m_ints = mints.ao_angular_momentum()
+    C = np.asarray(cc.ref.Ca_subset("AO", "ACTIVE"))
+    ref_mu = []
+    ref_m = []
+    for axis in range(3):
+        ref_mu.append(C.T @ np.asarray(dipole_ints[axis]) @ C)
+        ref_m.append(C.T @ (np.asarray(m_ints[axis])*-0.5) @ C)
+        assert np.allclose(ref_mu[axis],rtcc.mu[axis])
+        assert np.allclose(ref_m[axis]*1.0j,rtcc.m[axis])
+    ref_mu_tot = sum(ref_mu)/np.sqrt(3.0)
+
+    assert np.allclose(ref_mu_tot,rtcc.mu_tot)
+
+    rtcc = pycc.rtcc(cc, cclambda, ccdensity, V, magnetic = True, kick="Y")
+
+    mints = psi4.core.MintsHelper(cc.ref.basisset())
+    dipole_ints = mints.ao_dipole()
+    m_ints = mints.ao_angular_momentum()
+    C = np.asarray(cc.ref.Ca_subset("AO", "ACTIVE"))
+    ref_mu = []
+    ref_m = []
+    for axis in range(3):
+        ref_mu.append(C.T @ np.asarray(dipole_ints[axis]) @ C)
+        ref_m.append(C.T @ (np.asarray(m_ints[axis])*-0.5) @ C)
+        assert np.allclose(ref_mu[axis],rtcc.mu[axis])
+        assert np.allclose(ref_m[axis]*1.0j,rtcc.m[axis])
+
+    assert np.allclose(ref_mu[1],rtcc.mu_tot)