Re: [gmx-users] Calculating the average separation between two multi-atom groups

[Mark Abraham]

On 31/05/2012 11:24 AM, Andrew DeYoung wrote:

Hi,

I have a system in a slab geometry.  A surface exists at z = 0; many
hydrogens protrude from the surface, and these hydrogens are mostly (but not
precisely) immobile.  Above the surface, there is liquid, including the
anion BF4- (tetrahedral arrangement of fluorines around a central boron).
The liquid region is large, extending far from the surface, (far above the
surface, the liquid behaves like liquid in the bulk).

There is hydrogen-like bonding between the hydrogens protruding from the
surface and the fluorines in BF4-.  I would like to calculate the "average"
H-F distance, where H atoms protrude from the stationary surface and F atoms
exist on the BF4- ions.  But saying that I want to compute the "average" H-F
distance is very vague.  I can think of at least two possible, hopefully
reasonable, ways to formulate the problem:

(i) For the purpose of calculating the average H-F separation, only consider
fluorines on BF4- ions which are within a certain perpendicular distance z0
from the surface.  In other words, consider the BF4- ions which lie in the
region 0<  z<  z0 (where z0 is positive and very small compared to the z
dimension of the simulation box).  Then, using those BF4- ions, I calculate
the (time-averaged) H-F separation.

(ii) For the purpose of calculating the average H-F separation, only
consider fluorines when they are a certain small distance from any hydrogen.

Are (i) or (ii) these feasible?

For (i), I can think about using g_select to select BF4- ions which are a
distance of z0 or less from the surface at z = 0.  Maybe I would use a
selection like 'res_com of resname BF4 and z<  10' (where z0 = 10).  The
problem with this is that, I think, I would obtain an index file for each
simulation timestep.


Yes, you'd need to do some iteration over g_select and then another 
tool. GROMACS 5 will likely do that for you, but nirvana is some way off...



   So, I guess then if I have 200,000 simulation
timesteps, I would have to run g_bond 200,000 times!  (Or would g_dist be
appropriate here?)  Also, even my formulation in (i) is a little awkward;
fluorines at one edge of the xy dimension would be far from hydrogens
immobilized at the other side of the xy dimension, so I would get artifacts.


These tools *should* work correctly with PBC, but I would advise 
checking that they do, e.g. by making some very small subsets that you 
know cross boundaries and checking their stats agree with subsets that 
you know do not cross boundaries.



For (ii), it seems that g_hbond might be useful.  However, it does not seem
that fluorine is currently implemented as a hydrogen bond acceptor for use
in g_hbond.  I have never attempted to modify the Gromacs code and I am not
sure how easy this would be.  But if H-F is a hydrogen-like bond, then
(average) H-F bond length is what I am going after, I guess.

Do you know of any recipes with which to do this, or do you have any
suggestions?  Thanks so much!


You can match the trajectory to some .tpr that didn't create it. So if 
you generate a new topology that uses PO4 instead of BF4, while 
preserving the atom ordering over the whole system, you can fool g_hbond 
with a new .tpr without needing to touch any code. You don't even need 
sensible parameters for PO4, of course.


Alternatively, adding

((*top->atoms.atomname[n])[0] == 'F') ||

to the conditional that resembles

if ((bContact ||
 (((*top->atoms.atomname[n])[0] == 'O') ||
  (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N' &&
ISINGRP(datable[n])) {
datable[n] |= ACC; /* set the atom's acceptor flag in 
datable. */

add_acc(a,n,grp);
}

in search_acceptors() of src/gmxlib/gmx_hbond.c and recompiling is 
faster still - assuming you have no other atoms whose names start with 
"F". For safety, name the executable something that will remind you that 
it is not the standard one :-)


Mark
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] How GROMACS calculate the energy of hydrogen bond

[Acoot Brett]

Hi Mark,

It is confusing. As you know, for the same hydrogen bond in a protein, the 
related hydrogen bond angle and bond length can vary within a scope during the 
whole simulation process, however this small vibration of the hydrogen bond 
angle and length can lead to significant energy change, and correspondingly the 
energy of a hydrogen bond in simulation can be varied significantly. In 
comparison with hydrophobic effect, it would be too much is the energy of the 
hydrogen bond would be  not calculated  continuously.

Could you give some further clarification?

Cheers,

Acoot




 From: Mark Abraham 
To: Discussion list for GROMACS users  
Sent: Thursday, 31 May 2012 4:48 PM
Subject: Re: [gmx-users] How GROMACS calculate the energy of hydrogen bond
 

On 31/05/2012 4:42 PM, Acoot Brett wrote: 
Dear All, 
>The value of the energy of the hydrogen bond has relation with distance and 
>angle of the hydrogen bond related atoms. As for in the simulation process, 
>the distance and angle of the hydrogen bond related atoms may change 
>continuously. Will you please let me know based on which formula GROMACS 
>calculated the value of the energy of the hydrogen bonds?
There is no such formula used in MD force fields implemented in
GROMACS. The only non-bonded interactions are the ones you already
know about: electrostatics and VDW.
Observables like hydrogen bonds and the hydrophobic effect arise
from them.

Mark

-- 
gmx-users mailing list    gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] How GROMACS calculate the energy of hydrogen bond

[Mark Abraham]

On 31/05/2012 7:46 PM, Acoot Brett wrote:

Hi Mark,

It is confusing. As you know, for the same hydrogen bond in a protein, 
the related hydrogen bond angle and bond length can vary within a 
scope during the whole simulation process, however this small 
vibration of the hydrogen bond angle and length can lead to 
significant energy change, and correspondingly the energy of a 
hydrogen bond in simulation can be varied significantly. In comparison 
with hydrophobic effect, it would be too much is the energy of the 
hydrogen bond would be  not calculated  continuously.


It isn't, if the model physics isn't paramtrized to include it 
explicitly - which is the case for all the force fields in GROMACS.




Could you give some further clarification?


What are trying to do? Measuring "the strength of a hydrogen bond" 
requires you identify a state with and without it and a path between 
them over which you can integrate.


Mark



Cheers,

Acoot


*From:* Mark Abraham 
*To:* Discussion list for GROMACS users 
*Sent:* Thursday, 31 May 2012 4:48 PM
*Subject:* Re: [gmx-users] How GROMACS calculate the energy of 
hydrogen bond


On 31/05/2012 4:42 PM, Acoot Brett wrote:

Dear All,
The value of the energy of the hydrogen bond has relation with 
distance and angle of the hydrogen bond related atoms. As for in the 
simulation process, the distance and angle of the hydrogen bond 
related atoms may change continuously. Will you please let me know 
based on which formula GROMACS calculated the value of the energy of 
the hydrogen bonds?


There is no such formula used in MD force fields implemented in 
GROMACS. The only non-bonded interactions are the ones you already 
know about: electrostatics and VDW.
Observables like hydrogen bonds and the hydrophobic effect arise from 
them.


Mark

--
gmx-users mailing list gmx-users@gromacs.org 


http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!

Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org 
.

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





-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] OPLS specific problem in gromacs 4.5.4

[Ignacio Fernández Galván]

--- On Wed, 30/5/12, patrick wintrode  wrote:

> After creating a box using editconf and then solvating using genbox, I try
> to add ions and get the infamous "number of coordinates in .gro does not
> match number of coordinates in topology" error.

I've had this problem after solvating a PDB that already contained waters. The 
genbox command reported that N waters should be added, but since there were 
already n waters in the PDB, it added only N-n waters to the topology, while it 
should have added the whole lot of N molecules. Changing the number in the .top 
file solved this.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] PME nodes

[Ignacio Fernández Galván]

Hi all,

There must be something I don't fully understand, by running grompp on a 
system, I get this:

  Estimate for the relative computational load of the PME mesh part: 0.32

Good, that's approximately 1/3, or a 2:1 PP:PME ratio, which is the recommended 
value for a dodecahedral box. But then I run the dynamics with "mdrun_mpi -np 
8" (different cores in a single physical machine) and I get:

  Initializing Domain Decomposition on 8 nodes
  [...]
  Using 0 separate PME nodes

I would have expected at least 2 nodes (3:1, 0.25) to be used for PME, so 
there's obviously something wrong in my assumption.

Should I be looking somewhere in the output to find out why? Would it be better 
to try to get some dedicated PME node(s) (even in a single machine)?

Thanks,
Ignacio
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] PME nodes

[Peter C. Lai]

According to the manual, mdrun does not dedicate PME nodes unless -np > 11
You can manually specify dedicated PME nodes using -npme, but it is highly
system dependent on whether this will be faster on lowcore systems.

Also the estimate given by grompp may not be optimal during runtime. You'll
have to repeat runs on different node combinations or use g_tune_pme to  
discover the optimal PME:PP ratio for your particular system.

On 2012-05-31 04:11:15AM -0700, Ignacio Fernández Galván wrote:
> Hi all,
> 
> There must be something I don't fully understand, by running grompp on a 
> system, I get this:
> 
>   Estimate for the relative computational load of the PME mesh part: 0.32
> 
> Good, that's approximately 1/3, or a 2:1 PP:PME ratio, which is the 
> recommended value for a dodecahedral box. But then I run the dynamics with 
> "mdrun_mpi -np 8" (different cores in a single physical machine) and I get:
> 
>   Initializing Domain Decomposition on 8 nodes
>   [...]
>   Using 0 separate PME nodes
> 
> I would have expected at least 2 nodes (3:1, 0.25) to be used for PME, so 
> there's obviously something wrong in my assumption.
> 
> Should I be looking somewhere in the output to find out why? Would it be 
> better to try to get some dedicated PME node(s) (even in a single machine)?
> 
> Thanks,
> Ignacio
> -- 
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

-- 
==
Peter C. Lai| University of Alabama-Birmingham
Programmer/Analyst  | KAUL 752A
Genetics, Div. of Research  | 705 South 20th Street
p...@uab.edu| Birmingham AL 35294-4461
(205) 690-0808  |
==

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] PME nodes

[Mark Abraham]

On 31/05/2012 9:11 PM, Ignacio Fernández Galván wrote:

Hi all,

There must be something I don't fully understand, by running grompp on a 
system, I get this:

   Estimate for the relative computational load of the PME mesh part: 0.32

Good, that's approximately 1/3, or a 2:1 PP:PME ratio, which is the recommended value for 
a dodecahedral box. But then I run the dynamics with "mdrun_mpi -np 8" 
(different cores in a single physical machine) and I get:

   Initializing Domain Decomposition on 8 nodes
   [...]
   Using 0 separate PME nodes

I would have expected at least 2 nodes (3:1, 0.25) to be used for PME, so 
there's obviously something wrong in my assumption.

Should I be looking somewhere in the output to find out why? Would it be better 
to try to get some dedicated PME node(s) (even in a single machine)?


