testXCFunctional.C 3.44 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
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
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 <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
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.
21
// With a cube of side 1.0 and 32x32x32 points,
Francois Gygi committed
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 <iostream>
#include <vector>
#include "LDAFunctional.h"
#include "PBEFunctional.h"
#include "Timer.h"
#include <cassert>
#include <cmath>
Francois Gygi committed
32
#include <cstdlib>
Francois Gygi committed
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<vector<double> > rh;
  rh.resize(1);
  rh[0].resize(n3);
51

Francois Gygi committed
52 53 54
  XCFunctional *xcf_list[2];
  xcf_list[0] = new LDAFunctional(rh);
  xcf_list[1] = new PBEFunctional(rh);
55

Francois Gygi committed
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;
61

Francois Gygi committed
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];
    }
69

Francois Gygi committed
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;
75

Francois Gygi committed
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];
89

Francois Gygi committed
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();
106

Francois Gygi committed
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];
111

Francois Gygi committed
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] );
121

Francois Gygi committed
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;
}