next up previous contents index
Next: Index Up: eFF Programmers' Guide Previous: To access hidden features   Contents   Index

Subsections

Function prototypes and data structures


eff_access

This file defines accessor functions that should be used to get and set coordinates, energies, and forces of electrons and nuclei.

For example, these functions are called in eff_driver to print the energy of the system during dynamics runs, eff_bounds to wrap nuclear/electron positions to fit within the system bounds, eff_efield to apply electric field forces to the system, eff_properties to compute dipole moments and electron densities, eff_restrain to add restraint forces, and eff_dynamics and eff_minimize to perform dynamics and minimization.

There is some overhead associated with updating particle properties through function calls. In the minimization and dynamics routine, the overhead is reduced by copying particle coordinates and forces to a separate set of arrays, and performing multiple operations on the intermediate arrays before writing the positions and forces back to the main system. This preserves modularity at the cost of using additional memory. Such a policy makes it simple to use different minimization or dynamics routines, or to add extra degrees of freedom to the particles.

An exception is made in eff_update, where the underlying electron and nuclei data structures are accessed directly, bypassing the accessor functions.

The energies and forces are updated only after a call to UpdateEnergyForces() in eff_update. The nuclei indices range from 0 to NumNuclei() - 1, and the electron indices range from 0 to NumElectrons() - 1.


double GetTotalPE()

Returns the total potential energy of the system.


double GetNucleusPE(int i)

Returns the potential energy associated with nucleus i. Pairwise energies between two particles are split evenly between the two particles.


double GetElectronPE(int i)

Returns the potential energy associated with electron i. Pairwise energies between two particles are split evenly between the two particles.


void GetNuclearForce(int i, double *fx, double *fy, double *fz)

Gets the force (fx, fy, fz) acting on nucleus i.


void GetElectronForce(int i, double *fx, double *fy, double *fz, double *fr)

Gets the force (fx, fy, fz, fr ) acting on electron i.


void GetNuclearPosition(int i, double *x, double *y, double *z)

Gets the position (x, y, z) of nucleus i.


void SetNuclearPosition(int i, double x, double y, double z)

Sets the position of nucleus i to (x, y, z).


void GetElectronPosition(int i, double *x, double *y, double *z, double *r)

Gets the position (x, y, z) and size r of electron i.


void SetElectronPosition(int i, double x, double y, double z, double r)

Sets the position of electron i to x, y, z and its size to r.


double GetNuclearCharge(int i)

Returns the charge of nucleus i.


int GetElectronSpin(int i)

Returns the spin of electron i.


void SetElectronSpin(int i, int spin)

Sets the spin of electron i to spin.


int NumNuclei()

Returns the number of nuclei.


int NumElectrons()

Returns the number of electrons.


void AddNuclearForce(int i, double fx, double fy, double fz)

Adds a force (fx, fy, fz) to nucleus i.


void AddElectronForce(int i, double fx, double fy, double fz, double fr)

Adds a force (fx, fy, fz, fr) to nucleus i.


void AddNucleusPE(int i, double energy)

Adds energy to the potential energy of nucleus i.


void AddElectronPE(int i, double energy)

Adds energy to the potential energy of electron i.


void SetNucFX(int i, double val)

Sets the x-component of the force on nucleus i to val. The functions AddNuclearForce and AddElectronForce are preferred; this function is used only in eff_freeze to null certain degrees of freedom.


void SetNucFY(int i, double val)

Sets the y-component of the force on nucleus i to val.


void SetNucFZ(int i, double val)

Sets the z-component of the force on nucleus i to val.


void SetElecFX(int i, double val)

Sets the x-component of the force on electron i to val.


void SetElecFY(int i, double val)

Sets the y-component of the force on electron i to val.


void SetElecFZ(int i, double val)

Sets the z-component of the force on electron i to val.


void SetElecFR(int i, double val)

Sets the r-component of the force on electron i to val.


eff_bounds

This file contains functions that help to enforce rectangular periodic boundary conditions. The main function is ApplyPositionBounds, which wraps particles that have strayed outside the boundaries back to the interior of the unit cell.

Minimum image convention functions should properly be in this file, but are instead included as inline functions in eff_update for performance reasons. In eff_update, the functions BoundDX, BoundDY, and BoundDZ modify interatomic position vectors with the minimum image convention, so that the distance between two particles in a unit cell is always taken to be the shortest one subject to periodic boundaries.


void InitializeBounds()

Precomputes the box lengths and volume, and determines whether the individual dimensions of the box are periodic or not, based on the parameters params.x_bound[], params.y_bound[], params.z_bound[], and params.periodic.


int isWrapX()

Returns whether the X dimension is periodic or not. Used internally, and in eff_update to determine whether a minimum image correction needs to be applied.


int isWrapY()

Returns whether the Y dimension is periodic or not.


int isWrapZ()

Returns whether the Z dimension is periodic or not.


double GetVolume()

Returns the volume of the rectangular box.


void ApplyPositionBounds()

Enforces periodic boundaries by mapping particles that have strayed outside the unit cell back inside the unit cell. Used in eff_minimize and eff_dynamics.


eff_constants

Physical constants and conversion factors.


eff_create

In this version of eFF, the number of particles stays fixed throughout the simulation. The nuclei and electrons are created once at the start of simulation through two steps. First the Initialize function is called to dynamically allocate space for the particles. Then the functions AddNucleus and AddElectron are used to add the particles to the simulation space.

In this implementation, all of the functions in eff_create are called in eff_input, when the input .cfg file is read. The program counts the number of nuclei and electrons before adding them one at a time, so that exactly the amount of space needed is allocated.

It may be useful to add the ability to create and destroy particles during the simulation. Currently simple arrays hold the nuclear and electron coordinates. In a more advanced implementation, a structure such as a linked list could be used instead.


void Initialize(int maxnuclei, int maxelectrons)

Allocate enough space for maxnuclei nuclei and maxelectrons electrons.


int AddNucleus(double x, double y, double z, double q)

Add a nucleus at position (x, y, z) and charge q. If the number of nuclei exceeds the allocated space, the program will stop with an error.


int AddElectron(double x, double y, double z, int spin, double re)

Adds an electron at position (x, y, z) with spin spin and size re.


eff_driver

This file runs the main loop of the program. First, the main procedure allocates memory for the elements contained in the input file: for the nuclei and electrons via eff_create, for user-set individual masses and velocities via eff_initial, and for constraints and restraints via eff_freeze and eff_restrain. The program peeks into the input file to see how much space needs to be allocated for each section.

Then the program sets up the parameters via eff_param and loads in their default values.

Next, the program reads in values from the input file via eff_input, sets up the output file via eff_output, and prints out any information it can prior to actually running the simulation.

With the parameters read in, the program initializes the boundaries via eff_bounds, the electric field via eff_efield, and the Ewald solver via eff_ewald if periodic = true has been specified. When setting up Ewald, the program also may auto-set the cutoff parameters ewald_r_cutoff and ewald_k_cutoff to achieve the precision specified in ewald_log_precision.

Then the program branches off to either a single point calculation, using RunSinglePoint(); a minimization, using RunMinimize(); or a dynamics run using RunDynamics().

Finally, the main loop prints out elapsed time, updates the restart file as needed, and reclaims previously allocated memory.


void RunSinglePoint()

Runs a single point calculation by calling eff_update to compute an energy, and reporting it back via eff_output.


void RunMinimize()

In the minimization routine, the main loop is delegated to eff_minimize which in turn is driven through eff_shanno. The driver's role is limited to providing a callback function that outputs energies, positions, forces, etc. via eff_output every several iterations.


void MinimizeCallback(int num_iterations, int num_fevals, double f_val, double g_sq)

Callback function for minimizer, provided to eff_minimize, which updates it every print_every iterations with the progress of the minimization: number of iterations num_iterations, number of function evaluations num_fevals, the value of the function f_val, and the gradient squared g_sq. The function in turn prints an update out to the output .eff file, and goes through eff_output to prints out position, energies, forces, etc.


void RunDynamics()

In the dynamics routine, the main loop is actually present in eff_driver, and it takes care of initializing masses and velocities; outputting energies, positions, and forces via eff_output regularly; and using an adaptive step size if energy is not conserved. The primary calls are through eff_dynamics, which is in charge of allocating memory for dynamics, taking single time steps, and implementing a thermostat.


void AnnealSpins(double start_temperature, double end_temperature, int numiter)

This routine is commented out. When called, it randomly swaps opposite spins in a Metropolis Monte Carlo fashion while keeping all particle positions fixed. The temperature is varied from start_temperature to end_temperature over numiter iterations. We used this routine to anneal the spin configuration of lithium solid to a lower potential energy state.


eff_dynamics

This file implements a velocity Verlet algorithm to propagate nuclear and electron dynamics. The function Dynamics() performs a single iteration of dynamics using a timestep set by SetTimeStep(). The propagation of positions and velocities occurs in the functions UpdateDynamicsPositions() and UpdateDynamicsVelocities(). Two thermostats are provided as well, Nose-Hoover in the ApplyNoseHooverNuclei() routine, and Andersen as a modification of the UpdateDynamicsVelocities() function.

In this eFF implementation, particle velocities and masses are stored separately as arrays in eff_dynamics, and are not part of the global data structure that stores nuclear and electron coordinates. A separate set of accessor functions provides access to masses and velocities, as well as temperatures, kinetic energy, and the time step. There is also a function UpdateDynamicsVirial() that computes the kinetic energy virial contribution to pressure, as it is velocity dependent.

To improve modularity, all the positions and forces in the system are copied to one-dimensional arrays dyn_x and dyn_f, and the integrator acts on those arrays. This additional translation layer makes it easy to add new degrees of freedom, change the underlying data structure for representing nuclei and electrons, or to program different integrators. A similar approach is used in eff_minimize.


Data structures

