Dear Justin You are right about the cut-off. The vdw energy spike disappeared after I increased the cut-off in the rerun. But I still don't understand why? Can the cutoff error build up to several thousand kj/mol for 100AA protein.
Again I would like to emphasize NOT to use a xtc file in MM/PBSA type calculation. The error is just too large. best, dawei On Thu, Aug 11, 2011 at 9:03 AM, Justin A. Lemkul <jalem...@vt.edu> wrote: > > > Da-Wei Li wrote: > >> Dear Justin >> >> 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. >> >> 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. >> >> > Right, forgot about that. I still think the cutoffs are a problem, though. > > -Justin > > best, >> >> dawei >> >> >> On Thu, Aug 11, 2011 at 8:47 AM, Justin A. Lemkul <jalem...@vt.edu<mailto: >> jalem...@vt.edu>> wrote: >> >> >> >> Da-Wei Li wrote: >> >> Dear Mark and others >> >> I did more tests and thought that it might come from numerical >> error. The reasons are >> >> 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. >> >> 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. >> >> 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. >> >> >> 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. >> >> >> >> To Mark only >> >> Thanks. Here it is my script for >> >> rerun: mdrun -v -s pbsa.tpr -rerun coor.xtc -e rerun >> 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) >> >> rerun .mdp file: >> >> ********************************__**************************** >> >> ; Run parameters >> integrator = md ; leap-frog integrator >> nsteps = 50000000 ; 100 ns >> dt = 0.002 ; 2 fs >> ; Output control >> nstxout = 500000 ; save coordinates every 1000 ps >> nstvout = 500000 ; save velocities every 1000 ps >> nstxtcout = 5000 ; xtc compressed trajectory output >> every 1 ps >> nstenergy = 5000 ; save energies every 1 ps >> nstlog = 5000 ; update log file every 1 ps >> xtc_grps = Protein ; save protein part only >> ; Bond parameters >> continuation = yes ; Restarting after NPT >> constraint_algorithm = lincs ; holonomic constraints >> constraints = hbonds ; all bonds (even heavy atom-H bonds) >> constrained >> lincs_iter = 1 ; accuracy of LINCS >> lincs_order = 4 ; also related to accuracy >> ; Neighborsearching >> ns_type = grid ; search neighboring grid cels >> nstlist = 10 ; 20 fs >> rlist = 0.8 ; short-range neighborlist cutoff (in nm) >> rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm) >> rvdw = 1.0 ; short-range van der Waals cutoff (in nm) >> ; Electrostatics >> coulombtype = cut-off ; Particle Mesh Ewald for long-range >> electrostatics >> pme_order = 4 ; cubic interpolation >> fourierspacing = 0.12 ; grid spacing for FFT >> ; Temperature coupling is on >> tcoupl = no ; modified Berendsen thermostat >> tc-grps = System ; two coupling groups - more accurate >> tau_t = 0.1 ; time constant, in ps >> ref_t = 300 ; reference temperature, one for each >> group, in K >> ; Pressure coupling is on >> pcoupl = no ; Pressure coupling on in NPT >> pcoupltype = isotropic ; uniform scaling of box vectors >> tau_p = 2.0 ; time constant, in ps >> ref_p = 1.0 ; reference pressure, in bar >> compressibility = 4.5e-5 ; isothermal compressibility of >> water, bar^-1 >> ; Periodic boundary conditions >> pbc = no ; 3-D PBC >> ; Dispersion correction >> ;DispCorr = EnerPres ; account for cut-off vdW scheme >> DispCorr = no >> ; Velocity generation >> gen_vel = no ; Velocity generation is off >> >> >> >> ; IMPLICIT SOLVENT ALGORITHM >> implicit_solvent = GBSA >> gb_algorithm = OBC >> nstgbradii = 1 >> rgbradii = 0.8 >> gb_epsilon_solvent = 80 >> gb_saltconc = 0 >> gb_obc_alpha = 1 >> gb_obc_beta = 0.8 >> gb_obc_gamma = 4.85 >> gb_dielectric_offset = 0.009 >> sa_algorithm = Ace-approximation >> sa_surface_tension = 2.25936 >> >> ********************************__****************************** >> **__***************************** >> >> Thanks all. >> >> >> 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. >> >> -Justin >> >> -- ==============================**__========== >> >> Justin A. Lemkul >> Ph.D. Candidate >> ICTAS Doctoral Scholar >> MILES-IGERT Trainee >> Department of Biochemistry >> Virginia Tech >> Blacksburg, VA >> jalemkul[at]vt.edu <http://vt.edu> | (540) 231-9080 >> <tel:%28540%29%20231-9080> >> >> >> http://www.bevanlab.biochem.__**vt.edu/Pages/Personal/justin<http://vt.edu/Pages/Personal/justin> >> >> <http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin> >> > >> >> ==============================**__========== >> >> -- gmx-users mailing list gmx-users@gromacs.org >> <mailto:gmx-users@gromacs.org> >> >> >> http://lists.gromacs.org/__**mailman/listinfo/gmx-users<http://lists.gromacs.org/__mailman/listinfo/gmx-users> >> >> <http://lists.gromacs.org/**mailman/listinfo/gmx-users<http://lists.gromacs.org/mailman/listinfo/gmx-users> >> > >> Please search the archive at >> >> http://www.gromacs.org/__**Support/Mailing_Lists/Search<http://www.gromacs.org/__Support/Mailing_Lists/Search> >> >> <http://www.gromacs.org/**Support/Mailing_Lists/Search<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 >> <mailto:gmx-users-request@**gromacs.org<gmx-users-requ...@gromacs.org> >> >. >> >> Can't post? Read >> http://www.gromacs.org/__**Support/Mailing_Lists<http://www.gromacs.org/__Support/Mailing_Lists> >> >> <http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists> >> > >> >> >> > -- > ==============================**========== > > 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<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<http://lists.gromacs.org/mailman/listinfo/gmx-users> > Please search the archive at http://www.gromacs.org/** > Support/Mailing_Lists/Search<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<http://www.gromacs.org/Support/Mailing_Lists> >
-- 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
