diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 5e2cd1b..2055b3d 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -21,8 +21,8 @@ jobs: strategy: fail-fast: true matrix: - os: [macOS-latest, ubuntu-latest] - python-version: ['3.8', '3.9', '3.10'] + os: [macOS-latest] + python-version: ['3.11'] steps: - uses: actions/checkout@v3 @@ -35,7 +35,7 @@ jobs: ulimit -a - name: Create Psi4 Environment - uses: conda-incubator/setup-miniconda@v2 + uses: conda-incubator/setup-miniconda@v3 with: miniconda-version: "latest" activate-environment: p4env diff --git a/magpy/__init__.py b/magpy/__init__.py index 6113a74..9bd621d 100644 --- a/magpy/__init__.py +++ b/magpy/__init__.py @@ -13,6 +13,7 @@ from .hessian import Hessian from .apt import APT from .normal import normal +from .aat import AAT #from ._version import __version__ diff --git a/magpy/aat.py b/magpy/aat.py new file mode 100644 index 0000000..394daa0 --- /dev/null +++ b/magpy/aat.py @@ -0,0 +1,585 @@ +if __name__ == "__main__": + raise Exception("This file cannot be invoked on its own.") + +import psi4 +import magpy +import numpy as np +from .utils import * + + +class AAT(object): + + def __init__(self, molecule, charge=0, spin=1): + + # Ensure geometry remains fixed in space + molecule.fix_orientation(True) + molecule.fix_com(True) + molecule.update_geometry() + molecule.reinterpret_coordentry(False) + self.molecule = molecule + + self.charge = charge + self.spin = spin + + + def compute(self, method='HF', R_disp=0.0001, B_disp=0.0001, **kwargs): + + valid_methods = ['HF', 'CID'] + method = method.upper() + if method not in valid_methods: + raise Exception(f"{method:s} is not an allowed choice of method.") + self.method = method + + valid_normalizations = ['FULL', 'INTERMEDIATE'] + normalization = kwargs.pop('normalization', 'FULL').upper() + if normalization not in valid_normalizations: + raise Exception(f"{normalization:s} is not an allowed choice of normalization.") + self.normalization = normalization + + valid_orbitals = ['SPIN', 'SPATIAL'] + orbitals = kwargs.pop('orbitals', 'SPATIAL').upper() + if orbitals not in valid_orbitals: + raise Exception(f"{orbitals:s} is not an allowed choice of orbital representation.") + self.orbitals = orbitals + + # Extract kwargs + e_conv = kwargs.pop('e_conv', 1e-10) + r_conv = kwargs.pop('r_conv', 1e-10) + maxiter = kwargs.pop('maxiter', 400) + max_diis = kwargs.pop('max_diis', 8) + start_diis = kwargs.pop('start_diis', 1) + print_level = kwargs.pop('print_level', 0) + + mol = self.molecule + + # Compute the unperturbed HF wfn + H = magpy.Hamiltonian(mol) + scf0 = magpy.hfwfn(H, self.charge, self.spin) + scf0.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + if print_level > 0: + print("Psi4 SCF = ", self.run_psi4_scf(H.molecule)) + if method == 'CID': + if orbitals == 'SPATIAL': + ci0 = magpy.ciwfn(scf0, normalization=normalization) + else: + ci0 = magpy.ciwfn_so(scf0, normalization=normalization) + + # Magnetic field displacements + B_pos = [] + B_neg = [] + for B in range(3): + strength = np.zeros(3) + + # +B displacement + if print_level > 0: + print("B(%d)+ Displacement" % (B)) + strength[B] = B_disp + H = magpy.Hamiltonian(mol) + H.add_field(field='magnetic-dipole', strength=strength) + scf = magpy.hfwfn(H, self.charge, self.spin) + scf.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + scf.match_phase(scf0) + if method == 'HF': + B_pos.append(scf) + elif method == 'CID': + if orbitals == 'SPATIAL': + ci = magpy.ciwfn(scf, normalization=normalization) + else: + ci = magpy.ciwfn_so(scf, normalization=normalization) + ci.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + B_pos.append(ci) + + # -B displacement + if print_level > 0: + print("B(%d)- Displacement" % (B)) + strength[B] = -B_disp + H = magpy.Hamiltonian(mol) + H.add_field(field='magnetic-dipole', strength=strength) + scf = magpy.hfwfn(H, self.charge, self.spin) + scf.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + scf.match_phase(scf0) + if method == 'HF': + B_neg.append(scf) + elif method == 'CID': + if orbitals == 'SPATIAL': + ci = magpy.ciwfn(scf, normalization=normalization) + else: + ci = magpy.ciwfn_so(scf, normalization=normalization) + ci.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + B_neg.append(ci) + + # Atomic coordinate displacements + R_pos = [] + R_neg = [] + for R in range(3*mol.natom()): + + # +R displacement + if print_level > 0: + print("R(%d)+ Displacement" % (R)) + H = magpy.Hamiltonian(shift_geom(mol, R, R_disp)) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + scf = magpy.hfwfn(H, self.charge, self.spin) + scf.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + if print_level > 0: + print("Psi4 SCF = ", self.run_psi4_scf(H.molecule)) + scf.match_phase(scf0) + if method == 'HF': + R_pos.append(scf) + elif method == 'CID': + if orbitals == 'SPATIAL': + ci = magpy.ciwfn(scf, normalization=normalization) + else: + ci = magpy.ciwfn_so(scf, normalization=normalization) + ci.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + R_pos.append(ci) + + # -R displacement + if print_level > 0: + print("R(%d)- Displacement" % (R)) + H = magpy.Hamiltonian(shift_geom(mol, R, -R_disp)) + scf = magpy.hfwfn(H, self.charge, self.spin) + scf.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + if print_level > 0: + print("Psi4 SCF = ", self.run_psi4_scf(H.molecule)) + scf.match_phase(scf0) + if method == 'HF': + R_neg.append(scf) + elif method == 'CID': + if orbitals == 'SPATIAL': + ci = magpy.ciwfn(scf, normalization=normalization) + else: + ci = magpy.ciwfn_so(scf, normalization=normalization) + ci.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) + R_neg.append(ci) + + # Compute full MO overlap matrix for all combinations of perturbed MOs + S = [[[0 for k in range(4)] for j in range(3)] for i in range(3*mol.natom())] # list of overlap matrices + for R in range(3*mol.natom()): + if method == 'HF': + R_pos_C = R_pos[R].C + R_neg_C = R_neg[R].C + R_pos_H = R_pos[R].H.basisset + R_neg_H = R_neg[R].H.basisset + else: + R_pos_C = R_pos[R].hfwfn.C + R_neg_C = R_neg[R].hfwfn.C + R_pos_H = R_pos[R].hfwfn.H.basisset + R_neg_H = R_neg[R].hfwfn.H.basisset + + for B in range(3): + if method == 'HF': + B_pos_C = B_pos[B].C + B_neg_C = B_neg[B].C + B_pos_H = B_pos[B].H.basisset + B_neg_H = B_neg[B].H.basisset + else: + B_pos_C = B_pos[B].hfwfn.C + B_neg_C = B_neg[B].hfwfn.C + B_pos_H = B_pos[B].hfwfn.H.basisset + B_neg_H = B_neg[B].hfwfn.H.basisset + + S[R][B][0] = self.mo_overlap(R_pos_C, R_pos_H, B_pos_C, B_pos_H) + S[R][B][1] = self.mo_overlap(R_pos_C, R_pos_H, B_neg_C, B_neg_H) + S[R][B][2] = self.mo_overlap(R_neg_C, R_neg_H, B_pos_C, B_pos_H) + S[R][B][3] = self.mo_overlap(R_neg_C, R_neg_H, B_neg_C, B_neg_H) + + # Compute AAT components using finite-difference + if method == 'HF': + o = slice(0,scf0.ndocc) + elif method == 'CID': + o = slice(0,ci0.no) + no = ci0.no + nv = ci0.nv + + # + AAT_00 = np.zeros((3*mol.natom(), 3)) + for R in range(3*mol.natom()): + if method == 'HF': + C0_R_pos = C0_R_neg = 1.0 + else: + C0_R_pos = R_pos[R].C0 + C0_R_neg = R_neg[R].C0 + + for B in range(3): + if method == 'HF': + C0_B_pos = C0_B_neg = 1.0 + else: + C0_B_pos = B_pos[B].C0 + C0_B_neg = B_neg[B].C0 + + if method == 'HF': + pp = np.linalg.det(S[R][B][0][o,o]) + pm = np.linalg.det(S[R][B][1][o,o]) + mp = np.linalg.det(S[R][B][2][o,o]) + mm = np.linalg.det(S[R][B][3][o,o]) + AAT_00[R,B] = 2*(((pp - pm - mp + mm)/(4*R_disp*B_disp))).imag + else: + pp = self.det_overlap([0], [0], S[R][B][0], o, spins='AAAA') * C0_R_pos * C0_B_pos + pm = self.det_overlap([0], [0], S[R][B][1], o, spins='AAAA') * C0_R_pos * C0_B_neg + mp = self.det_overlap([0], [0], S[R][B][2], o, spins='AAAA') * C0_R_neg * C0_B_pos + mm = self.det_overlap([0], [0], S[R][B][3], o, spins='AAAA') * C0_R_neg * C0_B_neg + AAT_00[R,B] = (((pp - pm - mp + mm)/(4*R_disp*B_disp))).imag + + if method == 'HF': + return AAT_00 + + AAT_0D = np.zeros((3*mol.natom(), 3)) + AAT_D0 = np.zeros((3*mol.natom(), 3)) + for R in range(3*mol.natom()): + ci_R_pos = R_pos[R] + ci_R_neg = R_neg[R] + + for B in range(3): + ci_B_pos = B_pos[B] + ci_B_neg = B_neg[B] + + # + pp, pm, mp, mm = self.AAT_0D(ci_R_pos, ci_R_neg, ci_B_pos, ci_B_neg, S[R][B], o) + AAT_0D[R,B] = (((pp - pm - mp + mm)/(4*R_disp*B_disp))).imag + + #
+ pp, pm, mp, mm = self.AAT_D0(ci_R_pos, ci_R_neg, ci_B_pos, ci_B_neg, S[R][B], o) + AAT_D0[R,B] = (((pp - pm - mp + mm)/(4*R_disp*B_disp))).imag + + AAT_DD = np.zeros((3*mol.natom(), 3)) + for R in range(3*mol.natom()): + ci_R_pos = R_pos[R] + ci_R_neg = R_neg[R] + + for B in range(3): + ci_B_pos = B_pos[B] + ci_B_neg = B_neg[B] + + #
+ pp, pm, mp, mm = self.AAT_DD(ci_R_pos, ci_R_neg, ci_B_pos, ci_B_neg, S[R][B], o) + AAT_DD[R,B] = (((pp - pm - mp + mm)/(4*R_disp*B_disp))).imag + + return AAT_00, AAT_0D, AAT_D0, AAT_DD + + def mo_overlap(self, bra, bra_basis, ket, ket_basis): + """ + Compute the MO overlap matrix between two (possibly different) basis sets + + Parameters + ---------- + bra: MO coefficient matrix for the bra state (NumPy array) + bra_basis: Psi4 BasisSet object for the bra state + ket: MO coefficient matrix for the ket state (NumPy array) + ket_basis: Psi4 BasisSet object for the ket state + + Returns + ------- + S: MO-basis overlap matrix (NumPy array) + """ + # Sanity check + if (bra.shape[0] != ket.shape[0]) or (bra.shape[1] != ket.shape[1]): + raise Exception("Bra and Ket States do not have the same dimensions: (%d,%d) vs. (%d,%d)." % + (bra.shape[0], bra.shape[1], ket.shape[0], ket.shape[1])) + + # Get AO-basis overlap integrals + mints = psi4.core.MintsHelper(bra_basis) + if bra_basis == ket_basis: + S_ao = mints.ao_overlap().np + else: + S_ao = mints.ao_overlap(bra_basis, ket_basis).np + + # Transform to MO basis + S_mo = bra.T @ S_ao @ ket + + # Convert to spin orbitals + if self.orbitals == 'SPIN': + n = 2 * bra.shape[1] + S = np.zeros((n,n), dtype=S_mo.dtype) + for p in range(n): + for q in range(n): + S[p,q] = S_mo[p//2,q//2] * (p%2 == q%2) + return S + else: + return S_mo + + + + # Compute overlap between two determinants in (possibly) different bases + def det_overlap(self, bra_indices, ket_indices, S, o, spins='AAAA'): + """ + Compute the overlap between two Slater determinants (represented by strings of indices) + of equal length in (possibly) different basis sets using the determinant of their overlap. + + Parameters + ---------- + bra_indices: list of substitution indices + ket_indices: list of substitution indices + S: MO overlap between bra and ket bases (NumPy array) + o: Slice of S needed for determinant + spins: 'AAAA', 'AAAB', 'ABAA', or 'ABAB' (string) + """ + + if self.orbitals == 'SPIN': + S = S.copy() + + if len(bra_indices) == 4: # double excitation + i = bra_indices[0]; a = bra_indices[1] + j = bra_indices[2]; b = bra_indices[3] + S[[a,i],:] = S[[i,a],:] + S[[b,j],:] = S[[j,b],:] + + if len(ket_indices) == 4: # double excitation + i = ket_indices[0]; a = ket_indices[1] + j = ket_indices[2]; b = ket_indices[3] + S[:,[a,i]] = S[:,[i,a]] + S[:,[b,j]] = S[:,[j,b]] + + return np.linalg.det(S[o,o]) + + elif self.orbitals == 'SPATIAL': + S_alpha = S.copy() + S_beta = S.copy() + + if len(spins) != 4: + raise Exception(f"Excitations currently limited to doubles only: {spins:s}") + + bra_spin = spins[0] + spins[1] + ket_spin = spins[2] + spins[3] + + if len(bra_indices) == 4: # double excitation + i = bra_indices[0]; a = bra_indices[1] + j = bra_indices[2]; b = bra_indices[3] + if bra_spin == 'AA': + S_alpha[[a,i],:] = S_alpha[[i,a],:] + S_alpha[[b,j],:] = S_alpha[[j,b],:] + elif bra_spin == 'AB': + S_alpha[[a,i],:] = S_alpha[[i,a],:] + S_beta[[b,j],:] = S_beta[[j,b],:] + elif bra_spin == 'BB': + S_beta[[a,i],:] = S_beta[[i,a],:] + beta[[b,j],:] = S_beta[[j,b],:] + + if len(ket_indices) == 4: # double excitation + i = ket_indices[0]; a = ket_indices[1] + j = ket_indices[2]; b = ket_indices[3] + if ket_spin == 'AA': + S_alpha[:,[a,i]] = S_alpha[:,[i,a]] + S_alpha[:,[b,j]] = S_alpha[:,[j,b]] + elif ket_spin == 'AB': + S_alpha[:,[a,i]] = S_alpha[:,[i,a]] + S_beta[:,[b,j]] = S_beta[:,[j,b]] + elif bra_spin == 'BB': + S_beta[[a,i],:] = S_beta[[i,a],:] + beta[[b,j],:] = S_beta[[j,b],:] + + return np.linalg.det(S_alpha[o,o])*np.linalg.det(S_beta[o,o]) + else: + raise Exception("{orbitals:s} is not an allowed choice of orbital representation.") + + + def AAT_0D(self, ci_R_pos, ci_R_neg, ci_B_pos, ci_B_neg, S, o): + no = ci_R_pos.no + nv = ci_R_pos.nv + + pp = pm = mp = mm = 0.0 + if self.orbitals == 'SPATIAL': + for i in range(no): + for a in range(nv): + for j in range(no): + for b in range(nv): + det_AA = self.det_overlap([0], [i, a+no, j, b+no], S[0], o, spins='AAAA') + det_AB = self.det_overlap([0], [i, a+no, j, b+no], S[0], o, spins='ABAB') + pp += (0.5 * (ci_B_pos.C2[i,j,a,b] - ci_B_pos.C2[i,j,b,a]) * det_AA + ci_B_pos.C2[i,j,a,b] * det_AB) * ci_R_pos.C0 + + det_AA = self.det_overlap([0], [i, a+no, j, b+no], S[1], o, spins='AAAA') + det_AB = self.det_overlap([0], [i, a+no, j, b+no], S[1], o, spins='ABAB') + pm += (0.5 * (ci_B_neg.C2[i,j,a,b] - ci_B_neg.C2[i,j,b,a]) * det_AA + ci_B_neg.C2[i,j,a,b] * det_AB) * ci_R_pos.C0 + + det_AA = self.det_overlap([0], [i, a+no, j, b+no], S[2], o, spins='AAAA') + det_AB = self.det_overlap([0], [i, a+no, j, b+no], S[2], o, spins='ABAB') + mp += (0.5 * (ci_B_pos.C2[i,j,a,b] - ci_B_pos.C2[i,j,b,a]) * det_AA + ci_B_pos.C2[i,j,a,b] * det_AB) * ci_R_neg.C0 + + det_AA = self.det_overlap([0], [i, a+no, j, b+no], S[3], o, spins='AAAA') + det_AB = self.det_overlap([0], [i, a+no, j, b+no], S[3], o, spins='ABAB') + mm += (0.5 * (ci_B_neg.C2[i,j,a,b] - ci_B_neg.C2[i,j,b,a]) * det_AA + ci_B_neg.C2[i,j,a,b] * det_AB) * ci_R_neg.C0 + + elif self.orbitals == 'SPIN': + for i in range(0, no, 1): + for a in range(0, nv, 1): + for j in range(0, no, 1): + for b in range(0, nv, 1): + det = self.det_overlap([0], [i, a+no, j, b+no], S[0], o) + pp += 0.25 * ci_B_pos.C2[i,j,a,b] * det * ci_R_pos.C0 + det = self.det_overlap([0], [i, a+no, j, b+no], S[1], o) + pm += 0.25 * ci_B_neg.C2[i,j,a,b] * det * ci_R_pos.C0 + det = self.det_overlap([0], [i, a+no, j, b+no], S[2], o) + mp += 0.25 * ci_B_pos.C2[i,j,a,b] * det * ci_R_neg.C0 + det = self.det_overlap([0], [i, a+no, j, b+no], S[3], o) + mm += 0.25 * ci_B_neg.C2[i,j,a,b] * det * ci_R_neg.C0 + + return pp, pm, mp, mm + + def AAT_D0(self, ci_R_pos, ci_R_neg, ci_B_pos, ci_B_neg, S, o): + no = ci_R_pos.no + nv = ci_R_pos.nv + + pp = pm = mp = mm = 0.0 + if self.orbitals == 'SPATIAL': + for i in range(no): + for a in range(nv): + for j in range(no): + for b in range(nv): + det_AA = self.det_overlap([i, a+no, j, b+no], [0], S[0], o, spins='AAAA') + det_AB = self.det_overlap([i, a+no, j, b+no], [0], S[0], o, spins='ABAB') + pp += (0.5 * (ci_R_pos.C2[i,j,a,b] - ci_R_pos.C2[i,j,b,a]) * det_AA + ci_R_pos.C2[i,j,a,b] * det_AB) * ci_B_pos.C0 + + det_AA = self.det_overlap([i, a+no, j, b+no], [0], S[1], o, spins='AAAA') + det_AB = self.det_overlap([i, a+no, j, b+no], [0], S[1], o, spins='ABAB') + pm += (0.5 * (ci_R_pos.C2[i,j,a,b] - ci_R_pos.C2[i,j,b,a]) * det_AA + ci_R_pos.C2[i,j,a,b] * det_AB) * ci_B_neg.C0 + + det_AA = self.det_overlap([i, a+no, j, b+no], [0], S[2], o, spins='AAAA') + det_AB = self.det_overlap([i, a+no, j, b+no], [0], S[2], o, spins='ABAB') + mp += (0.5 * (ci_R_neg.C2[i,j,a,b] - ci_R_neg.C2[i,j,b,a]) * det_AA + ci_R_neg.C2[i,j,a,b] * det_AB) * ci_B_pos.C0 + + det_AA = self.det_overlap([i, a+no, j, b+no], [0], S[3], o, spins='AAAA') + det_AB = self.det_overlap([i, a+no, j, b+no], [0], S[3], o, spins='ABAB') + mm += (0.5 * (ci_R_neg.C2[i,j,a,b] - ci_R_neg.C2[i,j,b,a]) * det_AA + ci_R_neg.C2[i,j,a,b] * det_AB) * ci_B_neg.C0 + + elif self.orbitals == 'SPIN': + for i in range(0, no, 1): + for a in range(0, nv, 1): + for j in range(0, no, 1): + for b in range(0, nv, 1): + det = self.det_overlap([i, a+no, j, b+no], [0], S[0], o) + pp += 0.25 * ci_R_pos.C2[i,j,a,b] * det * ci_B_pos.C0 + det = self.det_overlap([i, a+no, j, b+no], [0], S[1], o) + pm += 0.25 * ci_R_pos.C2[i,j,a,b] * det * ci_B_neg.C0 + det = self.det_overlap([i, a+no, j, b+no], [0], S[2], o) + mp += 0.25 * ci_R_neg.C2[i,j,a,b] * det * ci_B_pos.C0 + det = self.det_overlap([i, a+no, j, b+no], [0], S[3], o) + mm += 0.25 * ci_R_neg.C2[i,j,a,b] * det * ci_B_neg.C0 + + return pp, pm, mp, mm + + def AAT_DD(self, ci_R_pos, ci_R_neg, ci_B_pos, ci_B_neg, S, o): + no = ci_R_pos.no + nv = ci_R_pos.nv + + pp = pm = mp = mm = 0.0 + if self.orbitals == 'SPATIAL': + for i in range(no): + for a in range(nv): + for j in range(no): + for b in range(nv): + for k in range(no): + for c in range(nv): + for l in range(no): + for d in range(nv): + ci_R = ci_R_pos; ci_B = ci_B_pos; disp = 0 + det_AA_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAA') + det_AA_BB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AABB') + det_AB_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAB') + det_AA_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAB') + det_AB_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAA') + pp += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_AA + pp += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_BB + pp += (1/2) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * ci_B.C2[k,l,c,d] *det_AA_AB + pp += (1/2) * ci_R.C2[i,j,a,b] * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AB_AA + pp += ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] *det_AB_AB + + ci_R = ci_R_pos; ci_B = ci_B_neg; disp = 1 + det_AA_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAA') + det_AA_BB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AABB') + det_AB_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAB') + det_AA_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAB') + det_AB_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAA') + pm += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_AA + pm += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_BB + pm += (1/2) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * ci_B.C2[k,l,c,d] *det_AA_AB + pm += (1/2) * ci_R.C2[i,j,a,b] * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AB_AA + pm += ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] *det_AB_AB + + ci_R = ci_R_neg; ci_B = ci_B_pos; disp = 2 + det_AA_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAA') + det_AA_BB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AABB') + det_AB_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAB') + det_AA_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAB') + det_AB_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAA') + mp += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_AA + mp += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_BB + mp += (1/2) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * ci_B.C2[k,l,c,d] *det_AA_AB + mp += (1/2) * ci_R.C2[i,j,a,b] * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AB_AA + mp += ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] *det_AB_AB + + ci_R = ci_R_neg; ci_B = ci_B_neg; disp = 3 + det_AA_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAA') + det_AA_BB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AABB') + det_AB_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAB') + det_AA_AB = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='AAAB') + det_AB_AA = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o, spins='ABAA') + mm += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_AA + mm += (1/8) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AA_BB + mm += (1/2) * (ci_R.C2[i,j,a,b] - ci_R.C2[i,j,b,a]) * ci_B.C2[k,l,c,d] *det_AA_AB + mm += (1/2) * ci_R.C2[i,j,a,b] * (ci_B.C2[k,l,c,d] - ci_B.C2[k,l,d,c]) *det_AB_AA + mm += ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] *det_AB_AB + + + elif self.orbitals == 'SPIN': + for i in range(0, no, 1): + for a in range(0, nv, 1): + for j in range(0, no, 1): + for b in range(0, nv, 1): + for k in range(0, no, 1): + for c in range(0, nv, 1): + for l in range(0, no, 1): + for d in range(0, nv, 1): + ci_R = ci_R_pos; ci_B = ci_B_pos; disp = 0 + det = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o) + pp += (1/16) * ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] * det + + ci_R = ci_R_pos; ci_B = ci_B_neg; disp = 1 + det = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o) + pm += (1/16) * ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] * det + + ci_R = ci_R_neg; ci_B = ci_B_pos; disp = 2 + det = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o) + mp += (1/16) * ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] * det + + ci_R = ci_R_neg; ci_B = ci_B_neg; disp = 3 + det = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[disp], o) + mm += (1/16) * ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] * det + + return pp, pm, mp, mm + + def nuclear(self): + """ + Computes the nuclear contribution to the atomic axial tensor (AAT). + + Parameters + ---------- + None + + Returns + ------- + aat_nuc: N*3 x 3 array of nuclear contributions to AAT + """ + + geom, mass, elem, Z, uniq = self.molecule.to_arrays() + natom = self.molecule.natom() + + AAT = np.zeros((natom*3,3)) + for M in range(natom): + for alpha in range(3): # atomic Cartesian coordinate + R = M*3 + alpha + for beta in range(3): # magnetic field coordinate + val = 0.0 + for gamma in range(3): # atomic Cartesian coordinate + AAT[R,beta] += (1/4) * levi([alpha, beta, gamma]) * geom[M, gamma] * Z[M] + + return AAT + + + def run_psi4_scf(self, molecule): + geom = molecule.create_psi4_string_from_molecule() + new_mol = psi4.geometry(geom) + new_mol.fix_orientation(True) + new_mol.fix_com(True) + new_mol.update_geometry() + + return psi4.energy('SCF') + diff --git a/magpy/aat_ci_so.py b/magpy/aat_ci_so.py index 64ec896..6c62138 100644 --- a/magpy/aat_ci_so.py +++ b/magpy/aat_ci_so.py @@ -207,9 +207,6 @@ def compute(self, R_disp, B_disp, e_conv=1e-10, r_conv=1e-10, maxiter=400, max_d det = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[R][B][disp], o) pp += (1/16) * ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] * det -# if abs(this.imag) > 1e-12: -# print("%16.13f %1d %1d %1d %1d <%1s%1s%1s%1s | %1s%1s%1s%1s> %1d %1d %1d %1d" % (this.imag, i,a,j,b, s[i%2],s[a%2],s[j%2],s[b%2],s[k%2],s[c%2],s[l%2],s[d%2], k,c,l,d)) - ci_R = ci_R_pos; ci_B = ci_B_neg; disp = 1 det = self.det_overlap([i, a+no, j, b+no], [k, c+no, l, d+no], S[R][B][disp], o) pm += (1/16) * ci_R.C2[i,j,a,b] * ci_B.C2[k,l,c,d] * det @@ -267,7 +264,7 @@ def mo_overlap(self, bra, bra_basis, ket, ket_basis): return S # Compute overlap between two determinants in (possibly) different bases - def det_overlap(self, bra_indices, ket_indices, S_inp, o): + def det_overlap(self, bra_indices, ket_indices, S, o): """ Compute the overlap between two Slater determinants (represented by strings of indices) of equal length in (possibly) different basis sets using the determinant of their overlap. @@ -280,7 +277,7 @@ def det_overlap(self, bra_indices, ket_indices, S_inp, o): o: Slice of S needed for determinant """ - S = S_inp.copy() + S = S.copy() if len(bra_indices) == 4: # double excitation i = bra_indices[0]; a = bra_indices[1] diff --git a/magpy/apt.py b/magpy/apt.py index f349d02..54b7301 100644 --- a/magpy/apt.py +++ b/magpy/apt.py @@ -8,13 +8,7 @@ class APT(object): - def __init__(self, molecule, charge=0, spin=1, method='HF'): - - valid_methods = ['HF', 'CID'] - method = method.upper() - if method not in valid_methods: - raise Exception(f"{method:s} is not an allowed choice of method.") - self.method = method + def __init__(self, molecule, charge=0, spin=1): # Ensure geometry remains fixed in space molecule.fix_orientation(True) @@ -29,7 +23,21 @@ def __init__(self, molecule, charge=0, spin=1, method='HF'): self.natom = self.geom.shape[0] - def compute(self, R_disp=0.001, F_disp=0.0001, e_conv=1e-10, r_conv=1e-10, maxiter=400, max_diis=8, start_diis=1, print_level=0): + def compute(self, method='HF', R_disp=0.001, F_disp=0.0001, **kwargs): + + valid_methods = ['HF', 'CID'] + method = method.upper() + if method not in valid_methods: + raise Exception(f"{method:s} is not an allowed choice of method.") + self.method = method + + # Extract kwargs + e_conv = kwargs.pop('e_conv', 1e-10) + r_conv = kwargs.pop('r_conv', 1e-10) + maxiter = kwargs.pop('maxiter', 400) + max_diis = kwargs.pop('max_diis', 8) + start_diis = kwargs.pop('start_diis', 1) + print_level = kwargs.pop('print_level', 0) params = [e_conv, r_conv, maxiter, max_diis, start_diis, print_level] @@ -64,37 +72,31 @@ def dipole(self, M, alpha, R_disp, F_disp, params): start_diis = params[4] print_level = params[5] - mu = np.zeros((3)) strength = np.eye(3) * F_disp for beta in range(3): H = magpy.Hamiltonian(shift_geom(self.molecule, M*3+alpha, R_disp)) H.add_field(field='electric-dipole', strength=strength[beta]) - scf = magpy.hfwfn(H, self.charge, self.spin, print_level) - escf, C = scf.solve_scf(e_conv, r_conv, maxiter, max_diis, start_diis) - -# if print_level > 0: -# print(f"{M:d}, {alpha:d}, {beta:d} ::: {R_disp:0.5f}; {F_disp:0.5f}") -# print(H.molecule.geometry().np) -# print(f"ESCF = {escf:18.15f}") + scf = magpy.hfwfn(H, self.charge, self.spin) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) if self.method == 'HF': E_pos = escf elif self.method == 'CID': ci = magpy.ciwfn(scf) - eci, C0, C2 = ci.solve_cid(e_conv, r_conv, maxiter, max_diis, start_diis) + eci, C0, C2 = ci.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) E_pos = eci + escf H = magpy.Hamiltonian(shift_geom(self.molecule, M*3+alpha, R_disp)) H.add_field(field='electric-dipole', strength=-1.0*strength[beta]) - scf = magpy.hfwfn(H, self.charge, self.spin, print_level) - escf, C = scf.solve_scf(e_conv, r_conv, maxiter, max_diis, start_diis) + scf = magpy.hfwfn(H, self.charge, self.spin) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) if self.method == 'HF': E_neg = escf elif self.method == 'CID': ci = magpy.ciwfn(scf) - eci, C0, C2 = ci.solve_cid(e_conv, r_conv, maxiter, max_diis, start_diis) + eci, C0, C2 = ci.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) E_neg = eci + escf mu[beta] = -(E_pos - E_neg)/(2 * F_disp) diff --git a/magpy/ciwfn.py b/magpy/ciwfn.py index edfc5e5..c5d5c96 100644 --- a/magpy/ciwfn.py +++ b/magpy/ciwfn.py @@ -14,8 +14,6 @@ def __init__(self, hfwfn, **kwargs): self.hfwfn = hfwfn - self.print_level = hfwfn.print_level - nfzc = hfwfn.H.basisset.n_frozen_core() valid_normalizations = ['FULL', 'INTERMEDIATE'] @@ -28,9 +26,6 @@ def __init__(self, hfwfn, **kwargs): no = self.no = hfwfn.ndocc - nfzc nv = self.nv = hfwfn.nbf - self.no - nfzc - if self.print_level > 0: - print("\nNMO = %d; NACT = %d; NO = %d; NV = %d" % (hfwfn.nbf, self.nt, self.no, self.nv)) - # Set up orbital subspace slices o = self.o = slice(0, no) v = self.v = slice(no, nt) @@ -43,13 +38,13 @@ def __init__(self, hfwfn, **kwargs): # If there are frozen core orbitals, build and add the frozen-core operator # (core contribution to Fock operator) to the one-electron Hamiltonian - efzc = 0 + self.efzc = 0 if nfzc > 0: C = self.hfwfn.C[:,:nfzc] # only core MOs Pc = contract('pi,qi->pq', C.conj(), C) ERI = self.hfwfn.H.ERI hc = h + 2.0 * contract('pqrs,pq->rs', ERI, Pc) - contract('pqrs,ps->qr', ERI, Pc) - efzc = contract('pq,pq->', (h+hc), Pc) + self.efzc = contract('pq,pq->', (h+hc), Pc) h = hc # Select active MOs @@ -73,15 +68,6 @@ def __init__(self, hfwfn, **kwargs): # Build MO-basis Fock matrix (diagonal for canonical MOs, but we don't assume them) F = self.F = self.h + contract('pmqm->pq', L[a,o,a,o]) - # SCF check - ESCF = efzc + 2.0 * contract('ii->',self.h[o,o]) + contract('ijij->', L[o,o,o,o]) - if self.print_level > 0: - print("ESCF (electronic) = ", ESCF) - print("ESCF (total) = ", ESCF+self.hfwfn.H.enuc) - print("HFWFN ESCF (electronic) = ", self.hfwfn.escf) - print("HFWFN ESCF (total) = ", self.hfwfn.escf + self.hfwfn.enuc) - self.E0 = self.hfwfn.escf + self.hfwfn.H.enuc - # Build orbital energy denominators eps_occ = np.diag(F)[o] eps_vir = np.diag(F)[v] @@ -90,19 +76,38 @@ def __init__(self, hfwfn, **kwargs): self.Dijab = Dijab - def solve_cid(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis=1): + def solve(self, **kwargs): + + # Extract kwargs + e_conv = kwargs.pop('e_conv', 1e-7) + r_conv = kwargs.pop('r_conv', 1e-7) + maxiter = kwargs.pop('maxiter', 100) + max_diis = kwargs.pop('max_diis', 8) + start_diis = kwargs.pop('start_diis', 1) + print_level = kwargs.pop('print_level', 0) + + if print_level > 0: + print("\nNMO = %d; NACT = %d; NO = %d; NV = %d" % (self.hfwfn.nbf, self.nt, self.no, self.nv)) o = self.o v = self.v no = self.no nv = self.nv - E0 = self.E0 F = self.F ERI = self.ERI L = self.L Dijab = self.Dijab - if self.print_level > 0: + # SCF check + ESCF = self.efzc + 2.0 * contract('ii->',self.h[o,o]) + contract('ijij->', L[o,o,o,o]) + E0 = self.hfwfn.escf + self.hfwfn.H.enuc + if print_level > 0: + print("\nESCF (electronic) = ", ESCF) + print("ESCF (total) = ", ESCF+self.hfwfn.H.enuc) + print("HFWFN ESCF (electronic) = ", self.hfwfn.escf) + print("HFWFN ESCF (total) = ", self.hfwfn.escf + self.hfwfn.enuc) + + if print_level > 0: print("\nSolving projected CID equations.") # initial guess amplitudes @@ -115,7 +120,7 @@ def solve_cid(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_dii # Setup DIIS object diis = DIIS(C2, max_diis) - if self.print_level > 0: + if print_level > 0: print("CID Iter %3d: CID Ecorr = %.15f dE = %.5E MP2" % (0, eci, -eci)) ediff = eci @@ -134,11 +139,11 @@ def solve_cid(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_dii ediff = eci - eci_last - if self.print_level > 0: + if print_level > 0: print('CID Iter %3d: CID Ecorr = %.15f dE = %.5E rms = %.5E' % (niter, eci, ediff, rms)) if ((abs(ediff) < e_conv) and (abs(rms) < r_conv)): - if self.print_level > 0: + if print_level > 0: print("\nCID Equations converged.") print("CID Correlation Energy = ", eci) print("CID Total Energy = ", eci + E0) @@ -147,7 +152,7 @@ def solve_cid(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_dii if self.normalization == 'FULL': C0, C2 = self.normalize(o, v, C2) norm = np.sqrt(C0*C0 + 2.0 * contract('ijab,ijab->', C2.conj(), C2) - contract('ijab,ijba', C2.conj(), C2)) - if self.print_level > 0: + if print_level > 0: print(f"Normalization check = {norm:18.12f}") self.C0 = C0 diff --git a/magpy/ciwfn_so.py b/magpy/ciwfn_so.py index 3d77170..aa50a11 100644 --- a/magpy/ciwfn_so.py +++ b/magpy/ciwfn_so.py @@ -14,8 +14,6 @@ def __init__(self, hfwfn, **kwargs): self.hfwfn = hfwfn - self.print_level = hfwfn.print_level - nfzc = hfwfn.H.basisset.n_frozen_core() valid_normalizations = ['FULL', 'INTERMEDIATE'] @@ -31,13 +29,13 @@ def __init__(self, hfwfn, **kwargs): # If there are frozen core orbitals, build and add the frozen-core operator # (core contribution to Fock operator) to the one-electron Hamiltonian - efzc = 0 + self.efzc = 0 if nfzc > 0: C = self.hfwfn.C[:,:nfzc] # only core MOs Pc = contract('pi,qi->pq', C.conj(), C) ERI = self.hfwfn.H.ERI hc = h + 2.0 * contract('pqrs,pq->rs', ERI, Pc) - contract('pqrs,ps->qr', ERI, Pc) - efzc = contract('pq,pq->', (h+hc), Pc) + self.efzc = contract('pq,pq->', (h+hc), Pc) h = hc # Select active MOs @@ -61,9 +59,6 @@ def __init__(self, hfwfn, **kwargs): no = self.no = 2*(hfwfn.ndocc - nfzc) nv = self.nv = self.nt - self.no - if self.print_level > 0: - print("\nNMO = %d; NACT = %d; NO = %d; NV = %d" % (hfwfn.nbf, self.nt, self.no, self.nv)) - # Set up orbital subspace slices o = self.o = slice(0, no) v = self.v = slice(no, nt) @@ -88,15 +83,6 @@ def __init__(self, hfwfn, **kwargs): # Build MO-basis Fock matrix (diagonal for canonical MOs, but we don't assume that) F = self.F = self.h + contract('pmqm->pq', ERI[a,o,a,o]) - # SCF check - ESCF = contract('ii->',self.h[o,o]) + 0.5 * contract('ijij->', ERI[o,o,o,o]) - if self.print_level > 0: - print("ESCF (electronic) = ", ESCF) - print("ESCF (total) = ", ESCF+self.hfwfn.H.enuc) - print("HFWFN ESCF (electronic) = ", self.hfwfn.escf) - print("HFWFN ESCF (total) = ", self.hfwfn.escf + self.hfwfn.enuc) - self.E0 = self.hfwfn.escf + self.hfwfn.H.enuc - # Build orbital energy denominators eps_occ = np.diag(F)[o] eps_vir = np.diag(F)[v] @@ -105,17 +91,39 @@ def __init__(self, hfwfn, **kwargs): self.Dijab = Dijab - def solve(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis=1): + def solve(self, **kwargs): + + # Extract kwargs + e_conv = kwargs.pop('e_conv', 1e-7) + r_conv = kwargs.pop('r_conv', 1e-7) + maxiter = kwargs.pop('maxiter', 100) + max_diis = kwargs.pop('max_diis', 8) + start_diis = kwargs.pop('start_diis', 1) + print_level = kwargs.pop('print_level', 0) + + if print_level > 0: + print("\nNMO = %d; NACT = %d; NO = %d; NV = %d" % (self.hfwfn.nbf, self.nt, self.no, self.nv)) o = self.o v = self.v no = self.no nv = self.nv - E0 = self.E0 F = self.F ERI = self.ERI Dijab = self.Dijab + # SCF check + ESCF = contract('ii->',self.h[o,o]) + 0.5 * contract('ijij->', ERI[o,o,o,o]) + E0 = self.hfwfn.escf + self.hfwfn.H.enuc + if print_level > 0: + print("ESCF (electronic) = ", ESCF) + print("ESCF (total) = ", ESCF+self.hfwfn.H.enuc) + print("HFWFN ESCF (electronic) = ", self.hfwfn.escf) + print("HFWFN ESCF (total) = ", self.hfwfn.escf + self.hfwfn.enuc) + + if print_level > 0: + print("\nSolving projected CID equations.") + # initial guess amplitudes -- intermediate normalization C0 = 1.0 C2 = ERI[o,o,v,v]/Dijab @@ -126,7 +134,7 @@ def solve(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis=1) # Setup DIIS object diis = DIIS(C2, max_diis) - if self.print_level > 0: + if print_level > 0: print("CID Iter %3d: CID Ecorr = %.15f dE = %.5E MP2" % (0, eci, -eci)) ediff = eci @@ -144,11 +152,11 @@ def solve(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis=1) ediff = eci - eci_last - if self.print_level > 0: + if print_level > 0: print('CID Iter %3d: CID Ecorr = %.15f dE = %.5E rms = %.5E' % (niter, eci, ediff, rms)) if ((abs(ediff) < e_conv) and (abs(rms) < r_conv)): - if self.print_level > 0: + if print_level > 0: print("\nCID Equations converged.") print("CID Correlation Energy = ", eci) print("CID Total Energy = ", eci + E0) @@ -157,7 +165,8 @@ def solve(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis=1) if self.normalization == 'FULL': C0, C2 = self.normalize(o, v, C2) norm = np.sqrt(C0*C0 + (1/4) * contract('ijab,ijab', C2.conj(), C2)) - print(f"Normalization check = {norm:18.12f}") + if print_level > 0: + print(f"Normalization check = {norm:18.12f}") self.C0 = C0 self.C2 = C2 diff --git a/magpy/hessian.py b/magpy/hessian.py index be5d2f6..e2d3c09 100644 --- a/magpy/hessian.py +++ b/magpy/hessian.py @@ -8,13 +8,7 @@ class Hessian(object): - def __init__(self, molecule, charge=0, spin=1, method='HF'): - - valid_methods = ['HF', 'CID'] - method = method.upper() - if method not in valid_methods: - raise Exception(f"{method:s} is not an allowed choice of method.") - self.method = method + def __init__(self, molecule, charge=0, spin=1): # Ensure geometry remains fixed in space molecule.fix_orientation(True) @@ -29,7 +23,21 @@ def __init__(self, molecule, charge=0, spin=1, method='HF'): self.natom = self.geom.shape[0] - def compute(self, disp=0.001, e_conv=1e-10, r_conv=1e-10, maxiter=400, max_diis=8, start_diis=1, print_level=0): + def compute(self, method='HF', disp=0.001, **kwargs): + + valid_methods = ['HF', 'CID'] + method = method.upper() + if method not in valid_methods: + raise Exception(f"{method:s} is not an allowed choice of method.") + self.method = method + + # Extract kwargs + e_conv = kwargs.pop('e_conv', 1e-10) + r_conv = kwargs.pop('r_conv', 1e-10) + maxiter = kwargs.pop('maxiter', 400) + max_diis = kwargs.pop('max_diis', 8) + start_diis = kwargs.pop('start_diis', 1) + print_level = kwargs.pop('print_level', 0) params = [e_conv, r_conv, maxiter, max_diis, start_diis, print_level] @@ -78,8 +86,8 @@ def energy(self, M1, alpha1, disp1, M2, alpha2, disp2, params): print_level = params[5] H = magpy.Hamiltonian(shift_geom(shift_geom(self.molecule, M1*3+alpha1, disp1), M2*3+alpha2, disp2)) - scf = magpy.hfwfn(H, self.charge, self.spin, print_level) - escf, C = scf.solve_scf(e_conv, r_conv, maxiter, max_diis, start_diis) + scf = magpy.hfwfn(H, self.charge, self.spin) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) if print_level > 0: print(f"{M1:d}, {alpha1:d}; {M2:d}, {alpha2:d} ::: {disp1:0.5f}; {disp2:0.5f}") print(H.molecule.geometry().np) @@ -89,5 +97,5 @@ def energy(self, M1, alpha1, disp1, M2, alpha2, disp2, params): return escf elif self.method == 'CID': ci = magpy.ciwfn(scf) - eci, C0, C2 = ci.solve_cid(e_conv, r_conv, maxiter, max_diis, start_diis) + eci, C0, C2 = ci.solve(e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) return eci + escf diff --git a/magpy/hfwfn.py b/magpy/hfwfn.py index b492cc5..85937c9 100644 --- a/magpy/hfwfn.py +++ b/magpy/hfwfn.py @@ -9,7 +9,7 @@ class hfwfn(object): - def __init__(self, H, charge=0, spin=1, print_level=0): + def __init__(self, H, charge=0, spin=1): # Keep the Hamiltonian (including molecule and basisset) self.H = H @@ -23,7 +23,6 @@ def __init__(self, H, charge=0, spin=1, print_level=0): # Determine number of orbitals self.nbf = H.basisset.nbf() - self.print_level = print_level def nelectron(self, charge): nelec = -charge @@ -32,9 +31,16 @@ def nelectron(self, charge): return nelec - def solve_scf(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis=1): - print_level = self.print_level + def solve(self, **kwargs): + + # Extract kwargs + e_conv = kwargs.pop('e_conv', 1e-7) + r_conv = kwargs.pop('r_conv', 1e-7) + maxiter = kwargs.pop('maxiter', 100) + max_diis = kwargs.pop('max_diis', 8) + start_diis = kwargs.pop('start_diis', 1) + print_level = kwargs.pop('print_level', 0) # Electronic Hamiltonian, including fields H = self.H diff --git a/magpy/normal.py b/magpy/normal.py index b8824e3..2f0cb8d 100644 --- a/magpy/normal.py +++ b/magpy/normal.py @@ -1,10 +1,25 @@ import psi4 import magpy import numpy as np -from .utils import AAT_nuc from opt_einsum import contract -def normal(molecule, method='HF', read_hessian=False, **kwargs): +def normal(molecule, method='HF', r_disp=0.001, f_disp=0.0001, b_disp=0.0001, **kwargs): + + valid_methods = ['HF', 'CID'] + method = method.upper() + if method not in valid_methods: + raise Exception(f"{method:s} is not an allowed choice of method.") + + # Extract kwargs + e_conv = kwargs.pop('e_conv', 1e-13) + r_conv = kwargs.pop('r_conv', 1e-13) + maxiter = kwargs.pop('maxiter', 400) + max_diis = kwargs.pop('max_diis', 8) + start_diis = kwargs.pop('start_diis', 1) + print_level = kwargs.pop('print_level', 0) + read_hessian = kwargs.pop('read_hessian', False) + if read_hessian == True: + fcm_file = kwargs.pop('fcm_file', 'fcm') # Physical constants and a few derived units _c = psi4.qcel.constants.get("speed of light in vacuum") # m/s @@ -32,17 +47,11 @@ def normal(molecule, method='HF', read_hessian=False, **kwargs): # Compute the Hessian [Eh/(a0^2)] if read_hessian is False: - hessian = magpy.Hessian(molecule, 0, 1, method) - disp = 0.001 - e_conv = 1e-13 - r_conv = 1e-13 - maxiter = 400 - max_diis = 8 - start_diis = 1 - print_level = 0 - H = hessian.compute(disp, e_conv, r_conv, maxiter, max_diis, start_diis, print_level) + hessian = magpy.Hessian(molecule) + H = hessian.compute(method, r_disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) else: - H = np.genfromtxt(kwargs.pop('file'), skip_header=1).reshape(3*molecule.natom(),3*molecule.natom()) + print("Using provided hessian...") + H = np.genfromtxt(fcm_file, skip_header=1).reshape(3*molecule.natom(),3*molecule.natom()) # Mass-weight the Hessian (Eh/(a0^2 m_e)) masses = np.array([(molecule.mass(i)*_u/_me) for i in range(molecule.natom())]) @@ -57,17 +66,12 @@ def normal(molecule, method='HF', read_hessian=False, **kwargs): S = np.flip(Lx, 1)[:,:(3*molecule.natom()-6)] freq = np.flip(np.sqrt(w))[:(3*molecule.natom()-6)] - for i in range(3*molecule.natom()-6): # Assuming non-linear molecules for now - print(f"{freq[i]*conv_freq_au2wavenumber:7.2f}") + #for i in range(3*molecule.natom()-6): # Assuming non-linear molecules for now + # print(f"{freq[i]*conv_freq_au2wavenumber:7.2f}") # Compute APTs and transform to normal mode basis - APT = magpy.APT(molecule, 0, 1, method) - r_disp = 0.001 - f_disp = 0.0001 - e_conv = 1e-13 - r_conv = 1e-13 - maxiter = 400 - P = APT.compute(r_disp, f_disp, e_conv, r_conv, maxiter) # 3N x 3 + APT = magpy.APT(molecule) + P = APT.compute(method, r_disp, f_disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) # (e a0)/(a0 sqrt(m_e)) P = P.T @ S # 3 x (3N-6) @@ -76,20 +80,18 @@ def normal(molecule, method='HF', read_hessian=False, **kwargs): for i in range(3*molecule.natom()-6): ir_intensities[i] = contract('j,j->', P[:,i], P[:,i]) - for i in range(3*molecule.natom()-6): # Assuming non-linear molecules for now - print(f"{freq[i]*conv_freq_au2wavenumber:7.2f} {ir_intensities[i]*conv_ir_au2kmmol:7.3f}") + #for i in range(3*molecule.natom()-6): # Assuming non-linear molecules for now + # print(f"{freq[i]*conv_freq_au2wavenumber:7.2f} {ir_intensities[i]*conv_ir_au2kmmol:7.3f}") # Compute AATs and transform to normal mode basis - r_disp = 0.0001 - b_disp = 0.0001 + r_disp = 0.0001 # need smaller displacement for AAT + AAT = magpy.AAT(molecule) if method == 'HF': - AAT = magpy.AAT_HF(molecule) - I = AAT.compute(r_disp, b_disp) # electronic contribution + I = AAT.compute('HF', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) elif method == 'CID': - AAT = magpy.AAT_CI(molecule) - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp) # electronic contribution + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) I = I_00 + I_DD - J = AAT_nuc(molecule) # nuclear contribution + J = AAT.nuclear() # nuclear contribution M = I + J # 3N x 3 M = M.T @ S # 3 x (3N-6) @@ -102,5 +104,7 @@ def normal(molecule, method='HF', read_hessian=False, **kwargs): print(" (cm-1) (km/mol) (esu**2 cm**2 10**44)") print("----------------------------------------------") for i in range(3*molecule.natom()-6): # Assuming non-linear molecules for now - print(f" {freq[i]*conv_freq_au2wavenumber:7.2f} {ir_intensities[i]*conv_ir_au2kmmol:8.3f} {rotatory_strengths[i] * conv_vcd_au2cgs:8.3f}") + print(f" {freq[i]*conv_freq_au2wavenumber:7.2f} {ir_intensities[i]*conv_ir_au2kmmol:8.3f} {rotatory_strengths[i]*conv_vcd_au2cgs:8.3f}") + + return freq*conv_freq_au2wavenumber, ir_intensities*conv_ir_au2kmmol, rotatory_strengths*conv_vcd_au2cgs diff --git a/magpy/tests/test_002_electric_dipole.py b/magpy/tests/test_002_electric_dipole.py index 4ddb804..aadec40 100644 --- a/magpy/tests/test_002_electric_dipole.py +++ b/magpy/tests/test_002_electric_dipole.py @@ -34,18 +34,15 @@ def test_ccsd_hf(): # Get MagPy energies and dipole moment H = magpy.Hamiltonian(mol) scf = magpy.hfwfn(H) - e_conv = 1e-13 - r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=1e-13, r_conv=1e-13) A = 0.0001 H.add_field(field='electric-dipole', strength=np.array([0.0, 0.0, A])) - escf_pos, C_pos = scf.solve_scf(e_conv, r_conv, maxiter) + escf_pos, C_pos = scf.solve(e_conv=1e-13, r_conv=1e-13) H.reset_V() H.add_field(field='electric-dipole', strength=np.array([0.0, 0.0, -A])) - escf_neg, C_neg = scf.solve_scf(e_conv, r_conv, maxiter) + escf_neg, C_neg = scf.solve(e_conv=1e-13, r_conv=1e-13) mu = -(escf_pos - escf_neg)/(2 * A) diff --git a/magpy/tests/test_003_AAT_HF.py b/magpy/tests/test_003_AAT_HF.py index 5873586..2f69111 100644 --- a/magpy/tests/test_003_AAT_HF.py +++ b/magpy/tests/test_003_AAT_HF.py @@ -21,11 +21,8 @@ def test_AAT_HF_H2O(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_HF(mol) - - r_disp = 0.0001 - b_disp = 0.0001 - I = AAT.compute(r_disp, b_disp) + AAT = magpy.AAT(mol) + I = AAT.compute() print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print(I) @@ -61,10 +58,8 @@ def test_AAT_HF_H2O2(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_HF(mol) - r_disp = 0.0001 - b_disp = 0.0001 - I = AAT.compute(r_disp, b_disp) + AAT = magpy.AAT(mol) + I = AAT.compute() print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print(I) @@ -105,10 +100,8 @@ def test_AAT_HF_etho(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_HF(mol) - r_disp = 0.0001 - b_disp = 0.0001 - I = AAT.compute(r_disp, b_disp) + AAT = magpy.AAT(mol) + I = AAT.compute() print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print(I) diff --git a/magpy/tests/test_004_CID.py b/magpy/tests/test_004_CID.py index 272867a..435a031 100644 --- a/magpy/tests/test_004_CID.py +++ b/magpy/tests/test_004_CID.py @@ -23,14 +23,13 @@ def test_CID_H2O_STO3G(): psi4.set_options({'freeze_core': 'false'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn(scf, normalization = 'intermediate') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -74.94207992819220 c4ci = -0.06865825074438 C0_ref = 1.0 @@ -39,7 +38,7 @@ def test_CID_H2O_STO3G(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn(scf, normalization = 'full') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -74.94207992819220 c4ci = -0.06865825074438 C0_ref = 0.978084764153443 @@ -49,14 +48,13 @@ def test_CID_H2O_STO3G(): psi4.set_options({'freeze_core': 'true'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn(scf, normalization = 'intermediate') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) C0_ref = 1.0 c4scf = -74.94207992819220 c4ci = -0.06859643905558 @@ -65,7 +63,7 @@ def test_CID_H2O_STO3G(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn(scf, normalization = 'full') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) C0_ref = 0.9780778514625926 c4scf = -74.94207992819220 c4ci = -0.06859643905558 @@ -90,14 +88,13 @@ def test_CID_H2O_CCPVDZ(): psi4.set_options({'freeze_core': 'false'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn(scf, normalization='intermediate') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) C0_ref = 1.0 c4scf = -75.98979581991861 c4ci = -0.21279410950205 @@ -106,7 +103,7 @@ def test_CID_H2O_CCPVDZ(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn(scf, normalization='full') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) C0_ref = 0.9712657483891116 c4scf = -75.98979581991861 c4ci = -0.21279410950205 @@ -116,14 +113,13 @@ def test_CID_H2O_CCPVDZ(): psi4.set_options({'freeze_core': 'true'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn(scf, normalization='intermediate') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) C0_ref = 1.0 c4scf = -75.98979581991861 c4ci = -0.21098966441656 @@ -132,7 +128,7 @@ def test_CID_H2O_CCPVDZ(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn(scf, normalization='full') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) C0_ref = 0.9712377677176858 c4scf = -75.98979581991861 c4ci = -0.21098966441656 @@ -165,10 +161,9 @@ def test_CID_H2O2_STO3G(): psi4.set_options({'freeze_core': 'false'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn(scf, normalization='intermediate') - eci, C0, C2 = cid.solve_cid(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) diff --git a/magpy/tests/test_005_AAT_CID.py b/magpy/tests/test_005_AAT_CID.py index 54437aa..a15291c 100644 --- a/magpy/tests/test_005_AAT_CID.py +++ b/magpy/tests/test_005_AAT_CID.py @@ -22,13 +22,13 @@ def test_AAT_CID_H2DIMER(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI(mol, 0, 1, 1, 'intermediate') + AAT = magpy.AAT(mol, 0, 1) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='intermediate') print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print("Hartree-Fock component:") print(I_00) @@ -120,13 +120,14 @@ def test_AAT_CID_H2DIMER_NORM(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI(mol, 0, 1, 1, 'full') + AAT = magpy.AAT(mol, 0, 1) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='full') + print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print("Hartree-Fock component:") print(I_00) @@ -218,13 +219,13 @@ def test_AAT_CID_H2O(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI(mol, 0, 1, 1, 'intermediate') + AAT = magpy.AAT(mol, 0, 1) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='intermediate') print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print("Hartree-Fock component:") print(I_00) @@ -274,15 +275,15 @@ def test_AAT_CID_H2O(): """) I_DD_ref = make_np_array(""" -[[-0. 0.00000000000005 -0.006571145180452] - [-0.000000000000113 0.000000000000231 0.000000000000325] - [ 0.036050121268829 -0.000000000000415 0. ] - [ 0.000000000000012 -0.000000000000017 0.001267365504146] - [ 0.000000000000069 0.000000000000001 -0.009999089093936] - [-0.020642645989297 0.016765832724568 -0.00000000000002 ] - [ 0.000000000000011 -0.000000000000008 0.001267365503931] - [-0.000000000000043 -0.000000000000052 0.009999089093942] - [-0.020642645989002 -0.016765832724794 0.00000000000002 ]] +[[ 0.000000000000076 0.000000000000007 -0.006571145180512] + [ 0.000000000002677 0.000000000001015 -0.000000000000175] + [ 0.03605012126928 -0.00000000000062 0.000000000000001] + [ 0.000000000000042 -0.000000000000022 0.001267365504038] + [ 0.000000000000008 0.000000000000042 -0.009999089093697] + [-0.020642645988524 0.016765832724844 -0.000000000000034] + [-0.000000000000056 -0.000000000000001 0.001267365504546] + [-0.000000000000143 0.000000000000071 0.009999089093786] + [-0.020642645988882 -0.016765832725199 0.000000000000003]] """) assert(np.max(np.abs(I_00_ref-I_00)) < 1e-9) @@ -304,13 +305,14 @@ def test_AAT_CID_H2O_NORM(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI(mol, 0, 1, 1, 'full') + AAT = magpy.AAT(mol, 0, 1) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='full') + print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print("Hartree-Fock component:") print(I_00) @@ -360,14 +362,14 @@ def test_AAT_CID_H2O_NORM(): """) I_DD_ref = make_np_array(""" -[[ 0.000000000000072 0.000000000000007 -0.006282596322971] - [ 0.000000000002559 0.00000000000097 -0.000000000000137] +[[ 0.000000000000072 0.000000000000007 -0.006282596323029] + [ 0.000000000002559 0.00000000000097 -0.000000000000168] [ 0.034467106277752 -0.000000000000593 0.000000000000001] - [ 0.00000000000004 -0.000000000000021 0.001211713581436] - [ 0.000000000000008 0.00000000000004 -0.00956001407078 ] + [ 0.00000000000004 -0.000000000000021 0.001211713581405] + [ 0.000000000000008 0.00000000000004 -0.009560014070841] [-0.019736196388032 0.016029619821933 -0.000000000000033] [-0.000000000000054 -0.000000000000001 0.001211713581891] - [-0.000000000000137 0.000000000000068 0.009560014070954] + [-0.000000000000137 0.000000000000068 0.009560014070926] [-0.019736196388375 -0.016029619822273 0.000000000000003]] """) diff --git a/magpy/tests/test_006_CID_SO.py b/magpy/tests/test_006_CID_SO.py index 79c6f5e..ec0a9ff 100644 --- a/magpy/tests/test_006_CID_SO.py +++ b/magpy/tests/test_006_CID_SO.py @@ -23,14 +23,13 @@ def test_CID_SO_H2O_STO3G(): psi4.set_options({'freeze_core': 'false'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn_so(scf, normalization='intermediate') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -74.94207992819220 c4ci = -0.06865825074438 C0_ref = 1.0 @@ -39,7 +38,7 @@ def test_CID_SO_H2O_STO3G(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn_so(scf, normalization='full') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -74.94207992819220 c4ci = -0.06865825074438 C0_ref = 0.9780847641534441 @@ -50,14 +49,13 @@ def test_CID_SO_H2O_STO3G(): # Test frozen core psi4.set_options({'freeze_core': 'true'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn_so(scf, normalization='intermediate') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -74.94207992819220 c4ci = -0.06859643905558 C0_ref = 1.0 @@ -66,7 +64,7 @@ def test_CID_SO_H2O_STO3G(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn_so(scf, normalization='full') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -74.94207992819220 c4ci = -0.06859643905558 C0_ref = 0.9780778514625926 @@ -91,14 +89,13 @@ def test_CID_SO_H2O_CCPVDZ(): psi4.set_options({'freeze_core': 'false'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn_so(scf, normalization='intermediate') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -75.98979581991861 c4ci = -0.21279410950205 C0_ref = 1.0 @@ -107,7 +104,7 @@ def test_CID_SO_H2O_CCPVDZ(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn_so(scf, normalization='full') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -75.98979581991861 c4ci = -0.21279410950205 C0_ref = 0.9712657483891116 @@ -117,14 +114,13 @@ def test_CID_SO_H2O_CCPVDZ(): psi4.set_options({'freeze_core': 'true'}) H = magpy.Hamiltonian(mol) - scf = magpy.hfwfn(H, 0, 1, 1) + scf = magpy.hfwfn(H, 0, 1) e_conv = 1e-13 r_conv = 1e-13 - maxiter = 100 - escf, C = scf.solve_scf(e_conv, r_conv, maxiter) + escf, C = scf.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) cid = magpy.ciwfn_so(scf, normalization='intermediate') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -75.98979581991861 c4ci = -0.21098966441656 C0_ref = 1.0 @@ -133,7 +129,7 @@ def test_CID_SO_H2O_CCPVDZ(): assert(abs(C0 - C0_ref) < 1e-11) cid = magpy.ciwfn_so(scf, normalization='full') - eci, C0, C2 = cid.solve(e_conv, r_conv, maxiter) + eci, C0, C2 = cid.solve(e_conv=e_conv, r_conv=r_conv, print_level=1) c4scf = -75.98979581991861 c4ci = -0.21098966441656 C0_ref = 0.9712377677176858 diff --git a/magpy/tests/test_007_AAT_CID_SO.py b/magpy/tests/test_007_AAT_CID_SO.py index 180ee12..43f1768 100644 --- a/magpy/tests/test_007_AAT_CID_SO.py +++ b/magpy/tests/test_007_AAT_CID_SO.py @@ -22,13 +22,13 @@ def test_AAT_CID_H2DIMER(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI_SO(mol, 0, 1, 1, 'intermediate') + AAT = magpy.AAT(mol, 0, 1) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='intermediate', orbitals='spin') print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print("Hartree-Fock component:") print(I_00) @@ -120,13 +120,13 @@ def test_AAT_CID_H2DIMER_NORM(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI_SO(mol, 0, 1, 1, 'full') + AAT = magpy.AAT(mol, 0, 1) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='full', orbitals='spin') print("\nElectronic Contribution to Atomic Axial Tensor (a.u.):") print("Hartree-Fock component:") print(I_00) @@ -204,6 +204,7 @@ def test_AAT_CID_H2DIMER_NORM(): assert(np.max(np.abs(I_D0_ref-I_D0)) < 1e-9) assert(np.max(np.abs(I_DD_ref-I_DD)) < 1e-9) +@pytest.mark.skip(reason="too long") def test_AAT_CID_H2O(): psi4.set_memory('2 GB') @@ -218,13 +219,14 @@ def test_AAT_CID_H2O(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI_SO(mol, 0, 1, 1, 'intermediate') + AAT = magpy.AAT(mol) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='intermediate', orbitals='spin') + print("Hartree-Fock component:") print(I_00) print("<0|D> Component\n") @@ -289,6 +291,7 @@ def test_AAT_CID_H2O(): assert(np.max(np.abs(I_D0_ref-I_D0)) < 1e-9) assert(np.max(np.abs(I_DD_ref-I_DD)) < 1e-9) +@pytest.mark.skip(reason="too long") def test_AAT_CID_H2O_NORM(): psi4.set_memory('2 GB') @@ -303,13 +306,13 @@ def test_AAT_CID_H2O_NORM(): rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) print(f" SCF Energy from Psi4: {rhf_e}") - AAT = magpy.AAT_CI_SO(mol, 0, 1, 1, 'full') + AAT = magpy.AAT(mol) r_disp = 0.0001 b_disp = 0.0001 e_conv = 1e-13 r_conv = 1e-13 - I_00, I_0D, I_D0, I_DD = AAT.compute(r_disp, b_disp, e_conv, r_conv) + I_00, I_0D, I_D0, I_DD = AAT.compute('CID', r_disp, b_disp, e_conv=e_conv, r_conv=r_conv, normalization='full', orbitals='spin') print("Hartree-Fock component:") print(I_00) print("<0|D> Component\n") diff --git a/magpy/tests/test_008_AAT_nuc.py b/magpy/tests/test_008_AAT_nuc.py index 22fd9a2..f9bd9b3 100644 --- a/magpy/tests/test_008_AAT_nuc.py +++ b/magpy/tests/test_008_AAT_nuc.py @@ -4,7 +4,6 @@ from ..data.molecules import * import numpy as np import os -from ..utils import AAT_nuc np.set_printoptions(precision=15, linewidth=200, threshold=200, suppress=True) @@ -25,9 +24,10 @@ def test_AAT_nuc(): no_com """) - #geom, mass, elem, Z, uniq = mol.to_arrays() - AAT = AAT_nuc(mol) - print(AAT) + AAT = magpy.AAT(mol) + I_nuc = AAT.nuclear() + print(I_nuc) + AAT_ref = np.array([ # Matches Dalton [ 0.00000000000000, -0.19050858260000, -2.63852838000000], [ 0.19050858260000, 0.00000000000000, 0.00000000000000], @@ -43,7 +43,7 @@ def test_AAT_nuc(): [ -0.42102591000000, 0.41162146750000, 0.00000000000000], ]) - assert(np.max(np.abs(AAT-AAT_ref)) < 1e-11) + assert(np.max(np.abs(I_nuc-AAT_ref)) < 1e-11) ## Ethylene Oxide mol = psi4.geometry(""" @@ -60,9 +60,10 @@ def test_AAT_nuc(): units bohr """) - #geom, mass, elem, Z, uniq = mol.to_arrays() - AAT = AAT_nuc(mol) - print(AAT) + AAT = magpy.AAT(mol) + I_nuc = AAT.nuclear() + print(I_nuc) + AAT_ref = np.array([ # Matches Dalton [ 0.00000000, 3.22387278, -0.00000000], [ -3.22387278, 0.00000000, 0.00000000], @@ -87,4 +88,4 @@ def test_AAT_nuc(): [ 0.59486586, 0.43724415, 0.00000000], ]) - assert(np.max(np.abs(AAT-AAT_ref)) < 1e-8) + assert(np.max(np.abs(I_nuc-AAT_ref)) < 1e-8) diff --git a/magpy/tests/test_009_hessian_HF.py b/magpy/tests/test_009_hessian_HF.py index 9865f5f..7e00b9b 100644 --- a/magpy/tests/test_009_hessian_HF.py +++ b/magpy/tests/test_009_hessian_HF.py @@ -28,7 +28,7 @@ def test_Hessian_H2O_STO3G(): units bohr """) - hessian = magpy.Hessian(mol, 0, 1, 'HF') + hessian = magpy.Hessian(mol, 0, 1) disp = 0.001 e_conv = 1e-14 r_conv = 1e-14 @@ -36,7 +36,7 @@ def test_Hessian_H2O_STO3G(): max_diis=8 start_diis=1 print_level=1 - hess = hessian.compute(disp, e_conv, r_conv, maxiter, max_diis, start_diis, print_level) + hess = hessian.compute('HF', disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) # CFOUR HF/STO-3G analytic Hessian hess_ref = np.array([ diff --git a/magpy/tests/test_010_hessian_CID.py b/magpy/tests/test_010_hessian_CID.py index 0441842..7e409f8 100644 --- a/magpy/tests/test_010_hessian_CID.py +++ b/magpy/tests/test_010_hessian_CID.py @@ -28,15 +28,15 @@ def test_Hessian_H2O_STO3G(): units bohr """) - hessian = magpy.Hessian(mol, 0, 1, 'CID') + hessian = magpy.Hessian(mol) disp = 0.001 - e_conv = 1e-14 - r_conv = 1e-14 + e_conv = 1e-13 + r_conv = 1e-13 maxiter = 400 max_diis=8 start_diis=1 print_level=1 - hess = hessian.compute(disp, e_conv, r_conv, maxiter, max_diis, start_diis, print_level) + hess = hessian.compute('CID', disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) # CFOUR CID/STO-3G finite-difference Hessian hess_ref = np.array([ diff --git a/magpy/tests/test_011_APT_HF.py b/magpy/tests/test_011_APT_HF.py index 98d34ea..3a72409 100644 --- a/magpy/tests/test_011_APT_HF.py +++ b/magpy/tests/test_011_APT_HF.py @@ -29,7 +29,7 @@ def test_APT_H2O_STO3G(): units bohr """) - apt = magpy.APT(mol, 0, 1, 'HF') + apt = magpy.APT(mol) R_disp = 0.001 F_disp = 0.0001 e_conv = 1e-14 @@ -38,7 +38,7 @@ def test_APT_H2O_STO3G(): max_diis=8 start_diis=1 print_level=1 - dipder = apt.compute(R_disp, F_disp, e_conv, r_conv, maxiter, max_diis, start_diis, print_level) + dipder = apt.compute('HF', R_disp, F_disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) print(dipder) # CFOUR HF/STO-3G analytic dipole derivatives @@ -54,5 +54,5 @@ def test_APT_H2O_STO3G(): [ 0. 0.2216402091 -0.0417659999]] """) - assert(np.max(np.abs(dipder-dipder_ref)) < 1e-6) + assert(np.max(np.abs(dipder-dipder_ref)) < 1e-5) diff --git a/magpy/tests/test_012_APT_CID.py b/magpy/tests/test_012_APT_CID.py index 0b1026a..2dd7a80 100644 --- a/magpy/tests/test_012_APT_CID.py +++ b/magpy/tests/test_012_APT_CID.py @@ -29,7 +29,7 @@ def test_APT_H2O_STO3G(): units bohr """) - apt = magpy.APT(mol, 0, 1, 'CID') + apt = magpy.APT(mol) R_disp = 0.001 F_disp = 0.0001 e_conv = 1e-14 @@ -38,7 +38,7 @@ def test_APT_H2O_STO3G(): max_diis=8 start_diis=1 print_level=1 - dipder = apt.compute(R_disp, F_disp, e_conv, r_conv, maxiter, max_diis, start_diis, print_level) + dipder = apt.compute('CID', R_disp, F_disp, e_conv=e_conv, r_conv=r_conv, maxiter=maxiter, max_diis=max_diis, start_diis=start_diis, print_level=print_level) # Gaussian09 CID/STO-3G finite-difference dipole derivatives dipder_ref = np.array([ @@ -53,5 +53,5 @@ def test_APT_H2O_STO3G(): [ 4.06138511E-13, -2.54923176E-01, -9.94990576E-02], ]) - assert(np.max(np.abs(dipder-dipder_ref)) < 1e-6) + assert(np.max(np.abs(dipder-dipder_ref)) < 1e-5) diff --git a/magpy/tests/test_013_VCD_HF.py b/magpy/tests/test_013_VCD_HF.py index b2c6fc1..2578d2b 100644 --- a/magpy/tests/test_013_VCD_HF.py +++ b/magpy/tests/test_013_VCD_HF.py @@ -7,7 +7,7 @@ np.set_printoptions(precision=10, linewidth=200, threshold=200, suppress=True) -def test_VCD_H2Dimer_STO3G(): +def test_VCD_H2O2_STO3GP(): psi4.set_memory('2 GB') psi4.set_output_file('output.dat', False) @@ -17,44 +17,7 @@ def test_VCD_H2Dimer_STO3G(): 'r_convergence': 1e-13}) psi4.set_options({'basis': 'STO-3G'}) - mol = psi4.geometry(moldict["(H2)_2"]) - magpy.normal(mol, 'HF', True, file="fcm_H2dimer_HF_STO3G.txt") - -def test_VCD_H2O2_STO3G(): - - psi4.set_memory('2 GB') - psi4.set_output_file('output.dat', False) - psi4.set_options({'scf_type': 'pk', - 'e_convergence': 1e-13, - 'd_convergence': 1e-13, - 'r_convergence': 1e-13}) - - psi4.set_options({'basis': 'STO-3G'}) - # HF/STO-3G Optimized geometry from CFOUR - mol = psi4.geometry(""" -H -1.838051419951917 1.472647809969243 0.806472638463773 -O -1.312852628446968 -0.129910361247786 -0.050815108044519 -O 1.312852628446968 0.129910361247786 -0.050815108044519 -H 1.838051419951917 -1.472647809969243 0.806472638463773 -no_com -no_reorient -symmetry c1 -units bohr - """) - - magpy.normal(mol, 'HF', True, file="fcm_H2O2_HF_STO3G.txt") - -def test_VCD_H2O2_STO3G_KP(): - - psi4.set_memory('2 GB') - psi4.set_output_file('output.dat', False) - psi4.set_options({'scf_type': 'pk', - 'e_convergence': 1e-13, - 'd_convergence': 1e-13, - 'r_convergence': 1e-13}) - - psi4.set_options({'basis': 'STO-3G'}) - # + # CFOUR HF/STO-3G optimized geometry mol = psi4.geometry(""" O 0.0000000000 1.3192641900 -0.0952542913 O -0.0000000000 -1.3192641900 -0.0952542913 @@ -66,6 +29,16 @@ def test_VCD_H2O2_STO3G_KP(): units bohr """) - magpy.normal(mol, 'HF') + freq, ir, vcd = magpy.normal(mol, 'HF') + + # CFOUR + freq_ref = np.array([4148.2766, 4140.8893, 1781.0495, 1589.6394, 1486.9510, 184.6327]) + # CFOUR + ir_ref = np.array([30.3023, 12.7136, 2.0726, 46.0798, 0.0168, 134.1975]) + # DALTON (*-1 because it appears they're rotatory strengths have the wrong sign) + vcd_ref = np.array([50.538, -53.528, 28.607, -14.650, 1.098, 101.378]) + assert(np.max(np.abs(freq-freq_ref)) < 0.1) + assert(np.max(np.abs(ir-ir_ref)) < 0.1) + assert(np.max(np.abs(vcd-vcd_ref)) < 0.1) diff --git a/magpy/tests/test_014_VCD_CI.py b/magpy/tests/test_014_VCD_CI.py index fed3c1b..f87059c 100644 --- a/magpy/tests/test_014_VCD_CI.py +++ b/magpy/tests/test_014_VCD_CI.py @@ -7,6 +7,7 @@ np.set_printoptions(precision=10, linewidth=200, threshold=200, suppress=True) +@pytest.mark.skip(reason="not ready") def test_VCD_H2Dimer_STO3G(): psi4.set_memory('2 GB') @@ -18,9 +19,10 @@ def test_VCD_H2Dimer_STO3G(): psi4.set_options({'basis': 'STO-3G'}) mol = psi4.geometry(moldict["(H2)_2"]) - magpy.normal(mol, 'CID', True, file="fcm_H2dimer_CID_STO3G.txt") + magpy.normal(mol, 'CID', read_hessian=True, fcm_file="fcm_H2dimer_CID_STO3G.txt") +@pytest.mark.skip(reason="not ready") def test_VCD_H2O2_STO3G(): psi4.set_memory('2 GB') @@ -43,8 +45,10 @@ def test_VCD_H2O2_STO3G(): units bohr """) - magpy.normal(mol, 'CID', True, file="fcm_H2O2_CID_STO3G.txt") + magpy.normal(mol, 'CID', read_hessian=True, fcm_file="fcm_H2O2_CID_STO3G.txt") + +@pytest.mark.skip(reason="not ready") def test_VCD_H2O2_STO3G_KP(): psi4.set_memory('2 GB') @@ -68,5 +72,3 @@ def test_VCD_H2O2_STO3G_KP(): """) magpy.normal(mol, 'CID') - - diff --git a/magpy/utils.py b/magpy/utils.py index f284573..08c362b 100644 --- a/magpy/utils.py +++ b/magpy/utils.py @@ -31,33 +31,6 @@ def levi(indexes): indexes[j], indexes[j + 1] = indexes[j + 1], indexes[j] return -1 * levi(indexes) -def AAT_nuc(molecule): - """ - Computes the nuclear contribution to the atomic axial tensor (AAT). - - Parameters - ---------- - molecule: Psi4 Molecule object - - Returns - ------- - aat_nuc: N*3 x 3 array of nuclear contributions to AAT - """ - geom, mass, elem, Z, uniq = molecule.to_arrays() - natom = len(Z) - - AAT = np.zeros((natom*3,3)) - for M in range(natom): - for alpha in range(3): # atomic Cartesian coordinate - R = M*3 + alpha - for beta in range(3): # magnetic field coordinate - val = 0.0 - for gamma in range(3): # atomic Cartesian coordinate - AAT[R,beta] += (1/4) * levi([alpha, beta, gamma]) * geom[M, gamma] * Z[M] - - return AAT - - def shift_geom(molecule, R, R_disp): """ Shift the R-th Cartesian coordinate of the given molecule geometry by