The arrays dyn_x, dyn_v, dyn_f, and dyn_prev_f store positions, velocities, forces, and previous forces respectively as one-dimensional arrays of length $ 3 N_{\mathrm{nuc}} + 4 N_{\mathrm{elec}}$, where $ N_{\mathrm{nuc}}$ is the number of nuclei and $ N_{\mathrm{elec}}$ is the number of electrons.

The degrees of freedom are stored as follows:

  nuc(1)  -- fx, fy, fz
  nuc(2)  -- fx, fy, fz
  ...
  elec(1) -- fx, fy, fz, fr
  elec(2) -- fx, fy, fz, fr
  ...
so that the nuclei occupy indices from $ 0 \ldots N_{\mathrm{nuc}} - 1$ and the electrons occupy indices from $ N_{\mathrm{nuc}} \ldots N_{\mathrm{nuc}} + N_{\mathrm{elec}} - 1$.

The functions GetDynamicsPositions() and GetDynamicsForces() translate positions and forces from the eFF system into the dynamics arrays, and the function SetDynamicsPositions() does the reverse operation.

The adaptive step size routine uses a spare set of arrays save_dyn_x, save_dyn_v, save_dyn_f, and save_dyn_prev_f to ``roll back'' an iteration when energy is not properly conserved. It does this using the functions SaveDynamics() and RevertDynamics().


void SetDynamicsLogfile(FILE *fp)

Sets the output logfile file handle to fp.


void AllocateDynamics()

Allocates space for the dynamics arrays discussed in the Data Structures section. Calls NumNuclei() and NumElectrons() via eff_access to determine how big the arrays should be.


void SaveDynamics()

Save the current dynamics arrays to a set of backup arrays (starting with save_), for possible restoration later. Used by the adaptive step size procedure in eff_driver.


void RevertDynamics()

Restores the current state to the dynamics state saved in SaveDynamics().


void InitializeDynamics()

Initializes the positions, velocities, masses, and forces for time zero. Must be called before Dynamics() can be valid.


void InitializeRandomVelocities(double temp)

Initializes the velocities of the nuclei to the given temperature temp using a Maxwell-Boltmann distribution. Random numbers are from eff_util. Random velocities can be applied to the electrons as well by uncommenting out a portion of the code. Calls NullCMMotion() to remove net translational drift.


void NullCMMotion()

Nulls out the net center-of-mass motion to prevent translational drift.


void Dynamics(enum ThermostatType s_thermostat, double temperature)

Performs a single time step of dynamics using the functions UpdateDynamicsPositions() and UpdateDynamicsVelocities(). Calls eff_update to get energies and forces, applies the Nose-Hoover or Andersen thermostat as needed, applies constraints and restraints and external electric fields, and updates the kinetic contribution to the pressure virial. Also moves positions and forces between the real system and the internal dynamics arrays.


void UpdateDynamicsPositions()

Propagates the particle positions with the velocity Verlet algorithm:

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

This routine is called by Dynamics() and should not be called separately.


void UpdateDynamicsVelocities(enum ThermostatType thermostat, double temperature)

Propagates the particle velocities with the velocity Verlet algorithm:

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

Also applies the Andersen thermostat if it is selected by the variable thermostat, which randomly replaces the velocity of the particle with one chosen from a Maxwell-Boltzmann distribution of velocities at the given temperature. This routine is called by Dynamics() and should not be called separately.


void ApplyNoseHoover(double temperature)

Applies the Nose-Hoover thermostat at the given 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)$  

Currently the function ApplyNoseHooverNuclei() is used instead.


void ApplyNoseHooverNuclei(double temperature)

Applies the Nose-Hoover thermostat to the nuclei only. Currently this is the version called by Dynamics(). This routine should not be called separately.


void GetDynamicsForces()

Copies nuclear and electron forces from the eFF system to the dynamics array dyn_f. Uses eff_access.


void GetDynamicsPositions()

Copies nuclear and electron positions from the eFF system to the dynamics array dyn_x.


void SetDynamicsPositions()

Copies nuclear and electron positions from the dynamics array dyn_x to the eFF system.


void UpdateDynamicsVirial()

Updates the pressure virial with the kinetic energy of the nuclei and electrons, placing separate contributions into rigid versus flexible pressures.


void SetTimeStep(double dt)

Sets the time step of the Verlet integrator to dt.


double GetTemperature()

Computes the temperature as $ \sum_{\mathrm{all}} \frac{1}{2} m v^{2} / \frac{3}{2} k_{B} N_{\mathrm{nuc}}$, where we sum over all degrees of freedom, but assume that the temperature is $ \frac{3}{2} k_{B}$ T per nucleus.


void SetNuclearVelocity(int i, double vx, double vy, double vz)

Sets the velocity of nucleus i to (vx, vy, vz). Used in eff_driver to set individual initial velocities.


void SetElectronVelocity(int i, double vx, double vy, double vz, double vr)

Sets the velocity of electron i to (vx, vy, vz, vr).


void GetNuclearVelocity(int i, double *vx, double *vy, double *vz)

Gets the velocity of nucleus i in (vx, vy, vz).


void GetElectronVelocity(int i, double *vx, double *vy, double *vz, double *vr)

Gets the velocity of electron i in (vx, vy, vz, vr).


void SetNuclearMass(int i, double mx, double my, double mz)

Sets the mass of nucleus i to (mx, my, mz). Used in eff_driver to set individual masses.


void SetElectronMass(int i, double mx, double my, double mz, double mr)

Sets the mass of electron i to (mx, my, mz, mr).


void GetNuclearMass(int i, double *mx, double *my, double *mz)

Gets the mass of nucleus i in (mx, my, mz). Not yet used.


void GetElectronMass(int i, double *mx, double *my, double *mz, double *mr)

Gets the mass of electron i in (mx, my, mz, mr).


double GetTotalKE()

Returns the total kinetic energy of the system $ \sum_{\mathrm{all}} \frac{1}{2} m v^{2}$.


double GetNucleusKE(int i)

Returns the kinetic energy of nucleus i:


double GetElectronKE(int i)

Returns the kinetic energy of electron i:

$\displaystyle T = \frac{1}{2}\left(m_{x} v_{x}^{2} + m_{y} v_{y}^{2} + m_{z} v_{z}^{2} + m_{r} v_{r}^{2}\right)
$


double GetElectronTranslationKE(int i)

Returns the translational kinetic energy of electron i:

$\displaystyle T = \frac{1}{2}\left(m_{x} v_{x}^{2} + m_{y} v_{y}^{2} + m_{z} v_{z}^{2}\right)
$


eff_efield

This file applies an electric field to the system that is constant in space, but may vary with time, either as a continuous sine wave or as a wave packet.


void SetupExternalField(double Ex0, double Ey0, double Ez0, double Efreq, double Epacket_duration)

Initializes the electric field with an amplitude (Ex0, Ey0, Ez0) in volts-per-centimeter, a frequency Efreq in Hertz, and a wave packet duration Epacket_duration in femtoseconds. The quantities are converted into the proper internal units using the constants V_PER_CM and T0_IN_FS.


void GetExternalField(double t, double *Ex, double *Ey, double *Ez)

This function outputs the electric field vector (Ex, Ey, Ez) at time t.

If Efreq and Epacket_duration are both zero, the electric field is constant over time. If the frequency is nonzero, but the wave packet duration still zero, it is a sine wave; otherwise if both the frequency and the wave packet duration are nonzero it is a wave packet.

This function is called by ApplyExternalField and is only used inside eff_efield.


void ApplyExternalField(double t)

Applies a uniform electric field to the system, taking the field from GetExternalField.


eff_eigen

This file computes the eigenvalues and eigenvectors of a matrix using Jacobi rotations. The routine was translated from an algorithm by H. Rutishauser in Numerische Mathematik 9, 1-10 (1966).


void jacobi(double *matrix, double eigenvalue[], double *eigenvector, int n)

Computes the eigenvalues and eigenvectors of a dimension n matrix. Used in eff_restrain to find least-squares lines and planes.


eff_erf

This file computes the function $ \mathrm{erf}(x)/x$ and its first and second derivatives if desired. These functions are used in eff_update and eff_ewald to compute the electrostatic interactions between Gaussian distributions of charge.

The error function is defined as

$\displaystyle \mathrm{erf}(x) = \frac{2}{\pi} \int_{0}^{x} e^{-t^{2}}
$

We compute the function using a Chebyshev polynomial expansion and bilinear mapping from J. Schonfelder, Math. Comp. 32(144):1232-40 (1978), and a routine adapted from the public domain SLATEC/FNLIB routine csevl.f by W. Fullerton, which uses the Clenshaw method to perform the expansion.

The user may select how many terms of the expansion to use using the constants ERF_TERMS1 (for $ x < 2$), ERF_TERMS2 (for $ x \geq 2$), and DERF_TERMS (for the first derivative, $ x < 2$).

The functions are replicated in eff_update as inline functions for performance reasons.


double erfoverx(double x)

Returns $ \mathrm{erf}(x)/x$ for any $ x>0$. Calls the Chebyshev expansion poly01() when $ x < 2$ and the expansion poly02() when $ x \geq 2$.


double erfoverx1(double x, double *df)

Returns $ \mathrm{erf}(x)/x$ and its first derivative for any $ x>0$. It requires less work to calculation the function and its derivative together than to compute them separately. Versus erfoverx(), the routine makes an added call to the built-in exponential function exp() when $ x>2$ and to the Chebyshev expasion poly1() when $ x \geq 2$.


double erfoverx2(double x, double *df, double *df2)

Returns $ \mathrm{erf}(x)/x$ and its first and second derivatives for any $ x>0$. Versus erfoverx1(), the routine makes an additional call to the polynomial poly2() when $ x < 2$. Currently not used, but could be pressed into service if analytic second derivatives were desired.


double poly01(double x)

Auxiliary polynomial $ P(x)$ (defined in the paper) for $ x < 2$.


double poly02(double x)

Auxiliary polynomial $ P(x)$ for $ x \geq 2$.


double poly1(double x)

Auxiliary polynomial $ P'(x)$ for $ x < 2$.


double poly2(double x)

