testLDAFunctional.C 2.46 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 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
// The xc energy must be -2.8105 a.u.
// dExc/da must be 0.911682

#include<iostream>
#include "LDAFunctional.h"
#include <cassert>
#include <cmath>
using namespace std;

int main(int argc, char **argv)
{
  // use: testxcf alat np
  if ( argc != 3 )
  {
    cout << " use: testLDAFunctional 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;
  double *rh = new double[n3];
  double *vxc = new double[n3];
  double *exc = new double[n3];
  double excsum = 0.0, dxcsum = 0.0;
48

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

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

Francois Gygi committed
74
  LDAFunctional xcf;
75

Francois Gygi committed
76 77 78 79
  int nspin = 1;
  xcf.rho = rh;
  xcf.vxc1 = vxc;
  xcf.exc = exc;
80

Francois Gygi committed
81
  xcf.setxc(n3,nspin);
82

Francois Gygi committed
83 84 85 86
  for ( int i = 0; i < n3; i++ )
    excsum += rh[i] * exc[i];
  for ( int i = 0; i < n3; i++ )
    dxcsum += rh[i] * ( exc[i] - vxc[i] );
87

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

Francois Gygi committed
90 91 92 93 94
  // Note: the energy variation is 3 * dExc/da * delta(a)
  cout << " dExc/da: " << dxcsum * omega / ( n3 * a ) << endl;

  return 0;
}