testLDAFunctional.C 2.36 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
// 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.
21
// With a cube of side 1.0 and 32x32x32 points,
Francois Gygi committed
22 23 24 25
// The xc energy must be -2.8105 a.u.
// dExc/da must be 0.911682

#include<iostream>
Francois Gygi committed
26
#include<vector>
Francois Gygi committed
27 28 29
#include "LDAFunctional.h"
#include <cassert>
#include <cmath>
Francois Gygi committed
30
#include <cstdlib>
Francois Gygi committed
31 32 33 34 35 36
using namespace std;

int main(int argc, char **argv)
{
  if ( argc != 3 )
  {
Francois Gygi committed
37 38
    cout << " use: testLDAFunctional alat n" << endl;
    return 1;
Francois Gygi committed
39 40 41 42 43
  }
  double a = atof(argv[1]);
  double omega = a*a*a;
  int n = atoi(argv[2]);
  int n3 = n*n*n;
Francois Gygi committed
44 45
  vector<vector<double> > rh(1);
  rh[0].resize(n3);
Francois Gygi committed
46
  double excsum = 0.0, dxcsum = 0.0;
47

Francois Gygi committed
48 49 50 51
  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;
52

Francois Gygi committed
53 54 55 56 57 58 59 60 61 62 63
  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 );
Francois Gygi committed
64 65
        rh[0][ii] = fac * exp( -r2 / (rc*rc) );
        sum += rh[0][ii];
Francois Gygi committed
66 67 68
      }
    }
  }
69
  sum = sum * omega / n3;
Francois Gygi committed
70
  // the density should be normalized
71 72
  cout << " Integrated density: " << sum << endl;

Francois Gygi committed
73 74
  LDAFunctional xcf(rh);
  xcf.setxc();
75

Francois Gygi committed
76
  for ( int i = 0; i < n3; i++ )
Francois Gygi committed
77
    excsum += xcf.rho[i] * xcf.exc[i];
Francois Gygi committed
78
  for ( int i = 0; i < n3; i++ )
Francois Gygi committed
79
    dxcsum += xcf.rho[i] * ( xcf.exc[i] - xcf.vxc1[i] );
80

Francois Gygi committed
81
  cout << " Total LDA xc energy: " << excsum * omega / n3 << endl;
82

Francois Gygi committed
83 84 85 86 87
  // Note: the energy variation is 3 * dExc/da * delta(a)
  cout << " dExc/da: " << dxcsum * omega / ( n3 * a ) << endl;

  return 0;
}