#!/usr/bin/env python
# Copyright 2014-2021 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Tianyu Zhu <zhutianyu1991@gmail.com>
#
"""
Spin-restricted MP2 quasiparticle energies
"""
import time
import numpy
import numpy as np
from scipy.optimize import newton
from pyscf import ao2mo, lib, scf
from pyscf.lib import logger
from pyscf.mp.mp2 import MP2, get_nocc, get_nmo, get_frozen_mask, _mo_energy_without_core, _mo_without_core, _mem_usage
from pyscf import __config__
from fcdmft.rpa.mp2.mol.qsmp2 import _make_eris
einsum = lib.einsum
[docs]
def kernel(mp2, mo_energy, mo_coeff, orbs):
"""MP2-corrected quasiparticle orbital energies
Parameters
----------
mp2 : MP2QP object
mo_energy : double array
molecular orbital energies
mo_coeff : double ndarray
molecular orbital coefficients
orbs : list
a list of orbital indices
Returns
-------
conv : bool
whether the quasiparticle equation is converged
mo_energy : double array
MP2-corrected quasiparticle orbital energies
"""
assert mp2.frozen is None
nmo = mp2.nmo
nocc = mp2.nocc
eta = mp2.eta
regularized = mp2.regularized
kappa = mp2.kappa
norm = mp2.norm
eia = mo_energy[:nocc, None] - mo_energy[None, nocc:]
mp2.eris = eris = mp2.ao2mo(mo_coeff)
conv = True
mf_mo_energy = np.array(mo_energy, copy=True)
mo_energy = np.zeros_like(mp2._scf.mo_energy)
mp2.z_factor = np.zeros_like(mp2._scf.mo_energy)
for p in orbs:
if mp2.linearized:
# linearized QP-MP2
ep = mf_mo_energy[p]
sigmaR = get_sigma_mp2(
nocc=nocc,
mo_energy=mf_mo_energy,
eris=eris,
p=p,
q=p,
freq=ep + 1j * eta,
eia=eia,
regularized=regularized,
kappa=kappa,
norm=norm,
).real
zn = get_mp2_sigma_derivative(
nocc=nocc,
mo_energy=mf_mo_energy,
eris=eris,
p=p,
q=p,
freq=ep + 1j * eta,
eia=eia,
regularized=regularized,
kappa=kappa,
norm=norm,
)
mp2.z_factor[p] = zn = 1.0 / (1.0 - zn)
mo_energy[p] = ep + zn * sigmaR
else:
# self-consistently solve QP equation
def quasiparticle(omega):
sigmaR = get_sigma_mp2(
nocc=nocc,
mo_energy=mf_mo_energy,
eris=eris,
p=p,
q=p,
freq=omega + 1j * eta,
eia=eia,
regularized=regularized,
kappa=kappa,
norm=norm,
).real
return omega - mf_mo_energy[p] - sigmaR
try:
mo_energy[p] = newton(quasiparticle, mf_mo_energy[p], tol=1e-5, maxiter=100)
except RuntimeError:
conv = False
mp2.mo_energy = mo_energy
with np.printoptions(threshold=nmo):
logger.debug(mp2, ' MP2QP mo_energy =\n%s', mo_energy)
return conv, mo_energy
[docs]
def get_sigma_mp2(nocc, mo_energy, eris, p, q, freq, eia, regularized=False, kappa=1.0, norm=2):
"""Compute the element of the MP2 self-energy in MO basis.
Parameters
----------
nocc : int
number of occupied orbitals.
mo_energy : double array
molecular orbital energies
eris : double ndarray
ERI in MO basis, nmo**4
p : int
row index of the self-energy matrix
q : int
column index of the self-energy matrix
freq : complex
frequency point
eia : double ndarray
occupied and virtual orbital energy difference
regularized : bool, optional
apply regularization, by default False
kappa : float, optional
kappa factor, by default 1.0
norm : int, optional
regularization norm, by default 2
Returns
-------
sigma : double
self-energy matrix element
"""
iapb = eris.ovnv[:, :, p, :]
iaqb = eris.ovnv[:, :, q, :]
iajp = eris.ovon[:, :, :, p]
iajq = eris.ovon[:, :, :, q]
ewa = freq - mo_energy[nocc:]
eiw = mo_energy[:nocc] - freq
e_iab = lib.direct_sum('ia+b->iab', eia, ewa)
e_iaj = -lib.direct_sum('ia+j->iaj', eia, eiw)
t_iab = iapb / e_iab
t_iaj = iajp / e_iaj
# apply regularization
if regularized:
t_iab = t_iab * (1 - np.exp(-np.abs(e_iab) * kappa)) ** norm
t_iaj = t_iaj * (1 - np.exp(-np.abs(e_iaj) * kappa)) ** norm
sigma = 2 * einsum('iab,iab', t_iab, iaqb.conj())
sigma -= einsum('iab,iba', t_iab, iaqb.conj())
sigma += 2 * einsum('iaj,iaj', t_iaj, iajq.conj())
sigma -= einsum('iaj,jai', t_iaj, iajq.conj())
return sigma
[docs]
def get_mp2_sigma_derivative(nocc, mo_energy, eris, p, q, freq, eia, regularized=False, kappa=1.0, norm=2):
"""Compute the element of the MP2 self-energy derivative in MO basis.
Parameters
----------
nocc : int
number of occupied orbitals.
mo_energy : double array
molecular orbital energies
eris : double ndarray
ERI in MO basis, nmo**4
p : int
row index of the self-energy matrix
q : int
column index of the self-energy matrix
freq : complex
frequency point
eia : double ndarray
occupied and virtual orbital energy difference
regularized : bool, optional
apply regularization, by default False
kappa : float, optional
kappa factor, by default 1.0
norm : int, optional
regularization norm, by default 2
Returns
-------
sigma : double
self-energy derivative matrix element
"""
iapb = eris.ovnv[:, :, p, :]
iaqb = eris.ovnv[:, :, q, :]
iajp = eris.ovon[:, :, :, p]
iajq = eris.ovon[:, :, :, q]
ewa = freq - mo_energy[nocc:]
eiw = mo_energy[:nocc] - freq
e_iab = lib.direct_sum('ia+b->iab', eia, ewa)
e_iaj = -lib.direct_sum('ia+j->iaj', eia, eiw)
t_iab = iapb / (e_iab**2)
t_iaj = iajp / (e_iaj**2)
# apply regularization
if regularized:
t_iab = t_iab * (1 - np.exp(-np.abs(e_iab) * kappa)) ** norm
t_iaj = t_iaj * (1 - np.exp(-np.abs(e_iaj) * kappa)) ** norm
derivative = -2 * einsum('iab,iab', t_iab, iaqb.conj())
derivative -= -einsum('iab,iba', t_iab, iaqb.conj())
derivative += -2 * einsum('iaj,iaj', t_iaj, iajq.conj())
derivative -= -einsum('iaj,jai', t_iaj, iajq.conj())
derivative = derivative.real
return derivative
[docs]
class MP2QP(MP2):
def __init__(self, mf, frozen=None, regularized=False, kappa=1.5, norm=2):
MP2.__init__(self, mf, frozen=frozen)
# options
self.eta = 5e-3
self.linearized = False # use linearized quasiparticle equation
self.regularized = regularized # use regularized MP2
self.kappa = kappa # regularization factor
self.norm = norm # regularization norm
# matrices
self.mo_energy = None # orbital energy
self.sigma = None # self-energy matrix
self.z_factor = None # Z-shot factor or quasiparticle renormalization factor
self.eris = None # two-electron integrals used in MP2
[docs]
def dump_flags(self):
log = logger.Logger(self.stdout, self.verbose)
log.info('')
log.info('******** %s ********', self.__class__)
log.info('method = %s', self.__class__.__name__)
nocc = self.nocc
nvir = self.nmo - nocc
log.info('MP2QP nocc = %d, nvir = %d', nocc, nvir)
log.info('use perturbative linearized QP eqn = %s', self.linearized)
if self.regularized is True:
log.info('regularized MP2 = %s', self.regularized)
log.info('kappa factor = %.2f', self.kappa)
log.info('regularized norm = %.2f', self.norm)
log.info('max_memory %d MB (current use %d MB)', self.max_memory, lib.current_memory()[0])
return self
@property
def nocc(self):
return self.get_nocc()
@nocc.setter
def nocc(self, n):
self._nocc = n
@property
def nmo(self):
return self.get_nmo()
@nmo.setter
def nmo(self, n):
self._nmo = n
get_nocc = get_nocc
get_nmo = get_nmo
get_frozen_mask = get_frozen_mask
[docs]
def kernel(self, mo_energy=None, mo_coeff=None, orbs=None):
"""
Args:
mo_energy : 1D array (nmo), mean-field mo energy
mo_coeff : 2D array (nmo, nmo), mean-field mo coefficient
orbs: list, orbital indices
Returns:
self.mo_energy : MP2 quasiparticle energy
"""
if mo_coeff is None:
mo_coeff = _mo_without_core(self, self._scf.mo_coeff)
if mo_energy is None:
mo_energy = _mo_energy_without_core(self, self._scf.mo_energy)
if orbs is None:
orbs = range(self.nmo)
cput0 = (time.process_time(), time.perf_counter())
self.dump_flags()
self.converged, self.mo_energy = kernel(self, mo_energy, mo_coeff, orbs=orbs)
logger.timer(self, 'MP2QP', *cput0)
return self.mo_energy
[docs]
def ao2mo(self, mo_coeff=None):
return _make_eris(self, mo_coeff, verbose=self.verbose)
[docs]
def get_graphical_solution(self, omega, orbs=None):
"""Get graphical solution of MP2.
Plot (omega, mo_energy_omega) and (omega, sigma), the intersection is the graphical solution.
Args:
omega (double array): frequency.
orbs (list, optional): index of desired orbitals. Defaults to None.
Returns:
mo_energy_omega (double ndarray): frequency - orbital energy.
sigma (double ndarray): dynamical diagonal self-energy.
"""
if orbs is None:
orbs = range(self.nmo)
sigma = np.zeros(shape=[len(orbs), len(omega)], dtype=np.double)
eia = self._scf.mo_energy[: self.nocc, None] - self._scf.mo_energy[None, self.nocc :]
if self.eris is None:
self.eris = self.ao2mo(self._scf.mo_coeff)
for p in orbs:
for w in range(len(omega)):
sigma[p, w] = get_sigma_mp2(
nocc=self.nocc,
mo_energy=self._scf.mo_energy,
eris=self.eris,
p=p,
q=p,
freq=omega[w] + 1j * self.eta,
eia=eia,
regularized=self.regularized,
kappa=self.kappa,
norm=self.norm,
).real
mo_energy_omega = lib.direct_sum('p+w->pw', -self._scf.mo_energy[orbs], omega)
return mo_energy_omega, sigma
if __name__ == '__main__':
from pyscf import gto, scf
mol = gto.Mole()
mol.verbose = 5
mol.atom = [[8, (0.0, 0.0, 0.0)], [1, (0.0, -0.7571, 0.5861)], [1, (0.0, 0.7571, 0.5861)]]
mol.basis = 'def2-tzvpp'
mol.build()
mf = scf.RHF(mol)
mf.kernel()
nocc = mol.nelectron // 2
mp2 = MP2QP(mf)
mp2.linearized = True
mp2.kernel(orbs=range(nocc - 1, nocc + 1))
mp2 = MP2QP(mf)
mp2.linearized = True
mp2.regularized = True
mp2.kappa = 1.5
mp2.kernel(orbs=range(nocc - 1, nocc + 1))