src/QMCHamiltonians/ECPotentialBuilder.cpp
/** reimplement simple table format used by NonLocalPPotential
 */
reimplement simple table format used by NonLocalPPotential
reimplement simple table format used by NonLocalPPotential
void ECPotentialBuilder::useSimpleTableFormat()
{
  SpeciesSet& Species(IonConfig.getSpeciesSet());
  int ng(Species.getTotalNum());
  int icharge(Species.addAttribute("charge"));
  for (int ig = 0; ig < ng; ig++)
  {
    std::vector<RealType> grid_temp, pp_temp;
    std::string species(Species.speciesName[ig]);
    std::string fname = species + ".psf";
    std::ifstream fin(fname.c_str(), std::ios_base::in);
    if (!fin)
    {
      ERRORMSG("Could not open file " << fname)
      exit(-1);
    }
    int npotentials;
    fin >> npotentials;
    int lmax      = -1;
    int numnonloc = 0;
    RealType rmax(0.0);
    app_log() << "  ECPotential for " << species << std::endl;
    NonLocalECPComponent* mynnloc = 0;
    typedef OneDimCubicSpline<RealType> CubicSplineFuncType;
    for (int ij = 0; ij < npotentials; ij++)
    {
      int angmom, npoints;
      fin >> angmom >> npoints;
      OneDimNumGridFunctor<RealType> inFunc;
      inFunc.put(npoints, fin);
      if (angmom < 0)
      {
        RealType zinv   = -1.0 / Species(icharge, ig);
        int ng          = npoints - 1;
        RealType rf     = 5.0;
        ng              = static_cast<int>(rf * 100) + 1; 
        GridType* agrid = new LinearGrid<RealType>;
        agrid->set(0, rf, ng);
        std::vector<RealType> pp_temp(ng);
        pp_temp[0] = 0.0;
        for (int j = 1; j < ng; j++)
        {
          RealType r((*agrid)[j]);
          pp_temp[j] = r * zinv * inFunc.splint(r);
        }
        pp_temp[ng - 1]          = 1.0;
        RadialPotentialType* app = new RadialPotentialType(agrid, pp_temp);
        app->spline();
        localPot[ig] = app;
        app_log() << "    LocalECP l=" << angmom << std::endl;
        app_log() << "      Linear grid=[0," << rf << "] npts=" << ng << std::endl;
        hasLocalPot = true; 
      }
      else
      {
        hasNonLocalPot = true; 
        if (mynnloc == 0)
          mynnloc = new NonLocalECPComponent;
        RealType rf     = inFunc.rmax();
        GridType* agrid = new LinearGrid<RealType>;
        int ng          = static_cast<int>(rf * 100) + 1;
        agrid->set(0.0, rf, ng);
        app_log() << "    NonLocalECP l=" << angmom << " rmax = " << rf << std::endl;
        app_log() << "      Linear grid=[0," << rf << "] npts=" << ng << std::endl;
        std::vector<RealType> pp_temp(ng);
        pp_temp[0] = inFunc(0);
        for (int j = 1; j < ng; j++)
        {
          pp_temp[j] = inFunc.splint((*agrid)[j]);
        }
        RadialPotentialType* app = new RadialPotentialType(agrid, pp_temp);
        app->spline();
        mynnloc->add(angmom, app);
        lmax = std::max(lmax, angmom);
        rmax = std::max(rmax, rf);
        numnonloc++;
      }
      if (mynnloc)
      {
        mynnloc->setLmax(lmax);
        mynnloc->setRmax(rmax);
        app_log() << "    Maximum cutoff of NonLocalECP " << rmax << std::endl;
      }
    }
    fin.close();
    if (mynnloc)
    {
      nonLocalPot[ig]   = mynnloc;
      int numsgridpts   = 0;
      std::string fname = species + ".sgr";
      std::ifstream fin(fname.c_str(), std::ios_base::in);
      if (!fin)
      {
        app_error() << "Could not open file " << fname << std::endl;
        exit(-1);
      }
      PosType xyz;
      RealType weight;
      while (fin >> xyz >> weight)
      {
        mynnloc->addknot(xyz, weight);
        numsgridpts++;
      }
      mynnloc->resize_warrays(numsgridpts, numnonloc, lmax);
    }
  } 
}
void ECPotentialBuilder :: useSimpleTableFormat ( ) { SpeciesSet & Species ( IonConfig . getSpeciesSet ( ) ) ; int ng ( Species . getTotalNum ( ) ) ; int icharge ( Species . addAttribute ( " " ) ) ; for ( int ig = 0 ; ig < ng ; ig ++ ) { std :: vector < RealType > grid_temp , pp_temp ; std :: string species ( Species . speciesName [ ig ] ) ; std :: string fname = species + " " ; std :: ifstream fin ( fname . c_str ( ) , std :: ios_base :: in ) ; if ( ! fin ) { ERRORMSG ( " " << fname ) exit ( -1 ) ; } int npotentials ; fin >> npotentials ; int lmax = -1 ; int numnonloc = 0 ; RealType rmax ( 0.0 ) ; app_log ( ) << " " << species << std :: endl ; NonLocalECPComponent * mynnloc = 0 ; typedef OneDimCubicSpline < RealType > CubicSplineFuncType ; for ( int ij = 0 ; ij < npotentials ; ij ++ ) { int angmom , npoints ; fin >> angmom >> npoints ; OneDimNumGridFunctor < RealType > inFunc ; inFunc . put ( npoints , fin ) ; if ( angmom < 0 ) { RealType zinv = -1.0 / Species ( icharge , ig ) ; int ng = npoints - 1 ; RealType rf = 5.0 ; ng = static_cast < int > ( rf * 100 ) + 1 ; GridType * agrid = new LinearGrid < RealType > ; agrid -> set ( 0 , rf , ng ) ; std :: vector < RealType > pp_temp ( ng ) ; pp_temp [ 0 ] = 0.0 ; for ( int j = 1 ; j < ng ; j ++ ) { RealType r ( ( * agrid ) [ j ] ) ; pp_temp [ j ] = r * zinv * inFunc . splint ( r ) ; } pp_temp [ ng - 1 ] = 1.0 ; RadialPotentialType * app = new RadialPotentialType ( agrid , pp_temp ) ; app -> spline ( ) ; localPot [ ig ] = app ; app_log ( ) << " " << angmom << std :: endl ; app_log ( ) << " " << rf << " " << ng << std :: endl ; hasLocalPot = true ; } else { hasNonLocalPot = true ; if ( mynnloc == 0 ) mynnloc = new NonLocalECPComponent ; RealType rf = inFunc . rmax ( ) ; GridType * agrid = new LinearGrid < RealType > ; int ng = static_cast < int > ( rf * 100 ) + 1 ; agrid -> set ( 0.0 , rf , ng ) ; app_log ( ) << " " << angmom << " " << rf << std :: endl ; app_log ( ) << " " << rf << " " << ng << std :: endl ; std :: vector < RealType > pp_temp ( ng ) ; pp_temp [ 0 ] = inFunc ( 0 ) ; for ( int j = 1 ; j < ng ; j ++ ) { pp_temp [ j ] = inFunc . splint ( ( * agrid ) [ j ] ) ; } RadialPotentialType * app = new RadialPotentialType ( agrid , pp_temp ) ; app -> spline ( ) ; mynnloc -> add ( angmom , app ) ; lmax = std :: max ( lmax , angmom ) ; rmax = std :: max ( rmax , rf ) ; numnonloc ++ ; } if ( mynnloc ) { mynnloc -> setLmax ( lmax ) ; mynnloc -> setRmax ( rmax ) ; app_log ( ) << " " << rmax << std :: endl ; } } fin . close ( ) ; if ( mynnloc ) { nonLocalPot [ ig ] = mynnloc ; int numsgridpts = 0 ; std :: string fname = species + " " ; std :: ifstream fin ( fname . c_str ( ) , std :: ios_base :: in ) ; if ( ! fin ) { app_error ( ) << " " << fname << std :: endl ; exit ( -1 ) ; } PosType xyz ; RealType weight ; while ( fin >> xyz >> weight ) { mynnloc -> addknot ( xyz , weight ) ; numsgridpts ++ ; } mynnloc -> resize_warrays ( numsgridpts , numnonloc , lmax ) ; } } }
reimplement simple table format used by NonLocalPPotential
reimplement simple table format used by NonLocalPPotential
// Read Number of potentials (local and non) for this atom //local potential, input is rescale by -r/z //use 1e-2 resolution //will create LocalECPotential //will create NonLocalECPotential //get the value //cout << "Spherical grid : " << numsgridpts << " points" << std::endl; //species