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