Generally mdrun does pretty well, given the constraints you've set for 
it. Here, you've implicitly let it choose (with mdrun -npme -1), and 
with fewer than a minimum number of nodes (10, in 4.5.5) it doesn't 
bother, since the book-keeping would be too costly. Otherwise, you can 
investigate the reasons for the choices mdrun made from the output in 
the .log file.


You can try mdrun -npme 2 or 3 if you like, but it's likely not faster 
or might even refuse to run. See manual 3.17, also.


Mark
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] (no subject)

[Subramaniam Boopathi]

how to assign charge and charge group number when new residue as being
added to the existing amino acid rtp file
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] (no subject)

[Mark Abraham]

On 31/05/2012 9:37 PM, Subramaniam Boopathi wrote:
how to assign charge and charge group number when new residue as being 
added to the existing amino acid rtp file


Read about the file format in the manual.

Mark
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] (no subject)

[Justin A. Lemkul]

On 5/31/12 7:37 AM, Subramaniam Boopathi wrote:

how to assign charge and charge group number when new residue as being added to
the existing amino acid rtp file



The details depend on the force field you're using.

http://www.gromacs.org/Documentation/How-tos/Parameterization

-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] About gomacs MPI installation

[vidhya sankar]

Dear justin
   Thank you for your previous reply
Now i am trying to run the parallel gromacs simulation (gromacs 4.5.5) first of 
all i have successfully installed debain package of  gromacs-openmpi 
for that i have configured and  compiled using the following command 
./configure  --enable-mpi --enable-double --prefix=/usr/local/mpigromacs  
--program-suffix=_mpi_d
then I issued the 
make 
make mdrun
make install
make install-mdrun
It compiles nicely and successfully
I have few doubt
1) should i need to install both  rpm of openmpi (coresponding to my OS) and 
Debain  of   gromacs-openmpi ? (compatible to my Linux OS)  Now i have both 
  Otherwise is it enough to install only Debain   gromacs-openmpi to install 
gromacs in paralleel

2) how to check Wheather Parallel installation of gromacs ? wheater is 
installed  parallely or not? 
I have installed on single  Hyper threading supported intel pentium i5 (dual 
core  processor,four thread)
it means i amusing thread based parallelism not mpi based parallelism (i know 
the performance may be poor, though later i will extend this to  clustering)

please tell me The Above my understandings is correct or not 
if ther are any  wrong things  in the compilation procedure Please give me some 
tips to rectify error?
I  am very grateful to your valuable reply


Also if the above procedure is correct
Can i run the gromacs parallel simulation using the following command
mdrun_mpi_d  -nt 8 -s topol.tpr

Don't i Need mpirun command ?


Thanks in Advance 

With regards 
S.vidhyasankar-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] About gomacs MPI installation

[Mark Abraham]

On 31/05/2012 10:17 PM, vidhya sankar wrote:

Dear justin
   Thank you for your previous reply
Now i am trying to run the parallel gromacs simulation (gromacs 4.5.5) 
first of all i have successfully installed debain package of  
gromacs-openmpi

for that i have configured and  compiled using the following command
./configure  --enable-mpi --enable-double 
--prefix=/usr/local/mpigromacs  --program-suffix=_mpi_d

then I issued the
make
make mdrun


As you will see in the 
http://www.gromacs.org/Documentation/Installation_Instructions, there is 
no need to do both "make mdrun" and "make" because the former does part 
of latter. Likewise for "make install-mdrun" and "make install". However 
it doesn't hurt.



make install
make install-mdrun
It compiles nicely and successfully
I have few doubt
1) should i need to install both  rpm of openmpi (coresponding to my 
OS) and Debain  of   gromacs-openmpi ? (compatible to my Linux OS) Now 
i have both


You've installed GROMACS from source already. You have no need of a 
binary distribution in an .rpm or .deb. However that binary probably is 
linked to FFTW3, which your source installation was not, so PME will 
likely be much faster. You should probably get rid of the source 
installation, and do it again following 
http://www.gromacs.org/Documentation/Installation_Instructions, or use 
the binary distribution.


  Otherwise is it enough to install only Debain   gromacs-openmpi to 
install gromacs in paralleel


Probably it is.



2) how to check Wheather Parallel installation of gromacs ? wheater is 
installed  parallely or not?


Something named gromacs-openmpi, or successfully configured --with-mpi 
is very likely to run in parallel with MPI. That's the point of the name!


I have installed on single  Hyper threading supported intel pentium i5 
(dual core  processor,four thread)
it means i amusing thread based parallelism not mpi based parallelism 
(i know the performance may be poor, though later i will extend this 
to  clustering)


No. See http://www.gromacs.org/Documentation/Installation_Instructions. 
Hyper-threading is useless for GROMACS, and you have configured with 
MPI, which prevents attempting to use threading of any sort.




please tell me The Above my understandings is correct or not
if ther are any  wrong things  in the compilation procedure Please 
give me some tips to rectify error?

I  am very grateful to your valuable reply


Also if the above procedure is correct
Can i run the gromacs parallel simulation using the following command
mdrun_mpi_d  -nt 8 -s topol.tpr


mdrun -nt is not available with an MPI-configured mdrun, and likely 
exits with an error.


Mark



Don't i Need mpirun command ?


Thanks in Advance

With regards
S.vidhyasankar






-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] How GROMACS calculate the energy of hydrogen bond

[lloyd riggs]

Dear All,

I have no clue what specifically you are trying, but I feal bad for all the 
physicist and quantum chemist whom have provided the software and continued to 
develop it.

Scanning in my free time, it seems a large amount of confusion on what people 
are trying to do stems from differences in what is taught textbook wise for 
things.

For instance a hydrogen bond to a physicist is an integration over space in 3 
dimensions including time and probabilities of occupied spaces (atom position 
variabilities reflected even more in proteins, ie the necessity of multiple MD 
runs with different starting conformations), Vs.  an organic chemist whom has 
cut offs, ie angles between two points and set distances between two atoms 
which generally reflect the means of calculated chemical energies within a 
range (say 80-90% which represent means, but usually from raw small molecules 
as determinants), Vs. Biologist whom have tables which either use a set 
distance and angle and little account of variability over time (ie a hydrogen 
bond equals 1.4 kCal/mol reflecting the absolute mean), conformations in amino 
acids, etc...

I think with gromacs it is very precise, as even the smallest energies between 
two interacting atoms is taken into account with accuracy reflected by the 
force fields used, and how they were derived.

Good luck, your going to start seeing more and more a flood of biologist.

Stephan Watkins

 Original-Nachricht 
> Datum: Thu, 31 May 2012 19:54:04 +1000
> Von: Mark Abraham 
> An: Discussion list for GROMACS users 
> Betreff: Re: [gmx-users] How GROMACS calculate the energy of hydrogen bond

> On 31/05/2012 7:46 PM, Acoot Brett wrote:
> > Hi Mark,
> >
> > It is confusing. As you know, for the same hydrogen bond in a protein, 
> > the related hydrogen bond angle and bond length can vary within a 
> > scope during the whole simulation process, however this small 
> > vibration of the hydrogen bond angle and length can lead to 
> > significant energy change, and correspondingly the energy of a 
> > hydrogen bond in simulation can be varied significantly. In comparison 
> > with hydrophobic effect, it would be too much is the energy of the 
> > hydrogen bond would be  not calculated  continuously.
> 
> It isn't, if the model physics isn't paramtrized to include it 
> explicitly - which is the case for all the force fields in GROMACS.
> 
> >
> > Could you give some further clarification?
> 
> What are trying to do? Measuring "the strength of a hydrogen bond" 
> requires you identify a state with and without it and a path between 
> them over which you can integrate.
> 
> Mark
> 
> >
> > Cheers,
> >
> > Acoot
> >
> > 
> > *From:* Mark Abraham 
> > *To:* Discussion list for GROMACS users 
> > *Sent:* Thursday, 31 May 2012 4:48 PM
> > *Subject:* Re: [gmx-users] How GROMACS calculate the energy of 
> > hydrogen bond
> >
> > On 31/05/2012 4:42 PM, Acoot Brett wrote:
> >> Dear All,
> >> The value of the energy of the hydrogen bond has relation with 
> >> distance and angle of the hydrogen bond related atoms. As for in the 
> >> simulation process, the distance and angle of the hydrogen bond 
> >> related atoms may change continuously. Will you please let me know 
> >> based on which formula GROMACS calculated the value of the energy of 
> >> the hydrogen bonds?
> >
> > There is no such formula used in MD force fields implemented in 
> > GROMACS. The only non-bonded interactions are the ones you already 
> > know about: electrostatics and VDW.
> > Observables like hydrogen bonds and the hydrophobic effect arise from 
> > them.
> >
> > Mark
> >
> > -- 
> > gmx-users mailing list gmx-users@gromacs.org 
> > 
> > http://lists.gromacs.org/mailman/listinfo/gmx-users
> > Please search the archive at 
> > http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> > Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-requ...@gromacs.org 
> > .
> > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> >
> >
> 

-- 
NEU: FreePhone 3-fach-Flat mit kostenlosem Smartphone!  

Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

RE: [gmx-users] REMD question

[Nathalia Garces]

Michael,

Thank you for your answer. On the other hand, I´m not implementing, I´m
using REMD.. I miss wrote it.

Nathalia 

-Mensaje original-
De: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] En
nombre de Michael Shirts
Enviado el: lunes, 28 de mayo de 2012 08:09 p.m.
Para: Discussion list for GROMACS users
Asunto: Re: [gmx-users] REMD question

Gromacs already supports replica exchange -- what particularly are you
implementing?

Equilibration of pressure is always a good idea -- even if you are running
NVT simulations, you want to get them to be at the equilibrium volume for
your system and temperature choice, which will require equilibration at
constant pressure.

On Mon, May 28, 2012 at 4:37 PM, Nathalia Garces 
wrote:
> Dear Gromacs Users,
>
>
>
> We are implementing REMD method in Gromacs in protein folding, in your 
> web page you give some steps that don´t mention any step about NPT 
> stabilization.  This step is necessary to run REMD simulations?
>
>
>
> Thank you in advance,
>
>
>
> Nathalia
>
>
>
>
> --
> gmx-users mailing list    gmx-users@gromacs.org 
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the www interface
or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] About gomacs MPI installation

[Justin A. Lemkul]

On 5/31/12 8:17 AM, vidhya sankar wrote:

Dear justin
Thank you for your previous reply
Now i am trying to run the parallel gromacs simulation (gromacs 4.5.5) first of
all i have successfully installed debain package of gromacs-openmpi
for that i have configured and compiled using the following command
./configure --enable-mpi --enable-double --prefix=/usr/local/mpigromacs
--program-suffix=_mpi_d
then I issued the
make
make mdrun
make install
make install-mdrun
It compiles nicely and successfully


