[gmx-users] How to recenter solvent around solute

[Dimitar Pachov]

Hi,

On Sat, May 7, 2011 at 9:45 PM, Justin A. Lemkul <jalemkul at vt.edu> wrote:


> Please clarify - do you wish to maintain the original triclinic
> representation (as -ur tric does) or do you wish to see the octahedral
> representation (as -pbc mol -ur compact gives)?  My answer was based on your
> request to "keep the original unit shape."
>
>
My original representation is truncated octahedron. Gromacs requires
triclinic vectors and does not care about unit cell shape at the moment the
md code starts. I wanted to keep MY original unit shape.


>>
> For complicated systems (or sometimes even for simple ones), it is quite
> common that multiple iterations of trjconv are necessary.  A suggested
> workflow is indeed documented:
>
>
> http://www.gromacs.org/Documentation/Terminology/Periodic_Boundary_Conditions#Suggested_trjconv_workflow
>
> The documentation must remain somewhat generic, as there are a number of
> specialty systems that can be considered.
>

I thought the guidelines given by the trjconv specific manual were
sufficient. I don't understand why people don't put this suggested workflow
in the trjconv manual. It is confusing otherwise. It never came to me that
one needed multiple usages of trjconv in order to achieve recentering around
solute.

I followed these instructions and got what I wanted. Basically, this was my
quick workflow:

=================================
psfinname="$prot-$shape-$flc"
topfile="$prot-$shape-$flc.top"
grofile="$shape-$fl-eq.gro"
inpfile="md.mdp"

name="run1"
newname="run1-recen"
indexfile="index-all.ndx"
tclfile="extr-frame-traj.tcl"

s="Protein-Nucleo"
indsolute=`cat $indexfile | awk -v site=$s '$1=="["{i++}$2==site{print
i-1}'`
cp $name.xtc $newname.xtc

#===> 1
intraj="$newname.xtc"
outtraj="traj/$newname-whole.xtc"
trjconv -f $intraj -s $name.tpr -n $indexfile -o $outtraj  -pbc whole > log1
2>&1 <<EOF1
0
EOF1
#===> 2
intraj="$outtraj"
outtraj="traj/$newname-cluster.xtc"
trjconv -f $intraj -s $name.tpr -n $indexfile -o $outtraj  -pbc cluster >
log2 2>&1 <<EOF2
$indsolute
0
EOF2
#===> 3
intraj="$outtraj"
outtraj="traj/$newname-nojump.xtc"
trjconv -f $intraj -n $indexfile -o $outtraj  -pbc nojump > log3 2>&1 <<EOF3
0
EOF3
#===> 4
intraj="$outtraj"
outtraj="traj/$newname-center.xtc"
trjconv -f $intraj -n $indexfile -o $outtraj  -center > log4 2>&1 <<EOF4
$indsolute
0
EOF4
#===> 5
## Get first frame from centered .xtc
fr=1
pdbfilename="`echo $outtraj | awk 'BEGIN{FS="."}{print $1}'`"
pdbfile="$pdbfilename-${fr}ps.pdb"
`which vmd` -dispdev text -e $tclfile -args $fr
${pathpsfin}/${psfinname}.psf $outtraj $pdbfilename > /dev/null

cat $pdbfile | sed 's/TIP3/SOL /' | sed 's/ OH2 SOL/ OW  SOL/' | sed 's/ H1
 SOL/ HW1 SOL/' | sed 's/ H2  SOL/ HW2 SOL/' | sed 's/ SOD SOD/ NA  NA /g' |
sed 's/ CLA CLA/ CL  CL /g' | sed 's/ HT1 MET/ H1  MET/' | sed 's/ HT2 MET/
H2  MET/' | sed 's/ HT3 MET/ H3  MET/' > $pdbfilename-${fr}ps-renamed.pdb

grompp -f $inpfile -c $pdbfilename-${fr}ps-renamed.pdb -n $indexfile -p
$topfile -o $newname-center.tpr > log5 2>&1
#===> 6
intraj="$outtraj"
outtraj="traj/$newname-compact.xtc"
trjconv -s $newname-center.tpr -f $intraj -n $indexfile -o $outtraj  -pbc
mol -ur compact  > log6 2>&1 <<EOF6
0
EOF6
=================================


However, it seems a long process, and if I have a large number of
trajectories, it might not be feasible. I actually do not know if all steps
are necessary. Particularly, I did step #5 just to get a new *tpr file from
the centered trajectory ... Is there a way to do that in a simpler way
within Gromacs?




On Sun, May 8, 2011 at 4:06 AM, Tsjerk Wassenaar <tsjerkw at gmail.com> wrote:

> Hi,
>
> Justin is very right to point at the workflow. In your case, the
> points 5 and 6 apply: centering, and putting things in some box. Maybe
> that doesn't work in one pass and you'll have to use trjconv twice.
>
> >> This does not work. Neither is the solute always surrounded by solvent
> for
> >> every frame (it goes out of the box) nor is the original unit cell
> >> (truncated octahedron) preserved.
>
> You have to think out of the box with PBC :) The solute is always
> surrounded by solvent, except where it touches its own image, in which
> case you should consider throwing away and redo the run in a larger
> box. There is no out of the box where atoms can go to. Furthermore,
> the unit cell does not have shape, it only has translational relations
> (lattice vectors).


It does not have a shape in Gromacs; otherwise it has.



> You can carve out any general shape you want, such
> as rectangular, truncated octahedron, rhombic dodecahedron. Note that
> only in special cases you can carve out regular instances of those
> (such as a cube for a rectangular brick).
> To understand this, here's a bit of an exercise:
>
> Draw three identical hexagons in contact with each other
> >From one of the hexagons, draw vectors to the other two, and draw the
> parallellogram spanned by these vectors.
> Place a hexagon at the farthest node, and not how the parallellogram
> has exactly the same content as any of the hexagons.
> Now turn the parallellogram into a rectangular brick by translating a
> part of it, and note how the brick still has the same contents as a
> hexagon.
> Feeling confident that any shape will do, find some tilings from
> Escher and try to see how you can fit that into a parallellogram,
> rectangle, or hexagon :) Of course the coolest would be to use those
> to define a unit cell to put your atoms in :p
>

I did not quite get why I needed all of this advice, but thanks for your
input anyway.

Best,
Dimitar




>
> Groetjes,
>
> Tsjerk
>
>
> --
> Tsjerk A. Wassenaar, Ph.D.
>
> post-doctoral researcher
> Molecular Dynamics Group
> * Groningen Institute for Biomolecular Research and Biotechnology
> * Zernike Institute for Advanced Materials
> University of Groningen
> The Netherlands
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110509/456c1925/attachment.html>