Re: [mdnalysis-discussion] Applying PBC to trajectory

[Oliver Beckstein]

Hi Andreas,

On 10 Sep, 2012, at 05:57, Andreas Rønnest wrote:

I simulate a DMPG bilayer in a water solution. For unknown reasons the system has a drift in the z axis, meaning the bilayer moves in the z axis (up/down). I use NAMD.

When using MDAnalysis I apply PBC manually by subtracting the center of mass of the DMPG bilayer from all atoms and then move the atoms sticking out of the new box, to the other end i.e. wrapping, by an if statement. This work for the most part, but if it was possible to write the corrected positions back into the universe.trajectory, it would help me immensely.

You can set all atom coordinates of an AtomGroup with the set_positions() [1] method or directly assign to AtomGroup.positions [2] – both methods are equivalent and it's pretty much a matter of taste. (If you're pressed for speed you can also directly modify the underlying Timestep._pos array [3] although the previous two methods are really recommended.)

(There are even more options, come to think: AtomGroup.translate() [4] should also work.)

The downside is that once you read the next timestep, the corrected positions are gone – but see below.

I hope my question is understandable, in short: "Is it possible to write new positions to the universe.trajectory?".

Yes: in oder to do something like

for ts in u.trajectory:

t = calculate_displacements(u)

# shift atoms

u.atoms.positions += t

# now analyze the shifted coordinates frame

analyze_coordinates(u.atoms.positions)

If you want to write a whole new trajectory, use a trajectory writer [5]  and do something like

W = MDAnalysis.Writer("remapped.dcd", numatoms=u.atoms.numberOfAtoms())

for ts in u.trajectory:

t = calculate_displacements(u)

u.atoms.positions += t

W.write(ts)

W.close()

Hope that helps,

Oliver

[1] http://packages.python.org/MDAnalysis/documentation_pages/core/AtomGroup.html#MDAnalysis.core.AtomGroup.AtomGroup.set_positions

[2] http://packages.python.org/MDAnalysis/documentation_pages/core/AtomGroup.html#MDAnalysis.core.AtomGroup.AtomGroup.positions

[3] http://packages.python.org/MDAnalysis/documentation_pages/coordinates/base.html#MDAnalysis.coordinates.base.Timestep._pos

[4] http://packages.python.org/MDAnalysis/documentation_pages/core/AtomGroup.html#MDAnalysis.core.AtomGroup.AtomGroup.translate

[5] http://packages.python.org/MDAnalysis/documentation_pages/coordinates/core.html#MDAnalysis.coordinates.core.writer

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst  * orbe...@gmail.com

[Andreas Rønnest]

Thanks a lot. I apologize for the late answer. I kinda got lost in this Google Groups thing.

I did as you said and the molecules are moved. Performance wise it's smarter for me to write a new dcd file, instead of fixing each frame and then analyzing. But now I've run into a problem with the dimensions/unitcell. Basically some of the values are moved around when writing the timesteps. The ._unitcell is correct, but dimensions are wrong. I can see other people have had the same problem, but also that it should have been fixed. I checked, and I should be using MDAnalysis 0.7.6.

Do you know what is going wrong?

my code: https://www.dropbox.com/s/3xjabamrr5omu5t/dcdpbc.py

[Oliver Beckstein]

Hi Andreas,

On 10 Oct, 2012, at 05:14, Andreas Rønnest wrote:

>
> I did as you said and the molecules are moved. Performance wise it's smarter for me to write a new dcd file, instead of fixing each frame and then analyzing. But now I've run into a problem with the dimensions/unitcell.
> Basically some of the values are moved around when writing the timesteps. The ._unitcell is correct, but dimensions are wrong.

Can you give an example? From what you're saying I don't understand what's wrong.

Best thing is to file a issue with a small example code snippet that shows the bug. Tell us what you would expect and what is actually seen. See http://code.google.com/p/mdanalysis/wiki/ReportingProblems

Once we have a well defined test case then it's pretty straightforward to find a solution.

> I can see other people have had the same problem,
> but also that it should have been fixed. I checked, and I should be using MDAnalysis 0.7.6.

