Charmm-gui lipid parameters appear different than those from pdb2gmx with the gromacs charmm port

GROMACS version: 2022.4
GROMACS modification: No

I am posting what appears to be a difference in parameters out of charmm-gui and gmx pdb2gmx with the gromacs port of the charmm ff in charmm36-jul2022.ff, which I downloaded from MacKerell Lab . I realize that this might not be the right forum to ask for help, but anyway I thought it might help somebody in the future to know that there are difference.

  1. Issue with proper dihedral around double bonds in acyl chains (e.g., POPC):
    This proper dihedral is defined in forcefield.itp from charmm-gui, but not in the standard gromacs port of charmm force field from the MacKerell website:
    CTL2 CEL1 CEL1 HEL1 9 1.800000e+02 2.510400e+01 2

  2. Issue with CHL1 dihedrals:
    This one is a quantitatively minor difference, but there is clearly a difference in charmm-gui and the distributed charm ff
    < phiA= 1.80000000e+02 cpA= 5.09611200e-01 phiB= 1.80000000e+02 cpB= 5.09611200e-01 mult=3


phiA= 1.80000000e+02 cpA= 5.09611000e-01 phiB= 1.80000000e+02 cpB= 5.09611000e-01 mult=3

  1. Issue with CHL1 LJ:
    There are more difference with CHL1 than this, but this is one example.
    Regular Charmm ffnonbonded.itp:
    CRL1 CTL1 1 0.338541512893 0.041840000000
    Charmm-gui forcefield.itp:
    CRL1 CTL1 1 3.57250385974e-01 1.15060000000e-01

It’s possible that I have made some mistakes here, but it’s worth noting that e.g., DPPC comes out with the same energies in both approaches given two separately generated topologies.

Thank you,
Chris.

It’s best to reach out to Alex MacKerell and Wonpil Im directly for any discrepancy of standard CHARMM files vs. CHARMM-GUI. What is in the GROMACS port for C36 comes directly from the CHARMM36 toppar files that Alex provides me. I am not really privy to the backend of CHARMM-GUI or what changes they may have made. That’s relevant to issue 1.

Issue 2 appears to be simply a rounding problem when compiling the .tpr file; the difference occurs at greater precision than our actual parameters.

Can you clarify issue 3? The values in ffnonbonded.itp are for 1-4 interactions but there is also an NBFIX for these two atom types when they are not in a 1-4 pair. Our nbfix.itp value is the same as the CHARMM-GUI parameters. What you list above from ffnonbonded.itp is for the 1-4 interaction.

Great. Thank you Justin. Good catch on #2. I hadn’t considered that. Regarding #3, that does indeed come up as a potential energy difference in the LJ1-4 term. I had figured that was derived from the parameters I posted initially, though possibly I am confused about how that is handled.

This is a test system with a single CHL1 molecule. I made 3 single-point energy evaluations (integrator = md; nsteps = 0; continuation = yes) with gmx 2022.4.

A) forcefield.itp and CHL1.itp from charmm-gui:

Energies (kJ/mol)
Bond U-B Proper Dih. LJ-14 Coulomb-14
4.64237e+01 1.61771e+02 1.01454e+02 3.70265e+01 -1.92795e+02
LJ (SR) Coulomb (SR) Coul. recip. Potential Kinetic En.
1.25963e+01 -6.75563e+01 1.04466e+01 1.09367e+02 1.30141e+02
Total Energy Temperature Pressure (bar)
2.39507e+02 1.42944e+02 1.20008e+00

B) forcefield.itp from charmm-gui and CHL1.itp from the output of pdb2gmx:

Energies (kJ/mol)
Bond U-B Proper Dih. LJ-14 Coulomb-14
4.64237e+01 1.61771e+02 1.01454e+02 3.70265e+01 -1.92795e+02
LJ (SR) Coulomb (SR) Coul. recip. Potential Kinetic En.
1.25963e+01 -6.75563e+01 1.04466e+01 1.09367e+02 1.30141e+02
Total Energy Temperature Pressure (bar)
2.39507e+02 1.42944e+02 1.20008e+00

C) forcefield.itp from charmm36-jul2022.ff and CHL1.itp from the output of pdb2gmx:

Energies (kJ/mol)
Bond U-B Proper Dih. LJ-14 Coulomb-14
4.64237e+01 1.61771e+02 1.01454e+02 3.51514e+01 -1.92795e+02
LJ (SR) Coulomb (SR) Coul. recip. Potential Kinetic En.
1.25963e+01 -6.75563e+01 1.04466e+01 1.07492e+02 1.30106e+02
Total Energy Temperature Pressure (bar)
2.37597e+02 1.42905e+02 1.12914e+00

(hopefully the alignment above works out. I see the difference in the LJ-14 term for A==B but !=C)

