Skip to content

Commit

Permalink
VV10 gradients (pyscf#1486)
Browse files Browse the repository at this point in the history
* VV10 derivatives without grids response

* implement grids response for VV10

* UKS grad for VV10 + some code cleanup
  • Loading branch information
kylebystrom authored Nov 14, 2022

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
1 parent 7862701 commit 654cc36
Showing 5 changed files with 277 additions and 9 deletions.
184 changes: 178 additions & 6 deletions pyscf/grad/rks.py
Original file line number Diff line number Diff line change
@@ -24,8 +24,11 @@
from pyscf import lib
from pyscf.lib import logger
from pyscf.grad import rhf as rhf_grad
from pyscf.dft import numint, radi, gen_grid
from pyscf.dft import numint, radi, gen_grid, xc_deriv
from pyscf import __config__
import ctypes

libdft = lib.load_library('libdft')


def get_veff(ks_grad, mol=None, dm=None):
@@ -45,11 +48,16 @@ def get_veff(ks_grad, mol=None, dm=None):
grids = ks_grad.grids
else:
grids = mf.grids
if mf.nlc != '':
if ks_grad.nlcgrids is not None:
nlcgrids = ks_grad.nlcgrids
else:
nlcgrids = mf.nlcgrids
if nlcgrids.coords is None:
nlcgrids.build(with_non0tab=True)
if grids.coords is None:
grids.build(with_non0tab=True)

if mf.nlc != '':
raise NotImplementedError
#enabling range-separated hybrids
omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, spin=mol.spin)

@@ -59,10 +67,24 @@ def get_veff(ks_grad, mol=None, dm=None):
exc, vxc = get_vxc_full_response(ni, mol, grids, mf.xc, dm,
max_memory=max_memory,
verbose=ks_grad.verbose)
if mf.nlc:
assert 'VV10' in mf.nlc.upper()
enlc, vnlc = get_vxc_full_response(ni, mol, nlcgrids,
mf.xc+'__'+mf.nlc, dm,
max_memory=max_memory,
verbose=ks_grad.verbose)
exc += enlc
vxc += vnlc
logger.debug1(ks_grad, 'sum(grids response) %s', exc.sum(axis=0))
else:
exc, vxc = get_vxc(ni, mol, grids, mf.xc, dm,
max_memory=max_memory, verbose=ks_grad.verbose)
if mf.nlc:
assert 'VV10' in mf.nlc.upper()
enlc, vnlc = get_vxc(ni, mol, nlcgrids, mf.xc+'__'+mf.nlc, dm,
max_memory=max_memory,
verbose=ks_grad.verbose)
vxc += vnlc
t0 = logger.timer(ks_grad, 'vxc', *t0)

