Skip to content

3.1 Water in a periodic vacuum box

These educational templates have not been executed and contain no fabricated numerical results. Inputs and filenames follow ABACUS 3.9.0 documentation; recheck every interface when changing version.

3.1.1 Question and declared model

Can a periodic electronic-structure calculation approximate one neutral water molecule well enough to compare forces and bond geometry? Start with a fixed, deliberately unoptimized molecule. The calculation is a gas-phase electronic model, not liquid water, a molecular crystal, a vibrational spectrum or a finite-temperature free energy. ABACUS still imposes periodicity; the surrounding vacuum reduces interactions with translated images rather than removing boundary conditions.

3.1.2 Physical interpretation

Use neutral closed-shell PBE with two hydrogen atoms and one oxygen atom. For common H and O pseudopotentials with one and six valence electrons respectively, the cell has eight valence electrons and four doubly occupied states. Inspect the actual UPF headers before relying on that count. Six bands provide a starting margin, not a converged excited-state basis. A single Gamma point samples the large molecular box. The molecule is polar, so its image interaction is not eliminated merely because the net charge is zero.

The force is the derivative of the same electronic energy used for comparing geometries:

\[\mathbf F_I=-\frac{\partial E}{\partial\mathbf R_I},\qquad \Delta E_L=E(L)-E(L_{\mathrm{ref}}).\]

Increasing the box changes both image separation and PW computational cost. Keep Cartesian molecular coordinates equivalent after recentering, so this test measures the boundary approximation rather than changing bond lengths.

Scientific schematic: Periodic images; Fixed O-H geometry; Vacuum size L; Check forces too. No measured data.

Open figure at full size

Original course illustration; conceptual geometry and curves, not calculated results.

3.1.3 Worked input and file changes

INPUT

INPUT_PARAMETERS
suffix water_box
calculation scf
basis_type pw
ks_solver cg
ntype 2
nbands 6
nspin 1
dft_functional PBE
pseudo_dir ./data/
ecutwfc 80
scf_thr 1e-9
scf_nmax 150
symmetry 0
cal_force 1
out_chg 1 10

STRU

ATOMIC_SPECIES
O 15.999 O_PBE_VERIFIED.upf
H 1.008 H_PBE_VERIFIED.upf

LATTICE_CONSTANT
1.889726125457828

LATTICE_VECTORS
16 0 0
0 16 0
0 0 16

ATOMIC_POSITIONS
Cartesian_angstrom
O
0.0
1
8.0000 8.0000 8.0000 m 1 1 1
H
0.0
2
8.7586 8.0000 8.5043 m 1 1 1
7.2414 8.0000 8.5043 m 1 1 1

KPT

K_POINTS
0
Gamma
1 1 1 0 0 0

The three blocks above are separate files. Replace the clearly marked pseudopotential filenames with local verified norm-conserving PBE files and retain their checksums. No numerical orbitals are needed for this PW example. The lattice constant is in Bohr, and multiplying it by the vectors gives a 16 Angstrom cubic cell; Cartesian_angstrom keeps atom positions explicitly in Angstrom. The H positions give a plausible starting geometry only. The 80 Ry wavefunction cutoff must be converged for the actual O and H files; it is not 80 eV. Archive the cell, units, valence counts and executable version before any run.

3.1.4 Run and inspect the evidence

Prepare independent 16, 20 and 24 Angstrom boxes, moving the whole molecule to each center without changing internal coordinates. First establish cutoff convergence at one large box, then repeat the box-size comparison at the accepted cutoff. A user may run the installation's documented executable in an authorized working directory; this lesson does not submit jobs. Read OUT.water_box/running_scf.log for SCF convergence, final energy and force units. With out_chg 1 10, version 3.9 writes OUT.water_box/SPIN1_CHG.cube. Inspect its grid and integral before drawing a density isosurface. After accepting the box, use the ionic-relaxation lesson to optimize coordinates without relaxing the vacuum cell.

3.1.5 Observable-specific convergence

Use an explicit teaching target of less than 1 meV per molecule between the largest two boxes and less than 0.005 eV/Angstrom change in any force component. These are proposed tolerances, not measured successes. Check more than energy because a tiny energy change can hide force errors. Translate the molecule inside a fixed box to detect grid-position sensitivity, and increase cutoff if forces depend materially on that translation. A neutral-cell density integral should match the pseudopotential valence-electron count within the integration and output-precision error. Record both SCF termination status and numerical changes; a nonconverged density fails even if adjacent energies look close.

3.1.6 Acceptance worksheet and provenance

For the vacuum worksheet, record box length in Angstrom, cutoff in Ry, final energy in eV and all nine force components in eV/Angstrom. Choose the largest accepted box as the reference. Subtract its molecular energy from each smaller-box energy; do not divide by three merely because there are three atoms when the tolerance is per molecule. For forces, calculate the maximum absolute component difference between equivalent atomic coordinates. Align atom identities before subtracting; the first hydrogen must not be exchanged with the second midway through the table.

Inspect the translation test separately. A rigid displacement inside the same periodic box preserves the physical molecule but changes its position relative to the numerical grid. Subtract translated and original force components after removing only the rigid coordinate shift. Archive the actual grid dimensions from each cube so a changed grid is not hidden behind an unchanged nominal cutoff. A vacuum test that changes both box and cutoff simultaneously cannot isolate the source of the error.

Box L (Angstrom) / cutoff (Ry) Energy difference (eV/molecule) Max force difference (eV/Angstrom) SCF and grid checked?
16 / accepted cutoff — — —
20 / accepted cutoff — — —
24 / accepted cutoff — — —

Keep the PP hashes, coordinate differences and cube header with this blank worksheet. These are planned comparisons, not calculated data.

3.1.7 Pitfalls and limits

Do not use cell-relax to optimize a molecule's vacuum: it changes an artificial boundary, not a physical lattice constant. Do not infer a dipole correction from a VASP keyword; only use an explicitly supported ABACUS option after checking its version and geometry. A one-point mesh does not demonstrate sufficient vacuum. Mixing PBE oxygen with another functional's hydrogen changes the Hamiltonian. Density lobes are visualization choices; their apparent extent depends on the chosen isovalue. A charged molecule requires a separate finite-size and electrostatic treatment beyond this neutral lesson.

3.1.8 Exercises with answer guidance

  1. Increase only the box length from 16 to 20 Angstrom while preserving the O-H vectors. Which coordinates must change? Guidance: recenter all atoms by the same translation; differences between their Cartesian coordinates must remain identical.
  2. Energy passes the box target but one force component changes by 0.012 eV/Angstrom. May the smaller box be used for geometry optimization? Guidance: no; force-sensitive use requires its own accepted image and cutoff errors. Report the failure rather than quoting the energy test alone.

3.1.9 Versioned references

The university-authored guide provides teaching context; older pages can predate 3.9. Use the linked official version for exact syntax. This is original instructional synthesis, and reference outputs are not presented as results of this course.