# script for mgpt t=0 eos with unrelaxed vacancy in bcc lattice:
# input for unrelaxed vacancy formation energy at constant atomic volume

echo screen

units		electron
atom_style	atomic

# Atomic volume for MGPT potential in a.u.
variable        atomic_vol equal 121.6

# Derive effective lattice volume from atomic volume for 249-site cell
variable        lat_vol equal ${atomic_vol}*249/250

# Derive lattice constant from lattice volume
variable        lattice_constant equal (${lat_vol}*2.0)^(1.0/3.0)

# Create bcc lattice with 5x5x5 unit cells (250 atoms)
lattice		bcc ${lattice_constant}
region		box block 0 5 0 5 0 5
create_box	1 box
create_atoms	1 box

# Remove central atom from bcc lattice to create vacancy
region vacancy sphere 2.5 2.5 2.5 0.1 units lattice
delete_atoms region vacancy

# Define potential for use in simulation
pair_style      mgpt

# Set parameters for potential:
#                  parameter files  atomic volume 
#pair_coeff      * * parmin potin    ${atomic_vol} 
pair_coeff      * * Ta6.8x.mgpt.parmin Ta6.8x.mgpt.potin  ${atomic_vol}

# Create velocities at 0 K
velocity	all create 0.0 87287

# Set neighbor list parameters
neighbor	0.1 bin
neigh_modify	every 1 delay 0 check yes

# Set up microcanonical integrator
fix		1 all nve

# Dump coordinates to file every 50 timesteps
dump		id all atom 50 dump.vac0-bcc

# Output thermodynamical data every 10 timesteps
thermo		10

# Set output quantities and output format
thermo_style custom step vol temp pe etotal press

## Example: Output floating point number with 5 digits exponential notation.
#thermo_modify format float %15.5e

# Run 0 timesteps
run		0
 
# Convert energy to rydbergs and pressure to gpa

variable natoms equal "count(all)"
variable voltot equal "vol"
variable atvol equal "v_voltot/v_natoms"
variable etot equal "2.0*pe"
variable etotry equal "v_etot/v_natoms"
variable ptot equal "press"
variable ptotgpa equal "v_ptot/1.0e+09"

print "number of atoms = ${natoms}"
print "atomic volume (a.u.) = ${atvol}"
print "total energy (ry/atom) = ${etotry}"
print "pressure (gpa) = ${ptotgpa}"
print "${natoms}   ${atvol}   ${etot}   ${ptotgpa}"
print "${atvol}   ${etotry}   ${ptotgpa}"
