Re: [gmx-users] Center to center distance in cylindrical micelles

[Antonio Baptista]

You didn't give us much info about your system, so I will assume that your 
box is a rectangular prism and that all the cylinders are necessarily 
oriented along the z axis (e.g., because they are "infinite"/periodic). 
Since the center of mass (COM) of each cylinder must be a point on its 
axis, you simply need to compute the COMs of each cylinder and then 
compute any intended distances between their axes using just the xy 
components.



On Fri, 17 Aug 2018, Shan Jayasinghe wrote:


Hi All,

Does anyone know how to calculate distance between two adjacent cylinders
in a  hexagonal array?

Thank you.


On Sun, Aug 12, 2018 at 8:01 AM Shan Jayasinghe <
shanjayasinghe2...@gmail.com> wrote:


Hi,

Actually, I want to calculate the lattice parameter (distance between two
cylinders) of hexagonal phase. Appreciate, if you could help me.

Thank you.




On Sat, Aug 11, 2018 at 11:11 PM Mark Abraham 
wrote:


Hi,

A center is a point. I don't see how you would define a central point of
an
infinite cylinder.

Mark

On Sat, Aug 11, 2018, 14:47 Shan Jayasinghe 


wrote:


Hi,

I'm sorry for not providing the detail description of my system.

Actually

it is a hexagonal phase. Therefore, four cylindrical micelles can be

seen

in a cell. I want to calculate the center to center distance of two
cylindrical micelles.

On Sat, Aug 11, 2018 at 8:35 PM Mark Abraham 
wrote:


Hi,

There is an axis of each cylinder but the two won't be parallel so

there

is

no unique distance to measure. So I suggest thinking carefully about

what

you really want.

Mark

On Sat, Aug 11, 2018, 06:58 Shan Jayasinghe <

shanjayasinghe2...@gmail.com>

wrote:


Hi,

I tried to used gmx distance. However, I don't understand how can I

define

the center of the circle face of the cylinder to the same in the

other

cylindrical micelle.

Can anyone help me?
Thank you.


On Sat, Aug 11, 2018 at 9:33 AM Mark Abraham <

mark.j.abra...@gmail.com



wrote:


Hi,

Have you looked at the different tools available and considered

what

might

be useful for you?

Mark

On Fri, Aug 10, 2018, 07:34 Shan Jayasinghe <

shanjayasinghe2...@gmail.com>

wrote:


Dear Gromacs Users,

I want to calculate the center to center distance of two

cylindrical

micelles in my simulation. What gmx command should I use to

calculate

the

distance? Can anyone help me?

Thank you.
--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List

before

posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit


https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users

or

send a mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List

before

posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit


https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users

or

send a mail to gmx-users-requ...@gromacs.org.




--
Best Regards
Shan Jayasinghe
--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users

or

send a mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.




--
Best Regards
Shan Jayasinghe
--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.




--
Best Regards
Shan Jayasinghe




--
Best Regards
Shan Jayasinghe
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to 

[gmx-users] Typo in table 3.1 of the reference manual?

[Antonio Baptista]

Dear all,

There is an inconsistency in the truncated octahedron entry in table 3.1 
of the reference manual (version 2018 and all previous ones I checked). 
The value indicated for the bc and ab angles is 71.53 degrees, but the 
vectors a, b and c indicated in the table give the value 70.53 degrees = 
acos(1/3). The indicated vectors correspond indeed to one of the possible 
sets of primitive vectors in the body-centered cubic lattice of which the 
truncated octahedron is the Wigner-Seitz primitive cell (namely the set 
with origin at a cube corner and pointing to three cube centers), so I 
belive that the typo is in the angles, not the vectors. Or am I missing 
something?


Best,
Antonio

--
António M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Molecular Simulation using Ionic liquid

[Antonio Baptista]

Hi Vidya,

You should try to make your model "fit" into the existing forcefield you 
will be using. Also, you should validate the model with known physical 
properties (e.g., density, self-diffusion, shear viscosity, isothermal 
compressibility). It's usually difficult to get them all good, of course, 
but playing with some (e.g., LJ) parameters may get you a reasonable 
compromise. An example of a parameterization of Bmim for a GROMOS ff: 
http://dx.doi.org/10.1021/jp061869s

And its use with a protein: http://dx.doi.org/10.1021/jp0766050

Best,
Antonio


On Tue, 25 Jul 2017, Mark Abraham wrote:


Hi,

How have others done similar work? Hint: extending a protein forcefield
means understanding how its parameters were developed. You can't just
assume that your favourite tool is the one for the job.

Mark

On Tue, 25 Jul 2017 08:37 Vidya Sundaram  wrote:


Hi,

I intend to simulate a protein of interest in the ionic liquid,1-Butyl-3-
methylimidazolium chloride (BmimCl). In order to generate the topology and
coordinate files of the IL, I tried optimizing the structure through
Gaussian09. Kindly clarify my doubts regarding the same.

1. Should the cation and anion structures be generated separately in
vacuum? or should it be generated by using water as solvent?

2. If it is to be generated in solvation model, I guess Gaussian needs
solvent descriptors (eps, epsinf and hydrogen bond acidity etc to be
specified when using SCF=SMD and Solvent= Generic). In that case, how do I
give these solvent descriptors separately for anion and cation?

Any information could really help.

Thanks in advance.

Vidya
--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Calculation of nematic order parameter using gromacs

[Antonio Baptista]

On Sat, 3 Jun 2017, nidhi sorout wrote:


Hello,

Thank you Antonio..

But my angle of interest is the angle between the molecular axis of protein
and the director. I am not able to understand here, from where I can choose
this 'director'?


In your case, what does 'director' mean exactly?



Nidhi

On 30 May 2017 8:55 p.m., "Antonio Baptista" <bapti...@itqb.unl.pt> wrote:


Hi Nidhi,

If I remember correctly (and your use of "p2" suggests so), that should be
the ensemble average of the 2nd-order Legendre polynomial of the angle
between the molecular axis and the membrane normal, right?

Although the order parameter computed by "gmx order" uses that same
definition, it takes each C-H bond of the aliphatic lipid tail, not the
overall molecular axis. So, "gmx order" is not what you want.

You can in principle compute what you need in two steps: (1) use "gmx
gangle" to compute the angle of interest for all molecules and all
snapshots (you will have to defined the vector of interest, say as the one
connecting the tail to the head); (2) do a small script to compute the
average from those data.

Best,
Antonio


On Tue, 30 May 2017, nidhi sorout wrote:

Dear All,

I want to calculate the nematic order parameter (p2) at each time step.
Is it possible to do this with "gmx order"?
If not than  please suggest the right way.

Thank you,
Nidhi
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support
/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.

--

Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support
/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Calculation of nematic order parameter using gromacs

[Antonio Baptista]

Hi Nidhi,

