Skip to content

MC-PDFT Hybrid Energy Decomposition #76

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Nov 1, 2024
Merged

MC-PDFT Hybrid Energy Decomposition #76

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
645d4ea
allow hybrid energy decomp and check for SS
matthew-hennefarth Oct 31, 2024
7611319
add test for sa hybrid decomp
matthew-hennefarth Oct 31, 2024
6bd5591
fix flake
matthew-hennefarth Oct 31, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 32 additions & 6 deletions pyscf/mcpdft/mcpdft.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,9 +206,13 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None,
E(MC-SCF) = E0 + E1 + E2c + Enc
E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc

Only compatible with pure translated or fully-translated fnals. If
mc.fcisolver.nroots > 1, lists are returned for everything except
the nuclear potential energy.
For hybrid functionals,

E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc + Enc

Where the Enc and EOTx/c terms are premultiplied by the hybrid factor. If
mc.fcisolver.nroots > 1, lists are returned for everything except the
nuclear potential energy.

Args:
mc : an instance of CASSCF or CASCI class
Expand Down Expand Up @@ -246,6 +250,8 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None,
EOTc = correlation part of translated functional
e_ncwfn : float or list of length nroots
E2ncc = <H> - E0 - E1 - E2c
If hybrid functional, this term is weighted appropriately. For pure
functionals, it is the full NC component
'''
if verbose is None: verbose = mc.verbose
log = logger.new_logger(mc, verbose)
Expand All @@ -271,8 +277,8 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None,
)

hyb_x, hyb_c = ot._numint.hybrid_coeff(ot.otxc)
if hyb_x > 1e-10 or hyb_c > 1e-10:
raise NotImplementedError("Decomp for hybrid PDFT fnals")
if abs(hyb_x - hyb_c) > 1e-11:
raise NotImplementedError("hybrid functionals with different exchange, correlations components")
if not isinstance(ot, transfnal):
raise NotImplementedError("Decomp for non-translated PDFT fnals")

Expand Down Expand Up @@ -309,11 +315,27 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None,
def _get_e_decomp(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None):
ncore = mc.ncore
ncas = mc.ncas
if ot is None: ot = mc.otfnal
if ot is None: ot = [mc.otfnal,]
if mo_coeff is None: mo_coeff = mc.mo_coeff
if ci is None: ci = mc.ci
if verbose is None: verbose = mc.verbose

if len(ot) == 1:
hyb_x, hyb_c = ot[0]._numint.hybrid_coeff(ot[0].otxc)

elif len(ot) == 2:
hyb_x, hyb_c = [
fnal._numint.hybrid_coeff(fnal.otxc)[idx] for idx, fnal in enumerate(ot)
]

else:
raise ValueError("ot must be length of 1 or 2")

if abs(hyb_x - hyb_c) > 1e-11:
raise NotImplementedError(
"hybrid functionals with different exchange, correlations components"
)

casdm1s = mc.make_one_casdm1s(ci, state=state)
casdm1 = casdm1s[0] + casdm1s[1]
casdm2 = mc.make_one_casdm2(ci, state=state)
Expand All @@ -334,6 +356,10 @@ def _get_e_decomp(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None):
max_memory=mc.max_memory)
for fnal in ot]
e_ncwfn = e_mcscf - e_nuc - e_1e - e_coul

if abs(hyb_x) > 1e-10:
e_ncwfn = hyb_x * e_ncwfn

return e_1e, e_coul, e_otxc, e_ncwfn

class _mcscf_env:
Expand Down
120 changes: 120 additions & 0 deletions pyscf/mcpdft/test/test_mcpdft.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,6 +256,28 @@ def test_decomposition_ss(self): # TODO
with self.subTest(objtype=objtype, symmetry=bool(ix), term="sanity"):
self.assertAlmostEqual(np.sum(test[:-1]), obj.e_tot, 9)

def test_decomposition_hybrid(self):
ref = [
1.0583544218,
-12.5375911135,
5.8093938665,
-1.6287258807,
-0.0619586538,
-0.5530763650,
]
terms = ["nuc", "core", "Coulomb", "OT(X)", "OT(C)", "WFN(XC)"]
for ix, mc in enumerate(mcp[0]):
mc_scf = mcpdft.CASSCF(mc, "tPBE0", 5, 2).run()
mc_ci = mcpdft.CASCI(mc, "tPBE0", 5, 2).run()
for obj, objtype in zip((mc_scf, mc_ci), ("CASSCF", "CASCI")):
test = obj.get_energy_decomposition()
for t, r, term in zip(test, ref, terms):
with self.subTest(objtype=objtype, symmetry=bool(ix), term=term):
self.assertAlmostEqual(t, r, delta=1e-5)
with self.subTest(objtype=objtype, symmetry=bool(ix), term="sanity"):
self.assertAlmostEqual(np.sum(test), obj.e_tot, 9)


def test_decomposition_sa(self):
ref_nuc = 1.0583544218
ref_states = np.array(
Expand Down Expand Up @@ -347,6 +369,104 @@ def test_decomposition_sa(self):
self.assertAlmostEqual(
np.sum(test[:-1]) + test_nuc, e_ref[state], 9
)

def test_decomposition_hybrid_sa(self):
ref_nuc = 1.0583544218
ref_states = np.array(
[
[
-12.5385413915,
5.8109724796,
-1.6290249990,
-0.0619850920,
-0.5531991067,
],
[
-12.1706553996,
5.5463231972,
-1.6201152933,
-0.0469559736,
-0.5485771470,
],
[
-12.1768195314,
5.5632261670,
-1.6164436229,
-0.0492571730,
-0.5471763843,
],
[
-12.1874226655,
5.5856701424,
-1.6111471613,
-0.0474456546,
-0.5422714244,
],
[
-12.1874226655,
5.5856701424,
-1.6111480360,
-0.0474456745,
-0.5422714244,
],
]
)
terms = ["core", "Coulomb", "OT(X)", "OT(C)", "WFN(XC)"]
for ix, (mc, ms) in enumerate(zip(mcp[1], [0, 1, "mixed"])):
s = bool(ix // 2)
mc_scf = mcpdft.CASSCF(mc, "tPBE0", 5, 2)
if ix == 0:
mc_scf = mc_scf.state_average(mc.weights)
else:
mc_scf = mc_scf.state_average_mix(mc.fcisolver.fcisolvers, mc.weights)
mc_scf.run(ci=mc.ci, mo_coeff=mc.mo_coeff)
objs = [
mc_scf,
]
objtypes = [
"CASSCF",
]
if ix != 1: # There is no CASCI equivalent to mcp[1][1]
mc_ci = mcpdft.CASCI(mc, "tPBE0", 5, 2)
mc_ci.fcisolver.nroots = (
5 - ix
) # Just check the A1 roots when symmetry is enabled
mc_ci.ci = mc.ci[: 5 - ix]
mc_ci.kernel()
objs.append(mc_ci)
objtypes.append("CASCI")
for obj, objtype in zip(objs, objtypes):
test = obj.get_energy_decomposition()
test_nuc, test_states = test[0], np.array(test[1:]).T
# Arrange states in ascending energy order
e_states = getattr(obj, "e_states", obj.e_tot)
idx = np.argsort(e_states)
test_states = test_states[idx, :]
e_ref = np.array(e_states)[idx]
with self.subTest(
objtype=objtype, symmetry=s, triplet_ms=ms, term="nuc"
):
self.assertAlmostEqual(test_nuc, ref_nuc, 9)
for state, (test, ref) in enumerate(zip(test_states, ref_states)):
for t, r, term in zip(test, ref, terms):
with self.subTest(
objtype=objtype,
symmetry=s,
triplet_ms=ms,
term=term,
state=state,
):
self.assertAlmostEqual(t, r, delta=1e-5)
with self.subTest(
objtype=objtype,
symmetry=s,
triplet_ms=ms,
term="sanity",
state=state,
):
self.assertAlmostEqual(
np.sum(test) + test_nuc, e_ref[state], 9
)

def test_energy_tot(self):
# Test both correctness and energy_tot function purity
Expand Down
Loading