next up previous contents index
Next: To add new features Up: eFF Programmers' Guide Previous: Recompiling the software   Contents   Index


To modify current features

Pauli potential

To modify the Pauli potential, change the inline function PauliElecElec() in eff_update, which computes the Pauli repulsion between two electrons. Be sure to update the analytic gradients as well as the energy. The function takes as input samespin, a flag which is 1 if the two spins are the same; rc, the distance between the two electrons; and re1 and re2, the sizes of the two electrons. The function outputs an energy, and forces frc, fre1, and re2 that are the negative gradients of energy with respect to distance and electron sizes. Arguments i and j are the indices of the electrons, but they are not used.

No other changes need to be made. As a note, the PauliElecElec() is called by UpdatePauliPeriodic() for Ewald-type calculations, and UpdateElecAndPauli() for all other calculations.

Kinetic energy

To modify the kinetic energy function, change the inline function UpdateKineticEnergy() in eff_update, which computes the kinetic energy of all the electrons.

The function does not take any arguments or return any values, but instead modifies directly the global data structure values elec[i].energy, the individual electron energy, and elec[i].fr, the individual electron radial force. If the kinetic energy depends on other coordinates, the forces or negative gradients relative to those coordinates need to be added too.

The function also calls AddSizeForceVirial() from eff_pressure to update the potential energy portion of the pressure virial, $ -r_{i} \cdot f_{i}$. Be sure to properly update this value. If the kinetic energy depends on other coordinates, the more general virial update function AddForceVirial() should be used instead.

Non-Ewald electrostatics

To modify the pairwise electrostatics functions, change any of three inline functions: ElecNucNuc(), for nucleus-nucleus repulsion; ElecNucElec(), for nucleus-electron attraction; and ElecElecElec(), for electron-electron repulsion. The arguments of these functions are similar to those in the Pauli potential described above. Beyond these functions, no other changes need to be made.

The current functions describe the interactions between Gaussian charge distributions, which are screened Coulomb potentials of the form:

$\displaystyle V = \frac{\mathrm{erf}(a x)}{x}

Functions such as ierfoverx1() compute the error function ratio and the first derivative. If one wishes to describe interactions involving other charge distributions, other auxiliary functions will be needed.

In our experience, the electrostatic energy has been the most difficult component of eFF to modify successfully. The Gaussian functions lack nuclear cusps and die off too quickly, but attempts to introduce these features by changing the electrostatic expressions have failed - there is a delicate cancellation of errors that is easily upset.

Ewald electrostatics

To modify the Ewald sum, change the real space sum RSpaceEnergy(), the reciprocal space sum KSpaceEnergy(), the self interaction SelfEnergy(), and the uniform background/charge interaction UniformChargeEnergy().

All of these functions take as input the arrays ewald_q, representing the charges, and ewald_alpha, representing the exponents of the Gaussian charges. They output an energy, and its gradients with respect to position and exponent, fx, fy, fz, and falpha.

If the goal is to use different kinds of charge distributions, new expressions for all of the above sums need to be derived. This includes an expression for the Fourier transform of the charge distribution, in order to compute the reciprocal space sum:

$\displaystyle E_{\mathrm{k-space}} = \frac{1}{2V} \int d\mathbf{k} \frac{4 \pi}{k^{2}} \left\vert\rho(\mathbf{k})\right\vert^{2}

In the current eFF, $ \rho(\mathbf{k})$ is a simple single Gaussian function, since the Fourier transform of a Gaussian is another Gaussian. The Fourier transform of other charge distributions may contain multiple terms, making the sum tedious to evaluate. In such cases, it may be useful to use other approaches to solving the Poisson equation, such as a Fast Fourier Transform on a reciprocal-space grid, or a multigrid or conjugate-gradient solver over a real-space grid. Such approaches have also been used to handle interactions more general than $ V=-1/r$, as in the Poisson-Boltzmann equation.

On the other hand, it may be difficult for grid-based algorithms to handle continuously varying electron sizes. The Ewald algorithm we have implemented, which handles charges of variable size, is slow partly because quantities such as $ e^{-k^{2}/4\alpha}$ are not precomputed, as they would be in a more traditional Ewald algorithm with fixed charge sizes. There may be a way to interpolate between precomputed values of this exponental function to make our algorithm more efficient.

As long as the interactions are all $ -1/r$, the pressure virial has the simple form $ \sum_{i} \mathbf{r_{i}} \cdot \mathbf{f_{i}} = U$, where $ U$ is the potential energy, and no further changes need to be made to EwaldVirial().

Pairwise cutoff function

To modify the pairwise cutoff function, change the inline functions cutoff() and dcutoff(), which are the cutoff function and its derivative respectively.

The function must have $ f(0)=1$ and $ f(1)=0$, and several of the higher order derivatives at the boundaries should be zero as well. The current cutoff function is a seventh-order spline; other lower-order splines are included in the code and commented out.

Integration algorithm

To modify the integration algorithm, change the functions UpdateDynamicsPositions() and UpdateDynamicsVelocities() in eff_dynamics. These two functions implement the two equations in the velocity Verlet algorithm:

$\displaystyle x_{i+1}$ $\displaystyle =$ $\displaystyle v_{i} \Delta t + \frac{1}{2 m} f_{i} \Delta t^{2}$  
$\displaystyle v_{i+1}$ $\displaystyle =$ $\displaystyle m \Delta t (f_{i} + f_{i+1})$  

The positions are stored in dyn_x, the velocities are stored in dyn_v, and the forces are stored in dyn_f. The program currently keeps one previous iteration of forces as well, dyn_prev_f.

To allocate additional arrays, modify AllocateDynamics(). Be sure also to duplicate these arrays for the functions SaveDynamics() and RevertDynamcs(), used by the adaptive step size procedure in eff_driver to roll back to a previous iteration when needed.


Currently two thermostats are implemented, Nosé-Hoover in ApplyNoseHooverNuclei(), which is called in Dynamics(), and the Andersen thermostat as a part of UpdateDynamicsVelocities().

The Nosé-Hoover thermostat could be improved by turning it into a chain thermostat. Currently, the temperature degrees of freedom are represented by a single damping factor $ \zeta$, which is modulated up and down to adjust the instantaneous temperature to match the desired temperature:

    $\displaystyle \Delta f_{i} = -\zeta m_{i} v_{i}$  
    $\displaystyle \Delta \zeta = \Delta t \left( \sum_{i} \frac{1}{2} m_{i} v_{i}^{2} - \frac{3}{2} N k_{B} T \right)$  

In a chain thermostat, there are several damping factors $ \zeta_{1} \ldots \zeta_{N}$. $ \zeta_{1}$ acts as a damping factor for the entire system, $ \zeta_{2}$ acts as a damping factor for $ \zeta_{1}$, $ \zeta_{3}$ for $ \zeta_{2}$ and so on.

To add such a thermostat, it would be best to make a separate function to update a set of thermostat global variables, and to call it at the same point the ApplyNoseHooverNuclei() function is called. The function would also apply the damping factor to the overall system, as it does now.

It may be useful to thermostat the electron degrees of freedom as well as the nuclear degrees of freedom.

next up previous contents index
Next: To add new features Up: eFF Programmers' Guide Previous: Recompiling the software   Contents   Index
Julius 2008-04-29