In [ ]:
import os, sys, time
import numpy as np
import mapp4py
from mapp4py import md
In [ ]:
from mapp4py import mpi
if mpi().rank!=0:
with open(os.devnull, 'w') as f:
sys.stdout = f;
Import the initial configuration. This will construct an mapp4py.md.atoms object.
It is just a bcc unit cell with two atoms.
In [ ]:
sim = md.atoms.import_cfg("configs/Fe_300K.cfg"); # sim's type is md.atoms
In [ ]:
sim.kB = 8.617330350e-5
sim.hP = 4.13566766225 * 0.1 * np.sqrt(1.60217656535/1.66053904020)
In [ ]:
sim *= [8, 8, 8]; ## sim.__imul__([8, 8, 8])
In [ ]:
# let us load the force field
sim.ff_eam_setfl("potentials/FeH-sina.eam.alloy");
to make our example a bit more interesting let us remove 10.0 % of atoms randomly
In [ ]:
np.random.seed(4534534);
frac = 0.1;
def remove_random(x, x_d, elem, x_dof, id):
if np.random.random() < frac:
return False;
# def remove_random():
# if np.random.random() < frac:
# return False;
sim.do(remove_random)
in order for the system not to start from zero temperature we need to set the velocities using gaussian distribution
In [ ]:
sim.create_temp(300.0, 846244)
In [ ]:
# 0.1 [T] (see above) it is about 1 fs
nvt=md.nvt(300.0,0.1);
# print thermodynamic variables every 10000 steps
nvt.ntally = 100;
# dump the snapshots every 10000 steps
nvt.export = md.export_cfg('dumps/dump',100,sort=True)
In [ ]:
def mu(p,T):
return -2.37+0.0011237850013293155*T+0.00004308665175*T*np.log(p)-0.000193889932875*T*np.log(T);
In [ ]:
muvt=md.muvt(mu(1.0e-3,300.0),300.0,0.1,'H',73108204);
muvt.nevery=100;
muvt.nattempts=40000;
muvt.ntally=1000;
muvt.export=md.export_cfg('dumps/dump',10000)
In [ ]:
sh=float('nan')
nst=md.nst([[0.0],[sh,0.0],[sh,sh,0.0]],300,0.1);
nst.ntally=1000;
nst.export=md.export_cfg('dumps/dump',1000);
In [ ]:
min_cg=md.min_cg(e_tol=1.0e-8);
min_cg.ntally=0;
min_cg.ls=mapp4py.ls_brent();
min_cg.export=md.export_cfg('dumps/dump',1000);
In [ ]:
min_lbfgs=md.min_cg(e_tol=1.0e-8);
min_lbfgs.ntally=0;
min_lbfgs.ls=mapp4py.ls_brent();
min_lbfgs.export=md.export_cfg('dumps/dump',1000);
In [ ]:
nvt.run(sim,100)
In [ ]:
nst.run(sim,100)
In [ ]:
sim.add_elem('H',1.007940)
muvt.run(sim,nsteps);