Dear Gromacs users, I really hope someone could give me some advise on my following problem. I am simulating liquid bulk hexane using a shell (or Drude oscillator) model that incorporates the electronic polarization in GROMACS 4.5.3. To do this, I use the force field parameters from the CHARMM drude-based polarizable force field developed by Vorobyov, Anisomov, and MacKerell (http://pubs.acs.org/doi/abs/10.1021/jp053182y).
For my simulation, I have a cubic box with a length of 3.2 nm consisting of 150 molecules, and I employ the following essential MD parameters: coulombtype = PME, rlist = rcoulomb = 1.5, pme_order = 6, fourierspacing = 0.1, vdw-type = switch, rvdw = 1.2, rvdw_switch = 1.0, dispcorr = no. The shell particles are optimized in every time step with settings: emtol = 0.1, niter = 100. The system is calculated for 5 ns. For the bulk properties at 298 K and 1 atm, I have calculated the density (654.7 g/cm3), the self-diffusion constant (3.2 x 10-5 cm2/s), surface tension (vapor-liquid hexane, 16.45 dyn/cm), and they seem to correspond quite fair with the experimental data (655 g/cm3, 4.0 x 10-5 cm2/s, and 17.9 dyn/cm, respectively). However, for the static dielectric constant (eps_0), I only obtained a value of around 1.18 (using both g_dipoles -f eq_npt.xtc -s eq_npt.tpr and g_current -f eq_npt.xtc -s eq_npt.tpr) In the paper, the authors didn't report the dielectric constant for hexane. However, they calculated this property for isobutane (1.81) and heptane (1.91). So, accordingly, the static dielectric constant for hexane should be located between those values. The experimentally extrapolated value for hexane at 298.15K is around 1.88. At the moment, I could not figure out why I couldn't get a positive result for eps_0, and also do not know what went wrong with my system. I am truly stuck. In the paper, the author did mention that for alkanes the high-frequency dependent dielectric constant (eps_inf) is more dominant. So, in order to calculate eps_0, they first calculated eps_inf, and aftewards used eps_inf as the input to calculate eps_0. (In GROMACS, if I understood correctly, g_dielectric is used to approximate eps_inf. But still, we need eps_0 as the input). The authors calculated eps_inf based on the dipole fluctuations of the simulation box associated with the movement of the shell particles, resulted from Langevin simulation on frozen nuclear configuration. However, they used extended Lagrangian formalism in their simulation, and this is not applicable in GROMACS. I don't know if I am missing something obvious. Again, I really hope if anyone could give me suggestion. I have also copied and pasted the top and itp files I used during my simulation below. I honestly apologize since they are quite long. Many thanks in advance. Regards, Aldi ITP File ; Polarizable model of heptane molecules ; Constructed by on a paper by Vorobroy, Anisimov, and MacKerell ; J. Phys. Chem. B, 2005, 109, 18988-18999 ; http://pubs.acs.org/doi/abs/10.1021/jp053182y [ defaults ] ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 2 yes 1 1 [ atomtypes ] ; m(C*) = 12.011 ; m(H*) = 1.008 ; m(D*) = 0.0 ; ; ; LJ parameters for all atom types were taken from the paper ; ; Optimised LJ Parameters ; type Rmin/2 eps ; CT2 2.010 0.056 ; CT3 2.040 0.078 ; HA2 1.340 0.035 ; HA3 1.340 0.024 ; ; eps(gmx) = eps(paper)*4.18400 ; sigma(gmx) = ((Rmin/2)*2)/(2^(1/6)) ; ; q(C_P) = q(C)-q(D*) ; q(C3) = -0.177, q(C2) = -0.156 ; q(D_C3) = -1.111429842, q(D_C2) = -0.9998924 ; q(C3_P) = q(C3) - q(D_C3) = -0.177 + 1.111429842 = 0.934429842 ; q(C2_P) = q(C2) - q(D_C2) = -0.156 + 0.9998924 = 0.8438924 ; name mass charge ptype sigma [nm] eps [kJ/mol] C3_P 12.011 0.934429842 A 0.363486677001 0.326352 C2_P 12.011 0.8438924 A 0.358141284692 0.234304 D_C3 0.0 -1.111429842 S 0.0000000000 0.000000 D_C2 0.0 -0.9998924 S 0.0000000000 0.000000 HA2 1.00800 0.078 A 0.238760856462 0.146440 HA3 1.00800 0.059 A 0.238760856462 0.100416 ; All internal parameters were taken from CHARMM27 force field implemented in GROMACS 4.5.3 [ nonbond_params ] ; sigma epsilon ; TARGET VALUES C3_P C3_P 1 0.3634867 0.3263520 C2_P C2_P 1 0.3581413 0.2343040 HA2 HA2 1 0.2387609 0.1464400 HA3 HA3 1 0.2387609 0.1004160 C3_P C2_P 1 0.3608140 0.2765241 C3_P HA2 1 0.3011238 0.2186115 C3_P HA3 1 0.3011238 0.1810275 C2_P HA2 1 0.2984511 0.1852336 C2_P HA3 1 0.2984511 0.1533880 HA2 HA3 1 0.2387609 0.1212638 [ bondtypes ] ; i j func b0 kb C3_P C2_P 1 0.1528 186188.0 HA3 C3_P 1 0.1111 269449.6 HA2 C2_P 1 0.1111 258571.2 C2_P C2_P 1 0.153 186188.0 [ angletypes ] ; i j k func th0 cth ub0 cub C3_P C2_P C2_P 5 115.00 485.344 0.2561 6694.4 HA3 C3_P C2_P 5 110.10 289.5328 0.2179 18853.104 HA3 C3_P HA3 5 108.40 297.064 0.1802 4518.72 HA2 C2_P C2_P 5 110.10 221.752 0.2179 18853.104 HA2 C2_P HA2 5 109.00 297.064 0.1802 4518.72 C3_P C2_P HA2 5 110.10 289.5328 0.2179 18853.104 C2_P C2_P C2_P 5 113.60 488.2728 0.2561 9338.688 [ dihedraltypes ] ; i j k l func phi0 cp mult ;CT2 CT2 CT2 CT2 9 0.00 0.6276 1 ;CT3 CT2 CT2 CT2 9 0.00 0.6276 1 C3_P C2_P C2_P C2_P 9 0 0.594128 5 C3_P C2_P C2_P C2_P 9 0 0.439320 4 C3_P C2_P C2_P C2_P 9 180 0.343088 3 C3_P C2_P C2_P C2_P 9 0 0.221752 2 C2_P C2_P C2_P C2_P 9 0 0.569024 5 C2_P C2_P C2_P C2_P 9 0 0.326352 4 C2_P C2_P C2_P C2_P 9 180 0.581576 3 C2_P C2_P C2_P C2_P 9 0 0.347272 2 HA2 C2_P C2_P HA2 9 0.00 0.81588 3 C3_P C2_P C2_P HA2 9 0.00 0.81588 3 HA2 C2_P C2_P C2_P 9 0.00 0.81588 3 HA3 C3_P C2_P HA2 9 0.00 0.66944 3 HA3 C3_P C2_P C2_P 9 0.00 0.66944 3 Top file ; Topology file for hexane #include "alkane_mod.itp" [ moleculetype ] ;Name nrexcl Hexane 2 [ atoms ] ;nr type rsnr resd atom cgnr charge mass 1 C3_P 1 HEX C1 1 0.934429842 12.011 2 D_C3 1 HEX D1 1 -1.111429842 0.0000 3 HA3 1 HEX H11 1 0.059 1.00800 4 HA3 1 HEX H12 1 0.059 1.00800 5 HA3 1 HEX H13 1 0.059 1.00800 6 C2_P 1 HEX C2 2 0.8438924 12.011 7 D_C2 1 HEX D2 2 -0.9998924 0.0000 8 HA2 1 HEX H21 2 0.078 1.00800 9 HA2 1 HEX H22 2 0.078 1.00800 10 C2_P 1 HEX C3 3 0.8438924 12.011 11 D_C2 1 HEX D3 3 -0.9998924 0.0000 12 HA2 1 HEX H31 3 0.078 1.00800 13 HA2 1 HEX H32 3 0.078 1.00800 14 C2_P 1 HEX C4 4 0.8438924 12.011 15 D_C2 1 HEX D4 4 -0.9998924 0.0000 16 HA2 1 HEX H41 4 0.078 1.00800 17 HA2 1 HEX H42 4 0.078 1.00800 18 C2_P 1 HEX C5 5 0.8438924 12.011 19 D_C2 1 HEX D5 5 -0.9998924 0.0000 20 HA2 1 HEX H51 5 0.078 1.00800 21 HA2 1 HEX H52 5 0.078 1.00800 22 C3_P 1 HEX C6 6 0.934429842 12.011 23 D_C3 1 HEX D6 6 -1.111429842 0.0000 24 HA3 1 HEX H61 6 0.059 1.00800 25 HA3 1 HEX H62 6 0.059 1.00800 26 HA3 1 HEX H63 6 0.059 1.00800 [ polarization ] ;i j type alpha(nm^3) 1 2 1 2.051e-03 6 7 1 1.660e-03 10 11 1 1.660e-03 14 15 1 1.660e-03 18 19 1 1.660e-03 22 23 1 2.051e-03 [ thole_polarization ] ;ai aj func a alpha(ai) alpha(aj) 1 2 6 7 2 2.6 2.051e-03 1.660e-03 1 2 10 11 2 2.6 2.051e-03 1.660e-03 6 7 10 11 2 2.6 1.660e-03 1.660e-03 6 7 14 15 2 2.6 1.660e-03 1.660e-03 10 11 14 15 2 2.6 1.660e-03 1.660e-03 10 11 18 19 2 2.6 1.660e-03 1.660e-03 14 15 18 19 2 2.6 1.660e-03 1.660e-03 14 15 22 23 2 2.6 1.660e-03 2.051e-03 18 19 22 23 2 2.6 1.660e-03 2.051e-03 [ bonds ] ;i j funct 1 3 1 1 4 1 1 5 1 1 6 1 6 8 1 6 9 1 6 10 1 10 12 1 10 13 1 10 14 1 14 16 1 14 17 1 14 18 1 18 20 1 18 21 1 18 22 1 22 24 1 22 25 1 22 26 1 ; drude-portion 1 2 5 6 7 5 10 11 5 14 15 5 18 19 5 22 23 5 [ angles ] ;i j k func 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 [ dihedrals ] ; i j k l funct 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 [ pairs ] ;i j funct 1 12 1 1 13 1 1 14 1 1 15 1 ; d-a (drude-atom pair interaction) 2 12 1 ; d-a 2 13 1 ; d-a 2 14 1 ; d-a 2 15 1 ; d-d (drude-drude pair interaction) ; 3 8 1 3 9 1 3 10 1 3 11 1 ; d-a 4 8 1 4 9 1 4 10 1 4 11 1 ; d-a 5 8 1 5 9 1 5 10 1 5 11 1 ; d-a 6 16 1 6 17 1 6 18 1 6 19 1 ; d-a 7 16 1 ; d-a 7 17 1 ; d-a 7 18 1 ; d-a 7 19 1 ; d-d (drude-drude pair interaction) ; 8 12 1 8 13 1 8 14 1 8 15 1 ; d-a 9 12 1 9 13 1 9 14 1 9 15 1 ; d-a 10 20 1 10 21 1 10 22 1 10 23 1 ; d-a 11 20 1 ; d-a 11 21 1 ; d-a 11 22 1 ; d-a 11 23 1 ; d-d (drude-drude interaction) ; 12 16 1 12 17 1 12 18 1 12 19 1 ; d-a 13 16 1 13 17 1 13 18 1 13 19 1 ; d-a 14 24 1 14 25 1 14 26 1 15 24 1 ; d-a 15 25 1 ; d-a 15 26 1 ; d-a 16 20 1 16 21 1 16 22 1 16 23 1 ; d-a 17 20 1 17 21 1 17 22 1 17 23 1 ; d-a 20 24 1 20 25 1 20 26 1 21 24 1 21 25 1 21 26 1 [ exclusions ] 1 3 4 5 6 7 8 9 10 11 2 3 4 5 6 7 8 9 10 11 ; drude 3 4 5 6 7 4 5 6 7 5 6 7 6 8 9 10 11 12 13 14 15 7 8 9 10 11 12 13 14 15 ; drude 8 9 10 11 9 10 11 10 12 13 14 15 16 17 18 19 11 12 13 14 15 16 17 18 19 ; drude 12 13 14 15 13 14 15 14 16 17 18 19 20 21 22 23 15 16 17 18 19 20 21 22 23 ; drude 16 17 18 19 17 18 19 18 20 21 22 23 24 25 26 19 20 21 22 23 24 25 26 20 21 22 23 21 22 23 22 24 25 26 23 24 25 26 [ system ] Drude-based hexane molecule [ molecules ] ;Name nr Hexane 150 -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
