Hi, Thank you. I tried again, this time doing 100 ps of NVT equlibration, followed by 500 ps of NVE. I have posted my results here: http://www.andrew.cmu.edu/user/adeyoung/july27/withnvt.pdf
But, they are not much better. For the NVT equilibration, I used v-rescale temperature coupling with tau_t = 2 ps. What kind of temperature coupling do you typically use or recommend? Thanks, Andrew Andrew DeYoung wrote: > Hi, > > I have been simulating 1000 SPC/E water molecules in the NVE ensemble (by > using tcoupl = no and pcoupl = no). In these simulations, I am considering > only Lennard-Jones interactions; I have "turned charges off" by setting the > partial charges on oxygen and hydrogen to zero in a local copy of spce.itp > that I include in my topology file. > > If you have time, I have posted a pdf file containing a summary of what I > have tried so far: http://www.andrew.cmu.edu/user/adeyoung/july27/july27.pdf > > At first, I tried using vdwtype = Cut-off. However, Gromacs gave me a note > message: > --- > NOTE: > You are using a cut-off for VdW interactions with NVE, for good energy > conservation use vdwtype = Shift (possibly with DispCorr) > --- > > I ignored this note message and ran the simulation anyway. Indeed, as > Gromacs warned me, the (total) energy conservation is very bad; the total > energy decreases seemingly significantly in what appears to be an almost > linear fashion. I have done energy minimization, equilibration, and > production dynamics, but for simplicity, I have drawn the equilibration and > production dynamics results in the same color in the pdf file that I have > posted. > > (As far as energy minimization, I did a series of two energy minimization > steps. First, I used integrator = l-bfgs with dt = 0.001 ps, nsteps = > 10000, emstep = 0.01 nm, and emtol = 10.0 kJ mol^-1 nm^-1. I got the > following results: > > Low-Memory BFGS Minimizer converged to Fmax < 10 in 0 steps > Potential Energy = -9.3022546e+02 > Maximum force = 5.5933684e-01 on atom 1123 > Norm of force = 1.9135195e-01 > > So, it seems that the energy of the configuration is quite low--the water > molecules are relatively far apart from one another (the simulation box is > quite large). > > Next, I used integrator = steep (passing the ending configuration from the > l-bfgs minimization to grompp using -t followed by the trajectory file) with > dt = 0.001 ps, nsteps = 10000, emstep = 0.01 nm, and emtol = 10.0 kJ mol^-1 > nm^-1. I got the following results: > > Steepest Descents converged to Fmax < 10 in 1 steps > Potential Energy = -1.2968502e+03 > Maximum force = 7.0817396e-02 on atom 820 > Norm of force = 2.6016856e-02 > > This configuration seems probably acceptable. So much for energy > minimization.) > > Secondly, I used vdwtype = Shift as the note message recommended. I tried > two different sets of values of rlist, rvdw_switch, rvdw, and rcoulomb. > But, in both cases, energy conservation is bad--almost as bad as when I used > vdwtype = Cut-off. > > If you have time, can you please look at the pdf file that I have posted > (http://www.andrew.cmu.edu/user/adeyoung/july27/july27.pdf)? Do you have > any suggestions of other vdw types or vdw parameters that I should try to > hopefully improve energy conservation? On page 1, I have displayed the > total energy results. On pages 2 and 3, I have pasted a typical .mdp file > that I used and also the results of other quantities from the simulations. > > One thing that I notice is, on the bottom of page 2, it appears that the > _potential_ energy conservation is good, whereas the kinetic energy > conservation is quite bad. Do you know, what sort of clue should this fact > give me as to what is going on with my system? > What seems most obvious to me is that your system is cooling down because your molecules are slowing down. Your gen_temp is set to 298 K, but your initial velocity distribution produces a temperature over 305 K, so you're not exploring the energy surface you think you are. Certainly in NVE temperature is not necessarily a conserved quantity (the total energy is), but you're not starting under the conditions you want. The temperature in all of these cases is only approaching the desired 298 K towards the end of the 500 ps timeframe you've shown. Try an NVT equilibration at 298 K first, then switch to NVE. The other thing to consider is whether or not this should even work. You've taken a condensed-phase water model that is largely held together by hydrogen bonding and transformed it into a semi-gas state and removed the most potent interactions that keep the molecules together. Force fields must balance LJ and electrostatic terms to reproduce some desired behavior. In this case, you're eliminating half of the equation (or more, depending on how the LJ and electrostatic terms are balanced to contribute to stability). I would suggest that if you have to hack this model under these conditions, start with NVT equilibration before moving to NVE. Even better would be to use a simpler system, like an ideal gas that only has LJ terms, and determine the proper protocol for obtaining a stable NVE ensemble. That way, you're not stabbing in the dark as to what the appropriate settings should be in the .mdp file. -Justin -- ======================================== Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ======================================== -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists