During the past ~4.5 years of my PhD, I have written an MD code, based largely on LAMMPS and NAMD, however reparallelized to suit my computational scenario better. One important thing to mention possibly for this discussion is that the code runs in SPDP mode to take advantage of the wide vector registers on xeon phi coprocessors. I am currently working with acetaminophen in water, however I cannot seem to rationalize the behavior I am observing with the temperature vs time data. I am using GAFF force field for acetaminophen and tip3p for water. The langevin thermostat is being used for temperature maintenance, and is being ran with a relaxation time of 2 ps. The shake algorithm is being used to constrain the motion of all bonded hydrogens. Currently, the system is a 2.75nm acetaminopehn nanocrystal embedded in a previously equilibrated box of acetaminopen + tip3p water. These are fairly large simulation I am trying to run, and consist of about 2600 acetaminophen molecules in about 30K water molecules for about 150K total atoms.  I have noticed that the system equilibrates at ~321 K, even though the set point is 320 K. I have figured out through testing that the set point offset is somehow related to the shake algorithm being applied to the solid nanocrystal. Currently,  I have the following test cases. I have more plots to show, but the website wont let me post more than 2. solid liquid shake : this is the case where shake is applied to all hydrogens in the system (solid acetaminophen nanocrystal + liquid acetaminophen +. The system was ran in the NPT ensemble with langevin thermostat (2 ps relaxation time) and berendsen barostat (1 ps relaxation time). The timestep was 1.0fs. As you can see from the first image, the system equilibrates at ~321 K, and is clearly above the 320 K set point.fluid-system shake: To test if there is a bug in my shake code that appears only when applied to acetaminophen, I ran an equilibration of a fluid mixture of acetaminophen and water. This system has no embedded solid sphere, and was generated using the insert-molecules and solvate command from gromacs. The following image shows the system equilibrating at 320 K correctly. Not sure what is so specific about the embedded sphere case thats making it bug out.one final thing of note is that I have tried varying the friction constant of the langevin thermostat. Increasing the friction accelerations temperature convergence to the steady state value, but does not affect the steady state value itself (i.e. 321 K)My specific questions are (1) Do others share any concern regarding the system equilibrating at a slightly higher (~0.8-1.0 K) temperature than the set point. I.E. Is it a general phenomenon for the shake algorithm to possible induce temperature set point offset? Is this kind of offset even considered strange? (2) Does anyone see something I am missing? I have been pounding my head on a desk for about a week now. Do others have thoughts as to why the shake algorithm, when applied to my embedded solid nanocrystal, causes the system to miss the set point, but not when the shake algorithm is applied to an acetaminophen + water fluid mixture? Any help is greatly appreciated. I am the only one in my research group who does MD, and cant ask these kind of very detailed questions to anyone else.
temperature set point maintenance with Langevin thermostat + shake algorithm
molecular dynamics;molecular mechanics;lammps