[gmx-users] strange vdw energy from rerun with GBSA model.

[Da-Wei Li]

Hello

Please see my response below.

On Thu, Aug 11, 2011 at 10:23 AM, Mark Abraham <Mark.Abraham at anu.edu.au>wrote:

>  On 11/08/2011 10:22 PM, 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.
>
>
> That doesn't make sense if the underlying frames are the same in your .xtc
> and .trr files. However, if you're talking about different numbers of
> frames, then you probably have transient periodicity artefacts.
>


They are the same but different precision. I used a mdp file without PBC. I
believe this is reasonable because I can fully reproduce the bond length and
angle energy if I use trr file but cannot reproduce if I use .xtc file.


>
> 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.
>
>
> Sounds like you have a non-equilibrium trajectory (from your EM?), in which
> case small deviations can have large effects.
>


No. It is from a stable NPT trajectory with explicit water models. The
forces are around 10000 but are not larger.


>
>
>
> 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.
>
>
> Several tens of kJ/mol error in the non-bonded energy looks OK compared to
> the magnitude of the non-bonded energy. Anyway, if your purpose is to
> compare energies computed from the same configurations, then you should use
> the same configurations for computing the energies, not approximations to
> those configurations.
>
>

That is why I say precesion in xtc file is not enough for bond length, angle
potential and vdw potential.


>
>
> 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
>
>
As Justin suggested, increase cutoff in rerun can remove the vdw energy
spike.




>
> 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)
>
>
> Here you are generating a nojump trajectory file, but you are not doing
> your rerun on it. That could explain some symptoms - you might not have
> whole molecules.
>
>
I am very sorry for this typo. I do use the nojump.xtc in my rerun.



>  Mark
>
>
>
> 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.
>
> dawei
>
>
>
>
>
>
>
>
>
>
>
> --
> gmx-users mailing list    gmx-users at 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-request at gromacs.org.
> Can't post? Read 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/a94b2aad/attachment.html>