Can you provide CHL1.itp that CHARMM-GUI gave you? There is an old version of cholesterol that has been replaced; we use the new one in the port, but that’s been the case for a while.

CHL1.itp from charmm-gui

;;
;; Generated by CHARMM-GUI FF-Converter
;;
;; Correspondance:
;; jul316@lehigh.edu or wonpil@lehigh.edu
;;
;; GROMACS topology file for CHL1
;;

[ moleculetype ]
; name nrexcl
CHL1 3

[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 CRL1 1 CHL1 C3 1 0.140000 12.0110 ; qtot 0.140
2 HGA1 1 CHL1 H3 2 0.090000 1.0080 ; qtot 0.230
3 OHL 1 CHL1 O3 3 -0.660000 15.9994 ; qtot -0.430
4 HOL 1 CHL1 H3’ 4 0.430000 1.0080 ; qtot -0.000
5 CRL2 1 CHL1 C4 5 -0.180000 12.0110 ; qtot -0.180
6 HGA2 1 CHL1 H4A 6 0.090000 1.0080 ; qtot -0.090
7 HGA2 1 CHL1 H4B 7 0.090000 1.0080 ; qtot -0.000
8 CEL1 1 CHL1 C5 8 0.000000 12.0110 ; qtot -0.000
9 CEL1 1 CHL1 C6 9 -0.150000 12.0110 ; qtot -0.150
10 HEL1 1 CHL1 H6 10 0.150000 1.0080 ; qtot -0.000
11 CRL2 1 CHL1 C7 11 -0.180000 12.0110 ; qtot -0.180
12 HGA2 1 CHL1 H7A 12 0.090000 1.0080 ; qtot -0.090
13 HGA2 1 CHL1 H7B 13 0.090000 1.0080 ; qtot -0.000
14 CRL1 1 CHL1 C8 14 -0.090000 12.0110 ; qtot -0.090
15 HGA1 1 CHL1 H8 15 0.090000 1.0080 ; qtot -0.000
16 CRL1 1 CHL1 C14 16 -0.090000 12.0110 ; qtot -0.090
17 HGA1 1 CHL1 H14 17 0.090000 1.0080 ; qtot -0.000
18 CRL2 1 CHL1 C15 18 -0.180000 12.0110 ; qtot -0.180
19 HGA2 1 CHL1 H15A 19 0.090000 1.0080 ; qtot -0.090
20 HGA2 1 CHL1 H15B 20 0.090000 1.0080 ; qtot -0.000
21 CRL2 1 CHL1 C16 21 -0.180000 12.0110 ; qtot -0.180
22 HGA2 1 CHL1 H16A 22 0.090000 1.0080 ; qtot -0.090
23 HGA2 1 CHL1 H16B 23 0.090000 1.0080 ; qtot -0.000
24 CRL1 1 CHL1 C17 24 -0.090000 12.0110 ; qtot -0.090
25 HGA1 1 CHL1 H17 25 0.090000 1.0080 ; qtot -0.000
26 CRL1 1 CHL1 C13 26 0.000000 12.0110 ; qtot -0.000
27 CTL3 1 CHL1 C18 27 -0.270000 12.0110 ; qtot -0.270
28 HAL3 1 CHL1 H18A 28 0.090000 1.0080 ; qtot -0.180
29 HAL3 1 CHL1 H18B 29 0.090000 1.0080 ; qtot -0.090
30 HAL3 1 CHL1 H18C 30 0.090000 1.0080 ; qtot -0.000
31 CRL2 1 CHL1 C12 31 -0.180000 12.0110 ; qtot -0.180
32 HGA2 1 CHL1 H12A 32 0.090000 1.0080 ; qtot -0.090
33 HGA2 1 CHL1 H12B 33 0.090000 1.0080 ; qtot -0.000
34 CRL2 1 CHL1 C11 34 -0.180000 12.0110 ; qtot -0.180
35 HGA2 1 CHL1 H11A 35 0.090000 1.0080 ; qtot -0.090
36 HGA2 1 CHL1 H11B 36 0.090000 1.0080 ; qtot -0.000
37 CRL1 1 CHL1 C9 37 -0.090000 12.0110 ; qtot -0.090
38 HGA1 1 CHL1 H9 38 0.090000 1.0080 ; qtot -0.000
39 CRL1 1 CHL1 C10 39 0.000000 12.0110 ; qtot -0.000
40 CTL3 1 CHL1 C19 40 -0.270000 12.0110 ; qtot -0.270
41 HAL3 1 CHL1 H19A 41 0.090000 1.0080 ; qtot -0.180
42 HAL3 1 CHL1 H19B 42 0.090000 1.0080 ; qtot -0.090
43 HAL3 1 CHL1 H19C 43 0.090000 1.0080 ; qtot -0.000
44 CRL2 1 CHL1 C1 44 -0.180000 12.0110 ; qtot -0.180
45 HGA2 1 CHL1 H1A 45 0.090000 1.0080 ; qtot -0.090
46 HGA2 1 CHL1 H1B 46 0.090000 1.0080 ; qtot -0.000
47 CRL2 1 CHL1 C2 47 -0.180000 12.0110 ; qtot -0.180
48 HGA2 1 CHL1 H2A 48 0.090000 1.0080 ; qtot -0.090
49 HGA2 1 CHL1 H2B 49 0.090000 1.0080 ; qtot -0.000
50 CTL1 1 CHL1 C20 50 -0.090000 12.0110 ; qtot -0.090
51 HAL1 1 CHL1 H20 51 0.090000 1.0080 ; qtot -0.000
52 CTL3 1 CHL1 C21 52 -0.270000 12.0110 ; qtot -0.270
53 HAL3 1 CHL1 H21A 53 0.090000 1.0080 ; qtot -0.180
54 HAL3 1 CHL1 H21B 54 0.090000 1.0080 ; qtot -0.090
55 HAL3 1 CHL1 H21C 55 0.090000 1.0080 ; qtot -0.000
56 CTL2 1 CHL1 C22 56 -0.180000 12.0110 ; qtot -0.180
57 HAL2 1 CHL1 H22A 57 0.090000 1.0080 ; qtot -0.090
58 HAL2 1 CHL1 H22B 58 0.090000 1.0080 ; qtot -0.000
59 CTL2 1 CHL1 C23 59 -0.180000 12.0110 ; qtot -0.180
60 HAL2 1 CHL1 H23A 60 0.090000 1.0080 ; qtot -0.090
61 HAL2 1 CHL1 H23B 61 0.090000 1.0080 ; qtot -0.000
62 CTL2 1 CHL1 C24 62 -0.180000 12.0110 ; qtot -0.180
63 HAL2 1 CHL1 H24A 63 0.090000 1.0080 ; qtot -0.090
64 HAL2 1 CHL1 H24B 64 0.090000 1.0080 ; qtot -0.000
65 CTL1 1 CHL1 C25 65 -0.090000 12.0110 ; qtot -0.090
66 HAL1 1 CHL1 H25 66 0.090000 1.0080 ; qtot -0.000
67 CTL3 1 CHL1 C26 67 -0.270000 12.0110 ; qtot -0.270
68 HAL3 1 CHL1 H26A 68 0.090000 1.0080 ; qtot -0.180
69 HAL3 1 CHL1 H26B 69 0.090000 1.0080 ; qtot -0.090
70 HAL3 1 CHL1 H26C 70 0.090000 1.0080 ; qtot -0.000
71 CTL3 1 CHL1 C27 71 -0.270000 12.0110 ; qtot -0.270
72 HAL3 1 CHL1 H27A 72 0.090000 1.0080 ; qtot -0.180
73 HAL3 1 CHL1 H27B 73 0.090000 1.0080 ; qtot -0.090
74 HAL3 1 CHL1 H27C 74 0.090000 1.0080 ; qtot -0.000

[ bonds ]
; ai aj funct b0 Kb
1 2 1
1 3 1
1 5 1
1 47 1
3 4 1
5 6 1
5 7 1
5 8 1
8 9 1
8 39 1
9 10 1
9 11 1
11 12 1
11 13 1
11 14 1
14 15 1
14 16 1
14 37 1
16 17 1
16 18 1
16 26 1
18 19 1
18 20 1
18 21 1
21 22 1
21 23 1
21 24 1
24 25 1
24 26 1
24 50 1
26 27 1
26 31 1
27 28 1
27 29 1
27 30 1
31 32 1
31 33 1
31 34 1
34 35 1
34 36 1
34 37 1
37 38 1
37 39 1
39 40 1
39 44 1
40 41 1
40 42 1
40 43 1
44 45 1
44 46 1
47 44 1
47 48 1
47 49 1
50 51 1
50 52 1
50 56 1
52 53 1
52 54 1
52 55 1
56 57 1
56 58 1
56 59 1
59 60 1
59 61 1
59 62 1
62 63 1
62 64 1
62 65 1
65 66 1
65 67 1
65 71 1
67 68 1
67 69 1
67 70 1
71 72 1
71 73 1
71 74 1

[ pairs ]
; ai aj funct c6 c12 or
; ai aj funct fudgeQQ q1 q2 c6 c12
1 9 1
1 39 1
1 45 1
1 46 1
2 4 1
2 6 1
2 7 1
2 8 1
2 44 1
2 48 1
2 49 1
3 6 1
3 7 1
3 8 1
3 44 1
3 48 1
3 49 1
4 5 1
4 47 1
5 10 1
5 11 1
5 37 1
5 40 1
5 44 1
5 48 1
5 49 1
6 9 1
6 39 1
6 47 1
7 9 1
7 39 1
7 47 1
8 12 1
8 13 1
8 14 1
8 34 1
8 38 1
8 41 1
8 42 1
8 43 1
8 45 1
8 46 1
8 47 1
9 15 1
9 16 1
9 37 1
9 40 1
9 44 1
10 12 1
10 13 1
10 14 1
10 39 1
11 17 1
11 18 1
11 26 1
11 34 1
11 38 1
11 39 1
12 15 1
12 16 1
12 37 1
13 15 1
13 16 1
13 37 1
14 19 1
14 20 1
14 21 1
14 24 1
14 27 1
14 31 1
14 35 1
14 36 1
14 40 1
14 44 1
15 17 1
15 18 1
15 26 1
15 34 1
15 38 1
15 39 1
16 22 1
16 23 1
16 25 1
16 28 1
16 29 1
16 30 1
16 32 1
16 33 1
16 34 1
16 38 1
16 39 1
16 50 1
17 19 1
17 20 1
17 21 1
17 24 1
17 27 1
17 31 1
17 37 1
18 25 1
18 27 1
18 31 1
18 37 1
18 50 1
19 22 1
19 23 1
19 24 1
19 26 1
20 22 1
20 23 1
20 24 1
20 26 1
21 27 1
21 31 1
21 51 1
21 52 1
21 56 1
22 25 1
22 26 1
22 50 1
23 25 1
23 26 1
23 50 1
24 28 1
24 29 1
24 30 1
24 32 1
24 33 1
24 34 1
24 53 1
24 54 1
24 55 1
24 57 1
24 58 1
24 59 1
25 27 1
25 31 1
25 51 1
25 52 1
25 56 1
26 35 1
26 36 1
26 37 1
26 51 1
26 52 1
26 56 1
27 32 1
27 33 1
27 34 1
27 50 1
28 31 1
29 31 1
30 31 1
31 38 1
31 39 1
31 50 1
32 35 1
32 36 1
32 37 1
33 35 1
33 36 1
33 37 1
34 40 1
34 44 1
35 38 1
35 39 1
36 38 1
36 39 1
37 41 1
37 42 1
37 43 1
37 45 1
37 46 1
37 47 1
38 40 1
38 44 1
39 48 1
39 49 1
40 45 1
40 46 1
40 47 1
41 44 1
42 44 1
43 44 1
45 48 1
45 49 1
46 48 1
46 49 1
50 60 1
50 61 1
50 62 1
51 53 1
51 54 1
51 55 1
51 57 1
51 58 1
51 59 1
52 57 1
52 58 1
52 59 1
53 56 1
54 56 1
55 56 1
56 63 1
56 64 1
56 65 1
57 60 1
57 61 1
57 62 1
58 60 1
58 61 1
58 62 1
59 66 1
59 67 1
59 71 1
60 63 1
60 64 1
60 65 1
61 63 1
61 64 1
61 65 1
62 68 1
62 69 1
62 70 1
62 72 1
62 73 1
62 74 1
63 66 1
63 67 1
63 71 1
64 66 1
64 67 1
64 71 1
66 68 1
66 69 1
66 70 1
66 72 1
66 73 1
66 74 1
67 72 1
67 73 1
67 74 1
68 71 1
69 71 1
70 71 1

[ angles ]
; ai aj ak funct th0 cth S0 Kub
2 1 3 5
2 1 5 5
2 1 47 5
3 1 5 5
3 1 47 5
5 1 47 5
1 3 4 5
1 5 6 5
1 5 7 5
1 5 8 5
6 5 7 5
6 5 8 5
7 5 8 5
5 8 9 5
5 8 39 5
9 8 39 5
8 9 10 5
8 9 11 5
10 9 11 5
9 11 12 5
9 11 13 5
9 11 14 5
12 11 13 5
12 11 14 5
13 11 14 5
11 14 15 5
11 14 16 5
11 14 37 5
15 14 16 5
15 14 37 5
16 14 37 5
14 16 17 5
14 16 18 5
14 16 26 5
17 16 18 5
17 16 26 5
18 16 26 5
16 18 19 5
16 18 20 5
16 18 21 5
19 18 20 5
19 18 21 5
20 18 21 5
18 21 22 5
18 21 23 5
18 21 24 5
22 21 23 5
22 21 24 5
23 21 24 5
21 24 25 5
21 24 26 5
21 24 50 5
25 24 26 5
25 24 50 5
26 24 50 5
16 26 24 5
16 26 27 5
16 26 31 5
24 26 27 5
24 26 31 5
27 26 31 5
26 27 28 5
26 27 29 5
26 27 30 5
28 27 29 5
28 27 30 5
29 27 30 5
26 31 32 5
26 31 33 5
26 31 34 5
32 31 33 5
32 31 34 5
33 31 34 5
31 34 35 5
31 34 36 5
31 34 37 5
35 34 36 5
35 34 37 5
36 34 37 5
14 37 34 5
14 37 38 5
14 37 39 5
34 37 38 5
34 37 39 5
38 37 39 5
8 39 37 5
8 39 40 5
8 39 44 5
37 39 40 5
37 39 44 5
40 39 44 5
39 40 41 5
39 40 42 5
39 40 43 5
41 40 42 5
41 40 43 5
42 40 43 5
39 44 45 5
39 44 46 5
39 44 47 5
45 44 46 5
45 44 47 5
46 44 47 5
1 47 44 5
1 47 48 5
1 47 49 5
44 47 48 5
44 47 49 5
48 47 49 5
24 50 51 5
24 50 52 5
24 50 56 5
51 50 52 5
51 50 56 5
52 50 56 5
50 52 53 5
50 52 54 5
50 52 55 5
53 52 54 5
53 52 55 5
54 52 55 5
50 56 57 5
50 56 58 5
50 56 59 5
57 56 58 5
57 56 59 5
58 56 59 5
56 59 60 5
56 59 61 5
56 59 62 5
60 59 61 5
60 59 62 5
61 59 62 5
59 62 63 5
59 62 64 5
59 62 65 5
63 62 64 5
63 62 65 5
64 62 65 5
62 65 66 5
62 65 67 5
62 65 71 5
66 65 67 5
66 65 71 5
67 65 71 5
65 67 68 5
65 67 69 5
65 67 70 5
68 67 69 5
68 67 70 5
69 67 70 5
65 71 72 5
65 71 73 5
65 71 74 5
72 71 73 5
72 71 74 5
73 71 74 5

[ dihedrals ]
; ai aj ak al funct phi0 cp mult
2 1 3 4 9
2 1 5 6 9
2 1 5 7 9
2 1 5 8 9
3 1 5 6 9
3 1 5 7 9
3 1 5 8 9
2 1 47 44 9
2 1 47 48 9
2 1 47 49 9
3 1 47 44 9
3 1 47 48 9
3 1 47 49 9
5 1 47 44 9
5 1 47 48 9
5 1 47 49 9
4 3 1 5 9
4 3 1 47 9
6 5 1 47 9
7 5 1 47 9
8 5 1 47 9
1 5 8 9 9
1 5 8 39 9
6 5 8 9 9
6 5 8 39 9
7 5 8 9 9
7 5 8 39 9
5 8 9 10 9
5 8 9 11 9
5 8 39 37 9
5 8 39 40 9
5 8 39 44 9
9 8 39 37 9
9 8 39 40 9
9 8 39 44 9
10 9 8 39 9
11 9 8 39 9
8 9 11 12 9
8 9 11 13 9
8 9 11 14 9
10 9 11 12 9
10 9 11 13 9
10 9 11 14 9
9 11 14 15 9
9 11 14 16 9
9 11 14 37 9
12 11 14 15 9
12 11 14 16 9
12 11 14 37 9
13 11 14 15 9
13 11 14 16 9
13 11 14 37 9
11 14 16 17 9
11 14 16 18 9
11 14 16 26 9
15 14 16 17 9
15 14 16 18 9
15 14 16 26 9
11 14 37 34 9
11 14 37 38 9
11 14 37 39 9
15 14 37 34 9
15 14 37 38 9
15 14 37 39 9
16 14 37 34 9
16 14 37 38 9
16 14 37 39 9
17 16 14 37 9
18 16 14 37 9
26 16 14 37 9
14 16 18 19 9
14 16 18 20 9
14 16 18 21 9
17 16 18 19 9
17 16 18 20 9
17 16 18 21 9
14 16 26 24 9
14 16 26 27 9
14 16 26 31 9
17 16 26 24 9
17 16 26 27 9
17 16 26 31 9
18 16 26 24 9
18 16 26 27 9
18 16 26 31 9
19 18 16 26 9
20 18 16 26 9
21 18 16 26 9
16 18 21 22 9
16 18 21 23 9
16 18 21 24 9
19 18 21 22 9
19 18 21 23 9
19 18 21 24 9
20 18 21 22 9
20 18 21 23 9
20 18 21 24 9
18 21 24 25 9
18 21 24 26 9
18 21 24 50 9
22 21 24 25 9
22 21 24 26 9
22 21 24 50 9
23 21 24 25 9
23 21 24 26 9
23 21 24 50 9
21 24 26 27 9
21 24 26 31 9
25 24 26 27 9
25 24 26 31 9
21 24 50 51 9
21 24 50 52 9
21 24 50 56 9
25 24 50 51 9
25 24 50 52 9
25 24 50 56 9
26 24 50 51 9
26 24 50 52 9
26 24 50 56 9
16 26 24 21 9
16 26 24 25 9
16 26 24 50 9
27 26 24 50 9
31 26 24 50 9
16 26 27 28 9
16 26 27 29 9
16 26 27 30 9
24 26 27 28 9
24 26 27 29 9
24 26 27 30 9
16 26 31 32 9
16 26 31 33 9
16 26 31 34 9
24 26 31 32 9
24 26 31 33 9
24 26 31 34 9
27 26 31 32 9
27 26 31 33 9
27 26 31 34 9
28 27 26 31 9
29 27 26 31 9
30 27 26 31 9
26 31 34 35 9
26 31 34 36 9
26 31 34 37 9
32 31 34 35 9
32 31 34 36 9
32 31 34 37 9
33 31 34 35 9
33 31 34 36 9
33 31 34 37 9
31 34 37 38 9
31 34 37 39 9
35 34 37 38 9
35 34 37 39 9
36 34 37 38 9
36 34 37 39 9
14 37 34 31 9
14 37 34 35 9
14 37 34 36 9
14 37 39 40 9
14 37 39 44 9
34 37 39 40 9
34 37 39 44 9
38 37 39 40 9
38 37 39 44 9
8 39 37 14 9
8 39 37 34 9
8 39 37 38 9
8 39 40 41 9
8 39 40 42 9
8 39 40 43 9
37 39 40 41 9
37 39 40 42 9
37 39 40 43 9
8 39 44 45 9
8 39 44 46 9
8 39 44 47 9
37 39 44 45 9
37 39 44 46 9
37 39 44 47 9
40 39 44 45 9
40 39 44 46 9
40 39 44 47 9
41 40 39 44 9
42 40 39 44 9
43 40 39 44 9
39 44 47 48 9
39 44 47 49 9
45 44 47 48 9
45 44 47 49 9
46 44 47 48 9
46 44 47 49 9
1 47 44 39 9
1 47 44 45 9
1 47 44 46 9
24 50 52 53 9
24 50 52 54 9
24 50 52 55 9
51 50 52 53 9
51 50 52 54 9
51 50 52 55 9
24 50 56 57 9
24 50 56 58 9
24 50 56 59 9
51 50 56 57 9
51 50 56 58 9
51 50 56 59 9
52 50 56 57 9
52 50 56 58 9
52 50 56 59 9
53 52 50 56 9
54 52 50 56 9
55 52 50 56 9
50 56 59 60 9
50 56 59 61 9
50 56 59 62 9
57 56 59 60 9
57 56 59 61 9
57 56 59 62 9
58 56 59 60 9
58 56 59 61 9
58 56 59 62 9
56 59 62 63 9
56 59 62 64 9
56 59 62 65 9
60 59 62 63 9
60 59 62 64 9
60 59 62 65 9
61 59 62 63 9
61 59 62 64 9
61 59 62 65 9
59 62 65 66 9
59 62 65 67 9
59 62 65 71 9
63 62 65 66 9
63 62 65 67 9
63 62 65 71 9
64 62 65 66 9
64 62 65 67 9
64 62 65 71 9
62 65 67 68 9
62 65 67 69 9
62 65 67 70 9
66 65 67 68 9
66 65 67 69 9
66 65 67 70 9
62 65 71 72 9
62 65 71 73 9
62 65 71 74 9
66 65 71 72 9
66 65 71 73 9
66 65 71 74 9
67 65 71 72 9
67 65 71 73 9
67 65 71 74 9
68 67 65 71 9
69 67 65 71 9
70 67 65 71 9

#ifdef POSRES
[ position_restraints ]
3 1 0.0 0.0 POSRES_FC_LIPID
#endif

#ifdef DIHRES
[ dihedral_restraints ]
#endif

And here’s the final possible combination. Seems that the difference tracks with the forcefield.itp not with the CHL1.itp file:

D) forcefield.itp from charmm36-jul2022.ff and CHL1.itp from charmm-gui:

