[gmx-users] strange vdw energy from rerun with GBSA model.
Da-Wei Li
lidawei at gmail.com
Thu Aug 11 15:40:24 CEST 2011
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 <jalemkul at 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 <jalemkul at vt.edu<mailto:
>> jalemkul at 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 at gromacs.org
>> <mailto:gmx-users at 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-request at gromacs.org
>> <mailto:gmx-users-request@**gromacs.org<gmx-users-request at 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 at 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-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<http://www.gromacs.org/Support/Mailing_Lists>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110811/939d61bb/attachment.html>
More information about the gromacs.org_gmx-users
mailing list