GROMACS version:2019
GROMACS modification: No
Here post your question
Dear Gromacs developers and users,
I working on a RNA-protein-ligand complex where the system is highly negatively charged and to neutralise the system I am using Na+ ions(~ 224 Na+ ions), after the system is well neutralised and minimized when i move to equilibration dynamics with whole of the complex including the ligand being restrained and use Ewald electrostatics to calculate the long range interactions it leads to artifacts (which shouldn’t have been the case as the system is well neutralised) with Na+ ions moving into regions with low dielectric ( near the ligand)
As suggested by Justin,
Tried increasing the box size which obviously is quite expensive to simulate, the similar sort of phenomenon is observed even when the box size is increased.( the previous dimensions were 2197 cubic nm which I increased to 3375 cubic nm)
These are the parameters I am using for equilibration
title = NVT equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
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
Thanking you in adavnce.
Kindly do provide your suggestions.
Regards,
Pallav
Attraction of ions to a highly charged species such as a large RNA is not unexpected. Can you provide some images of what you are defining as an artifact?
This is the print-out while generating the tpr file:
Command line:
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index_coupl.ndx -o nvt.tpr
Ignoring obsolete mdp entry 'title'
Setting the LD random seed to 769665142
Generated 100032 of the 100128 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 1
Generated 65937 of the 100128 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'Ion_chain_7'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_A'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_A2'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a2'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a3'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a4'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a5'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_a6'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_c'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_e'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_i'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'Protein_chain_k'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_v'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_x'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'RNA_chain_y'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'PAR'
turning H bonds into constraints...
Excluding 2 bonded neighbours molecule type 'SOL'
turning H bonds into constraints...
Excluding 2 bonded neighbours molecule type 'SOL'
Excluding 1 bonded neighbours molecule type 'NA'
turning H bonds into constraints...
Setting gen_seed to 1952425801
Velocities were taken from a Maxwell distribution at 300 K
Removing all charge groups because cutoff-scheme=Verlet
NOTE 1 [file topol.top, line 90]:
The bond in molecule-type RNA_chain_A2 between atoms 308 C2 and 309 O2
has an estimated oscillational period of 1.9e-02 ps, which is less than
10 times the time step of 2.0e-03 ps.
Maybe you forgot to change the constraints mdp option.
Number of degrees of freedom in T-Coupling group Protein_RNA_PAR_CYT_GUA_ADE_URA is 30360.79
Number of degrees of freedom in T-Coupling group Water_and_ions is 405585.22
Determining Verlet buffer for a tolerance of 0.005 kJ/mol/ps at 300 K
Calculated rlist for 1x1 atom pair-list as 1.234 nm, buffer size 0.034 nm
Set rlist, assuming 4x4 atom pair-list, to 1.200 nm, buffer size 0.000 nm
Note that mdrun will redetermine rlist based on the actual pair-list setup
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 84x84x84, spacing 0.155 0.155 0.155
Estimate for the relative computational load of the PME mesh part: 0.04
NOTE 2 [file nvt.mdp]:
This run will generate roughly 24550 Mb of data
Is this the final snapshot from the equilibration? How close are the ions to start (e.g. in em.gro)? Can you please address the question above regarding your force field and the nonbonded settings? If you’re using CHARMM, again those settings are wrong and will lead to incorrect behavior of the force field.
Yes, I have changed the non-bonded settings as advised which I had ignored previously and had attached the final equilibration snap in the reply to you 9 hours ago. I am using CHARMM36 force field, after changing the non-bonded settings the scenario has much improved with ions not attaching with the ligand as I was observing in previous simulations.
