# VARIABLES
variable        data_name      index    sys.npt.data
variable        settings_name  index    sys.pnipam.water.settings
variable        nsteps         index    10000
variable        avg_freq       index    10000
variable        coords_freq    index    10000
variable        vel_freq       index    10000
variable        thermo_freq    index    10000
variable        dump4avg       index    5000
variable        Tinit          index    250
variable        T0             index    250
variable        Tf             index    250
variable        P0             index    1.0
variable        Pf             index    1.0
variable        vseed1         index    8453
variable        vseed2         index    8892

#===========================================================
# GENERAL PROCEDURES
#===========================================================
units		real	# m* = mass/m, x* = x/sig, E* = E/eps
dimension	3	    # 3 dimensional simulation
newton		on	    # use Newton's 3rd law
boundary	p p p	# periodic boundary conditions 
atom_style	full    # molecular + charge

#===========================================================
# FORCE FIELD DEFINITION
#===========================================================
pair_style      lj/cut/coul/long 10   # inner outer_LJ outer_Coul (cutoff values, see LAMMPS Doc)
bond_style      harmonic
angle_style     harmonic
dihedral_style  fourier
kspace_style    pppm         1.0e-05
improper_style  class2                          # class2 type impropers
pair_modify     tail    yes mix geometric          # mix geometric shift yes
special_bonds   lj/coul 0.0 0.0 0.5
neigh_modify    delay   0 every 1 check yes

#===========================================================
# SETUP SIMULATIONS
#===========================================================
# READ IN COEFFICIENTS/COORDINATES/TOPOLOGY
read_data ${data_name}
include ${settings_name}

# SET RUN PARAMETERS
timestep	1.0		    # real-time units
run_style	verlet 		# Velocity-Verlet integrator

# SET OUTPUTS
thermo_style    custom step temp vol density etotal press
thermo_modify   format float %14.6f
thermo ${thermo_freq}

# DECLARE RELEVANT OUTPUT VARIABLES
variable        my_step 	equal   step
variable        my_temp 	equal   temp
variable        my_rho  	equal   density
variable        my_etot   	equal   etotal
variable        my_P    	equal   press
variable        my_vol  	equal   vol

# group atoms together of each polymer 
group polymer id 1:572  
compute rog polymer gyration 

fix rog_averages polymer ave/time 1 1 5000 c_rog[0] c_rog[1] c_rog[2] c_rog[3] file rad_gyr.npt.avg

fix averages all ave/time ${dump4avg} $(v_avg_freq/v_dump4avg) ${avg_freq} v_my_temp v_my_vol v_my_rho v_my_etot v_my_P file thermo.npt.avg 

# CREATE COORDINATE DUMPS 
dump crds all atom  ${coords_freq}  coords.npt.lammpstrj
dump vel  all custom ${vel_freq}    vels.npt.lammpstrj id type vx vy vz
dump_modify crds sort id
dump_modify crds scale no
dump_modify vel  sort id

#===========================================================
# RUN DYNAMICS
#===========================================================
#===========================================================
# define velocities...
velocity all create ${Tinit} ${vseed1} mom yes rot yes

# group atoms of the water type...  
group spce type 10 11  

# fix dynamics and rigid topologies with rattle 
fix frattle spce shake 0.0001 20 0 b 9 a 15                     # print out bond and angle statistics to make sure water is being held rigid 
fix npt_dynamics all npt temp ${T0} ${Tf} 100.0 iso ${P0} ${Pf} 1000.0

# fix abf
fix abf all colvars colvars.inp tstat npt_dynamics

print "About to kick off npt..." append print.statements 

run 	${nsteps}

unfix abf

# write_data sys.npt.r1.data
# print "Ran NPT step!" append print.statements

# unfix stuff... 
# unfix npt_dynamics 
# unfix frattle