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

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

structure       ./dimer-kffe-wb.psf
coordinates     ./dimer-kffe-wb.pdb
set temperature    800
set outputname     run1-kffe

firsttimestep      0


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

# Input
paraTypeCharmm	    on
parameters          ./par_all27_prot_lipid.inp 
temperature         $temperature


# Force-Field Parameters
exclude             scaled1-4
1-4scaling           1.0
cutoff               8.
switching           on
switchdist           7.
pairlistdist         9.5


# Integrator Parameters
timestep            2.0
rigidBonds          all
nonbondedFreq       1
fullElectFrequency  2  
stepspercycle       10


# Constant Temperature Control
langevin             on
langevinDamping      1.0
langevinTemp         $temperature
langevinHydrogen     no


# PBC
cellBasisVector1    40.0    0.    0.
cellBasisVector2    0.   40.0     0.
cellBasisVector3    0.    0.    40.0
cellOrigin          0.0 0.0 0.0

wrapAll             on

# PME (for full-system periodic electrostatics)
PME                 yes
PMEGridSizeX        48
PMEGridSizeY        48
PMEGridSizeZ        48


# Output
outputName          $outputname
restartfreq           1000
dcdfreq               1000
outputEnergies         1000
outputPressure         1000


tclBC on
tclBCScript {

  proc veclen2 {v1} {
    foreach {x1 y1 z1} $v1 { break }
    return [expr $x1*$x1 + $y1*$y1 + $z1*$z1]
  }

    proc calcforces {step unique R K} {
    if {$step % 20 == 0} {
    cleardrops
    }

    set R [expr 1.*$R]
    set R2 [expr $R*$R]
    set tol 2.0
    set r2 0.0
    #KFFE
    set massimo 160
    #KFFE
    set cut2 [expr ($R-$tol)*($R-$tol)]

    while {[nextatom]} {
      set atomix [getid]

      if { $atomix > $massimo } {
        dropatom
        continue
      }
      
      if { $atomix <= $massimo } {
      set rvec [getcoord]
      set r2 [veclen2 $rvec]
      }
      
      if { $r2 < $cut2 } {
        dropatom
        continue
      }

      if { ($r2 > $R2) && ($atomix <= $massimo)} {
        set r [expr sqrt($r2)]
        addenergy [expr $K*($r - $R)*($r - $R)]
        addforce [vecscale $rvec [expr -2.*$K*($r-$R)/$r]]
      }
    }
  }
}

tclBCArgs {20.0 2.0}

minimize 1000
reinitvels $temperature
run 1000000