Ewald in Sander

Right. The basic deal with the Ewald method is that the conditionally convergent Coulombic sum is transformed into a sum of two rapidly converging sums: the direct sum (a sum in real space) and a reciprocal sum (a sum in reciprocal space). "beta" is an arbitrary constant in the expressions which determines the relative rate of convergence of the two sums. If you tighten up the convergence of the direct sum (i.e. make it shorter ranged by increasing beta) this leads to loosening up the convergence of the reciprocal sum.

The total energy (assuming a neutral system) is independent of the "beta" so we can choose any one. In the PME, the reciprocal sum in *very* fast due to the interpolation of the charge grid and the FFTs used in the evaluation of the sum. So, in order to simplify and speed the calculation of the direct sum, we choose "beta" so that the direct sum becomes negligible outside the cutoff. This way, we only need to calculate the terms of the direct sum over our interacting pairs (the pairlist) and can calculate the direct sum at the same time we are doing the standard lennard jones terms.

I see an input that specifies what is considered "vanishing". What is not clear is the algorithm where, given this tolerance, the beta is calculated. Are there analytical expressions for this or is it a trial and error method?

It is a trial and error method. Basically, we start by assuming beta = 1 (i.e. direct sum is extremely short ranged) then make sure the direct sum value at the cutoff is less than the tolerance.

           erfc( beta * cutoff )
           ---------------------  < tolerance
                 cutoff

Once we find a "beta" that is less than the tolerance (lets assume this is 1.0), we perform a binary search (tolerance down to 2**-50) to get a beta that gives the above direct sum basically equal to the tolerance at the cutoff...

       beta_lo = 0
       beta_hi = 1.0
       do i = 1, 50
          beta = (beta_lo + beta_hi)/2
          if ( erfc( beta*cutoff ) / cutoff .ge. tolerance ) then
             beta_lo = beta
          else
             beta_hi = beta
          endif
       enddo

Secondly, obviously, the beta is system dependent. In the Cheatham et. al. paper on DNA, RNA and Ubiquitin simulation (recent JACS communication), there was no mention of what was the value of beta for 9A cutoff is. I am very interested in knowing what it is.

In the above implementation, the beta is cutoff size dependent. In the JACS 117, 4193 (1995) communication, for the DNA simulation, I used a cutoff of 9.0A (for the lennard jones) and specified a tolerance of 0.000001. This leads to a beta of 0.34883. With a cubic interpolation, this leads to an estimate RMS force error of 2-3x10-4.

The last question I have is, isn't the beta parameter a function of distribution of charges in case where exact Ewald summation is not used, such as the PME with the non-bonded cutoff? If that is the case, does the PME algorithm readjust the beta parameter every so often to account for this?

The energy is independent of "beta" in Ewald and in the PME assuming you are accurately interpolating the reciprocal sum. I don't understand what you mean by "a function of distribution of charges".

In the exact Ewald case, for efficiency, the "beta" is box size dependent. Assymptotic convergence for the two sums comes when

   beta ~= sqrt(pi) / box_length

   beta = pi / cutoff   ( pi/9.0 = 0.349)

1) The B-spline version of PME will not conserve momentum, and the net force for the whole system will not be zero, but a random number in the order of the RMS error in the force calculation. The authors recommended to remove this artifact by removing the average net force from each atom at each step of simulation. I wonder, if one always needs to set NTCM=1 when he wants to use PME, or it is still not enough?

In the current version, the energies are interpolated with the B-spline fit and the forces are calculated as the analytic derivative of this B-spline expression. This leads to very accurate energies but, as you point out, leads to incomplete force conservation. Therefore the net force for the whole system will not be zero unless the average net force is removed at each step. This is what is done (silently) within the PME code within AMBER and has been done this way since 4.1... This artifact could likely be avoided by interpolating the forces rather than the energies, and 5.1 may provide this option in the sander code.

