GROMACS version:2020.4

GROMACS modification: No

Greetings to all!

I am performing a 100 ns MD trajectory calculation for the complex of two proteins using CHARMM36 FF (2020). As a result, I get a model of a complex with a rather strongly distorted protein geometry (violation of planarity for planar groups, etc.), so coot is quite unhappy with it :). I use distance restrains for some important hydrogen bonds (“bond type 6” inside one protein chain and “pull code” for intermolecular interactions). How can the geometry of the model, obtained during the MD calculation, be optimized? Are there some tools for this in gromacs? Alternatively, do I need to set the simulation initially by imposing some additional restrictions to achieve good geometry of the amino acid residues and protein main chain? Thank you very much!

Hi,

I get a model of a complex with a rather strongly distorted protein geometry (violation of planarity for planar groups, etc.)

Which molecules/functonal groups have planar distortion? If those are standard residues, then it could be that you did not set your mdp parameter in line with CHARMM36 force field. Or in alternative you did not properly relax the starting structure in order to release possible steric clashes.

Best regards

Alessandra

Dear Alessandra! Thank you very much for the answer.

Well these are the standard residues like His, Trp, Tyr etc. What do you mean to set mdp parameter in line with CHARMM36 force field? My md.mdp looks like this:

title = Protein MD simulation

define = -DDISREA ; distance restrain

; Run parameters

integrator = md ; leap-frog integrator

nsteps = 50000000 ; 2 * 50000000 = 100000 ps (100 ns)

dt = 0.002 ; 2 fs

; Output control

nstenergy = 5000 ; save energies every 10.0 ps

nstlog = 5000 ; update log file every 10.0 ps

nstxout-compressed = 5000 ; save coordinates every 10.0 ps

; Bond parameters

continuation = yes ; continuing from NPT

constraint_algorithm = lincs ; holonomic constraints

constraints = h-bonds ; bonds to H are constrained

lincs-iter = 1 ; accuracy of LINCS

lincs-order = 4 ; also related to accuracy

; Neighbor searching and vdW

cutoff-scheme = Verlet

ns-type = grid ; search neighboring grid cells

nstlist = 20 ; largely irrelevant with Verlet

rlist = 1.2

vdwtype = cutoff

vdw-modifier = force-switch

rvdw-switch = 1.0

rvdw = 1.2 ; short-range van der Waals cutoff (in nm)

; Electrostatics

coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics

rcoulomb = 1.2

pme-order = 4 ; cubic interpolation

fourierspacing = 0.16 ; grid spacing for FFT

; Temperature coupling

tcoupl = V-rescale ; modified Berendsen thermostat

tc-grps = Protein Non-Protein ; two coupling groups - more accurate

tau-t = 0.1 0.1 ; time constant, in ps

ref-t = 310 310 ; reference temperature, one for each group, in K

; Pressure coupling

pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT

pcoupltype = isotropic ; uniform scaling of box vectors

tau-p = 2.0 ; time constant, in ps

ref-p = 1.0 ; reference pressure, in bar

compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1

; Periodic boundary conditions

pbc = xyz ; 3-D PBC

; Dispersion correction is not used for proteins with the C36 additive FF

DispCorr = no

; Velocity generation

gen-vel = no ; continuing from NPT equilibration

; Pull code for intermolecular distance restrains

pull = yes

pull-ngroups = 5

pull-ncoords = 3

pull-group1-name = a_496

pull-group2-name = a_1234

pull-group3-name = a_1235

pull-group4-name = b_3440

pull-group5-name = b_3443

pull-coord1-type = umbrella

pull-coord2-type = umbrella

pull-coord3-type = umbrella

pull-coord1-geometry = distance ; simple distance increase

pull-coord2-geometry = distance ; simple distance increase

pull-coord3-geometry = distance ; simple distance increase

pull-coord1-groups = 1 5

pull-coord2-groups = 2 4

pull-coord3-groups = 3 5

pull-coord1-dim = Y Y Y

pull-coord2-dim = Y Y Y

pull-coord3-dim = Y Y Y

pull-coord1-rate = 0.0 ; not pull, just distance restraint at a reference distance bw 2 groups

pull-coord2-rate = 0.0 ; not pull, just distance restraint at a reference distance bw 2 groups

pull-coord3-rate = 0.0 ; not pull, just distance restraint at a reference distance bw 2 groups

pull-coord1-k = 1000 ; kJ mol^-1 nm^-2

pull-coord2-k = 1000 ; kJ mol^-1 nm^-2

pull-coord3-k = 1000 ; kJ mol^-1 nm^-2

pull-coord1-init = 0.30

pull-coord2-init = 0.28

pull-coord3-init = 0.28

pull-nstxout = 0

pull-nstfout = 0

Is there something wrong with it that doesn’t match CHARMM36 FF?

Hi,

I guess you use the pull code between atoms, the applied potential is in “competition” with the force field potential, thus the use of not properly calibrate values for pulling could results in underised distortion of other interactions (e.i aromatic ring).

Best regards

Alessandra

Hi!

OK, but the geometry is far from ideal even for a model calculated without any additional constraints (pull code, distance restrains or position restrains). So, may be there is a way to correct the geometry of the resulting model in gromacs? Or the way to make an averaged model with good geometry for the entire trajectory or part of it?

Thanks!

Yours sincerely

Oleg

Hi,

There are several steps to consider:

- generation of the starting structure (e.i from experimental structure or from homology modelling, or).
- once you have a reasonable starting structure, you can energy minimize it (to assure that you are in one of the minima of the energy landscape and to eliminate steric clash)
- you can equilibrate the structure to get to the desidered temperature and pressure (using using position restrain for the macromolecule)
- then you can start to sample configuration

Did you see distorsion in the geometry after all those steps? Then I suggest to check carefully the starting structure.

Concerning the average structure, it is normal that an time-averaged structure show distortion in the amino-acid side chain geometry. For example the time-average of a rotating aromatic ring will be points on a line (not a ring).

Best regards

Alessandra

Thank you so much for your answers. I will analyze the structures in detail at all stages.

With best wishes

Oleg.