# 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

region		simRegion block -12 12 -3 3 -3 3
region		mdRegion block -8 8 -3 3 -3 3
boundary	f p p
create_box	1 simRegion
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
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

# define  layer
#               ID  group atc PhysicsType ParameterFile
fix             AtC internal   atc elastic     Ar_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
fix_modify  AtC  internal_atom_integrate off
fix iNVE internal nve

# 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

# 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.

# turn on multiscale
fix_modify      AtC control momentum ghost_flux

# new boundary conditions
fix_modify  AtC  fix velocity x rbc 0.00000004


#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  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_ghost_fluxFE 500 text
#dump            D1 all atom 100 bar1d_ghost_flux.dmp
timestep	5
thermo		100
run 		1000

# change nodes to fixed
fix_modify      AtC  unfix velocity x rbc
#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
