[WIP] MC-PDFT

examples/mcpdft/00-simple_mcpdft.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env/python
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+# 1. CASCI density
+
+mc0 = mcpdft.CASCI (mf, 'tPBE', 6, 8).run ()
+
+# 2. CASSCF density
+# Note that the MC-PDFT energy may not be lower, even though
+# E(CASSCF)<=E(CASCI).
+
+mc1 = mcpdft.CASSCF (mf, 'tPBE', 6, 8).run ()
+
+# 3. analyze () does the same thing as CASSCF analyze ()
+
+mc1.verbose = 4
+mc1.analyze ()
+
+# 4. Energy decomposition for additional analysis
+
+e_decomp = mc1.get_energy_decomposition (split_x_c=False)
+print ("e_nuc =",e_decomp[0])
+print ("e_1e =",e_decomp[1])
+print ("e_Coul =",e_decomp[2])
+print ("e_OT =",e_decomp[3])
+print ("e_ncwfn (not included in total energy) =",e_decomp[4])
+print ("e_PDFT - e_MCSCF =", mc1.e_tot - mc1.e_mcscf)
+print ("e_OT - e_ncwfn =", e_decomp[3] - e_decomp[4])
+e_decomp = mc1.get_energy_decomposition (split_x_c=True)
+print ("e_OT (x component) = ",e_decomp[3])
+print ("e_OT (c component) = ",e_decomp[4])
+
+

examples/mcpdft/01-different_functionals.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env/python
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+mc = mcpdft.CASCI (mf, 'tPBE', 6, 8).run ()
+
+# 1. Change the functional by setting mc.otxc
+
+mc.otxc = 'ftBLYP'
+mc.kernel ()
+
+# 2. Change the functional and compute energy in one line without
+# reoptimizing the wave function using compute_pdft_energy_
+
+mc.compute_pdft_energy_(otxc='ftPBE')
+
+#
+# The leading "t" or "ft" identifies either a "translated functional"
+# [JCTC 10, 3669 (2014)] or a "fully-translated functional"
+# [JCTC 11, 4077 (2015)]. It can be combined with a general PySCF
+# xc string containing any number of pure LDA or GGA functionals.
+# Meta-GGAs, built-in hybrid functionals ("B3LYP"), and range-separated
+# functionals are not supported.
+# 
+
+# 3. A translated user-defined compound functional
+
+mc.compute_pdft_energy_(otxc="t0.3*B88 + 0.7*SLATER,0.4*VWN5+0.6*LYP")
+
+# 4. A fully-translated functional consisting of "exchange" only
+
+mc.compute_pdft_energy_(otxc="ftPBE,")
+
+# 5. A fully-translated functional consisting of "correlation" only
+
+mc.compute_pdft_energy_(otxc="ft,PBE")
+
+# 6. The sum of 5 and 6 (look at the "Eot" output)
+
+mc.compute_pdft_energy_(otxc="ftPBE")
+
+

examples/mcpdft/02-hybrid_functionals.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env/python
+
+
+#
+# JPCL 11, 10158 (2020)
+# (but see #5 below)
+#
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+# 1. The only two predefined hybrid functionals as of writing
+
+mc = mcpdft.CASSCF (mf, 'tPBE0', 6, 8).run ()
+mc.compute_pdft_energy_(otxc='ftPBE0')
+
+# 2. Other predefined hybrid functionals are not supported
+
+try:
+    mc.compute_pdft_energy_(otxc='tB3LYP')
+except NotImplementedError as e:
+    print ("otxc='tB3LYP' results in NotImplementedError:")
+    print (str (e))
+
+# 3. Construct a custom hybrid functional by hand
+
+my_otxc = 't.8*B88 + .2*HF, .8*LYP + .2*HF'
+mc.compute_pdft_energy_(otxc=my_otxc)
+
+# 4. Construct the same custom functional using helper function
+
+my_otxc = 't' + mcpdft.hyb ('BLYP', .2)
+mc.compute_pdft_energy_(otxc=my_otxc)
+
+# 5. "lambda-MC-PDFT" of JCTC 16, 2274, 2020.
+
+my_otxc = 't' + mcpdft.hyb('PBE', .5, hyb_type='lambda')
+#       = 't0.5*PBE + 0.5*HF, 0.75*PBE + 0.5*HF'
+mc.compute_pdft_energy_(otxc=my_otxc)
+
+
+
+

