[gmx-users] Indexing problem when using genconf
Justin A. Lemkul
jalemkul at vt.edu
Tue May 31 05:22:45 CEST 2011
Ryan S Davis (rsdavis1) wrote:
> I wanted to copy a bilayer into a grid of 2x2x1 replicas. I used genconf and everything seemed to work fine exept that annoying feature
> that the command does not reorder the molecule types, so I end up with a .top file looking like this...
> 1 #include "martini_v2.1.itp"
> 2 #include "martini_v2.0_lipids.itp"
> 3 #include "martini_v2.0_cholesterol.itp"
> 5 [ system ]
> 6 CHOL
> 8 [ molecules ]
> 9 DPPC 832
> 10 CHOL 208
> 11 W 8320
> 12 DPPC 832
> 13 CHOL 208
> 14 W 8320
> 15 DPPC 832
> 16 CHOL 208
> 17 W 8320
> 18 DPPC 832
> 19 CHOL 208
> 20 W 8320
> Anyway, I run the simulation...no errors. I make an ndx file using make_ndx...indices look fine despite the repetitive order. HOWEVER, when I try to run commands such as
> trjconv with the index file as input, it reads all the way up to the first block of Waters and quits with the error
> Program trjconv, VERSION 4.0.7
> Source code file: gmx_trjconv.c, line: 1037
> Fatal error:
> Index 46593 is larger than the number of atoms in the trajectory file (46592)
> which I didnt expect, but makes perfect sense knowing that I specified in the .mdp file to not output water to the xtc file...
> xtc-grps = dppc chol
> Normally this isnt an issue because waters are typically last in the topology. But, I still need access to this data. How can I force the post-processing commands to read past the absent water blocks?
> The only options I see at the moments is to
> 1) scrap genconf, make new topology somehow, and rerun
> 2) reset to output water, and rerun
> 3) limit my analysis to the very sparse output from the .trr file
I see two viable options, one of which is to use your option 1 via a few
standard Unix tools to create a proper coordinate file, i.e.:
grep DPPC conf.gro > dppc
grep CHOL conf.gro > chol
grep W conf.gro > w
cat dppc chol w > system.gro
Add a title, number of atoms, and box vectors at the end, and you're done. Then
sum up the entries in [molecules] and the system should run just fine.
The second option is to just create a dummy system that has only DPPC and CHOL
molecules. Take your input coordinate file, strip out the waters, and renumber
with genconf. Remove all W blocks from the .top and create a new .tpr file.
This should now match the number of atoms in the trajectory and the order in
which they appear.
The latter option seems to be much more efficient and will not require you to
re-run any simulations, but I thought it would be worth mentioning the first
method, since I use it all the time in my heterogeneous membranes.
Justin A. Lemkul
ICTAS Doctoral Scholar
Department of Biochemistry
jalemkul[at]vt.edu | (540) 231-9080
More information about the gromacs.org_gmx-users