#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 2.07097278e-20
# 2.07e-20 J equals 5.00 kT (for T = 300 K)

shell cd ../

###----------INITIALIZATION----------###
units	si
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 7.28e-09

###-----------PAIR STYLE----------#
pair_style	lj/cut 2.390170685669979e-09

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

###----------NEIGHBOUR SETTINGS-----------#
neighbor	6.615e-09 multi
neigh_modify	every 1 delay 10 check yes one 512

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

bond_coeff	1 0.0503545233593595 2.028e-09  # 5.04E-02 J/m^2 equals 0.12 kT/A^2 (T = 300.0 K)

bond_coeff	2 0.0503545233593595 2.028e-09  # 5.04E-02 J/m^2 equals 0.12 kT/A^2 (T = 300.0 K)

bond_coeff	3 0.005725949892552008 6.014e-09  # 5.73E-03 J/m^2 equals 0.01 kT/A^2 (T = 300.0 K)

###----------PAIR SETTINGS-----------#
pair_coeff	1 1 2.07097278e-20 2.1294e-09 2.390170685669979e-09
pair_coeff	1 2 2.07097278e-20 2.1294e-09 2.390170685669979e-09
pair_coeff	1 3 2.07097278e-20 6.630435e-09 7.4424116512821575e-09
pair_coeff	2 2 2.07097278e-20 2.1294e-09 2.390170685669979e-09
pair_coeff	2 3 2.07097278e-20 6.630435e-09 7.4424116512821575e-09
pair_coeff	3 3 2.07097278e-20 1.1025000000000001e-08 1.2375144082610838e-08
pair_modify shift yes
pair_write 1 1 500 r 1e-10 100e-10 potential.txt LJ


###----------ANGLE INFO----------###
angle_style	harmonic
angle_coeff 1 2.2701047780769226e-21 180
# 2.27E-21 J/rad^2 equals 0.55 kT/rad^2 (T = 300.0 K)

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

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

###----------FIX COMMANDS BEFORE MINIMIZATION----------###
fix id_fix0 sb langevin 300 300 9.25550716157487e-14 1236 gjf vhalf
fix id_fix1 bb langevin 300 300 7.429614270957855e-14 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 -6.5125000000000006e-09 2.07097278e-20 1.1025000000000001e-08 1.2375144082610838e-08  # 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	5e-15
minimize 4.1419e-21                 4.1419e-21                 1000000                 1000000


###----------TIMESTEP----------###
timestep	3.3e-13 # Max stable dh without gjf half approx= 180fs.


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

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_si.in loop