# 3d Uniaxial Tension & Compression of a Ta crystal (extended Finnis-Sinclair potential)
# Sample is monitored using several parameters.

#(1) Initialization

echo both
units		metal
boundary	p p p

atom_style	atomic
neighbor	1.0 bin
neigh_modify	every 1 delay 2 check yes

#(2) Atom definition 
## Nominal, standard orientation
##lattice		bcc  3.303 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
## To orient the box differently, still need orthogonal axes.
## Separation between planes:Sqrt(2),Sqrt(6),Sqrt(3)(ao units) to use PBC properly.
lattice         bcc 3.303 orient x 1 -1 0 orient y 1 1 -2 orient z 1 1 1 spacing 1.41421356 2.44948974 1.73205081

#I want a roughly cubic cell
## I want Lx=30 ao~10 nm =Nx ax ==> 30 ao = Nx *Sqrt(2) ao ==> 
## Nx=30/Sqrt(2)~21,2 ==> use Nx=22.
## Ly=Lx=30 ao=Ny ay=Ny ao Sqrt(6) ==> Ny=30/Sqrt(6)~12.25 ==> use 14.
## Lz = Lx =30 ao=Nz az=Nz Sqrt(3) ==> Nz=30/Sqrt(3)~17,3 ==> use 18. 
# for convenience, specially when slicing, I prefer to center the box.
region		box block -11 11 -7 7 -9 9

create_box	1 box
create_atoms	1 box

mass		1 180.948

#(3) Potentials
pair_style	eam/fs
pair_coeff	* *  ta.efs Ta 

#(4) Define groups
# no need to define group since atoms will be deleted
# Using units box can include the radius in A, instead of ao.
region          void sphere 0. 0. 0. 15.0 units box
delete_atoms    region void

#If a void is present, it might help to minimize the configuration 
#minimize 1.0e-4 1.0e-6 100 1000

#(5) Initial velocities
velocity	all create 300.0 482748 dist gaussian

## various calculations. 
## Stress does not include thermal, which is wrong when there is flow of material.
## Caveats for cna: it assumes an spherical neighbor shell, radius R~1.207ao~4.0 in this case. 
##i)If material is expanded so one should use an increased shell for large strains: 1.3R~5.2
##ii)For uniaxial expansion, neighbor shells are no longer spherical, and should correct for that.

compute 1 all ke/atom
compute 2 all centro/atom ## Note: need to modify centro routine for bcc structures.
compute 3 all pe/atom
compute 4 all stress/atom pair
compute 5 all cna/atom 4.0 
## here I R=4 A is in box units, not lattice units. 
## This is confirmed by the fact that all atoms have cna=3==bcc except near the void.

## Here I can calculate the average Von Mises stress for the sample.
variable J2 equal sqrt(0.5*((pxx-pyy)^2+(pxx-pzz)^2+(pzz-pyy)^2+6*(pxy^2+pxz^2+pyz^2)))

thermo 100
thermo_style custom step temp ke pe etotal press pxx pyy pzz pxy pxz pyz ly lx lz v_J2
thermo_modify	lost warn norm yes

#(6) Equilibrate
## I'd be better to use aniso and drag to accelerate relaxation
##fix             1 all npt 300.0 300.0 0.5 xyz 0.0 0.0 0.5
##
## Using nph might work better than npt in these cases far from equil. 
fix              2 all  temp/rescale 1 300.0 300.0 1.0 1.0
## here it assumes that the initial pressure before the fix is ~0.1 GPa=1000.
fix 3 all nph aniso 1000.0 0.0 1000.0 0.0 1000.0 0.0 0.30 drag 1.50
run		2000

##unfix           1
unfix		    2
unfix		    3

#(7) Loading 
reset_timestep	0

fix	INTEGRATE_FIX  all nve

## need to dump individual velocities and stress to obtain local temp and stress maps. 
## Also get cna values, since csp not modified in this executable
dump OUT1 all custom 10000 dump_Ta.* x y z c_5 c_3 c_1 vx vy vz id c_4[1] c_4[2] c_4[3] c_2
dump OUT2 all custom 1000 Ta-notBCC.* x y z c_5 c_3 c_1 vx vy vz id c_4[1] c_4[2] c_4[3] c_2
## modify to get only non-bcc atoms (cna!=3)
dump_modify OUT2 thresh c_5 != 3 
restart	        25000 tmp.restart.*

thermo 100

## 2 possibilities for deformation: erate and scale
## erate: L(t) = L0 (1 + erate*t)
fix DEFORM_FIX   all deform 1 z erate 0.001
## this gives strain rate: 10^-3/10^-12 ps=10^9/s. 
timestep	0.002
run		75000
##simulating 1 ps deforms 10^-3. 
##If I want to deform 10%=0.1, I need 100 ps = 50000 steps, because step Dt=2 fs.  
## 15% deformation:150 ps ==>75000 steps.

write_restart restart.final-tension

## Now hold the load
unfix DEFORM_FIX
run 5000

## Now unload, compressing the sample.
fix DEFORM_FIX_1   all deform 1 z erate -0.001
## this gives strain rate: 10^-3/10^-12 ps=10^9/s. 
timestep	0.002
run		75000

write_restart restart.final-unload

### END of the run