A new LibXC interface and new MC-DCFT code (#70)

* new libxc interface * add cffi dependency * fix cffi dependency * bug fix * improve libxc interface - migrate all ctypes interfaces to cffi * bug fix * optimize performance for libxc interface * libxc interface enhancement * use new libxc interface for mcdcft * new mcdcft implementation * clean up * clean up * use Python 3.8 for automated test - trexio no longer publish wheels for cp37, which makes HDF5 backend missing for trexio * improve mcdcft __init__ * clean up and performance enhancement * add mcdcft preset * clean up, improve comments, and fix examples * fix edge case: natorb from molden file may not be square matrix * improve example and test; enhance chkfile handling * update libxc interface according to PySCF 2.7 (commit 410d960) * Revert "use Python 3.8 for automated test - trexio no longer publish wheels for cp37, which makes HDF5 backend missing for trexio" This reverts commit fe54aa1. * add comments to files adapted from PySCF core modules * implement LibXC interface using CFFI ABI mode * add example to set MC-DCFT functional parameters * new LibXC interface based on ctypes * resolve conflict * bug fix * clean up and enhance comments * enhance comment * minor enhancement
pyscf · Dec 4, 2024 · 7f130a6 · 7f130a6
1 parent 880b3d1
commit 7f130a6
Show file tree

Hide file tree

Showing 16 changed files with 3,550 additions and 407 deletions.
diff --git a/examples/dft2/24-reparameterize_xc_functional.py b/examples/dft2/24-reparameterize_xc_functional.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+
+'''
+Running DFT calculations with reparameterized functionals in Libxc.
+
+See also
+* Example 24-custom_xc_functional.py to customize XC functionals using the
+  functionals provided by Libxc library.
+'''
+
+from pyscf import gto
+from pyscf import dft, dft2
+import numpy as np
+
+# Patch the libxc module with the new implementation
+dft.libxc = dft2.libxc
+dft.numint.libxc = dft2.libxc
+dft.numint.LibXCMixin.libxc = dft2.libxc
+
+XC_ID_B97_2 = 410
+
+param_b97_1 = np.array([0.789518, 0.573805, 0.660975, 0.0, 0.0,
+                        0.0820011, 2.71681, -2.87103, 0.0, 0.0,
+                        0.955689, 0.788552, -5.47869, 0.0, 0.0,
+                        0.21
+])
+
+param_b97_3 = np.array([0.7334648, 0.292527, 3.338789, -10.51158, 10.60907,
+                        0.5623649, -1.32298, 6.359191, -7.464002, 1.827082,
+                        1.13383, -2.811967, 7.431302, -1.969342, -11.74423,
+                        2.692880E-01
+])
+
+mol = gto.M(
+    atom = '''
+    O  0.   0.       0.
+    H  0.   -0.757   0.587
+    H  0.   0.757    0.587 ''',
+    basis = 'ccpvdz')
+
+# Run normal B97-1
+print('Normal B97-1')
+mf = dft.RKS(mol, 'B97-1')
+e_b971 = mf.kernel()
+
+# Run normal B97-2
+print('\nNormal B97-2')
+mf.xc = 'B97-2'
+e_b972 = mf.kernel()
+
+# Run normal B97-3
+print('\nNormal B97-3')
+mf.xc = 'B97-3'
+e_b973 = mf.kernel()
+
+# Construct XC named `myfunctional` based on B97-2, but set its parameter to be B97-1
+print('\nReparameterized B97-2: will be the same as B97-1')
+dft2.libxc.register_custom_functional_('myfunctional', 'B97-2',
+                                       ext_params={XC_ID_B97_2: param_b97_1})
+mf.xc = 'myfunctional'
+e = mf.kernel()
+print('difference:', e - e_b971)
+
+# Update the parameter of `myfunctional` to be the same as B97-3 and rerun
+print('\nReparameterized B97-2: will be the same as B97-3')
+dft2.libxc.update_custom_functional_('myfunctional',
+                                     ext_params={XC_ID_B97_2: param_b97_3})
+e = mf.kernel()
+print('difference:', e - e_b973)
+print()
+
diff --git a/examples/mcdcft/01_h2_potential_curve.py b/examples/mcdcft/01_h2_potential_curve.py
@@ -1,22 +1,21 @@
-import pyscf
-from pyscf import scf, gto, mcscf, lib
+from pyscf import scf, gto, lib, mcscf, mcdcft
 import numpy as np
-#from pyscf.fci import csf_solver
-from pyscf.mcdcft import mcdcft
 
 '''
 This input file performs a potential energy scan of H2 using cBLYP (density
 coherence functional without reparameterization).
 '''
 
+preset = dict(args=dict(f=mcdcft.dcfnal.f_v1, negative_rho=True), xc_code='BLYP')
+mcdcft.dcfnal.register_dcfnal_('cBLYP', preset)
+
 def run(r, chkfile=None, mo_coeff=None):
-    r /= 2
-    mol = gto.M(atom=f'H  0 0 {r}; H 0 0 -{r}', basis='cc-pvtz', 
+    r *= 0.5
+    mol = gto.M(atom=f'H  0 0 {r}; H 0 0 -{r}', basis='cc-pvtz',
           symmetry=False, verbose=3, unit='angstrom')
     mf = scf.RHF(mol)
     mf.kernel()
-    mc = mcdcft.CASSCF(mf, 'BLYP', 2, 2, ot_name='cBLYP', 
-                       grids_level=6)
+    mc = mcdcft.CASSCF(mf, 'cBLYP', 2, 2, grids_level=6)
     #mc.fcisolver = csf_solver(mol, smult=1)
     mc.fix_spin_(ss=0)
     if mo_coeff is not None:
@@ -39,7 +38,6 @@ def scan():
         if mc and mc.converged:
             mo_coeff = mc.mo_coeff
             e_tot = mc.e_tot
-            mc.dump_mcdcft_chk(chkfile)
         else:
             print(f"Warning: bond distance {r:02.2f} not converged")
 

diff --git a/examples/mcdcft/02_h2_potential_curve_restart.py b/examples/mcdcft/02_h2_potential_curve_restart.py
@@ -1,34 +1,32 @@
-import pyscf
-from pyscf import scf, gto, mcscf, lib
+from pyscf import scf, gto, mcscf, lib, mcdcft
 import numpy as np
-#from pyscf.fci import csf_solver
-from pyscf.mcdcft import mcdcft
 
 '''
-This input file performs a potential energy scan of H2 using cBLYP-2 (density coherence functional with reparameterization). It will read the checkpoint files generated by example 01.
+This input file performs a potential energy scan of H2 using cBLYP-2 (density coherence functional
+with reparameterization). It will read the checkpoint files generated by example 01.
 '''
 
+dc_preset = dict(args=dict(f=mcdcft.dcfnal.f_v1, negative_rho=True), xc_code='0.92 * B88 , 1.59 * LYP')
+mcdcft.dcfnal.register_dcfnal_('cBLYP-2', dc_preset)
+
 def run(chkfile):
     mol = lib.chkfile.load_mol(chkfile)
     mol.verbose = 3
-    mf = scf.RHF(mol)
-    mc = mcdcft.CASSCF(mf, None, 2, 2, grids_level=6)
+    mc = mcdcft.DCFT_base(mol, grids_level=6)
     mc.load_mcdcft_chk(chkfile)
-    mc.recalculate_with_xc('0.92 * B88 , 1.59 * LYP', 
-                            ot_name='cBLYP-2', dump_chk=chkfile)
+    mc.recalculate_with_xc('cBLYP-2', dump_chk=chkfile)
     return mc
 
 
 def scan():
-    mo_coeff = None
     method = 'MC-DCFT'
 
     rrange = np.arange(0.5, 6.1, 0.1)
 
     for r in rrange:
         chkfile = f'{method}_{r:02.2f}.chk'
-        mc = run(chkfile)
-        
+        run(chkfile)
+
 
 if __name__ == '__main__':
     lib.num_threads(2)

diff --git a/examples/mcdcft/03_DC24_single_point.py b/examples/mcdcft/03_DC24_single_point.py
@@ -0,0 +1,14 @@
+from pyscf import scf, gto, mcdcft
+
+'''
+This input file performs a single point energy calculation of H2 with DC24.
+'''
+
+mol = gto.M(atom='H 0 0 0; H 0 0 0.74', basis='def2-tzvp',
+      symmetry=False, verbose=3, unit='angstrom')
+mf = scf.RHF(mol)
+mf.kernel()
+mc = mcdcft.CASSCF(mf, 'DC24', 2, 2, grids_level=(99, 590))
+mc.chkfile = 'H2_DC24.chk'
+mc.kernel()
+
diff --git a/examples/mcdcft/04_functional_parameter.py b/examples/mcdcft/04_functional_parameter.py
@@ -0,0 +1,41 @@
+from pyscf import scf, gto, mcdcft
+
+'''
+This input file performs a single point energy calculation of H2 with DC24.
+It demonstrates how users can define custom density coherence functionals,
+especially how functional parameters can be set.
+'''
+
+mol = gto.M(atom='H 0 0 0; H 0 0 0.74', basis='def2-tzvp',
+      symmetry=False, verbose=3, unit='angstrom')
+mf = scf.RHF(mol)
+mf.kernel()
+
+# Define a preset for the DC24 functional, which contains functional parameters.
+preset = {
+    # Display name of the new functional. It will be used in displayed messages and chkfiles.
+    'display_name': 'DC24,custom',
+    # The Kohn-Sham exchange correlation functional the DC functional is based on.
+    # Use the same notation as used in the KS module
+    'xc_code': 'HCTH407',
+    # Ratio of the MCSCF exchange-correlation energy
+    'hyb_x': 4.525671e-01,
+    # Functional parameters. We will set the functional parameters of HCTH407 (LibXC
+    #     functional 164) with the values listed here. The meaning of each parameter
+    #     can be obtained from the `xc-info` utility provided by LibXC.
+    'params': {164: [8.198942e-01, 4.106753e+00, -3.716774e+01, 1.100812e+02, -9.600026e+01,
+                     1.352989e+01, -6.881959e+01, 2.371350e+02, -3.433615e+02, 1.720927e+02,
+                     1.134169e+00, 1.148509e+01, -2.210990e+01, -1.006682e+02, 1.477906e+02]},
+}
+
+# Register the functional preset to the dcfnal module with identifier 'DC24_custom'
+mcdcft.dcfnal.register_dcfnal_('DC24_custom', preset)
+
+# Evaluate MC-DCFT energy of the custom functional
+mc = mcdcft.CASSCF(mf, 'DC24_custom', 2, 2, grids_level=(99, 590))
+e1 = mc.kernel()[0]
+
+# Compare with the results from pre-defined DC24
+e2 = mc.recalculate_with_xc('DC24')[0]
+print("Energy difference:", e1 - e2)
+
diff --git a/pyscf/dft2/__init__.py b/pyscf/dft2/__init__.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+# Copyright 2014-2022 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+from pyscf.dft2 import libxc