testXCFunctional.C 3.44 KB
 Francois Gygi committed May 20, 2004 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 May 20, 2004 15 16 17 18 19 20 // testXCFunctional.C // //////////////////////////////////////////////////////////////////////////////// // Test an XC 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 May 20, 2004 22 23 24 25 26 27 28 29 30 31 // The LDA xc energy must be -2.8105 a.u. // dExc/da must be 0.911682 #include #include #include "LDAFunctional.h" #include "PBEFunctional.h" #include "Timer.h" #include #include  Francois Gygi committed Aug 30, 2018 32 #include  Francois Gygi committed May 20, 2004 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 using namespace std; int main(int argc, char **argv) { // use: testxcf alat np if ( argc != 3 ) { cout << " use: testXCFunctional 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; vector > rh; rh.resize(1); rh[0].resize(n3);  Francois Gygi committed Oct 19, 2007 51   Francois Gygi committed May 20, 2004 52 53 54  XCFunctional *xcf_list[2]; xcf_list[0] = new LDAFunctional(rh); xcf_list[1] = new PBEFunctional(rh);  Francois Gygi committed Oct 19, 2007 55   Francois Gygi committed May 20, 2004 56 57 58 59 60  for ( int ixcf = 0; ixcf < 2; ixcf++ ) { Timer tm; XCFunctional *xcf = xcf_list[ixcf]; cout << endl << " Functional name: " << xcf->name() << endl;  Francois Gygi committed Oct 19, 2007 61   Francois Gygi committed May 20, 2004 62 63 64 65 66 67 68  double *grad_rho[3]; if ( xcf->isGGA() ) { grad_rho[0] = xcf->grad_rho[0]; grad_rho[1] = xcf->grad_rho[1]; grad_rho[2] = xcf->grad_rho[2]; }  Francois Gygi committed Oct 19, 2007 69   Francois Gygi committed May 20, 2004 70 71 72 73 74  const double rc = 0.1 * a; const double rc2 = rc * rc; const double pi = M_PI; const double fac = 1.0 / ( pow(pi,1.5) * rc*rc*rc ); double sum = 0.0;  Francois Gygi committed Oct 19, 2007 75   Francois Gygi committed May 20, 2004 76 77 78 79 80 81 82 83 84 85 86 87 88  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[0][ii] = fac * exp( -r2 / rc2 ); sum += rh[0][ii];  Francois Gygi committed Oct 19, 2007 89   Francois Gygi committed May 20, 2004 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105  if ( xcf->isGGA() ) { grad_rho[0][ii] = - rh[0][ii] * 2.0 * x / rc2; grad_rho[1][ii] = - rh[0][ii] * 2.0 * y / rc2; grad_rho[2][ii] = - rh[0][ii] * 2.0 * z / rc2; } } } } sum = sum * omega / n3; // the density should be normalized cout << " Integrated density: " << sum << endl; tm.start(); xcf->setxc(); tm.stop();  Francois Gygi committed Oct 19, 2007 106   Francois Gygi committed May 20, 2004 107 108 109 110  double *exc = xcf->exc; double excsum = 0.0; for ( int i = 0; i < n3; i++ ) excsum += rh[0][i] * exc[i];  Francois Gygi committed Oct 19, 2007 111   Francois Gygi committed May 20, 2004 112 113 114 115 116 117 118 119 120  cout << " Total xc energy: " << excsum * omega / n3 << endl; // if not a GGA, compute dExc/da if ( !xcf->isGGA() ) { double *vxc1 = xcf->vxc1; double dxcsum = 0.0; for ( int i = 0; i < n3; i++ ) dxcsum += rh[0][i] * ( exc[i] - vxc1[i] );  Francois Gygi committed Oct 19, 2007 121   Francois Gygi committed May 20, 2004 122 123 124 125 126 127 128  // Note: the energy variation is 3 * dExc/da * delta(a) cout << " dExc/da: " << dxcsum * omega / ( n3 * a ) << endl; } cout << " " << xcf->name() << " time: " << tm.real() << endl; } // ixcf return 0; }