[gmx-users] Simulations with Verlet scheme, virtual sites and heterogeneous systems

[Ivan Gladich]

Dear Gromacs users,
         I am experiencing a strange behavior using Verlet-cutoff scheme 
in combination with a water model with virtual sites.
I am simulating an ice slab with NE6 water model. NE6 water model is 
used for simulation of ice and it consists of 3 atom sites and 3 virtual 
sites.
My system has two vacuum interfaces as drawn below.

------------|------------|---------
vacuum  | ice         |  vacuum
                |               |
-------------------------------

My simulation run at 240 K and, at this temperature, my ice should not 
melt (for this water model Tm=289 K).
Using a group scheme , my ice remains stable over 10 ns

cutoff-scheme   = group
ns_type             = grid
nstlist                 = 1
rlist                   = 1
pbc                    = xyz
rcoulomb                 = 1.0
vdw-type                 = Cut-Off
rvdw                     = 1.0

*
**On the contrary*, using a Verlet cut-off scheme

cutoff-scheme   = Verlet
ns_type         = grid
nstlist         = 10
pbc             = xyz
vdw-type                 = Cut-Off
rvdw                     = 1.0

I start to observe that
1) sometimes  my ice suddenly (not in a realist way) melts after ~7 ns 
and my simulation crashes
OR
2) some particle are kicked out with high velocity and my run crashes 
with this message

-------------------------------------------------------
Program mdrun_4_7_gpu_s, VERSION 4.6.7
Source code file: /home/igladich/opt/gromacs-4.6.7/src/mdlib/pme.c, 
line: 851

Fatal error:
2 particles communicated to PME node 0 are more than 2/3 times the 
cut-off out of the domain decomposition cell of their charge group in 
dimension y.
This usually means that your system is not well equilibrated.
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------


I am uploading the temperature profile and the *.mdp file at

http://goo.gl/0zfqfe

So, my questions concern the Verlet scheme
1) Can  this scheme  properly work in the presence of an heterogeneous 
system as mine (with an interface on vacuum)?
2) Can the problem be the virtual sites in  my water models? Because I 
saw, running  grompp the following note...

NOTE 2 [file grompp-RUN-NE6-4_7.mdp]:
   There are 5760 non-linear virtual site constructions. Their contribution
   to the energy drift is approximated. In most cases this does not affect
   the energy drift significantly.

3) Is there a way to go around this issue?  With GPU+Verlet and I doing 
~110 ns/day while, with group+CPU 33 ns/day ....

Thanks for any possible help or suggestions
Sincerely
Ivan


-- 
------
Ivan Gladich, Ph.D.
Research Assistant
International School for Advanced Studies (SISSA)
Via Bonomea 265, I-34136, Trieste,
Italy

Room: 331
Tel: +390403787335
e-mail:igladich at sissa.it
web page:http://people.sissa.it/~igladich/