Unify MLWF and dipole print format

src/BOSampleStepper.cpp

@@ -290,16 +290,14 @@ void BOSampleStepper::step(int niter)
   {
     // MLWF can be computed at the gamma point only
     // There must be a single k-point, and it must be gamma
-    if ( wf.nkp() > 1 || ( wf.nkp()==1 && wf.kpoint(0) != D3vector(0,0,0) ) )
-    {
-      if ( onpe0 )
-      {
-        cout << " BOSampleStepper::step: MLWF can be computed at k=0 only"
-             << endl;
-        cout << " BOSampleStepper::step: cannot run" << endl;
-      }
-      return;
-    }
+
+    // Check that only the k=0 point is used
+    if ( wf.nkp()>1 || length(wf.kpoint(0)) != 0.0 )
+      throw runtime_error("BOSampleStepper: can compute MLWF at k=0 only");
+
+    // Check if nempty > 0
+    if ( wf.nempty() > 0 )
+      throw runtime_error("BOSampleStepper: cannot compute MLWF with nempty>0");
 
     for ( int isp_loc = 0; isp_loc < wf.nsp_loc(); ++isp_loc )
     {
@@ -1174,13 +1172,25 @@ void BOSampleStepper::step(int niter)
         if ( onpe0 )
         {
           D3vector idipole = atoms.dipole();
-          cout << setprecision(6);
-          cout << " <ionic_dipole> " << idipole
-               << " </ionic_dipole>" << endl;
-          cout << " <total_dipole> " << idipole + edipole_sum
-               << " </total_dipole>" << endl;
-          cout << " <total_dipole_length> " << length(idipole + edipole_sum)
-               << " </total_dipole_length>" << endl;
+          cout << "<dipole>" << endl;
+          cout << setprecision(8);
+          cout << " <dipole_ion>   "
+               << setw(14) << idipole.x
+               << setw(14) << idipole.y
+               << setw(14) << idipole.z
+               << " </dipole_ion>" << endl;
+          cout << " <dipole_el>    "
+               << setw(14) << edipole_sum.x
+               << setw(14) << edipole_sum.y
+               << setw(14) << edipole_sum.z
+               << " </dipole_el>" << endl;
+          D3vector dipole_total = idipole + edipole_sum;
+          cout << " <dipole_total> "
+               << setw(14) << dipole_total.x
+               << setw(14) << dipole_total.y
+               << setw(14) << dipole_total.z
+               << " </dipole_total>" << endl;
+          cout << "</dipole>" << endl;
         }
         tmap["mlwf"].stop();
       }

src/ComputeMLWFCmd.cpp

@@ -35,6 +35,10 @@ int ComputeMLWFCmd::action(int argc, char **argv)
   if ( wf.nkp()>1 || length(wf.kpoint(0)) != 0.0 )
     throw runtime_error("ComputeMLWFCmd: can only be used with 1 k-point");
 
+  // Check if nempty > 0
+  if ( wf.nempty() > 0 )
+    throw runtime_error("ComputeMLWFCmd: cannot run with nempty > 0");
+
   if ( onpe0 )
     cout << "<mlwfs>" << endl;
 
@@ -84,13 +88,25 @@ int ComputeMLWFCmd::action(int argc, char **argv)
   if ( onpe0 )
   {
     D3vector idipole = s->atoms.dipole();
-    cout << setprecision(6);
-    cout << "   <ionic_dipole> " << idipole
-         << " </ionic_dipole>" << endl;
-    cout << "   <total_dipole> " << idipole + edipole_sum
-         << " </total_dipole>" << endl;
-    cout << "   <total_dipole_length> " << length(idipole + edipole_sum)
-         << " </total_dipole_length>" << endl;
+    cout << "<dipole>" << endl;
+    cout << setprecision(8);
+    cout << " <dipole_ion>   "
+         << setw(14) << idipole.x
+         << setw(14) << idipole.y
+         << setw(14) << idipole.z
+         << " </dipole_ion>" << endl;
+    cout << " <dipole_el>    "
+         << setw(14) << edipole_sum.x
+         << setw(14) << edipole_sum.y
+         << setw(14) << edipole_sum.z
+         << " </dipole_el>" << endl;
+    D3vector dipole_total = idipole + edipole_sum;
+    cout << " <dipole_total> "
+         << setw(14) << dipole_total.x
+         << setw(14) << dipole_total.y
+         << setw(14) << dipole_total.z
+         << " </dipole_total>" << endl;
+    cout << "</dipole>" << endl;
   }
   return 0;
 }

