[gmx-users] How to use genbox to add ions or other molecules

[chris.neale at utoronto.ca]

genbox isn't too good at ensuring that the added ions/small molecules that it
adds are as far away from the protein as one would like. However, it can be
done. Here's how.

The following procedure will spread a molecule of interest into a pbc box
together with a protein and then solvate the entire system. The first script
echos instructions about how to use cutdown.pl to get rid of any molecules that
are too close to the protein.

STEP 1: make your molecule of interest (1 copy) and put it in ${MYTOP}/myfrag.gro

For example, here is an oxygen molecule with a bond length according to Fischer
and Lago:
GROtesk MACabre and Sinister
    2
    1FRAG   OOA    1   0.000   0.051   0.051
    1FRAG   OOB    2   0.102   0.051   0.051
   0.10200   0.10200   0.10200

* This script assumes myfrag.gro contains a single residue of type FRAG
  You can change that to fit your needs

STEP 2: Run the following script (I call it ./run1a)

#!/bin/sh
MOL=mymol
ED=/mydir/gromacs-3.3/exec/bin
MYTOP=/mydir/gromacs-3.3/exec/share/gromacs/top

${ED}/pdb2gmx -f ${MOL}.pdb -o ${MOL}_g.gro -p ${MOL}_g.top -ff oplsaa -n
${MOL}_g.ndx -q ${MOL}_g.pdb -water tip4p >& output.g_pdb2gmx

${ED}/editconf -c -f ${MOL}_g.gro -o ${MOL}_ge.gro -box 5.7 -bt dodecahedron >&
output.ge_editconf

${ED}/genbox -cp ${MOL}_ge.gro -ci ${MYTOP}/myfrag.gro -nmol 100 -seed ${SEED}
-o ${MOL}_geo.gro -p ${MOL}_g.top >& output.geo_genbox

${ED}/pdb2gmx -f ${MOL}_geo.gro -o ${MOL}_geo.gro -p ${MOL}_geo.top -ff oplsaa
-n ${MOL}_geo.ndx -q ${MOL}_geo.pdb -water tip4p >& output.geo_pdb2gmx

${ED}/grompp -f em.mdp -c ${MOL}_geo.gro -p ${MOL}_geo.top -o ${MOL}_geo.tpr >&
output.geo_grompp

echo "(1) g_mindist -f ${MOL}_geo.gro -s ${MOL}_geo.tpr -or mindistres.xvg"
echo "     - choose (12) FRAG then (1) Protein"
echo "(2) grep -v \@ mindistres.xvg | grep -v \# | sort -nr -k 2"
echo "(3) tail until finding the correct number (0.5 for 5A from protein)"
echo "(4) grep -v \@ mindistres.xvg | grep -v \# | sort -nr -k 2 | tail -X | awk
'{print $1}' | sort -n > list_to_remove.txt"
echo "     - then cp ${MOL}_geo.gro ${MOL}_geo_mod.gro"
echo "     - then use cutdown.pl to remove these FRAG residues from the
${MOL}_geo_mod.gro file"
echo "     - the residue number MUST BE INCREMENTED by the number or protein
residues (edit cutdown.pl before running it)"
echo "     - the 2nd line of the .gro file must be modified to represent the new
# of atoms (count it out yourself and change the number)"
echo "(5) run ./run1b "

###############################


STEP 3: Follow the instructions output at the end or the ./run1a script.
Here is cutdown.pl

#!/bin/bash
MOL=mymol
PROTRES=<insert number of protein residues here>
cp ${MOL}_geo_mod.gro temp.a
for i in $( cat ./list_to_remove.txt ); do
  j=$(expr $i + $PROTRES)
  grep -v " "$j"FRAG " temp.a > temp.b
  rm -f temp.a
  mv temp.b temp.a
done

###############################

STEP 3b: Check that temp.a differs from ${MOL}_geo_mod.gro only by having lost
those fragment inserts that were too close to the protein.

STEP 3c: Remember to change the number of residues at the top of temp.a

STEP 3d: mv temp.a ${MOL}_geo_mod.gro

STEP 4: Run the following script (I call it ./run1b)
#!/bin/sh
MOL=mymol
ED=/mydir/gromacs-3.3/exec/bin
MYTOP=/mydir/gromacs-3.3/exec/share/gromacs/top
SEED=12345

${ED}/pdb2gmx -f ${MOL}_geo_mod.gro -o ${MOL}_geo.gro -p ${MOL}_geo.top -ff
oplsaa -n ${MOL}_geo.ndx -q ${MOL}_geo.pdb -water tip4p >& output.geo_pdb2gmx

${ED}/genbox -cp ${MOL}_geo.gro -cs ${MYTOP}/tip4p.gro -o ${MOL}_geos.gro -seed
${SEED} -p ${MOL}_geo.top >& output.geos_genbox

###############################

Chris Neale.