From 618913623102f0f5e7ae0ad3f61266139dcf8655 Mon Sep 17 00:00:00 2001 From: "xiaojie.wu" Date: Sun, 12 Nov 2023 15:46:02 +0800 Subject: [PATCH] v0.6.6 --- examples/14-pcm_solvent.py | 4 - examples/dft_driver.py | 8 +- gpu4pyscf/__init__.py | 2 +- gpu4pyscf/scf/cphf.py | 4 +- gpu4pyscf/solvent/hessian/pcm.py | 25 +++- gpu4pyscf/solvent/pcm.py | 3 +- gpu4pyscf/solvent/tests/test_pcm_hessian.py | 142 +++++++------------- 7 files changed, 83 insertions(+), 105 deletions(-) diff --git a/examples/14-pcm_solvent.py b/examples/14-pcm_solvent.py index fef75200..8eaabbf8 100644 --- a/examples/14-pcm_solvent.py +++ b/examples/14-pcm_solvent.py @@ -44,10 +44,6 @@ hessobj = mf.Hessian() hess = hessobj.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): diff --git a/examples/dft_driver.py b/examples/dft_driver.py index 55ef556d..96b1718c 100644 --- a/examples/dft_driver.py +++ b/examples/dft_driver.py @@ -18,7 +18,6 @@ import argparse from pyscf import lib from gpu4pyscf.dft import rks -lib.num_threads(8) parser = argparse.ArgumentParser(description='Run DFT with GPU4PySCF for molecules') parser.add_argument("--input", type=str, default='benzene/coord') @@ -43,11 +42,14 @@ if args.solvent: mf_df = mf_df.PCM() - mf_df.lebedev_order = 29 - mf_df.method = args.solvent + mf_df.with_solvent.lebedev_order = 29 + mf_df.with_solvent.method = args.solvent + mf_df.with_solvent.eps = 78.3553 + mf_df.grids.atom_grid = (99,590) mf_df.direct_scf_tol = 1e-14 mf_df.direct_scf = 1e-14 +mf_df.conv_tol = 1e-12 e_tot = mf_df.kernel() scf_time = time.time() - start_time print(f'compute time for energy: {scf_time:.3f} s') diff --git a/gpu4pyscf/__init__.py b/gpu4pyscf/__init__.py index 681a2e88..550ac39c 100644 --- a/gpu4pyscf/__init__.py +++ b/gpu4pyscf/__init__.py @@ -1,2 +1,2 @@ from . import lib, grad, hessian, solvent, scf, dft -__version__ = '0.6.5' +__version__ = '0.6.6' diff --git a/gpu4pyscf/scf/cphf.py b/gpu4pyscf/scf/cphf.py index eb2b9da5..fdd45453 100644 --- a/gpu4pyscf/scf/cphf.py +++ b/gpu4pyscf/scf/cphf.py @@ -28,7 +28,7 @@ from gpu4pyscf.lib import logger def solve(fvind, mo_energy, mo_occ, h1, s1=None, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): + max_cycle=50, tol=1e-9, hermi=False, verbose=logger.WARN): ''' Args: fvind : function @@ -70,7 +70,7 @@ def vind_vo(mo1): # h1 shape is (:,nocc+nvir,nocc) def solve_withs1(fvind, mo_energy, mo_occ, h1, s1, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): + max_cycle=50, tol=1e-9, hermi=False, verbose=logger.WARN): '''For field dependent basis. First order overlap matrix is non-zero. The first order orbitals are set to C^1_{ij} = -1/2 S1 diff --git a/gpu4pyscf/solvent/hessian/pcm.py b/gpu4pyscf/solvent/hessian/pcm.py index f93aa598..5d2c4dac 100644 --- a/gpu4pyscf/solvent/hessian/pcm.py +++ b/gpu4pyscf/solvent/hessian/pcm.py @@ -57,10 +57,23 @@ def hess_nuc(pcmobj): 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]) + #de[ia,:] -= de_tmp.transpose([0,2,1]) + + + int2c2e_ip1ip2 = mol._add_suffix('int2c2e_ip1ip2') + v_ng_ip1ip2 = gto.mole.intor_cross(int2c2e_ip1ip2, fakemol, fakemol_nuc).reshape([3,3,-1,mol.natm]) + dv_g = numpy.einsum('n,xygn->gnxy', atom_charges, v_ng_ip1ip2) + dv_g = numpy.einsum('gnxy,g->gnxy', dv_g, q_sym) + + for ia in range(mol.natm): + p0, p1 = gridslice[ia] + de_tmp = numpy.sum(dv_g[p0:p1], axis=0) + 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] @@ -81,7 +94,7 @@ def hess_elec(pcmobj, dm, verbose=None): log = logger.new_logger(pcmobj, verbose) t1 = log.init_timer() pmol = pcmobj.mol.copy() - mol = pmol + mol = pmol.copy() coords = mol.atom_coords(unit='Bohr') def pcm_grad_scanner(mol): @@ -105,15 +118,18 @@ 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) + pcmobj.reset(pmol) return de def grad_qv(pcmobj, mo_coeff, mo_occ, atmlst=None, verbose=None): ''' + dv_solv / da slow version with finite difference ''' log = logger.new_logger(pcmobj, verbose) t1 = log.init_timer() - mol = pcmobj.mol.copy() + pmol = pcmobj.mol.copy() + mol = pmol.copy() if atmlst is None: atmlst = range(mol.natm) nao, nmo = mo_coeff.shape @@ -127,7 +143,7 @@ def pcm_vmat_scanner(mol): return v vmat = cupy.zeros([len(atmlst), 3, nao, nocc]) - eps = 1e-4 + eps = 1e-3 for i0, ia in enumerate(atmlst): for ix in range(3): dv = numpy.zeros_like(coords) @@ -145,6 +161,7 @@ def pcm_vmat_scanner(mol): grad_vmat = contract("iq,ip->pq", grad_vmat, mo_coeff) vmat[i0,ix] = grad_vmat t1 = log.timer_debug1('computing solvent grad veff', *t1) + pcmobj.reset(pmol) return vmat def make_hess_object(hess_method): diff --git a/gpu4pyscf/solvent/pcm.py b/gpu4pyscf/solvent/pcm.py index a310b1bd..e9634628 100644 --- a/gpu4pyscf/solvent/pcm.py +++ b/gpu4pyscf/solvent/pcm.py @@ -124,7 +124,7 @@ def gen_surface(mol, ng=302, vdw_scale=1.2): riJ = cupy.sum((atom_grid[:,None,:] - atom_coords[None,:,:])**2, axis=2)**0.5 diJ = (riJ - R_in_J) / R_sw_J diJ[:,ia] = 1.0 - diJ[diJ<1e-8] = 0.0 + diJ[diJ<1e-12] = 0.0 fiJ = switch_h(diJ) @@ -365,6 +365,7 @@ def Hessian(self, hess_method): def reset(self, mol=None): self.surface = None + self.intopt = None super().reset(mol) return self diff --git a/gpu4pyscf/solvent/tests/test_pcm_hessian.py b/gpu4pyscf/solvent/tests/test_pcm_hessian.py index 36e85c75..ff6074d4 100644 --- a/gpu4pyscf/solvent/tests/test_pcm_hessian.py +++ b/gpu4pyscf/solvent/tests/test_pcm_hessian.py @@ -28,111 +28,73 @@ def setUpModule(): H -0.7570000000 -0.0000000000 -0.4696000000 H 0.7570000000 0.0000000000 -0.4696000000 ''' - mol.basis = 'sto3g' + mol.basis = 'def2-tzvpp' mol.output = '/dev/null' mol.build(verbose=0) - epsilon = 35.9 + epsilon = 78.3553 lebedev_order = 29 - eps = 1e-4 + eps = 1e-3 xc = 'B3LYP' - tol = 1e-4 + tol = 1e-3 def tearDownModule(): global mol mol.stdout.close() del mol +def _check_hessian(method='C-PCM', ix=0, iy=0): + pmol = mol.copy() + pmol.build() + + mf = dft.rks.RKS(pmol, xc=xc).density_fit().PCM() + mf.with_solvent.method = method + mf.with_solvent.eps = epsilon + mf.with_solvent.lebedev_order = lebedev_order + mf.conv_tol = 1e-12 + mf.grids.atom_grid = (99,590) + mf.verbose = 0 + mf.kernel() + + g = mf.nuc_grad_method() + g.auxbasis_response = True + g.kernel() + g_scanner = g.as_scanner() + + coords = pmol.atom_coords() + v = np.zeros_like(coords) + v[ix,iy] = eps + pmol.set_geom_(coords + v, unit='Bohr') + pmol.build() + _, g0 = g_scanner(pmol) + + pmol.set_geom_(coords - v, unit='Bohr') + pmol.build() + _, g1 = g_scanner(pmol) + + h_fd = (g0 - g1)/2.0/eps + pmol.set_geom_(coords, unit='Bohr') + pmol.build() + + hobj = mf.Hessian() + hobj.set(auxbasis_response=2) + h = hobj.kernel() + + print(f"analytical Hessian H({ix},{iy})") + print(h[ix,:,iy,:]) + print(f"finite different Hessian H({ix},{iy})") + print(h_fd) + print('Norm of diff', np.linalg.norm(h[ix,:,iy,:] - h_fd)) + assert(np.linalg.norm(h[ix,:,iy,:] - h_fd) < tol) + class KnownValues(unittest.TestCase): def test_hess_cpcm(self): - pmol = mol.copy() - pmol.build() - - mf = dft.rks.RKS(pmol, xc=xc).density_fit().PCM() - mf.with_solvent.eps = epsilon - mf.with_solvent.lebedev_order = lebedev_order - mf.conv_tol = 1e-12 - mf.grids.atom_grid = (99,590) - mf.verbose = 0 - mf.kernel() - - g = mf.nuc_grad_method() - g.auxbasis_response = True - g.kernel() - g_scanner = g.as_scanner() - - ix = 0 - iy = 1 - coords = pmol.atom_coords() - v = np.zeros_like(coords) - v[ix,iy] = eps - pmol.set_geom_(coords + v, unit='Bohr') - pmol.build() - _, g0 = g_scanner(pmol) - - pmol.set_geom_(coords - v, unit='Bohr') - pmol.build() - _, g1 = g_scanner(pmol) - - h_fd = (g0 - g1)/2.0/eps - pmol.set_geom_(coords, unit='Bohr') - pmol.build() - - hobj = mf.Hessian() - hobj.set(auxbasis_response=2) - h = hobj.kernel() - - print(f"analytical Hessian H({ix},{iy})") - print(h[ix,:,iy,:]) - print(f"finite different Hessian H({ix},{iy})") - print(h_fd) - print('Norm of diff', np.linalg.norm(h[ix,:,iy,:] - h_fd)) - assert(np.linalg.norm(h[ix,:,iy,:] - h_fd) < tol) + _check_hessian(method='C-PCM', ix=0, iy=0) + _check_hessian(method='C-PCM', ix=0, iy=1) def test_hess_iefpcm(self): - pmol = mol.copy() - pmol.build() - - mf = dft.rks.RKS(pmol, xc=xc).density_fit().PCM() - mf.with_solvent.method = 'IEF-PCM' - mf.with_solvent.eps = epsilon - mf.with_solvent.lebedev_order = lebedev_order - mf.conv_tol = 1e-12 - mf.grids.atom_grid = (99,590) - mf.verbose = 0 - mf.kernel() - - g = mf.nuc_grad_method() - g.auxbasis_response = True - g.kernel() - g_scanner = g.as_scanner() - - ix = 0 - iy = 1 - coords = pmol.atom_coords() - v = np.zeros_like(coords) - v[ix,iy] = eps - pmol.set_geom_(coords + v, unit='Bohr') - pmol.build() - _, g0 = g_scanner(pmol) - - pmol.set_geom_(coords - v, unit='Bohr') - pmol.build() - _, g1 = g_scanner(pmol) - - h_fd = (g0 - g1)/2.0/eps - pmol.set_geom_(coords, unit='Bohr') - pmol.build() - - hobj = mf.Hessian() - hobj.set(auxbasis_response=2) - h = hobj.kernel() + _check_hessian(method='IEF-PCM', ix=0, iy=0) + _check_hessian(method='IEF-PCM', ix=0, iy=1) - print(f"analytical Hessian H({ix},{iy})") - print(h[ix,:,iy,:]) - print(f"finite different Hessian H({ix},{iy})") - print(h_fd) - print('Norm of diff', np.linalg.norm(h[ix,:,iy,:] - h_fd)) - assert(np.linalg.norm(h[ix,:,iy,:] - h_fd) < tol) if __name__ == "__main__": print("Full Tests for Hessian of PCMs") unittest.main() \ No newline at end of file