Greetings, I have been running simulations of 254 SPC/E water molecules using the OPLS force field. As Nilesh mentioned earlier today, I am using the NPT ensemble. I am setting the reference pressure at 1 bar, but when I do a sequence of minimization, equilibration, and dynamics steps, and then finally use g_energy to determine the average pressure, I find an average pressure of ~2 bar, rather than 1 bar.
Earlier today, Justin pointed out that this probably means that the system has not been equilibrated long enough, so I did another run with longer equilibration. However, I get a similar, puzzling result of ~2 bar. Finally, I tried changing the temperature and pressure coupling methods; this time, I get a seemingly more reasonable result for average pressure, but, as I will describe below and in my PDF file, still somewhat puzzling to me (I am new to the field of molecular dynamics). If you have time, could you please look at my PDF file at http://www.andrew.cmu.edu/user/adeyoung/may13.pdf that summarizes what I have tried? In case it is helpful, also, here is a text description of what I have tried: (i) T coupling = v-rescale; P coupling = parrinello-rahman. Reference pressure = 1 bar. 1 ns equilibration; 2 ns dynamics. Gromacs tells me that the average pressure is ~2.77 bar, but when I use g_energy to extract pressure as a function of time to an xvg list and then use software such as Mathematica or Matlab to compute the mean of the list, I find that the average pressure is ~1.24 bar. (ii) Use the same parameters as before, except equilibrate for 5 ns instead of 1 ns. Gromacs says that the average pressure is ~2.82 bar, but when I extract the pressure data and compute the mean of the list, I find that the average pressure is ~6.81 bar. (iii) Use the same parameters as in (ii), except use berendsen temperature coupling and berendsen pressure coupling. Gromacs says that the average pressure is ~1.00 bar, but when I extract the pressure data and compute the mean of the list, I find that the average pressure is ~3.19 bar. In my PDF file at http://www.andrew.cmu.edu/user/adeyoung/may13.pdf, I have plotted pressure vs time for the dynamics runs of (i), (ii), and (iii), and these plots show VERY large magnitude oscillations. If you have time, I have three questions about these results: (1) Does it seem reasonable that I obtain an average pressure closer to my reference pressure (1 bar) when I use berendsen temperature coupling and berendsen pressure coupling -- instead of v-rescale temperature coupling and parrinello-rahman pressure coupling? (2) When I use g_energy to extract pressure vs time, and then compute the mean of the list of pressures, why is my answer so different from Gromacs' "black box" calculation (that is, using g_energy to have Gromacs simply print the average pressure over the dynamics run) of the average pressure? For all of my runs, I am using a step size of 1 fs (=0.001 ps). However, when I extract the pressures (using the default settings), I get a pressure value for every 0.1 ps, not 0.001 ps. Could it be that the difference in average pressures is due to the fact that "under the hood" Gromacs is using the pressure data at every step (0.001 ps), instead of every 0.1 ps? But, even if this is true, I still wouldn't necessarily expect such a big disagreement between the two calculations. (3) When I plot pressure vs time (as I have done in my PDF file), why is there such a large magnitude of pressure fluctuation? For example, in my runs, the maximum pressure is on the order of 3000 bar, whereas the minimum pressure is on the order of -3000 bar. These pressures are unreasonably large in magnitude (despite the fact that the average pressure nevertheless turns out to be of the correct order of magnitude in the long run). Is this true? Also, is "negative pressure" unphysical? Or, does "negative pressure" correspond to "compression" and "positive pressure" corresponds to "expansion," or something like this? Thank you very much for your time. I truly appreciate it. Sincerely, Andrew DeYoung Carnegie Mellon University -- 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