[gmx-users] Indexing problem when using genconf

[Justin A. Lemkul]

Justin A. Lemkul wrote:
> 
> 
> 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"
>>   4   5 [ system ]
>>   6 CHOL
>>   7   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[29952] 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.
> 

...or just do this all in one step with tpbconv.  I always manage to quickly 
think of the hardest way to do things :)

-Justin

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================