//////////////////////////////////////////////////////////////////////////////// // // Copyright (c) 2008 The Regents of the University of California // // This file is part of Qbox // // Qbox is distributed under the terms of the GNU General Public License // as published by the Free Software Foundation, either version 2 of // the License, or (at your option) any later version. // See the file COPYING in the root directory of this distribution // or . // //////////////////////////////////////////////////////////////////////////////// // // 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. // With a cube of side 1.0 and 32x32x32 points, // 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; 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; 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]; } } } sum = sum * omega / n3; // the density should be normalized cout << " Integrated density: " << sum << endl; LDAFunctional xcf; int nspin = 1; xcf.rho = rh; xcf.vxc1 = vxc; xcf.exc = exc; xcf.setxc(n3,nspin); 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] ); cout << " Total LDA xc energy: " << excsum * omega / n3 << endl; // Note: the energy variation is 3 * dExc/da * delta(a) cout << " dExc/da: " << dxcsum * omega / ( n3 * a ) << endl; return 0; }