[Feature] Add Molpro as a backend via pymolpro

basisopt/api.py

@@ -149,6 +149,18 @@
     bo_logger.info("ORCA install dir at: %s", path)
 
 
+@register_backend
+def molpro(path: str):
+    """Tests pymolpro import and prepares to be used as calculation backend"""
+    try:
+        global _CURRENT_BACKEND
+        from basisopt.wrappers.molpro import MolproWrapper
+
+        _CURRENT_BACKEND = MolproWrapper()
+    except ImportError:
+        bo_logger.error("Molpro backend (using pymolpro) not found!")
+
+
 def run_calculation(
     evaluate: str = 'energy', mol: Molecule = None, params: dict[Any, Any] = {}
 ) -> int:

basisopt/wrappers/molpro.py

@@ -0,0 +1,154 @@
+# Wrappers for Molpro functionality using pymolpro (https://github.com/molpro/pymolpro)
+from pymolpro import Project
+
+from basisopt.bse_wrapper import internal_basis_converter
+from basisopt.containers import InternalBasis
+from basisopt.exceptions import FailedCalculation
+from basisopt.molecule import Molecule
+from basisopt.wrappers.wrapper import Wrapper, available
+
+
+class MolproWrapper(Wrapper):
+    """Wrapper for Molpro using pymolpro"""
+
+    def __init__(self):
+        super().__init__(name='Molpro')
+        self._method_strings = {
+            'hf': ['energy'],
+            'rhf': ['energy'],
+            'uhf': ['energy'],
+            'mp2': ['energy'],
+            'rmp2': ['energy'],
+            'ump2': ['energy'],
+            'ccsd': ['energy'],
+            'ccsd(t)': ['energy'],
+            'uccsd': ['energy'],
+            'uccsd(t)': ['energy'],
+            'rks': ['energy'],
+            'uks': ['energy'],
+        }
+
+    def convert_molecule(self, m: Molecule) -> str:
+        """Convert an internal Molecule object
+        to a Molpro geometry section and set charge/multiplicity
+        """
+        molstring = "geomtyp=xyz\n"
+        molstring += "geom={\n"
+        for i in range(m.natoms()):
+            molstring += m.get_line(i) + "\n"
+        molstring += "}\n"
+        molstring += f"set,charge={m.charge}\n"
+        spin = m.multiplicity - 1
+        molstring += f"set,spin={spin}\n"
+        return molstring
+
+    def _convert_params(self, method: str, **params) -> str:
+        """Helper function to pluck out the relevant
+        parameters for a given method.
+        Returns an empty string if the params are not found.
+        """
+        method_search = method + '-params'
+        param_options = params.get(method_search, '')
+        return param_options
+
+    def _command_string(self, method: str, **params) -> str:
+        """Helper function to turn an internal method
+        name into a Molpro command line
+        """
+        # Extract the relevant parameters for this method
+        method_options = self._convert_params(method, **params)
+        # Post-HF calculations need to specify the HF part too
+        # hence check for any user-supplied params for the HF part.
+        if method in ['mp2', 'ccsd', 'ccsd(t)']:
+            ref_options = self._convert_params('hf', **params)
+            command = f"{{hf{ref_options}}}\n{{{method}{method_options}}}"
+        elif method in ['rmp2', 'uccsd', 'uccsd(t)']:
+            ref_options = self._convert_params('rhf', **params)
+            command = f"{{rhf{ref_options}}}\n{{{method}{method_options}}}"
+        elif method in ['ump2']:
+            ref_options = self._convert_params('uhf', **params)
+            command = f"{{uhf{ref_options}}}\n{{{method}{method_options}}}"
+        # For DFT based methods, the functional must be specified separately
+        # to other params
+        elif method in ['rks', 'uks']:
+            if "functional" in params:
+                xcfun = params["functional"]
+                command = f"{{{method},{xcfun}{method_options}}}"
+            else:
+                raise KeyError("DFT functional not specified")
+        else:
+            command = f"{{{method}{method_options}}}"
+        return command + "\n"
+
+    def _convert_basis(self, basis: InternalBasis) -> str:
+        """Converts an InternalBasis to the basis string for Molpro
+        May be updated in future to enable auxiliary fitting sets
+
+        Arguments:
+            basis (InternalBasis): basis set to convert
+
+        Returns:
+            Molpro basis block string
+        """
+        molpro_basis = internal_basis_converter(basis, fmt="molpro").split("\n")
+        basis_str = ""
+        for line in molpro_basis[1:-1]:
+            basis_str += line + "\n"
+        return basis_str
+
+    def initialise(self, proj: Project, m: Molecule, tmp: str = "", **params):
+        """Initialises pymolpro before each calculation
+        - creates a pymolpro Project
+        - converts molecule
+        - sets options from globals
+        - converts basis set
+        """
+
+        # Handle options, molecule, basis
+        g_param_str = ""
+        glo_params = params.get("global-params", '')
+        if glo_params:
+            g_param_str = glo_params + "\n"
+        cmd = self._command_string(m.method, **params)
+        mol = self.convert_molecule(m)
+        basis = self._convert_basis(m.basis)
+
+        # Assemble into an input string and pass to the Project
+        molpro_input = g_param_str + mol + basis + cmd
+        print(molpro_input)
+        proj.write_input(molpro_input)
+
+    def _get_energy(self, proj: Project, meth: str) -> float:
+        """Helper function to retrieve the energy from
+        Molpro after a calculation"""
+        if meth in ['hf', 'rhf', 'uhf', 'rks', 'uks']:
+            energy = proj.energies()[0]
+        elif meth in ['mp2']:
+            energy = proj.energies()[-1]
+        elif meth in ['rmp2', 'ump2', 'ccsd']:
+            energy = proj.energies(method=f"{meth.upper()}")[-1]
+        elif meth in ['uccsd', 'uccsd(t)']:
+            energy = proj.energy(method=f"RHF-{meth.upper()}")
+        else:
+            energy = proj.energy(method=f"{meth.upper()}")
+        return energy
+
+    @available
+    def energy(self, mol, tmp="", **params):
+        name = "energy"
+        proj_name = f"{mol.name}-{mol.method}-" + name
+        p = Project(proj_name, location=tmp)
+        self.initialise(p, mol, tmp=tmp, **params)
+        p.run(wait=True)
+        if p.errors():
+            raise FailedCalculation
+        # Attempt to catch race cases where wait=True doesn't seem to be sufficient
+        p.wait()
+        # TODO use a helper function to extract the energy for a particular method
+        energy = self._get_energy(p, mol.method)
+        #        if mol.method == "hf":
+        #            energy = p.energy()
+        #        else:
+        #            energy = p.energy(method=f"{mol.method.upper()}")
+        p.clean()
+        return energy

