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!
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.
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.