If I remember correctly (and your use of "p2" suggests so), that should be 
the ensemble average of the 2nd-order Legendre polynomial of the angle 
between the molecular axis and the membrane normal, right?


Although the order parameter computed by "gmx order" uses that same 
definition, it takes each C-H bond of the aliphatic lipid tail, not the 
overall molecular axis. So, "gmx order" is not what you want.


You can in principle compute what you need in two steps: (1) use "gmx 
gangle" to compute the angle of interest for all molecules and all 
snapshots (you will have to defined the vector of interest, say as the one 
connecting the tail to the head); (2) do a small script to compute the 
average from those data.


Best,
Antonio


On Tue, 30 May 2017, nidhi sorout wrote:


Dear All,
I want to calculate the nematic order parameter (p2) at each time step.
Is it possible to do this with "gmx order"?
If not than  please suggest the right way.

Thank you,
Nidhi
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Order Parameter for HII phase

[Antonio Baptista]

Hi Mohsen,

I suggest you have a look at the QRB 1977 review by Seelig, that Tom 
already mentioned. Like for the planar case, they discuss how spectra of 
cylindrical lipid phases are related to the order parameter tensor (but I 
never looked into the details for that case). Anyway, I don't know if you 
want to compare your simulations with experimental data or are just 
looking for a convenient structural parameter. But, whatever the case is, 
reading the literature and thinking about the geometry of your system 
(Duliez's papers are a good train for that) should give you some hints for 
a relevant parameter for your system.


Once you have selected a parameter, you can compute it from stratch using 
other GROMACS tools besides g_order. In particular, you can use g_gangle 
to get several sorts of angles and then process them yourself. As an 
example, you can use this approach to compute SCD using the equation 
already mentioned by Tom, and then check if it agrees with the SCD given 
by g_order (I once did that, just to be sure).


Good luck! :)

Best,
Antonio


On Thu, 23 Feb 2017, Justin Lemkul wrote:




On 2/23/17 1:25 PM, Mohsen Ramezanpour wrote:

And I agree with Piggot. The paper by Chau is about on option in g_order.



Yes, and if memory serves (been a while since I used the program, so perhaps 
I am confusing something), this is what -o provides you so I thought it was 
relevant.  -od gives you the deuterium order parameters.  The documentation 
for gmx order is indeed very sparse.



Anyway, a general question:
Can we expect to find a published article for each/some module(s) in
Gromacs?


No.  Many of the tools are just implementations of simple algorithms used 
widely in simulations.  Some do have specific references and those are 
generally printed to the screen output in bold blocks of text.


With regards to deuterium order parameters, there is a well accepted and 
ubiquitously used equation, which Tom posted.



I mean how/where can we figure out the underlying algorithm and
comprehensive description of each analysis tool?



That's why GROMACS is open source :)  Efforts are certainly always made to 
provide references when possible.


-Justin


Cheers
Mohsen


On Thu, Feb 23, 2017 at 10:06 AM, Mohsen Ramezanpour <
ramezanpour.moh...@gmail.com> wrote:


Hi Justin, Piggot,

Thanks for your replies.
I agree with that. The problem is that the situation is straightforward
for bilayers as bilayers are usually in specific angles regarding the
magnetic field (e.g. they are perpendicular, as it is the case in
Vermeer study.)
For example the normal to the bilayer is z and we do not need be worry
about anything else.
However, in the case of HII phase, there are cylinders which are not
perfectly cylinder.
Even if they are perfect cylinders, I think I should use "membrane radial
axis" to mimic the local normal vector to the lipid surface.

Cheers
Mohsen


On Thu, Feb 23, 2017 at 8:08 AM, Piggot T.  wrote:


Hi,

I don't believe this Chau paper is for deuterium order parameters (as I
think you are most likely referring to) but is related to some of the 
other
order parameter options in gmx order. The Vermeer paper you mention gives 
a

good overview of deuterium order parameters and you can find more details
in the papers referenced therein (e.g. the Douliez et al. and the Seelig 
et

al. papers in particular).

For united-atom systems (as gmx order assumes), I believe the calculation
is performed in gmx order using the "SCD = (2/3) Sxx + (1/3) Syy" 
equation

(which is nicely derived in the 1998 Douliez paper:
http://aip.scitation.org/doi/pdf/10.1063/1.476823) . The calculation
requires an appropriate definition of the molecular frame of the system,
where the z-axis is the bilayer normal.

Note that if you are calculating order parameters for unsaturated carbons
(e.g. as in carbon-carbon double bonds), gmx order doesn't currently do
this correctly as it requires a different calculation with a different
molecular frame.

Cheers

Tom


From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se [
gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin
Lemkul [jalem...@vt.edu]
Sent: 23 February 2017 12:16
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Order Parameter for HII phase

On 2/22/17 9:51 PM, Mohsen Ramezanpour wrote:

Dear Gromacs users,

Unfortunately, I did not get any reply on this post.

I was wondering if anyone is aware of the original paper on g_order so

that

I can understand the underlying algorithm of the options?
I am not sure how to use it for HII phase order parameters.

I found this article but I am not sure if these are the ones, as they do
not describe g_order.

"A new order parameter for tetrahedral configurations" by Chau and

Hardwick

(1998)


This is the reference given at the end of the gmx order help description.

-Justin


and
"Acyl chain order parameter profiles in phospholipid bilayers:

Re: [gmx-users] gmx gangle

[Antonio Baptista]

Hi João,

I used gangle to compute several angle distributions relative to the 
z-axis membrane normal, but that was almost a year ago, so I don't recall 
the details. But I can tell you that, after many failed attempts, I gave 
up from using fancy stuff in the -group1 and -group2 options. Instead, I 
ended up writing a small bash script to process one atom pair at a time 
(just one group in the .ndx file per gangle run) and to collect everything 
at the end.


I understand that your problem is a bit different (your second vector is 
not an axis, but defined by a second atom pair), but maybe a similar 
"minimalist" approach could get the job done. Or not... Good luck!


Best,
Antonio

On Wed, 15 Feb 2017, João Henriques wrote:


Hello everyone,

It seems that my inquiries about gmx gangle and selections are either too
funky or too basic to justify getting a single answer from the community.
Still, I will try once more with an even simpler scenario, which will
hopefully attract an answer (even if it's just someone telling me to stop
reposting, despite the fact that it is not 100% a repost).

I have a trajectory where the water molecules are ordered by proximity to
the protein. I also have an index file with the N closest waters that I am
interested in analyzing. Now, my plan is to simply compute the O-H angle
distribution within the group composed by these N water molecules. For that
I type something along these lines:

gmx gangle -f xxx.xtc -s xxx.tpr -n xxx.ndx -g1 vector -group1 "group
"hydration_shell" and name OW HW1" -g2 vector -group2 "group
"hydration_shell" and name OW HW1" -oav -oall -oh

To my surprise, this command returns something completely different from
what I expected. Instead of the angular distribution function of all pairs
of O-H vectors from different water molecules, it outputs zero. Obviously,
this means that the command is computing the angle between the O-H vector
of each water molecule and itself. How can I modify my selection to produce
the desired result?

Thank you in advance,
Best regards,
João
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] GROMOS 54a7 parameters for rare amino acids (eg. beta-alanine)

[Antonio Baptista]

Hi Billy,

I know that people from the University of Minho, Portugal, have 
parameterized several non-proteinogenic amino acids for GROMOS 54A7, 
making the parameters freely available. At least two papers are published:


http://dx.doi.org/10.1021/jp4074587
http://dx.doi.org/10.1021/jp505400q

But you may also want to check with the authors, since other papers with 
more parameters had been submitted when the PhD student doing that work 
(Tarsila Castro) defended her thesis in April.


Hope this helps.

Best,
Antonio


--
António M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--

On Thu, 23 Jun 2016, Billy Williams-Noonan wrote:


Hi Gromacs Users,

 Does anyone know where I can find a set of download-able parameters for
rare amino acids?  Is there a repository of some kind?

Kind regards,

Billy

--
Billy Noonan*|*PhD Student*|*Bsci ( *Adv* ), IA Hon

*LinkedIn Profile

**|*   +61420 382 557

Monash Institute for Pharmaceutical Sciences ( *MIPS* )
Royal Parade, Parkville, 3052
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] PCA problems

