Coop coordination

CHANGELOG.md

@@ -2,7 +2,9 @@
 
 # 0.10.4 (Unrelease)
 ## New
+* Template to create B3LYP computations (#269)
 * Allow to compute both alphas/betas derivative couplings simultaneusly (#275)
+* Add nanoCAT dependency (#280)
 
 ## Changed
 * Do not remove the CP2K log files by default

nanoqm/workflows/schemas.py

@@ -303,7 +303,10 @@ def merge(d1: Dict[str, Any], d2: Dict[str, Any]) -> Dict[str, Any]:
 #: Input for a Crystal Orbital Overlap Population calculation
 dict_coop = {
     # List of the two elements to calculate the COOP for
-    "coop_elements": list}
+    "coop_elements": list,
+    # Coordination number to be considered for each of the two elements
+    Optional("elements_coordination", default=None): Or(None, list)
+}
 
 
 dict_merged_single_points = merge(dict_general_options, dict_single_points)

nanoqm/workflows/workflow_coop.py

@@ -1,4 +1,7 @@
 """Crystal Orbital Overlap Population calculation.
+The COOP is calculated between two selected elements.
+For each element, a specific coordination number can optionally be selected 
+by using the "elements_coordination" key in the yaml input.
 
 Index
 -----
@@ -10,12 +13,16 @@
 __all__ = ['workflow_crystal_orbital_overlap_population']
 
 import logging
-from typing import List, Tuple
+from typing import List, Tuple, Union
 
 import numpy as np
 
+from scm.plams import Molecule
+
 from qmflows.parsers.xyzParser import readXYZ
 
+from nanoCAT.recipes import coordination_number
+
 from ..common import (DictConfig, MolXYZ, h2ev,
                       number_spherical_functions_per_atom, retrieve_hdf5_data)
 from ..integrals.multipole_matrices import compute_matrix_multipole
@@ -44,11 +51,16 @@ def workflow_crystal_orbital_overlap_population(config: DictConfig):
     # Converting the xyz-file to a mol-file
     mol = readXYZ(config.path_traj_xyz)
 
+    if config.elements_coordination is not None:
+        geometry = Molecule(config.path_traj_xyz)
+    else:
+        geometry = None
+
     # Computing the indices of the atomic orbitals of the two selected
     # elements, and the overlap matrix that contains only elements related to
     # the two elements
     el_1_orbital_ind, el_2_orbital_ind, overlap_reduced = compute_overlap_and_atomic_orbitals(
-        mol, config)
+        mol, config, geometry)
 
     # Compute the crystal orbital overlap population between the two selected
     # elements
@@ -83,7 +95,8 @@ def get_eigenvalues_coefficients(config: DictConfig) -> Tuple[np.ndarray, np.nda
 
 
 def compute_overlap_and_atomic_orbitals(
-        mol: MolXYZ, config: DictConfig) -> Tuple[List[int], List[int], List[int]]:
+        mol: MolXYZ, config: DictConfig,
+        geometry: Union[Molecule, None]) -> Tuple[List[int], List[int], List[int]]:
     """Compute the indices of the atomic orbitals of the two selected elements.
 
     Computes the overlap matrix, containing only the elements related to those two elements.
@@ -102,8 +115,17 @@ def compute_overlap_and_atomic_orbitals(
     element_1 = config["coop_elements"][0]
     element_2 = config["coop_elements"][1]
 
-    element_1_index = [i for i, s in enumerate(mol) if element_1.lower() in s]
-    element_2_index = [i for i, s in enumerate(mol) if element_2.lower() in s]
+    if config["elements_coordination"] is None:
+        element_1_index = [i for i, s in enumerate(mol) if element_1.lower() in s]
+        element_2_index = [i for i, s in enumerate(mol) if element_2.lower() in s]
+
+    # For the two selected elements, only the indices corresponding to a given coordination number
+    else:
+        coord = coordination_number(geometry)
+        cn_1 = config["elements_coordination"][0]
+        cn_2 = config["elements_coordination"][1]
+        element_1_index = [i - 1 for i in coord[element_1][cn_1]]
+        element_2_index = [i - 1 for i in coord[element_2][cn_2]]
 
     # Making a list of the indices of the atomic orbitals for each of the two
     # elements

setup.py

@@ -151,7 +151,8 @@ def build_extensions(self):
         'h5py', 'mendeleev', 'more-itertools', 'noodles==0.3.3', 'numpy', 'pybind11>=2.2.4',
         'scipy', 'schema', 'pyyaml>=5.1',
         'plams@git+https://github.com/SCM-NV/PLAMS@master',
-        'qmflows@git+https://github.com/SCM-NV/qmflows@master'
+        'qmflows@git+https://github.com/SCM-NV/qmflows@master',
+        'nano-CAT@git+https://github.com/nlesc-nano/nano-CAT@master'
     ],
     cmdclass={'build_ext': BuildExt},
     ext_modules=[ext_pybind],

test/test_distribute.py

@@ -1,10 +1,12 @@
 """Test the distribution script."""
-from qmflows.type_hints import PathLike
-from pathlib import Path
-from subprocess import (PIPE, Popen)
 import fnmatch
-import shutil
 import os
+import re
+import shutil
+from pathlib import Path
+from subprocess import PIPE, Popen
+
+from qmflows.type_hints import PathLike
 
 
 def test_distribute(tmp_path: PathLike) -> None:
@@ -20,7 +22,8 @@ def call_distribute(tmp_path: PathLike, cmd: str) -> None:
     try:
         p = Popen(cmd, stdin=PIPE, stdout=PIPE, stderr=PIPE, shell=True)
         _, err = p.communicate()
-        if err:
+        error = re.search("error", err.decode(), re.IGNORECASE)
+        if error is not None:
             raise RuntimeError(err.decode())
         check_scripts()
     finally:

test/test_files/input_test_coop_coordination.yml

@@ -0,0 +1,26 @@
+workflow:
+  coop_calculation
+project_name: Cd33Se33
+
+active_space: [50, 50]
+path_hdf5: "test/test_files/Cd33Se33.hdf5"
+path_traj_xyz: "test/test_files/Cd33Se33.xyz"
+scratch_path: "/tmp/COOP"
+
+coop_elements: ["Cd", "Se"]
+elements_coordination: [4, 3]
+
+cp2k_general_settings:
+  basis:  "DZVP-MOLOPT-SR-GTH"
+  potential: "GTH-PBE"
+  cell_parameters: 20.0
+  periodic: none
+
+  cp2k_settings_main:
+    specific:
+      template: pbe_main
+
+  cp2k_settings_guess:
+    specific:
+      template:
+        pbe_guess

test/test_run_workflow.py

@@ -1,6 +1,7 @@
 """Test the CLI to run a workflow."""
 
 import os
+import re
 import shutil
 from pathlib import Path
 from subprocess import PIPE, Popen
@@ -37,7 +38,8 @@ def test_run_workflow(tmp_path):
     try:
         p = Popen(cmd, stdin=PIPE, stdout=PIPE, stderr=PIPE, shell=True)
         out, err = p.communicate()
-        if err:
+        error = re.search("error", err.decode(), re.IGNORECASE)
+        if error is not None:
             print("output: ", out)
             print("err: ", err)
             raise RuntimeError(err.decode())

test/test_workflow_coop.py

@@ -29,3 +29,21 @@ def test_workflow_coop(tmp_path: PathLike) -> None:
         print("scratch_path: ", tmp_path)
         print("Unexpected error:", sys.exc_info()[0])
         raise
+
+
+def test_workflow_coop_coordination(tmp_path: PathLike) -> None:
+    """Test the Crystal Orbital Overlap Population workflow."""
+    file_path = PATH_TEST / 'input_test_coop_coordination.yml'
+    config = process_input(file_path, 'coop_calculation')
+
+    # create scratch path
+    shutil.copy(config.path_hdf5, tmp_path)
+    config.path_hdf5 = join(tmp_path, "Cd33Se33.hdf5")
+    config.workdir = tmp_path
+    try:
+        workflow_crystal_orbital_overlap_population(config)
+        os.remove("COOP.txt")
+    except BaseException:
+        print("scratch_path: ", tmp_path)
+        print("Unexpected error:", sys.exc_info()[0])
+        raise