remove comment

examples/dft_driver.py

@@ -25,7 +25,7 @@
 parser.add_argument("--basis",    type=str,  default='def2-tzvpp')
 parser.add_argument("--auxbasis", type=str,  default='def2-tzvpp-jkfit')
 parser.add_argument("--xc",       type=str,  default='B3LYP')
-parser.add_argument("--solvent",  type=bool, default=False)
+parser.add_argument("--solvent",  type=str, default='')
 args = parser.parse_args()
 
 start_time = time.time()
@@ -44,11 +44,11 @@
 if args.solvent:
     mf_df = mf_df.PCM()
     mf_df.lebedev_order = 29
-    mf_df.method = 'IEF-PCM'
+    mf_df.method = args.solvent
 mf_df.grids.atom_grid = (99,590)
 mf_df.direct_scf_tol = 1e-14
 mf_df.direct_scf = 1e-14
-mf_df.kernel()
+e_tot = mf_df.kernel()
 scf_time = time.time() - start_time
 print(f'compute time for energy: {scf_time:.3f} s')
 
@@ -65,3 +65,6 @@
 h_dft = h.kernel()
 hess_time = time.time() - start_time
 print(f'compute time for hessian: {hess_time:.3f} s')
+
+import numpy
+numpy.savez('gpu4pyscf_out.npz', e_tot=e_tot, f=f, h_dft=h_dft)

gpu4pyscf/solvent/grad/pcm.py

@@ -176,11 +176,8 @@ def grad_elec(pcmobj, dm):
     '''
     if not pcmobj._intermediates:
         pcmobj.build()
-    mol = pcmobj.mol
 
     gridslice    = pcmobj.surface['gslice_by_atom']
-    grid_coords  = pcmobj.surface['grid_coords']
-    exponents    = pcmobj.surface['charge_exp']
     v_grids      = pcmobj._intermediates['v_grids']
     A            = pcmobj._intermediates['A']
     D            = pcmobj._intermediates['D']
@@ -212,32 +209,6 @@ def grad_elec(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)
-
-    # nuclei response
-    int2c2e_ip1 = mol._add_suffix('int2c2e_ip1')
-    v_ng_ip1 = gto.mole.intor_cross(int2c2e_ip1, fakemol_nuc, fakemol)
-    v_ng_ip1 = cupy.asarray(v_ng_ip1)
-
-    dv_g = contract('g,xng->nx', q_sym, v_ng_ip1)
-    de -= contract('nx,n->nx', dv_g, atom_charges)
-
-    # nuclei potential response
-    int2c2e_ip2 = mol._add_suffix('int2c2e_ip2')
-    v_ng_ip2 = gto.mole.intor_cross(int2c2e_ip2, fakemol_nuc, fakemol)
-    v_ng_ip2 = cupy.asarray(v_ng_ip2)
-
-    dv_g = contract('n,xng->gx', atom_charges, v_ng_ip2)
-    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)

gpu4pyscf/solvent/hessian/pcm.py

@@ -49,11 +49,8 @@ def hess_nuc(pcmobj):
     # 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)
+    dv_g = numpy.einsum('n,xyng->ngxy', atom_charges, v_ng_ip1ip2)
+    dv_g = numpy.einsum('ngxy,g->ngxy', dv_g, q_sym)
 
     de = numpy.zeros([mol.natm, mol.natm, 3, 3])
     for ia in range(mol.natm):
@@ -64,13 +61,10 @@ def hess_nuc(pcmobj):
 
     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):
@@ -82,11 +76,14 @@ def hess_nuc(pcmobj):
 def hess_elec(pcmobj, dm, verbose=None):
     '''
     slow version with finite difference
+    TODO: use analytical hess_nuc
     '''
     log = logger.new_logger(pcmobj, verbose)
     t1 = log.init_timer()
-    mol = pcmobj.mol.copy()
+    pmol = pcmobj.mol.copy()
+    mol = pmol
     coords = mol.atom_coords(unit='Bohr')
+
     def pcm_grad_scanner(mol):
         pcmobj.reset(mol)
         e, v = pcmobj._get_vind(dm)
@@ -108,7 +105,6 @@ 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):