From cda70b0fde6dd00a79cb84bec86c10afe1a1c47d Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 18:36:00 -0500 Subject: [PATCH 01/11] Functional version of evcurve workflow --- atomistics/workflows/evcurve/workflow.py | 63 +++++++++++++++++------- 1 file changed, 44 insertions(+), 19 deletions(-) diff --git a/atomistics/workflows/evcurve/workflow.py b/atomistics/workflows/evcurve/workflow.py index 89eb3b5e..339467b6 100644 --- a/atomistics/workflows/evcurve/workflow.py +++ b/atomistics/workflows/evcurve/workflow.py @@ -1,6 +1,7 @@ import numpy as np from ase.atoms import Atoms from collections import OrderedDict +from typing import Optional, List from atomistics.shared.output import OutputEnergyVolumeCurve from atomistics.workflows.evcurve.fit import EnergyVolumeFit @@ -103,6 +104,39 @@ def fit_ev_curve( ).fit_dict +def get_strains( + vol_range: Optional[float], + num_points: Optional[int], + strain_lst: Optional[List[float]] = None, +) -> np.ndarray: + if strain_lst is None: + return np.linspace( + -vol_range, + vol_range, + int(num_points), + ) + return strain_lst + + +def generate_structures_helper( + structure: Atoms, + vol_range: Optional[float], + num_points: Optional[int], + strain_lst: Optional[List[float]] = None, + axes: tuple[str, str, str] = ("x", "y", "z"), +) -> dict: + strain_lst = get_strains( + vol_range=vol_range, num_points=num_points, strain_lst=strain_lst + ) + return { + 1 + + np.round(strain, 7): _strain_axes( + structure=structure, axes=axes, volume_strain=strain + ) + for strain in strain_lst + } + + class EnergyVolumeCurveProperties: def __init__(self, fit_module): self._fit_module = fit_module @@ -165,8 +199,12 @@ def generate_structures(self) -> dict: (dict) """ self._structure_dict = OrderedDict( - generate_structure_dict( - structure=self.structure, axes=self.axes, strain_lst=self._get_strains() + generate_structures_helper( + structure=self.structure, + vol_range=self.vol_range, + num_points=self.num_points, + strain_lst=self.strains, + axes=self.axes, ) ) return {"calc_energy": self._structure_dict} @@ -213,21 +251,8 @@ def get_thermal_properties( ) def _get_strains(self): - strains = self.strains - if strains is None: - strains = np.linspace( - -self.vol_range, - self.vol_range, - int(self.num_points), - ) - return strains - - -def generate_structure_dict(structure, axes, strain_lst): - return { - 1 - + np.round(strain, 7): _strain_axes( - structure=structure, axes=axes, volume_strain=strain + return get_strains( + vol_range=self.vol_range, + num_points=self.num_points, + strain_lst=self.strains, ) - for strain in strain_lst - } From 66a2061fb3c9a748e9970ef101d12e2a206e370d Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 19:10:31 -0500 Subject: [PATCH 02/11] Test workflows in functional form --- .../workflows/elastic/elastic_moduli.py | 7 ++ atomistics/workflows/elastic/workflow.py | 2 +- atomistics/workflows/evcurve/workflow.py | 5 + tests/test_elastic_lammps.py | 118 ++++++++++++++++-- tests/test_evcurve_lammps.py | 57 +++++++++ 5 files changed, 175 insertions(+), 14 deletions(-) diff --git a/atomistics/workflows/elastic/elastic_moduli.py b/atomistics/workflows/elastic/elastic_moduli.py index 41c04da9..d5fb09c7 100644 --- a/atomistics/workflows/elastic/elastic_moduli.py +++ b/atomistics/workflows/elastic/elastic_moduli.py @@ -2,6 +2,8 @@ import numpy as np +from atomistics.shared.output import OutputElastic + def get_bulkmodul_voigt(elastic_matrix: np.ndarray) -> float: return ( @@ -214,3 +216,8 @@ def AVR(self) -> float: @cache def elastic_matrix_eigval(self): return get_elastic_matrix_eigval(elastic_matrix=self.elastic_matrix()) + + def to_dict(self, output_keys: tuple = OutputElastic.keys()) -> dict: + return OutputElastic( + **{k: getattr(self, k) for k in OutputElastic.keys()} + ).get(output_keys=output_keys) \ No newline at end of file diff --git a/atomistics/workflows/elastic/workflow.py b/atomistics/workflows/elastic/workflow.py index 8d8de8b5..aa6f7bdc 100644 --- a/atomistics/workflows/elastic/workflow.py +++ b/atomistics/workflows/elastic/workflow.py @@ -47,7 +47,7 @@ def generate_structures(self) -> dict: def analyse_structures( self, output_dict: dict, output_keys: tuple = OutputElastic.keys() - ): + ) -> dict: """ Args: diff --git a/atomistics/workflows/evcurve/workflow.py b/atomistics/workflows/evcurve/workflow.py index 339467b6..00f2e293 100644 --- a/atomistics/workflows/evcurve/workflow.py +++ b/atomistics/workflows/evcurve/workflow.py @@ -166,6 +166,11 @@ def fit_dict(self) -> dict: if k in self._fit_module.fit_dict.keys() } + def to_dict(self, output_keys: tuple = OutputEnergyVolumeCurve.keys()) -> dict: + return OutputEnergyVolumeCurve( + **{k: getattr(self, k) for k in OutputEnergyVolumeCurve.keys()} + ).get(output_keys=output_keys) + class EnergyVolumeCurveWorkflow(Workflow): def __init__( diff --git a/tests/test_elastic_lammps.py b/tests/test_elastic_lammps.py index 21618edb..44d5974a 100644 --- a/tests/test_elastic_lammps.py +++ b/tests/test_elastic_lammps.py @@ -5,6 +5,11 @@ import unittest from atomistics.workflows import ElasticMatrixWorkflow, optimize_positions_and_volume +from atomistics.workflows.elastic.workflow import ( + ElasticProperties, + analyse_structures_helper, + generate_structures_helper +) try: from atomistics.calculators import ( @@ -109,16 +114,103 @@ def test_calc_elastic(self): self.assertTrue(np.isclose(elastic_dict['youngsmodul_hill'], 101.45869947879392)) self.assertTrue(np.isclose(elastic_dict['poissonsratio_hill'], 0.2842453510798992)) self.assertTrue(np.isclose(elastic_dict['AVR'], 4.962492964955925)) - # self.assertTrue(np.isclose(elastic_dict['C_eigval'].eigenvalues, np.array( - # [235.12517572, 53.59208765, 53.59208765, 51.23853765, 51.23853765, 51.23853765] - # )).all()) - # self.assertTrue(np.isclose(elastic_dict['C_eigval'].eigenvectors, np.array( - # [ - # [-0.57735027, -0.62664396, 0.28808397, 0., 0.,0.], - # [-0.57735027, 0.76662983, 0.51758912, 0., 0., 0.], - # [-0.57735027, -0.13998587, -0.80567309, 0., 0., 0.], - # [0., 0., 0., 1., 0., 0.], - # [0., 0., 0., 0., 1., 0.], - # [0., 0., 0., 0., 0., 1.], - # ] - # )).all()) + + def test_calc_elastic_functions(self): + structure = bulk("Al", cubic=True) + df_pot_selected = get_potential_by_name( + potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', + resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), + ) + result_dict = evaluate_with_lammps( + task_dict={"optimize_positions_and_volume": structure}, + potential_dataframe=df_pot_selected, + ) + sym_dict, structure_dict = generate_structures_helper( + structure=result_dict["structure_with_optimized_positions_and_volume"], + eps_range=0.005, + num_of_point=5, + zero_strain_job_name="s_e_0", + sqrt_eta=True, + ) + result_dict = evaluate_with_lammps( + task_dict={"calc_energy": structure_dict}, + potential_dataframe=df_pot_selected, + ) + elastic_matrix, A2, strain_energy, ene0 = analyse_structures_helper( + output_dict=result_dict, + Lag_strain_list=sym_dict["Lag_strain_list"], + epss=sym_dict["epss"], + v0=sym_dict["v0"], + LC=sym_dict["LC"], + fit_order=2, + zero_strain_job_name="s_e_0", + ) + sym_dict["strain_energy"] = strain_energy + sym_dict["e0"] = ene0 + sym_dict["A2"] = A2 + elastic_dict = ElasticProperties(elastic_matrix=elastic_matrix).to_dict() + self.assertTrue(np.isclose(elastic_matrix, np.array( + [ + [114.10311701, 60.51102935, 60.51102935, 0., 0., 0.], + [60.51102935, 114.10311701, 60.51102935, 0., 0., 0.], + [60.51102935, 60.51102935, 114.10311701, 0., 0., 0.], + [0., 0., 0., 51.23853765, 0., 0.], + [0., 0., 0., 0., 51.23853765, 0.], + [0., 0., 0., 0., 0., 51.23853765] + ]) + ).all()) + self.assertEqual(sym_dict['SGN'], 225) + self.assertEqual(sym_dict['LC'], 'CI') + self.assertEqual(sym_dict['Lag_strain_list'], ['01', '08', '23']) + self.assertTrue(np.isclose(sym_dict['epss'], np.array([-0.005, -0.0025, 0., 0.0025, 0.005])).all()) + self.assertTrue(np.isclose(sym_dict["v0"], 66.43035441556098)) + self.assertTrue(np.isclose(sym_dict["e0"], -13.439999952735112)) + self.assertTrue(np.isclose(sym_dict['strain_energy'], np.array( + [ + [ + (-0.005, -13.436320248980278), + (-0.0025, -13.439079680886989), + (0.0, -13.439999952735112), + (0.0024999999999999996, -13.439084974614394), + (0.005, -13.436364320399795) + ], + [ + (-0.005, -13.43817471490433), + (-0.0025, -13.439544638502628), + (0.0, -13.439999952735112), + (0.0024999999999999996, -13.43954822781134), + (0.005, -13.438204192615181) + ], + [ + (-0.005, -13.436741954502294), + (-0.0025, -13.439195465714551), + (0.0, -13.439999952735112), + (0.0024999999999999996, -13.439213491269701), + (0.005, -13.436885713447486) + ] + ])).all() + ) + self.assertTrue(np.isclose(A2, np.array([2.20130388, 1.08985578, 1.91883479])).all()) + self.assertTrue(np.isclose(elastic_dict['bulkmodul_voigt'], 78.37505857279467)) + self.assertTrue(np.isclose(elastic_dict['shearmodul_voigt'], 41.46154012284969)) + self.assertTrue(np.isclose(elastic_dict['youngsmodul_voigt'], 105.73882997912072)) + self.assertTrue(np.isclose(elastic_dict['poissonsratio_voigt'], 0.2751435386362729)) + self.assertTrue(np.isclose(elastic_dict['elastic_matrix_inverse'], np.array( + [ + [0.01385733, -0.00480214, -0.00480214, 0., 0., 0.], + [-0.00480214, 0.01385733, -0.00480214, 0., 0., 0.], + [-0.00480214, -0.00480214, 0.01385733, 0., 0., 0.], + [0., 0., 0., 0.01951656, 0., 0.], + [0., 0., 0., 0., 0.01951656, 0.], + [0., 0., 0., 0., 0., 0.01951656] + ] + )).all()) + self.assertTrue(np.isclose(elastic_dict['bulkmodul_reuss'], 78.37505857279469)) + self.assertTrue(np.isclose(elastic_dict['shearmodul_reuss'], 37.54104251720356)) + self.assertTrue(np.isclose(elastic_dict['youngsmodul_reuss'], 97.11702764970639)) + self.assertTrue(np.isclose(elastic_dict['poissonsratio_reuss'], 0.29347803281170937)) + self.assertTrue(np.isclose(elastic_dict['bulkmodul_hill'], 78.37505857279467)) + self.assertTrue(np.isclose(elastic_dict['shearmodul_hill'], 39.501291320026624)) + self.assertTrue(np.isclose(elastic_dict['youngsmodul_hill'], 101.45869947879392)) + self.assertTrue(np.isclose(elastic_dict['poissonsratio_hill'], 0.2842453510798992)) + self.assertTrue(np.isclose(elastic_dict['AVR'], 4.962492964955925)) diff --git a/tests/test_evcurve_lammps.py b/tests/test_evcurve_lammps.py index a57da44d..d564894b 100644 --- a/tests/test_evcurve_lammps.py +++ b/tests/test_evcurve_lammps.py @@ -5,6 +5,14 @@ import unittest from atomistics.workflows import EnergyVolumeCurveWorkflow, optimize_positions_and_volume +from atomistics.workflows.evcurve.workflow import ( + EnergyVolumeCurveProperties, + fit_ev_curve_internal, + generate_structures_helper, + get_energy_lst, + get_volume_lst, + get_thermal_properties, +) try: @@ -58,3 +66,52 @@ def test_calc_evcurve(self): self.assertEqual(len(temperatures_ev), 2) self.assertEqual(len(volumes_ev), 2) self.assertTrue(volumes_ev[0] < volumes_ev[-1]) + + def test_calc_evcurve_functional(self): + structure = bulk("Al", cubic=True) + df_pot_selected = get_potential_by_name( + potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', + resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), + ) + result_dict = evaluate_with_lammps( + task_dict={"optimize_positions_and_volume": structure}, + potential_dataframe=df_pot_selected, + ) + structure_dict = generate_structures_helper( + structure=result_dict["structure_with_optimized_positions_and_volume"], + vol_range=0.05, + num_points=11, + strain_lst=None, + axes=('x', 'y', 'z'), + ) + result_dict = evaluate_with_lammps( + task_dict={"calc_energy": structure_dict}, + potential_dataframe=df_pot_selected, + ) + fit_dict = EnergyVolumeCurveProperties( + fit_module=fit_ev_curve_internal( + volume_lst=get_volume_lst(structure_dict=structure_dict), + energy_lst=get_energy_lst( + output_dict=result_dict, structure_dict=structure_dict + ), + fit_type='polynomial', + fit_order=3, + ) + ).to_dict() + thermal_properties_dict = get_thermal_properties( + fit_dict=fit_dict, + masses=structure.get_masses(), + t_min=1.0, + t_max=1500.0, + t_step=50.0, + temperatures=[100, 1000], + constant_volume=False, + output_keys=["temperatures", "volumes"], + ) + temperatures_ev, volumes_ev = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] + self.assertTrue(np.isclose(fit_dict['volume_eq'], 66.43019853103964)) + self.assertTrue(np.isclose(fit_dict['bulkmodul_eq'], 77.7250135953191)) + self.assertTrue(np.isclose(fit_dict['b_prime_eq'], 1.2795467367276832)) + self.assertEqual(len(temperatures_ev), 2) + self.assertEqual(len(volumes_ev), 2) + self.assertTrue(volumes_ev[0] < volumes_ev[-1]) From 69c6d24a4201d3aa3d3eb8813dfaf7d56cfd6451 Mon Sep 17 00:00:00 2001 From: pyiron-runner Date: Wed, 24 Apr 2024 00:11:51 +0000 Subject: [PATCH 03/11] Format black --- atomistics/workflows/elastic/elastic_moduli.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/atomistics/workflows/elastic/elastic_moduli.py b/atomistics/workflows/elastic/elastic_moduli.py index d5fb09c7..1350d162 100644 --- a/atomistics/workflows/elastic/elastic_moduli.py +++ b/atomistics/workflows/elastic/elastic_moduli.py @@ -218,6 +218,6 @@ def elastic_matrix_eigval(self): return get_elastic_matrix_eigval(elastic_matrix=self.elastic_matrix()) def to_dict(self, output_keys: tuple = OutputElastic.keys()) -> dict: - return OutputElastic( - **{k: getattr(self, k) for k in OutputElastic.keys()} - ).get(output_keys=output_keys) \ No newline at end of file + return OutputElastic(**{k: getattr(self, k) for k in OutputElastic.keys()}).get( + output_keys=output_keys + ) From fc70350c1cbc191552c1a2e6977b66b1e91c444f Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 19:35:50 -0500 Subject: [PATCH 04/11] Test workflows in functional form --- .../workflows/elastic/elastic_moduli.py | 6 +-- atomistics/workflows/elastic/helper.py | 26 +++++++++++++ atomistics/workflows/elastic/workflow.py | 17 ++------- atomistics/workflows/evcurve/workflow.py | 37 +++++++++++++------ tests/test_evcurve_lammps.py | 21 ++++------- 5 files changed, 65 insertions(+), 42 deletions(-) diff --git a/atomistics/workflows/elastic/elastic_moduli.py b/atomistics/workflows/elastic/elastic_moduli.py index d5fb09c7..1350d162 100644 --- a/atomistics/workflows/elastic/elastic_moduli.py +++ b/atomistics/workflows/elastic/elastic_moduli.py @@ -218,6 +218,6 @@ def elastic_matrix_eigval(self): return get_elastic_matrix_eigval(elastic_matrix=self.elastic_matrix()) def to_dict(self, output_keys: tuple = OutputElastic.keys()) -> dict: - return OutputElastic( - **{k: getattr(self, k) for k in OutputElastic.keys()} - ).get(output_keys=output_keys) \ No newline at end of file + return OutputElastic(**{k: getattr(self, k) for k in OutputElastic.keys()}).get( + output_keys=output_keys + ) diff --git a/atomistics/workflows/elastic/helper.py b/atomistics/workflows/elastic/helper.py index 2a20aff7..930796ac 100644 --- a/atomistics/workflows/elastic/helper.py +++ b/atomistics/workflows/elastic/helper.py @@ -4,6 +4,8 @@ import numpy as np import scipy.constants +from atomistics.shared.output import OutputElastic +from atomistics.workflows.elastic.elastic_moduli import ElasticProperties from atomistics.workflows.elastic.symmetry import ( symmetry_analysis, get_C_from_A2, @@ -89,6 +91,30 @@ def generate_structures_helper( def analyse_structures_helper( + output_dict: dict, + sym_dict: dict, + fit_order: int = 2, + zero_strain_job_name: str = "s_e_0", + output_keys: tuple = OutputElastic.keys(), +): + elastic_matrix, A2, strain_energy, ene0 = _get_elastic_matrix( + output_dict=output_dict, + Lag_strain_list=sym_dict["Lag_strain_list"], + epss=sym_dict["epss"], + v0=sym_dict["v0"], + LC=sym_dict["LC"], + fit_order=fit_order, + zero_strain_job_name=zero_strain_job_name, + ) + sym_dict["strain_energy"] = strain_energy + sym_dict["e0"] = ene0 + sym_dict["A2"] = A2 + return sym_dict, ElasticProperties(elastic_matrix=elastic_matrix).to_dict( + output_keys=output_keys + ) + + +def _get_elastic_matrix( output_dict: dict, Lag_strain_list: list[float], epss: float, diff --git a/atomistics/workflows/elastic/workflow.py b/atomistics/workflows/elastic/workflow.py index aa6f7bdc..d23a68e2 100644 --- a/atomistics/workflows/elastic/workflow.py +++ b/atomistics/workflows/elastic/workflow.py @@ -3,7 +3,6 @@ from atomistics.shared.output import OutputElastic from atomistics.workflows.interface import Workflow -from atomistics.workflows.elastic.elastic_moduli import ElasticProperties from atomistics.workflows.elastic.helper import ( generate_structures_helper, analyse_structures_helper, @@ -57,19 +56,11 @@ def analyse_structures( Returns: """ - elastic_matrix, A2, strain_energy, ene0 = analyse_structures_helper( + self._data, elastic_dict = analyse_structures_helper( output_dict=output_dict, - Lag_strain_list=self._data["Lag_strain_list"], - epss=self._data["epss"], - v0=self._data["v0"], - LC=self._data["LC"], + sym_dict=self._data, fit_order=self.fit_order, zero_strain_job_name=self.zero_strain_job_name, + output_keys=output_keys, ) - self._data["strain_energy"] = strain_energy - self._data["e0"] = ene0 - self._data["A2"] = A2 - elastic = ElasticProperties(elastic_matrix=elastic_matrix) - return OutputElastic( - **{k: getattr(elastic, k) for k in OutputElastic.keys()} - ).get(output_keys=output_keys) + return elastic_dict diff --git a/atomistics/workflows/evcurve/workflow.py b/atomistics/workflows/evcurve/workflow.py index 00f2e293..27c1bd1f 100644 --- a/atomistics/workflows/evcurve/workflow.py +++ b/atomistics/workflows/evcurve/workflow.py @@ -137,6 +137,25 @@ def generate_structures_helper( } +def analyse_structures_helper( + output_dict: dict, + structure_dict: dict, + fit_type: str = "polynomial", + fit_order: int = 3, + output_keys: tuple = OutputEnergyVolumeCurve.keys(), +) -> dict: + return EnergyVolumeCurveProperties( + fit_module=fit_ev_curve_internal( + volume_lst=get_volume_lst(structure_dict=structure_dict), + energy_lst=get_energy_lst( + output_dict=output_dict, structure_dict=structure_dict + ), + fit_type=fit_type, + fit_order=fit_order, + ) + ).to_dict(output_keys=output_keys) + + class EnergyVolumeCurveProperties: def __init__(self, fit_module): self._fit_module = fit_module @@ -217,19 +236,13 @@ def generate_structures(self) -> dict: def analyse_structures( self, output_dict: dict, output_keys: tuple = OutputEnergyVolumeCurve.keys() ) -> dict: - evcurve = EnergyVolumeCurveProperties( - fit_module=fit_ev_curve_internal( - volume_lst=get_volume_lst(structure_dict=self._structure_dict), - energy_lst=get_energy_lst( - output_dict=output_dict, structure_dict=self._structure_dict - ), - fit_type=self.fit_type, - fit_order=self.fit_order, - ) + self._fit_dict = analyse_structures_helper( + output_dict=output_dict, + structure_dict=self._structure_dict, + fit_type=self.fit_type, + fit_order=self.fit_order, + output_keys=output_keys, ) - self._fit_dict = OutputEnergyVolumeCurve( - **{k: getattr(evcurve, k) for k in OutputEnergyVolumeCurve.keys()} - ).get(output_keys=output_keys) return self.fit_dict def get_volume_lst(self) -> np.ndarray: diff --git a/tests/test_evcurve_lammps.py b/tests/test_evcurve_lammps.py index d564894b..b2075613 100644 --- a/tests/test_evcurve_lammps.py +++ b/tests/test_evcurve_lammps.py @@ -6,11 +6,8 @@ from atomistics.workflows import EnergyVolumeCurveWorkflow, optimize_positions_and_volume from atomistics.workflows.evcurve.workflow import ( - EnergyVolumeCurveProperties, - fit_ev_curve_internal, + analyse_structures_helper, generate_structures_helper, - get_energy_lst, - get_volume_lst, get_thermal_properties, ) @@ -88,16 +85,12 @@ def test_calc_evcurve_functional(self): task_dict={"calc_energy": structure_dict}, potential_dataframe=df_pot_selected, ) - fit_dict = EnergyVolumeCurveProperties( - fit_module=fit_ev_curve_internal( - volume_lst=get_volume_lst(structure_dict=structure_dict), - energy_lst=get_energy_lst( - output_dict=result_dict, structure_dict=structure_dict - ), - fit_type='polynomial', - fit_order=3, - ) - ).to_dict() + fit_dict = analyse_structures_helper( + output_dict=result_dict, + structure_dict=structure_dict, + fit_type="polynomial", + fit_order=3 + ) thermal_properties_dict = get_thermal_properties( fit_dict=fit_dict, masses=structure.get_masses(), From 27ed2523df1b60591410659f9f84048beb8b6eb8 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 19:35:56 -0500 Subject: [PATCH 05/11] fix test --- tests/test_elastic_lammps.py | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/tests/test_elastic_lammps.py b/tests/test_elastic_lammps.py index 44d5974a..83c844f7 100644 --- a/tests/test_elastic_lammps.py +++ b/tests/test_elastic_lammps.py @@ -6,7 +6,6 @@ from atomistics.workflows import ElasticMatrixWorkflow, optimize_positions_and_volume from atomistics.workflows.elastic.workflow import ( - ElasticProperties, analyse_structures_helper, generate_structures_helper ) @@ -136,20 +135,13 @@ def test_calc_elastic_functions(self): task_dict={"calc_energy": structure_dict}, potential_dataframe=df_pot_selected, ) - elastic_matrix, A2, strain_energy, ene0 = analyse_structures_helper( + sym_dict, elastic_dict = analyse_structures_helper( output_dict=result_dict, - Lag_strain_list=sym_dict["Lag_strain_list"], - epss=sym_dict["epss"], - v0=sym_dict["v0"], - LC=sym_dict["LC"], + sym_dict=sym_dict, fit_order=2, zero_strain_job_name="s_e_0", ) - sym_dict["strain_energy"] = strain_energy - sym_dict["e0"] = ene0 - sym_dict["A2"] = A2 - elastic_dict = ElasticProperties(elastic_matrix=elastic_matrix).to_dict() - self.assertTrue(np.isclose(elastic_matrix, np.array( + self.assertTrue(np.isclose(elastic_dict["elastic_matrix"], np.array( [ [114.10311701, 60.51102935, 60.51102935, 0., 0., 0.], [60.51102935, 114.10311701, 60.51102935, 0., 0., 0.], @@ -190,7 +182,7 @@ def test_calc_elastic_functions(self): ] ])).all() ) - self.assertTrue(np.isclose(A2, np.array([2.20130388, 1.08985578, 1.91883479])).all()) + self.assertTrue(np.isclose(sym_dict["A2"], np.array([2.20130388, 1.08985578, 1.91883479])).all()) self.assertTrue(np.isclose(elastic_dict['bulkmodul_voigt'], 78.37505857279467)) self.assertTrue(np.isclose(elastic_dict['shearmodul_voigt'], 41.46154012284969)) self.assertTrue(np.isclose(elastic_dict['youngsmodul_voigt'], 105.73882997912072)) From 0ef06ffae9230505125ec54026cd86c4e6679612 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 19:41:07 -0500 Subject: [PATCH 06/11] introduce evcurve helper module --- atomistics/workflows/evcurve/helper.py | 185 +++++++++++++++++++++++ atomistics/workflows/evcurve/workflow.py | 185 +---------------------- tests/test_evcurve_lammps.py | 4 +- 3 files changed, 193 insertions(+), 181 deletions(-) create mode 100644 atomistics/workflows/evcurve/helper.py diff --git a/atomistics/workflows/evcurve/helper.py b/atomistics/workflows/evcurve/helper.py new file mode 100644 index 00000000..89d843f7 --- /dev/null +++ b/atomistics/workflows/evcurve/helper.py @@ -0,0 +1,185 @@ +import numpy as np +from ase.atoms import Atoms +from typing import Optional, List + +from atomistics.shared.output import OutputEnergyVolumeCurve +from atomistics.workflows.evcurve.fit import EnergyVolumeFit + + +def _strain_axes( + structure: Atoms, volume_strain: float, axes: tuple[str, str, str] = ("x", "y", "z") +) -> Atoms: + """ + Strain box along given axes to achieve given *volumetric* strain. + + Returns a copy. + """ + axes = np.array([a in axes for a in ("x", "y", "z")]) + num_axes = sum(axes) + # formula calculates the strain along each axis to achieve the overall volumetric strain + # beware that: (1+e)**x - 1 != e**x + strains = axes * ((1 + volume_strain) ** (1.0 / num_axes) - 1) + return apply_strain(structure=structure, epsilon=strains, return_box=True) + + +def apply_strain( + structure: Atoms, epsilon: float, return_box: bool = False, mode: str = "linear" +) -> Atoms: + """ + Apply a given strain on the structure. It applies the matrix `F` in the manner: + + ``` + new_cell = F @ current_cell + ``` + + Args: + epsilon (float/list/ndarray): epsilon matrix. If a single number is set, the same + strain is applied in each direction. If a 3-dim vector is set, it will be + multiplied by a unit matrix. + return_box (bool): whether to return a box. If set to True, only the returned box will + have the desired strain and the original box will stay unchanged. + mode (str): `linear` or `lagrangian`. If `linear`, `F` is equal to the epsilon - 1. + If `lagrangian`, epsilon is given by `(F^T * F - 1) / 2`. It raises an error if + the strain is not symmetric (if the shear components are given). + """ + epsilon = np.array([epsilon]).flatten() + if len(epsilon) == 3 or len(epsilon) == 1: + epsilon = epsilon * np.eye(3) + epsilon = epsilon.reshape(3, 3) + if epsilon.min() < -1.0: + raise ValueError("Strain value too negative") + if return_box: + structure_copy = structure.copy() + else: + structure_copy = structure + cell = structure_copy.cell.copy() + if mode == "linear": + F = epsilon + np.eye(3) + elif mode == "lagrangian": + if not np.allclose(epsilon, epsilon.T): + raise ValueError("Strain must be symmetric if `mode = 'lagrangian'`") + E, V = np.linalg.eigh(2 * epsilon + np.eye(3)) + F = np.einsum("ik,k,jk->ij", V, np.sqrt(E), V) + else: + raise ValueError("mode must be `linear` or `lagrangian`") + cell = np.matmul(F, cell) + structure_copy.set_cell(cell, scale_atoms=True) + if return_box: + return structure_copy + + +def get_energy_lst(output_dict: dict, structure_dict: dict) -> list: + return [output_dict["energy"][k] for k in structure_dict.keys()] + + +def get_volume_lst(structure_dict: dict) -> list: + return [structure.get_volume() for structure in structure_dict.values()] + + +def fit_ev_curve_internal( + volume_lst: np.ndarray, energy_lst: np.ndarray, fit_type: str, fit_order: int +) -> EnergyVolumeFit: + fit_module = EnergyVolumeFit( + volume_lst=volume_lst, + energy_lst=energy_lst, + ) + fit_module.fit(fit_type=fit_type, fit_order=fit_order) + return fit_module + + +def fit_ev_curve( + volume_lst: np.ndarray, energy_lst: np.ndarray, fit_type: str, fit_order: int +) -> dict: + return fit_ev_curve_internal( + volume_lst=volume_lst, + energy_lst=energy_lst, + fit_type=fit_type, + fit_order=fit_order, + ).fit_dict + + +def get_strains( + vol_range: Optional[float], + num_points: Optional[int], + strain_lst: Optional[List[float]] = None, +) -> np.ndarray: + if strain_lst is None: + return np.linspace( + -vol_range, + vol_range, + int(num_points), + ) + return strain_lst + + +def generate_structures_helper( + structure: Atoms, + vol_range: Optional[float], + num_points: Optional[int], + strain_lst: Optional[List[float]] = None, + axes: tuple[str, str, str] = ("x", "y", "z"), +) -> dict: + strain_lst = get_strains( + vol_range=vol_range, num_points=num_points, strain_lst=strain_lst + ) + return { + 1 + + np.round(strain, 7): _strain_axes( + structure=structure, axes=axes, volume_strain=strain + ) + for strain in strain_lst + } + + +def analyse_structures_helper( + output_dict: dict, + structure_dict: dict, + fit_type: str = "polynomial", + fit_order: int = 3, + output_keys: tuple = OutputEnergyVolumeCurve.keys(), +) -> dict: + return EnergyVolumeCurveProperties( + fit_module=fit_ev_curve_internal( + volume_lst=get_volume_lst(structure_dict=structure_dict), + energy_lst=get_energy_lst( + output_dict=output_dict, structure_dict=structure_dict + ), + fit_type=fit_type, + fit_order=fit_order, + ) + ).to_dict(output_keys=output_keys) + + +class EnergyVolumeCurveProperties: + def __init__(self, fit_module): + self._fit_module = fit_module + + def volume_eq(self) -> float: + return self._fit_module.fit_dict["volume_eq"] + + def energy_eq(self) -> float: + return self._fit_module.fit_dict["energy_eq"] + + def bulkmodul_eq(self) -> float: + return self._fit_module.fit_dict["bulkmodul_eq"] + + def b_prime_eq(self) -> float: + return self._fit_module.fit_dict["b_prime_eq"] + + def volume(self) -> np.ndarray: + return self._fit_module.fit_dict["volume"] + + def energy(self) -> np.ndarray: + return self._fit_module.fit_dict["energy"] + + def fit_dict(self) -> dict: + return { + k: self._fit_module.fit_dict[k] + for k in ["fit_type", "least_square_error", "poly_fit", "fit_order"] + if k in self._fit_module.fit_dict.keys() + } + + def to_dict(self, output_keys: tuple = OutputEnergyVolumeCurve.keys()) -> dict: + return OutputEnergyVolumeCurve( + **{k: getattr(self, k) for k in OutputEnergyVolumeCurve.keys()} + ).get(output_keys=output_keys) diff --git a/atomistics/workflows/evcurve/workflow.py b/atomistics/workflows/evcurve/workflow.py index 27c1bd1f..c3c5a4f0 100644 --- a/atomistics/workflows/evcurve/workflow.py +++ b/atomistics/workflows/evcurve/workflow.py @@ -10,185 +10,12 @@ get_thermal_properties, OutputThermodynamic, ) - - -def _strain_axes( - structure: Atoms, volume_strain: float, axes: tuple[str, str, str] = ("x", "y", "z") -) -> Atoms: - """ - Strain box along given axes to achieve given *volumetric* strain. - - Returns a copy. - """ - axes = np.array([a in axes for a in ("x", "y", "z")]) - num_axes = sum(axes) - # formula calculates the strain along each axis to achieve the overall volumetric strain - # beware that: (1+e)**x - 1 != e**x - strains = axes * ((1 + volume_strain) ** (1.0 / num_axes) - 1) - return apply_strain(structure=structure, epsilon=strains, return_box=True) - - -def apply_strain( - structure: Atoms, epsilon: float, return_box: bool = False, mode: str = "linear" -) -> Atoms: - """ - Apply a given strain on the structure. It applies the matrix `F` in the manner: - - ``` - new_cell = F @ current_cell - ``` - - Args: - epsilon (float/list/ndarray): epsilon matrix. If a single number is set, the same - strain is applied in each direction. If a 3-dim vector is set, it will be - multiplied by a unit matrix. - return_box (bool): whether to return a box. If set to True, only the returned box will - have the desired strain and the original box will stay unchanged. - mode (str): `linear` or `lagrangian`. If `linear`, `F` is equal to the epsilon - 1. - If `lagrangian`, epsilon is given by `(F^T * F - 1) / 2`. It raises an error if - the strain is not symmetric (if the shear components are given). - """ - epsilon = np.array([epsilon]).flatten() - if len(epsilon) == 3 or len(epsilon) == 1: - epsilon = epsilon * np.eye(3) - epsilon = epsilon.reshape(3, 3) - if epsilon.min() < -1.0: - raise ValueError("Strain value too negative") - if return_box: - structure_copy = structure.copy() - else: - structure_copy = structure - cell = structure_copy.cell.copy() - if mode == "linear": - F = epsilon + np.eye(3) - elif mode == "lagrangian": - if not np.allclose(epsilon, epsilon.T): - raise ValueError("Strain must be symmetric if `mode = 'lagrangian'`") - E, V = np.linalg.eigh(2 * epsilon + np.eye(3)) - F = np.einsum("ik,k,jk->ij", V, np.sqrt(E), V) - else: - raise ValueError("mode must be `linear` or `lagrangian`") - cell = np.matmul(F, cell) - structure_copy.set_cell(cell, scale_atoms=True) - if return_box: - return structure_copy - - -def get_energy_lst(output_dict: dict, structure_dict: dict) -> list: - return [output_dict["energy"][k] for k in structure_dict.keys()] - - -def get_volume_lst(structure_dict: dict) -> list: - return [structure.get_volume() for structure in structure_dict.values()] - - -def fit_ev_curve_internal( - volume_lst: np.ndarray, energy_lst: np.ndarray, fit_type: str, fit_order: int -) -> EnergyVolumeFit: - fit_module = EnergyVolumeFit( - volume_lst=volume_lst, - energy_lst=energy_lst, - ) - fit_module.fit(fit_type=fit_type, fit_order=fit_order) - return fit_module - - -def fit_ev_curve( - volume_lst: np.ndarray, energy_lst: np.ndarray, fit_type: str, fit_order: int -) -> dict: - return fit_ev_curve_internal( - volume_lst=volume_lst, - energy_lst=energy_lst, - fit_type=fit_type, - fit_order=fit_order, - ).fit_dict - - -def get_strains( - vol_range: Optional[float], - num_points: Optional[int], - strain_lst: Optional[List[float]] = None, -) -> np.ndarray: - if strain_lst is None: - return np.linspace( - -vol_range, - vol_range, - int(num_points), - ) - return strain_lst - - -def generate_structures_helper( - structure: Atoms, - vol_range: Optional[float], - num_points: Optional[int], - strain_lst: Optional[List[float]] = None, - axes: tuple[str, str, str] = ("x", "y", "z"), -) -> dict: - strain_lst = get_strains( - vol_range=vol_range, num_points=num_points, strain_lst=strain_lst - ) - return { - 1 - + np.round(strain, 7): _strain_axes( - structure=structure, axes=axes, volume_strain=strain - ) - for strain in strain_lst - } - - -def analyse_structures_helper( - output_dict: dict, - structure_dict: dict, - fit_type: str = "polynomial", - fit_order: int = 3, - output_keys: tuple = OutputEnergyVolumeCurve.keys(), -) -> dict: - return EnergyVolumeCurveProperties( - fit_module=fit_ev_curve_internal( - volume_lst=get_volume_lst(structure_dict=structure_dict), - energy_lst=get_energy_lst( - output_dict=output_dict, structure_dict=structure_dict - ), - fit_type=fit_type, - fit_order=fit_order, - ) - ).to_dict(output_keys=output_keys) - - -class EnergyVolumeCurveProperties: - def __init__(self, fit_module): - self._fit_module = fit_module - - def volume_eq(self) -> float: - return self._fit_module.fit_dict["volume_eq"] - - def energy_eq(self) -> float: - return self._fit_module.fit_dict["energy_eq"] - - def bulkmodul_eq(self) -> float: - return self._fit_module.fit_dict["bulkmodul_eq"] - - def b_prime_eq(self) -> float: - return self._fit_module.fit_dict["b_prime_eq"] - - def volume(self) -> np.ndarray: - return self._fit_module.fit_dict["volume"] - - def energy(self) -> np.ndarray: - return self._fit_module.fit_dict["energy"] - - def fit_dict(self) -> dict: - return { - k: self._fit_module.fit_dict[k] - for k in ["fit_type", "least_square_error", "poly_fit", "fit_order"] - if k in self._fit_module.fit_dict.keys() - } - - def to_dict(self, output_keys: tuple = OutputEnergyVolumeCurve.keys()) -> dict: - return OutputEnergyVolumeCurve( - **{k: getattr(self, k) for k in OutputEnergyVolumeCurve.keys()} - ).get(output_keys=output_keys) +from atomistics.workflows.evcurve.helper import ( + get_strains, + get_volume_lst, + generate_structures_helper, + analyse_structures_helper +) class EnergyVolumeCurveWorkflow(Workflow): diff --git a/tests/test_evcurve_lammps.py b/tests/test_evcurve_lammps.py index b2075613..930c865a 100644 --- a/tests/test_evcurve_lammps.py +++ b/tests/test_evcurve_lammps.py @@ -5,10 +5,10 @@ import unittest from atomistics.workflows import EnergyVolumeCurveWorkflow, optimize_positions_and_volume -from atomistics.workflows.evcurve.workflow import ( +from atomistics.workflows.evcurve.debye import get_thermal_properties +from atomistics.workflows.evcurve.helper import ( analyse_structures_helper, generate_structures_helper, - get_thermal_properties, ) From 3e7ad510cc1513f377aff221729ff4780e373956 Mon Sep 17 00:00:00 2001 From: pyiron-runner Date: Wed, 24 Apr 2024 00:43:57 +0000 Subject: [PATCH 07/11] Format black --- atomistics/workflows/evcurve/workflow.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atomistics/workflows/evcurve/workflow.py b/atomistics/workflows/evcurve/workflow.py index c3c5a4f0..0cdeb58b 100644 --- a/atomistics/workflows/evcurve/workflow.py +++ b/atomistics/workflows/evcurve/workflow.py @@ -14,7 +14,7 @@ get_strains, get_volume_lst, generate_structures_helper, - analyse_structures_helper + analyse_structures_helper, ) From 36ef285ecaa3f0a2f924613d056d18b71f4ce121 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 19:45:52 -0500 Subject: [PATCH 08/11] fix quasi harmonic --- atomistics/workflows/quasiharmonic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/atomistics/workflows/quasiharmonic.py b/atomistics/workflows/quasiharmonic.py index 38f2453c..cc982555 100644 --- a/atomistics/workflows/quasiharmonic.py +++ b/atomistics/workflows/quasiharmonic.py @@ -2,8 +2,8 @@ import numpy as np from atomistics.shared.output import OutputThermodynamic -from atomistics.workflows.evcurve.workflow import ( - EnergyVolumeCurveWorkflow, +from atomistics.workflows.evcurve.workflow import EnergyVolumeCurveWorkflow +from atomistics.workflows.evcurve.helper import ( fit_ev_curve, _strain_axes, ) From d38b6c14427f6c555581dafc3d6af603107eb015 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 19:56:44 -0500 Subject: [PATCH 09/11] split tests --- tests/test_elastic_lammps.py | 97 -------------------- tests/test_elastic_lammps_functional.py | 117 ++++++++++++++++++++++++ tests/test_evcurve_lammps.py | 51 ----------- tests/test_evcurve_lammps_functional.py | 71 ++++++++++++++ 4 files changed, 188 insertions(+), 148 deletions(-) create mode 100644 tests/test_elastic_lammps_functional.py create mode 100644 tests/test_evcurve_lammps_functional.py diff --git a/tests/test_elastic_lammps.py b/tests/test_elastic_lammps.py index 83c844f7..a3cff96f 100644 --- a/tests/test_elastic_lammps.py +++ b/tests/test_elastic_lammps.py @@ -5,10 +5,6 @@ import unittest from atomistics.workflows import ElasticMatrixWorkflow, optimize_positions_and_volume -from atomistics.workflows.elastic.workflow import ( - analyse_structures_helper, - generate_structures_helper -) try: from atomistics.calculators import ( @@ -113,96 +109,3 @@ def test_calc_elastic(self): self.assertTrue(np.isclose(elastic_dict['youngsmodul_hill'], 101.45869947879392)) self.assertTrue(np.isclose(elastic_dict['poissonsratio_hill'], 0.2842453510798992)) self.assertTrue(np.isclose(elastic_dict['AVR'], 4.962492964955925)) - - def test_calc_elastic_functions(self): - structure = bulk("Al", cubic=True) - df_pot_selected = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', - resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), - ) - result_dict = evaluate_with_lammps( - task_dict={"optimize_positions_and_volume": structure}, - potential_dataframe=df_pot_selected, - ) - sym_dict, structure_dict = generate_structures_helper( - structure=result_dict["structure_with_optimized_positions_and_volume"], - eps_range=0.005, - num_of_point=5, - zero_strain_job_name="s_e_0", - sqrt_eta=True, - ) - result_dict = evaluate_with_lammps( - task_dict={"calc_energy": structure_dict}, - potential_dataframe=df_pot_selected, - ) - sym_dict, elastic_dict = analyse_structures_helper( - output_dict=result_dict, - sym_dict=sym_dict, - fit_order=2, - zero_strain_job_name="s_e_0", - ) - self.assertTrue(np.isclose(elastic_dict["elastic_matrix"], np.array( - [ - [114.10311701, 60.51102935, 60.51102935, 0., 0., 0.], - [60.51102935, 114.10311701, 60.51102935, 0., 0., 0.], - [60.51102935, 60.51102935, 114.10311701, 0., 0., 0.], - [0., 0., 0., 51.23853765, 0., 0.], - [0., 0., 0., 0., 51.23853765, 0.], - [0., 0., 0., 0., 0., 51.23853765] - ]) - ).all()) - self.assertEqual(sym_dict['SGN'], 225) - self.assertEqual(sym_dict['LC'], 'CI') - self.assertEqual(sym_dict['Lag_strain_list'], ['01', '08', '23']) - self.assertTrue(np.isclose(sym_dict['epss'], np.array([-0.005, -0.0025, 0., 0.0025, 0.005])).all()) - self.assertTrue(np.isclose(sym_dict["v0"], 66.43035441556098)) - self.assertTrue(np.isclose(sym_dict["e0"], -13.439999952735112)) - self.assertTrue(np.isclose(sym_dict['strain_energy'], np.array( - [ - [ - (-0.005, -13.436320248980278), - (-0.0025, -13.439079680886989), - (0.0, -13.439999952735112), - (0.0024999999999999996, -13.439084974614394), - (0.005, -13.436364320399795) - ], - [ - (-0.005, -13.43817471490433), - (-0.0025, -13.439544638502628), - (0.0, -13.439999952735112), - (0.0024999999999999996, -13.43954822781134), - (0.005, -13.438204192615181) - ], - [ - (-0.005, -13.436741954502294), - (-0.0025, -13.439195465714551), - (0.0, -13.439999952735112), - (0.0024999999999999996, -13.439213491269701), - (0.005, -13.436885713447486) - ] - ])).all() - ) - self.assertTrue(np.isclose(sym_dict["A2"], np.array([2.20130388, 1.08985578, 1.91883479])).all()) - self.assertTrue(np.isclose(elastic_dict['bulkmodul_voigt'], 78.37505857279467)) - self.assertTrue(np.isclose(elastic_dict['shearmodul_voigt'], 41.46154012284969)) - self.assertTrue(np.isclose(elastic_dict['youngsmodul_voigt'], 105.73882997912072)) - self.assertTrue(np.isclose(elastic_dict['poissonsratio_voigt'], 0.2751435386362729)) - self.assertTrue(np.isclose(elastic_dict['elastic_matrix_inverse'], np.array( - [ - [0.01385733, -0.00480214, -0.00480214, 0., 0., 0.], - [-0.00480214, 0.01385733, -0.00480214, 0., 0., 0.], - [-0.00480214, -0.00480214, 0.01385733, 0., 0., 0.], - [0., 0., 0., 0.01951656, 0., 0.], - [0., 0., 0., 0., 0.01951656, 0.], - [0., 0., 0., 0., 0., 0.01951656] - ] - )).all()) - self.assertTrue(np.isclose(elastic_dict['bulkmodul_reuss'], 78.37505857279469)) - self.assertTrue(np.isclose(elastic_dict['shearmodul_reuss'], 37.54104251720356)) - self.assertTrue(np.isclose(elastic_dict['youngsmodul_reuss'], 97.11702764970639)) - self.assertTrue(np.isclose(elastic_dict['poissonsratio_reuss'], 0.29347803281170937)) - self.assertTrue(np.isclose(elastic_dict['bulkmodul_hill'], 78.37505857279467)) - self.assertTrue(np.isclose(elastic_dict['shearmodul_hill'], 39.501291320026624)) - self.assertTrue(np.isclose(elastic_dict['youngsmodul_hill'], 101.45869947879392)) - self.assertTrue(np.isclose(elastic_dict['poissonsratio_hill'], 0.2842453510798992)) - self.assertTrue(np.isclose(elastic_dict['AVR'], 4.962492964955925)) diff --git a/tests/test_elastic_lammps_functional.py b/tests/test_elastic_lammps_functional.py new file mode 100644 index 00000000..3227df67 --- /dev/null +++ b/tests/test_elastic_lammps_functional.py @@ -0,0 +1,117 @@ +import os + +from ase.build import bulk +import numpy as np +import unittest + +from atomistics.workflows.elastic.workflow import ( + analyse_structures_helper, + generate_structures_helper +) + +try: + from atomistics.calculators import ( + evaluate_with_lammps, get_potential_by_name + ) + + skip_lammps_test = False +except ImportError: + skip_lammps_test = True + + +@unittest.skipIf( + skip_lammps_test, "LAMMPS is not installed, so the LAMMPS tests are skipped." +) +class TestElastic(unittest.TestCase): + def test_calc_elastic_functions(self): + structure = bulk("Al", cubic=True) + df_pot_selected = get_potential_by_name( + potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', + resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), + ) + result_dict = evaluate_with_lammps( + task_dict={"optimize_positions_and_volume": structure}, + potential_dataframe=df_pot_selected, + ) + sym_dict, structure_dict = generate_structures_helper( + structure=result_dict["structure_with_optimized_positions_and_volume"], + eps_range=0.005, + num_of_point=5, + zero_strain_job_name="s_e_0", + sqrt_eta=True, + ) + result_dict = evaluate_with_lammps( + task_dict={"calc_energy": structure_dict}, + potential_dataframe=df_pot_selected, + ) + sym_dict, elastic_dict = analyse_structures_helper( + output_dict=result_dict, + sym_dict=sym_dict, + fit_order=2, + zero_strain_job_name="s_e_0", + ) + self.assertTrue(np.isclose(elastic_dict["elastic_matrix"], np.array( + [ + [114.10311701, 60.51102935, 60.51102935, 0., 0., 0.], + [60.51102935, 114.10311701, 60.51102935, 0., 0., 0.], + [60.51102935, 60.51102935, 114.10311701, 0., 0., 0.], + [0., 0., 0., 51.23853765, 0., 0.], + [0., 0., 0., 0., 51.23853765, 0.], + [0., 0., 0., 0., 0., 51.23853765] + ]) + ).all()) + self.assertEqual(sym_dict['SGN'], 225) + self.assertEqual(sym_dict['LC'], 'CI') + self.assertEqual(sym_dict['Lag_strain_list'], ['01', '08', '23']) + self.assertTrue(np.isclose(sym_dict['epss'], np.array([-0.005, -0.0025, 0., 0.0025, 0.005])).all()) + self.assertTrue(np.isclose(sym_dict["v0"], 66.43035441556098)) + self.assertTrue(np.isclose(sym_dict["e0"], -13.439999952735112)) + self.assertTrue(np.isclose(sym_dict['strain_energy'], np.array( + [ + [ + (-0.005, -13.436320248980278), + (-0.0025, -13.439079680886989), + (0.0, -13.439999952735112), + (0.0024999999999999996, -13.439084974614394), + (0.005, -13.436364320399795) + ], + [ + (-0.005, -13.43817471490433), + (-0.0025, -13.439544638502628), + (0.0, -13.439999952735112), + (0.0024999999999999996, -13.43954822781134), + (0.005, -13.438204192615181) + ], + [ + (-0.005, -13.436741954502294), + (-0.0025, -13.439195465714551), + (0.0, -13.439999952735112), + (0.0024999999999999996, -13.439213491269701), + (0.005, -13.436885713447486) + ] + ])).all() + ) + self.assertTrue(np.isclose(sym_dict["A2"], np.array([2.20130388, 1.08985578, 1.91883479])).all()) + self.assertTrue(np.isclose(elastic_dict['bulkmodul_voigt'], 78.37505857279467)) + self.assertTrue(np.isclose(elastic_dict['shearmodul_voigt'], 41.46154012284969)) + self.assertTrue(np.isclose(elastic_dict['youngsmodul_voigt'], 105.73882997912072)) + self.assertTrue(np.isclose(elastic_dict['poissonsratio_voigt'], 0.2751435386362729)) + self.assertTrue(np.isclose(elastic_dict['elastic_matrix_inverse'], np.array( + [ + [0.01385733, -0.00480214, -0.00480214, 0., 0., 0.], + [-0.00480214, 0.01385733, -0.00480214, 0., 0., 0.], + [-0.00480214, -0.00480214, 0.01385733, 0., 0., 0.], + [0., 0., 0., 0.01951656, 0., 0.], + [0., 0., 0., 0., 0.01951656, 0.], + [0., 0., 0., 0., 0., 0.01951656] + ] + )).all()) + self.assertTrue(np.isclose(elastic_dict['bulkmodul_reuss'], 78.37505857279469)) + self.assertTrue(np.isclose(elastic_dict['shearmodul_reuss'], 37.54104251720356)) + self.assertTrue(np.isclose(elastic_dict['youngsmodul_reuss'], 97.11702764970639)) + self.assertTrue(np.isclose(elastic_dict['poissonsratio_reuss'], 0.29347803281170937)) + self.assertTrue(np.isclose(elastic_dict['bulkmodul_hill'], 78.37505857279467)) + self.assertTrue(np.isclose(elastic_dict['shearmodul_hill'], 39.501291320026624)) + self.assertTrue(np.isclose(elastic_dict['youngsmodul_hill'], 101.45869947879392)) + self.assertTrue(np.isclose(elastic_dict['poissonsratio_hill'], 0.2842453510798992)) + self.assertTrue(np.isclose(elastic_dict['AVR'], 4.962492964955925)) diff --git a/tests/test_evcurve_lammps.py b/tests/test_evcurve_lammps.py index 930c865a..2240f3a2 100644 --- a/tests/test_evcurve_lammps.py +++ b/tests/test_evcurve_lammps.py @@ -5,12 +5,6 @@ import unittest from atomistics.workflows import EnergyVolumeCurveWorkflow, optimize_positions_and_volume -from atomistics.workflows.evcurve.debye import get_thermal_properties -from atomistics.workflows.evcurve.helper import ( - analyse_structures_helper, - generate_structures_helper, -) - try: from atomistics.calculators import ( @@ -63,48 +57,3 @@ def test_calc_evcurve(self): self.assertEqual(len(temperatures_ev), 2) self.assertEqual(len(volumes_ev), 2) self.assertTrue(volumes_ev[0] < volumes_ev[-1]) - - def test_calc_evcurve_functional(self): - structure = bulk("Al", cubic=True) - df_pot_selected = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', - resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), - ) - result_dict = evaluate_with_lammps( - task_dict={"optimize_positions_and_volume": structure}, - potential_dataframe=df_pot_selected, - ) - structure_dict = generate_structures_helper( - structure=result_dict["structure_with_optimized_positions_and_volume"], - vol_range=0.05, - num_points=11, - strain_lst=None, - axes=('x', 'y', 'z'), - ) - result_dict = evaluate_with_lammps( - task_dict={"calc_energy": structure_dict}, - potential_dataframe=df_pot_selected, - ) - fit_dict = analyse_structures_helper( - output_dict=result_dict, - structure_dict=structure_dict, - fit_type="polynomial", - fit_order=3 - ) - thermal_properties_dict = get_thermal_properties( - fit_dict=fit_dict, - masses=structure.get_masses(), - t_min=1.0, - t_max=1500.0, - t_step=50.0, - temperatures=[100, 1000], - constant_volume=False, - output_keys=["temperatures", "volumes"], - ) - temperatures_ev, volumes_ev = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] - self.assertTrue(np.isclose(fit_dict['volume_eq'], 66.43019853103964)) - self.assertTrue(np.isclose(fit_dict['bulkmodul_eq'], 77.7250135953191)) - self.assertTrue(np.isclose(fit_dict['b_prime_eq'], 1.2795467367276832)) - self.assertEqual(len(temperatures_ev), 2) - self.assertEqual(len(volumes_ev), 2) - self.assertTrue(volumes_ev[0] < volumes_ev[-1]) diff --git a/tests/test_evcurve_lammps_functional.py b/tests/test_evcurve_lammps_functional.py new file mode 100644 index 00000000..59852ffc --- /dev/null +++ b/tests/test_evcurve_lammps_functional.py @@ -0,0 +1,71 @@ +import os + +from ase.build import bulk +import numpy as np +import unittest + +from atomistics.workflows.evcurve.debye import get_thermal_properties +from atomistics.workflows.evcurve.helper import ( + analyse_structures_helper, + generate_structures_helper, +) + + +try: + from atomistics.calculators import ( + evaluate_with_lammps, get_potential_by_name + ) + + skip_lammps_test = False +except ImportError: + skip_lammps_test = True + + +@unittest.skipIf( + skip_lammps_test, "LAMMPS is not installed, so the LAMMPS tests are skipped." +) +class TestEvCurve(unittest.TestCase): + def test_calc_evcurve_functional(self): + structure = bulk("Al", cubic=True) + df_pot_selected = get_potential_by_name( + potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', + resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), + ) + result_dict = evaluate_with_lammps( + task_dict={"optimize_positions_and_volume": structure}, + potential_dataframe=df_pot_selected, + ) + structure_dict = generate_structures_helper( + structure=result_dict["structure_with_optimized_positions_and_volume"], + vol_range=0.05, + num_points=11, + strain_lst=None, + axes=('x', 'y', 'z'), + ) + result_dict = evaluate_with_lammps( + task_dict={"calc_energy": structure_dict}, + potential_dataframe=df_pot_selected, + ) + fit_dict = analyse_structures_helper( + output_dict=result_dict, + structure_dict=structure_dict, + fit_type="polynomial", + fit_order=3 + ) + thermal_properties_dict = get_thermal_properties( + fit_dict=fit_dict, + masses=structure.get_masses(), + t_min=1.0, + t_max=1500.0, + t_step=50.0, + temperatures=[100, 1000], + constant_volume=False, + output_keys=["temperatures", "volumes"], + ) + temperatures_ev, volumes_ev = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] + self.assertTrue(np.isclose(fit_dict['volume_eq'], 66.43019853103964)) + self.assertTrue(np.isclose(fit_dict['bulkmodul_eq'], 77.7250135953191)) + self.assertTrue(np.isclose(fit_dict['b_prime_eq'], 1.2795467367276832)) + self.assertEqual(len(temperatures_ev), 2) + self.assertEqual(len(volumes_ev), 2) + self.assertTrue(volumes_ev[0] < volumes_ev[-1]) From f0e0e116a42f7bb9440de16fa18395777da5d54f Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 20:40:31 -0500 Subject: [PATCH 10/11] Test use assertAlmostEqual() --- tests/test_elastic_lammps_functional.py | 30 ++++++++++++------------- tests/test_evcurve_lammps_functional.py | 7 +++--- 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/tests/test_elastic_lammps_functional.py b/tests/test_elastic_lammps_functional.py index 3227df67..3e58e2bb 100644 --- a/tests/test_elastic_lammps_functional.py +++ b/tests/test_elastic_lammps_functional.py @@ -64,8 +64,8 @@ def test_calc_elastic_functions(self): self.assertEqual(sym_dict['LC'], 'CI') self.assertEqual(sym_dict['Lag_strain_list'], ['01', '08', '23']) self.assertTrue(np.isclose(sym_dict['epss'], np.array([-0.005, -0.0025, 0., 0.0025, 0.005])).all()) - self.assertTrue(np.isclose(sym_dict["v0"], 66.43035441556098)) - self.assertTrue(np.isclose(sym_dict["e0"], -13.439999952735112)) + self.assertAlmostEqual(sym_dict["v0"], 66.43035441556098) + self.assertAlmostEqual(sym_dict["e0"], -13.439999952735112) self.assertTrue(np.isclose(sym_dict['strain_energy'], np.array( [ [ @@ -92,10 +92,10 @@ def test_calc_elastic_functions(self): ])).all() ) self.assertTrue(np.isclose(sym_dict["A2"], np.array([2.20130388, 1.08985578, 1.91883479])).all()) - self.assertTrue(np.isclose(elastic_dict['bulkmodul_voigt'], 78.37505857279467)) - self.assertTrue(np.isclose(elastic_dict['shearmodul_voigt'], 41.46154012284969)) - self.assertTrue(np.isclose(elastic_dict['youngsmodul_voigt'], 105.73882997912072)) - self.assertTrue(np.isclose(elastic_dict['poissonsratio_voigt'], 0.2751435386362729)) + self.assertAlmostEqual(elastic_dict['bulkmodul_voigt'], 78.37505857279467) + self.assertAlmostEqual(elastic_dict['shearmodul_voigt'], 41.46154012284969) + self.assertAlmostEqual(elastic_dict['youngsmodul_voigt'], 105.73882997912072) + self.assertAlmostEqual(elastic_dict['poissonsratio_voigt'], 0.2751435386362729) self.assertTrue(np.isclose(elastic_dict['elastic_matrix_inverse'], np.array( [ [0.01385733, -0.00480214, -0.00480214, 0., 0., 0.], @@ -106,12 +106,12 @@ def test_calc_elastic_functions(self): [0., 0., 0., 0., 0., 0.01951656] ] )).all()) - self.assertTrue(np.isclose(elastic_dict['bulkmodul_reuss'], 78.37505857279469)) - self.assertTrue(np.isclose(elastic_dict['shearmodul_reuss'], 37.54104251720356)) - self.assertTrue(np.isclose(elastic_dict['youngsmodul_reuss'], 97.11702764970639)) - self.assertTrue(np.isclose(elastic_dict['poissonsratio_reuss'], 0.29347803281170937)) - self.assertTrue(np.isclose(elastic_dict['bulkmodul_hill'], 78.37505857279467)) - self.assertTrue(np.isclose(elastic_dict['shearmodul_hill'], 39.501291320026624)) - self.assertTrue(np.isclose(elastic_dict['youngsmodul_hill'], 101.45869947879392)) - self.assertTrue(np.isclose(elastic_dict['poissonsratio_hill'], 0.2842453510798992)) - self.assertTrue(np.isclose(elastic_dict['AVR'], 4.962492964955925)) + self.assertAlmostEqual(elastic_dict['bulkmodul_reuss'], 78.37505857279469) + self.assertAlmostEqual(elastic_dict['shearmodul_reuss'], 37.54104251720356) + self.assertAlmostEqual(elastic_dict['youngsmodul_reuss'], 97.11702764970639) + self.assertAlmostEqual(elastic_dict['poissonsratio_reuss'], 0.29347803281170937) + self.assertAlmostEqual(elastic_dict['bulkmodul_hill'], 78.37505857279467) + self.assertAlmostEqual(elastic_dict['shearmodul_hill'], 39.501291320026624) + self.assertAlmostEqual(elastic_dict['youngsmodul_hill'], 101.45869947879392) + self.assertAlmostEqual(elastic_dict['poissonsratio_hill'], 0.2842453510798992) + self.assertAlmostEqual(elastic_dict['AVR'], 4.962492964955925) diff --git a/tests/test_evcurve_lammps_functional.py b/tests/test_evcurve_lammps_functional.py index 59852ffc..f1d1e6b6 100644 --- a/tests/test_evcurve_lammps_functional.py +++ b/tests/test_evcurve_lammps_functional.py @@ -1,7 +1,6 @@ import os from ase.build import bulk -import numpy as np import unittest from atomistics.workflows.evcurve.debye import get_thermal_properties @@ -63,9 +62,9 @@ def test_calc_evcurve_functional(self): output_keys=["temperatures", "volumes"], ) temperatures_ev, volumes_ev = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] - self.assertTrue(np.isclose(fit_dict['volume_eq'], 66.43019853103964)) - self.assertTrue(np.isclose(fit_dict['bulkmodul_eq'], 77.7250135953191)) - self.assertTrue(np.isclose(fit_dict['b_prime_eq'], 1.2795467367276832)) + self.assertAlmostEqual(fit_dict['volume_eq'], 66.43019853103964) + self.assertAlmostEqual(fit_dict['bulkmodul_eq'], 77.7250135953191) + self.assertAlmostEqual(fit_dict['b_prime_eq'], 1.2795467367276832) self.assertEqual(len(temperatures_ev), 2) self.assertEqual(len(volumes_ev), 2) self.assertTrue(volumes_ev[0] < volumes_ev[-1]) From 54a9a6566ffad3e0ac6a96af2447173e90dd125d Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Tue, 23 Apr 2024 20:45:25 -0500 Subject: [PATCH 11/11] fix evcurve --- tests/test_evcurve_lammps_functional.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_evcurve_lammps_functional.py b/tests/test_evcurve_lammps_functional.py index f1d1e6b6..60d1231f 100644 --- a/tests/test_evcurve_lammps_functional.py +++ b/tests/test_evcurve_lammps_functional.py @@ -62,8 +62,8 @@ def test_calc_evcurve_functional(self): output_keys=["temperatures", "volumes"], ) temperatures_ev, volumes_ev = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] - self.assertAlmostEqual(fit_dict['volume_eq'], 66.43019853103964) - self.assertAlmostEqual(fit_dict['bulkmodul_eq'], 77.7250135953191) + self.assertAlmostEqual(fit_dict['volume_eq'], 66.43019790724685) + self.assertAlmostEqual(fit_dict['bulkmodul_eq'], 77.72501703646152) self.assertAlmostEqual(fit_dict['b_prime_eq'], 1.2795467367276832) self.assertEqual(len(temperatures_ev), 2) self.assertEqual(len(volumes_ev), 2)