CHARMM36 force field charge of ribonucleotide residues

User discussions
1

GROMACS version:2019
GROMACS modification: No
Here post your question
I am trying to simulate a protein-RNA complexed with a ligand where the ribonucleotide residues of which the total charges is -1 in the CHARMM36 force-field which makes the whole RNA protein complex highly negative (about -224) and to neutralize the system Sodium(Na+) was used as counter-ions and when the system is made to undergo dynamics with all of the biomolecule positionally restrained, these counter-ions which before equilibration were fine,these Na+ ions swims near the ligand and this phenomenon happens when Particle Mesh Ewald(PME) is used to calculate long range interactions but these same ions(Na+) don’t swim near the ligand when plain cut-off is used. Can there be some other way round where I can read these RNA residues as neutral residues such that I don’t need to add this many number of counter-ions which influence my ligand with that high number of positive charges. I have checked in the CHARMM force field such rtp entries for neutral residues are not provided and the residues which are there are only for the Nitrogenous bases.
These are the contents of my mdp file:

title                   = NVT equilibration 
define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 2500000     ; 2 * 2500000 = 5000 ps (5ns)
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Nonbonded settings 
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein  Non-Protein			; two coupling groups 
tau_t                   = 0.1 	   0.1			; time constant, in ps
ref_t                   = 300	   300   ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed

The RNA residues are as such attached here:

A1
A11628×806 75.8 KB

Thanking in advance.
Kindly do send in your suggestions.
Regards,
Pallav

2

It makes no physical sense to do this, nor is the force field parametrized to include neutral phosphate groups in the backbone of nucleic acids. The pKa is ~2 so under all relevant biological conditions, the phosphate groups carry a -1 charge.

The issue you’re probably facing is that the effective concentration of the ions is enormous for this kind of system. I’d wager that even neutralizing ions are in the high mM or near M concentration. You may want to increase the size of your box to dilute these ions to something a bit more reasonable (of course this is more expensive to simulate, but also probably more realistic).

Plain cutoffs are entirely inaccurate and should never be used for electrostatics. Whatever you observe with a plain cutoff is fictitious.

1 Like