Ewald

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)