<jupyter_start><jupyter_text>MD example - Lammps<jupyter_code># headers

# general modules
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import physical_constants

# pyiron modules
from pyiron_atomistics import Project
import pyiron_contrib

# conversion factors
KB = physical_constants['Boltzmann constant in eV/K'][0]
GPA_TO_BAR = 1e4
# project

pr = Project('md_test')
pr.remove_jobs_silently(recursive=True)
# basics

# structure specific
element = 'Al'
supercell = 3
cubic = True

# job specific
potential = '2008--Mendelev-M-I--Al--LAMMPS--ipr1'

# MD specific
temperature = 900. # K
pressure = 0.  # GPa
temperature_damping_timescale = 100.
pressure_damping_timescale = 1000.
timestep = 1.
n_steps = 5000
# create base structure

box = pr.create_ase_bulk(name=element, cubic=cubic).repeat(supercell)
n_atoms = box.get_number_of_atoms()
# create and run the NEB job

pr_prot = pr.create_group('protocol_md')  # create a new folder

md_ref = pr_prot.create_job(pr.job_type.Lammps, 'ref_md')
md_ref.structure = box
md_ref.potential = potential
md_ref.save()  # Don't forget this step!

# The same notebook should work for the other MD protocols:
# ProtoHarmMD: atoms behave as quantum-harmonic oscillators. Specify a spring or a force constant matrix!
# NOTE: ProtoHarmMD does not output pressure and volume, so some cells below may throw errors!
# ProtoConfinedMD: Regular MD, but atoms do not leave their lattice sites in the presence of a vacancy.
prot_md_job = pr_prot.create_job(pr.job_type.ProtoMD, 'prot_md_job')
prot_md_job.input.ref_job_full_path = md_ref.path
prot_md_job.input.structure = box
prot_md_job.input.temperature = temperature
prot_md_job.input.pressure = pressure  # not necessary for ProtoHarmMD!
prot_md_job.input.n_steps = n_steps
prot_md_job.input.temperature_damping_timescale = temperature_damping_timescale
prot_md_job.input.pressure_damping_timescale = pressure_damping_timescale  # not necessary for ProtoHarmMD!
prot_md_job.input.timestep = timestep

# prot_md_job.server.queue = 'cmti'
# prot_md_job.server.cores = 4
# prot_md_job.server.run_time = 3600

# set_output_whitelist sets how often an output of a particular vertex is stored in the archive.
# for example, here, the output 'volume' of vertex 'calc_static' is saved every 10 steps in the archive.
prot_md_job.set_output_whitelist(**
    {
    'barostat': {'pressure': 10},
    'verlet_velocities': {'energy_kin': 10},
    'calc_static': {'volume': 10}
    }
)  # ProtoHarmMD does not have 'barostat' and 'calc_static'!

prot_md_job.run()
# check what keys are in the output dictionary

prot_md_job.output.keys()
# run pyiron's lammps md

pr_pyir = pr.create_group('pyiron_md')

pyir_md_job = pr.create_job(pr.job_type.Lammps, 'pyir_md_job')
pyir_md_job.structure = box
pyir_md_job.potential = potential
pyir_md_job.calc_md(temperature=temperature, n_ionic_steps=n_steps, n_print=10, langevin=True)
# get lammps to use the berendsen barostat
pyir_md_job.input.control['fix'] = '1 all press/berendsen iso ' + str(pressure * GPA_TO_BAR) + \
                            ' ' + str(pressure * GPA_TO_BAR) + ' ' + str(pressure_damping_timescale)
pyir_md_job.run()
# main output protocol

protocol_volume = prot_md_job.graph.calc_static.archive.output.volume.data
protocol_pressure = prot_md_job.graph.barostat.archive.output.pressure.data
e_kin_md = np.array(prot_md_job.graph.verlet_velocities.archive.output.energy_kin.data)
protocol_temperature = (2. * prot_md_job.graph.verlet_velocities.archive.output.energy_kin.data) / (3 * KB * n_atoms)

onwards = 2
lattice_constant_protocol = np.cbrt(protocol_volume[onwards:]) / supercell
# main output lammps

lammps_volume = pyir_md_job.output.volume
lammps_pressure = [np.trace(p) for p in pyir_md_job.output.pressures]
lammps_temperature = pyir_md_job.output.temperature

lattice_constant_lammps = np.cbrt(lammps_volume[onwards:]) / supercell
# comparison plots

print('RESULTS\n')

plt.title('Lattice_constant')
plt.plot(lattice_constant_protocol[onwards:], label='protocol')
plt.plot(lattice_constant_lammps[onwards:], label='lammps')
print('Mean a protocol: ', np.mean(lattice_constant_protocol[onwards:]))
print('Mean a lammps: ', np.mean(lattice_constant_lammps[onwards:]))
print('\nStd a protocol: ', np.std(lattice_constant_protocol[onwards:]))
print('Std a lammps: ', np.std(lattice_constant_lammps[onwards:]))
plt.legend()
plt.show()

plt.title('Pressure')
plt.plot(protocol_pressure[onwards:], label='protocol')
plt.plot(lammps_pressure[onwards:], label='lammps')
print('Mean P protocol: ', np.mean(protocol_pressure[onwards:]))
print('Mean P lammps: ', np.mean(lammps_pressure[onwards:]))
print('\nStd P protocol: ', np.std(protocol_pressure[onwards:]))
print('Std P lammps: ', np.std(lammps_pressure[onwards:]))
plt.legend()
plt.show()

plt.title('Temperature')
plt.plot(protocol_temperature[onwards:], label='protocol')
plt.plot(lammps_temperature[onwards:], label='lammps')
print('Mean T protocol: ', np.mean(protocol_temperature[onwards:]))
print('Mean T lammps: ', np.mean(lammps_temperature[onwards:]))
print('\nStd T protocol: ', np.std(protocol_temperature[onwards:]))
print('Std T lammps: ', np.std(lammps_temperature[onwards:]))
plt.legend()
plt.show()<jupyter_output>RESULTS

Mean a protocol:  4.131081087250005
Mean a lammps:  4.134876908825429

Std a protocol:  0.005121686219488185
Std a lammps:  0.007008599297580372
{'fraction_non_alphanumeric': 0.11561228264890862, 'fraction_numerical': 0.021087680355160933, 'mean_word_length': 6.16158940397351, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '109977', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2098, 'n_tokens_neox': 1926, 'n_words': 407, 'path': 'notebooks/protocol_md_example.ipynb', 'repo_name': 'pyiron/pyiron_contrib', 'size': '5406'}