Energies (kJ/mol)
Bond U-B Proper Dih. LJ-14 Coulomb-14
4.64237e+01 1.61771e+02 1.01454e+02 3.51514e+01 -1.92795e+02
LJ (SR) Coulomb (SR) Coul. recip. Potential Kinetic En.
1.25963e+01 -6.75563e+01 1.04466e+01 1.07492e+02 1.30106e+02
Total Energy Temperature Pressure (bar)
2.37597e+02 1.42905e+02 1.12914e+00

PS: thanks a lot for your help tracking this down Justin!

Cross posted to the charmm-gui forum here: CHARMM-GUI

OK, I’m doing the comparison that I know how to do and what I use for checking the port, which is CHARMM vs. GROMACS energies. Everything looks good except for the LJ1-4 term, which is really weird because all the parameters involved look to be correctly implemented.

Thanks Justin! Do I understand correctly that you also checked POPC with your approach and didn’t find a difference with the charmm and gromacs energies in the proper dihedral term?

No problem with POPC. The energies are virtually identical; the single point energy is 11.722 kcal/mol in CHARMM and 11.723 kcal/mol in GROMACS.

Thanks Justin. I suppose for POPC that points to charmm-gui having an unexpected CTL2-CEL1-CEL1-HEL1 proper dihedral term. Sort of scary, since so many people use that to build lipid systems now.

