[gmx-users] Energy minimization problem with macromolecule and polarizable water model

[Steven Kirk]

Hello,

I recently obtained (from DvS) the .itp file for the SSWM4-DP 
polarizable water molecule, and a corresponding set of water box 
coordinates.

I use double precision executables in all of the calculations mentioned 
below.

I ran an energy minimization on this box of water, then a 500ps run 
(0.001ps timestep) at 300K with Berendsen pressure and temperature 
coupling (pressure time constant set to 10 ps), PME, EnerPres dispersion 
corrections, rlist=rcoulomb=rvdv=0.9. The shell particle was left massless.

Everything seemed to be OK in the run, total/kinetic/potential energies 
settled nicely with no extreme fluctuations, density ended up about 
998kg/m^3, no major temperature or box volume fluctuations.

So then I took the final water box and used genbox to solvate a charged 
macromolecule with my shiny new water model, then added Na+ ions with 
genion until charge balance was achieved.

I disabled the defaults line in the sswm4-dp.itp file, using instead 
OPLS-AA defaults and force field parameters for everything that wasn't 
water.

Seemingly, so far so good.

The problem arose when I tried to energy minimise my solvated 
macromolecule. Using the mdp settings:
cpp                 =  /lib/cpp
define              =
constraints         =  none
integrator          =  steep
nsteps              =  2000
;
;       Energy minimizing stuff
;
emtol               =  100
emstep              =  0.01
nstcomm             =  1
ns_type             =  grid
rlist               =  1.0
coulomb-type        =  pme
rcoulomb            =  1.0
rvdw                =  1.0
Tcoupl              =  no
Pcoupl              =  no
gen_vel             =  no
lincs_iter          = 4

the minimiser ran for about 100 steps, then converged to machine 
precision. Unfortunately the max force was still very large (E+06, on a 
  water hydrogen). Reducing emstep and changing to the l-bfgs minimiser 
(cg minimiser was disallowed because of presence of constraints, 
presumably from the water model) had very little effect on the max force.

In desperation I even tried a 10ps full MD run at 10K to see if this 
would ease the mysterious force problem a little, but got LINCS errors 
on the very first step.

I'm now wondering if the problem is due to (working backwards):
a) Incompatible settings in my production system (solvated macromolecule)
b) genbox somehow mangling my equilibrated water molecules
c) Some mistake in my original equilibration of the box of polarizable 
water.

The LINCS errors suggest a structural problem, but if that arose in step 
  c), surely it would be visible as LINCS errors in the water box 
equilibration run? This makes genbox my prime suspect.

I would be very grateful for anyone's opinions on where I have gone 
wrong, and what I can do to fix the problem.

Many thanks in advance,
Steve Kirk
-- 
Dr. Steven R. Kirk           <Steven.Kirk at hv.se, S.R.Kirk at physics.org>
Dept. of Technology, Mathematics & Computer Science  (P)+46 520 223215
University West                                      (F)+46 520 223299
P.O. Box 957 Trollhattan 461 29 SWEDEN       http://taconet.webhop.org