CHARMM Development Project Forums User Discussion & Questions Parameter Set Discussion H2PO4- toppars proposed. Feedback please!

Hey everybody! I've been looking quite a bit into modeling dihydrogen phosphate, H2PO4-, for the sake of ionizing a biomolecular system under physiologically relevant conditions. After reading posts and suggestions on this forum, VMD-L, and NAMD-L, I came up with the following strategy:
1) Create an H2PO4- fragment with VegaZZ ( http://www.vegazz.net/ ), and use UCSF Chimera ( http://www.cgl.ucsf.edu/chimera/ ) to create a .mol2 file.
2) Use SwissParam ( http://swissparam.ch/ ) to create CHARMM-compatible toppars.
3) Edit charges to match those published by Kaplus for H2PO4- ( The missing link between thermodynamics and structure in F1-ATPase ).
4) Refer to methylphosphate (RESI MP_1) topology and corresponding parameters to appropriately edit bond lengths and dihedrals ( top_all27_lipid.rtf and par_all27_lipid.prm in /toppar_history/c28a1_to_c31a1, the subdirectory of the latest-posted CHARMM toppars).

Below are the pasted topology and parameters that I ended up with. Could you please let me know how these look? I'm relatively new to MD and work mainly with VMD and NAMD while using CHARMM parameters. Any feedback would be really appreciated.

Topology:

MASS 201 PO4 30.973800
MASS 202 O2CM 15.999400
MASS 203 OR 15.999400
MASS 204 HOCO 1.007940

AUTOGENERATE ANGLES DIHE
DEFA FIRS NONE LAST NONE

RESI PH1 -1.000 ! H1
GROUP ! /
ATOM P PO4 1.2500 ! O1
ATOM O2 O2CM -0.7950 ! |
ATOM O1 OR -0.7200 ! O2==P==O3 (-)
ATOM H1 HOCO 0.3900 ! |
ATOM O3 O2CM -0.7950 ! O4
ATOM O4 OR -0.7200 ! /
ATOM H3 HOCO 0.3900 ! H3
BOND P O2
BOND P O1
BOND P O3
BOND P O4
BOND O1 H1
BOND O4 H3
!IMPH P O4 O1 O2
!IMPH P O2 O1 O3
IC O2 P O1 H1 1.65 98.26 25.66 116.99 1.00
IC O2 P O4 H3 1.65 95.82 154.31 119.06 0.99
IC O1 P O4 H3 1.63 104.37 -105.60 119.06 0.99
IC H1 O1 P O3 1.00 116.99 158.31 116.73 1.46
IC H1 O1 P O4 1.00 116.99 -72.56 104.37 1.65
IC O3 P O4 H3 1.46 115.75 24.12 119.06 0.99
!IC O4 O1 *P O2 0.00 0.00 120.00 0.00 0.00
!IC O2 O1 *P O3 0.00 0.00 -120.00 0.00 0.00

END

Parameters:

BONDS
PO4 O2CM 580.000 1.4800
PO4 OR 237.000 1.6370
OR HOCO 545.000 0.9600

ANGLES
O2CM PO4 OR 98.900 108.2300
O2CM PO4 O2CM 120.000 120.0000
OR PO4 OR 127.307 99.3110
PO4 OR HOCO 30.0 115.0 40.0 2.30

DIHEDRALS
O2CM PO4 OR HOCO 0.30 3 0.00 !X OHL PL X 0.30 3 0.00 ! terminal phosphate (from CHARMM27 lipid params)
OR PO4 OR HOCO 0.30 3 0.00 !X OHL PL X 0.30 3 0.00 ! terminal phosphate

!IMPROPER !There don't appear to be impropers in MP_1 topology, so I got rid of these SwissParam-generated ones
!PO4 OR OR O2CM 0.000 0 0.00
!PO4 O2CM OR O2CM 0.000 0 0.00

!Nonbonded below. Should be invisible in NAMD.
NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch -
cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5
PO4 0.000000 -0.585000 2.150000
O2CM 0.000000 -0.120000 1.700000
OR 0.000000 -0.152100 1.770000
HOCO 0.000000 -0.046000 0.224500

I tried these in NAMD, and the resulting H2PO4- appears to be well-behaved.

Thanks for the feedback,
JackB

Maybe you should have a look at the MacKerell lab site, and use the tools suggested there instead. You may also find some more recent parameter files as well.

The VMD/NAMD folks do not seem to be the best source of knowledge about the best way to develop CHARMM parameters, IMHO.

There are several issues with your strategy. More fundamentally, you're making it way to complicated. I personally would just make a copy of RESI MP_1, edit it into dihydrogen phosphate (with same atom types and charges on the new -OH group as on the existing one), and transfer 1 angle and 1 dihedral from RESI MP_1. Done. Simple is good.

Dear All:
I would like to use the phosphate ions H2PO4- and HPO4 2- in my MD simulations. According to your suggestions, I make a copy of MP_1 and MP_2, and edit them into H2PO4- and HPO4 2-, respectively.
RESI MP_2 -2.00 ! Methylphosphate, dianionic
GROUP !
ATOM P1 PL 1.10 !
ATOM O1 OSL -0.40 ! H11
ATOM O2 O2L -0.90 ! |
ATOM O3 O2L -0.90 ! H13--C1--H12
ATOM O4 O2L -0.90 ! |
GROUP ! O1
ATOM C1 CTL3 -0.27 ! |
ATOM H11 HAL3 0.09 ! (-) O4--P1(+)--O3 (-)
ATOM H12 HAL3 0.09 ! |
ATOM H13 HAL3 0.09 ! O2 (-)

