Repository navigation
BMD Compute Workflows
Authority. BMD Compute is authoritative for how BMD VASP calculations are generated and executed. BMDex supplies curated supporting data, evidence and tools. Scientific validation and adoption are separate human review.
This page explains why BMD Compute's workflows are built from the stages they contain. It does not define them.
This is the second of four pages:
- VASP — fundamentals: the density-functional physics and the core VASP settings.
- BMD Compute Workflows (this page): why calculations are sequences of stages and what passes between them.
-
BMD Compute INCARs:
the complete
INCARs of each workflow and what their settings mean. - BMD Compute Advanced Options: optional physics and method modifiers.
Each workflow section has two parts:
-
BMD workflow facts — what BMD Compute does, taken from the workflow
definitions and the generated
INCARs. These are established. - The physics — the general numerical and physical reasons why sequences like this are used in DFT practice. This is the wiki's explanation. Unless it says otherwise, it is not a statement of BMD Compute's own rationale.
In BMD Compute you choose a workflow by the result you want. Each workflow is a fixed sequence of stages, and each stage is a separate VASP run.
| Desired output | Stages, in order |
|---|---|
| Static Energy | PBE Static Energy |
| Relaxed Structure | PBE Geometry Optimisation → PBE Geometry Optimisation |
| Density of States | PBE Geometry Optimisation → HSE06 Static Energy → HSE06 Density of States |
| Band Structure | PBE Geometry Optimisation → HSE06 Static Energy → HSE06 Band Structure |
The complete INCAR for every stage is on
BMD Compute INCARs.
Depending on the structure, BMD Compute can also add treatments to these stages automatically: spin polarisation (every stage), spin–orbit coupling (every stage except Geometry Optimisation), a DFT-D3(BJ) van der Waals correction (PBE stages) and, for some oxides and fluorides of transition metals, a Hubbard U correction (DFT+U; PBE stages, never HSE06). See BMD Compute Advanced Options.
All stages of a workflow run one after another in a single cluster job. A stage counts as successful only if it converged; if a stage fails, or the job reaches its walltime, the workflow stops, later stages do not run, and no results are loaded. Running it again is a new submission that starts from the first stage: BMD Compute v1 does not continue a workflow from where it stopped.
A single VASP run answers one kind of question well, and different questions need settings that pull in different directions:
- Finding the structure means following forces and stresses downhill on the Born–Oppenheimer energy surface. It needs one self-consistent (SCF) solution per ionic step, and settings that behave well while the geometry changes.
- Establishing the electronic ground state of a fixed structure needs a well-converged basis and Brillouin-zone integration at that geometry.
- Spectral quantities, such as a density of states or a band structure, need Kohn–Sham eigenvalues on k-point sets chosen for that purpose.
Splitting these into stages lets each run use settings suited to its job.
Two kinds of information can, in principle, pass from one stage to the next:
- The structure: the lattice vectors and atomic positions at the end of a stage.
- The electronic state: the converged density and orbitals.
In VASP these correspond to files (a relaxed structure is written to
CONTCAR; the density can be written to CHGCAR and the orbitals to
WAVECAR), and a later run can read them.
What BMD Compute does:
- Each stage after the first starts from the final structure of the stage before it. In the Relaxed Structure workflow, Stage 2 starts from the structure produced by Stage 1; in the Density of States and Band Structure workflows, both HSE06 stages use the PBE-relaxed structure.
- No
WAVECARis carried between stages, and no stage in any of the four workflows writes one (LWAVE = Falsethroughout). - The stages that write
CHGCAR(LCHARG = True) are the PBE Static Energy stage and the HSE06 Static Energy stage. The Geometry Optimisation, HSE06 Density of States and HSE06 Band Structure stages do not. - No stage in the four workflows reads a previous stage's density or orbitals:
none sets
ISTARTorICHARG. The HSE06 Density of States and Band Structure stages are self-consistent, and take only the structure, the number of bands (NBANDS, 1.2 times that of the static stage) and the starting magnetic moments from the HSE06 Static Energy stage. ACHGCARmay be copied into their directory, but it is not read. - The only stages that do read a previous charge density are the PBE Density
of States and Band Structure stages available in Custom workflows. They run
non-self-consistently (
ICHARG = 11) from a PBE Static Energy stage that must come directly before them.
Previews and the inputs actually run. The inputs you preview before
submission are all generated from the structure you supplied. From Stage 2
onwards, the inputs actually run are regenerated from the previous stage's
output, so the structure, starting magnetic moments, k-points (and band path)
and NBANDS can differ from the preview. This is expected.
All k-point meshes are Γ-centred and generated automatically from a density of
points per reciprocal-lattice volume (pymatgen's automatic_density_by_vol),
so the mesh adapts to the cell. Each later stage regenerates its k-points for
the structure it starts from.
| Stage | Sampling |
|---|---|
| PBE Geometry Optimisation, PBE Static Energy, HSE06 Static Energy | Uniform mesh, reciprocal density 64 |
| HSE06 Density of States | Uniform mesh, reciprocal density 100 |
| HSE06 Band Structure | Uniform mesh at reciprocal density 64, plus a high-symmetry path at line density 40 with zero weight |
With SOC, symmetry is switched off, so the uniform part of the band-structure
KPOINTS lists every mesh point with equal weight.
- One stage: PBE Static Energy.
- It is a single-point calculation (
NSW = 0) at the structure you supply, using the tetrahedron method (ISMEAR = -5), and it writes the charge density, local potential, all-electron densities, electron localisation function and projected output.
Scientific result. The self-consistent PBE ground state of the given structure at fixed geometry: its total energy, electron density and electronic structure on the calculation's k-point mesh.
Why there are no earlier stages. The structure is taken as given. This is appropriate when the structure is already the one you want to study, for example an experimental structure or one relaxed consistently elsewhere. The energy is only as meaningful as the structure: the static energy of an unrelaxed structure is not an equilibrium energy, and energies should only be compared between calculations with consistent settings.
- Two consecutive PBE Geometry Optimisation stages.
- Both stages use identical
INCARsettings: variable-cell relaxation of positions, cell shape and volume (ISIF = 3) by conjugate gradients (IBRION = 2), up to 99 ionic steps, force criterion 0.01 eV/Å (EDIFFG = -0.01),ENCUT = 580.0. - Stage 2 starts from the structure produced by Stage 1.
Scientific result. The PBE equilibrium structure (lattice vectors and atomic positions): the structure at which the forces on all atoms and the stress on the cell vanish, to within the convergence criteria.
Why relax twice? Repeating a variable-cell relaxation from its own result is common DFT practice. One well-established reason is the basis-set effect called Pulay stress, explained below. Restarting can also help in other ways, for example when the first run starts far from equilibrium and uses up its ionic steps, or simply to confirm that the structure no longer changes. The wiki does not state BMD Compute's own rationale for the second stage beyond the sequence itself.
The plane waves in the basis are ei(k+G)·r, with
G the reciprocal-lattice vectors. At the start of a run, VASP selects the
G vectors whose kinetic energy ħ²|k+G|²/2m lies below ENCUT
for the initial cell, and keeps that set of plane waves for the whole run.
When a variable-cell relaxation changes the cell, the reciprocal lattice
changes with it, so the kinetic energies of the retained plane waves change.
If the cell contracts, the G vectors lengthen and the fixed set now
reaches above the nominal cutoff; if the cell expands, it falls short of it.
By the end of the run, the basis is no longer the one a fresh run at ENCUT
would use for the new cell.
In a complete basis this would not matter. In a truncated basis, the total energy depends on how complete the basis is, and the stress (the derivative of the energy with respect to strain) picks up a contribution from that dependence rather than from the physics. This spurious contribution is the Pulay stress. Its size depends on how strongly the total energy still changes with the cutoff. Because the energy decreases as the basis grows, the Pulay stress acts to shrink the cell, so relaxations at too low a cutoff tend to underestimate the equilibrium volume.
Two things reduce it:
-
A cutoff well above the POTCAR's recommended
ENMAX. The closer the basis is to complete, the less the energy depends on the cutoff, and so the smaller the Pulay stress. - Restarting the relaxation. A new run rebuilds the plane-wave set for the current cell at the nominal cutoff. If the first relaxation changed the volume substantially, the second relaxation corrects the structure with a basis consistent with the new cell. If the second relaxation barely changes the structure, the result is stable with respect to the basis; if it changes it noticeably, the first result was not reliable on its own.
Neither step removes basis-set error entirely; they reduce the inconsistency between the basis and the final cell.
- Three stages: PBE Geometry Optimisation → HSE06 Static Energy → HSE06 Density of States.
- Stage 1 has the same settings as each Geometry Optimisation stage in the Relaxed Structure workflow. There is one geometry optimisation stage in this workflow, not two.
- Stage 2, the HSE06 Static Energy stage, is a self-consistent single-point
hybrid-functional calculation (
NSW = 0). It is identical in the Density of States and Band Structure workflows. It writesCHGCARbut notWAVECAR. - Stage 3 is also a self-consistent single-point HSE06 calculation
(
NSW = 0), with the tetrahedron method (ISMEAR = -5), a fine energy grid for the DOS (NEDOS = 4001) and site- and orbital-projected output (LORBIT = 11).
Scientific result. The HSE06 electronic density of states,
g(E) = Σn ∫BZ δ(E − En(k)) dk / ΩBZ,
and its site- and orbital-projected components, evaluated at the PBE-relaxed structure. It shows how Kohn–Sham states are distributed in energy, and is used for example to read off band edges and gaps and to see which atoms and orbitals form the valence and conduction bands.
Why the PBE geometry optimisation comes first. Relaxation needs forces and stresses at every ionic step, and so an SCF solution at every step. With a hybrid functional each SCF solution is far more expensive than with PBE, because the exact-exchange operator couples pairs of orbitals across pairs of k-points. Relaxing with PBE gives the structure at a small fraction of the cost.
The result is therefore a hybrid-functional electronic structure evaluated at a PBE-relaxed geometry, not at the HSE06 equilibrium geometry. PBE tends to slightly overestimate the lattice constants of many solids, and band gaps depend on the lattice constant, so this is part of what the result means. When you report it, state both the functional used for the geometry (PBE) and the functional used for the electronic structure (HSE06).
Why there is an HSE06 Static Energy stage. This stage solves the hybrid-functional Kohn–Sham equations self-consistently at the fixed relaxed geometry. It establishes the HSE06 electronic ground state of that structure: the self-consistent orbitals, density and total energy at the hybrid level of theory, together with the full set of static outputs (charge density, local potential, all-electron densities, electron localisation function).
A hybrid Hamiltonian depends on the occupied orbitals, not only on the density, so a PBE ground state does not define an HSE06 one. The static stage is where the hybrid-functional ground state of the relaxed structure is established and recorded. In BMD Compute the final stage does not read the static stage's density or orbitals: it is self-consistent on its own mesh, starting from the same structure and taking its number of bands and starting moments from the static stage (see What passes between stages).
Why the final stage is specific to a density of states. A density of states is an integral over the whole Brillouin zone. It needs eigenvalues on a uniform k-point mesh and an integration method suited to that mesh, such as the tetrahedron method used in Stage 3. BMD Compute uses a denser uniform mesh for this stage than for the static stage (see k-point sampling).
- Three stages: PBE Geometry Optimisation → HSE06 Static Energy → HSE06 Band Structure.
- Stages 1 and 2 have the same settings as Stages 1 and 2 of the Density of States workflow.
- Stage 3 is a self-consistent single-point HSE06 calculation (
NSW = 0) with Gaussian smearing of very small width (ISMEAR = 0,SIGMA = 0.01).
Scientific result. The HSE06 band structure En(k): the Kohn–Sham eigenvalues of each band n along a path of high-symmetry directions in the Brillouin zone, at the PBE-relaxed structure. From it you can see whether the gap is direct or indirect, where the band extrema lie, and how the bands disperse (for example their curvature, which relates to effective masses).
Why the first two stages are needed. For the same reasons as in the Density of States workflow: the PBE stage establishes the geometry at affordable cost, and the HSE06 Static Energy stage establishes the self-consistent HSE06 electronic ground state of that structure. The same point about reporting both functionals applies.
Why the final stage differs from the density-of-states stage. A band structure is not an integral over the Brillouin zone; it is the set of eigenvalues at specific k-points along a path. Two consequences follow:
- Occupations. The tetrahedron method needs a uniform mesh, so it cannot be applied to points along a path. Gaussian smearing with a very small width is used instead, as in Stage 3.
- k-points and self-consistency. With a semilocal functional such as PBE, a converged density fixes the Kohn–Sham potential, and eigenvalues at any other k-points can be computed non-self-consistently in that fixed potential. With a hybrid functional this does not work, because the exact-exchange operator is built from the orbitals on the whole uniform mesh. The standard approach is a self-consistent run on a regular, weighted mesh with the path k-points added at zero weight, so that they receive eigenvalues from the self-consistent Hamiltonian without contributing to the density. Stage 3 follows this approach: it is self-consistent on a uniform mesh with the high-symmetry path added at zero weight (see k-point sampling). The manual walk-through in Calculating Band Structure in VASP shows the zero-weight procedure as a teaching example.