[gmx-users] Advice needed on choosing a VdW cutoff

[Steven Kirk]

Hello,

I would be very grateful for advice on the following system:

Consider a pair of spherical macromolecules of diameter ~2.5 nm, 
arranged along an x-axis so that their centres of mass are 4nm apart, 
then centered in a 14 x 10 x 10 nm periodic box, so that the minimum 
distance between either macromolecule's centre of mass and its closest 
simulation cell walls in the x,y, or z directions (standard orthogonal 
Cartesian axis set) is 5 nm.

(the macromolecules are silica particles with negatively charged surface 
sites and an equal number of Na+ counterions clustered around them, so 
each macromolecule-counterion 'bundle' is electrically neutral).

The box is then filled with TIP4P water, energy minimized, run with 
'all-bonds' position restraints at 300K for 100 ps, then the main MD run 
begins, reusing velocities from the last step of the position restraints 
to initialize the new run.

Here is the mdp file I am using for the main MD run:

title               =  Yo
cpp                 =  /lib/cpp
constraints         =  none
integrator          =  md
dt                  =  0.001    ; ps !
nsteps              =  1000000  ; total 1000 ps.
nstcomm             =  1
nstxout             =  10000
nstvout             =  10000
nstfout             =  0
nstlog              =  10000
nstenergy           =  10000
nstlist             =  10
ns_type             =  grid
coulombtype         = pme
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
; Berendsen temperature coupling is on in two groups
Tcoupl              =  berendsen
tc-grps             =  SNP      SOL      NA+
tau_t               =  0.1      0.1      0.1
ref_t               =  300      300      300
; Energy monitoring
energygrps          =  SNP  SOL   NA+
; Isotropic pressure coupling is not on
Pcoupl              =  no
Pcoupltype          = isotropic
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is off at 300 K.
gen_vel             =  no
gen_temp            =  300.0
gen_seed            =  173529

The problem is that when I run this simulation, the expected drift of 
the macromolecules towards each other does not occur. Assuming that I 
want to force every part of each macromolecule to 'see' every part of 
the other, this would suggest a value of rvdw of > 6.5 nm, but I have 
several worries about this:

1. I have never seen a recommended rvdw in this forum over 1.4 nm, in 
any model system
2. Should I use a standard, switched or other type of vdW cutoff?
3. Should I switch on long-range dispersion corrections (DispCorr = Ener) ?

My goal is to keep the periodic box for the advantages of PME, but 
somehow reassure myself that the macromolecules can 'see' each other via 
the vdW forces, so that they will drift together (the expected 
behaviour) over the course of the simulation.

I have trawled the mailing lists for advice on this topic - the only 
directly relevant post I could find involved the drifting apart of 
membranes over the course of a simulation, and if I remember correctly, 
the value of 'pme-order' was suggested as a culprit.

I use the default value of 'pme-order' in my simulations.

Can anyone please advise me as to what to do next? I do not want to 
abandon the investigation before I eliminate all possibilities that I am 
doing something stupid with the configuration of the vdW treatment. 
Maybe I do not have a long enough simulation, but the fact that I 
actually see *repulsion* suggests something more fundamental is wrong.

All helpful suggestions very gratefully recieved !

Regards,
Steve Kirk