Auxiliary polynomial $ P''(x)$ for $ x < 2$.


eff_ewald

This file computes the electrostatic interaction energy between a collection of periodic Gaussian charge densities:

$\displaystyle \rho_{i}(\mathbf{r}) = q_{i} \sum_{\mathbf{R_{L}}} e^{-\alpha_{i} \left(\mathbf{r} - \mathbf{r_{i}} - \mathbf{R_{L}}\right)^{2}}
$

where $ \mathbf{R_{L}}$ is a set of periodic lattice vectors.

The standard Ewald sum procedure operates on point charges. We have implemented a generalized procedure that works on charges of finite width. The key change we make is to define a threshold exponent $ \alpha_{\mathrm{max}}$. When $ \alpha < \alpha_{\mathrm{max}}$, i.e. the charges are diffuse, we sum the electrostatic energy entirely in reciprocal space:

$\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}
$

The periodicity of $ \rho$ in real space reduces this integral to a lattice sum which converges quickly in $ k$ when $ \alpha_{\mathrm{max}}$ is sufficiently small.

For $ \alpha \geq \alpha_{\mathrm{max}}$, we divide the charges into two portions:

    $\displaystyle \rho = \rho_{1} + \rho_{2}$  
    $\displaystyle \rho_{1} = N e^{-\alpha r^{2}} - N_{\mathrm{max}} e^{-\alpha_{\mathrm{max}} r^{2}}$  
    $\displaystyle \rho_{2} = N_{\mathrm{max}} e^{-\alpha_{\mathrm{max}} r^{2}}$  

We sum $ \rho_{1}$ over a real space lattice and sum $ \rho_{2}$ over a reciprocal space lattice. Since $ \rho_{1}$ is neutral for large $ r$ and localized to a small region of real space, the real space sum converges quickly in $ r$. Since $ \rho_{2}$ is diffuse for a small enough $ \alpha_{\mathrm{max}}$, the reciprocal space sum converges quickly in $ k$.

The routines EwaldRCutoff() and EwaldKCutoff() choose real and reciprocal space cutoffs to achieve a desired precision goal, guided by the known asymptotic behavior of the Gaussian distributions. We have chosen a default $ \alpha_{\mathrm{max}}$ that optimizes the overall CPU time, assuming that these automatic cutoffs are in use.

The final Ewald energy is

$\displaystyle E = E_{\mathrm{k-space}} + E_{\mathrm{r-space}} + E_{\mathrm{self}} + E_{\mathrm{uniform}}
$

The self energy is the energy of a Gaussian charge with itself - this quantity must be subtracted out, since in reality, charges do not repel themselves. The uniform energy represents the interaction of each Gaussian charge with the uniform neutralizing background charge.


Data structures

The program converts the eFF nuclei and electrons into Gaussian charge densities. Nuclei are represented as very sharply peaked Gaussians whose width is set by ewald_r_nuc. Hence all the arrays have a dimension equal to the total number of nuclei and electrons.

The arrays ewald_x, ewald_y, ewald_z contain the charge density positions; ewald_alpha contains the exponents; and ewald_q the charges.

The force arrays ewald_fx, ewald_fy, ewald_fz, ewald_falpha, store the negative derivatives of the energy relative to x, y, z, and alpha. The energy array ewald_energy stores the energy per particle.

The arrays ewald_S_real and ewald_S_imag represent the charge density in reciprocal space, which is a complex number. The arrays ewald_S_cap_real and ewald_S_cap_imag, are analogous, but for $ \alpha = \alpha_{\mathrm{max}}$.

The arrays lattice_kx, lattice_ky, lattice_kz, and lattice_k2 prestore reciprocal lattice points and the arrays lattice_rx, lattice_ry, and lattice_rz prestore real space lattice points. Those arrays are filled once at the start of the program.


void AllocateEwald()

Allocates space for the Ewald arrays discussed above.


void InitializeEwald(double Lx, double Ly, double Lz, double max_alpha, double max_r, double max_k, double nuc_alpha)

Sets the rectangular periodic box to have size (Lx, Ly, Lz); the variable $ \alpha_{\mathrm{max}}$ to be max_alpha, the real space and reciprocal space cutoffs to be max_r and max_k, and the exponent of the nuclei to be nuc_alpha.


void SetEwaldCharges(double nuc_alpha)

Converts the eFF nuclei and electrons into Gaussian charge distributions suitable for use in eff_ewald.


double TotalCharge()

Returns the total charge.


double TotalVolume()

Returns the total volume of the unit cell.


void AddEwaldEnergyForce(double *fx, double *fy, double *fz, double *falpha, double *energy)

Converts the energies and forces computed by the Ewald procedure into energies and forces that act on eFF nuclei and electrons, and adds them directly to the global data structure. It would be more appropriate to use eff_access instead. The Ewald arrays taken as input are ewald_fx, ewald_fy, ewald_fz, ewald_falpha, and ewald_energy.

Called several times in UpdateEwaldEnergy(). Note that the arrays such as ewald_energy are completely overwritten between calls. If this is not done, it is important to zero them, so that energies are not added twice to the eFF system.


void EwaldVirial(double *ewald_energy, double *ewald_falpha)

Computes the Ewald contribution to the pressure virial. It can be shown that for particles interacting through an $ -1/r$ potential, the virial is just the negative of the potential energy. This holds true if the charge density is uniformly scaled. If we assume the electron size stays rigid, we must make an additional correction that depends on ewald_falpha.


void UpdateEwaldEnergy()

Main Ewald procedure. Calls SetEwaldCharges() to get Gaussian charges from the eFF system, computes the sum of self, uniform, real space, and reciprocal space energies, and then calls AddEwaldEnergyForces() to apply the calculated energies and forces back to the eFF system. Also calls EwaldVirial() to update the pressure virial.


void InitializeKSpace()

Allocates and precomputes the points of the reciprocal space lattice, filling the arrays lattice_kx, lattice_ky, lattice_kz, and lattice_k2 (magnitude). Uses the k-space cutoff parameter ewald_max_k. The point (0, 0, 0) is excluded. This routine is run only once, inside InitializeEwald().


void InitializeRSpace()

Allocates and precomputes the points of the real space lattice, filling the arrays lattice_rx, lattice_ry, and lattice_rz. Uses the r-space cutoff parameter ewald_max_r. The point (0, 0, 0) is excluded. This routine is run only once, inside InitializeEwald().


void KSpaceEnergy(double *fx, double *fy, double *fz, double *falpha, double *energy)

Computes the reciprocal space energy:

$\displaystyle E_{\mathrm{k-space}} = \frac{2 \pi}{V} \sum_{k} \rho(\mathbf{k}) \rho_{\mathrm{max}}(\mathbf{k})
$

where
$\displaystyle \rho(\mathbf{k})$ $\displaystyle =$ $\displaystyle \sum_{i} q_{i} e^{-k^{2}/4 \alpha_{i}} e^{-i \mathbf{k} \cdot \mathbf{r_{i}}}$  
$\displaystyle \rho_{\mathrm{max}}(\mathbf{k})$ $\displaystyle =$ \begin{displaymath}\left\{
\begin{array}{ll}
\rho(\mathbf{k})\vert _{\alpha_{i}=...
...x}}\\
\rho(\mathbf{k}) & \mathrm{otherwise}
\end{array}\right.\end{displaymath}  

Places the energy in the variable energy. Also computes the negative gradients and places them in the variables fx, fy, fz, and falpha.


void RSpaceEnergy(double *fx, double *fy, double *fz, double *falpha, double *energy)

Computes the real space energy:

$\displaystyle E_{\mathrm{r-space}} = \sum_{\mathbf{R_{L}}} \sum_{i \neq j} E_{ij} \left(\mathbf{r_{ij}} - \mathbf{R_{L}}\right)
$

where
$\displaystyle E_{ij} =$   $\displaystyle \frac{1}{r_{ij}} \mathrm{erf} \sqrt{\frac{\alpha_{i} \alpha_{j}}{\alpha_{i} + \alpha_{j}}} r_{ij}$  
    $\displaystyle - \frac{1}{r_{ij}} \mathrm{erf} \sqrt{\frac{\alpha_{i} \alpha_{j}^{\mathrm{max}}}{\alpha_{i} + \alpha_{j}^{\mathrm{max}}}} r_{ij}$  

Only interactions between a charge and another charge with $ \alpha > \alpha_{\mathrm{max}}$ are considered. Places the energy in the variable energy. Also computes the negative gradients and places them in the variables fx, fy, fz, and falpha.


void SelfEnergy(double *falpha, double *energy)

Computes the self-energy:

\begin{displaymath}
E_{\mathrm{self}} = -\frac{1}{2} \sum_{i} q_{i}^{2} \frac{2}...
...pha_{i} + \alpha_{i}}} & \mathrm{otherwise}
\end{array}\right.
\end{displaymath}

Places the energy in the variable energy. Also computes the negative gradient with respect to the exponent and places it in the variable falpha.


void UniformChargeEnergy(double *falpha, double *energy)

Computes the interaction of each charge with a uniform neutralizing background:

$\displaystyle E_{\mathrm{uniform}} = \frac{1}{4} Q \sum_{i} \frac{\pi}{V} \left...
...\frac{1}{\alpha_{i}}\right) \mathrm{\ for\ \alpha_{i} > \alpha_{\mathrm{max}}}
$

Places the energy in the variable energy. Also computes the negative gradient with respect to the exponent and places it in the variable falpha.


double EwaldRCutoff(double precision, double alpha_cutoff, double min_alpha)

Computes the real-space lattice cutoff needed to achieve the negative log10 precision goal precision $ p$, assuming $ \alpha_{\mathrm{max}}$ equals alpha_cutoff and that we want to describe charges that have an exponent as small as $ \alpha_{\mathrm{min}} =$ min_alpha:

$\displaystyle r_{\mathrm{cutoff}} = \sqrt{\frac{(-\ln 10) \cdot p + 3}{\alpha_{\mathrm{max}} / \left(1 + \alpha_{\mathrm{max}} / \alpha_{\mathrm{min}}\right)}}
$


