Dear Justin<br><br>An implicit water simulaiton with this short cutoff is problematic but I think it is fine for rerun. I want to exactly repeat the original energies in the explicit water MD. <br><br>The manu say "neighbor list searching will be performed for every frame" with rerun option. So that I don't think this is the cause.<br>
<br>best,<br><br>dawei <br><br><div class="gmail_quote">On Thu, Aug 11, 2011 at 8:47 AM, Justin A. Lemkul <span dir="ltr"><<a href="mailto:jalemkul@vt.edu">jalemkul@vt.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<div><div></div><div class="h5"><br>
<br>
Da-Wei Li wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
Dear Mark and others<br>
<br>
I did more tests and thought that it might come from numerical error. The reasons are<br>
<br>
1. If I use .trr file instead of the low precision xtc file, things become better, ie, I get much less snapshots that has high energy.<br>
<br>
2. I supplied -pforce in my mdrun -rerun and found that the high vdw energy was usually caused by one pair of atoms, whose distance was just very near the clash zone, so that small error on the coordinates would cause large energy error. The force is always around 10000.<br>
<br>
3. Actually bond length and bond angle energies are also affected. I can fully reproduce these two energies if I use .trr file in my rerun but will get several tens of kj/mol error if I use .xtc file, for a protein of size of 100 AA.<br>
<br>
<br>
Now the question I still have is whether numerical error can be so large? The xtc file has a precision of 0.001 nm. Anyway, I will test more by using double precision Gromacs and define energy groups so that I can compare energy of protein directly between original MD and rerun.<br>
<br>
<br>
<br>
To Mark only<br>
<br>
Thanks. Here it is my script for<br>
<br>
rerun: mdrun -v -s pbsa.tpr -rerun coor.xtc -e rerun<br>
superpose: trjconv -s em.tpr -f coor.xtc -o nojump.xtc -pbc nojump (em.tpr is generated for energy minimization, protein is in the middle of the box)<br>
<br>
rerun .mdp file:<br>
<br>
******************************<u></u>****************************<br>
<br>
; Run parameters<br>
integrator = md ; leap-frog integrator<br>
nsteps = 50000000 ; 100 ns<br>
dt = 0.002 ; 2 fs<br>
; Output control<br>
nstxout = 500000 ; save coordinates every 1000 ps<br>
nstvout = 500000 ; save velocities every 1000 ps<br>
nstxtcout = 5000 ; xtc compressed trajectory output every 1 ps<br>
nstenergy = 5000 ; save energies every 1 ps<br>
nstlog = 5000 ; update log file every 1 ps<br>
xtc_grps = Protein ; save protein part only<br>
; Bond parameters<br>
continuation = yes ; Restarting after NPT<br>
constraint_algorithm = lincs ; holonomic constraints<br>
constraints = hbonds ; all bonds (even heavy atom-H bonds) constrained<br>
lincs_iter = 1 ; accuracy of LINCS<br>
lincs_order = 4 ; also related to accuracy<br>
; Neighborsearching<br>
ns_type = grid ; search neighboring grid cels<br>
nstlist = 10 ; 20 fs<br>
rlist = 0.8 ; short-range neighborlist cutoff (in nm)<br>
rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm)<br>
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)<br>
; Electrostatics<br>
coulombtype = cut-off ; Particle Mesh Ewald for long-range electrostatics<br>
pme_order = 4 ; cubic interpolation<br>
fourierspacing = 0.12 ; grid spacing for FFT<br>
; Temperature coupling is on<br>
tcoupl = no ; modified Berendsen thermostat<br>
tc-grps = System ; two coupling groups - more accurate<br>
tau_t = 0.1 ; time constant, in ps<br>
ref_t = 300 ; reference temperature, one for each group, in K<br>
; Pressure coupling is on<br>
pcoupl = no ; Pressure coupling on in NPT<br>
pcoupltype = isotropic ; uniform scaling of box vectors<br>
tau_p = 2.0 ; time constant, in ps<br>
ref_p = 1.0 ; reference pressure, in bar<br>
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1<br>
; Periodic boundary conditions<br>
pbc = no ; 3-D PBC<br>
; Dispersion correction<br>
;DispCorr = EnerPres ; account for cut-off vdW scheme<br>
DispCorr = no<br>
; Velocity generation<br>
gen_vel = no ; Velocity generation is off<br>
<br>
<br>
<br>
; IMPLICIT SOLVENT ALGORITHM<br>
implicit_solvent = GBSA<br>
gb_algorithm = OBC<br>
nstgbradii = 1<br>
rgbradii = 0.8<br>
gb_epsilon_solvent = 80<br>
gb_saltconc = 0<br>
gb_obc_alpha = 1<br>
gb_obc_beta = 0.8<br>
gb_obc_gamma = 4.85<br>
gb_dielectric_offset = 0.009<br>
sa_algorithm = Ace-approximation<br>
sa_surface_tension = 2.25936<br>
<br>
******************************<u></u>******************************<u></u>***************************<br>
<br>
Thanks all.<br>
<br>
</blockquote>
<br></div></div>
Using cutoffs this small may be the source of your problem. Proper implicit solvent calculations require longer cutoffs than would normally be used in explicit solvent MD. Try with longer (2.0 nm) or infinite cutoffs and a fixed neighbor list (nstlist = 0) and see if that smooths out the problem. What's likely happening now is that you've got interactions moving very quickly in and out of the very short cutoff, causing spikes in energy in between neighbor list updates.<br>
<br>
-Justin<br>
<br>
-- <br>
==============================<u></u>==========<br>
<br>
Justin A. Lemkul<br>
Ph.D. Candidate<br>
ICTAS Doctoral Scholar<br>
MILES-IGERT Trainee<br>
Department of Biochemistry<br>
Virginia Tech<br>
Blacksburg, VA<br>
jalemkul[at]<a href="http://vt.edu" target="_blank">vt.edu</a> | <a href="tel:%28540%29%20231-9080" value="+15402319080" target="_blank">(540) 231-9080</a><br>
<a href="http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin" target="_blank">http://www.bevanlab.biochem.<u></u>vt.edu/Pages/Personal/justin</a><br>
<br>
==============================<u></u>==========<div><div></div><div class="h5"><br>
-- <br>
gmx-users mailing list <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists/Search</a> before posting!<br>
Please don't post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<br>
Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists</a><br>
</div></div></blockquote></div><br>