RESI MP_2 -2.00 ! Methylphosphate, dianionic
GROUP !
ATOM P1 PL 1.10 !
ATOM O1 OSL -0.40 ! H11
ATOM O2 O2L -0.90 ! |
ATOM O3 O2L -0.90 ! H13--C1--H12
ATOM O4 O2L -0.90 ! |
GROUP ! O1
ATOM C1 CTL3 -0.27 ! |
ATOM H11 HAL3 0.09 ! (-) O4--P1(+)--O3 (-)
ATOM H12 HAL3 0.09 ! |
ATOM H13 HAL3 0.09 ! O2 (-)

RESI HP2 -1.00 ! H2PO4 -
GROUP !
! atom order for molvib
ATOM O1 OHL -0.77 !
ATOM P1 PL 1.50 !
ATOM O2 OHL -0.77 ! H1
ATOM O3 O2L -0.82 ! /
ATOM O4 O2L -0.82 ! O1
ATOM H1 HOL 0.34 ! |
ATOM H2 HOL 0.34 ! O4==P1==O3 (-)
! |
! O2
! \
! H2


RESI HP1 -2.00 ! HPO4 2-
GROUP !
ATOM P1 PL 1.10 !
ATOM O1 OHL -0.74 !
ATOM O2 O2L -0.90 !
ATOM O3 O2L -0.90 ! H1
ATOM O4 O2L -0.90 ! |
ATOM H1 HOL 0.34 ! O1
! |
! (-) O4--P1(+)--O3 (-)
! |
! O2 (-)

Due to both the O1 and O2 is OHL, The angle of O1-P1-O2 does not exist in the par_all27_lipid.prm file. So according the invariable angle in prm file we get the angle of O1-P1-O2 in H2PO4- is 102.54 degree. However, I am not completely sure about how to modify the charges of each atom after the editing. For H2PO4 -, I set the charges of P1, O3, O4, H2 same as in MP_1, and transfer the charge of O1 from -0.62 to -0.77 and the charge of O2 from -0.68 to -0.77 so that the total charge of dihydrogen phosphate is -1. For HPO4, I transfer Methyl to a hydrogen atom in MP_2, no angle in the system is changed. However I change the charge of O1 from -0.40 to -0.74 and the charge of hydrogen atom is the same as the hydrogen atom in H2PO4 - (0.34), the charge of the other atom does not change. Thus the total charge is -2.

Can anyone give a comment about whether this is a proper way to get what I need?

• Please learn to use the [code] tags (or the attachment function for larger files) to avoid messing up the formatting in your post.
• As indicated in my previous post, for HP1, it's more important to preserve the charges on the OH group (both oxygen and hydrogen) in MP_1 than the charge on the P (which is buried).
• The phosphate in CGenFF is actually an updated version of the one in the other parts of the additive force field. While it probably won't work well in a nucleic acid backbone, it's supposed to be more representative for a small isolated phosphate.
• For HP2 , one cannot just assign charges and parameters by analogy; optimization is required. This is because there do not exist OH charges and parameters for a double negative group, and they're expected to be qualitatively different from anything that is already present in the force field.

This is not a good English sentence and I have no idea what you mean, but it looks like the value 102.54 only occurs in IC tables. If you constructed an angle parameter based on a value from an IC table, then that's completely wrong. See this FAQ entry (and there's no guarantee whatsoever that the IC table geometry is even a relevant minimized structure). Better would be to simply copy the analogous parameter that contains the phosphate ester oxygen. Same for the missing dihedral parameter.

Dear Kenno:

Thanks for the prompt reply. I will try to be clearer for the following question.

1) "...For HP2, ...optimization is required...."

Would you kindly point out a protocol that I can follow to optimize the HP2 parameters? I've never done it before.

Sorry for my late reply, I was a bit tied up with the CGenFF release. By now, you surely must have Googled it already and found the tutorial and FAQ...

Another thing you could try is FFTK, but be aware that it minimizes of the rotation angle of the water molecules when doing the water interactions, which is incorrect because it tends to give rise to secondary interactions, most importantly hydrogen bonds. If you, as a user, can find a way to suppress the minimization of the angle but not the minimization of the distance (for example by doing a small part of the procedure outside of FFTK), you should be perfectly fine.

Apologies for my contribution so long after the initial discussion; I'm only stumbling upon this thread via a recent link from a related conversation on NAMD-L.



I'm one of the developers of ffTK and wanted to note that suppressing the rotational angle of the water during the QM optimization, as suggested by Kenno in the above quote, is relatively straightforward. Doing this requires only minor manipulation of the Gaussian input files generated by ffTK to remove the rotational degree of freedom.

• Setup the water interaction calculation in ffTK as usual
• Open the resulting *.gau file in a text editor
• In the Z-matrix portion of the molecular description, replace "dih" variable with whatever rotational angle you desire for the water molecule
• Below the Z-matrix portion where the optimized parameters are initialized, remove the line "dih 0.0"

After you make these modifications, you can load the input file into VMD using ffTK's "Load GAU Files" button to visually inspect the water positioning; specifically noting the now fixed rotational angle. If satisfied with the setup, run the Gaussian calculation and then apply the resulting LOG file to ffTK's charge optimization workflow as per usual.

That's it!

Regards,
Christopher Mayne

Thanks, Chris, for addressing this concern. We really appreciate your commitment to ffTK's CHARMM compatibility. BTW, AFAIK we don't have any policy against posting in old threads; anything goes as long as it's relevant enough to people who find the discussion through some search function or link (which your post certainly is).