Periodic Boundary Conditions

TeraChem can perform periodic electronic-structure calculations on three-dimensional (bulk) systems using periodic Gaussian basis functions and Ewald summation for the long-range Coulomb interaction. Both Hartree–Fock and DFT are supported, for single-point energies and analytic gradients.

Γ-point only

Periodic calculations are evaluated at the Γ point of the Brillouin zone — there is no k-point sampling. Converge properties with respect to system size by enlarging the simulation cell (supercells) rather than by sampling k-points. Only 3D periodicity is implemented; 1D (chains) and 2D (slabs) are not currently supported.

Activating a periodic calculation

Two things are required: the keyword periodic yes, and a unit cell defined on the comment line of the coordinate file.

periodic   yes

DFT additionally requires the periodic integral/grid engine, enabled with dftbox yes (see Methods below). Hartree–Fock works with periodic yes alone.

Specifying the unit cell

The unit cell is given on the second line (the comment line) of the XYZ coordinate file, in the format

3D  a  b  c  [alpha  beta  gamma]

where a b c are the cell-edge lengths in Ångström and alpha beta gamma are the cell angles in degrees. The angles are optional and default to 90°. The dimension token must be 3D.

For example, the conventional cubic cell of diamond (8 carbons):

diamond.xyz

8
3D 3.56679 3.56679 3.56679 90 90 90
C   0.0000000   0.0000000   0.0000000
C   0.0000000   1.7833950   1.7833950
C   1.7833950   1.7833950   0.0000000
C   1.7833950   0.0000000   1.7833950
C   2.6750925   0.8916975   2.6750925
C   0.8916975   0.8916975   0.8916975
C   0.8916975   2.6750925   2.6750925
C   2.6750925   2.6750925   0.8916975

A monoclinic cell simply supplies the non-trivial angle, e.g. a benzene crystal cell: 3D 5.5220 5.4396 7.6726 90 110.55 90.

Atomic positions are Cartesian (in Ångström) by default. To give them as fractional coordinates of the cell instead, add fractional_xyz yes to the input; TeraChem then converts them to Cartesian internally using the cell vectors.

Quick start

A minimal periodic Hartree–Fock single point:

periodic_hf.in

method        hf
run           energy
basis         6-31g
coordinates   benzene.xyz      # cell on the comment line: "3D a b c ..."
charge        0
spinmult      1
periodic      yes
guess         sad
precision     mixed
threall       1.0e-14
purify        no
end

Example 1: periodic Hartree–Fock energy

benzene_rhf_energy.in

method     hf
run        energy
basis      6-31g
coordinates benzene.xyz
charge     0
spinmult   1
periodic   yes
guess      sad
precision  mixed
threall    1.0e-14
purify     no
timings    yes
end

Example 2: periodic DFT gradient with dispersion

Periodic DFT requires dftbox yes. Empirical dispersion corrections (dftd) work in the periodic setting, including the three-body Axilrod–Teller–Muto term (dftd_atm):

benzene_blyp_d3_gradient.in

method      blyp
run         gradient
basis       6-31g
coordinates benzene.xyz
dftbox      yes
dftgrid     3
dftd        d3-bj
dftd_atm    yes
dftd_print  yes
charge      0
spinmult    1
periodic    yes
guess       sad
precision   mixed
threall     1.0e-14
purify      no
end

Methods and run types

Setting	Notes
`method hf`	Periodic Hartree–Fock. No `dftbox` needed.
`method <functional>` + `dftbox yes`	Periodic DFT. The `dftbox` engine provides the periodic integrals and grid.
`run energy`	Periodic single-point energy.
`run gradient`	Analytic nuclear gradient (enables geometry optimization and Born–Oppenheimer MD on the periodic energy).
`dftd <model>`	Empirical dispersion (e.g. `d3-bj`), including `dftd_atm` three-body term.

Integration grids (DFT)

Periodic DFT exchange–correlation can be integrated on either grid:

