Re: [AMBER] [Sender Not Verified] Re: Issues with simulating 5'-phosphorylated RNA strands

From: Michael Bruist via AMBER <>
Date: Mon, 19 Jun 2023 12:53:04 -0400


I have run across similar problems and came up with two different methods
that seem to have stabilized the phosphates.

Method 1

Dissociate the problematic hydrogen, giving the terminal phosphate a charge
of –2.

Since the pKa of that group is close to 7, this is not an abnormal choice.
The atom type of the dissociated hydrogen was changed and charged on P and
OP1, OP2, and OP3 were altered, all in leap.

source leaprc.RNA.OL3

source leaprc.water.tip3p

loadOff terminalphos.lib

rna=loadPDB ../../TLS1_HPO3_old/TLS1_clipped.chimera.pdb

remove rna.1 rna.1.36 # removes the HOP from residue

set rna.1.35 type "O2" # changed atom type from OH
to O2 (OP3)

set rna.1.35 charge -0.9534 # changes charge of oxygen with
changed type

set rna.1.1 charge 1.1662 # changes charge on phosphorus

set rna.1.2 charge -0.9534 # changes charge on OP2

set rna.1.3 charge -0.9534 # changes charge on OP1

check RNA

Method 2

Increase the force constraint on the P - OP3 - HO bond angle.

In all cases where I have been able to see images before the crash, this
angle has been less than 90°. I increased it approximately 3x, to give it
the same value as the O - P - P bonds, so it is not completely abnormal.

With this value, a 100 nsec run of a 4-way junction molecule with four 5’
phosphates gave one P -O -H angle that dropped to 95°.

These modifications were made using parmed on the parm and crd file created
in leap with the problematic phosphate(s).

Note parmed will only recognize an angle listed in one direction:

printAngles :1.P :1.OP3 :1.HP3 # This is recognized

printAngles :1.HP3 :1.OP3 :1.P # gives blank information with no

Get started:

parmed -p ../build3/TLS1_clipped_K_tmp3.parm7 -c

printAngles :1.P :1.OP3 :1.HP3 # provides the current
properties of this angle

setAngle :1.P :1.OP3 :1.HP3 140.0 108.5 # 140 is the new spring
constant; 108.5 is the original equilibrium angle

printAngles :1.P :1.OP3 :1.HP3 # The new value was accepted

setAngle :11.P :11.OP3 :11.HP3 140.0 108.5 # The other three phosphates
are modified.

setAngle :21.P :21.OP3 :21.HP3 140.0 108.5

setAngle :31.P :31.OP3 :31.HP3 140.0 108.5

printDetails .OP3 # Various ways to give the program
a chance to say something is wrong.



outparm TLS1_short_HPO3_k140.parm7 # This prints out
the new parm file. Be sure to give it a distinctive name. This file will
work with previous coordinate files because no atoms were added or deleted.


On Thu, Mar 23, 2023 at 10:57 AM Harvey Dale via AMBER <>

> Dear David and Thomas,
> Thanks very much for the advice. As both of you have described, setting a
> small vdW radius (0.3) and epsilon (0.05) for the phosphate hydrogen
> appears to have resolved the problems with the instabilities I was
> observing; I must confess I have not yet tried stiffening the P-O-H angle,
> but suspect this too will yield the same result.
> Will certainly look to validate this approach as I run more simulations
> down the line.
> Thanks again,
> Harvey
> ________________________________________________
> Dr Harvey J. A. Dale
> Research Fellow
> Royal Commission for the Exhibition of 1851
> John Henry Coates Research Fellow
> Emmanuel College | Cambridge | CB2 3AP |
> Postdoctoral Researcher
> MRC Lab of Molecular Biology | Cambridge | CB2 0QH
> _______________________________________________
> -----Original Message-----
> From: David A Case via AMBER <>
> Sent: Wednesday, March 22, 2023 3:07 PM
> To: Steinbrecher, Thomas <>; AMBER Mailing
> List <>
> Subject: Re: [AMBER] [Sender Not Verified] Re: Issues with simulating
> 5'-phosphorylated RNA strands
> On Wed, Mar 22, 2023, Steinbrecher, Thomas via AMBER wrote:
> >
> >my experience here may be outdated, but the stability issues you
> >mention are (or were some years ago) typical for MD simulations of
> >protonated phosphates. As you describe, the imbalance between
> >electrostatics and vdW repulsion is the issue here. A normal hydroxyl H
> >does not need a LJ center, because it is sufficiently protected by both
> >the O-LJ volume plus, and this is key, the relatively stiff angle
> >bending in e.g. C-C-O groups, so the H is kept away from the heavy
> >atoms three bonds away. In contrast, a phosphate has O-P-O angle terms
> >which are softer plus relatively high charges on both H and O, so the
> >collapse you describe between atoms 1 and 4 in a O-P-O-H substructure
> >can occur. The ad hoc solution of giving the H non-zero vdW parameters
> >will make the simulation stable (most likely) but in principle you'd be
> >on the hook to show that this still describes the phosphate group well
> >in the environment you simulate it in. The cowardly solution would be
> >to claim that the phosphate has pK around 7 anyway and in the presence
> >of all those nucleic acid associated metal cations, you may as well
> simulate it deprotonated.
> Thomas' post seems spot-on to me. One other thing that I remember doing
> is to increase the force constant for the P-O-H angle in the monophosphate:
> this will help keep the proton from finding a neighboring hydrogen. I'd
> double the force constant -- having a slightly-too-stiff angle is not
> really going to do anything bad to your sampling.
> Let us know if either of both of these ideas (non-zero LJ [try the value
> for hydrogens bonded to nitrogen] or stiffer angle) help. We should update
> the parameters so that others don't see the problems you are having.
> ...thanks for the detailed reports....dac
> _______________________________________________
> AMBER mailing list
> _______________________________________________
> AMBER mailing list
AMBER mailing list
Received on Mon Jun 19 2023 - 10:00:03 PDT
Custom Search