added Thomas-Fermi density mixing...EJB

Nwpw/nwpwlib/Control/Control2.cpp

@@ -503,6 +503,10 @@ Control2::Control2(const int np0, const std::string rtdbstring)
    if (rtdbjson["nwpw"]["scf_alpha"].is_number_float())
       pscf_alpha = rtdbjson["nwpw"]["scf_alpha"];
 
+   pscf_beta = 0.25;
+   if (rtdbjson["nwpw"]["scf_beta"].is_number_float())
+      pscf_beta = rtdbjson["nwpw"]["scf_beta"];
+
    pkerker_g0 = 0.0;
    if (rtdbjson["nwpw"]["kerker_g0"].is_number_float())
       pkerker_g0 = rtdbjson["nwpw"]["kerker_g0"];

Nwpw/nwpwlib/Control/Control2.hpp

@@ -19,7 +19,7 @@ class Control2 {
 
    std::string myrtdbstring, xcstring;
    double punita[9], ptolerances[3], pscaling[2];
-   double ptime_step, pfake_mass, pscf_alpha, pecut, pwcut, prcut;
+   double ptime_step, pfake_mass, pscf_alpha, pscf_beta, pecut, pwcut, prcut;
    double pkerker_g0,pfractional_kT,pfractional_temperature,pfractional_alpha;
    double pbo_time_step;
    double ptotal_charge;
@@ -186,6 +186,7 @@ class Control2 {
    double Sprecondition() { return psprecondition;}
 
    double scf_alpha() { return pscf_alpha; }
+   double scf_beta() { return pscf_beta; }
    double kerker_g0() { return pkerker_g0; }
    double fractional_kT() { return pfractional_kT; }
    double fractional_temperature() { return pfractional_temperature; }

Nwpw/nwpwlib/parse/parse_pwdft.cpp

@@ -1561,12 +1561,13 @@ static json parse_nwpw(json nwpwjson, int *curptr,
        if (mystring_contains(line, "johnson-pulay")) nwpwjson["scf_algorithm"] = 2;
        if (mystring_contains(line, "diis"))          nwpwjson["scf_algorithm"] = 2;
        if (mystring_contains(line, "anderson"))      nwpwjson["scf_algorithm"] = 3;
-       if (mystring_contains(line, "TF"))            nwpwjson["scf_algorithm"] = 4;
+       if (mystring_contains(line, "thomas-fermi"))  nwpwjson["scf_algorithm"] = 4;
        if (mystring_contains(line, "alpha")) 
           nwpwjson["scf_alpha"] = std::stod(mystring_trim(mystring_split(line, "alpha")[1]));
+       if (mystring_contains(line, "beta")) 
+          nwpwjson["scf_beta"] = std::stod(mystring_trim(mystring_split(line, "beta")[1]));
        if (mystring_contains(line, "kerker")) 
           nwpwjson["kerker_g0"] = std::stod(mystring_trim(mystring_split(line, "kerker")[1]));
-
        if (mystring_contains(line, " iterations")) 
           nwpwjson["ks_maxit_orb"] = std::stoi(mystring_trim(mystring_split(line, " iterations")[1]));
        if (mystring_contains(line, "outer_iterations")) 

Nwpw/nwpwlib/utilities/nwpw_scf_mixing.hpp

@@ -57,12 +57,13 @@ class nwpw_scf_mixing : public nwpw_kerker {
     *        nwpw_scf_mixing::nwpw_scf_mixing       *
     *                                               *
     *************************************************/
-   nwpw_scf_mixing(PGrid *mygrid0, const double g0, const int algorithm0, const double alpha0, const int max_m0, 
+   nwpw_scf_mixing(PGrid *mygrid0, const double g0, const int algorithm0, const double alpha0, const double beta0, const int max_m0, 
                    const int ispin0, const int nsize0, double *rho_in) : nwpw_kerker(mygrid0, g0)
    {
       parall = mygrid0->d3db::parall;
       algorithm = algorithm0;
       alpha = alpha0;
+      beta  = beta0;
       max_m = max_m0;
       n2ft3d = nsize0;
       nsize  = nsize0*ispin0;
@@ -101,6 +102,7 @@ class nwpw_scf_mixing : public nwpw_kerker {
       /* local Thomas-Fermi mixing */
       else if (algorithm==4)
       {
+         rho_list = new (std::nothrow) double[nsize*3]();
          reset_mix(rho_in);
       }
       std::memcpy(rho_list, rho_in, nsize*sizeof(double));
@@ -512,14 +514,34 @@ class nwpw_scf_mixing : public nwpw_kerker {
       /* local Thomas Fermi mixing */
       if (algorithm==4)
       {
+         const double twothirds = 2.0/3.0;
          double *rr = rho_list;
-         double *ss = rho_list + nsize;
-         double *tt = rho_list + 2*nsize;
-         double *ff = rho_list + 3*nsize;
+         double *ff = rho_list+nsize;
+         double *tf = rho_list+2*nsize;
+         std::memcpy(ff,rr,nsize*sizeof(double));
+
+         // compute the  residual ff = vout - vold 
+         DAXPY_PWDFT(nsize,mrone,vout,one,ff,one);
+         DSCAL_PWDFT(nsize,mrone,ff,one);
+
+         // scf_error = sqrt(<ff|ff>)
+         double scf_error = DDOT_PWDFT(nsize,ff,one,ff,one);
+         *scf_error0 = std::sqrt(parall->SumAll(1,scf_error))/((double) nsize);
+         //scf_error = std::sqrt(scf_error);
 
-         std::memcpy(ss,rr,nsize*sizeof(double));
-         std::memcpy(rr,vnew,nsize*sizeof(double));
-         std::memcpy(tt,vout,nsize*sizeof(double));
+         // Apply Kerker filter to smooth residual ff
+         for (auto ms=0; ms<ispin; ++ms)
+            kerker_G(ff + ms*n2ft3d);
+
+         // Apply TF mixing
+         // tf = ff/(1+alpha*rho(n-1)**(2/3))= ff/(1+alpha*rr**
+         // rho(n) = rho(n-1) + beta*tf
+         for (auto i=0; i<nsize; ++i)
+            tf[i] = ff[i]/(1.0 + alpha*std::pow(rr[i],twothirds));
+
+         std::memcpy(vnew,rr,nsize*sizeof(double)); // 
+         DAXPY_PWDFT(nsize,beta,tf,one,vnew,one);  // vnew(n) = rho(n-1) + beta*tf
+         std::memcpy(rr,vnew,nsize*sizeof(double)); //vm=vnew
       }
    }
 };

Nwpw/pspw/minimizer/cgsd_bybminimize.cpp

@@ -13,7 +13,7 @@
 #include "Pneb.hpp"
 #include "util_date.hpp"
 #include "util_linesearch.hpp"
-#include "nwpw_scf_mixing.hpp"
+//#include "nwpw_scf_mixing.hpp"
 
 namespace pwdft {
 
@@ -81,12 +81,6 @@ double cgsd_bybminimize(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
    //**** bybminimizer ****
    //**********************
 
-   nwpw_scf_mixing scfmix(mygrid,g0,
-                          scf_algorithm,scf_alpha,diis_histories,
-                          mygrid->ispin,mygrid->n2ft3d,vall_out);
-
-
-
    /* iniitialize blocked cg */
 
    // Making an extra call to electron.run and energy

Nwpw/pspw/minimizer/cgsd_energy.cpp

@@ -92,6 +92,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
    double tole = control.tolerances(0);
    double tolc = control.tolerances(1);
    double scf_alpha = control.scf_alpha();
+   double scf_beta  = control.scf_beta();
    double kerker_g0 = control.kerker_g0();
    int diis_histories = control.diis_histories();
    int scf_algorithm = control.scf_algorithm();
@@ -343,7 +344,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
       double scf_error;
 
       nwpw_scf_mixing scfmix(mygrid,kerker_g0,
-                             scf_algorithm,scf_alpha,diis_histories,
+                             scf_algorithm,scf_alpha,scf_beta,diis_histories,
                              mygrid->ispin,mygrid->n2ft3d,mymolecule.rho1);
 
       while ((icount < (it_out*it_in)) && (!converged))

Nwpw/pspw/minimizer/pspw_bomd.cpp

@@ -406,10 +406,13 @@ int pspw_bomd(MPI_Comm comm_world0,std::string &rtdbstring,std::ostream &coutput
            if (control.scf_algorithm()==2) coutput << "     SCF algorithm        = Johnson-Pulay mixing"
                                                    << " (" << Ifmt(3) <<  control.diis_histories() << " histories)\n";
            if (control.scf_algorithm()==3) coutput << "     SCF algorithm        = Anderson mixing\n";
-           if (control.scf_algorithm()==4) coutput << "     SCF algorithm        = TF mixing\n";
+           if (control.scf_algorithm()==4) coutput << "     SCF algorithm        = Thomas-Fermi mixing\n";
            if (control.minimizer()==5) coutput << "     SCF mixing type      = potential\n";
            if (control.minimizer()==8) coutput << "     SCF mixing type      = density\n";
-           coutput << "     SCF mixing parameter = " << control.scf_alpha() << std::endl;
+           if (control.scf_algorithm()==4)
+              coutput << "     SCF mixing parameters: alpha=" << control.scf_alpha() << " beta=" << control.scf_beta() << std::endl;
+           else
+              coutput << "     SCF mixing parameter = " << control.scf_alpha() << std::endl;
            coutput << "     Kerker damping       = " << control.kerker_g0() << std::endl;
            coutput << std::endl;
            if (control.fractional())

Nwpw/pspw/minimizer/pspw_geovib.cpp

@@ -345,35 +345,37 @@ int pspw_geovib(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream &cou
               << Ffmt(12,8) << myewald.rsalpha() << " rs =" << Ffmt(12,8) << myewald.rs()
               << ")" << std::endl;
 
-      if (flag > 0) {
-        coutput << std::endl;
-        coutput << " technical parameters:\n";
-        if (control.nolagrange()) coutput << "      diabling Lagrange multiplier" << std::endl;
-        if (control.io_buffer()) coutput << "      using io buffer " << std::endl;
-        coutput << "      fixed step: time step =" << Ffmt(12, 2)
-                << control.time_step() << "  ficticious mass =" << Ffmt(12, 2)
-                << control.fake_mass() << std::endl;
-        coutput << "      tolerance =" << Efmt(12, 3) << control.tolerances(0)
-                << " (energy) " << Efmt(12, 3) << control.tolerances(1)
-                << " (density) " << Efmt(12, 3) << control.tolerances(2)
-                << " (ion)\n";
-        coutput << "      max iterations = " << Ifmt(10)
-                << control.loop(0) * control.loop(1) << " (" << Ifmt(5)
-                << control.loop(0) << " inner " << Ifmt(5) << control.loop(1)
-                << " outer)\n";
-        if (control.minimizer() == 1)
-          coutput << "      minimizer = Grassmann conjugate gradient\n";
-        if (control.minimizer() == 2)
-          coutput << "      minimizer = Grassmann lmbfgs\n";
-        if (control.minimizer() == 4)
-          coutput << "      minimizer = Stiefel conjugate gradient\n";
-        if (control.minimizer() == 5)
-          coutput << "      minimizer = scf (potential)\n";
-        if (control.minimizer() == 7)
-          coutput << "      minimizer = Stiefel lmbfgs\n";
-        if (control.minimizer() == 8)
-          coutput << "      minimizer = scf (density)\n";
-        if ((control.minimizer() == 5) || (control.minimizer() == 8)) {
+      if (flag > 0) 
+      {
+         coutput << std::endl;
+         coutput << " technical parameters:\n";
+         if (control.nolagrange()) coutput << "      diabling Lagrange multiplier" << std::endl;
+         if (control.io_buffer()) coutput << "      using io buffer " << std::endl;
+         coutput << "      fixed step: time step =" << Ffmt(12, 2)
+                 << control.time_step() << "  ficticious mass =" << Ffmt(12, 2)
+                 << control.fake_mass() << std::endl;
+         coutput << "      tolerance =" << Efmt(12, 3) << control.tolerances(0)
+                 << " (energy) " << Efmt(12, 3) << control.tolerances(1)
+                 << " (density) " << Efmt(12, 3) << control.tolerances(2)
+                 << " (ion)\n";
+         coutput << "      max iterations = " << Ifmt(10)
+                 << control.loop(0) * control.loop(1) << " (" << Ifmt(5)
+                 << control.loop(0) << " inner " << Ifmt(5) << control.loop(1)
+                 << " outer)\n";
+         if (control.minimizer() == 1)
+            coutput << "      minimizer = Grassmann conjugate gradient\n";
+         if (control.minimizer() == 2)
+            coutput << "      minimizer = Grassmann lmbfgs\n";
+         if (control.minimizer() == 4)
+            coutput << "      minimizer = Stiefel conjugate gradient\n";
+         if (control.minimizer() == 5)
+            coutput << "      minimizer = scf (potential)\n";
+         if (control.minimizer() == 7)
+            coutput << "      minimizer = Stiefel lmbfgs\n";
+         if (control.minimizer() == 8)
+            coutput << "      minimizer = scf (density)\n";
+         if ((control.minimizer() == 5) || (control.minimizer() == 8)) 
+         {
             coutput << std::endl;
             coutput << " Kohn-Sham scf parameters:\n";
             coutput << "     Kohn-Sham algorithm  = conjugate gradient\n";
@@ -384,10 +386,13 @@ int pspw_geovib(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream &cou
             if (control.scf_algorithm()==2) coutput << "     SCF algorithm        = Johnson-Pulay mixing"
                                                     << " (" << Ifmt(3) <<  control.diis_histories() << " histories)\n";
             if (control.scf_algorithm()==3) coutput << "     SCF algorithm        = Anderson mixing\n";
-            if (control.scf_algorithm()==4) coutput << "     SCF algorithm        = TF mixing\n";
+            if (control.scf_algorithm()==4) coutput << "     SCF algorithm        = Thomas-Fermi mixing\n";
             if (control.minimizer()==5) coutput << "     SCF mixing type      = potential\n";
             if (control.minimizer()==8) coutput << "     SCF mixing type      = density\n";
-            coutput << "     SCF mixing parameter = " << control.scf_alpha() << std::endl;
+            if (control.scf_algorithm()==4)
+               coutput << "     SCF mixing parameters: alpha=" << control.scf_alpha() << " beta=" << control.scf_beta() << std::endl;
+            else
+               coutput << "     SCF mixing parameter = " << control.scf_alpha() << std::endl;
             coutput << "     Kerker damping       = " << control.kerker_g0() << std::endl;
             coutput << std::endl;
             if (control.fractional())

Nwpw/pspw/minimizer/pspw_minimizer.cpp

@@ -322,8 +322,8 @@ int pspw_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream &
       {
          coutput << std::endl;
          coutput << " technical parameters:\n";
-        if (control.nolagrange()) coutput << "      disabling Lagrange multiplier " << std::endl;
-        if (control.io_buffer()) coutput << "      using io buffer " << std::endl;
+         if (control.nolagrange()) coutput << "      disabling Lagrange multiplier " << std::endl;
+         if (control.io_buffer()) coutput << "      using io buffer " << std::endl;
          coutput << "      fixed step: time step =" << Ffmt(12, 2)
                  << control.time_step() << "  ficticious mass =" << Ffmt(12, 2)
                  << control.fake_mass() << std::endl;
@@ -353,10 +353,13 @@ int pspw_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream &
             if (control.scf_algorithm()==2) coutput << "     SCF algorithm        = Johnson-Pulay mixing"
                                                     << " (" << Ifmt(3) <<  control.diis_histories() << " histories)\n";
             if (control.scf_algorithm()==3) coutput << "     SCF algorithm        = Anderson mixing\n";
-            if (control.scf_algorithm()==4) coutput << "     SCF algorithm        = TF mixing\n";
+            if (control.scf_algorithm()==4) coutput << "     SCF algorithm        = Thomas-Fermi mixing\n";
             if (control.minimizer()==5) coutput << "     SCF mixing type      = potential\n";
             if (control.minimizer()==8) coutput << "     SCF mixing type      = density\n";
-            coutput << "     SCF mixing parameter = " << control.scf_alpha() << std::endl;
+            if (control.scf_algorithm()==4) 
+               coutput << "     SCF mixing parameters: alpha=" << control.scf_alpha() << " beta=" << control.scf_beta() << std::endl;
+            else
+               coutput << "     SCF mixing parameter = " << control.scf_alpha() << std::endl;
             coutput << "     Kerker damping       = " << control.kerker_g0() << std::endl;
             coutput << std::endl;
             if (control.fractional())