# Computes elastic waves propagating in and out of a finite temperature region
#AtC ThermoElastic Coupling
echo both

units		real
atom_style	atomic

# create domain
#lattice	type reduced density rho* = 4*(sigma/a)^3, where N=4 for fcc, s = 3.405 A (Wagner) and a = 5.25 A (Ashcroft & Mermin, p. 70)
lattice         fcc 5.2582305 origin 0.25 0.25 0.25

# NOTE following 3 lines added for restart
boundary        f p p
pair_style      lj/cut 13.
read_data       temp.init
region		simRegion block -12 12 -3 3 -3 3
# create atoms, NOTE commented out for restart
#region		mdRegion block -8 8 -3 3 -3 3
#boundary	f p p
#region		mdBox block -9 9 -3 3 -3 3
#create_box	1 mdBox
#create_atoms	1 region mdRegion
#mass		1 39.95

# specify interal/ghost atoms
region		mdInternal block -6 6 -3 3 -3 3
region		leftghost block -8 -6 -3 3 -3 3
region		rightghost block 6 8 -3 3 -3 3
group		internal region mdInternal
group		Lghost region leftghost
group		Rghost region rightghost
group		ghosts union Lghost Rghost

# velocities have Vcm = 0, NOTE next for lines commented out for restart
#velocity	internal create 40. 87287 mom yes loop geom
#pair_style	lj/cut 13.
##pair_coeff	1 1 0.010323166 3.405 13.
#pair_coeff  	1 1 .2381 3.405 13.

neighbor	5. bin
neigh_modify	every 10 delay 0 check no
#write_restart  tinit.dat

# zero initial momentum
fix             AtC internal   atc elastic     Ar_thermo_elastic.mat
#               ID  part keywords    nx ny nz region
fix_modify      AtC mesh create 12  1  1  simRegion f p p
fix_modify      AtC mesh create_faceset obndy box -6.0 6.0 -INF INF -INF  INF outward
fix_modify      AtC internal_quadrature off
fix_modify	AtC fix displacement x all 0.
fix_modify      AtC fix displacement y all 0.
fix_modify      AtC fix displacement z all 0.
fix_modify      AtC fix velocity x all 0.
fix_modify      AtC fix velocity y all 0.
fix_modify      AtC fix velocity z all 0.
fix_modify      AtC control momentum glc_velocity
#fix_modify      AtC  output bar1d_thermo_elastic_initializeFE 1 text #binary
timestep        0
thermo          1
run             1
unfix AtC

# define  layer
#               ID  group atc PhysicsType ParameterFile
fix             AtC internal   atc thermo_elastic     Ar_thermo_elastic.mat
fix_modify	AtC boundary ghosts

#               ID  part keywords    nx ny nz region
fix_modify      AtC mesh create 12  1  1  simRegion f p p
fix_modify      AtC mesh create_faceset obndy box -6.0 6.0 -INF INF -INF  INF outward
fix_modify      AtC internal_quadrature off

# initial conditions 
fix_modify  AtC  initial displacement x all 0.0
fix_modify  AtC  fix displacement y all 0.0
fix_modify  AtC  fix displacement z all 0.0
fix_modify  AtC  initial velocity x all 0.0
fix_modify  AtC  fix velocity y all 0.0
fix_modify  AtC  fix velocity z all 0.0
fix_modify  AtC  fix temperature all 20.

# set node sets and bcs
#           ID  mesh create_nodeset tag xmin xmax ymin ymax zmin zmax
fix_modify  AtC mesh create_nodeset lbc -12.1  -11.9   -INF INF  -INF INF
fix_modify  AtC mesh create_nodeset rbc  11.9   12.1   -INF INF  -INF INF
#fix_modify  AtC  fix velocity x rbc 0.
#fix_modify  AtC  fix displacement x rbc 0.
fix_modify  AtC  fix displacement x lbc 0.
fix_modify  AtC  fix velocity x lbc 0.

# rescale thermostat for initial atomic temperatures
fix_modify      AtC  control thermal rescale 10
fix_modify      AtC control momentum ghost_flux
fix_modify      AtC  output bar1d_thermo_elasticFE 100 text #binary
timestep        5
thermo          100
run             1000

# free all nodes
#fix_modify  AtC  unfix displacement x all
#fix_modify  AtC  unfix velocity x all
fix_modify  AtC  unfix temperature all

# new boundary conditions
fix_modify  AtC  fix velocity x rbc 0.00000004
fix_modify  AtC fix temperature lbc 20.


#fix_modify      AtC  output follow_ex.fe 50
#fix_modify	AtC  localized_lambda on
#fix_modify	AtC  momentum control glc_velocity
#fix_modify	AtC  momentum control flux faceset obndy
#fix_modify     AtC control momentum flux
#fix_modify	AtC control momentum ghost_flux
fix_modify	AtC control thermal flux
#fix_modify      AtC  filter scale 1000.0

# run to extension
compute		myTemp internal temp
compute		atomStress internal stress/atom NULL
compute		avgStress internal reduce sum c_atomStress[1] c_atomStress[2] c_atomStress[3]
variable	myPres equal -(c_avgStress[1]+c_avgStress[2]+c_avgStress[3])/(3*vol)
thermo_style	custom step c_myTemp v_myPres pe
fix_modify      AtC  output bar1d_thermo_elasticFE 500 text
#dump            D1 all custom 100 bar1d_thermo_elastic.dmp id type x y z vx vy vz
timestep	5
thermo		100
run 		1000

# change nodes to fixed
fix_modify      AtC  unfix velocity x rbc 0.
#fix_modify	AtC  fix displacement x rbc 0.0002

#fix_modify      AtC  output bar1d_fluxFE 500 text binary

# run to equilibrium
timestep        5
thermo		100
run 		3000