if abs(hyb) < 1e-10 and abs(alpha) < 1e-10:
@@ -109,7 +131,29 @@ def get_vxc(ni, mol, grids, xc_code, dms, relativity=0, hermi=1,
_gga_grad_sum_(vmat[idm], mol, ao, wv, mask, ao_loc)

elif xctype == 'NLC':
raise NotImplementedError('NLC')
nlc_pars = ni.nlc_coeff(xc_code)
ao_deriv = 2
vvrho = []
for ao, mask, weight, coords \
in ni.block_loop(mol, grids, nao, ao_deriv, max_memory):
vvrho.append([make_rho(idm, ao[:4], mask, 'GGA')
for idm in range(nset)])

vv_vxc = []
for idm in range(nset):
rho = numpy.hstack([r[idm] for r in vvrho])
vxc = numint._vv10nlc(rho, grids.coords, rho, grids.weights,
grids.coords, nlc_pars)[1]
vv_vxc.append(xc_deriv.transform_vxc(rho, vxc, 'GGA', spin=0))

p1 = 0
for ao, mask, weight, coords \
in ni.block_loop(mol, grids, nao, ao_deriv, max_memory):
p0, p1 = p1, p1 + weight.size
for idm in range(nset):
wv = vv_vxc[idm][:,p0:p1] * weight
wv[0] *= .5 # *.5 because vmat + vmat.T at the end
_gga_grad_sum_(vmat[idm], mol, ao, wv, mask, ao_loc)

elif xctype == 'MGGA':
ao_deriv = 2
@@ -182,6 +226,69 @@ def _tau_grad_dot_(vmat, mol, ao, wv, mask, ao_loc, dR1_on_bra=True):
aow = numint._scale_ao(ao[3], wv, aow)
_d1_dot_(vmat, mol, [ao[6], ao[8], ao[9]], aow, mask, ao_loc, True)

def _vv10nlc_grad(rho, coords, vvrho, vvweight, vvcoords, nlc_pars):
# VV10 gradient term from Vydrov and Van Voorhis 2010 eq. 25-26
# https://doi.org/10.1063/1.3521275
thresh=1e-8

#output
exc=numpy.zeros((rho[0,:].size,3))

#outer grid needs threshing
threshind=rho[0,:]>=thresh
coords=coords[threshind]
R=rho[0,:][threshind]
Gx=rho[1,:][threshind]
Gy=rho[2,:][threshind]
Gz=rho[3,:][threshind]
G=Gx**2.+Gy**2.+Gz**2.

#inner grid needs threshing
innerthreshind=vvrho[0,:]>=thresh
vvcoords=vvcoords[innerthreshind]
vvweight=vvweight[innerthreshind]
Rp=vvrho[0,:][innerthreshind]
RpW=Rp*vvweight
Gxp=vvrho[1,:][innerthreshind]
Gyp=vvrho[2,:][innerthreshind]
Gzp=vvrho[3,:][innerthreshind]
Gp=Gxp**2.+Gyp**2.+Gzp**2.

#constants and parameters
Pi=numpy.pi
Pi43=4.*Pi/3.
Bvv, Cvv = nlc_pars
Kvv=Bvv*1.5*Pi*((9.*Pi)**(-1./6.))
Beta=((3./(Bvv*Bvv))**(0.75))/32.

#inner grid
W0p=Gp/(Rp*Rp)
W0p=Cvv*W0p*W0p
W0p=(W0p+Pi43*Rp)**0.5
Kp=Kvv*(Rp**(1./6.))

#outer grid
W0tmp=G/(R**2)
W0tmp=Cvv*W0tmp*W0tmp
W0=(W0tmp+Pi43*R)**0.5
K=Kvv*(R**(1./6.))

vvcoords = numpy.asarray(vvcoords, order='C')
coords = numpy.asarray(coords, order='C')
F = numpy.empty(R.shape +(3,), order='C')
libdft.VXC_vv10nlc_grad(F.ctypes.data_as(ctypes.c_void_p),
vvcoords.ctypes.data_as(ctypes.c_void_p),
coords.ctypes.data_as(ctypes.c_void_p),
W0p.ctypes.data_as(ctypes.c_void_p),
W0.ctypes.data_as(ctypes.c_void_p),
K.ctypes.data_as(ctypes.c_void_p),
Kp.ctypes.data_as(ctypes.c_void_p),
RpW.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(vvcoords.shape[0]),
ctypes.c_int(coords.shape[0]))
#exc is multiplied by Rho later
exc[threshind] = F
return exc, Beta

def get_vxc_full_response(ni, mol, grids, xc_code, dms, relativity=0, hermi=1,
max_memory=2000, verbose=None):
@@ -236,7 +343,54 @@ def get_vxc_full_response(ni, mol, grids, xc_code, dms, relativity=0, hermi=1,
rho = vxc = wv = None

elif xctype == 'NLC':
raise NotImplementedError('NLC')
nlc_pars = ni.nlc_coeff(xc_code)
ao_deriv = 2
vvrho = []
vvcoords = []
vvweights = []
for atm_id, (coords, weight) in enumerate(grids_noresponse_cc(grids)):
mask = gen_grid.make_mask(mol, coords)
ao = ni.eval_ao(mol, coords, deriv=ao_deriv, non0tab=mask,
cutoff=grids.cutoff)
vvrho.append(make_rho(0, ao[:4], mask, 'GGA'))
vvcoords.append(coords)
vvweights.append(weight)

vv_vxc = []
vvcoords_flat = numpy.vstack(vvcoords)
vvweights_flat = numpy.concatenate(vvweights)
vvrho_flat = numpy.hstack(vvrho)
for atm_id, (coords, weight, weight1) in enumerate(grids_response_cc(grids)):
rho = vvrho[atm_id]
mask = gen_grid.make_mask(mol, coords)
ao = ni.eval_ao(mol, coords, deriv=ao_deriv, non0tab=mask,
cutoff=grids.cutoff)

exc, vxc = numint._vv10nlc(rho, coords, vvrho_flat, vvweights_flat,
vvcoords_flat, nlc_pars)
vv_vxc = xc_deriv.transform_vxc(rho, vxc, 'GGA', spin=0)
wv = vv_vxc * weight
wv[0] *= .5
vtmp = numpy.zeros((3,nao,nao))
_gga_grad_sum_(vtmp, mol, ao, wv, mask, ao_loc)
vmat += vtmp

vvrho_sub = numpy.hstack(
[r for i, r in enumerate(vvrho) if i != atm_id])
vvcoords_sub = numpy.vstack(
[r for i, r in enumerate(vvcoords) if i != atm_id])
vvweights_sub = numpy.concatenate(
[r for i, r in enumerate(vvweights) if i != atm_id])
egrad, Beta = _vv10nlc_grad(rho, coords, vvrho_sub,
vvweights_sub, vvcoords_sub,
nlc_pars)
# account for factor of 2 in double integration
exc -= 0.5 * Beta
# response of weights
excsum += 2 * numpy.einsum('r,r,nxr->nx', exc, rho[0], weight1)
# response of grids coordinates
excsum[atm_id] += 2 * numpy.einsum('xij,ji->x', vtmp, dms)
excsum[atm_id] += numpy.einsum('r,rx->x', rho[0]*weight, egrad)

elif xctype == 'MGGA':
ao_deriv = 2
@@ -378,6 +532,23 @@ def get_du(ia, ib): # JCP 98, 5612 (1993); (B10)
yield coords, w0, w1


def grids_noresponse_cc(grids):
# same as above but without the response, for nlc grids response routine
assert grids.becke_scheme == gen_grid.original_becke
mol = grids.mol
atom_grids_tab = grids.gen_atomic_grids(mol, grids.atom_grid,
grids.radi_method,
grids.level, grids.prune)
coords_all, weights_all = gen_grid.get_partition(mol, atom_grids_tab,
grids.radii_adjust,
grids.atomic_radii,
grids.becke_scheme,
concat=False)
natm = mol.natm
for ia in range(natm):
yield coords_all[ia], weights_all[ia]


class Gradients(rhf_grad.Gradients):

# This parameter has no effects for HF gradients. Add this attribute so that
@@ -387,10 +558,11 @@ class Gradients(rhf_grad.Gradients):
def __init__(self, mf):
rhf_grad.Gradients.__init__(self, mf)
self.grids = None
self.nlcgrids = None
# This parameter has no effects for HF gradients. Add this attribute so that
# the kernel function can be reused in the DFT gradients code.
self.grid_response = False
self._keys = self._keys.union(['grid_response', 'grids'])
self._keys = self._keys.union(['grid_response', 'grids', 'nlcgrids'])

def dump_flags(self, verbose=None):
rhf_grad.Gradients.dump_flags(self, verbose)
21 changes: 21 additions & 0 deletions pyscf/grad/test/test_rks.py
Original file line number Diff line number Diff line change
@@ -246,6 +246,27 @@ def test_finite_diff_rks_grad_gga(self):
e2 = mf_scanner(mol1.set_geom_('O 0. 0. -.001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))
self.assertAlmostEqual(g[0,2], (e1-e2)/2e-3*lib.param.BOHR, 5)

def test_finite_diff_rks_grad_nlc(self):
#[[ 2.91036539e-16 1.22693574e-15 2.45978284e-02]
# [ 2.83888198e-17 2.66388957e-02 -1.23039325e-02]
# [ 1.17327811e-16 -2.66388957e-02 -1.23039325e-02]]
mf = mol.RKS()
mf.set(xc='VV10', nlc='VV10', conv_tol=1e-12)
mf.nlcgrids.level = 1
mf.kernel()
g = mf.nuc_grad_method().set().kernel()
self.assertAlmostEqual(lib.fp(g), -0.049431714073528615, 6)

mf.nlcgrids.level = 0
mf.kernel()
g = mf.nuc_grad_method().set(grid_response=True).kernel()

mol1 = mol.copy()
mf_scanner = mf.as_scanner()
e1 = mf_scanner(mol1.set_geom_('O 0. 0. 0.001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))
e2 = mf_scanner(mol1.set_geom_('O 0. 0. -.001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))
self.assertAlmostEqual(g[0,2], (e1-e2)/2e-3*lib.param.BOHR, 5)

def test_finite_diff_rks_grad_mgga(self):
mf = mol.RKS().run(xc='m06l', conv_tol=1e-12)
g = mf.nuc_grad_method().set(grid_response=True).kernel()
21 changes: 21 additions & 0 deletions pyscf/grad/test/test_uks.py
Original file line number Diff line number Diff line change
@@ -122,6 +122,27 @@ def test_finite_diff_uks_grad_gga(self):
e2 = mf_scanner(mol1.set_geom_('O 0. 0. -.0001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))
self.assertAlmostEqual(g[0,2], (e1-e2)/2e-4*lib.param.BOHR, 6)

def test_finite_diff_uks_grad_nlc(self):
#[[ 3.19690405e-16 -9.39540337e-16 5.09520937e-02]
# [-5.47247224e-17 6.32050882e-02 -2.54793086e-02]
# [ 6.75779981e-17 -6.32050882e-02 -2.54793086e-02]]
mf = mol.UKS()
mf.set(xc='VV10', conv_tol=1e-12, nlc='VV10')
mf.nlcgrids.level = 1
mf.kernel()
g = mf.nuc_grad_method().kernel()
self.assertAlmostEqual(lib.fp(g), -0.11368788988328639, 6)

mf.nlcgrids.level = 0
mf.kernel()
g = mf.Gradients().set(grid_response=True).kernel()

mol1 = mol.copy()
mf_scanner = mf.as_scanner()
e1 = mf_scanner(mol1.set_geom_('O 0. 0. 0.0001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))
e2 = mf_scanner(mol1.set_geom_('O 0. 0. -.0001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))
self.assertAlmostEqual(g[0,2], (e1-e2)/2e-4*lib.param.BOHR, 6)

def test_finite_diff_uks_grad_mgga(self):
mf = mol.UKS().run(xc='m06l', conv_tol=1e-12)
g = mf.nuc_grad_method().set(grid_response=True).kernel()
25 changes: 22 additions & 3 deletions pyscf/grad/uks.py
Original file line number Diff line number Diff line change
@@ -45,11 +45,16 @@ def get_veff(ks_grad, mol=None, dm=None):
grids = ks_grad.grids
else:
grids = mf.grids
if mf.nlc != '':
if ks_grad.nlcgrids is not None:
nlcgrids = ks_grad.nlcgrids
else:
nlcgrids = mf.nlcgrids
if nlcgrids.coords is None:
nlcgrids.build(with_non0tab=True)
if grids.coords is None:
grids.build(with_non0tab=True)

if mf.nlc != '':
raise NotImplementedError
#enabling range-separated hybrids
omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, spin=mol.spin)

@@ -60,9 +65,22 @@ def get_veff(ks_grad, mol=None, dm=None):
max_memory=max_memory,
verbose=ks_grad.verbose)
logger.debug1(ks_grad, 'sum(grids response) %s', exc.sum(axis=0))
if mf.nlc:
assert 'VV10' in mf.nlc.upper()
enlc, vnlc = rks_grad.get_vxc_full_response(
ni, mol, nlcgrids, mf.xc+'__'+mf.nlc, dm[0]+dm[1],
max_memory=max_memory, verbose=ks_grad.verbose)
exc += enlc
vxc += vnlc
else:
exc, vxc = get_vxc(ni, mol, grids, mf.xc, dm,
max_memory=max_memory, verbose=ks_grad.verbose)
if mf.nlc:
assert 'VV10' in mf.nlc.upper()
enlc, vnlc = rks_grad.get_vxc(ni, mol, nlcgrids, mf.xc+'__'+mf.nlc,
dm[0]+dm[1], max_memory=max_memory,
verbose=ks_grad.verbose)
vxc += vnlc
t0 = logger.timer(ks_grad, 'vxc', *t0)

if abs(hyb) < 1e-10:
@@ -237,8 +255,9 @@ class Gradients(uhf_grad.Gradients):
def __init__(self, mf):
uhf_grad.Gradients.__init__(self, mf)
self.grids = None
self.nlcgrids = None
self.grid_response = False
self._keys = self._keys.union(['grid_response', 'grids'])
self._keys = self._keys.union(['grid_response', 'grids', 'nlcgrids'])

def dump_flags(self, verbose=None):
uhf_grad.Gradients.dump_flags(self, verbose)
35 changes: 35 additions & 0 deletions pyscf/lib/dft/nr_numint.c
Original file line number Diff line number Diff line change
@@ -278,3 +278,38 @@ void VXC_vv10nlc(double *Fvec, double *Uvec, double *Wvec,
}
}
}

void VXC_vv10nlc_grad(double *Fvec, double *vvcoords, double *coords,
double *W0p, double *W0, double *K, double *Kp, double *RpW,
int vvngrids, int ngrids)
{
#pragma omp parallel
{
double DX, DY, DZ, R2;
double gp, g, gt, T, Q, FX, FY, FZ;
int i, j;
#pragma omp for schedule(static)
for (i = 0; i < ngrids; i++) {
FX = 0;
FY = 0;
FZ = 0;
for (j = 0; j < vvngrids; j++) {
DX = vvcoords[j*3+0] - coords[i*3+0];
DY = vvcoords[j*3+1] - coords[i*3+1];
DZ = vvcoords[j*3+2] - coords[i*3+2];
R2 = DX*DX + DY*DY + DZ*DZ;
gp = R2*W0p[j] + Kp[j];
g = R2*W0[i] + K[i];
gt = g + gp;
T = RpW[j] / (g*gp*gt);
Q = T * (W0[i]/g + W0p[j]/gp + (W0[i]+W0p[j])/gt);
FX += Q * DX;
FY += Q * DY;
FZ += Q * DZ;
}
Fvec[i*3+0] = FX * -3;
Fvec[i*3+1] = FY * -3;
Fvec[i*3+2] = FZ * -3;
}
}
}

0 comments on commit 654cc36

Please sign in to comment.