double EwaldKCutoff(double precision, double alpha_cutoff)

Computes the reciprocal-space lattice cutoff needed to achieve the negative log10 precision goal precision $ p$, assuming $ \alpha_{\mathrm{max}}$ equals alpha_cutoff:

$\displaystyle k_{\mathrm{cutoff}} = \sqrt{4 \alpha_{\mathrm{max}} \left((-\ln 10) \cdot p + 5\right)}
$


eff_freeze

This file keeps track of and applies coordinate constraints. When a coordinate is constrained, the forces corresponding to that coordinate are set to zero at every minimization or dynamics step via the function ApplyFreezes().

New constraints are added via functions such as AddFreezeNucX(); these functions are called in eff_input when in the input file the character # appears after a nuclear or electron coordinate.


Data structures

The constraints themselves are dynamically allocated and stored in a linked list declared via eff_listint. Thus the number of constraints that can be added is limited only by the system's memory. Separate lists are used to store constraints on different coordinates: freeze_nuc_x, freeze_nuc_y, freeze_nuc_z, freeze_elec_x, freeze_elec_y, freeze_elec_z, freeze_elec_r.


void InitializeFreezes()

Initializes the linked lists containing the constraints. Uses eff_listint.


void ClearFreezes()

Deallocates the constraint lists via eff_listint.


void AddFreezeNucX(int i)

Adds a constraint to the x-coordinate of nucleus i.


void AddFreezeNucY(int i)

Adds a constraint to the y-coordinate of nucleus i.


void AddFreezeNucZ(int i)

Adds a constraint to the z-coordinate of nucleus i.


void AddFreezeElecX(int i)

Adds a constraint to the x-coordinate of electron i.


void AddFreezeElecY(int i)

Adds a constraint to the y-coordinate of electron i.


void AddFreezeElecZ(int i)

Adds a constraint to the z-coordinate of electron i.


void AddFreezeElecR(int i)

Adds a constraint to the r-coordinate of electron i.


void ApplyFreezes()

Applies all of constraints by setting to zero the forces on any coordinates present in the constraint lists. Called in eff_dynamics and eff_minimize.


eff_global

Contains global data structures shared by many source files: nucleus data, electron data, and parameters.

To keep the program modular, the nucleus and electron data should only be accessed through the routines in eff_access and eff_create. However two source files break this guideline and access the global data structures nuc[ ] and elec[ ] directly: eff_update and eff_ewald. If changes are made to the nucleus and electron data structures, routines in those source files should be updated as well.

Every user-defined parameter is stored as a global variable accessible to all routines in the program. There have been some efforts to restrict the usage of these parameters to eff_driver but they have been only somewhat successful.

Nucleus data

  typedef struct
  {
    double x, y, z, q;
    double fx, fy, fz, energy;
  } NUCLEUS;

The variables x, y, z, and q are the nuclei position and charge. The variables fx, fy, fz, and energy are the nuclear forces and single particle energy.

An extern variable of type (NUCLEUS *) is defined in eff_global, which is allocated as an array with NumNuclei() entries in eff_create. Hence any source file which includes eff_global has access to variables of the form nuc[i].*.

Electron data

  typedef struct
  {
    double x, y, z, r;
    int spin;
    double fx, fy, fz, fr, energy;
  } ELECTRON;

The variables x, y, z, and r are the electron positions and size. The variable spin is the electron spin, and can be either plus or minus one. The variables fx, fy, fz, fr, and energy are the electron forces and single particle energy.

An extern variable of type (ELECTRON *) is defined in eff_global, which is allocated as an array with NumElectrons() entries in eff_create. Hence any source file which includes eff_global has access to variables of the form elec[i].*.

Parameters

Tables 6.1 and 6.2 lists the variables are defined inside the typedef PARAMS. In some cases, enumerated types are used and their entries are hidden here; however, their definitions in the header file are self-explanatory.

An extern variable of type PARAMS is defined in eff_global. Hence any source file which includes eff_global has access to variables of the form params.*.


Table 6.1: Parameter global data structure, part I
Parameter Description Default
Cutoffs    
double taper_cutoff General pair cutoff 10000.0
Boundaries    
double x_bound[2] Rectangular box xmin, xmax -5000, 5000
double y_bound[2] Rectangular box ymin, ymax -5000, 5000
double z_bound[2] Rectangular box zmin, zmax -5000, 5000
periodic Ewald bounds, min image, none? none
Ewald    
double ewald_log_precision -log 10 desired Ewald precision -6.0
double ewald_re_cutoff Real/reciprocal space switchover 3.54
int ewald_autoset Automatically set r and k-space cutoffs? yes
double ewald_r_cutoff Manually set r-space cutoff 7.0
double ewald_k_cutoff Manually set k-space cutoff 8.0
double ewald_nuc_r Radius of nucleus in Ewald 1e-10
double ewald_max_re Maximum electron size Ewald can describe 4.5
Minimization    
min Conjugate gradient or Newton-Raphson? CG
min_freeze Freeze nuclei and/or electrons? no
Dynamics    
thermostat None, Nose-Hoover, or Anderson none
double andersen_coupling Andersen collision rate 0.1
double nose_hoover_coupling Nose-Hoover effective mass 1.0
double start_temperature Initial T for rand. velocities and thermostat 0
double end_temperature End temperature ramp for thermostat 0
double dt Time step for integrator 0.005
adaptive_step_size On or off? off
double adaptive_energy Energy conservation goal 0.0001
int adaptive_num_tries Number of subdivisions to attempt 5



Table 6.2: Parameter global data structure, part II.
Updating    
int print_every Output information every print_every steps 50
int num_steps Total number of dynamics or minimization steps to take 10000
Action    
calc Single point, dynamics, or minimization single point
Output    
output_position None, every step, or only at end all
output_velocity None, every step, or only at end all
output_energy_forces None, every step, or only at end none
output_restart None, every step, or only at end all
output_restraints None, every step, or only at end all
Masses    
double electron_mass Dynamics mass of electron to use 1.0
Random numbers    
unsigned int rand_seed Pseudorandom number seed 10000
External field    
double Ex, Ey, Ez Amplitude of electric field 0, 0, 0
double Efreq Frequency of electric field (0 for constant) 0
double Epacket_duration Length of wave packet (0 for continuous) 0
Electron size limit    
int size_limit Size limit is on or off? on
double size_limit_stiffness Constant in harmonic restraining potential 1.0



eff_initial

This file keeps track of individual initial velocities defined in the nuc_velocities and elec_velocities section of the .cfg input file. It does the same for the individual masses defined in the nuc_masses and elec_masses section as well. This file performs bookkeeping only, and does not actually set the initial velocities or masses.

Accessor functions are provided to add and retrieve the entries. Entries are added in eff_input by functions such as ProcessNucVelocities() when the input .cfg file is read in.

Entries are retrieved in eff_driver by the RunDynamics() function, which then calls eff_dynamics to initialize the velocities and masses. Note that this is done after the default velocities and masses have been set.


Data structures

The initial velocities and masses are stored in fixed arrays, unlike the constraints in eff_freezes or the restraints in eff_restraints, which are stored in linked lists. Hence the number of initial velocities and masses must be known ahead of time, so that the arrays can be allocated.

The arrays used are nuc_v_array, which contains components of (vx, vy, vz); elec_v_array, which contains components of (vx, vy, vz, vr); nuc_m_array, which contains components of (mx, my, mz); and elec_m_array, which contains components of (mx, my, mz, mr).


void AllocateInitialVelocities(int max_nuc_v, int max_elec_v)

Allocates the initial velocity arrays with max_nuc_v nuclear velocities and max_elec_v electron velocities.


void AllocateInitialMasses(int max_nuc_m, int max_elec_m)

Allocates the individual mass arrays with max_nuc_m nuclear masses and max_elec_m electron masses.


int GetNumNucV()

Returns the number of nuclear velocities stored.


int GetNumElecV()

Returns the number of electron velocities stored.


int GetNumNucM()

Returns the number of nuclear masses stored.


int GetNumElecM()

Returns the number of electron masses stored.


void AddNucV(int i, double vx, double vy, double vz)

Adds a nuclear velocity (vx, vy, vz) corresponding to nucleus i.


void AddElecV(int i, double vx, double vy, double vz, double vr)

Adds an electron velocity (vx, vy, vz, vr) corresponding to electron i.


void GetNucV(int idx, int *i, double *vx, double *vy, double *vz)

Gets stored entry idx, filling up the variables i, the nucleus index, and (vx, vy, vz), the nucleus velocity.


void GetElecV(int idx, int *i, double *vx, double *vy, double *vz, double *vr)

Gets stored entry idx, filling up the variables i, the electron index, and (vx, vy, vz, vr), the electron velocity.


void AddNucM(int i, double mx, double my, double mz)

Adds a nuclear mass (mx, my, mz) corresponding to nucleus i.


void AddElecM(int i, double mx, double my, double mz, double mr)

Adds an electron mass (mx, my, mz, mr) corresponding to electron i.


void GetNucM(int idx, int *i, double *mx, double *my, double *mz)

Gets stored entry idx, filling up the variables i, the nucleus index, and (mx, my, mz), the nucleus mass.


void GetElecM(int idx, int *i, double *mx, double *my, double *mz, double *mr)

Gets stored entry idx, filling up the variables i, the electron index, and (mx, my, mz, mr), the electron mass.


eff_input

Reads in the .cfg input file, which contains calculation parameters, nuclear positions and charges, electron positions and sizes and spins, nuclear and electron initial velocities, nuclear and electron individual masses, and restraints.

Once the data is read, updates the eFF system using functions in eff_params, eff_create, eff_freezes eff_initial, and eff_restraints.

Also contains auxiliary functions to match strings, strip newlines, and parse strings.


int CountLines(FILE *fp, char *header)

Returns the number of lines under a given header, such as nuclei. Used to allocate array sizes before the data has been read, required in eff_create and eff_initial.


