Skip to content
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

Jump to bottom

GAMESS EFP Support #256

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

GAMESS EFP Support #256

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
16b90b8
support GAMESS EFP-EFP energy, empty QM molecule workaround
zachglick Jun 27, 2020
6114dba
lint and remove comment
zachglick Jun 27, 2020
b957f46
makefp support
zachglick Jun 28, 2020
1291e9e
fix formatting
zachglick Jun 28, 2020
8adce37
fix formatting
zachglick Jun 28, 2020
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
22 changes: 21 additions & 1 deletion qcengine/programs/gamess/germinate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,18 @@
from typing import Any, Dict

from ...exceptions import InputError


def muster_modelchem(method: str, derint: int) -> Dict[str, Any]:
"""Converts the QC method into GAMESS keywords."""

method = method.lower()
if method.endswith("-makefp"):
method = method[:-7]
run_makefp = True
else:
run_makefp = False

opts = {}

runtyp = {
Expand All @@ -14,7 +22,13 @@ def muster_modelchem(method: str, derint: int) -> Dict[str, Any]:
#'properties': 'prop',
}[derint]

opts["contrl__runtyp"] = runtyp
if run_makefp and runtyp != "energy":
raise InputError(f"GAMESS MAKEFP must be run with energy, not {runtyp}")

if run_makefp:
opts["contrl__runtyp"] = "makefp"
else:
opts["contrl__runtyp"] = runtyp

if method == "gamess":
pass
Expand All @@ -35,4 +49,10 @@ def muster_modelchem(method: str, derint: int) -> Dict[str, Any]:
elif method == "ccsd(t)":
opts["contrl__cctyp"] = "ccsd(t)"

elif method == "efp":
opts["contrl__coord"] = "fragonly"

else:
raise InputError(f"Unrecognized method type '{method}'.")

