[gmx-users] problem with the energy minimization of a system composed of two NA+ ions and explicit water

[Justin A. Lemkul]

Ozge Engin wrote:
> Hi all,
> 
> First of all, thank you for your quick reply Justin. Here, I want to 
> state the problem in a clearer way.
> 
> I have been trying to calculate the PMF between two NA+ ions in explicit 
> water. In order to do that I first pulled one of the NA+ ions 
> continuously with respect to the other until the farthest point on the 
> PMF curve was reached. After that I split the continuous pulling 
> trajectory from the points that are in correspondence with the 
> constraint points of the PMF curve. These were done via *genera_start_4* 
> script.
> 

You don't have a PMF curve.  I expect you're referring either to the pullf.xvg 
or pullx.xvg files generated from the COM pulling output.  PMF is calculated 
later, through umbrella sampling.

> I took each starting conformation, that was obtained from continuous 
> pulling, solvated with water, and then energy minimized. Then, I 
> prepared the systems for constraint pulling MD production runs. These 

Weren't these systems already solvated in the above pulling?  If you add water 
differently to each system, you may get unreliable results.  Use 
water-containing configurations from the original SMD trajectory (above) as your 
input configurations.

> were done via *input* script. Interestingly, systems (ion+water) at some 
> of the constraint points could not been energy minimized, it reached to 
> the machine epsilon, but not to the desired criteria. The force per atom 
> was too high, around the order of ^4. Moreover, this issue was not 
> distance-dependent, which means that the energy minimization worked for 
> smaller constraint points, but not for some of the larger constraint 
> points. When I looked at the problematic trajectories, I saw that bonds 
> of a water molecule near the box boundary was longer. I thought that 
> this may be due to genbox command. Because of the small radius of NA+ 

Perhaps.  See the above comment.  It's better to not re-solvate after the fact. 
  Was the initial pull actually conducted in vacuo?

> ion it may not check the overlaps among the atoms. In order to test 
> whether it is the case I increased the LJ parameters of the ion. This 
> worked for some of the points, but not for others. (Note that the 
> diameter is still smaller than the constraint distance between the two 
> ions.) 
> 
> I use ffG45a3 force field, and the scripts given above. I tried this 
> with different versions of Gromacs: 4.0.7, 4.0.5, and faced with similar 
> problems.
> 

Is the total pull distance less than half the box vector along which you are 
pulling?  If not, you will be getting bizarre behavior from the pull code 
related to PBC.

Also realize that when you say you are doing constraint pulling, it will require 
some trickery to calculate the PMF after your simulations are done.  What you 
want is umbrella sampling.  You can get the PMF from constraint forces, but 
there is no automated tool to do this in Gromacs (see the list archive for a few 
posts on how you could approach that).  I'd say make your life easy and use 
umbrella sampling for the actual data collection :)

-Justin

> Kind Regards
> 
> -- 
> Ozge Engin
> ★☆
> 

-- 
========================================

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

========================================