If you have installed a Debian package, why are you re-compiling from source? 
Now you may have duplicate executables on your system.



I have few doubt
1) should i need to install both rpm of openmpi (coresponding to my OS) and
Debain of gromacs-openmpi ? (compatible to my Linux OS) Now i have both
Otherwise is it enough to install only Debain gromacs-openmpi to install gromacs
in paralleel



I would assume the Debian package contains everything you need; refer to its 
documentation.  Either you need a pre-built package or you need to compile from 
source, not both, unless the gromacs-openmpi package contains only mdrun, in 
which case you will need the other tools.



2) how to check Wheather Parallel installation of gromacs ? wheater is installed
parallely or not?


If you have mdrun_mpi (or mdrun_mpi_d from source) and it executes in parallel 
via mpirun, then it is installed and working.



I have installed on single Hyper threading supported intel pentium i5 (dual core
processor,four thread)
it means i amusing thread based parallelism not mpi based parallelism (i know
the performance may be poor, though later i will extend this to clustering)


Use of threading and MPI are mutually exclusive.  If you invoke --enable-mpi, 
you are not using threading.  If your architecture supports threading, there is 
no need to compile with MPI support.  Configuration will detect necessary 
options and the resulting mdrun executable will have an argument called -nt, 
which allows you to execute the process using threads.




please tell me The Above my understandings is correct or not
if ther are any wrong things in the compilation procedure Please give me some
tips to rectify error?
I am very grateful to your valuable reply


Also if the above procedure is correct
Can i run the gromacs parallel simulation using the following command
mdrun_mpi_d -nt 8 -s topol.tpr

Don't i Need mpirun command ?



The above command is incorrect.  If compiled with MPI, there should be no -nt 
flag.  If you invoke MPI support when configuring, you need mpirun:


mpirun -np X mdrun_mpi -s topol.tpr

All of these issues are discussed in the manual and on the Gromacs website.

-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] surface tension, Long range LJ correction, TIP4P/2005

[MD]

Hi All,

I really need to know how to apply long range LJ correction to calculate 
surface tension of TIP4P/2005 water. I can get 65 dyn even i use vdw cut-off = 
1.4 nm, but from the reference people can get 69 dyn.
I included the LJ  Long range LJ correction using the following .mdp,
please note that i used: DispCorr  = EnerPres, which means i included the long 
range LJ correction for energy and pressure, but why i can get only 65 dyn for 
surface tension. I can see Disper. corr. in the .log file, but i saw
Table routines are used for coulomb: TRUE
Table routines are used for vdw: FALSE
I spend two weeks on this ,but still failed to know the reason, because only 
myself do MD in my department. Can anyone help??

title= Yo
cpp  = /usr/bin/cpp
include  =
define   =
integrator   = md
tinit= 0
dt   = 0.001
nsteps   = 500
init_step= 0
comm-mode= Linear
nstcomm  = 1
comm-grps=
bd-fric  = 0
ld-seed  = 1993
niter= 20
nstxout  = 5000
nstvout  = 8000
nstfout  = 8000
nstcheckpoint= 1000
nstlog   = 5000
nstenergy= 5000
nstxtcout= 500
xtc-precision= 1000
xtc-grps =
energygrps   =
nstlist  = 5
ns_type  = grid
pbc  = xyz
rlist= 1.4
domain-decomposition = no
coulombtype  = PME
rcoulomb-switch  = 0
rcoulomb = 1.4
epsilon-r= 1
vdw-type = Cut-off
rvdw-switch  = 0
rvdw = 1.4
DispCorr = EnerPres
table-extension  = 1
fourierspacing   = 0.12
fourier_nx   = 0
fourier_ny   = 0
fourier_nz   = 0
pme_order= 4
ewald_rtol   = 1e-05
ewald_geometry   = 3d
epsilon_surface  = 0
optimize_fft = no
gb_algorithm = Still
nstgbradii   = 1
rgbradii = 2
gb_saltconc  = 0
implicit_solvent = No
Tcoupl   = v-rescale
tc-grps  = System
tau_t= 0.1
ref_t= 300
Pcoupl   = no
Pcoupltype   = isotropic
tau_p= 1
compressibility  = 4.5e-5
ref_p= 1.0
andersen_seed= 815131
annealing= no
annealing_npoints=
annealing_time   =
annealing_temp   =
gen_vel  = yes
gen_temp = 300
gen_seed = 1993
constraints  = none
constraint-algorithm = Lincs
unconstrained-start  = no
Shake-SOR= no
shake-tol= 1e-04
lincs-order  = 4
lincs-iter   = 1
lincs-warnangle  = 30
morse= no











-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] surface tension, Long range LJ correction, TIP4P/2005

[Mark Abraham]

On 1/06/2012 12:18 AM, MD wrote:

Hi All,

I really need to know how to apply long range LJ correction to 
calculate surface tension of TIP4P/2005 water. I can get 65 dyn even i 
use vdw cut-off = 1.4 nm, but from the reference people can get 69 dyn.

I included the LJ  Long range LJ correction using the following .mdp,
please note that i used: DispCorr  = EnerPres, which means i included 
the long range LJ correction for energy and pressure, but why i can 
get only 65 dyn for surface tension. I can see Disper. corr. in the 
.log file, but i saw

Table routines are used for coulomb: TRUE
Table routines are used for vdw: FALSE


That's normal for your choices of coulombtype and vdw-type.

I spend two weeks on this ,but still failed to know the reason, 
because only myself do MD in my department. Can anyone help??


Yes. You can, by finding that 69 dyn reference and reading it :-) 
Nothing else is worthwhile. I'm now going to stop giving you the same 
advice each time you ask the same question.


Mark



title= Yo
cpp& nbsp; = /usr/bin/cpp
include  =
define   =
integrator   = md
tinit= 0
dt   = 0.001
nsteps   = 500
init_step= 0< br>comm-mode= Linear
nstcomm  = 1
comm-grps=
bd-fric  = 0
ld-seed  = 1993
niter= 20
nstxout  = 5000
nstvout  = 8000
nstfout &nb sp;   = 8000
nstcheckpoint= 1000
nstlog   = 5000
nstenergy= 5000
nstxtcout= 500
xtc-precision= 1000
xtc-grps =
energygrps   =
nstlist  = 5
ns_type & nbsp;= grid
pbc  = xyz
rlist= 1.4
domain-decomposition = no
coulombtype  = PME
rcoulomb-switch  = 0
rcoulomb = 1.4
epsilon-r= 1
vdw-type = Cut-off
rvdw-switch  = 0*rvdw = 1.4
DispCorr = EnerPres
table-extension  = 1
fourierspacing   = 0.12
fourier_nx   = 0
fourier_ny   = 0
fourier_nz   = 0
pme_order= 4
ewald_rtol   = 1e-05
ewald_geometry &nb sp;= 3d
epsilon_surface  = 0
optimize_fft = no
gb_algorithm = Still
nstgbradii   = 1
rgbradii = 2
gb_saltconc  = 0
implicit_solvent = No
Tcoupl   = v-rescale
tc-grps  = System
tau_t &nbs p;   = 0.1
ref_t= 300
Pcoupl   = no
Pcoupltype   = isotropic
tau_p= 1
compressibility  = 4.5e-5
ref_p= 1.0
andersen_seed= 815131
annealing= no
annealing_npoints =
annealing_time   =
annealing_temp   =
gen_vel  = yes
gen_temp = 300
gen_seed = 1993
constraints  = none
constraint-algorithm = Lincs
unconstrained-start  = no
Shake-SOR= no
shake-tol= 1e-04< br>lincs-order  = 4
lincs-iter   = 1
lincs-warnangle  = 30
morse= no











*
*



*


-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: Calculating the average separation between two multi-atom groups

[Andrew DeYoung]

Mark,

Thank you SO much for your help.  I will try these suggestions.  You're
awesome.

Andrew

> Hi,
>
> I have a system in a slab geometry.  A surface exists at z = 0; many
> hydrogens protrude from the surface, and these hydrogens are mostly (but
not
> precisely) immobile.  Above the surface, there is liquid, including the
> anion BF4- (tetrahedral arrangement of fluorines around a central boron).
> The liquid region is large, extending far from the surface, (far above the
> surface, the liquid behaves like liquid in the bulk).
>
> There is hydrogen-like bonding between the hydrogens protruding from the
> surface and the fluorines in BF4-.  I would like to calculate the
"average"
> H-F distance, where H atoms protrude from the stationary surface and F
atoms
> exist on the BF4- ions.  But saying that I want to compute the "average"
H-F
> distance is very vague.  I can think of at least two possible, hopefully
> reasonable, ways to formulate the problem:
>
> (i) For the purpose of calculating the average H-F separation, only
consider
> fluorines on BF4- ions which are within a certain perpendicular distance
z0
> from the surface.  In other words, consider the BF4- ions which lie in the
> region 0<  z<  z0 (where z0 is positive and very small compared to the z
> dimension of the simulation box).  Then, using those BF4- ions, I
calculate
> the (time-averaged) H-F separation.
>
> (ii) For the purpose of calculating the average H-F separation, only
> consider fluorines when they are a certain small distance from any
hydrogen.
>
> Are (i) or (ii) these feasible?
>
> For (i), I can think about using g_select to select BF4- ions which are a
> distance of z0 or less from the surface at z = 0.  Maybe I would use a
> selection like 'res_com of resname BF4 and z<  10' (where z0 = 10).  The
> problem with this is that, I think, I would obtain an index file for each
> simulation timestep.

Yes, you'd need to do some iteration over g_select and then another 
tool. GROMACS 5 will likely do that for you, but nirvana is some way off...

>So, I guess then if I have 200,000 simulation
> timesteps, I would have to run g_bond 200,000 times!  (Or would g_dist be
> appropriate here?)  Also, even my formulation in (i) is a little awkward;
> fluorines at one edge of the xy dimension would be far from hydrogens
> immobilized at the other side of the xy dimension, so I would get
artifacts.

These tools *should* work correctly with PBC, but I would advise 
checking that they do, e.g. by making some very small subsets that you 
know cross boundaries and checking their stats agree with subsets that 
you know do not cross boundaries.

> For (ii), it seems that g_hbond might be useful.  However, it does not
seem
> that fluorine is currently implemented as a hydrogen bond acceptor for use
> in g_hbond.  I have never attempted to modify the Gromacs code and I am
not
> sure how easy this would be.  But if H-F is a hydrogen-like bond, then
> (average) H-F bond length is what I am going after, I guess.
>
> Do you know of any recipes with which to do this, or do you have any
> suggestions?  Thanks so much!

You can match the trajectory to some .tpr that didn't create it. So if 
you generate a new topology that uses PO4 instead of BF4, while 
preserving the atom ordering over the whole system, you can fool g_hbond 
with a new .tpr without needing to touch any code. You don't even need 
sensible parameters for PO4, of course.

