#Input script for LAMMPS, simulating a rigid finite-sized particle tethered to a surface with a semi-flexible polymer linker.

# fb_radius: 1.014e-09 m (fixed bead on surface)
# sb_radius: 1.014e-09 m (tether polymer beads)
# bb_radius: 5e-09 m (tethered nanoparticle)

###----------VARIABLES: MAIN----------###
variable quote string "'" 
variable lj_eps equal 1

shell cd ../

###----------INITIALIZATION----------###
units	lj
boundary	p p f
dimension	3
atom_style	hybrid angle sphere
log	35030_log.log 

#restart 20560164 RESTARTS/35030_restart*.restart 

###-----------COMM INFO----------#
# Make sure comm distance > cutoff + neighbor bin
comm_modify cutoff 3.6

###-----------PAIR STYLE----------#
pair_style	lj/cut 1.18

###---------READ DATAFILE---------#
read_data	35030_df_lj.df
echo screen

###----------NEIGHBOUR SETTINGS-----------#
neighbor	3.26 multi
neigh_modify	every 1 delay 10 check yes one 512

###----------BOND INFO-----------#
bond_style	harmonic

bond_coeff	1 10 1
bond_coeff	2 10 1
bond_coeff	3 1.14 2.965

###----------PAIR SETTINGS-----------#
pair_coeff	1 1 1 1.05 1.18
pair_coeff	1 2 1 1.05 1.18
pair_coeff	1 3 1 3.27 3.67
pair_coeff	2 2 1 1.05 1.18
pair_coeff	2 3 1 3.27 3.67
pair_coeff	3 3 1 5.43 6.10

pair_modify shift yes

###----------ANGLE INFO----------###
angle_style	harmonic
angle_coeff 1 0.11 180

###----------GROUP INFO----------###
group fb type 1
group sb type 2
group bb type 3

###----------VELOCITIES----------###
velocity all create 0.2 1236
velocity fb zero linear

###----------FIX COMMANDS BEFORE MINIMIZATION----------###
fix id_fix0 sb langevin 0.2 0.2 0.0053 1236 gjf vhalf
fix id_fix1 bb langevin 0.2 0.2 0.0042 1236 omega yes gjf vhalf
fix top sb wall/reflect zhi EDGE
fix lowwall sb wall/reflect zlo EDGE
#fix bblow bb wall/lj93 zlo -5.4 1 5.44 6.1 # Disable wall potential; first get gjf vhalf + sphere to work
fix	newton0 sb nve
fix	newton_sphere2 bb nve/sphere
velocity fb zero linear
fix zeroforce fb setforce 0 0 0

###----------COMPUTES----------###
compute ke1 all ke/atom
compute pebond1 all pe/atom bond
compute pepair1 all pe/atom pair

###----------MINIMIZATION----------###
timestep	0.000286
minimize 0.2                 0.2                 1000000                 1000000

###----------TIMESTEP----------###
timestep	0.01716

##-------THERMO OUTPUT--------###
thermo_style custom step temp etotal ke pe
thermo 822407

###----------RUNS & DUMP ----------###
variable maxRunsPerDump equal 2000000000
variable totalRuns equal 21000000
variable noSubRuns equal ceil(${totalRuns}/${maxRunsPerDump})
variable remainingRuns equal ${totalRuns}

variable i loop ${noSubRuns}

print ${noSubRuns}

label loop
dump	1 all custom 8224065 35030_dump$i.dump id type x y z vx vy vz fx fy fz omegax omegay omegaz c_ke1 c_pepair1 c_pebond1 id mass
dump_modify 1 sort id pad 12 delay 0
if "${remainingRuns}>${maxRunsPerDump}" then &
 "run ${maxRunsPerDump}" &
 "variable remainingRuns equal ${remainingRuns}-${maxRunsPerDump}" &
elif ${remainingRuns}<=0 &
 "undump 1" &
 "next i" &
else &
 "run ${remainingRuns}" &
 "variable remainingRuns equal 0" 
undump 1
next i
jump 35030_input_lj.in loop