testLDAFunctional.C 2.46 KB
 Francois Gygi committed Oct 02, 2003 1 2 //////////////////////////////////////////////////////////////////////////////// //  Francois Gygi committed Aug 13, 2008 3 4 5 6 // Copyright (c) 2008 The Regents of the University of California // // This file is part of Qbox //  Francois Gygi committed Sep 08, 2008 7 8 // Qbox is distributed under the terms of the GNU General Public License // as published by the Free Software Foundation, either version 2 of  Francois Gygi committed Aug 13, 2008 9 10 11 12 13 14 // the License, or (at your option) any later version. // See the file COPYING in the root directory of this distribution // or . // //////////////////////////////////////////////////////////////////////////////// //  Francois Gygi committed Oct 02, 2003 15 16 17 18 19 20 // testLDAFunctional.C // //////////////////////////////////////////////////////////////////////////////// // Test the LDA functional by computing the xc energy of a gaussian // of width 0.1 a.u. in a cube of side 1.0 a.u.  Francois Gygi committed Oct 19, 2007 21 // With a cube of side 1.0 and 32x32x32 points,  Francois Gygi committed Oct 02, 2003 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 // The xc energy must be -2.8105 a.u. // dExc/da must be 0.911682 #include #include "LDAFunctional.h" #include #include using namespace std; int main(int argc, char **argv) { // use: testxcf alat np if ( argc != 3 ) { cout << " use: testLDAFunctional alat np" << endl; return 0; } assert(argc==3); double a = atof(argv[1]); double omega = a*a*a; int n = atoi(argv[2]); int n3 = n*n*n; double *rh = new double[n3]; double *vxc = new double[n3]; double *exc = new double[n3]; double excsum = 0.0, dxcsum = 0.0;  Francois Gygi committed Oct 19, 2007 48   Francois Gygi committed Oct 02, 2003 49 50 51 52  double rc = 0.1 * a; double pi = 4.0 * atan(1.0); double fac = 1.0 / ( pow(pi,1.5) * rc*rc*rc ); double sum = 0.0;  Francois Gygi committed Oct 19, 2007 53   Francois Gygi committed Oct 02, 2003 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69  for ( int i = 0; i < n; i++ ) { double x = ( i * a ) / n - a/2; for ( int j = 0; j < n; j++ ) { double y = ( j * a ) / n - a/2; for ( int k = 0; k < n; k++ ) { double z = ( k * a ) / n - a/2; double r2 = x*x + y*y + z*z; int ii = i + n * ( j + n * k ); rh[ii] = fac * exp( -r2 / (rc*rc) ); sum += rh[ii]; } } }  Francois Gygi committed Oct 19, 2007 70  sum = sum * omega / n3;  Francois Gygi committed Oct 02, 2003 71  // the density should be normalized  Francois Gygi committed Oct 19, 2007 72 73  cout << " Integrated density: " << sum << endl;  Francois Gygi committed Oct 02, 2003 74  LDAFunctional xcf;  Francois Gygi committed Oct 19, 2007 75   Francois Gygi committed Oct 02, 2003 76 77 78 79  int nspin = 1; xcf.rho = rh; xcf.vxc1 = vxc; xcf.exc = exc;  Francois Gygi committed Oct 19, 2007 80   Francois Gygi committed Oct 02, 2003 81  xcf.setxc(n3,nspin);  Francois Gygi committed Oct 19, 2007 82   Francois Gygi committed Oct 02, 2003 83 84 85 86  for ( int i = 0; i < n3; i++ ) excsum += rh[i] * exc[i]; for ( int i = 0; i < n3; i++ ) dxcsum += rh[i] * ( exc[i] - vxc[i] );  Francois Gygi committed Oct 19, 2007 87   Francois Gygi committed Oct 02, 2003 88  cout << " Total LDA xc energy: " << excsum * omega / n3 << endl;  Francois Gygi committed Oct 19, 2007 89   Francois Gygi committed Oct 02, 2003 90 91 92 93 94  // Note: the energy variation is 3 * dExc/da * delta(a) cout << " dExc/da: " << dxcsum * omega / ( n3 * a ) << endl; return 0; }