more byb updates...EJB

Nwpw/nwpwlib/D3dB/Pneb.cpp

@@ -2773,6 +2773,23 @@ void Pneb::g_project_out_filled_below(double *psi, const int ms, const int k, do
 }
 
 
+/**************************************
+ *                                    *
+ *  Pneb::g_project_out_filled_above  *
+ *                                    *
+ **************************************/
+void Pneb::g_project_out_filled_above(double *psi, const int ms, const int k, double *Horb) 
+{
+   int ishift = ms*neq[0]*2*PGrid::npack(1);
+   for (auto ka=k+1; ka<neq[ms]; ++ka)
+   {
+     int indx = 2*PGrid::npack(1)*ka + ishift;
+     double w = -PGrid::cc_pack_dot(1,psi+indx,Horb);
+     PGrid::cc_pack_daxpy(1,w,psi+indx,Horb);
+   }
+}
+
+
 /*********************************
  *                               *
  *  Pneb::g_project_out_virtual  *

Nwpw/nwpwlib/D3dB/Pneb.hpp

@@ -242,6 +242,7 @@ class Pneb : public PGrid, public d1db {
    void g_project_out_filled(double *, const int, double *);
    void g_project_out_virtual(const int, const int *, const int, double *,  double *);
    void g_project_out_filled_below(double *, const int, const int, double *);
+   void g_project_out_filled_above(double *, const int, const int, double *);
 
    void g_norm(double *);
 

Nwpw/nwpwlib/utilities/nwpw_kerker.hpp

@@ -87,7 +87,9 @@ class nwpw_kerker {
       if (dokerker)
       {
          mygrid->rr_SMul(scal1,v,tmpv);
+         mygrid->r_zero_ends(tmpv);
          mygrid->rc_pfft3f(0,tmpv);
+         mygrid->c_pack(0,tmpv);
 
          int k1=0;
          int k2=0;
@@ -98,8 +100,10 @@ class nwpw_kerker {
             k1 += 2;
             k2 += 2;
          }
+         mygrid->c_unpack(0,tmpv);
          mygrid->cr_pfft3b(0,tmpv);
          mygrid->rr_copy(tmpv,v);
+         mygrid->r_zero_ends(v);
       }
    }
 };

Nwpw/nwpwlib/utilities/nwpw_scf_mixing.hpp

@@ -54,7 +54,8 @@ class nwpw_scf_mixing : public nwpw_kerker {
       n2ft3d = nsize0;
       nsize  = nsize0*ispin0;
       ispin  = ispin0;
-      std::cout << "nwpw_scf_mixing algorithm=" << algorithm << std::endl;
+
+      // std::cout << "nwpw_scf_mixing algorithm=" << algorithm << std::endl;
       /* simple mixing */
       if (algorithm==0)
       {
@@ -79,7 +80,7 @@ class nwpw_scf_mixing : public nwpw_kerker {
       /* Anderson mixing */
       else if (algorithm==3)
       {
-         rho_list = new (std::nothrow) double[nsize*3]();
+         rho_list = new (std::nothrow) double[nsize*4]();
          nwpw_scf_mixing_reset(rho_in);
       }
 
@@ -144,8 +145,12 @@ class nwpw_scf_mixing : public nwpw_kerker {
          *scf_error0= scf_error;
 
          // kerker stuff here
-
-         std::memcpy(rr,vnew,nsize*sizeof(double));
+         for (auto ms=0; ms<ispin; ++ms)
+            kerker_G(ff + ms*n2ft3d);
+
+         std::memcpy(vnew,rr,nsize*sizeof(double));
+         DAXPY_PWDFT(nsize,alpha,ff,one,vnew,one);  //vnew = vnew+alpha*f
+         std::memcpy(rr,vnew,nsize*sizeof(double)); //vm=vnew
       }
 
       /* Broyden mixing */
@@ -378,33 +383,84 @@ class nwpw_scf_mixing : public nwpw_kerker {
       /* Anderson density mixing */
       if (algorithm==3)
       {
-         double *rr = rho_list;
-         double *ff = rho_list + nsize;
-         double *tt = rho_list + 2*nsize;
+         if ((m==2) || (m==1) || (deltae>0.0))
+         {
+            double *rr = rho_list;
+            double *ff = rho_list + nsize;
+            double *tt = rho_list + 2*nsize;
 
-         std::memcpy(tt,vout,nsize*sizeof(double));
-         std::memcpy(ff,rr,nsize*sizeof(double));
-         std::memcpy(rr,vnew,nsize*sizeof(double));
+            // ff=vout-vm
+            std::memcpy(ff,vout,nsize*sizeof(double));
+            DAXPY_PWDFT(nsize,mrone,rr,one,ff,one);
 
+            double scf_error = DDOT_PWDFT(nsize,ff,one,ff,one);
+            parall->SumAll(1,scf_error);
+            scf_error = std::sqrt(scf_error);
+            *scf_error0 = scf_error;
+
+            for (auto ms=0; ms<ispin; ++ms)
+               kerker_G(ff + ms*n2ft3d);
+
+            std::memcpy(vnew,rr,nsize*sizeof(double));
+            DAXPY_PWDFT(nsize,alpha,ff,one,vnew,one);
+
+            std::memcpy(tt,vout,nsize*sizeof(double));
+            std::memcpy(ff,rr,nsize*sizeof(double));
+            std::memcpy(rr,vnew,nsize*sizeof(double));
+         }
+         else
+         {
+            double *rr = rho_list;           // vm 
+            double *ss = rho_list + nsize;   // vm1
+            double *tt = rho_list + 2*nsize; // vout1
+            double *ff = rho_list + 3*nsize; //kk
 
-         // ff=vout-vm
-         std::memcpy(ff,vout,nsize*sizeof(double));
-         DAXPY_PWDFT(nsize,mrone,rr,one,ff,one);
-         double scf_error = DDOT_PWDFT(nsize,ff,one,ff,one);
-         parall->SumAll(1,scf_error);
-         scf_error = std::sqrt(scf_error);
-         *scf_error0 = scf_error;
+            std::memcpy(ff,vout,nsize*sizeof(double));
 
-         //do ms=1,ispin
-         //   shift = (ms-1)*n2ft3d
-         //   call kerker_G(dbl_mb(ff_ptr+shift))
-         //end do
+            DAXPY_PWDFT(nsize,mrone,rr,one,ff,one);  // ff_ptr = vout - rr_ptr
+            DAXPY_PWDFT(nsize,mrone,ss,one,tt,one);  // tt_ptr = vout1 - vm1
 
-         std::memcpy(vnew,rr,nsize*sizeof(double));
-         DAXPY_PWDFT(nsize,alpha,ff,one,vnew,one);
-         std::memcpy(tt,vout,nsize*sizeof(double));
-         std::memcpy(ff,rr,nsize*sizeof(double));
-         std::memcpy(rr,vnew,nsize*sizeof(double));
+            double scf_error = DDOT_PWDFT(nsize,ff,one,ff,one);
+            parall->SumAll(1,scf_error);
+            scf_error = std::sqrt(scf_error);
+            *scf_error0 = scf_error;
+
+            for (auto ms=0; ms<ispin; ++ms)
+            {
+               int shift = ms*n2ft3d;
+
+               kerker_G(ff + shift);
+               kerker_G(tt + shift);
+
+               // generate beta 
+               double p00 = DDOT_PWDFT(n2ft3d,ff+shift,one,ff+shift,one); // p00 = <ff_ptr|ff_ptr>
+               double p01 = DDOT_PWDFT(n2ft3d,ff+shift,one,tt+shift,one); // p01 = <ff_ptr|tt_ptr>
+               double p11 = DDOT_PWDFT(n2ft3d,tt+shift,one,tt+shift,one); // p11 = <tt_ptr|tt_ptr>
+               parall->SumAll(1,p00);
+               parall->SumAll(1,p01);
+               parall->SumAll(1,p11);
+               double r00 = p00-2.00*p01+p11;
+               beta = (p00-p01)/r00;
+
+               if ((r00<0.0) || (beta<(-1.5))) beta = 1.0e-3;
+               if (beta>1.5)  beta = 1.5;
+
+               double ombeta = 1.0-beta;
+
+               //*** vnew = (1-beta)*vm + beta*vm1 + alpha*((1-beta)*fm + beta*fm1) ***
+               std::memcpy(vnew+shift,rr+shift,n2ft3d*sizeof(double));
+
+               DSCAL_PWDFT(n2ft3d,ombeta,vnew+shift,one);
+               DAXPY_PWDFT(n2ft3d,beta,ss+shift,one,vnew+shift,one);
+               DSCAL_PWDFT(n2ft3d,ombeta,ff+shift,one);
+               DAXPY_PWDFT(n2ft3d,beta,tt+shift,one,ff+shift,one);
+               DAXPY_PWDFT(n2ft3d,alpha,ff+shift,one,vnew+shift,one);
+            }
+            std::memcpy(ss,rr,  nsize*sizeof(double));
+            std::memcpy(rr,vnew,nsize*sizeof(double));
+            std::memcpy(tt,vout,nsize*sizeof(double));
+         }
+         ++m; 
       }
 
       /* local Thomas Fermi mixing */

Nwpw/pspw/lib/electron/psi_H.cpp

@@ -471,7 +471,7 @@ void psi_H_orb(Pneb *mygrid, Kinetic_Operator *myke, Pseudopotential *mypsp,
 
    /* apply k-space operators */
    myke->ke_orb(orb,Horb);
- 
+
    /* apply non-local PSP  - Expensive */
    mypsp->v_nonlocal_orb(orb, Horb);
 

Nwpw/pspw/lib/molecule/Molecule.cpp

@@ -293,7 +293,7 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
    double *g  = mygrid->c_pack_allocate(1);
    double *t  = mygrid->c_pack_allocate(1);
 
-   bool precondition = true;
+   bool precondition = false;
    bool done = false;
    double error0 = 0.0;
    double e0 = 0.0;
@@ -306,19 +306,20 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
    while (!done)
    {
       ++it;
+      error0 = std::abs(e0-eold);
       eold = e0;
 
       //calculate residual (steepest descent) direction for a single band
       psi_get_gradient(orb, vall+ms*n2ft3d, g);
       e0 = mygrid->cc_pack_dot(1,orb,g);
 
       e0 = -e0;
+      //std::cout << "it =" << it << " e0=" << e0 << std::endl;
 
 
-      double percent_error=0.0;
+      double percent_error = 0.0;
       if(error0>1.0e-11) percent_error = std::abs(e0-eold)/error0;
 
-
       precondition = (std::abs(e0-eold)>(sp*maxerror));
 
       done = ((it > maxit_orb) || (std::abs(e0-eold)<maxerror));
@@ -339,7 +340,6 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
       mygrid->cc_pack_copy(1,r1,t);
 
 
-      std::cout << "lmbda_r0=" << lmbda_r0 << std::endl;
       if (it>1) 
          if (!std::isnan(lmbda_r0)) 
             mygrid->cc_pack_daxpy(1,(lmbda_r1/lmbda_r0),t0,t);
@@ -359,8 +359,6 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
          de0 = mygrid->cc_pack_dot(1,t,t);
          de0 = 1.0/std::sqrt(de0);
          if (std::isnan(de0)) de0=0.0;
-         std::cout << "de0=" << de0 << std::endl;
-      std::cout << "A theta=" << theta << std::endl;
          mygrid->c_pack_SMul(1,de0,t);
          de0 = mygrid->cc_pack_dot(1,t,g);
 
@@ -377,7 +375,6 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
       de0 = -2.0*de0;
 
       psi_linesearch_update(e0,de0,&theta,vall+ms*n2ft3d,orb,t);
-      std::cout << "B theta=" << theta << std::endl;
 
       done = ((it > maxit_orb) ||  (std::abs(e0-eold) < maxerror));
       //done = true;
@@ -388,10 +385,11 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
    mygrid->c_pack_deallocate(g);
    mygrid->c_pack_deallocate(t);
 
-   *error_out = std::abs(e0-eold);
    e0         = -e0;
+   *error_out = (e0-eold);
 
    bool lprint = (mygrid->d3db::parall->is_master());
+   lprint = false;
    if (lprint) coutput << std::setw(12) << "orbital" << std::setw(4) << k+1
            << " current e=" << std::setw(10) << std::scientific << std::setprecision(3) << e0
            << " (error=" << std::setw(9) << std::scientific << std::setprecision(3) << (*error_out) << ")"
@@ -414,17 +412,19 @@ double Molecule::psi_KS_update(const int maxit_orb, const double maxerror,
                                double *error_out, std::ostream &coutput)
 
 {        
+   double esum = 0.0;
 
    for (auto ms=0; ms<ispin; ++ms)
    {
       int ishift = ms*neq[0]*2*mygrid->PGrid::npack(1);
-      for (auto i=0; i<neq[ms]; ++i)
+      for (auto i=neq[ms]-1; i>=0; --i)
       {
          int indx = 2*mygrid->PGrid::npack(1)*i + ishift;
          double *orb = psi + indx;
 
          // orthogonalize to lower orbitals
-         mygrid->g_project_out_filled_below(psi1, ms, i, orb);
+         //mygrid->g_project_out_filled_below(psi1, ms, i, orb);
+         mygrid->g_project_out_filled_above(psi1, ms, i, orb);
 
          // normalize
          double norm = mygrid->cc_pack_dot(1,orb,orb);
@@ -433,8 +433,11 @@ double Molecule::psi_KS_update(const int maxit_orb, const double maxerror,
 
          double e0 = psi_KS_update_orb(ms, i, maxit_orb, maxerror, perror, vall, orb,
                                        error_out, coutput);
+         esum += e0;
       }
    }
+
+   return esum;
 }        
 
 /********************************************

Nwpw/pspw/minimizer/cgsd_bybminimize2.cpp

@@ -55,53 +55,15 @@ double cgsd_bybminimize2(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
 
 
    // ion-ion energy 
-   double eion = mymolecule.eion();
-
- //  int sd_it = 10
- //  int cg_it = 10
- //   if (set_iterations) then
- //       it_in = iterations
- //       sd_it = 2 
- //       cg_it = 1 
- //     else
- //       it_in = control_it_in()*control_it_out()
- //       sd_it = 10
- //       cg_it = 10
- //     end if
-
-//      int maxit_orbs = control_ks_maxit_orbs() 
-      //precondition = control_precondition()
-      //ispin = control_ispin()
-      double deltav_old = 10.0;
-      double deltav     = 0.0;
-
- //     stalled       = false;
-  //    std::vector<double> deltae_history;
-   //   deltae_history.push_back(0.0);
-  //    deltae_history.push_back(0.0);
-   //   deltae_history.push_back(0.0);
-    //  deltae_history.push_back(0.0);
-      int stalled_count     = 0;
-      int sd_count          = 0;
-
- //     vall_in = 
-  //    vall_out = 
- //     vall_junk = 
- //     rho_in = 
+   //double eion = mymolecule.eion();
 
    //**********************
    //**** bybminimizer ****
    //**********************
-
-   //nwpw_scf_mixing scfmix(mygrid,g0,
-   //                       scf_algorithm,scf_alpha,diis_histories,
-   //                       mygrid->ispin,mygrid->n2ft3d,vall_out);
-
    double e0 = mymolecule.psi_KS_update(ks_it_in,ks_deltae,perror,vall,
                                         ispin, neq, mymolecule.psi1,
                                         deltae, std::cout);
 
-
    /* iniitialize blocked cg */
 
    // Making an extra call to electron.run and energy