src/QMCHamiltonians/ECPComponentBuilder.cpp

//////////////////////////////////////////////////////////////////////////////////////
// This file is distributed under the University of Illinois/NCSA Open Source License.
// See LICENSE file in top directory for details.
//
// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
//
// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
//                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////


#include "ECPComponentBuilder.h"
#include "Numerics/GaussianTimesRN.h"
#include "Numerics/Quadrature.h"
#include "QMCHamiltonians/NonLocalECPComponent.h"
#include "Numerics/Transform2GridFunctor.h"
#include "Utilities/IteratorUtility.h"
#include "Utilities/SimpleParser.h"
#include "Message/CommOperators.h"
#include <cmath>
#include "Utilities/qmc_common.h"


namespace qmcplusplus
{
ECPComponentBuilder::ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule)
    : MPIObjectBase(c),
      NumNonLocal(0),
      Lmax(0),
      Nrule(nrule),
      Srule(8),
      AtomicNumber(0),
      Zeff(0),
      RcutMax(-1),
      Species(aname),
      grid_global(0)
{
  angMon["s"] = 0;
  angMon["p"] = 1;
  angMon["d"] = 2;
  angMon["f"] = 3;
  angMon["g"] = 4;
  angMon["h"] = 5;
  angMon["i"] = 6;
  angMon["j"] = 7;
  angMon["k"] = 8;
  angMon["0"] = 0;
  angMon["1"] = 1;
  angMon["2"] = 2;
  angMon["3"] = 3;
  angMon["4"] = 4;
  angMon["5"] = 5;
  angMon["6"] = 6;
  angMon["7"] = 7;
  angMon["8"] = 8;
}

bool ECPComponentBuilder::parse(const std::string& fname, xmlNodePtr cur)
{
  const XMLAttrString cutoff_str(cur, "cutoff");
  if (!cutoff_str.empty())
    RcutMax = std::stod(cutoff_str);

  return read_pp_file(fname);
}

int ReadFileBuffer::get_file_length(std::ifstream* f) const
{
  f->seekg(0, std::ios::end);
  int len = f->tellg();
  f->seekg(0, std::ios::beg);
  return len;
}

void ReadFileBuffer::reset()
{
  if (is_open)
  {
    if (myComm == NULL || myComm->rank() == 0)
    {
      delete fin;
      fin = NULL;
    }
    delete[] cbuffer;
    cbuffer = NULL;
    is_open = false;
    length  = 0;
  }
}

bool ReadFileBuffer::open_file(const std::string& fname)
{
  reset();

  if (myComm == NULL || myComm->rank() == 0)
  {
    fin = new std::ifstream(fname.c_str());
    if (fin->is_open())
      is_open = true;
  }
  if (myComm)
    myComm->bcast(is_open);
  return is_open;
}

bool ReadFileBuffer::read_contents()
{
  if (!is_open)
    return false;

  if (myComm == NULL || myComm->rank() == 0)
  {
    length = get_file_length(fin);
  }
  if (myComm)
    myComm->bcast(length);

  cbuffer         = new char[length + 1];
  cbuffer[length] = '\0';

  if (myComm == NULL || myComm->rank() == 0)
    fin->read(cbuffer, length);

  if (myComm != NULL)
    myComm->bcast(cbuffer, length);

  return true;
}


bool ECPComponentBuilder::read_pp_file(const std::string& fname)
{
  ReadFileBuffer buf(myComm);
  bool okay = buf.open_file(fname);
  if (!okay)
    myComm->barrier_and_abort("ECPComponentBuilder::read_pp_file  Missing PP file " + fname + "\n");

  okay = buf.read_contents();
  if (!okay)
    myComm->barrier_and_abort("ECPComponentBuilder::read_pp_file Unable to read PP file " + fname + "\n");

  xmlDocPtr m_doc = xmlReadMemory(buf.contents(), buf.length, NULL, NULL, 0);

  if (m_doc == NULL)
  {
    xmlFreeDoc(m_doc);
    myComm->barrier_and_abort("ECPComponentBuilder::read_pp_file xml file " + fname + " is invalid");
  }
  // Check the document is of the right kind
  xmlNodePtr cur = xmlDocGetRootElement(m_doc);
  if (cur == NULL)
  {
    xmlFreeDoc(m_doc);
    myComm->barrier_and_abort("Empty document");
  }
  bool success = put(cur);
  xmlFreeDoc(m_doc);
  return success;
}

