testSlaterDet.C 5.13 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 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
// testSlaterDet.C
//
////////////////////////////////////////////////////////////////////////////////

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

#include <iostream>
#include <fstream>
#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: testSlaterDet a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut nst kx ky kz npr npc
Francois Gygi committed
41 42 43 44 45 46 47
    if ( argc != 17 )
    {
      cout <<
      "use: testSlaterDet a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut nst kx ky kz npr npc"
      << endl;
      return 1;
    }
Francois Gygi committed
48
    double err;
Francois Gygi committed
49

Francois Gygi committed
50 51 52 53 54 55 56
    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);
    double ecut = atof(argv[10]);
    int nst = atoi(argv[11]);
    D3vector kpoint(atof(argv[12]),atof(argv[13]),atof(argv[14]));
57

Francois Gygi committed
58 59
    int npr = atoi(argv[15]);
    int npc = atoi(argv[16]);
60

Francois Gygi committed
61
    Timer tm;
62

63
    Context ctxt(MPI_COMM_WORLD,npr,npc);
64

Francois Gygi committed
65 66
    SlaterDet sd(ctxt,kpoint);
    cout << sd.context();
67

Francois Gygi committed
68 69 70 71 72 73 74
    sd.resize(cell,cell,ecut,nst);
    if ( ctxt.myproc() == 0 )
    {
      cout << " kpoint: " << sd.basis().kpoint() << endl;
      cout << " ngw:    " << sd.basis().size() << endl;
      cout << " nst:    " << sd.nst() << endl;
    }
75

Francois Gygi committed
76 77 78
    err = sd.ortho_error();
    if ( ctxt.myproc() == 0 )
      cout << " Initial orthogonality error before rand " << err << endl;
79

Francois Gygi committed
80 81 82 83 84 85
    sd.randomize(0.2/sd.basis().size());
    if ( ctxt.myproc() == 0 )
      cout << " sd.randomize: done" << endl;
    err = sd.ortho_error();
    if ( ctxt.myproc() == 0 )
      cout << " Orthogonality error after rand " << err << endl;
86

Francois Gygi committed
87 88 89 90 91 92 93 94 95 96 97
    tm.reset();
    tm.start();
    sd.gram();
    tm.stop();
    cout << " Gram: CPU/Real: " << tm.cpu() << " / " << tm.real() << endl;
    err = sd.ortho_error();
    if ( ctxt.myproc() == 0 )
      cout << " Gram orthogonality error " << err << endl;

    SlaterDet sdm(ctxt,kpoint);
    sdm.resize(cell,cell,ecut,nst);
98

Francois Gygi committed
99 100
    sdm.c() = sd.c();
    sd.randomize(0.1/sd.basis().size());
101

Francois Gygi committed
102 103 104 105 106 107 108 109
    tm.reset();
    tm.start();
    sd.riccati(sdm);
    tm.stop();
    cout << " riccati: CPU/Real: " << tm.cpu() << " / " << tm.real() << endl;
    err = sd.ortho_error();
    if ( ctxt.myproc() == 0 )
      cout << " Riccati orthogonality error " << err << endl;
110

Francois Gygi committed
111 112 113 114 115 116 117 118 119 120
    sd.riccati(sdm);
    err = sd.ortho_error();
    if ( ctxt.myproc() == 0 )
    {
      cout << " Riccati orthogonality error " << err << endl;
      cout << " sd.total size (MB): " << setprecision(4)
           << sd.memsize() / 1048576.0 << endl;
      cout << " sd.local size (MB): " << setprecision(4)
           << sd.localmemsize() / 1048576.0 << endl;
    }
121

122
    //cout << " ekin: " << setprecision(8) << sd.ekin(occ) << endl;
123

Francois Gygi committed
124 125 126
    // compute charge density in real space
    FourierTransform ft(sd.basis(),
      2*sd.basis().np(0), 2*sd.basis().np(1), 2*sd.basis().np(2));
127

Francois Gygi committed
128 129
    vector<complex<double> > f(ft.np012loc());
    vector<double> rho(ft.np012loc());
130

Francois Gygi committed
131 132 133
    Timer tmrho;
    tmrho.reset();
    tmrho.start();
134

Francois Gygi committed
135 136
    cout << " compute_density..." << endl;
    double weight = 1.0;
137
    sd.compute_density(ft,weight,&rho[0]);
138

Francois Gygi committed
139
    tmrho.stop();
140
    cout << " compute_density: CPU/Real: "
Francois Gygi committed
141
         << tmrho.cpu() << " / " << tmrho.real() << endl;
142

Francois Gygi committed
143 144 145 146 147 148 149 150 151
    // 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;
    MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,sd.context().comm());
    sum = tsum;
#endif
152
    cout << ctxt.mype() << ": "
Francois Gygi committed
153
         << " rho: " << sum * cell.volume() / ft.np012() << endl;
154

Francois Gygi committed
155 156 157 158 159 160 161 162 163
    SlaterDet sdp(sd);
    //SlaterDet sdp(ctxt,kpoint);
    //sdp.resize(cell,cell,ecut,nst);
    //sdp = sd;
    sdp.randomize(0.001);
    sdp.gram();
    err = sdp.ortho_error();
    if ( ctxt.myproc() == 0 )
      cout << " Gram orthogonality error " << err << endl;
164

Francois Gygi committed
165 166 167 168 169
    Timer tmv;
    tmv.reset();
    tmv.start();
    sd.rs_mul_add(ft,&rho[0],sdp);
    tmv.stop();
170
    cout << " rs_mul_add: CPU/Real: "
Francois Gygi committed
171
         << tmv.cpu() << " / " << tmv.real() << endl;
172

Francois Gygi committed
173 174 175 176 177 178
#if 0
    ofstream outfile("sd.dat");
    tm.reset();
    tm.start();
    outfile << sd;
    tm.stop();
179
    cout << " write: CPU/Real: "
Francois Gygi committed
180 181
         << tm.cpu() << " / " << tm.real() << endl;
#endif
182

Francois Gygi committed
183 184
    for ( TimerMap::iterator i = sd.tmap.begin(); i != sd.tmap.end(); i++ )
      cout << ctxt.mype() << ": "
185
           << setw(15) << (*i).first
Francois Gygi committed
186 187 188 189 190 191
           << " : " << setprecision(3) << (*i).second.real() << endl;
  }
#if USE_MPI
  MPI_Finalize();
#endif
}