[gmx-users] Re: Problem with trjconv and centering bilayer.

[Ioannis Beis]

Quoting gmx-users-request at gromacs.org:


> Message: 2
> Date: Mon, 16 Jan 2012 13:38:11 -0500
> From: "Justin A. Lemkul" <jalemkul at vt.edu>
> Subject: Re: [gmx-users] Problem with trjconv and centering bilayer.
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <4F146E93.4090102 at vt.edu>
> Content-Type: text/plain; charset=UTF-8; format=flowed
>
>
>
> Ioannis Beis wrote:
>> Dear gromacs users,
>>
>> I am trying to center the trajectory of a bilayer in the rectangular
>> simulation box in the frame of my effort to calculate the bilayer
>> thickness with g_dist. According to the visualization, the upper layer
>> of the membrane lies on the lowest part of the box and the lower part of
>> the membrane on the highest part of the box, with the water in the
>> middle. There should not be anything wrong with drifting along the
>> z-axis, since the initial structure was at the very bottom of the box,
>> so even a small drift downwards -due to the pbc- would place the lower
>> part of the membrane on the top of the box.
>>
>> I tried issuing:
>>
>> trjconv -f my_trajectory.xtc -s my_initial_structure.gro -o
>> my_trajectory_no_jump.xtc -pbc nojump
>>
>> with 0 (for system) and subsequently
>>
>> trjconv -f my_trajectory_no_jump.xtc -s DOPC_1DOG.tpr -o
>> my_trajectory_no_jump_center.xtc -center -boxcenter tric -pbc mol
>>
>> with 1 (other) for centering and 0 (system) for output.
>>
>> I used this kind of commands some time ago and they worked fine. Now
>> strangely the bilayer remains around the position that I described in
>> the beginning. Is it possible that the program treats the bilayer
>> separated as it looks in vmd and considers the geometrical center in the
>> middle of the water instead of the middle of the bilayer?
>>
>
> That's exactly what it does.
>
>> I tried the -trans flag to translate the bilayer along z-axis near the
>> center and repeated the steps, but the new trajectory wasn't even
>> visible in vmd. In addition, the .gro files generated by trjconv are
>> apparently trajectory files instead of structure files. I don't feel
>> confident about using editconf for operations like translation of a
>> single initial structure, because as far as I understand editconf is a
>> tool meant mostly for setting up systems; thus I was sceptical about
>> messing a structure with editconf that I would later on use together
>> with an .xtc file as input for trjconv.  However, trjconv doesn't seem
>> to generate structure output files. I don't understand the meaning of a
>> .gro file as trajectory file since there are already at least 3 file
>> formats for trajectories.
>>
>
> Generally speaking, a trajectory is just a series of coordinates, so you can
> output it into a number of formats, including .gro and .pdb, among   
> others.  You
> get a chain of coordinate files that may (or may not) end up being useful in
> some applications.
>
>> So how can I bring my bilayer's trajectory in the center of the unit
>> cell for reliable thickness calculation without drifts? Is it possible
>> at the same time to have clear visualization without disturbances?
>> trjconv with -pbc mol still gives rise to lines in the visualization,
>> apparently as a result of atom jumps. Sadly trajectory files cannot be
>> inspected and visualization is quite handy for certain types of feedback
>> in various data analysis-related tasks, so it would be nice if the
>> trajectories used for analysis also look proper in vmd.
>>
>
> I can think of two approaches, the first of which I have used so it   
> should work ;)
>
> 1. Provide a custom index group specifying only a single lipid atom   
> from the end
> of a hydrocarbon chain and center on it.  Therefore, its geometric   
> center has to
> be the center of the box, and it should bring the rest of the membrane to the
> (visual) center of the unit cell.
>
> 2. Calculate the distance with the trajectory you have now, and   
> subtract it from
> the z-length (assuming the membrane plane is x-y) of the box (stored  
>  in the .edr
> file).

First of all thank you very much for the reply. To make sure there is  
not something that I misunderstand about how trjconv works, I would  
like to mention what I did and how I anticipate it.

I made the index file with a last carbon from a random acyl chain and issued:

trjconv -f my_init_traj.xtc -o towards_center_traj.xtc -pbc mol  
-center -boxcenter tric -n tip_atom.ndx

with centering at the index file atom and then:

trjconv -f towards_center_traj.xtc -s my_binary.tpr -o  
centered_traj.xtc -pbc mol -center -boxcenter tric

with centering at the middle of the bilayer. The final .xtc output is  
to be used as input in the g_dist.

I assume that .tpr is used only for the masses of the atoms (so we are  
talking about a mass-weighted geometrical center, which isn't  
mentioned in the manual). Therefore it is mandatory in the second  
command, but not in the first. Is this going to give the correct result?

A relevant concern is whether the -pbc flag is used properly. I just  
used the mol value assuming that is appropriate for this context (is  
it necessary?). Does it make any difference in the results if I use  
-pbc nojump or no -pbc at all in the first command? To be honest, I do  
not understand well the pbc treatment of trjconv. For example it is  
not clear who is translated where exactly and certain options imply  
modification of the size of the box. Is there a clear description  
somewhere? The manual is confusing. In most cases one uses the  
original trajectory, but it is important to know what does what in pbc  
for certain tasks.

Thank you in advance!

Ioannis


>
> -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
>
> ========================================
>