examples/mcpdft/03-metaGGA_functionals.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env/python
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+# The translation of Meta-GGAs and hybrid-Meta-GGAs [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121]
+
+# Translated-Meta-GGA
+mc = mcpdft.CASCI(mf, 'tM06L', 6, 8).run ()
+
+# Hybrid-Translated-Meta-GGA
+tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average')
+mc = mcpdft.CASCI(mf, tM06L0, 6, 8).run ()
+
+# MC23: meta-hybrid on-top functional [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121]
+
+# State-Specific
+mc = mcpdft.CASCI(mf, 'MC23', 6, 8)
+mc.kernel()
+
+# State-average
+nroots=2
+mc = mcpdft.CASCI(mf, 'MC23', 6, 8)
+mc.fcisolver.nroots=nroots
+mc.kernel()[0]

examples/mcpdft/11-grid_scheme.py

@@ -0,0 +1,40 @@
+#!/usr/bin/env python
+
+from pyscf import gto, scf, dft
+from pyscf import mcpdft
+
+# See also pyscf/examples/dft/11-grid-scheme.py
+
+mol = gto.M(
+    verbose = 0,
+    atom = '''
+    o    0    0.       0.
+    h    0    -0.757   0.587
+    h    0    0.757    0.587''',
+    basis = '6-31g')
+mf = scf.RHF (mol).run ()
+mc = mcpdft.CASSCF(mf, 'tLDA', 4, 4).run()
+print('Default grid setup.  E = %.12f' % mc.e_tot)
+
+# See pyscf/dft/radi.py for more radial grid schemes
+mc.grids.radi_method = dft.mura_knowles
+print('radi_method = mura_knowles.  E = %.12f' % mc.kernel()[0])
+
+# All in one command:
+mc = mcpdft.CASSCF (mf, 'tLDA', 4, 4,
+                    grids_attr={'radi_method': dft.gauss_chebyshev})
+print('radi_method = gauss_chebyshev.  E = %.12f' % mc.kernel ()[0])
+
+# Or inline with an already-built mc object:
+e = mc.compute_pdft_energy_(grids_attr={'radi_method': dft.delley})[0]
+print('radi_method = delley.  E = %.12f' % e)
+
+# All grids attributes can be addressed in any of the ways above
+# There is also a shortcut to address grids.level:
+
+mc = mcpdft.CASSCF(mf, 'tLDA', 4, 4, grids_level=4).run()
+print('grids.level = 4.  E = %.12f' % mc.e_tot)
+
+e = mc.compute_pdft_energy_(grids_level=2)[0]
+print('grids.level = 2.  E = %.12f' % e)
+

examples/mcpdft/15-state_average.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+'''
+State average
+
+This works the same way as mcscf.state_average_, but you
+must use the method attribute (mc.state_average, mc.state_average_)
+instead of the function call.
+'''
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M(
+    atom = [
+        ['O', ( 0., 0.    , 0.   )],
+        ['H', ( 0., -0.757, 0.587)],
+        ['H', ( 0., 0.757 , 0.587)],],
+    basis = '6-31g',
+    symmetry = 1)
+
+mf = scf.RHF(mol)
+mf.kernel()
+
+mc = mcpdft.CASSCF(mf, 'tPBE', 4, 4).state_average_([.64,.36]).run (verbose=4)
+
+print ("Average MC-PDFT energy =", mc.e_tot)
+print ("E_PDFT-E_CASSCF for state 0 =", mc.e_states[0]-mc.e_mcscf[0])
+print ("E_PDFT-E_CASSCF for state 1 =", mc.e_states[1]-mc.e_mcscf[1])
+
+

examples/mcpdft/16-multi_state.py

