efp: start pylibefp interface

qcengine/programs/base.py

@@ -7,6 +7,7 @@
 from ..exceptions import InputError, ResourceError
 from .cfour import CFOURHarness
 from .dftd3 import DFTD3Harness
+from .pylibefp import PylibEFPHarness
 from .entos import EntosHarness
 from .gamess import GAMESSHarness
 from .molpro import MolproHarness
@@ -101,3 +102,4 @@ def list_available_programs() -> Set[str]:
 register_program(NWChemHarness())
 register_program(CFOURHarness())
 register_program(EntosHarness())
+register_program(PylibEFPHarness())

qcengine/programs/pylibefp.py

@@ -0,0 +1,165 @@
+"""
+Calls the PylibEFP interface to LibEFP.
+"""
+import pprint
+from typing import Dict
+
+import qcelemental as qcel
+from qcelemental.models import Provenance, Result
+from qcelemental.util import safe_version, which_import
+
+#from ..exceptions import InputError, RandomError, ResourceError, UnknownError
+from .model import ProgramHarness
+
+pp = pprint.PrettyPrinter(width=120, compact=True, indent=1)
+
+
+class PylibEFPHarness(ProgramHarness):
+
+    _defaults = {
+        "name": "PylibEFP",
+        "scratch": False,
+        "thread_safe": False,
+        "thread_parallel": False,  # can be but not the way Psi usually builds it
+        "node_parallel": False,
+        "managed_memory": False,
+    }
+    version_cache: Dict[str, str] = {}
+
+    class Config(ProgramHarness.Config):
+        pass
+
+    @staticmethod
+    def found(raise_error: bool = False) -> bool:
+        return which_import('pylibefp',
+                            return_bool=True,
+                            raise_error=raise_error,
+                            raise_msg='Please install via `conda install pylibefp -c psi4`.')
+
+    def get_version(self) -> str:
+        self.found(raise_error=True)
+
+        which_prog = which_import('pylibefp')
+        if which_prog not in self.version_cache:
+            import pylibefp
+            self.version_cache[which_prog] = safe_version(pylibefp.__version__)
+
+        candidate_version = self.version_cache[which_prog]
+
+        if "undef" in candidate_version:
+            raise TypeError(
+                "Using custom build without tags. Please pull git tags with `git pull origin master --tags`.")
+
+        return candidate_version
+
+    def compute(self, input_model: 'ResultInput', config: 'JobConfig') -> 'Result':
+        """
+        Runs PylibEFP in API mode
+        """
+        self.found(raise_error=True)
+
+        #        if parse_version(self.get_version()) < parse_version("1.2"):
+        #            raise ResourceError("Psi4 version '{}' not understood.".format(self.get_version()))
+
+        # Setup the job
+        input_data = input_model.dict(encoding="json")
+        #        input_data["nthreads"] = config.ncores
+        #        input_data["memory"] = int(config.memory * 1024 * 1024 * 1024 * 0.95)  # Memory in bytes
+        input_data["success"] = False
+        input_data["return_output"] = True
+
+        import pylibefp
+        efpobj = pylibefp.core.efp()
+
+        efp_extras = input_model.molecule.extras['efp_molecule']['extras']
+        efpobj.add_potential(efp_extras['fragment_files'])
+        efpobj.add_fragment(efp_extras['fragment_files'])
+        for ifr, (hint_type, geom_hint) in enumerate(zip(efp_extras['hint_types'], efp_extras['geom_hints'])):
+            print('SEt_frag_coordinates', ifr, hint_type, geom_hint)
+            efpobj.set_frag_coordinates(ifr, hint_type, geom_hint)
+        efpobj.prepare()
+
+        # print efp geom in [A]
+        print(efpobj.banner())
+        print(efpobj.geometry_summary(units_to_bohr=qcel.constants.bohr2angstroms))
+
+        if input_model.model.method != 'efpefp':
+            raise InputError
+
+        # set keywords
+        efpobj.set_opts(input_model.keywords)
+        #efpobj.set_opts(efpopts, label='psi', append='psi')
+
+        if input_model.driver == 'energy':
+            do_gradient = False
+        elif input_model.driver == 'gradient':
+            do_gradient = True
+        else:
+            raise InputError
+
+        # compute and report
+        efpobj.compute(do_gradient=do_gradient)
+        print(efpobj.energy_summary(label='psi'))
+
+        ene = efpobj.get_energy(label='psi')
+
+        pp.pprint(ene)
+        print('<<< get_opts():  ', efpobj.get_opts(), '>>>')
+        #print('<<< summary():   ', efpobj.summary(), '>>>')
+        print('<<< get_energy():', ene, '>>>')
+        print('<<< get_atoms(): ', efpobj.get_atoms(), '>>>')
+        print(efpobj.energy_summary())
+        print(efpobj.geometry_summary(units_to_bohr=qcel.constants.bohr2angstroms))
+        print(efpobj.geometry_summary(units_to_bohr=1.0))
+
+        ###### psi4 proc
+        #def run_efp(name, **kwargs):
+        #        try:
+        #            efpobj = efp_molecule.EFP
+        #        except AttributeError:
+        #            raise ValidationError("""Method 'efp' not available without EFP fragments in molecule""")
+        #
+        #        core.set_variable('EFP ELST ENERGY', ene['electrostatic'] + ene['charge_penetration'] + ene['electrostatic_point_charges'])
+        #        core.set_variable('EFP IND ENERGY', ene['polarization'])
+        #        core.set_variable('EFP DISP ENERGY', ene['dispersion'])
+        #        core.set_variable('EFP EXCH ENERGY', ene['exchange_repulsion'])
+        #        core.set_variable('EFP TOTAL ENERGY', ene['total'])
+        #        core.set_variable('CURRENT ENERGY', ene['total'])
+        #
+        #        if do_gradient:
+        #            core.print_out(efpobj.gradient_summary())
+        #
+        #            core.set_variable('EFP TORQUE', torq)
+        #
+        #                output_data = input_data
+
+        if input_model.driver == 'energy':
+            retres = ene['total']
+
+
+#        elif input_model.driver == 'gradient':
+#            torq = efpobj.get_gradient()
+#            #torq = core.Matrix.from_array(np.asarray(torq).reshape(-1, 6))
+#            retres = torq
+
+        output_data = {
+            'schema_name': 'qcschema_output',
+            'schema_version': 1,
+            # 'extras': {
+            #     'outfiles': outfiles,
+            # },
+            'properties': {},
+            'provenance': Provenance(creator="PylibEFP", version=self.get_version(), routine="pylibefp"),
+            'return_result': retres,
+            # 'stdout': stdout,
+        }
+
+        #        # got to even out who needs plump/flat/Decimal/float/ndarray/list
+        #        # Decimal --> str preserves precision
+        #        output_data['extras']['qcvars'] = {
+        #            k.upper(): str(v) if isinstance(v, Decimal) else v
+        #            for k, v in qcel.util.unnp(qcvars, flat=True).items()
+        #        }
+
+        output_data['success'] = True
+        return Result(**{**input_model.dict(), **output_data})

