diff --git a/README.md b/README.md index e108a23..952fb6f 100644 --- a/README.md +++ b/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 diff --git a/embasi/embedding.py b/embasi/embedding.py index c6d0acf..3fbd06e 100644 --- a/embasi/embedding.py +++ b/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 diff --git a/embasi/workflows/extrapolation.py b/embasi/workflows/extrapolation.py new file mode 100644 index 0000000..e60d5da --- /dev/null +++ b/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) diff --git a/examples/FHIaims/basis_extrapolation/extrapolation_example.py b/examples/FHIaims/basis_extrapolation/extrapolation_example.py new file mode 100644 index 0000000..4f914ef --- /dev/null +++ b/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) diff --git a/tests/FHI-aims/projection_embedding/serial/basis_extrapolation/test_extrapolation.py b/tests/FHI-aims/projection_embedding/serial/basis_extrapolation/test_extrapolation.py new file mode 100644 index 0000000..dfe4927 --- /dev/null +++ b/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)) + diff --git a/tests/FHI-aims/projection_embedding/serial/basis_extrapolation/test_extrapolation/test_basis_extrapolation_serial.csv b/tests/FHI-aims/projection_embedding/serial/basis_extrapolation/test_extrapolation/test_basis_extrapolation_serial.csv new file mode 100644 index 0000000..b4e64e5 --- /dev/null +++ b/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 diff --git a/tests/fhiaims_extrapolate_test_generator.py b/tests/fhiaims_extrapolate_test_generator.py new file mode 100644 index 0000000..b621419 --- /dev/null +++ b/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, + } + + + + diff --git a/tests/fhiaims_meoh_test_generator.py b/tests/fhiaims_meoh_test_generator.py index b1da9a4..2bc760a 100644 --- a/tests/fhiaims_meoh_test_generator.py +++ b/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