ExternalPotential.C 7.34 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2016 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ExternalPotential.C
//
////////////////////////////////////////////////////////////////////////////////

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

24 25 26
#include "Basis.h"
#include "ExternalPotential.h"
#include "FourierTransform.h"
27
#include "Function3d.h"
28 29
#include "Base64Transcoder.h"

30
////////////////////////////////////////////////////////////////////////////////
31 32
bool abs_compare(const double &a, const double &b)
{
33 34
  return (abs(a) < abs(b));
}
35

36
////////////////////////////////////////////////////////////////////////////////
37 38
void ExternalPotential::update(const ChargeDensity& cd)
{
39 40 41 42 43
  const Context* ctxt = s_.wf.spincontext();
  bool onpe0 = ctxt->onpe0();
  int nprow = ctxt->nprow();
  int myrow = ctxt->myrow();
  int mycol = ctxt->mycol();
44
  MPI_Comm vcomm = cd.vcomm();
45

46
  Timer tm_read_vext;
47 48
  double time, tmin, tmax;

49 50
  // In cube mode and xml mode, the external potential is
  // read by processors on row 0 and stored in vext_read
51 52 53
  vector<double> vext_read, vext_read_loc;

  tm_read_vext.start();
54
  if ( fmt_ == "cube" )
55
  {
56 57
    // read cube file, n_'s are determined by cube file
    if ( myrow == 0 )
58
    {
59
      ifstream vfile(filename_.c_str());
60 61 62 63 64 65 66 67 68 69
      if (!vfile)
      {
        if (mycol == 0)
          cout << "  ExternalPotential::update: file not found: "
               << filename_ << endl;
        ctxt->abort(1);
      }
      string tmpstr;
      for (int i = 0; i < 2; i++)
        getline(vfile, tmpstr);  // skip comments
70 71
      int natom;
      vfile >> natom;
72 73
      getline(vfile, tmpstr);
      for (int i = 0; i < 3; i++)
74 75
      {
        vfile >> n_[i];
76
        getline(vfile, tmpstr);
77
      }
78 79
      for (int i = 0; i < natom; i++)
        getline(vfile, tmpstr);  // skip atom coordinates
80
      vext_read.resize(n_[0] * n_[1] * n_[2]);
81 82 83
      for (int nx = 0; nx < n_[0]; nx++)
        for (int ny = 0; ny < n_[1]; ny++)
          for (int nz = 0; nz < n_[2]; nz++)
84 85
          {
            const int ir = nx + ny * n_[0] + nz * n_[0] * n_[1];
86
            vfile >> vext_read[ir];
87
          }
88
      vfile.close();
89
    }
90
    MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
91
  }
92
  else if ( fmt_ == "xml" )
93
  {
94 95
    if (myrow == 0)
    {
96 97 98 99 100 101
      Function3d f;
      f.read(filename_);
      vext_read = f.val;
      n_[0] = f.nx;
      n_[1] = f.ny;
      n_[2] = f.nz;
102
    }
103
    MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
104
  }
105
  tm_read_vext.stop();
mahe committed
106 107 108 109 110
  // at this point, all processes have correct n_ regardless of io mode

  // now construct 2 Fourier transforms ft1 and ft2
  // ft1 is used to transform vext_read_loc to vext_g (G space)
  // ft2 is used to transform vext_g to vext_r_ (R space)
111
  // the whole process is a Fourier interpolation/extrapolation
112
  // of the external potential on the charge density grid
113

mahe committed
114 115
  // create a Basis with largest possible ecut that is compatible with
  // the external potential grid from file
116
  const UnitCell& cell = cd.vbasis()->cell();
117 118 119
  int n0max = (n_[0]-2)/2;
  int n1max = (n_[1]-2)/2;
  int n2max = (n_[2]-2)/2;
120 121 122
  double ecut0 = 0.5 * norm2(cell.b(0)) * n0max*n0max;
  double ecut1 = 0.5 * norm2(cell.b(1)) * n1max*n1max;
  double ecut2 = 0.5 * norm2(cell.b(2)) * n2max*n2max;
123
  ecut_ = min(min(ecut0,ecut1),ecut2);
124
  Basis basis(vcomm,D3vector(0,0,0));
125
  basis.resize(cell,cell,ecut_);
126

127
  FourierTransform *vft = cd.vft();
128

129 130 131
  FourierTransform ft1(basis,n_[0],n_[1],n_[2]);
  vext_read_loc.resize(ft1.np012loc());
  vector<complex<double> > vext_g(basis.localsize());
132

133 134 135
  // check that the basis fits in the vft grid
  //assert(basis.fits_in_grid(vft->np0(),vft->np1(),vft->np2()));

136 137 138
  FourierTransform ft2(basis,vft->np0(),vft->np1(),vft->np2());
  vext_r_.resize(ft2.np012loc());

139 140 141 142 143 144
  // xml mode or cube mode: processors on row 0 scatter
  // vext to other rows
  vector<int> scounts(nprow,0);
  vector<int> sdispls(nprow,0);
  int displ = 0;
  for ( int iproc = 0; iproc < nprow; iproc++ )
145
  {
146 147 148
    sdispls[iproc] = displ;
    scounts[iproc] = ft1.np012loc(iproc);
    displ += ft1.np012loc(iproc);
149
  }
150 151
  MPI_Scatterv(&vext_read[0],&scounts[0],&sdispls[0],MPI_DOUBLE,
               &vext_read_loc[0],ft1.np012loc(),MPI_DOUBLE,0,vcomm);
mahe committed
152

153 154 155 156 157 158
  // now vext_read_loc on all processors contains the correct portion of vext
  // Fourier forward transform vext_read_loc to vext_g
  vector<complex<double> > tmp_r(ft1.np012loc());
  for ( int ir = 0; ir < tmp_r.size(); ir++ )
    tmp_r[ir] = complex<double>(vext_read_loc[ir],0);
  ft1.forward(&tmp_r[0],&vext_g[0]);
mahe committed
159
  // Fourier backward transform vext_g to vext_r_
160 161
  tmp_r.resize(ft2.np012loc());
  ft2.backward(&vext_g[0],&tmp_r[0]);
162 163
  for ( int i = 0; i < vext_r_.size(); i++ )
    vext_r_[i] = real(tmp_r[i]);
164 165 166 167
  if ( onpe0 )
  {
    cout << "  ExternalPotential::update: read external potential from file:  "
         << filename_ << endl;
168 169
    cout << "  ExternalPotential::update: grid size n0 = "
         << n_[0] << ", n1 = " << n_[1] << ", n2 = " << n_[2] << endl;
170 171 172 173
    cout << "  ExternalPotential::update: ecut:  " << ecut_ << endl;
    if ( amplitude_ != 0 )
      cout << "  ExternalPotential::update: amplitude:  " << amplitude_ << endl;
  }
174 175 176 177 178 179
  if ( amplitude_ == 0.0 )
  {
    // If amplitude_ = 0.0, use following scheme to get an amplitude.
    // Empirically, an absolute magnitude of 1.0E-3 ~ 1.0E-5 hartree for Vext
    // would be suitable. Here the amplitude is set to scale the
    // magnitude of vext to be 1.0E-3 hartree
180
    if (vext_read_loc.size() > 0)
181 182
      magnitude_ = abs(*max_element(vext_read_loc.begin(),
                       vext_read_loc.end(), abs_compare));
183 184
    ctxt->dmax('C',1,1,&magnitude_,1);
    MPI_Bcast(&magnitude_,1,MPI_DOUBLE,0,vcomm);
185
    amplitude_ = 1.0E-3 / magnitude_;
186
    if ( onpe0 )
187 188 189
      cout << "  ExternalPotential::update: amplitude automatically"
           << " determined to be " << amplitude_
           << " (max abs value of vext = " << magnitude_ << ")" << endl;
190 191 192 193 194 195 196 197 198
  }

  time = tm_read_vext.real();
  tmin = time;
  tmax = time;
  ctxt->dmin(1,1,&tmin,1);
  ctxt->dmax(1,1,&tmax,1);
  if ( onpe0 )
  {
199
    cout << "  ExternalPotential::update: vext read time "
200
         << "min: " << tmin << " max: " << tmax << endl;
201
  }
202
}
203

204
////////////////////////////////////////////////////////////////////////////////
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
double ExternalPotential::compute_eext(const ChargeDensity& cd)
{
  // Eext =  integral ( rhor * vext_r )
  double eext = 0.0;
  double omega = s_.wf.cell().volume();
  const int np012loc = cd.vft()->np012loc();
  const int np012 = cd.vft()->np012();
  const std::vector<std::vector<double> > &rhor = cd.rhor;
  for ( int ispin = 0; ispin < s_.wf.nspin(); ispin++ )
  {
    for ( int ir = 0; ir < np012loc; ir++ )
    {
      assert( rhor[ispin].size() == np012loc );
      assert( vext_r_.size() == np012loc );
      eext += rhor[ispin][ir] * v(ir);
    }
  }
  double eext_sum = 0.0;
  MPI_Allreduce(&eext,&eext_sum,1,MPI_DOUBLE,MPI_SUM,cd.vbasis()->comm());
  eext_sum =  eext_sum * omega / np012;
  return eext_sum;
}