void ProcessInputFile(FILE *fp)

Reads in the .cfg file line-by-line. Determines which section each line belongs to, and sends it to the appropriate processing function, e.g. ProcessParam() or ProcessNucleus(). Also removes the newline character.


void ProcessNucleus(char *str)

From a string, reads in the nucleus position and charge, and adds it to the eFF system using AddNucleus() in eff_create. Also looks for the character # appended to any coordinate and calls functions in eff_freeze if necessary.


void ProcessElectron(char *str)

From a string, reads in the electron position, spin, and radius, and adds it to the eFF system using AddElectron() in eff_create. Also looks for the character # appended to any coordinate and calls functions in eff_freeze if necessary.


void ProcessNucVelocities(char *str)

From a string, reads in an initial nuclear velocity, and adds it to the eFF system via eff_initial.


void ProcessElecVelocities(char *str)

From a string, reads in an initial electron velocity, and adds it to the eFF system via eff_initial.


void ProcessNucMasses(char *str)

From a string, reads in an initial nuclear mass, and adds it to the eFF system via eff_initial.


void ProcessElecMasses(char *str)

From a string, reads in an initial electron mass, and adds it to the eFF system via eff_initial.


void ProcessParam(char *str)

From a string, reads in a parameter and its arguments, and sends it to ExecuteParameter() in eff_ param. The source file eff_param contains a system to dynamically add commands and function pointers, and to match parameters to their associated functions.


void ProcessRestraints(char *str)

From a string, reads in a restraint, and adds it to the eFF system via eff_restraints.


void NucElecIndex(char *str, int *index)

Helper function for ProcessRestraints(). Recognizes the labels used to refer to nuclei and electrons, and converts them to a numerical form the program can use. In restraints, nuclei are referenced by their number (e.g. 5), and electrons by a number with the prefix e (e.g. e5).


int string_match(char *str1, char *str2)

Returns 1 if str1 and str2 match, disregarding case.


void StripNewline(char *str)

Removes the terminating newline character from str, if present.


void ParseString(char *buffer, int *argc, char **argv)

Tokenizes the input string buffer, recognizing both the tab character and space as separators, and fills argc, the number of substrings, and argv *, an array of pointers to the substrings. Both argc and argv must be defined and properly allocated before the function call.


eff_listint

This file contains functions to handle linked lists of integers, which the program uses in situations where a list is needed, but the number of entries is not known in advance.

Accessor functions are provided to create, adds elements to, and destroy the lists. Reading elements from the linked list is accomplished by accessing the underlying data structure directly, and not through accessor functions.

The integer lists are used in eff_update to keep track of which particles are within different spatial bins, and in eff_freeze to keep track of coordinate constraints.


Data structure

The file defines two structure types, LIST_INT for the list itself, and LIST_INT_ELEMENT for the elements of the linked list.

LIST_INT contains pointers to the head and tail elements of the list:

  typedef struct
  {
    LIST_INT_ELEMENT *head;
    LIST_INT_ELEMENT *tail;
  } LIST_INT;
LIST_INT_ELEMENT contains the list element itself, as well as a pointer to the next element of the list:
  typedef struct LIST_INT_ELEMENT
  {
    struct LIST_INT_ELEMENT *next;
    int element;
  } LIST_INT_ELEMENT;
To read elements of a list that has already been created, the following loop can be used in any file that includes eff_listint:
  LIST_INT integer_list;
  ... (initialize and fill with elements) ...

  LIST_INT_ELEMENT *ptr;
  for (ptr = integer_list->head; ptr != 0; ptr = ptr->next)
  {
    printf("Element %i\n", integer_list->element);
  }


LIST_INT *new_list_int()

Creates a new integer list, e.g.
  LIST_INT integer_list = new_list_int();


void delete_list_int(LIST_INT *list)

Deallocates a list. Before calling this function, call clear_list_int to deallocate the individual elements of the list.


void add_list_int(LIST_INT *list, int newelement)

Adds an element to the tail of the list, e.g.
  add_list_int(integer_list, 4);


void clear_list_int(LIST_INT *list)

Deallocates all the individual list items; the list is empty afterwards.


eff_listptr

This file performs the functions as eff_listint, but with elements of type (void *) instead of integers. Since the data structures and prototype functions are identical in all but element type, they will not be repeated here. Void pointers require a type cast function, e.g. (double *) ptr, before they are considered valid pointers and can be accessed.

The pointer lists are used in eff_restrain to keep track of restraints that appear in the restraints section of the .cfg file.


eff_minimize

This file interfaces with the Shanno minimization code in eff_shanno, initializing position and gradient arrays, updating gradient arrays with values from the eFF system, and updating positions with values from the minimizer.


Data structures

As in eff_dynamics, all the nuclear and electronic degrees of freedom are mapped to a single one-dimensional position array min_x:

  nuc(1)  -- x, y, z
  nuc(2)  -- x, y, z
  ...
  elec(1) -- x, y, z, log r
  elec(2) -- x, y, z, log r
  ...
Thus the array has dimension $ 3 N_{\mathrm{nuc}} + 4 N_{\mathrm{elec}}$, where $ N_{\mathrm{nuc}}$ is the number of nuclei and $ N_{\mathrm{elec}}$ is the number of electrons.

A key difference from the variable mapping versus eff_dynamics is that $ \log r$ is stored instead of r. This ensures that $ r>0$ even though the variables range unconstrained from $ --\infty$ to $ \infty$ during the minimization. It also improves the efficiency of the minimizer when small electrons are involved, for example core electrons around high $ Z$ nuclei.

The array min_g stores the gradients of the total potential energy with respect to the variables in min_f.

Finally, an array min_temp is defined which serves as temporary storage for the Shanno minimizer. Its size depends on the minimization scheme selected. If min_x and min_g have dimension $ N$, the storage array has size $ 5 N + 2$ for the conjugate gradient optimizer, and size $ N (N + 7) / 2$ for the Newton-Raphson optimizer.


void AllocateMinimize(enum MinimizeType minmethod)

Allocates space for the min_x, min_g, and min_temp arrays. The space set aside for min_temp depends on the minimization algorithm minmethod selected.


int MinStorageRequirements(enum MinimizeType minmethod)

Computes the storage requirements for the min_temp array, using the formula given in the Data Structures section, which is a function of the minimization algorithm minmethod selected.


void CalcContractFG(double *x, double *f, double *g)

A callback function that gives the scalar function value f and the function gradient vector g at a position vectorx. In eFF, the f is the potential energy of the system, g are the forces on the nuclei and electrons, and x are the positions of the nuclei and electrons.

The function operates using the helper functions GetMinimizeForces() to fetch eFF system forces, SetMinimizePositions() to set eFF system positions, and UpdateEnergyForces() from eff_update to compute the energy of the system and the forces on the particles.

The function also applies restraints, constraints, and external electric fields. Note that the restraint energy and forces are included explicitly in the minimization; however those energy components are not included in the sum GetTotalPE() in eff_update.


enum MinimizeResult Minimize(enum MinimizeType s_minmethod, double eps, double acc, int max_numsteps, int print_every, void (*MinCallback)(int, int, double, double))

Interface function to the Shanno algorithm in eff_shanno. Passes along the method type, either conjugate gradient or Newton, an epsilon value eps, the desired accuracy acc, the maximum number of steps to take max_numsteps, how often to print out energy and position information (print_every), and the callback function MinCallback which is called every print_every steps with information on the iteration, function value, and gradient.


void GetMinimizeForces(double *forces)

Maps eFF system forces into the array forces. Used by CalcContractFG() to determine the energy f and gradients g at a given position x.


void GetInitialPosition(double *positions)

Maps the current eFF position into the array positions. Used by Minimize() to start the Shanno minimizer off at the current system position.


void SetMinimizePositions(double *positions)

Maps positions from the array positions into the eFF system. Used by CalcContractFG() to determine the energy f and gradients g at a given position x. The position x is assigned by the Shanno algorithm.


eff_output

This file provides functions for writing to output .eff files and restart .cfg.restart files.

Writing to .eff files is done via calls to individual functions such as OutputEnergyForces() and OutputVelocities(). This arrangement simplifies switching off and on different outputs, for example, by the parameters output_energy_forces and output_velocity.

For modularity, all of the writing to the .eff file should be done through eff_output. This rule is violated in a few places: (1) in eff_driver, minimization and dynamics iteration updates are written directly, as are adaptive step size updates, (2) in eff_restrain, information on the different restraints is written out directly, and (3) in eff_update, out of bounds warnings are written out directly.

Writing to the .cfg.restart file is done via a single call to OutputRestartFiles(). If the input file format is changed, this function should be changed too.


void OutputEnergyForces(FILE *out_fp)

Outputs individual particle energies and forces. See section [*].


void OutputPositions(FILE *out_fp)

Outputs individual particle positions. See section [*].


void OutputVelocities(FILE *out_fp)

Outputs individual particle velocities. See section [*].


void OutputEnergy(FILE *out_fp)

Outputs the dipole moment, total energy, and pressure at the end of the calculation. See section [*].


void OutputTimeElapsed(FILE *out_fp)

Outputs the elapsed time at the end of the calculation. See section [*].


void OutputRestartFile(char *restart_filname)

Outputs the restart file, given the filename restart_filename. The parameters are printed out with a call to OutputParams(), and all the other quantities are retrieved via accessor functions: nuclear and electron coordinates from eff_access, nuclear and electron velocities from eff_dynamics, individually specified masses from eff_initial, and restraints with a call to OutputRestartRestraints() in eff_restraints.


void OutputFixedHeader(FILE *out_fp)

Outputs the header at the start of the .eff file, which contains fixed quantities relating to the nuclei and electrons. See section [*].


void OutputDynamicsMasses(FILE *out_fp)

Outputs individually specified masses into the .eff file.


void OutputParams(char *tag, FILE *out_fp)

