diff --git a/openms/lib/boson.py b/openms/lib/boson.py index a691c42..e8ce175 100644 --- a/openms/lib/boson.py +++ b/openms/lib/boson.py @@ -24,13 +24,20 @@ from pyscf.lib import logger from openms.lib.ov_blocks import one_e_blocks, block_diag from openms.lib.ov_blocks import two_e_blocks, two_e_blocks_full -from openms.mqed import qedhf, scqedhf, vtqedhf def get_bosonic_Ham(nmodes, nboson_states, omega, za, Fa): r"""Construct Bosonic Hamiltonian in different representation after integrating out the electronic DOF. + TODO: Finish constructing the Hb matrix in different representation: + + - Fock/CS representation. + - FIXME: squeezed coherent state :math:`S(F)D(z)\ket{0}` + - FIXME: Displaced squeezed state :math:`D(z)S(F)\ket{0}` + - FIXME: Squeezed LF :math:`S(F) D(\hat{f})\ket{0}` + - FIXME: LF squeezed :math:`D(\hat{f})S(F) \ket{0}` + Bosonic Hamiltonain after integrating out the electronic DOF: .. math:: @@ -55,12 +62,11 @@ def get_bosonic_Ham(nmodes, nboson_states, omega, za, Fa): .. math:: - H_{VSQ} = &\omega (\cosh(2r)\omega_\alpha b^\dagger_\alpha b_\alpha + \sinh^2(r) + H_{VSQ} = & \omega_\alpha[\cosh(2r) b^\dagger_\alpha b_\alpha + \sinh^2(r)] -\frac{1}{2}\sinh(2r)[b^2_\alpha + b^{\dagger 2}_\alpha] \\ & + e^{-r} \sqrt{\frac{\omega_\alpha}{2}}\langle\Delta\lambda_\alpha\cdot\boldsymbol{D} - \rangle (b^\dagger_\alpha + b_\alpha) - + TBA. + \rangle (b^\dagger_\alpha + b_\alpha). """ boson_size = sum(nboson_states) @@ -149,25 +155,25 @@ def get_quadrupole_ao(mol, add_nuc_dipole=True, origin_shift=None): Diagonal elements correspond to: - quadrupole_ao[0] = :math:`\bm{Q}_{XX}` + quadrupole_ao[0] = :math:`\boldsymbol{Q}_{XX}` - quadrupole_ao[4] = :math:`\bm{Q}_{YY}` + quadrupole_ao[4] = :math:`\boldsymbol{Q}_{YY}` - quadrupole_ao[8] = :math:`\bm{Q}_{ZZ}` + quadrupole_ao[8] = :math:`\boldsymbol{Q}_{ZZ}` The lower- and upper-triangles are related as: - quadrupole_ao[1] = :math:`\bm{Q}_{YX}` + quadrupole_ao[1] = :math:`\boldsymbol{Q}_{YX}` - quadrupole_ao[2] = :math:`\bm{Q}_{ZX}` + quadrupole_ao[2] = :math:`\boldsymbol{Q}_{ZX}` - quadrupole_ao[3] = :math:`(\bm{Q}_{YX})^T` + quadrupole_ao[3] = :math:`(\boldsymbol{Q}_{YX})^T` - quadrupole_ao[5] = :math:`\bm{Q}_{ZY}` + quadrupole_ao[5] = :math:`\boldsymbol{Q}_{ZY}` - quadrupole_ao[6] = :math:`(\bm{Q}_{ZX})^T` + quadrupole_ao[6] = :math:`(\boldsymbol{Q}_{ZX})^T` - quadrupole_ao[7] = :math:`(\bm{Q}_{ZY})^T` + quadrupole_ao[7] = :math:`(\boldsymbol{Q}_{ZY})^T` Parameters ---------- @@ -234,14 +240,14 @@ class Boson(object): mol : :class:`MoleBase ` PySCF molecule object. omega : float, :class:`list[floats] `, :class:`~numpy.ndarray` - Bosonic frequencies, :math:`\om_\al`. + Bosonic frequencies, :math:`\omega_\alpha`. vec : list, :class:`~numpy.ndarray` - Bosonic mode vectors, :math:`\bm{e}_\al`. + Bosonic mode vectors, :math:`\boldsymbol{e}_\alpha`. Keyword Arguments ----------------- gfac : float, :class:`list[floats] `, :class:`~numpy.ndarray` - Coupling constants, :math:`\la_\al`, **optional**. + Coupling constants, :math:`\lambda_\alpha`, **optional**. nboson_states : int, :class:`list[ints] `, :class:`~numpy.ndarray` Number states for each mode, **optional** ``default = [1] * nmodes``. @@ -253,11 +259,11 @@ class Boson(object): cavity_type : str Type of cavity environment. omega : :class:`~numpy.ndarray` - Bosonic frequencies, :math:`\om_\al`. + Bosonic frequencies, :math:`\omega_\alpha`. vec : :class:`~numpy.ndarray` - Bosonic mode vectors, :math:`\bm{e}_\alpha`. + Bosonic mode vectors, :math:`\boldsymbol{e}_\alpha`. gfac : :class:`~numpy.ndarray` - Coupling constants, :math:`\la_\al`. + Coupling constants, :math:`\lambda_\alpha`. nmodes : int Number of boson modes, :code:`len(omega)`. nboson_states : :class:`~numpy.ndarray` @@ -297,16 +303,16 @@ class Boson(object): Controls how OEI component is computed, **optional** ``default = True``. couplings : :class:`~numpy.ndarray` - Coupling terms, :math:`\la_\al`, + Coupling terms, :math:`\lambda_\alpha`, copy of :attr:`gfac`. couplings_bilinear : :class:`~numpy.ndarray` Bilinear coupling terms, - :math:`\sqrt{\frac{\om_\al}{2}} \la_\al`. + :math:`\sqrt{\frac{\omega_\alpha}{2}} \lambda_\alpha`. couplings_self : :class:`~numpy.ndarray` Dipole self-energy coupling terms, - :math:`\frac{1}{2} \la_\al^2`. + :math:`\frac{1}{2} \lambda_\alpha^2`. couplings_var : :class:`~numpy.ndarray` - Variational parameters, :math:`\bm{f}`. + Variational parameters, :math:`\boldsymbol{f}`. - ``Default = 0`` for each mode. - For SC-QED-HF, ``= 1`` for each ``mode``. @@ -316,7 +322,7 @@ class Boson(object): Difference between :math:`1.0` and :attr:`couplings_var`. ``Default = 1`` for each ``mode``. optimize_varf : bool - Optimize variational parameters, :math:`\bm{f}`. + Optimize variational parameters, :math:`\boldsymbol{f}`. ``Default = False``. Raises @@ -397,6 +403,16 @@ def __init__( self.omega = omega.copy() self.nmodes = omega.size + # Cavity mode photon occupation numbers + if isinstance(nboson_states, list): + if len(nboson_states) != len(omega): + raise ValueError( + "nboson_state must be an integer or a list with the same as the length of 'omega'" + ) + self.nboson_states = nboson_states + else: + self.nboson_states = [nboson_states for i in range(self.nmodes)] + # Molecular coupling strengths self.gfac = None @@ -509,6 +525,8 @@ def __init__( self.optimize_varf = True self.couplings_var = 0.5 * numpy.ones(self.nmodes) + self.add_dse = kwargs.get("add_dse", True) + self.squeezed_cs = kwargs.get("squeezed_cs", False) self.optimize_vsq = kwargs.get("optimize_vsq", False) self.squeezed_var = numpy.zeros(self.nmodes) @@ -516,6 +534,8 @@ def __init__( self.optimize_vsq = False self.squeezed_var = kwargs["squeezed_var"] + # self.polarizations = numpy.zeros((3, self.nmodes), dtype=float) #replaced by vec + # Pre-factors of PF Hamiltonian components from coupling strengths self.couplings = numpy.asarray(self.gfac, dtype=float) self.couplings_bilinear = numpy.zeros(self.nmodes, dtype=float) @@ -553,6 +573,8 @@ def __init__( # self.shift = shift # ------------------------------------ + # TODO: reorganize the following + # get molecular attributes self.spin = mol.spin self.intor = mol.intor self.symmetry = mol.symmetry @@ -569,11 +591,42 @@ def __init__( self.energy_nuc = mol.energy_nuc self.nelec = mol.nelec # ------------------------------------ + if self.verbose > 3: + #self.print_summary() + self.dump_flags() # --------------------- # General boson methods # --------------------- + def get_boson_occ(self): + r"""Return photon ``mode`` density matrix. + + Parameters + ---------- + mode : int + Index for ``mode`` stored in the object. + + Returns + ------- + :class:`~numpy.ndarray` + photon mode density matrix + """ + + nocc = numpy.zeros(self.nmodes) + for mode in range(self.nmodes): + mdim = self.nboson_states[mode] + idx = 0 + if mode > 0: + idx = sum(self.nboson_states[:mode]) + bc = self.boson_coeff[idx:idx+mdim, 0].copy() + + bc = bc * bc + n0 = numpy.array(range(self.nboson_states[mode])) + nocc[mode] = numpy.dot(bc, n0) + print('bc =', bc, '\n', n0) + return nocc + def get_boson_dm(self, mode=None): r"""Return photon ``mode`` density matrix. @@ -595,7 +648,7 @@ def get_boson_dm(self, mode=None): mdim = self.nboson_states[mode] idx = 0 if mode > 0: - idx = sum(self.nboson_states[:mode]) + mode + idx = sum(self.nboson_states[:mode]) # + mode bc = self.boson_coeff[idx:idx+mdim, 0].copy() return numpy.outer(numpy.conj(bc), bc) @@ -609,15 +662,19 @@ def update_boson_coeff(self, eref, dm): Store eigenvectors in :attr:`boson_coeff`, update photonic energy, :attr:`e_boson`. + Note: self.z_alpha can be different frm :math:`Tr[gD]` since self.z_alpha can + be a fixed input variable. + .. math:: - \hat{H}_\tt{photon} - &= \sum_\al \hat{H}^\al_\tt{photon} \\ - \hat{H}^\al_\tt{photon} - &= \om_\al \cb{\al} \ab{\al} - + \sqrt{\frac{\om_\al}{2}} \bm{e}_\al + \hat{H}_\text{photon} + &= \sum_\alpha \hat{H}^\alpha_\text{photon} \\ + \hat{H}^\alpha_\text{photon} + &= \omega_\alpha \hat{b}^{\dagger}_{\alpha} \hat{b}_{\alpha} + + \sqrt{\frac{\omega_\alpha}{2}} \boldsymbol{e}_\alpha \cdot \sum_{\mu\nu} \rho_{\mu\nu} - \cdot \mel*{\mu}{\la_\al \cdot \hat{D}}{\nu} (\cb{\al} + \ab{\al}) - where :math:`\bm{\rho}` is the provided electronic density matrix, ``dm``. + \cdot \bra{\mu}{\lambda_\alpha \cdot \hat{D}}\ket{\nu} (\hat{b}^{\dagger}_{\alpha} + \hat{b}_{\alpha}) + + where :math:`\boldsymbol{\rho}` is the provided electronic density matrix, ``dm``. -------------------------------------------------------------------- @@ -625,22 +682,23 @@ def update_boson_coeff(self, eref, dm): states as: .. math:: - \cb{} \ab{} \ket{m} &= m \ket{m} \\ - \cb{} \ket{m} &= \sqrt{m+1} \ket{m+1} \\ - \ab{} \ket{m} &= \sqrt{m} \ket{m-1} + \hat{b}^{\dagger} \hat{b} \ket{m} &= m \ket{m} \\ + \hat{b}^{\dagger} \ket{m} &= \sqrt{m+1} \ket{m+1} \\ + \hat{b} \ket{m} &= \sqrt{m} \ket{m-1} + In the Fock/particle number basis, the off-diagonal bilinear interaction terms are: .. math:: - \mel*{n}{\cb{} + \ab{}}{m} - &= \mel*{n}{\cb{}}{m} + \mel*{n}{\ab{}}{m} \\ - &= \mel*{n}{\sqrt{m+1}}{m+1} + \mel*{n}{\sqrt{m}}{m-1} \\ - &= \d{m+1}{n} \sqrt{m+1} + \d{m-1}{n} \sqrt{m} + \bra{n}{\hat{b}^{\dagger} + \hat{b}}\ket{m} + &= \bra{n}{\hat{b}^{\dagger}}\ket{m} + \bra{n}{\hat{b}}\ket{m} \\ + &= \bra{n}{\sqrt{m+1}}\ket{m+1} + \bra{n}{\sqrt{m}}\ket{m-1} \\ + &= \delta_{m+1, n} \sqrt{m+1} + \delta_{m-1, n} \sqrt{m} .. note:: In coherent-state representation, photonic Hamiltonian has only diagonal terms. - :math:`\hat{H}_{\tt{photon}}=\sum_{\al}\om_{\al}\cb{\al}\ab{\al}` + :math:`\hat{H}_{\text{photon}}=\sum_{\alpha}\omega_{\alpha}\hat{b}^{\dagger}_{\alpha}\hat{b}_{\alpha}` In the VSQ representation: (TBA) @@ -755,7 +813,7 @@ def update_mean_field(self, mf, **kwargs): Keyword Arguments ----------------- couplings_var : float, :class:`list[floats] `, :class:`~numpy.ndarray` - Variational parameter values, :math:`f_\al`. + Variational parameter values, :math:`f_\alpha`. **Optional** optimize_varf : bool Whether to optimize variational parameters. @@ -870,7 +928,8 @@ def dump_flags(self): if self.origin_shift is not None: logger.info(self, 'origin_shift = %.3f %.3f %.3f', *(self.origin_shift)) logger.info(self, 'optimize_varf = %s', self.optimize_varf) - logger.info(self, 'use_cs = %s', self.use_cs) + logger.info(self, 'use_cs = %s', self.use_cs) + logger.info(self, 'add_dse = %s', self.add_dse) logger.info(self, 'nmodes = %d\n', self.nmodes) for a in range(self.nmodes): logger.info(self, '------ cavity mode #%s ------', (a + 1)) @@ -893,7 +952,7 @@ def dump_flags(self): # ---------------------- def update_cs(self): - r"""Template method to update :math:`z_\al` values in CS representation.""" + r"""Template method to update :math:`z_\alpha` values in CS representation.""" raise NotImplementedError("Subclasses must implement this method.") def get_gmat_so(self): @@ -932,6 +991,8 @@ def get_dipole_ao(self): r"""Template method to construct dipole moment matrix in AO basis.""" raise NotImplementedError("Subclasses must implement this method.") + def get_geb_ao_1der(self, mode): + raise NotImplementedError def construct_g_dse_JK(self): r""" @@ -983,14 +1044,14 @@ def update_cs(self, dm): .. warning:: Relate to the working equations, found in the theory, the values of :math:`z_\alpha` computed here differ by - a factor of :math:`-\frac{1}{\sqrt{2\om_\al}}`. + a factor of :math:`-\frac{1}{\sqrt{2\omega_\alpha}}`. """ if self.use_cs: self.z_alpha = lib.einsum("pq, Xpq ->X", dm, self.gmat) - def add_oei_ao(self, dm, s1e=None, residue=False): + def add_oei_ao(self, dm, s1e=None, residue=False, compute_grad=False): r"""Compute QED-RHF boson-mediated 1e- integrals. return DSE-mediated oei.. This is universal for bare HF or QED-HF. @@ -1013,24 +1074,26 @@ def add_oei_ao(self, dm, s1e=None, residue=False): from the DSE term of the PF Hamiltonian: .. math:: - \bm{h}_{\tt{DSE}} &= \bm{h}_{\tt{bare}} - \frac{1}{2} - \sum_\al \sum_{\mu\nu} \rho_{\mu\nu} - \cdot \bm{\tilde{q}}^\al_{\mu\nu} \\ - \bm{\tilde{q}}^\al_{\mu\nu} &= \la_{\al}^2 - \cdot \bm{Q}^\al_{\mu\nu} - where :math:`\bm{Q}^\al` is the polarized quadrupole moment matrix + \boldsymbol{h}_{\tt{DSE}} &= \boldsymbol{h}_{\tt{bare}} - \frac{1}{2} + \sum_\alpha \sum_{\mu\nu} \rho_{\mu\nu} + \cdot \boldsymbol{\tilde{q}}^\alpha_{\mu\nu} \\ + \boldsymbol{\tilde{q}}^\alpha_{\mu\nu} &= \lambda_{\alpha}^2 + \cdot \boldsymbol{Q}^\alpha_{\mu\nu} + + where :math:`\boldsymbol{Q}^\alpha` is the polarized quadrupole moment matrix in the AO basis. In CS representation, when :attr:`use_cs` ``= True``, additional DSE-mediated contribution to the OEI is included: .. math:: - \bm{h}_{\tt{CS}} &= \bm{h}_{\tt{DSE}} - \bm{\tilde{d}}^\al - \sum_\al \sum_{\mu\nu} \rho_{\mu\nu} - \cdot \mel*{\mu}{\hat{D}}{\nu} \\ - \bm{\tilde{d}}^\al_{\mu\nu} &= \la_\al - \cdot \bm{D}^\al_{\mu\nu} - where :math:`\bm{D}^\al` is the polarized dipole moment matrix in + \boldsymbol{h}_{\tt{CS}} &= \boldsymbol{h}_{\tt{DSE}} - \boldsymbol{\tilde{d}}^\alpha + \sum_\alpha \sum_{\mu\nu} \rho_{\mu\nu} + \cdot \bra{\mu}{\hat{D}}\ket{\nu} \\ + \boldsymbol{\tilde{d}}^\alpha_{\mu\nu} &= \lambda_\alpha + \cdot \boldsymbol{D}^\alpha_{\mu\nu} + + where :math:`\boldsymbol{D}^\alpha` is the polarized dipole moment matrix in the AO basis. Also in CS representation, :math:`E_{\tt{QED-HF}}` includes DSE @@ -1038,31 +1101,33 @@ def add_oei_ao(self, dm, s1e=None, residue=False): .. math:: E_{\tt{QED-HF}} = E_{\tt{HF}} - + \frac{1}{2} \sum_\al + + \frac{1}{2} \sum_\alpha \langle - [ \la_\al \cdot - (\hat{D} - \ev*{\hat{D}}_{\mu\nu}) + [ \lambda_\alpha \cdot + (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu}) ]^2 \rangle + that can also be included via the OEI, since: .. math:: \langle - [ \la_\al - \cdot (\hat{D} - \ev*{\hat{D}}_{\mu\nu})]^2 + [ \lambda_\alpha + \cdot (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu})]^2 \rangle - &= \langle [ \la_\al - \cdot (\hat{D} - \ev*{\hat{D}}_{\mu\nu})]^2 + &= \langle [ \lambda_\alpha + \cdot (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu})]^2 \rangle - \cdot \frac{\Tr[\bm{S} \cdot \bm{\rho}]} + \cdot \frac{\text{Tr}[\boldsymbol{S} \cdot \boldsymbol{\rho}]} {N_{\tt{elec}}} \\ &= \tt{Tr} \left[ - \frac{[ \la_\al \cdot - (\hat{D} - \ev*{\hat{D}}_{\mu\nu})]^2} - {N_{\tt{elec}}} \bm{S} \cdot \bm{\rho} + \frac{[ \lambda_\alpha \cdot + (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu})]^2} + {N_{\tt{elec}}} \boldsymbol{S} \cdot \boldsymbol{\rho} \right] - where :math:`\bm{\rho}` is the provided electronic density matrix, + + where :math:`\boldsymbol{\rho}` is the provided electronic density matrix, ``dm``. Parameters @@ -1095,24 +1160,41 @@ def add_oei_ao(self, dm, s1e=None, residue=False): gvar2 = numpy.ones(self.nmodes) if residue: - gvar2 = self.couplings_res**2 # element-wise + if self.add_dse: + # (1-f) ^2 + gvar2 = self.couplings_res**2 # element-wise + else: + # f^2 - 2f = f * (f-2) + gvar2 = self.couplings_var * (self.couplings_var - 2.0) + dg = - 2.00 * self.couplings_res # 2(f-1) oei = - lib.einsum("Xpq, X->pq", self.gmat, gvar2 * self.z_alpha) - if self.complete_basis: - oei -= 0.5 * lib.einsum("Xpq, X->pq", self.q_lambda_ao, gvar2) - else: - s_eval, s_evec = scipy.linalg.eigh(s1e) - idx = s_eval > 1e-15 - s_inv = numpy.dot(s_evec[:, idx] / s_eval[idx], s_evec[:, idx].conj().T) + if compute_grad: + oei_grad = - lib.einsum("Xpq, X->pq", self.gmat, dg * self.z_alpha) + + # the following part comes from the DSE + if self.add_dse: + if self.complete_basis: + tmp_term = - self.q_lambda_ao + # oei -= 0.5 * lib.einsum("Xpq, X->pq", self.q_lambda_ao, gvar2) + else: + s_eval, s_evec = scipy.linalg.eigh(s1e) + idx = s_eval > 1e-15 + s_inv = numpy.dot(s_evec[:, idx] / s_eval[idx], s_evec[:, idx].conj().T) + tmp_term = lib.einsum("Xpr, rs, Xsq-> Xpq", self.gmat, s_inv, self.gmat) - tmp_term = lib.einsum("Xpr, rs, Xsq-> Xpq", self.gmat, s_inv, self.gmat) oei += 0.5 * lib.einsum("X, Xpq-> pq", gvar2, tmp_term) + if compute_grad: + oei_grad += 0.5 * lib.einsum("X, Xpq-> pq", dg, tmp_term) del (tmp_term) + # FIXME: this will cause problem to VT/VSQ gradient with respect to f as e_boson is not function of f # DSE + boson energy (beyond |0> state) z2s = (0.5 * numpy.sum(self.z_alpha**2 * gvar2) + self.e_boson)* s1e/self._mol.nelectron # z2s = (0.5 * numpy.sum(self.z_alpha**2 * gvar2))* s1e/self._mol.nelectron oei += z2s + if compute_grad: + oei_grad += 0.5 * numpy.sum(self.z_alpha**2 * dg) * s1e/self._mol.nelectron # bilinear term # off-diaognal (photonic) block @@ -1120,7 +1202,7 @@ def add_oei_ao(self, dm, s1e=None, residue=False): for imode in range(self.nmodes): shift = numpy.eye(self.gmat.shape[1]) * self.z_alpha[imode] #/ mol.nelectron shift = s1e * self.z_alpha[imode] - gtmp = (self.gmat[imode] - shift) * self.couplings_res[imode] + gtmp = (self.gmat[imode] - shift) gtmp *= numpy.sqrt(0.5*self.omega[imode]) ci = self.boson_coeff[idx : idx + self.nboson_states[imode], 0] @@ -1129,9 +1211,15 @@ def add_oei_ao(self, dm, s1e=None, residue=False): h_od = numpy.diag(numpy.sqrt(numpy.arange(1, mdim)), k = 1) \ + numpy.diag(numpy.sqrt(numpy.arange(1, mdim)), k = -1) ph_exp_val = numpy.sum(h_od * pdm) - oei += gtmp * ph_exp_val + + oei += ph_exp_val * gtmp * self.couplings_res[imode] + if compute_grad: + # d(1-f) /df = -1 + oei_grad -= ph_exp_val * gtmp idx += self.nboson_states[imode] + if compute_grad: + return oei, oei_grad return oei @@ -1147,24 +1235,26 @@ def get_dse_hcore(self, dm=None, s1e=None, residue=False): from the DSE term of the PF Hamiltonian: .. math:: - \bm{h}_{\tt{DSE}} &= \bm{h}_{\tt{bare}} - \frac{1}{2} - \sum_\al \sum_{\mu\nu} \rho_{\mu\nu} - \cdot \bm{\tilde{q}}^\al_{\mu\nu} \\ - \bm{\tilde{q}}^\al_{\mu\nu} &= \la_{\al}^2 - \cdot \bm{Q}^\al_{\mu\nu} - where :math:`\bm{Q}^\al` is the polarized quadrupole moment matrix + \boldsymbol{h}_{\tt{DSE}} &= \boldsymbol{h}_{\tt{bare}} - \frac{1}{2} + \sum_\alpha \sum_{\mu\nu} \rho_{\mu\nu} + \cdot \boldsymbol{\tilde{q}}^\alpha_{\mu\nu} \\ + \boldsymbol{\tilde{q}}^\alpha_{\mu\nu} &= \lambda_{\alpha}^2 + \cdot \boldsymbol{Q}^\alpha_{\mu\nu} + + where :math:`\boldsymbol{Q}^\alpha` is the polarized quadrupole moment matrix in the AO basis. In CS representation, when :attr:`use_cs` ``= True``, additional DSE-mediated contribution to the OEI is included: .. math:: - \bm{h}_{\tt{CS}} &= \bm{h}_{\tt{DSE}} - \bm{\tilde{d}}^\al - \sum_\al \sum_{\mu\nu} \rho_{\mu\nu} - \cdot \mel*{\mu}{\hat{D}}{\nu} \\ - \bm{\tilde{d}}^\al_{\mu\nu} &= \la_\al - \cdot \bm{D}^\al_{\mu\nu} - where :math:`\bm{D}^\al` is the polarized dipole moment matrix in + \boldsymbol{h}_{\tt{CS}} &= \boldsymbol{h}_{\tt{DSE}} - \boldsymbol{\tilde{d}}^\alpha + \sum_\alpha \sum_{\mu\nu} \rho_{\mu\nu} + \cdot \bra{\mu}{\hat{D}}\ket{\nu} \\ + \boldsymbol{\tilde{d}}^\alpha_{\mu\nu} &= \lambda_\alpha + \cdot \boldsymbol{D}^\alpha_{\mu\nu} + + where :math:`\boldsymbol{D}^\alpha` is the polarized dipole moment matrix in the AO basis. Also in CS representation, :math:`E_{\tt{QED-HF}}` includes DSE @@ -1172,31 +1262,33 @@ def get_dse_hcore(self, dm=None, s1e=None, residue=False): .. math:: E_{\tt{QED-HF}} = E_{\tt{HF}} - + \frac{1}{2} \sum_\al + + \frac{1}{2} \sum_\alpha \langle - [ \la_\al \cdot - (\hat{D} - \ev*{\hat{D}}_{\mu\nu}) + [ \lambda_\alpha \cdot + (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu}) ]^2 \rangle + that can also be included via the OEI, since: .. math:: \langle - [ \la_\al - \cdot (\hat{D} - \ev*{\hat{D}}_{\mu\nu})]^2 + [ \lambda_\alpha + \cdot (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu})]^2 \rangle - &= \langle [ \la_\al - \cdot (\hat{D} - \ev*{\hat{D}}_{\mu\nu})]^2 + &= \langle [ \lambda_\alpha + \cdot (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu})]^2 \rangle - \cdot \frac{\Tr[\bm{S} \cdot \bm{\rho}]} + \cdot \frac{\text{Tr}[\boldsymbol{S} \cdot \boldsymbol{\rho}]} {N_{\tt{elec}}} \\ &= \tt{Tr} \left[ - \frac{[ \la_\al \cdot - (\hat{D} - \ev*{\hat{D}}_{\mu\nu})]^2} - {N_{\tt{elec}}} \bm{S} \cdot \bm{\rho} + \frac{[ \lambda_\alpha \cdot + (\hat{D} - \langle{\hat{D}}\rangle_{\mu\nu})]^2} + {N_{\tt{elec}}} \boldsymbol{S} \cdot \boldsymbol{\rho} \right] - where :math:`\bm{\rho}` is the provided electronic density matrix, + + where :math:`\boldsymbol{\rho}` is the provided electronic density matrix, ``dm``. Parameters @@ -1287,14 +1379,14 @@ def get_dse_jk( set to an array of zeros when not in CS representation. .. math:: - J^\tt{DSE}_{\mu\nu} &= \sum_{\la\si} \rho_{\la\si} - \sum_\al (\mu\nu|\la\si)^\al - = \sum_{\la\si} \rho_{\la\si} - \sum_\al g^\al_{\mu\nu} g^\al_{\la\si} \\ - K^\tt{DSE}_{\mu\nu} &= \sum_{\la\si} \rho_{\la\si} - \sum_\al (\mu\si|\la\nu)^\al - = \sum_{\la\si} \rho_{\la\si} - \sum_\al g^\al_{\mu\si} g^\al_{\la\nu} + J^\tt{DSE}_{\mu\nu} &= \sum_{\lambda\sigma} \rho_{\lambda\sigma} + \sum_\alpha (\mu\nu|\lambda\sigma)^\alpha + = \sum_{\lambda\sigma} \rho_{\lambda\sigma} + \sum_\alpha g^\alpha_{\mu\nu} g^\alpha_{\lambda\sigma} \\ + K^\tt{DSE}_{\mu\nu} &= \sum_{\lambda\sigma} \rho_{\lambda\sigma} + \sum_\alpha (\mu\sigma|\lambda\nu)^\alpha + = \sum_{\lambda\sigma} \rho_{\lambda\sigma} + \sum_\alpha g^\alpha_{\mu\sigma} g^\alpha_{\lambda\nu} Parameters ---------- @@ -1360,9 +1452,9 @@ def get_polarized_quadrupole_ao(self, mode): photon ``mode`` and quadrupole moment matrix :attr:`quadrupole_ao`. .. math:: - \bm{Q}_\al = \bm{e}_\al - \cdot \mel*{\mu}{\hat{Q}}{\nu} - \cdot \bm{e}_\al + \boldsymbol{Q}_\alpha = \boldsymbol{e}_\alpha + \cdot \bra{\mu}{\hat{Q}}\ket{\nu} + \cdot \boldsymbol{e}_\alpha Parameters ---------- @@ -1383,8 +1475,9 @@ def get_q_lambda_ao(self): r"""Compute dipole self-energy matrix for all modes in AO basis. .. math:: - \bm{\tilde{q}} = \sum_\al \cdot - \mel*{\mu}{(\bm{\la}_\al\cdot\hat{D})^2}{\nu} + \boldsymbol{\tilde{q}} = \sum_\alpha \cdot + \bra{\mu}{(\boldsymbol{\lambda}_\alpha\cdot\hat{D})^2}\ket{\nu} + Function only runs if :attr:`q_lambda_ao` is not ``None``. Returns @@ -1445,8 +1538,8 @@ def get_gmat_ao(self): r"""Compute interaction matrix for all modes in AO basis. .. math:: - \bm{\tilde{d}} = \sum_\al \cdot - \mel*{\mu}{\bm{\la}_\al\cdot\hat{D}}{\nu} + \boldsymbol{\tilde{d}} = \sum_\alpha \cdot + \mel*{\mu}{\boldsymbol{\la}_\alpha\cdot\hat{D}}{\nu} Function only runs if :attr:`gmat` is not ``None``. @@ -1474,8 +1567,9 @@ def get_gmatao(self): TODO: get_gmat_ao or this one will be deprecated .. math:: - \bm{\tilde{d}} = \sum_\al \cdot - \mel*{\mu}{\bm{\la}_\al\cdot\hat{D}}{\nu} + \boldsymbol{\tilde{d}} = \sum_\alpha \cdot + \bra{\mu}{\boldsymbol{\lambda}_\alpha\cdot\hat{D}}\ket{\nu} + Function only runs if :attr:`gmat` is not ``None``. @@ -1499,13 +1593,15 @@ def get_gmatao(self): def get_geb_ao(self, mode): r"""Return bilinear interaction term of ``mode`` in :class:`Boson`. + Gets the bilinear interaction term in the AO basis, i.e., the g in g(b+b^+). + The bilinear interaction term of ``mode`` is :attr:`gmat[mode] `, scaled by a frequency-dependent term: .. math:: - \sqrt{\frac{\om_\al}{2}} \bm{\tilde{d}}^\al - = \sqrt{\frac{\om_\al}{2}} - \cdot \mel*{\mu}{\bm{\la}_\al \cdot \hat{D}}{\nu} + \sqrt{\frac{\omega_\alpha}{2}} \boldsymbol{\tilde{d}}^\alpha + = \sqrt{\frac{\omega_\alpha}{2}} + \cdot \bra{\mu}{\boldsymbol{\lambda}_\alpha \cdot \hat{D}}\ket{\nu} Parameters ---------- @@ -1517,8 +1613,12 @@ def get_geb_ao(self, mode): :class:`~numpy.ndarray` Interaction matrix of mode scaled by frequency term. """ - return (numpy.sqrt(0.5 * self.omega[mode]) * self.gmat[mode]) + logger.debug(self, " construct bilinear interation term in AO") + g_eb = self.get_polarized_dipole_ao(mode) + logger.debug(self, f" Norm of gao without w {numpy.linalg.norm(g_eb)}") + g_eb *= self.couplings_bilinear[mode] + return g_eb def get_bdag_plus_b_sq_expval(self, mode): @@ -1553,7 +1653,7 @@ def get_mos(self): """ mf = self._mf - self.nmo = 2 * self.nao + self.nmo = 2 * self._mol.nao_nr() if mf.mo_coeff is None: mf.kernel() if mf.mo_occ.ndim == 1: @@ -1784,6 +1884,7 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) + self.vmat = None # def relax(self): raise NotImplementedError("Method not implemented!") @@ -1814,10 +1915,13 @@ def relax(self): hf.diis_space = 10 mf = hf.run() + qed = Photon(mol, mf=mf, omega=omega, vec=vec, gfac=gfac) + qed.kernel() + for i in range(nmodes): - g_eb = hf.qed.get_geb_ao(i) + g_eb = qed.get_geb_ao(i) print("e-ph coupling matrix of mode ", i, "is \n", g_eb) # get I - I = hf.qed.get_I() - F = hf.qed.g_fock() + I = qed.get_I() + F = qed.g_fock() diff --git a/openms/mqed/qedhf.py b/openms/mqed/qedhf.py index 149228b..a2e7443 100644 --- a/openms/mqed/qedhf.py +++ b/openms/mqed/qedhf.py @@ -15,20 +15,6 @@ # Ilia Mazin # -import numpy - -from pyscf.dft import rks -from pyscf import lib -from pyscf import scf -from pyscf.scf import hf -from pyscf.lib import logger - -import openms -from openms.lib import boson -from openms import __config__ - -TIGHT_GRAD_CONV_TOL = getattr(__config__, "TIGHT_GRAD_CONV_TOL", True) - r""" Theoretical background of QEDHF ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -39,7 +25,7 @@ .. math:: \ket{\Psi} = &\sum_n c_n \prod_\alpha e^{z_\alpha b^\dagger_\alpha - z^*_\alpha b_\alpha } \ket{HF}\otimes{n_\alpha} \\ - = &\sum_n c_n U(\mathbf{z}) \ket{HF}\otimes{n_\alhpa}. + = &\sum_n c_n U(\mathbf{z}) \ket{HF}\otimes{n_\alpha}. where :math:`z_\alpha=-\frac{\lambda_\alpha\cdot\langle\boldsymbol{D}\rangle}{\sqrt{2\omega_\alpha}}` denotes the photon displacement due to the coupling with electrons. @@ -60,10 +46,26 @@ .. math:: - E_{QEDHF}= E_{HF} + \frac{1}{2}\langle \boldsymbol{lambda}\cdot [\boldsymbol{D}-\langle \boldsymbol{D}\rangle)]^2\rangle, + E_{QEDHF}= E_{HF} + \frac{1}{2}\langle \boldsymbol{\lambda}\cdot [\boldsymbol{D}-\langle \boldsymbol{D}\rangle)]^2\rangle, """ +import numpy + +from pyscf.dft import rks +from pyscf import lib +from pyscf import scf +from pyscf.scf import hf +from pyscf.lib import logger + +import openms +from openms.lib import boson +from openms import __config__ + +TIGHT_GRAD_CONV_TOL = getattr(__config__, "TIGHT_GRAD_CONV_TOL", True) + + + def kernel( mf, conv_tol=1e-10, conv_tol_grad=None, dump_chk=True, dm0=None, callback=None, conv_check=True, **kwargs): @@ -511,7 +513,7 @@ def get_hcore(self, mol=None, dm=None): return (self.bare_h1e + self.oei) - def get_dse_jk(self, dm, residue=False): + def get_dse_jk(self, dm, residue=False, compute_grad=False): r"""DSE-mediated JK terms To be deprecated, use boson.get_dse_jk instead in the future @@ -521,6 +523,13 @@ def get_dse_jk(self, dm, residue=False): effective DSE-mediated eri is: I' = lib.einsum("Xpq, Xrs->pqrs", gmat, gmat) + + The residual DSE are + + .. math :: + H_{DSE} &= \frac{(1-f)^2}{2}(\boldsymbol{\lambda}\cdot\boldsymbol{D})^2 \textit{(with bare DSE)} \\ + &= \frac{f^2-2f}{2}(\boldsymbol{\lambda}\cdot\boldsymbol{D})^2 \textit{(without bare DSE)} \\ + """ dm_shape = dm.shape @@ -529,21 +538,40 @@ def get_dse_jk(self, dm, residue=False): n_dm = dm.shape[0] logger.debug(self, "No. of dm is %d", n_dm) - vj_dse = numpy.zeros((n_dm,nao,nao)) - vk_dse = numpy.zeros((n_dm,nao,nao)) + grad_vjk = numpy.zeros((n_dm, nao, nao)) + vj_dse = numpy.zeros((n_dm, nao, nao)) + vk_dse = numpy.zeros((n_dm, nao, nao)) gtmp = self.qed.gmat - if residue: gtmp = self.qed.gmat * self.qed.couplings_res + gvar2 = numpy.ones(self.qed.nmodes) + if residue: + if self.qed.add_dse: + gvar2 = self.qed.couplings_res**2 # element-wise + else: + gvar2 = self.qed.couplings_var * (self.qed.couplings_var - 2.0) + dg = - 2.00 * self.qed.couplings_res # 2(f-1) + for i in range(n_dm): # DSE-medaited J - scaled_mu = lib.einsum("pq, Xpq ->X", dm[i], gtmp)# <\lambada * D> - vj_dse[i] += lib.einsum("Xpq, X->pq", gtmp, scaled_mu) + scaled_mu = lib.einsum("pq, Xpq ->X", dm[i], gtmp) # <\lambada * D> + vj_dse[i] += lib.einsum("Xpq, X, X->pq", gtmp, scaled_mu, gvar2) + if compute_grad: + # FIXME: this gradient should be mode specific + grad_vjk[i] += lib.einsum("Xpq, X, X->pq", gtmp, scaled_mu, dg) # DSE-mediated K - vk_dse[i] += lib.einsum("Xpr, Xqs, rs -> pq", gtmp, gtmp, dm[i]) + # vk_dse[i] += lib.einsum("Xpr, Xqs, rs -> pq", gtmp, gtmp, dm[i]) + gdm = lib.einsum("Xqs, rs -> Xqr", gtmp, dm[i]) + vk_dse += lib.einsum("Xpr, Xqr, X-> pq", gtmp, gdm, gvar2) + if compute_grad: + # FIXME: this gradient should be mode specific + grad_vjk -= 0.5 * lib.einsum("Xpr, Xqr, X-> pq", gtmp, gdm, dg) vj = vj_dse.reshape(dm_shape) vk = vk_dse.reshape(dm_shape) + if compute_grad: + grad_vjk = grad_vjk.reshape(dm_shape) + return vj, vk, grad_vjk return vj, vk def get_jk(self, mol=None, dm=None, hermi=1, with_j=True, with_k=True, omega=None): @@ -561,9 +589,10 @@ def get_jk(self, mol=None, dm=None, hermi=1, with_j=True, with_k=True, omega=Non else: vj, vk = RHF.get_jk(self, mol, dm, hermi, with_j, with_k, omega) - vj_dse, vk_dse = self.get_dse_jk(dm) - vj += vj_dse - vk += vk_dse + if self.qed.add_dse: + vj_dse, vk_dse = self.get_dse_jk(dm) + vj += vj_dse + vk += vk_dse return vj, vk @@ -716,7 +745,7 @@ def scf(self, dm0=None, **kwargs): module. Initial density matrix guess moved here from kernel function, to allow for computing initial QED quantities that require electronic density matrix, such - as the :math:`\bm{z}_{\al}` values for calculations + as the :math:`\bm{z}_{\alpha}` values for calculations being performed with the photonic wavefunction in the coherent-state (CS) represenation. @@ -777,6 +806,15 @@ def post_kernel(self, envs): logger.note(self, f"{breakline}\n") return self + def make_boson_nocc_org(self, nfock=10): + r"""boson occ in original fock representation""" + import math + + _, _, rho_b = self.make_rdm1_org(self.mo_coeff, self.mo_occ, nfock=nfock) + nocc = numpy.dot(numpy.diagonal(rho_b), numpy.array(range(nfock))) + print("nocc =", nocc) + return nocc + def make_rdm1_org(self, mo_coeff, mo_occ, nfock=2, **kwargs): r""" """ @@ -802,7 +840,7 @@ def make_rdm1_org(self, mo_coeff, mo_occ, nfock=2, **kwargs): rho_tot[m, :, n, :] = rho_tmp rho_e = numpy.einsum("mpmq->pq", rho_tot) - rho_b = numpy.einsum("mpnp->mn", rho_tot) / numpy.trace(rho) + rho_b = numpy.einsum("mpnp->mn", rho_tot) / numpy.trace(rho_e) return rho_tot, rho_e, rho_b diff --git a/openms/mqed/scqedhf.py b/openms/mqed/scqedhf.py index 3be416b..be3570d 100644 --- a/openms/mqed/scqedhf.py +++ b/openms/mqed/scqedhf.py @@ -15,6 +15,30 @@ # Ilia Mazin # +r""" +Theoretical background of SC/VT-QEDHF methods +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +SC-QEDHF module for solving the QED Hamiltonian. The kernel is also used for VT-QEDHF. + +The WF ansatz for SC-QEDHF is: + +.. math:: + + \ket{\Phi} = e^{} \ket{\Phi_0} + +The corresponding HF Energy within SC-QEDHF formalism becomes: + +.. math:: + + E = \sum_{pq} h_{pq} D_{pq} G_{pq} + \cdots + + + +""" + + + import time import numpy from scipy import linalg @@ -38,27 +62,6 @@ FORCE_PIVOTED_CHOLESKY = getattr(__config__, "FORCE_PIVOTED_CHOLESKY", False) LINEAR_DEP_TRIGGER = getattr(__config__, "LINEAR_DEP_TRIGGER", 1e-10) -r""" -Theoretical background of SC/VT-QEDHF methods -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -SC-QEDHF module for solving the QED Hamiltonian. The kernel is also used for VT-QEDHF. - -The WF ansatz for SC-QEDHF is: - -.. math:: - - \ket{\Phi} = e^{} \ket{\Phi_0} - -The corresponding HF Energy within SC-QEDHF formalism becomes: - -.. math:: - - E = \sum_{pq} h_{pq} D_{pq} G_{pq} + \cdots - - - -""" def unitary_transform(U, A): r"U^T A U" @@ -963,6 +966,66 @@ def scf(self, dm0=None, **kwargs): return self.e_tot kernel = lib.alias(scf, alias_name='kernel') + + def get_jk_chols(self, mol=None, dm=None, hermi=1, with_j=True, with_k=True, omega=None): + r"""get jk matrix in the in DO and with chols + + .. math:: + + I_{ijkl} = I^e_{ijkl} + g_{ij}g_{kl} + + where :math:`I^e` is the pure electronci two-body integral. + the latter term counts for the photon-mediated correlations + + .. math:: + + J_{uv} &= \sum_{ls} D_{ls}(uv|ls) \\ + K_{uv} &= \sum_{ls} D_{ls}(us|lv) + + Hence, the photon-mediated part of JK is: + + .. math:: + + J^p_{uv} &= \sum_{ls} D_{ls}(uv|ls) = \sum_{ls} D_{ls} g_{uv} g_{ls} \\ + K^p_{uv} &= \sum_{ls} D_{ls} g_{us} g_{lv} + + """ + # Note the incore version, which initializes an _eri array in memory. + if mol is None: + mol = self.mol + if dm is None: + dm = self.make_rdm1() + dm = self.get_dm_do(dm, self.ao2dipole) + + vj = vk = None + + # replace the following with chols, without using eri_DO + fc_factor = self.FC_factor(self.eta, imode, onebody=False) + fc_factor *= (1.0 * self.eri_DO - 0.5 * self.eri_DO.transpose(0, 3, 2, 1)) + # G_{pqrs} * [I_{pqrs} - 0.5 * I_{ps rq} ] * rho_{rs} + # = L_{\gamma, pq} L_{\gamma, rs} * \rho_{rs} + # - L_{\gamma, pq} L_{\gamma, rs} * \rho_{rq} + + # vj = G_{pqrs} I_{pqrs} * \rho_{rs} + # = G_{pqrs} L_{\gamma, pq} L_{\gamma, rs} * \rho_{rs} + + # vk = G_{psrq} I_{psrq} * \rho_{rs} + # = G_{psrq} L_{\gamma, ps} L_{\gamma, rq} * \rho_{rs} + + # the above is the same as the following, using the following if we + # want to separate it into vj and vk + + # add dressing factor to two-body integrals (todo) + for imode in range(self.qed.nmodes): + #U = self.ao2dipole[imode] + factor = self.FC_factor(self.eta, imode, onebody=False) + eri_tmp = self.eri_DO * factor + vj, vk = hf.dot_eri_dm(eri_tmp, dm, hermi, with_j, with_k) + del eri_tmp + + return vj, vk + + # =============================== # below are deprecated functions # =============================== diff --git a/openms/mqed/vtqedhf.py b/openms/mqed/vtqedhf.py index b8022dc..c9d9b67 100644 --- a/openms/mqed/vtqedhf.py +++ b/openms/mqed/vtqedhf.py @@ -14,25 +14,6 @@ # Author: Yu Zhang # -import numpy -from scipy import linalg -import warnings - -from pyscf import lib -from pyscf.lib import logger - -import openms -from openms.mqed import scqedhf -from openms import __config__ - -TIGHT_GRAD_CONV_TOL = getattr(__config__, "TIGHT_GRAD_CONV_TOL", True) -LINEAR_DEP_THRESHOLD = getattr(__config__, "LINEAR_DEP_THRESHOLD", 1e-8) -CHOLESKY_THRESHOLD = getattr(__config__, "CHOLESKY_THRESHOLD", 1e-10) -FORCE_PIVOTED_CHOLESKY = getattr(__config__, "FORCE_PIVOTED_CHOLESKY", False) -LINEAR_DEP_TRIGGER = getattr(__config__, "LINEAR_DEP_TRIGGER", 1e-10) - - - r""" Theoretical background of VT-QEDHF methods @@ -75,6 +56,23 @@ """ +import numpy +from scipy import linalg +import warnings + +from pyscf import lib +from pyscf.lib import logger + +import openms +from openms.mqed import scqedhf +from openms import __config__ + +TIGHT_GRAD_CONV_TOL = getattr(__config__, "TIGHT_GRAD_CONV_TOL", True) +LINEAR_DEP_THRESHOLD = getattr(__config__, "LINEAR_DEP_THRESHOLD", 1e-8) +CHOLESKY_THRESHOLD = getattr(__config__, "CHOLESKY_THRESHOLD", 1e-10) +FORCE_PIVOTED_CHOLESKY = getattr(__config__, "FORCE_PIVOTED_CHOLESKY", False) +LINEAR_DEP_TRIGGER = getattr(__config__, "LINEAR_DEP_TRIGGER", 1e-10) + def newton_update(var, gradient, precond0): r"""TBA""" @@ -123,6 +121,7 @@ def __init__(self, mol, **kwargs): self.qed.update_couplings() self.vhf_dse = None + self.grad_vhf_dse = None def get_hcore(self, mol=None, dm=None, dress=True): @@ -130,9 +129,10 @@ def get_hcore(self, mol=None, dm=None, dress=True): h1e = super().get_hcore(mol, dm, dress) if dm is not None: - + # this is zero for scqedhf # self.oei = self.qed.get_dse_hcore(dm, residue=True) # this is zero for scqedhf - self.oei = self.qed.add_oei_ao(dm, residue=True) # this is zero for scqedhf + self.oei, self.grad_oei = self.qed.add_oei_ao(dm, residue=True, compute_grad=True) + return h1e + self.oei return h1e @@ -152,7 +152,8 @@ def get_veff(self, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1): # two-electron part (residue, to be tested!) # vj_dse, vk_dse = self.qed.get_dse_jk(dm, residue=True) - vj_dse, vk_dse = self.get_dse_jk(dm, residue=True) + vj_dse, vk_dse, self.grad_vhf_dse = self.get_dse_jk(dm, residue=True, compute_grad=True) + self.vhf_dse = vj_dse - 0.5 * vk_dse vhf += self.vhf_dse @@ -334,11 +335,17 @@ def grad_var_params(self, dm_do, g_DO, dm=None): if abs(1.0 - self.qed.couplings_var[0]) > 1.0e-4 and self.vhf_dse is not None: # only works for nmode == 1 - self.dse_tot = self.energy_elec(dm, self.oei, self.vhf_dse)[0] - self.vlf_grad[0] -= self.dse_tot * 2.0 / (1.0 - self.qed.couplings_var[0]) - + ## E_DSE = (1-f)^2 * original E_DSE, + ## so the gradient = -2(1-f) * original E_DSE = - E_DSE*2/(1-f) + # old code + # self.dse_tot = self.energy_elec(dm, self.oei, self.vhf_dse)[0] + # self.vlf_grad[0] -= self.dse_tot * 2.0 / (1.0 - self.qed.couplings_var[0]) + + # Now, we compute the gradient from add_oei and get_dse_jk + self.vlf_grad[0] += self.energy_elec(dm, self.grad_oei, self.grad_vhf_dse)[0] return self + def norm_var_params(self): var_norm = linalg.norm(self.eta_grad) / numpy.sqrt(self.eta.size) var_norm += linalg.norm(self.vlf_grad) / numpy.sqrt(self.vlf_grad.size) @@ -382,14 +389,14 @@ def set_var_params(self, params): fsize += etasize varsize = self.qed.couplings_var.size - if params.size > fsize: + if params.size > fsize and self.qed.optimize_varf: self.qed.couplings_var = params[fsize : fsize + varsize].reshape( self.vlf_grad.shape ) self.qed.update_couplings() fsize += varsize - if params.size > fsize: + if params.size > fsize and self.qed.optimize_vsq: self.qed.squeezed_var = params[fsize :].reshape( self.vsq_grad.shape ) @@ -419,6 +426,16 @@ def set_params(self, params, fock_shape=None): ) return f + def make_boson_nocc_org(self, nfock=10): + r"""boson occ in original fock representation""" + import math + + _, _, rho_b = self.make_rdm1_org(self.mo_coeff, self.mo_occ, nfock=nfock) + nocc = numpy.dot(numpy.diagonal(rho_b), numpy.array(range(nfock))) + print("nocc =", nocc) + return nocc + + def make_rdm1_org(self, mo_coeff, mo_occ, nfock=2, **kwargs): r"""One-particle density matrix in original AO-Fock representation @@ -457,7 +474,8 @@ def make_rdm1_org(self, mo_coeff, mo_occ, nfock=2, **kwargs): z0 = numpy.exp(-0.5 * zalpha ** 2) zm = z0 * zalpha ** m * numpy.sqrt(math.factorial(m)) - zn = z0 * (-zalpha) ** n * numpy.sqrt(math.factorial(n)) + # zn = z0 * (-zalpha) ** n * numpy.sqrt(math.factorial(n)) + zn = z0 * (zalpha) ** n * numpy.sqrt(math.factorial(n)) rho_tmp = rho_DO * numpy.outer(zm, zn) # back to AO @@ -465,7 +483,7 @@ def make_rdm1_org(self, mo_coeff, mo_occ, nfock=2, **kwargs): rho_tot[m, :, n, :] = rho_tmp rho_e = numpy.einsum("mpmq->pq", rho_tot) - rho_b = numpy.einsum("mpnp->mn", rho_tot) / numpy.trace(rho) + rho_b = numpy.einsum("mpnp->mn", rho_tot) / numpy.trace(rho_e) return rho_tot, rho_e, rho_b