####-----Micro PMB(bond-base)-PD  simulation of Ti ----####

#### Initialization ##############################################
dimension       3
units           si
boundary        s s s
atom_style      peri
atom_modify     map array
neighbor        0.00085 bin

#### simulation parameters ########################################
variable   lat_const equal 0.0003 						        # m
variable   horizon equal ${lat_const}*3.00001
variable   volume equal (${lat_const}^3)

#### material properties ##########################################
#metal (Ti)
variable   E  equal 11.0e10 						        	# Pa, Young's modulus
variable   v  equal 0.33 							        # Poisson's ratio
variable   K  equal ((${E}*6*(1-4*${v}))/(3.14*${horizon}^4*(1-2*${v})*(1+${v}))) 	# Bulk Modulus
variable   C  equal (((1+${v})/(1-4*${v}))*${K})					# micro Modulus
variable   S  equal 0.011                                                               # s00
			 				
#### geometry ######################################################
read_data    atom.data
region       metal 	block EDGE EDGE      EDGE      EDGE    EDGE EDGE  units box
region       top 	block EDGE EDGE      0.0024    0.003   EDGE EDGE  units box
region       bottom     block EDGE EDGE     -0.003    -0.0024  EDGE EDGE  units box

group        topface    region top 
group        botface    region bottom

### initialize ######################################################
set group all        density 4357
set group all        volume ${volume}
velocity all set 0.0 0.0 0.0 sum no units box

### force field #####################################################
pair_style      peri/pmb

pair_coeff * *  ${C}   ${horizon}  ${S}  0.3

#### Set Bonds ######################################################
fix	10      all nve
run	2
unfix	10
timestep        1e-8
#####################################################################
fix   20 topface   move linear NULL 0.1 NULL units box
fix   30 botface   setforce  0.   0.  0.

### fixes ###########################################################
group mid          subtract all topface botface
fix   40           mid nve

### post-processing #################################################

##### position ###############################################
compute         dy         topface	reduce ave y
compute         dyb   	   botface      reduce ave y

##### energy #################################################
compute         peratom    all 		        pe/atom
compute         pe 	   all 		        reduce sum c_peratom

#### force ###################################################
compute         Yforce   mid                    reduce sum fy
compute 	Yforce1  topface 	        reduce sum fy
compute 	Yforce2  botface 	        reduce sum fy

### stress-strain ############################################
variable     len      equal  "ly" 
variable     L0       equal   ${len}
print       "Initial  Length, L0: ${L0}"
variable     strain   equal  "(ly - v_L0)/v_L0"
variable     Ny     equal  "v_strain"

compute      stpa     mid    stress/atom NULL
compute      st       mid    reduce sum  c_stpa[2]
variable     Sy       equal  (c_st/(0.8*vol))
  
### output #########################################################
thermo          1000
thermo_style    custom step c_dy c_dyb c_Yforce c_Yforce1  c_Yforce2 c_pe 
thermo_modify   norm no

dump   1        all custom 1000 Pd.dump1. id type xs ys zs 
dump   2        all       custom   1000 pd.dump2. id type x y z vx vy vz fx fy fz
dump_modify     1   format float %.3e first no

fix def1 all print 1000 "${Ny} ${Sy}" file deformation.txt title "Strain sigmayy" screen no

### run ############################################################
run            100000
print         "All done"