Merge branch 'master' into small-batch-grids

.github/workflows/unittest.yml

@@ -23,7 +23,7 @@ jobs:
       run: |
         pip3 config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simple
         python3 -m pip install --upgrade pip
-        pip3 install flake8 pytest coverage gpu4pyscf-libxc-cuda11x pytest-cov pyscf-dispersion
+        pip3 install flake8 pytest coverage pytest-cov pyscf-dispersion
         pip3 install pyscf --upgrade
         git config --global core.compression 9
     - name: Build GPU4PySCF

gpu4pyscf/drivers/dft_driver.py

@@ -25,6 +25,7 @@
 
 from pyscf import lib, gto
 from pyscf import dft, scf
+from pyscf.hessian import thermo
 from pyscf.lib import logger
 
 def warmup(atom=None):
@@ -88,6 +89,7 @@ def run_dft(mol_name, config, charge=None, spin=0):
     with_solvent   = config.get('with_solvent',   False)
     with_grad      = config.get('with_grad',      True)
     with_hess      = config.get('with_hess',      True)
+    with_thermo    = config.get('with_thermo',    False) 
     save_density   = config.get('save_density',   False)
     input_dir      = config.get('input_dir',      './')
 
@@ -133,11 +135,17 @@ def run_dft(mol_name, config, charge=None, spin=0):
         mf = mf.to_gpu()
 
     mf.chkfile = None
-    if with_solvent:
+    if with_solvent and config['solvent']['method'].endswith(('PCM', 'pcm')):
         mf = mf.PCM()
         mf.with_solvent.lebedev_order = 29
         mf.with_solvent.method = config['solvent']['method'].replace('PCM','-PCM')
         mf.with_solvent.eps = config['solvent']['eps']
+
+    if with_solvent and config['solvent']['method'].endswith(('smd', 'SMD')):
+        mf = mf.SMD()
+        mf.with_solvent.lebedev_order = 29
+        mf.with_solvent.method = 'SMD'
+        mf.with_solvent.solvent = config['solvent']['solvent']
 
     mf.direct_scf_tol = direct_scf_tol
     mf.chkfile = None
@@ -150,25 +158,27 @@ def run_dft(mol_name, config, charge=None, spin=0):
     scf_time = time.time() - start_time
     print(f'compute time for energy: {scf_time:.3f} s')
 
-    e1     = mf.scf_summary.get('e1',         0.0)
-    e_coul = mf.scf_summary.get('coul',       0.0)
-    e_xc   = mf.scf_summary.get('exc',        0.0)
-    e_disp = mf.scf_summary.get('dispersion', 0.0)
+    e1        = mf.scf_summary.get('e1',         0.0)
+    e_coul    = mf.scf_summary.get('coul',       0.0)
+    e_xc      = mf.scf_summary.get('exc',        0.0)
+    e_disp    = mf.scf_summary.get('dispersion', 0.0)
+    e_solvent = mf.scf_summary.get('e_solvent',  0.0)
 
     data_file = mol_name[:-4] + '_pyscf.h5'
     import h5py
     h5f = h5py.File(f'{local_dir}/{data_file}', 'w')
-    h5f.create_dataset('e_tot',    data=e_tot)
-    h5f.create_dataset('e1',       data=e1)
-    h5f.create_dataset('e_coul',   data=e_coul)
-    h5f.create_dataset('e_xc',     data=e_xc)
-    h5f.create_dataset('e_disp',   data=e_disp)
-    h5f.create_dataset('scf_time', data=scf_time)
+    h5f.create_dataset('e_tot',     data=e_tot)
+    h5f.create_dataset('e1',        data=e1)
+    h5f.create_dataset('e_coul',    data=e_coul)
+    h5f.create_dataset('e_xc',      data=e_xc)
+    h5f.create_dataset('e_disp',    data=e_disp)
+    h5f.create_dataset('e_solvent', data=e_solvent)
+    h5f.create_dataset('scf_time',  data=scf_time)
 
     dm = mf.make_rdm1()
     if isinstance(dm, cupy.ndarray): dm = dm.get()
     h5f.create_dataset('dm',       data=dm)
-    
+
     if save_density and xc.lower() != 'hf':
         weights = mf.grids.weights
         coords = mf.grids.coords
