Merge pull request #239 from SCM-NV/template

 ENH: Add a template for cell optimization calculation

src/qmflows/__init__.py

@@ -13,7 +13,7 @@
 from .packages import (
     adf, cp2k, cp2k_mm, dftb, orca, run, PackageWrapper)
 
-from .templates import (freq, geometry, singlepoint, ts, md)
+from .templates import (freq, geometry, singlepoint, ts, md, cell_opt)
 from .settings import Settings
 from .examples import (example_H2O2_TS, example_freqs, example_generic_constraints,
                        example_partial_geometry_opt)
@@ -25,5 +25,6 @@
     'Angle', 'Dihedral', 'Distance', 'Settings',
     'adf', 'cp2k', 'cp2k_mm', 'dftb', 'orca', 'run', 'PackageWrapper',
     'example_H2O2_TS', 'example_freqs', 'example_generic_constraints',
-    'example_partial_geometry_opt', 'freq', 'geometry', 'singlepoint', 'ts', 'md',
+    'example_partial_geometry_opt',
+    'freq', 'geometry', 'singlepoint', 'ts', 'md', 'cell_opt',
     'find_first_job', 'select_max', 'select_min']

src/qmflows/packages/cp2k_mm.py

@@ -16,6 +16,7 @@
 from os.path import join, abspath
 from typing import Union, Any, ClassVar, Dict, Type
 
+import numpy as np
 from scm import plams
 
 from .packages import Result, parse_output_warnings, load_properties
@@ -200,13 +201,29 @@ def _parse_periodic(s: Settings, key: str, value: Any, mol: plams.Molecule) -> N
         if value.lower() == 'none':
             ewald.ewald_type = value
         elif not ewald.ewald_type:
-            ewald.ewald_type = 'EWALD'
+            ewald.ewald_type = 'SPME'
+
+    @staticmethod
+    def _parse_gmax(s: Settings, key: str, value: Any, mol: plams.Molecule) -> None:
+        """Set the keyword for the ``gmax`` parameter."""
+        ar = np.asarray(value)
+        if not issubclass(ar.dtype.type, np.str_):
+            ar = ar.astype(np.int64, copy=False, casting='same_kind')
+
+        ewald = s.specific.cp2k.force_eval.mm.poisson.ewald
+        if ar.ndim == 0:
+            ewald.gmax = str(ar.item())
+        elif ar.ndim == 1:
+            ewald.gmax = " ".join(str(i) for i in ar)
+        else:
+            raise RuntimeError(f"gmax:{value!r}\nformat not recognized")
 
 
 CP2KMM.SPECIAL_FUNCS = {
     'psf': CP2KMM._parse_psf,
     'prm': CP2KMM._parse_prm,
-    'periodic': CP2KMM._parse_periodic
+    'periodic': CP2KMM._parse_periodic,
+    'gmax': CP2KMM._parse_gmax,
 }
 
 for k in CP2K_KEYS_ALIAS:

src/qmflows/templates/__init__.py

@@ -1,5 +1,5 @@
 """Various Settings templates used throughout QMFlows."""
 
-from .templates import (freq, geometry, singlepoint, ts, md)
+from .templates import (freq, geometry, singlepoint, ts, md, cell_opt)
 
-__all__ = ['freq', 'geometry', 'singlepoint', 'ts', 'md']
+__all__ = ['freq', 'geometry', 'singlepoint', 'ts', 'md', 'cell_opt']

src/qmflows/templates/templates.py

@@ -1,6 +1,6 @@
 """Default options to call some quantum packages."""
 
-__all__ = ['freq', 'geometry', 'singlepoint', 'ts', 'md']
+__all__ = ['freq', 'geometry', 'singlepoint', 'ts', 'md', 'cell_opt']
 
 import yaml
 
