Ru/H Part 3: Initial reference data MD simulation + single-point replays¶
Important
First read Ru/H introduction and follow all parts in order.
To follow along,
Download
03_Ru_H2_gas_snapshots.ipynb(see also: how to install Jupyterlab in AMS)
Initial imports¶
import scm.plams as plams
from scm.params import ResultsImporter
from scm.plams import Settings, AMSJob, log, Molecule, packmol_on_slab
from pathlib import Path
import matplotlib.pyplot as plt
# common_ru_h.py must exist in the current working directory
from common_ru_h import (
rotation,
dft_settings,
QEKPointsConfig,
m3gnet_up_settings,
replay_settings,
slice_slab,
check_installation,
)
Initialize PLAMS working directory¶
old_ref_dir = "reference_data_2"
check_installation(old_ref_dir)
new_ref_dir = "reference_data_3"
ri = ResultsImporter.from_yaml(old_ref_dir)
plams.init()
Current AMS version: 2024.102
05-31 14:23:05 m3gnet is installed: M3GNet ML Backend v[0.2.4] - build:0 [06668e0a45ce742d8f66ff23484b8a1e]
05-31 14:23:06 qe is installed: Quantum ESPRESSO (AMSPIPE) v[7.1] - build:115 [777d72eb480fe4d632a003cc62e9c1cb]
PLAMS working folder: /path/plams_workdir.003
Construct the Ru(10-10)/H2(gas) interface¶
For details about the construction of the slab, see the previous notebook.
optimized_bulk = ri.job_collection["hcp_lattopt_Ru_dft"].molecule
slab_100 = slice_slab(
optimized_bulk, miller=(1, 0, 0), thickness=7.0, cell_z=15, ref_atom=0
)
slab_100 = slab_100.supercell(3, 2, 1)
for at in slab_100:
at.properties = Settings() # remove details about supercell generation
plams.plot_molecule(slab_100, rotation=rotation)
plt.title("Ru(10-10)");
Now use the packmol_slab function to add hydrogen molecules:
from scm.plams import packmol_on_slab
h2_mol = plams.from_smiles("[HH]")
density = 0.3 # approximate density of the gas phase in g/cm^3
slab_100_H2_gas_raw = packmol_on_slab(slab_100, h2_mol, density=density)
slab_100_H2_gas = plams.preoptimize(
slab_100_H2_gas_raw, settings=m3gnet_up_settings(), maxiterations=100
)
plams.plot_molecule(slab_100_H2_gas, rotation=rotation, radii=0.8)
Now let’s run some MD:
mdjob = plams.AMSNVTJob(
settings=m3gnet_up_settings(),
name="md_ru10-10_h2_m3gnet",
temperature=501,
nsteps=10000,
molecule=slab_100_H2_gas,
samplingfreq=100,
)
mdjob.run();
[31.05|14:23:35] JOB md_ru10-10_h2_m3gnet STARTED
[31.05|14:23:35] JOB md_ru10-10_h2_m3gnet RUNNING
[31.05|15:08:27] JOB md_ru10-10_h2_m3gnet FINISHED
[31.05|15:08:28] JOB md_ru10-10_h2_m3gnet SUCCESSFUL
from scm.params import ResultsImporter
ri = ResultsImporter.from_yaml(old_ref_dir)
settings = dft_settings(QEKPointsConfig(3, 3, 1))
settings += replay_settings(mdjob.results.rkfpath(), frames=[10, 30, 50, 70, 90])
dft_replay_job = plams.AMSJob(
settings=settings,
name="snapshots_from_md_ru10-10_h2_dft",
)
dft_replay_job.run(watch=True)
[31.05|15:08:28] JOB snapshots_from_md_ru10-10_h2_dft STARTED
[31.05|15:08:28] JOB snapshots_from_md_ru10-10_h2_dft RUNNING
[31.05|15:08:29] snapshots_from_md_ru10-10_h2_dft: AMS 2024.102 RunTime: May31-2024 15:08:29 ShM Nodes: 1 Procs: 1
[31.05|15:08:29] snapshots_from_md_ru10-10_h2_dft: Starting trajectory replay ...
[31.05|15:08:29] snapshots_from_md_ru10-10_h2_dft: Replaying frame #1/5 (#10 in original trajectory)
[31.05|15:08:29] snapshots_from_md_ru10-10_h2_dft: NOTE: a single QE.Label is assigned to atoms of different species.
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: AMS Pseudopotentials Finder
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: ---------------------------
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: PP library: /home/user/.scm/packages/AMS2024.1.packages/qe-_5399rm6/content/upf_files
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: * Ru
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: Path: GGA/PBE/SR/SSSP_Efficiency_v1.3.0/UPFs
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: File: Ru_ONCV_PBE-1.0.oncvpsp.upf
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: * H
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: Path: GGA/PBE/SR/SSSP_Efficiency_v1.3.0/UPFs
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: File: H.pbe-rrkjus_psl.1.0.0.UPF
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: recommended ecutwfc: 60.0 Ry
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: recommended ecutrho: 480.0 Ry
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: recommended ratio: 8.0x
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: NOTE: System%ecutwfc (40.0 Ry) is below the recommended threshold (60.0 Ry)
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: NOTE: System%ecutrho not specified in input. Using recommended ecutrho to ecutwfc ratio: ecutrho = 8.0 * ecutwfc = 320.0 Ry
[31.05|15:08:30] snapshots_from_md_ru10-10_h2_dft: Starting QuantumEspresso worker ... (see output file for progress information)
[31.05|15:14:14] snapshots_from_md_ru10-10_h2_dft: Replaying frame #2/5 (#30 in original trajectory)
[31.05|15:20:39] snapshots_from_md_ru10-10_h2_dft: Replaying frame #3/5 (#50 in original trajectory)
[31.05|15:28:07] snapshots_from_md_ru10-10_h2_dft: Replaying frame #4/5 (#70 in original trajectory)
[31.05|15:37:02] snapshots_from_md_ru10-10_h2_dft: Replaying frame #5/5 (#90 in original trajectory)
[31.05|15:45:36] snapshots_from_md_ru10-10_h2_dft: Trajectory replay complete.
[31.05|15:46:36] snapshots_from_md_ru10-10_h2_dft: Number of QE evaluations: 5
[31.05|15:46:36] snapshots_from_md_ru10-10_h2_dft: NORMAL TERMINATION
[31.05|15:46:36] JOB snapshots_from_md_ru10-10_h2_dft FINISHED
[31.05|15:46:36] JOB snapshots_from_md_ru10-10_h2_dft SUCCESSFUL
<scm.plams.interfaces.adfsuite.ams.AMSResults at 0x14c005cbf7c0>
ri.add_trajectory_singlepoints(dft_replay_job, properties=["energy", "forces"])
ri.store(new_ref_dir)
['reference_data_3/job_collection.yaml',
'reference_data_3/results_importer_settings.yaml',
'reference_data_3/training_set.yaml']
Complete Python code¶
#!/usr/bin/env amspython
# coding: utf-8
# ## Initial imports
import scm.plams as plams
from scm.params import ResultsImporter
from scm.plams import Settings, AMSJob, log, Molecule, packmol_on_slab
from pathlib import Path
import matplotlib.pyplot as plt
# common_ru_h.py must exist in the current working directory
from common_ru_h import (
rotation,
dft_settings,
QEKPointsConfig,
m3gnet_up_settings,
replay_settings,
slice_slab,
check_installation,
)
# ## Initialize PLAMS working directory
old_ref_dir = "reference_data_2"
check_installation(old_ref_dir)
new_ref_dir = "reference_data_3"
ri = ResultsImporter.from_yaml(old_ref_dir)
plams.init()
# ## Construct the Ru(10-10)/H2(gas) interface
#
# For details about the construction of the slab, see the previous notebook.
optimized_bulk = ri.job_collection["hcp_lattopt_Ru_dft"].molecule
slab_100 = slice_slab(optimized_bulk, miller=(1, 0, 0), thickness=7.0, cell_z=15, ref_atom=0)
slab_100 = slab_100.supercell(3, 2, 1)
for at in slab_100:
at.properties = Settings() # remove details about supercell generation
plams.plot_molecule(slab_100, rotation=rotation)
plt.title("Ru(10-10)")
# Now use the ``packmol_slab`` function to add hydrogen molecules:
from scm.plams import packmol_on_slab
h2_mol = plams.from_smiles("[HH]")
density = 0.3 # approximate density of the gas phase in g/cm^3
slab_100_H2_gas_raw = packmol_on_slab(slab_100, h2_mol, density=density)
slab_100_H2_gas = plams.preoptimize(slab_100_H2_gas_raw, settings=m3gnet_up_settings(), maxiterations=100)
plams.plot_molecule(slab_100_H2_gas, rotation=rotation, radii=0.8)
# Now let's run some MD:
mdjob = plams.AMSNVTJob(
settings=m3gnet_up_settings(),
name="md_ru10-10_h2_m3gnet",
temperature=501,
nsteps=10000,
molecule=slab_100_H2_gas,
samplingfreq=100,
)
mdjob.run()
from scm.params import ResultsImporter
ri = ResultsImporter.from_yaml(old_ref_dir)
settings = dft_settings(QEKPointsConfig(3, 3, 1))
settings += replay_settings(mdjob.results.rkfpath(), frames=[10, 30, 50, 70, 90])
dft_replay_job = plams.AMSJob(
settings=settings,
name="snapshots_from_md_ru10-10_h2_dft",
)
dft_replay_job.run(watch=True)
ri.add_trajectory_singlepoints(dft_replay_job, properties=["energy", "forces"])
ri.store(new_ref_dir)