Feature: add coulomb potential calculation

source/module_cell/atom_pseudo.cpp

@@ -240,6 +240,9 @@ void Atom_pseudo::bcast_atom_pseudo(void)
 // == end of pseudo_vl ==
 
 // == pseudo ==
+	if(nbeta == 0)
+		return;
+
 	if(GlobalV::MY_RANK != 0)
 	{
 		delete[] lll;

source/module_cell/atom_spec.cpp

@@ -130,6 +130,7 @@ void Atom::bcast_atom(void)
     Parallel_Common::bcast_int( nw );
     Parallel_Common::bcast_int( stapos_wf );
     Parallel_Common::bcast_string( label );
+    Parallel_Common::bcast_bool(coulomb_potential);
     if(GlobalV::MY_RANK!=0)
     {
         delete[] l_nchi;

source/module_cell/atom_spec.h

@@ -45,6 +45,9 @@ class Atom
 	double* angle1;//spin angle, added by zhengdy-soc
 	double* angle2;
     ModuleBase::Vector3<double> *m_loc_;
+    // Coulomb potential v(r) = z/r
+    // It is a local potentail, and has no non-local potential parts.
+    bool coulomb_potential = false;
 
     void print_Atom(std::ofstream &ofs);
     void update_force(ModuleBase::matrix &fcs);

source/module_cell/pseudo.cpp

@@ -30,7 +30,8 @@ void pseudo::set_pseudo(const Pseudopot_upf& upf)
 	this->set_pseudo_atom(upf);
 	this->set_pseudo_vl(upf);
 
-
+	if (nbeta == 0)
+		return;
 	delete[] lll;
 	lll = new int[nbeta];
 	assert(lll != 0);
@@ -218,7 +219,6 @@ void pseudo::set_pseudo_atom(const Pseudopot_upf &upf)
 	delete[] rho_atc;
 	rho_atc = new double[mesh];
 	assert(rho_atc != 0);
-	ModuleBase::GlobalFunc::ZEROS(rho_atc, mesh);
 
 	for (int i = 0;i < nchi;i++)
 	{
@@ -292,7 +292,10 @@ void pseudo::set_pseudo_vl(const Pseudopot_upf &upf)
 	vloc_at = new double[mesh];
 	assert(vloc_at != 0);
 
-	for (int i = 0;i < mesh;i++)
+    if (upf.coulomb_potential)
+        return;
+
+    for (int i = 0;i < mesh;i++)
 	{
 		vloc_at[i] = upf.vloc[i];
 	}

source/module_cell/pseudo.h

@@ -1,89 +1,89 @@
 #ifndef PSEUDO_H
 #define PSEUDO_H
 
-#include "read_pp.h"
 #include "../module_io/output.h"
+#include "read_pp.h"
 
 //-----------------------------------------
 // read in pseudopotentials
 // mohan update 2021-05-01
 //-----------------------------------------
 class pseudo
 {
-	public:
-      pseudo();
-      ~pseudo();
-
-      //<PP_HEADER>
-      bool has_so;         // if .true. includes spin-orbit
-      int nv;              // UPF file version number
-      std::string psd;     // Element label
-      std::string pp_type; // Pseudo type ( NC or US )
-      bool tvanp;          // .true. if Ultrasoft
-      bool nlcc;           // Non linear core corrections(bool)
-      std::string xc_func; // Exch-Corr type
-      int zv;              // z valence
-      double etotps;       // total energy
-      double ecutwfc;      // suggested cut-off for wfc
-      double ecutrho;      // suggested cut-off for rho
-      int lmax;            // maximum angular momentum component
-      int mesh;            // number of point in the radial mesh
-      int nchi;            // nwfc,number of wavefunctions
-      int nbeta;           // number of projectors
-      int nqlc;            // number of angular momenta in Q
-      int kkbeta;          // kkbeta, point where the beta are zero
-
-      std::string* els = nullptr; // els[nchi]
-      int* lchi = nullptr;        // lchi[nchi]
-      double* oc = nullptr;       // oc[nchi]
-
-      double* jjj = nullptr;  // total angual momentum, jjj[nbeta]
-      double* jchi = nullptr; // jchi(nwfc), added by zhengdy-soc
-      int* nn = nullptr;
-
-      // Local pseudopotentials
-      double* vloc_at = nullptr; // [mesh], local potential( = pseudopot_upf.vloc )
-
-      // <PP_MESH>
-      double* r = nullptr;   // radial logaritmic mesh, r[0:mesh-1]
-      double* rab = nullptr; // derivative of the radial mesh, rab[0:mesh-1]
-
-      //<PP_NLCC>
-      double* rho_atc = nullptr; // radial core charge density, rho_atc[0:mesh-1]
-
-      //<PP_RHOATOM>
-      double* rho_at = nullptr; // radial atomic charge density, rho_at[0:mesh-1]
-
-      // <PP_PSWFC>
-      ModuleBase::matrix chi; // radial atomic orbitals, chi(nchi, mesh)
-
-      // other
-      int msh;     // number of points up to rcut
-      double rcut; // cut-off radius
-
-      // <PP_BETA>
-      int* lll = nullptr; // lll(nbeta), angular momentum of the beta function
-
-      // <PP_DIJ>
-      ModuleBase::matrix dion;  // dion(nbeta,nbeta)
-      ModuleBase::matrix betar; // (nbeta, mesh), radial beta_{mu} functions
-
-      // other
-      int nh; // number of beta functions per atomic type
-
-      // uspp
-      ModuleBase::realArray qfuncl; // qfuncl(2*lmax+1,nbeta*(nbeta+1)/2,mesh) Q_{mu,nu}(|r|) function for |r|> r_L
-      ModuleBase::matrix qqq;       // qqq(nbeta,nbeta) q_{mu,nu}
-
-      void set_pseudo_h(const Pseudopot_upf& upf);
-      void set_pseudo_atom(const Pseudopot_upf& upf);
-      void set_pseudo_vl(const Pseudopot_upf& upf);
-      void set_pseudo(const Pseudopot_upf& upf);
-
-      void print_pseudo_h(std::ofstream& ofs);
-      void print_pseudo_atom(std::ofstream& ofs);
-      void print_pseudo_vl(std::ofstream& ofs);
-      void print_pseudo(std::ofstream& ofs);
+  public:
+    pseudo();
+    ~pseudo();
+
+    //<PP_HEADER>
+    bool has_so = false;  // if .true. includes spin-orbit
+    int nv = 0;           // UPF file version number
+    std::string psd;      // Element label
+    std::string pp_type;  // Pseudo type ( NC or US )
+    bool tvanp = false;   // .true. if Ultrasoft
+    bool nlcc = false;    // Non linear core corrections(bool)
+    std::string xc_func;  // Exch-Corr type
+    int zv = 0;           // z valence
+    double etotps = 0.0;  // total energy
+    double ecutwfc = 0.0; // suggested cut-off for wfc
+    double ecutrho = 0.0; // suggested cut-off for rho
+    int lmax = 0;         // maximum angular momentum component
+    int mesh = 0;         // number of point in the radial mesh
+    int nchi = 0;         // nwfc,number of wavefunctions
+    int nbeta = 0;        // number of projectors
+    int nqlc = 0;         // number of angular momenta in Q
+    int kkbeta = 0;       // kkbeta, point where the beta are zero
+
+    std::string* els = nullptr; // els[nchi]
+    int* lchi = nullptr;        // lchi[nchi]
+    double* oc = nullptr;       // oc[nchi]
+
+    double* jjj = nullptr;  // total angual momentum, jjj[nbeta]
+    double* jchi = nullptr; // jchi(nwfc), added by zhengdy-soc
+    int* nn = nullptr;
+
+    // Local pseudopotentials
+    double* vloc_at = nullptr; // [mesh], local potential( = pseudopot_upf.vloc )
+
+    // <PP_MESH>
+    double* r = nullptr;   // radial logaritmic mesh, r[0:mesh-1]
+    double* rab = nullptr; // derivative of the radial mesh, rab[0:mesh-1]
+
+    //<PP_NLCC>
+    double* rho_atc = nullptr; // radial core charge density, rho_atc[0:mesh-1]
+
+    //<PP_RHOATOM>
+    double* rho_at = nullptr; // radial atomic charge density, rho_at[0:mesh-1]
+
+    // <PP_PSWFC>
+    ModuleBase::matrix chi; // radial atomic orbitals, chi(nchi, mesh)
+
+    // other
+    int msh = 0;       // number of points up to rcut
+    double rcut = 0.0; // cut-off radius
+
+    // <PP_BETA>
+    int* lll = nullptr; // lll(nbeta), angular momentum of the beta function
+
+    // <PP_DIJ>
+    ModuleBase::matrix dion;  // dion(nbeta,nbeta)
+    ModuleBase::matrix betar; // (nbeta, mesh), radial beta_{mu} functions
+
+    // other
+    int nh = 0; // number of beta functions per atomic type
+
+    // uspp
+    ModuleBase::realArray qfuncl; // qfuncl(2*lmax+1,nbeta*(nbeta+1)/2,mesh) Q_{mu,nu}(|r|) function for |r|> r_L
+    ModuleBase::matrix qqq;       // qqq(nbeta,nbeta) q_{mu,nu}
+
+    void set_pseudo_h(const Pseudopot_upf& upf);
+    void set_pseudo_atom(const Pseudopot_upf& upf);
+    void set_pseudo_vl(const Pseudopot_upf& upf);
+    void set_pseudo(const Pseudopot_upf& upf);
+
+    void print_pseudo_h(std::ofstream& ofs);
+    void print_pseudo_atom(std::ofstream& ofs);
+    void print_pseudo_vl(std::ofstream& ofs);
+    void print_pseudo(std::ofstream& ofs);
 };
 
 #endif // PSEUDO_H

source/module_cell/read_atoms.cpp

@@ -65,6 +65,10 @@ int UnitCell::read_atom_species(std::ifstream &ifa, std::ofstream &ofs_running)
 					{
 						pseudo_type[i] = one_string;
 					}
+					else if (one_string == "1/r")
+					{
+						atoms[i].coulomb_potential = true;
+					}
 					else
 					{
 						GlobalV::ofs_warning << "unrecongnized pseudopotential type: " << one_string << ", check your STRU file." << std::endl;

source/module_cell/read_cell_pseudopots.cpp

@@ -19,6 +19,7 @@ void UnitCell::read_cell_pseudopots(const std::string &pp_dir, std::ofstream &lo
 	for (int i = 0;i < ntype;i++)
 	{
 		Pseudopot_upf upf;
+        upf.coulomb_potential=this->atoms[i].coulomb_potential;
 
 		// mohan update 2010-09-12	
 		int error = 0;
@@ -39,11 +40,13 @@ void UnitCell::read_cell_pseudopots(const std::string &pp_dir, std::ofstream &lo
                 // average pseudopotential if needed
                 error_ap = upf.average_p(GlobalV::soc_lambda); // added by zhengdy 2020-10-20
             }
+            this->atoms[i].coulomb_potential = upf.coulomb_potential;
         }
 
 #ifdef __MPI
 		Parallel_Common::bcast_int(error);
 		Parallel_Common::bcast_int(error_ap);
+        Parallel_Common::bcast_bool(atoms[i].coulomb_potential);
 #endif
 
 		if(error_ap) 

source/module_cell/read_pp.h

@@ -35,7 +35,7 @@ class Pseudopot_upf
     double etotps;            // total energy
     double ecutwfc;           // suggested cut-off for wfc
     double ecutrho;           // suggested cut-off for rho
-    int nv;                   // UPF file version number
+    int nv = 0;               // UPF file version number
     int lmax;                 // maximum angular momentum component in beta
     int lmax_rho;             // maximum angular momentum component in rho (should be 2*lmax)
     int nwfc;                 // number of wavefunctions
@@ -61,6 +61,7 @@ class Pseudopot_upf
     double* rab = nullptr;           // rab(mesh) dr(x)/dx (x=linear grid)
     double* rho_atc = nullptr;       // rho_atc(mesh) atomic core charge "Nonlinear Core Correction"
     double* vloc = nullptr;          // vloc(mesh) local atomic potential
+    bool coulomb_potential = false;  // coulomb potentail : z/r
     ModuleBase::matrix chi;          // chi(nwfc,mesh) atomic wavefcts
     double* rho_at = nullptr;        // rho_at(mesh) atomic charge
     int* lll = nullptr;              // lll(nbeta):angular momentum of projector i

source/module_cell/read_pp_upf100.cpp

@@ -81,13 +81,16 @@ int Pseudopot_upf::read_pseudo_upf(std::ifstream &ifs)
 		ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_NLCC>");
 	}
 
-	// Search for Local potential
-	ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_LOCAL>");
-	//---------------------
-	// call member function
-	//---------------------
-	read_pseudo_local(ifs);
-	ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_LOCAL>");
+    if (!this->coulomb_potential)
+    {
+        // Search for Local potential
+        ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_LOCAL>");
+        //---------------------
+        // call member function
+        //---------------------
+        read_pseudo_local(ifs);
+        ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_LOCAL>");
+    }
 
 	// Search for Nonlocal potential
 	ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_NONLOCAL>");
@@ -150,6 +153,11 @@ void Pseudopot_upf::read_pseudo_header(std::ifstream &ifs)
 	{
 		this->tvanp = false;
 	}
+	else if(pp_type == "1/r")
+	{
+		this->tvanp = false;
+		this->coulomb_potential = true;
+	}
 	else
 	{
 		// A bug here!!! can't quit together.
@@ -202,6 +210,12 @@ void Pseudopot_upf::read_pseudo_header(std::ifstream &ifs)
 	{
 		ifs >> els[i] >> this->lchi[i] >> this->oc[i];
 	}
+	if (this->coulomb_potential)
+	{
+		this->nbeta = 0;
+        this->lmax = 0;
+        this->lloc = 0;
+	}
 	return;
 }
 

source/module_cell/read_pp_upf201.cpp

@@ -40,22 +40,25 @@ int Pseudopot_upf::read_pseudo_upf201(std::ifstream &ifs)
     //--------------------------------------
     //-              PP_LOCAL              -
     //--------------------------------------
-    if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_LOCAL", true, false))
+    if (!this->coulomb_potential)
     {
-        ifs.ignore(150, '>'); // skip type, size, columns and so on.
-    }
-    else
-    {
-        ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_LOCAL>");
-    }
-    delete[] vloc;
-    this->vloc = new double[mesh];
-    ModuleBase::GlobalFunc::ZEROS(vloc, mesh);
-    for (int ir = 0; ir < mesh; ir++)
-    {
-        ifs >> this->vloc[ir];
+        if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_LOCAL", true, false))
+        {
+            ifs.ignore(150, '>'); // skip type, size, columns and so on.
+        }
+        else
+        {
+            ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_LOCAL>");
+        }
+        delete[] vloc;
+        this->vloc = new double[mesh];
+        ModuleBase::GlobalFunc::ZEROS(vloc, mesh);
+        for (int ir = 0; ir < mesh; ir++)
+        {
+            ifs >> this->vloc[ir];
+        }
+        ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_LOCAL>");
     }
-    ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_LOCAL>");
 
     //--------------------------------------
     //-            PP_NONLOCAL             -
@@ -251,8 +254,7 @@ void Pseudopot_upf::read_pseudo_upf201_header(std::ifstream& ifs)
         {
             if (val[ip] == "T" || val[ip] == "TRUE" || val[ip] == "True" || val[ip] == "true")
             {
-                ModuleBase::WARNING_QUIT("Pseudopot_upf::read_pseudo_upf201_header",
-                                         "COULOMB POTENTIAL IS NOT SUPPORTED");
+                this->coulomb_potential = true;
             }
         }
         else if (name[ip] == "has_so")
@@ -337,6 +339,12 @@ void Pseudopot_upf::read_pseudo_upf201_header(std::ifstream& ifs)
             ModuleBase::WARNING("PP_HEADRER reading", warningstr);
         }
     }
+    if (this->coulomb_potential)
+    {
+        this->nbeta = 0;
+        this->lmax = 0;
+        this->lloc = 0;
+    }
 }
 
 void Pseudopot_upf::read_pseudo_upf201_mesh(std::ifstream& ifs)
@@ -422,6 +430,10 @@ void Pseudopot_upf::read_pseudo_upf201_mesh(std::ifstream& ifs)
 void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs)
 {
     ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_NONLOCAL>");
+    if (nbeta == 0)
+    {
+        return;
+    }
     std::string word;
     std::string name[50];
     std::string val[50];