Parameters

Actions

calc = single_pt

Options: single_pt, minimize, or dynamics. Single_pt computes the energy of a single configuration of nuclei and electrons, without changing any of the coordinates. See also output_energy_forces, which controls the output of energies and forces on individual particles.

Minimize varies the nuclear and electron coordinates to minimize the total energy of the system, using gradient descent algorithms that find the nearest local minimum. See also min, which specifies the kind of descent algorithm to use, min_freeze, which permits the freezing of nuclei or electrons during minimization, and num_steps and print_every. Constraints and restraints are also applicable during minimization.

Dynamics propagates the nuclear and electron coordinates via semiclassical Ehrenfest dynamics using a velocity Verlet algorithm, making it possible to explore the quantum dynamics of excited systems. See also dt for specifying the time step; thermostat, start_temperature, and end_temperature for performing NVT dynamics; electron_mass for setting the dynamic electron mass; and num_steps and print_every. Constraints and restraints are also applicable during dynamics.

Updating

print_every = 100

The program will output geometry and energy information into the .eff and .cfg.restart files every print_every iterations. See also output_position, output_velocity, output_energy_forces, output_restart, and output_restraints, which set which quantities are output at those points.

num_steps = 10000

Sets the number of dynamics or minimization steps.

Minimization

min = conjugate_gradient

Options: conjugate_gradient or newton. Both methods are gradient descent methods, using information about the total energy and the its gradient to approach the nearest local minimum over several consecutive iterations. The conjugate gradient algorithm requires $O(N)$ memory, and is the method of choice for most systems. The quasi-Newton algorithm finds the minimum in fewer steps but requires $O(N^{2})$ memory. It may be preferred for small systems with up to hundreds of atoms and electrons.

min_freeze = none

Options: none, nuclei, or electrons. Freezes all nuclei or all electron coordinates during minimization. This is a basic option, and there are more selective ways to freeze particles available. To freeze individual nuclear or electron coordinates, use constraints instead, specified by # appended to a given coordinate in the @nuclei or @electron section in the input file. To fix relative distances or angles, to fix particles relative to each other, or to force them into a plane or line, use restraints, specified in a separate @restraints section.

Dynamics

thermostat = none

Options: none, andersen, or nose-hoover. By default, the program performs NVE dynamics with a velocity Verlet algorithm, keeping the number of particles N and the volume V fixed while conserving energy E. The option is also available to perform NVT dynamics, so that the temperature T rather than the energy is kept fixed. This is done by applying a thermostat.

The Andersen thermostat uses a stochastic method to create the correct ensemble of velocities. At each time step, a particle has a chance, related to the parameter andersen_coupling, of having its velocity reset to a vector picked out of a Maxwell-Boltzmann distribution at the desired temperature. This simulates random collisions with an external heat bath.

The Nose-Hoover thermostat is a deterministic method where the particles are slowed down or sped up by a velocity dependent force. The force is proportional to the difference between the system's instantaneous temperature and the desired temperature, and also to the coupling parameter nose_hoover_coupling. This method is derived from an extended Lagrangian which contains an extra degree of freedom corresponding to the external heat reservoir.

andersen_coupling = 0.1

Units: $\mathrm{fs}^{-1}$. Andersen thermostat coupling parameter $\alpha$, corresponding to a collision rate. The probability that a particle will collide with the heat bath and have its temperature reset is set to be $P_{\mathrm{collide}} = \alpha \cdot dt$ for a differential time interval $dt$. Hence the probability that some collision will occur over a time $t$ is $1 - e^{-\alpha t}$. In the program, we assume that the time step $\Delta t$ is small enough that we can take $P_{\mathrm{collide}} = \alpha \cdot \Delta t$.

If $\alpha$ is too low, the thermostat is inefficient, but if $\alpha$ is too high, random collisions dominate and the dynamics is destroyed.

nose_hoover_coupling = 1

Units: $\mathrm{amu} \cdot \mathrm{bohr}^{2}$. Nosé-Hoover thermostat coupling parameter Q, corresponding to an effective mass for the Nosé-Hoover parameter $\zeta$. In the modified equations of motion, we have

$$m \ddot{r} = \vec F - \zeta m \dot {r}$$