qcengine/programs/tests/test_pylibefp.py

@@ -0,0 +1,62 @@
+from qcelemental.testing import compare, compare_recursive, compare_values
+
+import qcengine as qcng
+from qcengine.testing import using_pylibefp
+
+b2a = 0.529177
+a2b = 1.0 / b2a
+
+
+@using_pylibefp
+def test_total_1a():
+
+    resi = {
+        'molecule': {
+            'geometry': [0, 0, 0],  # dummy
+            'symbols': ['He'],  # dummy
+            'extras': {
+                'efp_molecule': {
+                    'geometry': [],
+                    'symbols': [],
+                    'extras': {
+                        'fragment_files': ['h2o', 'nh3'],
+                        'hint_types': ['xyzabc', 'xyzabc'],
+                        'geom_hints': [
+                            [0.0 * a2b, 0.0 * a2b, 0.0 * a2b, 1.0, 2.0, 3.0],
+                            [5.0 * a2b, 0.0 * a2b, 0.0 * a2b, 5.0, 2.0, 8.0],
+                        ]
+                    }
+                }
+            }
+        },
+        'driver': 'energy',
+        'model': {
+            'method': 'efpefp',
+        },
+        'keywords': {
+            'elec': True,
+            'elec_damp': 'screen',
+            'xr': True,
+            'pol': True,
+            'disp': True,
+            'disp_damp': 'tt',
+        }
+    }
+
+    res = qcng.compute(resi, 'pylibefp', raise_error=True, return_dict=False)
+
+    atol = 1.e-6
+    assert compare_values(0.0001922903, res.return_result, atol=atol)
+
+
+#    expected_ene = blank_ene()
+#    expected_ene['elec'] = expected_ene['electrostatic'] = 0.0002900482
+#    expected_ene['xr'] = expected_ene['exchange_repulsion'] = 0.0000134716
+#    expected_ene['pol'] = expected_ene['polarization'] = 0.0002777238 - expected_ene['electrostatic']
+#    expected_ene['disp'] = expected_ene['dispersion'] = -0.0000989033
+#    expected_ene['total'] = 0.0001922903
+#    assert compare(2, asdf.get_frag_count(), sys._getframe().f_code.co_name + ': nfrag')
+#    assert compare_values(0.0, asdf.get_frag_charge(1), sys._getframe().f_code.co_name + ': f_chg', atol=1.e-6)
+#    assert compare(1, asdf.get_frag_multiplicity(1), sys._getframe().f_code.co_name + ': f_mult')
+#    assert compare('NH3', asdf.get_frag_name(1), sys._getframe().f_code.co_name + ': f_name')
+#    assert compare_recursive(expected_ene, ene, sys._getframe().f_code.co_name + ': ene', atol=1.e-6)

qcengine/testing.py

@@ -131,6 +131,7 @@ def get_job(self):
     "cfour": which('xcfour', return_bool=True),
     "gamess": which('rungms', return_bool=True),
     "nwchem": which('nwchem', return_bool=True),
+    "pylibefp": which_import("pylibefp", return_bool=True),
 }
 
 
@@ -158,6 +159,7 @@ def _build_pytest_skip(program):
 using_cfour = _build_pytest_skip("cfour")
 using_gamess = _build_pytest_skip("gamess")
 using_nwchem = _build_pytest_skip("nwchem")
+using_pylibefp = _build_pytest_skip("pylibefp")
 
 using_dftd3_321 = pytest.mark.skipif(is_program_new_enough("dftd3", "3.2.1") is False,
                                      reason='DFTD3 does not include 3.2.1 features. Update package and add to PATH')