Is there an Issue or email that you're referring to in particular? I quickly checked the CHANGELOG http://code.google.com/p/mdanalysis/source/browse/package/CHANGELOG?name=develop and don't see any mentioning.

> Do you know what is going wrong?

Not without a more detailed issue report :-).

>
> my code: https://www.dropbox.com/s/3xjabamrr5omu5t/dcdpbc.py

That looks ok to me, you're not doing anything to the unitcell.

Oliver

[Andreas Rønnest]

> Can you give an example? From what you're saying I don't understand what's wrong. 

Yes, of course. With the supplied code I've made this test dcd file. Asking for u.trajectory.ts.dimensions yields this result:

array([ 62.2541275 ,  86.92596436,  90.        ,  61.25470734,

        90.        ,  90.        ], dtype=float32)

Asking for u.trajectory.ts._unitcell yields this result:

array([ 62.2541275 ,  61.25470734,  86.92596436,  90.        ,

        90.        ,  90.        ], dtype=float32)

They are supposed to be the same, and in the original (which is too big to supply) both are equal to the _unitcell array seen above.

test.dcd: https://www.dropbox.com/s/gp5mi3dsgxn2lzn/test.dcd

> Is there an Issue or email that you're referring to in particular? I quickly checked the CHANGELOG http://code.google.com/p/mdanalysis/source/browse/package/CHANGELOG?name=develop and don't see any mentioning. 

I was referring to this one:

https://groups.google.com/forum/?fromgroups=#!topic/mdnalysis-discussion/ezMsY5O1yuE

> Not without a more detailed issue report :-).

What, are you telling me you don't know how to read minds?

> That looks ok to me, you're not doing anything to the unitcell. 

Good.

For now I've worked around the problem with the line:

xbox, zbox, bla, ybox = u.trajectory.ts.dimensions[:4]

But it seems the wrong dimensions, pbc wise, messes with center of mass time series. I'm not certain of this, as I couldn't quite understand whether it was dimensions or _unitcell that was used internally for pbc.

Anyways, thanks again for your help and patience.

Andreas

[Oliver Beckstein]

Hi Andreas,

can you please open an Issue in the tracker – this seems to be a real problem. Can you please also supply a PSF or PDB so that I can load your test dcd? If the files are big and you don't want to attach them to the Issue report then just post the drop box links there.

On 15 Oct, 2012, at 04:52, Andreas Rønnest wrote:

> > Can you give an example? From what you're saying I don't understand what's wrong.
>
> Yes, of course. With the supplied code I've made this test dcd file. Asking for u.trajectory.ts.dimensions yields this result:
> array([ 62.2541275 , 86.92596436, 90. , 61.25470734,
> 90. , 90. ], dtype=float32)
>
> Asking for u.trajectory.ts._unitcell yields this result:
> array([ 62.2541275 , 61.25470734, 86.92596436, 90. ,
> 90. , 90. ], dtype=float32)
>
> They are supposed to be the same, and in the original (which is too big to supply) both are equal to the _unitcell array seen above.
>
> test.dcd: https://www.dropbox.com/s/gp5mi3dsgxn2lzn/test.dcd
>
>
> > Is there an Issue or email that you're referring to in particular? I quickly checked the CHANGELOG http://code.google.com/p/mdanalysis/source/browse/package/CHANGELOG?name=develop and don't see any mentioning.
>
> I was referring to this one:
> https://groups.google.com/forum/?fromgroups=#!topic/mdnalysis-discussion/ezMsY5O1yuE

In the light of what you're saying, it seems I misunderstood the real problem in that thread but as there was no follow-up and no issue report I forgot about this.

Note to everyone in general: If a problem with MDAnalysis bothers you then file an issue and provide the appropriate information – if you don't do this then the developers will assume that it's not that bad and spend their time doing other things – which then might lead to real bugs being overlooked (as might have happened in this case).

>
> > Not without a more detailed issue report :-).
>
> What, are you telling me you don't know how to read minds?

Well, the card you're holding up right now is .....

the Ace of Spades.

Right?

