[gmx-users] Re: Problems with umbrella sampling

[Thomas Schlesier]

Hi,
Berk is right, for the actual umbrella sampling 'pull=umbrella' is 
needed and 'pull_dim = Y Y Y' is also better (you want your umbrella 
potential in all directions).
To generate the window structures 'pull_dim= Y Y Y' is also better (with 
'Y N N' you need to be sure that the two ions are parallel to the x-axis).
But the biggest problem is 'pull_start=yes' AND 'pull_init1=0.1':
With 'pull_start=yes' gromacs calculates the distance between the two 
ions and uses this as 'pull_init', then with the second thing you 
increase this this value then. So the spring doesn't lie on top of I- 
but 0.1nm behind it, and so there should be an initial force.
The other thing where i'm not sure about is the 'pull_vec1'. I left it 
blank (default then was 0.0 0.0 0.0) because i thought with pull_dim=Y Y 
Y  and pull_geometry=distance gromacs calculates the vector so, that the 
pulling is along the vector connecting the two ions. Don't know if a 
entry in 'pull_vec1' changes something.
To pull the two ions together i should be sufficient to change the sign 
of the 'pull_rate1'.
Greetings
Thomas





Hi,

I think you want pull_dim = Y Y Y
and for umbrella sampling you need pull=umbrella

Berk

Date: Thu, 11 Mar 2010 17:37:34 -0800
From: jrcasey20 at gmail.com
To: gmx-users at gromacs.org
Subject: [gmx-users] Problems with umbrella sampling

Hello,

I am interested in creating a PMF for Na+ and I- in the presence of THF 
(tetrahydrofuran).  I am having a couple of issues though.

First, I have followed the umbrella sampling tutorial on the gromacs 
website very careful, and I think I am doing everything right.  I 
minimized the energy of my system, did a short equilibrium run, and then 
pulled Na+ and I- away from each other at a specified rate.  This 
generated starting configurations, which I used in my umbrella sampling 
(see attached file ).  Once I had many umbrella sampling windows, I used 
g_wham to recombine them.  First I looked at the histograms.  I have 
good overlap, but I notice some very odd features:


1)  I see substantial shifting towards lower distances between Na+ and 
I-, no matter what the initial distances they start separated at.  I 
would expect this if I am near an equilibrium distance, but I see it for 
all distances.

2)  My original pull simulation (where I specificed a rate) started Na+ 
and I- 0.27 nm apart.  All of my starting configurations are 0.27 nm 
apart or further.  If I run some of the umbrella sampling simulations at 
0.27 nm, I see that the Na+ and I- reach distances as low as 0.03 nm! 
This is completely unphysical (both in that my spring constant is strong 
enough to supposedly prevent such large drifts and because the two 
molecules should not want to be on top of each other).


To check these issues out, I tried to do a couple of things:

I want to pull the Na+ and I- closer together than 0.27 nm, and see if I 
get the steep increase in energy, like I expect.  I am having problems 
with this part though.  I have not been able to pull my molecules closer 
together when starting at 0.27 nm, so I tried to start them closer 
together and then pull them apart by using the following in my .mdp file:


pull = constraint
pull_geometry = distance
pull_dim = Y N N
pull_start = yes
pull_ngroups = 1
pull_group0 = Na+
pull_group1 = I-
pull_vec1 = 1 0 0
pull_init1 = 0.1 (I think this means that I- starts 0.1 nm away from Na+?)

pull_rate1 = 0.005

When I do this, I still get my molecules being pulled apart from 0.27 nm 
to further distances apart.  I also tried using no pull_init1, but using 
pull_vec1 = -1 0 0, but that didn't work either.  Whenever I tried to 
use pull_geometry = position, I got segmentation faults...


Has anyone else experienced these issues?  Is there a way to bring Na+ 
and I- closer?  I have checked older posts, but the suggestions I keep 
seeing (i.e. use pull_init1) aren't working for me.

Thank you!

Jenny