#-----INITIALIZATION-----
clear
units		real
atom_style	full
timestep 0.1        #Timestep is 0.1 fs
boundary p p p
#-----Force Fields-----
pair_style	lj/cut/coul/long 2.5 2.5
pair_modify	tail yes
bond_style	harmonic
angle_style	harmonic
improper_style cvff
kspace_style	pppm 1.0e-4
read_data	cyclopentanol.txt  group cyclopentanol extra/bond/types 19 extra/atom/types 19 extra/angle/types 37 extra/improper/types 12
replicate	5 6 4
read_data	cpme.txt add append toff 16 boff 16 aoff 31 ioff 10 
variable dt equal 0.1
variable Tdamp equal 100*${dt}
variable Pdamp equal 100*${dt}
#-----Equilibriation-----
neighbor 2.0 bin
neigh_modify delay 1 every 1 check no

#----Run1 for equilibration-----
velocity cyclopentanol create 300 12345678 dist uniform
velocity cpme create 300 12345678 dist uniform
run 0
velocity cyclopentanol scale 500
velocity cpme scale 500
minimize	1e-20 1e-20 10000 10000
fix 1 all nve
thermo 100
thermo_style custom step temp press vol pe ke etotal
dump  2 all atom 100  dump_nve_cyclopentanol_equilibriation_run.atom
dump 3 all cfg 1000 dump.equi_*.cfg mass type xs ys zs
dump_modify 3 element C1 C2 C3 C4 C5 O6 H7 H8 H9 H10 H11 H12 H13 H14 H15 H16 C17 O18 C19 C20 C21 C22 C23 H24 H25 H26 H27 H28 H29 H30 H31 H32 H33 H34 H35
run 100000           
unfix 1
undump 3
undump 2
print "NVE Equilibration run completed"
write_restart nve_equi_all.dat 

#------Run2 for equilibration-------
fix 44 all npt temp 300 300 ${Tdamp} iso 1.0 1.0 ${Pdamp} pchain 10        #Starting pressure is determined from previous NVE run
fix 4 all press/berendsen iso 1.0 1.0 10000
#compute 5 all temp
#compute 6 all pressure
thermo 500
thermo_style custom step temp press vol pe ke etotal

#----Output Style----------
variable        Nevery equal 20                      #Specify the distance (in timesteps) between samples for computing ensemble aver$
                                                     #  (Interval between samples in a block)
variable        Nrepeat equal 50                     #Specify the number of samples per output  of thermodynamic averages
                                                     #  (Samples per block output)
variable        Nfreq equal ${Nevery}*${Nrepeat}     #Specify the dump interval (in timesteps)
                                                     #  (Some people call these blocks)
variable        PotentialEnergy equal pe          #Call the intermolecular potential energy $PotentialEnergy
variable        Pressure equal press                 #Call the pressure $Pressure
variable        Temperature equal temp               #Call the instantaneous temperature $Temperature
#special_bond coul 1e-50 1e-50 1e-50
compute myrdf all rdf 150 cutoff 2.5
fix             12 all ave/time ${Nevery} ${Nrepeat} ${Nfreq} v_Temperature v_PotentialEnergy v_Pressure file ave2_all.out format %.8g  
fix		13 all ave/time ${Nevery} ${Nrepeat} ${Nfreq} c_myrdf[*] file tmp1_all.rdf mode vector
#fix		colvr cyclopentanol colvars colvars.in
dump  14 all atom 50000  dump_npt_cyclopentanol_production_run.atom                         ## output final configuration into dump_NEW.atom
dump 15 all cfg 1000000 dump.prod_*.cfg mass type xs ys zs
dump_modify 15 element C1 C2 C3 C4 C5 O6 H7 H8 H9 H10 H11 H12 H13 H14 H15 H16 C17 O18 C19 C20 C21 C22 C23 H24 H25 H26 H27 H28 H29 H30 H31 H32 H33 H34 H35
                                                     #  (MD step, temperature, internal energy, pressure)    
run 1000000            #Total production time 1 ns
undump 14
undump 15
unfix 4
unfix 44
unfix 12
unfix 13
write_restart final_cyclopentanol.restart
write_data final_cyclopentanol.dat
print "The system is equilibrated and stable along with completion of Run4"
write_restart final_cyclopentanol.restart
#----PMF calculation----
#fix		colvr cyclopentanol colvars colvars.in input out.colvars.state output out2
#run 100