GAMESS EFP Support

[zachglick]

Description

Adds support for GAMESS's MAKEFP method, which performs a single-point energy calculation on a monomer, then generates an effective fragment potential (EFP) from the wavefunction of that monomer. A MAKEFP calculation is specified with QCEngine by appending -makefp to the name of the method (i.e. hf-makefp runs MAKEFP at the Hartree-Fock level). The EFP file generated by GAMESS is stored in extras->outfiles->gamess.efp.

Also adds support for GAMESS calculations that make use of EFPs. To run an EFP calculation, the user must provide the entire GAMESS $EFRAG block as a string (this block contains fragment positions as well as the fragment potential(s) generated with MAKEFP). This string is expected in Molecule->extras->efp.

Some things that would help with this (relevant to MolSSI/QCElemental#124):

• Support for "empty" QCElemental Molecule. (Since EFP fragments are are input through Molecule.extras, a Molecule used in EFP calculations will have no actual QM-treated atoms.)
• Agree on a standard for specifying EFP fragments and EFP potentials, and if/how they should be validated. With this PR, the user will be expected to handle this on their own by providing the entire GAMESS-ready $EFRAG block.

An example script demonstrating how both of these new features can be used together:

GAMESS EFP Demo

from qcelemental.testing import compare_values
from qcelemental.models import Molecule, AtomicInput

import qcengine as qcng

if __name__ == "__main__":

  # H2O monomer
  mol = Molecule.from_data("""
  # R=0.958 A=104.5
  H                  0.000000000000     1.431430901356     0.984293362719
  O                  0.000000000000     0.000000000000    -0.124038860300
  H                  0.000000000000    -1.431430901356     0.984293362719
  units au
  """)

  # run GAMESS MAKEFP on a HF wavefunction through QCEngine
  # this is a normal HF calculation, with some additional post-processing
  inp = AtomicInput(
      molecule=mol,
      driver="energy",
      model={"method": "hf-makefp", "basis": "accd"},
      keywords={'contrl__ispher' : '1'},
  )
  res = qcng.compute(inp, program="gamess", raise_error=True, return_dict=True)
  assert res["success"] is True

  # read the .efp file generated by GAMESS for H2O
  fragment_potential = res["extras"]["outfiles"]["gamess.efp"]

  # remove the first two lines of the efp file
  fragment_potential = '\n'.join(fragment_potential.split('\n')[2:])

  # next, get the geometry of two water fragments
  # I've manually written in GAMESS syntax for now
  # in the future, users will be able to specify fragment geometry with a qcel mol
  fragment_geometry = """ $efrag
coord=cart
FRAGNAME=FRAGNAME
A01H    0.000000000000    1.423241232739    0.978661909286      
A02O    0.000000000000    0.000000000000   -0.123329194776      
A03H    0.000000000000   -1.423241232739    0.978661909286      
FRAGNAME=FRAGNAME
A01H    0.000000000000    0.000000000000   -3.9   
A02O    0.000000000000    0.000000000000   -5.5
A03H   -1.3               0.000000000000   -7.1
$end
  """
  # the entire EFP part of a GAMESS input file
  efp_command = fragment_geometry.rstrip() + '\n ' + fragment_potential.lstrip()

  # now we'll make a dummy molecule to hold the efp command
  # the symbols and geometry don't matter at all 
  dummy_mol = Molecule(
      symbols=mol.symbols,
      geometry=mol.geometry,
      extras={"efp" : efp_command}
  )

  # evaluate the efp energy of the H2O dimer using our calculated EFP potential
  inp2 = AtomicInput(
      molecule=dummy_mol,
      driver="energy",
      model={"method": "efp"},
      keywords={},
  )
  res2 = qcng.compute(inp2, program="gamess", raise_error=True, return_dict=True)
  assert res2["success"] is True

  efp_ref = -0.0068308265
  efp_ene = res2["return_result"]
  atol = 1.0e-8

  assert compare_values(efp_ref, efp_ene, atol=atol)
  print(f'computed interaction energy is {627.509*efp_ene:0.5f} kcal/mol')

Changelog description

Support for generating GAMESS EFP and calculations using EFPs

Status

• Code base linted
• Ready to go

[codecov]

Codecov Report

Merging #256 into master will increase coverage by 0.02%.
The diff coverage is n/a.

[loriab]

looking good

[loriab]

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).

[zachglick]

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.

[loriab]

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

[zachglick]

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?

[loriab]

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.

[loriab]

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

[zachglick]

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