Alternatively, adding

((*top->atoms.atomname[n])[0] == 'F') ||

to the conditional that resembles

 if ((bContact ||
  (((*top->atoms.atomname[n])[0] == 'O') ||
   (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N' &&
 ISINGRP(datable[n])) {
 datable[n] |= ACC; /* set the atom's acceptor flag in 
datable. */
 add_acc(a,n,grp);
 }

in search_acceptors() of src/gmxlib/gmx_hbond.c and recompiling is 
faster still - assuming you have no other atoms whose names start with 
"F". For safety, name the executable something that will remind you that 
it is not the standard one :-)

Mark

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re : How GROMACS calculate the energy of hydrogen bond

[ABEL Stephane 175950]

Hi, 

Here, my 2 cents worth. You can also estimate (roughly !!) the HB energy 
between two groups (say NH---CO) by using the Kabsch and Sander function  
described  in Kabsch, W.; Sander, C.  Biopolymers 1983, 22, 2577−2637).

Quoting:" E = qlq2(1/r(0N) + l/r(CH) - l/r(OH) - l/r(CN))*f with q1 = 0.42e and 
q 2 = 0.20e, e being the unit electron charge and r(AB) the interatomic 
distance from A to B. In chemical units, r is in angstroms,
the dimensional factor f = 332, and E is in kcal/mol. A good H bond has about 
-3 kcal/mol binding energy. We choose a generous cutoff to allow for bifurcated 
H bonds and errors in coordinates and assign an H bond
between C=O of residue i and N-H of residue j if E is less than the cutoff, 
i.e., “Hbond(ij)=: [E < -0.5kcal/mole].”

To obtain  HB energy value E, you need only the distance between the donnor and 
acceptor groups. 

HTH 

Stephane 




--

Message: 2
Date: Thu, 31 May 2012 15:00:05 +0200
From: "lloyd riggs" 
Subject: Re: [gmx-users] How GROMACS calculate the energy of hydrogen
bond
To: Discussion list for GROMACS users 
Message-ID: <20120531130005.302...@gmx.net>
Content-Type: text/plain; charset="utf-8"

Dear All,

I have no clue what specifically you are trying, but I feal bad for all the 
physicist and quantum chemist whom have provided the software and continued to 
develop it.

Scanning in my free time, it seems a large amount of confusion on what people 
are trying to do stems from differences in what is taught textbook wise for 
things.

For instance a hydrogen bond to a physicist is an integration over space in 3 
dimensions including time and probabilities of occupied spaces (atom position 
variabilities reflected even more in proteins, ie the necessity of multiple MD 
runs with different starting conformations), Vs.  an organic chemist whom has 
cut offs, ie angles between two points and set distances between two atoms 
which generally reflect the means of calculated chemical energies within a 
range (say 80-90% which represent means, but usually from raw small molecules 
as determinants), Vs. Biologist whom have tables which either use a set 
distance and angle and little account of variability over time (ie a hydrogen 
bond equals 1.4 kCal/mol reflecting the absolute mean), conformations in amino 
acids, etc...

I think with gromacs it is very precise, as even the smallest energies between 
two interacting atoms is taken into account with accuracy reflected by the 
force fields used, and how they were derived.

Good luck, your going to start seeing more and more a flood of biologist.

Stephan Watkins

 Original-Nachricht 
> Datum: Thu, 31 May 2012 19:54:04 +1000
> Von: Mark Abraham 
> An: Discussion list for GROMACS users 
> Betreff: Re: [gmx-users] How GROMACS calculate the energy of hydrogen bond

> On 31/05/2012 7:46 PM, Acoot Brett wrote:
> > Hi Mark,
> >
> > It is confusing. As you know, for the same hydrogen bond in a protein,
> > the related hydrogen bond angle and bond length can vary within a
> > scope during the whole simulation process, however this small
> > vibration of the hydrogen bond angle and length can lead to
> > significant energy change, and correspondingly the energy of a
> > hydrogen bond in simulation can be varied significantly. In comparison
> > with hydrophobic effect, it would be too much is the energy of the
> > hydrogen bond would be  not calculated  continuously.
>
> It isn't, if the model physics isn't paramtrized to include it
> explicitly - which is the case for all the force fields in GROMACS.
>
> >
> > Could you give some further clarification?
>
> What are trying to do? Measuring "the strength of a hydrogen bond"
> requires you identify a state with and without it and a path between
> them over which you can integrate.
>
> Mark
>
> >
> > Cheers,
> >
> > Acoot
> >
> > 
> > *From:* Mark Abraham 
> > *To:* Discussion list for GROMACS users 
> > *Sent:* Thursday, 31 May 2012 4:48 PM
> > *Subject:* Re: [gmx-users] How GROMACS calculate the energy of
> > hydrogen bond
> >
> > On 31/05/2012 4:42 PM, Acoot Brett wrote:
> >> Dear All,
> >> The value of the energy of the hydrogen bond has relation with
> >> distance and angle of the hydrogen bond related atoms. As for in the
> >> simulation process, the distance and angle of the hydrogen bond
> >> related atoms may change continuously. Will you please let me know
> >> based on which formula GROMACS calculated the value of the energy of
> >> the hydrogen bonds?
> >
> > There is no such formula used in MD force fields implemented in
> > GROMACS. The only non-bonded interactions are the ones you already
> > know about: electrostatics and VDW.
> > Observables like hydrogen bonds and the hydrophobic effect arise from
> > them.
> >
> > Mark
> >
> > --
> > gmx-users mailing list gmx-users@gromacs.org
> > 
> > http://lists.gro

[gmx-users] About usage of node in mpi calculation

[vidhya sankar]

Dear justin ,
  Very very thanks for your previous patience reply
When i run the mpi calculation using more than one node  (each node have 16 
processor)  which option do i need to use in the following command  may i use 
-npme option?

mpirun  -np 19  mdrun_mpi_d  -s toplo.tpr

Thanks in Advance
S.Vidhyasankar
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] About usage of node in mpi calculation

[Justin A. Lemkul]

On 5/31/12 12:00 PM, vidhya sankar wrote:

Dear justin ,
Very very thanks for your previous patience reply
When i run the mpi calculation using more than one node (each node have 16
processor) which option do i need to use in the following command may i use
-npme option?
mpirun -np 19 mdrun_mpi_d -s toplo.tpr



The above command will not work on 16 processors, since it specifies 19. 
Assuming that is a typo, there's nothing necessarily wrong with that command at 
all.  The -npme option is only necessary when you need to override whatever 
mdrun detects as being optimal (which is usually fine and does not require 
intervention).


-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] (no subject)

[Marc Charendoff]

http://livehistorytours.com/wp-content/uploads/wapple-architect/images/postHeader/loveit.php?public138.jpg-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] boundaries in PMF

[Rebeca García Fandiño]

Hi,
I am trying to calculate a PMF for an ion along a channel. Everything went OK, 
but when I used g_wham I got a profile with strange dimensions for the x-axis.
What are the boundaries g_wham is using for calculating the units of x-axis?

I have used:
g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist 
file_histo_output.xvg -unit kCal -cycl yes

(version 4.0.7)

and the units I got in the x-axis are determined by the boundaries:

"Determined boundaries to 0.35 to 2.603290 "

Which are these units? nm?
 
The z coordinate for my ion explores at least 5 nm!!

I am a bit confused about that. 

Any help will be appreciated.
Thanks a lot for your help in advance!

Best wishes,

Dr.Rebeca Garcia
Santiago de Compostela University
Spain
  -- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] boundaries in PMF

[Justin A. Lemkul]

On 5/31/12 1:20 PM, Rebeca García Fandiño wrote:

Hi,
I am trying to calculate a PMF for an ion along a channel. Everything went OK,
but when I used g_wham I got a profile with strange dimensions for the x-axis.
What are the boundaries g_wham is using for calculating the units of x-axis?



The values are based on the restraint distances along the reaction coordinate.


I have used:
g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist
file_histo_output.xvg -unit kCal -cycl yes

(version 4.0.7)

and the units I got in the x-axis are determined by the boundaries:

"Determined boundaries to 0.35 to 2.603290 "

Which are these units? nm?



Yes.


The z coordinate for my ion explores at least 5 nm!!

I am a bit confused about that.



The exact output depends on how you set up the umbrella sampling (in the .mdp 
file).  The range of values corresponds to whatever the distances are that are 
sampled in the various windows.  Perhaps there is a sign issue here?  Do you 
have some restraints that are at negative displacement and others at positive? 
Did you set up the .mdp files properly to account for this behavior?


-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

RE: [gmx-users] boundaries in PMF

[Rebeca García Fandiño]

Thanks a lot for your quick answer.
The mdp file I have used is copied below. What is strange is that when I look 
at the *gro files for the different windows (100 windows in total), i. e:

window 1:  8704Na  Na56458   5.134   5.085   5.824
window 50:  8704Na  Na56458   5.053   5.081   3.459
window 100:  8704Na  Na56458   4.990   5.042   0.951

you can see that the z-coordinate goes from 0.951 to 5.824 nm 

I should have a total distance in the x-axis of about 5 nm, and however, the 
boundaries calculated by g_wham are

 "Determined boundaries to 0.35 to 2.603290 "

Can you see anything in the mdp file which is causing this behaviour...?

Thanks again for your help!




MDP FILE USED:
title   = Umbrella pulling simulation
define  = -DPOSRES_B
define  =
; Run parameters
integrator  = md
dt  = 0.002
tinit   = 0
nsteps  = 50   ; 1 ns
nstcomm = 10
; Output parameters
nstxout = 5000 ; every 10 ps
nstvout = 5000
nstfout = 5000
nstxtcout   = 5000  ; every 10 ps
nstenergy   = 5000
; Bond parameters
constraint_algorithm= lincs
constraints = all-bonds
continuation= yes
; Single-range cutoff scheme
nstlist = 5
ns_type = grid
rlist   = 1.4
rcoulomb= 1.4
rvdw= 1.4
; PME electrostatics parameters
coulombtype = PME
fourierspacing  = 0.12
fourier_nx  = 0
fourier_ny  = 0
fourier_nz  = 0
pme_order   = 4
ewald_rtol  = 1e-5
optimize_fft= yes
;Berendsen temperature coupling is on
Tcoupl  = V-rescale
tau_t   = 0.1 0.10.1
tc-grps = MOL dop WAT_Cl_Na
ref_t   = 300 300300
; Pressure coupling
Pcoupl   = Parrinello-Rahman
Pcoupltype   = Semiisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p= 1.0 1.0
compressibility  = 4.6E-5 4.6E-5
ref-p= 1.0 1.0
; Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Long-range dispersion correction
DispCorr= EnerPres
; Pull code
pull= umbrella
pull_geometry   = distance
pull_dim= N N Y
pull_start  = yes
pull_ngroups= 1
pull_group0 = MOL
pull_group1 = Na
pull_init1  = 0
pull_rate1  = 0.0
pull_k1 = 1000  ; kJ mol^-1 nm^-2
pull_nstxout= 1000  ; every 2 ps
pull_nstfout= 1000  ; every 2 ps









