v0.6.6

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):

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')

gpu4pyscf/__init__.py

@@ -1,2 +1,2 @@
 from . import lib, grad, hessian, solvent, scf, dft
-__version__ = '0.6.5'
+__version__ = '0.6.6'

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

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):

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
 

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()