From cdfc191bff8f2ea274160da62e8c406c9866f787 Mon Sep 17 00:00:00 2001 From: hjung Date: Wed, 5 Jun 2024 12:08:51 +0000 Subject: [PATCH 1/5] espresso.py modified --- wfl/calculators/espresso.py | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/wfl/calculators/espresso.py b/wfl/calculators/espresso.py index 80e8b454..8826a120 100644 --- a/wfl/calculators/espresso.py +++ b/wfl/calculators/espresso.py @@ -17,7 +17,8 @@ from ase.io.espresso import kspacing_to_grid from .wfl_fileio_calculator import WFLFileIOCalculator -from .utils import save_results, parse_genericfileio_profile_argv +from wfl.utils.save_calc_results import save_calc_results +from .utils import parse_genericfileio_profile_argv # NOMAD compatible, see https://nomad-lab.eu/prod/rae/gui/uploads _default_keep_files = ["*.pwo"] @@ -77,20 +78,10 @@ def __init__(self, keep_files="default", rundir_prefix="run_QE_", raise ValueError("Cannot specify both calculator_exec and profile") if " -in " in calculator_exec: raise ValueError("calculator_exec should not include espresso command line arguments such as ' -in PREFIX.pwi'") + if "pseudo_dir" not in kwargs_command: + raise ValueError(f"calculator_exec also requires pseudo_dir to create EspressoProfile") + kwargs_command["profile"] = EspressoProfile(calculator_exec, kwargs_command.pop("pseudo_dir")) - # newer syntax, but pass binary without a keyword (which changed from "argv" to "exc" - # to "binary" over time), assuming it's first argument - argv = shlex.split(calculator_exec) - try: - kwargs_command["profile"] = EspressoProfile(argv=argv) - except TypeError: - binary, parallel_info = parse_genericfileio_profile_argv(argv) - # argument names keep changing (e.g. pseudo_path -> pseudo_dir), just pass first two as positional - # and hope order doesn't change - if "pseudo_dir" not in kwargs_command: - raise ValueError(f"calculator_exec also requires pseudo_dir to create EspressoProfile") - kwargs_command["profile"] = EspressoProfile(binary, kwargs_command.pop("pseudo_dir"), - parallel_info=parallel_info) elif "profile" not in kwargs_command: raise ValueError("EspressoProfile is defined but neither calculator_exec nor profile was specified") From 60e2cc1f711fdec73d1e4de06d483becb47f7fd3 Mon Sep 17 00:00:00 2001 From: hjung Date: Fri, 26 Jul 2024 15:07:28 +0000 Subject: [PATCH 2/5] md modified --- wfl/generate/md/__init__.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/wfl/generate/md/__init__.py b/wfl/generate/md/__init__.py index 2d59b6d9..c3bbc68d 100644 --- a/wfl/generate/md/__init__.py +++ b/wfl/generate/md/__init__.py @@ -19,7 +19,7 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBerendsen", temperature=None, temperature_tau=None, - pressure=None, pressure_tau=None, compressibility_au=None, compressibility_fd_displ=0.01, + pressure=None, pressure_tau=None, compressibility_fd_displ=0.01, logger=None, logfile=None, traj_step_interval=1, skip_failures=True, results_prefix='last_op__md_', verbose=False, update_config_type="append", traj_select_during_func=lambda at: True, traj_select_after_func=None, abort_check=None, rng=None, _autopara_per_item_info=None): @@ -53,10 +53,12 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere pressure_tau: float, default None time scale for Berendsen constant P volume rescaling (fs) ignored if pressure is None, defaults to 3*temperature_tau - compressibility_au: float, default None - compressibility, if available, for NPTBerendsen compressibility_fd_displ: float, default 0.01 finite difference in strain to use when computing compressibility for NPTBerendsen + logger: default None + ase.md.MDLogger - derived class that records user customized properties over MD simulation + logfile: default None + name of logfile used by the logger traj_step_interval: int, default 1 interval between trajectory snapshots skip_failures: bool, default True @@ -95,7 +97,9 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere all_trajs = [] - if verbose: + if logger is not None: + assert logfile is not None, "logfile should be provided." + elif verbose: logfile = '-' else: logfile = None @@ -130,7 +134,10 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere rng = _autopara_per_item_info[at_i].get("rng") at.calc = calculator - if pressure is not None and compressibility_au is None: +# if logger is not None: +# logger = logger(at, logfile, mode="w") + compressibility = None + if pressure is not None: pressure = sample_pressure(pressure, at, rng=rng) at.info['MD_pressure_GPa'] = pressure # convert to ASE internal units @@ -144,15 +151,18 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere Em = at.get_potential_energy() at.set_cell(c0, scale_atoms=True) d2E_dF2 = (Ep + Em - 2.0 * E0) / (compressibility_fd_displ ** 2) - compressibility_au = at.get_volume() / d2E_dF2 + compressibility = at.get_volume() / d2E_dF2 if temperature is not None: # set initial temperature assert rng is not None MaxwellBoltzmannDistribution(at, temperature_K=temperature[0]['T_i'], force_temp=True, communicator=None, rng=rng) Stationary(at, preserve_temperature=True) - - stage_kwargs = {'timestep': dt * fs, 'logfile': logfile} + + if logger is not None: + stage_kwargs = {'timestep': dt * fs, 'logfile': None} + else: + stage_kwargs = {'timestep': dt * fs, "logfile": logfile} if temperature_tau is None: # NVE @@ -171,7 +181,7 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere if pressure is not None: md_constructor = NPTBerendsen stage_kwargs['pressure_au'] = pressure - stage_kwargs['compressibility_au'] = compressibility_au + stage_kwargs['compressibility_au'] = compressibility stage_kwargs['taut'] = temperature_tau * fs stage_kwargs['taup'] = pressure_tau * fs if pressure_tau is not None else temperature_tau * fs * 3 else: @@ -211,6 +221,7 @@ def process_step(interval): if not first_step_of_later_stage and cur_step % interval == 0: at.info['MD_time_fs'] = cur_step * dt at.info['MD_step'] = cur_step + at.info["MD_current_temperature"] = at.get_temperature() at_save = at_copy_save_calc_results(at, prefix=results_prefix) if traj_select_during_func(at): @@ -235,6 +246,8 @@ def process_step(interval): md = md_constructor(at, **stage_kwargs) md.attach(process_step, 1, traj_step_interval) + if logger is not None: + md.attach(logger(md, at, logfile, mode="w")) if stage_i > 0: first_step_of_later_stage = True From 76874358633505b6c374e999c235fce886d9caf3 Mon Sep 17 00:00:00 2001 From: hjung Date: Tue, 30 Jul 2024 14:49:17 +0000 Subject: [PATCH 3/5] MD logger patch --- wfl/generate/md/__init__.py | 52 +++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/wfl/generate/md/__init__.py b/wfl/generate/md/__init__.py index c3bbc68d..124efcc3 100644 --- a/wfl/generate/md/__init__.py +++ b/wfl/generate/md/__init__.py @@ -6,6 +6,7 @@ from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary from ase.md.verlet import VelocityVerlet from ase.md.langevin import Langevin +from ase.md.logger import MDLogger from ase.units import GPa, fs from wfl.autoparallelize import autoparallelize, autoparallelize_docstring @@ -19,10 +20,10 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBerendsen", temperature=None, temperature_tau=None, - pressure=None, pressure_tau=None, compressibility_fd_displ=0.01, logger=None, logfile=None, + pressure=None, pressure_tau=None, compressibility_au=None, compressibility_fd_displ=0.01, traj_step_interval=1, skip_failures=True, results_prefix='last_op__md_', verbose=False, update_config_type="append", - traj_select_during_func=lambda at: True, traj_select_after_func=None, abort_check=None, rng=None, - _autopara_per_item_info=None): + traj_select_during_func=lambda at: True, traj_select_after_func=None, abort_check=None, + logger_kwargs=None, logger_interval=None, rng=None, _autopara_per_item_info=None): """runs an MD trajectory with aggresive, not necessarily physical, integrators for sampling configs. By default calculator properties for each frame stored in keys prefixed with "last_op__md_", which may be overwritten by next operation. @@ -53,12 +54,10 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere pressure_tau: float, default None time scale for Berendsen constant P volume rescaling (fs) ignored if pressure is None, defaults to 3*temperature_tau + compressibility_au: float, default None + compressibility, if available, for NPTBerendsen compressibility_fd_displ: float, default 0.01 finite difference in strain to use when computing compressibility for NPTBerendsen - logger: default None - ase.md.MDLogger - derived class that records user customized properties over MD simulation - logfile: default None - name of logfile used by the logger traj_step_interval: int, default 1 interval between trajectory snapshots skip_failures: bool, default True @@ -83,6 +82,11 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere checks the MD snapshots and aborts the simulation on some condition. rng: numpy.random.Generator, default None random number generator to use (needed for pressure sampling, initial temperature, or Langevin dynamics) + logger_kwargs: dict, default None + kwargs to MDLogger to attach to each MD run, including "logfile" as string to which + config number will be appended + logger_interval: int, default None + interval for logger _autopara_per_item_info: dict INTERNALLY used by autoparallelization framework to make runs reproducible (see wfl.autoparallelize.autoparallelize() docs) @@ -97,13 +101,15 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere all_trajs = [] - if logger is not None: - assert logfile is not None, "logfile should be provided." - elif verbose: + if verbose: logfile = '-' else: logfile = None + if logger_kwargs is not None: + logger_constructor = logger_kwargs.pop("logger", MDLogger) + logger_logfile = logger_kwargs["logfile"] + if temperature_tau is None and (temperature is not None and not isinstance(temperature, (float, int))): raise RuntimeError(f'NVE (temperature_tau is None) can only accept temperature=float for initial T, got {type(temperature)}') @@ -132,12 +138,10 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere # get rng from autopara_per_item info if available ("rng" arg that was passed in was # already used by autoparallelization framework to set "rng" key in per-item dict) rng = _autopara_per_item_info[at_i].get("rng") + item_i = _autopara_per_item_info[at_i].get("item_i") at.calc = calculator -# if logger is not None: -# logger = logger(at, logfile, mode="w") - compressibility = None - if pressure is not None: + if pressure is not None and compressibility_au is None: pressure = sample_pressure(pressure, at, rng=rng) at.info['MD_pressure_GPa'] = pressure # convert to ASE internal units @@ -151,18 +155,15 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere Em = at.get_potential_energy() at.set_cell(c0, scale_atoms=True) d2E_dF2 = (Ep + Em - 2.0 * E0) / (compressibility_fd_displ ** 2) - compressibility = at.get_volume() / d2E_dF2 + compressibility_au = at.get_volume() / d2E_dF2 if temperature is not None: # set initial temperature assert rng is not None MaxwellBoltzmannDistribution(at, temperature_K=temperature[0]['T_i'], force_temp=True, communicator=None, rng=rng) Stationary(at, preserve_temperature=True) - - if logger is not None: - stage_kwargs = {'timestep': dt * fs, 'logfile': None} - else: - stage_kwargs = {'timestep': dt * fs, "logfile": logfile} + + stage_kwargs = {'timestep': dt * fs, 'logfile': logfile} if temperature_tau is None: # NVE @@ -181,7 +182,7 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere if pressure is not None: md_constructor = NPTBerendsen stage_kwargs['pressure_au'] = pressure - stage_kwargs['compressibility_au'] = compressibility + stage_kwargs['compressibility_au'] = compressibility_au stage_kwargs['taut'] = temperature_tau * fs stage_kwargs['taup'] = pressure_tau * fs if pressure_tau is not None else temperature_tau * fs * 3 else: @@ -245,9 +246,14 @@ def process_step(interval): at.info['MD_temperature_K'] = stage_kwargs['temperature_K'] md = md_constructor(at, **stage_kwargs) + md.attach(process_step, 1, traj_step_interval) - if logger is not None: - md.attach(logger(md, at, logfile, mode="w")) + if logger_kwargs is not None: + logger_kwargs["logfile"] = f"{logger_logfile}.item_{item_i}" + logger_kwargs["dyn"] = md + logger_kwargs["atoms"] = at + logger = logger_constructor(**logger_kwargs) + md.attach(logger, logger_interval) if stage_i > 0: first_step_of_later_stage = True From 707072a96b806def25fdccfc25252d5131082ee6 Mon Sep 17 00:00:00 2001 From: hjung Date: Tue, 30 Jul 2024 14:55:18 +0000 Subject: [PATCH 4/5] documentation modified --- wfl/generate/md/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wfl/generate/md/__init__.py b/wfl/generate/md/__init__.py index 124efcc3..91402e69 100644 --- a/wfl/generate/md/__init__.py +++ b/wfl/generate/md/__init__.py @@ -84,7 +84,7 @@ def _sample_autopara_wrappable(atoms, calculator, steps, dt, integrator="NVTBere random number generator to use (needed for pressure sampling, initial temperature, or Langevin dynamics) logger_kwargs: dict, default None kwargs to MDLogger to attach to each MD run, including "logfile" as string to which - config number will be appended + config number will be appended. User defined ase.md.MDLogger derived class can be provided with "logger" as key. logger_interval: int, default None interval for logger _autopara_per_item_info: dict From 5bc586cf8dccdc958fe88264102700c697db3508 Mon Sep 17 00:00:00 2001 From: hjung Date: Tue, 30 Jul 2024 15:56:41 +0000 Subject: [PATCH 5/5] pytest for MD logger added --- tests/test_md.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tests/test_md.py b/tests/test_md.py index 340bdc1d..fcbb99ed 100644 --- a/tests/test_md.py +++ b/tests/test_md.py @@ -2,12 +2,15 @@ import pytest import numpy as np +import os +from pathlib import Path from ase import Atoms import ase.io from ase.build import bulk from ase.calculators.emt import EMT from ase.units import fs +from ase.md.logger import MDLogger from wfl.autoparallelize import autoparainfo from wfl.generate import md @@ -216,3 +219,30 @@ def test_md_abort_function(cu_slab): rng=np.random.default_rng(1)) assert len(list(atoms_traj)) < 501 + + +def test_md_attach_logger(cu_slab): + + calc = EMT() + autopara_info = autoparainfo.AutoparaInfo(num_python_subprocesses=2, num_inputs_per_python_subprocess=1, skip_failed=False) + + inputs = ConfigSet([cu_slab, cu_slab]) + outputs = OutputSpec() + + logger_kwargs = { + "logger" : MDLogger, + "logfile" : "test_log", + } + + atoms_traj = md.md(inputs, outputs, calculator=calc, integrator="Langevin", steps=300, dt=1.0, + temperature=500.0, temperature_tau=100/fs, logger_kwargs=logger_kwargs, logger_interval=1, + rng=np.random.default_rng(1), autopara_info=autopara_info,) + + atoms_traj = list(atoms_traj) + atoms_final = atoms_traj[-1] + + workdir = Path(os.getcwd()) + + assert len(atoms_traj) == 602 + assert all([Path(workdir / "test_log.item_0").is_file(), Path(workdir / "test_log.item_1").is_file()]) +