> Date: Thu, 31 May 2012 13:25:26 -0400
> From: jalem...@vt.edu
> To: gmx-users@gromacs.org
> Subject: Re: [gmx-users] boundaries in PMF
> 
> 
> 
> On 5/31/12 1:20 PM, Rebeca García Fandiño wrote:
> > Hi,
> > I am trying to calculate a PMF for an ion along a channel. Everything went 
> > OK,
> > but when I used g_wham I got a profile with strange dimensions for the 
> > x-axis.
> > What are the boundaries g_wham is using for calculating the units of x-axis?
> >
> 
> The values are based on the restraint distances along the reaction coordinate.
> 
> > I have used:
> > g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist
> > file_histo_output.xvg -unit kCal -cycl yes
> >
> > (version 4.0.7)
> >
> > and the units I got in the x-axis are determined by the boundaries:
> >
> > "Determined boundaries to 0.35 to 2.603290 "
> >
> > Which are these units? nm?
> >
> 
> Yes.
> 
> > The z coordinate for my ion explores at least 5 nm!!
> >
> > I am a bit confused about that.
> >
> 
> The exact output depends on how you set up the umbrella sampling (in the .mdp 
> file).  The range of values corresponds to whatever the distances are that 
> are 
> sampled in the various windows.  Perhaps there is a sign issue here?  Do you 
> have some restraints that are at negative displacement and others at 
> positive? 
> Did you set up the .mdp files properly to account for this behavior?
> 
> -Justin
> 
> -- 
> 
> 
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> -- 
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
  -- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read ht

[gmx-users] MD Training

[edward.de...@gmail.com]

Dear all,
does anyone know of any MD or gromacs training course in the coming months? 
All the best
 g-

--
Sent from my Nokia phone
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] boundaries in PMF

[Justin A. Lemkul]

On 5/31/12 1:38 PM, Rebeca García Fandiño wrote:

Thanks a lot for your quick answer.
The mdp file I have used is copied below. What is strange is that when I look at
the *gro files for the different windows (100 windows in total), i. e:

window 1: 8704Na Na56458 5.134 5.085 5.824
window 50: 8704Na Na56458 5.053 5.081 3.459
window 100: 8704Na Na56458 4.990 5.042 0.951

you can see that the z-coordinate goes from 0.951 to 5.824 nm

I should have a total distance in the x-axis of about 5 nm, and however, the
boundaries calculated by g_wham are

"Determined boundaries to 0.35 to 2.603290 "

Can you see anything in the mdp file which is causing this behaviour...?



