working commit

applications/allElectron/carbon/main.cc

@@ -21,7 +21,7 @@ const double kb = 3.166811429e-06;
 //#define meshFileName "../../../data/mesh/singleAtom.inp"
 
 //dft header
-#include "../../../src/dft/dft.cc"
+#include "../../../src2/dft/dft.cc"
 
 //Atom locations
 void getAtomicLocations(dealii::Table<2,double>& atoms, std::map<unsigned int, std::string>& initialGuessFiles){

applications/allElectron/carbon/temp

@@ -0,0 +1,97 @@
+[100%] Built target main
+[ 50%] Built target main
+[100%] Run main with Release configuration
+number of MPI processes: 1
+Refinement levels executed: 9
+Smallest element size: 3.91e-02
+ Number of global active cells: 512
+number of elements: 512
+number of degrees of freedom: 41265
+reading initial guess for rho
+initial total charge: 5.9999937449e+00 
+numAtoms: 1
+Atom core (6) located with node id 37036 in processor 0 and added 
+massVector norm: 82106.3
+reading initial guess for PSI
+poisson solve: initial residual:9.897861e+00, current residual:9.013215e-12, nsteps:1454, tolerance criterion:9.897861e-12
+
+
+Begin SCF Iteration:1
+poisson solve: initial residual:6.000385e+00, current residual:5.548787e-12, nsteps:1450, tolerance criterion:6.000385e-12
+bUp: 1.8334349407e+05
+bLow: 0
+a0: -10
+ 0 l2: 2.2144004761e+00 linf: 1.5729174242e-01 
+ 1 l2: 2.2479967764e+00 linf: 1.9821693451e-01 
+ 2 l2: 2.2507174743e+00 linf: 2.4823420842e-01 
+ 3 l2: 2.2507174743e+00 linf: 2.4823420842e-01 
+ 4 l2: 2.2507174743e+00 linf: 2.4823420842e-01 
+Total time for only chebyshev filter: 1.25842mins
+ 0 l2: 1.3920183398e+00 linf: 3.5493186433e-02 
+ 1 l2: 4.4307214332e-03 linf: 1.1297295180e-04 
+ 2 l2: 1.4594461733e-11 linf: 6.3949892641e-13 
+ 3 l2: 1.4594461745e-11 linf: 6.3949906741e-13 
+ 4 l2: 1.4594461722e-11 linf: 6.3949863441e-13 
+eigen value 0: -9.9468289835e+00
+eigen value 1: -4.9876707353e-01
+eigen value 2: -7.4697089272e-02
+eigen value 3: -7.4697089272e-02
+eigen value 4: -7.4697089272e-02
+Total time for chebyshev filter: 1.25915mins
+Fermi energy Residual: 6.21724893790087662637e-15 
+Fermi energy: -7.57946224788558914343e-02 
+poisson solve: initial residual:6.000385e+00, current residual:5.806451e-12, nsteps:1459, tolerance criterion:6.000385e-12
+partialOccupancy 0: 1.00000000000000000000e+00 
+partialOccupancy 1: 1.00000000000000000000e+00 
+partialOccupancy 2: 3.33333333335232628869e-01 
+partialOccupancy 3: 3.33333333333678427657e-01 
+partialOccupancy 4: 3.33333333331091941076e-01 
+Total energy: -3.71716096143394452156e+01 
+Band energy: -2.10405862925250097817e+01 
+Kinetic energy: 3.74341842907320341283e+01 
+Exchange energy: -4.35485728308278297050e+00 
+Correlation energy: -3.71312296237106986840e-01 
+Electrostatic energy: -6.98796243257515925507e+01 
+Repulsive energy: 0.00000000000000000000e+00 
+SCF iteration: 1 complete
+
+
+Begin SCF Iteration:2
+Mixing Scheme: iter:2, norm:1.769735e-04
+poisson solve: initial residual:6.000385e+00, current residual:5.994737e-12, nsteps:1458, tolerance criterion:6.000385e-12
+bUp: 1.8334349415e+05
+bLow: -0.0746971
+a0: -9.94683
+ 0 l2: 1.0000000000e+00 linf: 2.5497642841e-02 
+ 1 l2: 1.0000000000e+00 linf: 3.0147484616e-02 
+ 2 l2: 1.0000000000e+00 linf: 4.5711959505e-02 
+ 3 l2: 1.0000000000e+00 linf: 4.6756352302e-02 
+ 4 l2: 1.0000000000e+00 linf: 4.5698780592e-02 
+Total time for only chebyshev filter: 1.24024mins
+ 0 l2: 1.0023233770e+00 linf: 2.5557062659e-02 
+ 1 l2: 1.0466293154e-05 linf: 2.6686803551e-07 
+ 2 l2: 1.0590253907e-10 linf: 2.7172682491e-12 
+ 3 l2: 5.3163123015e-11 linf: 1.3701052649e-12 
+ 4 l2: 4.2536019674e-10 linf: 1.0864862631e-11 
+eigen value 0: -9.9483901631e+00
+eigen value 1: -5.0112644342e-01
+eigen value 2: -1.9892759180e-01
+eigen value 3: -1.9892759151e-01
+eigen value 4: -1.9892759126e-01
+Total time for chebyshev filter: 1.24093mins
+Fermi energy Residual: 8.88178419700125232339e-16 
+Fermi energy: -2.00025124730548176100e-01 
+poisson solve: initial residual:6.000385e+00, current residual:5.552124e-12, nsteps:1465, tolerance criterion:6.000385e-12
+partialOccupancy 0: 1.00000000000000000000e+00 
+partialOccupancy 1: 1.00000000000000000000e+00 
+partialOccupancy 2: 3.33333372348329382007e-01 
+partialOccupancy 3: 3.33333331883818750896e-01 
+partialOccupancy 4: 3.33333295767851256475e-01 
+Total energy: -3.74227688801870144175e+01 
+Band energy: -2.12968883960854711290e+01 
+Kinetic energy: 3.71902845108251653983e+01 
+Exchange energy: -4.35554543191728260609e+00 
+Correlation energy: -3.71353313545927321560e-01 
+Electrostatic energy: -6.98861546455489701657e+01 
+Repulsive energy: 0.00000000000000000000e+00 
+SCF iteration: 2 complete

include2/dft.h

@@ -0,0 +1,102 @@
+#ifndef dft_H_
+#define dft_H_
+#include <iostream>
+#include <iomanip> 
+#include <numeric>
+#include <sstream> 
+#include "headers.h"
+#include "poisson.h"
+#include "eigen.h"
+//include alglib
+//#include "/nfs/mcfs_home/rudraa/Public/alglib/cpp/src/interpolation.h"
+//#include "/nfs/mcfs_home/rudraa/Public/libxc/libxc-2.2.0/installDir/include/xc.h"
+#include "/opt/software/numerics/alglib/cpp/src/interpolation.h"
+#include "/opt/software/numerics/libxc-2.2.2/installDir/include/xc.h"
+
+//std::cout << std::setprecision(18) << std::scientific;
+
+//Initialize Namespace
+using namespace dealii;
+//blas-lapack routines
+extern "C"{
+ void dgemv_(char* TRANS, const int* M, const int* N, double* alpha, double* A, const int* LDA, double* X, const int* INCX, double* beta, double* C, const int* INCY);
+ void dgesv_( int* n, int* nrhs, double* a, int* lda, int* ipiv, double* b, int* ldb, int* info );
+ void dscal_(int *n, double *alpha, double *x, int *incx);
+ void daxpy_(int *n, double *alpha, double *x, int *incx, double *y, int *incy);
+ void dgemm_(char* transA, char* transB, int *m, int *n, int *k, double *alpha, double *A, int *lda, double *B, int *ldb, double *beta, double *C, int *ldc);
+ void dsyevd_(char* jobz, char* uplo, int* n, double* A, int *lda, double* w, double* work, int* lwork, int* iwork, int* liwork, int* info);
+}
+xc_func_type funcX, funcC;
+
+//Define dft class
+class dftClass{
+ friend class poissonClass;
+ friend class eigenClass; 
+ public:
+ dftClass();
+ void run();
+ Table<2,double> atomLocations;
+ std::map<unsigned int, std::string> initialGuessFiles;
+
+ private:
+ void mesh();
+ void init();
+ void initRho();
+ double totalCharge();
+ void locateAtomCoreNodes();
+ double mixing_simple();
+ double mixing_anderson();
+ void compute_energy();
+ void compute_fermienergy();
+ double repulsiveEnergy();
+ void compute_rhoOut();
+ void readPSIRadialValues(std::vector<std::vector<std::vector<double> > >& singleAtomPSI);
+ void readPSI();
+
+ //FE data structres
+ parallel::distributed::Triangulation<3> triangulation;
+ FE_Q<3> FE;
+ DoFHandler<3> dofHandler;
+ MatrixFree<3,double> matrix_free_data, matrix_free_dataGauss;
+
+ //parallel objects
+ MPI_Comm mpi_communicator;
+ const unsigned int n_mpi_processes;
+ const unsigned int this_mpi_process;
+ IndexSet locally_owned_dofs;
+ IndexSet locally_relevant_dofs;
+
+ poissonClass poisson;
+ eigenClass eigen;
+ ConstraintMatrix constraintsNone;
+ std::vector<double> eigenValues;
+ std::vector<parallel::distributed::Vector<double>*> eigenVectors;
+
+ //parallel message stream
+ ConditionalOStream pcout; 
+
+ //compute-time logger
+ TimerOutput computing_timer;
+
+ //dft related objects
+ std::map<dealii::CellId, std::vector<double> > *rhoInValues, *rhoOutValues;
+ std::vector<std::map<dealii::CellId,std::vector<double> >*> rhoInVals, rhoOutVals;
+
+ //map of atom node number and atomic weight
+ std::map<unsigned int, double> atoms; 
+
+ //fermi energy
+ double fermiEnergy;
+
+ //chebyshev filter variables and functions
+ double bUp, bLow, a0;
+ vectorType vChebyshev, v0Chebyshev, fChebyshev;
+ std::vector<parallel::distributed::Vector<double>*> PSI, tempPSI, tempPSI2, tempPSI3;
+ void chebyshevSolver();
+ double upperBound();
+ void gramSchmidt(std::vector<vectorType*>& X);
+ void chebyshevFilter(std::vector<vectorType*>& X, unsigned int m, double a, double b, double a0); 
+ void rayleighRitz(std::vector<vectorType*>& X);
+};
+
+#endif

include2/eigen.h

@@ -0,0 +1,55 @@
+#ifndef eigen_H_
+#define eigen_H_
+#include "headers.h"
+
+//Define eigenClass class
+class eigenClass
+{
+ friend class dftClass;
+public:
+ eigenClass(dftClass* _dftPtr);
+ void HX(const std::vector<vectorType*> &src, std::vector<vectorType*> &dst);
+ void XHX(const std::vector<vectorType*> &src); 
+ void vmult (vectorType &dst, const vectorType &src) const;
+ private:
+ void implementHX(const dealii::MatrixFree<3,double> &data,
+ std::vector<vectorType*> &dst, 
+ const std::vector<vectorType*> &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const;
+ void init ();
+ void computeMassVector();
+ void computeVEffectiveRHS(std::map<dealii::CellId,std::vector<double> >* rhoValues, const vectorType& phi); 
+ void MX (const dealii::MatrixFree<3,double> &data,
+ vectorType &dst,
+ const vectorType &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const;
+ void computeVEffective(std::map<dealii::CellId,std::vector<double> >* rhoValues, const vectorType& phi);
+
+ //pointer to dft class
+ dftClass* dftPtr;
+
+ //FE data structres
+ dealii::FE_Q<3> FE;
+
+ //constraints
+ dealii::ConstraintMatrix constraintsNone; 
+
+ //data structures
+ vectorType massVector, vEffective, rhsVeff;
+ std::vector<double> XHXValue;
+ std::vector<double>* XHXValuePtr;
+ std::vector<vectorType*> HXvalue;
+
+ //parallel objects
+ MPI_Comm mpi_communicator;
+ const unsigned int n_mpi_processes;
+ const unsigned int this_mpi_process;
+ dealii::ConditionalOStream pcout;
+
+ //compute-time logger
+ dealii::TimerOutput computing_timer;
+ //mutex thread for managing multi-thread writing to XHXvalue
+ mutable dealii::Threads::Mutex assembler_lock;
+};
+
+#endif

include2/headers.h

@@ -0,0 +1,51 @@
+#ifndef headers_H_
+#define headers_H_
+
+//Include all deal.II header file
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/table.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+//Include generic C++ headers
+#include <fstream>
+#include <iostream>
+
+#endif

include2/poisson.h

@@ -0,0 +1,58 @@
+#ifndef poisson_H_
+#define poisson_H_
+#include "headers.h"
+
+typedef double dataType;
+typedef dealii::parallel::distributed::Vector<double> vectorType;
+
+class dftClass;
+
+//Define poisson class
+class poissonClass
+{
+ friend class dftClass; 
+public:
+ poissonClass(dftClass* _dftPtr);
+ void computeLocalJacobians();
+ void vmult(vectorType &dst, const vectorType &src) const;
+ void precondition_Jacobi(vectorType& dst, const vectorType& src, const double omega) const;
+ void unsubscribe (const char *identifier=0) const{}; //function needed to mimic SparseMatrix for Jacobi Preconditioning
+ bool operator!= (double val) const {return true;}; //function needed to mimic SparseMatrix
+private:
+ void init ();
+ void computeRHS(std::map<dealii::CellId,std::vector<double> >* rhoValues);
+ void solve(vectorType& phi, std::map<dealii::CellId,std::vector<double> >* rhoValues=0);
+ void AX(const dealii::MatrixFree<3,double> &data,
+ vectorType &dst, 
+ const vectorType &src,
+ const std::pair<unsigned int,unsigned int> &cell_range) const;
+ //pointer to dft class
+ dftClass* dftPtr;
+
+ //FE data structres
+ dealii::FE_Q<3> FE;
+ std::map<dealii::CellId,std::vector<double> > localJacobians;
+ std::map<dealii::CellId,std::vector<double> >* localJacobiansPtr; //this ptr created to circumvent problem with const definition of vmult and Ax
+ //constraints
+ dealii::ConstraintMatrix constraintsNone, constraintsZero, constraints1byR;
+ std::map<dealii::types::global_dof_index, double> valuesZero, values1byR;
+
+ //data structures
+ vectorType rhs, Ax;
+ vectorType jacobianDiagonal;
+ vectorType phiTotRhoIn, phiTotRhoOut, phiExt;
+ vectorType rhs2;
+ double jacobianDiagonalValue;
+ double relaxation; //relaxation parameter for Jacobi Preconditioning
+
+ //parallel objects
+ MPI_Comm mpi_communicator;
+ const unsigned int n_mpi_processes;
+ const unsigned int this_mpi_process;
+ dealii::ConditionalOStream pcout;
+
+ //compute-time logger
+ dealii::TimerOutput computing_timer;
+};
+
+#endif

src/dft/mesh.cc

@@ -85,7 +85,7 @@ void dftClass::mesh(){
  }
  }
 
-
+ /*
  cell = triangulation.begin_active(), end_cell = triangulation.end();
  double hmin=L, Lmin=L;
  for ( ; cell != end_cell; ++cell){
@@ -99,7 +99,7 @@ void dftClass::mesh(){
  }
  }
  triangulation.execute_coarsening_and_refinement();
-
+ */
  /*
  //define mesh parameters
  double L=20;