diff --git a/examples/14-pcm_solvent.py b/examples/14-pcm_solvent.py index 6e17fc17..fef75200 100644 --- a/examples/14-pcm_solvent.py +++ b/examples/14-pcm_solvent.py @@ -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]) + diff --git a/examples/15-chelpg.py b/examples/15-chelpg.py index f102741b..473d500d 100644 --- a/examples/15-chelpg.py +++ b/examples/15-chelpg.py @@ -17,7 +17,6 @@ from gpu4pyscf.dft import rks from gpu4pyscf.qmmm import chelpg - mol = gto.Mole() mol.verbose = 0 mol.output = None @@ -30,7 +29,7 @@ mol.unit = 'B' mol.build() mol.verbose = 6 - + xc = 'b3lyp' mf = rks.RKS(mol, xc=xc) mf.grids.level = 5 diff --git a/gpu4pyscf/grad/rks.py b/gpu4pyscf/grad/rks.py index 39345219..82b0befd 100644 --- a/gpu4pyscf/grad/rks.py +++ b/gpu4pyscf/grad/rks.py @@ -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) diff --git a/gpu4pyscf/scf/hf.py b/gpu4pyscf/scf/hf.py index 8e1a9855..8bcfd65e 100644 --- a/gpu4pyscf/scf/hf.py +++ b/gpu4pyscf/scf/hf.py @@ -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 diff --git a/gpu4pyscf/solvent/_attach_solvent.py b/gpu4pyscf/solvent/_attach_solvent.py index 52f090e1..23e08af2 100644 --- a/gpu4pyscf/solvent/_attach_solvent.py +++ b/gpu4pyscf/solvent/_attach_solvent.py @@ -13,6 +13,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import cupy from pyscf import lib from pyscf.lib import logger from pyscf.solvent._attach_solvent import _Solvation @@ -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) diff --git a/gpu4pyscf/solvent/grad/pcm.py b/gpu4pyscf/solvent/grad/pcm.py index 99518394..f3b782f7 100644 --- a/gpu4pyscf/solvent/grad/pcm.py +++ b/gpu4pyscf/solvent/grad/pcm.py @@ -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) @@ -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() @@ -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) @@ -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) @@ -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 diff --git a/gpu4pyscf/solvent/hessian/pcm.py b/gpu4pyscf/solvent/hessian/pcm.py index 2ad8e137..c068f538 100644 --- a/gpu4pyscf/solvent/hessian/pcm.py +++ b/gpu4pyscf/solvent/hessian/pcm.py @@ -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 ''' @@ -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 @@ -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): @@ -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) @@ -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): diff --git a/gpu4pyscf/solvent/pcm.py b/gpu4pyscf/solvent/pcm.py index 2b348740..a310b1bd 100644 --- a/gpu4pyscf/solvent/pcm.py +++ b/gpu4pyscf/solvent/pcm.py @@ -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]) @@ -302,24 +302,26 @@ 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): ''' @@ -327,20 +329,18 @@ def _get_v(self, dms): ''' 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): @@ -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 diff --git a/gpu4pyscf/solvent/tests/test_pcm_hessian.py b/gpu4pyscf/solvent/tests/test_pcm_hessian.py index 08cb1cc0..756ca285 100644 --- a/gpu4pyscf/solvent/tests/test_pcm_hessian.py +++ b/gpu4pyscf/solvent/tests/test_pcm_hessian.py @@ -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