No, but you need to be careful how you're interpreting what you get.  The 
z-coordinate of your ion is not what's relevant; its displacement along the 
z-axis from the reference group is what matters.  So unless your reference group 
is centered at z=0 (which it shouldn't, based on the way Gromacs builds boxes), 
you won't have a range equivalent to the z-coordinate.  The output of grompp 
states what the restraint distance is in all cases, based on the coordinates and 
dimensions chosen for the restraint.


-Justin


Thanks again for your help!




MDP FILE USED:
title = Umbrella pulling simulation
define = -DPOSRES_B
define =
; Run parameters
integrator = md
dt = 0.002
tinit = 0
nsteps = 50 ; 1 ns
nstcomm = 10
; Output parameters
nstxout = 5000 ; every 10 ps
nstvout = 5000
nstfout = 5000
nstxtcout = 5000 ; every 10 ps
nstenergy = 5000
; Bond parameters
constraint_algorithm = lincs
constraints = all-bonds
continuation = yes
; Single-range cutoff scheme
nstlist = 5
ns_type = grid
rlist = 1.4
rcoulomb = 1.4
rvdw = 1.4
; PME electrostatics parameters
coulombtype = PME
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
;Berendsen temperature coupling is on
Tcoupl = V-rescale
tau_t = 0.1 0.1 0.1
tc-grps = MOL dop WAT_Cl_Na
ref_t = 300 300 300
; Pressure coupling
Pcoupl = Parrinello-Rahman
Pcoupltype = Semiisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1.0 1.0
compressibility = 4.6E-5 4.6E-5
ref-p = 1.0 1.0
; Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Long-range dispersion correction
DispCorr = EnerPres
; Pull code
pull = umbrella
pull_geometry = distance
pull_dim = N N Y
pull_start = yes
pull_ngroups = 1
pull_group0 = MOL
pull_group1 = Na
pull_init1 = 0
pull_rate1 = 0.0
pull_k1 = 1000 ; kJ mol^-1 nm^-2
pull_nstxout = 1000 ; every 2 ps
pull_nstfout = 1000 ; every 2 ps









 > Date: Thu, 31 May 2012 13:25:26 -0400
 > From: jalem...@vt.edu
 > To: gmx-users@gromacs.org
 > Subject: Re: [gmx-users] boundaries in PMF
 >
 >
 >
 > On 5/31/12 1:20 PM, Rebeca García Fandiño wrote:
 > > Hi,
 > > I am trying to calculate a PMF for an ion along a channel. Everything went 
OK,
 > > but when I used g_wham I got a profile with strange dimensions for the 
x-axis.
 > > What are the boundaries g_wham is using for calculating the units of 
x-axis?
 > >
 >
 > The values are based on the restraint distances along the reaction 
coordinate.
 >
 > > I have used:
 > > g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist
 > > file_histo_output.xvg -unit kCal -cycl yes
 > >
 > > (version 4.0.7)
 > >
 > > and the units I got in the x-axis are determined by the boundaries:
 > >
 > > "Determined boundaries to 0.35 to 2.603290 "
 > >
 > > Which are these units? nm?
 > >
 >
 > Yes.
 >
 > > The z coordinate for my ion explores at least 5 nm!!
 > >
 > > I am a bit confused about that.
 > >
 >
 > The exact output depends on how you set up the umbrella sampling (in the .mdp
 > file). The range of values corresponds to whatever the distances are that are
 > sampled in the various windows. Perhaps there is a sign issue here? Do you
 > have some restraints that are at negative displacement and others at 
positive?
 > Did you set up the .mdp files properly to account for this behavior?
 >
 > -Justin
 >
 > --
 > 
 >
 > Justin A. Lemkul, Ph.D.
 > Research Scientist
 > Department of Biochemistry
 > Virginia Tech
 > Blacksburg, VA
 > jalemkul[at]vt.edu | (540) 231-9080
 > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
 >
 > 
 > --
 > gmx-users mailing list gmx-users@gromacs.org
 > http://lists.gromacs.org/mailman/listinfo/gmx-users
 > Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
 > Please don't post (un)subscribe requests to the list. Use the
 > www interface or send it to gmx-users-requ...@gromacs.org.
 > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists


--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://ww

Re: [gmx-users] mutational analyses: Cystine and indels

[Frederico Moraes Ferreira]

Hi List,
Has anybody have something else to add to the Peter's comments?
Basically, I'm still unsure about modeling deletions.
Any help is appreciated.

Em 30-05-2012 17:37, Peter C. Lai escreveu:

You might be able to use MODELLER for generating the helix deletions since
it is alignment-based. If you use automodel.make(exit_stage=2) it will
generate coordinates based on the specified sequence-spatial alignment but
will stop there (otherwise it will try to run simulated annealing in vacuum).
So what it will do is overlay the sideschains in your new sequence over the
backbone of your beginning structure.

You can then parameterize and refine with gromacs. Distance/angle restraints
will be your friend here (set up distance restraints between N-H to keep the
helix intact while you run energy minimzation and short relaxation MD).

You are correct in that point indels that do not remove or add a full turn
to the helix will be pretty drastic, so you will obviously have to find some
ways to do validation, particularly when assessing if vicinal CYS will form
disulfide bridge or not (which will also obviously do things to helix
geometry) but I think these are hypotheses still well suited for MD to probe.

On 2012-05-30 03:10:39PM -0300, Frederico Moras Ferreira wrote:

Hi Gromacs list,
I would like to study two different mutations in a protein molecule,
perhaps using MD.
One of them is a insert deletion and the other is a mutation to Cys
beside an existing Cys.
The problem is that most of the molecular modeling programs, including
Rosetta, do not handle with indels nor with cystine residues.
Such deletion is located in a helix so as the molecular modeling program
should be capable of moving a entire helix and rotate one or both parts
it. That's not something trivial!
The other mutation sounds worse. Rosetta can't even deal with mutation
to Cys. Assuming that such mutation is not going to disrupt the protein
folding, it probably will form a cystine residue. Not only because of
the distance between cysteines, but also because the pH of the
crystallization conditions.
My question is: has anybody managed such kinds of mutations? or perhaps,
could someone shade some light on these problems?
All the best,
Fred

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

RE: [gmx-users] boundaries in PMF

[Rebeca García Fandiño]

Thanks again, and sorry for insisting, but I am not sure of understanding it 
totally.
So, the boundaries g_wham calculates are not related to the dimensions of my 
channel? Would be any way to convert these units into the position of the ion 
in the channel in each case?
Thanks again a lot for your help.
Best wishes,
Rebeca.

> Date: Thu, 31 May 2012 13:45:18 -0400
> From: jalem...@vt.edu
> To: gmx-users@gromacs.org
> Subject: Re: [gmx-users] boundaries in PMF
> 
> 
> 
> On 5/31/12 1:38 PM, Rebeca García Fandiño wrote:
> > Thanks a lot for your quick answer.
> > The mdp file I have used is copied below. What is strange is that when I 
> > look at
> > the *gro files for the different windows (100 windows in total), i. e:
> >
> > window 1: 8704Na Na56458 5.134 5.085 5.824
> > window 50: 8704Na Na56458 5.053 5.081 3.459
> > window 100: 8704Na Na56458 4.990 5.042 0.951
> >
> > you can see that the z-coordinate goes from 0.951 to 5.824 nm
> >
> > I should have a total distance in the x-axis of about 5 nm, and however, the
> > boundaries calculated by g_wham are
> >
> > "Determined boundaries to 0.35 to 2.603290 "
> >
> > Can you see anything in the mdp file which is causing this behaviour...?
> >
> 
> No, but you need to be careful how you're interpreting what you get.  The 
> z-coordinate of your ion is not what's relevant; its displacement along the 
> z-axis from the reference group is what matters.  So unless your reference 
> group 
> is centered at z=0 (which it shouldn't, based on the way Gromacs builds 
> boxes), 
> you won't have a range equivalent to the z-coordinate.  The output of grompp 
> states what the restraint distance is in all cases, based on the coordinates 
> and 
> dimensions chosen for the restraint.
> 
> -Justin
> 
> > Thanks again for your help!
> >
> >
> >
> >
> > MDP FILE USED:
> > title = Umbrella pulling simulation
> > define = -DPOSRES_B
> > define =
> > ; Run parameters
> > integrator = md
> > dt = 0.002
> > tinit = 0
> > nsteps = 50 ; 1 ns
> > nstcomm = 10
> > ; Output parameters
> > nstxout = 5000 ; every 10 ps
> > nstvout = 5000
> > nstfout = 5000
> > nstxtcout = 5000 ; every 10 ps
> > nstenergy = 5000
> > ; Bond parameters
> > constraint_algorithm = lincs
> > constraints = all-bonds
> > continuation = yes
> > ; Single-range cutoff scheme
> > nstlist = 5
> > ns_type = grid
> > rlist = 1.4
> > rcoulomb = 1.4
> > rvdw = 1.4
> > ; PME electrostatics parameters
> > coulombtype = PME
> > fourierspacing = 0.12
> > fourier_nx = 0
> > fourier_ny = 0
> > fourier_nz = 0
> > pme_order = 4
> > ewald_rtol = 1e-5
> > optimize_fft = yes
> > ;Berendsen temperature coupling is on
> > Tcoupl = V-rescale
> > tau_t = 0.1 0.1 0.1
> > tc-grps = MOL dop WAT_Cl_Na
> > ref_t = 300 300 300
> > ; Pressure coupling
> > Pcoupl = Parrinello-Rahman
> > Pcoupltype = Semiisotropic
> > ; Time constant (ps), compressibility (1/bar) and reference P (bar)
> > tau-p = 1.0 1.0
> > compressibility = 4.6E-5 4.6E-5
> > ref-p = 1.0 1.0
> > ; Generate velocities is off
> > gen_vel = no
> > ; Periodic boundary conditions are on in all directions
> > pbc = xyz
> > ; Long-range dispersion correction
> > DispCorr = EnerPres
> > ; Pull code
> > pull = umbrella
> > pull_geometry = distance
> > pull_dim = N N Y
> > pull_start = yes
> > pull_ngroups = 1
> > pull_group0 = MOL
> > pull_group1 = Na
> > pull_init1 = 0
> > pull_rate1 = 0.0
> > pull_k1 = 1000 ; kJ mol^-1 nm^-2
> > pull_nstxout = 1000 ; every 2 ps
> > pull_nstfout = 1000 ; every 2 ps
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >  > Date: Thu, 31 May 2012 13:25:26 -0400
> >  > From: jalem...@vt.edu
> >  > To: gmx-users@gromacs.org
> >  > Subject: Re: [gmx-users] boundaries in PMF
> >  >
> >  >
> >  >
> >  > On 5/31/12 1:20 PM, Rebeca García Fandiño wrote:
> >  > > Hi,
> >  > > I am trying to calculate a PMF for an ion along a channel. Everything 
> > went OK,
> >  > > but when I used g_wham I got a profile with strange dimensions for the 
> > x-axis.
> >  > > What are the boundaries g_wham is using for calculating the units of 
> > x-axis?
> >  > >
> >  >
> >  > The values are based on the restraint distances along the reaction 
> > coordinate.
> >  >
> >  > > I have used:
> >  > > g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist
> >  > > file_histo_output.xvg -unit kCal -cycl yes
> >  > >
> >  > > (version 4.0.7)
> >  > >
> >  > > and the units I got in the x-axis are determined by the boundaries:
> >  > >
> >  > > "Determined boundaries to 0.35 to 2.603290 "
> >  > >
> >  > > Which are these units? nm?
> >  > >
> >  >
> >  > Yes.
> >  >
> >  > > The z coordinate for my ion explores at least 5 nm!!
> >  > >
> >  > > I am a bit confused about that.
> >  > >
> >  >
> >  > The exact output depends on how you set up the umbrella sampling (in the 
> > .mdp
> >  > file). The range of values corresponds to whatever the distances are 
> > that are
> >  > sampled in the various windows. Perhaps there is a sign issue here? Do 
> > you
> >  > have so

Re: [gmx-users] boundaries in PMF

[Justin A. Lemkul]

On 5/31/12 1:58 PM, Rebeca García Fandiño wrote:

Thanks again, and sorry for insisting, but I am not sure of understanding it
totally.
So, the boundaries g_wham calculates are not related to the dimensions of my
channel? Would be any way to convert these units into the position of the ion in
the channel in each case?


The pullx.xvg files contain helpful information in this case - the COM position 
of the reference group and restrained group along the restrained dimension(s). 
This should be what you're looking for, correct?


The WHAM algorithm doesn't care about the absolute position of anything, it's 
the COM distance between the restrained and reference groups that matter.  It's 
that figure that should motivate how the windows were constructed and then analyzed.


-Justin


Thanks again a lot for your help.
Best wishes,
Rebeca.

 > Date: Thu, 31 May 2012 13:45:18 -0400
 > From: jalem...@vt.edu
 > To: gmx-users@gromacs.org
 > Subject: Re: [gmx-users] boundaries in PMF
 >
 >
 >
 > On 5/31/12 1:38 PM, Rebeca García Fandiño wrote:
 > > Thanks a lot for your quick answer.
 > > The mdp file I have used is copied below. What is strange is that when I
look at
 > > the *gro files for the different windows (100 windows in total), i. e:
 > >
 > > window 1: 8704Na Na56458 5.134 5.085 5.824
 > > window 50: 8704Na Na56458 5.053 5.081 3.459
 > > window 100: 8704Na Na56458 4.990 5.042 0.951
 > >
 > > you can see that the z-coordinate goes from 0.951 to 5.824 nm
 > >
 > > I should have a total distance in the x-axis of about 5 nm, and however, 
the
 > > boundaries calculated by g_wham are
 > >
 > > "Determined boundaries to 0.35 to 2.603290 "
 > >
 > > Can you see anything in the mdp file which is causing this behaviour...?
 > >
 >
 > No, but you need to be careful how you're interpreting what you get. The
 > z-coordinate of your ion is not what's relevant; its displacement along the
 > z-axis from the reference group is what matters. So unless your reference 
group
 > is centered at z=0 (which it shouldn't, based on the way Gromacs builds 
boxes),
 > you won't have a range equivalent to the z-coordinate. The output of grompp
 > states what the restraint distance is in all cases, based on the coordinates 
and
 > dimensions chosen for the restraint.
 >
 > -Justin
 >
 > > Thanks again for your help!
 > >
 > >
 > >
 > >
 > > MDP FILE USED:
 > > title = Umbrella pulling simulation
 > > define = -DPOSRES_B
 > > define =
 > > ; Run parameters
 > > integrator = md
 > > dt = 0.002
 > > tinit = 0
 > > nsteps = 50 ; 1 ns
 > > nstcomm = 10
 > > ; Output parameters
 > > nstxout = 5000 ; every 10 ps
 > > nstvout = 5000
 > > nstfout = 5000
 > > nstxtcout = 5000 ; every 10 ps
 > > nstenergy = 5000
 > > ; Bond parameters
 > > constraint_algorithm = lincs
 > > constraints = all-bonds
 > > continuation = yes
 > > ; Single-range cutoff scheme
 > > nstlist = 5
 > > ns_type = grid
 > > rlist = 1.4
 > > rcoulomb = 1.4
 > > rvdw = 1.4
 > > ; PME electrostatics parameters
 > > coulombtype = PME
 > > fourierspacing = 0.12
 > > fourier_nx = 0
 > > fourier_ny = 0
 > > fourier_nz = 0
 > > pme_order = 4
 > > ewald_rtol = 1e-5
 > > optimize_fft = yes
 > > ;Berendsen temperature coupling is on
 > > Tcoupl = V-rescale
 > > tau_t = 0.1 0.1 0.1
 > > tc-grps = MOL dop WAT_Cl_Na
 > > ref_t = 300 300 300
 > > ; Pressure coupling
 > > Pcoupl = Parrinello-Rahman
 > > Pcoupltype = Semiisotropic
 > > ; Time constant (ps), compressibility (1/bar) and reference P (bar)
 > > tau-p = 1.0 1.0
 > > compressibility = 4.6E-5 4.6E-5
 > > ref-p = 1.0 1.0
 > > ; Generate velocities is off
 > > gen_vel = no
 > > ; Periodic boundary conditions are on in all directions
 > > pbc = xyz
 > > ; Long-range dispersion correction
 > > DispCorr = EnerPres
 > > ; Pull code
 > > pull = umbrella
 > > pull_geometry = distance
 > > pull_dim = N N Y
 > > pull_start = yes
 > > pull_ngroups = 1
 > > pull_group0 = MOL
 > > pull_group1 = Na
 > > pull_init1 = 0
 > > pull_rate1 = 0.0
 > > pull_k1 = 1000 ; kJ mol^-1 nm^-2
 > > pull_nstxout = 1000 ; every 2 ps
 > > pull_nstfout = 1000 ; every 2 ps
 > >
 > >
 > >
 > >
 > >
 > >
 > >
 > >
 > >
 > > > Date: Thu, 31 May 2012 13:25:26 -0400
 > > > From: jalem...@vt.edu
 > > > To: gmx-users@gromacs.org
 > > > Subject: Re: [gmx-users] boundaries in PMF
 > > >
 > > >
 > > >
 > > > On 5/31/12 1:20 PM, Rebeca García Fandiño wrote:
 > > > > Hi,
 > > > > I am trying to calculate a PMF for an ion along a channel. Everything
went OK,
 > > > > but when I used g_wham I got a profile with strange dimensions for the
x-axis.
 > > > > What are the boundaries g_wham is using for calculating the units of
x-axis?
 > > > >
 > > >
 > > > The values are based on the restraint distances along the reaction
coordinate.
 > > >
 > > > > I have used:
 > > > > g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist
 > > > > file_histo_output.xvg -unit kCal -cycl yes
 > > > >
 > > > > (version 4.0.7)
 > > > >
 > > > > and the units I got in the x-axis

[gmx-users] Re: boundaries in PMF

[Thomas Schlesier]

Where is the center of mass of reference group (MOL) located?

It seems that the COM is near the middle of the ion channel. Since you 
use 'pull_geometry=distance', g_wham will look only for the distance 
between 'MOL' and 'Na' and that leads to problem.
If the com of 'MOL' sits in the center of the channel (which is around 
5nm long), one has 2.5nm in both directions. Since g_wham uses only the 
distance you get the PMF for half of the channel, but with the data of 
both parts.
If the channel would be symmetric and the com of 'MOL' would lie exactly 
in the middle of the channel, this could be ok. But if even one of both 
assumptions is wrong, the results would be errorous.


A better approach would be to use 'pull_geometry=direction', since the 
you define a vector along which the windows lie and do not have the 
problem that the distance is in both directions (along positive and 
negative vector) the same.
Only problem could be that g_wham supports 'pull_geometry=direction' 
only from gromacs 4.5.x (don't know this, since instead of umbrella 
smapling i use thermodynamik integration, where one uses constraints 
(instead of restraints) and integrates the constraint-force; but the 
conceptual problem with 'distance/direction' is the same).


Another approach (with 'pull_geometry=distance') would be to use a 
reference group which is just outside of the channel (like going some 
steps away from the channel, along the vector which goes through the 
channel). If there is only water, it would be bad, because then the 
reference group would be move away.
But then on could use the entry and exit of the channel as a reference 
point for two simulations. In the case that the entry is the reference 
group, the PMF would be ill defined near the entry, but from the 
simulation with the exit as reference you would get the right PMF for 
the entry region and vice versa.


Greetings
Thomas


Am 31.05.2012 19:39, schrieb gmx-users-requ...@gromacs.org:

Thanks a lot for your quick answer. The mdp file I have used is copied
below. What is strange is that when I look at the *gro files for the
different windows (100 windows in total), i. e: window 1: 8704Na Na56458
5.134 5.085 5.824 window 50: 8704Na Na56458 5.053 5.081 3.459 window
100: 8704Na Na56458 4.990 5.042 0.951 you can see that the z-coordinate
goes from 0.951 to 5.824 nm I should have a total distance in the x-axis
of about 5 nm, and however, the boundaries calculated by g_wham are
"Determined boundaries to 0.35 to 2.603290 " Can you see anything in
the mdp file which is causing this behaviour...? Thanks again for your
help! MDP FILE USED: title = Umbrella pulling simulation define =
-DPOSRES_B define = ; Run parameters integrator = md dt = 0.002 tinit =
0 nsteps = 50 ; 1 ns nstcomm = 10 ; Output parameters nstxout = 5000
; every 10 ps nstvout = 5000 nstfout = 5000 nstxtcout = 5000 ; every 10
ps nstenergy = 5000 ; Bond parameters constraint_algorithm = lincs
constraints = all-bonds continuation = yes ; Single-range cutoff scheme
nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb = 1.4 rvdw = 1.4 ; PME
electrostatics parameters coulombtype = PME fourierspacing = 0.12
fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol =
1e-5 optimize_fft = yes ;Berendsen temperature coupling is on Tcoupl =
V-rescale tau_t = 0.1 0.1 0.1 tc-grps = MOL dop WAT_Cl_Na ref_t = 300
300 300 ; Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype =
Semiisotropic ; Time constant (ps), compressibility (1/bar) and
reference P (bar) tau-p = 1.0 1.0 compressibility = 4.6E-5 4.6E-5 ref-p
= 1.0 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary
conditions are on in all directions pbc = xyz ; Long-range dispersion
correction DispCorr = EnerPres ; Pull code pull = umbrella pull_geometry
= distance pull_dim = N N Y pull_start = yes pull_ngroups = 1
pull_group0 = MOL pull_group1 = Na pull_init1 = 0 pull_rate1 = 0.0
pull_k1 = 1000 ; kJ mol^-1 nm^-2 pull_nstxout = 1000 ; every 2 ps
pull_nstfout = 1000 ; every 2 ps

>  Date: Thu, 31 May 2012 13:25:26 -0400
>  From:jalem...@vt.edu
>  To:gmx-users@gromacs.org
>  Subject: Re: [gmx-users] boundaries in PMF
>
>
>
>  On 5/31/12 1:20 PM, Rebeca Garc?a Fandi?o wrote:

>  >  Hi,
>  >  I am trying to calculate a PMF for an ion along a channel. Everything 
went OK,
>  >  but when I used g_wham I got a profile with strange dimensions for the 
x-axis.
>  >  What are the boundaries g_wham is using for calculating the units of 
x-axis?
>  >

>
>  The values are based on the restraint distances along the reaction 
coordinate.
>

>  >  I have used:
>  >  g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist
>  >  file_histo_output.xvg -unit kCal -cycl yes
>  >
>  >  (version 4.0.7)
>  >
>  >  and the units I got in the x-axis are determined by the boundaries:
>  >
>  >  "Determined boundaries to 0.35 to 2.603290"
>  >
>  >  Which are these units? nm?
>  >

>
>  Yes.
>

>  >  The z coordinate for my ion explores

RE: [gmx-users] Re: boundaries in PMF

[Rebeca García Fandiño]

Hi,
the center of mass of my channel is at the middle of the ion channel, and it is 
a symmetric system, so I suppose these results would be OK. Anyway, I will 
check the options you propose.
Thanks a lot for all your comments!!
Best wishes,
Rebeca.

> Date: Thu, 31 May 2012 20:08:26 +0200
> From: schl...@uni-mainz.de
> To: gmx-users@gromacs.org
> Subject: [gmx-users] Re: boundaries in PMF
> 
> Where is the center of mass of reference group (MOL) located?
> 
> It seems that the COM is near the middle of the ion channel. Since you 
> use 'pull_geometry=distance', g_wham will look only for the distance 
> between 'MOL' and 'Na' and that leads to problem.
> If the com of 'MOL' sits in the center of the channel (which is around 
> 5nm long), one has 2.5nm in both directions. Since g_wham uses only the 
> distance you get the PMF for half of the channel, but with the data of 
> both parts.
> If the channel would be symmetric and the com of 'MOL' would lie exactly 
> in the middle of the channel, this could be ok. But if even one of both 
> assumptions is wrong, the results would be errorous.
> 
> A better approach would be to use 'pull_geometry=direction', since the 
> you define a vector along which the windows lie and do not have the 
> problem that the distance is in both directions (along positive and 
> negative vector) the same.
> Only problem could be that g_wham supports 'pull_geometry=direction' 
> only from gromacs 4.5.x (don't know this, since instead of umbrella 
> smapling i use thermodynamik integration, where one uses constraints 
> (instead of restraints) and integrates the constraint-force; but the 
> conceptual problem with 'distance/direction' is the same).
> 
> Another approach (with 'pull_geometry=distance') would be to use a 
> reference group which is just outside of the channel (like going some 
> steps away from the channel, along the vector which goes through the 
> channel). If there is only water, it would be bad, because then the 
> reference group would be move away.
> But then on could use the entry and exit of the channel as a reference 
> point for two simulations. In the case that the entry is the reference 
> group, the PMF would be ill defined near the entry, but from the 
> simulation with the exit as reference you would get the right PMF for 
> the entry region and vice versa.
> 
> Greetings
> Thomas
> 
> 
> Am 31.05.2012 19:39, schrieb gmx-users-requ...@gromacs.org:
> > Thanks a lot for your quick answer. The mdp file I have used is copied
> > below. What is strange is that when I look at the *gro files for the
> > different windows (100 windows in total), i. e: window 1: 8704Na Na56458
> > 5.134 5.085 5.824 window 50: 8704Na Na56458 5.053 5.081 3.459 window
> > 100: 8704Na Na56458 4.990 5.042 0.951 you can see that the z-coordinate
> > goes from 0.951 to 5.824 nm I should have a total distance in the x-axis
> > of about 5 nm, and however, the boundaries calculated by g_wham are
> > "Determined boundaries to 0.35 to 2.603290 " Can you see anything in
> > the mdp file which is causing this behaviour...? Thanks again for your
> > help! MDP FILE USED: title = Umbrella pulling simulation define =
> > -DPOSRES_B define = ; Run parameters integrator = md dt = 0.002 tinit =
> > 0 nsteps = 50 ; 1 ns nstcomm = 10 ; Output parameters nstxout = 5000
> > ; every 10 ps nstvout = 5000 nstfout = 5000 nstxtcout = 5000 ; every 10
> > ps nstenergy = 5000 ; Bond parameters constraint_algorithm = lincs
> > constraints = all-bonds continuation = yes ; Single-range cutoff scheme
> > nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb = 1.4 rvdw = 1.4 ; PME
> > electrostatics parameters coulombtype = PME fourierspacing = 0.12
> > fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol =
> > 1e-5 optimize_fft = yes ;Berendsen temperature coupling is on Tcoupl =
> > V-rescale tau_t = 0.1 0.1 0.1 tc-grps = MOL dop WAT_Cl_Na ref_t = 300
> > 300 300 ; Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype =
> > Semiisotropic ; Time constant (ps), compressibility (1/bar) and
> > reference P (bar) tau-p = 1.0 1.0 compressibility = 4.6E-5 4.6E-5 ref-p
> > = 1.0 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary
> > conditions are on in all directions pbc = xyz ; Long-range dispersion
> > correction DispCorr = EnerPres ; Pull code pull = umbrella pull_geometry
> > = distance pull_dim = N N Y pull_start = yes pull_ngroups = 1
> > pull_group0 = MOL pull_group1 = Na pull_init1 = 0 pull_rate1 = 0.0
> > pull_k1 = 1000 ; kJ mol^-1 nm^-2 pull_nstxout = 1000 ; every 2 ps
> > pull_nstfout = 1000 ; every 2 ps
> >> >  Date: Thu, 31 May 2012 13:25:26 -0400
> >> >  From:jalem...@vt.edu
> >> >  To:gmx-users@gromacs.org
> >> >  Subject: Re: [gmx-users] boundaries in PMF
> >> >
> >> >
> >> >
> >> >  On 5/31/12 1:20 PM, Rebeca Garc?a Fandi?o wrote:
> >>> >  >  Hi,
> >>> >  >  I am trying to calculate a PMF for an ion along a channel. 
> >>> > Everything went OK,
> >>> >  >  but when I u

Re: [gmx-users] Re: boundaries in PMF

[Justin A. Lemkul]

On 5/31/12 3:37 PM, Rebeca García Fandiño wrote:

Hi,
the center of mass of my channel is at the middle of the ion channel, and it is
a symmetric system, so I suppose these results would be OK. Anyway, I will check
the options you propose.


If you are sampling regions above and below the channel/membrane, then the 
"distance" geometry is not appropriate; you need either "direction" or (perhaps 
the most flexible option) "position."  There are a number of useful discussions 
on such topics in the list archive and in my umbrella sampling tutorial for you 
to consider.  You can likely create a series of .tpr files from new .mdp files 
with correct options to run the analysis.


I suspect this is why g_wham is finding a range of values only equal to half of 
your expected reaction coordinate; it is considering only the positive 
displacement along the reaction coordinate.


-Justin


Thanks a lot for all your comments!!
Best wishes,
Rebeca.

 > Date: Thu, 31 May 2012 20:08:26 +0200
 > From: schl...@uni-mainz.de
 > To: gmx-users@gromacs.org
 > Subject: [gmx-users] Re: boundaries in PMF
 >
 > Where is the center of mass of reference group (MOL) located?
 >
 > It seems that the COM is near the middle of the ion channel. Since you
 > use 'pull_geometry=distance', g_wham will look only for the distance
 > between 'MOL' and 'Na' and that leads to problem.
 > If the com of 'MOL' sits in the center of the channel (which is around
 > 5nm long), one has 2.5nm in both directions. Since g_wham uses only the
 > distance you get the PMF for half of the channel, but with the data of
 > both parts.
 > If the channel would be symmetric and the com of 'MOL' would lie exactly
 > in the middle of the channel, this could be ok. But if even one of both
 > assumptions is wrong, the results would be errorous.
 >
 > A better approach would be to use 'pull_geometry=direction', since the
 > you define a vector along which the windows lie and do not have the
 > problem that the distance is in both directions (along positive and
 > negative vector) the same.
 > Only problem could be that g_wham supports 'pull_geometry=direction'
 > only from gromacs 4.5.x (don't know this, since instead of umbrella
 > smapling i use thermodynamik integration, where one uses constraints
 > (instead of restraints) and integrates the constraint-force; but the
 > conceptual problem with 'distance/direction' is the same).
 >
 > Another approach (with 'pull_geometry=distance') would be to use a
 > reference group which is just outside of the channel (like going some
 > steps away from the channel, along the vector which goes through the
 > channel). If there is only water, it would be bad, because then the
 > reference group would be move away.
 > But then on could use the entry and exit of the channel as a reference
 > point for two simulations. In the case that the entry is the reference
 > group, the PMF would be ill defined near the entry, but from the
 > simulation with the exit as reference you would get the right PMF for
 > the entry region and vice versa.
 >
 > Greetings
 > Thomas
 >
 >
 > Am 31.05.2012 19:39, schrieb gmx-users-requ...@gromacs.org:
 > > Thanks a lot for your quick answer. The mdp file I have used is copied
 > > below. What is strange is that when I look at the *gro files for the
 > > different windows (100 windows in total), i. e: window 1: 8704Na Na56458
 > > 5.134 5.085 5.824 window 50: 8704Na Na56458 5.053 5.081 3.459 window
 > > 100: 8704Na Na56458 4.990 5.042 0.951 you can see that the z-coordinate
 > > goes from 0.951 to 5.824 nm I should have a total distance in the x-axis
 > > of about 5 nm, and however, the boundaries calculated by g_wham are
 > > "Determined boundaries to 0.35 to 2.603290 " Can you see anything in
 > > the mdp file which is causing this behaviour...? Thanks again for your
 > > help! MDP FILE USED: title = Umbrella pulling simulation define =
 > > -DPOSRES_B define = ; Run parameters integrator = md dt = 0.002 tinit =
 > > 0 nsteps = 50 ; 1 ns nstcomm = 10 ; Output parameters nstxout = 5000
 > > ; every 10 ps nstvout = 5000 nstfout = 5000 nstxtcout = 5000 ; every 10
 > > ps nstenergy = 5000 ; Bond parameters constraint_algorithm = lincs
 > > constraints = all-bonds continuation = yes ; Single-range cutoff scheme
 > > nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb = 1.4 rvdw = 1.4 ; PME
 > > electrostatics parameters coulombtype = PME fourierspacing = 0.12
 > > fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol =
 > > 1e-5 optimize_fft = yes ;Berendsen temperature coupling is on Tcoupl =
 > > V-rescale tau_t = 0.1 0.1 0.1 tc-grps = MOL dop WAT_Cl_Na ref_t = 300
 > > 300 300 ; Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype =
 > > Semiisotropic ; Time constant (ps), compressibility (1/bar) and
 > > reference P (bar) tau-p = 1.0 1.0 compressibility = 4.6E-5 4.6E-5 ref-p
 > > = 1.0 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary
 > > conditions are on in all

Re: [gmx-users] Re: boundaries in PMF

[Justin A. Lemkul]

On 5/31/12 3:51 PM, Rebeca García Fandiño wrote:

OK,
however, just one point about your last comment:

 > I suspect this is why g_wham is finding a range of values only equal to half 
of
 > your expected reaction coordinate; it is considering only the positive
 > displacement along the reaction coordinate.

It seems like all the channel is explored, not only one half. If g_wham was only
considering the positive displacement I should see a profile for just one half
of the channel, shouldn´t I? and I can see a profile typical for the entire 
channel.



Then perhaps it's just a small g_wham output bug (boundary values).  You didn't 
mention that your profile looked normal ;)


