update examples

examples/00-h2o.py

@@ -36,12 +36,12 @@
     atom=atom,                         # water molecule
     basis='def2-tzvpp',                # basis set
     output='./pyscf.log',              # save log file
-    verbose=6                         # control the level of print info
+    verbose=6                          # control the level of print info
     )
 
 mf_GPU = rks.RKS(                      # restricted Kohn-Sham DFT
     mol,                               # pyscf.gto.object
-    xc='b3lyp'                        # xc funtionals, such as pbe0, wb97m-v, tpss,
+    xc='b3lyp'                         # xc funtionals, such as pbe0, wb97m-v, tpss,
     ).density_fit()                    # density fitting
 
 mf_GPU.grids.atom_grid = (99,590)      # (99,590) lebedev grids, (75,302) is often enough
@@ -51,7 +51,7 @@
 
 # Compute Energy
 e_dft = mf_GPU.kernel()
-print(f"total energy = {e_dft}") # -76.26736519501688
+print(f"total energy = {e_dft}")       # -76.46668196729536
 
 # Compute Gradient
 g = mf_GPU.nuc_grad_method()

examples/02-h2o_geomopt.py

@@ -43,4 +43,4 @@ def callback(envs):
 mol_eq = optimize(mf_GPU, maxsteps=20, callback=callback)
 print("Optimized coordinate:")
 print(mol_eq.atom_coords())
-print('geometry optimization took', time.time() - start_time, 's')
+print('Geometry optimization took', time.time() - start_time, 's')

examples/04-h2o_esp.py

@@ -21,6 +21,7 @@
 import numpy as np
 from pyscf import gto
 from gpu4pyscf.dft import rks
