Skip to content

Commit

Permalink
fixed bugs for pcm hessian
Browse files Browse the repository at this point in the history
  • Loading branch information
wxj6000 committed Nov 10, 2023
1 parent fe3d206 commit eb5d5eb
Show file tree
Hide file tree
Showing 9 changed files with 131 additions and 51 deletions.
22 changes: 11 additions & 11 deletions examples/14-pcm_solvent.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,24 +34,24 @@
mf.small_rho_cutoff = 1e-10
mf.with_solvent.lebedev_order = 29 # 302 Lebedev grids
mf.with_solvent.method = 'C-PCM'

mf.with_solvent.eps = 78.3553
mf.kernel()

g = mf.nuc_grad_method()
g.auxbasis_response = True
f = g.kernel()
gradobj = mf.nuc_grad_method()
gradobj.auxbasis_response = True
f = gradobj.kernel()

hessobj = mf.Hessian()
hess = hessobj.kernel()

h = mf.Hessian()
hess = h.kernel()
print(hess[0,0])
print(hess[1,0])
print(hess[2,0])

# mass weighted hessian
mass = [15.99491, 1.00783, 1.00783]
for i in range(3):
for j in range(3):
hess[i,j] = hess[i,j]/np.sqrt(mass[i]*mass[j])
n = hess.shape[0]
#hess = hess.transpose([0,2,1,3])#.reshape(3*n,3*n)
print(hess[0,0])
print(hess[1,0])
print(hess[2,0])


3 changes: 1 addition & 2 deletions examples/15-chelpg.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@
from gpu4pyscf.dft import rks
from gpu4pyscf.qmmm import chelpg


mol = gto.Mole()
mol.verbose = 0
mol.output = None
Expand All @@ -30,7 +29,7 @@
mol.unit = 'B'
mol.build()
mol.verbose = 6

xc = 'b3lyp'
mf = rks.RKS(mol, xc=xc)
mf.grids.level = 5
Expand Down
2 changes: 0 additions & 2 deletions gpu4pyscf/grad/rks.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,9 +28,7 @@
from gpu4pyscf.dft import numint, xc_deriv, rks
from gpu4pyscf.dft.numint import _GDFTOpt, AO_THRESHOLD
from gpu4pyscf.lib.cupy_helper import contract, get_avail_mem, add_sparse, tag_array, load_library

from gpu4pyscf.lib import logger

from pyscf import __config__