It’s not unexpected; it’s there. I forgot about wildcards.

CHARMM:

X    CEL1 CEL1 X        0.4500  1   180.00 ! 2-butene, adm jr., 4/04
X    CEL1 CEL1 X        8.5000  2   180.00 !

GROMACS:

       X     CEL1     CEL1        X     9   180.000000     1.882800     1 ; 2-butene, adm jr., 4/04
       X     CEL1     CEL1        X     9   180.000000    35.564000     2 ;

CHARMM-GUI is just providing it as an explicit parameter rather than the wildcard version.

Quick summary: For POPC, I think the presence of a single (extra) explicit parameter for CTL2-CEL1-CEL1-HEL1 dihedrals in the charmm-gui forcefield.itp:
CTL2 CEL1 CEL1 HEL1 9 1.800000e+02 2.510400e+01 2

is clobbering and replacing the 2x wildcard parameters from charmm36-jul2022.ff/ffbonded.itp (that are in fact also in the charmm-gui forcefield.itp but not used due to the explicit parameter)
X CEL1 CEL1 X 9 1.800000e+02 1.882800e+00 1
X CEL1 CEL1 X 9 1.800000e+02 3.556400e+01 2

Indeed, commenting out that CTL2-CEL1-CEL1-HEL1 dihedral parameter in the forcefield.itp file from charmm-gui leads to the same energies as charmm36-jul2022.ff/ffbonded.itp

