Abrupt temperature change during heating

GROMACS version:2023.1
GROMACS modification: Yes/No

I am using the annealing code to warm up my system:

tcoupl                  = V-rescale
tc_grps                 = SOLU SOLV
tau_t                   = 1.0 1.0
ref_t                   = 303 303
;
pcoupl                  = Parrinello-Rahman
pcoupltype              = isotropic
tau_p                   = 5.0
compressibility         = 4.5e-5
ref_p                   = 1.0
;
constraints             = h-bonds
constraint_algorithm    = LINCS
continuation            = yes
;
nstcomm                 = 100
comm_mode               = linear
comm_grps               = SOLU SOLV
;
; Simulated annealing
annealing	= single single
annealing_npoints	= 9 9
annealing_time 	= 0   500 1000 1500 2000 3000 5000 9500 10000  0   500 1000 1500 2000 3000 5000 9500 10000
annealing_temp	= 303 303 453  453  503  503  553  553  303    303 303 453  453  503  503  553  553  303

There were no errors in either the grompp or mdrun steps. But when I checked the temperature, I found that when the system heated up above 500K, a strange temperature curve appeared:

When these peculiar temperature changes occur, multiple glycan branches (the system is a glycosylated protein, and the protein backbone is constrained by posres) appear in the trajectory to simultaneously swing in one direction.

Is this phenomenon reasonable? If it doesn’t, how can I improve my simulation? Thank you!

Updates:

I refer to this mdp parameter: http://www.mdtutorials.com/gmx/membrane_protein/Files/anneal_npt.mdp

I think it may be caused by the large volume change caused by the heating under the NPT ensemble. Why does this official simulation take the approach of heating up under NPT?

Indeed, heating in an NPT ensemble may/will cause the box size to change. In general, I think that is what you want. You would get a very high pressure if you would heat it under NVT conditions. But heating with Parrinello-Rahman is probably not a great idea, since you can get oscillations. I would recommend using c-rescale instead. For most purposes, there’s no reason to use anything else.

On the other hand, you are heating the system quite a lot. If there’s water in the system it would start boiling, which might not be great in combination with NPT.

You also describe that some glycan branches swing simultaneously in the same direction. This makes me think of possible issues with center of mass motion removal. I think it would be interesting to see if there would be any difference if you only used comm_grps = system. In combination, I think you should also test using only one temperature coupling group.

Thanks for replying!

I checked the details of the trajectory: the moment when glycan swing in one direction is exactly the moment when the whole system undergoes a dramatic volume expansion.

I’m actually trying to reproduce a simulation done by someone else, where they observed that a particular glycan would undergo a conformational change during annealing. But they don’t provide specific warming information. Although 550K is too high for the protein system, I did not observe the reported conformational changes at low temperatures. That’s why I try to warm up more than 550K.

I think in this case I should try to use the NVT ensemble for warming while put constraints on the protein backbone.

It might be that the conformational changes cause the volume changes rather than the other way around. I.e., that what you observe are transitions upon heating.