Note that this force artifact is different than the "flying ice cube" phenomenon which relates to growing center of mass motion (see Harvey et al., J Comp Chem 19, 726-740 (1998)). With energy conservation, center of mass motion will grow as a random walk which can lead to displacement of the coordinates over a long simulation. More catastrophic is when there is an energy drain in the system (with temperature coupling applied to bring the slowly dropping temperature back up to target); in this case, the COM motion can grow almost exponentially, such that by after ~1ns, ~50K of the motion is in COM translation. The major artifact of this problem (growth of COM translation) can be eliminated by periodically removing the COM translation during dynamics. This can be done in 5.0 through the NTCM/NSCM flags or in 4.1 through the application of appropriate patches. Although there has previously been arguments on this list on this point (particularly with respect to calculation of transport properties), I recommend that with the current sander code that COM motion (translation) should be removed in periodic boundary simulations with Ewald at every restart or every 50 ps or so. The error of removing this velocity is small when this is done (and definitely critical to prevent major artifactual behavior from the growth of the center of mass motion), and errors in the calculated transport properties are likely smaller than that error resulting from the constant coordinate rescaling in constant pressure calculations.

2) In the same article it is said that it is not possible to apply the dispersion approximation ( metioned in this article) to the old Amber force field due to the existance of 10-12 interactions. I want to make sure if this statement is still corrent in the current implementation of AMBER 5.0. What would happen if one use PME with parm91.dat?

First off, the dispersion PME approximation is not included in any released version of AMBER. Part of the reason for this is pragmatic in that with the AMBER (Lorentz-Bertholet) combining rules (for vdw radii), to perform the PME on the r**-6 would require seven separate PME evaluations for the dispersion terms. This makes it prohibitively expensive. Given that the vdw correction is small (and may only be necessary for crystal simulations) at reasonable cutoffs (say in the 9-14 angstrom range), the vdw are calculated as normal within the PME code (except that an ATOM based cutoff is applied rather than a residue based as in standard sander).

Now if you did happen to have the PME with PME dispersion corrections, it would not be wise to use this with the Weiner et al. force field or any with 10-12 interactions since no correction is applied to the r**-10 and the code may not be set up for proper operation in this case.

Since the dispersion PME is not included, and the PME only applied to the electrostatic interactions, PME can be applied with any of the standard force fields.

Tom Cheatham

    Cutoff, extra skin =  8.0 1.0 whereas largest inscribed sphere has radius 9.0

The reason the program is exiting is that the CUT + SKIN > SPHERE. "SPHERE" is the largest incribed sphere possible for your unit cell. "CUT" is the chosen cutoff and "SKIN" is the buffer for the pairlist that is defined in the source code (src/gibbs/ewald/unitcell.h) and set to 1.0 ansgtrom. In other words with gibbs and PME the MAXIMUM cutoff you can currently use is equal to the largest possible incribed sphere in that unit cell minus 1.0.

To get around this problem, decrease CUT slightly or increase the size of the unit cell slightly (build a bigger box). Remember too that when running constant pressure dynamics, the box size fluctuates, so allow for this.

Tom Cheatham

It defines both the border between the direct and reciprocal space summation for the electrostatic interactions, and the cutoff for the considered VDW interactions.

Is there a PME algorithm also for the VDW interactions?

Although an algorithm for performing PME on the VDW interactions has been described [in the smooth PME paper of Essman et al. J Chem Phys 103, 8577-8593 (1995)] this was never included in the released versions of the code. Part of the reason for this is that to do a PME treatment on the van der Waals (only the r**-6 part) requires 7 PME evaluations due to the use of the Lorentz-Bertholet combining rules in AMBER. Although in principle we might be able to perform some kind of multiple time step update to the long range vdw (since this energy/force may slowly vary) or somehow impose geometric combining rules (which will alter the force field), support for this hasn't been put into the standard code nor extensively investigated in the published literature. Another reason for not including this option is that the standardly applied cutoffs in the 9 angstrom range may be sufficent (although this point is arguable in crystal simulations). Tom Darden has investigated this issue in detail. [Now one reason we might want to do a long range vdw correction is that then we can put the direct space cutoff really small, say ~5-6 angstroms, and then get some speed up (although accuracy may suffer in the reciprocal space calculation).]

So, the cutoff specifies both the point where the vdw interactions are truncated and the point where the direct space pairlist is built to. In fact, the same pairlist is used for both the vdw and electrostatics (direct space) in AMBER currently.

