Commit 7a008115 by Francois Gygi

Update spectrum command

Add calculation of eigenvalues and eigenfunctions.
Modify command syntax to include optional energy range,
optional gaussian width.
Changed parameter de to const value (0.01).
parent 6ba92e50
......@@ -492,17 +492,6 @@ JDWavefunctionStepper.o: D3vector.h UnitCell.h
KpointCmd.o: UserInterface.h D3vector.h Sample.h AtomSet.h Context.h blacs.h
KpointCmd.o: Atom.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
KpointCmd.o: ExtForceSet.h Wavefunction.h Control.h
LBFGS_IonicStepper.o: LBFGS_IonicStepper.h IonicStepper.h Sample.h AtomSet.h
LBFGS_IonicStepper.o: Context.h blacs.h Atom.h D3vector.h UnitCell.h
LBFGS_IonicStepper.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
LBFGS_IonicStepper.o: Wavefunction.h Control.h Species.h
LBFGS_IonicStepper.o: IonicStepper.h Sample.h AtomSet.h Context.h blacs.h
LBFGS_IonicStepper.o: Atom.h D3vector.h UnitCell.h D3tensor.h blas.h
LBFGS_IonicStepper.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
LBFGS_IonicStepper.o: Species.h
LBFGS_hessian.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
LBFGS_hessian.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
LBFGS_hessian.o: Wavefunction.h Control.h
LDAFunctional.o: LDAFunctional.h XCFunctional.h
LDAFunctional.o: XCFunctional.h
LineMinimizer.o: LineMinimizer.h
......@@ -739,10 +728,13 @@ SpeciesHandler.o: SpeciesHandler.h StructureHandler.h Species.h StrX.h
SpeciesHandler.o: StructureHandler.h
SpeciesReader.o: Species.h SpeciesReader.h StructuredDocumentHandler.h StrX.h
SpeciesReader.o: StructureHandler.h SpeciesHandler.h
SpectrumCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
SpectrumCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
SpectrumCmd.o: ExtForceSet.h Wavefunction.h Control.h MLWFTransform.h
SpectrumCmd.o: BasisMapping.h
SpectrumCmd.o: SpectrumCmd.h UserInterface.h Context.h blacs.h
SpectrumCmd.o: ChargeDensity.h Timer.h D3vector.h EnergyFunctional.h
SpectrumCmd.o: StructureFactor.h ElectricEnthalpy.h Matrix.h Wavefunction.h
SpectrumCmd.o: UnitCell.h SlaterDet.h Basis.h Sample.h AtomSet.h Atom.h
SpectrumCmd.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Control.h
SpectrumCmd.o: MLWFTransform.h BasisMapping.h
SpectrumCmd.o: UserInterface.h
StatusCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
StatusCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
StatusCmd.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h Timer.h
......@@ -850,11 +842,6 @@ qb.o: ThWidth.h Vext.h ExternalPotential.h WfDiag.h WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
release.o: release.h
sinft.o: sinft.h
spectrumCmd.o: SpectrumCmd.h UserInterface.h Sample.h AtomSet.h Context.h
spectrumCmd.o: blacs.h Atom.h D3vector.h UnitCell.h D3tensor.h blas.h
spectrumCmd.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
spectrumCmd.o: MLWFTransform.h BasisMapping.h SlaterDet.h Basis.h Matrix.h
spectrumCmd.o: Timer.h
spline.o: spline.h
testAndersonMixer.o: AndersonMixer.h
testBase64Transcoder.o: Base64Transcoder.h
......@@ -895,7 +882,6 @@ testWavefunction.o: SlaterDet.h Basis.h Matrix.h Timer.h
testXCFunctional.o: LDAFunctional.h XCFunctional.h PBEFunctional.h Timer.h
testXMLGFPreprocessor.o: Context.h blacs.h Matrix.h XMLGFPreprocessor.h
test_fftw.o: Timer.h readTSC.h
test_sym.o: Basis.h D3vector.h UnitCell.h
test_vext.o: Function3d.h D3vector.h
testjacobi.o: Timer.h Context.h blacs.h Matrix.h jacobi.h
testjade.o: Timer.h Context.h blacs.h Matrix.h jade.h
......
......@@ -21,7 +21,9 @@
#include <fstream>
#include <cassert>
#include "Context.h"
#include "SlaterDet.h"
#include "ChargeDensity.h"
#include "EnergyFunctional.h"
#include "MLWFTransform.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
......@@ -40,40 +42,91 @@ int SpectrumCmd::action(int argc, char **argv)
return 1;
}
if ( !( argc == 4 || argc == 6 ) )
if ( !( argc == 2 || argc == 3 || argc == 5 ) )
{
if ( ui->onpe0() )
{
cout << " use: spectrum [emin emax] de width filename" << endl;
cout << " use: spectrum [emin emax] [width] filename" << endl;
}
return 1;
}
// Compute eigenvalues using the current wave function wf
Wavefunction dwf(wf);
ChargeDensity cd(wf);
cd.update_density();
EnergyFunctional ef(*s,cd);
const bool compute_stress = false;
ef.update_vhxc(compute_stress);
const bool compute_forces = false;
const bool compute_hpsi = true;
valarray<double> sigma_eks;
vector<vector<double> > fion;
ef.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
const bool compute_eigvec = true;
wf.diag(dwf,compute_eigvec);
if ( ui->onpe0() )
{
cout << "<eigenset>" << endl;
// print eigenvalues
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
const int nst = wf.sd(ispin,ikp)->nst();
const double eVolt = 2.0 * 13.6058;
cout << " <eigenvalues spin=\"" << ispin
<< "\" kpoint=\""
<< setprecision(8)
<< wf.sd(ispin,ikp)->kpoint()
<< "\" weight=\""
<< setprecision(8)
<< wf.weight(ikp)
<< "\" n=\"" << nst << "\">" << endl;
for ( int i = 0; i < nst; i++ )
{
cout << setw(12) << setprecision(5)
<< wf.sd(ispin,ikp)->eig(i)*eVolt;
if ( i%5 == 4 ) cout << endl;
}
if ( nst%5 != 0 ) cout << endl;
cout << " </eigenvalues>" << endl;
}
}
cout << "</eigenset>" << endl;
}
const double eVolt = 2.0 * 13.6058;
const UnitCell& cell = wf.cell();
// emin, emax: bounds of plot in eV
// de: energy spacing of plot values in eV
// width: gaussian width of convolution in eV
double emin = 0.0, emax = 0.0, de = 0.01, width = 0.05, erange = 0.0;
// de: energy spacing of plot values in eV (fixed at 0.01)
// width: gaussian width of convolution in eV (default 0.05)
const double de = 0.01;
double emin = 0.0, emax = 0.0, width = 0.05, erange = 0.0;
const char *spfilename;
// spectrum width
if ( argc == 4 )
// spectrum file
if ( argc == 2 )
{
spfilename = argv[1];
}
// spectrum width file
if ( argc == 3 )
{
de = atof(argv[1]);
width = atof(argv[2]);
spfilename = argv[3];
width = atof(argv[1]);
spfilename = argv[2];
}
// spectrum emin emax de width
if ( argc == 6 )
// spectrum emin emax de width file
if ( argc == 5 )
{
emin = atof(argv[1]);
emax = atof(argv[2]);
de = atof(argv[3]);
width = atof(argv[4]);
spfilename = argv[5];
width = atof(argv[3]);
spfilename = argv[4];
erange = emax - emin + 3 * width;
if ( emax <= emin )
{
......
......@@ -24,8 +24,7 @@
#include <string>
#include "UserInterface.h"
#include "Sample.h"
#include "MLWFTransform.h"
class Sample;
class SpectrumCmd : public Cmd
{
......@@ -41,19 +40,20 @@ class SpectrumCmd : public Cmd
{
return
"\n spectrum\n\n"
" syntax: spectrum de width filename\n\n"
" spectrum [emin emax] de width filename\n\n"
" The spectrum command computes the dipole transition strength\n"
" between occupied and empty orbitals, following a calculation of\n"
" Kohn-Sham eigenvalues and eigenfunctions. The corresponding absorption\n"
" syntax: spectrum filename\n"
" spectrum width filename\n"
" spectrum emin emax width filename\n\n"
" The spectrum command computes the dipole transition strengths\n"
" between occupied and empty orbitals. It computes Kohn-Sham eigenvalues\n"
" and eigenfunctions of the current wave function using the the current\n"
" value of the xc variable. The corresponding absorption\n"
" spectrum is written on an output file after convolution with a\n"
" gaussian function.\n"
" emin, emax: energy range (optional)\n"
" de : energy spacing of values in the output file\n"
" width : width of the gaussian function used in the convolution\n"
" width : width of the gaussian (optional) (default 0.05 eV)\n"
" filename : output file name\n"
" If the emin, emax parameters are not given, the energy range is defined\n"
" by all possible transitions between occupied and empty orbitals,\n"
" If emin and emax are not given, the energy range includes\n"
" all possible transitions between occupied and empty orbitals.\n"
" All energy parameters must be given in eV.\n\n";
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment