# Drude for GROMACS **URL:** \ **Category:** User discussions\ **Created:** [January 12, 2021, 5:27pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430 "2021-01-12T17:27:24Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 12, 2021, 5:27pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/1 "2021-01-12T17:27:24Z")
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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 13, 2021, 12:05pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/2 "2021-01-13T12:05:09Z")
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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 13, 2021, 2:39pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/3 "2021-01-13T14:39:32Z")
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. ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 13, 2021, 3:05pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/4 "2021-01-13T15:05:55Z")
> [@bellafrizzo](#): > > 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. 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? ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 13, 2021, 3:28pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/5 "2021-01-13T15:28:28Z")
So I have to transpose bonded and non bonded parameters from CHARMM drude\_toppar\_2019.tgz to GROMACS drude-2013f\_2018\_10.tgz? ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 13, 2021, 3:44pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/6 "2021-01-13T15:44:37Z")
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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 13, 2021, 3:56pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/7 "2021-01-13T15:56:17Z")
Ok, Thank you very much for your time! ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 25, 2021, 4:01pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/8 "2021-01-25T16:01:02Z")
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? ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 25, 2021, 8:46pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/9 "2021-01-25T20:46:45Z")
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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 26, 2021, 10:55am UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/10 "2021-01-26T10:55:36Z")
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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 26, 2021, 2:02pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/11 "2021-01-26T14:02:04Z")
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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 26, 2021, 3:15pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/12 "2021-01-26T15:15:50Z")
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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 27, 2021, 1:57pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/13 "2021-01-27T13:57:27Z")
> [@bellafrizzo](#): > > P1 04 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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 28, 2021, 1:52pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/14 "2021-01-28T13:52:28Z")
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. ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 29, 2021, 1:50pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/15 "2021-01-29T13:50:22Z")
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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 29, 2021, 2:12pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/16 "2021-01-29T14:12:35Z")
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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 29, 2021, 2:23pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/17 "2021-01-29T14:23:41Z")
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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 29, 2021, 2:37pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/18 "2021-01-29T14:37:41Z")
Ok, I have simply equilibrated a hydrogen phosphate molecule in CHARMM and then _pdb2gmx_ed 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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 29, 2021, 3:23pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/19 "2021-01-29T15:23:30Z")
> [@bellafrizzo](#): > > [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 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. ---
**Author:** ![bellafrizzo](https://avatars.discourse-cdn.com/v4/letter/b/b487fb/32.png) [@bellafrizzo](https://gromacs.bioexcel.eu/u/bellafrizzo)\ **Post date:** [January 29, 2021, 4:05pm UTC](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430/20 "2021-01-29T16:05:38Z")
Thank you for everything. Have a nice weekend Matteo [Next page](https://gromacs.bioexcel.eu/t/drude-for-gromacs/1430.md?page=2)