Drude for GROMACS

User discussions
1

Hello there,

I’m simulating a ternary system water/tetradecane/phosphate ions both with standard Gromacs and Drude one. At a glance I’ve realized phosphate ions in Drude simulation tend to fluctuate around their position (like being constrained) unlike moving around as in the standard code, despite parameters being comparable. Below mdp parameters and topology for Drude, hoping someone help me to clarify this

Thank you in advance for your time
Matteo

; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.0005
nsteps = 200000
; 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 = 20
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
; 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 = Berendsen
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 = 298.15
gen-seed = -1
; OPTIONS FOR BONDS
;constraints = none ; can also be h-bonds, not all-bonds
;continuation = no
; 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

;
; File ‘topol.top’ was generated
; By user: frigerio (5421)
; On host: b-an01.hpc2n.umu.se
; At date: Tue Dec 8 10:24:01 2020

;
; This is a standalone topology file
;
; Created by:
; :-) GROMACS - gmx pdb2gmx, 2016-dev (-:
;
; Executable: /hpc2n/eb/software/MPI/GCC/6.4.0-2.28/impi/2017.3.196/GROMACS/2016.x-drude-20180214-g3f7439a/bin/gmx
; Data prefix: /hpc2n/eb/software/MPI/GCC/6.4.0-2.28/impi/2017.3.196/GROMACS/2016.x-drude-20180214-g3f7439a
; Command line:
; gmx pdb2gmx -f mhp.pdb -o mhp.gro
; Force field was read from current directory or a relative path - path added.
;

; Include forcefield parameters
#include “./drude-2013f_2018_10.ff/forcefield.itp”

[ moleculetype ]
; Name nrexcl
TETD 3

[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 TETD rtp TETD q 0.0
1 CD33A 1 TETD C1 1 2.30823 11.611 ; qtot 2.308
2 DRUD 1 TETD DC1 2 -2.48523 0 ; qtot -0.177
3 HDA3A 1 TETD H1 3 0.059 1.008 ; qtot -0.118
4 HDA3A 1 TETD H2 4 0.059 1.008 ; qtot -0.059
5 HDA3A 1 TETD H3 5 0.059 1.008 ; qtot 0
6 CD32A 1 TETD C2 6 2.07983 11.611 ; qtot 2.08
7 DRUD 1 TETD DC2 7 -2.23583 0 ; qtot -0.156
8 HDA2A 1 TETD H4 8 0.078 1.008 ; qtot -0.078
9 HDA2A 1 TETD H5 9 0.078 1.008 ; qtot 0
10 CD32A 1 TETD C3 10 2.07983 11.611 ; qtot 2.08
11 DRUD 1 TETD DC3 11 -2.23583 0 ; qtot -0.156
12 HDA2A 1 TETD H6 12 0.078 1.008 ; qtot -0.078
13 HDA2A 1 TETD H7 13 0.078 1.008 ; qtot 0
14 CD32A 1 TETD C4 14 2.07983 11.611 ; qtot 2.08
15 DRUD 1 TETD DC4 15 -2.23583 0 ; qtot -0.156
16 HDA2A 1 TETD H8 16 0.078 1.008 ; qtot -0.078
17 HDA2A 1 TETD H9 17 0.078 1.008 ; qtot 0
18 CD32A 1 TETD C5 18 2.07983 11.611 ; qtot 2.08
19 DRUD 1 TETD DC5 19 -2.23583 0 ; qtot -0.156
20 HDA2A 1 TETD H10 20 0.078 1.008 ; qtot -0.078
21 HDA2A 1 TETD H11 21 0.078 1.008 ; qtot 0
22 CD32A 1 TETD C6 22 2.07983 11.611 ; qtot 2.08
23 DRUD 1 TETD DC6 23 -2.23583 0 ; qtot -0.156
24 HDA2A 1 TETD H12 24 0.078 1.008 ; qtot -0.078
25 HDA2A 1 TETD H13 25 0.078 1.008 ; qtot 0
26 CD32A 1 TETD C7 26 2.07983 11.611 ; qtot 2.08
27 DRUD 1 TETD DC7 27 -2.23583 0 ; qtot -0.156
28 HDA2A 1 TETD H14 28 0.078 1.008 ; qtot -0.078
29 HDA2A 1 TETD H15 29 0.078 1.008 ; qtot 0
30 CD32A 1 TETD C8 30 2.07983 11.611 ; qtot 2.08
31 DRUD 1 TETD DC8 31 -2.23583 0 ; qtot -0.156
32 HDA2A 1 TETD H16 32 0.078 1.008 ; qtot -0.078
33 HDA2A 1 TETD H17 33 0.078 1.008 ; qtot 0
34 CD32A 1 TETD C9 34 2.07983 11.611 ; qtot 2.08
35 DRUD 1 TETD DC9 35 -2.23583 0 ; qtot -0.156
36 HDA2A 1 TETD H18 36 0.078 1.008 ; qtot -0.078
37 HDA2A 1 TETD H19 37 0.078 1.008 ; qtot 0
38 CD32A 1 TETD C10 38 2.07983 11.611 ; qtot 2.08
39 DRUD 1 TETD DC10 39 -2.23583 0 ; qtot -0.156
40 HDA2A 1 TETD H20 40 0.078 1.008 ; qtot -0.078
41 HDA2A 1 TETD H21 41 0.078 1.008 ; qtot 0
42 CD32A 1 TETD C11 42 2.07983 11.611 ; qtot 2.08
43 DRUD 1 TETD DC11 43 -2.23583 0 ; qtot -0.156
44 HDA2A 1 TETD H22 44 0.078 1.008 ; qtot -0.078
45 HDA2A 1 TETD H23 45 0.078 1.008 ; qtot 0
46 CD32A 1 TETD C12 46 2.07983 11.611 ; qtot 2.08
47 DRUD 1 TETD DC12 47 -2.23583 0 ; qtot -0.156
48 HDA2A 1 TETD H24 48 0.078 1.008 ; qtot -0.078
49 HDA2A 1 TETD H25 49 0.078 1.008 ; qtot 0
50 CD32A 1 TETD C13 50 2.07983 11.611 ; qtot 2.08
51 DRUD 1 TETD DC13 51 -2.23583 0 ; qtot -0.156
52 HDA2A 1 TETD H26 52 0.078 1.008 ; qtot -0.078
53 HDA2A 1 TETD H27 53 0.078 1.008 ; qtot 0
54 CD33A 1 TETD C14 54 2.30823 11.611 ; qtot 2.308
55 DRUD 1 TETD DC14 55 -2.48523 0 ; qtot -0.177
56 HDA3A 1 TETD H28 56 0.059 1.008 ; qtot -0.118
57 HDA3A 1 TETD H29 57 0.059 1.008 ; qtot -0.059
58 HDA3A 1 TETD H30 58 0.059 1.008 ; qtot 0

[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
1 6 1
6 7 1
6 8 1
6 9 1
6 10 1
10 11 1
10 12 1
10 13 1
10 14 1
14 15 1
14 16 1
14 17 1
14 18 1
18 19 1
18 20 1
18 21 1
18 22 1
22 23 1
22 24 1
22 25 1
22 26 1
26 27 1
26 28 1
26 29 1
26 30 1
30 31 1
30 32 1
30 33 1
30 34 1
34 35 1
34 36 1
34 37 1
34 38 1
38 39 1
38 40 1
38 41 1
38 42 1
42 43 1
42 44 1
42 45 1
42 46 1
46 47 1
46 48 1
46 49 1
46 50 1
50 51 1
50 52 1
50 53 1
50 54 1
54 55 1
54 56 1
54 57 1
54 58 1

[ pairs ]
; ai aj funct c0 c1 c2 c3
1 12 1
1 13 1
1 14 1
1 15 1
2 15 1
3 10 1
3 11 1
4 10 1
4 11 1
5 10 1
5 11 1
6 16 1
6 17 1
6 18 1
6 19 1
7 19 1
8 14 1
8 15 1
9 14 1
9 15 1
10 20 1
10 21 1
10 22 1
10 23 1
11 23 1
12 2 1
12 18 1
12 19 1
13 2 1
13 18 1
13 19 1
14 2 1
14 24 1
14 25 1
14 26 1
14 27 1
15 27 1
16 7 1
16 22 1
16 23 1
17 7 1
17 22 1
17 23 1
18 7 1
18 28 1
18 29 1
18 30 1
18 31 1
19 31 1
20 11 1
20 26 1
20 27 1
21 11 1
21 26 1
21 27 1
22 11 1
22 32 1
22 33 1
22 34 1
22 35 1
23 35 1
24 15 1
24 30 1
24 31 1
25 15 1
25 30 1
25 31 1
26 15 1
26 36 1
26 37 1
26 38 1
26 39 1
27 39 1
28 19 1
28 34 1
28 35 1
29 19 1
29 34 1
29 35 1
30 19 1
30 40 1
30 41 1
30 42 1
30 43 1
31 43 1
32 23 1
32 38 1
32 39 1
33 23 1
33 38 1
33 39 1
34 23 1
34 44 1
34 45 1
34 46 1
34 47 1
35 47 1
36 27 1
36 42 1
36 43 1
37 27 1
37 42 1
37 43 1
38 27 1
38 48 1
38 49 1
38 50 1
38 51 1
39 51 1
40 31 1
40 46 1
40 47 1
41 31 1
41 46 1
41 47 1
42 31 1
42 52 1
42 53 1
42 54 1
42 55 1
43 55 1
44 35 1
44 50 1
44 51 1
45 35 1
45 50 1
45 51 1
46 35 1
46 56 1
46 57 1
46 58 1
48 39 1
48 54 1
48 55 1
49 39 1
49 54 1
49 55 1
50 39 1
52 43 1
53 43 1
54 43 1
56 47 1
57 47 1
58 47 1

[ exclusions ]
; i excluded from i
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
2 1 3 4 5 6 7 8 9 10 11 12 13 14 15
3 1 2 4 5 6 7 8 9 10 11
4 1 2 3 5 6 7 8 9 10 11
5 1 2 3 4 6 7 8 9 10 11
6 1 2 3 4 5 7 8 9 10 11 12 13 14 15 16 17 18 19
7 1 2 3 4 5 6 8 9 10 11 12 13 14 15 16 17 18 19
8 1 2 3 4 5 6 7 9 10 11 12 13 14 15
9 1 2 3 4 5 6 7 8 10 11 12 13 14 15
10 1 2 3 4 5 6 7 8 9 11 12 13 14 15 16 17 18 19 20 21 22 23
11 1 2 6 7 8 9 10 12 13 14 15 16 17 18 19 20 21 22 23
12 1 6 7 8 9 10 11 13 14 15 16 17 18 19
13 1 6 7 8 9 10 11 12 14 15 16 17 18 19
14 1 6 7 8 9 10 11 12 13 15 16 17 18 19 20 21 22 23 24 25 26 27
15 2 6 7 10 11 12 13 14 16 17 18 19 20 21 22 23 24 25 26 27
16 6 10 11 12 13 14 15 17 18 19 20 21 22 23
17 6 10 11 12 13 14 15 16 18 19 20 21 22 23
18 6 10 11 12 13 14 15 16 17 19 20 21 22 23 24 25 26 27 28 29 30 31
19 7 10 11 14 15 16 17 18 20 21 22 23 24 25 26 27 28 29 30 31
20 10 14 15 16 17 18 19 21 22 23 24 25 26 27
21 10 14 15 16 17 18 19 20 22 23 24 25 26 27
22 10 14 15 16 17 18 19 20 21 23 24 25 26 27 28 29 30 31 32 33 34 35
23 11 14 15 18 19 20 21 22 24 25 26 27 28 29 30 31 32 33 34 35
24 14 18 19 20 21 22 23 25 26 27 28 29 30 31
25 14 18 19 20 21 22 23 24 26 27 28 29 30 31
26 14 18 19 20 21 22 23 24 25 27 28 29 30 31 32 33 34 35 36 37 38 39
27 15 18 19 22 23 24 25 26 28 29 30 31 32 33 34 35 36 37 38 39
28 18 22 23 24 25 26 27 29 30 31 32 33 34 35
29 18 22 23 24 25 26 27 28 30 31 32 33 34 35
30 18 22 23 24 25 26 27 28 29 31 32 33 34 35 36 37 38 39 40 41 42 43
31 19 22 23 26 27 28 29 30 32 33 34 35 36 37 38 39 40 41 42 43
32 22 26 27 28 29 30 31 33 34 35 36 37 38 39
33 22 26 27 28 29 30 31 32 34 35 36 37 38 39
34 22 26 27 28 29 30 31 32 33 35 36 37 38 39 40 41 42 43 44 45 46 47
35 23 26 27 30 31 32 33 34 36 37 38 39 40 41 42 43 44 45 46 47
36 26 30 31 32 33 34 35 37 38 39 40 41 42 43
37 26 30 31 32 33 34 35 36 38 39 40 41 42 43
38 26 30 31 32 33 34 35 36 37 39 40 41 42 43 44 45 46 47 48 49 50 51
39 27 30 31 34 35 36 37 38 40 41 42 43 44 45 46 47 48 49 50 51
40 30 34 35 36 37 38 39 41 42 43 44 45 46 47
41 30 34 35 36 37 38 39 40 42 43 44 45 46 47
42 30 34 35 36 37 38 39 40 41 43 44 45 46 47 48 49 50 51 52 53 54 55
43 31 34 35 38 39 40 41 42 44 45 46 47 48 49 50 51 52 53 54 55
44 34 38 39 40 41 42 43 45 46 47 48 49 50 51
45 34 38 39 40 41 42 43 44 46 47 48 49 50 51
46 34 38 39 40 41 42 43 44 45 47 48 49 50 51 52 53 54 55 56 57 58
47 35 38 39 42 43 44 45 46 48 49 50 51 52 53 54 55 56 57 58
48 38 42 43 44 45 46 47 49 50 51 52 53 54 55
49 38 42 43 44 45 46 47 48 50 51 52 53 54 55
50 38 42 43 44 45 46 47 48 49 51 52 53 54 55 56 57 58
51 39 42 43 46 47 48 49 50 52 53 54 55 56 57 58
52 42 46 47 48 49 50 51 53 54 55 56 57 58
53 42 46 47 48 49 50 51 52 54 55 56 57 58
54 42 46 47 48 49 50 51 52 53 55 56 57 58
55 43 46 47 50 51 52 53 54 56 57 58
56 46 50 51 52 53 54 55 57 58
57 46 50 51 52 53 54 55 56 58
58 46 50 51 52 53 54 55 56 57

[ angles ]
; ai aj ak funct c0 c1 c2 c3
3 1 4 5
3 1 5 5
3 1 6 5
4 1 5 5
4 1 6 5
5 1 6 5
1 6 8 5
1 6 9 5
1 6 10 5
8 6 9 5
8 6 10 5
9 6 10 5
6 10 12 5
6 10 13 5
6 10 14 5
12 10 13 5
12 10 14 5
13 10 14 5
10 14 16 5
10 14 17 5
10 14 18 5
16 14 17 5
16 14 18 5
17 14 18 5
14 18 20 5
14 18 21 5
14 18 22 5
20 18 21 5
20 18 22 5
21 18 22 5
18 22 24 5
18 22 25 5
18 22 26 5
24 22 25 5
24 22 26 5
25 22 26 5
22 26 28 5
22 26 29 5
22 26 30 5
28 26 29 5
28 26 30 5
29 26 30 5
26 30 32 5
26 30 33 5
26 30 34 5
32 30 33 5
32 30 34 5
33 30 34 5
30 34 36 5
30 34 37 5
30 34 38 5
36 34 37 5
36 34 38 5
37 34 38 5
34 38 40 5
34 38 41 5
34 38 42 5
40 38 41 5
40 38 42 5
41 38 42 5
38 42 44 5
38 42 45 5
38 42 46 5
44 42 45 5
44 42 46 5
45 42 46 5
42 46 48 5
42 46 49 5
42 46 50 5
48 46 49 5
48 46 50 5
49 46 50 5
46 50 52 5
46 50 53 5
46 50 54 5
52 50 53 5
52 50 54 5
53 50 54 5
50 54 56 5
50 54 57 5
50 54 58 5
56 54 57 5
56 54 58 5
57 54 58 5

[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
3 1 6 8 9
3 1 6 9 9
3 1 6 10 9
4 1 6 8 9
4 1 6 9 9
4 1 6 10 9
5 1 6 8 9
5 1 6 9 9
5 1 6 10 9
1 6 10 12 9
1 6 10 13 9
1 6 10 14 9
8 6 10 12 9
8 6 10 13 9
8 6 10 14 9
9 6 10 12 9
9 6 10 13 9
9 6 10 14 9
6 10 14 16 9
6 10 14 17 9
6 10 14 18 9
12 10 14 16 9
12 10 14 17 9
12 10 14 18 9
13 10 14 16 9
13 10 14 17 9
13 10 14 18 9
10 14 18 20 9
10 14 18 21 9
10 14 18 22 9
16 14 18 20 9
16 14 18 21 9
16 14 18 22 9
17 14 18 20 9
17 14 18 21 9
17 14 18 22 9
14 18 22 24 9
14 18 22 25 9
14 18 22 26 9
20 18 22 24 9
20 18 22 25 9
20 18 22 26 9
21 18 22 24 9
21 18 22 25 9
21 18 22 26 9
18 22 26 28 9
18 22 26 29 9
18 22 26 30 9
24 22 26 28 9
24 22 26 29 9
24 22 26 30 9
25 22 26 28 9
25 22 26 29 9
25 22 26 30 9
22 26 30 32 9
22 26 30 33 9
22 26 30 34 9
28 26 30 32 9
28 26 30 33 9
28 26 30 34 9
29 26 30 32 9
29 26 30 33 9
29 26 30 34 9
26 30 34 36 9
26 30 34 37 9
26 30 34 38 9
32 30 34 36 9
32 30 34 37 9
32 30 34 38 9
33 30 34 36 9
33 30 34 37 9
33 30 34 38 9
30 34 38 40 9
30 34 38 41 9
30 34 38 42 9
36 34 38 40 9
36 34 38 41 9
36 34 38 42 9
37 34 38 40 9
37 34 38 41 9
37 34 38 42 9
34 38 42 44 9
34 38 42 45 9
34 38 42 46 9
40 38 42 44 9
40 38 42 45 9
40 38 42 46 9
41 38 42 44 9
41 38 42 45 9
41 38 42 46 9
38 42 46 48 9
38 42 46 49 9
38 42 46 50 9
44 42 46 48 9
44 42 46 49 9
44 42 46 50 9
45 42 46 48 9
45 42 46 49 9
45 42 46 50 9
42 46 50 52 9
42 46 50 53 9
42 46 50 54 9
48 46 50 52 9
48 46 50 53 9
48 46 50 54 9
49 46 50 52 9
49 46 50 53 9
49 46 50 54 9
46 50 54 56 9
46 50 54 57 9
46 50 54 58 9
52 50 54 56 9
52 50 54 57 9
52 50 54 58 9
53 50 54 56 9
53 50 54 57 9
53 50 54 58 9

[ thole_polarization ]
; ai aj ak al funct c0 c1 c2 c3
1 2 6 7 2 0.002051 0.001660 1.3000 1.3000
1 2 10 11 2 0.002051 0.001660 1.3000 1.3000
6 7 10 11 2 0.001660 0.001660 1.3000 1.3000
6 7 14 15 2 0.001660 0.001660 1.3000 1.3000
10 11 14 15 2 0.001660 0.001660 1.3000 1.3000
10 11 18 19 2 0.001660 0.001660 1.3000 1.3000
14 15 18 19 2 0.001660 0.001660 1.3000 1.3000
14 15 22 23 2 0.001660 0.001660 1.3000 1.3000
18 19 22 23 2 0.001660 0.001660 1.3000 1.3000
18 19 26 27 2 0.001660 0.001660 1.3000 1.3000
22 23 26 27 2 0.001660 0.001660 1.3000 1.3000
22 23 30 31 2 0.001660 0.001660 1.3000 1.3000
26 27 30 31 2 0.001660 0.001660 1.3000 1.3000
26 27 34 35 2 0.001660 0.001660 1.3000 1.3000
30 31 34 35 2 0.001660 0.001660 1.3000 1.3000
30 31 38 39 2 0.001660 0.001660 1.3000 1.3000
34 35 38 39 2 0.001660 0.001660 1.3000 1.3000
34 35 42 43 2 0.001660 0.001660 1.3000 1.3000
38 39 42 43 2 0.001660 0.001660 1.3000 1.3000
38 39 46 47 2 0.001660 0.001660 1.3000 1.3000
42 43 46 47 2 0.001660 0.001660 1.3000 1.3000
42 43 50 51 2 0.001660 0.001660 1.3000 1.3000
46 47 50 51 2 0.001660 0.001660 1.3000 1.3000
46 47 54 55 2 0.001660 0.002051 1.3000 1.3000
50 51 54 55 2 0.001660 0.002051 1.3000 1.3000

[ moleculetype ]
; Name nrexcl
MHP 3

[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 MHP rtp MHP q -2.0
1 OD31D 1 MHP O1 1 1.97173 15.599 ; qtot 1.972
2 DRUD 1 MHP DO1 2 -1.97173 0 ; qtot 0
3 LPDNA1 1 MHP LPA 3 -0.345 0 ; qtot -0.345
4 LPDNA1 1 MHP LPB 4 -0.345 0 ; qtot -0.69
5 PD1AN 1 MHP P1 5 3.66551 30.574 ; qtot 2.976
6 DRUD 1 MHP DP1 6 -1.93551 0 ; qtot 1.04
7 OD2C2C 1 MHP O2 7 0.60529 15.5994 ; qtot 1.645
8 DRUD 1 MHP DO2 8 -1.69229 0 ; qtot -0.047
9 OD2C2C 1 MHP O3 9 0.60529 15.5994 ; qtot 0.5583
10 DRUD 1 MHP DO3 10 -1.69229 0 ; qtot -1.134
11 OD2C2C 1 MHP O4 11 0.60529 15.5994 ; qtot -0.5287
12 DRUD 1 MHP DO4 12 -1.69229 0 ; qtot -2.221
13 HDP1A 1 MHP H 13 0.221 1.008 ; qtot -2

[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 5 1
1 13 1
5 6 1
5 7 1
5 9 1
7 8 1
9 10 1
11 12 1

[ pairs ]
ai aj funct c0 c1 c2 c3
7 13 1
9 13 1
13 8 1
13 10 1

[ exclusions ]
i excluded from i
1 2 3 4 5 6 7 8 9 10 13
2 1 3 4 5 6 7 8 9 10 13
3 4
4 3
5 1 2 3 4 6 7 8 9 10 13
6 1 2 3 4 5 7 8 9 10 13
7 1 2 3 4 5 6 8 9 10 13
8 1 2 3 4 5 6 7 9 10 13
9 1 2 3 4 5 6 7 8 10 13
10 1 2 3 4 5 6 7 8 9 13
11 12
12 11
13 1 2 3 4 5 6 7 9

[ angles ]
; ai aj ak funct c0 c1 c2 c3
5 1 13 5
1 5 7 5
1 5 9 5
7 5 9 5

[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
13 1 5 7 9
13 1 5 9 9

[ thole_polarization ]
; ai aj ak al funct c0 c1 c2 c3
1 2 5 6 2 0.001291 0.001244 1.0910 1.6780
1 2 7 8 2 0.001291 0.000951 1.0910 1.0830
1 2 9 10 2 0.001291 0.000951 1.0910 1.0830
5 6 7 8 2 0.001244 0.000951 1.6780 1.0830
5 6 9 10 2 0.001244 0.000951 1.6780 1.0830
7 8 9 10 2 0.000951 0.000951 1.0830 1.0830

[ virtual_sites3 ]
; ai aj ak al funct c0 c1
3 1 5 13 3 110.00 0.035
4 1 5 13 3 110.00 0.035

; Include Position restraint file
#ifdef POSRES
#include “posre.itp”
#endif

; Include water topology
#include “./drude-2013f_2018_10.ff/swm4-ndp.itp”

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif

; Include topology for ions
#include “./drude-2013f_2018_10.ff/ions.itp”

[ system ]
; Name
Title in water

[ molecules ]
; Compound #mols
TETD 158
MHP 5
SOL 1434
NA 10

2

How did you generate these topologies?

Slow diffusion is possible if the phosphate ions are perhaps strongly attracted to something nearby. But there are no obvious restraints at work here, so nothing is preventing them from moving.

The choice of no constraints and dt = 0.5 fs is odd. Why have you done this? Our Drude model constrains bonds to H and uses dt = 1 fs.

3

Hi Justin,

How did you generate these topologies?

Phosphate ion is from Kognole’s ‘Balanced polarizable Drude force field parameters for molecular anions: phosphates, sulfates, sulfamates, and oxides’ I have adapted for GROMACS. Tetradecane instead is a slight modification of parameters for Pentadecane in the .rtp file.

Slow diffusion is possible if the phosphate ions are perhaps strongly attracted to something nearby. But there are no obvious restraints at work here, so nothing is preventing them from moving.

I said so due to the fact that that simulations with CHARMM additive forcefield highlight fast diffusion and the generation of a phosphate cluster over time.

The choice of no constraints and dt = 0.5 fs is odd. Why have you done this? Our Drude model constrains bonds to H and uses dt = 1 fs.

Due to this:
The bond in molecule-type MHP between atoms 1 O1 and 13 H has an estimated oscillational period of 9.1e-03 ps, which is less than 10 times the time step of 1.0e-03 ps.
Maybe you forgot to change the constraints mdp option.

But I guess I can constraint bonds to H and solve it.

4

OK, that’s why I asked - there have been considerable changes to the force field as a result of this work, including new LJ parameters. You’re hacking topologies into the old 2013f force field, but the new phosphates, sulfates, etc. are dependent upon a new, developmental parameter set that has not been adapted for GROMACS yet. Have you made modifications to ffnonbonded.itp and ffbonded.itp to update the relevant atom types that have been changed?

5

So I have to transpose bonded and non bonded parameters from CHARMM drude_toppar_2019.tgz to GROMACS drude-2013f_2018_10.tgz?

6

No, you need the parameters from the Kognole paper (and SI). These are not yet publicly available in CHARMM format as they are part of our developmental repository, but everything you need is described in the paper.

7

Ok, Thank you very much for your time!

8

Hi Justin,

I have modified the non-bonded parameters in agreement with drude_toppar_2019.tgz, but my simulation still suffers from slow diffusion (for an equal long simulation using additive charmm I have reasonable diffusion coefficients). From the Kognole paper only anionic oxygen parameters are explicitly cited, but maybe also other parameters are affected. Should I look at the interaction energies in the SI and retrieve parameters from those?

9

You can check the interaction energies but that’s not something easily done in GROMACS (quite straightforward in CHARMM), but I’d suggest you just inquire with Alex’s group to see if they’ve looked at diffusion, etc. at all in their parametrization. It may not have been a relevant quantity for them, especially with the goal of using these molecules as building blocks of biomolecules, not necessarily used individually.

10

Thank you for the advice.
One more thing: whenever I try to use anisotropic polarization in my topology I have a segmentation fault during mdrun.

Steepest Descents:
Tolerance (Fmax) = 1.00000e+03
Number of steps = -1
Step= 0, Dmax= 1.0e-02 nm, Epot= -nan Fmax= 5.76182e+06, atom= 63
Segmentation fault

11

I can’t provide any useful comment without seeing the relevant topology section. Conversion of the anisotropy topology entries is rather complicated, so if it’s failing, you’ve probably made a mistake somewhere along the way, but again I can’t say.

12

Yeah sorry, so here the .rtp entry and to follow the top file.

[ MHP ]
[ atoms ]
; name type charge cgnr alpha thole
O1 OD31D 1.97173 0 0.001291 1.091
DO1 DRUD -1.97173 1
P1 PD1AN 3.66551 2 0.001244 1.678
DP1 DRUD -1.93551 3
O2 OD2C2C 0.60529 4 0.000951 1.083
DO2 DRUD -1.69229 5
O3 OD2C2C 0.60529 6 0.000951 1.083
DO3 DRUD -1.69229 7
O4 OD2C2C 0.60529 8 0.000951 1.083
DO4 DRUD -1.69229 9
LPA LPDNA1 -0.345 10
LPB LPDNA1 -0.345 11
H HDP1A 0.221 12

[ bonds ]
O1 DO1
P1 DP1
O2 DO2
O3 DO3
O4 DO4
P1 O1
P1 O2
P1 O3
P1 04
O1 H

[ anisotropic_polarization ]
; drude ai aj ak al a11 a22 a33
DO1 O1 P1 LPA LPB 0.76473 1.16239 1.07288

[ virtual_sites3 ]
; 3fad
; site from func theta d
LPA O1 P1 H 3 110.00 0.035
LPB O1 P1 H 3 110.00 0.035

;
; File ‘topol.top’ was generated
; By user: frigerio (5421)
; On host: b-an01.hpc2n.umu.se
; At date: Tue Jan 26 11:47:17 2021

;
; This is a standalone topology file
;
; Created by:
; :-) GROMACS - gmx pdb2gmx, 2016-dev (-:
;
; Executable: /hpc2n/eb/software/MPI/GCC/6.4.0-2.28/OpenMPI/2.1.1/GROMACS/2016.x-drude-20180214-g3f7439a/bin/gmx
; Data prefix: /hpc2n/eb/software/MPI/GCC/6.4.0-2.28/OpenMPI/2.1.1/GROMACS/2016.x-drude-20180214-g3f7439a
; Command line:
; gmx pdb2gmx -f input.gro -o mhp.gro
; Force field was read from current directory or a relative path - path added.
;

; Include forcefield parameters
#include “./drude-2013f_2018_10.ff/forcefield.itp”

[ moleculetype ]
; Name nrexcl
MHP 3

[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 MHP rtp MHP q -2.0
1 OD31D 1 MHP O1 1 1.97173 15.599 ; qtot 1.972
2 DRUD 1 MHP DO1 2 -1.97173 0 ; qtot 0
3 LPDNA1 1 MHP LPA 3 -0.345 0 ; qtot -0.345
4 LPDNA1 1 MHP LPB 4 -0.345 0 ; qtot -0.69
5 PD1AN 1 MHP P1 5 3.66551 30.574 ; qtot 2.976
6 DRUD 1 MHP DP1 6 -1.93551 0 ; qtot 1.04
7 OD2C2C 1 MHP O2 7 0.60529 15.5994 ; qtot 1.645
8 DRUD 1 MHP DO2 8 -1.69229 0 ; qtot -0.047
9 OD2C2C 1 MHP O3 9 0.60529 15.5994 ; qtot 0.5583
10 DRUD 1 MHP DO3 10 -1.69229 0 ; qtot -1.134
11 OD2C2C 1 MHP O4 11 0.60529 15.5994 ; qtot -0.5287
12 DRUD 1 MHP DO4 12 -1.69229 0 ; qtot -2.221
13 HDP1A 1 MHP H 13 0.221 1.008 ; qtot -2

[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 5 1
1 13 1
5 6 1
5 7 1
5 9 1
7 8 1
9 10 1
11 12 1

[ pairs ]
; ai aj funct c0 c1 c2 c3
7 13 1
9 13 1
13 8 1
13 10 1

[ exclusions ]
; i excluded from i
1 2 3 4 5 6 7 8 9 10 13
2 1 3 4 5 6 7 8 9 10 13
3 4
4 3
5 1 2 3 4 6 7 8 9 10 13
6 1 2 3 4 5 7 8 9 10 13
7 1 2 3 4 5 6 8 9 10 13
8 1 2 3 4 5 6 7 9 10 13
9 1 2 3 4 5 6 7 8 10 13
10 1 2 3 4 5 6 7 8 9 13
11 12
12 11
13 1 2 3 4 5 6 7 9

[ angles ]
; ai aj ak funct c0 c1 c2 c3
5 1 13 5
1 5 7 5
1 5 9 5
7 5 9 5

[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
13 1 5 7 9
13 1 5 9 9

[ anisotropic_polarization ]
; ai aj ak al am funct c0 c1 c2
2 1 5 3 4 2 0.76473 1.16239 1.07288

[ thole_polarization ]
; ai aj ak al funct c0 c1 c2 c3
1 2 5 6 2 0.001291 0.001244 1.0910 1.6780
1 2 7 8 2 0.001291 0.000951 1.0910 1.0830
1 2 9 10 2 0.001291 0.000951 1.0910 1.0830
5 6 7 8 2 0.001244 0.000951 1.6780 1.0830
5 6 9 10 2 0.001244 0.000951 1.6780 1.0830
7 8 9 10 2 0.000951 0.000951 1.0830 1.0830

[ virtual_sites3 ]
; ai aj ak al funct c0 c1
3 1 5 13 3 110.00 0.035
4 1 5 13 3 110.00 0.035

; Include Position restraint file
#ifdef POSRES
#include “posre.itp”
#endif

; Include water topology
#include “./drude-2013f_2018_10.ff/swm4-ndp.itp”

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif

; Include topology for ions
#include “./drude-2013f_2018_10.ff/ions.itp”

[ system ]
; Name
Title

[ molecules ]
; Compound #mols
MHP 8
SOL 4111
NA 16

13

Here’s your problem. You’ve got a 0 (zero) instead of O (letter) in your .rtp entry, so the whole topology is wrong, because you’re missing the P1-O4 bond and all subsequent related interactions. Anisotropy looks fine, probably just random chance that it led to some instability in conjunction with the incorrect bonded structure.

14

Thank you Justin, but unfortunately the error remains.
It has to do with anisotropic polarization as once I comment out those lines mdrun goes fine.
Maybe a formatting error, cause once again mdrun cannot compute potential energy.

15

If you can send me a .tpr file, I can try to look into this, but I’m not sure how quickly I’ll be able to diagnose it. The problem should not be in the anisotropy; I have tested that very thoroughly and everything matches CHARMM forces and energies. Likely the anisotropy is exposing something else that is wrong in the topology.

16

Yeah, no problem.
I’ll send it at jalemkul@vt.edu as gromacs forum does not allow to upload it here.
Take your time and thank again for your support.

Cheers
Matteo

17

Your molecules do not have proper coordinates. The lone pairs are on top of one another, which leads to an undefined force. You need to have sensible coordinates to start. I also recommend restraining all real atoms for the first round of energy minimization, then minimizing the entire system without restraints if needed.

18

Ok, I have simply equilibrated a hydrogen phosphate molecule in CHARMM and then pdb2gmxed it in Drude Gromacs. How can I do that? Maybe with Drude Prepper?
I don’t think it is possible the latter cause the hydrogen phosphate topology is not included yet in drude charmm version

19

I see it now. You’re defining the lone pairs to be in the same place, and also using the incorrect geometry. Tetrahedral lone pairs have to use the 3out geometry, because they also effectively rely on a relative orientation that is governed by a dihedral (as defined in the CHARMM topology).

Here is the correct .rtp entry:

[ HP_2 ]
  [ atoms ]
  ;  name   type     charge cgnr     alpha    thole
       P1   PD1A    3.66551   0   0.001244   1.6780
      DP1   DRUD   -1.93551   1
       O2 OD2C2B    0.60529   2   0.000951   1.0830
      DO2   DRUD   -1.69229   3
       O3 OD2C2B    0.60529   4   0.000951   1.0830
      DO3   DRUD   -1.69229   5
       O4 OD2C2B    0.60529   6   0.000951   1.0830
      DO4   DRUD   -1.69229   7
       O1  OD31D    1.97173   8   0.001291   1.0910
      DO1   DRUD   -1.97173   9
      LPA LPDNA1   -0.34500  10
      LPB LPDNA1   -0.34500  11
        H  HDP1A    0.22100  12
  [ bonds ]
       P1   DP1
       O2   DO2
       O3   DO3
       O4   DO4
       O1   DO1
       P1    O1
       P1    O2
       P1    O3
       P1    O4
       O1     H
  [ virtual_sites3 ]
  ; 3out
  ;  site            from   func            a            b            c
      LPA    O1    P1     H      4 -0.0893945185 -0.0094584498 -2.4599354697
      LPB    O1    P1     H      4 -0.0893945185 -0.0094584498 2.4599354697

I should note that I pulled this from our internal (development) version of the force field, which has anisotropy commented out. There is a comment to “check on this” so please inquire with Alex MacKerell’s group about whether the O1 anisotropy should indeed be added to this .rtp entry and thus the topology.

20

Thank you for everything.

Have a nice weekend
Matteo