bool ECPComponentBuilder::put(xmlNodePtr cur)
{
  //  int nk=0;
  //vector<RealType> kpts;
  std::vector<xmlNodePtr> semiPtr;
  cur = cur->children;
  while (cur != NULL)
  {
    std::string cname((const char*)cur->name);
    if (cname == "header")
    {
      Zeff         = std::stoi(XMLAttrString{cur, "zval"});
      AtomicNumber = std::stoi(XMLAttrString{cur, "atomic-number"});
    }
    else if (cname == "grid")
    {
      //capture the global grid
      grid_global = createGrid(cur);
    }
    else if (cname == "L2")
    {
      buildL2(cur);
    }
    else if (cname == "semilocal")
    {
      semiPtr.push_back(cur); //save the pointer
    }
    else if (cname == "local")
    {
      buildLocal(cur);
    }
    cur = cur->next;
  }
  if (semiPtr.size())
  {
    if (!pp_nonloc)
      pp_nonloc = std::make_unique<NonLocalECPComponent>();
    if (pp_so == 0)
      pp_so = std::make_unique<SOECPComponent>();
    if (pp_loc)
    {
      for (int i = 0; i < semiPtr.size(); i++)
        addSemiLocal(semiPtr[i]);
    }
    else
      buildSemiLocalAndLocal(semiPtr);
  }
  if (pp_nonloc)
  {
    SetQuadratureRule(Nrule);
    app_log() << "    Non-local pseudopotential parameters" << std::endl;
    pp_nonloc->print(app_log());
    app_log() << "    Maximum cutoff radius " << pp_nonloc->Rmax << std::endl;
  }
  return true;
}

void ECPComponentBuilder::printECPTable()
{
  if (!qmc_common.io_node || qmc_common.mpi_groups > 1)
    return;

  char fname[12];
  sprintf(fname, "%s.pp.dat", Species.c_str());
  std::ofstream fout(fname);
  fout.setf(std::ios::scientific, std::ios::floatfield);
  fout.precision(12);
  int nl      = pp_nonloc ? pp_nonloc->nlpp_m.size() : 0;
  RealType d  = 1.7e-2;
  RealType rt = 0.13 * d;
  if (nl)
  {
    fout << "#  Lmax = " << Lmax + 1 << " nonlocal L channels" << nl << std::endl;
    fout << "#  Units = bohr hartree " << std::endl;
    fout << "#  r  -r*V/zeff   Vnl ... " << std::endl;
    while (rt < 5)
    {
      fout << rt << std::setw(25) << pp_loc->splint(rt);
      for (int l = 0; l < nl; l++)
        fout << std::setw(25) << pp_nonloc->nlpp_m[l]->splint(rt);
      fout << std::endl;
      rt += d;
    }
  }
  else
  {
    fout << "#  Units = bohr hartree " << std::endl;
    fout << "#  r  -r*V/zeff " << std::endl;
    while (rt < 5)
    {
      fout << rt << std::setw(25) << pp_loc->splint(rt) << std::endl;
      ;
      rt += d;
    }
  }
}

void ECPComponentBuilder::SetQuadratureRule(int rule)
{
  app_log() << "  Quadrature Nrule: " << rule << std::endl;
  Quadrature3D<RealType> myRule(rule);
  pp_nonloc->sgridxyz_m    = myRule.xyz_m;
  pp_nonloc->sgridweight_m = myRule.weight_m;
  // Allocate storage for wave function ratios
  pp_nonloc->resize_warrays(myRule.nk, NumNonLocal, Lmax);
  if (pp_so)
  { //added here bc must have nonlocal terms to have SO contributions
    pp_so->sgridxyz_m    = myRule.xyz_m;
    pp_so->sgridweight_m = myRule.weight_m;
    pp_so->resize_warrays(myRule.nk, NumSO, Srule);
  }
}

} // namespace qmcplusplus