Merge pull request #278 from SCM-NV/openshell

compute derivative couplings for multiple spin states

CHANGELOG.md

@@ -3,18 +3,16 @@
 # 0.10.3 (Unreleased)
 ## New
 * Template to create B3LYP computations (#269)
+* Add support for derivative couplings for system with more than one spin state (#275)
 
 ## Fixed
 * Fix distribution error (#272)
-
-
-# 0.10.2
-## Fixed
+* Fix molecular orbital error [in qmflows](https://github.com/SCM-NV/qmflows/pull/213) (#270)
 * multivalue settings issue (#260)
 * CP2K executable (#264)
 
 # 0.10.1
-## Added
+## New
 * Keywords to print eigenvalues and eigenvectors (#248)
 
 ## Fixed

nanoqm/integrals/nonAdiabaticCoupling.py

@@ -83,7 +83,8 @@ def calculate_couplings_levine(dt: float, w_jk: Matrix,
     Garrett A. Meek and Benjamin G. Levine.
     dx.doi.org/10.1021/jz5009449 | J. Phys. Chem. Lett. 2014, 5, 2351−2356
 
-    .. NOTE::
+    Notes
+    -----
         In numpy sinc is defined as sin(pi * x) / (pi * x)
 
 
@@ -155,10 +156,10 @@ def calculate_couplings_levine(dt: float, w_jk: Matrix,
 
 def correct_phases(overlaps: Tensor3D, mtx_phases: Matrix) -> np.ndarray:
     """Correct the phases for all the overlaps."""
-    nOverlaps = overlaps.shape[0]  # total number of overlap matrices
+    noverlaps = overlaps.shape[0]  # total number of overlap matrices
     dim = overlaps.shape[1]  # Size of the square matrix
 
-    for k in range(nOverlaps):
+    for k in range(noverlaps):
         # Extract phases
         phases_t0, phases_t1 = mtx_phases[k: k + 2]
         phases_t0 = phases_t0.reshape(dim, 1)

nanoqm/schedule/components.py

@@ -63,8 +63,7 @@ def calculate_mos(config: DictConfig) -> List[str]:
 
     Returns
     -------
-    list
-        path to nodes in the HDF5 file to MO energies and MO coefficients.
+        paths to the datasets in the HDF5 file containging both the MO energies and MO coefficients
 
     """
     # Read Cell parameters file
@@ -82,7 +81,7 @@ def calculate_mos(config: DictConfig) -> List[str]:
 
     # orbital type is either an empty string for restricted calculation
     # or alpha/beta for unrestricted calculations
-    orbital_type = config.orbital_type
+    orbitals_type = config.orbitals_type
 
     for j, gs in enumerate(config.geometries):
 
@@ -96,7 +95,7 @@ def calculate_mos(config: DictConfig) -> List[str]:
 
         # Path where the MOs will be store in the HDF5
         root = join(config.project_name, f'point_{k}',
-                    config.package_name, 'mo', orbital_type)
+                    config.package_name, 'mo', orbitals_type)
         dict_input["node_MOs"] = [
             join(
                 root, 'eigenvalues'), join(
@@ -132,7 +131,7 @@ def calculate_mos(config: DictConfig) -> List[str]:
 
             # Store the computation
             if config["compute_orbitals"]:
-                orbitals.append(store_MOs(config, dict_input, promise_qm))
+                orbitals.append(store_molecular_orbitals(config, dict_input, promise_qm))
             else:
                 orbitals.append(None)
             energies.append(store_enery(config, dict_input, promise_qm))
@@ -143,7 +142,7 @@ def calculate_mos(config: DictConfig) -> List[str]:
 
 
 @schedule
-def store_MOs(
+def store_molecular_orbitals(
         config: DictConfig, dict_input: DefaultDict[str, Any], promise_qm: PromisedObject) -> str:
     """Store the MOs in the HDF5.
 
@@ -180,7 +179,7 @@ def save_orbitals_in_hdf5(mos: OrbitalType, config: DictConfig, job_name: str) -
 
 
 def dump_orbitals_to_hdf5(
-        data: InfoMO, config: DictConfig, job_name: str, orbital_type: str = "") -> None:
+        data: InfoMO, config: DictConfig, job_name: str, orbitals_type: str = "") -> None:
     """Store the result in HDF5 format.
 
     Parameters
@@ -191,15 +190,14 @@ def dump_orbitals_to_hdf5(
         Dictionary with the job configuration
     job_name
         Name of the current job
-    orbital_type
+    orbitals_type
         Either an empty string for MO coming from a restricted job or alpha/beta
         for unrestricted MO calculation
     """
-    es = join("cp2k", "mo", orbital_type, "eigenvalues")
-    css = join("cp2k", "mo", orbital_type, "coefficients")
+    root = join("cp2k", "mo", orbitals_type)
 
-    for path, array in zip((es, css), (data.eigenvalues, data.eigenvectors)):
-        path_property = join(config.project_name, job_name, path)
+    for name, array in zip(("eigenvalues", "coefficients"), (data.eigenvalues, data.eigenvectors)):
+        path_property = join(config.project_name, job_name, root, name)
         store_arrays_in_hdf5(config.path_hdf5, path_property, array)
 
 
@@ -303,7 +301,7 @@ def create_point_folder(
     return folders
 
 
-def split_file_geometries(path_xyz: PathLike) -> Sequence[MolXYZ]:
+def split_file_geometries(path_xyz: PathLike) -> Sequence[str]:
     """Read a set of molecular geometries in xyz format."""
     # Read Cartesian Coordinates
     with open(path_xyz) as f:
@@ -318,9 +316,9 @@ def create_file_names(work_dir: PathLike, i: int) -> JobFiles:
     file_xyz = join(work_dir, f'coordinates_{i}.xyz')
     file_inp = join(work_dir, f'point_{i}.inp')
     file_out = join(work_dir, f'point_{i}.out')
-    file_MO = join(work_dir, f'mo_coeff_{i}.out')
+    file_mo = join(work_dir, f'mo_coeff_{i}.out')
 
-    return JobFiles(file_xyz, file_inp, file_out, file_MO)
+    return JobFiles(file_xyz, file_inp, file_out, file_mo)
 
 
 def adjust_cell_parameters(

nanoqm/schedule/scheduleCoupling.py

@@ -91,8 +91,9 @@ def lazy_couplings(
         fixed_phase_overlaps = correct_phases(overlaps, mtx_phases)
 
     # Write the overlaps in text format
-    logger.debug("Writing down the overlaps in ascii format")
-    write_overlaps_in_ascii(fixed_phase_overlaps)
+    if config.write_overlaps:
+        logger.debug("Writing down the overlaps in ascii format")
+        write_overlaps_in_ascii(fixed_phase_overlaps)
 
     # Compute the couplings using either the levine method
     # or the 3Points approximation
@@ -112,7 +113,7 @@ def lazy_couplings(
 
 def compute_the_fixed_phase_overlaps(
         paths_overlaps: List[str], path_hdf5: PathLike, project_name: str,
-        enumerate_from: int, nHOMO: int) -> Tuple[List[np.ndarray], np.ndarray]:
+        enumerate_from: int, nHOMO: int) -> Tuple[np.ndarray, np.ndarray]:
     """Fix the phase of the overlaps.
 
     First track the unavoided crossings between Molecular orbitals and
@@ -187,7 +188,7 @@ def calculate_couplings(config: DictConfig, i: int, fixed_phase_overlaps: Tensor
 
     # Path were the couplinp is store
     k = i + config.enumerate_from
-    path = join(config.project_name, f'coupling_{k}')
+    path = join(config.project_name, config.orbitals_type, f'coupling_{k}')
 
     # Skip the computation if the coupling is already done
     if is_data_in_hdf5(config.path_hdf5, path):
@@ -380,7 +381,7 @@ def single_machine_overlaps(
 
 def create_overlap_path(config: DictConfig, i: int) -> str:
     """Create the path inside the HDF5 where the overlap is going to be store."""
-    root = join(config.project_name, 'overlaps_{}'.format(
+    root = join(config.project_name, config.orbitals_type, 'overlaps_{}'.format(
         i + config.enumerate_from))
     return join(root, 'mtx_sji_t0')
 

nanoqm/workflows/initialization.py

@@ -23,7 +23,7 @@
 from os.path import join
 from pathlib import Path
 from subprocess import PIPE, Popen
-from typing import List, Sequence, Union
+from typing import List, Union
 
 import numpy as np
 import pkg_resources
@@ -147,7 +147,7 @@ def read_swaps(path_hdf5: PathLike, project_name: str) -> Matrix:
         raise RuntimeError(msg)
 
 
-def split_trajectory(path: PathLike, nblocks: int, pathOut: PathLike) -> Sequence[PathLike]:
+def split_trajectory(path: PathLike, nblocks: int, pathOut: PathLike) -> List[str]:
     """Split an XYZ trajectory in n Block and write them in a given path.
 
     Parameters

nanoqm/workflows/schemas.py

@@ -171,7 +171,7 @@ def merge(d1: Dict[str, Any], d2: Dict[str, Any]) -> Dict[str, Any]:
 
     # Empty string for restricted calculation or either alpha/beta
     # for unrestricted calculation
-    Optional("orbital_type", default=""): any_lambda(("", "alpha", "beta"))
+    Optional("orbitals_type", default=""): any_lambda(("", "alphas", "betas"))
 }
 
 #: Dict with input options to run a derivate coupling workflow

nanoqm/workflows/workflow_coupling.py

@@ -34,7 +34,7 @@
 logger = logging.getLogger(__name__)
 
 
-def workflow_derivative_couplings(config: DictConfig) -> Tuple[List[PathLike], List[PathLike]]:
+def workflow_derivative_couplings(config: DictConfig) -> Tuple[List[str], List[str]]:
     """Compute the derivative couplings for a molecular dynamic trajectory."""
     # Dictionary containing the general configuration
     config.update(initialize(config))

test/test_coupling.py

@@ -1,19 +1,35 @@
 """Test the derivative coupling calculation."""
 import os
 import shutil
+from typing import Sequence
 
 import numpy as np
 
-from nanoqm.common import is_data_in_hdf5, retrieve_hdf5_data
+from nanoqm.common import DictConfig, is_data_in_hdf5, retrieve_hdf5_data
 from nanoqm.workflows.input_validation import process_input
 from nanoqm.workflows.workflow_coupling import workflow_derivative_couplings
 
 from .utilsTest import PATH_TEST, remove_files
 
 
 def test_fast_couplings(tmp_path):
+    """Check the derivative couplings workflow"""
+    run_derivative_coupling(tmp_path, 'input_fast_test_derivative_couplings.yml')
+
+
+def test_unrestricted_alphas(tmp_path):
+    """Test the derivative coupling for the alphas spin orbitals."""
+    run_derivative_coupling(tmp_path, 'input_couplings_alphas.yml')
+
+
+def test_unrestricted_betas(tmp_path):
+    """Test the derivative coupling for the alphas spin orbitals."""
+    run_derivative_coupling(tmp_path, 'input_couplings_betas.yml')
+
+
+def run_derivative_coupling(tmp_path: str, input_file: str) -> None:
     """Check that the couplings run."""
-    path_input = PATH_TEST / 'input_fast_test_derivative_couplings.yml'
+    path_input = PATH_TEST / input_file
     config = process_input(path_input, 'derivative_couplings')
     config["scratch_path"] = tmp_path
     tmp_hdf5 = os.path.join(tmp_path, 'fast_couplings.hdf5')
@@ -27,10 +43,10 @@ def test_fast_couplings(tmp_path):
         remove_files()
 
 
-def check_couplings(config: dict, tmp_hdf5: str) -> None:
+def check_couplings(config: DictConfig, tmp_hdf5: str) -> None:
     """Check that the couplings have meaningful values."""
     def create_paths(keyword: str) -> list:
-        return [f'{config.project_name}/{keyword}_{x}'
+        return [os.path.join(config.project_name, config.orbitals_type, f'{keyword}_{x}')
                 for x in range(len(config.geometries) - 1)]
     overlaps = create_paths('overlaps')
     couplings = create_paths('coupling')
@@ -51,7 +67,7 @@ def create_paths(keyword: str) -> list:
     assert (not np.all(np.isnan(tensor_couplings)))
 
 
-def check_hamiltonians(hamiltonians: str) -> None:
+def check_hamiltonians(hamiltonians: Sequence[str]) -> None:
     """Check that the hamiltonians were written correctly."""
     energies = np.stack([np.diag(np.loadtxt(ts[1])) for ts in hamiltonians])
     couplings = np.stack([np.loadtxt(ts[0]) for ts in hamiltonians])