$$\zeta' = \frac{1}{Q}\left(\sum_{i} \frac{1}{2} m \dot {r}^{2} - \frac{3}{2} N k_{B} T\right)$$

Thus $\zeta$ acts as a drag force on the particles, or equivalently a time-rescaling factor, while Q determines the response of $\zeta$ to the difference between the actual and desired system temperature.

If Q is too large, the thermostat is inefficient, but if Q is too small, the temperature will oscillate too much, and the momentum fluctuations will not couple properly to the heat bath.

start_temperature = 0

Units: Kelvin. Sets the initial velocities of the nuclei to a Maxwell-Boltzmann velocity distribution with the desired temperature. Entries in the @velocities section take precedence over this option, so that if the velocity of a specific nucleus is given, it will overwrite any random initialization generated by this option.

When the Andersen or Nosé-Hoover thermostats are selected, also sets the target temperature for the thermostat.

end_temperature = 0

Units: Kelvin. This option is only applied when a thermostat is present, so that the program is attempting the regulate the temperature throughout the simulation. If the value is the same as the one given in start_temperature, there is no effect. If the value is explicitly set to be different than start_temperature, the option produces a linear temperature ramp that stretches from the start to the end of the simulation.

dt = 0.005

Units: femtoseconds. Time step for dynamics. The program uses atomic units internally, but with the modification that the unit mass is taken to be one atomic mass unit, or one-twelfth the mass of a carbon atom. Thus the internal time unit (1 atomic time unit, or 1 atu) is

$$dt = \frac{\hbar}{E_{h}} \cdot \sqrt{\frac{\mathrm{amu}}{m_{e}}} = 1.03275 \times 10^{-15}\ \mathrm{seconds}$$

which is close to 1 fs. We bring this up because there are two similar time units in play: (1) The time step for dynamics is input in femtoseconds, and the conversion factor is applied internally. (2) Velocities are specified in units of bohr/atu and not bohr/fs.

The minimum time step is determined by the fastest oscillations in the system, usually the vibrations of the core electrons surrounding the heaviest nuclei. Typically time steps of 0.01 to 0.02 fs are reasonable, but for studying Auger processes, a very short time step of 0.001 fs is needed. It is possible to freeze core electrons relative to nuclei by specifying a tether in the @restraints section, which makes longer time steps possible.

electron_mass = 1

Units: atomic mass units. The dynamic electron mass $m_{e}$ relates the forces on the electron wave packets to the rate of change of their translational and radial velocities:

$$\mathbf{{\dot p}_{x}} = -\nabla_{x}V, \ \mathbf{p_{x}} = m_{e} \mathbf{\dot{x}}$$

$$\dot{p}_{s} = -\partial V/\partial s, \ p_{s} = (3 m_{e} /4) \dot{s}$$

When the wave packet is in a harmonic potential, $m_{e}$ should be the true mass of the electron, 1/1822.89 amu = 0.000548579 amu.

In more complex potentials, however, the dynamic electron mass may be much higher. It can be estimated in semiconductors using band theory, or in metals using many-body theory. In eFF, $m_{e}$ is a fixed parameter that we set to reproduce some characteristic known time-scale of electron motion in the system we are studying. In general, the time step of the dynamics should be scaled to $m_{e}^{-1/2}$.

A large $m_{e}$ may be used to simulate adiabatic ground-state potentials, as in Car-Parinello dynamics. In some special cases, such as the motion of electrons near a core-hole, we have found it necessary to use a large $m_{e}$ to reproduce experimental core-hole lifetimes. Smaller $m_{e}$ on the order of the true electron mass is appropriate for studying the response of molecules to high electric fields, proton stopping in solids, and Auger ionization.

The dynamic electron mass $m_{e}$ is separate from the electron mass present in the kinetic energy term in Schrödinger's equation, and changing $m_{e}$ does not affect electron sizes or bonding in the ground state.

adaptive_step_size = false

Options: false or true. Shrinks the step size to make sure energy is conserved. If adaptive_step_size = true and thermostat = none, the program checks to make sure that the total energy does not change by more than adaptive_energy between iterations. If the energy change exceeds adaptive_energy, the program rolls back to the previous iteration, scales the step size dt by a factor of 0.5, and tries the iteration again. The process is repeated until the energy difference critereon is satisfied. After the iteration is taken, the step size is reset to its original value.