@@ -228,9 +238,24 @@ def run_dft(mol_name, config, charge=None, spin=0):
             start_time = time.time()
             h = mf.Hessian()
             h.auxbasis_response = 2
-            h_dft = h.kernel().transpose([0,2,1,3]).reshape([3*natm, 3*natm])
+            _h_dft = h.kernel()
+            h_dft = _h_dft.transpose([0,2,1,3]).reshape([3*natm, 3*natm])
             hess_time = time.time() - start_time
             print(f'compute time for hessian: {hess_time:.3f} s')
+
+            if with_thermo:
+                # harmonic analysis
+                start_time = time.time()
+                normal_mode = thermo.harmonic_analysis(mol, _h_dft)
+
+                thermo_dat = thermo.thermo(
+                    mf,                            # GPU4PySCF object
+                    normal_mode['freq_au'],
+                    298.15,                            # room temperature
+                    101325)                            # standard atmosphere
+                thermo_time = time.time() - start_time
+                print(f'compute time for harmonic analysis: {thermo_time:.3f} s')
+
         except Exception as exc:
             print(traceback.format_exc())
             print(exc)
@@ -240,6 +265,26 @@ def run_dft(mol_name, config, charge=None, spin=0):
         h5f.create_dataset('hess', data=h_dft)
         h5f.create_dataset('hess_time', data=hess_time)
 
+        if with_thermo: 
+            h5f.create_dataset('freq_au',         data=normal_mode['freq_au'])
+            h5f.create_dataset('freq_wavenumber', data=normal_mode['freq_wavenumber'])
+            h5f.create_dataset('E_tot',           data=thermo_dat['E_tot'][0])
+            h5f.create_dataset('H_tot',           data=thermo_dat['H_tot'][0])
+            h5f.create_dataset('G_tot',           data=thermo_dat['G_tot'][0])
+            h5f.create_dataset('E_elec',          data=thermo_dat['E_elec'][0])
+            h5f.create_dataset('E_trans',         data=thermo_dat['E_trans'][0])
+            h5f.create_dataset('E_rot',           data=thermo_dat['E_rot'][0])
+            h5f.create_dataset('E_vib',           data=thermo_dat['E_vib'][0])
+            h5f.create_dataset('E_0K',            data=thermo_dat['E_0K'][0])
+            h5f.create_dataset('H_elec',          data=thermo_dat['H_elec'][0])
+            h5f.create_dataset('H_trans',         data=thermo_dat['H_trans'][0])
+            h5f.create_dataset('H_rot',           data=thermo_dat['H_rot'][0])
+            h5f.create_dataset('H_vib',           data=thermo_dat['H_vib'][0])
+            h5f.create_dataset('G_elec',          data=thermo_dat['G_elec'][0])
+            h5f.create_dataset('G_trans',         data=thermo_dat['G_trans'][0])
+            h5f.create_dataset('G_rot',           data=thermo_dat['G_rot'][0])
+            h5f.create_dataset('G_vib',           data=thermo_dat['G_vib'][0])
+
     h5f.close()
 
     # copy the files to destination folder

gpu4pyscf/drivers/dft_sample.json

@@ -12,11 +12,13 @@
     ],
     "basis": "def2-tzvpp",
     "xc": "b3lyp",
-    "verbose": 6,
-    "with_solvent": false,
+    "verbose": 4,
+    "with_solvent": true,
+    "with_thermo": true,
     "solvent": {
         "eps": 78.3553,
-        "method": "CPCM"
+        "solvent": "water",
+        "method": "SMD"
     },
     "with_gpu": true,
     "with_df": true,

gpu4pyscf/solvent/_attach_solvent.py

@@ -104,7 +104,7 @@ def energy_elec(self, dm=None, h1e=None, vhf=None):
         if isinstance(e_solvent, cupy.ndarray):
             e_solvent = e_solvent.get()[()]
         e_tot += e_solvent
-        self.scf_summary['e_solvent'] = vhf.e_solvent.real
+        self.scf_summary['e_solvent'] = e_solvent
 
         if (hasattr(self.with_solvent, 'method') and self.with_solvent.method.upper() == 'SMD'):
             if self.with_solvent.e_cds is None: