using unwrap transformation

[Azadeh Alavi]

Dear all

I have a problem using unwrap module. The error I get is this:

 transform = mda.transformations.wrap.unwrap(ag)
AttributeError: module 'MDAnalysis.transformations' has no attribute 'wrap'

I am also not sure if I import it correctly:

from MDAnalysis import transformations

I also have another question. If this module works, then for the calculation of distance, the periodic images are taking into account automatically? Because now I am using the "distance_array()" module.

Thank you very much for your help.

Best regards,

Azadeh

[Oliver Beckstein]

Hi Azadeh,

The transformations.wrap() is currently only in the development version and will be in upcoming 1.0. 

However, if you use distance_array() then you should just supply the box argument with the unit cell because the function will then automatically take PBC into account (see docs for detail). You wouldn’t need transformations for that. 

Oliver 

Am 15.05.2020 um 06:25 schrieb Azadeh Alavi <azade...@gmail.com>:



--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/a4333ecb-8a2a-416c-97f4-89d073796f00%40googlegroups.com.

[Azadeh Alavi]

Hi Oliver

Thanks a lot for your prompt reply and the explanations. Sorry, I did not realize that it is under development. Is it possible that I use it in this mode?

Actually, besides distance, I also need to calculate a vector from one point of a molecule to another, therefore, the PBC should be taken into account. Is there any tool I can use for this purpose?

My exact problem is that there are two similar molecules in my system. For each molecule, I should calculate a vector from one point of the molecule to the other point, and then find two points on these two vectors which are at closest distance from each other (for this I consider 100 points on these two vectors). And then use these two points along with two other points on the two vectors to calculate a dihedral angle.

I would also need to calculate the angle between these two vectors.

Another point is that I want to calculate the com distance of these two molecules. I do know how I should change the command below so that the PBC is taken into account. I have also used the pbc flag, but it did not change anything.

com_distance = distance_array(sel1.center_of_geometry(), sel2.center_of_geometry(),box=u.dimensions)

Thank you very much for your help.

Best,

Azadeh

[Oliver Beckstein]

On May 16, 2020, at 4:23 AM, Azadeh Alavi <azade...@gmail.com> wrote:

Hi Oliver

Thanks a lot for your prompt reply and the explanations. Sorry, I did not realize that it is under development. Is it possible that I use it in this mode?

Yes, you just have to install from source; see https://www.mdanalysis.org/UserGuide/contributing_code.html#building-mdanalysis for help (scroll up/down for context).

Actually, besides distance, I also need to calculate a vector from one point of a molecule to another, therefore, the PBC should be taken into account. Is there any tool I can use for this purpose?

I don’t think that there’s an explicit function for that – something like (non-existent) `lib.calc_difference_vectors()`. Raise an issue, should be easy enough to add.

My exact problem is that there are two similar molecules in my system. For each molecule, I should calculate a vector from one point of the molecule to the other point, and then find two points on these two vectors which are at closest distance from each other (for this I consider 100 points on these two vectors). And then use these two points along with two other points on the two vectors to calculate a dihedral angle.

lib.distances.calc_dihedrals()  https://www.mdanalysis.org/mdanalysis/documentation_pages/lib/distances.html#MDAnalysis.lib.distances.calc_dihedrals

I would also need to calculate the angle between these two vectors.

lib.distances.calc_angles() https://www.mdanalysis.org/mdanalysis/documentation_pages/lib/distances.html#MDAnalysis.lib.distances.calc_angles

Another point is that I want to calculate the com distance of these two molecules. I do know how I should change the command below so that the PBC is taken into account. I have also used the pbc flag, but it did not change anything.

com_distance = distance_array(sel1.center_of_geometry(), sel2.center_of_geometry(),box=u.dimensions)

Possibly need to use options for center_of_geometry():

sel1.center_of_geometry(pbc=True, unwrap=True)

To be honest, I am not 100% sure on the options – try it out on an example where you know what the answer should be. Let us know what you expect to see and what you get instead. 

Oliver

Thank you very much for your help.

Best,

Azadeh

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/4da8ca1b-c5ca-4f79-9a54-dfab8f132513%40googlegroups.com.

--

Oliver Beckstein

orbe...@gmail.com

[Azadeh Alavi]

Dear Oliver

Thanks a lot again for your reply and the information.

Sorry, I forgot to say that I have also used the unwrap option. When I set both of these parameters to True, the error is:

ValueError: both 'pbc' and 'unwrap' can not be set to true

And when I set pbc to False:

MDAnalysis.exceptions.NoDataError: AtomGroup.unwrap() not available; this requires Bonds

When I transform the trajectory using the gromacs trjconv pbc tool, the result is the figure on the left. But when I use the normal trajectory, the result is the right one (please see the attachment).

Oh yes, I am already using calc_dihedrals(), however since as an input I give the points on the two vectors that I talked about, I was not sure if it would give a correct result.

Thank you very much.

Best regards,

Azadeh

[Lily Wang]

Dear Azadeh,

Sorry, I forgot to say that I have also used the unwrap option. When I set both of these parameters to True, the error is:

ValueError: both 'pbc' and 'unwrap' can not be set to true