[Antonio Baptista]

On Wed, 8 Jun 2016, Tsjerk Wassenaar wrote:


Hey :)




 That usually gives a fitted ensemble that more closely retains the
original RMSD values between all pairs of structures.



This should read: ... a fitted ensemble of which the sum of the traces of
all pairwise inner product matrices is closer to minimal.
The pairwise RMSDs (and their sum) remain what they are.


Yep, I didn't phrase it the best way. Thanks for clarifying... :)



Cheers,

Tsjerk

--
Tsjerk A. Wassenaar, Ph.D.
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] PCA problems

[Antonio Baptista]

Hi James,

If your molecule shows some flexiblility, I would suggest using as a 
reference the structure of your original ensemble that produces the fitted 
ensemble with the lowest sum of RMSD values to that structure (or their 
square). That usually gives a fitted ensemble that more closely retains 
the original RMSD values between all pairs of structures.


If your molecule is well-structured, it wouldn't make much difference 
which reference you use, as Tsjerk said.


Best,
Antonio


On Wed, 8 Jun 2016, jkrie...@mrc-lmb.cam.ac.uk wrote:


ok thanks Tsjerk. I think that makes sense now.

Best wishes
James


Hi James,

That's silly! Ambiguous means that the same structure can have multiple
solutions in a fit. The fit to a single reference structure (with more
than
three atoms) is never ambiguous. Can never, by definition!

Now if you have two reference structures at hand, and they have (quite)
different structures, then fitting on one may give a different ensemble
from fitting on the other. The fit is not consistent, and the
inconsistency
is worse for flexible molecules. Different ensembles will mean different
correlations, thus giving different principal components.

Progressive fitting does not solve the problem. In fact, progressive
fitting _is_ ambiguous. Let's say we have a series of conformations ABCAC.
Then we  fit C once to B, which was fitted to A, and later we have C
fitted
to A, which was fitted to the previous C. Note that in practice the
situation will be much worse as we can approach a certain configuration
from many sides. Using B as reference will yield a fit that is different
from using A as reference, so the structure C will have two different
orientations in the resulting ensemble. Hence, the fit is ambiguous.

For structured proteins, the difference will not matter much. However, in
long trajectories there may be an added contribution (drift) of the
orientation.

Hope it helps,

Tsjerk

On Wed, Jun 8, 2016 at 5:33 PM,  wrote:


Thanks Tsjerk,

