From 4a2c3dffbd6644159ca6b60e7c6ef42de151ace4 Mon Sep 17 00:00:00 2001 From: Andrey Sobolev Date: Wed, 3 Jul 2024 20:34:05 +0200 Subject: [PATCH] Adding MD input set to FHI-aims (#3896) * Added MD input set and generator for aims * Added MD input set and generator for aims * Conditionally add velocities to site properties * pre-commit auto-fixes --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- src/pymatgen/io/aims/inputs.py | 19 +++- src/pymatgen/io/aims/sets/core.py | 101 +++++++++++++++++- .../md-si/control.in.gz | Bin 0 -> 1216 bytes .../md-si/geometry.in.gz | Bin 0 -> 106 bytes .../md-si/parameters.json | 1 + tests/io/aims/input_files/geometry.in.si.gz | Bin 222 -> 237 bytes .../io/aims/input_files/geometry.in.si.ref.gz | Bin 243 -> 275 bytes tests/io/aims/input_files/si_ref.json.gz | Bin 570 -> 581 bytes tests/io/aims/test_sets/test_md_generator.py | 48 +++++++++ 9 files changed, 165 insertions(+), 4 deletions(-) create mode 100644 tests/io/aims/aims_input_generator_ref/md-si/control.in.gz create mode 100644 tests/io/aims/aims_input_generator_ref/md-si/geometry.in.gz create mode 100644 tests/io/aims/aims_input_generator_ref/md-si/parameters.json create mode 100644 tests/io/aims/test_sets/test_md_generator.py diff --git a/src/pymatgen/io/aims/inputs.py b/src/pymatgen/io/aims/inputs.py index 342c01b40d0..e03d0355f75 100644 --- a/src/pymatgen/io/aims/inputs.py +++ b/src/pymatgen/io/aims/inputs.py @@ -60,7 +60,7 @@ def from_str(cls, contents: str) -> Self: ] species, coords, is_frac, lattice_vectors = [], [], [], [] - charges_dct, moments_dct = {}, {} + charges_dct, moments_dct, velocities_dct = {}, {}, {} for line in content_lines: inp = line.split() @@ -74,6 +74,8 @@ def from_str(cls, contents: str) -> Self: moments_dct[len(coords) - 1] = float(inp[1]) if inp[0] == "initial_charge": charges_dct[len(coords) - 1] = float(inp[1]) + if inp[0] == "velocity": + velocities_dct[len(coords) - 1] = [float(x) for x in inp[1:]] charge = np.zeros(len(coords)) for key, val in charges_dct.items(): @@ -83,6 +85,10 @@ def from_str(cls, contents: str) -> Self: for key, val in moments_dct.items(): magmom[key] = val + velocity: list[None | list[float]] = [None for _ in coords] + for key, v in velocities_dct.items(): + velocity[key] = v + if len(lattice_vectors) == 3: lattice = Lattice(lattice_vectors) for cc in range(len(coords)): @@ -96,6 +102,9 @@ def from_str(cls, contents: str) -> Self: raise ValueError("Incorrect number of lattice vectors passed.") site_props = {"magmom": magmom, "charge": charge} + if velocities_dct: + site_props["velocity"] = velocity + if lattice is None: structure = Molecule(species, coords, np.sum(charge), site_properties=site_props) else: @@ -137,13 +146,17 @@ def from_structure(cls, structure: Structure | Molecule) -> Self: charges = structure.site_properties.get("charge", np.zeros(len(structure.species))) magmoms = structure.site_properties.get("magmom", np.zeros(len(structure.species))) - for species, coord, charge, magmom in zip(structure.species, structure.cart_coords, charges, magmoms): + velocities = structure.site_properties.get("velocity", [None for _ in structure.species]) + for species, coord, charge, magmom, v in zip( + structure.species, structure.cart_coords, charges, magmoms, velocities + ): content_lines.append(f"atom {coord[0]: .12e} {coord[1]: .12e} {coord[2]: .12e} {species}") if charge != 0: content_lines.append(f" initial_charge {charge:.12e}") if magmom != 0: content_lines.append(f" initial_moment {magmom:.12e}") - + if v is not None and any(v_i != 0.0 for v_i in v): + content_lines.append(f" velocity {' '.join([f'{v_i:.12e}' for v_i in v])}") return cls(_content="\n".join(content_lines), _structure=structure) @property diff --git a/src/pymatgen/io/aims/sets/core.py b/src/pymatgen/io/aims/sets/core.py index db689c60701..c858b2f1236 100644 --- a/src/pymatgen/io/aims/sets/core.py +++ b/src/pymatgen/io/aims/sets/core.py @@ -2,7 +2,7 @@ from __future__ import annotations -from dataclasses import dataclass +from dataclasses import dataclass, field from typing import TYPE_CHECKING from pymatgen.core import Structure @@ -14,6 +14,13 @@ from pymatgen.core import Molecule +_valid_dynamics: dict[str, tuple[str, ...]] = { + "nve": ("", "4th_order", "damped"), + "nvt": ("andersen", "berendsen", "parrinello", "nose-hoover"), + "gle": ("thermostat",), +} + + @dataclass class StaticSetGenerator(AimsInputGenerator): """Common class for ground-state generators. @@ -97,3 +104,95 @@ def get_parameter_updates(self, structure: Structure | Molecule, prev_parameters dict: The updated for the parameters for the output section of FHI-aims """ return {"use_pimd_wrapper": (self.host, self.port)} + + +@dataclass +class MDSetGenerator(AimsInputGenerator): + """ + A class for generating FHI-aims input sets for molecular dynamics calculations. + + Parameters + ---------- + ensemble + Molecular dynamics ensemble to run. Options include `nvt`, `nve`, and `gle`. + Default: `nve` + ensemble_specs + A dictionary containing the specifications of the molecular dynamics ensemble. + Valid keys are `type` (the ensemble type, valid types are defined in + `_valid_dynamics` dict), and `parameter` - the control parameter for the thermostat + (not used for `nve` and `nve_4th_order`). + temp + Thermostat temperature. Default: None + time + Simulation time (in picoseconds). Negative value stands for indefinite run. + Default: 5 ps + time_step + The time step (in picoseconds) for the simulation. default: 1 fs + **kwargs + Other keyword arguments that will be passed to :obj:`AimsInputGenerator`. + """ + + calc_type: str = "md" + ensemble: str = "nve" + ensemble_specs: dict[str, Any] = field(default_factory=dict) + temp: float | None = None + time: float = 5.0 + time_step: float = 0.001 + init_velocities: bool = True + + def get_parameter_updates(self, structure: Structure | Molecule, prev_parameters: dict[str, Any]) -> dict: + """Get the parameter updates for the calculation. + + Parameters + ---------- + structure (Structure or Molecule): + The structure to calculate the bands for + prev_parameters (Dict[str, Any]): + The previous parameters + + Returns + ------- + dict + A dictionary of updates to apply. + """ + updates: dict[str, Any] = {"MD_run": [self.time], "MD_time_step": self.time_step} + + # check for ensemble type validity + default_ensemble_types = {"nve": "", "nvt": "parrinello", "gle": "thermostat"} + if self.ensemble not in _valid_dynamics: + raise ValueError(f"Ensemble {self.ensemble} not valid") + ensemble_type = self.ensemble_specs.get("type", default_ensemble_types[self.ensemble]) + if ensemble_type not in _valid_dynamics[self.ensemble]: + raise ValueError( + f"Type {ensemble_type} is not valid for {self.ensemble} ensemble. " + f"Valid types are: {' ,'.join(_valid_dynamics[self.ensemble])}" + ) + ensemble_name = f"{self.ensemble.upper()}_{ensemble_type}" if ensemble_type else self.ensemble.upper() + updates["MD_run"].append(ensemble_name) + + # add temperature + if self.ensemble == "nve": + if not self.init_velocities and "velocity" not in structure.site_properties: + raise ValueError("Velocities must be initialized for NVE ensemble") + else: + if self.temp is None: + raise ValueError(f"Temperature must be set for {ensemble_name} ensemble") + updates["MD_run"].append(self.temp) + + # check for ensemble control parameter + ensemble_parameter = self.ensemble_specs.get("parameter", None) + if ensemble_name not in ("NVE", "NVE_4th_order"): + if ensemble_parameter is None: + raise ValueError(f"Ensemble {ensemble_name} parameter is not defined") + updates["MD_run"].append(ensemble_parameter) + + # ...and put everything in the string + updates["MD_run"] = " ".join(map(str, updates["MD_run"])) + + # initialize velocities + if self.init_velocities: + if self.temp is None: + raise ValueError("Temperature must be set for velocity initialisation") + updates["MD_MB_init"] = self.temp + + return updates diff --git a/tests/io/aims/aims_input_generator_ref/md-si/control.in.gz b/tests/io/aims/aims_input_generator_ref/md-si/control.in.gz new file mode 100644 index 0000000000000000000000000000000000000000..d50a761f109b7bd5d82724f018d06b1cc2f65f14 GIT binary patch literal 1216 zcmV;x1V8&9iwFpbjd^AO17mM)baHQOE@^H6wN}e++d2^4`zr>{rU7hLk|jTidlyNA zpk1^mkRn@wmMEJIB?=^!HvRg0NZF2G>D{>ULe`v_;hEuZhFAYy_8I&t=>MK}DF#L; z?&PzyPQGQpovvk(|JH_+Ex6)l@Z~NyUGqohal|6}`uII>xiPX4N@UaqYd`3bLIM1h~}QQi*Y~7R3|s=7o%~!L438 zjx8?fr*vy_Yhhs|oRiI}U$fQLjS)s6tc7so9SgQ%qn!l}-4R#bl!RJaA<8AU(t3}U z^&uuv61|TW@1r^0lqD9fsOp-+N!@s&+qik|px2tu{i{7Kj&r(0GyLPh6K zoo+1Ef+M|rsa~koNVE9WipI+4e$_seLWx>5?u6x6@hl8f6tX(n#S#p^jp$)Sp!&{p zMM^9nICt;^NL~|fTY-~&=n;xPI#g8kf~65*vDUH>A?cT{b0WVd5_~t$gWQl=DZS)s z2csZ-1Mt}zc(~~bB|3Y+J5B>sb9*TF2Y=3@d6o^Xf`-c1+&Ck|o<_RCiq7S&KI}9T zmhG+LKVlXI?NF0>DcT%%x~hVL7;l{!UP`WVXuU3#CX{WgH*k$gm#vbhQdkFsH!Bo} z2_0=4Oh`o~MLFodOW&}4z4<k z$N`$p;|Kvw(z9rmC2*WAqO<5UnWsdv*?Dvp&1OW4**VV%OQY$W;^dVTs$Q6F`}bBS z;t;~6t2cC^yIp1a)-bH`QFAryOXF{xnLg%cu7qzho37|MTG?P9=qndx6H;~-YI~@6 z7y>AYuH{YfGT2CRNDe&3>@3vw(B)%29A&&LL#ELf`nb3K-ZHxP6UOMW>s$POe+`oc zn!5HK)c4bxd`s^q18+1W2BQ!IsK$&US-09a4%MSoh#TR3DH-Tg(pY!BnMNOhJA>6d z04^o&HpYJN#hlR%Ml)_g;gxJYd+{W^|5BX{?>$3uXe~n;`?!pL?#IXn8fUW@tpUU5 z(YO5ZWX7V2|LRZA`lU3utWilROZG(xb{l7ZUzB3>PPbY}6mF3dlnX~s3^X1DE*0fl zJmM*v?a|n>xH_V-IFMwtq6JG2WGR>(N@B2x*fa|)Fwf1KBqnUW*l~bKw8vo^QMrp~ zkC??)PLF6*G(V!|jVC9vmS$|;Zfa?p&zrj0n`i86{f5)J_oF z;HZJV{oL%d)`QU>qr+$Yzl4gH$0Q0kyuPgrxe<;uU;`%%4mn;Vc{x)0-$gWxpKu<@ zpMbv*S`LmZbmY#1VG17hexew>L(i~iyZtix-V^_@43z@4ZED?A6=`2(40vzi literal 222 zcmV<403rV$iwFoW(_Cc$17~G#ZDn+Fc`j*gE^}!Bosh8#0x=AR_j`)KEjoDA);smU zK}69(9PR;*QlkcJ3%w}#_EtB&x|kuum;e70k}v+NEQif@H_~J~IMEerIHR}7;Zcw) zI*RjEK+SlBTSx8;d_h>Uw`8Cpd%J2T?_llj6uh(O7zWE`J&eA6{tV*{=4Eew-299Y=? z3hd+CfyzP3Qs$-3Kne969UgYU&M;Gr)dg$rJ%C8WpTZ|009106dWYo;wo%PTj;pp? z`qTfaV%w(Hs;xdU9Nu&{_fGan5zIx7m`qpM@PUp&-ine*rVF7K+Q3e_&}){-OzlNo ZuhiwFozUqNL6|7T@yZDn+Fc`j*gE^}!va%E-!t&qVA!Y~Ym?|X{Cj^diN zRNcU<=s`UA08(b0g>98C3ckG^DCnl>!QApU-xr#A{-=s!SS(j#W9t@5ay7BP!ogOQ zgE)p(@agOeswz=h7N=`%xTHofwcHBpLddhVT0;Rz9Z|U}Hm1>8a{D$0mwcN;gPgER0v# zuJ*qdTxa7CSgn*wl{yDYo_Uzp%#7i**?Pg*|fdu zvb&MwigjH8zxd!+VYG1e_^o-1bJIA{xBw{ToaEy9Ah=WNeiZ$sfvgDcq$ zZmh`P`MU{y>BvJY%wL$@hg^O`6ZcrYneY(sdRD%A|M!_c*xD0pmUn7HchjM26?t02 zamise;+e004f&WPO=d>Pk6A0uJPZA--SKGYeP_j5BkjE^!n{P0lM*0WQCiffK+7Fu zyOAj%!X*L(00A*2Ko z3}Ct7_Dt&lVwzJ@69g&5hyYYA!Glq~4g+3-92QkU38^o{dGaf?*0jPpDZ)>5nnpV5 z2IssVr&YzXt`%HH(DahaQlW$%56sKK5CskAy{BN5INvha8B0u8u;h<(EpU_S>+#n= zP@Q)*n#!Mbdu?%9wiwFo|UO{C5|8r?ya%E;NYIARH0Ns^MkDD+MhVS_m5$7%P2OoA%+f&u^ z9=p}baosEu#>kjx(2vPyPo)wI3ol82GxHS0PDezL*OLP=rt`Ac()bJZBp7zfnMn2L*+li)__JxGCP zdz6+iAFjo7^Aa46bN5xY8nh=O1SD)RCu9n7=Z8bh+#c&D>-87Q#co>gnpu`oAvx&ek4h zi@b9ix*HD_tH`%C9F`nrBb<5jX~>^R(r9KBeVDc4!sF>D?SThNZ#yg28gcKY2(uD6 zB^f}ppcIz~BmWJO!^mKW@(ck200dsi6FxFI%`wSJQkDQB;)~Vos}n42F(86PSy(kd zpwAtIxnb5Igp}kY2Qc%vBa_mBIZY`k34#>nhyYah;7QfJbOWA&6ygF?LduD_&OU`! zt5#@J3HK3Q=8>*4gLB!B%c|l@*9u+*(D;$*sHlY;7v{7Wf}rBOw-gK#*Lx;AVUg)_ z7X5Ix0(YstoqzrX)n!+MsqC&hN+tc>xIpIXb46SOx-~n9dhF(I1O8v92t`)D`!_}- zcl`@C2MMM)rvSbGW?@*&a>VYdXb?k;{?YlG?q-jhJf4K^D%IE%O<&i$-{f`o3w6F^ Iiunlu04+TjasU7T diff --git a/tests/io/aims/test_sets/test_md_generator.py b/tests/io/aims/test_sets/test_md_generator.py new file mode 100644 index 00000000000..13c6f8f6859 --- /dev/null +++ b/tests/io/aims/test_sets/test_md_generator.py @@ -0,0 +1,48 @@ +"""Tests the MD input set generator""" + +from __future__ import annotations + +from pathlib import Path + +import pytest +from pymatgen.io.aims.sets.core import MDSetGenerator +from pymatgen.util.testing.aims import Si, compare_files + +module_dir = Path(__file__).resolve().parents[1] +species_dir = module_dir / "species_directory" +ref_path = (module_dir / "aims_input_generator_ref").resolve() + + +def test_si_md(tmp_path): + # default behaviour + parameters = { + "species_dir": str(species_dir / "light"), + "k_grid": [2, 2, 2], + } + test_name = "md-si" + # First do the exceptions + with pytest.raises(ValueError, match="Ensemble something not valid"): + MDSetGenerator(ensemble="something").get_input_set(Si) + with pytest.raises(ValueError, match="Type parrinello is not valid for nve ensemble"): + MDSetGenerator(ensemble="nve", ensemble_specs={"type": "parrinello"}).get_input_set(Si) + with pytest.raises(ValueError, match="Velocities must be initialized"): + MDSetGenerator(ensemble="nve", init_velocities=False).get_input_set(Si) + with pytest.raises(ValueError, match="Temperature must be set"): + MDSetGenerator(ensemble="nve").get_input_set(Si) + with pytest.raises(ValueError, match="Temperature must be set"): + MDSetGenerator(ensemble="nvt").get_input_set(Si) + with pytest.raises(ValueError, match="parameter is not defined"): + MDSetGenerator(ensemble="nve", ensemble_specs={"type": "damped"}, temp=300).get_input_set(Si) + # then do the actual input set + generator = MDSetGenerator( + ensemble="nvt", + ensemble_specs={"type": "parrinello", "parameter": 0.4}, + temp=300, + time=10.0, + time_step=0.002, + user_params=parameters, + ) + input_set = generator.get_input_set(Si) + input_set.write_input(tmp_path / test_name) + + return compare_files(test_name, tmp_path, ref_path)