AMBER Archive (2003)

Subject: Re: Bug found in sander - incorrect conversion factor causes randomized velocities to be too cold

From: David A. Case (case_at_scripps.edu)
Date: Mon Jan 20 2003 - 01:26:16 CST

On Sat, Jan 18, 2003, John D. Chodera wrote:
>
> I hate to trouble you, but I was wondering if there was some way to get
> reasonable energy conservation for constant-energy (ntt=0) simulations
> without solvent when cutoffs are used in AMBER7. We're trying to use GBSA
> and want to compute free energy differences in the canonical ensemble, and
> I'm a bit hesitant to use the Berendsen thermostat because many texts note
> it doesn't sample from a proper canonical ensemble. We were therefore
> thinking of randomizing velocities every so often or using an
> Andersen-like thermostat. Unfortunately, energy conservation with
> reasonably finite cutoffs seems to be rather poor. Is there a
> force-switch function available in the nonbonded calculations? The switch
> eedmeth=3 seems to enable some sort of switch, but I'm afraid I couldn't
> really tell what was going on from a cursory reading o the PME code in
> which nonbonded calculations now seem to be embedded.
>

There is no force switch currently implemented. The eedmeth=3 is related
to PME, and is not active for GBSA simulations. Given the complexity of
the long-range behavior of the GB potential, it is not clear how best to
implement a smooth cutoff, which in any event, leads to the usual problems
with cutoffs (even if it would conserve total energy).

Since the potential falls off rapidly with distance, in many cases a cutoff
off of 15-20 Ang. will suffice; you can also play with nrespa (try up to
4), but we recognize the the GB method is not yet as efficient as we would
like.

The Berendsen thermostat has a worse reputation than it really deserves,
although for a small number of degrees of freedom (which you might have with
implicit solvent) it can be less effective than other methods. You might
want to look at:

%A T. Morishita
%T Fluctuation formulas in molecular-dynamics simulations with the weak
coupling heat bath
%J J. Chem. Phys.
%V 113
%P 2976
%D 2000

Vrand is indeed roughly like an Andersen thermostat (depending on details
of implementation. Note that, in its current form in Amber 7, the kinetic
energy on the step where it is invoked is not good if SHAKE is being used.
This has no effect on the trajectory (and hence on any averages over
coordinates or over potential energies), but it will influence what the system
thinks the average temperature is, and the average and fluctuations in
total energy.

..good luck....dac 

-- 

==================================================================
David A. Case                     |  e-mail:      case_at_scripps.edu
Dept. of Molecular Biology, TPC15 |  fax:          +1-858-784-8896
The Scripps Research Institute    |  phone:        +1-858-784-9768
10550 N. Torrey Pines Rd.         |  home page:                   
La Jolla CA 92037  USA            |    http://www.scripps.edu/case
==================================================================