#3. Step: No constraints
#Minimization and heating with constraints
# NO PME (using electrostatic cutoff switch)
#SHAKE for protein, SETTLE for waters, 2fs timesteps
# Using psf prepared by Gernot on 12.08.2010, where patching of GLUP was followed by AUTO ANGL DIHE

set temperature 300 ;# target temperature
              set propsteps 48000 ;# propagation steps
seed            28885

# initial config
coordinates     eclc_g11_inwaterbox.pdb #Just for setup, WILL BE OVERWRITTEN BY bincoordinates
bincoordinates  mini_and_heating_afterstep2.coor # Actually used coordinates
# firsttimestep 1
# binvelocities xxx.vel

# output params - this block specifies the output
outputname      mini_heating_propagation_afterstep3 ;# The file name (without extension!), to which final coordinates and velocities are written #(appended extensions are *.coor or *.vel)
dcdfreq 	100 ;# Frequency of writing structural snapshots to dcd trajectory file
binaryoutput    yes ;# is binary format used for saving last structure?

outputenergies 100;
# Interval, in integration steps, of writing energies to stdout
outputtiming 100 ;# Interval of writing timing information (speed, memory allocation, etc) to stdout
restartname mini_heating_propagation_step3 ;# The file name (without extension), which holds the restart structure and velocities (appended extensions are *.coor or *.vel)
restartfreq 100000 ;# Interval between writing out the restart
# coordinates and velocities
binaryrestart yes ;# are restart files in binary format?
DCDfile mini_heating_propagation_step3_noconstr.dcd ;# Trajectory filename (binary file)

# integrator params
timestep        2.0 ;

### structure and force field params
structure       eclc_g11_inwaterbox.xplor.psf ;
paraTypeCharmm  on ;
parameters      par_all22_prot.inp ;
exclude         scaled1-4 ;# Exclude/scale local (along the sequence) 
1-4scaling      1.0 ;# Scale factor for (i,i+4) EL interactions

###Periodic Boundary Conditions
# this block defines periodic boundary conditions......
cellBasisVector1 110.0 0.0 0.0
# Direction of the x basis vector for a unit cell
cellBasisVector2 0.0 110.0 0.0
# Direction of the y basis vector for a unit cell
cellBasisVector3 0.0 0.0 110.0
# Direction of the z basis vector for a unit cell
cellOrigin 0.0 0.0 0.0
# Position of the unit cell center
wrapWater on  	# Are water molecules translated back to the unit cell (purely cosmetic option, has no effect on simulations
wrapAll off     # Are other molecules such as proteins translated back 

###Electrostatics and Van der Waals
dielectric	1.0 ;# Value of the dielectric constant
switching       on ;# Energy Switch VdW interactions and partition EL into local and non-local contributions
switchdist      10.0 ;# distance at which to activate switching function for van der Waals (and electrostatic in sph. cut.)
cutoff          12.0    ;# In PME, cutoff dictates the separation between long and short range forces and doesn't simply turn off interactions.
pairlistdist    13.5    ;# CUTNB in charmm
stepspercycle   18	;# INBFRQ in charmm, Frequency of updating Verlet list (in integration steps) Def. 20
#fullElectfrequency 8 ;# number of timesteps between full electrostatic evaluations Acceptable Values. positive integer factor of stepspercycle
fullElectfrequency 3 ;# number of timesteps between full electrostatic evaluations Acceptable Values. positive integer factor of stepspercycle
#margin          0.0 ;# Extra distance used in selecting patches
#nonbondedFreq 1      ;# timesteps between nonbonded evaluation. Positive integer factor of fullElectFrequency.

### Particle Mesh Ewald
PME off
#PME on ;# Use PME for electrostatic calculation - USe only factors of 2,3,5
#PMEGridSizeX 100
#PMEGridSizeY 100
#PMEGridSizeZ 100
##ischrgd=1, ;# NAMD doesnt force neutralization of char_charge
#PMETolerance 0.000001 ;# Def 1E-6
#PMEInterpOrder 4 ;# Def 4 (cubic+1)

### SETTLE
rigidBonds all ;# Use SHAKE on solute (and SETTLE on water according to useSettle)
useSettle on; #Apply SETTLE (implemented for water)
rigidTolerance 0.00001 ;# Desired accuracy in maintaining SHAKEed bond lengths
rigidIterations 100 ;# Maximum number of SHAKE iterations 

### Temperature control: Langevin
langevin on ;#Langevin dynamics active. 
langevinTemp $temperature ;#This temperature will be roughly maintained across the selected atoms through the addition of friction and random forces.
langevinDamping 2. ;# damping coefficient of 2/ps
langevinHydrogen no ;# don't couple bath to hydrogens

###Pressure control: "extendedSystem" $NAME.xsc holds the periodic cell dimensions to continue after a constant 
#pressure simulation. If present, cell parameters will be ignored.
#langevinPiston on
#useFlexibleCell no ;# to prevent the cell from becoming highly extended in one dimension
#useGroupPressure yes ;# needed for rigid bonds
#useConstantArea no ;# no for water box, maybe for membrane
#langevinPistonTarget 1.01325 ;# pressure in bar -> 1 atm
#langevinPistonPeriod 100. ;# oscillation period around 100 fs------ damping 50 --- period 400  
#langevinPistonDecay 50. ;# oscillation decay time of 50 fs -------decay 200 ---- damping 200
#langevinPistonTemp $temperature ;# coupled to heat bath

temperature 0 ;  # initialize velocities randomly, USE ONLY WHEN STARTING FROM SCRATCH

minimize 180; # AT LEAST 600 worked for this script

output minimized

#Heating

set IHTFRQ 900
set TEMINC 5

for {set i 5} {$i <= $temperature} {incr i $TEMINC } {
langevinTemp $i
reinitvels $i
run $IHTFRQ
}

output equilibrated

#Propagation

langevinTemp $temperature
run $propsteps