Source code for pyscf.adc.radc

# Copyright 2014-2022 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: Samragni Banerjee <samragnibanerjee4@gmail.com>
#         Alexander Sokolov <alexander.y.sokolov@gmail.com>
#

'''
Restricted algebraic diagrammatic construction
'''
import numpy as np
import pyscf.ao2mo as ao2mo
from pyscf import lib
from pyscf.lib import logger
from pyscf.adc import radc_ao2mo
from pyscf.adc import dfadc
from pyscf.adc import radc_amplitudes
from pyscf import __config__
from pyscf import df
from pyscf import symm


# Excited-state kernel
[docs] def kernel(adc, nroots=1, guess=None, eris=None, verbose=None): adc.method = adc.method.lower() if adc.method not in ("adc(2)", "adc(2)-x", "adc(3)"): raise NotImplementedError(adc.method) cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(adc.stdout, adc.verbose) if adc.verbose >= logger.WARN: adc.check_sanity() adc.dump_flags() if eris is None: eris = adc.transform_integrals() imds = adc.get_imds(eris) matvec, diag = adc.gen_matvec(imds, eris) guess = adc.get_init_guess(nroots, diag, ascending=True) conv, adc.E, U = lib.linalg_helper.davidson_nosym1( lambda xs : [matvec(x) for x in xs], guess, diag, nroots=nroots, verbose=log, tol=adc.conv_tol, max_cycle=adc.max_cycle, max_space=adc.max_space, tol_residual=adc.tol_residual) adc.U = np.array(U).T.copy() if adc.compute_properties: adc.P,adc.X = adc.get_properties(nroots) header = ("\n*************************************************************" "\n ADC calculation summary" "\n*************************************************************") logger.info(adc, header) for n in range(nroots): print_string = ('%s root %d | Energy (Eh) = %14.10f | Energy (eV) = %12.8f ' % (adc.method, n, adc.E[n], adc.E[n]*27.2114)) if adc.compute_properties: print_string += ("| Spec factors = %10.8f " % adc.P[n]) print_string += ("| conv = %s" % conv[n]) logger.info(adc, print_string) log.timer('ADC', *cput0) return adc.E, adc.U, adc.P, adc.X
[docs] class RADC(lib.StreamObject): '''Ground state calculations Attributes: verbose : int Print level. Default value equals to :class:`Mole.verbose` max_memory : float or int Allowed memory in MB. Default value equals to :class:`Mole.max_memory` incore_complete : bool Avoid all I/O. Default is False. method : string nth-order ADC method. Options are : ADC(2), ADC(2)-X, ADC(3). Default is ADC(2). >>> mol = gto.M(atom = 'H 0 0 0; F 0 0 1.1', basis = 'ccpvdz') >>> mf = scf.RHF(mol).run() >>> myadc = adc.RADC(mf).run() Saved results e_corr : float MPn correlation correction e_tot : float Total energy (HF + correlation) t1, t2 : T amplitudes t1[i,a], t2[i,j,a,b] (i,j in occ, a,b in virt) ''' incore_complete = getattr(__config__, 'adc_radc_RADC_incore_complete', False) async_io = getattr(__config__, 'adc_radc_RADC_async_io', True) blkmin = getattr(__config__, 'adc_radc_RADC_blkmin', 4) memorymin = getattr(__config__, 'adc_radc_RADC_memorymin', 2000) _keys = set(( 'tol_residual','conv_tol', 'e_corr', 'method', 'mo_coeff', 'mol', 'mo_energy', 'incore_complete', 'scf_energy', 'e_tot', 't1', 'frozen', 'chkfile', 'max_space', 't2', 'mo_occ', 'max_cycle', 'imds', 'method', 'method_type', 'with_df', 'compute_properties', 'approx_trans_moments', 'evec_print_tol', 'spec_factor_print_tol', 'ncvs', 'E', 'U', 'P', 'X', )) def __init__(self, mf, frozen=0, mo_coeff=None, mo_occ=None): from pyscf import gto if 'dft' in str(mf.__module__): raise NotImplementedError('DFT reference for UADC') if mo_coeff is None: mo_coeff = mf.mo_coeff if mo_occ is None: mo_occ = mf.mo_occ self.mol = mf.mol self._scf = mf self.verbose = self.mol.verbose self.stdout = self.mol.stdout self.max_memory = mf.max_memory self.max_space = getattr(__config__, 'adc_radc_RADC_max_space', 12) self.max_cycle = getattr(__config__, 'adc_radc_RADC_max_cycle', 50) self.conv_tol = getattr(__config__, 'adc_radc_RADC_conv_tol', 1e-12) self.tol_residual = getattr(__config__, 'adc_radc_RADC_tol_res', 1e-6) self.scf_energy = mf.e_tot self.frozen = frozen self.incore_complete = self.incore_complete or self.mol.incore_anyway self.mo_coeff = mo_coeff self.mo_occ = mo_occ self.e_corr = None self.t1 = None self.t2 = None self.imds = lambda:None self._nocc = mf.mol.nelectron//2 self._nmo = mo_coeff.shape[1] self._nvir = self._nmo - self._nocc self.mo_energy = mf.mo_energy self.chkfile = mf.chkfile self.method = "adc(2)" self.method_type = "ip" self.with_df = None self.compute_properties = True self.approx_trans_moments = False self.evec_print_tol = 0.1 self.spec_factor_print_tol = 0.1 self.ncvs = None self.E = None self.U = None self.P = None self.X = None compute_amplitudes = radc_amplitudes.compute_amplitudes compute_energy = radc_amplitudes.compute_energy transform_integrals = radc_ao2mo.transform_integrals_incore
[docs] def dump_flags(self, verbose=None): logger.info(self, '') logger.info(self, '******** %s ********', self.__class__) logger.info(self, 'max_space = %d', self.max_space) logger.info(self, 'max_cycle = %d', self.max_cycle) logger.info(self, 'conv_tol = %s', self.conv_tol) logger.info(self, 'max_memory %d MB (current use %d MB)', self.max_memory, lib.current_memory()[0]) return self
[docs] def dump_flags_gs(self, verbose=None): logger.info(self, '') logger.info(self, '******** %s ********', self.__class__) logger.info(self, 'max_memory %d MB (current use %d MB)', self.max_memory, lib.current_memory()[0]) return self
[docs] def kernel_gs(self): assert (self.mo_coeff is not None) assert (self.mo_occ is not None) self.method = self.method.lower() if self.method not in ("adc(2)", "adc(2)-x", "adc(3)"): raise NotImplementedError(self.method) if self.verbose >= logger.WARN: self.check_sanity() self.dump_flags_gs() nmo = self._nmo nao = self.mo_coeff.shape[0] nmo_pair = nmo * (nmo+1) // 2 nao_pair = nao * (nao+1) // 2 mem_incore = (max(nao_pair**2, nmo**4) + nmo_pair**2) * 8/1e6 mem_now = lib.current_memory()[0] if getattr(self, 'with_df', None) or getattr(self._scf, 'with_df', None): if getattr(self, 'with_df', None): self.with_df = self.with_df else: self.with_df = self._scf.with_df def df_transform(): return radc_ao2mo.transform_integrals_df(self) self.transform_integrals = df_transform elif (self._scf._eri is None or (mem_incore+mem_now >= self.max_memory and not self.incore_complete)): def outcore_transform(): return radc_ao2mo.transform_integrals_outcore(self) self.transform_integrals = outcore_transform eris = self.transform_integrals() self.e_corr, self.t1, self.t2 = radc_amplitudes.compute_amplitudes_energy( self, eris=eris, verbose=self.verbose) self._finalize() return self.e_corr, self.t1, self.t2
[docs] def kernel(self, nroots=1, guess=None, eris=None): assert (self.mo_coeff is not None) assert (self.mo_occ is not None) self.method = self.method.lower() if self.method not in ("adc(2)", "adc(2)-x", "adc(3)"): raise NotImplementedError(self.method) if self.verbose >= logger.WARN: self.check_sanity() self.dump_flags_gs() nmo = self._nmo nao = self.mo_coeff.shape[0] nmo_pair = nmo * (nmo+1) // 2 nao_pair = nao * (nao+1) // 2 mem_incore = (max(nao_pair**2, nmo**4) + nmo_pair**2) * 8/1e6 mem_now = lib.current_memory()[0] if getattr(self, 'with_df', None) or getattr(self._scf, 'with_df', None): if getattr(self, 'with_df', None): self.with_df = self.with_df else: self.with_df = self._scf.with_df def df_transform(): return radc_ao2mo.transform_integrals_df(self) self.transform_integrals = df_transform elif (self._scf._eri is None or (mem_incore+mem_now >= self.max_memory and not self.incore_complete)): def outcore_transform(): return radc_ao2mo.transform_integrals_outcore(self) self.transform_integrals = outcore_transform eris = self.transform_integrals() self.e_corr, self.t1, self.t2 = radc_amplitudes.compute_amplitudes_energy( self, eris=eris, verbose=self.verbose) self._finalize() self.method_type = self.method_type.lower() if (self.method_type == "ea"): e_exc, v_exc, spec_fac, x, adc_es = self.ea_adc(nroots=nroots, guess=guess, eris=eris) elif(self.method_type == "ip"): if not isinstance(self.ncvs, type(None)) and self.ncvs > 0: e_exc, v_exc, spec_fac, x, adc_es = self.ip_cvs_adc( nroots=nroots, guess=guess, eris=eris) else: e_exc, v_exc, spec_fac, x, adc_es = self.ip_adc( nroots=nroots, guess=guess, eris=eris) else: raise NotImplementedError(self.method_type) self._adc_es = adc_es return e_exc, v_exc, spec_fac, x
def _finalize(self): '''Hook for dumping results and clearing up the object.''' logger.note(self, 'E_corr = %.8f', self.e_corr) return self
[docs] def ea_adc(self, nroots=1, guess=None, eris=None): from pyscf.adc import radc_ea adc_es = radc_ea.RADCEA(self) e_exc, v_exc, spec_fac, x = adc_es.kernel(nroots, guess, eris) return e_exc, v_exc, spec_fac, x, adc_es
[docs] def ip_adc(self, nroots=1, guess=None, eris=None): from pyscf.adc import radc_ip adc_es = radc_ip.RADCIP(self) e_exc, v_exc, spec_fac, x = adc_es.kernel(nroots, guess, eris) return e_exc, v_exc, spec_fac, x, adc_es
[docs] def ip_cvs_adc(self, nroots=1, guess=None, eris=None): from pyscf.adc import radc_ip_cvs adc_es = radc_ip_cvs.RADCIPCVS(self) e_exc, v_exc, spec_fac, x = adc_es.kernel(nroots, guess, eris) return e_exc, v_exc, spec_fac, x, adc_es
[docs] def density_fit(self, auxbasis=None, with_df=None): if with_df is None: self.with_df = df.DF(self._scf.mol) self.with_df.max_memory = self.max_memory self.with_df.stdout = self.stdout self.with_df.verbose = self.verbose if auxbasis is None: self.with_df.auxbasis = self._scf.with_df.auxbasis else: self.with_df.auxbasis = auxbasis else: self.with_df = with_df return self
[docs] def analyze(self): self._adc_es.analyze()
[docs] def compute_dyson_mo(self): return self._adc_es.compute_dyson_mo()
[docs] def make_rdm1(self): return self._adc_es.make_rdm1()
if __name__ == '__main__': from pyscf import scf from pyscf import gto from pyscf import adc r = 1.098 mol = gto.Mole() mol.atom = [ ['N', (0., 0. , -r/2 )], ['N', (0., 0. , r/2)],] mol.basis = {'N':'aug-cc-pvdz'} mol.verbose = 0 mol.build() mf = scf.RHF(mol) mf.conv_tol = 1e-12 mf.kernel() myadc = adc.ADC(mf) ecorr, t_amp1, t_amp2 = myadc.kernel_gs() print(ecorr - -0.32201692360512324) myadcip = adc.radc_ip.RADCIP(myadc) e,v,p = kernel(myadcip,nroots=3) print("ADC(2) IP energies") print (e[0] - 0.5434389910483670) print (e[1] - 0.6240296243595950) print (e[2] - 0.6240296243595956) print("ADC(2) IP spectroscopic factors") print (p[0] - 1.7688097076459075) print (p[1] - 1.8192921131700284) print (p[2] - 1.8192921131700293) myadcea = adc.radc_ea.RADCEA(myadc) e,v,p = kernel(myadcea,nroots=3) print("ADC(2) EA energies") print (e[0] - 0.0961781923822576) print (e[1] - 0.1258326916409743) print (e[2] - 0.1380779405750178) print("ADC(2) EA spectroscopic factors") print (p[0] - 1.9832854445007961) print (p[1] - 1.9634368668786559) print (p[2] - 1.9783719593912672) myadc = adc.ADC(mf) myadc.method = "adc(3)" ecorr, t_amp1, t_amp2 = myadc.kernel_gs() print(ecorr - -0.31694173142858517) myadcip = adc.radc_ip.RADCIP(myadc) e,v,p = kernel(myadcip,nroots=3) print("ADC(3) IP energies") print (e[0] - 0.5667526829981027) print (e[1] - 0.6099995170092525) print (e[2] - 0.6099995170092529) print("ADC(3) IP spectroscopic factors") print (p[0] - 1.8173191958988848) print (p[1] - 1.8429224413853840) print (p[2] - 1.8429224413853851) myadcea = adc.radc_ea.RADCEA(myadc) e,v,p = kernel(myadcea,nroots=3) print("ADC(3) EA energies") print (e[0] - 0.0936790850738445) print (e[1] - 0.0983654552141278) print (e[2] - 0.1295709313652367) print("ADC(3) EA spectroscopic factors") print (p[0] - 1.8324175318668088) print (p[1] - 1.9840991060607487) print (p[2] - 1.9638550014980212) myadc.method = "adc(2)-x" e,v,p = myadc.kernel(nroots=4) print("ADC(2)-x IP energies") print (e[0] - 0.5405255360673724) print (e[1] - 0.6208026698756577) print (e[2] - 0.6208026698756582) print (e[3] - 0.6465332771967947) myadc.method_type = "ea" e,v,p = myadc.kernel(nroots=4) print("ADC(2)-x EA energies") print (e[0] - 0.0953065329985665) print (e[1] - 0.1238833070823509) print (e[2] - 0.1365693811939308) print (e[3] - 0.1365693811939316)