ComputeMLWFCmd.C 2.75 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
// ComputeMLWFCmd.C:
//
////////////////////////////////////////////////////////////////////////////////

#include "ComputeMLWFCmd.h"
#include<iostream>
#include "Context.h"
#include "SlaterDet.h"
using namespace std;

Francois Gygi committed
25
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
26 27 28
int ComputeMLWFCmd::action(int argc, char **argv)
{
  Wavefunction& wf = s->wf;
29

30
  // Check that only the k=0 point is used
31
  if ( wf.nkp()>1 || length(wf.kpoint(0)) != 0.0 )
32 33 34 35 36 37 38 39 40
  {
    if ( ui->onpe0() )
    {
      cout << " ComputeMLWFCmd::action: compute_mlwf can only be used at\n"
           << " the Gamma point (k=0)" << endl;
    }
    return 1;
  }

Francois Gygi committed
41
  if ( ui->onpe0() )
42 43 44 45
    cout << "<mlwfs>" << endl;

  D3vector edipole_sum;
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
Francois Gygi committed
46
  {
47 48 49 50 51 52 53 54
    SlaterDet& sd = *(wf.sd(ispin,0));
    MLWFTransform* mlwft = new MLWFTransform(sd);

    mlwft->update();
    mlwft->compute_transform();
    mlwft->apply_transform(sd);

    if ( ui->onpe0() )
Francois Gygi committed
55
    {
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
      cout << " <mlwfset spin=\"" << ispin
           << "\" size=\"" << sd.nst() << "\">" << endl;
      for ( int i = 0; i < sd.nst(); i++ )
      {
        D3vector ctr = mlwft->center(i);
        double sp = mlwft->spread(i);
        cout.setf(ios::fixed, ios::floatfield);
        cout.setf(ios::right, ios::adjustfield);
        cout << "   <mlwf center=\"" << setprecision(6)
             << setw(12) << ctr.x
             << setw(12) << ctr.y
             << setw(12) << ctr.z
             << " \" spread=\" " << sp << " \"/>"
             << endl;
      }
      D3vector edipole = mlwft->dipole();
      cout << " <electronic_dipole spin=\"" << ispin << "\"> " << edipole
           << " </electronic_dipole>" << endl;
      cout << " </mlwfset>" << endl;
      edipole_sum += edipole;
Francois Gygi committed
76
    }
77 78 79 80 81
    delete mlwft;
  }

  if ( ui->onpe0() )
  {
Francois Gygi committed
82
    D3vector idipole = s->atoms.dipole();
83
    cout << "   <ionic_dipole> " << idipole
Francois Gygi committed
84
         << " </ionic_dipole>" << endl;
85
    cout << "   <total_dipole> " << idipole + edipole_sum
Francois Gygi committed
86
         << " </total_dipole>" << endl;
87
    cout << "   <total_dipole_length> " << length(idipole + edipole_sum)
Francois Gygi committed
88
         << " </total_dipole_length>" << endl;
89
    cout << "</mlwfs>" << endl;
Francois Gygi committed
90
  }
Francois Gygi committed
91
  return 0;
Francois Gygi committed
92
}