testChargeDensity.C 4.08 KB
Newer Older
Francois Gygi committed
1 2 3 4 5 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 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155
////////////////////////////////////////////////////////////////////////////////
//
// testChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: testChargeDensity.C,v 1.1 2003-01-25 01:23:31 fgygi Exp $

#include "Context.h"
#include "Wavefunction.h"
#include "ChargeDensity.h"
#include "SlaterDet.h"
#include "FourierTransform.h"
#include "Timer.h"

#include <iostream>
#include <iomanip>
#include <cassert>
using namespace std;

#ifdef USE_MPI
#include <mpi.h>
#endif

int main(int argc, char **argv)
{
#if USE_MPI
  MPI_Init(&argc,&argv);
#endif
  {
    // use: 
    // testChargeDensity a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut nel nempty nspin nkp
    assert(argc==15);
    D3vector a(atof(argv[1]),atof(argv[2]),atof(argv[3]));
    D3vector b(atof(argv[4]),atof(argv[5]),atof(argv[6]));
    D3vector c(atof(argv[7]),atof(argv[8]),atof(argv[9]));
    UnitCell cell(a,b,c);
    cout << " volume: " << cell.volume() << endl;
    double ecut = atof(argv[10]);
    int nel = atoi(argv[11]);
    int nempty = atoi(argv[12]);
    int nspin = atoi(argv[13]);
    int nkp = atoi(argv[14]);
    
    Timer tm;
    
    Context ctxt;
    Wavefunction wf(ctxt);
    
    tm.reset(); tm.start();
    wf.resize(cell,cell,ecut);
    tm.stop();
    cout << " wf.resize: CPU/Real: " 
         << tm.cpu() << " / " << tm.real() << endl;

    tm.reset(); tm.start();
    wf.set_nel(nel);
    tm.stop();
    cout << " wf.set_nel: CPU/Real: " 
         << tm.cpu() << " / " << tm.real() << endl;
         
    tm.reset(); tm.start();
    wf.set_nspin(nspin);
    tm.stop();
    cout << " wf.set_nspin: CPU/Real: " 
         << tm.cpu() << " / " << tm.real() << endl;
    
    for ( int ikp = 0; ikp < nkp-1; ikp++ )
    {
      wf.add_kpoint(D3vector((0.5*(ikp+1))/(nkp-1),0,0),1.0);
    }
    
    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
    {
      for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
      {
        if ( wf.sd(ispin,ikp) != 0 && wf.sdcontext(ispin,ikp)->active() )
        {
        cout << "wf.sd(ispin=" << ispin << ",ikp=" << ikp << "): "
             << wf.sd(ispin,ikp)->c().m() << "x"
             << wf.sd(ispin,ikp)->c().n() << endl;
        cout << ctxt.mype() << ":"
             << " sdcontext[" << ispin << "][" << ikp << "]: " 
             << wf.sd(ispin,ikp)->context();
        }
      }
    }

    tm.reset();
    tm.start();
    wf.randomize(0.1);
    tm.stop();
    cout << " wf.randomize: CPU/Real: " 
         << tm.cpu() << " / " << tm.real() << endl;

    tm.reset();
    tm.start();
    wf.gram();
    cout << " wf.gram: CPU/Real: " 
         << tm.cpu() << " / " << tm.real() << endl;
         
    wf.update_occ();
        
    // compute charge density in real space
    Timer tmrho;
    tmrho.reset();
    tmrho.start();
    cout << " ChargeDensity::ctor..." << endl;
    ChargeDensity cd(wf);
    
    cout << "done" << endl;
    
    tmrho.stop();
    cout << " ChargeDensity::ctor: CPU/Real: " 
         << tmrho.cpu() << " / " << tmrho.real() << endl;
         
    tmrho.reset();
    tmrho.start();
    cout << " ChargeDensity::update_density..." << endl;
    cd.update_density();
    tmrho.stop();
    cout << " ChargeDensity::update_density: CPU/Real: " 
         << tmrho.cpu() << " / " << tmrho.real() << endl;
         
    // print the first few Fourier coefficients of the charge density
    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
    {
      for ( int i = 0; i < cd.vbasis()->localsize(); i++ )
      {
        cout << ctxt.mype() << ": rho(ispin=" << ispin << ",ig=" << i << " (" 
        << cd.vbasis()->idx(3*i) << " "
        << cd.vbasis()->idx(3*i+1) << " " << cd.vbasis()->idx(3*i+2) << ") "
        << cd.rhog[ispin][i] << endl;
      }
    }
         
#if 0
    // integral of rho in r space
    double sum = 0.0;
    for ( int i = 0; i < rho.size(); i++ )
      sum += rho[i];
#if USE_MPI
    double tsum;
    int mycol = sd.context().mycol();
    MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,sd.context().comm());
    sum = tsum;
#endif
    cout << ctxt.mype() << ": " << " rho: " << sum / ft.np012() << endl;

#endif

  }
#if USE_MPI
  MPI_Finalize();
#endif
}