[gmx-users] Dielectric constant for alkanes based on shell polarizable model
aldi asmadi
aldi.asmadi at gmail.com
Tue Feb 22 14:50:37 CET 2011
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
More information about the gromacs.org_gmx-users
mailing list