From 4542d745f168fba458301bc9926a843880d5f8b0 Mon Sep 17 00:00:00 2001 From: sovietdevil <2558390548@qq.com> Date: Fri, 26 Jul 2024 15:54:52 +0200 Subject: [PATCH 1/7] Add ShiftML version 2 example, and correct the variable name in readme --- README.md | 2 +- src/shiftml/ase/calculator.py | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index fd60cec..3ad8126 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ from shiftml.ase import ShiftML frame = bulk("C", "diamond", a=3.566) calculator = ShiftML("ShiftML1.0") -cs_iso = calc.get_cs_iso(frame) +cs_iso = calculator.get_cs_iso(frame) print(cs_iso) diff --git a/src/shiftml/ase/calculator.py b/src/shiftml/ase/calculator.py index 791f223..c2ad35d 100644 --- a/src/shiftml/ase/calculator.py +++ b/src/shiftml/ase/calculator.py @@ -14,16 +14,19 @@ url_resolve = { "ShiftML1.0": "https://tinyurl.com/3xwec68f", "ShiftML1.1": "https://tinyurl.com/53ymkhvd", + "ShiftML2.0": "https://tinyurl.com/9v8ppnru", } resolve_outputs = { "ShiftML1.0": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, "ShiftML1.1": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, + "ShiftML2.0": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, } resolve_fitted_species = { "ShiftML1.0": set([1, 6, 7, 8, 16]), "ShiftML1.1": set([1, 6, 7, 8, 16]), + "ShiftML2.0": set([1, 6, 7, 8, 9, 11, 12, 15, 16, 17, 19, 20]), } From 56735307016f8f3326a23fcff0fa06027d4ad743 Mon Sep 17 00:00:00 2001 From: sovietdevil <2558390548@qq.com> Date: Mon, 29 Jul 2024 10:42:07 +0200 Subject: [PATCH 2/7] Add mean and standard deviation functions to the ShiftML2.0 calculator. --- src/shiftml/ase/calculator.py | 93 ++++++++++++++++--- tests/test_ase.py | 162 +++++++++++++++++++++++++++++++++- 2 files changed, 240 insertions(+), 15 deletions(-) diff --git a/src/shiftml/ase/calculator.py b/src/shiftml/ase/calculator.py index c2ad35d..02e3627 100644 --- a/src/shiftml/ase/calculator.py +++ b/src/shiftml/ase/calculator.py @@ -14,13 +14,17 @@ url_resolve = { "ShiftML1.0": "https://tinyurl.com/3xwec68f", "ShiftML1.1": "https://tinyurl.com/53ymkhvd", - "ShiftML2.0": "https://tinyurl.com/9v8ppnru", + "ShiftML2.0": "https://tinyurl.com/2mp8emsd", } resolve_outputs = { "ShiftML1.0": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, "ShiftML1.1": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, - "ShiftML2.0": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, + "ShiftML2.0": { + "mtt::cs_iso_mean": ModelOutput(quantity="", unit="ppm", per_atom=True), + "mtt::cs_iso_std": ModelOutput(quantity="", unit="ppm", per_atom=True), + "mtt::cs_iso_ensemble": ModelOutput(quantity="", unit="ppm", per_atom=True), + }, } resolve_fitted_species = { @@ -153,25 +157,86 @@ def __init__(self, model_version, force_download=False): raise e super().__init__(model_file) + self.model_version = model_version def get_cs_iso(self, atoms): """ Compute the shielding values for the given atoms object """ + if self.model_version == "ShiftML1.0" or self.model_version == "ShiftML1.1": + assert ( + "mtt::cs_iso" in self.outputs.keys() + ), "model does not support chemical shielding prediction" + + if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): + raise ValueError( + f"Model is fitted only for the following atomic numbers:\ + {self.fitted_species}. The atomic numbers in the atoms object are:\ + {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ + with only the fitted species." + ) - assert ( - "mtt::cs_iso" in self.outputs.keys() - ), "model does not support chemical shielding prediction" + out = self.run_model(atoms, self.outputs) + cs_iso = out["mtt::cs_iso"].block(0).values.detach().numpy() - if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): - raise ValueError( - f"Model is fitted only for the following atomic numbers:\ - {self.fitted_species}. The atomic numbers in the atoms object are:\ - {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ - with only the fitted species." - ) + elif self.model_version == "ShiftML2.0": + assert ( + "mtt::cs_iso_mean" in self.outputs.keys() + ), "model does not support chemical shielding prediction" + + if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): + raise ValueError( + f"Model is fitted only for the following atomic numbers:\ + {self.fitted_species}. The atomic numbers in the atoms object are:\ + {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ + with only the fitted species." + ) - out = self.run_model(atoms, self.outputs) - cs_iso = out["mtt::cs_iso"].block(0).values.detach().numpy() + out = self.run_model(atoms, self.outputs) + cs_iso = out["mtt::cs_iso_mean"].block(0).values.detach().numpy() return cs_iso + + def get_cs_iso_std(self, atoms): + if self.model_version == "ShiftML2.0": + assert ( + "mtt::cs_iso_std" in self.outputs.keys() + ), "model does not support chemical shielding prediction" + + if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): + raise ValueError( + f"Model is fitted only for the following atomic numbers:\ + {self.fitted_species}. The atomic numbers in the atoms object are:\ + {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ + with only the fitted species." + ) + + out = self.run_model(atoms, self.outputs) + cs_iso_std = out["mtt::cs_iso_std"].block(0).values.detach().numpy() + else: + raise RuntimeError("Version not supporting uncertainty quantification.") + + return cs_iso_std + + def get_cs_iso_ensemble(self, atoms): + if self.model_version == "ShiftML2.0": + assert ( + "mtt::cs_iso_ensemble" in self.outputs.keys() + ), "model does not support chemical shielding prediction" + + if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): + raise ValueError( + f"Model is fitted only for the following atomic numbers:\ + {self.fitted_species}. The atomic numbers in the atoms object are:\ + {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ + with only the fitted species." + ) + + out = self.run_model(atoms, self.outputs) + cs_iso_ensemble = ( + out["mtt::cs_iso_ensemble"].block(0).values.detach().numpy() + ) + else: + raise RuntimeError("Version not supporting uncertainty quantification.") + + return cs_iso_ensemble diff --git a/tests/test_ase.py b/tests/test_ase.py index 69f8ca3..482a944 100644 --- a/tests/test_ase.py +++ b/tests/test_ase.py @@ -7,6 +7,144 @@ from shiftml.ase import ShiftML expected_output = np.array([137.5415, 137.5415]) +expected_ensemble_v2 = np.array( + [ + [ + 114.8194808959961, + 113.47244262695312, + 117.47064208984375, + 115.61190795898438, + 130.88909912109375, + 131.42332458496094, + 120.96844482421875, + 115.95867919921875, + 116.98180389404297, + 135.3658447265625, + 120.45016479492188, + 123.00967407226562, + 137.23724365234375, + 129.23104858398438, + 131.00619506835938, + 130.82601928710938, + 121.90162658691406, + 120.66400909423828, + 109.59469604492188, + 118.66798400878906, + 126.18386840820312, + 124.9156494140625, + 120.90362548828125, + 106.26658630371094, + 128.32107543945312, + 125.82593536376953, + 121.3394775390625, + 127.37902069091797, + 122.92572784423828, + 126.26400756835938, + 112.87037658691406, + 112.48919677734375, + 126.00082397460938, + 109.98661804199219, + 110.7204818725586, + 107.30191040039062, + 113.85182189941406, + 110.24645233154297, + 133.27935791015625, + 126.40534973144531, + 133.42047119140625, + 112.2728271484375, + 126.27506256103516, + 117.58969116210938, + 119.17208099365234, + 121.65959167480469, + 115.62092590332031, + 118.12762451171875, + 119.478271484375, + 137.32974243164062, + 120.26103210449219, + 118.25013732910156, + 121.78120422363281, + 125.66693115234375, + 112.0889892578125, + 115.92691802978516, + 121.31621551513672, + 118.76759338378906, + 126.86924743652344, + 129.01571655273438, + 109.53144073486328, + 110.71353149414062, + 125.9607925415039, + 108.36444091796875, + ], + [ + 114.8194808959961, + 113.47244262695312, + 117.47064208984375, + 115.61190795898438, + 130.88909912109375, + 131.42332458496094, + 120.96844482421875, + 115.95867919921875, + 116.98180389404297, + 135.3658447265625, + 120.45016479492188, + 123.00967407226562, + 137.23724365234375, + 129.23104858398438, + 131.00619506835938, + 130.82601928710938, + 121.90162658691406, + 120.66400909423828, + 109.59469604492188, + 118.66798400878906, + 126.18386840820312, + 124.9156494140625, + 120.90362548828125, + 106.26658630371094, + 128.32107543945312, + 125.82593536376953, + 121.3394775390625, + 127.37902069091797, + 122.92572784423828, + 126.26400756835938, + 112.87037658691406, + 112.48919677734375, + 126.00082397460938, + 109.98661804199219, + 110.7204818725586, + 107.30191040039062, + 113.85182189941406, + 110.24645233154297, + 133.27935791015625, + 126.40534973144531, + 133.42047119140625, + 112.2728271484375, + 126.27506256103516, + 117.58969116210938, + 119.17208099365234, + 121.65959167480469, + 115.62092590332031, + 118.12762451171875, + 119.478271484375, + 137.32974243164062, + 120.26103210449219, + 118.25013732910156, + 121.78120422363281, + 125.66693115234375, + 112.0889892578125, + 115.92691802978516, + 121.31621551513672, + 118.76759338378906, + 126.86924743652344, + 129.01571655273438, + 109.53144073486328, + 110.71353149414062, + 125.9607925415039, + 108.36444091796875, + ], + ] +) +expected_mean_v2 = np.array([120.85137, 120.85137]) +expected_std_v2 = np.array([7.7993703, 7.7993703]) def test_shiftml1_regression(): @@ -62,7 +200,7 @@ def test_shiftml1_size_extensivity_test(): def test_shftml1_fail_invalid_species(): - """Test ShiftML1.o for non-fitted species""" + """Test ShiftML1.0 for non-fitted species""" frame = bulk("Si", "diamond", a=3.566) model = ShiftML("ShiftML1.0") @@ -73,3 +211,25 @@ def test_shftml1_fail_invalid_species(): assert "Model is fitted only for the following atomic numbers:" in str( exc_info.value ) + + +def test_shiftml2_regression_mean(): + """Regression test for the ShiftML2.0 model.""" + + frame = bulk("C", "diamond", a=3.566) + model = ShiftML("ShiftML2.0", force_download=True) + out_mean = model.get_cs_iso(frame) + out_std = model.get_cs_iso_std(frame) + out_ensemble = model.get_cs_iso_ensemble(frame) + + assert np.allclose( + out_mean.flatten(), expected_mean_v2 + ), "ShiftML2 failed regression mean test" + + assert np.allclose( + out_std.flatten(), expected_std_v2 + ), "ShiftML2 failed regression variance test" + + assert np.allclose( + out_ensemble.flatten(), expected_ensemble_v2 + ), "ShiftML2 failed regression ensemble test" From b91754f47b10b1140d6d81aa8633d6c7350529f4 Mon Sep 17 00:00:00 2001 From: sovietdevil <2558390548@qq.com> Date: Mon, 29 Jul 2024 11:23:49 +0200 Subject: [PATCH 3/7] Unify formats of different ShiftML versions, correct mistakes and simplify calculator functions. --- src/shiftml/ase/calculator.py | 97 ++++++++++++----------------------- tests/test_ase.py | 2 +- 2 files changed, 33 insertions(+), 66 deletions(-) diff --git a/src/shiftml/ase/calculator.py b/src/shiftml/ase/calculator.py index 02e3627..e091ff2 100644 --- a/src/shiftml/ase/calculator.py +++ b/src/shiftml/ase/calculator.py @@ -14,14 +14,14 @@ url_resolve = { "ShiftML1.0": "https://tinyurl.com/3xwec68f", "ShiftML1.1": "https://tinyurl.com/53ymkhvd", - "ShiftML2.0": "https://tinyurl.com/2mp8emsd", + "ShiftML2.0": "https://tinyurl.com/bdcp647w", } resolve_outputs = { "ShiftML1.0": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, "ShiftML1.1": {"mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True)}, "ShiftML2.0": { - "mtt::cs_iso_mean": ModelOutput(quantity="", unit="ppm", per_atom=True), + "mtt::cs_iso": ModelOutput(quantity="", unit="ppm", per_atom=True), "mtt::cs_iso_std": ModelOutput(quantity="", unit="ppm", per_atom=True), "mtt::cs_iso_ensemble": ModelOutput(quantity="", unit="ppm", per_atom=True), }, @@ -34,6 +34,16 @@ } +def is_fitted_on(atoms, fitted_species): + if not set(atoms.get_atomic_numbers()).issubset(fitted_species): + raise ValueError( + f"Model is fitted only for the following atomic numbers:\ + {fitted_species}. The atomic numbers in the atoms object are:\ + {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ + with only the fitted species." + ) + + class ShiftML(MetatensorCalculator): """ ShiftML calculator for ASE @@ -163,80 +173,37 @@ def get_cs_iso(self, atoms): """ Compute the shielding values for the given atoms object """ - if self.model_version == "ShiftML1.0" or self.model_version == "ShiftML1.1": - assert ( - "mtt::cs_iso" in self.outputs.keys() - ), "model does not support chemical shielding prediction" - - if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): - raise ValueError( - f"Model is fitted only for the following atomic numbers:\ - {self.fitted_species}. The atomic numbers in the atoms object are:\ - {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ - with only the fitted species." - ) - - out = self.run_model(atoms, self.outputs) - cs_iso = out["mtt::cs_iso"].block(0).values.detach().numpy() + assert ( + "mtt::cs_iso" in self.outputs.keys() + ), "model does not support chemical shielding prediction" - elif self.model_version == "ShiftML2.0": - assert ( - "mtt::cs_iso_mean" in self.outputs.keys() - ), "model does not support chemical shielding prediction" - - if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): - raise ValueError( - f"Model is fitted only for the following atomic numbers:\ - {self.fitted_species}. The atomic numbers in the atoms object are:\ - {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ - with only the fitted species." - ) + is_fitted_on(atoms, self.fitted_species) - out = self.run_model(atoms, self.outputs) - cs_iso = out["mtt::cs_iso_mean"].block(0).values.detach().numpy() + out = self.run_model(atoms, self.outputs) + cs_iso = out["mtt::cs_iso"].block(0).values.detach().numpy() return cs_iso def get_cs_iso_std(self, atoms): - if self.model_version == "ShiftML2.0": - assert ( - "mtt::cs_iso_std" in self.outputs.keys() - ), "model does not support chemical shielding prediction" - - if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): - raise ValueError( - f"Model is fitted only for the following atomic numbers:\ - {self.fitted_species}. The atomic numbers in the atoms object are:\ - {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ - with only the fitted species." - ) + assert ( + "mtt::cs_iso_std" in self.outputs.keys() + ), "model does not support chemical shielding prediction" - out = self.run_model(atoms, self.outputs) - cs_iso_std = out["mtt::cs_iso_std"].block(0).values.detach().numpy() - else: - raise RuntimeError("Version not supporting uncertainty quantification.") + is_fitted_on(atoms, self.fitted_species) + + out = self.run_model(atoms, self.outputs) + cs_iso_std = out["mtt::cs_iso_std"].block(0).values.detach().numpy() return cs_iso_std def get_cs_iso_ensemble(self, atoms): - if self.model_version == "ShiftML2.0": - assert ( - "mtt::cs_iso_ensemble" in self.outputs.keys() - ), "model does not support chemical shielding prediction" - - if not set(atoms.get_atomic_numbers()).issubset(self.fitted_species): - raise ValueError( - f"Model is fitted only for the following atomic numbers:\ - {self.fitted_species}. The atomic numbers in the atoms object are:\ - {set(atoms.get_atomic_numbers())}. Please provide an atoms object\ - with only the fitted species." - ) + assert ( + "mtt::cs_iso_ensemble" in self.outputs.keys() + ), "model does not support chemical shielding prediction" - out = self.run_model(atoms, self.outputs) - cs_iso_ensemble = ( - out["mtt::cs_iso_ensemble"].block(0).values.detach().numpy() - ) - else: - raise RuntimeError("Version not supporting uncertainty quantification.") + is_fitted_on(atoms, self.fitted_species) + + out = self.run_model(atoms, self.outputs) + cs_iso_ensemble = out["mtt::cs_iso_ensemble"].block(0).values.detach().numpy() return cs_iso_ensemble diff --git a/tests/test_ase.py b/tests/test_ase.py index 482a944..06c18e7 100644 --- a/tests/test_ase.py +++ b/tests/test_ase.py @@ -231,5 +231,5 @@ def test_shiftml2_regression_mean(): ), "ShiftML2 failed regression variance test" assert np.allclose( - out_ensemble.flatten(), expected_ensemble_v2 + out_ensemble.flatten(), expected_ensemble_v2.flatten() ), "ShiftML2 failed regression ensemble test" From a64ea167b681ed85427d18e7bc364e832a622cc6 Mon Sep 17 00:00:00 2001 From: sovietdevil <2558390548@qq.com> Date: Mon, 29 Jul 2024 11:41:38 +0200 Subject: [PATCH 4/7] Enlarge tolerance of regression check to avoid the error raised from numeric calculations. --- tests/test_ase.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_ase.py b/tests/test_ase.py index 06c18e7..db397fd 100644 --- a/tests/test_ase.py +++ b/tests/test_ase.py @@ -223,13 +223,13 @@ def test_shiftml2_regression_mean(): out_ensemble = model.get_cs_iso_ensemble(frame) assert np.allclose( - out_mean.flatten(), expected_mean_v2 + out_mean.flatten(), expected_mean_v2, atol=0.01 ), "ShiftML2 failed regression mean test" assert np.allclose( - out_std.flatten(), expected_std_v2 + out_std.flatten(), expected_std_v2, atol=0.01 ), "ShiftML2 failed regression variance test" assert np.allclose( - out_ensemble.flatten(), expected_ensemble_v2.flatten() + out_ensemble.flatten(), expected_ensemble_v2.flatten(), atol=0.01 ), "ShiftML2 failed regression ensemble test" From 07a84bbae6a521ffb0be02a383c74bd944e6bf75 Mon Sep 17 00:00:00 2001 From: sovietdevil <2558390548@qq.com> Date: Tue, 30 Jul 2024 11:12:39 +0200 Subject: [PATCH 5/7] Add structure prediction functions based on the ShiftML model --- src/shiftml/ase/calculator.py | 3 +- src/shiftml/csp/structure_pred.py | 102 ++++++++++++++++++++++++++++++ 2 files changed, 103 insertions(+), 2 deletions(-) create mode 100644 src/shiftml/csp/structure_pred.py diff --git a/src/shiftml/ase/calculator.py b/src/shiftml/ase/calculator.py index e091ff2..49eee95 100644 --- a/src/shiftml/ase/calculator.py +++ b/src/shiftml/ase/calculator.py @@ -13,7 +13,7 @@ url_resolve = { "ShiftML1.0": "https://tinyurl.com/3xwec68f", - "ShiftML1.1": "https://tinyurl.com/53ymkhvd", + "ShiftML1.1": "https://tinyurl.com/f237evr3", "ShiftML2.0": "https://tinyurl.com/bdcp647w", } @@ -167,7 +167,6 @@ def __init__(self, model_version, force_download=False): raise e super().__init__(model_file) - self.model_version = model_version def get_cs_iso(self, atoms): """ diff --git a/src/shiftml/csp/structure_pred.py b/src/shiftml/csp/structure_pred.py new file mode 100644 index 0000000..3a086f1 --- /dev/null +++ b/src/shiftml/csp/structure_pred.py @@ -0,0 +1,102 @@ +import numpy as np +import pandas as pd +from sklearn.metrics import root_mean_squared_error + +from shiftml.ase import ShiftML + + +def load_experimental(filename, skiprows=1): + """ + Function to import experimental .csv file and + convert to atomic indexes and chemical shieldings + + Parameters + ---------- + filname : str + name of the .csv file with path + skiprows : int + number of the rows to skip, default 1 + + Outputs + ------- + list_atom : np.array(List[str]) + list of the atomic indexes, in string array + list_cs: np.array(List[str]) + chemical sheilding for each corresponding index, in float array + """ + exp_data = pd.read_csv(filename, skiprows=skiprows) + atom_label = exp_data["Atom Label"] + atom_cs = exp_data["^(1)Hdelta(ppm)"] + list_atom = atom_label.to_list() + list_cs = atom_cs.array + return list_atom, list_cs + + +def extract_from_array(array, num_atoms, list_atom): + """ + Function to extract regression array from a given + symbol list of atoms and data exported from the model + """ + num_molecules = int(len(array) / num_atoms) + data_fit = array.reshape((-1, num_molecules))[:, 0] + X = [] + for atom_string in list_atom: + label_list = atom_string.split(",") + if len(label_list) == 1: + label = int(label_list[0]) + X.append(data_fit[label - 1]) + else: + X.append( + sum([data_fit[int(label_str) - 1] for label_str in label_list]) + / len(label_list) + ) + return np.array(X) + + +def structure_prediction( + model_version, frames, list_atom, list_cs, GIPAW_avail=True, cs_sym="CS" +): + """ + Function to select the suitable structures + based on a set of candidate structures, + given rmse of the linear regression results + + Parameters + ---------- + model_version : str + The version of the ShiftML model to use. Supported versions are + "ShiftML1.0", "ShiftML1.1", and "ShiftML2.0". + frames : List[ase.Atoms] + A list of candidate structures. + list_atom : np.array(List[str]) + An array of atom symbols included in the structure. + list_cs: np.array(List[float]) + An array of chemical shielding values corresponding to the atom symbols. + """ + calculator = ShiftML(model_version) + number_list = list_atom[-1].split(",") + num_atoms = float(number_list[-1]) + rmse_rec1 = np.array([]) + rmse_rec2 = np.array([]) + + for frame in frames: + Y = list_cs + atom_label = frame.get_atomic_numbers() == 1 + array = calculator.get_cs_iso(frame).ravel()[atom_label] + X = extract_from_array(array, num_atoms, list_atom) + slope = -1 + intercept = np.mean(Y) - slope * np.mean(X) + rmse = root_mean_squared_error(slope * X + intercept, Y) + rmse_rec1 = np.append(rmse_rec1, rmse) + if GIPAW_avail: + array = frame[atom_label].arrays[cs_sym].ravel() + X = extract_from_array(array, num_atoms, list_atom) + slope = -1 + intercept = np.mean(Y) - slope * np.mean(X) + rmse = root_mean_squared_error(slope * X + intercept, Y) + rmse_rec2 = np.append(rmse_rec2, rmse) + rmse_rec = (rmse_rec1, rmse_rec2) + else: + rmse_rec = rmse_rec1 + + return rmse_rec From 9f2f7a99f617e9897a18e7274fac332b0d0bf2e4 Mon Sep 17 00:00:00 2001 From: sovietdevil <2558390548@qq.com> Date: Tue, 30 Jul 2024 13:10:41 +0200 Subject: [PATCH 6/7] changes debugging level of ase calculator --- src/shiftml/ase/calculator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shiftml/ase/calculator.py b/src/shiftml/ase/calculator.py index 49eee95..99b1f2d 100644 --- a/src/shiftml/ase/calculator.py +++ b/src/shiftml/ase/calculator.py @@ -8,7 +8,7 @@ # For now we set the logging level to DEBUG logformat = "%(asctime)s - %(levelname)s - %(message)s" -logging.basicConfig(level=logging.DEBUG, format=logformat) +logging.basicConfig(level=logging.INFO, format=logformat) url_resolve = { From 52936976f4e08ed2304724b001bd14759dbe66ef Mon Sep 17 00:00:00 2001 From: sovietdevil <2558390548@qq.com> Date: Tue, 30 Jul 2024 14:40:29 +0200 Subject: [PATCH 7/7] A notebook example for the chemical structure identification based on the predicted chemical shielding. --- examples/test_pred.ipynb | 244 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 examples/test_pred.ipynb diff --git a/examples/test_pred.ipynb b/examples/test_pred.ipynb new file mode 100644 index 0000000..1469e05 --- /dev/null +++ b/examples/test_pred.ipynb @@ -0,0 +1,244 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "from sklearn.metrics import root_mean_squared_error\n", + "from ase.io import read\n", + "from matplotlib import pyplot as plt\n", + "\n", + "from shiftml.ase import ShiftML\n", + "import urllib.request\n", + "\n", + "\n", + "def extract_from_array(array, num_atoms, list_atom):\n", + " \"\"\"\n", + " Function to extract regression array from a given\n", + " symbol list of atoms and data exported from the model\n", + " \"\"\"\n", + " num_molecules = int(len(array) / num_atoms)\n", + " data_fit = array.reshape((-1, num_molecules))[:, 0]\n", + " X = []\n", + " for atom_string in list_atom:\n", + " label_list = atom_string.split(\",\")\n", + " if len(label_list) == 1:\n", + " label = int(label_list[0])\n", + " X.append(data_fit[label - 1])\n", + " else:\n", + " X.append(\n", + " sum([data_fit[int(label_str) - 1] for label_str in label_list])\n", + " / len(label_list)\n", + " )\n", + " return np.array(X)\n", + "\n", + "\n", + "def structure_prediction(\n", + " model_version, frames, list_atom, list_cs, GIPAW_avail=True, cs_sym=\"CS\"\n", + "):\n", + " \"\"\"\n", + " Function to select the suitable structures\n", + " based on a set of candidate structures,\n", + " given rmse of the linear regression results\n", + "\n", + " Parameters\n", + " ----------\n", + " model_version : str\n", + " The version of the ShiftML model to use. Supported versions are\n", + " \"ShiftML1.0\", \"ShiftML1.1\", and \"ShiftML2.0\".\n", + " frames : List[ase.Atoms]\n", + " A list of candidate structures.\n", + " list_atom : np.array(List[str])\n", + " An array of atom symbols included in the structure.\n", + " list_cs: np.array(List[float])\n", + " An array of chemical shielding values corresponding to the atom symbols.\n", + " \"\"\"\n", + " calculator = ShiftML(model_version)\n", + " number_list = list_atom[-1].split(\",\")\n", + " num_atoms = float(number_list[-1])\n", + " rmse_rec1 = np.array([])\n", + " rmse_rec2 = np.array([])\n", + "\n", + " for frame in frames:\n", + " Y = list_cs\n", + " atom_label = frame.get_atomic_numbers() == 1\n", + " array = calculator.get_cs_iso(frame).ravel()[atom_label]\n", + " X = extract_from_array(array, num_atoms, list_atom)\n", + " slope = -1\n", + " intercept = np.mean(Y) - slope * np.mean(X)\n", + " rmse = root_mean_squared_error(slope * X + intercept, Y)\n", + " rmse_rec1 = np.append(rmse_rec1, rmse)\n", + " if GIPAW_avail:\n", + " array = frame[atom_label].arrays[cs_sym].ravel()\n", + " X = extract_from_array(array, num_atoms, list_atom)\n", + " slope = -1\n", + " intercept = np.mean(Y) - slope * np.mean(X)\n", + " rmse = root_mean_squared_error(slope * X + intercept, Y)\n", + " rmse_rec2 = np.append(rmse_rec2, rmse)\n", + " rmse_rec = (rmse_rec1, rmse_rec2)\n", + " else:\n", + " rmse_rec = rmse_rec1\n", + "\n", + " return rmse_rec\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "# atoms list defines groups of atoms to average over\n", + "\n", + "list_atom = ['1', '2', '3', '4', '5', '6', '7', '8', '9', '10',\n", + " '11,12,13', '14', '15', '16', '17', '18', '19,20,21']\n", + "\n", + "list_cs = np.array([3.76, 3.78, 5.63, 3.32, 3.49, 3.06, 2.91, 3.38, 2.56, 2.12, 1.04, 8.01, 8.01,\n", + " 8.01, 8.01, 8.01, 3.78])" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "urllib.request.urlretrieve(\"https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-06972-x/MediaObjects/41467_2018_6972_MOESM8_ESM.zip\", \"cocaine.zip\")\n", + "import zipfile\n", + "with zipfile.ZipFile(\"cocaine.zip\",\"r\") as zip_ref:\n", + " zip_ref.extractall(\".\")\n", + "\n", + "frames = read(\"./Supplementary_Dataset_6/cocaine_QuantumEspresso.xyz\", \":\")" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2024-07-30 13:26:15,484 - INFO - rascaline version: 0.1.0.dev558\n", + "2024-07-30 13:26:15,485 - INFO - rascaline-torch is installed, importing rascaline-torch\n", + "2024-07-30 13:26:15,485 - INFO - Found model version in url_resolve\n", + "2024-07-30 13:26:15,486 - INFO - Resolving model version to model files at url: https://tinyurl.com/f237evr3\n", + "2024-07-30 13:26:15,486 - INFO - Found ShiftML1.1 in cache, and importing it from here: /Users/zhangyuxuan/Library/Caches/shiftml/ShiftML1.1\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+qElEQVR4nO3deXiU5aHH/d8kZIMsyJaFBAJioYigQolA2TSyaEGMFNwKoqWnFs4B87ZFrBKXHqKe1mItLW+tFL1akcWAsXJQS4nGglK2o75UWpESliSA1oSEJTjzvH8MGTIkgUxmMvc8M9/Pdc0V5sn9zNyZDJPfc68Oy7IsAQAAGBJlugIAACCyEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGNXOdAVawuVy6ciRI0pKSpLD4TBdHQAA0AKWZenEiRPKyMhQVFTz7R+2CCNHjhxRVlaW6WoAAIBWOHjwoDIzM5v9vi3CSFJSkiT3D5OcnGy4NgAAoCWqq6uVlZXl+TveHFuEkfqumeTkZMIIAAA2c6khFgxgBQAARhFGAACAUYQRAABglC3GjLSE0+nU2bNnTVcj4kVHR6tdu3ZMwQYAtFhYhJGamhodOnRIlmWZrgoktW/fXunp6YqNjTVdFQCADdg+jDidTh06dEjt27dX165duSI3yLIs1dXV6dixY9q/f7+uuOKKiy5yAwCAFAZh5OzZs7IsS127dlVCQoLp6kS8hIQExcTE6MCBA6qrq1N8fLzpKgEAQlzYXLbSIhI6aA0BAPjC9i0jgeJ0SqWlUnm5lJ4ujRwpRUebrhWAcMbnDuBGGJFUVCTNmycdOnT+WGam9OyzUl6euXoBCF987gDnRXx7elGRNHWq9weCJB0+7D5eVGSmXgDCF587gDefw8i7776rSZMmKSMjQw6HQ+vXr79o+aKiIt14443q2rWrkpOTNWzYML355putrW9AOZ3uK5OmZgTXH5s/312urVRUVGjevHnq06eP4uPjlZqaqhEjRug3v/mNTp48KUnKzs7WkiVLPOdkZ2fL4XDI4XCoQ4cOuvbaa7VmzRqvxz116pQ6deqkLl266MyZM57jy5YtU1JSkr766ivPsZqaGsXExGjMmDFej1FSUiKHw6F9+/YF/gcHIlQofO4AocbnMFJbW6tBgwZp6dKlLSr/7rvv6sYbb9SGDRu0Y8cOjR07VpMmTdKuXbt8rmyglZY2vjJpyLKkgwfd5drCZ599pmuuuUZvvfWWFi9erF27dmnr1q368Y9/rD/96U/685//3Oy5jz/+uMrLy7Vr1y594xvf0PTp07VlyxbP91999VVdeeWV6tevn1dgHDt2rGpqarR9+3bPsdLSUqWlpemDDz7Q6dOnPcc3b96sHj166PLLLw/sDw5EMNOfO0Ao8nnMyMSJEzVx4sQWl294RS9Jixcv1muvvabXX39d11xzja9PH1Dl5YEt56sf/OAHateunbZv364OHTp4jvfu3Vu33HLLRRdxS0pKUlpamtLS0rR06VL94Q9/0Ouvv67hw4dLkl544QXdfffdsixLL7zwgqZPny5J6tu3r9LT01VSUqLrrrtOkrsF5JZbbtFf/vIXvf/++54WkpKSEo0dO7ZtfvgwwOBDtIbpzx0gFAV9zIjL5dKJEyfUqVOnZsucOXNG1dXVXre2kJ4e2HK++Pzzz/XWW29pzpw5XkGkoZZOV27Xrp1iYmJUV1cnSdq3b5+2bt2qadOmadq0aSotLdWBAwc85ceOHavNmzd77m/evFljxozR6NGjPcdPnTqlDz74gDDSjKIiKTtbGjtWuvNO99fsbPv19TudUkmJtHKl+ytdA23P5OcOEKqCHkZ+9rOfqaamRtOmTWu2TGFhoVJSUjy3rKysNqnLyJHu0evN/c13OKSsLHe5QPv0009lWZb69u3rdbxLly5KTExUYmKiFixYcMnHqaurU2FhoaqqqnT99ddLkpYvX66JEyfqsssuU6dOnTR+/Hj9/ve/95wzduxY/fWvf9VXX32lEydOaNeuXRo9erRGjRqlkpISSdLWrVt15swZwkgTwmXwYbgEKrsx+bkDhKqghpGXX35Zjz32mFavXq1u3bo1W27hwoWqqqry3A4ePNgm9YmOdk+jkxp/MNTfX7IkuE3v27Zt0+7du3XllVd6DTy90IIFC5SYmKj27dvrqaee0pNPPqmbb75ZTqdTL774ou6++25P2bvvvlsrVqyQy+WSJI0ZM0a1tbX629/+ptLSUn3ta19T165dNXr0aM+4kZKSEvXu3Vs9evRo85/ZTsJl8GG4BCo7CsXPHcC0oIWRV155Rd/97ne1evVq5ebmXrRsXFyckpOTvW5tJS9PWrtWysjwPp6Z6T7eVvP9+/TpI4fDob1793od7927t/r06XPJpe1/9KMfaffu3Tp06JD+/e9/e1pR3nzzTR0+fFjTp09Xu3bt1K5dO91+++06cOCANm3a5HnuzMxMbd68WZs3b9bo0aMlSRkZGcrKytKWLVu0efNmT0sLzguHwYfhEqjszNTnDhCqghJGVq5cqVmzZmnlypW6+eab2/S5LMuS6+RJn25TJpzUxztOKlanFaMzemP9ae37/9zHfXkcX3YN7ty5s2688Ub96le/Um1trc8/Z5cuXdSnTx+lpaV5jS154YUXdPvtt2v37t1et9tvv10vvPCCp9zYsWNVUlKikpISrym9o0aN0v/+7/9q27ZtdNE0IRwGH4ZDoAoHeXnSnj3n72/YIO3fTxBBZPJ5Nk1NTY0+/fRTz/39+/dr9+7d6tSpk3r06KGFCxfq8OHDeumllyS5u2ZmzpypZ599Vjk5OaqoqJDk3lAtJSUlQD/GedapU9p77eBWnbu7fvjGAunTi5ZsWt+dO+Ro377F5X/9619rxIgRGjJkiB599FENHDhQUVFR+tvf/qZPPvlEgwf79nMcO3ZMr7/+uoqLizVgwACv782YMUO33nqrvvjiC3Xq1Eljx47VnDlzdPbsWU/LiCSNHj1ac+fOVV1dHWGkCeEw+DAcAlW4aNgVM2oUXTOIXD63jGzfvl3XXHONZ1pufn6+rrnmGi1atEiSVF5errKyMk/53/72t/rqq680Z84cpaene27z5s0L0I9gX5dffrl27dql3NxcLVy4UIMGDdKQIUP03HPP6Yc//KGeeOIJnx7vpZdeUocOHXTDDTc0+t4NN9yghIQE/eEPf5Dkbhk5deqU+vTpo9TUVE+50aNH68SJE54pwPAWDoMPwyFQAQgvDsuXvgVDqqurlZKSoqqqqkbjR06fPq39+/erV69eio+Pl2VZsk6dMlJPR0ICuwer8e8k3NQP/pS8x13U/+pDvc/f6XTPmjl8uOlxIw6HO3Dt38+VelurrZUSE93/rqmRmpnlD9jWxf5+NxR2G+U5HA6fukoAX9UPPvyv/3L/Qa+XmemeBRHKQUQ6P5tj6lR38GgqUDGbA0AwRfxGeUBr2H3wIbM5AISSsGsZAYLF7oMP8/Kk3Fypfhz5hg3SuHH2+zkA2B8tI0AEs3ugAhAeCCMAAMAowggAADCKMAIAAIxiAGs9l1M6sEWqqZQSU6Wew6UoOtABAGhrtIxI0p5iackA6cVvSa/e5/66ZID7uEEOh0Pr169v9vslJSVyOBz68ssvPcfWr1+vPn36KDo6WvPnz2/zOgIA4C/CyJ5iafUMqfqI9/HqcvfxNgwkx44d0/33368ePXooLi5OaWlpGj9+vP7617+26Pzhw4ervLzca4+f//iP/9DUqVN18OBBPfHEE7rnnns0ZcqURuc6HA45HA69//77XsfPnDmjzp07y+FwqKSkxKv8xYIRAACtFdndNC6ntHGBpKZWxLckOaSND0r9bm6TLpvbbrtNdXV1evHFF9W7d29VVlZq06ZN+vzzz1t0fmxsrNLS0jz3a2pqdPToUY0fP14ZF65m1YSsrCz9/ve/13XXXec5tm7dOiUmJuqLL77w/QcCAKAVIrtl5MCWxi0iXiyp+rC7XIB9+eWXKi0t1VNPPaWxY8eqZ8+eGjp0qBYuXKjJkyd7yh0/fly33nqr2rdvryuuuELFxedbahp205SUlCgpKUmSdP3118vhcGjMmDF68cUX9dprr3laQhq2dsycOVOvvPKKTjXYy2f58uWaOXNmwH9eAACaE9lhpKYysOV8kJiYqMTERK1fv15nzpxpttxjjz2madOm6cMPP9RNN92ku+66q8lWi+HDh2vv3r2SpFdffVXl5eUqLi7WtGnTNGHCBJWXl6u8vFzDhw/3nDN48GBlZ2fr1VdflSSVlZXp3Xff1Xe+850A/7QAADQvssNIYmpgy/mgXbt2WrFihV588UV17NhRI0aM0EMPPaQPP/zQq9w999yjO+64Q3369NHixYtVU1Ojbdu2NXq82NhYdevWTZLUqVMnpaWlKTk5WQkJCZ7xKGlpaYqNjfU6795779Xy5cslSStWrNBNN92krl27BvznBQCgOZEdRnoOl5IzJDmaKeCQkru7y7WB2267TUeOHFFxcbEmTJigkpISXXvttVqxYoWnzMCBAz3/7tChg5KTk3X06NGA1eHuu+/W1q1b9dlnn2nFihW69957A/bYAAC0RGSHkahoacJT5+5cGEjO3Z/wZJuuNxIfH68bb7xRjzzyiLZs2aJ77rlHBQUFnu/HxMR418rhkMvlCtjzd+7cWd/61rd033336fTp05o4cWLAHhsAgJaI7DAiSf0nS9NekpLSvI8nZ7iP95/c9HltVZ3+/VVbWxuwx4uNjZXT6bxomXvvvVclJSWaMWOGotkpDQAQZGE3tdeyLFkNZoe0SHauNOsd6Rdfc9+f9kep12h3i8jJky1+GEdCghyO5rp8vH3++ef69re/rXvvvVcDBw5UUlKStm/frqefflq33HKLb/W/iOzsbL355pvau3evOnfurJSUlEatLRMmTNCxY8eUnJx80cfav3+/du/e7XXsiiuuUIcOHQJWXwBA5Am/MHLqlPZeO7iVZ59bm2PNj1p1dt+dO+Ro375FZRMTE5WTk6Nf/OIX2rdvn86ePausrCzNnj1bDz30UKuevymzZ89WSUmJhgwZopqaGm3evFljxozxKuNwONSlS5dLPlZ+fn6jY6WlpfrmN78ZqOoCACJQ2IURu4iLi1NhYaEKCwubLWNZjRdja7j0+5gxY7zKdOzYsdE5Xbt21VtvvdWix77Y41ysPAAA/gi7MOJISFDfnTuMPTcAAMHidEqlpVJ5uZSeLo0cKdlx6F/4hRGHo8VdJQAA2FVRkTRvnnTo0PljmZnSs89KeXnm6tUazKYBAMBmioqkqVO9g4gkHT7sPl5UZKZerUUYAQDARpxOd4tIU0P56o/Nn+8uZxeEEQAAbKS0tHGLSEOWJR086C5nF2ETRpjtETr4XQBA2ykvD2y5UGD7MFK/YmhdXZ3hmqDeyXMLxV24uBoAwH/p6YEtFwpsP5umXbt2at++vY4dO6aYmBhFRdk+X9mWZVk6efKkjh49qo4dO7K0PAC0gZEj3bNmDh9uetyIw+H+/siRwa9ba9k+jDgcDqWnp2v//v06cOBA0J7XsqQzZ9wDhKKjpbg49xsA7kXT0tLSLl0QAOCz6Gj39N2pU91/dxoGkvq/Q0uW2Gu9EduHEcm9GdwVV1wRtK6at96SFi+WKirOH0tLkx56SBo3LihVCFkxMTG0iABAG8vLk9aulf7rv9wtJPUyM91BxG7rjIRFGJGkqKgoxcfHt/nz1M/tvrBprKxMmjLF/eaw25sAAGA/eXlSbq6UkuK+v2GD+4LYjteDDLDwQTjO7QYA2FfD4DFqlD2DiEQY8Uk4zu0GAMA0wogPwnFuNwAApoXNmJFgCMe53UAkC5cdT+2M3wEkWkZ8Uj+3u7kpvA6HlJVlr7ndQKQqKpKys6WxY6U773R/zc623wZjdsbvAPUIIz6on9stNQ4kdp3bDUSicNvx1I74HaAhwoiP6ud2Z2R4H8/MZFovYAfMijOP3wEuRBhphbw8ac+e8/c3bJD27yeIAHbArDjz+B3gQoSRVgqXud1ApGFWnHn8DnAhwgiAiMKsOPP4HeBChBEAEYVZceYF8nfgdEolJdLKle6vjDOxJ8IIgIjCrDjzAvU7YGpw+CCMGEKaB8xhVpx5/v4OmBocXggjBpDmAfOYFWdea38HTA0OPz6HkXfffVeTJk1SRkaGHA6H1q9ff8lzSkpKdO211youLk59+vTRihUrWlHV8ECaB0IHs+LMa83vgKnB4cfnMFJbW6tBgwZp6dKlLSq/f/9+3XzzzRo7dqx2796t+fPn67vf/a7efPNNnytrd6R5APAfU4PDj88b5U2cOFETJ05scflly5apV69e+vnPfy5J+vrXv6733ntPv/jFLzR+/Hhfn97WfEnzY8YErVoAYCtMDQ4/bT5mZOvWrcrNzfU6Nn78eG3durXZc86cOaPq6mqvWzggzQOA/5ieHX7aPIxUVFQoNTXV61hqaqqqq6t16tSpJs8pLCxUSkqK55aVldXW1QwK0jwA+I/p2eEnJGfTLFy4UFVVVZ7bwYMHTVcpIEjzABAYTM8OLz6PGfFVWlqaKisrvY5VVlYqOTlZCQkJTZ4TFxenuLi4tq5a0NWn+alT3cGj4UBW0jxgP06ne4xXebm7RXPkSP7/BlNenpSbK6WkuO9v2CCNG8fvwI7avGVk2LBh2rRpk9ext99+W8OGDWvrpw5JpHkgPLBeUGhgenZ48DmM1NTUaPfu3dq9e7ck99Td3bt3q6ysTJK7i2XGjBme8t///vf12Wef6cc//rE++eQT/frXv9bq1av1wAMPBOYnsCEWWwoTLqdG9yzV7QPWKqqsVHIxJztSsF4QEFg+d9Ns375dY8eO9dzPz8+XJM2cOVMrVqxQeXm5J5hIUq9evfTGG2/ogQce0LPPPqvMzEz97ne/i7hpvRcizdvcnmIlbFigknuOuO+vkpScIU14Suo/2WjV0LYutV6Qw+FeL+iWW/h/DbSUz2FkzJgxspr6X3hOU6urjhkzRrt27fL1qcKby6nRPbcoPalSUWWpUt/hUhSfXLawp1haPUMOXfD/oLpcWj1DmvYSgSSMsV4QEHhtPoAVTeCq2r5cTmnjAkmWGk+KsiQ5pI0PSv1uJlyGKdYLAgIvJKf2hrX6q+qaI97H66+q9xSbqRda5sAWqfrIRQpYUvVhdzmEJdYLAgKPMBJMl7yqlvuqmoGQoaum8tJlfCkH22G9ICDwCCPBxFW1/SWmXrqML+VgOwFd/ZMZWYAkwkhwcVVtfz2Hu8f3NNG25eaQkru7yyFsBWS9oD3FSvh/B6jknm9p5W33KWHVt6QlA+iqRUQijAQTV9X2FxXtHmispjrbzt2f8CSDVyOAX+sFMXYM8EIYCSauqsND/8nStJdkJaZ5H0/OYFpvhGnVekGMHQMaIYwEE1fV4aP/ZJ2atc1z99Rta6X5HxFEcGmMHQsrTqdUUiKtXOn+6iRDtgphJNi4qg4fDUKjK5NF69BCjB0LG+xPFDiEERO4qgYiF2PHwgL7EwUWYcQUrqqByMTYMdu71P5Eknt/IrpsWo4wAgDBxNgx2/NlfyK0DGEEAIKNsWO2xv5EgcdGebAtp9N95VFe7t4HZORItmyHjfSfrFOpY9ThuSxJ7rFjCVdeT4uIDbA/UeDRMgJbYhQ7wgJjx2yJ/YkCjzAC22EUOwAPA/v7BHR/IkgijMBmGMUOwMPg/j4B2Z8IHoQR2Aqj2AFICon9ffzanwheCCOwFUaxAwil/X1atT8RGiGMwFYYxQ4goPv7uJzS/lLpo7Xur2xQaARTe2Er9aPYDx9uetyIw+H+PqPYW8jl1OieW5SeVKmoslSpLzM6YAOB2t9nT7G7haVhsEnOcC9Kx1ovQUUYga3Uj2KfOtUdPBoGEkax+2hPsRI2LFDJPec+iFeJD2LYQyD29zk35sTTrVOvfsyJXRafC5MLCrppYDuMYg+AEBj8B7Sav/v7NBhz0lhwx5w4nVJJibRypfurTzMBDc4mCjTCCGwpEKPY/foQsLMQGvwHtIq/+/sEcsyJH/xavDHMLigIIzYWsX9Mz/FnFHtRkdS7l1OP3lOq4sVr9eg9perdyxkZC6aFyAcx4Bd/9vcJ1JgTP/i1eGMYXlAwZsSmiorci381fCNnZrrHU9BNcXFFRdIfHy7We7ctUFbK+T/KB6syNP/hpyRNDu/XMAQ+iIGAaO3+PoEYc+KHSy3e6HC4F2+85ZZmLrJ8uaDoZY/R/LSM2BDLobee0yn97zPFWvPtGeqe7P2fuXtyudZ8e4Y2/qI4vFuZDH8QAwHVmv19/B1z4ie/F28MwwsKwojNsBx6A63Yk6L0XacWfcPdvBl1wedQlMP9Aj485EGVvhvGL6DhD2KEjojt6m0w5qTx/4MWjDnxk9+LN4bhBQVhxGZYDv2cVo4id362RVkpRxoFkXpRDks9Ug7L+VkYj5fwd/AfwkLE73xdP+YkyXvMidWSMSd+8nvxxjC8oCCM2AzLocuvUeTpiS1rtmxpOdvyZ/BfAEXslblhdPW6FX0yWV9fus1zf8If1ip7yUcq+qRt3//1izdeuONvPYdDysq6yOKNYXhBQRixmYhfDt3PUeR9h7Ss2bKl5Wyt/2SdmnX+g/jUbWul+R8FLYhE/JW5IXT1utUHsoOHzv/BLi0broOHon0LZK3oLq5fvFFqHEhavHhjiFxQBAphpLVa8QYMBL8Ttd35OS01utdwnYzJkMtq+gV0WQ6djOmu6F72ad70S2sG/wUAV+bm0NUbwEDmx6JjAVm80fAFRSARRlrD4Kp3AUnUdubvKPKoaLW/9Sk5HGoUSFyWQw6H1P5WezVv2g1X5mbR1RugQBaARccCsXijqQuKQCOM+CoEVr2L6OXQAzGKvP9kOaa9JF0wcM2RkuE+bsOrCjvhytysiO/qVQACWQAXHfNn8cZwQhjxRQiteheQRG1HgRpF3kTzpsOmzZt2w5W5WRHf1asABDJWMQ44wogvQuwNGJGJOpCjyMOkedNuuDI3K+K7ehWAQBaGi46ZRhjxBW/A0BBmo8gjDVfm5oVUV6+ByQB+B7IwXHTMNMKIL3gDho4wGkUeabgyDw15edKej88HgfdfKdX+fc7gBhGDkwHqA9mFLXAtCmRhuOiYaYQRX/AGDC10s9hWoK7MWTTND3uKlfi780Eg5+/fUvRzwQkC9c8fCpMBduw4f39dUQvH3oXhomOmEUZ8wRsQCBh/B2GzaJofTAeBEJoMEN3gr+CIET60yNFdHFCEEV/xBgQCprWDsFk0zQ+hEARCbDJAq9FdHDCEkdbgDQgERisGL7Jomp9CIQiE02QAuosDop3pCtgWb0DAP3uKlbBhgUruOfeHcZXcLYwTnrposPdl0bQxYwJa4/AQCkGAyQC4AC0jAILPjzELIbVomqE9qvwSCkGAyQC4QKvCyNKlS5Wdna34+Hjl5ORo27ZtFy2/ZMkS9e3bVwkJCcrKytIDDzyg06dPt6rCAGzOzzELIbNomsFpqX4JhSDAZABcwOcwsmrVKuXn56ugoEA7d+7UoEGDNH78eB09erTJ8i+//LIefPBBFRQU6O9//7teeOEFrVq1Sg899JDflYd/mBYJI/wcsxASi6aZno3ij1AJAkwGcLNj61ob8DmMPPPMM5o9e7ZmzZql/v37a9myZWrfvr2WL1/eZPktW7ZoxIgRuvPOO5Wdna1x48bpjjvuuGRrCtoW0yJhjJ9jFowvmhYKs1H8FSpBINInA9i1da0N+BRG6urqtGPHDuXm5p5/gKgo5ebmauvWrU2eM3z4cO3YscMTPj777DNt2LBBN910U7PPc+bMGVVXV3vdEDhMi4RRARizYHQ581CYjRIIoRIEInUyQIi0roVKC7lPYeT48eNyOp1KTfX+kEhNTVVFRUWT59x55516/PHH9c1vflMxMTG6/PLLNWbMmIt20xQWFiolJcVzy8rK8qWauAimRcK4AI1ZMLZzdSjMRgmUSA0CpoVI61ootZC3+WyakpISLV68WL/+9a+1c+dOFRUV6Y033tATTzzR7DkLFy5UVVWV53bw4MG2rmbE8GVaJNAmAjhmIdpxvr99THapoh1BSNGhMBsF9hYCrWuh1kLu0zojXbp0UXR0tCorvRN/ZWWl0tLSmjznkUce0Xe+8x1997vflSRdddVVqq2t1fe+9z395Cc/UVRU4zwUFxenuLg4X6qGFgqpaZGIXPVjFjb8WI6aBm+25Ax3EGlJV0Er1ynxW33LTnW5PFexXhzu7zMtFc0x3Lp2qRZyh8PdQn7LLcHbsNKnlpHY2FgNHjxYmzZt8hxzuVzatGmThg0b1uQ5J0+ebBQ4os/9dFZTrwTaVMPpjlENripH9yxVVIOryjafFgn4M2bBZH97qMxGgX0Zbl0LxRZyn1dgzc/P18yZMzVkyBANHTpUS5YsUW1trWbNmiVJmjFjhrp3767CwkJJ0qRJk/TMM8/ommuuUU5Ojj799FM98sgjmjRpkieUIHjqp0UOTSrWkvELlJVy/sP8YFWG5r/5lP5WM7ltp0UC9VozZuGS/e0Od397v5vbLhAEomUHkctw61ootpD7HEamT5+uY8eOadGiRaqoqNDVV1+tjRs3ega1lpWVebWEPPzww3I4HHr44Yd1+PBhde3aVZMmTdJ///d/B+6nQItFR0trHi/W0H/N0IX/Cbonl2vNt2doW/ZLio7mwxQhypf+9l5tmKr7T9ap1DHq8Jx7gP2p29Yq4crraRHBpdW3rq2eIUsOObw+i9u+dS1kFg5soFV708ydO1dz585t8nslJSXeT9CunQoKClRQUNCap0KguZy67t8LZDkaX1VGOSxZcui6fz8oudrwqvIcp9PdDFhe7n7TjxwZvP5J2Fgg+9tdTo3uuUXpSZWKKkuV+vo4o4TZKGgtg61r9S3khw83PW7E4XB/P5gt5GyUF2nOXVU2P6kyOFeVRUXuAVQN+y0zM92LWbX51EzYW6D6200NgAXqGWpdq184cOpUd/BoGEiCsnBgE9goL9KEwBoJoTalzBiXU9pfKn201v01lFfsbIaRBZMCsU5JiCw4BZhqXTO6cGATCCORxvAobhZdO2dPsXvZ5xe/Jb16n/urzZaBNrZgkr+zWUJkwSnANGMLBzaBMBJpDO/YGYpTyoLu3FV5o0GYNroqN9665c/eKiGw4BQQKhp2xYwaZW7cHmEk0hheIyEUp5QFVYOr8sbscVUeMq1brV2nJAS6KgF4I4xEIoM7dobilLKgCoOr8pBq3WpNfzvLuQMhhzASqQzt2Fk/pezCrd/rORxSVlZwp5QFVRhcldu+dctwVyWAxggjkczAKO76KWVS40BiakpZUIXBVbntW7dYzh0IOYQRBF2oTSkLqjC4Kg+L1i2DXZUAGiOMwIhQmlIWVA2uyhsHEntclYdN65ahrsqwEgZr5SA0EEZgTKhMKQu6c1flSrLvVXnYtG6xnHvrhcFaOQgdhBHAhP6TpTnnr8p1l/2uyiO2dQuetXKsC2aGWTZaKwehhTBiZy6nRvcs1e0D1iqqjCZS22l4Fd7TnlflEdu6FcnOrZVjNTn813KvlhPia+Ug9LBRnl2xyRcAE0Jks02EF1pG7IhNvgAY4qpu2Ro4LS0HSIQR+2GTLwAGffhZy9bAaWk5QCKM2E8ILSduZPt4AEZ9cmq4DlZlyGU13VHjshwqq+quT06F7lo5CD2EEbsJkeXEjW0fD8CotIxozdvoXivnwkBSf3/+xieVlsFoZrQcYcRuQmA5cePbxwMwZuRI6W81k/XtNS/pyAnvtXIOVWfo22te0vbayaG9Ai9CDmHEbgwvJx4y28cDMKJ+Bd51n0xW/1+fXytnwh/WqvcvP9K6TybbYwVehBTCiN0Y3uQrpLaPB2BE/Qq8qWnnP2dKy4Yro3u0vVbgRcggjNiRwU2+bL99PICAyMuTduw4f39dESvwovVY9Myu+k/WqdQx6vBcliT3Jl8JV17f5qt42n77eAABE93gcnbECLpm0Hq0jNiZgU2+wmL7eABASCGMwCdhs308ACBkEEbgs7DZPh4AEBIYM4JWycuTcnOllBT3/Q0bpHHjaBEBAPiOMIJWY/t4s5xO9xTq8nL3gOGRI/kdALAnwghgQ0VF7sXnGq75kpnpHs9DNxkAu2HMCGAzLMcPINwQRgAbYTl+AOGIMALYCMvxw4vLKe0vlT5a6/7qIoXCnhgzAthIw2X2oxxOjeyxRelJlSo/karSsuFyWdGNyiFM7SmW9b8L5DhxxHPISsqQY+JTbbolBMKMy6nRPd2fI1FlqVLf4CygeSHCCGAj9cvs39qvWM9OWKCslPN/iA5WZWjexqe07pPJLMcf7vYUy1o9Q5ZleS0+aFWXS6tnyNHGe1QhTOwpVsKGBSq559znyCq59zibEPxASzcNYCMjR0rfHVGstdNmqHvyEa/vdU8u19ppMzT7m8Usxx/OXE6dXLdAlmUp6oJVkKMclixLOrnuQbpscHF7it3Btcb7c0TnAq32FAe1OoQRwEaiHU49O3GBpKb/EEnSkgkPKtrBH6Jw5dy/Re3PHmn0+68X5bDU/uxhOfdvCW7FYB8up7TR/TnS+G10biT8xuAGWsIIYCcHWvaHSAeC9IfI5dTonqW6fcBaRZUxgDIY9m6vDGg5RKADW6TqIxcpYEnVQfwcEWNGAHupaeEfmJaW80cI9TdHkvKaVPUPYDlEoFD6HDmHlhHAThJTA1uutUKsvzmSRPceroNVGXJZTTePuSyHyqq6K7r38CDXDLYRKp8jDRBGADvpOdzd+tBET6+bQ0ru7i7XVkKwvzmSjBwVrcf/9pQkNQok9fd/uv1JjRzFRkVoRih8jlyAMAIY0nCV1HdLW7hqalS0uxtETUWBc/cnPNm26wSEYH9zJImOlibmT9a317ykIyfSvL53qDpD317zkiY8MJlNE9G8UPgcubBKQXsmAB5FRdLXG3ToT5woZWe3cF+Z/pOlaS/JSvT+Q6TkDCkY60uEYH9zpMnLk+766WTdsG6b59iEP6zVqKKPdNdPJ7NZIi7N9OfIBQgjQJDVb3R3+LD3cZ82uus/Wadmnf9DdOq2tdL8j4LzARKC/c2RKC9P+tv281eu+UuGa99n0QQRtJzJz5ELEEaAIAroRncNmlBdmUFcwjkE+5sjVXSDT/ARI0TXDHxn6nPkwmoYeVYgQoXFRnch2N8MwN5aFUaWLl2q7OxsxcfHKycnR9u2bbto+S+//FJz5sxRenq64uLi9LWvfU0bNmxoVYUBO2vpBnYhv9FdiPU3A7A3nxc9W7VqlfLz87Vs2TLl5ORoyZIlGj9+vPbu3atu3bo1Kl9XV6cbb7xR3bp109q1a9W9e3cdOHBAHTt2DET9AVtp6QZ2ttjorv9knUodow7PZUly9zcnXHk9LSIAfOZzGHnmmWc0e/ZszZo1S5K0bNkyvfHGG1q+fLkefPDBRuWXL1+uL774Qlu2bFFMTIwkKTs7279aAzY1cqSUmdl48Go9h8P9fdtsdBci/c0A7M2nbpq6ujrt2LFDubm55x8gKkq5ubnaunVrk+cUFxdr2LBhmjNnjlJTUzVgwAAtXrxYzouM0Dtz5oyqq6u9bkA4iI6Wnn3W/W/HBcMt6u8vWcJARACRxacwcvz4cTmdTqWmek/ZS01NVUVFRZPnfPbZZ1q7dq2cTqc2bNigRx55RD//+c/105/+tNnnKSwsVEpKiueWlZXlSzWBkJaXJ61d27grJjPTfZypmQAiTZvPpnG5XOrWrZt++9vfavDgwZo+fbp+8pOfaNmyZc2es3DhQlVVVXluBw8ebOtqAkGVlyft2HH+/roiaf9+ggiAyOTTmJEuXbooOjpalZXeKytWVlYqLS2tyXPS09MVExOj6Abtzl//+tdVUVGhuro6xcbGNjonLi5OcXFxvlQNsB3WiAAAN59aRmJjYzV48GBt2rTJc8zlcmnTpk0aNmxYk+eMGDFCn376qVwul+fYP/7xD6WnpzcZRAAAQGTxuZsmPz9fzz//vF588UX9/e9/1/3336/a2lrP7JoZM2Zo4cKFnvL333+/vvjiC82bN0//+Mc/9MYbb2jx4sWaM2dO4H4KAABgWz5P7Z0+fbqOHTumRYsWqaKiQldffbU2btzoGdRaVlamqKjzGScrK0tvvvmmHnjgAQ0cOFDdu3fXvHnztGDBgsD9FAAAwLZ8DiOSNHfuXM2dO7fJ75WUlDQ6NmzYML3//vuteSoAABDm2JsGAAAY1aqWEUCS5HJqdM8tSk+qVFRZqtSXFTgBAL4jjKB19hQrYcMCldxzxH1/ldybpE14ik3SAAA+oZsGvttTLK2eIUfNEe/j1eXS6hnu77eEy6nRPUt1+4C1iiorlVzNbxEAAAhftIzANy6ntHGBJEuORt+0JDmkjQ9K/W6+eJcNLSsAgHNoGYFvDmyRqo9cpIAlVR92l2tOoFpWAABhgTAC39RUXrrMxcpdsmVF7pYVumwAIGIQRuCbxNRLl7lYuUC0rAAAwgphBL7pOdw9tqOJdg03h5Tc3V2uKf62rAAAwg5hBL6JinYPMlVTHS3n7k94svnBq/62rAAAwg5hBL7rP1ma9pKsxDTv48kZ0rSXLj4bxt+WFQBA2CGMoHX6T9apWds8d0/dtlaa/9Glp+X627ICAAg7hBG0XoPA4Mr0YSl4f1pWAABhh0XPYEb/yTqVOkYdnsuS5G5ZSbjyelpEACAC0TICc1rbsgIACCuEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgVORO7XU53Zux1VS6lx7vyWwOAABMiMwwsqfYvY19w91jkzPcK4Oy4BYAAEEVed00e4ql1TMab2NfXe4+vqfYTL0AAIhQkRVGXE53i4isJr557tjGB93lAABAUERWN82BLY1bRLxYUvVh7Xj3Df27W85FH6qu+qRuOvfvd/YeVWxye5+q4u/5oVAHzrf3+aFQB86P7PNDoQ6cf/58kyKrZaSmskXFYk8fa+OKAACAepEVRhJTW1SsLr5rG1cEAADUi6ww0nO4e9aMHE1+25JDpxPS9e8uQ4JbLwAAIlhkhZGoaPf0XUkXBhLr3P291zzEeiMAAARRZIURyb2OyLSXpKQ0r8NnEtL04fBf6ljmeEMVAwAgMkXWbJp6/SdLvcdIT2ZJknaOfF5fpH6TFhEAAAyIvJaReg2Cx5ddv0EQAQDAkMgNIwAAICQQRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEa1KowsXbpU2dnZio+PV05OjrZt29ai81555RU5HA5NmTKlNU8LAADCkM9hZNWqVcrPz1dBQYF27typQYMGafz48Tp69OhFz/vXv/6lH/7whxo5cmSrKwsAAMKPz2HkmWee0ezZszVr1iz1799fy5YtU/v27bV8+fJmz3E6nbrrrrv02GOPqXfv3n5VGAAAhBefwkhdXZ127Nih3Nzc8w8QFaXc3Fxt3bq12fMef/xxdevWTffdd1+LnufMmTOqrq72ugEAgPDkUxg5fvy4nE6nUlNTvY6npqaqoqKiyXPee+89vfDCC3r++edb/DyFhYVKSUnx3LKysnypJgAAsJE2nU1z4sQJfec739Hzzz+vLl26tPi8hQsXqqqqynM7ePBgG9YSAACY1M6Xwl26dFF0dLQqKyu9jldWViotLa1R+X379ulf//qXJk2a5DnmcrncT9yunfbu3avLL7+80XlxcXGKi4vzpWoAAMCmfGoZiY2N1eDBg7Vp0ybPMZfLpU2bNmnYsGGNyvfr108fffSRdu/e7blNnjxZY8eO1e7du+l+AQAAvrWMSFJ+fr5mzpypIUOGaOjQoVqyZIlqa2s1a9YsSdKMGTPUvXt3FRYWKj4+XgMGDPA6v2PHjpLU6DgAAIhMPoeR6dOn69ixY1q0aJEqKip09dVXa+PGjZ5BrWVlZYqKCv2FXZ1OKfrcvz/aHqN+35Cioy96CgAAaAM+hxFJmjt3rubOndvk90pKSi567ooVK1rzlAFVVCQ9+P9I/5jpvv/Q9y9T+07x+sHCExp54xmzlQMAIMKEfhNGgBUVSVOnSocPex8/fjRKjz+QotK3GTgLAEAwRVQYcTqlefMky2rim5ZDkvSbJ5PkdAa3XgAARLKICiOlpdKhQxcpYDl0rCJaH++ICVqdAACIdBEVRsrLW1bu82OMZAUAIFgiKoykp7esXOeu9NMAABAsERVGRo6UMjMlh6OZAg5LXdOcGjD4bFDrBQBAJIuoMBIdLT37rPvfjQKJwz2q9f4HT7DeCAAAQRRRYUSS8vKktWsbd9l0TXVp0S+qWGcEAIAga9WiZ3aXlyfljpb0nPv+o8/+WwO/GUOLCAAABkRcy0i96AY/+ZXXnCWIAABgSMSGEQAAEBoIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCqVWFk6dKlys7OVnx8vHJycrRt27Zmyz7//PMaOXKkLrvsMl122WXKzc29aHkAABBZfA4jq1atUn5+vgoKCrRz504NGjRI48eP19GjR5ssX1JSojvuuEObN2/W1q1blZWVpXHjxunw4cN+Vx4AANifz2HkmWee0ezZszVr1iz1799fy5YtU/v27bV8+fImy//xj3/UD37wA1199dXq16+ffve738nlcmnTpk1+Vx4AANifT2Gkrq5OO3bsUG5u7vkHiIpSbm6utm7d2qLHOHnypM6ePatOnTo1W+bMmTOqrq72ugEAgPDkUxg5fvy4nE6nUlNTvY6npqaqoqKiRY+xYMECZWRkeAWaCxUWFiolJcVzy8rK8qWaAADARoI6m+bJJ5/UK6+8onXr1ik+Pr7ZcgsXLlRVVZXndvDgwSDWEgAABFM7Xwp36dJF0dHRqqys9DpeWVmptLS0i577s5/9TE8++aT+/Oc/a+DAgRctGxcXp7i4OF+qBgAAbMqnlpHY2FgNHjzYa/Bp/WDUYcOGNXve008/rSeeeEIbN27UkCFDWl9bAAAQdnxqGZGk/Px8zZw5U0OGDNHQoUO1ZMkS1dbWatasWZKkGTNmqHv37iosLJQkPfXUU1q0aJFefvllZWdne8aWJCYmKjExMYA/CgAAsCOfw8j06dN17NgxLVq0SBUVFbr66qu1ceNGz6DWsrIyRUWdb3D5zW9+o7q6Ok2dOtXrcQoKCvToo4/6V3sAAGB7PocRSZo7d67mzp3b5PdKSkq87v/rX/9qzVMAAIAIwd40AADAqFa1jISb0X27qcNlHXw6p/bftdJGc+eHQh04397nh0IdOD+yzw+FOnD++fNNomUEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEa1KowsXbpU2dnZio+PV05OjrZt23bR8mvWrFG/fv0UHx+vq666Shs2bGhVZQEAQPjxOYysWrVK+fn5Kigo0M6dOzVo0CCNHz9eR48ebbL8li1bdMcdd+i+++7Trl27NGXKFE2ZMkUff/yx35UHAAD2187XE5555hnNnj1bs2bNkiQtW7ZMb7zxhpYvX64HH3ywUflnn31WEyZM0I9+9CNJ0hNPPKG3335bv/rVr7Rs2TI/q996LpdLrq8ckqSv/v2FvnKd9un8r6pOGj0/FOrA+fY+PxTqwPmRfX4o1IHzz5/vcrl8OjeQHJZlWS0tXFdXp/bt22vt2rWaMmWK5/jMmTP15Zdf6rXXXmt0To8ePZSfn6/58+d7jhUUFGj9+vX6v//7vyaf58yZMzpz5oznfnV1tbKyslRVVaXk5OSWVveiqj47qCM3jQvIYwEAYHcZG95SSu+sgD5mdXW1UlJSLvn326dumuPHj8vpdCo1NdXreGpqqioqKpo8p6KiwqfyklRYWKiUlBTPLSsrsC+OJHVIaR/wxwQAwK5M/l30uZsmGBYuXKj8/HzP/fqWkUCKuuwyXfHX9wL6mAAA2FXUZZcZe26fwkiXLl0UHR2tyspKr+OVlZVKS0tr8py0tDSfyktSXFyc4uLifKmaz6KiohTVuXObPgcAALg0n7ppYmNjNXjwYG3atMlzzOVyadOmTRo2bFiT5wwbNsyrvCS9/fbbzZYHAACRxedumvz8fM2cOVNDhgzR0KFDtWTJEtXW1npm18yYMUPdu3dXYWGhJGnevHkaPXq0fv7zn+vmm2/WK6+8ou3bt+u3v/1tYH8SAABgSz6HkenTp+vYsWNatGiRKioqdPXVV2vjxo2eQaplZWWKijrf4DJ8+HC9/PLLevjhh/XQQw/piiuu0Pr16zVgwIDA/RQAAMC2fJraa0pLpwYBAIDQ0SZTewEAAAKNMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwyufl4E2oXyS2urracE0AAEBL1f/dvtRi77YIIydOnJAkZWVlGa4JAADw1YkTJ5SSktLs922xN43L5dKRI0eUlJQkh8MRsMetrq5WVlaWDh48yJ43rcRr6B9eP//xGvqH189/vIbNsyxLJ06cUEZGhtcmuheyRctIVFSUMjMz2+zxk5OTeQP5idfQP7x+/uM19A+vn/94DZt2sRaRegxgBQAARhFGAACAUREdRuLi4lRQUKC4uDjTVbEtXkP/8Pr5j9fQP7x+/uM19J8tBrACAIDwFdEtIwAAwDzCCAAAMIowAgAAjCKMAAAAoyI6jCxdulTZ2dmKj49XTk6Otm3bZrpKtvDoo4/K4XB43fr162e6WiHt3Xff1aRJk5SRkSGHw6H169d7fd+yLC1atEjp6elKSEhQbm6u/vnPf5qpbIi61Gt4zz33NHpfTpgwwUxlQ1BhYaG+8Y1vKCkpSd26ddOUKVO0d+9erzKnT5/WnDlz1LlzZyUmJuq2225TZWWloRqHlpa8fmPGjGn0Hvz+979vqMb2ErFhZNWqVcrPz1dBQYF27typQYMGafz48Tp69KjpqtnClVdeqfLycs/tvffeM12lkFZbW6tBgwZp6dKlTX7/6aef1i9/+UstW7ZMH3zwgTp06KDx48fr9OnTQa5p6LrUayhJEyZM8Hpfrly5Mog1DG3vvPOO5syZo/fff19vv/22zp49q3Hjxqm2ttZT5oEHHtDrr7+uNWvW6J133tGRI0eUl5dnsNahoyWvnyTNnj3b6z349NNPG6qxzVgRaujQodacOXM8951Op5WRkWEVFhYarJU9FBQUWIMGDTJdDduSZK1bt85z3+VyWWlpadb//M//eI59+eWXVlxcnLVy5UoDNQx9F76GlmVZM2fOtG655RYj9bGjo0ePWpKsd955x7Is93suJibGWrNmjafM3//+d0uStXXrVlPVDFkXvn6WZVmjR4+25s2bZ65SNhaRLSN1dXXasWOHcnNzPceioqKUm5urrVu3GqyZffzzn/9URkaGevfurbvuuktlZWWmq2Rb+/fvV0VFhdf7MSUlRTk5ObwffVRSUqJu3bqpb9++uv/++/X555+brlLIqqqqkiR16tRJkrRjxw6dPXvW633Yr18/9ejRg/dhEy58/er98Y9/VJcuXTRgwAAtXLhQJ0+eNFE927HFRnmBdvz4cTmdTqWmpnodT01N1SeffGKoVvaRk5OjFStWqG/fviovL9djjz2mkSNH6uOPP1ZSUpLp6tlORUWFJDX5fqz/Hi5twoQJysvLU69evbRv3z499NBDmjhxorZu3aro6GjT1QspLpdL8+fP14gRIzRgwABJ7vdhbGysOnbs6FWW92FjTb1+knTnnXeqZ8+eysjI0IcffqgFCxZo7969KioqMlhbe4jIMAL/TJw40fPvgQMHKicnRz179tTq1at13333GawZItntt9/u+fdVV12lgQMH6vLLL1dJSYluuOEGgzULPXPmzNHHH3/MWK9Wau71+973vuf591VXXaX09HTdcMMN2rdvny6//PJgV9NWIrKbpkuXLoqOjm40SryyslJpaWmGamVfHTt21Ne+9jV9+umnpqtiS/XvOd6PgdW7d2916dKF9+UF5s6dqz/96U/avHmzMjMzPcfT0tJUV1enL7/80qs870Nvzb1+TcnJyZEk3oMtEJFhJDY2VoMHD9amTZs8x1wulzZt2qRhw4YZrJk91dTUaN++fUpPTzddFVvq1auX0tLSvN6P1dXV+uCDD3g/+uHQoUP6/PPPeV+eY1mW5s6dq3Xr1ukvf/mLevXq5fX9wYMHKyYmxut9uHfvXpWVlfE+1KVfv6bs3r1bkngPtkDEdtPk5+dr5syZGjJkiIYOHaolS5aotrZWs2bNMl21kPfDH/5QkyZNUs+ePXXkyBEVFBQoOjpad9xxh+mqhayamhqvq6P9+/dr9+7d6tSpk3r06KH58+frpz/9qa644gr16tVLjzzyiDIyMjRlyhRzlQ4xF3sNO3XqpMcee0y33Xab0tLStG/fPv34xz9Wnz59NH78eIO1Dh1z5szRyy+/rNdee01JSUmecSApKSlKSEhQSkqK7rvvPuXn56tTp05KTk7Wf/7nf2rYsGG67rrrDNfevEu9fvv27dPLL7+sm266SZ07d9aHH36oBx54QKNGjdLAgQMN194GTE/nMem5556zevToYcXGxlpDhw613n//fdNVsoXp06db6enpVmxsrNW9e3dr+vTp1qeffmq6WiFt8+bNlqRGt5kzZ1qW5Z7e+8gjj1ipqalWXFycdcMNN1h79+41W+kQc7HX8OTJk9a4ceOsrl27WjExMVbPnj2t2bNnWxUVFaarHTKaeu0kWb///e89ZU6dOmX94Ac/sC677DKrffv21q233mqVl5ebq3QIudTrV1ZWZo0aNcrq1KmTFRcXZ/Xp08f60Y9+ZFVVVZmtuE04LMuyghl+AAAAGorIMSMAACB0EEYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAY9f8Dpx0NhnWwP8oAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "rmse1, rmse2 = structure_prediction(\"ShiftML1.1\", frames, list_atom, list_cs)\n", + "sequence = np.array(range(len(rmse1)))\n", + "plt.stem(sequence, rmse2, 'b', label='GIPAW')\n", + "plt.stem(sequence, rmse1, 'tab:orange', label=\"ShiftML\")\n", + "plt.fill_between(sequence, 0.33+0.16, 0.33-0.16, alpha=0.3)\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "# atoms list defines groups of atoms to average over\n", + "\n", + "list_atom = ['1', '2', '3', '4', '5', '6', '7', '8', '9', '10,11',\n", + " '12', '13', '14', '15', '16,17', '18,19', '20', '21,22', '23,24,25,26,27,28,29,30,31']\n", + "\n", + "list_cs = np.array([6.92, 8.69, 9.01, 8.47, 15.37, 7.73, 9.64, 2.90, 1.78, 1.88, \n", + " 1.8, 1.6, 0.44, 1.54, 1.88, 0.8, 1, 1.74, 0.73])" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "urllib.request.urlretrieve(\"https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-06972-x/MediaObjects/41467_2018_6972_MOESM6_ESM.zip\", \"AZD.zip\")\n", + "\n", + "import zipfile\n", + "with zipfile.ZipFile(\"AZD.zip\",\"r\") as zip_ref:\n", + " zip_ref.extractall(\".\")\n", + "\n", + "frames = read(\"./Supplementary_Dataset_4/AZD_QuantumEspresso.xyz\", \":\")" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2024-07-30 13:34:45,537 - INFO - rascaline version: 0.1.0.dev558\n", + "2024-07-30 13:34:45,539 - INFO - rascaline-torch is installed, importing rascaline-torch\n", + "2024-07-30 13:34:45,539 - INFO - Found model version in url_resolve\n", + "2024-07-30 13:34:45,540 - INFO - Resolving model version to model files at url: https://tinyurl.com/f237evr3\n", + "2024-07-30 13:34:45,542 - INFO - Found ShiftML1.1 in cache, and importing it from here: /Users/zhangyuxuan/Library/Caches/shiftml/ShiftML1.1\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAykElEQVR4nO3df1iUVf7/8deAMqAC/iBAFEWzMjMRNf2ga2Jhaq5mZlrZYlZen0o/aWxt0g+pbZNsV5cuM/3oatrVulpmZmtruWykJX3MH+zWt7ItLUkBdU0QVNCZ+f4xMjoByuDAmWGej+uaK+bMued+30vpa899zrktDofDIQAAAEOCTBcAAAACG2EEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFHNTBdQF3a7XQcPHlR4eLgsFovpcgAAQB04HA4dP35ccXFxCgqqffzDL8LIwYMHFR8fb7oMAABQDwUFBerYsWOtn/tFGAkPD5fkvJiIiAjD1QAAgLooLS1VfHy86+/x2vhFGKm6NRMREUEYAQDAz1xsigUTWAEAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABG+cWmZ7gAu036YZtUViy1ipE6D5SCgk1XBQBAnRFG/NmXG6RNj0ulB8+1RcRJI+ZKPcaYqwsAAA9wm8ZffblBeiPNPYhIUmmhs/3LDWbqAgDAQ4QRf2S3OUdE5Kjhw7Ntm2Y5+wEA4OMII/7oh23VR0TcOKTSA85+AAD4OMKIPyor9m4/AAAMIoz4o1Yx3u0HAIBBhBF/1Hmgc9WMLLV0sEgRHZz9AADwcYQRfxQU7Fy+K6l6IDn7fsQL7DcCAPALhBF/1WOMNOE1KTzWvT0iztnOPiMAAD/Bpmf+rMcYqWuK9EK88/2ktdLlNzAiAu9gd18AjYQw4u/O/8uBvyzgLezuC6ARcZsGgDt29wXQyAgjAM5hd18ABngcRrZs2aLRo0crLi5OFotF69evr/Oxn3zyiZo1a6bevXt7eloAjYHdfQEY4HEYKS8vV2JiohYuXOjRcceOHVNaWppuvPFGT0+JQGC3Sfu2Sp+vdf6T/+dtBrv7AjDA4wmsI0eO1MiRIz0+0QMPPKC77rpLwcHBHo2mIAAwWdJ3sLsvAAMaZc7Iq6++qr179yozM7NO/SsqKlRaWur2QhPFZEnfwu6+AAxo8DDy73//W7NmzdLrr7+uZs3qNhCTlZWlyMhI1ys+Pr6Bq4QRTJb0PezuC8CABg0jNptNd911l5599lldeeWVdT4uIyNDJSUlrldBQUEDVgljmCzpm9jdF7g45rl5VYNuenb8+HHt2LFDu3fv1vTp0yVJdrtdDodDzZo10wcffKAbbrih2nFWq1VWq7UhS4MvYLKk72J3X6B2zHPzugYNIxEREfr888/d2l555RX94x//0Nq1a9WlS5eGPD18HZMlfRu7+wLVVc1z+/nt5ap5bowe1ovHYaSsrEzffvut6/2+ffuUn5+vtm3bqlOnTsrIyNCBAwf02muvKSgoSD179nQ7Pjo6WqGhodXaEYCqJkuWFqrmeSMW5+dMlgTgCy46z83inOfWfRTh3UMezxnZsWOHkpKSlJSUJElKT09XUlKSZs+eLUkqLCzU/v37vVslmiYmSwLwJ8xzazAej4ykpKTI4agpFTqtWLHigsc/88wzeuaZZzw9LZqqqsmSf/uNdLzwXHtEnDOIMNwJwFcwz63B8NRemMdkSQD+gHluDYYH5cE3MFkSgK9jU8AGQxgBAKAumOfWYAgjAOAP2GTLN7ApYINgzggA+Do22fItzHPzOkZGAMCX8TBJ38Q8N68ijACAr+JhkggQhBEA8FVssoUAQRgBAF/FJlsIEIQRAPBVbLKFAEEYAQBfxSZbCBCEEQDwVWyyhQBBGAEAX8YmWwgAbHoGAL6OTbbQxBFGAMAfsMkWGoLd5lwaXlbsnAht6N8twggAAIHIhx4zwJwRAAACjY89ZoAwAgBAIPHBxwwQRgAACCQ++JgBwggAAIHEBx8zQBgBACCQ+OBjBlhNA6Bp85Gli4DPqHrMQGmhap43YnF+3oiPGSCMAGi6fGjpIuAzqh4z8EaanI8VOD+QmHnMALdpADRNPrZ0EfApPvaYAcIIgKbHB5cuAj6nxxhp2vZz7yetlWZ+bmTUkDACoOnxwaWLgE/ykccMEEYAND0+uHQRQO0IIwBqZDvvDsaWre7vfZ4PLl0EUDvCCIBq1q2Tru5x7v3IkVJCgrPdL1QtXaxaGVCNRYro0KhLFwHUjjACwM26ddL48dKBA+7tBw442/0ikFQtXZRUPZCYWboIoHaEEQAuNps0Y4bkqGERSlXbzJl+csvGx5YuAqgdYQSAy9at0o8/1v65wyEVFDj7+QUfWroIoHaEEQAuhYXe7ecTfGTpIoDaEUYAuLRv791+AFAXhBEALoMHSx07SpZaFqFYLFJ8vLMfAHiLx2Fky5YtGj16tOLi4mSxWLR+/foL9l+3bp2GDRumyy67TBEREUpOTtb7779f33oBNKDgYOmll5w//zyQVL3Pznb2AwBv8TiMlJeXKzExUQsXLqxT/y1btmjYsGF67733tHPnTg0dOlSjR4/W7t27PS4WQMMbN05au7b6rZiOHZ3t48aZqQtA09XM0wNGjhypkSNH1rl/dna22/s5c+bonXfe0bvvvqukpCRPTw+gEYwbJ6UOkbTA+f7tddKNIxkRAarYbFLVfw5btkqDUvjv41I0+pwRu92u48ePq23btrX2qaioUGlpqdsLQOMKPu9Ph0GD+IMWqOL3OxT7oEYPI3/4wx9UVlamCRMm1NonKytLkZGRrld8fHwjVggAQM2axA7FPqhRw8iqVav07LPP6o033lB0dHSt/TIyMlRSUuJ6FRQUNGKVAABU16R2KPYxHs8Zqa/Vq1fr/vvv15tvvqnU1NQL9rVarbJarY1UGeBldpv0wzbn4+lbxbDRFtBEeLJDcUpKo5XVJDRKGPnLX/6ie++9V6tXr9aoUaMa45SAGV9ukDY9LpUePNcWEed8aBtbkAN+rUnuUOwjPL5NU1ZWpvz8fOXn50uS9u3bp/z8fO3fv1+S8xZLWlqaq/+qVauUlpamefPmacCAASoqKlJRUZFKSkq8cwWAr/hyg/RGmhznBxFJjtJC6Y005+cA/BY7FDccj8PIjh07lJSU5FqWm56erqSkJM2ePVuSVFhY6AomkrRkyRKdOXNG06ZNU/v27V2vGTNmeOkSAB9gt0mbHpdDjhoeWO+QQ5I2zXL2A+CX2KG44Xh8myYlJUWOmmbvnLVixQq397m5uZ6eAvA/P2yTSg9WCyJVLHJIpQec/brwJxXgj6p2KB4/nh2KvY1n0wBeYC8t9mo/AL6JHYobBmEE8IJ/7Y3xaj8AvmvcOGnnznPv314n7dtHELkUhBHAC74+OVAFJXGyO2q+UWN3WLS/pIO+PjmwkStDU3H+3hVbtrKXhWnsUOxdhBHAC2LjgjVj01xJqhZIqt7P3PSCYuP4EwueY/txNHWEEcALBg+WPisbo9vffE0Hj8e6ffZjaZxuf/M17Sgfwyx7eIztxxEIGm0HVqApOzfLfow2701R6Szn85RGvL5Wf993g+yOYK1dy1AuPHOx7cctFuf247fc4kf/brFDMWpAGAG8pGqW/ePp5/5g3bp/oOI6BCs7m8lt8FyT236cHYpRC27TAF7ELHt4U5PafvzsDsX62Q7FYodiiDACeB2z7OEtTWb78bM7FEs1bZh5to0digMaYQQAfFST2X787A7FtTtvh2IEJMIIAPioqonRkp9vP15Wx52H69oPTQ5hBAB8WJPYfrxVHXcerms/NDmEET/HroxA0+f3E6M7D3SumrnAoyQV0cHZDwGJMOLH2JURCBx+PTE6KNi5fFdS9UBy9v2IF9hvJIARRvwUuzIC8Cs9xkgTXpPC3XcoVkScs519RgIaYcQPXWxXRsm5KyO3bAD4lB5jpGnbz72ftFaa+TlBBIQRf+TJrowA4FPOvxXDVvA4izDih5rUrowAAGN8ZREEYcQPNZldGQEAxvjSIgjCiB9qMrsynsdX0jkABAJfWwRBGPFDTWZXxrN8KZ0DQFPni4sgCCN+qknsyijfS+cA0NT54iIIwogf8/ddGX0xnQNAU+eLiyAII37On3dl9MV0DgBNnS8ugiCMwBhfTOcA0NT54iIIwgiM8cV0jqaHlVqAO19cBEEYgTG+mM7RtLBSC6iZry2CIIzAGF9M52g6WKnlmxip8h2+tAiCMAKjfC2do2lgpZZvYqTK9/jKIgjCCIzzpXSOpoGVWr6HkSpcCGEEPsFX0jmaBlZq+RZGqnAxhBEATQ4rtXwLI1W4GMIIgCaHlVq+hZEqXAxhBECTw0ot38JIFS6GMAKgSWKllu9gpAoX43EY2bJli0aPHq24uDhZLBatX7/+osfk5uaqT58+slqt6tatm1asWFGPUgHAM6zU8g2MVOFiPA4j5eXlSkxM1MKFC+vUf9++fRo1apSGDh2q/Px8zZw5U/fff7/ef/99j4sFAE+xUss3MFKFC2nm6QEjR47UyJEj69x/8eLF6tKli+bNmydJuvrqq/Xxxx/rj3/8o4YPH+7p6QEAfmrcOCl1iKQFzvdvr5NuHElARCPMGcnLy1Nqaqpb2/Dhw5WXl1frMRUVFSotLXV7AQD8HyNVqEmDh5GioiLFxMS4tcXExKi0tFQnT56s8ZisrCxFRka6XvHx8Q1dJgAAMMQnV9NkZGSopKTE9SooKDBdEgAAaCAezxnxVGxsrIqLi93aiouLFRERobCwsBqPsVqtslqtDV0aAADwAQ0+MpKcnKycnBy3ts2bNys5ObmhTw0AAPyAx2GkrKxM+fn5ys/Pl+Rcupufn6/9+/dLct5iSUtLc/V/4IEHtHfvXv3mN7/R119/rVdeeUVvvPGGHnnkEe9cAQAA8Gseh5EdO3YoKSlJSUlJkqT09HQlJSVp9uzZkqTCwkJXMJGkLl26aOPGjdq8ebMSExM1b948/elPf2JZLwAAkFSPOSMpKSly1PQc6LNq2l01JSVFu3fv9vRUAAAgAPjkahoAABA4CCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjGpmugAAADxls9l0+vRpY+evqKxQcKv4cz+fCjZWy6W41Oto3ry5goMv/doJIwAAv+FwOFRUVKRjx44ZrcNutyto0Dznz0cKFXTUP280eOM6WrdurdjYWFkslnrXQRgBAPiNqiASHR2tFi1aXNJfgJfCdsam4KMVzp/bJii4mX+OjFzKdTgcDp04cUKHDh2SJLVv377eddQrjCxcuFC///3vVVRUpMTERC1YsED9+/evtX92drYWLVqk/fv3KyoqSuPHj1dWVpZCQ0PrXTgAILDYbDZXEGnXrp3ZWs7YFNzMGYRsoaH+HUYu4TrCwsIkSYcOHVJ0dHS9b9l4PB6zZs0apaenKzMzU7t27VJiYqKGDx/uSkY/t2rVKs2aNUuZmZn66quvtGzZMq1Zs0ZPPPFEvQoGAASmqjkiLVq0MFwJzlf1+7iUOTweh5H58+dr6tSpmjJlinr06KHFixerRYsWWr58eY39t23bpkGDBumuu+5SQkKCbrrpJt15553avn17vYsGAAQuU7dmUDNv/D48CiOVlZXauXOnUlNTz31BUJBSU1OVl5dX4zEDBw7Uzp07XeFj7969eu+993TzzTdfQtkAANSfzSbl5kp/+Yvznzab6YoCm0dzRo4cOSKbzaaYmBi39piYGH399dc1HnPXXXfpyJEj+sUvfiGHw6EzZ87ogQceuOBtmoqKClVUVLjel5aWelImAAC1WrdOmjFD+vHHc20dO0ovvSSNG2eurkDW4GuRcnNzNWfOHL3yyivatWuX1q1bp40bN+q5556r9ZisrCxFRka6XvHx8Q1dJgAgAKxbJ40f7x5EJOnAAWf7unUNe/6ioiLNmDFD3bp1U2hoqGJiYjRo0CAtWrRIJ06ckCQlJCQoOzvbdUxCQoIsFossFotatmypPn366M0333T73pMnT6pt27aKiopy+z/zixcvVnh4uM6cOeNqKysrU/PmzZWSkuL2Hbkf5cpisei7777z/oVfhEdhJCoqSsHBwSouLnZrLy4uVmxsbI3HPP300/rVr36l+++/X9dee61uvfVWzZkzR1lZWbLb7TUek5GRoZKSEteroKDAkzIBAKjGZnOOiDgc1T+raps5s+Fu2ezdu1dJSUn64IMPNGfOHO3evVt5eXn6zW9+o7/+9a/6+9//Xuuxv/3tb1VYWKjdu3fruuuu08SJE7Vt2zbX52+99ZauueYade/eXevXr3e1Dx06VGVlZdqxY4erbevWrYqNjdX//d//6dSpU6723NxcderUSZdffrl3L7wOPAojISEh6tu3r3JyclxtdrtdOTk5Sk5OrvGYEydOKCjI/TRVS38cNf0bIclqtSoiIsLtBQDApdi6tfqIyPkcDqmgwNmvITz00ENq1qyZduzYoQkTJujqq69W165ddcstt2jjxo0aPXp0rceGh4crNjZWV155pRYuXKiwsDC9++67rs+XLVumu+++W3fffbeWLVvmar/qqqvUvn175ebmutpyc3N1yy23qEuXLvr0/z51tX/00UcaOnSody+6jjy+TZOenq6lS5dq5cqV+uqrr/Tggw+qvLxcU6ZMkSSlpaUpIyPD1X/06NFatGiRVq9erX379mnz5s16+umnNXr0aK9sIQsAQF0UFnq3nyf+85//6IMPPtC0adPUsmXLGvvUdVVKs2bN1Lx5c1VWVkqSvvvuO+Xl5WnChAmaMGGCtm7dqh9++MHVf+jQofrwww9d7z/88EOlpKRoyJAhrpBy8uQpbd++3VgY8XjTs4kTJ+rw4cOaPXu2ioqK1Lt3b23atMk1qXX//v1uIyFPPfWULBaLnnrqKR04cECXXXaZRo8ereeff957VwEAwEXUdYPQS9hItFbffvutHA6HrrrqKrf2qKgo162SadOmae7cuRf8nsrKSs2bN08lJSW64YYbJEnLly/XyJEj1aZNG0nS8OHD9eqrr+qZZ56R5AwjM2fO1JkzZ3Ty5Ent3r1bQ4YM0enTp7Vo0SLpwXHK2/kvVVRU+E8YkaTp06dr+vTpNX52/lCQ5ExwmZmZyszMrM+pAADwisGDnatmDhyoed6IxeL8fPDgxqtp+/btstvtmjRpktvE0597/PHH9dRTT+nUqVNq1aqVXnjhBY0aNUo2m00rV67USy+95Op7991369FHH9Xs2bMVFBSklJQUlZeX67PPPtNPP/2kK6+8UpdddpmGDBmiKVOm6NSpCuXm7VTXrl3VqVOnxrjsang2DQAgIAQHO5fvjh/vDB7nB5KqOyTZ2c5+3tatWzdZLBbt2bPHrb1r166Szm2rXpvHHntM99xzj1q1aqWYmBjXLZ33339fBw4c0MSJE93622w25eTkaNiwYerWrZs6duyoDz/8UD/99JOGDBkiSYqLi1N8fLy27finPtz2mYammBkVkRphaS8AAL5i3Dhp7VopLs69vWNHZ3tD7TPSrl07DRs2TC+//LLKy8s9Pj4qKkrdunWr9nTcZcuW6Y477lB+fr7b64477nCbyDp06FDl5uYqNzfXbUnv4F8M1t8+/ETb8/9ftaW+jYmREQCA33I4HHKcPOnRMWNHSDfslGJig+SQRevXOzTsBruCgyX7ibp9h/2MTUEOh0dbob/yyisaNGiQ+vXrp2eeeUa9evVSUFCQPvvsM3399dfq27evR9dx+PBhvfvuu9qwYYN69uzp9llaWppuvfVWHT16VG3bttXQoUM1bdo0nT592jUyIknXX3+9Hn54uiorTxNGAACoD8fJk9rTx7O/xKvkV80lfVz6th7HX7VxuSxhdX/6/OWXX67du3drzpw5ysjI0I8//iir1aoePXro0Ucf1UMPPeTR+V977TW1bNlSN954Y7XPbrzxRoWFhen111/Xww8/rKFDh+rkyZPq3r272y7q119/vY6XleuqyxPUviFm7tYRYQQAgEbSvn17LViwQAsWLKi1z/fff3/B91V+/etf69e//nWNn4WEhOinn35yvU9ISKhxb6/OnTvLcWCXJMnk43kIIwAAv2UJC9NVu3Y2+nltZ2yyHP93o5+3qSKMAAD8lsVikaVFi0Y/r+OMTZayus8XwYWxmgYAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARrG0FwAQeOw26YdtUlmx1CpG6jxQCmqAJ+ShThgZAQAEli83SNk9pZW/lN66z/nP7J7OdoMsFovWr19f6+e5ubmyWCw6duyYq239+vXq1q2bgoODNXPmzAavsaEQRgAAgePLDdIbaVLpQff20kJnewMGksOHD+vBBx9Up06dZLVaFRsbq+HDh+uTTz6p0/EDBw5UYWGhIiMjXW3//d//rfHjx6ugoEDPPfec7rnnHo0dO7basRaLRRaLRZ9++qlbe0VFhdpdM1SWDn2U+1GuW/8LBSNv4zYNACAw2G3SpsclVX9Gi7PNIm2aJXUf1SC3bG677TZVVlZq5cqV6tq1q4qLi5WTk6P//Oc/dTo+JCREsbGxrvdlZWU6dOiQhg8frri4uIseHx8fr1dffVX/9V//5Wpbv369WrUM09FjJZ5fkBcxMgIACAw/bKs+IuLGIZUecPbzsmPHjmnr1q2aO3euhg4dqs6dO6t///7KyMjQmDFjXP2OHDmiW2+9VS1atNAVV1yhDRvOjdScf5smNzdX4eHhkqQbbrhBFotFKSkpWrlypd555x3XSEhubq7r+MmTJ2v16tU6efKkq+3VFa9q8u2jvX69niKMAAACQ1mxd/t5oFWrVmrVqpXWr1+vioqKWvs9++yzmjBhgv71r3/p5ptv1qRJk3T06NFq/QYOHKg9e/ZIkt566y0VFhZqw4YNmjBhgkaMGKHCwkIVFhZq4MCBrmP69u2rhIQEvfXWW5Kk/fv3a+vWrfrVbaO8fLWeI4wAAAJDqxjv9vNAs2bNtGLFCq1cuVKtW7fWoEGD9MQTT+hf//qXW7977rlHd955p7p166Y5c+aorKxM27dvr/Z9ISEhio6OliS1bdtWsbGxioiIUFhYmGs+SmxsrEJCQtyOu/fee7V8+XJJ0ooVKzRy5Ehd1q6N16/XU4QRAEBg6DxQioiTVNvTdi1SRAdnvwZw22236eDBg9qwYYNGjBih3Nxc9enTRytWrHD16dWrl+vnli1bKiIiQocOHfJaDXfffbfy8vK0d+9erVixQvfcc4/XvvtSEEYAAIEhKFgaMffsm58HkrPvR7zQoPuNhIaGatiwYXr66ae1bds23XPPPcrMzHR93rx5c/eqLBbZ7Xavnb9du3b65S9/qfvuu0+nTp3SyBEjvfbdl4IwAgAIHD3GSBNek8Jj3dsj4pztPcbUfFxDldOjh8rLy732fSEhIbLZbBfsc++99yo3N1dpaWkKDvaNjd5Y2gsA8FsOh0OO81aH1ElCqjTlI+mPVzrfT/iz1GWIc0TkxIk6fYX9jE1BDocsltpu+bj7z3/+o9tvv1333nuvevXqpfDwcO3YsUMvvviibrnlFs/qv4CEhAS9//772rNnj9q1a6fIyMhqoy0jRozQ4cOHFRERccHv2rdvn/Lz893arrjiCrVs2dJr9VYhjAAA/Jbj5Ent6dO3nkef3ZvjzcfqdfRVG5fLEhZap76tWrXSgAED9Mc//lHfffedTp8+rfj4eE2dOlVPPPFEvc5fk6lTpyo3N1f9+vVTWVmZPvzwQ6WkpLj1sVgsioqKkiTZztQ+ipKenl6tbevWrfrFL37htXqrEEYAAGhgVqtVWVlZysrKqrWPw1F9M7bzt35PSUlx69O6detqx1x22WX64IMP6vTdru+JDJfjwC7ZonvVqX9DIIwAAPyWJSxMV+3a2ejntZ2xyXL8341+3qaKMAIA8FsWi0WWFi0a/byOMzZZyuo2XwQXF7BhxGaTtm6VCgul9u2lwYMlH5lUDABAQAnIMLJunTRjhvTjj+faOnaUXnpJGjfOXF0AAASigNtnZN06afx49yAiSQcOONvXrTNTFwAAgSqgwojN5hwRqWmScFXbzJnOfgAA39TYKz1wYd74fQRUGNm6tfqIyPkcDqmgwNkPAOBbqjbvOlHHjcnQOKp+Hz/fXM0TATVnpLDQu/0AAI0nODhYrVu3dj04rkWLFnXeAdXbbGdsCj7jHBGwnTql4Gb+uQLiUq7D4XDoxIkTOnTokFq3bn1JW8sHVBhp3967/QAAjSs21vlMGW8+ybY+7Ha7gkoPO38+blVQkH/eaPDGdbRu3dr1e6mvgAojgwc7V80cOFDzvBGLxfn54MGNXxsA4OIsFovat2+v6OhonT592lgdJ0pOqMXfJjp//tUWtYhs/L1OvOFSr6N58+ZeedhevcLIwoUL9fvf/15FRUVKTEzUggUL1L9//1r7Hzt2TE8++aTWrVuno0ePqnPnzsrOztbNN99c78LrIzjYuXx3/Hhn8Dg/kFSN9GVns98IAPi64OBgo0+ctZ20KbSswPlziFWhoXV7Ro2v8ZXr8Hg8Zs2aNUpPT1dmZqZ27dqlxMREDR8+vNYhs8rKSg0bNkzff/+91q5dqz179mjp0qXq0KHDJRdfH+PGSWvXSnFx7u0dOzrb2WcEAIDG5fHIyPz58zV16lRNmTJFkrR48WJt3LhRy5cv16xZs6r1X758uY4ePapt27a5ZtomJCRcWtWXaNw4KTVViox0vn/vPemmmxgRAQDABI9GRiorK7Vz506lpqae+4KgIKWmpiovL6/GYzZs2KDk5GRNmzZNMTEx6tmzp+bMmSPbBTbzqKioUGlpqdvL284PHtdfTxABAMAUj8LIkSNHZLPZFBMT49YeExOjoqKiGo/Zu3ev1q5dK5vNpvfee09PP/205s2bp9/97ne1nicrK0uRkZGuV3x8vCdlAgAAP9Lga5Hsdruio6O1ZMkS9e3bVxMnTtSTTz6pxYsX13pMRkaGSkpKXK+CgoKGLhMAABji0ZyRqKgoBQcHq7i42K29uLi41jXG7du3r7b05+qrr1ZRUZEqKysVEhJS7Rir1Sqr1epJaQAAwE95NDISEhKivn37Kicnx9Vmt9uVk5Oj5OTkGo8ZNGiQvv32W9ntdlfbN998o/bt29cYRAAAQGDx+DZNenq6li5dqpUrV+qrr77Sgw8+qPLyctfqmrS0NGVkZLj6P/jggzp69KhmzJihb775Rhs3btScOXM0bdo0710FAADwWx4v7Z04caIOHz6s2bNnq6ioSL1799amTZtck1r379/vtp1sfHy83n//fT3yyCPq1auXOnTooBkzZujxxx/33lUAAAC/Va8dWKdPn67p06fX+Flubm61tuTkZH366af1ORUAAGji/PPJPgAAoMkgjAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAo5qZLsAX/OOrYoW1MF1F/VSWntDNZ3/+aM8hhUT454U0leuQms61cB2+hevwLU3xOkxiZAQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABgVuE/ttds0pPM2tQ8vVrsjITrRsZ8UFGy6KgAAAk5ghpEvNyjsvceVe89B5/tPpVNhsdqT9KQOdxxutjYAAAJM4N2m+XKD9EaaLGUH3ZqtJ4vVa9vDuuzH9w0VBgBAYKpXGFm4cKESEhIUGhqqAQMGaPv27XU6bvXq1bJYLBo7dmx9Tnvp7DZp0+OSHLL87COLHJKkq3bPcfYDAACNwuMwsmbNGqWnpyszM1O7du1SYmKihg8frkOHDl3wuO+//16PPvqoBg8eXO9iL9kP26TSg7V+bJFDoScL1ebIjkYsCgCAwOZxGJk/f76mTp2qKVOmqEePHlq8eLFatGih5cuX13qMzWbTpEmT9Oyzz6pr166XVPAlKSuuU7eQU4cbuBAAAFDFozBSWVmpnTt3KjU19dwXBAUpNTVVeXl5tR7329/+VtHR0brvvvvqdJ6KigqVlpa6vbyiVUydulWGXuad8wEAgIvyKIwcOXJENptNMTHuf6nHxMSoqKioxmM+/vhjLVu2TEuXLq3zebKyshQZGel6xcfHe1Jm7ToPlCLipGozRpwcsuhUWHv9FNXPO+cDAAAX1aCraY4fP65f/epXWrp0qaKioup8XEZGhkpKSlyvgoIC7xQUFCyNmCup+hTWqvd7kp5gvxEAABqRR/uMREVFKTg4WMXF7nMviouLFRsbW63/d999p++//16jR492tdntdueJmzXTnj17dPnll1c7zmq1ymq1elJa3fUYI014TY73fiNLWaGruSIsVnuSnmCfEQAAGplHIyMhISHq27evcnJyXG12u105OTlKTk6u1r979+76/PPPlZ+f73qNGTNGQ4cOVX5+vvduv3iqxxidnHJuOfKn1y3Vx6P+QRABAMAAj3dgTU9P1+TJk9WvXz/1799f2dnZKi8v15QpUyRJaWlp6tChg7KyshQaGqqePXu6Hd+6dWtJqtbe6M67FXO03XUK4dYMAABGeBxGJk6cqMOHD2v27NkqKipS7969tWnTJtek1v379ysoKPA2dgUAAPVTr2fTTJ8+XdOnT6/xs9zc3Aseu2LFivqcEgAANFEMYQAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCqXmFk4cKFSkhIUGhoqAYMGKDt27fX2nfp0qUaPHiw2rRpozZt2ig1NfWC/QEAQGDxOIysWbNG6enpyszM1K5du5SYmKjhw4fr0KFDNfbPzc3VnXfeqQ8//FB5eXmKj4/XTTfdpAMHDlxy8QAAwP95HEbmz5+vqVOnasqUKerRo4cWL16sFi1aaPny5TX2//Of/6yHHnpIvXv3Vvfu3fWnP/1JdrtdOTk5l1w8AADwfx6FkcrKSu3cuVOpqannviAoSKmpqcrLy6vTd5w4cUKnT59W27Zta+1TUVGh0tJStxcAAGiaPAojR44ckc1mU0xMjFt7TEyMioqK6vQdjz/+uOLi4twCzc9lZWUpMjLS9YqPj/ekTAAA4EcadTXNCy+8oNWrV+vtt99WaGhorf0yMjJUUlLiehUUFDRilQAAoDE186RzVFSUgoODVVxc7NZeXFys2NjYCx77hz/8QS+88IL+/ve/q1evXhfsa7VaZbVaPSkNAAD4KY9GRkJCQtS3b1+3yadVk1GTk5NrPe7FF1/Uc889p02bNqlfv371rxYAADQ5Ho2MSFJ6eromT56sfv36qX///srOzlZ5ebmmTJkiSUpLS1OHDh2UlZUlSZo7d65mz56tVatWKSEhwTW3pFWrVmrVqpUXLwUAAPgjj8PIxIkTdfjwYc2ePVtFRUXq3bu3Nm3a5JrUun//fgUFnRtwWbRokSorKzV+/Hi378nMzNQzzzxzadUDAAC/53EYkaTp06dr+vTpNX6Wm5vr9v7777+vzykAAECAqFcYaWqGXBWtlm1ami6jXsp/Kpc2OX/mOnxDU7kWrsO3cB2+pSleh0k8KA8AABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABG1SuMLFy4UAkJCQoNDdWAAQO0ffv2C/Z/88031b17d4WGhuraa6/Ve++9V69iAQBA0+NxGFmzZo3S09OVmZmpXbt2KTExUcOHD9ehQ4dq7L9t2zbdeeeduu+++7R7926NHTtWY8eO1RdffHHJxQMAAP/XzNMD5s+fr6lTp2rKlCmSpMWLF2vjxo1avny5Zs2aVa3/Sy+9pBEjRuixxx6TJD333HPavHmzXn75ZS1evPgSy68/u90u+xmLJOnMT0d1xn7KWC2X4kzJCa7DxzSVa+E6fAvX4Vua4nXY7XZjdVgcDoejrp0rKyvVokULrV27VmPHjnW1T548WceOHdM777xT7ZhOnTopPT1dM2fOdLVlZmZq/fr1+uc//1njeSoqKlRRUeF6X1paqvj4eJWUlCgiIqKu5V5Qyd4CHbz5Jq98FwAA/i7uvQ8U2TXeq99ZWlqqyMjIi/797dFtmiNHjshmsykmJsatPSYmRkVFRTUeU1RU5FF/ScrKylJkZKTrFR/v3f9xJKllZAuvfycAAP7K5N+LHt+maQwZGRlKT093va8aGfGmoDZtdMUnH3v1OwEA8FdBbdoYO7dHYSQqKkrBwcEqLi52ay8uLlZsbGyNx8TGxnrUX5KsVqusVqsnpXksKChIQe3aNeg5AADAxXl0myYkJER9+/ZVTk6Oq81utysnJ0fJyck1HpOcnOzWX5I2b95ca38AABBYPL5Nk56ersmTJ6tfv37q37+/srOzVV5e7lpdk5aWpg4dOigrK0uSNGPGDA0ZMkTz5s3TqFGjtHr1au3YsUNLlizx7pUAAAC/5HEYmThxog4fPqzZs2erqKhIvXv31qZNm1yTVPfv36+goHMDLgMHDtSqVav01FNP6YknntAVV1yh9evXq2fPnt67CgAA4Lc8WtprSl2XBgEAAN/RIEt7AQAAvI0wAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADDK4+3gTajaJLa0tNRwJQAAoK6q/t6+2GbvfhFGjh8/LkmKj483XAkAAPDU8ePHFRkZWevnfvFsGrvdroMHDyo8PFwWi8Vr31taWqr4+HgVFBTwzBsfwO/D9/A78S38PnwLv4+LczgcOn78uOLi4tweovtzfjEyEhQUpI4dOzbY90dERPAvkg/h9+F7+J34Fn4fvoXfx4VdaESkChNYAQCAUYQRAABgVECHEavVqszMTFmtVtOlQPw+fBG/E9/C78O38PvwHr+YwAoAAJqugB4ZAQAA5hFGAACAUYQRAABgFGEEAAAYFdBhZOHChUpISFBoaKgGDBig7du3my4pIGVlZem6665TeHi4oqOjNXbsWO3Zs8d0WTjrhRdekMVi0cyZM02XErAOHDigu+++W+3atVNYWJiuvfZa7dixw3RZActms+npp59Wly5dFBYWpssvv1zPPffcRZ+/gtoFbBhZs2aN0tPTlZmZqV27dikxMVHDhw/XoUOHTJcWcD766CNNmzZNn376qTZv3qzTp0/rpptuUnl5uenSAt5nn32m//3f/1WvXr1MlxKwfvrpJw0aNEjNmzfX3/72N3355ZeaN2+e2rRpY7q0gDV37lwtWrRIL7/8sr766ivNnTtXL774ohYsWGC6NL8VsEt7BwwYoOuuu04vv/yyJOfzb+Lj4/U///M/mjVrluHqAtvhw4cVHR2tjz76SNdff73pcgJWWVmZ+vTpo1deeUW/+93v1Lt3b2VnZ5suK+DMmjVLn3zyibZu3Wq6FJz1y1/+UjExMVq2bJmr7bbbblNYWJhef/11g5X5r4AcGamsrNTOnTuVmprqagsKClJqaqry8vIMVgZJKikpkSS1bdvWcCWBbdq0aRo1apTbfydofBs2bFC/fv10++23Kzo6WklJSVq6dKnpsgLawIEDlZOTo2+++UaS9M9//lMff/yxRo4cabgy/+UXD8rztiNHjshmsykmJsatPSYmRl9//bWhqiA5R6hmzpypQYMGqWfPnqbLCVirV6/Wrl279Nlnn5kuJeDt3btXixYtUnp6up544gl99tlnevjhhxUSEqLJkyebLi8gzZo1S6WlperevbuCg4Nls9n0/PPPa9KkSaZL81sBGUbgu6ZNm6YvvvhCH3/8selSAlZBQYFmzJihzZs3KzQ01HQ5Ac9ut6tfv36aM2eOJCkpKUlffPGFFi9eTBgx5I033tCf//xnrVq1Stdcc43y8/M1c+ZMxcXF8Tupp4AMI1FRUQoODlZxcbFbe3FxsWJjYw1VhenTp+uvf/2rtmzZoo4dO5ouJ2Dt3LlThw4dUp8+fVxtNptNW7Zs0csvv6yKigoFBwcbrDCwtG/fXj169HBru/rqq/XWW28ZqgiPPfaYZs2apTvuuEOSdO211+qHH35QVlYWYaSeAnLOSEhIiPr27aucnBxXm91uV05OjpKTkw1WFpgcDoemT5+ut99+W//4xz/UpUsX0yUFtBtvvFGff/658vPzXa9+/fpp0qRJys/PJ4g0skGDBlVb6v7NN9+oc+fOhirCiRMnFBTk/tdncHCw7Ha7oYr8X0COjEhSenq6Jk+erH79+ql///7Kzs5WeXm5pkyZYrq0gDNt2jStWrVK77zzjsLDw1VUVCRJioyMVFhYmOHqAk94eHi1+TotW7ZUu3btmMdjwCOPPKKBAwdqzpw5mjBhgrZv364lS5ZoyZIlpksLWKNHj9bzzz+vTp066ZprrtHu3bs1f/583XvvvaZL81+OALZgwQJHp06dHCEhIY7+/fs7Pv30U9MlBSRJNb5effVV06XhrCFDhjhmzJhhuoyA9e677zp69uzpsFqtju7duzuWLFliuqSAVlpa6pgxY4ajU6dOjtDQUEfXrl0dTz75pKOiosJ0aX4rYPcZAQAAviEg54wAAADfQRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABg1P8HoEQyJHxqw2sAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "rmse1, rmse2 = structure_prediction(\"ShiftML1.1\", frames, list_atom, list_cs)\n", + "sequence = np.array(range(len(rmse1)))\n", + "plt.stem(sequence, rmse2, 'b', label='GIPAW')\n", + "plt.stem(sequence, rmse1, 'tab:orange', label=\"ShiftML\")\n", + "plt.fill_between(sequence, 0.33+0.16, 0.33-0.16, alpha=0.3)\n", + "plt.legend()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "DevShiftML", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}