Outputs the calculation parameters. This function does double duty, it is used by OutputRestartFile() to write to .cfg.restart files, and by eff_driver to write to .eff files. The tag prepends a prefix to each line, and is set to [params] when writing to .eff files.

The parameters are read from the global data structure params.*. Some helper functions are used to translate enumerated data types into strings for output. Their purpose is clear, and so their names are listed below without further elaboration: PeriodicTypeToString, MinimizeTypeToString, MinimizeFreezeTypeToString, ThermostatTypeToString, ActionTypeToString, OutputTypeToString, AdaptiveStepSizeTypeToString, and BoolToString.


void PrintAnalyticForces()

This function prints out to the standard output all the nuclear and electron forces. It is used for debugging purposes only, and is present in commented-out form in eff_driver.


void PrintNumericalForces()

This function prints out to the standard output all the nuclear and electron forces, evaluated using finite-difference numerical derivatives. It is used for debugging purposes only, and is present in commented-out form in eff_driver.


eff_params

This file parses calculation parameters that appear in the params section of the .cfg input file, and sets the corresponding variable in the global data structure params.*. It also sets default values for parameters via the function SetDefaultParameters().

The program provides a general scheme for adding parameters. The function AddParameters( ) associates each parameter name to a function which parses the parameter string. Multiple parameters can map onto a single parsing function, which can then contain switch statements to handle the multiple parameters.

Often parameters that are in the same logical category are parsed together. For example, the parameters x_bound, y_bound, and z_bound are parsed by ParamXYZBound, and the parameters ewald_log_precision and ewald_re_cutoff are parsed by ParamEwald. This division makes to code easier to understand than having a single global parsing function which handles all the parameters.

All of the parameters are added under the function InitializeParameters(), called in eff_driver. Once this is done, the function ExecuteParameter is called in eff_input as the parameters appear in the .cfg input file.


Data structures

The function AddParameters( ) associates each parameter name with a function pointer via a structure of type PARAM_FUNCTIONS:
  typedef struct
  {
    char name[255];
    void (*fp)(int, char **);
  } PARAM_FUNCTIONS;
All of the name/function combinations are stored in a fixed-length array param_functions[ ]. Currently the array size is set at 100; this value should be increased as needed.

The parser functions take as arguments the integer argc, which are the number of words present on the parameter string, and the array of pointers argv, which point to the individual words present on the parameter string. For example, the string

  x_bound = 5
will have argc equal to 3, and
  argv[0] = "x_bound"
  argv[1] = "="
  argv{2] = "5"


void InitializeParameters()

Initializes the name/function parameter array params_functions[ ], and associates all of the parameter names with their parsing functions by calling AddParameter().


void AddParameter(char *name, void (*fp)(int, char **))

Associates the parameter name with the parsing function fp. Called repeatedly in InitializeParameters().


void ExecuteParameter(char *name, int argc, char **argv)

Executes the parsing function matched to name, with the argument count argc and the argument list argv. Called by eff_input.


void SetDefaultParameters()

Sets the parameters to the default values listed in Table 6.1.


Parsing functions

Table 6.3 lists the parsing functions and the parameter names they handle. The purpose of the parsing functions is clear, and so no further elaboration is given.


Table 6.3: Parameter parsing functions and parameter names.
Parsing function Parameters
ParamCutoff taper_cutoff
ParamXYZBound x_bound, y_bound, z_bound
ParamPeriodic periodic
ParamEwald ewald_log_precision, ewald_re_cutoff
  ewald_autoset, ewald_r_cutoff, ewald_k_cutoff
  ewald_nuc_r, ewald_max_re
ParamMin min
ParamMinFreeze min_freeze
ParamThermostat thermostat
ParamNoseHooverCoupling nose_hoover_coupling
ParamAndersenCoupling andersen_coupling
ParamTemperature start_temperature, end_temperature
ParamAdaptiveStepSize adaptive_step_size
ParamAdaptiveEnergy adaptive_energy
ParamAdapativeNumTries adaptive_num_tries
ParamDt dt
ParamPrintEvery print_every
ParamNumSteps num_steps
ParamCalc calc
ParamOutput output_position, output_velocity,
  output_energy_forces, output_restart,
  output_restraints
ParamElectronMass electron_mass
ParamRandSeed rand_seed
ParamEField e_field
ParamEFieldFreq e_field_freq
ParamEFieldPacketDuration e_field_packet_duration
ParamSizeLimit size_limit
ParamSizeLimitStiffness size_limit_stiffness



void ReadValue(char *str, double *val)

A helper function for the parser functions. Translates the string str to a double val, with the provision that the strings -infinity and infinity are translated into the global constant INFINITY, defined in eff_global as 1000000.


eff_pressure

This file computes the instantaneous pressure of a periodic system using the virial expression
$\displaystyle P$ $\displaystyle =$ $\displaystyle \frac{dE}{dV}$  
  $\displaystyle =$ $\displaystyle \frac{N k_{B} T}{V} + \frac{1}{3 V} \sum_{i} \mathbf{r_{i}} \cdot \mathbf{f_{i}}$  
  $\displaystyle =$ $\displaystyle \frac{1}{3 V} \left(2 \sum_{i} \frac{1}{2} m_{i} v_{i}^{2} + \sum_{i} \mathbf{r_{i}} \cdot \mathbf{f_{i}}\right)$  
  $\displaystyle =$ $\displaystyle \frac{1}{3 V} \left(\mathbf{ke\_virial} + \mathbf{pe\_virial}\right)$  

where N is the number of particles, T is the instantaneous temperature, V is the volume, $ \mathbf{r_{i}}$ and $ \mathbf{f_{i}}$ are the particle positions and forces, and $ m_{i}$ and $ v_{i}$ are the particle masses and velocities.

The virial sum is built up in pieces throughout the calculation. Each time a force on a particle is calculated in eff_update, a function call of the form AddForceVirial() is made to update the potential energy virial. Each time a velocity is computed, an analogous call is made in eff_dynamics to AddKineticEnergyVirial() to update the kinetic energy virial. The potential energy virial for electrostatic sums has a particularly simple form, and is handled through AddPotentialEnergyVirial().

Two separate pressures are computed, one assuming that electron sizes change upon compression, and one that assumes that they do not. The program keeps track of these virial sums separately, and when the virial add functions are called, the argument electron_type specifies which sum to update. For example, when electron_type is rigid, no radial force contribution is included, as opposed to when electron_type is flexible.


Data structures

The program tracks the virial sums in the four variables rigid_ke_virial, rigid_pe_virial, flexible_ke_virial, and flexible_pe_virial.

The variable electron_type can be set to the values rigid, flexible, or both. If the virial add functions are called with the parameter rigid, only the rigid virial is updated; with flexible, only the flexible virial is updated; and with both, both of the virials are updated.


void ClearVirial()

Resets the virial sums to zero.


void AddForceVirial(double rx, double ry, double rz, double rr, double fx, double fy, double fz, double fr, enum ElectronType electron_type)

Adds to the potential energy virial the term $ r_{x} f_{x} + r_{y} f_{y} + r_{z} f_{z} + r_{r} f_{r}$. The parameter electron_type determines which virial sums are updated.


void AddSizeForceVirial(double r, double f, enum ElectronType electron_type)

Adds to the potential energy virial the term $ r \cdot f$.


void AddRadialForceVirial(double r, double f, enum ElectronType electron_type)

Adds to the potential energy virial the term $ r \cdot f$. The same as the AddSizeForceVirial().


void AddPotentialEnergyVirial(double energy, enum ElectronType electron_type)

Adds to the potential energy virial the value energy. Used in computing the potential energy virial for Coulomb systems, where a simplification occurs: $ \sum_{i} \mathbf{r_{i}} \cdot \mathbf{f_{i}} = U$, where $ U$ is the potential energy.


void AddKineticEnergyVirial(double energy, enum ElectronType electron_type)

Adds to the kinetic energy virial the value $ 2 \cdot \mathbf{energy}$.


double GetRigidPEPressure(double volume)

Returns the rigid potential energy component of pressure, which equals $ \frac{1}{3V} \cdot \mathbf{rigid\_pe\_virial}$. Units are in GPa.


double GetFlexiblePEPressure(double volume)

Returns the flexible potential energy component of pressure, which equals $ \frac{1}{3V} \cdot \mathbf{flexible\_pe\_virial}$. Units are in GPa.


double GetRigidKEPressure(double volume)

Returns the rigid kinetic energy component of pressure, which equals $ \frac{1}{3V} \cdot \mathbf{rigid\_ke\_virial}$. Units are in GPa.


double GetFlexibleKEPressure(double volume)

Returns the flexible kinetic energy component of pressure, which equals $ \frac{1}{3V} \cdot \mathbf{flexible\_ke\_virial}$. Units are in GPa.


eff_properties

Computes properties of the system at the end of the calculation.


void CalcDipole(double *dipole_x, double *dipole_y, double *dipole_z)

Computes the dipole moment of the system (dipole_x, dipole_y, dipole_z) in atomic units.


double ElectronDensity(double x, double y, double z)

Returns the electron density at position (x, y, z). Currently not used in the code.


eff_restrain

This file track of and applies geometry restraints, for example keeping two atoms a fixed distance apart from each other. The restraints are enforced via the function ApplyRestraints(), which acts by applying harmonic restraining potentials to sets of atoms. The function is called at every minimization or dynamics step.

These potentials have a minimum when the geometric condition is satisfied, and have steep walls to push the system toward the desired configuration. However, there may be cases where large energy barriers are present and opposing forces resist the applied restraint, and in the end the restraint may not be fully satisfied.

The function calls and data structures are similar to those in eff_freeze. New restraints are added via functions such as AddDistanceRestraint(); these functions are called in eff_input when lines in the @restraints section are read.

The tether restraint is a special kind of restraint, in that it averages and evenly distributes forces on a given a set of particles. In the case of dynamics, the forces are scaled by particle masses so that the accelerations are equal for all particle in the set, and the initial velocities are averaged over the set as well.