Isn't the progressive fit supposed to rotate everything back into the
same
orientation without having to worry about inferring that orientation
from
a reference structure that doesn't align well? Each configuration should
in theory align well to its predecessor all the way back to the starting
structure (which is what I'd usually take as a reference anyway).

The original note I was thinking of says as follows:

'''Before a PCA, all structures should be superimposed onto a common
reference
structure. This can be problematic for very flexible systems such as
peptides,
where the fit may be ambiguous, leading to artificial structural
transitions. In
certain cases, such problems may be alleviated by using a progressive
fit,
where
each structure is superimposed onto the previous one. It is also
important
to note
that when results of different PCAs are to be compared with each other,
then
each individual PCA should be based on the same reference structure used
for
superposition.'''

Please could you explain further what it is I have misunderstood.

Also would you say a progressive fit is a bad idea for more structured
proteins?

Many thanks
James


Hi James,

'Spurious alignment' is the dependence of the resulting ensemble on

the

reference structure. Unfortunately, that's not solved by a progressive
fit.
Rather, in a progressive fit, the same configuration can have multiple
orientations, based on the previous structures, which is also

problematic

when you're trying to understand spatial correlations between atoms

within

their reference frame.

Cheers,

Tsjerk

On Wed, Jun 8, 2016 at 9:27 AM,  wrote:


Dear Teresa,

That sounds like a periodic boundary issue to me. It could be fixed

by

using a tpr instead of a gro as the gmx covar manual says "All
structures
are fitted to the structure in the structure file. When this is not a
run
input file periodicity will not be taken into account." Alternatively

if

you don't have a tpr you could use gmx trjconv first with -pbc whole

or

-pbc nojump.

I also remember reading (I think it was in the Hayward and de Groot
review
2008) that fitting peptides to a reference structure can cause

spurious

alignments. I don't know if this is also related to what you're

seeing

but
it might be worth using gmx trjconv again with-fit progressive then

use

-nofit in gmx covar.

Best wishes
James


Dear GROMACS community

I am trying to complete a PCA analysis of my peptide adsorbed to a
surface. However when I use :

gmx covar -s trajectory.gro  -f md_golp_vacuo.xtc

and select the protein for both the least squares fit and

covariance

calculation, followed by


gmx anaeig -s trajectory.gro -f md_golp_vacuo.trr -filt filter1.gro
-first 1 -last 1 -skip 100

and I select the peptide for the least squares and covariance
calculation

My peptide is now broken up into pieces. Is this right?



Best
Teresa
--
Gromacs Users mailing list

* Please search the archive 

Re: [gmx-users] Correct method to do Cartesian PCA

[Antonio Baptista]

On Fri, 8 Jan 2016, Bin Liu wrote:


Hi Tsjerk,

I replicated their settings in their papers. The system size I used is
larger than what they had. Could you elaborate on the reference used for
fitting? Thank you so much.


About the importance of the reference structure used for fitting in 
Cartesian PCA, you may want to look at this: 
http://dx.doi.org/10.1021/jp902991u


Best,

--
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--




Cheers,

Bin

Message: 3
Date: Fri, 8 Jan 2016 18:29:26 +0100
From: Tsjerk Wassenaar 
To: Discussion list for GROMACS users 
Subject: Re: [gmx-users] Correct method to do Cartesian PCA
Message-ID:
   
Content-Type: text/plain; charset=UTF-8

Hi Bin,

Did you use the same force field, simulation length, other parameters? The
system size might play a role too. And for trialanine the reference used
for fitting may matter too. But first check all the other
similarities/differences.

Cheers,

Tsjerk
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Force-field bias???

[Antonio Baptista]

On Wed, 5 Aug 2015, Justin Lemkul wrote:




On 8/5/15 10:32 AM, Smith, Micholas D. wrote:
OPLS-AA 2001 (the one implimented in gromacs) can be a little slow in 
folding

(in both my experience, and from   doi:10.1016/j.bpj.2009.04.061 ). That is
not necessarily a bad thing if the peptide is suppose to take a while to
fold. If you just want to get to the folded structure, and you don't care
about the pathway, you could also use an enhanced sampling technique.

An alternative force-field that seems to handle beta structures pretty 
well,

but can sometime accelerate the folding is the older GROMOS53A force-field;
though since it is united atom, I have found that sometime the contacts 
don't

quite match with NMR data.



Beware, it's well known that 53A6 systematically under-stabilizes helices 
such that you get inappropriate unfolding.  It samples the unfolded state 
quite well, but again, this may be application-specific.


That under-stabilization of helices in 53A6 seems to have been mostly 
solved by the few changes that led to 54A7, which is able to reproduce 
both the thermodynamics and the kinetics of helix formation in a cyclic 
peptide model (doi: 10.1021/ct400529k). So, I would say that 54A7 is 
largely a fixed 53A6 in that respect (plus a few changes on choline 
groups and ions).


Of course, the fact that the kinetics seems good is not necessarily a good
thing if you just want thermodynamics (relative populations), because
sampling may turn out to be too slow for your system... or not...

Best,
Antonio

--
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--



-Justin


Force-field choice is a non-trival, but largely difficult problem to deal
with. In general, though, I tend to avoid the older Amber force-fields when
simulating beta-rich peptides as have had a bias towards helical structures
(see the references within the paper I noted above).

Good Luck,

Micholas

=== Micholas Dean Smith, PhD. Post-doctoral Research
Associate University of Tennessee/Oak Ridge National Laboratory Center for
Molecular Biophysics

 From:
gromacs.org_gmx-users-boun...@maillist.sys.kth.se
gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of SAPNA 
BORAH

sapnauser...@gmail.com Sent: Wednesday, August 05, 2015 10:20 AM To:
gmx-us...@gromacs.org Subject: Re: [gmx-users] Force-field bias???

Dear all,

Thank you for your answers which definitely clears some of my doubts. I am
trying to unfold a lectin tetramer protein that comprises largely of beta
strands, and I have used the OPLS ff, so far I have not been successful, 
can

OPLS be a problem at all??

Thanks..

Sapna Mayuri Borah c/o Dr. A. N. Jha Research student Tezpur University,
India

On Wed, Aug 5, 2015 at 5:44 PM, Smith, Micholas D. smit...@ornl.gov 
wrote:



You should also note that force-field bias may be system dependent. So if
you are unsure of the dynamics from the onset, you should likely test a
variety of force-fields and/or try to find a comparison to NMR (or other)
experimental methods.

-Micholas

=== Micholas Dean Smith, PhD. Post-doctoral Research
Associate University of Tennessee/Oak Ridge National Laboratory Center for
Molecular Biophysics

 From:
gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of Justin
Lemkul jalem...@vt.edu Sent: Wednesday, August 05, 2015 7:31 AM To:
gmx-us...@gromacs.org Subject: Re: [gmx-users] Force-field bias???

On 8/5/15 6:24 AM, SAPNA BORAH wrote:

Dear all,

I have a general query about unfolding simulations. Is there a bias
among force-fields selected for unfolding, i.e. is it possible that
unfolding

of

the protein may be different with a change in force fields applied?
Literature has given some support on this, however, I am not sure how

this

may happen.



Yes, this is true.  It usually comes down to errors in the balance between
intra-protein and protein-water interactions.  Incorrect balance favors
either the folded or unfolded state.  No force field is perfect, but some
are better than others and there's lots of literature on this topic
regarding which are the most successful.

-Justin

-- ==

Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences School of Pharmacy Health Sciences
Facility II, Room 629 University of Maryland, Baltimore 20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

== -- Gromacs Users 

Re: [gmx-users] definition of eigenvector in gromacs

[Antonio Baptista]

Tsjerk, that's one of the coolest explanations of PCA I've ever read! :)

Just one additional comment for Brett: Tsjerk's example assumes implicitly 
that your hands are actually moving parallel to the desk on average. On my 
desk, the meal is usually on top of the pile of papers that mysteriously 
keeps growing to the left of my keyboard. So, my hands don't move parallel 
to the desk, but along an inclined plane that raises from the keyboard to 
the meal position and is slightly tilted towards me (yeah, I know I should 
sit straight...). Therefore, the average position of my hands and their 
direction of major extent (eigenvector 1) should be marked on that plane, 
not on the desk, and the second direction (eigenvector 2) is a line also 
on that plane drawn perpendicularly to the first. As in Tsjerk's example, 
if I project every position of my hands onto this resulting PCA plane I 
will probably capture most of their motion, without having to consider the 
third direction perpendicular to the plane.


My point is that the subspace obtained from PCA is often different from 
the one that seems more natural to us (which was perhaps not obvious in 
Tsjerk's example). Although my desk is the natural reference frame when 
I'm working, my hands' motion can be more easily described in terms of 
positions on the plane I've just described, since I roughly need just 2 
coordinates instead of the original 3. That's the aim of PCA -- to 
describe the original motion using fewer coordinates.


Incidentally, this also shows that a messy desk is scientifically more 
fruitful than a tidy one... ;)


Cheers,
Antonio

On Mon, 9 Mar 2015, Tsjerk Wassenaar wrote:


Hi Brett,

Let's say you're sitting at your _desk_ writing that paper with a deadline
yesterday and you put a quick _meal_ next to you, wondering why on _earth_
you keep up with this. Your hands are moving between the meal and the
keyboard. You notice that the average position of your hands is somewhere
between the two and mark the mean position on your desk. Then you draw a
line through it that corresponds to the major extent of the motion of your
hands, and write 'eigenvector 1' along it. You add a line through the
average position, perfectly perpendicular to the first, and write
'eigenvector 2' along it. Now you can project every position of your hands
onto your desk, giving it an 'eigenvector 1' coordinate (or score) and an
'eigenvector 2' coordinate (or score). You notice that it's only part of
the total motion, as you neglect the height, which will be a third line,
perpendicular to the desk.

You can look at it a bit differently and say that your desk is the subspace
of your real space, spanned by the two perpendicular vectors, which
together describe most of your hand's motion.

I hope this makes some sense :)