More details:

These are the relevant atoms in both charmm-gui and charmm36-jul2022.ff versions of POPC
atom (134):
atom[59]={name=“C28”}
atom[62]={name=“C29”}
atom[64]={name=“C210”}
atom[65]={name=“H101”}
atom[66]={name=“C211”}

#full list of PDIHS with 62/64 (or 64/62) in middle from charmm-gui:
192 type=318 (PDIHS) 59 62 64 65
193 type=319 (PDIHS) 59 62 64 66
194 type=320 (PDIHS) 59 62 64 66
195 type=321 (PDIHS) 63 62 64 65
196 type=318 (PDIHS) 63 62 64 66

#relevant charmm-gui dihedral definitions:
functype[318]=PDIHS, phiA= 1.80000000e+02, cpA= 2.51040000e+01, phiB= 1.80000000e+02, cpB= 2.51040000e+01, mult=2
functype[319]=PDIHS, phiA= 1.80000000e+02, cpA= 1.88280000e+00, phiB= 1.80000000e+02, cpB= 1.88280000e+00, mult=1
functype[320]=PDIHS, phiA= 1.80000000e+02, cpA= 3.55640000e+01, phiB= 1.80000000e+02, cpB= 3.55640000e+01, mult=2
functype[321]=PDIHS, phiA= 1.80000000e+02, cpA= 4.18400000e+00, phiB= 1.80000000e+02, cpB= 4.18400000e+00, mult=2