doc/src/_tutorials/multimol.rst

@@ -5,7 +5,7 @@
 Multiple-molecule basis optimization
 ====================================
 
-This example demonstrates how BasisOpt can be used to optimize a double-zeta molecule-optimized basis set. In this case five molecules (water, methane, methanol, formaldehyde and oxygen) are used for the optimization, resulting in a single set of exponents for hydrogen, carbon and oxygen. Please note that as the exponents are tightly coupled this optimization can take a long time to run. The full example script can be found in ``examples/multi-molecule/multi-molecule.py``.
+This example demonstrates how BasisOpt can be used to optimize a double-zeta molecule-optimized basis set. In this case five molecules (water, methane, methanol, formaldehyde and oxygen) are used for the optimization, resulting in a single set of exponents for hydrogen, carbon and oxygen. Please note that as the exponents are tightly coupled this optimization can take a long time to run. The full example script can be found in ``examples/multi-molecule/multi_molecule.py``.
 
 This example also shows the use of the logger functionality in BasisOpt.
 
@@ -18,19 +18,19 @@ In this example we use the Psi4 program to carry out the electronic structure ca
 
 	mb = MolecularBasis(name="double")
 	list_of_mols = ['water', 'methane', 'methanol', 'formaldehyde', 'oxygen']
-	mol_objs = [
-    	bo.molecule.Molecule.from_xyz(mol+'.xyz', name=mol)
-    	for mol in list_of_mols
-		]
+	mol_objs = [bo.molecule.Molecule.from_xyz(mol + '.xyz', name=mol) for mol in list_of_mols]
 	mb = MolecularBasis(name="double", molecules=mol_objs)
 