-Justin

--


Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

RE: [gmx-users] Re: boundaries in PMF

[Rebeca García Fandiño]

Yes,
my profile seems normal, the only problem are the units in the x-axis, because 
I expected them to be in the range of the dimensions of the channel.
I will try if I see differences with the other options you proposed, anyway.
Thanks a lot for your help!!
Best wishes, 
Rebeca.

> Date: Thu, 31 May 2012 15:52:48 -0400
> From: jalem...@vt.edu
> To: gmx-users@gromacs.org
> Subject: Re: [gmx-users] Re: boundaries in PMF
> 
> 
> 
> On 5/31/12 3:51 PM, Rebeca García Fandiño wrote:
> > OK,
> > however, just one point about your last comment:
> >
> >  > I suspect this is why g_wham is finding a range of values only equal to 
> > half of
> >  > your expected reaction coordinate; it is considering only the positive
> >  > displacement along the reaction coordinate.
> >
> > It seems like all the channel is explored, not only one half. If g_wham was 
> > only
> > considering the positive displacement I should see a profile for just one 
> > half
> > of the channel, shouldn´t I? and I can see a profile typical for the entire 
> > channel.
> >
> 
> Then perhaps it's just a small g_wham output bug (boundary values).  You 
> didn't 
> mention that your profile looked normal ;)
> 
> -Justin
> 
> -- 
> 
> 
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> -- 
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
  -- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

