[gmx-users] problem with simulation of phosphatidic acid, possibly in my topology?

[himanshu khandelia]

Hi,

I am trying to run a simulation of a POPC bilayer mixed with some
mono-anionic Phosphatidic acid (POPA) , where the choline group is
replaced by a hydrogen atom. However, the energy of my simulation box
diverges, and I am trying to fix the problem. I also tried running a
simulation of PA solvated in water and a single sodium ion, but even
that does not work. I will try to provide as much detail as possible,
I hope someone can help point out some obvious error which I have not
been able to debug over 3 days. In the simulation of a single PA
molecule in water, MDRUN keeps complaining that the distance between
atoms H1 and O4 is very high (> 4 nm or so). The problem probably lies
in my topology, possibly in the bonded interactions?

Thank you very very much for all suggestions.

The headgroup of phosphatidic acid in my topology looks something like this:

           H1
            |
 (-) O4  O2
         \ /
          P3 (+)
         / \
   (-)O5   O6
            |
           C7-----etc..... (rest of lipid)

The same region in the POPC topology looks something like this:

           C6-----etc (choline group)
            |
(-) O9  O7
         \ /
          P8 (+)
         / \
 (-)O10  O11
            |
           C12---------etc .... (rest of lipid)


Thus, only one new atom type, a H-atom was introduced. I replaced the
C6 of POPC by H1 on POPA. O7, P8, O9 were replaced by O2, P3, O4 and
so on (so atom indices differ by a factor of 5)

In the following, I have copy-pasted the topology of POPA and POPC for
the above regions. The topology for POPC is the one on Tieleman's
webpage, of course. My topology is directly derived from it.

###############################################################################################
###### FOR POPA #######
 [ atoms ]
;   nr    type   resnr  residu    atom    cgnr        charge          mass
     1     LHO       1    PAPA     H1        1          0.400
1.0080 ; this is the new atom type, I defined it in the lipid.itp file
and also described non-bonded interactions and pairtypes for this
atom.
     2     LOS       1    PAPA     O2        1         -0.600   15.9994  ;
     3     LP        1    PAPA     P3        1          1.500   30.9738
     4     LOM       1    PAPA     O4        1         -0.800   15.9994
     5     LOM       1    PAPA     O5        1         -0.800   15.9994
     6     LOS       1    PAPA     O6        1         -0.700
15.9994	; qtot:-1
     7     LC2       1    PAPA     C7        2          0.400   14.0270	

 [bonds ]
;  ai    aj funct
       1       2       1 0.10000E+00 0.31380E+06
       2       3       1 0.16100E+00 0.25100E+06
       3       4       1 0.14800E+00 0.37660E+06
       3       5       1 0.14800E+00 0.37660E+06
       3       6       1 0.16100E+00 0.25100E+06
       6       7       1 0.14300E+00 0.25100E+06
       ....etc
[ pairs ]
;  ai    aj funct
   1     4      1
   1     5      1
   1     6      1
   2     7      1
   3     8      1
   4     7      1
.... etc ....

[ angles ]
;  ai    aj    ak funct
      1        2       3       1 0.10960E+03 0.39750E+03
      2        3       4       1 0.10960E+03 0.39750E+03
      2        3       5       1 0.10960E+03 0.39750E+03
      2        3       6       1 0.10300E+03 0.39750E+03
      3        6       7       1 0.12000E+03 0.39750E+03
      4        3       5       1 0.12000E+03 0.58580E+03
      4        3       6       1 0.10960E+03 0.39750E+03
      5        3       6       1 0.10960E+03 0.39750E+03

[ dihedrals ]
;  ai    aj    ak    al funct   phi0     cp     mult
    1     2     3     6     1    0.0     1.05   3
    1     2     3     6     1    0.0     3.14   2
    2     3     6     7     1    0.0     1.05   3
    2     3     6     7     1    0.0     3.14   2
#####################################################################################

#####################################################################################
###### FOR POPC #######

 [ atoms ]
;   nr    type   resnr  residu    atom    cgnr        charge          mass

; 1 to 5 are from the choline group

     6     LC2       1    POPC      C6       1         0.4000
14.0270	; qtot:1.0
     7     LOS       1    POPC      O7       1        -0.800
15.9994	; qtot:0.54
     8     LP        1    POPC      P8       1         1.700
30.9738	; qtot:2.3
     9     LOM       1    POPC      O9       1        -0.800
15.9994	; qtot:1.5
    10     LOM       1    POPC     O10       1        -0.800
15.9994	; qtot:0.7
    11     LOS       1    POPC     O11       1        -0.700   15.9994	; qtot:0
    12     LC2       1    POPC     C12       2         0.400
14.0270	; qtot:0.08

[ bonds ]
;  ai    aj funct
       6       7       1 0.14300E+00 0.25100E+06
       7       8       1 0.16100E+00 0.25100E+06
       8       9       1 0.14800E+00 0.37660E+06
       8      10       1 0.14800E+00 0.37660E+06
       8      11       1 0.16100E+00 0.25100E+06
      11      12       1 0.14300E+00 0.25100E+06

[ pairs ]
;  ai    aj funct
    4     7     1
    5     8     1
    6     9     1
    6    10     1
    6    11     1
    7    12     1

[ angles ]
;  ai    aj    ak funct
        6       7       8       1 0.12000E+03 0.39750E+03
       7       8       9       1 0.10960E+03 0.39750E+03
       7       8      10       1 0.10960E+03 0.39750E+03
       7       8      11       1 0.10300E+03 0.39750E+03
       8      11      12       1 0.12000E+03 0.39750E+03
       9       8      10       1 0.12000E+03 0.58580E+03
       9       8      11       1 0.10960E+03 0.39750E+03
      10       8      11       1 0.10960E+03 0.39750E+03


[ dihedrals ]
;  ai    aj    ak    al funct   phi0     cp     mult
    6     7     8    11     1    0.0     1.05   3
    6     7     8    11     1    0.0     3.14   2
    7     8    11    12     1    0.0     1.05   3
    7     8    11    12     1    0.0     3.14   2
###############################################################################################