#full list with 62/64 (or 64/62) in middle from charmm36-jul2022.ff:
192 type=349 (PDIHS) 59 62 64 65
193 type=350 (PDIHS) 59 62 64 65
194 type=349 (PDIHS) 59 62 64 66
195 type=350 (PDIHS) 59 62 64 66
196 type=351 (PDIHS) 63 62 64 65
197 type=349 (PDIHS) 63 62 64 66
198 type=350 (PDIHS) 63 62 64 66

#relevant charmm36-jul2022.ff dihedral definitions:
functype[349]=PDIHS, phiA= 1.80000000e+02, cpA= 1.88280000e+00, phiB= 1.80000000e+02, cpB= 1.88280000e+00, mult=1
functype[350]=PDIHS, phiA= 1.80000000e+02, cpA= 3.55640000e+01, phiB= 1.80000000e+02, cpB= 3.55640000e+01, mult=2
functype[351]=PDIHS, phiA= 1.80000000e+02, cpA= 4.18400000e+00, phiB= 1.80000000e+02, cpB= 4.18400000e+00, mult=2

my conclusions from the above data:

  1. The CTL2-CEL1-CEL1-CTL2 dihedral (C28-C29-C210-C211 = 59-62-64-66) is the same for both cases (2x terms)
  2. The HEL1-CEL1-CEL1-CEL1 dihedral (H91-C29-C210-H101 = 63-62-64-65) is the same for both cases (1x term)
  3. Charmm-gui is adding phiA= 1.80000000e+02, cpA= 2.51040000e+01, mult=2 to the CTL2-CEL1-CEL1-HEL1 dihedrals (C28-C29-C210-H101 = 59-62-64-65 and H91-C29-C210-C211 = 63-62-64-66) but charmm36-jul2022.ff is not (this is the extra dihedral potential I was referring to earlier)
  4. charmm36-jul2022.ff is adding phiA= 1.80000000e+02, cpA= 1.88280000e+00, mult=1 and phiA= 1.80000000e+02, cpA= 3.55640000e+01, mult=2 to the CTL2-CEL1-CEL1-HEL1 dihedrals (C28-C29-C210-H101 = 59-62-64-65 and H91-C29-C210-C211 = 63-62-64-66) but charmm-gui is not (I did not appreciate this part earlier)