src/ElectricEnthalpy.cpp

@@ -645,58 +645,60 @@ void ElectricEnthalpy::print(ostream& os) const
   // print MLWF centers if pol_type_ == MLWF or MLWF_REF or MLWF_REF_Q
   if ( pol_type_ == mlwf || pol_type_ == mlwf_ref || pol_type_ == mlwf_ref_q )
   {
+    os << "<mlwfs>" << endl;
     int nst = sd_.nst();
     os << " <mlwfset spin=\"0\" size=\"" << nst << "\">" << endl;
     for ( int i = 0; i < nst; i++ )
     {
       if ( pol_type_ == mlwf )
       {
-        os << "   <mlwf> <center> " << setprecision(6)
-           << setw(12) << mlwfc_[i].x << " "
-           << setw(12) << mlwfc_[i].y << " "
-           << setw(12) << mlwfc_[i].z
+        os << " <mlwf> <center> " << setprecision(8)
+           << setw(14) << mlwfc_[i].x << " "
+           << setw(14) << mlwfc_[i].y << " "
+           << setw(14) << mlwfc_[i].z
            << " </center> </mlwf>"
            << endl;
       }
       else
       {
-        os << "   <mlwf> <center> " << setprecision(6)
-           << setw(12) << mlwfc_[i].x + correction_[i].x << " "
-           << setw(12) << mlwfc_[i].y + correction_[i].y << " "
-           << setw(12) << mlwfc_[i].z + correction_[i].z
+        os << " <mlwf> <center> " << setprecision(8)
+           << setw(14) << mlwfc_[i].x + correction_[i].x << " "
+           << setw(14) << mlwfc_[i].y + correction_[i].y << " "
+           << setw(14) << mlwfc_[i].z + correction_[i].z
            << " </center> </mlwf>"
            << endl;
 
         if ( compute_quadrupole_ )
-          os << "   <quad>"
-             << setw(12) << quad_[i][0] << " "
-             << setw(12) << quad_[i][4] << " "
-             << setw(12) << quad_[i][8] << " "
-             << setw(12) << quad_[i][1] << " "
-             << setw(12) << quad_[i][2] << " "
-             << setw(12) << quad_[i][5]
+          os << " <quad>"
+             << setw(14) << quad_[i][0] << " "
+             << setw(14) << quad_[i][4] << " "
+             << setw(14) << quad_[i][8] << " "
+             << setw(14) << quad_[i][1] << " "
+             << setw(14) << quad_[i][2] << " "
+             << setw(14) << quad_[i][5]
              << " </quad>" << endl;
       }
     }
     os << " </mlwfset>" << endl;
+    os << "</mlwfs>" << endl;
   }
 
   // print dipole
-  os << setprecision(10) << fixed << right;
-  os << " <dipole>\n";
-  os << "   <dipole_ion>   "
+  os << setprecision(8) << fixed << right;
+  os << "<dipole>\n";
+  os << " <dipole_ion>   "
      << setw(14) << dipole_ion_.x << " "
      << setw(14) << dipole_ion_.y << " "
      << setw(14) << dipole_ion_.z << " </dipole_ion>\n";
-  os << "   <dipole_el>    "
+  os << " <dipole_el>    "
      << setw(14) << dipole_el_.x << " "
      << setw(14) << dipole_el_.y << " "
      << setw(14) << dipole_el_.z << " </dipole_el>\n";
-  os << "   <dipole_total> "
+  os << " <dipole_total> "
      << setw(14) << dipole_total_.x << " "
      << setw(14) << dipole_total_.y << " "
      << setw(14) << dipole_total_.z << " </dipole_total>\n";
