#load force field

#include "loadff.inc"

#load md system

loadmsf file="step7_two_phases.msf"

#load MD core and setup

loadsetup file="step9_meshing.mds"

#feed in trj

trjin file="step9_meshing.trj"

#remove files if already exists

if(-f zrgy.txt) rm zrgy.txt

if(-f aout.txt) rm aout.txt

if(-f cout.txt) rm cout.txt

#define anions, cations

defcls anions  select="resn PF*"

defcls cations select="resn BM*"

#tell # of anions, cations

tell nres select="anions"

nresa=$nres

tell nres select="cations"

nresc=$nres

ifrm=0

while($ifrm<$ntotframes){

  trjreadnextframe

  mdtime=$curmdtime

  if($mdtime<0) goto done #double check whether trjreadnextframe successfully or not

  #print progress for every 2ps

  tmod=$mdtime%2

  if($tmod==0) echo -hs "current MD time: $mdtime ps"

  #calculate RMSD for all heavy atoms comparing to initial structure

  calc zrgy select="segn ILC"

  echo "$mdtime $value" >> zrgy.txt

  #report # of anions moved into the water already

  naout=0

  ires=0

  while($ires<$nresa){

    tell center select="ires $ires in anions"

    if(abs($zce)>24) naout++

    ires++

  }

  echo "$mdtime $naout" >> aout.txt

  #report # of cations moved into the water already

  ncout=0

  ires=0

  while($ires<$nresc){

    tell center select="ires $ires in cations"

    if(abs($zce)>24) ncout++

    ires++

  }

  echo "$mdtime $ncout" >> cout.txt

  ifrm++

}

stop

:done

echo -hs "something went wrong when reading trjactory data"

abort