Tom Cheatham

> I would like to use the PME method for my simulation. The parameters of
> my box are  :
> 65.30 73.07 62.27 90.0000000  90.0000000  90.0000000
>
>
> For PME, and in order to run efficiently, the size of the charge grid in
> each dimension should be a  "product of powers of 2, 3, and 5" and the
> size must be choice in order to have a grid spacing of ~1.0 angstrom .
> few choice for PME parameters:
>                     75.00  75.00 75.00  90.00  90.00  90.00
>                     60.00  75.00 60.00  90.00  90.00  90.00
>
> So the best Choice seem to be
>                                        75.00  75.00 75.00  90.00  90.00
> 90.00
>
> Because I am a new user and I not sure to understand all the influence
> of these parameter on the calculation could you please help my to make
> the good choice?

First of all, leave the box size alone. I would probably use 64, 72 and 64 for my charge grid dimensions. The density doesn't need to be exactly 1.0 . Close to it is good enough.

On the other hand, the code will run ok with any charge grid size, but be incredibly slow. It needs to be a multiple of powers of 2,3 and 5 due to the FFTs. These use a recursive method to speed the transformation of the grid,and the recursion depends on peeling off factors of 2,3 or 5 from the array size. So for example, if a grid dimension is divisible by 13, i.e. you specify 65 for the first dimension, the FFT routine will split 65 into 5 x 13. The factor of 5 is good and it will start work. However it then will encounter the 13, which it doesn't have a fast routine for, and go into slow default mode, to do the fourier transform of an array of size 65, which is quite slow.

On the other hand, 64 = 2 x 2 x 2 x 2 x 2 x 2, and 72 = 2 x 2 x 2 x 3 x 3 and the recursion can work all the way down, and so the whole array can be done fast.

Tom Darden

There are bound to be various opinions on this, however my opinion is that dual temperature coupling is not desired and in general unnecessary with PME since (*in proper usage*) PME conserves energy.

Dual temperature coupling, as I understand it, is applied to avoid the hot solvent/cold solute problem which is a direct manifestation of localized energy sources/sinks which lead to higher/lower temperatures in different parts of the system. The largest source of disparity in the temperature of the solute and solvent is likely due to group-based truncation and the phenomenon of waters crossing back and forth the cutoff. This effectively heats the water up such that there is a discernable difference in the effective temperature of the solute (which is cooler and therefore moves less) and the solvent (even though the entire system is effectively at 300K). Another source of disparity comes from the Berendsen pressure coupling which is dissipative (since whenever the pressure is below zero or above the target pressure-- which is common-- the system is relaxed by contracting or expanding the box), and which may lead to effectively different rates of energy drain from the solvent and solute degrees of freedom. [Other sources of incomplete energy conservation are larger time steps and SHAKE tolerances that are too high.]

The hack to fix this is to separately couple the solute and solvent to a temperature bath which hides the problem. In fact, this may not fully solve the problem since there still may be effective temperature gradients (perhaps comparing the water near the surface of the solute to the bulk solvent), although I have not tested this myself. A way around this is to separately couple each degree of freedom or apply Langevin dynamics (however these approaches are not currently implemented in AMBER) which remove localized temperature gradients.

Now, why not use dual temperature coupling with PME? My worry with dual temperature coupling is that since there are fewer degrees of freedom in the solute (typically) this leads to effectively tighter coupling and less temperature fluctuation in the solute (unless the coupling time constant is made sufficiently long) than is expected. The tighter coupling could prevent sampling in principle by limiting larger kinetic energy fluctuations. While this may not affect equilibrium properties, it may change the dynamics and/or time sacle of barrier crossing. However, again I have not fully investigated this.

With standard cutoff methods, applying the dual cutoff is definately preferred over the hot solvent/cold solute problem. With PME however, and with the major energy drains plugged, the disparity in the solute/solvent temperature does not appear in my experience in multi-nanosecond length simulation, even with constant pressure and Berendsen pressure coupling. Therefore, I do not apply dual temperature coupling with PME. However, this does not mean it could be an issue in longer simulations...

Tom Cheatham