-  os << " </dipole>\n";
+  os << "</dipole>\n";
 
   if ( compute_quadrupole_ )
   {
@@ -721,30 +723,30 @@ void ElectricEnthalpy::print(ostream& os) const
     D3tensor q_traceless = q_total;
     q_traceless.traceless();
 
-    os << " <quadrupole> " << endl;
-    os << "   <quadrupole_ion> " << endl
+    os << "<quadrupole> " << endl;
+    os << " <quadrupole_ion> " << endl
        << q_ion
-       << "   </quadrupole_ion>" << endl;
-    os << "   <quadrupole_el> " << endl
+       << " </quadrupole_ion>" << endl;
+    os << " <quadrupole_el> " << endl
        << q_mlwf
-       << "   </quadrupole_el>" << endl;
-    os << "   <quadrupole_total> " << endl
+       << " </quadrupole_el>" << endl;
+    os << " <quadrupole_total> " << endl
        << q_total
-       << "   </quadrupole_total>" << endl;
-    os << "   <traceless_quadrupole> " << endl
+       << " </quadrupole_total>" << endl;
+    os << " <traceless_quadrupole> " << endl
        << q_traceless
-       << "   </traceless_quadrupole>" << endl;
+       << " </traceless_quadrupole>" << endl;
     char uplo = 'u';
     D3vector eigval;
     D3tensor eigvec;
     q_traceless.syev(uplo, eigval, eigvec);
-    os << "   <traceless_quadrupole_eigval> " << endl
-       << "    " << eigval << endl
-       << "   </traceless_quadrupole_eigval>" << endl;
-    os << "   <traceless_quadrupole_eigvec> " << endl
+    os << " <traceless_quadrupole_eigval> " << endl
+       << "  " << eigval << endl
+       << " </traceless_quadrupole_eigval>" << endl;
+    os << " <traceless_quadrupole_eigvec> " << endl
        << eigvec
-       << "   </traceless_quadrupole_eigvec>" << endl;
-    os << " </quadrupole> " << endl;
+       << " </traceless_quadrupole_eigvec>" << endl;
+    os << "</quadrupole> " << endl;
   }
 
 }

src/MLWFTransform.cpp

@@ -407,38 +407,26 @@ void MLWFTransform::apply_transform(SlaterDet& sd)
 ////////////////////////////////////////////////////////////////////////////////
 void MLWFTransform::print(ostream& os) const
 {
-  double total_spread[3] = { 0.0, 0.0, 0.0 };
   for ( int i = 0; i < sd_.nst(); i++ )
   {
     D3vector ctr = center(i);
-    double spi[3];
-    for (int j=0; j<3; j++)
-    {
-      spi[j] = spread2(i,j);
-      total_spread[j] += spi[j];
-    }
 
     os.setf(ios::fixed, ios::floatfield);
     os.setf(ios::right, ios::adjustfield);
     os << "   <mlwf>\n"
-       << "     <center>  " << setprecision(6)
-       << setw(12) << ctr.x
-       << setw(12) << ctr.y
-       << setw(12) << ctr.z
+       << "     <center>  " << setprecision(8)
+       << setw(14) << ctr.x
+       << setw(14) << ctr.y
+       << setw(14) << ctr.z
        << " </center>\n"
        << "     <spread2> "
-       << setw(12) << spi[0]
-       << setw(12) << spi[1]
-       << setw(12) << spi[2]
+       << setw(14) << spread2(i,0)
+       << setw(14) << spread2(i,1)
+       << setw(14) << spread2(i,2)
        << " </spread2>\n"
        << "   </mlwf>"
        << endl;
   }
-  D3vector edipole = dipole();
-  os << "   <e_dipole>  ";
-  for ( int j = 0; j < 3; j++ )
-    os << setprecision(6) << setw(12) << edipole[j];
-  os << " </e_dipole>" << endl;
 }
 ////////////////////////////////////////////////////////////////////////////////
 ostream& operator<<(ostream& os, MLWFTransform& mlwft)