RE: [gmx-users] Re: boundaries in PMF

[Rebeca García Fandiño]

OK,
however, just one point about your last comment:

> I suspect this is why g_wham is finding a range of values only equal to half 
> of 
> your expected reaction coordinate; it is considering only the positive 
> displacement along the reaction coordinate.

It seems like all the channel is explored, not only one half. If g_wham was 
only considering the positive displacement I should see a profile for just one 
half of the channel, shouldn´t I? and I can see a profile typical for the 
entire channel.

Rebeca.

> Date: Thu, 31 May 2012 15:42:06 -0400
> From: jalem...@vt.edu
> To: gmx-users@gromacs.org
> Subject: Re: [gmx-users] Re: boundaries in PMF
> 
> 
> 
> On 5/31/12 3:37 PM, Rebeca García Fandiño wrote:
> > Hi,
> > the center of mass of my channel is at the middle of the ion channel, and 
> > it is
> > a symmetric system, so I suppose these results would be OK. Anyway, I will 
> > check
> > the options you propose.
> 
> If you are sampling regions above and below the channel/membrane, then the 
> "distance" geometry is not appropriate; you need either "direction" or 
> (perhaps 
> the most flexible option) "position."  There are a number of useful 
> discussions 
> on such topics in the list archive and in my umbrella sampling tutorial for 
> you 
> to consider.  You can likely create a series of .tpr files from new .mdp 
> files 
> with correct options to run the analysis.
> 
> I suspect this is why g_wham is finding a range of values only equal to half 
> of 
> your expected reaction coordinate; it is considering only the positive 
> displacement along the reaction coordinate.
> 
> -Justin
> 
> > Thanks a lot for all your comments!!
> > Best wishes,
> > Rebeca.
> >
> >  > Date: Thu, 31 May 2012 20:08:26 +0200
> >  > From: schl...@uni-mainz.de
> >  > To: gmx-users@gromacs.org
> >  > Subject: [gmx-users] Re: boundaries in PMF
> >  >
> >  > Where is the center of mass of reference group (MOL) located?
> >  >
> >  > It seems that the COM is near the middle of the ion channel. Since you
> >  > use 'pull_geometry=distance', g_wham will look only for the distance
> >  > between 'MOL' and 'Na' and that leads to problem.
> >  > If the com of 'MOL' sits in the center of the channel (which is around
> >  > 5nm long), one has 2.5nm in both directions. Since g_wham uses only the
> >  > distance you get the PMF for half of the channel, but with the data of
> >  > both parts.
> >  > If the channel would be symmetric and the com of 'MOL' would lie exactly
> >  > in the middle of the channel, this could be ok. But if even one of both
> >  > assumptions is wrong, the results would be errorous.
> >  >
> >  > A better approach would be to use 'pull_geometry=direction', since the
> >  > you define a vector along which the windows lie and do not have the
> >  > problem that the distance is in both directions (along positive and
> >  > negative vector) the same.
> >  > Only problem could be that g_wham supports 'pull_geometry=direction'
> >  > only from gromacs 4.5.x (don't know this, since instead of umbrella
> >  > smapling i use thermodynamik integration, where one uses constraints
> >  > (instead of restraints) and integrates the constraint-force; but the
> >  > conceptual problem with 'distance/direction' is the same).
> >  >
> >  > Another approach (with 'pull_geometry=distance') would be to use a
> >  > reference group which is just outside of the channel (like going some
> >  > steps away from the channel, along the vector which goes through the
> >  > channel). If there is only water, it would be bad, because then the
> >  > reference group would be move away.
> >  > But then on could use the entry and exit of the channel as a reference
> >  > point for two simulations. In the case that the entry is the reference
> >  > group, the PMF would be ill defined near the entry, but from the
> >  > simulation with the exit as reference you would get the right PMF for
> >  > the entry region and vice versa.
> >  >
> >  > Greetings
> >  > Thomas
> >  >
> >  >
> >  > Am 31.05.2012 19:39, schrieb gmx-users-requ...@gromacs.org:
> >  > > Thanks a lot for your quick answer. The mdp file I have used is copied
> >  > > below. What is strange is that when I look at the *gro files for the
> >  > > different windows (100 windows in total), i. e: window 1: 8704Na 
> > Na56458
> >  > > 5.134 5.085 5.824 window 50: 8704Na Na56458 5.053 5.081 3.459 window
> >  > > 100: 8704Na Na56458 4.990 5.042 0.951 you can see that the z-coordinate
> >  > > goes from 0.951 to 5.824 nm I should have a total distance in the 
> > x-axis
> >  > > of about 5 nm, and however, the boundaries calculated by g_wham are
> >  > > "Determined boundaries to 0.35 to 2.603290 " Can you see anything 
> > in
> >  > > the mdp file which is causing this behaviour...? Thanks again for your
> >  > > help! MDP FILE USED: title = Umbrella pulling simulation define =
> >  > > -DPOSRES_B define = ; Run parameters integrator = md dt = 0.002 tinit =
> >  > > 0 nsteps = 5