adaptive_energy = 0.0001

Units: Hartrees. Maximum allowed difference in energy per iteration when adaptive_step_size = true.

Electric field

e_field = 0 0 0

Units: Volts per centimeter. Electric field vector $\mathbf{E_{0}} = [E_{x}, E_{y}, E_{z}]$ to apply. Affects minimization and dynamics calculations. By default the field is constant over time, but it may be modified by the parameters e_field_freq and e_field_packet_duration below to be a continuous sine wave or a wave packet.

e_field_freq = 0

Units: Hertz. Electric frequency $f$. If e_field_freq is zero, the field is assumed to have constant amplitude. Otherwise, the applied field is

$$\mathbf{E} = \mathbf{E_{0}} \sin \left(2 \pi f t\right)$$

if the e_field_packet_duration is zero.

e_field_packet_duration = 0

Units: femtoseconds. Electric field packet duration $t_{\mathrm{packet}}$. If both e_field_freq and e_field_packet_duration are nonzero, the applied field is

$$\mathbf{E} = \mathbf{E_{0}} \sin \left(2 \pi f t\right) \exp\left(-\frac{(t - 3 t_{\mathrm{packet}})^{2}}{2 t_{\mathrm{packet}}^{2}}\right)$$

so that the standard deviation of the wave packet is $t_{\mathrm{packet}}$, and the wave packet begins at time zero three standard deviations away from the maximum.

Electron size limit

size_limit = true

Options: true or false. For periodic systems, applies a harmonic restraint to limit the electron size, using the potential $V = \frac{1}{2} k_{s} {s}^{2}$ for $s>L_{\mathrm{box}}/2$, where $L_{\mathrm{box}}$ is the smallest dimension of the rectangular bounding box. If this potential is not applied, electrons once ionized can expand without limit, causing errors in the computed pressure. See also size_limit_stiffness.

size_limit_stiffness = 1

Units: hartree per $\mathit{bohr}^{2}$. Harmonic stiffness parameter $k_{s}$ for the size_limit restraint.

Output

output_position = all

Options: none, all, or end. Output positions of nuclei and electrons at intervals defined by the print_every parameter. Positions are in bohr.

output_velocity = all

Options: none, all, or end. Output velocities of nuclei and electrons at intervals defined by the print_every parameter. Velocities are in bohr per atu ($\sim$ 1.03 fs).

output_energy_forces = none

Options: none, all, or end. Output energy and forces on individual nuclei and electrons at intervals defined by the print_every parameter. Energies are in hartrees, forces are in hartrees per bohr.

output_restart = all

Options: none, all, or end. Updates the restart file at intervals defined by the print_every parameter. Each update overwrites the old restart file, so that only the last point calculated is kept in the restart file.

output_restraints = all

Options: none, all, or end. Output value of restraints at intervals defined by the print_every parameter.

Random numbers

rand_seed = 10000

Specifies a random number seed used to initialize the master pseudorandom number generator. If the seed is the same between two runs, the same random numbers will be generated, useful for reproducing the results of a prior run. Conversely, if different random initial conditions are desired for a series of runs, it is essential to use a different seed for each calculation.

Ewald parameters

ewald_re_cutoff = 3.54

Units: bohr. The Ewald algorithm operates by separating the long-range electrostatic sum into real space and reciprocal space sums. The parameter ewald_re_cutoff acts as the threshold $s_{\mathrm{cutoff}}$ between real and reciprocal sums. Charges with $s < s_{\mathrm{cutoff}}$ are summed over real space with multiple replicas, while charges with a larger width are summed over reciprocal space, leaving a equal and opposite combination of charges to be summed over real space.

ewald_autoset = true

Options: true or false. Automatically estimates the real space cutoff ewald_r_cutoff and reciprocal space cutoff ewald_k_cutoff needed to achieve the precision specified by the variable ewald_log_precision. Depends on the variables ewald_re_cutoff and ewald_max_re as well.

The cutoffs are set by the program to minimize CPU time, assuming a certain ratio of computational expense between real space and reciprocal space sums. The formulas are

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

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

