######################################################
## INPUT AND OUTPUT FILES                           ##
### heating with restrain 20
######################################################
set input       min2
set output       heat
set inputname    $input
set outputname   $output
bincoordinates   $inputname.restart.coor
# COM : there will be no velocities from min2 simulation.. no dynamics
# COM : just minimization with constraints
#binvelocities    $inputname.restart.vel
extendedSystem   $inputname.restart.xsc
set restrain      restrain.pdb
coordinates      ref.pdb
structure        ref.psf
parameters       par_all27_prot_lipid.prm
#parameters       cholesterol.prm

######################################################
## INPUT SETTINGS                                   ##
######################################################
# COM : check the first time step for subsequent simulations
firsttimestep 0  
set temperature    0
set temperature1    300
set cons  1
#set start 1

######################################################
## JOB DESCRIPTION                                  ##
######################################################

#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input
paraTypeCharmm	    on
temperature         $temperature
binaryoutput        yes

# Force-Field Parameters
exclude             1-4
1-4scaling          1.0
cutoff              10.
switching           on
switchdist          8.5
pairlistdist        12.0

# Integrator Parameters
# COM : timestep is 1.0 during heating step
timestep            1.0  ;# 1fs/step
rigidBonds          all  ;# needed for 2fs steps
nonbondedFreq       1
fullElectFrequency  2  
stepspercycle       10

# for periodic boundary conditions
if {0} {
cellBasisVector1  60.0 0.0 0.0  
cellBasisVector2  0.0 77.0 0.0
cellBasisVector3  0.0 0.0 171.0
cellOrigin 0.0 0.0 0.0        
} 
wrapAll on
wrapWater on
DCDUnitCell     yes
margin 15

PME		yes 
#PMEGridSizeX 110 # grid-pts along cellBasisVector1.
#PMEGridSizeY 115 # along cellBasisVector2.
#PMEGridSizeZ 107 # along cellBasisVector3
PMEGridSpacing 1.0


# Output
outputName          $outputname

restartfreq         500 ;# 500steps = every 1ps
dcdfreq             500
xstFreq             500
outputEnergies      500
outputPressure      500
#velDCDfile					vel.dcd
#velDCDfreq          500
#############################################################
## EXTRA PARAMETERS                                        ##
#############################################################
langevin            on    ;# do langevin dynamics
langevinDamping     5     ;# damping coefficient (gamma) of 5/ps
langevinTemp $temperature
langevinHydrogen    off    ;# don't couple langevin bath to hydrogens

# Constant Pressure Control (variable volume)
# COM : heat in the NPT ensemble with flexible cell
 if {1} {
 useGroupPressure      yes ;# needed for 2fs steps
 useFlexibleCell       yes  ;# no for water box, yes for membrane
 useConstantArea       no  ;# no for water box, yes for membrane

 langevinPiston        on
 langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
 langevinPistonPeriod  200.
 langevinPistonDecay   50.
# langevinPistonTemp    $temperature
# COM : langevin piston temperature should be positive > 0
 langevinPistonTemp    1.0
 }


if {$cons == 1} {
        constraints on
	consexp 2
	consref   ref.pdb
	conskfile $restrain
	conskcol B
#	selectConstraints on
#	selectConstrX off
#	selectConstrY off
#	selectConstrZ on
} else {
        constraints off
}



#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################

# COM : put the constraint to 20 here.. heat it slowly to 300 K.
set ck 20
constraintScaling $ck

# COM: run dynamics at zero temperature for 500 steps
run 500
set freq 700
set teminc 5

for {set i 5} {$i <= $temperature1} {incr i $teminc} {
# Constant Temperature Control
langevinTemp $i
# Constant Pressure Control piston temperature
langevinPistonTemp $i
#COM : added this step to reinitialze velocities to desired temperature
reinitvels $i
run $freq
# 42 ps
}

# output equilibrated structures