Basic input options¶
Now we will move to a more complicated (and realistic) input file. We will obtain the ground state of the nitrogen atom. We will introduce several basic input variables and will give a more detailed description of the output for this example.
mkdir -p 2-basic_input_options
Generating the input file¶
This sample input file lets us obtain the ground state of the nitrogen atom, within the LDA approximation, in a closed-shell (unpolarized) configuration (as explained below, you need an auxiliary .xyz input). Note that this is not the correct ground state of the nitrogen atom! However, it will permit us to describe some of the most important input variables:
%%writefile 2-basic_input_options/inp
stdout = "stdout_gs.txt"
stderr = "stderr_gs.txt"
CalculationMode = gs
UnitsOutput = eV_Angstrom
Nitrogen_mass = 14.0
%Species
'N' | species_pseudo | set | standard | lmax | 1 | lloc | 0 | mass | Nitrogen_mass
%
XYZCoordinates = 'N.xyz'
# 2 extra states needed for chebyshev filtering (1 because of degeneracy and 1 for better convergence)
ExtraStates = 2
%Occupations
2 | 1 | 1 | 1
%
BoxShape = sphere
Radius = 5.0*angstrom
Spacing = 0.18*angstrom
Writing 2-basic_input_options/inp
We have introduced here several new variables:
UnitsOutput = eV_Angstrom: Two different unit systems may be used for output: the usual atomic units (which is the default, and the ones used internally in the code); and the system in which the Ångström is substituted for the atomic unit of length, and the electronvolt is substituted for the atomic unit of energy. You can find a more detailed description of units in Octopus in the Units page of the manual.The following entry in the input file is not a variable that Octopus will read directly, but rather illustrates the possibility of writing “user-defined” values and expressions to simplify the input file. In this case, we define the nitrogen mass (
Nitrogen_mass = 14.0) (note that in this case, as an exception, the value is expected to be in the so-called “atomic mass units”, rather than in “atomic units”). This definition may be used elsewhere in the input file.The
Speciesblock should contain the list of species that are present in the system to be studied. In this case we have only one species: nitrogen. The first field is a string that defines the name of the species, “N” in this case. The second field defines the type of species, in this casespecies_pseudo. Then a list of parameters follows. The parameters are specified by a first field with the parameter name and the field that follows with the value of the parameter. Some parameters are specific to a certain species while others are accepted by all species. In our examplesetinstructs Octopus to use a pseudopotential for nitrogen from thestandardset. This happens to be a Troullier-Martins pseudopotential defined in theN.psffile found in the directoryshare/octopus/pseudopotentials/PSF. Then come maximumlmax- component of the pseudopotential to consider in the calculation, and thelloc- component to consider as local. Generally, you want to set the maximum “l” to the highest available in the pseudopotential and the local “l” equal to the maximum “l”. Finally, the mass of the species can also be modified from the default values by settingmassparameter.XYZCoordinates= ‘N.xyz’: The geometry of the molecule (in this case, a single atom in the grid origin) is described in this case in a file with the well knownXYZformat. The file for this outrageously simple case is given by:
%%writefile 2-basic_input_options/N.xyz
1
This is a comment line
N 0 0 0
Writing 2-basic_input_options/N.xyz
ExtraStates= 2: By default, Octopus performs spin-unpolarized calculations (restricted closed-shell, in Hartree-Fock terminology). It then places two electrons in each orbital. The number of orbitals, or Kohn-Sham states, is then calculated by counting the number of valence electrons present in the system, and dividing by two. In this case, since we have five valence electrons, the code would use three orbitals. However, we know beforehand that the HOMO orbital has a three-fold degeneracy, and as a consequence we need to put each one of the three p electrons in a different orbital. In addition to that, the chebyshev_filter eigensolver requires an additional empty state for better convergence. We therefore need two more orbitals, which we get with this line in the input file.%Occupationsblock: Generally, the occupations of the Kohn-Sham orbitals are automatically decided by the code, filling the lowest-energy orbitals. However, if we have degeneracies in the LUMO as in this case, the user may want to accommodate the electrons in a certain predefined way. In this example, the obvious way to fill the orbitals of the nitrogen atom is to put two electrons in the first and deepest orbital (the s orbital), and then one electron on each of the second, third and fourth orbitals (the p orbitals, which should be degenerate).BoxShape = sphere: This is the choice of the shape of the simulation box, which in this case is set to be a sphere (other possible choices areminimum,species_pseudo, orset).Radius = 5.0*angstrom: The radius of the sphere that defines the simulation box.Spacing = 0.18*angstrom: As you should know, Octopus works in a real-space regular cubic mesh. This variable defines the spacing between points, a key numerical parameter, in some ways equivalent to the energy cutoff in plane-wave calculations.
Output¶
Once you have constructed the input file and created the .xyz file, you may unleash Octopus on it. Lets now go over some of the sections of the output.
!cd 2-basic_input_options && octopus
Species¶
****************************** Species *******************************
Species 'N'
type : pseudopotential
file : '/app/share/octopus/pseudopotentials/PSF/N.psf'
file format : PSF
valence charge : 5.0
atomic number : 7
form on file : semilocal
orbital origin : calculated
lmax : 1
llocal : 0
projectors per l : 1
total projectors : 1
application form : kleinman-bylander
orbitals : 16
bound orbitals : 16
**********************************************************************
Here the code searches for the needed pseudopotential files, and informs the user about its success or failure. In this case, only the N.psf file is required. Once that file has been processed, some information about it is written to the output. One of the most important pieces of information to be found here is the valence charge, which tells us how many electrons from this species will be considered in the calculation.
Grid¶
******************************** Grid ********************************
Simulation Box:
Type = sphere
Radius [A] = 5.000
Main mesh:
Spacing [A] = ( 0.180, 0.180, 0.180) volume/point [A^3] = 0.00583
# inner mesh = 89727
# total mesh = 127183
Grid Cutoff [eV] = 1160.586810 Grid Cutoff [Ry] = 85.301565
**********************************************************************
This step is about the construction of the mesh. As requested in the input file, a sphere of radius 5 Å is used, which contains a cubic regular real-space grid with spacing 0.18 Å. This implies 89727 points (inner mesh = 89727). For the sake of comparison with plane-wave-based codes, this is more or less equivalent to a plane-wave calculation that imposes a density cutoff of 1160.595 eV = 42.6 Hartree (except that in this case there is no artificial periodic repetition of the system).
Mixing¶
Input: [MixField = potential] (what to mix during SCF cycles)
Input: [MixingScheme = broyden]
During the self-consistent procedure one has to use a mixing scheme to help convergence. One can mix either the density or the potential, and there are several mixing schemes available.
Eigensolver¶
**************************** Eigensolver *****************************
Input: [Eigensolver = chebyshev_filter]
Input: [SubspaceDiagonalization = standard]
**********************************************************************
Here we see that the eigensolver used will be chebyshev_filter.
LCAO¶
After some output you should see something like:
Info: Performing initial LCAO calculation with 4 orbitals.
Info: Getting Hamiltonian matrix elements.
Eigenvalues [eV]
#st Spin Eigenvalue Occupation
1 -- -17.398903 2.000000
2 -- -6.429122 1.000000
3 -- -6.429122 1.000000
4 -- -6.429122 1.000000
Generating random wavefunctions for states 5 and above
Orthogonalizing wavefunctions.
Info: Ground-state restart information will be written to '/builds/octopus-code/octopus/doc/jupyter_tutorials/1-basics/2-basic_input_options/restart//gs'.
This is the first step of a ground-state calculation: obtaining a reasonably good starting density and Kohn-Sham orbitals to feed in the self-consistent (SCF) procedure. For this purpose, Octopus performs an initial calculation restricted to the basis set of atomic orbitals ( Linear Combination of Atomic Orbitals, LCAO). The resulting eigenvalues of this calculation are written to standard output.
Wavefunction kind¶
Info: SCF using real wavefunctions.
Very often one can work with real wave-functions. This is particularly helpful as calculations with real wave-functions are much faster than with complex ones. However, if a magnetic field is present, if the system is periodic, or if spin-orbit coupling is present, complex wave-functions are mandatory. But don’t worry: the program is able to figure out by itself what to use.
SCF¶
*********************** SCF CYCLE ITER # 1 ************************
etot = -2.54181207E+02 abs_ev = 2.19E-01 rel_ev = 4.04E-03
ediff = 9.34E+00 abs_dens = 2.99E-01 rel_dens = 5.98E-02
Matrix vector products: 25
Converged eigenvectors: 0
# State Eigenvalue [eV] Occupation Error
1 -17.420851 2.000000 ( 4.4E-06)
2 -6.487569 1.000000 ( 1.2E-04)
3 -6.487569 1.000000 ( 1.2E-04)
4 -6.487569 1.000000 ( 1.2E-04)
5 0.946413 0.000000 ( 2.4E-01)
Density of states:
-----------------------------------------%----------------------------
-----------------------------------------%----------------------------
-----------------------------------------%----------------------------
-----------------------------------------%----------------------------
-----------------------------------------%----------------------------
-----------------------------------------%----------------------------
-----------------------------------------%----------------------------
%----------------------------------------%---------------------------%
%----------------------------------------%---------------------------%
%----------------------------------------%---------------------------%
^
Elapsed time for SCF step 1: 0.68
**********************************************************************
Now the SCF cycle starts. For every step, Octopus outputs several pieces of information:
The values
abs_densandrel_densare to monitor the absolute and relative convergence of the density, whilerel_evandabs_evare two alternative measures of the convergence, based on measuring the difference between input and output eigenvalues. The SCF procedure, by default, is stopped whenrel_densis smaller than \(10^{-5}\). This may be altered with the appropriate input variables (see in the manual the variables ConvAbsDens, ConvRelDens, ConvAbsEv and ConvRelEv).The line
Matrix vector products: 25tells us that the Hamiltonian was applied 25 times. This gives us an idea of the computational cost.The line
Converged eigenvectors: 0tells us that upon completion of the diagonalization procedure, none of the orbitals met the required precision criterion for the wavefunctions. In a following example, we will modify this criterion in the input file.The list of eigenvalues is then printed, along with their errors: how much they deviate from “exact” eigenvalues of the current Hamiltonian. This number is the so-called “residue”.
When the calculation is finished, Octopus prints a summary of the run, showing the number of iterations required and the total run time:
Info: Writing states. 2026/02/07 at 03:59:51
Info: Finished writing states. 2026/02/07 at 03:59:51
Info: SCF converged in 14 iterations
Info: Number of matrix-vector products: 350
Info: Finished writing information Ground-state to '/builds/octopus-code/octopus/doc/jupyter_tutorials/1-basics/2-basic_input_options/restart//gs'.
Calculation ended on 2026/02/07 at 03:59:51
Walltime: 05.492s
Octopus emitted 5 warnings.
You can now take a look at the file static/info that will hold a summary of the calculation.
******************************** Grid ********************************
Simulation Box:
Type = sphere
Radius [A] = 5.000
Main mesh:
Spacing [A] = ( 0.180, 0.180, 0.180) volume/point [A^3] = 0.00583
# inner mesh = 89727
# total mesh = 127183
Grid Cutoff [eV] = 1160.586810 Grid Cutoff [Ry] = 85.301565
**********************************************************************
***************************** Symmetries *****************************
Symmetry elements : (i) (Cinf) (sigma)
Symmetry group : Kh
**********************************************************************
**************************** Theory Level ****************************
Input: [TheoryLevel = kohn_sham]
Exchange-correlation:
Exchange
Slater exchange (LDA)
[1] P. A. M. Dirac., Math. Proc. Cambridge Philos. Soc. 26, 376 (1930)
[2] F. Bloch., Z. Phys. 57, 545 (1929)
Correlation
Perdew & Zunger (Modified) (LDA)
[1] J. P. Perdew and A. Zunger., Phys. Rev. B 23, 5048 (1981)
**********************************************************************
SCF converged in 14 iterations
Some of the states are not fully converged!
Eigenvalues [eV]
#st Spin Eigenvalue Occupation
1 -- -18.282871 2.000000
2 -- -7.302321 1.000000
3 -- -7.302321 1.000000
4 -- -7.302321 1.000000
5 -- 0.552496 0.000000
Energy [eV]:
Total = -262.24113162
Free = -262.24113162
-----------
Ion-ion = 0.00000000
Eigenvalues = -58.47270531
Hartree = 221.90846241
Int[n*v_xc] = -76.98672965
Exchange = -51.56021520
Correlation = -7.28647835
vanderWaals = 0.00000000
Delta XC = 0.00000000
Entropy = 4.15888308
-TS = -0.00000000
Photon ex. = 0.00000000
Kinetic = 170.64410402
External = -595.94700450
Non-local = -69.80010537
Int[n*v_E] = 0.00000000
Dipole: [A] [Debye]
<x> = -1.84733E-10 -8.87309E-10
<y> = -3.10102E-11 -1.48948E-10
<z> = -1.06390E-10 -5.11013E-10
Convergence:
abs_energy = 8.44785428E-05 ( 0.00000000E+00) [eV]
rel_energy = 3.22140921E-07 ( 0.00000000E+00)
abs_dens = 1.26515020E-06 ( 0.00000000E+00)
rel_dens = 2.53030040E-07 ( 1.00000000E-06)
abs_evsum = 6.21385844E-05 ( 0.00000000E+00) [eV]
rel_evsum = 1.06269552E-06 ( 0.00000000E+00)
Forces on the ions [eV/A]
Ion x y z
1 N 2.89490183E-09 6.05662767E-10 1.59470546E-09
----------------------------------------------------------
Max abs force 2.89490183E-09 6.05662767E-10 1.59470546E-09
Total force 2.89490183E-09 6.05662767E-10 1.59470546E-09
Total torque 0.00000000E+00 0.00000000E+00 0.00000000E+00
Restarting¶
Any ground-state calculation may be restarted later (to refine it if it did not converge properly, or with any other purpose), provided that the contents of the restart directory are preserved. You can try this now, just by running Octopus again. You will notice that Octopus did not give any warning about missing restart information.
Caution
Octopus will overwrite files written during the first call, so you loose data such as stdout and stderr from the first call.
# extract runtime of first call to compare with restart (stdout_gs.txt will be overwritten)
runtime_first_call, = !grep -am 1 "Walltime" 2-basic_input_options/stdout_gs.txt
!cd 2-basic_input_options && octopus
Info: Ground-state restart information will be read from '/builds/octopus-code/octopus/doc/jupyter_tutorials/1-basics/2-basic_input_options/restart//gs'.
The second run is shorter as it uses the restart files
print(runtime_first_call)
!grep -am 1 "Walltime" 2-basic_input_options/stdout_gs.txt
Walltime: 05.492s
Walltime: 02.946s
This is useful if you change slightly the parameters of the simulation (for example the XC functional or the convergence criteria). If you change the grid parameters Octopus will not be able to restart from the previous calculation. If you do not want Octopus to try to restart a calculation, you can set the variable FromScratch.
In case you were wondering what the restart information looks like, you can have a look at the contents of the restart directory. This is where the files needed to restart a calculation are stored. It may contain several sub-directories depending on the calculations previously performed. In this case, it just contains one:
!ls 2-basic_input_options/restart/
gs
!ls 2-basic_input_options/restart/gs
density df_0103.obf dv_0103.obf indices.obf vhxc wfns
density.obf df_0104.obf dv_0104.obf mixing vhxc.obf
df_0101.obf dv_0101.obf f_old_01.obf occs vin_old_01.obf
df_0102.obf dv_0102.obf grid states wavefunctions.bp
Octopus stores each individual state in a different binary (yet platform-independent) file. Some other useful quantities, like the density, are also stored in binary form. The other files are text files that contain diverse control information. Which files exactly are generated in this folder also might depend on the installation options of Octopus. It is unlikely that you will ever have to work directly with these files, but you may take a look around if you are curious.
Tutorial Validation Checks¶
from postopus import Run, open_textfile
import numpy as np
run = Run("2-basic_input_options")
logfile = open_textfile(run / "stdout_gs.txt")
total_energy = run.scf.info.get_total_energy()
assert(run.scf.info.scf_is_converged())
np.testing.assert_allclose(total_energy, -262.241131, atol=1e-4)
np.testing.assert_equal(int(logfile.grepfield('inner mesh', 4)), 89727)
np.testing.assert_allclose(float(logfile.grepfield('Grid Cutoff', 4)), 1160.586810, atol=1e-4)