Feature : add paw_nl_force (not tested yet) (#3279)

* update version of libpaw_interface

* Feature : add paw_nl_force

* update ut for paw

---------

Co-authored-by: wenfei-li <liwenfei@gmail.com>
Co-authored-by: Qianrui Liu <76200646+Qianruipku@users.noreply.github.com>

Author: 3 people

deps/libpaw_interface

@@ -19,6 +19,7 @@ class Paw_Atom
 
     //pass <psi|ptilde> from outside and saves it
     void set_ca(std::vector<std::complex<double>> & ca_in, const double weight_in);
+    void set_dca(std::vector<std::vector<std::complex<double>>> & dca_in, const double weight_in);
 
     void init_rhoij(); //set rhoij according to occupation number in xml file
 
@@ -50,6 +51,7 @@ class Paw_Atom
 
     std::vector<std::complex<double>> ca; //coefficients <psi|ptilde> for a given psi
     std::vector<std::vector<double>> rhoij; //on-site density matrix, upper triangular
+    std::vector<std::vector<std::complex<double>>> dca; // d/dR_I <psi|ptilde>
 
     std::vector<std::vector<double>> dij; //nonlocal pseudopotential strength
     std::vector<double> sij; //<phi|phi> - <phitilde|phitilde>

source/module_cell/module_paw/paw_atom.h

@@ -210,6 +210,21 @@ void Paw_Cell::set_paw_k(
     {
         gnorm[ipw] = std::sqrt(kpg[ipw][0]*kpg[ipw][0] + kpg[ipw][1]*kpg[ipw][1] + kpg[ipw][2]*kpg[ipw][2]) * tpiba;
     }
+
+    std::complex<double> i_cplx = (0.0,1.0);
+    // ikpg : i (k+G)
+    //if(GlobalV::CAL_FORCE || GlobalV::CAL_STRESS)
+    {
+        ikpg.resize(npw);
+        for(int ipw = 0; ipw < npw; ipw ++)
+        {
+            ikpg[ipw].resize(3);
+            for(int i = 0; i < 3; i ++)
+            {
+                ikpg[ipw][i] = kpg[ipw][i] * tpiba * i_cplx;
+            }
+        }
+    }
 }
 
 void Paw_Cell::set_isk(const int nk, const int * isk_in)
@@ -665,4 +680,81 @@ void Paw_Cell::paw_nl_psi(const int mode, const std::complex<double> * psi, std:
             }
         }
     }
+}
+
+void Paw_Cell::paw_nl_force(const std::complex<double> * psi, const double * epsilon, const double * weight, const int nbands , double * force)
+{
+    ModuleBase::TITLE("Paw_Cell","paw_nl_force");
+
+    for(int iband = 0; iband < nbands; iband ++)
+    {
+        for(int iat = 0; iat < nat; iat ++)
+        {
+            // ca : <ptilde(G)|psi(G)>
+            // = \sum_G [\int f(r)r^2j_l(r)dr] * [(-i)^l] * [ylm(\hat{G})] * [exp(-GR_I)] *psi(G)
+            // = \sum_ipw ptilde * fact * ylm * sk * psi (in the code below)
+            // This is what is called 'becp' in nonlocal pp
+            std::vector<std::complex<double>> ca;
+            std::vector<std::vector<std::complex<double>>> dca;
+
+            const int it = atom_type[iat];
+            const int nproj = paw_element_list[it].get_mstates();
+            const int proj_start = start_iprj_ia[iat];
+
+            ca.resize(nproj);
+            dca.resize(3);
+            for(int i = 0; i < 3; i ++)
+            {
+                dca[i].resize(nproj);
+            }
+
+            for(int iproj = 0; iproj < nproj; iproj ++)
+            {
+                ca[iproj] = 0.0;
+                dca[0][iproj] = 0.0;
+                dca[1][iproj] = 0.0;
+                dca[2][iproj] = 0.0;
+                
+                // consider use blas subroutine for this part later
+                for(int ipw = 0; ipw < npw; ipw ++)
+                {
+                    std::complex<double> overlp = psi[ipw] * std::conj(vkb[iproj+proj_start][ipw]);
+                    ca[iproj] += overlp;
+                    dca[0][iproj] += overlp * ikpg[ipw][0];
+                    dca[1][iproj] += overlp * ikpg[ipw][1];
+                    dca[2][iproj] += overlp * ikpg[ipw][2];
+                }
+            }
+
+#ifdef __MPI
+            Parallel_Reduce::reduce_pool(ca.data(), nproj);
+            Parallel_Reduce::reduce_pool(dca[0].data(), nproj);
+            Parallel_Reduce::reduce_pool(dca[1].data(), nproj);
+            Parallel_Reduce::reduce_pool(dca[2].data(), nproj);
+#endif
+            // sum_ij (D_ij - epsilon_n O_ij) ca_j
+            std::vector<std::complex<double>> v_ca;
+            v_ca.resize(nproj);
+
+            for(int iproj = 0; iproj < nproj; iproj ++)
+            {
+                v_ca[iproj] = 0.0;
+                for(int jproj = 0; jproj < nproj; jproj ++)
+                {
+                    double coeff = paw_atom_list[iat].get_dij()[current_spin][iproj*nproj+jproj] -
+                        paw_atom_list[iat].get_sij()[iproj*nproj+jproj] * epsilon[iband];
+                    v_ca[iproj] += coeff * ca[jproj];
+                }
+            }
+
+            // force += conjg(v_ca[iproj]) * d_ca[iproj]
+            // \sum_i ptilde_{iproj}(G) v_ca[iproj]
+            for(int iproj = 0; iproj < nproj; iproj ++)
+            {
+                force[iat*3] += (std::conj(v_ca[iproj]) * dca[0][iproj]).real();
+                force[iat*3+1] += (std::conj(v_ca[iproj]) * dca[1][iproj]).real();
+                force[iat*3+2] += (std::conj(v_ca[iproj]) * dca[2][iproj]).real();
+            }
+        }
+    }
 }

source/module_cell/module_paw/paw_cell.cpp

@@ -162,6 +162,9 @@ class Paw_Cell
     // I'd rather also calculate it once and save it
     std::vector<double> gnorm;
 
+    // i(k+G), used in force calculation
+    std::vector<std::vector<std::complex<double>>> ikpg;
+
     void set_ylm(const int npw_in, const double ** kpg);
 
     std::vector<int> isk;
@@ -176,6 +179,9 @@ class Paw_Cell
     // mode = 0 : V_{NL}|psi>, mode = 1 : (S+I)|psi>
     void paw_nl_psi(const int mode, const std::complex<double> * psi, std::complex<double> * vnlpsi);
 
+    
+    void paw_nl_force(const std::complex<double> * psi, const double * epsilon, const double * weight, const int nbands , double * force);
+
     // set by providing dij explicitly
     void set_dij(const int iat, double** dij_in){paw_atom_list[iat].set_dij(dij_in);}
     // set by extracting dij from libpaw_interface

source/module_cell/module_paw/paw_cell.h

@@ -45,3 +45,8 @@ AddTest(
 )
 add_dependencies(Test_Paw4 libpaw_interface)
 
+AddTest(
+  TARGET Test_Paw5
+  LIBS ${math_libs} base device
+  SOURCES test_paw5.cpp ../paw_element.cpp ../paw_sphbes.cpp ../paw_cell.cpp ../paw_atom.cpp
+)