#!/usr/bin/env python
# Copyright 2014-2020 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: Qiming Sun <osirpt.sun@gmail.com>
#
import ctypes
from functools import reduce
import numpy
from pyscf import lib
from pyscf.lib import logger
from pyscf.ao2mo import _ao2mo
from pyscf.mcscf.casci import CASCI
from pyscf import df
[docs]
def density_fit(casscf, auxbasis=None, with_df=None):
'''Generate DF-CASSCF for given CASSCF object. It is done by overwriting
three CASSCF member functions:
* casscf.ao2mo which generates MO integrals
* casscf.get_veff which generate JK from core density matrix
* casscf.get_jk which
Args:
casscf : an CASSCF object
Kwargs:
auxbasis : str or basis dict
Same format to the input attribute mol.basis. If auxbasis is
None, auxiliary basis based on AO basis (if possible) or
even-tempered Gaussian basis will be used.
Returns:
An CASSCF object with a modified J, K matrix constructor which uses density
fitting integrals to compute J and K
Examples:
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvdz', verbose=0)
>>> mf = scf.RHF(mol)
>>> mf.scf()
>>> mc = DFCASSCF(mf, 4, 4)
-100.05994191570504
'''
if with_df is None:
if (getattr(casscf._scf, 'with_df', None) and
(auxbasis is None or auxbasis == casscf._scf.with_df.auxbasis)):
with_df = casscf._scf.with_df
else:
with_df = df.DF(casscf.mol)
with_df.max_memory = casscf.max_memory
with_df.stdout = casscf.stdout
with_df.verbose = casscf.verbose
with_df.auxbasis = auxbasis
if isinstance(casscf, _DFCAS):
if casscf.with_df is None:
casscf.with_df = with_df
elif getattr(casscf.with_df, 'auxbasis', None) != auxbasis:
#logger.warn(casscf, 'DF might have been initialized twice.')
casscf = casscf.copy()
casscf.with_df = with_df
return casscf
if isinstance(casscf, CASCI):
cls = _DFCASCI
else:
cls = _DFCASSCF
return lib.set_class(cls(casscf, with_df), (cls, casscf.__class__))
# A tag to label the derived MCSCF class
class _DFCAS:
__name_mixin__ = 'DF'
_keys = set(['with_df'])
def __init__(self, mc, with_df):
self.__dict__.update(mc.__dict__)
#self.grad_update_dep = 0
self.with_df = with_df
def undo_df(self):
'''Remove the DFCASCI/DFCASSCF Mixin'''
obj = lib.view(self, lib.drop_class(self.__class__, _DFCAS))
del obj.with_df
return obj
def dump_flags(self, verbose=None):
super().dump_flags(verbose)
logger.info(self, 'DFCASCI/DFCASSCF: density fitting for JK matrix '
'and 2e integral transformation')
return self
def reset(self, mol=None):
self.with_df.reset(mol)
return super().reset(mol)
# Modify get_veff for JK matrix of core density because get_h1eff calls
# self.get_veff to generate core JK
def get_veff(self, mol=None, dm=None, hermi=1):
if dm is None:
mocore = self.mo_coeff[:,:self.ncore]
dm = numpy.dot(mocore, mocore.T) * 2
vj, vk = self.get_jk(mol, dm, hermi)
return vj - vk * .5
# only approximate jk for self.update_jk_in_ah
def get_jk(self, mol, dm, hermi=1, with_j=True, with_k=True, omega=None):
if self.with_df:
return self.with_df.get_jk(dm, hermi,
with_j=with_j, with_k=with_k, omega=omega)
else:
return super().get_jk(mol, dm, hermi,
with_j=with_j, with_k=with_k, omega=omega)
def _exact_paaa(self, mo, u, out=None):
if self.with_df:
nmo = mo.shape[1]
ncore = self.ncore
ncas = self.ncas
nocc = ncore + ncas
mo1 = numpy.dot(mo, u)
mo1_cas = mo1[:,ncore:nocc]
paaa = self.with_df.ao2mo([mo1, mo1_cas, mo1_cas, mo1_cas], compact=False)
return paaa.reshape(nmo,ncas,ncas,ncas)
else:
return super()._exact_paaa(mol, u, out)
def nuc_grad_method(self):
raise NotImplementedError
def _state_average_nuc_grad_method (self, state=None):
raise NotImplementedError
class _DFCASCI(_DFCAS):
def get_h2eff(self, mo_coeff=None):
if self.with_df:
ncore = self.ncore
nocc = ncore + self.ncas
if mo_coeff is None:
mo_coeff = self.mo_coeff[:,ncore:nocc]
elif mo_coeff.shape[1] != self.ncas:
mo_coeff = mo_coeff[:,ncore:nocc]
return self.with_df.ao2mo(mo_coeff)
else:
return super(_DFCAS, self).get_h2eff(mo_coeff)
class _DFCASSCF(_DFCAS):
get_h2eff = _DFCASCI.get_h2eff
def ao2mo(self, mo_coeff=None):
if self.with_df:
return _ERIS(self, mo_coeff, self.with_df)
else:
return super(_DFCAS, self).ao2mo(mo_coeff)
def nuc_grad_method(self):
from pyscf.df.grad import casscf
return casscf.Gradients (self)
def _state_average_nuc_grad_method (self, state=None):
from pyscf.df.grad import sacasscf
return sacasscf.Gradients (self, state=state)
[docs]
def approx_hessian(casscf, auxbasis=None, with_df=None):
'''Approximate the orbital hessian with density fitting integrals
Note this function has no effects if the input casscf object is DF-CASSCF.
It only modifies the orbital hessian of normal CASSCF object.
Args:
casscf : an CASSCF object
Kwargs:
auxbasis : str or basis dict
Same format to the input attribute mol.basis.
The default basis 'weigend+etb' means weigend-coulomb-fit basis
for light elements and even-tempered basis for heavy elements.
Returns:
A CASSCF object with approximated JK contraction for orbital hessian
Examples:
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvdz', verbose=0)
>>> mf = scf.RHF(mol)
>>> mf.scf()
>>> mc = mcscf.approx_hessian(mcscf.CASSCF(mf, 4, 4))
-100.06458716530391
'''
if isinstance(casscf, CASCI):
return casscf # because CASCI does not need orbital optimization
if getattr(casscf, 'with_df', None):
return casscf
if with_df is None:
if (getattr(casscf._scf, 'with_df', None) and
(auxbasis is None or auxbasis == casscf._scf.with_df.auxbasis)):
with_df = casscf._scf.with_df
else:
with_df = df.DF(casscf.mol)
with_df.max_memory = casscf.max_memory
with_df.stdout = casscf.stdout
with_df.verbose = casscf.verbose
if auxbasis is not None:
with_df.auxbasis = auxbasis
dfcas = _DFHessianCASSCF(casscf, with_df)
return lib.set_class(dfcas, (_DFHessianCASSCF, casscf.__class__))
class _DFHessianCASSCF:
__name_mixin__ = 'DFHessian'
_keys = set(['with_df'])
def __init__(self, mc, with_df):
self.__dict__.update(mc.__dict__)
#self.grad_update_dep = 0
self.with_df = with_df
def undo_df(self):
'''Remove the DFHessianCASSCF Mixin'''
obj = lib.view(self, lib.drop_class(self.__class__, _DFHessianCASSCF))
del obj.with_df
return obj
def dump_flags(self, verbose=None):
super().dump_flags(verbose)
logger.info(self, 'CASSCF: density fitting for orbital hessian')
def reset(self, mol=None):
self.with_df.reset(mol)
return super().reset(mol)
def ao2mo(self, mo_coeff):
# the exact integral transformation
eris = super().ao2mo(mo_coeff)
log = logger.Logger(self.stdout, self.verbose)
# Add the approximate diagonal term for orbital hessian
t1 = t0 = (logger.process_clock(), logger.perf_counter())
mo = numpy.asarray(mo_coeff, order='F')
nao, nmo = mo.shape
ncore = self.ncore
eris.j_pc = numpy.zeros((nmo,ncore))
k_cp = numpy.zeros((ncore,nmo))
fmmm = _ao2mo.libao2mo.AO2MOmmm_nr_s2_iltj
fdrv = _ao2mo.libao2mo.AO2MOnr_e2_drv
ftrans = _ao2mo.libao2mo.AO2MOtranse2_nr_s2
max_memory = self.max_memory - lib.current_memory()[0]
blksize = max(4, int(min(self.with_df.blockdim, max_memory*.3e6/8/nmo**2)))
bufs1 = numpy.empty((blksize,nmo,nmo))
for eri1 in self.with_df.loop(blksize):
naux = eri1.shape[0]
buf = bufs1[:naux]
fdrv(ftrans, fmmm,
buf.ctypes.data_as(ctypes.c_void_p),
eri1.ctypes.data_as(ctypes.c_void_p),
mo.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(naux), ctypes.c_int(nao),
(ctypes.c_int*4)(0, nmo, 0, nmo),
ctypes.c_void_p(0), ctypes.c_int(0))
bufd = numpy.einsum('kii->ki', buf)
eris.j_pc += numpy.einsum('ki,kj->ij', bufd, bufd[:,:ncore])
k_cp += numpy.einsum('kij,kij->ij', buf[:,:ncore], buf[:,:ncore])
t1 = log.timer_debug1('j_pc and k_pc', *t1)
eris.k_pc = k_cp.T.copy()
log.timer('ao2mo density fit part', *t0)
return eris
def get_jk(self, mol, dm, hermi=1, with_j=True, with_k=True, omega=None):
if self.with_df:
return self.with_df.get_jk(dm, hermi,
with_j=with_j, with_k=with_k, omega=omega)
else:
return super().get_jk(mol, dm, hermi,
with_j=with_j, with_k=with_k, omega=omega)
class _ERIS:
def __init__(self, casscf, mo, with_df):
log = logger.Logger(casscf.stdout, casscf.verbose)
mol = casscf.mol
nao, nmo = mo.shape
ncore = casscf.ncore
ncas = casscf.ncas
nocc = ncore + ncas
naoaux = with_df.get_naoaux()
mem_incore, mem_outcore, mem_basic = _mem_usage(ncore, ncas, nmo)
mem_now = lib.current_memory()[0]
max_memory = max(3000, casscf.max_memory*.9-mem_now)
if max_memory < mem_basic:
log.warn('Calculation needs %d MB memory, over CASSCF.max_memory (%d MB) limit',
(mem_basic+mem_now)/.9, casscf.max_memory)
t1 = t0 = (logger.process_clock(), logger.perf_counter())
self.feri = lib.H5TmpFile()
self.ppaa = self.feri.create_dataset('ppaa', (nmo,nmo,ncas,ncas), 'f8')
self.papa = self.feri.create_dataset('papa', (nmo,ncas,nmo,ncas), 'f8')
self.j_pc = numpy.zeros((nmo,ncore))
k_cp = numpy.zeros((ncore,nmo))
mo = numpy.asarray(mo, order='F')
fxpp = lib.H5TmpFile()
blksize = max(4, int(min(with_df.blockdim, (max_memory*.95e6/8-naoaux*nmo*ncas)/3/nmo**2)))
bufpa = numpy.empty((naoaux,nmo,ncas))
bufs1 = numpy.empty((blksize,nmo,nmo))
fmmm = _ao2mo.libao2mo.AO2MOmmm_nr_s2_iltj
fdrv = _ao2mo.libao2mo.AO2MOnr_e2_drv
ftrans = _ao2mo.libao2mo.AO2MOtranse2_nr_s2
fxpp_keys = []
b0 = 0
for k, eri1 in enumerate(with_df.loop(blksize)):
naux = eri1.shape[0]
bufpp = bufs1[:naux]
fdrv(ftrans, fmmm,
bufpp.ctypes.data_as(ctypes.c_void_p),
eri1.ctypes.data_as(ctypes.c_void_p),
mo.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(naux), ctypes.c_int(nao),
(ctypes.c_int*4)(0, nmo, 0, nmo),
ctypes.c_void_p(0), ctypes.c_int(0))
fxpp_keys.append([str(k), b0, b0+naux])
fxpp[str(k)] = bufpp.transpose(1,2,0)
bufpa[b0:b0+naux] = bufpp[:,:,ncore:nocc]
bufd = numpy.einsum('kii->ki', bufpp)
self.j_pc += numpy.einsum('ki,kj->ij', bufd, bufd[:,:ncore])
k_cp += numpy.einsum('kij,kij->ij', bufpp[:,:ncore], bufpp[:,:ncore])
b0 += naux
t1 = log.timer_debug1('j_pc and k_pc', *t1)
self.k_pc = k_cp.T.copy()
bufs1 = bufpp = None
t1 = log.timer('density fitting ao2mo pass1', *t0)
mem_now = lib.current_memory()[0]
nblk = int(max(8, min(nmo, ((max_memory-mem_now)*1e6/8-bufpa.size)/(ncas**2*nmo))))
bufs1 = numpy.empty((nblk,ncas,nmo,ncas))
dgemm = lib.numpy_helper._dgemm
for p0, p1 in prange(0, nmo, nblk):
#tmp = numpy.dot(bufpa[:,p0:p1].reshape(naoaux,-1).T,
# bufpa.reshape(naoaux,-1))
tmp = bufs1[:p1-p0]
dgemm('T', 'N', (p1-p0)*ncas, nmo*ncas, naoaux,
bufpa.reshape(naoaux,-1), bufpa.reshape(naoaux,-1),
tmp.reshape(-1,nmo*ncas), 1, 0, p0*ncas, 0, 0)
self.papa[p0:p1] = tmp.reshape(p1-p0,ncas,nmo,ncas)
bufaa = bufpa[:,ncore:nocc,:].copy().reshape(-1,ncas**2)
bufs1 = bufpa = None
t1 = log.timer('density fitting papa pass2', *t1)
mem_now = lib.current_memory()[0]
nblk = int(max(8, min(nmo, (max_memory-mem_now)*1e6/8/(nmo*naoaux+ncas**2*nmo))))
bufs1 = numpy.empty((nblk,nmo,naoaux))
bufs2 = numpy.empty((nblk,nmo,ncas,ncas))
for p0, p1 in prange(0, nmo, nblk):
nrow = p1 - p0
buf = bufs1[:nrow]
tmp = bufs2[:nrow].reshape(-1,ncas**2)
for key, col0, col1 in fxpp_keys:
buf[:nrow,:,col0:col1] = fxpp[key][p0:p1]
lib.dot(buf.reshape(-1,naoaux), bufaa, 1, tmp)
self.ppaa[p0:p1] = tmp.reshape(p1-p0,nmo,ncas,ncas)
bufs1 = bufs2 = buf = None
t1 = log.timer('density fitting ppaa pass2', *t1)
self.feri.flush()
dm_core = numpy.dot(mo[:,:ncore], mo[:,:ncore].T)
vj, vk = casscf.get_jk(mol, dm_core)
self.vhf_c = reduce(numpy.dot, (mo.T, vj*2-vk, mo))
t0 = log.timer('density fitting ao2mo', *t0)
def _mem_usage(ncore, ncas, nmo):
outcore = basic = ncas**2*nmo**2*2 * 8/1e6
incore = outcore + (ncore+ncas)*nmo**3*4/1e6
return incore, outcore, basic
[docs]
def prange(start, end, step):
for i in range(start, end, step):
yield i, min(i+step, end)
if __name__ == '__main__':
from pyscf import gto
from pyscf import scf
from pyscf import mcscf
from pyscf.mcscf import addons
mol = gto.Mole()
mol.atom = [
['O', ( 0., 0. , 0. )],
['H', ( 0., -0.757, 0.587)],
['H', ( 0., 0.757 , 0.587)],]
mol.basis = {'H': 'cc-pvdz',
'O': 'cc-pvdz',}
mol.build()
m = scf.RHF(mol)
ehf = m.scf()
mc = approx_hessian(mcscf.CASSCF(m, 6, 4))
mc.verbose = 4
mo = addons.sort_mo(mc, m.mo_coeff, (3,4,6,7,8,9), 1)
emc = mc.kernel(mo)[0]
print(ehf, emc, emc-ehf)
#-76.0267656731 -76.0873922924 -0.0606266193028
print(emc - -76.0873923174, emc - -76.0926176464)
mc = approx_hessian(mcscf.CASSCF(m, 6, (3,1)))
mc.verbose = 4
emc = mc.mc2step(mo)[0]
print(emc - -75.7155632535814)
mf = scf.density_fit(m)
mf.kernel()
#mc = density_fit(mcscf.CASSCF(mf, 6, 4))
#mc = mcscf.CASSCF(mf, 6, 4)
mc = mcscf.DFCASSCF(mf, 6, 4)
mc.verbose = 4
mo = addons.sort_mo(mc, mc.mo_coeff, (3,4,6,7,8,9), 1)
emc = mc.kernel(mo)[0]
print(emc, 'ref = -76.0917567904955', emc - -76.0917567904955)
mc.with_dep4 = True
mc.max_cycle_micro = 10
emc = mc.kernel(mo)[0]
print(emc, 'ref = -76.0917567904955', emc - -76.0917567904955)
#mc = density_fit(mcscf.CASCI(mf, 6, 4))
#mc = mcscf.CASCI(mf, 6, 4)
mc = mcscf.DFCASCI(mf, 6, 4)
mo = addons.sort_mo(mc, mc.mo_coeff, (3,4,6,7,8,9), 1)
emc = mc.kernel(mo)[0]
print(emc, 'ref = -76.0476686258461', emc - -76.0476686258461)