Cheers,

Tsjerk





On Mon, Mar 9, 2015 at 2:33 AM, Brett brettliu...@163.com wrote:


Dear All,

Recently I have read a gromacs related article which contain a sentence
Trajectories pojected on the first two eigenvectors. Will you please tell
me the meaning of eigenvector and Trajectories pojected on the first two
eigenvectors? Will you please also explain it to me using common-used
words simple as desk, meal, earth, etc?

Brett


--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.





--
Tsjerk A. Wassenaar, Ph.D.
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.



--
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--

--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] PCA projections using g_anaeig

[Antonio Baptista]

Hi James,

Yes, but all you need are the 1D projections, since they are the 
coordinates in the PCA reference frame. So, just combine them in any way 
you want, as Tsjerk indicated. Also, keep in mind that there is actually 
no transformation of the conformation space when you do PCA, just a 
rotation of the coordinate axes.


But beware of 2D projections alone -- two different 3D distributions can 
give *exactly* the same 1-vs-2, 1-vs-3 and 2-vs-3 projections! So, if you 
want to get a better understanding of your distribution (as seems to be 
the case) try also to look at least at the 3D projection 1-vs-2-vs-3. For 
more dimensions you cannot simply look, of course... :) And be skeptical 
about the basins you see in 2D PCA maps, since each of them typically have 
many disparate conformations lumped together. You may want to look at 
http://dx.doi.org/10.1021/jp902991u


Best,
Antonio


On Wed, 28 Jan 2015, James Starlight wrote:


Thx, Tsjerk!

So normally gromacs make 2D projections  exactly only on 1 and 2 components
(even if 100 modes have been calculated for instance) doest it?

2015-01-28 16:13 GMT+01:00 Tsjerk Wassenaar tsje...@gmail.com:


Hi James,

If you have the 1D projections as 1.xvg and 3.xvg:

paste 1.xvg 3.xvg | awk '/^[^;@#]/{print $2, $4}'  2d.dat

Cheers,

Tsjerk
On Jan 28, 2015 3:25 PM, James Starlight jmsstarli...@gmail.com wrote:

 Dear Gromacs users!

 I'd like to ask about possibility to make projection of the MD trajectory
 onto the subset of selected PC calculated by g-covar


 e.g using

 g_anaeig –v eigenvec.trr –s md.tpr –f md.xtc –first 1 –last 3  -2d


 produce projection of the md.xtc onto the 1 and 2 components although 1-3
 modes has been calculated. How it would be possible to make the same
graph
 using 1 and 3rd components instead?


 Thanks for help,

 James
 --
 Gromacs Users mailing list

 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] constant pH simulations in GROMACS 4.5

[Antonio Baptista]

Hi,