@@ -348,3 +348,59 @@
             project: cp2k
             run_type: MD
 """, Loader=yaml.FullLoader))
+
+#: Templates for cell optimization calculations.
+cell_opt = Settings(yaml.load("""
+specific:
+    adf: Null
+
+    ams: Null
+
+    dftb: Null
+
+    cp2k: Null
+
+    orca: Null
+
+    cp2k_mm:
+        force_eval:
+            stress_tensor: analytical
+            method: FIST
+            mm:
+                print:
+                    ff_info low:
+                        spline_data: .FALSE.
+                        spline_info: .FALSE.
+                forcefield:
+                    ei_scale14: 1.0
+                    vdw_scale14: 1.0
+                    ignore_missing_critical_params: ''
+                    do_nonbonded: ''
+                    shift_cutoff: .TRUE.
+                    spline:
+                        r0_nb: 0.25
+                poisson:
+                    periodic: xyz
+                    ewald:
+                        ewald_type: spme
+                        o_spline: 4
+            subsys:
+                cell:
+                    periodic: xyz
+                topology:
+                    coord_file_format: 'OFF'
+                    center_coordinates:
+                        center_point: 0.0 0.0 0.0
+        motion:
+            cell_opt:
+                type: direct_cell_opt
+                optimizer: lbfgs
+                max_iter: 500
+            print:
+                cell low:
+                    filename: ''
+        global:
+            print_level: low
+            project: cp2k
+            run_type: cell_opt
+""", Loader=yaml.FullLoader))

test/test_cp2k_mm_mock.py

@@ -8,7 +8,7 @@
 from pytest_mock import MockFixture
 from scm.plams import Molecule
 
-from qmflows import Settings, cp2k_mm, singlepoint, geometry, freq, md
+from qmflows import Settings, cp2k_mm, singlepoint, geometry, freq, md, cell_opt
 from qmflows.utils import InitRestart
 from qmflows.packages.cp2k_mm import CP2KMM_Result
 from qmflows.test_utils import get_mm_settings, PATH, PATH_MOLECULES
@@ -139,3 +139,47 @@ def test_c2pk_md_mock(mocker) -> None:
     assertion.eq(result.status, 'successful')
 
     assertion.isfile(result.results['cp2k-1_1000.restart'])
+
+
+def test_c2pk_cell_opt_mock(mocker) -> None:
+    """Mock a call to CP2K."""
+    mol = Molecule(PATH / 'cspbbr3_3d.xyz')
+
+    s = Settings()
+    s.specific.cp2k += cell_opt.specific.cp2k_mm.copy()
+    s.specific.cp2k.motion.cell_opt.max_iter = 10
+    s.specific.cp2k.motion.print['forces low'].filename = ''
+
+    s.gmax = [22, 22, 22]
+    s.cell_parameters = [25.452, 35.995, 24.452]
+    s.charge = {
+        'param': 'charge',
+        'Cs': 0.2,
+        'Pb': 0.4,
+        'Br': -0.2,
+    }
+    s.lennard_jones = {
+        'param': ('sigma', 'epsilon'),
+        'unit': ('nm', 'kjmol'),
+        'Cs Cs': (0.585, 1),
+        'Cs Pb': (0.510, 1),
+        'Br Se': (0.385, 1),
+        'Pb Pb': (0.598, 1),
+        'Br Pb': (0.290, 1),
+        'Br Br': (0.426, 1),
+    }
+
+    job = cp2k_mm(s, mol)
+    run_mocked = mock_runner(mocker, settings=s, jobname="cp2k_mm_cell_opt")
+    result = run_mocked(job)
+    assertion.eq(result.status, 'successful')
+
+    ref_volume = np.load(PATH / 'volume.npy')
+    ref_coordinates = np.load(PATH / 'coordinates.npy')
+    ref_forces = np.load(PATH / 'forces.npy')
+    ref_lattice = np.load(PATH / 'lattice.npy')
+
+    np.testing.assert_allclose(result.volume, ref_volume)
+    np.testing.assert_allclose(result.coordinates, ref_coordinates)
+    np.testing.assert_allclose(result.forces, ref_forces)
+    np.testing.assert_allclose(result.lattice, ref_lattice)

test/test_files/cp2k_mm_special_keyword.yaml

@@ -50,11 +50,12 @@ cp2k:
         parmtype: CHM
       poisson:
         ewald:
-          ewald_type: none
-        periodic: none
+          ewald_type: SPME
+          gmax: '22 22 22'
+        periodic: xyz
     subsys:
       cell:
-        periodic: none
+        periodic: xyz
       kind C2O3:
         element: C
       kind C331: