Basis extrapolation

README.md

@@ -2,7 +2,7 @@
 
 EmbASI is a python-based wrapper for QM/QM embedding simulations using the ASI (Atomic Simulation Interface). The wrapper handles communication between QM code library calls through C-based callbacks, and imports/exports relevant matices such as density matrices and hamiltonians required for embedding schemes such as Projection-Based Embedding (PbE). Atomic information is communicated through Atomic Simulation Environment (ASE) atoms objects.
 
-# Supported QM Pakage(s)
+# Supported QM Package(s)
 
 - FHI-aims
 

embasi/embedding.py

@@ -321,7 +321,7 @@ def run(self):
 
 
 class ProjectionEmbedding(EmbeddingBase):
-    """Implementation of Projeciton-Based Embedding with ASI communication
+    """Implementation of Projection-Based Embedding with ASI communication
 
     A class controlling the interaction between two subsystems calculated at
     two levels of theory (low-level and high-level) with the

embasi/workflows/extrapolation.py

@@ -0,0 +1,116 @@
+import os, typing
+from embasi.embedding import ProjectionEmbedding
+
+#os.environ["ASI_LIB_PATH"] = "/home/dchen/Software/FHIaims/_build_lib/libaims.250711.scalapack.mpi.so"
+
+class Extrapolation:
+    """
+
+    Provides an estimated value based on the assumption that the
+    existing trend would continue. Allows us to predict energy at
+    the CBS Limit by using results from smaller basis sets, while
+    maintaining a high degree of accuracy as well as a more efficient
+    method which requires less computing power as well as less time.
+
+    Parameters
+    ----------
+    file1: str
+        Name of file which contains the first basis set
+    file2: str
+        Name of file which contains the second basis set
+    path: str
+        Name of directory which contain the basis set
+    atom: ASE Atom
+        Object of an atom, which contains information about the atom
+    embed_mask: list[int] or int
+        Assigns either the first in atoms to region 1, or an index of
+        int values 1 and 2 to each embedding layer. WARNING: The atoms
+        object will be reordered such that embedding layer 1 appear
+        first
+    calc_ll: ASE FileIOCalculator
+        Calculator object for layer 1
+    calc_hl: ASE FileIOCalculator
+        Calculator object for layer 2
+    asi_path: str
+        Name of directory where ASI (Atomic Simulation Interface) is installed
+    projection1_param: dict
+        Additional parameters for the first projection
+    projection2_param: dict
+        Additional parameters for the second projection
+    d: float
+        Value used for the formula E(∞)
+    alpha: float
+        Value used for the formula E(∞)
+
+    """
+    def __init__(self, file1, file2, path, atom, embed_mask, calc_ll, calc_hl, asi_path, projection1_param = {} , projection2_param = {} , d=2.85, alpha = 4.49):
+        self.asi_path:str = asi_path
+        os.environ["ASI_LIB_PATH"] = self.asi_path
+        self.file1:str= file1
+        self.file2:str = file2
+        self.path:str = path
+        self.atom:object = atom
+        self.embed_mask: typing.List[int] = embed_mask
+        self.calc_ll = calc_ll
+        self.calc_hl = calc_hl
+        self.mu_val: float = 1.e+6
+        self.options: typing.List[str] = [file1,file2]
+        self.results = []
+        self.projection1_param = projection1_param
+        self.projection2_param = projection2_param
+        self.d = d
+        self.alpha = alpha
+
+    def checkInParam(self, item: str, default: str, cycle: int) -> str:
+        if cycle == 0:
+            app_Dict = self.projection1_param.items()
+        else:
+            app_Dict = self.projection2_param.items()
+
+        for key, val in app_Dict:
+            if key == item:
+                return val
+
+        return default
+
+    @property
+    def extrapolate(self) -> float:
+        try:
+            conv1 = int(self.file1)
+            conv2 = int(self.file2)
+
+            if conv1 > conv2:
+                pass
+            else:
+                raise ValueError(
+                    "File1 should be bigger than File2."
+                )
+        except:
+            raise TypeError(
+                "File1: int\nFile2: int"
+            )
+
+        for index, (item, projection_parameters) in enumerate(zip(self.options, [self.projection1_param, self.projection2_param])):
+            os.environ["AIMS_SPECIES_DIR"] = f"{self.path}{item}Z"
+
+            projection = ProjectionEmbedding(
+                self.atom,
+                embed_mask=self.embed_mask,
+                calc_base_hl=self.calc_hl,
+                calc_base_ll=self.calc_ll,
+                mu_val=self.mu_val,
+                total_energy_corr= self.checkInParam("total_energy_corr","1storder", index),
+                localisation = self.checkInParam("localisation","SPADE", index),
+                projection = self.checkInParam("projection","level-shift", index)
+            )
+
+            for key, val in projection_parameters.items():
+                setattr(projection, key, val)
+
+            projection.run()
+
+            energy = projection.DFT_AinB_total_energy
+
+            self.results.append(energy)
+
+        return (self.results[0]*(int(self.file1)+self.d)**self.alpha - self.results[1]*(int(self.file2)+self.d)**self.alpha)/((int(self.file1) + self.d)**self.alpha-(int(self.file2)+self.d)**self.alpha)

examples/FHIaims/basis_extrapolation/extrapolation_example.py

@@ -0,0 +1,74 @@
+from embasi.workflows.extrapolation import Extrapolation
+from ase.data.s22 import s26,create_s22_system
+from ase.calculators.aims import Aims, AimsProfile
+'''
+A example of running an FHIaims QM/MM embedding simulation of the extrapolation
+of a methanol and a methanol dimer to calculate the end dissociation energy.
+'''
+
+
+# os.environ["ASI_LIB_PATH"] = "/home/dchen/Software/FHIaims/_build_lib/libaims.250711.scalapack.mpi.so"
+# os.environ['AIMS_SPECIES_DIR'] = "/home/dchen/Software/FHIaims/species_defaults/NAO-VCC-nZ/NAO-VCC-
+# The 'AIMS_SPECIES_DIR' and 'ASI_LIB_PATH' are passed into the class to allow for easier retrieval.
+
+calc_ll = Aims(
+    xc='PBE',
+    profile=AimsProfile(command="asi-doesnt-need-command"),
+    KS_method="parallel",
+    RI_method="LVL",
+    collect_eigenvectors=True,
+    density_update_method='density_matrix', # for DM export
+    atomic_solver_xc="PBE",
+    compute_kinetic=True,
+    override_initial_charge_check=True,
+    override_illconditioning=True
+)
+
+calc_hl = Aims(
+    xc='PBE',
+    profile=AimsProfile(command="asi-doesnt-need-command"),
+    KS_method="parallel",
+    RI_method="LVL",
+    collect_eigenvectors=True,
+    density_update_method='density_matrix',  # for DM export
+    atomic_solver_xc="PBE",
+    compute_kinetic=True,
+    override_initial_charge_check=True,
+    override_illconditioning=True
+)
+
+# Import dimer from s26 test set
+methanol_dimer_idx = s26[22]
+methanol_dimer = create_s22_system(methanol_dimer_idx)
+
+# Set-up calculator parameters (similar to FHIaims Calculator for
+# ASE) for low-level and high-level calculations. Below are the
+# absolute minimum parameters required for normal operation.
+
+test = Extrapolation(
+    file1= "3",
+    file2 = "2",
+    path = "/home/dchen/Software/FHIaims/species_defaults/NAO-VCC-nZ/NAO-VCC-",
+    asi_path = "/home/dchen/Software/FHIaims/_build_lib/libaims.250711.scalapack.mpi.so",
+    atom = methanol_dimer,
+    embed_mask = [2,1,2,2,2,1,2,1,2,2,2,1],
+    calc_ll = calc_ll,
+    calc_hl = calc_hl,
+)
+
+test2 = Extrapolation(
+    file1= "3",
+    file2 = "2",
+    path = "/home/dchen/Software/FHIaims/species_defaults/NAO-VCC-nZ/NAO-VCC-",
+    asi_path = "/home/dchen/Software/FHIaims/_build_lib/libaims.250711.scalapack.mpi.so",
+    atom = methanol_dimer[:6],
+    embed_mask = [2,1,2,2,2,1],
+    calc_ll = calc_ll,
+    calc_hl = calc_hl
+)
+
+energy1 = test.extrapolate  # Returns the extrapolated energy for methanol
+energy2 = test2.extrapolate  # Returns the extrapolated energy for the methanol dimer
+
+# Calculate dimer dissociation energy in the usual way
+print(energy1 - 2 * energy2)

