[gmx-developers] [Fwd: coulomb SR in 3.3cvs]

[Berk Hess]

Michel Cuendet wrote:

>
> Hi guys,
>
> Are there any news about that issue (see below my original message to 
> this list)?
>
> I've been trying to figure out what could in theory cause a relaxation 
> of the coulomb SR energy in 3.3.a lower than in 3.2.1. It could be 
> something with the constraints (I constrain all bonds involving 
> hydrogens) ?
>
> By minimizing with 321 and subsequently with 33a, I see the same kind 
> of energy difference as reported below. But there is no big apprarent 
> movement. RMSD's between the two conformations are of the order of 
> 0.0005 nm for the protein, and 0.0015 for the water, uniformly spread 
> over atom numbers except for a few molecules jumping the cell boundary...
>
> I would really like to know if 3.3.a is usable or not...
>
> Thanks to anyone commenting on that issue,
>
> Michel
>
>
> -------- Original Message --------
> Subject: coulomb SR in 3.3cvs
> Date: Mon, 11 Jul 2005 17:34:42 +0200
> From: Michel Cuendet <michel.cuendet at epfl.ch>
> Reply-To: michel.cuendet at epfl.ch
> To: gmx-users at gromacs.org
> References: <20050711100013.A17C62509A at xray.bmc.uu.se>
>
> Hi,
>
> I compiled the 3.3 CVS version of 27.05.2005 on itanium2 with the intel
> 8.1 compiler. (for those who don't know, the cvs version brings a
> terrific speedup on itaniums). Tests were performed with both versions
> in double precision, and were reproduced on an athlon platform.
>
> I compared the different energies obtained with 3.2.1 and 3.3, by taking
> the same trr (from a 3.2.1 equilibration, 85000 atoms, PME) and running
> it with both versions for about 12ps. From the tables below, one sees
> that there is a gap between Coulomb (SR) of 3450 KJ/mol, and between
> potential energies of  3280 KJ/mol.
>
> A close inspection of the energy files reveals that the version 3.3
> reproduces the same energy as 3.2.1 at the first step, but the system
> relaxes very quickly (0.5ps) to a lower energy with 3.3.
>
> Rerunning with 3.3 the trajectory produced by 3.2.1 surprisingly gives
> the correct energies.
>
> Any explaination for that ?
>
Maybe it can be explained by a change in the way Ekin and the temperature
are determined.

I changed the kinetic energy calculation in the CVS from being based on
the average of the n-1/2 and n+1/2 velocities to n-1/2 and n+1/2 Ekin.
With the old method the discretization error in Ekin (always causing
Ekin to be too small) was 4 times larger than with the new code.
The error causes the temperature to appear as too low, so the temperature
coupling scales it up too much, which effectively results in an ensemble
at a too high temperature.

The error in Ekin and T goes with dt^2, so making the time step smaller
should make the discrepancy between 3.2.1 and cvs smaller.

With water at 300 K in 3.2.1 the error in T is 1 K for dt=2 fs and 4 K 
for dt=4 fs.
In the CVS version the errors are a factor 4 smaller.

For your 85000 atom system this error is on the order of magnitude of
the error in Coulomb SR you have reported.

Berk.