#############################################################
## MINIMIZACION: Prueba 1L2Y                               ##
#############################################################

# Archivo de parámetors para realizar dinámica molecular con namd
# NOTE: No incluir acentos

#############################################################
## PARÁMETROS AJUSTABLES                                   ##
#############################################################

## Declara el nombre del archivo .pdb y psf (sin extenciones)
set inputname      tc5b_wb
## Lee los archivos de coordenadas (pdb) y topologia (psf)
structure          $inputname.psf
coordinates        $inputname.pdb
## set restart        min_eq

## Nombre de los archivos de salida
outputName         ${inputname}_min
set temperature    0

# Continua una simulacion desde un archivo restart 0 = NO, 1 = SI
# En este caso no usamos archivos restart
if {0} {
binCoordinates     $restart.restart.coor
# Si se incluyen velocidades se debe remover temperatura
# Pues esta esta ya dada por la velocidad de los atomos del archivo restart
binVelocities      $restart.restart.vel
extendedSystem     $restart.restart.xsc
} 

binaryoutput       off
firsttimestep      0

#############################################################
## PARAMETROS DE SIMULACION                                ##
#############################################################

# Archivos de parametros de Charmm
paraTypeCharmm      on
# Dirección del archivo de parametros
parameters          ./par_all27_prot_lipid.prm ;

## NOTE: No usar temperatura para velocidades iniciales
## si se ha especificado un archivo .vel!
## Considera que al minimizar, este parametro no es usado
temperature         $temperature
 
## Condiciones periódicas de frontera
## NOTE: No especificar cell basis si tambien se especifica 
## un archivo .xsc!
if {1} {
cellBasisVector1    61.51300048828125   0.0 0.0
cellBasisVector2    0.0 63.739999771118164  0.0
cellBasisVector3    0.0 0.0 51.224998474121094
cellOrigin  19.30879020690918   25.521724700927734  5.089336395263672
}

wrapWater           on
wrapAll             on

## Parámetros de campo de fuerza
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.
switching           on
switchdist          10.
pairlistdist        13.5

## Parámetros de integracion
timestep            1.0    # t de integracion en fs
rigidBonds          water  # all/ none /water
stepspercycle       10  # pasos por ciclo de eval.
nonbondedFreq       1   # Freq de eval. para pares no enlazados
fullElectFrequency  2   # Freq de eval. para inter. elec.

## PME (for full-system periodic electrostatics)
## El tamaño del grid determina parcilamente la exactitud 
## y eficiencia del PME. Por velocidad, PMEGridSize_ deberia 
## tener dimensiones cuyos factores sean enteros pequeños (2, 3 and 5).
if {1} {
PME                 yes
PMEGridSizeX        64 ;# 2**6
PMEGridSizeY        64 ;# 2**6
PMEGridSizeZ        64 ;# 2**6
}


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

# Minimization de 10000 pasos
minimize            10000
exit
#############################################################
## CALENTAMIENTO NVT: Prueba 1L2Y                          ##
#############################################################

## Escribe aquí los comentarios necesarios para describir el trabajo

#############################################################
## PARÁMETROS AJUSTABLES                                   ##
#############################################################

## Declara el nombre del archivo .pdb y psf (sin extenciones)
set inputname      tc5b_wb
## Lee los archivos de coordenadas (pdb) y topologia (psf)
structure          $inputname.psf
coordinates        $inputname.pdb
## Ahora si se utiliza un archivo restart, el de la minimizacion
set restart        ${inputname}_min

## Nombre de los archivos de salida
outputName         ${inputname}_heat
set temperature    0

# Continua una simulacion desde un archivo restart 0 = NO, 1 = SI
if {1} {
binCoordinates     $restart.restart.coor
## Si se incluyen velocidades se debe remover temperatura
## Pues esta esta ya dada por la velocidad de los atomos del archivo restart
## No se utilizan velocidades pues el restart de la minimizacion no las genera
## Pero si usamos el archivo .xsc con los datos de las dimensiones de la caja
# binVelocities      $restart.restart.vel
extendedSystem     $restart.restart.xsc
} 

binaryoutput       off
firsttimestep      0

#############################################################
## PARAMETROS DE SIMULACION                                ##
#############################################################

# Archivos de parametros de Charmm
paraTypeCharmm      on
parameters          ./par_all27_prot_lipid.prm

