1/** \page Outputs Output variables
2
3\image html tutorials/Orbitals.jpg
4
5The previous tutorial introduced a basic calculation of a water molecule,
6and visualized one possible output: the electron density.
7This tutorial explores other possible outputs from the same calculation.
8It also discusses how to access the binary output files from Octave/MATLAB or Python/NumPy,
9which you can skip till later if/when you need to do custom postprocessing on JDFTx outputs.
10
11Here's the same input file water.in as before,
12with the introductory comments removed,
13and with more output options requested:
14
15    lattice Cubic 10  #Alternate to explicitly specifying lattice vectors
16    elec-cutoff 20 100
17    ion-species GBRV/h_pbe.uspp
18    ion-species GBRV/o_pbe.uspp
19    coords-type cartesian
20    ion O   0.00 0.00  0.00  0
21    ion H   0.00 1.13 +1.45  1
22    ion H   0.00 1.13 -1.45  1
23
24    dump-name water.$VAR              #Filename pattern for outputs
25
26    #Outputs every few electronic steps:
27    dump Electronic State        #The first parameter is when to dump, the remaining are what
28    dump-interval Electronic 5   #How often to dump; this requests every 5 electronic steps
29
30    #Output at the end:
31    dump End Ecomponents ElecDensity EigStats Vscloc Dtot RealSpaceWfns
32
33Here we invoke the \ref CommandDump command multiple times to request
34different quantities to be output at various stages of the calculation.
35The \ref CommandDumpInterval controls the interval at each type
36of calculation stage where output can be retrieved.
37
38Now run
39
40    jdftx -i water.in | tee water.out
41
42and examine water.out.
43Every 5 electronic steps, the calculation outputs State, which is
44the collection of all variables necessary to continue the calculation.
45In this case, that consists of wfns, the electronic wavefunctions (more precisely, Kohn-Sham orbitals)
46as can be seen in the repeated lines 'Dumping water.wfns ... done.'
47
48After completing electronic minimization, the calculation dumps
49
50Filename           | Description
51-------------------|------------------------------------
52water.Ecomponents  | Ecomponents (energy components)
53water.eigStats     | EigStats (eigenvalue statistics)
54water.n            | ElecDensity (electron density)
55water.d_tot        | Dtot (which is the net electrostatic potential)
56water.Vscloc       | Vscloc (which is the self-consistent Kohn-Sham potential)
57water.wfns_0_x.rs  | with x=0,1,2,3 for RealSpaceWfns (real-space wave functions)
58
59Examine water.Ecomponents.
60It contains the energy components printed in exactly the same format as the output file.
61The energy components reported in this case are
62
63Component    | Description
64-------------|------------------------------------
65Eewald       | %Ewald: nuclear-nuclear %Coulomb energy
66EH           | Hartree: electron-electron mean-field %Coulomb energy
67Eloc         | Electron-nuclear interaction (local part of pseudopotential)
68Enl          | Electron-nuclear interaction (nonlocal part of pseudopotential)
69Exc          | Exchange-Correlation energy
70Exc_core     | Exchange-Correlation energy subtraction for core electrons
71KE           | Kinetic energy of electrons
72Etot         | Total energy
73
74Next look at water.eigStats (same content also reproduced in output file).
75It contains statistics of the Kohn-Sham eigenvalues
76such as the minimum and maximum, and the highest-occupied orbital (HOMO) energies.
77Some of the energies, such as the lowest-unoccupied orbital (LUMO) energy
78and electron chemical potential (mu) are reported as unavailable
79since the present calculation did not include any unoccupied states.
80
81The remaining outputs are binary files containing the values of
82the corresponding quantities on a uniform grid,
83in this case on a 40 x 40 x 40 grid for the wavefunctions and
84a 48 x 48 x 48 grid for the remainder (see the lines starting with
85'Chosen fftbox' in the initialization section of the output file).
86
87Within the [Octave](http://www.octave.org) or MATLAB commandline,
88we can examine the binary electron density file which contains
89the electron density n(r) at each grid point
90(in electrons/bohr<sup>3</sup> units) using:
91
92     fp = fopen('water.n', 'rb');  %# read as a binary file
93     n = fread(fp, 'double');      %# stored as double-precision floating point numbers
94     %# n = swapbytes(n);          %# uncomment on big-endian machines; data is always stored litte-endian
95     fclose(fp);
96
97     S = [ 48 48 48 ];             %# grid dimensions (see output file)
98     V = 1000;                     %# unit cell volume (see output file)
99     dV = V / prod(S);
100     Nelectrons = dV * sum(n)
101
102     n = reshape(n,S)(:,:,1);      %# convert to 3D and extract yz slice
103     n = circshift(n, S(1:2)/2);   %# shift molecule from origin to box center
104     contour(n)
105
106This should print Nelectrons = 8, the number of valence electrons in the calculation
107and then present a contour plot of the yz slice of the electron density.
108Note that Octave/MATLAB differ from C++ in the array ordering convention,
109so the first array dimension is z, the second is y and the final one is x.
110The circshift prevents the tearing of the molecule across the boundaries.
111(The big-endian fix only applies if you are working on a PowerPC architecture,
112such as the Bluegene supercomputer and very old Macs;
113this is very unlikely.)
114
115Alternately, if you prefer Python/NumPy, you can achieve exactly the same result with:
116
117     import numpy as np
118     import matplotlib.pyplot as plt
119
120     n = np.fromfile('water.n',   # read from binary file
121               dtype=np.float64)  # as 64-bit (double-precision) floats
122     #n = swapbytes(n)            # uncomment on big-endian machines; data is always stored litte-endian
123
124     S = [48,48,48]               # grid dimensions (see output file)
125     V = 1000.0                   # unit cell volume (see output file)
126     dV = V / np.prod(S)
127     print "Nelectrons = ", dV * np.sum(n)
128
129     n = n.reshape(S)[0,:,:]        # convert to 3D and extract yz slice
130     n = np.roll(n, S[1]/2, axis=0) # shift molecule from origin to box center (along y)
131     n = np.roll(n, S[2]/2, axis=1) # shift molecule from origin to box center (along z)
132     plt.contour(n)
133     plt.show()
134
135The only difference here is that the default data order in NumPy is the same as C++.
136
137The above scripts show you how you can directly access the
138binary output data for custom post-processing and plotting.
139However, for simple visualization tasks, the createXSF script should suffice,
140as we demonstrated with the electron density at the end of the [previous tutorial](FirstCalc.html).
141Now visualize the Kohn-Sham orbitals using:
142
143    createXSF water.out water-x.xsf water.wfns_0_x.rs
144
145for each x in 0,1,2,3, and view each of the generated XSF files
146using VESTA to get the plots shown below and at the start of this tutorial.
147(You will need to change the boundary setting as in the previous tutorial.)
148
149The orbitals are sorted by their corresponding energy eigenvalue.
150In VESTA, the yellow and cyan correspond to different signs of the plotted quantity.
151Note that the lowest orbital has the same sign thoughout,
152whereas the remaining three have symmetric positive and negative regions
153separated by a nodal plane (where the orbital is zero).
154Two of the three nodal planes are exact planes constrained by reflection symmetry
155of the molecule, while the third is only approximately a plane. (Which one is that?)
156
157\image html tutorials/Orbitals.jpg
158
159*/
160