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)35 int 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