‘pbc’ and ‘unwrap’ are mutually exclusive because they kind of do opposite things. ‘pbc’ accounts for periodic images when calculating distances, whereas ‘unwrap’ transforms your trajectory into the one image, typically going over the cell boundary. The docs have a nice image of this. 

And when I set pbc to False:

MDAnalysis.exceptions.NoDataError: AtomGroup.unwrap() not available; this requires Bonds

Your topology file needs to be in a format that contains bond information  or you will need to guess which atoms are bonded based on their distance from each other. Given that you are using GROMACS, the TPR format will do that (an important note for that format is that atom ID numbers and residue resid numbers are indexed from 0, if you are using those to create your vectors). Alternatively, you can choose to guess bonds. For example:

u = mda.Universe(my_top_file.gro, my_traj.xtc, guess_bonds=True)

However, this will depend on MDAnalysis correctly guessing the elements of your atoms, which is not a given.

Oh yes, I am already using calc_dihedrals(), however since as an input I give the points on the two vectors that I talked about, I was not sure if it would give a correct result.

You can pass a box into calc_dihedrals(a, b, c, d, box=u.dimensions) to take periodic images into account.

Hope that helps.

Cheers,

Lily

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/e1a64c59-ca00-4423-8c93-d1a3fa31d941%40googlegroups.com.
<bitmap.png>

[Azadeh Alavi]

Dear Lily

Thank you very much for your reply and the explanations.

Actually, as it is written in the link you have sent, the coordinates are not read from the TPR file. So I am wondering how one can use it in the universe. 

So I did this:

u = mda.Universe(gro_file.gro, traj_file.xtc, guess_bonds=True, vdwradii={'W':0.263, 'G':0.263, 'D':0.263, 'R':0.263 ...}

I had to give the vdwradii flag, since it was complaining. Of course I am not sure if I did it correctly.

And then gave the pbc=True for center of geometry. This is what I should do, right?

dist_com = distance_array(sel1.center_of_geometry(pbc=True), sel2.center_of_geometry(pbc=True),box=box)

The result is however the same as before.

Perhaps the bonds are not properly considered since when I gave the unwrap=True in the above command, the error was as follows. I am actually doing a coarse-grain simulation. Perhaps the only way is to transform the trajectory via gromacs...

 dist_com = distance_array(sel1.center_of_geometry(unwrap=True), sel2.center_of_geometry(unwrap=True),box=box)
 File "/home/aalaviza/.local/lib64/python3.6/site-packages/MDAnalysis/lib/util.py", line 1825, in wrapper                          
   return groupmethod(group, *args, **kwargs)                                                                                          
 File "/home/aalaviza/.local/lib64/python3.6/site-packages/MDAnalysis/core/groups.py", line 360, in wrapped   
   return function(group, *args, **kwargs)          

File "/home/aalaviza/.local/lib64/python3.6/site-packages/MDAnalysis/core/groups.py", line 880, in center_of_geometry
   return self.center(None, pbc=pbc, compound=compound, unwrap=unwrap)                                        
 File "/home/aalaviza/.local/lib64/python3.6/site-packages/MDAnalysis/lib/util.py", line 1825, in wrapper                             
   return groupmethod(group, *args, **kwargs)                                                                                          
 File "/home/aalaviza/.local/lib64/python3.6/site-packages/MDAnalysis/core/groups.py", line 360, in wrapped   
   return function(group, *args, **kwargs)                                                                                             
 File "/home/aalaviza/.local/lib64/python3.6/site-packages/MDAnalysis/core/groups.py", line 756, in center                                                                                                                                                                     
   coords = atoms.unwrap(compound=comp, reference=None, inplace=False)                                        
 File "/home/aalaviza/.local/lib64/python3.6/site-packages/MDAnalysis/core/groups.py", line 1615, in unwrap                                                                                                                                                                    
   positions = mdamath.make_whole(unique_atoms, inplace=False)                                                                         
 File "MDAnalysis/lib/_cutil.pyx", line 310, in MDAnalysis.lib._cutil.make_whole                                                                                                                                                                                               
ValueError: AtomGroup was not contiguous from bonds, process failed              

Thanks a lot.

Best,

Azadeh

[Oliver Beckstein]

Hi Azadeh,

Do you have a TPR file? Then try

u = mda.Universe(TPR, XTC)

This will read atoms and bonds from the TPR and coordinates from the XTC.

Oliver

[Azadeh Alavi]

Hi Oliver

Thanks a lot again for your favor.

I did it and gave the TPR and XTC files. The resid numbers were completely different !

Unfortunately, the result is the same as before and the pbc is not taken into account...

Best,

Azadeh

[Oliver Beckstein]

Hi Azadeh,

You might have to raise an issue in the issue tracker and provide detailed information on what goes wrong with the PBC handling. I’m out of simple ideas.

The residue numbers in a TPR are different but that’s a different issue.

Oliver

[Azadeh Alavi]

Hi Oliver

I am so sorry for my late reply and also because I was wrong. The center of mass distance was correctly calculated under the PBC when using the TPR file and guess_bonds in universe.

Many thanks to you and Lily, and of course to all the MDAnalysis community.

Best,

Azadeh