#!/usr/bin/env python
# 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: Qiming Sun <osirpt.sun@gmail.com>
#
'''
Generalized Kohn-Sham
'''
import numpy
import scipy.linalg
from pyscf import lib
from pyscf.lib import logger
from pyscf.scf import ghf
from pyscf.dft import rks
from pyscf.dft.numint2c import NumInt2C
[docs]
def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
'''Coulomb + XC functional
.. note::
This function will change the ks object.
Args:
ks : an instance of :class:`RKS`
XC functional are controlled by ks.xc attribute. Attribute
ks.grids might be initialized.
dm : ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
dm_last : ndarray or a list of ndarrays or 0
The density matrix baseline. If not 0, this function computes the
increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last : ndarray or a list of ndarrays or 0
The reference Vxc potential matrix.
hermi : int
Whether J, K matrix is hermitian
| 0 : no hermitian or symmetric
| 1 : hermitian
| 2 : anti-hermitian
Returns:
matrix Veff = J + Vxc. Veff can be a list matrices, if the input
dm is a list of density matrices.
'''
if mol is None: mol = ks.mol
if dm is None: dm = ks.make_rdm1()
ks.initialize_grids(mol, dm)
t0 = (logger.process_clock(), logger.perf_counter())
ground_state = isinstance(dm, numpy.ndarray) and dm.ndim == 2
if hermi == 2: # because rho = 0
n, exc, vxc = 0, 0, 0
else:
max_memory = ks.max_memory - lib.current_memory()[0]
ni = ks._numint
n, exc, vxc = ni.get_vxc(mol, ks.grids, ks.xc, dm,
hermi=hermi, max_memory=max_memory)
logger.debug(ks, 'nelec by numeric integration = %s', n)
if ks.nlc or ni.libxc.is_nlc(ks.xc):
if ni.libxc.is_nlc(ks.xc):
xc = ks.xc
else:
assert ni.libxc.is_nlc(ks.nlc)
xc = ks.nlc
n, enlc, vnlc = ni.nr_nlc_vxc(mol, ks.nlcgrids, xc, dm,
hermi=hermi, max_memory=max_memory)
exc += enlc
vxc += vnlc
logger.debug(ks, 'nelec with nlc grids = %s', n)
t0 = logger.timer(ks, 'vxc', *t0)
if not ni.libxc.is_hybrid_xc(ks.xc):
vk = None
if (ks._eri is None and ks.direct_scf and
getattr(vhf_last, 'vj', None) is not None):
ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
vj = ks.get_j(mol, ddm, hermi)
vj += vhf_last.vj
else:
vj = ks.get_j(mol, dm, hermi)
vxc += vj
else:
omega, alpha, hyb = ni.rsh_and_hybrid_coeff(ks.xc, spin=mol.spin)
if (ks._eri is None and ks.direct_scf and
getattr(vhf_last, 'vk', None) is not None):
ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
vj, vk = ks.get_jk(mol, ddm, hermi)
vk *= hyb
if omega != 0:
vklr = ks.get_k(mol, ddm, hermi, omega=omega)
vklr *= (alpha - hyb)
vk += vklr
vj += vhf_last.vj
vk += vhf_last.vk
else:
vj, vk = ks.get_jk(mol, dm, hermi)
vk *= hyb
if omega != 0:
vklr = ks.get_k(mol, dm, hermi, omega=omega)
vklr *= (alpha - hyb)
vk += vklr
vxc += vj - vk
if ground_state:
exc -= numpy.einsum('ij,ji', dm, vk).real * .5
if ground_state:
ecoul = numpy.einsum('ij,ji', dm, vj).real * .5
else:
ecoul = None
vxc = lib.tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk)
return vxc
energy_elec = rks.energy_elec
[docs]
class GKS(rks.KohnShamDFT, ghf.GHF):
'''Generalized Kohn-Sham'''
get_veff = get_veff
energy_elec = energy_elec
def __init__(self, mol, xc='LDA,VWN'):
ghf.GHF.__init__(self, mol)
rks.KohnShamDFT.__init__(self, xc)
self._numint = NumInt2C()
[docs]
def dump_flags(self, verbose=None):
ghf.GHF.dump_flags(self, verbose)
rks.KohnShamDFT.dump_flags(self, verbose)
logger.info(self, 'collinear = %s', self._numint.collinear)
if self._numint.collinear[0] == 'm':
logger.info(self, 'mcfun spin_samples = %s', self._numint.spin_samples)
logger.info(self, 'mcfun collinear_thrd = %s', self._numint.collinear_thrd)
logger.info(self, 'mcfun collinear_samples = %s', self._numint.collinear_samples)
return self
@property
def collinear(self):
return self._numint.collinear
@collinear.setter
def collinear(self, val):
self._numint.collinear = val
@property
def spin_samples(self):
return self._numint.spin_samples
@spin_samples.setter
def spin_samples(self, val):
self._numint.spin_samples = val
[docs]
def nuc_grad_method(self):
raise NotImplementedError
[docs]
def to_hf(self):
'''Convert to GHF object.'''
return self._transfer_attrs_(self.mol.GHF())
if __name__ == '__main__':
from pyscf import gto
mol = gto.Mole()
mol.verbose = 3
mol.atom = 'H 0 0 0; H 0 0 1; O .5 .6 .2'
mol.basis = 'ccpvdz'
mol.build()
mf = GKS(mol)
mf.xc = 'b3lyp'
mf.kernel()
dm = mf.init_guess_by_1e(mol)
dm = dm + 0j
nao = mol.nao_nr()
numpy.random.seed(12)
dm[:nao,nao:] = numpy.random.random((nao,nao)) * .1j
dm[nao:,:nao] = dm[:nao,nao:].T.conj()
mf.kernel(dm)
mf.canonicalize(mf.mo_coeff, mf.mo_occ)
mf.analyze()
print(mf.spin_square())
print(mf.e_tot - -76.2760115704274)