Thanks everyone for your reply. I see your point that it could just be some error cancellation from constraints and parameters.
I was previously using gromacs-5.1.4, and changed to gromacs-2018, and the drift seems to be much smaller! I get something similar to what Mark had posted -0.000025 kJ/mol/ps/atom for 1000 SPCE waters with settles at 298K for single precision. With double precision, there is no noticeable drift for a short test run of only 200 ps. Any idea why this is? Is there some default setting that changed? Mark, regarding your earlier questions: My application is to calculate very precise and accurate dynamic properties (diffusion, viscosity) of NaCl in water. There are a number of works in literature where people use NVT simulations to calculate these properties, but I wanted to check if the thermostat affected the dynamics any by running an NVE simulation. Also I want to better understand how NVE simulations work and all the controlling factors. >From figure 3.5 in the reference manual, I think 10^-5 drift (kJ/mol/ps) is acceptable. Can you help me better understand what verlet-buffer-tolerance does? My understanding from the documentation is that it determines the width of the buffer outside of the cutoff for the neighborlist, but I am not clear why energy can be added or lost by manipulating this parameter. Also, I am not sure how to judge what is an acceptable width of fluctuation NVE simulations? Does it correspond to random error (from the paper <https://www.biorxiv.org/content/biorxiv/early/2016/10/24/083055.full.pdf>Mark referenced)? Thank you! On Tue, Mar 20, 2018 at 1:39 AM, Mark Abraham <mark.j.abra...@gmail.com> wrote: > Hi, > > On Mon, Mar 19, 2018 at 10:08 PM Jo <jojo412...@gmail.com> wrote: > > > Hello, > > > > Thank you for your response. I have found a way to conserve energy of my > > box of SPC/E water by removing 'settles' from my topology file. > > > > Thus your water is no longer rigid. > > > > Previously, the total energy of the system was consistently increasing or > > decreasing by a significant amount resulting in no energy conservation. > > > > You're doing a numerical computation with finite-precision forces, > positions and velocities and a discrete fixed time step. You are guaranteed > to have both random and systematic errors compared to exact arithmetic. > Maybe you can find a parameter combination that will appear to lack > systematic errors, but that doesn't mean errors are absent, or the > combination is "better," as you have no doubt already read about in the > article I shared. > > > > I have defined 'bonds' and 'angles' with length and force constant for a > > flexible SPC/E water, with no 'settles' and no 'constraints'. For these > > parameters, the total energy converges to approximately the correct > energy > > and fluctuates withing 40 kJ/mol for a timestep of 1 fs, which is good. > > > > Why is it good? :-) > > > > However, I don't understand why or how removing the constraints allowed > for > > energy conservation. > > > It is implementing different equations of motion. If the only observable of > interest is the potential energy, and the multiple sources of numerical > inaccuracy are cancelling suitable for you, and you can afford to generate > enough samples at that short time step, then go right ahead :-) > > > > Does simply defining 'bonds' and 'angles' actually > > force the atoms to to constrained as a rigid molecule. > > > No, see the documentation about bonded interactions... > > > > If not, why does > > adding some sort of constraint make the system loose energy conservation? > > > > It is implementing different equations of motion, and they have different > numerical properties. Whether those properties are sufficient for any > purpose is an open question. I note that you ignored my earlier series of > questions on this topic, which were intended to help you learn ;-) > > How can I go about constraining the atoms as a rigid water molecule so that > > I don't loose energy conservation? > > > > As you can infer from Figure 3.5 in the reference manual, you can choose > the drift from the verlet buffer to add enough energy to offset the drift > from SETTLE. But obviously that would be a pair of compensating systematic > errors - so you won't learn anything from merely observing > that the systematic drift is zero. > > Mark > > > > Thank you in advance, > > > > Jo > > > > On Fri, Mar 16, 2018 at 4:51 PM, Mark Abraham <mark.j.abra...@gmail.com> > > wrote: > > > > > Hi, > > > > > > On Fri, Mar 16, 2018 at 7:42 PM Jo <jojo412...@gmail.com> wrote: > > > > > > > Thank you for your reply! > > > > > > > > I am attempting to conserve energy in an NVE run of 1000 SPCE water > - I > > > > have tried a number of different verlet-buffer-tolerances (0.001 to > > > 5e-5), > > > > sometimes the run output file suggests a specific verlet > > > buffer-tolerance. > > > > However I am still experiencing ~600 kJ/mol shift per ns. > > > > > > > > > You can't be getting that over that whole range. If you're doing NVE > with > > > tolerance 5e-5 then the drift is negative (from SETTLE, e.g. I just > > > observed -0.000069 kJ/mol/ps/atom on 1728 tip3p waters, GROMACS 2018 > > mixed > > > precision on a GPU, with SETTLE at 300K and 2fs timestep with PME and > > other > > > settings at defaults). > > > > > > I have tried > > > > double precision which does not seem to make a difference. I also > use > > a > > > > potential modifier (potential-shift) to ensure the potential reaches > 0. > > > > > > > > > You can never have an unshifted potential with the Verlet scheme. > > > > > > > > > > I > > > > have tried turning off the charges (to see if the ewald parameters > are > > > the > > > > source of the problem), but the same massive energy shift occurs - > > which > > > > leads me to believe the ewald parameters are not the source of the > > > problem. > > > > > > > > I am using 'settle' to constrain the water molecule. I am suspicious > > > that > > > > this could be the cause of the energy shift. Although, looking at > some > > of > > > > the previous posts on gromacs email forums, it appears that 'settle' > > has > > > > not been a problem for energy conservation previously. > > > > > > > > I have also tried different versions of Gromacs (5.1.4 and 2018), but > > the > > > > energy shift still occurs. > > > > > > > > Can you recommend any other parameters to change? Or any other > > > strategies > > > > to go about conserving energy for NVE? > > > > > > > > > First, given that you can't have zero drift in a numerical simulation > > with > > > a finite time step (and particularly not with constraints), have you > > > decided what level you regard as acceptable? For what application? Have > > you > > > considered the thoughts at https://www.biorxiv.org/node/23099.full? > > > > > > Mark > > > -- > > > Gromacs Users mailing list > > > > > > * Please search the archive at http://www.gromacs.org/ > > > Support/Mailing_Lists/GMX-Users_List before posting! > > > > > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > > > > > * For (un)subscribe requests visit > > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > > > send a mail to gmx-users-requ...@gromacs.org. > > > > > -- > > Gromacs Users mailing list > > > > * Please search the archive at > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > > posting! > > > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > > > * For (un)subscribe requests visit > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > > send a mail to gmx-users-requ...@gromacs.org. > > > -- > Gromacs Users mailing list > > * Please search the archive at http://www.gromacs.org/ > Support/Mailing_Lists/GMX-Users_List before posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.