testLDAFunctional.C 2.36 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 // The xc energy must be -2.8105 a.u. // dExc/da must be 0.911682 #include  Francois Gygi committed Aug 30, 2018 26 #include  Francois Gygi committed Oct 02, 2003 27 28 29 #include "LDAFunctional.h" #include #include  Francois Gygi committed Aug 30, 2018 30 #include  Francois Gygi committed Oct 02, 2003 31 32 33 34 35 36 using namespace std; int main(int argc, char **argv) { if ( argc != 3 ) {  Francois Gygi committed Aug 30, 2018 37 38  cout << " use: testLDAFunctional alat n" << endl; return 1;  Francois Gygi committed Oct 02, 2003 39 40 41 42 43  } double a = atof(argv[1]); double omega = a*a*a; int n = atoi(argv[2]); int n3 = n*n*n;  Francois Gygi committed Aug 30, 2018 44 45  vector > rh(1); rh[0].resize(n3);  Francois Gygi committed Oct 02, 2003 46  double excsum = 0.0, dxcsum = 0.0;  Francois Gygi committed Oct 19, 2007 47   Francois Gygi committed Oct 02, 2003 48 49 50 51  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 52   Francois Gygi committed Oct 02, 2003 53 54 55 56 57 58 59 60 61 62 63  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 );  Francois Gygi committed Aug 30, 2018 64 65  rh[0][ii] = fac * exp( -r2 / (rc*rc) ); sum += rh[0][ii];  Francois Gygi committed Oct 02, 2003 66 67 68  } } }  Francois Gygi committed Oct 19, 2007 69  sum = sum * omega / n3;  Francois Gygi committed Oct 02, 2003 70  // the density should be normalized  Francois Gygi committed Oct 19, 2007 71 72  cout << " Integrated density: " << sum << endl;  Francois Gygi committed Aug 30, 2018 73 74  LDAFunctional xcf(rh); xcf.setxc();  Francois Gygi committed Oct 19, 2007 75   Francois Gygi committed Oct 02, 2003 76  for ( int i = 0; i < n3; i++ )  Francois Gygi committed Aug 30, 2018 77  excsum += xcf.rho[i] * xcf.exc[i];  Francois Gygi committed Oct 02, 2003 78  for ( int i = 0; i < n3; i++ )  Francois Gygi committed Aug 30, 2018 79  dxcsum += xcf.rho[i] * ( xcf.exc[i] - xcf.vxc1[i] );  Francois Gygi committed Oct 19, 2007 80   Francois Gygi committed Oct 02, 2003 81  cout << " Total LDA xc energy: " << excsum * omega / n3 << endl;  Francois Gygi committed Oct 19, 2007 82   Francois Gygi committed Oct 02, 2003 83 84 85 86 87  // Note: the energy variation is 3 * dExc/da * delta(a) cout << " dExc/da: " << dxcsum * omega / ( n3 * a ) << endl; return 0; }