Our PB/MC/MD approach to constant-pH simulations is implemented up to 
version 4.0.7 of GROMACS, but we don't currently have a usable 
distribution -- after using some collaborators as guinea pigs, we came to 
the conclusion that the current implementation is unintelligible, since no 
one besides those who wrote the scripts and are familiar with all the 
interfaced programs (GROMACS+MEAD+PETIT+etc) could actually use it... :-( 
We have already wrote a new and cleaner implementation which is being 
tested, so I hope we could finally make it available soon.


Best,
Antonio

--
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--

On Thu, 22 Jan 2015, Gerrit Groenhof wrote:


Hi,

Our lamda-dynamics based approach works only in gromacs 3.3.3 at the moment.

I don't know the status of the PB/MC/MD constant pH approach by Victorio 
Baptista and co-workers. Maybe that works in combination with gromacs 4.5, but 
he can tell you that.

Best,

Gerrit



Message: 6
Date: Thu, 22 Jan 2015 19:29:45 +0530
From: Sahithya S Iyer sah2...@gmail.com
To: gmx-us...@gromacs.org
Subject: [gmx-users] constant pH simulations in GROMACS 4.5
Message-ID:
cacprdiqc9ufnbogcym0n0xbvy7pfndd4dz0fvzn1-gevtcb...@mail.gmail.com
Content-Type: text/plain; charset=ISO-8859-1

Hi users
Have standard methods been developed in
Gromacs
4.5 version to perform constant pH simulations ?
--
Sahithya



Our methods so far only works in gromacs 3.3.3.

I don't know the status of the constant pH approach by victorio baptista and 
co-workers. Maybe that works in combination with gromacs 4.5

Best,

Gerrit



--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Pressure Question

[Antonio Baptista]

On Thu, 6 Nov 2014, Téletchéa Stéphane wrote:


Le 06/11/2014 06:16, Antonio Baptista a écrit :


In particular, the virial-based instantaneous pressure (call it P') 
computed in simulations has its ensemble average equal to the thermodynamic 
pressure P (check any good book on molecular simulation). But, as others 
already pointed out, this P' is well-known to show extremelly large 
fluctuations, meaning that its average computed from the simulation has 
usually a very large statistical spread. In other words, although the 
ensemble average of P' is strictly equal to P, its simulation average is a 
random variable that often shows large deviations from P (especially for 
short simulations). To get an idea of what is an acceptable error for the 
average of P', you may look at its distribution histogram in the NPT 
simulation. 


Dear Antonio,

Sorry if my message sound aggressive when I talked about totally 
irrevelevant, I will clarify my thoughts.


No problem, Stéphane. I was just trying to avoid propagating the wrong 
idea that a parameter is irrelevant in an ensemble where its value is not 
explicitly imposed, an idea I saw stated before in several discussions. 
Anyway, it seems I misunderstood you, since you say you were actually 
making the same point... :)


Best,
Antonio



From a theoretical point of view, you are right, each ensemble is accessible.

From a biological point of view, though, the concept of fixing the volume is 
less reasonable:
we live at constant pressure and temperature, and also at tighly controlled 
pH, and salt concentrations.


The volume varies though, as you feel it when the weather is getting hot or 
cold.


My point was exactly what your are telling in a more formal way than me:
this P' is well-known to show extremely large fluctuations

Well, digging a bit more on my feeling, I also found opposite arguments on 
the AMBER mailing list,

like here: http://archive.ambermd.org/201103/0431.html

So I'll got back again on my research and adjust my mind on the actual 
bleeding edge simulations

taking into account all the recent code and force fields progresses.

Best,

Stéphane

--
Team Protein Design In Silico
UFIP, UMR 6286 CNRS,
UFR Sciences et Techniques,
2, rue de la Houssinière, Bât. 25,
44322 Nantes cedex 03, France
Tél : +33 251 125 636
Fax : +33 251 125 632
http://www.ufip.univ-nantes.fr/ - http://www.steletch.org

--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!


* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Pressure Question

[Antonio Baptista]

Well, it is definitely *not* totally irrelevant to talk about pressure 
when doing simulations in the NVT (canonical) or NVE (microcanonical) 
ensembles. Pressure, like temperature, volume, energy, numbers of 
particles, etc, is a thermodynamic property which is *always* defined for 
*any* system in equilibrium. These parameters characterize the 
thermodynamic state, regardless of the ensemble that you choose to 
describe the system in terms of statistical mechanics. And the microscopic 
counterparts of these thermodynamic parameters are defined in such a way 
that their ensemble average must necessarily equal the thermodynamic 
value, regardless of ensemble (although their fluctuations are 
ensemble-dependent, becoming asymptotically identical only for macroscopic 
systems).


In particular, the virial-based instantaneous pressure (call it P') 
computed in simulations has its ensemble average equal to the 
thermodynamic pressure P (check any good book on molecular simulation). 
But, as others already pointed out, this P' is well-known to show 
extremelly large fluctuations, meaning that its average computed from the 
simulation has usually a very large statistical spread. In other words, 
although the ensemble average of P' is strictly equal to P, its simulation 
average is a random variable that often shows large deviations from P 
(especially for short simulations). To get an idea of what is an 
acceptable error for the average of P', you may look at its distribution 
histogram in the NPT simulation.


As for the equilibration of the system, the only thing that matters is 
which termodynamic state you are aiming at and what is the best way to get 
there. For example, if you choose NVT but happen to start with a volume 
which is a bit too large (eg, because the parameterized model acquires a 
higher density than the true experimental value that you assumed when 
setting the box size), you may get into trouble because the system may 
then want to separate into two phases but, being unable to do so in a 
small simulation box, ends up in a weird metastable state (eg, if you take 
an amount of water into a syringe, seal the tip, and then further pull the 
piston, you will get an empty region that is actually filled with water 
vapor, because having only liquid water filling that volume at that 
temperature is not thermodynamically stable). So, it is usually a good 
idea equilibrate in NPT, because the system finds its proper density at 
some temperature and pressure, whose regions of interest you usually know 
for the system you are studying. Once a good volume is found for some 
relevant P and T values, you can do the same for the energy: use NVT and 
let the system find a good energy for that T value, moving then to NVE. 
So, if for some reason you really want to do a simulation in the NVE 
ensemble, the suggested sequential procedure NPT  NVT  NVE sounds 
reasonable to me. Actually, I have a simpler suggestion: just run an NPT 
simulation, look at the 2D distribution histogram of V and E values, 
choose one representative snapshot that is in the central region of that 
distribution, and use that snapshot to run NVE. In any case, unless you 
have some experimental indication of the material density (N/V) *and* of 
the energy density (E/V) of your system, which would be extremelly 
unusual, you will have to follow some kind of approach similar to this. Of 
course, we may also ask why you think you need an NVE simulation, but that 
is an entirely different question...


Best,
Antonio


--
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--

On Wed, 5 Nov 2014, Johnny Lu wrote:


Well... I finally have to accept that I really need NVE. Some forum posts
suggested to run NPT, then NVT, and finally NVE.

For ideal gas law, if we fix the number of molecules and the temperature,
then fixing the volume at different value means changing the pressure.
So, I guess pressure is still meaningful, if I have to decide the correct
volume to run the NVT simulation.

For a large number of ideal gas molecules, averages of all thermodynamic
variables should be defined, regardless of which three of them is fixed or
which ensemble do I choose to represent the state.
For a small number of molecules, the distribution of the thermodynamic
variables should still have a definition (so, the average and the
distribution of pressure is still meaningful).

For a protein simulation, it still has a partition function, and equation
of state.

I mainly want to know how much error in pressure (deviation from 1 bar) is
acceptable.

On Wed, Nov 5, 2014 at 2:35 PM, Téletchéa Stéphane 
stephane.teletc...@univ-nantes.fr wrote:


Le 04/11/2014 18:00, Johnny Lu 

Re: [gmx-users] Use NVT to mimic NVE

[Antonio Baptista]

Then, run true NVE, for the reasons we already pointed.

On Sat, 11 Oct 2014, Johnny Lu wrote:


For dynamics with correct rate and correct fluctuation.

On Sat, Oct 11, 2014 at 8:09 AM, Dr. Vitaly Chaban vvcha...@gmail.com
wrote:


what is it needed for?



Dr. Vitaly V. Chaban

Виталий Витальевич ЧАБАН


On Fri, Oct 10, 2014 at 7:41 PM, Johnny Lu johnny.lu...@gmail.com wrote:
 Hi.

 Is it a good idea to mimic NVE by a NVT simulation with a large
temperature
 coupling time constant, to reduce the effect of the thermostat ?

 If I use V-rescale thermostat, what artifacts will the simulation get if
I
 use a large coupling time (say, 500 ps) and single precision gromacs ?

 Thank you.
 --
 Gromacs Users mailing list

 * Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.
--
Gromacs Users mailing list

* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.


--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Principal axes of inertia and g_principal

[Antonio Baptista]

Hi Nicola,

On Thu, 11 Sep 2014, nicola staffolani wrote:


​Dear GMX community,

regarding the program g_principal and acknowledged the bug reported here
http://permalink.gmane.org/gmane.science.biology.gromacs.user/66719​
​, I would like to understand what the meaning of the output of this
program is.
So, left aside moi.dat, the output files are axis1.dat, axis2.dat 
axis3.dat.
Supposing that I have already corrected the output files (or reinstalled
GROMACS having fixed the bug in the source file as explained always here
http://t110399.science-biology-gromacs-user.biotalk.us/g-principal-bug-or-very-bad-choice-of-filenames-t110399.html),
in each of the files, for each time instant there are three coordinates: x,
y and z: what do they mean? My interpretation is that if I connect the
point, whose coordinates are those printed out in the output files, with
the origin of the Cartesian axes, then I get a line coinciding with the
corresponding (principal) axis: is this interpretation correct?


Yes, that is correct. In other words, they are the coordinates of the 
point P=O+v, where O is the origin and v is the vector of unit norm that 
points along that principal axis. Of course, for graphical representation 
purposes, you would displace each of those 3 principal axis vectors to 
make them point away from the center of mass of your system, since the 
whole idea is to get a reference frame which is natural to that system 
(e.g., see Goldstein or any other book on rigid body mechanics).




To get the principal axes, I have to diagonalize the moment of inertia
symmetric 3×3 matrix; one way to do it is by solving the 3rd order equation
for the eigenvalues: is this the way by which GROMACS finds then the
eigenvectors of the moment of inertia matrix, from which it calculates the
aforementioned coordinates?


As I pointed in the links you mentioned, I didn't look very deep into the 
GROMACS code, but the vectors that g_principal writes (after the 
correction) are exactly the ones obtained from the diagonalization you 
describe. But check it yourself for a simple case. My rule of thumb: never 
trust an analysis tool until you check it yourself... :)


Best,
Antonio



Thank you in advance,

Nicola   ​

--
Nicola Staffolani PhD
Biophysics  Nanoscience Centre http://www.unitus.it/biophysics
Università della Tuscia
Largo dell'Università s.n.c., I-01100 Viterbo
email: n.staffol...@unitus.it
tel.: +39 0761 35 70 27; fax: +39 0761 35 71 36
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

--
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] isothermal-isometric ensemble on gromacs

[Antonio Baptista]

Well, the terminology is an old one, but not very common. An isometric 
ensemble is generally one in which the extension in the mechanical work 
term is kept constant. Thus, it may mean that the volume is kept constant 
in a volume/pressure system, that the area is kept constant in a 
area/surface tension system, that the length is kept constant in a 
length/linear tension system, etc. It is perhaps more common in the last 
case. Anyway, in the paper you mentioned it is probably just an unusual 
way of saying constant-volume, as Justin pointed.


Cheers,
Antonio


On Tue, 3 Jun 2014, Justin Lemkul wrote:




On 6/3/14, 10:05 AM, Carlos Navarrro Retamal wrote:


Dear gromacs users,
I’m trying to replicate the following paper 

http://www.ncbi.nlm.nih.gov/pubmed/?term=Desiccation+Induced+Structural+Alterations+in+a+66-Amino+Acid+Fragment+of+an+Anhydrobiotic+Nematode+Late+Embryogenesis+Abundant+(LEA)+Protein.
According to the methodology, they used gromacs 3.3.3 software and they 

performed several MD simulations using the isothermal-isometric ensemble.
This is the first time i’ve heard about this ensemble. Which kind of 
properties does this ensemble considerer? An even more importantly, how can i 
set up this ensemble in a MD simulations on gromacs?


From briefly looking at the methods, it sounds like they put all of their 
systems in boxes of the same size (hence isometric), and if those boxes were 
invariant (as is implied), then it's a constant volume ensemble, so it's NVT 
in 
this case.  The language is a bit nonstandard here, but that's what I'd 
assume 
they used.


-Justin

--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!


* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] g_principal -- bug or very bad choice of filenames

