From 54a37246ae548420bb0c3fc5ba48df7fd3a52d24 Mon Sep 17 00:00:00 2001 From: Chaoqun Zhang Date: Tue, 20 Jan 2026 19:45:09 -0500 Subject: [PATCH 1/3] add GHF-based pprpa and gradient for rpprpa --- examples/grad/05-rpprpa_grad.py | 95 ++ examples/pprpa/06-gpprpa_excitation_energy.py | 48 + pyscf/grad/rpprpa.py | 879 ++++++++++++++++++ pyscf/pprpa/gpprpa_davidson.py | 418 +++++++++ pyscf/pprpa/rpprpa_davidson.py | 14 +- pyscf/pprpa/rpprpa_direct.py | 18 +- 6 files changed, 1458 insertions(+), 14 deletions(-) create mode 100644 examples/grad/05-rpprpa_grad.py create mode 100644 examples/pprpa/06-gpprpa_excitation_energy.py create mode 100644 pyscf/grad/rpprpa.py create mode 100644 pyscf/pprpa/gpprpa_davidson.py diff --git a/examples/grad/05-rpprpa_grad.py b/examples/grad/05-rpprpa_grad.py new file mode 100644 index 000000000..3525e0536 --- /dev/null +++ b/examples/grad/05-rpprpa_grad.py @@ -0,0 +1,95 @@ +r""" + This is an example to calculate the molecular geometrical gradient of the ppRPA energy. +""" +from pyscf import gto, scf, dft +from pyscf.pprpa import rpprpa_davidson, rpprpa_direct +from pyscf.grad import rpprpa +import numpy as np + +mult = "s" +istate = 0 +nocc = 5 +nvir = 10 + +def mfobj(dx): + mol = gto.Mole() + mol.verbose = 0 + mol.atom = [["O", (0.0, 0.0, 0.0)], + ["H", (0.0, 1.0, 1.0+dx)], + ["H", (0.0, -1.0, 1.0)]] + mol.basis = "cc-pvdz" + mol.charge = 0 + mol.unit = "Bohr" + mol.build() + + # density fitting is required to ensure the numerical and analytical gradients are the same + # mf = scf.RHF(mol).density_fit() + + # DFT grid response not implemented in pprpa + # Higher grid level for better accuracy + mf = dft.RKS(mol, xc="b3lyp").density_fit() + mf.grids.level = 7 + mf.conv_tol = 1e-12 + return mf + +def pprpaobj(mf, nocc_act, nvir_act): + mp = rpprpa_davidson.RppRPADavidson(mf, nocc_act, nvir_act, channel="pp", nroot=3) + mp.mu = 0.0 + return mp + + + + +dx = 0.001 +mf1 = mfobj(dx) +mf2 = mfobj(-dx) + +e1 = mf1.kernel() +# f1 = mf1.get_fock(dm=dm0_hf) +h1 = mf1.get_hcore() + +e2 = mf2.kernel() +# f2 = mf2.get_fock(dm=dm0_hf) +h2 = mf2.get_hcore() + +# dfock = (f1 - f2)/(2.0*dx) +dh = (h1 - h2)/(2.0*dx) +e_hf = (e1 - e2)/dx/2.0 + +mp1 = pprpaobj(mf1, nocc, nvir) +mp1.kernel(mult) +mp1.analyze() +if mult == "s": + e1 += mp1.exci_s[istate] if mp1.channel == "pp" else -mp1.exci_s[istate] +else: + e1 += mp1.exci_t[istate] if mp1.channel == "pp" else -mp1.exci_t[istate] +mp2 = pprpaobj(mf2, nocc, nvir) +mp2.kernel(mult) +mp2.analyze() +if mult == "s": + e2 += mp2.exci_s[istate] if mp2.channel == "pp" else -mp2.exci_s[istate] +else: + e2 += mp2.exci_t[istate] if mp2.channel == "pp" else -mp2.exci_t[istate] + +e_rpa = (e1 - e2)/dx/2.0 + +mf = mfobj(0.0) +mf.kernel() +if mult == "s": + oo_dim = nocc * (nocc + 1) // 2 +else: + oo_dim = nocc * (nocc - 1) // 2 +mp = pprpaobj(mf, nocc, nvir) +mp.kernel(mult) +mp.analyze() + +xy = mp.xy_s[istate] if mult == "s" else mp.xy_t[istate] +from lib_pprpa.grad import pprpa +mpg = mp.Gradients(mult,istate) +mpg.kernel() + +print("analytical (Total): ", mpg.de) +print("numerical (Total): ", e_rpa) +print("diff: ", e_rpa - mpg.de[1,2]) +# print(e1,e2) +print(mpg.de[0]+mpg.de[1]+mpg.de[2]) # Should be zero diff --git a/examples/pprpa/06-gpprpa_excitation_energy.py b/examples/pprpa/06-gpprpa_excitation_energy.py new file mode 100644 index 000000000..ed3f9a091 --- /dev/null +++ b/examples/pprpa/06-gpprpa_excitation_energy.py @@ -0,0 +1,48 @@ +from pyscf import gto, dft + +# For ppRPA excitation energy of N-electron system in particle-particle channel +# mean field is (N-2)-electron + +mol = gto.Mole() +mol.verbose = 4 +mol.atom = [ + ["O", (0.0, 0.0, 0.0)], + ["H", (0.0, -0.7571, 0.5861)], + ["H", (0.0, 0.7571, 0.5861)]] +mol.basis = 'def2-svp' +# create a (N-2)-electron system for charged-neutral H2O +mol.charge = 2 +mol.build() + +# =====> Part I. Restricted and its real-valued generalized ppRPA <===== +# restricted KS-DFT calculation as starting point of RppRPA +rmf = dft.RKS(mol) +rmf.xc = "b3lyp" +rmf.kernel() +from pyscf.pprpa.rpprpa_davidson import RppRPADavidson +pp = RppRPADavidson(rmf, nocc_act=3, nvir_act=None, nroot=1) +# solve for singlet states +pp.kernel("s") +# solve for triplet states +pp.kernel("t") +pp.analyze() + +# Davidson algorithm for GKS, equivalent to the above RKS case +from pyscf.pprpa.gpprpa_davidson import GppRPADavidson +gmf = rmf.to_gks() +pp = GppRPADavidson(gmf, nocc_act=6, nvir_act=None, nroot=4) +# solve for singlet and triplet states +pp.kernel() +pp.analyze() + + + +# =====> Part II. Complex-valued generalized ppRPA <===== +gmf = dft.GKS(mol).x2c1e() +gmf.xc = "b3lyp" +gmf.kernel() +pp = GppRPADavidson(gmf, nocc_act=6, nvir_act=None, nroot=4) +# solve for "singlet" and "triplet" states +# triplets are no longer degenerate with SOC +pp.kernel() +pp.analyze() \ No newline at end of file diff --git a/pyscf/grad/rpprpa.py b/pyscf/grad/rpprpa.py new file mode 100644 index 000000000..b3148ebfa --- /dev/null +++ b/pyscf/grad/rpprpa.py @@ -0,0 +1,879 @@ +import numpy as np +from functools import reduce +from pyscf import lib, scf +from pyscf.df.df_jk import _DFHF +from pyscf.lib import logger +from pyscf.df.grad import rhf as rhf_grad +from pyscf.grad import rks as rks_grad + + +def grad_elec(pprpa_grad, xy, mult, atmlst=None): + mf = pprpa_grad.base._scf + pprpa = pprpa_grad.base + mol = mf.mol + mf_grad = mf.nuc_grad_method() + if atmlst is None: + atmlst = range(mol.natm) + assert mult in ['t', 's'], 'mult = {}. is not valid in grad_elec'.format(mult) + + nocc_all = mf.mol.nelectron // 2 + nocc = pprpa.nocc_act + nvir = pprpa.nvir_act + nfrozen_occ = nocc_all - nocc + + hcore_deriv = mf_grad.hcore_generator(mol) + s1 = mf_grad.get_ovlp(mol) + dm0, i_int = make_rdm1_relaxed_rhf_pprpa( + pprpa, mf, xy=xy, mult=mult, cphf_max_cycle=pprpa_grad.cphf_max_cycle, cphf_conv_tol=pprpa_grad.cphf_conv_tol + ) + dm0 = np.einsum('pi,ij,qj->pq', mf.mo_coeff, dm0, mf.mo_coeff, optimize=True) + pprpa_grad.rdm1e = dm0 + dm0_hf = mf.make_rdm1() + i_int = np.einsum('pi,ij,qj->pq', mf.mo_coeff, i_int, mf.mo_coeff, optimize=True) + i_int -= mf_grad.make_rdm1e(mf.mo_energy, mf.mo_coeff, mf.mo_occ) + + occ_y_mat, vir_x_mat = get_xy_full(xy, pprpa.oo_dim, mult) + coeff_occ = mf.mo_coeff[:, nfrozen_occ : nfrozen_occ + nocc] + coeff_vir = mf.mo_coeff[:, nfrozen_occ + nocc : nfrozen_occ + nocc + nvir] + xy_ao = np.einsum('pi,ij,qj->pq', coeff_vir, vir_x_mat, coeff_vir, optimize=True) + np.einsum( + 'pi,ij,qj->pq', coeff_occ, occ_y_mat, coeff_occ, optimize=True + ) + + aux_response = False + if isinstance(mf, _DFHF): + # aux_response is Ture by default in DFHF + # To my opinion, aux_response should always be True for DFHF + aux_response = mf_grad.auxbasis_response + else: + logger.warn(pprpa, + 'The analytical gradient of the ppRPA must be used with the density\n\ + fitting mean-field method. The calculation will proceed but the analytical\n\ + gradient is no longer exact (does NOT agree with numerical gradients).') + + if not hasattr(mf, 'xc'): # RHF + vj, vk = mf_grad.get_jk(mol, (dm0_hf, dm0, xy_ao), hermi=0) + vhf = np.zeros_like(vj) + + vhf[:2] = vj[:2] - 0.5 * vk[:2] + vhf[2] = vk[2] + if aux_response: + vhf_aux = np.zeros_like(vj.aux) + vhf_aux[:2, :2] = vj.aux[:2, :2] - 0.5 * vk.aux[:2, :2] + if mult == 's': + vhf_aux[2, 2] = vk.aux[2, 2] + else: + vhf_aux[2, 2] = -vk.aux[2, 2] + vhf = lib.tag_array(vhf, aux=vhf_aux) + + aoslices = mol.aoslice_by_atom() + de = np.zeros((len(atmlst), 3)) + for k, ia in enumerate(atmlst): + p0, p1 = aoslices[ia, 2:] + h1ao = hcore_deriv(ia) + h1ao[:,p0:p1] += vhf[0,:,p0:p1] + h1ao[:,:,p0:p1] += vhf[0,:,p0:p1].transpose(0,2,1) + de[k] += np.einsum('xij,ij->x', h1ao, dm0+dm0_hf) + # nabla was applied on bra in s1, *2 for the contributions of nabla|ket> + de[k] += np.einsum('xij,ij->x', vhf[1, :, p0:p1], dm0_hf[p0:p1, :]) * 2 + de[k] += np.einsum('xij,ij->x', vhf[2, :, p0:p1], xy_ao[p0:p1, :]) * 2 + + de[k] += np.einsum('xij,ji->x', s1[:, p0:p1], i_int[:, p0:p1]) * 2 + + if aux_response: + de[k] += vhf.aux[0, 1, ia] + 0.5*vhf.aux[0, 0, ia] + de[k] += vhf.aux[1, 0, ia] + 0.5*vhf.aux[0, 0, ia] + de[k] += vhf.aux[2, 2, ia] + else: # RKS + # The grid response by default is not included. + # Even if grid_response is set to True, the response is not complete. + # It will include the response of the Vxc but NOT the fxc. + # For benchmarking, one can use high grid level to avoid this error. + # mf_grad.grid_response = True + + vj, vk = mf_grad.get_jk(mol, xy_ao, hermi=0) + vhf = vk + if aux_response: + vxc, vjk = get_veff_df_rks(mf_grad, mol, (dm0_hf, dm0)) + if mult == 's': + vhf_aux = vk.aux[0, 0] + else: + vhf_aux = -vk.aux[0, 0] + vhf = lib.tag_array(vhf, aux=vhf_aux) + else: + vxc, vjk = get_veff_rks(mf_grad, mol, (dm0_hf, dm0)) + + vjk[1] += _contract_xc_kernel(mf, mf.xc, dm0, None, False, False, True)[0][1:]*0.5 + + aoslices = mol.aoslice_by_atom() + de = np.zeros((len(atmlst), 3)) + for k, ia in enumerate(atmlst): + p0, p1 = aoslices[ia, 2:] + h1ao = hcore_deriv(ia) + h1ao[:, p0:p1] += vxc[0, :, p0:p1] + vjk[0, :, p0:p1] + h1ao[:, :, p0:p1] += vxc[0, :, p0:p1].transpose(0, 2, 1) + vjk[0, :, p0:p1].transpose(0, 2, 1) + de[k] += np.einsum('xij,ij->x', h1ao, dm0 + dm0_hf) + # nabla was applied on bra in s1, *2 for the contributions of nabla|ket> + de[k] += np.einsum('xij,ij->x', vjk[1, :, p0:p1], dm0_hf[p0:p1, :]) * 2 + de[k] += np.einsum('xij,ij->x', vhf[:, p0:p1], xy_ao[p0:p1, :]) * 2 + + de[k] += np.einsum('xij,ji->x', s1[:, p0:p1], i_int[:, p0:p1]) * 2 + + if aux_response: + de[k] += vjk.aux[0, 1, ia] + 0.5*vjk.aux[0, 0, ia] + de[k] += vjk.aux[1, 0, ia] + 0.5*vjk.aux[0, 0, ia] + de[k] += vhf.aux[ia] + if mf_grad.grid_response: + de[k] += vxc.exc1_grid[0, ia] + + return de + + +def get_xy_full(xy, oo_dim, mult='t'): + """Expand the lower triangular xy matrix to the full matrix.""" + vv_dim = len(xy) - oo_dim + if mult == 't': + ndim_v = round(0.5 * (np.sqrt(8 * vv_dim + 1) + 1)) if vv_dim > 0 else 0 + ndim_o = round(0.5 * (np.sqrt(8 * oo_dim + 1) + 1)) if oo_dim > 0 else 0 + occ_y_mat = np.zeros((ndim_o, ndim_o), dtype=xy.dtype) + vir_x_mat = np.zeros((ndim_v, ndim_v), dtype=xy.dtype) + occ_y_mat[np.tril_indices(ndim_o, -1)] = xy[:oo_dim] + vir_x_mat[np.tril_indices(ndim_v, -1)] = xy[oo_dim:] + occ_y_mat = occ_y_mat - occ_y_mat.T + vir_x_mat = vir_x_mat - vir_x_mat.T + else: + ndim_v = round(0.5 * (np.sqrt(8 * vv_dim + 1) - 1)) if vv_dim > 0 else 0 + ndim_o = round(0.5 * (np.sqrt(8 * oo_dim + 1) - 1)) if oo_dim > 0 else 0 + occ_y_mat = np.zeros((ndim_o, ndim_o), dtype=xy.dtype) + vir_x_mat = np.zeros((ndim_v, ndim_v), dtype=xy.dtype) + occ_y_mat[np.tril_indices(ndim_o)] = xy[:oo_dim] + vir_x_mat[np.tril_indices(ndim_v)] = xy[oo_dim:] + occ_y_mat = occ_y_mat + occ_y_mat.T + vir_x_mat = vir_x_mat + vir_x_mat.T + np.fill_diagonal(occ_y_mat, 1.0 / np.sqrt(2.0) * occ_y_mat.diagonal()) + np.fill_diagonal(vir_x_mat, 1.0 / np.sqrt(2.0) * vir_x_mat.diagonal()) + + return occ_y_mat, vir_x_mat + + +def make_rdm1_unrelaxed_from_xy_full(occ_y_mat, vir_x_mat, diag=True): + """Make unrelaxed one-particle density matrix from the full X and Y matrices.""" + if diag: + di = -np.einsum('ij,ij->i', occ_y_mat.conj(), occ_y_mat) + da = np.einsum('ab,ab->a', vir_x_mat.conj(), vir_x_mat) + # combine the two parts + den = np.concatenate((di, da)) + else: + print('Warning: non-diagonal 1e-RDM is not well-defined for pp-RPA with both DEA and DIP blocks.') + den_v = np.einsum('ac,bc->ba', vir_x_mat.conj(), vir_x_mat) + den_o = -np.einsum('ik,jk->ij', occ_y_mat.conj(), occ_y_mat) + den = np.zeros((den_v.shape[0] + den_o.shape[0], den_v.shape[1] + den_o.shape[1]), dtype=den_v.dtype) + den[: den_o.shape[0], : den_o.shape[1]] = den_o + den[den_o.shape[0] :, den_o.shape[1] :] = den_v + return den + + +def make_rdm1_relaxed_rhf_pprpa(pprpa, mf, xy=None, mult='t', istate=0, cphf_max_cycle=20, cphf_conv_tol=1.0e-8): + r"""Calculate relaxed density matrix (and the I intermediates) + for given pprpa and mean-field object. + Args: + pprpa: a pprpa object. + mf: a mean-field RHF/RKS object. + Returns: + den_relaxed: the relaxed one-particle density matrix (nmo_full, nmo_full) + i_int: the I intermediates (nmo_full, nmo_full) + Both are in the MO basis. + """ + assert mult in ['t', 's'], 'mult = {}. is not valid in make_rdm1_relaxed_pprpa'.format(mult) + + if xy is None: + if mult == 's': + xy = pprpa.xy_s[istate] + else: + xy = pprpa.xy_t[istate] + nocc_all = mf.mol.nelectron // 2 + nvir_all = mf.mol.nao - nocc_all + nocc = pprpa.nocc_act + nvir = pprpa.nvir_act + nfrozen_occ = nocc_all - nocc + nfrozen_vir = nvir_all - nvir + if mult == 's': + oo_dim = (nocc + 1) * nocc // 2 + else: + oo_dim = (nocc - 1) * nocc // 2 + + # create slices + slice_p = choose_slice('p', nfrozen_occ, nocc, nvir, nfrozen_vir) # all active + slice_i = choose_slice('i', nfrozen_occ, nocc, nvir, nfrozen_vir) # active occupied + slice_a = choose_slice('a', nfrozen_occ, nocc, nvir, nfrozen_vir) # active virtual + slice_ip = choose_slice('ip', nfrozen_occ, nocc, nvir, nfrozen_vir) # frozen occupied + slice_ap = choose_slice('ap', nfrozen_occ, nocc, nvir, nfrozen_vir) # frozen virtual + slice_I = choose_slice('I', nfrozen_occ, nocc, nvir, nfrozen_vir) # all occupied + slice_A = choose_slice('A', nfrozen_occ, nocc, nvir, nfrozen_vir) # all virtual + + orbA = mf.mo_coeff[:, slice_A] + orbI = mf.mo_coeff[:, slice_I] + orbp = mf.mo_coeff[:, slice_p] + orbi = mf.mo_coeff[:, slice_i] + orba = mf.mo_coeff[:, slice_a] + occ_y_mat, vir_x_mat = get_xy_full(xy, oo_dim, mult) + if nfrozen_occ > 0 or nfrozen_vir > 0: + mo_ene_full = mf.mo_energy + Lpq_full = pprpa.ao2mo(full_mo=True) + else: + mo_ene_full = pprpa.mo_energy + if pprpa.Lpq is not None: + Lpq_full = pprpa.Lpq + elif pprpa.Lpi is not None and pprpa.Lpa is not None: + Lpq_full = np.concatenate((pprpa.Lpi, pprpa.Lpa), axis=2) + else: + raise RuntimeError('Lpq or Lpi/Lpa is required in make_rdm1_relaxed_rhf_pprpa') + + # set singlet=None, generate function for CPHF type response kernel + vresp = mf.gen_response(singlet=None, hermi=1) + den_u = make_rdm1_unrelaxed_from_xy_full(occ_y_mat, vir_x_mat) + den_u_ao = np.einsum('pi,i,qi->pq', orbp, den_u, orbp, optimize=True) + veff_den_u = reduce(np.dot, (mf.mo_coeff.T, vresp(den_u_ao) * 2, mf.mo_coeff)) + + # calculate I' first + i_prime = np.zeros((len(mo_ene_full), len(mo_ene_full)), dtype=occ_y_mat.dtype) + # I' active-active block + i_prime[slice_p, slice_p] += contraction_2rdm_Lpq( + occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, 'p', 'p' + ) + i_prime[slice_a, slice_i] += veff_den_u[slice_a, slice_i] + for p in choose_range('p', nfrozen_occ, nocc, nvir, nfrozen_vir): + i_prime[p, p] += mo_ene_full[p] * den_u[p - nfrozen_occ] + + if nfrozen_vir > 0: + # I' frozen virtual-active block + i_prime[slice_ap, slice_p] += contraction_2rdm_Lpq( + occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, 'ap', 'p' + ) + i_prime[slice_ap, slice_i] += veff_den_u[slice_ap, slice_i] + if nfrozen_occ > 0: + # I' frozen occupied-active block + i_prime[slice_ip, slice_p] += contraction_2rdm_Lpq( + occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, 'ip', 'p' + ) + # I' all virtual-frozen occupied block + i_prime[slice_A, slice_ip] += veff_den_u[slice_A, slice_ip] + + # calculate I'' next + i_prime_prime = np.zeros_like(i_prime) + # I'' active virtual-all occupied block + i_prime_prime[slice_a, slice_I] = i_prime[slice_a, slice_I] - i_prime[slice_I, slice_a].T + # I'' = I' blocks + i_prime_prime[slice_A, slice_a] = i_prime[slice_A, slice_a] + i_prime_prime[slice_I, slice_i] = i_prime[slice_I, slice_i] + i_prime_prime[slice_ap, slice_I] = i_prime[slice_ap, slice_I] + + d_prime = np.zeros_like(i_prime_prime) + threshold = 1.0e-6 + # D' all occupied-active occupied block + for i in choose_range('I', nfrozen_occ, nocc, nvir, nfrozen_vir): + for j in choose_range('i', nfrozen_occ, nocc, nvir, nfrozen_vir): + denorm = mo_ene_full[j] - mo_ene_full[i] + factor = 1.0 / denorm if abs(denorm) >= threshold else 0.0 + d_prime[i, j] = factor * i_prime_prime[i, j] + + # D' all virtual-active virtual block + for a in choose_range('A', nfrozen_occ, nocc, nvir, nfrozen_vir): + for b in choose_range('a', nfrozen_occ, nocc, nvir, nfrozen_vir): + denorm = mo_ene_full[b] - mo_ene_full[a] + factor = 1.0 / denorm if abs(denorm) >= threshold else 0.0 + d_prime[a, b] = factor * i_prime_prime[a, b] + + x_int = i_prime_prime[slice_A, slice_I].copy() + d_ao = reduce(np.dot, (orbI, d_prime[slice_I, slice_i], orbi.T)) + d_ao += reduce(np.dot, (orbA, d_prime[slice_A, slice_a], orba.T)) + d_ao += d_ao.T + x_int += reduce(np.dot, (orbA.T, vresp(d_ao) * 2, orbI)) + + def fvind(x): + dm = reduce(np.dot, (orbA, x.reshape(nvir + nfrozen_vir, nocc + nfrozen_occ) * 2, orbI.T)) + v1ao = vresp(dm + dm.T) + return reduce(np.dot, (orbA.T, v1ao, orbI)).ravel() + + from pyscf.scf import cphf + + d_prime[slice_A, slice_I] = cphf.solve( + fvind, mo_ene_full, mf.mo_occ, x_int, max_cycle=cphf_max_cycle, tol=cphf_conv_tol + )[0].reshape(nvir + nfrozen_vir, nocc + nfrozen_occ) + + i_int = -np.einsum('qp,p->qp', d_prime, mo_ene_full) + # I all occupied-all occupied block + dp_ao = reduce(np.dot, (mf.mo_coeff, d_prime, mf.mo_coeff.T)) + veff_dp_II = reduce(np.dot, (orbI.T, vresp(dp_ao + dp_ao.T), orbI)) + i_int[slice_I, slice_I] -= 0.5 * veff_den_u[slice_I, slice_I] + i_int[slice_I, slice_I] -= veff_dp_II + # I active virtual-all occupied block + i_int[slice_I, slice_a] -= i_prime[slice_I, slice_a] + + # I active-active block extra term + for i in choose_range('p', nfrozen_occ, nocc, nvir, nfrozen_vir): + for j in choose_range('p', nfrozen_occ, nocc, nvir, nfrozen_vir): + denorm = mo_ene_full[j] - mo_ene_full[i] + if abs(denorm) < threshold: + i_int[i, j] -= 0.5 * i_prime[i, j] + + den_relaxed = d_prime + # active-active block + for p in choose_range('p', nfrozen_occ, nocc, nvir, nfrozen_vir): + den_relaxed[p, p] += 0.5 * den_u[p - nfrozen_occ] + den_relaxed = den_relaxed + den_relaxed.T + i_int = i_int + i_int.T + + return den_relaxed, i_int + + +def choose_slice(label, nfrozen_occ, nocc, nvir, nfrozen_vir): + """Choose the slice for the given label. + "i" for active occupied orbitals; + "a" for active virtual orbitals; + "p" for all active orbitals; + "ip" for frozen occupied orbitals; + "ap" for frozen virtual orbitals; + "I" for all occupied orbitals; + "A" for all virtual orbitals; + "P" for all orbitals; + + In the energy ordering, frozen occ -> active occ -> active vir -> frozen vir. + """ + if label == 'i': + return slice(nfrozen_occ, nfrozen_occ + nocc) + elif label == 'a': + return slice(nfrozen_occ + nocc, nfrozen_occ + nocc + nvir) + elif label == 'p': + return slice(nfrozen_occ, nfrozen_occ + nocc + nvir) + elif label == 'ip': + return slice(0, nfrozen_occ) + elif label == 'ap': + return slice(nfrozen_occ + nocc + nvir, nfrozen_occ + nocc + nvir + nfrozen_vir) + elif label == 'I': + return slice(0, nfrozen_occ + nocc) + elif label == 'A': + return slice(nfrozen_occ + nocc, nfrozen_occ + nocc + nvir + nfrozen_vir) + elif label == 'P': + return slice(0, nfrozen_occ + nocc + nvir + nfrozen_vir) + else: + raise ValueError('label = {}. is not valid in choose_slice'.format(label)) + + +def choose_range(label, nfrozen_occ, nocc, nvir, nfrozen_vir): + """Choose the range list for the given label. + "i" for active occupied orbitals; + "a" for active virtual orbitals; + "p" for all active orbitals; + "ip" for frozen occupied orbitals; + "ap" for frozen virtual orbitals; + "I" for all occupied orbitals; + "A" for all virtual orbitals; + "P" for all orbitals; + + In the energy ordering, frozen occ -> active occ -> active vir -> frozen vir. + """ + if label == 'i': + return range(nfrozen_occ, nfrozen_occ + nocc) + elif label == 'a': + return range(nfrozen_occ + nocc, nfrozen_occ + nocc + nvir) + elif label == 'p': + return range(nfrozen_occ, nfrozen_occ + nocc + nvir) + elif label == 'ip': + return range(0, nfrozen_occ) + elif label == 'ap': + return range(nfrozen_occ + nocc + nvir, nfrozen_occ + nocc + nvir + nfrozen_vir) + elif label == 'I': + return range(0, nfrozen_occ + nocc) + elif label == 'A': + return range(nfrozen_occ + nocc, nfrozen_occ + nocc + nvir + nfrozen_vir) + elif label == 'P': + return range(0, nfrozen_occ + nocc + nvir + nfrozen_vir) + else: + raise ValueError('label = {}. is not valid in choose_slice'.format(label)) + + +def contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, label2): + r""" + Contraction in the form of (anti-symmetrized or symmetrized) + I_{tp} = \frac{1}{2} \sum_{qrs} \Gamma_{pq,rs} \langle tq||rs \rangle + = \frac{1}{2} \sum_{qrsP} XY_{pq}^* XY_{rs} + (L^P_{tr} L^P_{qs} \pm L^P_{ts} L^P_{qr}) + = \sum_{qrsP} XY_{pq}^* XY_{rs} L^P_{tr} L^P_{qs} + I_{ti} = \sum_{jklP} Y_{ij}^* Y_{kl} L^P_{tk} L^P_{jl} + + \sum_{jcdP} Y_{ij}^* X_{cd} L^P_{tc} L^P_{jd} + I_{ta} = \sum_{bklP} X_{ab}^* Y_{kl} L^P_{tk} L^P_{bl} + + \sum_{bcdP} X_{ab}^* X_{cd} L^P_{tc} L^P_{bd} + Args: + occ_y_mat: coefficients for occupied orbitals Y + vir_x_mat: coefficients for virtual orbitals X + Lpq_full: full Lpq matrix including frozen orbitals. + nocc: number of active occupied orbitals. + nvir: number of active virtual orbitals. + nfrozen_occ: number of frozen occupied orbitals. + nfrozen_vir: number of frozen virtual orbitals. + label1: label for the first index t + label2: label for the second index p + Returns: + out: contracted intermediates. + """ + # qrs are all active space indices. + slice1 = choose_slice(label1, nfrozen_occ, nocc, nvir, nfrozen_vir) + slice_i = choose_slice('i', nfrozen_occ, nocc, nvir, nfrozen_vir) + slice_a = choose_slice('a', nfrozen_occ, nocc, nvir, nfrozen_vir) + naux = Lpq_full.shape[0] + + # Special cases for TDA + if label1 == 'p': + if nocc == 0: + label1 = 'a' + elif nvir == 0: + label1 = 'i' + if label2 == 'p': + if nocc == 0: + label2 = 'a' + elif nvir == 0: + label2 = 'i' + + if label1 == 'i': + n1 = nocc + elif label1 == 'a': + n1 = nvir + elif label1 == 'ip': + n1 = nfrozen_occ + elif label1 == 'ap': + n1 = nfrozen_vir + else: + n1 = nocc + nvir + if label2 == 'i': + # Slow but more readable version + # out = np.einsum("ij,kl,Ptk,Pjl->ti", occ_y_mat.conj(), occ_y_mat, + # Lpq_full[:,slice1,slice_i], + # Lpq_full[:,slice_i,slice_i], + # optimize=True) + # out+= np.einsum("ij,cd,Ptc,Pjd->ti", occ_y_mat.conj(), vir_x_mat, + # Lpq_full[:,slice1,slice_a], + # Lpq_full[:,slice_i,slice_a], + # optimize=True) + L1i = np.ascontiguousarray(Lpq_full[:, slice1, slice_i]).reshape(-1, nocc) # (Pt,k) + Lij = np.ascontiguousarray(Lpq_full[:, slice_i, slice_i]).reshape(-1, nocc).conj() # (P,j,l)*->(Pl,j) + L1i = np.matmul(L1i, occ_y_mat).reshape(naux, n1, nocc) # (Pt,k)(k,l)->(P,t,l) + L1i = L1i.transpose(1, 0, 2).reshape(n1, -1) # (P,t,l)->(t,Pl) + tmp = np.matmul(L1i, Lij) # (t,Pl)(Pl,j)->(t,j) + out = np.matmul(tmp, occ_y_mat.T.conj()) # (t,j)(j,i) -> (t,i) + + if nvir > 0: + L1a = np.ascontiguousarray(Lpq_full[:, slice1, slice_a]).reshape(-1, nvir) # (Pt,c) + Lai = np.ascontiguousarray(Lpq_full[:, slice_a, slice_i]).reshape(-1, nocc).conj() # (P,d,j)*->(Pd,j) + L1a = np.matmul(L1a, vir_x_mat).reshape(naux, n1, nvir) # (Pt,c)(c,d)->(P,t,d) + L1a = L1a.transpose(1, 0, 2).reshape(n1, -1) # (P,t,d)->(t,Pd) + tmp = np.matmul(L1a, Lai) # (t,Pd)(Pd,j)->(t,j) + out += np.matmul(tmp, occ_y_mat.T.conj()) # (t,j)(j,i) -> (t,i) + elif label2 == 'a': + # Slow but more readable version + # out = np.einsum("ab,cd,Ptc,Pbd->ta", vir_x_mat.conj(), vir_x_mat, + # Lpq_full[:,slice1,slice_a], + # Lpq_full[:,slice_a,slice_a], + # optimize=True) + # out+= np.einsum("ab,kl,Ptk,Pbl->ta", vir_x_mat.conj(), occ_y_mat, + # Lpq_full[:,slice1,slice_i], + # Lpq_full[:,slice_a,slice_i], + # optimize=True) + L1a = np.ascontiguousarray(Lpq_full[:, slice1, slice_a]).reshape(-1, nvir) # (Pt,c) + Lab = np.ascontiguousarray(Lpq_full[:, slice_a, slice_a]).reshape(-1, nvir).conj() # (P,b,d)*->(Pd,b) + L1a = np.matmul(L1a, vir_x_mat).reshape(naux, n1, nvir) # (Pt,c)(c,d)->(P,t,d) + L1a = L1a.transpose(1, 0, 2).reshape(n1, -1) # (P,t,d)->(t,Pd) + tmp = np.matmul(L1a, Lab) # (t,Pd)(Pd,b)->(t,b) + out = np.matmul(tmp, vir_x_mat.T.conj()) # (t,b)(b,a) -> (t,a) + + if nocc > 0: + L1i = np.ascontiguousarray(Lpq_full[:, slice1, slice_i]).reshape(-1, nocc) # (Pt,k) + Lia = np.ascontiguousarray(Lpq_full[:, slice_i, slice_a]).reshape(-1, nvir).conj() # (P,l,b)*->(Pl,b) + L1i = np.matmul(L1i, occ_y_mat).reshape(naux, n1, nocc) # (Pt,k)(k,l)->(P,t,l) + L1i = L1i.transpose(1, 0, 2).reshape(n1, -1) # (P,t,l)->(t,Pl) + tmp = np.matmul(L1i, Lia) # (t,Pl)(Pl,b)->(t,b) + out += np.matmul(tmp, vir_x_mat.T.conj()) # (t,b)(b,a) -> (t,a) + elif label2 == 'p': + # slow (more copies) but more readable version + # out = np.concatenate(( + # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "i"), + # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "a")), axis=1) + out = np.zeros((n1, nocc + nvir), dtype=occ_y_mat.dtype) + L1i = np.ascontiguousarray(Lpq_full[:, slice1, slice_i]).reshape(-1, nocc) # (Pt,k) + Lia = np.ascontiguousarray(Lpq_full[:, slice_i, slice_a]).reshape(-1, nvir).conj() # (P,l,b)*->(Pl,b) + Lij = np.ascontiguousarray(Lpq_full[:, slice_i, slice_i]).reshape(-1, nocc).conj() # (P,j,l)*->(Pl,j) + L1i = np.matmul(L1i, occ_y_mat).reshape(naux, n1, nocc) # (Pt,k)(k,l)->(P,t,l) + L1i = L1i.transpose(1, 0, 2).reshape(n1, -1) # (P,t,l)->(t,Pl) + tmp = np.matmul(L1i, Lia) # (t,Pl)(Pl,b)->(t,b) + out[:, nocc:] = np.matmul(tmp, vir_x_mat.T.conj()) # (t,b)(b,a) -> (t,a) + tmp = np.matmul(L1i, Lij) # (t,Pl)(Pl,j)->(t,j) + out[:, :nocc] = np.matmul(tmp, occ_y_mat.T.conj()) # (t,j)(j,i) -> (t,i) + + L1a = np.ascontiguousarray(Lpq_full[:, slice1, slice_a]).reshape(-1, nvir) # (Pt,c) + Lab = np.ascontiguousarray(Lpq_full[:, slice_a, slice_a]).reshape(-1, nvir).conj() # (P,b,d)*->(Pd,b) + Lai = np.ascontiguousarray(Lpq_full[:, slice_a, slice_i]).reshape(-1, nocc).conj() # (P,d,j)*->(Pd,j) + L1a = np.matmul(L1a, vir_x_mat).reshape(naux, n1, nvir) # (Pt,c)(c,d)->(P,t,d) + L1a = L1a.transpose(1, 0, 2).reshape(n1, -1) # (P,t,d)->(t,Pd) + tmp = np.matmul(L1a, Lab) # (t,Pd)(Pd,b)->(t,b) + out[:, nocc:] += np.matmul(tmp, vir_x_mat.T.conj()) # (t,b)(b,a) -> (t,a) + tmp = np.matmul(L1a, Lai) # (t,Pd)(Pd,j)->(t,j) + out[:, :nocc] += np.matmul(tmp, occ_y_mat.T.conj()) # (t,j)(j,i) -> (t,i) + else: + raise ValueError('label2 = {}. is not valid in contraction_2rdm_Lpq'.format(label2)) + + return out + + +''' +The following functions are modified from pyscf functions +to accommodate auxbasis_response in pprpa. +''' +def get_veff_df_rks(ks_grad, mol=None, dm=None): + """Coulomb + XC functional response + Modified from pyscf.df.grad.rks.get_veff + """ + if mol is None: + mol = ks_grad.mol + if dm is None: + dm = ks_grad.base.make_rdm1() + t0 = (logger.process_clock(), logger.perf_counter()) + + mf = ks_grad.base + ni = mf._numint + grids, nlcgrids = rks_grad._initialize_grids(ks_grad) + + mem_now = lib.current_memory()[0] + max_memory = max(2000, ks_grad.max_memory * 0.9 - mem_now) + if ks_grad.grid_response: + exc = [] + vxc = [] + for dmi in dm: + exci, vxci = rks_grad.get_vxc_full_response( + ni, mol, grids, mf.xc, dmi, max_memory=max_memory, verbose=ks_grad.verbose + ) + exc.append(exci) + vxc.append(vxci) + exc = np.asarray(exc) + vxc = np.asarray(vxc) + if mf.do_nlc(): + if ni.libxc.is_nlc(mf.xc): + xc = mf.xc + else: + xc = mf.nlc + enlc, vnlc = rks_grad.get_nlc_vxc_full_response( + ni, mol, nlcgrids, xc, dm, max_memory=max_memory, verbose=ks_grad.verbose + ) + exc += enlc + vxc += vnlc + logger.debug1(ks_grad, 'sum(grids response) %s', exc.sum(axis=0)) + else: + exc, vxc = rks_grad.get_vxc(ni, mol, grids, mf.xc, dm, max_memory=max_memory, verbose=ks_grad.verbose) + if mf.do_nlc(): + if ni.libxc.is_nlc(mf.xc): + xc = mf.xc + else: + xc = mf.nlc + enlc, vnlc = rks_grad.get_nlc_vxc(ni, mol, nlcgrids, xc, dm, max_memory=max_memory, verbose=ks_grad.verbose) + vxc += vnlc + t0 = logger.timer(ks_grad, 'vxc', *t0) + + vjk = np.zeros_like(vxc) + if not ni.libxc.is_hybrid_xc(mf.xc): + vj = ks_grad.get_j(mol, dm) + vjk += vj + if ks_grad.auxbasis_response: + e1_aux = vj.aux + else: + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, spin=mol.spin) + vj, vk = ks_grad.get_jk(mol, dm) + if ks_grad.auxbasis_response: + vk.aux *= hyb + vk[:] *= hyb # Don't erase the .aux tags! + if omega != 0: # For range separated Coulomb operator + # TODO: replaced with vk_sr which is numerically more stable for + # inv(int2c2e) + vk_lr = ks_grad.get_k(mol, dm, omega=omega) + vk[:] += vk_lr * (alpha - hyb) + if ks_grad.auxbasis_response: + vk.aux[:] += vk_lr.aux * (alpha - hyb) + vjk += vj - vk * 0.5 + if ks_grad.auxbasis_response: + e1_aux = vj.aux - vk.aux * 0.5 + + if ks_grad.auxbasis_response: + vjk = lib.tag_array(vjk, aux=e1_aux) + if ks_grad.grid_response: + vxc = lib.tag_array(vxc, exc1_grid=exc) + return vxc, vjk + + +def get_veff_rks(ks_grad, mol=None, dm=None): + """ + First order derivative of DFT effective potential matrix (wrt electron coordinates) + + Args: + ks_grad : grad.uhf.Gradients or grad.uks.Gradients object + """ + if mol is None: + mol = ks_grad.mol + if dm is None: + dm = ks_grad.base.make_rdm1() + t0 = (logger.process_clock(), logger.perf_counter()) + + mf = ks_grad.base + ni = mf._numint + grids, nlcgrids = rks_grad._initialize_grids(ks_grad) + + mem_now = lib.current_memory()[0] + max_memory = max(2000, ks_grad.max_memory * 0.9 - mem_now) + if ks_grad.grid_response: + exc, vxc = rks_grad.get_vxc_full_response( + ni, mol, grids, mf.xc, dm, max_memory=max_memory, verbose=ks_grad.verbose + ) + if mf.do_nlc(): + if ni.libxc.is_nlc(mf.xc): + xc = mf.xc + else: + xc = mf.nlc + enlc, vnlc = rks_grad.get_nlc_vxc_full_response( + ni, mol, nlcgrids, xc, dm, max_memory=max_memory, verbose=ks_grad.verbose + ) + exc += enlc + vxc += vnlc + logger.debug1(ks_grad, 'sum(grids response) %s', exc.sum(axis=0)) + else: + exc, vxc = rks_grad.get_vxc(ni, mol, grids, mf.xc, dm, max_memory=max_memory, verbose=ks_grad.verbose) + if mf.do_nlc(): + if ni.libxc.is_nlc(mf.xc): + xc = mf.xc + else: + xc = mf.nlc + enlc, vnlc = rks_grad.get_nlc_vxc(ni, mol, nlcgrids, xc, dm, max_memory=max_memory, verbose=ks_grad.verbose) + vxc += vnlc + t0 = logger.timer(ks_grad, 'vxc', *t0) + + vjk = np.zeros_like(vxc) + if not ni.libxc.is_hybrid_xc(mf.xc): + vj = ks_grad.get_j(mol, dm) + vjk += vj + else: + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, spin=mol.spin) + vj, vk = ks_grad.get_jk(mol, dm) + vk *= hyb + if omega != 0: + vk += ks_grad.get_k(mol, dm, omega=omega) * (alpha - hyb) + vjk += vj - vk * 0.5 + + vxc = lib.tag_array(vxc, exc1_grid=exc) + return vxc, vjk + + +def _contract_xc_kernel(mf, xc_code, dmvo, dmoo=None, with_vxc=True, with_kxc=True, singlet=True, max_memory=2000): + from pyscf import lib + from pyscf.lib import logger + from pyscf.grad.tdrks import _lda_eval_mat_, _gga_eval_mat_, _mgga_eval_mat_ + + mol = mf.mol + grids = mf.grids + + ni = mf._numint + xctype = ni._xc_type(xc_code) + + mo_coeff = mf.mo_coeff + mo_occ = mf.mo_occ + nao, nmo = mo_coeff.shape + shls_slice = (0, mol.nbas) + ao_loc = mol.ao_loc_nr() + + # dmvo ~ reduce(np.dot, (orbv, Xai, orbo.T)) + dmvo = (dmvo + dmvo.T) * 0.5 # because K_{ia,jb} == K_{ia,bj} + + f1vo = np.zeros((4, nao, nao)) # 0th-order, d/dx, d/dy, d/dz + deriv = 2 + if dmoo is not None: + f1oo = np.zeros((4, nao, nao)) + else: + f1oo = None + if with_vxc: + v1ao = np.zeros((4, nao, nao)) + else: + v1ao = None + if with_kxc: + k1ao = np.zeros((4, nao, nao)) + deriv = 3 + else: + k1ao = None + + if xctype == 'HF': + return f1vo, f1oo, v1ao, k1ao + elif xctype == 'LDA': + fmat_, ao_deriv = _lda_eval_mat_, 1 + elif xctype == 'GGA': + fmat_, ao_deriv = _gga_eval_mat_, 2 + elif xctype == 'MGGA': + fmat_, ao_deriv = _mgga_eval_mat_, 2 + logger.warn(mf, 'PPRPA-MGGA Gradients may be inaccurate due to grids response') + else: + raise NotImplementedError(f'td-rks for functional {xc_code}') + + if singlet: + for ao, mask, weight, coords in ni.block_loop(mol, grids, nao, ao_deriv, max_memory): + if xctype == 'LDA': + ao0 = ao[0] + else: + ao0 = ao + rho = ni.eval_rho2(mol, ao0, mo_coeff, mo_occ, mask, xctype, with_lapl=False) + vxc, fxc, kxc = ni.eval_xc_eff(xc_code, rho, deriv, xctype=xctype)[1:] + + rho1 = ni.eval_rho(mol, ao0, dmvo, mask, xctype, hermi=1, with_lapl=False) * 2 # *2 for alpha + beta + if xctype == 'LDA': + rho1 = rho1[np.newaxis] + wv = np.einsum('yg,xyg,g->xg', rho1, fxc, weight) + fmat_(mol, f1vo, ao, wv, mask, shls_slice, ao_loc) + + if dmoo is not None: + rho2 = ni.eval_rho(mol, ao0, dmoo, mask, xctype, hermi=1, with_lapl=False) * 2 + if xctype == 'LDA': + rho2 = rho2[np.newaxis] + wv = np.einsum('yg,xyg,g->xg', rho2, fxc, weight) + fmat_(mol, f1oo, ao, wv, mask, shls_slice, ao_loc) + if with_vxc: + fmat_(mol, v1ao, ao, vxc * weight, mask, shls_slice, ao_loc) + if with_kxc: + wv = np.einsum('yg,zg,xyzg,g->xg', rho1, rho1, kxc, weight) + fmat_(mol, k1ao, ao, wv, mask, shls_slice, ao_loc) + else: + for ao, mask, weight, coords in ni.block_loop(mol, grids, nao, ao_deriv, max_memory): + if xctype == 'LDA': + ao0 = ao[0] + else: + ao0 = ao + rho = ni.eval_rho2(mol, ao0, mo_coeff, mo_occ, mask, xctype, with_lapl=False) + rho *= 0.5 + rho = np.repeat(rho[np.newaxis], 2, axis=0) + vxc, fxc, kxc = ni.eval_xc_eff(xc_code, rho, deriv, xctype=xctype)[1:] + # fxc_t couples triplet excitation amplitudes + # 1/2 int (tia - tIA) fxc (tjb - tJB) = tia fxc_t tjb + fxc_t = fxc[:, :, 0] - fxc[:, :, 1] + fxc_t = fxc_t[0] - fxc_t[1] + + rho1 = ni.eval_rho(mol, ao0, dmvo, mask, xctype, hermi=1, with_lapl=False) + if xctype == 'LDA': + rho1 = rho1[np.newaxis] + wv = np.einsum('yg,xyg,g->xg', rho1, fxc_t, weight) + fmat_(mol, f1vo, ao, wv, mask, shls_slice, ao_loc) + + if dmoo is not None: + # fxc_s == 2 * fxc of spin restricted xc kernel + # provides f1oo to couple the interaction between first order MO + # and density response of tddft amplitudes, which is described by dmoo + fxc_s = fxc[0, :, 0] + fxc[0, :, 1] + rho2 = ni.eval_rho(mol, ao0, dmoo, mask, xctype, hermi=1, with_lapl=False) + if xctype == 'LDA': + rho2 = rho2[np.newaxis] + wv = np.einsum('yg,xyg,g->xg', rho2, fxc_s, weight) + fmat_(mol, f1oo, ao, wv, mask, shls_slice, ao_loc) + if with_vxc: + vxc = vxc[0] + fmat_(mol, v1ao, ao, vxc * weight, mask, shls_slice, ao_loc) + if with_kxc: + # kxc in terms of the triplet coupling + # 1/2 int (tia - tIA) kxc (tjb - tJB) = tia kxc_t tjb + kxc = kxc[0, :, 0] - kxc[0, :, 1] + kxc = kxc[:, :, 0] - kxc[:, :, 1] + wv = np.einsum('yg,zg,xyzg,g->xg', rho1, rho1, kxc, weight) + fmat_(mol, k1ao, ao, wv, mask, shls_slice, ao_loc) + + f1vo[1:] *= -1 + if f1oo is not None: + f1oo[1:] *= -1 + if v1ao is not None: + v1ao[1:] *= -1 + if k1ao is not None: + k1ao[1:] *= -1 + + return f1vo, f1oo, v1ao, k1ao + + +class Gradients(rhf_grad.Gradients): + cphf_max_cycle = 20 + cphf_conv_tol = 1e-8 + + _keys = { + 'cphf_max_cycle', + 'cphf_conv_tol', + 'mol', + 'base', + 'chkfile', + 'state', + 'atmlst', + 'de', + } + + def __init__(self, pprpa, mult='t', state=0): + assert isinstance(pprpa._scf, scf.hf.RHF) + self.base = pprpa + self.mol = pprpa._scf.mol + self.state = state + self.verbose = self.mol.verbose + self.mult = mult + + self.rdm1e = None + self.atmlst = None + self.de = None + + def dump_flags(self, verbose=None): + log = logger.new_logger(self, verbose) + log.info('\n') + log.info('******** %s gradients for %s ********', self.base.__class__, self.base._scf.__class__) + log.info('cphf_conv_tol = %g', self.cphf_conv_tol) + log.info('cphf_max_cycle = %d', self.cphf_max_cycle) + log.info('State ID = %d', self.state) + log.info('max_memory %d MB (current use %d MB)', self.base._scf.max_memory, lib.current_memory()[0]) + log.info('\n') + return self + + def grad_elec(self, xy, mult, atmlst): + return grad_elec(self, xy, mult, atmlst) + + def _finalize(self): + if self.verbose >= logger.NOTE: + logger.note( + self, '--------- %s gradients for state %d ----------', self.base.__class__.__name__, self.state + ) + self._write(self.mol, self.de, self.atmlst) + logger.note(self, '----------------------------------------------') + + def kernel(self, xy=None, state=None, mult=None, atmlst=None): + if mult is None: + mult = self.mult + if xy is None: + if state is None: + state = self.state + else: + self.state = state + if mult == 't': + xy = self.base.xy_t[state] + else: + xy = self.base.xy_s[state] + if atmlst is None: + atmlst = self.atmlst + else: + self.atmlst = atmlst + + if self.verbose >= logger.INFO: + self.dump_flags() + + de = self.grad_elec(xy, mult, atmlst) + self.de = de = de + self.grad_nuc(atmlst=atmlst) + if self.mol.symmetry: + self.de = self.symmetrize(self.de, atmlst) + self._finalize() + return self.de + + +Grad = Gradients + +from pyscf.pprpa.rpprpa_direct import RppRPADirect +from pyscf.pprpa.rpprpa_davidson import RppRPADavidson + +RppRPADirect.Gradients = RppRPADavidson.Gradients = lib.class_as_method(Gradients) \ No newline at end of file diff --git a/pyscf/pprpa/gpprpa_davidson.py b/pyscf/pprpa/gpprpa_davidson.py new file mode 100644 index 000000000..1032a03c7 --- /dev/null +++ b/pyscf/pprpa/gpprpa_davidson.py @@ -0,0 +1,418 @@ +import numpy, scipy +from pyscf import df, dft, pbc, scf +from pyscf.ao2mo._ao2mo import nr_e2 +from pyscf.lib import logger, current_memory +from .rpprpa_davidson import RppRPADavidson, pprpa_get_trial_vector, pprpa_expand_space +from .rpprpa_direct import pprpa_print_a_pair, inner_product, pprpa_orthonormalize_eigenvector + + + +def kernel(pprpa): + # initialize trial vector and product matrix + data_type = pprpa._scf.mo_coeff.dtype + tri_vec = numpy.zeros(shape=[pprpa.max_vec, pprpa.full_dim], dtype=data_type) + tri_vec_sig = numpy.zeros(shape=[pprpa.max_vec], dtype=data_type) + if pprpa.channel == "pp": + ntri = min(pprpa.nroot * 4, pprpa.vv_dim) + else: + ntri = min(pprpa.nroot * 4, pprpa.oo_dim) + tri_vec[:ntri], tri_vec_sig[:ntri] = pprpa_get_trial_vector(pprpa, ntri) + mv_prod = numpy.zeros_like(tri_vec) + + iter = 0 + nprod = 0 # number of contracted vectors + while iter < pprpa.max_iter: + logger.info( + pprpa, "\nppRPA Davidson %d-th iteration, ntri= %d , nprod= %d .", + iter + 1, ntri, nprod) + mv_prod[nprod:ntri] = pprpa_contraction(pprpa, tri_vec[nprod:ntri]) + nprod = ntri + + # get ppRPA matrix and metric matrix in subspace + m_tilde = numpy.dot(tri_vec[:ntri].conj(), mv_prod[:ntri].T) + w_tilde = numpy.zeros_like(m_tilde) + for i in range(ntri): + if inner_product(tri_vec[i].conj(), tri_vec[i], pprpa.oo_dim).real > 0: + w_tilde[i, i] = 1 + else: + w_tilde[i, i] = -1 + + # diagonalize subspace matrix + if data_type == numpy.double: + alphar, _, beta, _, v_tri, _, _ = scipy.linalg.lapack.dggev( + m_tilde, w_tilde, compute_vl=0) + elif data_type == numpy.complex128: + alphar, beta, _, v_tri, _, _ = scipy.linalg.lapack.zggev( + m_tilde, w_tilde, compute_vl=0) + e_tri = (alphar / beta).real + v_tri = v_tri.T # Fortran matrix to Python order + + if pprpa.channel == "pp": + # sort eigenvalues and eigenvectors by ascending order + idx = e_tri.argsort() + e_tri = e_tri[idx] + v_tri = v_tri[idx, :] + + # re-order all states by signs, first hh then pp + sig = numpy.zeros(shape=[ntri], dtype=int) + for i in range(ntri): + if numpy.sum(v_tri[i].conj() * tri_vec_sig[:ntri] * v_tri[i]).real > 0: + sig[i] = 1 + else: + sig[i] = -1 + + hh_index = numpy.where(sig < 0)[0] + pp_index = numpy.where(sig > 0)[0] + e_tri_hh = e_tri[hh_index] + e_tri_pp = e_tri[pp_index] + e_tri[:len(hh_index)] = e_tri_hh + e_tri[len(hh_index):] = e_tri_pp + v_tri_hh = v_tri[hh_index] + v_tri_pp = v_tri[pp_index] + v_tri[:len(hh_index)] = v_tri_hh + v_tri[len(hh_index):] = v_tri_pp + + # get only two-electron addition energy + first_state=len(hh_index) + pprpa.exci = e_tri[first_state:first_state+pprpa.nroot] + else: + # sort eigenvalues and eigenvectors by descending order + idx = e_tri.argsort()[::-1] + e_tri = e_tri[idx] + v_tri = v_tri[idx, :] + + # re-order all states by signs, first pp then hh + sig = numpy.zeros(shape=[ntri], dtype=int) + for i in range(ntri): + if numpy.sum(v_tri[i].conj() * tri_vec_sig[:ntri] * v_tri[i]).real > 0: + sig[i] = 1 + else: + sig[i] = -1 + + hh_index = numpy.where(sig < 0)[0] + pp_index = numpy.where(sig > 0)[0] + e_tri_hh = e_tri[hh_index] + e_tri_pp = e_tri[pp_index] + e_tri[:len(pp_index)] = e_tri_pp + e_tri[len(pp_index):] = e_tri_hh + v_tri_hh = v_tri[hh_index] + v_tri_pp = v_tri[pp_index] + v_tri[:len(pp_index)] = v_tri_pp + v_tri[len(pp_index):] = v_tri_hh + + # get only two-electron removal energy + first_state=len(pp_index) + pprpa.exci = e_tri[first_state:first_state+pprpa.nroot] + + ntri_old = ntri + conv, ntri = pprpa_expand_space( + pprpa=pprpa, first_state=first_state, tri_vec=tri_vec, + tri_vec_sig=tri_vec_sig, mv_prod=mv_prod, v_tri=v_tri) + logger.info(pprpa, "add %d new trial vectors.", ntri - ntri_old) + + iter += 1 + if conv is True: + break + + assert conv is True, "ppRPA Davidson is not converged!" + logger.info( + pprpa, "\nppRPA Davidson converged in %d iterations, final subspace size = %d", + iter, nprod) + + pprpa_orthonormalize_eigenvector(pprpa.multi, pprpa.nocc_act, pprpa.exci, pprpa.xy) + + return + + +def pprpa_contraction(pprpa, tri_vec): + """GppRPA contraction. + + Args: + pprpa (GppRPA_Davidson): GppRPA_Davidson object. + tri_vec (double/complex ndarray): trial vector. + + Returns: + mv_prod (double/complex ndarray): product between ppRPA matrix and trial vectors. + """ + # for convenience, use XX as XX_act in this function + nocc, nvir, nmo = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act + mo_energy = pprpa.mo_energy_act + Lpi = pprpa.Lpi + Lpa = pprpa.Lpa + naux = Lpi.shape[0] + + ntri = tri_vec.shape[0] + mv_prod = numpy.zeros(shape=[ntri, pprpa.full_dim], dtype=tri_vec.dtype) + + tri_row_o, tri_col_o = numpy.tril_indices(nocc, -1) + tri_row_v, tri_col_v = numpy.tril_indices(nvir, -1) + + for ivec in range(ntri): + z_oo = numpy.zeros(shape=[nocc, nocc], dtype=tri_vec.dtype) + z_vv = numpy.zeros(shape=[nvir, nvir], dtype=tri_vec.dtype) + z_oo[tri_row_o, tri_col_o] = tri_vec[ivec][:pprpa.oo_dim] + z_oo[numpy.diag_indices(nocc)] *= 1.0 / numpy.sqrt(2) + z_oo = numpy.ascontiguousarray(z_oo.T) + z_vv[tri_row_v, tri_col_v] = tri_vec[ivec][pprpa.oo_dim:] + z_vv[numpy.diag_indices(nvir)] *= 1.0 / numpy.sqrt(2) + z_vv = numpy.ascontiguousarray(z_vv.T) + + # Lpqz_{L,pr} = \sum_s Lpq_{L,ps} z_{rs} + Lpq_z = numpy.zeros(shape=[naux * nmo, nmo], dtype=tri_vec.dtype) + Lpq_z[:, :nocc] = Lpi.reshape(naux * nmo, nocc) @ z_oo + Lpq_z[:, nocc:] = Lpa.reshape(naux * nmo, nvir) @ z_vv + + # transpose and reshape for faster multiplication + Lpq_z = Lpq_z.reshape(naux, nmo, nmo).transpose(1, 0, 2).conj() + Lpq_z = Lpq_z.reshape(nmo, naux * nmo) + + # MV_{pq} = \sum_{Lr} Lpq_{L,pr} Lpqz_{L,qr} - Lpq_{L,qr} Lpqz_{L,pr} + # -MV_{qp}* = - \sum_{Lr} Lpq_{L,rp} Lpqz_{L,qr}^* + Lpq_{L,rq} Lpqz_{L,pr}^* + mv_prod_full = numpy.zeros(shape=[nmo, nmo], dtype=tri_vec.dtype) + mv_prod_full[:nocc, :nocc] = Lpq_z[:nocc] @ Lpi.reshape(naux * nmo, nocc) + mv_prod_full[nocc:, nocc:] = Lpq_z[nocc:] @ Lpa.reshape(naux * nmo, nvir) + mv_prod_full -= mv_prod_full.T + mv_prod_full[numpy.diag_indices(nmo)] *= 1.0 / numpy.sqrt(2) + mv_prod[ivec][: pprpa.oo_dim] =\ + -mv_prod_full[:nocc, :nocc][tri_row_o, tri_col_o].conj() + mv_prod[ivec][pprpa.oo_dim:] = \ + -mv_prod_full[nocc:, nocc:][tri_row_v, tri_col_v].conj() + + orb_sum_oo = (mo_energy[None, : nocc] + mo_energy[: nocc, None]) - 2.0 * pprpa.mu + orb_sum_oo = orb_sum_oo[tri_row_o, tri_col_o] + orb_sum_vv = (mo_energy[None, nocc:] + mo_energy[nocc:, None]) - 2.0 * pprpa.mu + orb_sum_vv = orb_sum_vv[tri_row_v, tri_col_v] + for ivec in range(ntri): + oz_oo = -orb_sum_oo * tri_vec[ivec][:pprpa.oo_dim] + mv_prod[ivec][:pprpa.oo_dim] += oz_oo + oz_vv = orb_sum_vv * tri_vec[ivec][pprpa.oo_dim:] + mv_prod[ivec][pprpa.oo_dim:] += oz_vv + + return mv_prod + + +def ao2mo(pprpa, full_mo=False): + """Get three-center density-fitting matrix in MO active space. + + Args: + pprpa: GppRPA object. + + Returns: + Lpq (complex/double ndarray): three-center DF matrix in MO active space. + """ + mf = pprpa._scf + mo_coeff = mf.mo_coeff + nocc = pprpa.nocc + + nao = pprpa._scf.mol.nao_nr() + mo = numpy.asarray(mo_coeff, order='F') + if full_mo: + nocc_act, nvir_act, nmo_act = nocc, mo.shape[1] - nocc, mo.shape[1] + ijslice = (0, mo.shape[1], 0, mo.shape[1]) + else: + nocc_act, nvir_act, nmo_act = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act + ijslice = (nocc-nocc_act, nocc+nvir_act, nocc-nocc_act, nocc+nvir_act) + + if isinstance(mf, (scf.ghf.GHF, dft.gks.GKS)): + # molecule + if getattr(mf, 'with_df', None): + pprpa.with_df = mf.with_df + else: + pprpa.with_df = df.DF(mf.mol) + if pprpa.auxbasis is not None: + pprpa.with_df.auxbasis = pprpa.auxbasis + else: + pprpa.with_df.auxbasis = df.make_auxbasis( + mf.mol, mp2fit=True) + pprpa._keys.update(['with_df']) + + naux = pprpa.with_df.get_naoaux() + mem_incore = (2*nao**2*naux) * 8/1e6 + mem_incore += (nmo_act**2*naux) * 16/1e6 + mem_now = current_memory()[0] + + mo_a = mf.mo_coeff[:nao, :] + mo_b = mf.mo_coeff[nao:, :] + + Lpq = None + if (mem_incore + mem_now < pprpa.max_memory) or pprpa.mol.incore_anyway: + if mo_a.dtype == numpy.double: # real coefficients, no SOC + Lpq = nr_e2(pprpa.with_df._cderi, mo_a, ijslice, aosym='s2') + Lpq += nr_e2(pprpa.with_df._cderi, mo_b, ijslice, aosym='s2') + else: + Lpq = numpy.zeros((naux, nmo_act * nmo_act), dtype=numpy.complex128) + ijslice = (ijslice[0], ijslice[1], mo_a.shape[1] + ijslice[2], mo_a.shape[1] + ijslice[3]) + mo_tmp = numpy.asarray(numpy.hstack((mo_a.real, mo_a.real)), order='F') + Lpq += nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + mo_tmp = numpy.asarray(numpy.hstack((mo_a.imag, mo_a.imag)), order='F') + Lpq += nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + mo_tmp = numpy.asarray(numpy.hstack((mo_b.real, mo_b.real)), order='F') + Lpq += nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + mo_tmp = numpy.asarray(numpy.hstack((mo_b.imag, mo_b.imag)), order='F') + Lpq += nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + mo_tmp = numpy.asarray(numpy.hstack((mo_a.real, mo_a.imag)), order='F') + Lpq += 1j * nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + mo_tmp = numpy.asarray(numpy.hstack((mo_a.imag, -mo_a.real)), order='F') + Lpq += 1j * nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + mo_tmp = numpy.asarray(numpy.hstack((mo_b.real, mo_b.imag)), order='F') + Lpq += 1j * nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + mo_tmp = numpy.asarray(numpy.hstack((mo_b.imag, -mo_b.real)), order='F') + Lpq += 1j * nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') + + return Lpq.reshape(naux, nmo_act, nmo_act) + else: + logger.warn(pprpa, 'Memory may not be enough!') + raise NotImplementedError + elif isinstance(mf, (pbc.scf.ghf.GHF, pbc.dft.gks.GKS)): + raise NotImplementedError("GppRPA for periodic systems is not implemented yet.") + else: + raise NotImplementedError("GppRPA only supports GHF/GKS reference now.") + + +def complex_matrix_norm(matrix): + """Calculate the Complex norm of each matrix element. + + Args: + matrix (complex ndarray): input matrix. + + Returns: + out (double ndarray): norm square of each element. + """ + out = numpy.zeros(matrix.shape, dtype=numpy.double) + out = numpy.power(numpy.abs(matrix.real), 2) + numpy.power(numpy.abs(matrix.imag), 2) + return out + + +# analysis functions +def pprpa_davidson_print_eigenvector(pprpa, exci0, exci, xy): + """Print dominant components of an eigenvector. + + Args: + pprpa (RppRPADavidson): ppRPA object. + exci0 (double): lowest eigenvalue. + exci (double array): ppRPA eigenvalue. + xy (double ndarray): ppRPA eigenvector. + """ + nocc = pprpa.nocc + nocc_act, nvir_act = pprpa.nocc_act, pprpa.nvir_act + oo_dim = int((nocc_act - 1) * nocc_act / 2) + logger.info(pprpa, "\n print GppRPA excitations: triplet\n") + + tri_row_o, tri_col_o = numpy.tril_indices(nocc_act, -1) + tri_row_v, tri_col_v = numpy.tril_indices(nvir_act, -1) + + nroot = len(exci) + au2ev = 27.211386 + if pprpa.channel == "pp": + for iroot in range(nroot): + logger.info(pprpa, "#%-d excitation: exci= %-12.4f eV 2e= %-12.4f eV", + iroot + 1, (exci[iroot] - exci0) * au2ev, exci[iroot] * au2ev) + if nocc_act > 0: + full = numpy.zeros(shape=[nocc_act, nocc_act], dtype=xy.dtype) + full[tri_row_o, tri_col_o] = xy[iroot][:oo_dim] + full = complex_matrix_norm(full) + pairs = numpy.argwhere(full > pprpa.print_thresh) + for i, j in pairs: + pprpa_print_a_pair( + pprpa, is_pp=False, p=i, q=j, percentage=full[i, j]) + + full = numpy.zeros(shape=[nvir_act, nvir_act], dtype=xy.dtype) + full[tri_row_v, tri_col_v] = xy[iroot][oo_dim:] + full = complex_matrix_norm(full) + pairs = numpy.argwhere(full > pprpa.print_thresh) + for a, b in pairs: + pprpa_print_a_pair( + pprpa, is_pp=True, p=a+nocc, q=b+nocc, percentage=full[a, b]) + logger.info(pprpa, "") + else: + for iroot in range(nroot): + logger.info(pprpa, "#%-d de-excitation: exci= %-12.4f eV 2e= %-12.4f eV", + iroot + 1, (exci[iroot] - exci0) * au2ev, exci[iroot] * au2ev) + full = numpy.zeros(shape=[nocc_act, nocc_act], dtype=xy.dtype) + full[tri_row_o, tri_col_o] = xy[iroot][:oo_dim] + full = complex_matrix_norm(full) + pairs = numpy.argwhere(full > pprpa.print_thresh) + for i, j in pairs: + pprpa_print_a_pair( + pprpa, is_pp=False, p=i, q=j, percentage=full[i, j]) + + if nvir_act > 0: + full = numpy.zeros(shape=[nvir_act, nvir_act], dtype=xy.dtype) + full[tri_row_v, tri_col_v] = xy[iroot][oo_dim:] + full = complex_matrix_norm(full) + pairs = numpy.argwhere(full > pprpa.print_thresh) + for a, b in pairs: + pprpa_print_a_pair( + pprpa, is_pp=True, p=a+nocc, q=b+nocc, percentage=full[a, b]) + logger.info(pprpa, "") + + return + + +def analyze_pprpa_davidson(pprpa): + """Analyze ppRPA (Davidson algorithm) excited states. + + Args: + pprpa (RppRPADavidson): GppRPA object. + """ + logger.info(pprpa, "\nanalyze GppRPA results.") + + assert pprpa.exci is not None + if pprpa.channel == "pp": + exci0 = pprpa.exci[0] + else: + exci0 = pprpa.exci[0] + pprpa_davidson_print_eigenvector( + pprpa, exci0=exci0, exci=pprpa.exci, xy=pprpa.xy) + return + + +class GppRPADavidson(RppRPADavidson): + def __init__( + self, mf, nocc_act=None, nvir_act=None, channel="pp", nroot=5, + max_vec=200, max_iter=100, residue_thresh=1.0e-7, print_thresh=0.1, + auxbasis=None): + super().__init__( + mf, nocc_act, nvir_act, channel, nroot, max_vec, max_iter, + residue_thresh, print_thresh, auxbasis) + # Most of the GPPRPA formulations are just the complex-valued + # extention of those in triplet RPPRA, since they both use + # the antisymmetric two-electron integrals. + self.multi = "t" + + analyze = analyze_pprpa_davidson + ao2mo = ao2mo + + def check_memory(self): + """Check required memory. + """ + # intermediate in contraction; mv_prod, tri_vec, xy + mem = (3 * self.max_vec * self.full_dim) * 8 / 1.0e6 + if self._scf.mo_coeff.dtype == numpy.complex128: + mem *= 2 + if mem < 1000: + logger.info(self, "GppRPA needs at least %d MB memory.", mem) + else: + logger.info(self, "GppRPA needs at least %.1f GB memory.", mem/1.0e3) + return + + def kernel(self): + """Run GppRPA davidson diagonalization. + """ + self.check_parameter() + self.check_memory() + + cput0 = (logger.process_clock(), logger.perf_counter()) + if self.Lpi is None or self.Lpa is None: + Lpq = self.ao2mo() + self.Lpi = numpy.ascontiguousarray(Lpq[:, :, :self.nocc_act]) + self.Lpa = numpy.ascontiguousarray(Lpq[:, :, self.nocc_act:]) + logger.timer(self, "GppRPA integral transformation:", *cput0) + self.dump_flags() + kernel(pprpa=self) + logger.timer(self, "GppRPA Davidson:", *cput0) + + self.exci_t = self.exci.copy() + self.xy_t = self.xy.copy() + return + \ No newline at end of file diff --git a/pyscf/pprpa/rpprpa_davidson.py b/pyscf/pprpa/rpprpa_davidson.py index d0a438195..98f5e45ff 100644 --- a/pyscf/pprpa/rpprpa_davidson.py +++ b/pyscf/pprpa/rpprpa_davidson.py @@ -327,10 +327,10 @@ def pprpa_expand_space( Args: pprpa (ppRPA_Davidson): ppRPA_Davidson object. first_state (int): index of first particle-particle or hole-hole state. - tri_vec (double ndarray): trial vector. + tri_vec (double/complex ndarray): trial vector. tri_vec_sig (int array): signature of trial vector. - mv_prod (double ndarray): product matrix of ppRPA matrix and trial vector. - v_tri (double ndarray): eigenvector of subspace matrix. + mv_prod (double/complex ndarray): product matrix of ppRPA matrix and trial vector. + v_tri (double/complex ndarray): eigenvector of subspace matrix. Returns: conv (bool): if Davidson algorithm is converged. @@ -394,7 +394,7 @@ def pprpa_expand_space( for ivec in range(ntri): # compute product between new vector and old vector - inp = -inner_product(residue[iroot], tri_vec[ivec], pprpa.oo_dim) + inp = -inner_product(tri_vec[ivec].conj(), residue[iroot], pprpa.oo_dim) # eliminate parallel part if tri_vec_sig[ivec] < 0: inp = -inp @@ -402,10 +402,10 @@ def pprpa_expand_space( # add a new trial vector if len(residue[iroot][abs(residue[iroot]) > residue_thresh]) > 0: - assert ntri < max_vec, ( + assert ntri <= max_vec, ( "ppRPA Davidson expansion failed! ntri %d exceeds max_vec %d!" % (ntri, max_vec)) - inp = inner_product(residue[iroot], residue[iroot], pprpa.oo_dim) + inp = inner_product(residue[iroot].conj(), residue[iroot], pprpa.oo_dim) tri_vec_sig[ntri] = 1 if inp > 0 else -1 tri_vec[ntri] = residue[iroot] / numpy.sqrt(abs(inp)) ntri = ntri + 1 @@ -663,7 +663,7 @@ def check_memory(self): return def kernel(self, multi): - """Run ppRPA direct diagonalization. + """Run ppRPA davidson diagonalization. Args: multi (char): multiplicity. diff --git a/pyscf/pprpa/rpprpa_direct.py b/pyscf/pprpa/rpprpa_direct.py index e222245db..dee8c0864 100644 --- a/pyscf/pprpa/rpprpa_direct.py +++ b/pyscf/pprpa/rpprpa_direct.py @@ -284,7 +284,7 @@ def get_chemical_potential(nocc, mo_energy): return mu -def ao2mo(pprpa): +def ao2mo(pprpa, full_mo=False): """Get three-center density-fitting matrix in MO active space. Args: @@ -296,11 +296,15 @@ def ao2mo(pprpa): mf = pprpa._scf mo_coeff = mf.mo_coeff nocc = pprpa.nocc - nocc_act, nvir_act, nmo_act = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act nao = mo_coeff.shape[0] mo = numpy.asarray(mo_coeff, order='F') - ijslice = (nocc-nocc_act, nocc+nvir_act, nocc-nocc_act, nocc+nvir_act) + if full_mo: + nocc_act, nvir_act, nmo_act = nocc, mo.shape[1] - nocc, mo.shape[1] + ijslice = (0, mo.shape[1], 0, mo.shape[1]) + else: + nocc_act, nvir_act, nmo_act = pprpa.nocc_act, pprpa.nvir_act, pprpa.nmo_act + ijslice = (nocc-nocc_act, nocc+nvir_act, nocc-nocc_act, nocc+nvir_act) if isinstance(mf, (scf.rhf.RHF, dft.rks.RKS)): # molecule @@ -367,7 +371,7 @@ def pprpa_orthonormalize_eigenvector(multi, nocc, exci, xy): multi (string): multiplicity. nocc (int): number of occupied orbitals. exci (double array): ppRPA eigenvalue. - xy (double ndarray): ppRPA eigenvector. + xy (double/complex ndarray): ppRPA eigenvector. """ nroot = xy.shape[0] @@ -379,18 +383,18 @@ def pprpa_orthonormalize_eigenvector(multi, nocc, exci, xy): # determine the vector is pp or hh sig = numpy.zeros(shape=[nroot], dtype=numpy.double) for i in range(nroot): - sig[i] = 1 if inner_product(xy[i], xy[i], oo_dim) > 0 else -1 + sig[i] = 1 if inner_product(xy[i].conj(), xy[i], oo_dim).real > 0 else -1 # eliminate parallel component for i in range(nroot): for j in range(i): if abs(exci[i] - exci[j]) < 1.0e-7: - inp = inner_product(xy[i], xy[j], oo_dim) + inp = inner_product(xy[j].conj(), xy[i], oo_dim) xy[i] -= sig[j] * xy[j] * inp # normalize for i in range(nroot): - inp = inner_product(xy[i], xy[i], oo_dim) + inp = inner_product(xy[i].conj(), xy[i], oo_dim).real inp = numpy.sqrt(abs(inp)) xy[i] /= inp From a7e97ca7e32a77f3b095f4b180a705ad1f5a00f3 Mon Sep 17 00:00:00 2001 From: Chaoqun Zhang Date: Wed, 21 Jan 2026 16:13:25 -0500 Subject: [PATCH 2/3] fix code style --- pyscf/grad/rpprpa.py | 11 ++++++----- pyscf/pprpa/gpprpa_davidson.py | 9 ++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/pyscf/grad/rpprpa.py b/pyscf/grad/rpprpa.py index b3148ebfa..dfc1c3f0c 100644 --- a/pyscf/grad/rpprpa.py +++ b/pyscf/grad/rpprpa.py @@ -45,7 +45,7 @@ def grad_elec(pprpa_grad, xy, mult, atmlst=None): # To my opinion, aux_response should always be True for DFHF aux_response = mf_grad.auxbasis_response else: - logger.warn(pprpa, + logger.warn(pprpa, 'The analytical gradient of the ppRPA must be used with the density\n\ fitting mean-field method. The calculation will proceed but the analytical\n\ gradient is no longer exact (does NOT agree with numerical gradients).') @@ -101,7 +101,7 @@ def grad_elec(pprpa_grad, xy, mult, atmlst=None): vhf = lib.tag_array(vhf, aux=vhf_aux) else: vxc, vjk = get_veff_rks(mf_grad, mol, (dm0_hf, dm0)) - + vjk[1] += _contract_xc_kernel(mf, mf.xc, dm0, None, False, False, True)[0][1:]*0.5 aoslices = mol.aoslice_by_atom() @@ -421,7 +421,7 @@ def contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ slice_i = choose_slice('i', nfrozen_occ, nocc, nvir, nfrozen_vir) slice_a = choose_slice('a', nfrozen_occ, nocc, nvir, nfrozen_vir) naux = Lpq_full.shape[0] - + # Special cases for TDA if label1 == 'p': if nocc == 0: @@ -496,7 +496,8 @@ def contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ # slow (more copies) but more readable version # out = np.concatenate(( # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "i"), - # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "a")), axis=1) + # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "a")), + # axis=1) out = np.zeros((n1, nocc + nvir), dtype=occ_y_mat.dtype) L1i = np.ascontiguousarray(Lpq_full[:, slice1, slice_i]).reshape(-1, nocc) # (Pt,k) Lia = np.ascontiguousarray(Lpq_full[:, slice_i, slice_a]).reshape(-1, nvir).conj() # (P,l,b)*->(Pl,b) @@ -876,4 +877,4 @@ def kernel(self, xy=None, state=None, mult=None, atmlst=None): from pyscf.pprpa.rpprpa_direct import RppRPADirect from pyscf.pprpa.rpprpa_davidson import RppRPADavidson -RppRPADirect.Gradients = RppRPADavidson.Gradients = lib.class_as_method(Gradients) \ No newline at end of file +RppRPADirect.Gradients = RppRPADavidson.Gradients = lib.class_as_method(Gradients) diff --git a/pyscf/pprpa/gpprpa_davidson.py b/pyscf/pprpa/gpprpa_davidson.py index 1032a03c7..ac677b8de 100644 --- a/pyscf/pprpa/gpprpa_davidson.py +++ b/pyscf/pprpa/gpprpa_davidson.py @@ -258,7 +258,7 @@ def ao2mo(pprpa, full_mo=False): Lpq += 1j * nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') mo_tmp = numpy.asarray(numpy.hstack((mo_b.imag, -mo_b.real)), order='F') Lpq += 1j * nr_e2(pprpa.with_df._cderi, mo_tmp, ijslice, aosym='s2') - + return Lpq.reshape(naux, nmo_act, nmo_act) else: logger.warn(pprpa, 'Memory may not be enough!') @@ -267,7 +267,7 @@ def ao2mo(pprpa, full_mo=False): raise NotImplementedError("GppRPA for periodic systems is not implemented yet.") else: raise NotImplementedError("GppRPA only supports GHF/GKS reference now.") - + def complex_matrix_norm(matrix): """Calculate the Complex norm of each matrix element. @@ -375,8 +375,8 @@ def __init__( super().__init__( mf, nocc_act, nvir_act, channel, nroot, max_vec, max_iter, residue_thresh, print_thresh, auxbasis) - # Most of the GPPRPA formulations are just the complex-valued - # extention of those in triplet RPPRA, since they both use + # Most of the GPPRPA formulations are just the complex-valued + # extention of those in triplet RPPRA, since they both use # the antisymmetric two-electron integrals. self.multi = "t" @@ -415,4 +415,3 @@ def kernel(self): self.exci_t = self.exci.copy() self.xy_t = self.xy.copy() return - \ No newline at end of file From a4506297a6878e17e2ea459eac0b62ec4593ec8d Mon Sep 17 00:00:00 2001 From: Chaoqun Zhang Date: Sun, 25 Jan 2026 17:58:26 -0500 Subject: [PATCH 3/3] use pyscf.data.nist HARTREE2EV to keep consistency, fix code style --- pyscf/grad/rpprpa.py | 2 +- pyscf/pprpa/gpprpa_davidson.py | 4 ++-- pyscf/pprpa/rpprpa_davidson.py | 3 ++- pyscf/pprpa/rpprpa_direct.py | 3 ++- pyscf/pprpa/upprpa_direct.py | 3 ++- 5 files changed, 9 insertions(+), 6 deletions(-) diff --git a/pyscf/grad/rpprpa.py b/pyscf/grad/rpprpa.py index dfc1c3f0c..6f47e6e10 100644 --- a/pyscf/grad/rpprpa.py +++ b/pyscf/grad/rpprpa.py @@ -496,7 +496,7 @@ def contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ # slow (more copies) but more readable version # out = np.concatenate(( # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "i"), - # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "a")), + # contraction_2rdm_Lpq(occ_y_mat, vir_x_mat, Lpq_full, nocc, nvir, nfrozen_occ, nfrozen_vir, label1, "a")), # axis=1) out = np.zeros((n1, nocc + nvir), dtype=occ_y_mat.dtype) L1i = np.ascontiguousarray(Lpq_full[:, slice1, slice_i]).reshape(-1, nocc) # (Pt,k) diff --git a/pyscf/pprpa/gpprpa_davidson.py b/pyscf/pprpa/gpprpa_davidson.py index ac677b8de..ebdf833c4 100644 --- a/pyscf/pprpa/gpprpa_davidson.py +++ b/pyscf/pprpa/gpprpa_davidson.py @@ -4,7 +4,7 @@ from pyscf.lib import logger, current_memory from .rpprpa_davidson import RppRPADavidson, pprpa_get_trial_vector, pprpa_expand_space from .rpprpa_direct import pprpa_print_a_pair, inner_product, pprpa_orthonormalize_eigenvector - +from pyscf.data.nist import HARTREE2EV def kernel(pprpa): @@ -302,7 +302,7 @@ def pprpa_davidson_print_eigenvector(pprpa, exci0, exci, xy): tri_row_v, tri_col_v = numpy.tril_indices(nvir_act, -1) nroot = len(exci) - au2ev = 27.211386 + au2ev = HARTREE2EV if pprpa.channel == "pp": for iroot in range(nroot): logger.info(pprpa, "#%-d excitation: exci= %-12.4f eV 2e= %-12.4f eV", diff --git a/pyscf/pprpa/rpprpa_davidson.py b/pyscf/pprpa/rpprpa_davidson.py index 98f5e45ff..ade769be2 100644 --- a/pyscf/pprpa/rpprpa_davidson.py +++ b/pyscf/pprpa/rpprpa_davidson.py @@ -22,6 +22,7 @@ from pyscf.lib import einsum, logger, StreamObject from pyscf.mp.mp2 import get_nocc, get_nmo +from pyscf.data.nist import HARTREE2EV from pyscf.pprpa.rpprpa_direct import ao2mo, inner_product, ij2index, \ get_chemical_potential, pprpa_orthonormalize_eigenvector, pprpa_print_a_pair @@ -440,7 +441,7 @@ def pprpa_davidson_print_eigenvector(pprpa, multi, exci0, exci, xy): tri_row_v, tri_col_v = numpy.tril_indices(nvir_act, is_singlet-1) nroot = len(exci) - au2ev = 27.211386 + au2ev = HARTREE2EV if pprpa.channel == "pp": for iroot in range(nroot): logger.info(pprpa, "#%-d %s excitation: exci= %-12.4f eV 2e= %-12.4f eV", diff --git a/pyscf/pprpa/rpprpa_direct.py b/pyscf/pprpa/rpprpa_direct.py index dee8c0864..3240bcd54 100644 --- a/pyscf/pprpa/rpprpa_direct.py +++ b/pyscf/pprpa/rpprpa_direct.py @@ -26,6 +26,7 @@ from pyscf.mp.mp2 import get_nocc, get_nmo from pyscf.pbc.df.fft_ao2mo import _format_kpts from pyscf.scf.hf import KohnShamDFT +from pyscf.data.nist import HARTREE2EV def diagonalize_pprpa_singlet(nocc, mo_energy, Lpq, mu=None): @@ -441,7 +442,7 @@ def pprpa_print_direct_eigenvector(pprpa, multi, exci0, exci, xy): tri_row_o, tri_col_o = numpy.tril_indices(nocc_act, is_singlet-1) tri_row_v, tri_col_v = numpy.tril_indices(nvir_act, is_singlet-1) - au2ev = 27.211386 + au2ev = HARTREE2EV for istate in range(min(pprpa.hh_state, oo_dim)): logger.info(pprpa, "#%-d %s de-excitation: exci= %-12.4f eV 2e= %-12.4f eV", diff --git a/pyscf/pprpa/upprpa_direct.py b/pyscf/pprpa/upprpa_direct.py index fdcb7d8d6..d2cbf38a1 100644 --- a/pyscf/pprpa/upprpa_direct.py +++ b/pyscf/pprpa/upprpa_direct.py @@ -24,6 +24,7 @@ from pyscf.mp.ump2 import get_nocc, get_nmo from pyscf.pbc.df.fft_ao2mo import _format_kpts from pyscf.ao2mo._ao2mo import nr_e2 +from pyscf.data.nist import HARTREE2EV from pyscf.scf.hf import KohnShamDFT from pyscf.pprpa.rpprpa_direct import inner_product, ij2index, \ pprpa_orthonormalize_eigenvector, \ @@ -343,7 +344,7 @@ def pprpa_print_direct_eigenvector( tri_row_o, tri_col_o = numpy.tril_indices(nocc, -1) tri_row_v, tri_col_v = numpy.tril_indices(nvir, -1) - au2ev = 27.211386 + au2ev = HARTREE2EV # =====================> two-electron removal <====================== for istate in range(min(hh_state, oo_dim)):