+from gpu4pyscf.gto.int3c1e import int1e_grids
 
 atom ='''
 O       0.0000000000    -0.0000000000     0.1174000000
@@ -33,10 +34,8 @@
 mf.kernel()
 dm = mf.make_rdm1()  # compute one-electron density matrix
 
-# Use default mesh grids
-coords = mf.grids.coords.get()
+# Use default Lebedev grids
+coords = mf.grids.coords
 
-# The efficiency can be improved if needed
-from pyscf import df
-fakemol = gto.fakemol_for_charges(coords)
-v = np.einsum('ijp,ij->p', df.incore.aux_e2(mol, fakemol), dm)
+# Calculate electrostatic potential
+v = int1e_grids(mol, coords, dm=dm) # performing 'ijp,ij->p' efficiently

examples/05-h2o_multipole_moment.py

@@ -32,10 +32,10 @@
 mf.kernel()
 dm = mf.make_rdm1()
 
-dip = mf.dip_moment(unit='DEBYE', dm=dm.get())
+dip = mf.dip_moment(unit='DEBYE', dm=dm)
 print('dipole moment:')
 print(dip)
 
-quad = mf.quad_moment(unit='DEBYE-ANG', dm=dm.get())
+quad = mf.quad_moment(unit='DEBYE-ANG', dm=dm)
 print('quadrupole moment:')
 print(quad)

examples/14-pcm_solvent.py

@@ -31,9 +31,9 @@
 mf = rks.RKS(mol, xc='HYB_GGA_XC_B3LYP').density_fit()
 mf = mf.PCM()
 mf.grids.atom_grid = (99,590)
-mf.with_solvent.lebedev_order = 29 # 302 Lebedev grids
-mf.with_solvent.method = 'IEF-PCM'
-mf.with_solvent.eps = 78.3553
+mf.with_solvent.lebedev_order = 29  # 302 Lebedev grids
+mf.with_solvent.method = 'IEF-PCM'   # Can be C-PCM, SS(V)PE, COSMO
+mf.with_solvent.eps = 78.3553        # Dielectric constant
 mf.kernel()
 
 gradobj = mf.nuc_grad_method()

examples/15-chelpg.py

@@ -32,18 +32,18 @@
 mol.basis = '631g'
 mol.unit = 'B'
 mol.build()
-mol.verbose = 6
+mol.verbose = 4
 
 xc = 'b3lyp'
 mf = rks.RKS(mol, xc=xc)
 mf.grids.level = 5
 mf.kernel()
 q = chelpg.eval_chelpg_layer_gpu(mf)
-print('partial charge with CHELPG, using modified Bondi radii')
+print('Partial charge with CHELPG, using modified Bondi radii')
 print(q) # [ 0.04402311  0.11333945 -0.25767919  0.10031663]
 
 # Customize the radii used for calculating CHELPG charges
 from pyscf.data import radii
 q = chelpg.eval_chelpg_layer_gpu(mf, Rvdw=radii.UFF)
-print('partial charge with CHELPG, using UFF radii')
+print('Partial charge with CHELPG, using UFF radii')
 print(q)

examples/16-smd_solvent.py

@@ -28,16 +28,14 @@
 
 mol = pyscf.M(atom=atom, basis='def2-tzvpp', verbose=1)
 mf = dft.rks.RKS(mol, xc='HYB_GGA_XC_B3LYP').density_fit()
-mf = mf.SMD()
 mf.grids.atom_grid = (99,590)
-mf.with_solvent.lebedev_order = 29 # 302 Lebedev grids
-mf.with_solvent.method = 'SMD'
-mf.with_solvent.solvent = 'water'
-e_tot = mf.kernel()
-print('total energy with SMD:', e_tot)
+e_gas = mf.kernel()
+print('total energy in gas phase:', e_gas)
 
-gradobj = mf.nuc_grad_method()
-f = gradobj.kernel()
+mf = mf.SMD()   # Add SMD model to the mean-field object
+mf.with_solvent.lebedev_order = 29 # 302 Lebedev grids,
+mf.with_solvent.solvent = 'water' # Has to be a string, lookup the solvent name from https://comp.chem.umn.edu/solvation/mnsddb.pdf
+e_smd = mf.kernel()
+print('total energy in water:', e_smd)
 
-hessobj = mf.Hessian()
-h = hessobj.kernel()
+print('Solvation free energy:', e_smd - e_gas)

examples/19-unrestricted_dft.py

@@ -49,14 +49,3 @@
 
 hobj_with_pcm = mf_with_pcm.Hessian()
 h = hobj_with_pcm.kernel()
-
-# SCF, gradient, and Hessian for DF-UKS with IEF-PCM
-mf_with_smd = mf.SMD()
-mf_with_smd.with_solvent.solvent = 'water'
-mf_with_smd.kernel()
-
-gobj_with_smd = mf_with_smd.nuc_grad_method()
-g = gobj_with_smd.kernel()
-
-hobj_with_smd = mf_with_smd.Hessian()
-h = hobj_with_smd.kernel()

examples/20-dfmp2.py

@@ -35,7 +35,18 @@
 e_corr, t2 = ptobj.kernel()
 e_mp2 = e_hf + e_corr
 
+# It prints out MP2 energies, those energies are assessible in the PT object.
+print('MP2 correlation energy:', ptobj.emp2)
+print('SCS MP2 correlation energy:', ptobj.emp2_scs)
+print('Total energy with SCS MP2:', ptobj.e_tot_scs)
+
+print('----- frozen core --------')
+
 # frozen core
 ptobj.frozen = [0]
 e_corr, t2 = ptobj.kernel()
 e_mp2 = e_hf + e_corr
+
+print('MP2 correlation energy:', ptobj.emp2)
+print('SCS MP2 correlation energy:', ptobj.emp2_scs)
+print('Total energy with SCS MP2:', ptobj.e_tot_scs)

examples/22-resp_charge.py

@@ -42,11 +42,11 @@
 print(q0)
 
 # RESP charge // first stage fitting
-q1 = esp.resp_solve(mol, dm)    
+q1 = esp.resp_solve(mol, dm)
 
-# Add constraint: fix those charges in the second stage 
+# Add constraint: fix those charges in the second stage
 # q2[4] = q1[4]
-# q2[5] = q1[5] 
+# q2[5] = q1[5]
 # q2[6] = q1[6]
 # q2[7] = q1[7]
 sum_constraints = []
@@ -58,7 +58,7 @@
 equal_constraints = [[1,2,3]]
 
 # RESP charge // second stage fitting
-q2 = esp.resp_solve(mol, dm, resp_a=1e-3, 
+q2 = esp.resp_solve(mol, dm, resp_a=1e-3,
                     sum_constraints=sum_constraints,
                     equal_constraints=equal_constraints)
 print('Fitted RESP charge')

examples/24-cp_bsse.py → examples/25-cp_bsse.py

@@ -21,15 +21,15 @@
 from gpu4pyscf.dft import rks
 
 atom_A = [
-('O', (0.000000, 0.000000, 0.000000)),
-('H', (0.000000, 0.757160, 0.586260)),
-('H', (0.000000, -0.757160, 0.586260))
+    ('O', (0.000000, 0.000000, 0.000000)),
+    ('H', (0.000000, 0.757160, 0.586260)),
+    ('H', (0.000000, -0.757160, 0.586260))
 ]
 
 atom_B = [
-('O', (0.000000, 0.000000, 2.913530)),
-('H', (0.000000, 0.757160, 3.499790)),
-('H', (0.000000, -0.757160, 3.499790))
+    ('O', (0.000000, 0.000000, 2.913530)),
+    ('H', (0.000000, 0.757160, 3.499790)),
+    ('H', (0.000000, -0.757160, 3.499790))
 ]
 
 atom_AB = atom_A + atom_B
@@ -51,7 +51,7 @@
 mol_B_ghost.build()
 
 def solve_dft(mol, xc='b3lyp'):
-    mf = rks.RKS(mol, xc='b3lyp').density_fit()
+    mf = rks.RKS(mol, xc=xc).density_fit()
     mf.grids.atom_grid = (99,590)
     return mf.kernel()