[Antonio Baptista]

On Fri, 7 Feb 2014, Justin Lemkul wrote:




On 2/7/14, 6:08 PM, Antonio Baptista wrote:

Dear all,

This is a follow-up to an old thread on g_principal, which continues (as of
version 4.6.5) to suffer from what I would call a bug or, at least, a very 
bad
and misleading choice of output file names. This message is primarily a 
further
warning to users, but I also hope that it promotes the solution of the 
problem

(I explicitly indicate below the small required code changes).

The program g_principal calculates the three principal axes of inertia for 
a

group of atoms and generates three output files with the default names
axis{1,2,3}.dat. Although the content of these files is not explained, 
their
names naturally suggest that axisN.dat contains the xyz components of the 
Nth
principal axis (presumably in the conventional major-to-minor order). 
Indeed,
each of these files contains (besides the time) 3 real values per frame, 
and

this interpretation also yields three vectors which are orthonormal, as one
would expect for the new frame defining the principal axes.

However, this natural interpretation is wrong, and not because of a 
different

axes order. As correctly identified by Chris Neale in the message included
below, the file axis1.dat contains the x components of the 1st, 2nd and 3rd
axes, the file axis2.dat their y components, and the file axis3.dat their z
components (he checked with VMD, we did it with in-house programs in C and
Octave). I think this a rather convoluted and extremely misleading way to 
output

the principal axes.

My guess is that this problem derives from a simple confusion between the 
form

of the orthogonal matrix obtained when solving the eigenvalue problem --
depending on the adopted algorithm, the vectors defining the new basis (the
principal axes) can be either the columns or the rows of this matrix. I 
didn't

look beyond gmx_principal.c, but the problem can be easily solved there by
replacing the current lines

 fprintf(axis1, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][XX], axes[YY][XX], axes[ZZ][XX]);
 fprintf(axis2, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][YY], axes[YY][YY], axes[ZZ][YY]);
 fprintf(axis3, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][ZZ], axes[YY][ZZ], axes[ZZ][ZZ]);

with

 fprintf(axis1, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][XX], axes[XX][YY], axes[XX][ZZ]);
 fprintf(axis2, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[YY][XX], axes[YY][YY], axes[YY][ZZ]);
 fprintf(axis3, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[ZZ][XX], axes[ZZ][YY], axes[ZZ][ZZ]);

With this change, the file names would make perfect sense, with axisN.dat 
simply
containing the components of the Nth principal axis. (Note that both the 
by-row
and by-column readings give orthonormal vectors, because this is an 
orthogonal

matrix.)

The file moi.dat also produced by g_principal is fine, containing the 
moments of
inertia along the principal axes in the proper order (lowest to highest, 
since
the inertia _around_ those axes increases as one goes from the major to the 
minor).


I believe that, as it stands now, g_principal is misleading many users into 
the
wrong interpretation of its output. Maybe some developer wants to have a 
look at

this issue and introduce my suggested fix.



Please file an issue on redmine.gromacs.org with all of the above 
information. Thanks for the thorough report!


More than a weak later, http://redmine.gromacs.org/login keeps telling me 
that Your account was created and is now pending administrator approval. 
Is this normal? Or maybe I missed some email message...


Best,
Antonio




-Justin

--
==

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!


* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.



--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] g_principal -- bug or very bad choice of filenames

[Antonio Baptista]

Szilárd already activated it a couple of hours ago (just noticed he 
replied to me only). The issue is already submitted also. Thanks anyway.


Best,
Antonio

On Wed, 19 Feb 2014, Mark Abraham wrote:


It shows as active to me...

Mark


On Tue, Feb 18, 2014 at 11:46 PM, Antonio Baptista bapti...@itqb.unl.ptwrote:


On Fri, 7 Feb 2014, Justin Lemkul wrote:




On 2/7/14, 6:08 PM, Antonio Baptista wrote:


Dear all,

This is a follow-up to an old thread on g_principal, which continues (as
of
version 4.6.5) to suffer from what I would call a bug or, at least, a
very bad
and misleading choice of output file names. This message is primarily a
further
warning to users, but I also hope that it promotes the solution of the
problem
(I explicitly indicate below the small required code changes).

The program g_principal calculates the three principal axes of inertia
for a
group of atoms and generates three output files with the default names
axis{1,2,3}.dat. Although the content of these files is not explained,
their
names naturally suggest that axisN.dat contains the xyz components of
the Nth
principal axis (presumably in the conventional major-to-minor order).
Indeed,
each of these files contains (besides the time) 3 real values per frame,
and
this interpretation also yields three vectors which are orthonormal, as
one
would expect for the new frame defining the principal axes.

However, this natural interpretation is wrong, and not because of a
different
axes order. As correctly identified by Chris Neale in the message
included
below, the file axis1.dat contains the x components of the 1st, 2nd and
3rd
axes, the file axis2.dat their y components, and the file axis3.dat
their z
components (he checked with VMD, we did it with in-house programs in C
and
Octave). I think this a rather convoluted and extremely misleading way
to output
the principal axes.