Developers wishing to add new restraints should note that for the tether restraint to work properly, it must be applied last in the function ApplyRestraints().

Restraints can act on both nuclei and electrons together. In the input file, for example, a distance restraint setting the distance between nucleus 2 and electron 2 to 1.0 bohr is written as

  @restraints
  distance 2 e2 1.0
Internally, the program translates the indices ``2'' and ``e2'' into particle indices 2 and $ 2 + N_{\mathrm{nuclei}}$ via the function NucElecIndex() in eff_input. Helper functions in eff_restrain then take these particle indices as arguments. Examples include AddNucElecForce(), GetNucElecPosition(), and GetNucElecForce().


Data structures

The restraints are dynamically allocated and stored in a linked list declared via eff_listptr. Thus the number of restraints that can be added is limited only by the system's memory. Separate lists are used to store different types of restraints: restraint_coordinate, restraint_distance, restraint_angle, restraint_dihedral, restraint_line, restraint_plane, and restraint_tether.

Data structures for the individual restraints containing particle indices and the restraint value. For example, the data structure for the dihedral restraint is as follows:

  typedef struct
  {
     int idx1, idx2, idx3, idx4;
     double angle;
  } RES_DIHEDRAL;

In some cases, the restraint contains a variable number of particles. In these cases - restraint_line, restraint_plane, and restraint_tether - the data structure is defined as follows:

  typedef struct
  {
    int *indices;
    int numindices;
  } RES_TETHER;
containing a pointer to an integer array indices, where the size of the array is numindices.

The linked list stores void pointers, which must be cast to the appropriate type using an explicit cast. For example, to apply the planar restraints, the following code is used:

  LIST_PTR_ELEMENT *ptr;   // ptr->element has type (void *)
  for (ptr = restraint_plane->head; ptr != 0; ptr = ptr->next)
  {
    RES_PLANE *element = (RES_PLANE *) ptr->element;
    restraints_energy += ApplyPlaneRestraint(element->indices, element->numindices, 0);
  }

The harmonic restraints have spring constants K_COORD, K_DISTANCE, K_ANGLE, K_DIHEDRAL, K_LINE, and K_PLANE.


void InitializeRestraints()

Initializes empty linked lists for all the restraint types. Called once in eff_driver when the program starts.


void ClearRestraints()

Clears all the restraint lists. Called once in eff_driver when the program ends.


void AddCoordinateRestraint(int idx, double x, double y, double z)

Adds a restraint forcing particle idx to coordinate (x, y, z).


void AddDistanceRestraint(int idx1, int idx2, double dist)

Adds a restraint forcing the distance between particles idx1 and idx2 to dist.


void AddAngleRestraint(int idx1, int idx2, int idx3, double angle)

Adds a restraint forcing the angle between particles idx1, idx2, and idx3 to angle.


void AddDihedralRestraint(int idx1, int idx2, int idx3, int idx4, double angle)

Adds a restraint forcing the dihedral between particles idx1, idx2, idx3, and idx4 to angle.


void AddLineRestraint(int *indices, int numindices)

Adds a restraint forcing the numindices particles listed in indices onto a line. There must be at least three particles for this restraint to be valid.


void AddPlaneRestraint(int *indices, int numindices)

Adds a restraint forcing the numindices particles listed in indices onto a plane. There must be at least four particles for this restraint to be valid.


void AddTetherRestraint(int *indices, int numindices)

Adds a restraint that averages and evenly distributes the forces on the numindices particles listed in indices.


double ApplyRestraints()

Applies all of the added restraints to the eFF system. Does this by applying functions such as ApplyCoordinateRestraint() and ApplyDistanceRestraint() to each of the elements of the restraint lists. This function is called once per minimization or dynamics iteration.


void ApplyTetherVelocityRestraints()

Executes the function ApplyTetherVelocityRestraint() over all the stored tethers. This function is called once in the entire dynamics run, right after the initial velocities have been set.


void AddsNucElecForce(int idx, double fx, double fy, double fz)

Adds to particle idx the force (fx, fy, fz). The particle indices are discussed in the Data Structures section above.


void GetNucElecPosition(int idx, double *x, double *y, double *z)

Gets the position of particle idx and places it in (x, y, z).


void GetNucElecForce(int idx, double *fx, double *fy, double *fz)

Gets the force on particle idx and places it in (fx, fy, fz).


void GetNucElecMass(int idx, double *mx, double *my, double *mz)

Gets the mass of particle idx and places it in (mx, my, mz).


void GetNucElecVelocity(int idx, double *vx, double *vy, double *vz)

Gets the velocity of particle idx and places it in (vx, vy, vz).


void SetNucElecVelocity(int idx, double vx, double vy, double vz)

Sets the velocity of particle idx to (vx, vy, vz).


void NucElecString(char *str, int idx)

Translates the particle index idx to a string. For example, if the number of nuclei is 5, the index 2 would be translated to ``2'', while the index 5 would be translated to ``e1''.


double ApplyCoordinateRestraint(int idx, double x0, double y0, double z0, double *dist)

Applies a coordinate restraint to particle idx, forcing it to have position (x0, y0, z0). Returns the restraint energy. If dist is not NULL, places into dist the distance between the particle and (x0, y0, z0).


double ApplyDistanceRestraint(int idx1, int idx2, double dist0, double *dist)

Applies a distance restraint to particles idx1 and idx2, forcing them to be separated by distance dist0. Returns the restraint energy. If dist is not NULL, places into dist the distance between the two particles.


double ApplyAngleRestraint(int idx1, int idx2, int idx3, double angle0, double *angle)

Applies an angle restraint to particles idx1, idx2, and idx3 forcing them to subtend the angle angle0. Returns the restraint energy. If angle is not NULL, places into angle the angle between the three particles.


double ApplyDihedralRestraint(int idx1, int idx2, int idx3, int idx4, double d0, double *angle)

Applies an dihedral restraint to particles idx1, idx2, idx3, and idx4 forcing them to subtend the dihedral d0. Returns the restraint energy. If angle is not NULL, places into angle the dihderal between the four particles.


double ApplyLineRestraint(int *indices, int numindices, double *dists)

Applies a line restraint to the numindices particles listed in indices, forcing them to lie on a line. If dists is not NULL, places into the array the distance between the particles and the least-squares line, measured orthogonally.


double ApplyPlaneRestraint(int *indices, int numindices, double *dists)

Applies a plane restraint to the numindices particles listed in indices, forcing them to lie on a plane. If dists is not NULL, places into the array the distance between the particles and the least-squares plane, measured orthogonally.


double ApplyTetherRestraint(int *indices, int numindices, double *dists)

Applies a tether restraint to the numindices particles listed in indices, averaging and evenly distributing the forces on the particles. If a dynamics calculation is being performed, scales the average velocities by the mass of each particle, so that the accelerations are the same; also see ApplyTetherRestraintVelocity() used to set the initial velocities. If dists is not NULL, places zeros into the array (not used).


double ApplyTetherRestraintVelocity(int *indices, int numindices)

Applies a tether restraint to the numindices particles listed in indices, averaging and evenly distributes the velocities on the particles. This only needs to be done once, at the beginning of the dynamics simulation.


void OutputRestraints(FILE *out_fp)

Outputs value of restraints to the .eff output file, handle out_fp.


void OutputRestartRestraints(FILE *restart_fp)

Outputs value of restraints to the .cfg.restart restart file, handle restart_fp.


eff_shanno

This file contains the unconstrained quasi-Newton/conjugate gradient optimization subroutine conmin from D. F. Shanno and K. H. Phua, Algorithm 500. Minimization of unconstrained multivariate functions [E4]. ACM Transactions of Mathematical Software, 2(1):87-94 (1976). The subroutine is free for public use.

We have converted the original Fortran subroutine into a C function, and modified the array indices so that they start at 0 instead of 1.

The subroutine can operate in two modes, a conjugate gradient mode requiring $ 7 n$ words of temporary storage, or a quasi-Newton BFGS (Broyden-Fletcher-Goldfarb-Shanno) mode requiring $ n^{2}/2 + 11 n / 2$ words of temporary storage.

The conmin subroutine is called by the Minimize( ) function in eff_minimize.


int conmin(int n, double *x, double *f, double *g, int *ifun, int *iter, double eps, int mxfun, double *w, int iout, int mdim, double acc, int nmeth, void (*calcfg)(double *, double *, double *), void (*update_func)(int, int, double, double))

From the comments provided in the original conmin algorithm, the arguments are listed in Table 6.4 and the return value error codes are listed in Table 6.5.


Table 6.4: Arguments for Shanno's conmin function in eff_shanno.
Argument Purpose
n Number of variables in the function to be minimized
x Vector containing the initial guess.
  Upon exiting, holds the final vector.
f Upon exiting, holds the final (minimum) function value.
g Upon exiting, holds the final gradient.
ifun Upon exiting, holds the number of function/gradient evaluations.
iter Upon exiting, holds the number of iterations.
eps User-supplied convergence parameter.
mxfun Maximum number of evaluations allowed.
w Temporary storage vector.
iout If 1, enables output.
mdim Dimension of w.
acc User-supplied machine accuracy.
nmeth 1 for BFGS, 0 for CG.
calcfg Callback func(x, f, g) calculates function at x, and returns
  its value in f and the gradient in g.
update_func Callback output(num_iterations, num_fevals, f_val, g_sq).



Table 6.5: Error codes for Shanno's conmin function in eff_shanno.
Return value Meaning
0 Algorithm has converged.
1 Maximum number of function evaluations used.
2 Linear search fails, is gradient wrong?
3 Search vector not a descent direction, roundoff error?



double dmin(double a, double b)

Returns the minimum of $ a$ and $ b$. Used in conmin.


double dmax(double a, double b)

Returns the maximum of $ a$ and $ b$. Used in conmin.


eff_timing

This file implements a stopwatch used to determine how long a calculation takes to complete.

Two timers are available, CPUTimeElapsed(), which has a resolution of less than a second but which overflows after 36 minutes, and TotalTimeElapsed(), which has a resolution of one second but should not overflow. Both values are printed out by eff_output at the end of the calculation.


