Failing to implement triclinic_to_orthorhombic and nojump

[sanjiv856]

I am new to this and still trying to figure out both MDAnalysis and LiPyphilic (and MD analysis in general). I ran an MD with a peptide in presence of POPG and POPE, the system was built using CHARMM-GUI, ran MD with Gromacs (2021.3). Doing analysis using MDAnalysis(v2.0.2) and lipyphilic (v0.10.0) on Ubuntu (20.03).

The problem I am facing: After transformations from MDAnalysis (i.e. transformations.unwrap(ag), transformations.center_in_box(protein, center = 'mass'), transformations.wrap(ag, compound='fragments'), transformations.fit_rot_trans(protein, reference, plane='xy', weights="mass"), the system still appears broken with lipids jumping around.

I can not figure out a way to use triclinic_to_orthorhombic and nojump from LipPyphilic together, which I believe would help.

Bringing files in:
u = mda.Universe(working_tpr, working_trr) u.transfer_to_memory(step=10)

Defining variables:
ag = u.atoms # selecting all atoms protein = u.select_atoms('protein') not_protein = u.select_atoms('not protein') ref_u = u.copy() reference = ref_u.select_atoms("protein")

Loading new universe before each transformation trial:
`u = mda.Universe(working_tpr, working_trr)
u.transfer_to_memory(step=10)

workflow = (triclinic_to_orthorhombic(ag), nojump(ag=ag, nojump_x=True, nojump_y=True, nojump_z=False))
u.trajectory.add_transformations(*workflow)`

The error I am getting:
ValueError: nojump requires an orthorhombic box. Please use the on-the-fly transformation :class:lipyphilic.transformations.triclinic_to_orthorhombic before calling nojump
or
ValueError: Can't add transformations again. Please create a new Universe object

Please suggest.

Regards,

Sanjiv Kumar

[p-j-smith]

Hi Sanjiv,

Welcome to the lipyphilic disussions!

The error that you get:

ValueError: Can't add transformations again. Please create a new Universe object

occurs when you try to add two sets of transformations to an MDAnalysis universe. This is an MDAnalysis error that is raised, not a lipyphilic one.

To use the lipyphilic transformations you mentioned, this should work for you:

import MDAnalysis as mda
from lipyphilic.transformations import triclinic_to_orthorhombic, nojump

u = mda.Universe(tpr, trr)
u.transfer_to_memory(step=10)
ag = u.select_atoms("all")
workflow = (
    triclinic_to_orthorhombic(ag),
    nojump(ag=ag, nojump_x=True, nojump_y=True, nojump_z=True),
)
u.trajectory.add_transformations(*workflow)

If you then try to add any other transformations later on (u.trajectory.add_transformations(*workflow)), you will get the error you mentioned.

However, I'm not sure triclinic_to_orthorhombic and no_jump are what you're after. triclinic_to_orthorhombic is only useful if you have a triclinic system, rather than the much more common orthorhombic box. And no_jump is really only useful if you want to look at the diffusion of your proteins or lipids.

If you want to look at e.g. the structure of your protein over time, or specific interactions between your protein and lipids, then your original approach is nearly correct. However, by applying the rotation, you should be aware that you will then not be able to wrap or unwrap your molecules later on, or take into account periodic boundary conditions. This is because there is no (simple) way to apply periodic boundary conditions to a system which has undergone an arbitrary rotation.

To align your trajectory to a centered protein, you'll first need to create a reference structure of your centered protein by following the MDAnalysis userguide:

import MDAnalysis as mda
import MDanalysis.transformations as trans

u = mda.Universe(tpr, trr)
protein = u.select_atoms("protein")
non_protein = u.select_atoms("not protein")

transforms = [
    trans.unwrap(protein),
    trans.center_in_box(protein, wrap=True),
    trans.wrap(non_protein),
]
u.trajectory.add_transformations(*transforms)

protein.write("reference-protein-structure.pdb")

Once you've created this reference structure, you'll need to re-load your original universe (as only a single set of transformations can be applied), then apply the rotation and transformation:

import MDAnalysis as mda
import MDanalysis.transformations as trans

ref_u = mda.Universe("reference-protein-structure.pdb")
u = mda.Universe(tpr, trr)
protein = u.select_atoms("protein")

transforms = [
    trans.unwrap(non_protein),
    trans.fit_rot_trans(protein, ref_u, plane="xy", weights="mass"),
]
u.trajectory.add_transformations(*transforms)

Then the trajectory will be aligned to the centered reference structure, and all molecules will be whole because of the unwrap you applied before the transformation. However, the molecules will not be in the primary unit cell due to the rotation that was applied, and again, you should not wrap the molecules after the rotation as they will be wrapped incorrectly. If you need to apply periodic boundaries at all in your analysis, then you should apply transformations to center the protein and unwrap the molecules but not align protein structure.

Let me know if you have any more questions about the lipyphilic transformations, and if the MDAnalysis transformatios I've suggested don't do what you're after, it's probably best to ask a question on the MDAnalysis forum.

I hope that helps!

Cheers,
Paul

P.S there may be some typos in the code as I wrote it directly here and haven't tested it

[sanjiv856]

Dear Paul,

Thanks a lot for the suggestions, very insightful. Your codes ran perfectly, without a typo.

Lipyphilic.transformations

The lipyphilic transformation you suggested:

import MDAnalysis as mda
from lipyphilic.transformations import triclinic_to_orthorhombic, nojump
u = mda.Universe(tpr, trr)
u.transfer_to_memory(step=10)
ag = u.select_atoms("all")
workflow = (
triclinic_to_orthorhombic(ag),
nojump(ag=ag, nojump_x=True, nojump_y=True, nojump_z=True),
)
u.trajectory.add_transformations(*workflow)

though I see that we are applying triclinic_to_orthorhombic(ag) before nojump , leads to this error:

ValueError: nojump requires an orthorhombic box. Please use the on-the-fly transformation :class: lipyphilic.transformations.triclinic_to_orthorhombic before calling nojump

I am using Box type – HEXAGONAL while creating the system in CHARMM-GUI, could this be an issue?

MDanalysis.transformations

The MDanalysis.transformations you suggested above, though the membrane is not breaking as it was with my previous transformations, the overall structure of the system (peptide and lipids) looks stretched as if PBC's are not applied, image attached.

Along with the standard analysis i.e., RMSD, RMSF, Rg, SASA etc. I wish to use almost all of the useful functionalities in lipyphilic i.e., lipyphilic.lib.area_per_lipid.AreaPerLipid. If I understand correctly, some of these analyses are not affected by PBC not being applied.

Again, thanks a lot for the help, for the MDanalysis.transformations, I will write to the forum you suggested.

With warm regards,

Sanjiv Kumar

[p-j-smith]

Hi @sanjiv856 sorry for taking over a year to respond! I imagine you have moved on from this now but just wanted to reply to say thank you for raising the issue about using nojump with triclinic boxes - it does indeed raise an error whereas it should be a warning.

And about the transformations - you're right that most of the functions in lipyphilic don't require molecules to be whole (and in fact many of them will wrap atoms back in to the primary cell). But unwrapping is useful for visualisation purposes, and nojump unwrapping is required for calculating diffusion coefficients

[sanjiv856]

Dear Paul (@p-j-smith) ,

Thank you for taking time to look at it. I still am working on those lines and experimenting LiPyphilic to get the results for some of the analysis (area per lipid and membrane thickness), it is extremely helpful.

I realised that if I do not use nojump, the transformation goes through and does the analysis well:
u.trajectory.add_transformations(triclinic_to_orthorhombic(ag=ag))

I was just worried if not including nojump is problematic for calculating area per lipid and membrane thickness.

I am still experimenting with LiPyphilic to calculate lateral diffusion of lipids using lipyphilic.lib.lateral_diffusion.MSD, lateral diffusion coefficient, and Lipid order parameters (SCD), but worried if not including nojump may be an issue.

It would be helpful for a newbie like me, if it could be added somewhere in documentation, preferably with every tool, if the analysis is affected by PBC, or nojump, unwrap or any such other transformation.

I know it should have been other way round but because of LiPhyphilic I stated learning more about MDAnalysis and Python in general.

Thanks again for reverting and this wonderful tool.

Regards,

Sanjiv

[p-j-smith]

I'm happy to hear lipyphilic is useful for you! The only analysis you should need to use nojump for is lipyphilic.lib.lateral_diffusion.MSD - all the other ones can be used without any transformation so long as your membrane is not broken across the z boundary (lipids being split across x or y is okay).

It would be helpful for a newbie like me, if it could be added somewhere in documentation, preferably with every tool, if the analysis is affected by PBC, or nojump, unwrap or any such other transformation.

That's a nice idea - I'll make an issue to add this to the docs - thanks for the suggestion!

I know it should have been other way round but because of LiPhyphilic I stated learning more about MDAnalysis and Python in general.

That's wonderful 😊 I'm so happy to have been your gateway into Python!

[lucyraven]

Hi!
I hope it's okay to jump onto this thread as my problem is included in that of Sanjiv above.

I am trying to apply the triclinic_to_ortherhombic + nojump transformations to my universe, and get the error:
ValueError: nojump requires an orthorhombic box. Please use the on-the-fly transformation :class:lipyphilic.transformations.triclinic_to_orthorhombic before calling nojump

However, let me run you through my process:
u = mda.Universe(gro, xtc)
protein = u.select_atoms('protein')
syst = u.select_atoms('all')

u.dimensions
Out: array([137.2366 , 137.2366 , 137.23651 , 59.999977, 59.999977, 90. ], dtype=float32)

u.trajectory.add_transformations(triclinic_to_orthorhombic(protein), nojump(protein, nojump_x=True, nojump_y=True, nojump_z=True))

Out: ValueError: nojump requires an orthorhombic box. Please use the on-the-fly transformation :class:lipyphilic.transformations.triclinic_to_orthorhombic before calling nojump

BUT, if I do:
u.trajectory.add_transformations(triclinic_to_orthorhombic(protein))

u.dimensions
​Out: array([137.2366 , 137.2366 , 137.23651, 90. , 90. , 90. ], dtype=float32)

To me, that sounds as if triclinic_to_orthorhombic does the right thing, but nojump doesn't recognize this. Also, I get the same problem when using the code snippet on lipyphilic transformation provided by Poul too. Can anyone help me understand whats going on?

Kindly,
Lucy

[sanjiv856]

Dear Lucy,

I hope someone more experienced than me answers your query.

I tried writing to the MDAnalysis forum, but got no response there.

Got few suggestions from GROMACS forum. Here is how I proceeded.

And tried MDAnalysis as follows, you may need to tweak the code a bit with proper file names:

 Trail 2:  Working, better then trial 1 , Failed for Membrane Only system 

u = mda.Universe(working_tpr, working_trr)

prot = u.select_atoms("protein")

lipid = u.select_atoms('resname POPE POPG and name P or protein')

reference = u.select_atoms("protein")

ag = u.atoms

transformations = (mda.transformations.unwrap(ag),
                    mda.transformations.center_in_box(lipid, center='mass'),
#                   mda.transformations.center_in_box(prot, center='mass'),
#                    mda.transformations.center_in_box(lipid, center='geometry'),
                   mda.transformations.wrap(ag, compound='fragments'),
#                    mda.transformations.fit_rot_trans(prot, reference)
                  )

u.trajectory.add_transformations(*transformations)

# nv.show_mdanalysis(u)
# Writing trejectrory
# n_atoms=len(u.atoms)
# .atoms.n_atoms

with mda.Writer(xtc_pbc_corrected_trial2, n_atoms=u.atoms.n_atoms) as W:
    for ts in u.trajectory:
        W.write(u)

Again, my python skills are not as good, just started learning. Please note how I selected lipids which did the trick.

Try multiple combinations in transformation.

Hope that gives you some direction until someone more experienced sheds some light.

Regards,

Sanjiv

[p-j-smith]

Hi @lucyraven sorry for the very slow reply! And thanks for reporting this - nojumo shouldn't raise an error here so I will fix that. It should just give a warning that you need to also use the triclinic_to_orthorhombic transformation

[lucyraven]

Dear Sanjiv, Thanks for taking time to reply! Lucy Kate Thøgersen, PhD Senior Specialist Raven Biosciences

…

________________________________ From: Sanjiv Kumar ***@***.***> Sent: Tuesday, May 31, 2022 14:18 To: p-j-smith/lipyphilic ***@***.***> Cc: Lucy Ladefoged Thøgersen ***@***.***>; Comment ***@***.***> Subject: Re: [p-j-smith/lipyphilic] Failing to implement triclinic_to_orthorhombic and nojump (Discussion #94) Dear Lucy, I hope someone more experienced than me answers your query. I tried writing to the MDAnalysis forum<https://groups.google.com/g/mdnalysis-discussion/c/u_-2D0AmM14/m/LfGjXFJRAQAJ>, but got no response there. Got few suggestions from GROMACS forum. Here is how I proceed<https://gromacs.bioexcel.eu/t/removing-pbc-for-peptide-lipid-system/3691>. And tried MDAnalysis as follows, you may need to tweak the code a bit with proper file names: Trail 2: Working, better then trial 1 , Failed for Membrane Only system u = mda.Universe(working_tpr, working_trr) prot = u.select_atoms("protein") lipid = u.select_atoms('resname POPE POPG and name P or protein') reference = u.select_atoms("protein") ag = u.atoms transformations = (mda.transformations.unwrap(ag), mda.transformations.center_in_box(lipid, center='mass'), # mda.transformations.center_in_box(prot, center='mass'), # mda.transformations.center_in_box(lipid, center='geometry'), mda.transformations.wrap(ag, compound='fragments'), # mda.transformations.fit_rot_trans(prot, reference) ) u.trajectory.add_transformations(*transformations) # nv.show_mdanalysis(u) # Writing trejectrory # n_atoms=len(u.atoms) # .atoms.n_atoms with mda.Writer(xtc_pbc_corrected_trial2, n_atoms=u.atoms.n_atoms) as W: for ts in u.trajectory: W.write(u) Again, my python skills are not as good, just started learning. Please note how I selected lipids which did the trick. Try multiple combinations in transformation. Hope that gives you some direction until someone more experienced sheds some light. Regards, Sanjiv — Reply to this email directly, view it on GitHub<#94 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AXUXHFJ7LDY5EPC5FOGIQXLVMX7R3ANCNFSM5OBJZVQQ>. You are receiving this because you commented.Message ID: ***@***.***>

[p-j-smith]

Hi @sanjiv856 and @lucyraven, there are a few issues with the nojump transformation (#105, #113), and now that there is an implementation in MDAnalysis I will remove the one in lipyphilic and recommend people to use MDAnalysis.transformations.NoJump instead. Thanks for reporting these issues, and sorry for any hassle!