Thanks,
Oliver

>
> > That looks ok to me, you're not doing anything to the unitcell.
> Good.
>
> For now I've worked around the problem with the line:
> xbox, zbox, bla, ybox = u.trajectory.ts.dimensions[:4]
>
> But it seems the wrong dimensions, pbc wise, messes with center of mass time series. I'm not certain of this, as I couldn't quite understand whether it was dimensions or _unitcell that was used internally for pbc.
>
> Anyways, thanks again for your help and patience.
>
> Andreas
>

[Naveen Michaud-Agrawal]

On Mon, Oct 15, 2012 at 7:52 AM, Andreas Rønnest <dan...@gmail.com> wrote:

> Can you give an example? From what you're saying I don't understand what's wrong. 

Yes, of course. With the supplied code I've made this test dcd file. Asking for u.trajectory.ts.dimensions yields this result:

array([ 62.2541275 ,  86.92596436,  90.        ,  61.25470734,

        90.        ,  90.        ], dtype=float32)

Asking for u.trajectory.ts._unitcell yields this result:

array([ 62.2541275 ,  61.25470734,  86.92596436,  90.        ,

        90.        ,  90.        ], dtype=float32)

They are supposed to be the same, and in the original (which is too big to supply) both are equal to the _unitcell array seen above.

Hi Andreas,

It looks like the timestep._unitcell (which should be [A,C,alpha,B,beta,gamma] for DCD files) is actually storing [A, B, C, alpha, beta, gamma]. It looks like the DCD writer isn't creating the correct underlying order for the unitcell array (as defined by the CHARMM docs). How are you generating your test.dcd?

Naveen

[Naveen Michaud-Agrawal]

Oli,

I think I've found the error (although I've been away from the code so long that I might be off):

in DCDWriter.write_next_timestep, the code converts the dimensions of the passed in timestep object to a unitcell by calling base.Writer.convert_dimensions_to_unitcell(), which returns a unitcell structured as [A,B,C,alpha,beta,gamma], which is then passed directly into the low level DCD writing routines which expects [A,alpha, B, beta, gamma, C].

I've attached a patch since I'm not able to get log in to my google code account at the moment.

Naveen

--
-----------------------------------
Naveen Michaud-Agrawal

[Oliver Beckstein]

Naveen,

Thanks, the patch seems to fix the problem that the DCDWriter writes the wrong CHARMM unitcell.

I added a test case (which only passes with your patch) and pushed everything.

Andreas, can you try the latest develop branch from git?

Thanks,
Oliver

p.s.: apologies if you should get two emails with the above contents... mail issues.

On 15 Oct, 2012, at 14:47, Naveen Michaud-Agrawal wrote:

> I think I've found the error (although I've been away from the code so long that I might be off):
>
> in DCDWriter.write_next_timestep, the code converts the dimensions of the passed in timestep object to a unitcell by calling base.Writer.convert_dimensions_to_unitcell(), which returns a unitcell structured as [A,B,C,alpha,beta,gamma], which is then passed directly into the low level DCD writing routines which expects [A,alpha, B, beta, gamma, C].
>
> I've attached a patch since I'm not able to get log in to my google code account at the moment.
>

[Andreas Rønnest]

You guys work fast. I'll try it out tomorrow and report back.

Thank you again for your help.

Btw. I was holding my student card. Close though.

Andreas

[Andreas Rønnest]

Okay, I couldn't wait.

It works. This is the output from u.trajectory.ts.dimensions:

array([ 62.2541275 ,  61.25470734,  86.92596436,  90.        ,

        90.        ,  90.        ], dtype=float32)

Thank you very much.

Andreas

[Oliver Beckstein]

Thanks for the feedback, I'm glad that we eliminated this bug.
Oliver

On 16 Oct, 2012, at 10:25, Andreas Rønnest wrote:

> Okay, I couldn't wait.
>
> It works. This is the output from u.trajectory.ts.dimensions:
>
> array([ 62.2541275 , 61.25470734, 86.92596436, 90. ,
> 90. , 90. ], dtype=float32)
>
> Thank you very much.