tests/FHI-aims/projection_embedding/serial/basis_extrapolation/test_extrapolation.py

@@ -0,0 +1,17 @@
+import pytest
+import os
+import sys
+from pytest_regressions.data_regression import DataRegressionFixture
+from pytest_regressions.num_regression import NumericRegressionFixture
+
+sys.path.append(os.getcwd() + "/../../../")
+print(sys.path)
+from fhiaims_extrapolate_test_generator import FHIaims_projection_embedding_extrapolate_test
+
+@pytest.mark.mpi
+def test_basis_extrapolation_serial(num_regression: NumericRegressionFixture, tmp_path):
+    test = FHIaims_projection_embedding_extrapolate_test(test_dir=tmp_path)
+    test.run_test()
+    test_vals = test.output_ref_values()
+    num_regression.check(test_vals, default_tolerance=dict(atol=1e-6))
+

tests/FHI-aims/projection_embedding/serial/basis_extrapolation/test_extrapolation/test_basis_extrapolation_serial.csv

@@ -0,0 +1,2 @@
+,Methanol, Methanol Dimer, Dissociation energy
+0,-6300.0137534067026,-3149.8960985437543,-0.2215563191935091

tests/fhiaims_extrapolate_test_generator.py

@@ -0,0 +1,132 @@
+from ase.data.s22 import s22, s26, create_s22_system
+from ase.calculators.aims import Aims, AimsProfile
+from embasi.workflows.extrapolation import Extrapolation
+import os
+
+
+class FHIaims_projection_embedding_extrapolate_test():
+    def __init__(self, localisation="SPADE", projection="level-shift", total_energy_corr="1storder", parallel=False, gc=False,
+                 calc_ll_param_dict=None, calc_hl_param_dict=None, test_dir=None, basis_dir1="3", basis_dir2="2"):
+        # Set-up calculator parameters (similar to FHIaims Calculator for
+        # ASE) for low-level and high-level calculations. Below are the
+        # absolute minimum parameters required for normal operation.
+
+        if calc_hl_param_dict is None:
+            calc_hl_param_dict = {}
+        if calc_ll_param_dict is None:
+            calc_ll_param_dict = {}
+
+        self.calc_ll = Aims(
+            xc='PBE',
+            profile=AimsProfile(command="asi-doesnt-need-command"),
+            KS_method="parallel",
+            RI_method="LVL",
+            collect_eigenvectors=True,
+            density_update_method='density_matrix',  # for DM export
+            atomic_solver_xc="PBE",
+            compute_kinetic=True,
+            override_initial_charge_check=True,
+        )
+
+        for key, val in calc_ll_param_dict.items():
+            self.calc_ll.parameters[key] = val
+
+        self.calc_hl = Aims(
+            xc='PBE',
+            profile=AimsProfile(command="asi-doesnt-need-command"),
+            KS_method="parallel",
+            RI_method="LVL",
+            collect_eigenvectors=True,
+            density_update_method='density_matrix',  # for DM export
+            atomic_solver_xc="PBE",
+            compute_kinetic=True,
+            override_initial_charge_check=True,
+        )
+        for key, val in calc_hl_param_dict.items():
+            self.calc_hl.parameters[key] = val
+
+        self.localisation = localisation
+        self.projection = projection
+        self.total_energy_corr = total_energy_corr
+        self.parallel = False
+        self.gc = False
+        self.test_dir = test_dir
+        self.basis_dir1 = basis_dir1
+        self.basis_dir2 = basis_dir2
+        #os.environ["AIMS_SPECIES_DIR"] = os.environ["AIMS_ROOT_DIR"] + "/species_defaults/" + basis_dir
+        self.energy = None
+        self.energy2 = None
+        self.diff_energy = None
+
+    def run_test(self):
+
+        # Creates a S22 system s26[22], with 22 being the ID from ASE
+        methanol_dimer = create_s22_system(s26[22])
+        energy = Extrapolation(
+            file1=self.basis_dir1,
+            file2=self.basis_dir2,
+            path=os.environ["AIMS_ROOT_DIR"] + "/species_defaults/NAO-VCC-nZ/NAO-VCC-",
+            asi_path=os.environ["AIMS_LIB_PATH"],
+            atom=methanol_dimer,
+            embed_mask=[2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 1],
+            calc_ll=self.calc_ll,
+            calc_hl=self.calc_hl,
+            projection1_param={
+                "localisation": self.localisation,
+                "projection": self.projection,
+                "total_energy_corr": self.total_energy_corr,
+                "parallel" : self.parallel,
+                "gc" : self.gc,
+            },
+            projection2_param={
+                "localisation": self.localisation,
+                "projection": self.projection,
+                "total_energy_corr": self.total_energy_corr,
+                "parallel": self.parallel,
+                "gc": self.gc,
+            }
+        )
+
+        energy2 = Extrapolation(
+            file1=self.basis_dir1,
+            file2=self.basis_dir2,
+            path=os.environ["AIMS_ROOT_DIR"] + "/species_defaults/NAO-VCC-nZ/NAO-VCC-",
+            asi_path=os.environ["AIMS_LIB_PATH"],
+            atom=methanol_dimer[:6],
+            embed_mask=[2,1,2,2,2,1],
+            calc_ll=self.calc_ll,
+            calc_hl=self.calc_hl,
+            projection1_param={
+                "localisation": self.localisation,
+                "projection": self.projection,
+                "total_energy_corr": self.total_energy_corr,
+                "parallel": self.parallel,
+                "gc": self.gc,
+            },
+            projection2_param={
+                "localisation": self.localisation,
+                "projection": self.projection,
+                "total_energy_corr": self.total_energy_corr,
+                "parallel": self.parallel,
+                "gc": self.gc,
+            }
+        )
+
+        energy = energy.extrapolate  # Finds the monomer energy of the methanol
+        self.energy = energy  # Stores the energy
+
+        energy2 = energy2.extrapolate  # Finds the monomer energy of the methanol dimer
+        self.energy2 = energy2
+
+        self.diff_energy = energy - 2 * energy2 # Calculates the dissociation energy: diss_energy = energy2 - 2(energy1)
+
+    def output_ref_values(self):
+        return {
+            "Methanol Dimer": self.energy2,
+            "Methanol": self.energy,
+            "Dissociation Energy": self.diff_energy,
+        }
+
+
+
+

tests/fhiaims_meoh_test_generator.py

@@ -24,7 +24,7 @@ def __init__(self, localisation="SPADE", projection="huzinaga", parallel=False,
                        )
 
         for key, val in calc_ll_param_dict.items():
-            self.calculator_ll.parameters[key]=val
+            self.calc_ll.parameters[key]=val
 
         self.calc_hl = Aims(xc='PBE', profile=AimsProfile(command="asi-doesnt-need-command"),
                        KS_method="parallel",
@@ -37,7 +37,7 @@ def __init__(self, localisation="SPADE", projection="huzinaga", parallel=False,
                        )
 
         for key, val in calc_hl_param_dict.items():
-            self.calculator_hl.parameters[key]=val
+            self.calc_hl.parameters[key]=val
 
         self.localisation = localisation
         self.projection = projection