-We also set some calculation parameters that will be used in the optimization strategy. Specifically, we select the :math:`\omega` B97X-D density functional and turn off the use of density fitting in Psi4.
+We also set some calculation parameters that will be used in the optimization strategy. Specifically, we select the :math:`\omega` B97X-D density functional, adjust the DFT integration grid, and turn off the use of density fitting in Psi4.
 
 .. code-block:: python
 
 	params = {
     	'functional': "wb97x-d",
-    	'scf_type': "pk"
+    	'scf_type': "pk",
+		'dft_spherical_points': "974",
+		'dft_radial_points': "175",
+		'dft_pruning_scheme': "none"
 	}
 
 Optimization strategy and first run
@@ -68,7 +68,7 @@ Within these iterations we also use the logging capability in BasisOpt to produc
 	    bo_logger.info("Starting consistency iteration %d", counter+1)
 	    mb.optimize()
 	    e_opt.append(strategy.last_objective)
-	    e_diff = strategy.delta_objective
+	    e_diff = abs(strategy.last_objective - e_opt[counter])
 	    bo_logger.info("Objective function difference from previous iteration: %f\n", e_diff)
 	    counter += 1
 

doc/src/molpro.rst

@@ -0,0 +1,80 @@
+.. basisopt molpro wrapper file
+
+.. _`sec:molpro`:
+
+=========================
+Using Molpro as a backend
+=========================
+
+The Molpro backend for BasisOpt uses the pymolpro Python library. Both the commercial 
+`Molpro quantum chemistry package <https://www.molpro.net>`_ and the free `pymolpro library <https://molpro.github.io/pymolpro/index.html#>`_ 
+must be available to use the backend. The backend can then specified in the same manner as for other quantum chemistry packages.
+
+.. code-block:: python 
+
+	bo.set_backend('molpro')
+	bo.set_tmp_dir('/tmp/')
+
+This guide is intended to highlight how to specify certain types of Molpro calculation via BasisOpt, including passing Molpro options and
+directives.
+
+Density Functional Theory calculations
+======================================
+
+In Molpro, Kohn-Sham DFT calculations are requested using either the ``rks`` or ``uks`` commands, for spin-restricted and spin-unrestricted, respectively. The exchange-correlation functional must also be specified. To maintain consistency with other BasisOpt backends, the choice of
+functional is specified by setting the ``functional`` params. For example, a PBE0 calculation using the spin-restricted Kohn-Sham program
+can be setup as:
+
+.. code-block:: python 
+
+	params = {'functional': "pbe0"}
+	strategy = Strategy()
+	mb.setup(method='rks', strategy=strategy, params=params)
+
+Specifying options and directives
+=================================
+
+A `params` block is used to pass options and directives to Molpro on a per-method basis. This is done using a *method*-params key and, as
+an example, requesting that no orbitals are treated with a frozen-core approximation in an MP2 calculation can be specified in the 
+following way:
+
+.. code-block:: python 
+
+	params = {'mp2-params': ";core,0"}
+	strategy = Strategy()
+	mb.setup(method='mp2', strategy=strategy, params=params)
+
+Note that the ``;`` and ``,`` separators used in Molpro input files must be included as part of the params values.
+
+Specifying global parameters and thresholds
+===========================================
+
+Molpro global parameters and thresholds are also specified in the params dictionary, using the global-params key. For example,
+to set a global energy convergence threshold:
+
+.. code-block:: python 
+
+	params = {'global-params': "gthresh,energy=1.d-12"}
+
+Post-HF methods
+===============
+
+When running a post-Hartree-Fock calculation in Molpro, the reference wavefunction must also be specified. In BasisOpt, the reference
+calculation is somewhat *hard-coded* and will be determined from the post-HF method selected. Currently, requesting a method of `mp2`, 
+`ccsd` or `ccsd(t)` will use a `hf` reference. The methods `rmp2`, `uccsd` and `uccsd(t)` will use an `rhf` reference, and `ump2` uses
+the spin-unrestricted `uhf` reference.
+
+It is possible to pass Molpro options and/or directives to both the reference and the correlated calculation. For example, to increase
+the accuracy of the HF reference and also request that no orbitals are treated with a frozen-core approximation in an MP2 calculation:
+
+.. code-block:: python 
+
+	params = {'hf-params': ";accu,14", 'mp2-params': ";core,0"}
+	strategy = Strategy()
+	mb.setup(method='mp2', strategy=strategy, params=params)
+
+Parameters for an RHF reference would use an ``rhf-params`` key, while a UHF reference would require ``uhf-params``.
+
+
+.. toctree::
+   :hidden:

