#!/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>
#
'''
Analytical nuclear gradients for domain decomposition COSMO
See also
[1] Fast Domain Decomposition Algorithm for Continuum Solvation Models: Energy and First Derivatives.
F. Lipparini, B. Stamm, E. Cances, Y. Maday, B. Mennucci
J. Chem. Theory Comput., 9, 3637-3648 (2013)
http://dx.doi.org/10.1021/ct400280b
[2] Quantum, classical, and hybrid QM/MM calculations in solution: General implementation of the ddCOSMO linear scaling strategy.
F. Lipparini, G. Scalmani, L. Lagardere, B. Stamm, E. Cances, Y. Maday, J.-P.Piquemal, M. J. Frisch, B. Mennucci
J. Chem. Phys., 141, 184108 (2014)
http://dx.doi.org/10.1063/1.4901304
''' # noqa: E501
import numpy
from pyscf import lib
from pyscf.lib import logger
from pyscf import gto
from pyscf import df
from pyscf.dft import gen_grid, numint
from pyscf.symm import sph
from pyscf.solvent import ddcosmo
from pyscf.solvent._attach_solvent import _Solvation
from pyscf.grad import rhf as rhf_grad
from pyscf.grad import rks as rks_grad
from pyscf.grad import tdrhf as tdrhf_grad # noqa
# TODO: Define attribute grad_method.base to point to the class of the 0th
# order calculation for all gradients class. Then this function can be
# extended and used as the general interface to initialize solvent gradients.
[docs]
def make_grad_object(grad_method):
'''For grad_method in vacuum, add nuclear gradients of solvent pcmobj'''
# Zeroth order method object must be a solvation-enabled method
assert isinstance(grad_method.base, _Solvation)
if grad_method.base.with_solvent.frozen:
raise RuntimeError('Frozen solvent model is not avialbe for energy gradients')
name = (grad_method.base.with_solvent.__class__.__name__
+ grad_method.__class__.__name__)
return lib.set_class(WithSolventGrad(grad_method),
(WithSolventGrad, grad_method.__class__), name)
[docs]
class WithSolventGrad:
_keys = set(['de_solvent', 'de_solute'])
def __init__(self, grad_method):
self.__dict__.update(grad_method.__dict__)
self.de_solvent = None
self.de_solute = None
[docs]
def undo_solvent(self):
cls = self.__class__
name_mixin = self.base.with_solvent.__class__.__name__
obj = lib.view(self, lib.drop_class(cls, WithSolventGrad, name_mixin))
del obj.de_solvent
del obj.de_solute
return obj
[docs]
def kernel(self, *args, dm=None, **kwargs):
if dm is None:
dm = self.base.make_rdm1(ao_repr=True)
self.de_solvent = kernel(self.base.with_solvent, dm)
self.de_solute = super().kernel(*args, **kwargs)
self.de = self.de_solute + self.de_solvent
if self.verbose >= logger.NOTE:
logger.note(self, '--------------- %s (+%s) gradients ---------------',
self.base.__class__.__name__,
self.base.with_solvent.__class__.__name__)
rhf_grad._write(self, self.mol, self.de, self.atmlst)
logger.note(self, '----------------------------------------------')
return self.de
def _finalize(self):
# disable _finalize. It is called in grad_method.kernel method
# where self.de was not yet initialized.
pass
[docs]
def kernel(pcmobj, dm, verbose=None):
mol = pcmobj.mol
natm = mol.natm
lmax = pcmobj.lmax
if pcmobj.grids.coords is None:
pcmobj.grids.build(with_non0tab=True)
if not (isinstance(dm, numpy.ndarray) and dm.ndim == 2):
# UHF density matrix
dm = dm[0] + dm[1]
r_vdw = ddcosmo.get_atomic_radii(pcmobj)
coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order)
ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, lmax, True))
fi = ddcosmo.make_fi(pcmobj, r_vdw)
ui = 1 - fi
ui[ui<0] = 0
cached_pol = ddcosmo.cache_fake_multipoles(pcmobj.grids, r_vdw, lmax)
nlm = (lmax+1)**2
L0 = ddcosmo.make_L(pcmobj, r_vdw, ylm_1sph, fi)
L0 = L0.reshape(natm*nlm,-1)
L1 = make_L1(pcmobj, r_vdw, ylm_1sph, fi)
phi0 = ddcosmo.make_phi(pcmobj, dm, r_vdw, ui, ylm_1sph)
phi1 = make_phi1(pcmobj, dm, r_vdw, ui, ylm_1sph)
L0_X = numpy.linalg.solve(L0, phi0.ravel()).reshape(natm,-1)
psi0, vmat, L0_S = \
ddcosmo.make_psi_vmat(pcmobj, dm, r_vdw, ui, ylm_1sph,
cached_pol, L0_X, L0)
e_psi1 = make_e_psi1(pcmobj, dm, r_vdw, ui, ylm_1sph,
cached_pol, L0_X, L0)
dielectric = pcmobj.eps
if dielectric > 0:
f_epsilon = (dielectric-1.)/dielectric
else:
f_epsilon = 1
de = .5 * f_epsilon * e_psi1
de+= .5 * f_epsilon * numpy.einsum('jx,azjx->az', L0_S, phi1)
de-= .5 * f_epsilon * numpy.einsum('aziljm,il,jm->az', L1, L0_S, L0_X)
return de
[docs]
def make_L1(pcmobj, r_vdw, ylm_1sph, fi):
# See JCTC, 9, 3637, Eq (18)
mol = pcmobj.mol
natm = mol.natm
lmax = pcmobj.lmax
eta = pcmobj.eta
nlm = (lmax+1)**2
coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order)
ngrid_1sph = weights_1sph.size
atom_coords = mol.atom_coords()
ylm_1sph = ylm_1sph.reshape(nlm,ngrid_1sph)
Lmat = numpy.zeros((natm,3,natm,nlm,natm,nlm))
fi1 = make_fi1(pcmobj, pcmobj.get_atomic_radii())
for ja in range(natm):
part_weights = weights_1sph.copy()
part_weights[fi[ja]>1] /= fi[ja,fi[ja]>1]
part_weights1 = numpy.zeros((natm,3,ngrid_1sph))
tmp = part_weights[fi[ja]>1] / fi[ja,fi[ja]>1]
part_weights1[:,:,fi[ja]>1] = -tmp * fi1[:,:,ja,fi[ja]>1]
for ka in ddcosmo.atoms_with_vdw_overlap(ja, atom_coords, r_vdw):
vjk = r_vdw[ja] * coords_1sph + atom_coords[ja] - atom_coords[ka]
rv = lib.norm(vjk, axis=1)
tjk = rv / r_vdw[ka]
wjk0 = pcmobj.regularize_xt(tjk, eta)
wjk1 = regularize_xt1(tjk, eta)
sjk = vjk.T / rv
wjk1 = 1./r_vdw[ka] * wjk1 * sjk
wjk01 = wjk0 * part_weights1
wjk0 *= part_weights
wjk1 *= part_weights
pol0 = sph.multipoles(vjk, lmax)
pol1 = multipoles1(vjk, lmax)
p1 = 0
for l in range(lmax+1):
fac = 4*numpy.pi/(l*2+1) / r_vdw[ka]**(l+1)
p0, p1 = p1, p1 + (l*2+1)
a = numpy.einsum('xn,zn,mn->zxm', ylm_1sph, wjk1, pol0[l])
a+= numpy.einsum('xn,n,zmn->zxm', ylm_1sph, wjk0, pol1[l])
Lmat[ja,:,ja,:,ka,p0:p1] += -fac * a
Lmat[ka,:,ja,:,ka,p0:p1] -= -fac * a
a = numpy.einsum('xn,azn,mn->azxm', ylm_1sph, wjk01, pol0[l])
Lmat[:,:,ja,:,ka,p0:p1] += -fac * a
return Lmat
[docs]
def multipoles1(r, lmax, reorder_dipole=True):
ngrid = r.shape[0]
xs = numpy.ones((lmax+1,ngrid))
ys = numpy.ones((lmax+1,ngrid))
zs = numpy.ones((lmax+1,ngrid))
for i in range(1,lmax+1):
xs[i] = xs[i-1] * r[:,0]
ys[i] = ys[i-1] * r[:,1]
zs[i] = zs[i-1] * r[:,2]
ylms = []
for l in range(lmax+1):
nd = (l+1)*(l+2)//2
c = numpy.empty((nd,3,ngrid))
k = 0
for lx in reversed(range(0, l+1)):
for ly in reversed(range(0, l-lx+1)):
lz = l - lx - ly
c[k,0] = lx * xs[lx-1] * ys[ly] * zs[lz]
c[k,1] = ly * xs[lx] * ys[ly-1] * zs[lz]
c[k,2] = lz * xs[lx] * ys[ly] * zs[lz-1]
k += 1
ylm = gto.cart2sph(l, c.reshape(nd,3*ngrid).T)
ylm = ylm.reshape(3,ngrid,l*2+1).transpose(0,2,1)
ylms.append(ylm)
# when call libcint, p functions are ordered as px,py,pz
# reorder px,py,pz to p(-1),p(0),p(1)
if (not reorder_dipole) and lmax >= 1:
ylms[1] = ylms[1][:,[1,2,0]]
return ylms
[docs]
def regularize_xt1(t, eta):
xt = numpy.zeros_like(t)
# no response if grids are inside the cavity
# inner = t <= 1-eta
# xt[inner] = 0
on_shell = (1-eta < t) & (t < 1)
ti = t[on_shell]
xt[on_shell] = -30./eta**5 * (1-ti)**2 * (1-eta-ti)**2
return xt
[docs]
def make_fi1(pcmobj, r_vdw):
coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order)
mol = pcmobj.mol
eta = pcmobj.eta
natm = mol.natm
atom_coords = mol.atom_coords()
ngrid_1sph = coords_1sph.shape[0]
fi1 = numpy.zeros((natm,3,natm,ngrid_1sph))
for ia in range(natm):
for ja in ddcosmo.atoms_with_vdw_overlap(ia, atom_coords, r_vdw):
v = r_vdw[ia]*coords_1sph + atom_coords[ia] - atom_coords[ja]
rv = lib.norm(v, axis=1)
t = rv / r_vdw[ja]
xt1 = regularize_xt1(t, eta)
s_ij = v.T / rv
xt1 = 1./r_vdw[ja] * xt1 * s_ij
fi1[ia,:,ia] += xt1
fi1[ja,:,ia] -= xt1
fi = ddcosmo.make_fi(pcmobj, r_vdw)
fi1[:,:,fi<1e-20] = 0
return fi1
[docs]
def make_phi1(pcmobj, dm, r_vdw, ui, ylm_1sph):
mol = pcmobj.mol
natm = mol.natm
if not (isinstance(dm, numpy.ndarray) and dm.ndim == 2):
dm = dm[0] + dm[1]
tril_dm = lib.pack_tril(dm+dm.T)
nao = dm.shape[0]
diagidx = numpy.arange(nao)
diagidx = diagidx*(diagidx+1)//2 + diagidx
tril_dm[diagidx] *= .5
atom_coords = mol.atom_coords()
atom_charges = mol.atom_charges()
coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order)
#extern_point_idx = ui > 0
fi1 = make_fi1(pcmobj, pcmobj.get_atomic_radii())
fi1[:,:,ui==0] = 0
ui1 = -fi1
ngrid_1sph = weights_1sph.size
v_phi0 = numpy.empty((natm,ngrid_1sph))
for ia in range(natm):
cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph
d_rs = atom_coords.reshape(-1,1,3) - cav_coords
v_phi0[ia] = numpy.einsum('z,zp->p', atom_charges, 1./lib.norm(d_rs,axis=2))
phi1 = -numpy.einsum('n,ln,azjn,jn->azjl', weights_1sph, ylm_1sph, ui1, v_phi0)
for ia in range(natm):
cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph
for ja in range(natm):
rs = atom_coords[ja] - cav_coords
d_rs = lib.norm(rs, axis=1)
v_phi = atom_charges[ja] * numpy.einsum('px,p->px', rs, 1./d_rs**3)
tmp = numpy.einsum('n,ln,n,nx->xl', weights_1sph, ylm_1sph, ui[ia], v_phi)
phi1[ja,:,ia] += tmp # response of the other atoms
phi1[ia,:,ia] -= tmp # response of cavity grids
int3c2e = mol._add_suffix('int3c2e')
int3c2e_ip1 = mol._add_suffix('int3c2e_ip1')
aoslices = mol.aoslice_by_atom()
for ia in range(natm):
cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph
#fakemol = gto.fakemol_for_charges(cav_coords[ui[ia]>0])
fakemol = gto.fakemol_for_charges(cav_coords)
v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1')
v_phi = numpy.einsum('ij,ijk->k', dm, v_nj)
phi1[:,:,ia] += numpy.einsum('n,ln,azn,n->azl', weights_1sph, ylm_1sph, ui1[:,:,ia], v_phi)
v_e1_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1, comp=3, aosym='s1')
phi1_e2_nj = numpy.einsum('ij,xijr->xr', dm, v_e1_nj)
phi1_e2_nj += numpy.einsum('ji,xijr->xr', dm, v_e1_nj)
phi1[ia,:,ia] += numpy.einsum('n,ln,n,xn->xl', weights_1sph, ylm_1sph, ui[ia], phi1_e2_nj)
for ja in range(natm):
shl0, shl1, p0, p1 = aoslices[ja]
phi1_nj = numpy.einsum('ij,xijr->xr', dm[p0:p1 ], v_e1_nj[:,p0:p1])
phi1_nj += numpy.einsum('ji,xijr->xr', dm[:,p0:p1], v_e1_nj[:,p0:p1])
phi1[ja,:,ia] -= numpy.einsum('n,ln,n,xn->xl', weights_1sph, ylm_1sph, ui[ia], phi1_nj)
return phi1
[docs]
def make_e_psi1(pcmobj, dm, r_vdw, ui, ylm_1sph, cached_pol, Xvec, L):
mol = pcmobj.mol
natm = mol.natm
lmax = pcmobj.lmax
grids = pcmobj.grids
if not (isinstance(dm, numpy.ndarray) and dm.ndim == 2):
dm = dm[0] + dm[1]
ni = numint.NumInt()
make_rho, nset, nao = ni._gen_rho_evaluator(mol, dm)
den = numpy.empty((4,grids.weights.size))
ao_loc = mol.ao_loc_nr()
vmat = numpy.zeros((3,nao,nao))
psi1 = numpy.zeros((natm,3))
i1 = 0
for ia, (coords, weight, weight1) in enumerate(rks_grad.grids_response_cc(grids)):
i0, i1 = i1, i1 + weight.size
ao = ni.eval_ao(mol, coords, deriv=1)
mask = gen_grid.make_mask(mol, coords)
den[:,i0:i1] = make_rho(0, ao, mask, 'GGA')
fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)]
eta_nj = 0
p1 = 0
for l in range(lmax+1):
fac = 4*numpy.pi/(l*2+1)
p0, p1 = p1, p1 + (l*2+1)
eta_nj += fac * numpy.einsum('mn,m->n', fak_pol[l], Xvec[ia,p0:p1])
psi1 -= numpy.einsum('n,n,zxn->zx', den[0,i0:i1], eta_nj, weight1)
psi1[ia] -= numpy.einsum('xn,n,n->x', den[1:4,i0:i1], eta_nj, weight)
vtmp = numpy.zeros((3,nao,nao))
aow = numpy.einsum('pi,p->pi', ao[0], weight*eta_nj)
rks_grad._d1_dot_(vtmp, mol, ao[1:4], aow, mask, ao_loc, True)
vmat += vtmp
aoslices = mol.aoslice_by_atom()
for ia in range(natm):
shl0, shl1, p0, p1 = aoslices[ia]
psi1[ia] += numpy.einsum('xij,ij->x', vmat[:,p0:p1], dm[p0:p1]) * 2
return psi1
if __name__ == '__main__':
from pyscf import scf
mol = gto.M(atom='H 0 0 0; H 0 1 1.2; H 1. .1 0; H .5 .5 1', unit='B')
mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol))
mf.kernel()
de = mf.nuc_grad_method().kernel()
de_cosmo = kernel(mf.with_solvent, mf.make_rdm1())
dm1 = mf.make_rdm1()
mol = gto.M(atom='H 0 0 -0.001; H 0 1 1.2; H 1. .1 0; H .5 .5 1', unit='B')
mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol))
e1 = mf.kernel()
e1_cosmo = mf.with_solvent.energy(dm1)
mol = gto.M(atom='H 0 0 0.001; H 0 1 1.2; H 1. .1 0; H .5 .5 1', unit='B')
mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol))
e2 = mf.kernel()
e2_cosmo = mf.with_solvent.energy(dm1)
print(abs((e2-e1)/0.002 - de[0,2]).max())
print(abs((e2_cosmo-e1_cosmo)/0.002 - de_cosmo[0,2]).max())