The forcefield.itp from charm-guy does indeed have the 2x X-CEL1-CEL1-X dihedral parameters that match charmm36-jul2022.ff:
X CEL1 CEL1 X 9 1.800000e+02 1.882800e+00 1
X CEL1 CEL1 X 9 1.800000e+02 3.556400e+01 2

However, it also has an explicit parameter that differs from charmm36-jul2022.ff:
CTL2 CEL1 CEL1 HEL1 9 1.800000e+02 2.510400e+01 2

And another that is the same as charmm36-jul2022.ff:
HEL1 CEL1 CEL1 HEL1 9 1.800000e+02 4.184000e+00 2

Here’s those matching wildcard parameters from charmm36-jul2022.ff/ffbonded.itp:
X CEL1 CEL1 X 9 180.000000 1.882800 1 ; 2-butene, adm jr., 4/04
X CEL1 CEL1 X 9 180.000000 35.564000 2 ;

And the matching specific parameter in charmm36-jul2022.ff/ffbonded.itp that behaves the same as charmm-gui:
HEL1 CEL1 CEL1 HEL1 9 180.000000 4.184000 2 ; 2-butene, adm jr., 8/98 update

But charmm36-jul2022.ff/ffbonded.itp doesn’t have a specific parameter for CTL2-CEL1-CEL1-HEL1 dihedrals.

