echo		screen
units		real
atom_style	full
boundary p p p

# Variable definitions

variable	project		index	"polymer"	# project name
variable	source		index	.	# data directory
variable	params		index	.	# parameter directory
variable	temperature	index	303		# system temperature
variable	tdamp		index	50		# temperature damping
variable    	pdamp           index   500             # pressure damping
variable	dielectric	index	1.0		# medium dielectric
variable	kappa		index	4		# electrostatics kappa
variable	cutoff		index	12.0		# standard cutoff
variable	charge_cutoff	index	12.0		# charge cutoff
variable	precision	index	0.00001		# kspace precision
variable	lseed		index	723853		# langevin seed
variable	vseed		index	486234		# velocity seed
variable	trun		index	40000000	# run time
variable	frestart	index	0		# 0: equil, 1: restart
variable	dtrestart	index	10000		# delta restart time
variable	dtdump		index	10000		# delta dump time
variable	dtthermo	index	10000		# delta thermo time
variable	timestep	index	0.5		# integration time step
variable	restart		index	${params}/${project}.restart

if "${frestart} != 0" then &
"variable	data		index	${restart}" &
else &
"variable	data		index	${params}/${project}.data" &


# Interaction potential definition

pair_style lj/class2/coul/cut ${cutoff} ${charge_cutoff}
bond_style	class2
angle_style	class2
dihedral_style	class2
improper_style	class2
pair_modify	mix sixthpower tail yes
special_bonds	lj/coul 0 0 1

if "${frestart} != 0" then "read_restart ${data}" else "read_data ${data}"
include		${params}/${project}.params

neighbor 2.0 bin
neigh_modify  delay 0 every 1 check yes

# Integration conditions (check)

timestep	${timestep}
dielectric	${dielectric}
fix		mom all momentum 100 linear 1 1 1 angular

if "${frestart} != 0" then "jump SELF simulate"

group liion type 1
group otfion molecule 31:60
group backbone molecule 1:30
group ion union liion otfion


# Minimization
reset_timestep 0 
thermo 10 
thermo_style custom step temp pe etotal lx ly lz vol density press pxx pyy pzz

fix min all box/relax aniso 0.0 vmax 0.001
min_style sd 1e-5
minimize 1e-25 1e-25 5000 10000 

min_style cg 1e-5
minimize 1e-25 1e-25 5000 10000
unfix min

# specify initial velocity of atoms
velocity	all create ${temperature} ${vseed} &
		dist gaussian rot yes mom yes sum yes

reset_timestep 0
thermo		${dtthermo}
dump           1 all custom 5000 equilibration.lammpstrj id type x y z

#room temp and pressure
fix npt1 all npt temp ${temperature} ${temperature} ${tdamp} aniso 1.0 1.0 ${pdamp}  
run 1000000
unfix npt1

# constant room temp
fix nvt_eq all nvt temp ${temperature} ${temperature} ${tdamp}
run 4000000
unfix nvt_eq

#room temp and pressure
fix npt2 all npt temp ${temperature} ${temperature} ${tdamp} aniso 1.0 1.0 ${pdamp}  
run 5000000
unfix npt2

undump 1
write_restart	${project}.restart

# Simulation

label		simulate

# Integrator
reset_timestep 0

pair_style lj/class2/coul/long ${cutoff} ${charge_cutoff}
bond_style	class2
angle_style	class2
dihedral_style	class2
improper_style	class2
pair_modify	mix sixthpower tail yes
special_bonds	lj/coul 0 0 1

include		${params}/${project}.params
kspace_style pppm/cg ${precision}

fix             nvt3 all nvt temp ${temperature} ${temperature} ${tdamp}
fix            3 all efield 5e-8 0.0 0.0
fix_modify 3 energy yes
fix_modify 3 virial yes

run ${trun}