return opts
98 changes: 97 additions & 1 deletion qcengine/programs/gamess/harvester.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,10 @@ def harvest(p4Mol, gamessout: str, **largs) -> Tuple[PreservingDict, Molecule, l
if outMol:
outqcvar["NUCLEAR REPULSION ENERGY"] = outMol.nuclear_repulsion_energy()
if p4Mol:
if abs(outMol.nuclear_repulsion_energy() - p4Mol.nuclear_repulsion_energy()) > 1.0e-3:
# temporary hack until qcel lets us do EFP with an 'empty' QM molecule
if p4Mol.extras is not None and "efp" in p4Mol.extras:
outMol = p4Mol
elif abs(outMol.nuclear_repulsion_energy() - p4Mol.nuclear_repulsion_energy()) > 1.0e-3:
raise ValueError(
"""gamess outfile (NRE: %f) inconsistent with Psi4 input (NRE: %f)."""
% (outMol.nuclear_repulsion_energy(), p4Mol.nuclear_repulsion_energy())
Expand Down Expand Up @@ -319,6 +322,95 @@ def harvest_outfile_pass(outtext):
logger.debug("matched dft xc")
qcvar["DFT XC ENERGY"] = mobj.group(1)

# Process EFP
mobj = re.search(
r"^\s+"
+ r"ELECTROSTATIC ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"REPULSION ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"POLARIZATION ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"--------------------------------------"
+ r"\s*"
r"^\s+" + r"FINAL EFP ENERGY =" + r"\s+" + NUMBER + r"\s*",
outtext,
zachglick marked this conversation as resolved.
Show resolved Hide resolved
re.MULTILINE,
)
if mobj:
logger.debug("matched efp")
print("matched efp")
qcvar["EFP ELECTROSTATIC ENERGY"] = mobj.group(1)
qcvar["EFP REPULSION ENERGY"] = mobj.group(2)
qcvar["EFP POLARIZATION ENERGY"] = mobj.group(3)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nice for now, may have to revisit in future. i'm against a 3rd labels set. https://github.com/loriab/pylibefp/blob/master/pylibefp/wrapper.py#L961-L966

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Noted, that's going to be a headache later on

qcvar["EFP DISPERSION ENERGY"] = 0.0
qcvar["EFP CHARGE TRANSFER ENERGY"] = 0.0
qcvar["EFP TOTAL ENERGY"] = mobj.group(4)
qcvar_coord = Molecule(validate=False, symbols=[], geometry=[])

mobj = re.search(
r"^\s+"
+ r"ELECTROSTATIC ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"REPULSION ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"POLARIZATION ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"TOTAL DISPERSION ENERGY\(E6\+E7\+E8\) ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"E7 DISPERSION ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"E6 DISPERSION ENERGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"CHARGE TRANSFER ENRGY ="
+ r"\s+"
+ NUMBER
+ r"\s*"
+ r"^\s+"
+ r"--------------------------------------"
+ r"\s*"
r"^\s+" + r"FINAL EFP ENERGY =" + r"\s+" + NUMBER + r"\s*",
outtext,
re.MULTILINE,
)
if mobj:
logger.debug("matched efp")
print("matched efp")
qcvar["EFP ELECTROSTATIC ENERGY"] = mobj.group(1)
qcvar["EFP REPULSION ENERGY"] = mobj.group(2)
qcvar["EFP POLARIZATION ENERGY"] = mobj.group(3)
qcvar["EFP DISPERSION ENERGY"] = mobj.group(4)
qcvar["EFP CHARGE TRANSFER ENERGY"] = mobj.group(7)
qcvar["EFP TOTAL ENERGY"] = mobj.group(8)
qcvar_coord = Molecule(validate=False, symbols=[], geometry=[])

# Process Geometry
mobj = re.search(
# fmt: off
Expand Down Expand Up @@ -445,6 +537,10 @@ def harvest_outfile_pass(outtext):
qcvar["CURRENT REFERENCE ENERGY"] = qcvar["DFT TOTAL ENERGY"]
qcvar["CURRENT ENERGY"] = qcvar["DFT TOTAL ENERGY"]

if "EFP TOTAL ENERGY" in qcvar:
qcvar["CURRENT REFERENCE ENERGY"] = qcvar["EFP TOTAL ENERGY"]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

unless something needs it, maybe drop the current ref line. it's questionable.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had to add that line since CURRENT REFERENCE ENERGY is expected here:

QCEngine/qcengine/programs/gamess/runner.py

Line 155 in 4525eb3

retres = qcvars[f"CURRENT {input_model.driver.upper()}"]

What exactly is questionable about it? Would setting CURRENT ENERGY make more sense? Does CURRENT REFERENCE ENERGY imply a full QM calculation?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you could get away with CURRENT ENERGY only that'd be better, but if CURRENT REFERENCE ENERGY needed, fine. It's questionable b/c it's an interaction energy like SAPT rather than a total energy.

qcvar["CURRENT ENERGY"] = qcvar["EFP TOTAL ENERGY"]

if "FCI TOTAL ENERGY" in qcvar: # and 'FCI CORRELATION ENERGY' in qcvar:
qcvar["CURRENT ENERGY"] = qcvar["FCI TOTAL ENERGY"]

Expand Down
21 changes: 18 additions & 3 deletions qcengine/programs/gamess/runner.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ def compute(self, input_data: AtomicInput, config: "TaskConfig") -> AtomicResult
def build_input(
self, input_model: AtomicInput, config: "TaskConfig", template: Optional[str] = None
) -> Dict[str, Any]:
gamessrec = {"infiles": {}, "scratch_directory": config.scratch_directory}
gamessrec = {"infiles": {}, "outfiles": [], "scratch_directory": config.scratch_directory}

opts = copy.deepcopy(input_model.keywords)

Expand All @@ -107,7 +107,18 @@ def build_input(
# Handle conversion from schema (flat key/value) keywords into local format
optcmd = format_keywords(opts)

gamessrec["infiles"]["gamess.inp"] = optcmd + molcmd
# Currently passing in the entire GAMESS-formatted EFP command (i.e. the $EFRAG section)
# In the future, should be passed in as actual data fragments, coordinates, etc.)
if input_model.molecule.extras is not None and "efp" in input_model.molecule.extras:
efpcmd = input_model.molecule.extras["efp"]
gamessrec["infiles"]["gamess.inp"] = optcmd + efpcmd + molcmd
else:
gamessrec["infiles"]["gamess.inp"] = optcmd + molcmd
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what's in molecule.extras["efp"] again? I don't recall it from MolSSI/QCElemental#124 (comment) (though that needn't be a strict blueprint).

Copy link
Contributor Author

@zachglick zachglick Jun 28, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose I should rename that variable from molecule.extras["efp"] to molecule.extras["efp_molecule"] in order to be more consistent with the blueprint. In the blueprint, molecule.extras["efp_molecule"] is a schema-like dictionary defining all of the fragment locations and potentials. In this PR, molecule.extras["efp_molecule"] is instead expected to be an equivalent GAMESS-ready string, which is directly inserted into the input file.

This is a short-tem solution. Eventually, we need an input-parsing function that takes in your efp_molecule and converts it to the GAMESS-ready string. I held off on doing this right away because I wasn't sure how final the blueprint was.

I added an example script to the PR description that demonstrates how you'd run MAKEFP and then get an EFP interaction energy. This will give you a better idea of what the efp_molecule string looks like in this PR.


# Store the generated efp if running makefp
if opts["contrl__runtyp"] == "makefp":
gamessrec["outfiles"].append("gamess.efp")

gamessrec["command"] = [which("rungms"), "gamess"] # rungms JOB VERNO NCPUS >& JOB.log &

return gamessrec
Expand All @@ -133,7 +144,11 @@ def build_input(
def execute(self, inputs, extra_outfiles=None, extra_commands=None, scratch_name=None, timeout=None):

success, dexe = execute(
inputs["command"], inputs["infiles"], [], scratch_messy=False, scratch_directory=inputs["scratch_directory"]
inputs["command"],
inputs["infiles"],
inputs["outfiles"],
scratch_messy=False,
scratch_directory=inputs["scratch_directory"],
)
return success, dexe

Expand Down