testVWN.C 1.89 KB
 Francois Gygi committed Aug 22, 2015 1 2 3 4 //////////////////////////////////////////////////////////////////////////////// // // testVWN.C //  Francois Gygi committed Mar 05, 2016 5 // test the Vosko-Wilk-Nusair functional  Francois Gygi committed Aug 22, 2015 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 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 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 // //////////////////////////////////////////////////////////////////////////////// #include #include "VWNFunctional.h" #include "LDAFunctional.h" #include using namespace std; // c1 = (3.D0/(4.D0*pi))**third const double c1 = 0.6203504908994001; const double c3 = -0.610887057711; double rho_of_rs(double rs) { return pow(c1/rs,3.0); } int main() { const int n = 500; vector > rho; rho.resize(1); rho[0].resize(n); LDAFunctional f(rho); double rs[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.5, 10.0, 15.0, 20.0, 50.0, 100.0 }; // print table for given rs values (Table 5 in paper) for ( int i = 0; i < 13; i++ ) rho[0][i] = rho_of_rs(rs[i]); f.setxc(); for ( int i = 0; i < 13; i++ ) { double vx = c3 / rs[i]; double ex = 0.75 * vx; cout << "rs=" << rs[i] << " " << f.rho[i] << " " << f.exc[i] << " " << f.vxc1[i] << endl; } cout << endl; #if 0 const double drh = 0.01; for ( int i = 0; i < rho[0].size(); i++ ) { double rh = drh * i; rho[0][i] = rh; } f.setxc(); for ( int i = 1; i < rho[0].size()-1; i++ ) { cout << f.rho[i] << " " << f.exc[i] << " " << f.vxc1[i] << " " << (f.rho[i+1]*f.exc[i+1]-f.rho[i-1]*f.exc[i-1])/(2*drh) << endl; } #endif // spin-polarized test rho.resize(2); rho[0].resize(500); rho[1].resize(500); for ( int i = 0; i < 13; i++ ) { rho[0][i] = rho_of_rs(rs[i]); rho[1][i] = 0.0; } LDAFunctional fs(rho); fs.setxc(); for ( int i = 0; i < 13; i++ ) { const double c4 = -0.769669463118; double vx = c4 / rs[i]; double ex = 0.75 * vx; cout << "rs=" << rs[i] << " " << fs.rho_up[i] << " " << fs.exc[i] << " " << fs.vxc1_up[i] << endl; } cout << endl; }