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

# Minimization of a water cube on graphene monolayer. Water cube is 1000 atoms. Graphene is 5nmx5nm


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

structure          ./pore(5,2).psf
coordinates	   ./pore(5,2)_m.coor

set itemp	    0
set ftemp	    300
set outputname      pore(5,2)_i

firsttimestep      0


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

# Input
paraTypeCharmm	    on
parameters          ./graphene_par_all27_prot_na.prm
temperature         $itemp


# Force-Field Parameters
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.0
switching           on
switchdist          10.0
pairlistdist        14.0

# Heating parameters
reassignFreq 1 # Number of steps between reassignment velocities
reassignIncr 0.001 # Increment used to adjust temperature during temperature reassignment
reassignHold $ftemp # The value of temperature to be kept after heating is completed 

langevin on
langevinTemp $ftemp
# LangevinPiston on
# LangevinPistonTarget 1.01325
# LangevinPistonPeriod 200
# LangevinPistonDecay 100
# LangevinPistonTemp $ftemp
# useFlexibleCell off

# Integrator Parameters
timestep            1.0  ;
rigidBonds          none  ; 
nonbondedFreq       1
fullElectFrequency  1  
stepspercycle       5

# Periodic Boundary Conditions
# minmax results {-2.7060000896453857 -1.8209999799728394 -0.27799999713897705} {50.89799880981445 52.16999816894531 0.2549999952316284}
#center 24.105510711669922 25.098291397094727 0.0010484461672604084
cellBasisVector1     105.0820088119507    0    0.0
cellBasisVector2     0.0  107.8139978942871   0.0
cellBasisVector3     0.0    0  62.39100137329102
cellOrigin           44.24100041389465 5.218999862670898 -0.2154998779296875

wrapAll             on


# PME (for full-system periodic electrostatics)
PME                 yes
PMEGridSpacing      1.228

#manual grid definition
PMEGridSizeX        108
PMEGridSizeY        108
PMEGridSizeZ        64


# Output
outputName          $outputname

restartfreq         500     ;# 500steps = every 1ps
dcdfreq             250
xstFreq             250
outputEnergies      100

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

binaryrestart no
binaryoutput no

fixedAtoms on
fixedAtomsFile ./pore(5,2).pdb
fixedAtomsCol B

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

tclForces		on

tclForcesScript {
  set cellLengthZ 58.80699983215332 ;# The length of the unit cell in the z-direction, should be the same as the value in "cellBasisVector3"
  set LowerBoundary -12.5
  set UpperBoundary 12.5
  set force {0 0 0.1}
  set watIdList {}
  for {set i 577} {$i<1320} {incr i 3} {
    lappend watIdList $i
    addatom $i
  }
  proc calcforces {} {
    global cellLengthZ LowerBoundary UpperBoundary force watIdList
    loadcoords coorList
    foreach i $watIdList {
      set z [lindex $coorList($i) 2]
      set z [expr $z-round($z/$cellLengthZ)*$cellLengthZ]  ;# Translate this coordinate into its corresponding value in the original unit cell
      if {$z>$UpperBoundary || $z<$LowerBoundary} {
        addforce $i $force
      }
    }
  }
}


run 300000