# needs description
#AtC Thermal 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

# create atoms
region		simRegion block -12 12 -3 3 -3 3
region		atomRegion block -9 9 -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
#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

# define  layer
#               ID  group atc PhysicsType ParameterFile
fix             AtC internal   atc elastic     Ar_elastic.mat
#fix_modify	AtC boundary Lghost
#fix_modify	AtC boundary Rghost
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

# initial conditions 
fix_modify  AtC  initial displacement x all 0.0
fix_modify  AtC  initial displacement y all 0.0
fix_modify  AtC  initial displacement z all 0.0
fix_modify  AtC  initial velocity x all 0.0
fix_modify  AtC  initial velocity y all 0.0
fix_modify  AtC  initial 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.00000004
#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.


#fix_modify      AtC  output follow_ex.fe 50
fix_modify     AtC  internal_quadrature off
#fix_modify	AtC  control lumped_lambda_solve on
#fix_modify	AtC  momentum control glc_velocity
#fix_modify	AtC  momentum control flux faceset obndy
fix_modify     AtC control momentum flux interpolate
#fix_modify      AtC  filter scale 1000.0

# run to extension
compute		myTemp internal temp
compute		atomStress internal stress/atom
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_fluxFE 10 text
timestep	5
thermo		100
run 		1000

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

fix_modify      AtC  output bar1d_fluxFE 500 text

# run to equilibrium
timestep        5
thermo		100
run 		10000