# add counter ions
#
# author: mengen
#
# usage: ./modar step3_add_ions.inp anion=POT cation=CLA ionc=0.1 method=hybrid
# anion default is K+
# cation default is CL-
# concentration default is 0.1 mol/ml
#
#load force field
include "loadff.inc"
#anion and cation name
if(! -var anton) anion="POT"
if(! -var cation) cation="CLA"
if(! -var method) {
method="hybrid"
echo -hs "no adding ions method specified"
echo -hs "the default method hydrid will be used"
echo -hs "or please make a choice bellow"
echo -hs " 1 MC"
echo -hs " 2 best size"
echo -hs " 3 hybrid"
echo -hsn " which one:"
ichoice=gettext()
if($ichoice==1) method="mc"
if($ichoice==2) method="bs"
}
echo -hs " using method: $method"
#concentration of ions, default value is 0.1 mol/ml
if(! -var ionc) {
ionc=0.10
echo -hs "| no ionc (ion concentration) specified, "
echo -hs " default value 0.1 mol/ml will be used"
echo -hsn " or tell the concentration here:"
newionc=gettext()
if($newionc>0.001) ionc=$newionc
}
echo -hs "| ion concentration is $ionc mol/ml"
if($ionc>100 || $ionc<0) {
echo -hs "| ** invalid ion concentration $ionc mol/ml"
}
#show arguments
echo -hs "| anion is $anion"
echo -hs "| cation is $cation"
#load solute
loadmsf file="step1_single.msf"
#load PBC crystal setup
include "pbc.inc"
tell msf_info
#calc ions accessible volume
tell boxsize
tell volume gridsize=0.5 rprobe=1.4 select="all"
accvol=$boxvol-$molvol
#compute number of ions to be added
nion=$ionc*6.021*0.0001*$accvol
npos= int($nion-$netcharge/2) + 1
nneg=$npos+$netcharge
#show number of ions to add
echo -hs "| number of $anion need to add: $npos"
echo -hs "| number of $cation need to add: $nneg"
#buildionacc gridsize=2.0 ion=CLA outpdbfile=iongrids.pdb
#nonbond setup
nonbond type=pme nblcutoff=12.0 nbcutoff=10.0 swcutoff=8.0 eps=80.0 \
beta=0.34 ftx=$fftx fty=$ffty ftz=$fftz bsorder=6
savepdb fmt=pdb file="step3_with_ions.pdb"
#build mdcore
build mdcore timestep=0.002 nstep=100 nprint=1 \
trjfile="test.dcd" frqtrj=100 trjfmt=charmm \
rstfile="test.rst" frqrst=100 \
randseed=314159 randtype=old \
frqresort=20000
nnadd=0
npadd=0
#init system random seed for MC move try
randseed iseed=20000 sys=true
echo -hs "| adding ions ..."
echo -hs "| it may take a few minitues, please wait a while"
nptot=$npos
nntot=$nneg
#adding ions
if($method=="hybrid") goto hybrid
else if($method=="bs") goto bestsite
:mconly
echo -hs "| adding ions using MC method"
while($npos>0 || $nneg>0) {
baddneg=1
if($npos>$nneg) baddneg=0
else if($npos<$nneg) baddneg=1
else {
#calc the voltage for all accessible grids to determin to add positive or negative first
addion testonly=true ion=$cation gridsize=1.0
if(abs($vmin)>$vmax) baddneg=0
}
if($baddneg==1) {
nnadd=$nnadd+1
echo -hs "| adding a $cation, $nnadd/$nntot"
addion method=MC ion=$cation nMCtry=100 gridsize=1.0
nneg=$nneg-1
}
else {
npadd=$npadd+1
echo -hs "| adding a $anion, $npadd/$nptot"
addion method=MC ion=$anion nMCtry=100 gridsize=1.0
npos=$npos-1
}
}
goto doneaddions
:bestsite
echo -hs "| adding ions using best site method"
while($npos>0 || $nneg>0) {
baddneg=1
addion testonly=true ion=$cation gridsize=1.0
echo -hs "| vmin,vmax,vavg,vsig: $vmin, $vmax, $vavg,$vsig"
if($npos>$nneg) baddneg=0
else if($npos<$nneg) baddneg=1
else{
#calc the voltage for all accessible grids to determin to add positive or negative first
if(abs($vmin)>$vmax) baddneg=0
}
if($baddneg==1) {
nnadd=$nnadd+1
echo -hs "| adding a $cation, $nnadd/$nntot"
addion method=bs ion=$cation nMCtry=100 gridsize=1.0
nneg=$nneg-1
}
else{
npadd=$npadd+1
echo -hs "| adding a $anion, $npadd/$nptot"
addion method=bs ion=$anion nMCtry=100 gridsize=1.0
npos=$npos-1
}
}
goto doneaddions
:hybrid
echo -hs "| adding ions using hybrid method"
while($npos>0 || $nneg>0){
baddneg=1
addion testonly=true ion=$cation gridsize=1.0
echo -hs "| vmin,vmax,vavg,vsig: $vmin, $vmax, $vavg,$vsig"
method=mc
if($npos>$nneg) baddneg=0
else if($npos<$nneg) baddneg=1
else {
#calc the voltage for all accessible grids to determin to add positive or negative first
if(abs($vmin)>$vmax) baddneg=0
}
if($baddneg==0) {
if($vmax-$vavg>0.018) method=bs
}
else{
if($vavg-$vmin>0.018) method=bs
}
if($baddneg==1) {
nnadd=$nnadd+1
echo -hs "| adding a $cation, $nnadd/$nntot using $method method"
addion method=$method ion=$cation nMCtry=100 gridsize=1.0
nneg=$nneg-1
}
else{
npadd=$npadd+1
echo -hs "| adding a $anion, $npadd/$nptot using $method method"
addion method=$method ion=$anion nMCtry=100 gridsize=1.0
npos=$npos-1
}
}
:doneaddions
#save solute with ions
savemsf file="step3_with_ions.msf"
savepdb fmt=pdb file="step3_with_ions.pdb"
echo -hs "| output files:"
echo -hs "| step3_with_ions.msf"
echo -hs "| step3_with_ions.pdb"
echo -hs "| "
echo -hs "| it's ready to go next step to put all stuff together:"
echo -hs "| ./modar -nobak step4_merge_all.inp"
|
