1 //////////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2008 The Regents of the University of California 4 // 5 // This file is part of Qbox 6 // 7 // Qbox is distributed under the terms of the GNU General Public License 8 // as published by the Free Software Foundation, either version 2 of 9 // the License, or (at your option) any later version. 10 // See the file COPYING in the root directory of this distribution 11 // or <http://www.gnu.org/licenses/>. 12 // 13 //////////////////////////////////////////////////////////////////////////////// 14 // 15 // testXCFunctional.cpp 16 // 17 //////////////////////////////////////////////////////////////////////////////// 18 19 // Test an XC functional by computing the xc energy of a gaussian 20 // of width 0.1 a.u. in a cube of side 1.0 a.u. 21 // With a cube of side 1.0 and 32x32x32 points, 22 // The LDA xc energy must be -2.8105 a.u. 23 // dExc/da must be 0.911682 24 25 #include <iostream> 26 #include <vector> 27 #include "LDAFunctional.h" 28 #include "PBEFunctional.h" 29 #include "Timer.h" 30 #include <cassert> 31 #include <cmath> 32 #include <cstdlib> 33 using namespace std; 34 main(int argc,char ** argv)35int main(int argc, char **argv) 36 { 37 // use: testxcf alat np 38 if ( argc != 3 ) 39 { 40 cout << " use: testXCFunctional alat np" << endl; 41 return 0; 42 } 43 assert(argc==3); 44 double a = atof(argv[1]); 45 double omega = a*a*a; 46 int n = atoi(argv[2]); 47 int n3 = n*n*n; 48 vector<vector<double> > rh; 49 rh.resize(1); 50 rh[0].resize(n3); 51 52 XCFunctional *xcf_list[2]; 53 xcf_list[0] = new LDAFunctional(rh); 54 xcf_list[1] = new PBEFunctional(rh); 55 56 for ( int ixcf = 0; ixcf < 2; ixcf++ ) 57 { 58 Timer tm; 59 XCFunctional *xcf = xcf_list[ixcf]; 60 cout << endl << " Functional name: " << xcf->name() << endl; 61 62 double *grad_rho[3]; 63 if ( xcf->isGGA() ) 64 { 65 grad_rho[0] = xcf->grad_rho[0]; 66 grad_rho[1] = xcf->grad_rho[1]; 67 grad_rho[2] = xcf->grad_rho[2]; 68 } 69 70 const double rc = 0.1 * a; 71 const double rc2 = rc * rc; 72 const double pi = M_PI; 73 const double fac = 1.0 / ( pow(pi,1.5) * rc*rc*rc ); 74 double sum = 0.0; 75 76 for ( int i = 0; i < n; i++ ) 77 { 78 double x = ( i * a ) / n - a/2; 79 for ( int j = 0; j < n; j++ ) 80 { 81 double y = ( j * a ) / n - a/2; 82 for ( int k = 0; k < n; k++ ) 83 { 84 double z = ( k * a ) / n - a/2; 85 double r2 = x*x + y*y + z*z; 86 int ii = i + n * ( j + n * k ); 87 rh[0][ii] = fac * exp( -r2 / rc2 ); 88 sum += rh[0][ii]; 89 90 if ( xcf->isGGA() ) 91 { 92 grad_rho[0][ii] = - rh[0][ii] * 2.0 * x / rc2; 93 grad_rho[1][ii] = - rh[0][ii] * 2.0 * y / rc2; 94 grad_rho[2][ii] = - rh[0][ii] * 2.0 * z / rc2; 95 } 96 } 97 } 98 } 99 sum = sum * omega / n3; 100 // the density should be normalized 101 cout << " Integrated density: " << sum << endl; 102 103 tm.start(); 104 xcf->setxc(); 105 tm.stop(); 106 107 double *exc = xcf->exc; 108 double excsum = 0.0; 109 for ( int i = 0; i < n3; i++ ) 110 excsum += rh[0][i] * exc[i]; 111 112 cout << " Total xc energy: " << excsum * omega / n3 << endl; 113 114 // if not a GGA, compute dExc/da 115 if ( !xcf->isGGA() ) 116 { 117 double *vxc1 = xcf->vxc1; 118 double dxcsum = 0.0; 119 for ( int i = 0; i < n3; i++ ) 120 dxcsum += rh[0][i] * ( exc[i] - vxc1[i] ); 121 122 // Note: the energy variation is 3 * dExc/da * delta(a) 123 cout << " dExc/da: " << dxcsum * omega / ( n3 * a ) << endl; 124 } 125 cout << " " << xcf->name() << " time: " << tm.real() << endl; 126 } // ixcf 127 return 0; 128 } 129