## NOTE: No usar temperatura para velocidades iniciales
## si se ha especificado un archivo .vel!
temperature         $temperature
 
## Condiciones periódicas de frontera
## NOTE: No especificar cell basis si tambien se especifica 
## un archivo .xsc! En este caso ponemos "0" pues el archivo 
## restart contiene también las dim del sistema
if {0} { 
cellBasisVector1    61.51300048828125   0.0 0.0
cellBasisVector2    0.0 63.739999771118164  0.0
cellBasisVector3    0.0 0.0 51.224998474121094
cellOrigin  19.30879020690918   25.521724700927734  5.089336395263672
}

wrapWater           on
wrapAll             on

## Parámetros de campo de fuerza
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.
switching           on
switchdist          10.
pairlistdist        13.5

## Parámetros de integracion
timestep            1.0    # t de integracion en fs
rigidBonds          water  # all/ none /water
stepspercycle       10  # pasos por ciclo de eval.
nonbondedFreq       1   # Freq de eval. para pares no enlazados
fullElectFrequency  2   # Freq de eval. para inter. elec.

## PME (for full-system periodic electrostatics)
## El tamaño del grid determina parcilamente la exactitud 
## y eficiencia del PME. Por velocidad, PMEGridSize_ deberia 
## tener dimensiones cuyos factores sean enteros pequeños (2, 3 and 5).
if {1} {
PME                 yes
PMEGridSizeX        64 ;# 2**6
PMEGridSizeY        64 ;# 2**6
PMEGridSizeZ        64 ;# 2**6
}

# Termostato
langevin            on    ;# do langevin dynamics
langevinDamping     5     ;# damping coefficient (gamma) of 5/ps
langevinTemp        $temperature
langevinHydrogen    no    ;# don't couple langevin bath to hydrogens

## Barsostato: Constant Pressure Control (variable volume)
## Si estos parametros se ejecutan ( = "1") se realizara una simulacion NPT
## De otro modo, si = "0", se realiza una simulacion NVT
# Para SA usaremos NVT
if {0} {
useGroupPressure      yes ;# needed for 2fs steps
useFlexibleCell       no  ;# 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  100.0
langevinPistonDecay   50.0
langevinPistonTemp    $temperature
}

## Intervalo de guardado de las coordenadas y archivos
restartfreq         10000     ;# 500steps = every 1ps
dcdfreq             10000
xstFreq             10000
outputEnergies      10000
outputPressure      10000

#############################################################
## SCRIPT DE EJECUCION: SIMULATED ANNEALING                ##
#############################################################

##   Asigna velocidades a la temperatura dada
reinitvels          $temperature

## Ejecuta por 10000 pasos = 10 ps
## run 10000

##  temperatura mínima de inicio
set temp_min 0
##  tempratura máxima
set temp_max 400
##  temperatura final de la DM
set temp_final 310

## ------------------------
##  PARTE 1: Incrementar en 50 intervalos de 0 a 400K durante 100 ps
### Número  total de pasos, intervalos, incremento y pasos por intervalo para dcada ciclo
set stps 100000
set intv 50
set incremento_temp [expr ($temp_max - $temp_min)/$intv]
set stps_by_intv [expr $stps/$intv]

for { set temperature $temp_min } { $temperature <= $temp_max } { incr temperature $incremento_temp } {
    
    reassignTemp $temperature
    langevinTemp $temperature
    run $stps_by_intv
}

# -------------------------
#   PARTE 2: Mantiene la temperatura a 400K por 100 ps
run 100000

# -------------------------
#   PARTE 3: Disminuye en 50 intervalos la Temp de 400K a 310K en 0.2 ns
#   temperature == 400
set temp_final 310

set stps 200000
set intv 50
set decremento_temp [expr ($temp_final - $temp_max)/$intv]
set stps_by_intv [expr $stps/$intv]

for { set temperature $temp_max } { $temperature >= $temp_final } { incr temperature $decremento_temp } {
    
    reassignTemp $temperature
    langevinTemp $temperature
    run $stps_by_intv
}

# -------------------------
#   PARTE 4: Templado - > Mantiene a 310K por 600 000 fs por 0.6 ns
run 600000

# -----
# Tiempo total de Simulated annealing : 1 ns