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