@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+
+'''
+multi-state
+
+A "quasidegenerate" extension of MC-PDFT for states close in energy. Currently
+only "compressed multi-state" (CMS-) [JCTC 16, 7444 (2020)] and "extended multi-
+state" (XMS-) [Faraday Discuss 224, 348-372 (2020)] is supported.
+'''
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M(
+    atom = [
+        ['Li', ( 0., 0.    , 0.   )],
+        ['H', ( 0., 0., 3)]
+    ], basis = 'sto-3g',
+    symmetry = 0 # symmetry enforcement is not recommended for MS-PDFT
+    )
+
+mf = scf.RHF(mol)
+mf.kernel()
+
+mc = mcpdft.CASSCF(mf, 'tpbe', 2, 2)
+mc.fix_spin_(ss=0) # often necessary!
+mc_sa = mc.state_average ([.5, .5]).run ()
+
+# For CMS-PDFT
+mc_cms = mc.multi_state([.5, .5], "cms").run(verbose=4)
+
+# For XMS-PDFT
+mc_xms = mc.multi_state([.5, .5], "xms").run(verbose=4)
+
+print ('{:>21s} {:>12s} {:>12s}'.format ('state 0','state 1', 'gap'))
+fmt_str = '{:>9s} {:12.9f} {:12.9f} {:12.9f}'
+print (fmt_str.format ('CASSCF', mc_sa.e_mcscf[0], mc_sa.e_mcscf[1],
+    mc_sa.e_mcscf[1]-mc_sa.e_mcscf[0]))
+print (fmt_str.format ('MC-PDFT', mc_sa.e_states[0], mc_sa.e_states[1],
+    mc_sa.e_states[1]-mc_sa.e_states[0]))
+print (fmt_str.format ('CMS-PDFT', mc_cms.e_states[0], mc_cms.e_states[1],
+    mc_cms.e_states[1]-mc_cms.e_states[0]))
+print (fmt_str.format ('XMS-PDFT', mc_xms.e_states[0], mc_xms.e_states[1],
+    mc_xms.e_states[1]-mc_xms.e_states[0]))

examples/mcpdft/41-state_average.py

@@ -0,0 +1,80 @@
+#!/usr/bin/env python
+
+'''
+State average over states of different spins and/or spatial symmetry
+
+In MC-PDFT, you must use the method attribute rather than the function:
+
+NO  : mc = mcscf.state_average_mix (mc, [solver1, solver2], weights)
+NO  : mcscf.state_average_mix_(mc, [solver1, solver2], weights)
+YES : mc = mc.state_average_mix ([solver1, solver2], weights)
+YES : mc.state_average_mix_([solver1, solver2], weights)
+'''
+
+import numpy as np
+from pyscf import gto, scf, mcpdft, fci
+
+r = 1.8
+mol = gto.Mole()
+mol.atom = [
+    ['C', ( 0., 0.    , -r/2   )],
+    ['C', ( 0., 0.    ,  r/2)],]
+mol.basis = 'cc-pvdz'
+mol.unit = 'B'
+mol.symmetry = True
+mol.build()
+mf = scf.RHF(mol)
+mf.irrep_nelec = {'A1g': 4, 'E1gx': 0, 'E1gy': 0, 'A1u': 4,
+                  'E1uy': 2, 'E1ux': 2, 'E2gx': 0, 'E2gy': 0, 'E2uy': 0, 'E2ux': 0}
+ehf = mf.kernel()
+#mf.analyze()
+
+#
+# state-average over 1 triplet + 2 singlets
+# Note direct_spin1 solver is called here because the CI solver will take
+# spin-mix solution as initial guess which may break the spin symmetry
+# required by direct_spin0 solver
+#
+weights = np.ones(3)/3
+solver1 = fci.direct_spin1_symm.FCI(mol)
+solver1.spin = 2
+solver1 = fci.addons.fix_spin(solver1, shift=.2, ss=2)
+solver1.nroots = 1
+solver2 = fci.direct_spin0_symm.FCI(mol)
+solver2.spin = 0
+solver2.nroots = 2
+
+mc = mcpdft.CASSCF(mf, 'tPBE', 8, 8)
+mc.state_average_mix_([solver1, solver2], weights)
+
+# Mute warning msgs
+mc.check_sanity = lambda *args: None
+
+mc.verbose = 4
+mc.kernel()
+
+
+#
+# Example 2: Mix FCI wavefunctions with different symmetry irreps
+#
+mol = gto.Mole()
+mol.build(atom='''
+ O     0.    0.000    0.1174
+ H     0.    0.757   -0.4696
+ H     0.   -0.757   -0.4696
+''', symmetry=True, basis='631g')
+
+mf = scf.RHF(mol).run()
+
+weights = [.5, .5]
+solver1 = fci.direct_spin1_symm.FCI(mol)
+solver1.wfnsym= 'A1'
+solver1.spin = 0
+solver2 = fci.direct_spin1_symm.FCI(mol)
+solver2.wfnsym= 'A2'
+solver2.spin = 0
+
+mc = mcpdft.CASSCF(mf, 'tPBE', 4, 4)
+mc.state_average_mix_([solver1, solver2], weights)
+mc.verbose = 4
+mc.kernel()