void StartTimer()

Starts the stopwatch timer.


void StopTimer()

Stops the stopwatch timer.


double CPUTimeElapsed()

CPU time elapsed in between the time StartTimer() and StopTimer() is called.


double TotalTimeElapsed()

User time elapsed in between the time StartTimer() and StopTimer() is called.


eff_update

This file computes the energy and forces associated with an eFF system. In eFF, the overall energy is a sum of a Hartree product kinetic energy, a Hartree product electrostatic energy, and an antisymmetrization (Pauli) correction:

$\displaystyle E = E_{\mathrm{ke}} + E_{\mathrm{nuc \cdot nuc}} + E_{\mathrm{nuc \cdot elec}} + E_{\mathrm{elec \cdot elec}} + E_{\mathrm{Pauli}}
$

which can be broken down further as follows:
$\displaystyle E_{\mathrm{ke}}$ $\displaystyle =$ $\displaystyle \sum_{i} \frac{3}{2} \ \frac{1}{s_{i}^{2}}$  
$\displaystyle E_{\mathrm{nuc \cdot nuc}}$ $\displaystyle =$ $\displaystyle \sum_{i<j} \frac{Z_{i} Z_{j}}{R_{ij}}$  
$\displaystyle E_{\mathrm{nuc(i)\cdot elec(j)}}$ $\displaystyle =$ $\displaystyle -\sum_{i,j} \frac{Z_{i}}{R_{ij}}\ \mathrm{Erf}\left(\frac{\sqrt{2} \, R_{ij}}{s_{j}}\right)$  
$\displaystyle E_{\mathrm{elec(i) \cdot elec(j)}}$ $\displaystyle =$ $\displaystyle \sum_{i<j} \frac{1}{r_{ij}}\ \mathrm{Erf}\left( \frac{\sqrt{2} \, r_{ij}}{\sqrt{s_{i}^{2} + s_{j}^{2}}}\right)$  
$\displaystyle E_{\mathrm{Pauli}}$ $\displaystyle \sum_{\sigma_{i} = \sigma_{j}} E(\uparrow \uparrow)_{ij} + \sum_{\sigma_{i} \neq \sigma_{j}}
E(\uparrow \downarrow)_{ij}$  

where $ E(\uparrow \uparrow)$ and $ E(\uparrow \downarrow)$ are the Pauli potential functions:
$\displaystyle E(\uparrow \uparrow)_{ij}$ $\displaystyle =$ $\displaystyle \left(\frac{S_{ij}^{2}}{1-S_{ij}^{2}} + (1 - \rho) \frac{S_{ij}^{2}}{1 + S_{ij}^{2}}\right) \Delta T_{ij}$  
$\displaystyle E(\uparrow \downarrow)_{ij}$ $\displaystyle =$ $\displaystyle \frac{\rho S_{ij}^{2}}{1 + S_{ij}^{2}} \Delta T_{ij}$  

where $ \Delta T$ is a measure of the kinetic energy change upon antisymmetrization, and $ S$ is the overlap between two wave packets:
    $\displaystyle \Delta T_{ij} = \frac{3}{2}\left(\frac{1}{\bar{s}_{1}^{2}}+\frac{...
...2}+\bar{s}_{2}^{2})-2 \bar{r}_{12}^{2})}{(\bar{s}_{1}^{2}+\bar{s}_{2}^{2})^{2}}$  
    $\displaystyle S_{ij} = \left(\frac{2}{\bar{s}_{i}/\bar{s}_{j}+\bar{s}_{j}/\bar{s}_{i}}\right)^{3/2} \exp(-\bar{r}_{ij}^{2}/(\bar{s}_{i}^{2}+\bar{s}_{j}^{2}))$  

where $ \rho = -0.2$, $ \bar{x}_{ij} = x_{ij} \cdot 1.125$, and $ \bar{s}_{i} = s_{i} \cdot 0.9$.

The main function UpdateEnergyForces() updates the energy and forces, calling other functions in eff_update to compute individual energy components, and them summing them all together. The final results can be retrieved using functions from eff_access such as GetTotalPE(), GetNuclearForce(), and GetElectronForce().

For systems whose dimensions are larger than the pairwise cutoff taper_cutoff, eFF saves on computational cost by partitioning space into rectangular bins smaller than the cutoff length, and only considering interactions between nearest-neighbor bins. The routine AssignParticleBins() assigns the particles to the bins, and a routine inside UpdateEnergyForces() calls UpdateElecAndPauli() only over nearest-neighbor bins.

When such systems have periodic boundaries, the unit cell is so large that each particle can only ``see'' particles within the same unit cell. We treat the interactions using the minimum image convention, and wrap around space so that bins on one edge of the system can interact with bins on the opposite edge of the system. Functions such as LeftStencil() and RightStencil() help manage this process.

To apply the minimum image convention, we use functions such BoundDX(), which convert interatomic distances within a periodic cell to a minimum image distance. In smaller periodic systems, each particle may see multiple images of itself and particles from other cells, and a more sophisticated algorithm is needed to sum long range forces. We sum electrostatic forces using the Ewald algorithm via UpdateEwaldEnergy() in eff_ewald. The Pauli potential is shorter range, and we assume the minimum image convention holds - we perform this sum using UpdatePauliPeriodic(). The kinetic energy in all cases is updated with the same function UpdateKineticEnergy().

Many functions in this source file are declared inline to improve performance, such as ierfoverx1(), which evaluates $ \mathrm{erf}(x)/x$. All gradients are computed analytically.


Data structures

We partition space into rectangular bins, and each bin $ (n_{x}, n_{y}, n_{z})$ has an index

$\displaystyle i = n_{x} + n_{y} \cdot \mathrm{numbins\_x} + n_{z} \cdot \mathrm{numbins\_x} \cdot \mathrm{numbins\_y}
$

To keep track of which particles are in which bin, we associate with each bin index $ i$ a list of nuclei indices nuc_bin[i] and a list of electron indices elec_bin[i]. The variables nuc_bin and elec_bin are arrays of pointers to integer lists, and have type (LIST_INT **).


void IntitializeEnergyUpdater()

Initializes the rectangular bins. The number of bins along a given direction is calculated as follows:
  numbins_x = (int) (length_x / params.taper_cutoff);
  if (numbins_x == 0) numbins_x = 1;
Hence the program tries to break space down into segments the size of the cutoff, but if fractional units are left over, the program stretches out the bin slightly so that all the bins fit evenly into space. Each bin is guaranteed to be at least the cutoff length in size along every dimension.

The program also determines whether the minimum image convention needs to be applied along the x, y, and z dimensions, based on the user-set parameters, and sets the flags minimage_x, minimage_y, and minimage_z accordingly.


int LeftStencil(int numbins, int pbc)

The left range of bins that can interact with a bin, given the number of bins numbins along a direction and whether or not periodic boundary conditions are present along that direction (pbc). Usually equal to -1, so that a nearest neighbor bin can interact with a central bin, but some cases need to be handled specially.


int RightStencil(int numbins, int pbc)

Similar to LeftStencil, but represents the right range of bins, and is usually equal to +1.


void AssignParticleBins()

Assigns particles to the rectangular bins. Clears the lists, then for each particle determine which bin contains it, then adds the particles to the bins. If a particle doesn't belong to any bin, prints out a warning. This function is called every time UpdateEnergyForces() is called.


void UpdateEnergyForces()

Updates the total potential energy, and the forces on individual particles. Begins by computing the kinetic energy, then branches into two cases.

If Ewald summation is being used, computes electrostatics via UpdateEwaldEnergy() in eff_ewald, and Pauli repulsion via UpdatePauliPeriodic().

Otherwise, breaks space down into bins using AssignParticleBins(), and call UpdateElecAndPauli() over all nearest-neighbor bins that are expected to interact with each other, accounting for periodic boundary conditions if present.


void UpdateKineticEnergy()

Computes the total kinetic energy and associated radial forces.


void UpdateElecAndPauli(int bin1, int bin2, int wrap_x, int wrap_y, int wrap_z)

Computes the total electrostatic interaction and Pauli repulsion between particles in bin1 and bin2. The flags wrap_x, wrap_y, and wrap_z control whether or not the minimum image convention is applied in the x, y, and z dimensions.


void UpdatePauliPeriodic()

Computes the Pauli repulsion for a periodic system using the minimum image convention.


void UpdateWidthRestraint()

Applies a restraint potential $ V = \frac{1}{2} k_{s} s^{2}$ to electrons that have a radius greater than half the size of the smallest rectangular box dimension.


eff_util

This file is a collection of miscellaneous utility functions to report errors, generate random numbers, convert strings from uppercase to lowercase, get information on elements, and handle file extensions.


void error(const char *format, ...)

Reports an error and exits the program. Called with the same arguments as printf:
  error("Error: 1 plus 1 is %i\n", 3);


int rand_int(int max)

Returns an integer from 0 to $ \mathrm{max} - 1$.


double rand_uni()

Return a double from 0 to 1 (will never be exactly 1), distributed uniformly.


double rand_gauss()

Returns a double from a Gaussian distribution, mean 0, standard deviation 1.


void lcase(char *str)

Converts the string str to a lower case string.


double GetMass(int Z)

Returns the atomic mass in amu of the element with nuclear charge Z.


char *GetType(int Z)

Returns the element symbol (H, He, Li, etc.) of the element with nuclear charge Z.


int GetZ(char *str)

Returns the nuclear charge Z of the element with element symbol str.


void GetBasename(char *str, char *basename)

Returns the basename of a filename basename.extension.


void GetExtension(char *str, char *extension)

Returns the extension of a filename basename.extension.


void rand_seed(unsigned int seed)

Seeds the pseudorandom generator with seed.


double min(double x, double y)

Returns the minimum of x and y.


next up previous contents index
Next: Index Up: eFF Programmers' Guide Previous: To access hidden features   Contents   Index
Julius 2008-04-29