MIN_BLK_SIZE = getattr(__config__, 'min_grid_blksize', 128*128)
Expand Down
1 change: 0 additions & 1 deletion gpu4pyscf/scf/hf.py
Original file line number Diff line number Diff line change
Expand Up @@ -467,7 +467,6 @@ def _gen_rhf_response(mf, mo_coeff=None, mo_occ=None,
orbital hessian or CPHF will be generated. If singlet is boolean,
it is used in TDDFT response kernel.
'''

if mo_coeff is None: mo_coeff = mf.mo_coeff
if mo_occ is None: mo_occ = mf.mo_occ
mol = mf.mol
Expand Down
2 changes: 1 addition & 1 deletion gpu4pyscf/solvent/_attach_solvent.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.

import cupy
from pyscf import lib
from pyscf.lib import logger
from pyscf.solvent._attach_solvent import _Solvation
Expand Down Expand Up @@ -107,7 +108,6 @@ def gen_response(self, *args, **kwargs):
is_uhf = isinstance(self, scf.uhf.UHF)
# singlet=None is orbital hessian or CPHF type response function
singlet = kwargs.get('singlet', True)

singlet = singlet or singlet is None
def vind_with_solvent(dm1):
v = vind(dm1)
Expand Down
42 changes: 38 additions & 4 deletions gpu4pyscf/solvent/grad/pcm.py
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,36 @@ def get_dD_dS(surface, dF, with_S=True, with_D=False):

return dD, dS, dSii

def grad_kernel(pcmobj, dm):
def grad_nuc(pcmobj):
if not pcmobj._intermediates:
pcmobj.build()
mol = pcmobj.mol
q_sym = pcmobj._intermediates['q_sym'].get()
gridslice = pcmobj.surface['gslice_by_atom']
grid_coords = pcmobj.surface['grid_coords'].get()
exponents = pcmobj.surface['charge_exp'].get()

atom_coords = mol.atom_coords(unit='B')
atom_charges = numpy.asarray(mol.atom_charges(), dtype=numpy.float64)
fakemol_nuc = gto.fakemol_for_charges(atom_coords)
fakemol = gto.fakemol_for_charges(grid_coords, expnt=exponents**2)

int2c2e_ip1 = mol._add_suffix('int2c2e_ip1')

v_ng_ip1 = gto.mole.intor_cross(int2c2e_ip1, fakemol_nuc, fakemol)

dv_g = numpy.einsum('g,xng->nx', q_sym, v_ng_ip1)
de = -numpy.einsum('nx,n->nx', dv_g, atom_charges)

v_ng_ip1 = gto.mole.intor_cross(int2c2e_ip1, fakemol, fakemol_nuc)

dv_g = numpy.einsum('n,xgn->gx', atom_charges, v_ng_ip1)
dv_g = numpy.einsum('gx,g->gx', dv_g, q_sym)

de -= numpy.asarray([numpy.sum(dv_g[p0:p1], axis=0) for p0,p1 in gridslice])
return de

def grad_elec(pcmobj, dm):
'''
dE = 0.5*v* d(K^-1 R) *v + q*dv
v^T* d(K^-1 R)v = v^T*K^-1(dR - dK K^-1R)v = v^T K^-1(dR - dK q)
Expand All @@ -163,7 +192,6 @@ def grad_kernel(pcmobj, dm):
vK_1 = cupy.linalg.solve(K.T, v_grids)

# ----------------- potential response -----------------------
atom_coords = mol.atom_coords(unit='B')

intopt = pcmobj.intopt
intopt.clear()
Expand All @@ -184,6 +212,10 @@ def grad_kernel(pcmobj, dm):
dvj= 2.0 * cupy.asarray([cupy.sum(dvj[:,p0:p1], axis=1) for p0,p1 in aoslice[:,2:]])
de = dq + dvj

#de += grad_nuc(mol, pcmobj.surface, q_sym)

'''
atom_coords = mol.atom_coords(unit='B')
atom_charges = cupy.asarray(mol.atom_charges(), dtype=numpy.float64)
fakemol_nuc = gto.fakemol_for_charges(atom_coords)
fakemol = gto.fakemol_for_charges(grid_coords.get(), expnt=exponents.get()**2)
Expand All @@ -205,7 +237,7 @@ def grad_kernel(pcmobj, dm):
dv_g = contract('gx,g->gx', dv_g, q_sym)
de -= cupy.asarray([cupy.sum(dv_g[p0:p1], axis=0) for p0,p1 in gridslice])

'''
## --------------- response from stiffness matrices ----------------
gridslice = pcmobj.surface['gslice_by_atom']
dF, dA = get_dF_dA(pcmobj.surface)
Expand Down Expand Up @@ -284,7 +316,9 @@ def kernel(self, *args, dm=None, atmlst=None, **kwargs):
if dm is None:
dm = self.base.make_rdm1(ao_repr=True)

self.de_solvent = grad_kernel(self.base.with_solvent, dm)
self.de_solvent = grad_elec(self.base.with_solvent, dm)
self.de_solvent+= grad_nuc(self.base.with_solvent)

self.de_solute = grad_method_class.kernel(self, *args, **kwargs)
self.de = self.de_solute + self.de_solvent

Expand Down
62 changes: 58 additions & 4 deletions gpu4pyscf/solvent/hessian/pcm.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,14 +25,61 @@
from pyscf import gto, df
from pyscf.grad import rhf as rhf_grad
from gpu4pyscf.solvent.pcm import PI, switch_h
from gpu4pyscf.solvent.grad.pcm import grad_switch_h, get_dF_dA, get_dD_dS, grad_kernel
from gpu4pyscf.solvent.grad.pcm import grad_switch_h, get_dF_dA, get_dD_dS, grad_elec, grad_nuc
from gpu4pyscf.df import int3c2e
from gpu4pyscf.lib.cupy_helper import contract
from gpu4pyscf.lib import logger

libdft = lib.load_library('libdft')

def hess_kernel(pcmobj, dm, verbose=None):
def hess_nuc(pcmobj):
if not pcmobj._intermediates:
pcmobj.build()
mol = pcmobj.mol
q_sym = pcmobj._intermediates['q_sym'].get()
gridslice = pcmobj.surface['gslice_by_atom']
grid_coords = pcmobj.surface['grid_coords'].get()
exponents = pcmobj.surface['charge_exp'].get()

atom_coords = mol.atom_coords(unit='B')
atom_charges = numpy.asarray(mol.atom_charges(), dtype=numpy.float64)
fakemol_nuc = gto.fakemol_for_charges(atom_coords)
fakemol = gto.fakemol_for_charges(grid_coords, expnt=exponents**2)

# nuclei potential response
int2c2e_ip1ip2 = mol._add_suffix('int2c2e_ip1ip2')
v_ng_ip1ip2 = gto.mole.intor_cross(int2c2e_ip1ip2, fakemol_nuc, fakemol).reshape([3,3,mol.natm,-1])
charge3 = numpy.tile(atom_charges, [3,1])

q3 = numpy.tile(q_sym, [3,1])
dv_g = numpy.einsum('xn,xyng->ngxy', charge3, v_ng_ip1ip2)
dv_g = numpy.einsum('ngxy,yg->ngxy', dv_g, q3)

de = numpy.zeros([mol.natm, mol.natm, 3, 3])
for ia in range(mol.natm):
p0, p1 = gridslice[ia]
de_tmp = numpy.sum(dv_g[:,p0:p1], axis=1)
de[:,ia] -= de_tmp
de[ia,:] -= de_tmp.transpose([0,2,1])

int2c2e_ipip1 = mol._add_suffix('int2c2e_ipip1')
v_ng_ipip1 = gto.mole.intor_cross(int2c2e_ipip1, fakemol_nuc, fakemol).reshape([3,3,mol.natm,-1])

dv_g = numpy.einsum('g,xyng->nxy', q_sym, v_ng_ipip1)
for ia in range(mol.natm):
de[ia,ia] -= dv_g[ia] * atom_charges[ia]

v_ng_ipip1 = gto.mole.intor_cross(int2c2e_ipip1, fakemol, fakemol_nuc).reshape([3,3,-1,mol.natm])

dv_g = numpy.einsum('n,xygn->gxy', atom_charges, v_ng_ipip1)
dv_g = numpy.einsum('g,gxy->gxy', q_sym, dv_g)
for ia in range(mol.natm):
p0, p1 = gridslice[ia]
de[ia,ia] -= numpy.sum(dv_g[p0:p1], axis=0)

return de

def hess_elec(pcmobj, dm, verbose=None):
'''
slow version with finite difference
'''
Expand All @@ -43,7 +90,8 @@ def hess_kernel(pcmobj, dm, verbose=None):
def pcm_grad_scanner(mol):
pcmobj.reset(mol)
e, v = pcmobj._get_vind(dm)
return grad_kernel(pcmobj, dm)
#return grad_elec(pcmobj, dm)
return grad_nuc(pcmobj) + grad_elec(pcmobj, dm)

de = numpy.zeros([mol.natm, mol.natm, 3, 3])
eps = 1e-3
Expand All @@ -60,6 +108,7 @@ def pcm_grad_scanner(mol):
g1 = pcm_grad_scanner(mol)
de[ia,:,ix] = (g0 - g1)/2.0/eps
t1 = log.timer_debug1('solvent energy', *t1)

return de

def grad_qv(pcmobj, mo_coeff, mo_occ, atmlst=None, verbose=None):
Expand Down Expand Up @@ -94,6 +143,7 @@ def pcm_vmat_scanner(mol):
mol.set_geom_(coords - dv, unit='Bohr')
mol.build()
vmat1 = pcm_vmat_scanner(mol)

grad_vmat = (vmat0 - vmat1)/2.0/eps
grad_vmat = contract("ij,jq->iq", grad_vmat, mocc)
grad_vmat = contract("iq,ip->pq", grad_vmat, mo_coeff)
Expand All @@ -117,9 +167,13 @@ def kernel(self, *args, dm=None, atmlst=None, **kwargs):
dm = kwargs.pop('dm', None)
if dm is None:
dm = self.base.make_rdm1(ao_repr=True)
self.de_solvent = hess_kernel(self.base.with_solvent, dm, verbose=self.verbose)
is_equilibrium = self.base.with_solvent.equilibrium_solvation
self.base.with_solvent.equilibrium_solvation = True
self.de_solvent = hess_elec(self.base.with_solvent, dm, verbose=self.verbose)
#self.de_solvent+= hess_nuc(self.base.with_solvent)
self.de_solute = hess_method_class.kernel(self, *args, **kwargs)
self.de = self.de_solute + self.de_solvent
self.base.with_solvent.equilibrium_solvation = is_equilibrium
return self.de

def make_h1(self, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None):
Expand Down
46 changes: 21 additions & 25 deletions gpu4pyscf/solvent/pcm.py
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ def gen_surface(mol, ng=302, vdw_scale=1.2):
swf = cupy.prod(fiJ, axis=1)
idx = w*swf > 1e-16

p0, p1 = p1, p1+sum(idx)
p0, p1 = p1, p1+sum(idx).get()
gslice_by_atom.append([p0,p1])
grid_coords.append(atom_grid[idx,:3])
weights.append(w[idx])
Expand Down Expand Up @@ -302,45 +302,45 @@ def _get_vind(self, dms):

K = self._intermediates['K']
R = self._intermediates['R']
v_grids = self._get_v(dms)
v_grids_e = self._get_v(dms)
v_grids = self.v_grids_n - v_grids_e

b = cupy.dot(R, v_grids.T)
q = cupy.linalg.solve(K, b.T)
q = cupy.linalg.solve(K, b).T

vK_1 = cupy.linalg.solve(K.T, v_grids)
q_sym = (q + cupy.dot(R.T, vK_1))/2.0
vK_1 = cupy.linalg.solve(K.T, v_grids.T)
qt = cupy.dot(R.T, vK_1).T
q_sym = (q + qt)/2.0

vmat = self._get_vmat(q_sym)
epcm = 0.5 * cupy.dot(q_sym, v_grids)
epcm = 0.5 * cupy.dot(v_grids[0], q_sym[0])

self._intermediates['K'] = K
self._intermediates['R'] = R
self._intermediates['q'] = q
self._intermediates['q_sym'] = q_sym
self._intermediates['v_grids'] = v_grids
self._intermediates['q'] = q[0]
self._intermediates['q_sym'] = q_sym[0]
self._intermediates['v_grids'] = v_grids[0]

return epcm, vmat
return epcm, vmat[0]

def _get_v(self, dms):
'''
return electrostatic potential on surface
'''
nset = dms.shape[0]
ngrids = self.surface['grid_coords'].shape[0]
v_grids = cupy.empty([nset, ngrids])
v_grids_e = cupy.empty([nset, ngrids])
for i in range(nset):
v_grids_e = 2.0*int3c2e.get_j_int3c2e_pass1(self.intopt, dms[0])
v_grids[i] = self.v_grids_n - v_grids_e
return v_grids
v_grids_e[i] = 2.0*int3c2e.get_j_int3c2e_pass1(self.intopt, dms[i])
return v_grids_e

def _get_vmat(self, q):
assert q.ndim == 2
nset = q.shape[0]
nao = self.mol.nao
vmat = cupy.empty([nset, nao, nao])
for i in range(nset):
vmat[i] = -int3c2e.get_j_int3c2e_pass2(self.intopt, q[i])
if nset == 1:
return vmat[0]
return vmat

def nuc_grad_method(self, grad_method):
Expand Down Expand Up @@ -377,20 +377,16 @@ def _B_dot_x(self, dms):

K = self._intermediates['K']
R = self._intermediates['R']
v_grids = self._get_v(dms)
v_grids = -self._get_v(dms)

b = cupy.dot(R, v_grids.T)
q = cupy.linalg.solve(K, b.T)
q = cupy.linalg.solve(K, b).T

vK_1 = cupy.linalg.solve(K.T, v_grids.T)
q_sym = (q + cupy.dot(R.T, vK_1))/2.0
qt = cupy.dot(R.T, vK_1).T
q_sym = (q + qt)/2.0

vmat = self._get_vmat(q_sym)

self._intermediates['K'] = K
self._intermediates['R'] = R
self._intermediates['q'] = q
self._intermediates['q_sym'] = q_sym
self._intermediates['v_grids'] = v_grids

return vmat

2 changes: 1 addition & 1 deletion gpu4pyscf/solvent/tests/test_pcm_hessian.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ def setUpModule():
H 0.7570000000 0.0000000000 -0.4696000000
'''
mol.basis = 'sto3g'
#mol.output = '/dev/null'
mol.output = '/dev/null'
mol.build(verbose=0)
epsilon = 35.9
lebedev_order = 29
Expand Down

0 comments on commit eb5d5eb

Please sign in to comment.