diff --git a/.ci_support/environment-lammps.yml b/.ci_support/environment-lammps.yml index 41e98976..2f29837a 100644 --- a/.ci_support/environment-lammps.yml +++ b/.ci_support/environment-lammps.yml @@ -8,3 +8,4 @@ dependencies: - iprpy-data =2023.07.25 - dynaphopy =1.18.0 - lammpsparser =0.0.1 +- pyscal3 = 3.3.1 diff --git a/src/atomistics/calculators/lammps/commands.py b/src/atomistics/calculators/lammps/commands.py index 4f16146e..bf2b0eb3 100644 --- a/src/atomistics/calculators/lammps/commands.py +++ b/src/atomistics/calculators/lammps/commands.py @@ -20,7 +20,7 @@ LAMMPS_TIMESTEP = "timestep {{timestep}}" -LAMMPS_VELOCITY = "velocity all create $(2 * {{ temp }}) {{seed}} dist {{dist}}" +LAMMPS_VELOCITY = "velocity all create $({{velocity_rescale_factor}} * {{ temp }}) {{seed}} dist {{dist}}" LAMMPS_ENSEMBLE_NPT = "fix ensemble all npt temp {{Tstart}} {{Tstop}} {{Tdamp}} iso {{Pstart}} {{Pstop}} {{Pdamp}}" diff --git a/src/atomistics/calculators/lammps/helpers.py b/src/atomistics/calculators/lammps/helpers.py index 5472d43d..9a44568c 100644 --- a/src/atomistics/calculators/lammps/helpers.py +++ b/src/atomistics/calculators/lammps/helpers.py @@ -102,6 +102,7 @@ def lammps_thermal_expansion_loop( Pdamp: float = 1.0, seed: int = 4928459, dist: str = "gaussian", + velocity_rescale_factor: float = 2.0, lmp=None, output_keys=OutputThermalExpansion.keys(), **kwargs, @@ -115,6 +116,7 @@ def lammps_thermal_expansion_loop( timestep=timestep, seed=seed, dist=dist, + velocity_rescale_factor=velocity_rescale_factor, ), lmp=lmp, **kwargs, diff --git a/src/atomistics/calculators/lammps/libcalculator.py b/src/atomistics/calculators/lammps/libcalculator.py index 42b2e81f..f28094fe 100644 --- a/src/atomistics/calculators/lammps/libcalculator.py +++ b/src/atomistics/calculators/lammps/libcalculator.py @@ -1,6 +1,6 @@ from __future__ import annotations -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Optional import numpy as np import pandas @@ -49,11 +49,12 @@ def optimize_positions_and_volume_with_lammpslib( maxiter: int = 100000, maxeval: int = 10000000, thermo: int = 10, + vmax: Optional[float] = None, lmp=None, **kwargs, ) -> Atoms: template_str = ( - LAMMPS_MINIMIZE_VOLUME + _get_vmax_command(vmax=vmax) + "\n" + LAMMPS_THERMO_STYLE + "\n" @@ -154,12 +155,12 @@ def calc_molecular_dynamics_nvt_with_lammpslib( timestep: float = 0.001, seed: int = 4928459, dist: str = "gaussian", - disable_initial_velocity: bool = False, + velocity_rescale_factor: Optional[float] = 2.0, lmp=None, output_keys=OutputMolecularDynamics.keys(), **kwargs, ) -> dict: - if not disable_initial_velocity: + if velocity_rescale_factor is not None: init_str = ( LAMMPS_THERMO_STYLE + "\n" @@ -180,6 +181,7 @@ def calc_molecular_dynamics_nvt_with_lammpslib( timestep=timestep, seed=seed, dist=dist, + velocity_rescale_factor=velocity_rescale_factor, ) else: init_str = ( @@ -232,12 +234,17 @@ def calc_molecular_dynamics_npt_with_lammpslib( Pdamp: float = 1.0, seed: int = 4928459, dist: str = "gaussian", - disable_initial_velocity: bool = False, + velocity_rescale_factor: Optional[float] = 2.0, + couple_xyz: bool = False, lmp=None, output_keys=OutputMolecularDynamics.keys(), **kwargs, ) -> dict: - if not disable_initial_velocity: + if couple_xyz: + LAMMPS_ENSEMBLE_NPT_XYZ = LAMMPS_ENSEMBLE_NPT + " couple xyz" + else: + LAMMPS_ENSEMBLE_NPT_XYZ = LAMMPS_ENSEMBLE_NPT + if velocity_rescale_factor is not None: init_str = ( LAMMPS_THERMO_STYLE + "\n" @@ -247,7 +254,7 @@ def calc_molecular_dynamics_npt_with_lammpslib( + "\n" + LAMMPS_VELOCITY + "\n" - + LAMMPS_ENSEMBLE_NPT + + LAMMPS_ENSEMBLE_NPT_XYZ ) input_template = Template(init_str).render( thermo=thermo, @@ -261,6 +268,7 @@ def calc_molecular_dynamics_npt_with_lammpslib( timestep=timestep, seed=seed, dist=dist, + velocity_rescale_factor=velocity_rescale_factor, ) else: init_str = ( @@ -270,7 +278,7 @@ def calc_molecular_dynamics_npt_with_lammpslib( + "\n" + LAMMPS_THERMO + "\n" - + LAMMPS_ENSEMBLE_NPT + + LAMMPS_ENSEMBLE_NPT_XYZ ) input_template = Template(init_str).render( thermo=thermo, @@ -314,12 +322,12 @@ def calc_molecular_dynamics_nph_with_lammpslib( Pdamp: float = 1.0, seed: int = 4928459, dist: str = "gaussian", - disable_initial_velocity: bool = False, + velocity_rescale_factor: Optional[float] = 2.0, lmp=None, output_keys=OutputMolecularDynamics.keys(), **kwargs, ) -> dict: - if not disable_initial_velocity: + if velocity_rescale_factor is not None: init_str = ( LAMMPS_THERMO_STYLE + "\n" @@ -340,6 +348,7 @@ def calc_molecular_dynamics_nph_with_lammpslib( timestep=timestep, seed=seed, dist=dist, + velocity_rescale_factor=velocity_rescale_factor, ) else: init_str = ( @@ -389,12 +398,12 @@ def calc_molecular_dynamics_langevin_with_lammpslib( Tdamp: float = 0.1, seed: int = 4928459, dist: str = "gaussian", - disable_initial_velocity: bool = False, + velocity_rescale_factor: Optional[float] = 2.0, lmp=None, output_keys=OutputMolecularDynamics.keys(), **kwargs, ): - if not disable_initial_velocity: + if velocity_rescale_factor is not None: init_str = ( LAMMPS_THERMO_STYLE + "\n" @@ -417,6 +426,7 @@ def calc_molecular_dynamics_langevin_with_lammpslib( timestep=timestep, seed=seed, dist=dist, + velocity_rescale_factor=velocity_rescale_factor, ) else: init_str = ( @@ -473,6 +483,7 @@ def calc_molecular_dynamics_thermal_expansion_with_lammpslib( Pdamp: float = 1.0, seed: int = 4928459, dist: str = "gaussian", + couple_xyz: bool = False, lmp=None, output_keys=OutputThermalExpansion.keys(), **kwargs, @@ -487,7 +498,11 @@ def calc_molecular_dynamics_thermal_expansion_with_lammpslib( + LAMMPS_VELOCITY + "\n" ) - run_str = LAMMPS_ENSEMBLE_NPT + "\n" + LAMMPS_RUN + if couple_xyz: + LAMMPS_ENSEMBLE_NPT_XYZ = LAMMPS_ENSEMBLE_NPT + " couple xyz" + else: + LAMMPS_ENSEMBLE_NPT_XYZ = LAMMPS_ENSEMBLE_NPT + run_str = LAMMPS_ENSEMBLE_NPT_XYZ + "\n" + LAMMPS_RUN temperature_lst = np.arange(Tstart, Tstop + Tstep, Tstep).tolist() return lammps_thermal_expansion_loop( structure=structure, @@ -593,3 +608,13 @@ def evaluate_with_lammpslib( ) lmp.close() return results_dict + + +def _get_vmax_command(vmax: Optional[float]) -> str: + if vmax is not None: + if isinstance(vmax, float): + return LAMMPS_MINIMIZE_VOLUME + " vmax {vmax}".format(vmax=vmax) + else: + raise TypeError("vmax must be a float.") + else: + return LAMMPS_MINIMIZE_VOLUME diff --git a/src/atomistics/calculators/lammps/melting.py b/src/atomistics/calculators/lammps/melting.py new file mode 100644 index 00000000..4824c928 --- /dev/null +++ b/src/atomistics/calculators/lammps/melting.py @@ -0,0 +1,294 @@ +import numpy as np +import operator +import random +from ase.build import bulk +from ase.data import reference_states, atomic_numbers +from structuretoolkit.analyse import ( + get_adaptive_cna_descriptors, + get_diamond_structure_descriptors, +) +from atomistics.shared.output import OutputMolecularDynamics +from atomistics.calculators.lammps import ( + optimize_positions_and_volume_with_lammpslib, + calc_molecular_dynamics_npt_with_lammpslib, +) + + +def _check_diamond(structure): + """ + Utility function to check if the structure is fcc, bcc, hcp or diamond + + Args: + structure (pyiron_atomistics.structure.atoms.Atoms): Atomistic Structure object to check + + Returns: + bool: true if diamond else false + """ + cna_dict = get_adaptive_cna_descriptors( + structure=structure, mode="total", ovito_compatibility=True + ) + dia_dict = get_diamond_structure_descriptors( + structure=structure, mode="total", ovito_compatibility=True + ) + return ( + cna_dict["CommonNeighborAnalysis.counts.OTHER"] + > dia_dict["IdentifyDiamond.counts.OTHER"] + ) + + +def _analyse_structure(structure, mode="total", diamond=False): + """ + Use either common neighbor analysis or the diamond structure detector + + Args: + structure (pyiron_atomistics.structure.atoms.Atoms): The structure to analyze. + mode ("total"/"numeric"/"str"): Controls the style and level + of detail of the output. + - total : return number of atoms belonging to each structure + - numeric : return a per atom list of numbers- 0 for unknown, + 1 fcc, 2 hcp, 3 bcc and 4 icosa + - str : return a per atom string of sructures + diamond (bool): Flag to either use the diamond structure detector or + the common neighbor analysis. + + Returns: + (depends on `mode`) + """ + if not diamond: + return get_adaptive_cna_descriptors( + structure=structure, mode=mode, ovito_compatibility=True + ) + else: + return get_diamond_structure_descriptors( + structure=structure, mode=mode, ovito_compatibility=True + ) + + +def _analyse_minimized_structure(structure): + """ + + Args: + ham (GenericJob): + + Returns: + + """ + diamond_flag = _check_diamond(structure=structure) + final_structure_dict = _analyse_structure( + structure=structure, mode="total", diamond=diamond_flag + ) + key_max = max(final_structure_dict.items(), key=operator.itemgetter(1))[0] + number_of_atoms = len(structure) + distribution_initial = final_structure_dict[key_max] / number_of_atoms + distribution_initial_half = distribution_initial / 2 + return ( + structure, + key_max, + number_of_atoms, + distribution_initial_half, + final_structure_dict, + ) + + +def _next_calc(structure, potential, temperature, seed, run_time_steps=10000): + """ + Calculate NPT ensemble at a given temperature using the job defined in the project parameters: + - job_type: Type of Simulation code to be used + - project: Project object used to create the job + - potential: Interatomic Potential + - queue (optional): HPC Job queue to be used + + Args: + structure (pyiron_atomistics.structure.atoms.Atoms): Atomistic Structure object to be set to the job as input sturcture + temperature (float): Temperature of the Molecular dynamics calculation + run_time_steps (int): Number of Molecular dynamics steps + + Returns: + Final Atomistic Structure object + """ + output_md_dict = calc_molecular_dynamics_npt_with_lammpslib( + structure=structure, + potential_dataframe=potential, + Tstart=temperature, + Tstop=temperature, + Tdamp=0.1, + run=run_time_steps, + thermo=1000, + timestep=0.001, + Pstart=0.0, + Pstop=0.0, + Pdamp=1.0, + seed=seed, + dist="gaussian", + velocity_rescale_factor=1.0, + lmp=None, + output_keys=OutputMolecularDynamics.keys(), + ) + structure_md = structure.copy() + structure_md.set_positions(output_md_dict["positions"][-1]) + structure_md.set_cell(output_md_dict["cell"][-1]) + return structure_md + + +def _next_step_funct( + number_of_atoms, + key_max, + structure_left, + structure_right, + potential, + temperature_left, + temperature_right, + distribution_initial_half, + structure_after_minimization, + run_time_steps, + crystalstructure, + seed, +): + """ + + Args: + number_of_atoms: + key_max: + structure_left: + structure_right: + temperature_left: + temperature_right: + distribution_initial_half: + structure_after_minimization: + run_time_steps: + + Returns: + + """ + structure_left_dict = _analyse_structure( + structure=structure_left, + mode="total", + diamond=crystalstructure.lower() == "diamond", + ) + structure_right_dict = _analyse_structure( + structure=structure_right, + mode="total", + diamond=crystalstructure.lower() == "diamond", + ) + temperature_diff = temperature_right - temperature_left + if ( + structure_left_dict[key_max] / number_of_atoms > distribution_initial_half + and structure_right_dict[key_max] / number_of_atoms > distribution_initial_half + ): + structure_left = structure_right.copy() + temperature_left = temperature_right + temperature_right += temperature_diff + structure_right = _next_calc( + structure=structure_after_minimization, + temperature=temperature_right, + potential=potential, + seed=seed, + run_time_steps=run_time_steps, + ) + elif ( + structure_left_dict[key_max] / number_of_atoms + > distribution_initial_half + > structure_right_dict[key_max] / number_of_atoms + ): + temperature_diff /= 2 + temperature_left += temperature_diff + structure_left = _next_calc( + structure=structure_after_minimization, + temperature=temperature_left, + potential=potential, + seed=seed, + run_time_steps=run_time_steps, + ) + elif ( + structure_left_dict[key_max] / number_of_atoms < distribution_initial_half + and structure_right_dict[key_max] / number_of_atoms < distribution_initial_half + ): + temperature_diff /= 2 + temperature_right = temperature_left + temperature_left -= temperature_diff + structure_right = structure_left.copy() + structure_left = _next_calc( + structure=structure_after_minimization, + temperature=temperature_left, + potential=potential, + seed=seed, + run_time_steps=run_time_steps, + ) + else: + raise ValueError("We should never reach this point!") + return structure_left, structure_right, temperature_left, temperature_right + + +def estimate_melting_temperature( + element, + potential, + strain_run_time_steps=1000, + temperature_left=0, + temperature_right=1000, + number_of_atoms=8000, + seed=None, +): + if seed is None: + seed = random.randint(0, 99999) + crystalstructure = reference_states[atomic_numbers[element]]["symmetry"] + if crystalstructure == "hcp": + basis = bulk(name=element, orthorhombic=True) + else: + basis = bulk(name=element, cubic=True) + basis_lst = [basis.repeat([i, i, i]) for i in range(5, 30)] + basis = basis_lst[ + np.argmin([np.abs(len(b) - number_of_atoms / 2) for b in basis_lst]) + ] + + structure_opt = optimize_positions_and_volume_with_lammpslib( + structure=basis, + potential_dataframe=potential, + min_style="cg", + etol=0.0, + ftol=0.0001, + maxiter=100000, + maxeval=10000000, + thermo=10, + lmp=None, + ) + + ( + structure_after_minimization, + key_max, + number_of_atoms, + distribution_initial_half, + _, + ) = _analyse_minimized_structure(structure=structure_opt) + + structure_left = structure_after_minimization + structure_right = _next_calc( + structure=structure_after_minimization, + temperature=temperature_right, + seed=seed, + potential=potential, + run_time_steps=strain_run_time_steps, + ) + temperature_step = temperature_right - temperature_left + + while temperature_step > 10: + ( + structure_left, + structure_right, + temperature_left, + temperature_right, + ) = _next_step_funct( + number_of_atoms=number_of_atoms, + key_max=key_max, + structure_left=structure_left, + structure_right=structure_right, + potential=potential, + temperature_left=temperature_left, + temperature_right=temperature_right, + distribution_initial_half=distribution_initial_half, + structure_after_minimization=structure_after_minimization, + run_time_steps=strain_run_time_steps, + seed=seed, + crystalstructure=crystalstructure, + ) + temperature_step = temperature_right - temperature_left + return int(round(temperature_left)) diff --git a/tests/test_lammpslib_md.py b/tests/test_lammpslib_md.py index 94041dd1..0fabe2c4 100644 --- a/tests/test_lammpslib_md.py +++ b/tests/test_lammpslib_md.py @@ -71,7 +71,7 @@ def test_lammps_md_nvt_all_no_velocity(self): seed=4928459, dist="gaussian", lmp=None, - disable_initial_velocity=True, + velocity_rescale_factor=None, ) self.assertEqual(result_dict["positions"].shape, (10, 32, 3)) self.assertEqual(result_dict["velocities"].shape, (10, 32, 3)) @@ -163,7 +163,7 @@ def test_lammps_md_npt_all_no_velocity(self): seed=4928459, dist="gaussian", lmp=None, - disable_initial_velocity=True, + velocity_rescale_factor=None, ) self.assertEqual(result_dict["positions"].shape, (10, 32, 3)) self.assertEqual(result_dict["velocities"].shape, (10, 32, 3)) diff --git a/tests/test_melting_lammps.py b/tests/test_melting_lammps.py new file mode 100644 index 00000000..eb113982 --- /dev/null +++ b/tests/test_melting_lammps.py @@ -0,0 +1,33 @@ +import os + +import unittest + + +try: + from atomistics.calculators import get_potential_by_name + from atomistics.calculators.lammps.melting import estimate_melting_temperature + + 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 TestLammpsMelting(unittest.TestCase): + def test_estimate_melting_temperature(self): + potential = get_potential_by_name( + potential_name="1999--Mishin-Y--Al--LAMMPS--ipr1", + resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), + ) + melting_temp = estimate_melting_temperature( + element="Al", + potential=potential, + strain_run_time_steps=1000, + temperature_left=0, + temperature_right=1000, + number_of_atoms=8000, + seed=None, + ) + self.assertIn(melting_temp, [1008, 1023])