In my distant memory, I think that wildcard parameters may only be used when specific parameters don’t exist (or perhaps ordering is relevant for that general rule). This would explain why the forcefield.itp from charmm-gui is yielding different parameters for the CTL2-CEL1-CEL1-HEL1 dihedral and gives different energies for POPC in the proper dihedral term.

Dear Justin:

Jumin Lee is looking into the POPC and DPPEE dihedral parameters (see the charmm-gui forum link in this thread).

Regarding CHL1, this is what he said on the charmm-gui forum (pasted below). Perhaps it will help track down the difference between charmm and gromacs energies for CHL1.

For “CRL1 CTL1”, what CHARMM-GUI provides is correct.
If the NBFIX parameter present, it should be applied to 1-4 non-bonded interaction as well.
The forcefield.itp from CHARMM-GUI does not contain “CRL1 CTL1” pairtype in [ pairtypes ], so the pairtype for “CRL1 CTL1” is automatically generated using the parameter in [ nonbond_params ].

Thank you,
Chris.

This is really interesting, and something I did not even know about CHARMM. We will need to fix the port, as I have confirmed that removing the NBFIX pairs from [pairtypes] leads to the correct energies. This should be easy to do in the conversion process, but I will confirm with @awacha on that.

Great news. Thank you Justin.

Hi Chris,

thanks for finding this about the NBFIX! If I understand you correctly, NBFIX takes precedence even in 1-4 LJ interactions, so when [pairtypes] are generated, \epsilon and \sigma values should be taken from there, instead of combining them from those given at a per-atomtype basis in the CHARMM prm files. I’ve updated the conversion program, so it’s now Justin’s turn to verify if the parameters are correct.

On the dihedral parameters I think that the Gromacs port is consistent with the CHARMM toppar files. I can confirm that no CTL2-CEL1-CEL1-HEL1 dihedral parameters (nor HEL1-CEL1-CEL1-CTL2 for that matter) exists in the most recent version (2022), this set of atom types is only matched by the X-CEL1-CEL1-X wildcard entry.

Thanks again,

Andras

Thank you Andras for the quick fix!

I have run through my battery of tests and the latest port from Andras’ software is correct for these odd cases involving NBFIXes and special 1-4 LJ. Cholesterol energies now match between CHARMM and GROMACS at the same level of precision as everything else we test for. The new port with the corrected files is available on the MacKerell Lab webpage.