examples/multi-molecule/multi-molecule.py → examples/multi-molecule/multi_molecule.py

@@ -11,7 +11,13 @@
 mol_objs = [bo.molecule.Molecule.from_xyz(mol + '.xyz', name=mol) for mol in list_of_mols]
 mb = MolecularBasis(name="double", molecules=mol_objs)
 
-params = {'functional': "wb97x-d", 'scf_type': "pk"}
+params = {
+    'functional': "wb97x-d",
+    'scf_type': "pk",
+    'dft_spherical_points': "974",
+    'dft_radial_points': "175",
+    'dft_pruning_scheme': "none",
+}
 
 strategy = Strategy()
 strategy.params = params
@@ -31,7 +37,7 @@
     bo_logger.info("Starting consistency iteration %d", counter + 1)
     mb.optimize()
     e_opt.append(strategy.last_objective)
-    e_diff = strategy.delta_objective
+    e_diff = abs(strategy.last_objective - e_opt[counter])
     bo_logger.info("Objective function difference from previous iteration: %f\n", e_diff)
     counter += 1
 

examples/multi-molecule/opt_basis.txt

@@ -2,45 +2,45 @@ spherical
 basis={
 !
 ! hydrogen             (4s,1p) -> [4s,1p]
-s, H , 1.797959e+01, 2.758758e+00, 6.245358e-01, 1.517951e-01
+s, H , 1.798531e+01, 2.759596e+00, 6.247444e-01, 1.517825e-01
 c, 1.1, 1.000000e+00
 c, 2.2, 1.000000e+00
 c, 3.3, 1.000000e+00
 c, 4.4, 1.000000e+00
-p, H , 8.438481e-01
+p, H , 8.446020e-01
 c, 1.1, 1.000000e+00
 !
 ! carbon               (7s,4p,1d) -> [7s,4p,1d]
-s, C , 2.958169e+03, 4.440114e+02, 1.011358e+02, 2.862702e+01, 9.097819e+00, 3.047830e+00, 4.275333e-01
+s, C , 2.958325e+03, 4.440341e+02, 1.011409e+02, 2.862842e+01, 9.098225e+00, 3.047955e+00, 4.275739e-01
 c, 1.1, 1.000000e+00
 c, 2.2, 1.000000e+00
 c, 3.3, 1.000000e+00
 c, 4.4, 1.000000e+00
 c, 5.5, 1.000000e+00
 c, 6.6, 1.000000e+00
 c, 7.7, 1.000000e+00
-p, C , 1.337480e+01, 3.008887e+00, 8.123437e-01, 2.404405e-01
+p, C , 1.337360e+01, 3.008598e+00, 8.122578e-01, 2.403973e-01
 c, 1.1, 1.000000e+00
 c, 2.2, 1.000000e+00
 c, 3.3, 1.000000e+00
 c, 4.4, 1.000000e+00
-d, C , 8.322163e-01
+d, C , 8.324954e-01
 c, 1.1, 1.000000e+00
 !
 ! oxygen               (7s,4p,1d) -> [7s,4p,1d]
-s, O , 2.309354e+03, 3.473414e+02, 7.891545e+01, 2.185828e+01, 6.678970e+00, 1.004907e+00, 2.891048e-01
+s, O , 2.309352e+03, 3.473404e+02, 7.891515e+01, 2.185820e+01, 6.678948e+00, 1.004889e+00, 2.890962e-01
 c, 1.1, 1.000000e+00
 c, 2.2, 1.000000e+00
 c, 3.3, 1.000000e+00
 c, 4.4, 1.000000e+00
 c, 5.5, 1.000000e+00
 c, 6.6, 1.000000e+00
 c, 7.7, 1.000000e+00
-p, O , 1.745149e+01, 3.818731e+00, 1.035449e+00, 2.602673e-01
+p, O , 1.745163e+01, 3.818766e+00, 1.035460e+00, 2.602755e-01
 c, 1.1, 1.000000e+00
 c, 2.2, 1.000000e+00
 c, 3.3, 1.000000e+00
 c, 4.4, 1.000000e+00
-d, O , 9.894508e-01
+d, O , 9.895363e-01
 c, 1.1, 1.000000e+00
 }