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

# Minimization and Equilibration of 
# Ubiquitin in generalized Born implicit solvent


#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

#structure          ionized.psf

parmfile comb.prmtop
coordinates        comb.pdb
#ambercoor comb.inpcrd
set inputname out_minimize_explicit
bincoordinates $inputname.coor
binvelocities $inputname.vel
extendedSystem $inputname.xsc

set temperature    310
set outputname     out_thermalize

firsttimestep      0


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

# Input
amber on

rigidBonds          water  ;#Let the hydrogens wiggle during minimization.
useSettle           on
rigidTolerance      1.0e-8
cutoff              9.0
pairlistdist        11.0
switching           off
exclude             scaled1-4
readexclusions      yes
1-4scaling          0.8333333
scnb                2.0
zeromomentum        on
ljcorrection        on
watermodel          tip4



#temperature         $temperature


# Force-Field Parameters

# Integrator Parameters
timestep            1.0  ;# 1fs/step
nonbondedFreq       1
fullElectFrequency  1
stepspercycle       10


# Constant Temperature Control
langevin            on    ;# do langevin dynamics
langevinDamping     1     ;# damping coefficient (gamma) of 1/ps
langevinTemp        $temperature
langevinHydrogen    off    ;# don't couple langevin bath to hydrogens

#d is recovered from the last line of the crd file. 
#I've increased it from 120.07 to 123 because I was getting very high pressures during minimization. 
#set d 123.0

#cellBasisVector1 $d                   0.0                                 0.0
#cellBasisVector2 [expr {-$d/3.0}]     [expr {2.0/3.0*sqrt(2.0)*$d}]       0.0
#cellBasisVector3 [expr {-$d/3.0}]     [expr {-sqrt(2.0)/3*$d}]            [expr {-sqrt(6.0)/3.0*$d}]

wrapWater on
wrapNearest on

PME yes
PMEGridSpacing 1.0

useGroupPressure yes
useFlexibleCell no
useConstantArea no

langevinPiston on 
langevinPistonTarget 1.01325
langevinPistonPeriod 100.0
langevinPistonDecay  50.0
langevinPistonTemp   $temperature


# Output
outputName          $outputname

restartfreq         500     ;# 500steps = every 500fs
dcdfreq             100
xstFreq             100
outputEnergies      100
outputPressure      100


#constraints on
#consref ionized.pdb
#conskfile ionized.pdb
#conskcol B

colvars on
colvarsConfig ends.col


#############################################################
## EXTRA PARAMETERS                                        ##
#############################################################


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

# Minimization
run 10000
#reinitvels          $temperature

#run 5000 
#5 ps thermalization.

#set constraintScale 0.5

#for {set i 0} {$i < 5} {incr i} {
#    constraintScaling $constraintScale
#    run 1000
#    set constraintScale [expr {$constraintScale / 2.0}]
#}