Gromacs Drude version

User discussions
1

Hello there,

I successfully simulated a ternary system water/solvent/phosphate with a low concentration of the latter. I just wanted to do the same with a increased amount of phosphate ions (let’s say 10 times). Minimization and nvt equilibrations are fine. NPT one gives a segmentation fault, namely

step 28643: EM did not converge in 100 iterations, RMS force 10.914
step 28644: EM did not converge in 100 iterations, RMS force 130.329
step 28645: EM did not converge in 100 iterations, RMS force 435.571
step 28646: EM did not converge in 100 iterations, RMS force 834.122
step 28647: EM did not converge in 100 iterations, RMS force 1377.708
step 28648: EM did not converge in 100 iterations, RMS force 2349.287
step 28649: EM did not converge in 100 iterations, RMS force 4753.038
step 28692: EM did not converge in 100 iterations, RMS force 3.740
step 28693: EM did not converge in 100 iterations, RMS force 32.670
step 28694: EM did not converge in 100 iterations, RMS force 213.652
step 28695: EM did not converge in 100 iterations, RMS force 414.540
step 28696: EM did not converge in 100 iterations, RMS force 840.783
step 28697: EM did not converge in 100 iterations, RMS force 1445.265
step 28698: EM did not converge in 100 iterations, RMS force 2244.515
step 28699: EM did not converge in 100 iterations, RMS force 37.925

I post here the mdp. file I used:

; RUN CONTROL PARAMETERS
;define = -DPOSRES
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.001
nsteps = 10000000
comm-mode = linear
nstcomm = 100
comm-grps = System
; OUTPUT CONTROL OPTIONS
nstxout = 1000
nstvout = 1000
nstfout = 1000
; Output frequency for energies to log file and energy file
nstlog = 100
nstcalcenergy = 1
nstenergy = 100
; NEIGHBORSEARCHING PARAMETERS
cutoff-scheme = verlet
nstlist = 10
ns-type = Grid
pbc = xyz
rlist = 1.2
; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype = PME
rcoulomb = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = EnerPres
pme_order = 4
fourierspacing = 0.1
; TEMPERATURE COUPLING
tcoupl = V-rescale
; Groups to couple separately
tc-grps = System
tau-t = 0.1
ref-t = 300
; PRESSURE COUPLING IS NOT YET SUPPORTED
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau-p = 1.0
ref-p = 1.0
compressibility = 4.5e-5
refcoord-scaling = all
; GENERATE VELOCITIES FOR STARTUP RUN, OTHERWISE DO
; NOT GENERATE NEW VELOCITIES IF PREVIOUSLY EQUILIBRATED
gen-vel = no
gen-temp = 300
gen-seed = -1
; OPTIONS FOR BONDS
constraints = h-bonds ; can also be h-bonds, not all-bonds
continuation = yes
; DRUDE PARAMETERS
drude = yes
drude-mode = SCF
;drude-rhyp = yes
drude-khyp = 16736000.0
drude-r = 0.02
drude-pow = 4
emtol = 1.0
niter = 100
drude-hyper = yes

Thank you in advance
Matteo

2

Did you build, minimize, and successfully equilibrate with an additive force field first? Likely your initial conditions are simply too far from equilibrium to work.

3

Hi Justin,

You mean equilibrate through the additive ff and then like pdb2gmx with the drude ff?

4

Yes, this is the conventional approach to preparing a polarizable system. It is easy to cause instability in an additive simulation, and extraordinarily easy to break a polarizable one, if any element of it is far from equilibrium.