BisectionCmd.h 3.75 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
////////////////////////////////////////////////////////////////////////////////
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// BisectionCmd.h:
//
////////////////////////////////////////////////////////////////////////////////

#ifndef BISECTIONCMD_H
#define BISECTIONCMD_H

#include <iostream>
#include "UserInterface.h"
#include "Sample.h"
#include <cstdlib>
#include <bitset>
#include "Bisection.h"
28
#include "Timer.h"
Francois Gygi committed
29 30 31 32 33 34 35 36 37

class BisectionCmd : public Cmd
{
  public:

  Sample *s;

  BisectionCmd(Sample *sample) : s(sample) {};

38 39
  const char *name(void) const { return "bisection"; }
  const char *help_msg(void) const
Francois Gygi committed
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
  {
    return
    "\n bisection\n\n"
    " syntax: bisection lx ly lz threshold\n\n"
    "   The bisection command performs recursive bisection on the\n"
    "   wavefunctions.\n\n";
  }

  int action(int argc, char **argv)
  {
    if ( argc != 5 )
    {
      if ( ui->onpe0() )
      {
        cout << " use: bisection lx ly lz threshold" << endl;
      }
      return 1;
    }

59 60
    Timer tm;

Francois Gygi committed
61 62 63 64 65 66 67
    Wavefunction &wf=s->wf;
    double epsilon=atof(argv[4]);
    int nLevels[3];
    nLevels[0]=atoi(argv[1]);
    nLevels[1]=atoi(argv[2]);
    nLevels[2]=atoi(argv[3]);

68 69 70 71 72 73 74 75 76
    if ( epsilon < 0.0 )
    {
      if ( ui->onpe0() )
      {
        cout << " BisectionCmd: threshold must be non-negative" << endl;
      }
      return 1;
    }

77 78 79 80 81 82 83 84 85
    if ( wf.nkp() > 1 )
    {
      if ( ui->onpe0() )
      {
        cout << " BisectionCmd: only implemented for k=0" << endl;
      }
      return 1;
    }

86 87 88 89 90 91 92 93 94 95 96
    if ( nLevels[0] < 0 || nLevels[0] > 5 ||
         nLevels[1] < 0 || nLevels[1] > 5 ||
         nLevels[2] < 0 || nLevels[2] > 5 )
    {
      if ( ui->onpe0() )
      {
        cout << " BisectionCmd: levels must be in [0,5]" << endl;
      }
      return 1;
    }

97
    tm.reset();
98
    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
Francois Gygi committed
99
    {
100 101 102 103
      SlaterDet& sd = *wf.sd(0,0);
      Bisection bisection(sd,nLevels);
      const int maxsweep = 20;
      const double tol = 1.e-8;
104
      tm.start();
105 106 107
      bisection.compute_transform(sd,maxsweep,tol);
      bisection.compute_localization(epsilon);
      bisection.forward(sd);
108
      tm.stop();
109 110 111
      vector<long int> localization = bisection.localization();

      if ( ui->onpe0() )
Francois Gygi committed
112
      {
113 114
        cout << " BisectionCmd: lx=" << nLevels[0]
             << " ly=" << nLevels[1]
Francois Gygi committed
115 116
             << " lz=" << nLevels[2]
             << " threshold=" << epsilon << endl;
117 118 119 120 121 122 123 124 125 126 127 128 129 130

        // print localization vectors and overlaps
        int sum = 0;
        const int nst = localization.size();
        for ( int i = 0; i < nst; i++ )
        {
          int count = 0;
          for ( int j = 0; j < nst; j++ )
          {
            if ( bisection.overlap(i,j) )
              count++;
          }
          cout << "localization[" << i << "]: "
               << localization[i] << " "
131 132 133
               << bitset<30>(localization[i])
               << " size: " << bisection.size(i)
               << " overlaps: " << count << endl;
134 135
          sum += count;
        }
136 137 138 139
        cout << " Bisection: total size:    ispin=" << ispin
             << ": " << bisection.total_size() << endl;
        cout << " Bisection: pair fraction: ispin=" << ispin
             << ": " << bisection.pair_fraction() << endl;
Francois Gygi committed
140 141
      }
    }
142 143
    if ( ui->onpe0() )
      cout << " Bisection: time=" << tm.real() << endl;
Francois Gygi committed
144 145 146 147
    return 0;
  }
};
#endif