Source code for dpgen2.exploration.task.lmp.lmp_input
import random
from typing import (
List,
Optional,
)
import dpdata
import numpy as np
import scipy.constants as pc
from packaging.version import (
Version,
)
from dpgen2.constants import (
lmp_traj_name,
)
def _sample_sphere():
while True:
vv = np.array([np.random.normal(), np.random.normal(), np.random.normal()])
vn = np.linalg.norm(vv)
if vn < 0.2:
continue
return vv / vn
[docs]
def make_lmp_input(
conf_file: str,
ensemble: str,
graphs: List[str],
nsteps: int,
dt: float,
neidelay: Optional[int],
trj_freq: int,
mass_map: List[float],
temp: float,
tau_t: float = 0.1,
pres: Optional[float] = None,
tau_p: float = 0.5,
use_clusters: bool = False,
relative_f_epsilon: Optional[float] = None,
relative_v_epsilon: Optional[float] = None,
pka_e: Optional[float] = None,
ele_temp_f: Optional[float] = None,
ele_temp_a: Optional[float] = None,
nopbc: bool = False,
max_seed: int = 1000000,
deepmd_version="2.0",
trj_seperate_files=True,
):
if (ele_temp_f is not None or ele_temp_a is not None) and Version(
deepmd_version
) < Version("1"):
raise RuntimeError(
"the electron temperature is only supported by deepmd-kit >= 1.0.0, please upgrade your deepmd-kit"
)
if ele_temp_f is not None and ele_temp_a is not None:
raise RuntimeError(
"the frame style ele_temp and atom style ele_temp should not be set at the same time"
)
if "npt" in ensemble and pres is None:
raise RuntimeError("the pressre should be provided for npt ensemble")
ret = "variable NSTEPS equal %d\n" % nsteps
ret += "variable THERMO_FREQ equal %d\n" % trj_freq
ret += "variable DUMP_FREQ equal %d\n" % trj_freq
ret += "variable TEMP equal %f\n" % temp
if ele_temp_f is not None:
ret += "variable ELE_TEMP equal %f\n" % ele_temp_f
if ele_temp_a is not None:
ret += "variable ELE_TEMP equal %f\n" % ele_temp_a
if pres is not None:
ret += "variable PRES equal %f\n" % pres
ret += "variable TAU_T equal %f\n" % tau_t
if pres is not None:
ret += "variable TAU_P equal %f\n" % tau_p
ret += "\n"
ret += "units metal\n"
if nopbc:
ret += "boundary f f f\n"
else:
ret += "boundary p p p\n"
ret += "atom_style atomic\n"
ret += "\n"
ret += "neighbor 1.0 bin\n"
if neidelay is not None:
ret += "neigh_modify delay %d\n" % neidelay
ret += "\n"
ret += "box tilt large\n"
ret += (
'if "${restart} > 0" then "read_restart dpgen.restart.*" else "read_data %s"\n'
% conf_file
)
ret += "change_box all triclinic\n"
for jj in range(len(mass_map)):
ret += "mass %d %f\n" % (jj + 1, mass_map[jj])
graph_list = ""
for ii in graphs:
graph_list += ii + " "
if Version(deepmd_version) < Version("1"):
# 0.x
ret += "pair_style deepmd %s ${THERMO_FREQ} model_devi.out\n" % graph_list
else:
# 1.x
keywords = ""
if use_clusters:
keywords += "atomic "
if relative_f_epsilon is not None:
keywords += "relative %s " % relative_f_epsilon
if relative_v_epsilon is not None:
keywords += "relative_v %s " % relative_v_epsilon
if ele_temp_f is not None:
keywords += "fparam ${ELE_TEMP}"
if ele_temp_a is not None:
keywords += "aparam ${ELE_TEMP}"
ret += (
"pair_style deepmd %s out_freq ${THERMO_FREQ} out_file model_devi.out %s\n"
% (graph_list, keywords)
)
ret += "pair_coeff * *\n"
ret += "\n"
ret += "thermo_style custom step temp pe ke etotal press vol lx ly lz xy xz yz\n"
ret += "thermo ${THERMO_FREQ}\n"
if trj_seperate_files:
ret += "dump 1 all custom ${DUMP_FREQ} traj/*.lammpstrj id type x y z fx fy fz\n"
else:
ret += (
"dump 1 all custom ${DUMP_FREQ} %s id type x y z fx fy fz\n"
% lmp_traj_name
)
ret += "restart 10000 dpgen.restart\n"
ret += "\n"
if pka_e is None:
ret += 'if "${restart} == 0" then "velocity all create ${TEMP} %d"' % (
random.randrange(max_seed - 1) + 1
)
else:
sys = dpdata.System(conf_file, fmt="lammps/lmp")
sys_data = sys.data
pka_mass = mass_map[sys_data["atom_types"][0] - 1]
pka_vn = (
pka_e
* pc.electron_volt
/ (0.5 * pka_mass * 1e-3 / pc.Avogadro * (pc.angstrom / pc.pico) ** 2)
) # type: ignore
pka_vn = np.sqrt(pka_vn)
print(pka_vn)
pka_vec = _sample_sphere()
pka_vec *= pka_vn
ret += "group first id 1\n"
ret += 'if "${restart} == 0" then "velocity first set %f %f %f"\n' % (
pka_vec[0],
pka_vec[1],
pka_vec[2],
)
ret += "fix 2 all momentum 1 linear 1 1 1\n"
ret += "\n"
if ensemble.split("-")[0] == "npt":
assert pres is not None
if nopbc:
raise RuntimeError("ensemble %s is conflicting with nopbc" % ensemble)
if ensemble == "npt" or ensemble == "npt-i" or ensemble == "npt-iso":
ret += "fix 1 all npt temp ${TEMP} ${TEMP} ${TAU_T} iso ${PRES} ${PRES} ${TAU_P}\n"
elif ensemble == "npt-a" or ensemble == "npt-aniso":
ret += "fix 1 all npt temp ${TEMP} ${TEMP} ${TAU_T} aniso ${PRES} ${PRES} ${TAU_P}\n"
elif ensemble == "npt-t" or ensemble == "npt-tri":
ret += "fix 1 all npt temp ${TEMP} ${TEMP} ${TAU_T} tri ${PRES} ${PRES} ${TAU_P}\n"
elif ensemble == "nvt":
ret += "fix 1 all nvt temp ${TEMP} ${TEMP} ${TAU_T}\n"
elif ensemble == "nve":
ret += "fix 1 all nve\n"
else:
raise RuntimeError("unknown emsemble " + ensemble)
if nopbc:
ret += "velocity all zero linear\n"
ret += "fix fm all momentum 1 linear 1 1 1\n"
ret += "\n"
ret += "timestep %f\n" % dt
ret += "run ${NSTEPS} upto\n"
return ret