Molecular (atom-centered) grid — the default, selected with the usual dftgrid N keyword. This is the grid used in Examples above; its construction is documented under Integration Grids.
Uniform real-space grid — request a regular grid spanning the cell with either an explicit size or a target density:

periodic_uniform_grid_size     100,100,100   # nx,ny,nz grid points
# or
periodic_uniform_grid_density  40000         # points per unit cell volume

diamond_blyp_uniform_grid.in

method     blyp
run        energy
basis      6-31g
coordinates diamond.xyz
dftbox     yes
periodic   yes
periodic_uniform_grid_size 100,100,100
charge     0
spinmult   1
guess      sad
precision  double
threall    1.0e-14
purify     no
end

Range separation and Ewald summation

The Coulomb interaction is split into a short-ranged real-space part and a long-ranged reciprocal-space part controlled by the Ewald parameter \(\omega\).

By default (\(\omega\) not specified) TeraChem chooses \(\omega\) automatically to balance the computational cost of the two lattice sums:

\[ \omega = \frac{2^{1/6}\,\sqrt{\pi}}{V^{1/3}}, \]

where \(V\) is the unit-cell volume (in atomic units, bohr³), giving \(\omega\) in bohr⁻¹. This comes from approximating the number of lattice points inside the real-space cutoff (\(r_\text{real}\approx 1/\omega\)) as \(r_\text{real}^3/V\) and inside the reciprocal-space cutoff (\(r_\text{recip}\approx 2\omega\)) as \(r_\text{recip}^3 V/(8\pi^3)\), then minimizing the weighted sum \(w\,n_\text{real} + n_\text{recip}\) over \(\omega\) (with the real-space term weighted \(w=2\)), which yields \(\omega^6 = w\pi^3/V^2\).

A larger \(\omega\) shifts work into the reciprocal-space sum; a smaller \(\omega\) into the real-space sum. Override the automatic value with periodic_omega only for accuracy/timing studies:

periodic_omega 0.1

Lattice-sum cutoffs

The summations over periodic images are truncated using drop tolerances (smaller = tighter = more images retained). These rarely need changing; the defaults are 1.0e-14.

Keyword	Default	Description
`periodic_orbital_cutoff`	1.0e-14	Tolerance for including periodic images of the Gaussian basis functions
`periodic_ewald_cutoff_real`	1.0e-14	Tolerance for the real-space Ewald lattice sum
`periodic_ewald_cutoff_reciprocal`	1.0e-14	Tolerance for the reciprocal-space Ewald lattice sum
`periodic_cutoff_all`	—	Convenience keyword: sets all three cutoffs above to the same value

Summary of keywords

Keyword	Type	Default	Description
`periodic`	bool	no	Master switch — activate 3D periodic boundary conditions
`fractional_xyz`	bool	no	Interpret atomic coordinates as fractional cell coordinates
`dftbox`	bool	no	Required for periodic DFT (periodic integrals and grid)
`periodic_omega`	float	automatic	Ewald range-separation parameter \(\omega\) (bohr⁻¹); see formula above
`periodic_uniform_grid_size`	3×int	—	Uniform real-space grid dimensions `nx,ny,nz` (DFT)
`periodic_uniform_grid_density`	float	—	Uniform real-space grid density (points per cell volume; DFT)
`periodic_orbital_cutoff`	float	1.0e-14	Basis-function image cutoff tolerance
`periodic_ewald_cutoff_real`	float	1.0e-14	Real-space Ewald cutoff tolerance
`periodic_ewald_cutoff_reciprocal`	float	1.0e-14	Reciprocal-space Ewald cutoff tolerance
`periodic_cutoff_all`	float	—	Set all three image/Ewald cutoffs at once

The unit cell itself is not an input keyword — it is read from the comment line of the coordinate file (see Specifying the unit cell).

Limitations

3D periodicity only — 1D and 2D systems are not supported.
Γ-point only — no k-point sampling; use supercells to approach the bulk limit.
Periodic DFT requires dftbox yes.