testChargeDensity.C 4.47 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6 7 8 9 10 11 12 13 14
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License 
// as published by the Free Software Foundation, either version 2 of 
// 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
// testChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
18
// $Id: testChargeDensity.C,v 1.3 2008-08-13 06:39:43 fgygi Exp $
Francois Gygi committed
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

#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
  {
42
    // use:
Francois Gygi committed
43 44 45 46 47 48 49 50 51 52 53 54
    // 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]);
55

Francois Gygi committed
56
    Timer tm;
57

Francois Gygi committed
58 59
    Context ctxt;
    Wavefunction wf(ctxt);
60

Francois Gygi committed
61 62 63
    tm.reset(); tm.start();
    wf.resize(cell,cell,ecut);
    tm.stop();
64
    cout << " wf.resize: CPU/Real: "
Francois Gygi committed
65 66 67 68 69
         << tm.cpu() << " / " << tm.real() << endl;

    tm.reset(); tm.start();
    wf.set_nel(nel);
    tm.stop();
70
    cout << " wf.set_nel: CPU/Real: "
Francois Gygi committed
71
         << tm.cpu() << " / " << tm.real() << endl;
72

Francois Gygi committed
73 74 75
    tm.reset(); tm.start();
    wf.set_nspin(nspin);
    tm.stop();
76
    cout << " wf.set_nspin: CPU/Real: "
Francois Gygi committed
77
         << tm.cpu() << " / " << tm.real() << endl;
78

Francois Gygi committed
79 80 81 82
    for ( int ikp = 0; ikp < nkp-1; ikp++ )
    {
      wf.add_kpoint(D3vector((0.5*(ikp+1))/(nkp-1),0,0),1.0);
    }
83

Francois Gygi committed
84 85 86 87 88 89 90 91 92 93
    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() << ":"
94
             << " sdcontext[" << ispin << "][" << ikp << "]: "
Francois Gygi committed
95 96 97 98 99 100 101 102 103
             << wf.sd(ispin,ikp)->context();
        }
      }
    }

    tm.reset();
    tm.start();
    wf.randomize(0.1);
    tm.stop();
104
    cout << " wf.randomize: CPU/Real: "
Francois Gygi committed
105 106 107 108 109
         << tm.cpu() << " / " << tm.real() << endl;

    tm.reset();
    tm.start();
    wf.gram();
110
    cout << " wf.gram: CPU/Real: "
Francois Gygi committed
111
         << tm.cpu() << " / " << tm.real() << endl;
112

Francois Gygi committed
113
    wf.update_occ();
114

Francois Gygi committed
115 116 117 118 119 120
    // compute charge density in real space
    Timer tmrho;
    tmrho.reset();
    tmrho.start();
    cout << " ChargeDensity::ctor..." << endl;
    ChargeDensity cd(wf);
121

Francois Gygi committed
122
    cout << "done" << endl;
123

Francois Gygi committed
124
    tmrho.stop();
125
    cout << " ChargeDensity::ctor: CPU/Real: "
Francois Gygi committed
126
         << tmrho.cpu() << " / " << tmrho.real() << endl;
127

Francois Gygi committed
128 129 130 131 132
    tmrho.reset();
    tmrho.start();
    cout << " ChargeDensity::update_density..." << endl;
    cd.update_density();
    tmrho.stop();
133
    cout << " ChargeDensity::update_density: CPU/Real: "
Francois Gygi committed
134
         << tmrho.cpu() << " / " << tmrho.real() << endl;
135

Francois Gygi committed
136 137 138 139 140
    // 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++ )
      {
141
        cout << ctxt.mype() << ": rho(ispin=" << ispin << ",ig=" << i << " ("
Francois Gygi committed
142 143 144 145 146
        << cd.vbasis()->idx(3*i) << " "
        << cd.vbasis()->idx(3*i+1) << " " << cd.vbasis()->idx(3*i+2) << ") "
        << cd.rhog[ispin][i] << endl;
      }
    }
147

Francois Gygi committed
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
#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
}