- Pauli potential
- Kinetic energy
- Non-Ewald electrostatics
- Ewald electrostatics
- Pairwise cutoff function
- Integration algorithm
- Thermostat

Pauli potential

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

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,
. 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

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

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

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:

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 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 , the pressure virial has the simple form
, where 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 and , 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:

The positions are stored in

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.

Thermostat

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 , which is modulated up and down to adjust the instantaneous temperature to match the desired temperature:

In a chain thermostat, there are several damping factors . acts as a damping factor for the entire system, acts as a damping factor for , for 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.