My guess is that this problem derives from a simple confusion between
the form
of the orthogonal matrix obtained when solving the eigenvalue problem --
depending on the adopted algorithm, the vectors defining the new basis
(the
principal axes) can be either the columns or the rows of this matrix. I
didn't
look beyond gmx_principal.c, but the problem can be easily solved there
by
replacing the current lines

 fprintf(axis1, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][XX], axes[YY][XX], axes[ZZ][XX]);
 fprintf(axis2, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][YY], axes[YY][YY], axes[ZZ][YY]);
 fprintf(axis3, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][ZZ], axes[YY][ZZ], axes[ZZ][ZZ]);

with

 fprintf(axis1, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[XX][XX], axes[XX][YY], axes[XX][ZZ]);
 fprintf(axis2, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[YY][XX], axes[YY][YY], axes[YY][ZZ]);
 fprintf(axis3, %15.10f %15.10f  %15.10f  %15.10f\n, t,
axes[ZZ][XX], axes[ZZ][YY], axes[ZZ][ZZ]);

With this change, the file names would make perfect sense, with
axisN.dat simply
containing the components of the Nth principal axis. (Note that both the
by-row
and by-column readings give orthonormal vectors, because this is an
orthogonal
matrix.)

The file moi.dat also produced by g_principal is fine, containing the
moments of
inertia along the principal axes in the proper order (lowest to highest,
since
the inertia _around_ those axes increases as one goes from the major to
the minor).

I believe that, as it stands now, g_principal is misleading many users
into the
wrong interpretation of its output. Maybe some developer wants to have a
look at
this issue and introduce my suggested fix.



Please file an issue on redmine.gromacs.org with all of the above
information. Thanks for the thorough report!



More than a weak later, http://redmine.gromacs.org/login keeps telling me
that Your account was created and is now pending administrator approval.
Is this normal? Or maybe I missed some email message...

Best,
Antonio





-Justin

--
==

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/
Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.

 --

Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/
Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe

[gmx-users] g_principal -- bug or very bad choice of filenames

[Antonio Baptista]

Dear all,

This is a follow-up to an old thread on g_principal, which continues 
(as of version 4.6.5) to suffer from what I would call a bug or, at 
least, a very bad and misleading choice of output file names. This 
message is primarily a further warning to users, but I also hope that it 
promotes the solution of the problem (I explicitly indicate below the 
small required code changes).


The program g_principal calculates the three principal axes of inertia 
for a group of atoms and generates three output files with the default 
names axis{1,2,3}.dat. Although the content of these files is not 
explained, their names naturally suggest that axisN.dat contains the xyz 
components of the Nth principal axis (presumably in the conventional 
major-to-minor order). Indeed, each of these files contains (besides the 
time) 3 real values per frame, and this interpretation also yields three 
vectors which are orthonormal, as one would expect for the new frame 
defining the principal axes.


However, this natural interpretation is wrong, and not because of a 
different axes order. As correctly identified by Chris Neale in the 
message included below, the file axis1.dat contains the x components of 
the 1st, 2nd and 3rd axes, the file axis2.dat their y components, and 
the file axis3.dat their z components (he checked with VMD, we did it 
with in-house programs in C and Octave). I think this a rather 
convoluted and extremely misleading way to output the principal axes.


My guess is that this problem derives from a simple confusion between 
the form of the orthogonal matrix obtained when solving the eigenvalue 
problem -- depending on the adopted algorithm, the vectors defining the 
new basis (the principal axes) can be either the columns or the rows of 
this matrix. I didn't look beyond gmx_principal.c, but the problem can 
be easily solved there by replacing the current lines


fprintf(axis1, %15.10f %15.10f  %15.10f  %15.10f\n, t, 
axes[XX][XX], axes[YY][XX], axes[ZZ][XX]);
fprintf(axis2, %15.10f %15.10f  %15.10f  %15.10f\n, t, 
axes[XX][YY], axes[YY][YY], axes[ZZ][YY]);
fprintf(axis3, %15.10f %15.10f  %15.10f  %15.10f\n, t, 
axes[XX][ZZ], axes[YY][ZZ], axes[ZZ][ZZ]);


with

fprintf(axis1, %15.10f %15.10f  %15.10f  %15.10f\n, t, 
axes[XX][XX], axes[XX][YY], axes[XX][ZZ]);
fprintf(axis2, %15.10f %15.10f  %15.10f  %15.10f\n, t, 
axes[YY][XX], axes[YY][YY], axes[YY][ZZ]);
fprintf(axis3, %15.10f %15.10f  %15.10f  %15.10f\n, t, 
axes[ZZ][XX], axes[ZZ][YY], axes[ZZ][ZZ]);


With this change, the file names would make perfect sense, with 
axisN.dat simply containing the components of the Nth principal axis. 
(Note that both the by-row and by-column readings give orthonormal 
vectors, because this is an orthogonal matrix.)


The file moi.dat also produced by g_principal is fine, containing the 
moments of inertia along the principal axes in the proper order (lowest 
to highest, since the inertia _around_ those axes increases as one goes 
from the major to the minor).


I believe that, as it stands now, g_principal is misleading many users 
into the wrong interpretation of its output. Maybe some developer wants 
to have a look at this issue and introduce my suggested fix.



Best,
Antonio

--
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de 
Lisboa

Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619 email: bapti...@itqb.unl.pt
fax:   +351-214411277 WWW:   http://www.itqb.unl.pt/~baptista
--


Nov 02, 2010; 4:46pm Chris Neale
Principal axis (g_principal)

Dear Sarah,

I just ran into the same thing so I'm posting this as a very late
response to your querry.

I think that g_principal prints the transpose of the principal axes
matrix, which is the exact opposite of what one would expect.

(1) First, I define a system that is longest in Z, second longest in
Y, and shortest in X:

cneale@iram:~$ cat a.gro
a
12
 1CL- CL1   0.000   0.000   0.000
 2CL- CL2   0.000   1.000   0.000
 3CL- CL3   0.000   2.000   0.000
 4CL- CL4   1.000   0.000   0.000
 5CL- CL5   1.000   1.000   0.000
 6CL- CL6   1.000   2.000   0.000
 7CL- CL7   0.000   0.000  10.000
 8CL- CL8   0.000   1.000  10.000
 9CL- CL9   0.000   2.000  10.000
10CL- CL   10   1.000   0.000  10.000
11CL- CL   11   1.000   1.000  10.000
12CL- CL   12   1.000   2.000  10.000
50 50 50
cneale@iram:~$ cat a.top
#include ffoplsaa.itp
#include ions.itp


[ system ]
; name
grid


[ molecules ]
; namenumber
CL-  12



(2) Next I calculate the principal axes:


$ echo 0 |g_principal -f a.gro -s a.tpr /dev/null 21; cat axis1.dat
axis2.dat axis3.dat
0.000.00 0.00 1.00