where $\alpha_{\mathrm{cutoff}} = 2 / s_{\mathrm{cutoff}}^{2}$ ( $s_{\mathrm{cutoff}} =$ ewald_re_cutoff), $\alpha_{\mathrm{min}} = 2 / s_{\mathrm{max}}^{2}$ ( $s_{\mathrm{max}} =$ ewald_max_re), and $p =$ ewald_log_precision.

If either ewald_r_cutoff or ewald_k_cutoff are set manually, the program assumes the user wishes to override the automatic settings, and the autoset feature is turned off.

ewald_log_precision = -6

Desired precision of the Ewald sum, log base 10 of the energy in Hartrees. Used in the ewald_autoset procedure.

ewald_max_re = 4.5

Units: bohr. Largest electron the Ewald sum should be able to describe. Used in the ewald_autoset procedure.

ewald_r_cutoff = 7

Units: bohr. Real space cutoff for replica sum. Setting this will turn ewald_autoset to false.

ewald_k_cutoff = 8

Units: $\mathit{bohr}^{-1}$. Reciprocal space cutoff for lattice sum. Setting this will turn ewald_autoset to false.

ewald_nuc_r = 1e-10

Units: bohr. Radius of nuclei, taken to be very sharply localized Gaussian charges.

Periodic boundary conditions

periodic = none

Options: true, minimage_x, minimage_y, minimage_z, minimage_xy, minimage_xz, minimage_xz, minimage_xyz, false. eFF provides methods for handling periodic boundaries in two cases: (1) small unit cells, where each particle is expected to interact with multiple copies of itself and (2) large unit cells, where each particle ``sees'' every other particle at most once (minimum image convention).

For case (1), periodic = true computes electrostatics using an Ewald sum, which uses a combined real and reciprocal space sum to efficiently sum periodic replicas. Pauli repulsion interactions are shorter range and are summed using a minimum image convention. When this option is used, the pairwise cutoff splines are turned off, and the parameter taper_cutoff has no effect.

For case (2), the periodic options starting with minimage_ compute both electrostatic interactions and Pauli repulsion using a minimum image convention. For this method to be valid, the pairwise cutoff taper_cutoff must be smaller than $L_\mathrm{box}/2$, where $L_{\mathrm{box}}$ is the minimum dimension of the rectangular unit cell. 1D and 2D and 3D periodic bounds can be treated in this way. For example, periodic = minimage_xy creates a system that is periodic in the $x$ and $y$ directions.

bound_x = -10000 10000

bound_y = -10000 10000

bound_z = -10000 10000

Units: bohr. Rectangular boundaries of the system. In a dynamics simulation, if periodic boundaries are enabled for one of the dimensions, particles exiting on one side will enter through the opposite side. Otherwise, particles that leave the boundaries will have no force or energy associated with them, and will continue indefinitely in a straight line at constant velocity.

Pairwise cutoff

taper_cutoff = 1000

Units: bohr. All pairwise interactions are multiplied by a taper function which is 1 at zero distance and gradually levels off to 0 at distance $d_{\mathrm{cutoff}} = \textbf{taper\_cutoff}$. This enables the overall computational expense to be linear scaling with the number of particles, and makes calculations with hundreds of thousands of electrons practical.

The function is the seventh-order spline

$$f = 20 x^{7}-70 x^{6} + 84 x^{5} - 35 x^{4} + 1$$

where $x = d / d_{\mathrm{cutoff}}$, and it is defined such that the first and second and third derivatives are zero at the endpoints. It differs from traditional molecular dynamics cutoffs in that it is operative over the entire range $d < d_{\mathrm{cutoff}}$ rather than a narrow range $d_{\mathrm{inner}} < d < d_{\mathrm{outer}}$. This minimizes discontinuities in the energy function or its derivatives.

A $d_\mathrm{cutoff}$ of 20 bohr is considered an aggressive cutoff with probable errors of $\sim$ 2 kcal/mol/atom, while 50 bohr is a mild cutoff with probable errors on the order of $< 0.01$ kcal/mol/atom. If taper_cutoff is larger than half of the smallest box dimension, and minimum image periodic boundaries are applied, the program will terminate with an error.