Merge branch 'develop'

src/AndersonMixer.cpp

@@ -96,13 +96,17 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
       const int kmi = ( k_ - i + nmax_ ) % ( nmax_ + 1 );
       assert(kmi>=0);
       assert(kmi<nmax_+1);
-      // cout << "k=" << k_ << " i=" << i << " kmi=" << kmi << endl;
+#ifdef DEBUG
+      cout << "k=" << k_ << " i=" << i << " kmi=" << kmi << endl;
+#endif
       for ( int j = 0; j <= i; j++ )
       {
         const int kmj = ( k_ - j + nmax_ ) % ( nmax_ + 1 );
         assert(kmj>=0);
         assert(kmj<nmax_+1);
-        // cout << "k=" << k_ << " j=" << j << " kmj=" << kmj << endl;
+#ifdef DEBUG
+        cout << "k=" << k_ << " j=" << j << " kmj=" << kmj << endl;
+#endif
         double sum = 0.0;
         for ( int l = 0; l < m_; l++ )
           sum += (f_[k_][l] - f_[kmi][l]) * (f_[k_][l] - f_[kmj][l]);
@@ -141,10 +145,12 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
         int info;
         dsyev(&jobz,&uplo,&n_,&a[0],&n_,&w[0],&work[0],&lwork,&info);
 
+#ifdef DEBUG
         cout << "AndersonMixer: eigenvalues: ";
         for ( int i = 0; i < n_; i++ )
           cout << w[i] << "  ";
         cout << endl;
+#endif
         if ( info != 0 )
         {
           cerr << " AndersonMixer: Error in dsyev" << endl;
@@ -229,21 +235,25 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
           }
           norm_ok = theta_sum <= 1.0;
 #endif
-          // cout << " tp = " << tikhonov_parameter
-          //      << " AndersonMixer: theta = ";
-          // for ( int i = 0; i < theta.size(); i++ )
-          //   cout << theta[i] << " ";
-          // cout << endl;
+#ifdef DEBUG
+          cout << " tp = " << tikhonov_parameter
+               << " AndersonMixer: theta = ";
+          for ( int i = 0; i < theta.size(); i++ )
+            cout << theta[i] << " ";
+          cout << endl;
+#endif
 
           tikhonov_parameter *= 2.0;
           iter++;
         }
       }
 
+#ifdef DEBUG
       cout << " AndersonMixer: theta = ";
       for ( int i = 0; i < theta.size(); i++ )
         cout << theta[i] << " ";
       cout << endl;
+#endif
 
     }
 

src/AtomSet.cpp

@@ -19,6 +19,7 @@
 #include "AtomSet.h"
 #include "Species.h"
 #include "NameOf.h"
+#include "MPIdata.h"
 #include "sampling.h"
 #include <iostream>
 #include <algorithm>
@@ -47,7 +48,7 @@ bool AtomSet::addSpecies(Species* sp, string name)
   if ( s != 0 )
   {
     // species is already defined: substitute with new definition
-    if ( ctxt_.onpe0() )
+    if ( MPIdata::onpe0() )
     {
       cout << " AtomSet::addSpecies: species " << name
            << " is already defined" << endl;
@@ -75,7 +76,7 @@ bool AtomSet::addSpecies(Species* sp, string name)
     atom_list.resize(atom_list.size()+1);
   }
 
-  if ( ctxt_.onpe0() )
+  if ( MPIdata::onpe0() )
   {
     cout << endl << " species " << sp->name() << ":" << endl;
     sp->info(cout);
@@ -91,7 +92,7 @@ bool AtomSet::addAtom(Atom *a)
   if ( findAtom(a->name()) )
   {
     // this name is already in the atom list, reject atom definition
-    if ( ctxt_.onpe0() )
+    if ( MPIdata::onpe0() )
       cout << " AtomSet:addAtom: atom " << a->name()
            << " is already defined" << endl;
     return false;
@@ -103,7 +104,7 @@ bool AtomSet::addAtom(Atom *a)
   if ( !s )
   {
     // species not found, cannot define atom
-    if ( ctxt_.onpe0() )
+    if ( MPIdata::onpe0() )
       cout << " AtomSet:addAtom: species " << spname
            << " is undefined" << endl;
     return false;
@@ -164,7 +165,7 @@ bool AtomSet::delAtom(string name)
   }
 
   // this name was not found in the atom list
-  if ( ctxt_.onpe0() )
+  if ( MPIdata::onpe0() )
     cout << " AtomSet:delAtom: no such atom: " << name << endl;
   return false;
 }
@@ -230,7 +231,7 @@ void AtomSet::listAtoms(void) const
   {
     for ( int ia = 0; ia < atom_list[is].size(); ia++ )
     {
-      if ( ctxt_.onpe0() )
+      if ( MPIdata::onpe0() )
         cout << *atom_list[is][ia];
     }
   }
@@ -241,7 +242,7 @@ void AtomSet::listSpecies(void) const
 {
   for ( int is = 0; is < species_list.size(); is++ )
   {
-    if ( ctxt_.onpe0() )
+    if ( MPIdata::onpe0() )
     {
       cout << endl << " species " << spname[is] << ":" << endl;
       species_list[is]->info(cout);
@@ -345,14 +346,7 @@ void AtomSet::sync_positions(vector<vector<double> >& tau)
   {
     int m = tau[is].size();
     double* p = &tau[is][0];
-    if ( ctxt_.onpe0() )
-    {
-      ctxt_.dbcast_send(m,1,p,m);
-    }
-    else
-    {
-      ctxt_.dbcast_recv(m,1,p,m,0,0);
-    }
+    MPI_Bcast(p,m,MPI_DOUBLE,0,MPIdata::comm());
   }
 }
 
@@ -400,14 +394,7 @@ void AtomSet::sync_velocities(vector<vector<double> >& vel)
   {
     int m = vel[is].size();
     double* p = &vel[is][0];
-    if ( ctxt_.onpe0() )
-    {
-      ctxt_.dbcast_send(m,1,p,m);
-    }
-    else
-    {
-      ctxt_.dbcast_recv(m,1,p,m,0,0);
-    }
+    MPI_Bcast(p,m,MPI_DOUBLE,0,MPIdata::comm());
   }
 }
 
@@ -764,28 +751,21 @@ void AtomSet::sync_cell(void)
 ////////////////////////////////////////////////////////////////////////////////
 void AtomSet::sync_cell(UnitCell& cell)
 {
-  if ( ctxt_.onpe0() )
-  {
-    double sbuf[9];
-    for ( int i = 0; i < 9; i++ )
-      sbuf[i] = cell.amat(i);
-    ctxt_.dbcast_send(9,1,sbuf,9);
-  }
-  else
-  {
-    double rbuf[9];
-    ctxt_.dbcast_recv(9,1,rbuf,9,0,0);
-    D3vector a0(rbuf[0],rbuf[1],rbuf[2]);
-    D3vector a1(rbuf[3],rbuf[4],rbuf[5]);
-    D3vector a2(rbuf[6],rbuf[7],rbuf[8]);
-    cell.set(a0,a1,a2);
-  }
+  double buf[9];
+  for ( int i = 0; i < 9; i++ )
+    buf[i] = cell.amat(i);
+  MPI_Bcast(buf,9,MPI_DOUBLE,0,MPIdata::comm());
+
+  D3vector a0(buf[0],buf[1],buf[2]);
+  D3vector a1(buf[3],buf[4],buf[5]);
+  D3vector a2(buf[6],buf[7],buf[8]);
+  cell.set(a0,a1,a2);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 ostream& operator << ( ostream &os, const AtomSet &as )
 {
-  if ( as.context().onpe0() )
+  if ( MPIdata::onpe0() )
   {
     os << "<atomset>\n";
     os << as.cell();

src/AtomSet.h

@@ -19,7 +19,6 @@
 #ifndef ATOMSET_H
 #define ATOMSET_H
 
-#include "Context.h"
 #include "Atom.h"
 #include "UnitCell.h"
 #include "D3tensor.h"
@@ -34,8 +33,6 @@ class AtomSet
 {
   private:
 
-  const Context& ctxt_;
-
   int nel_;
   std::map<std::string,int> na_;  // na_[sp_name]: num. of at. of spec. sp_name
   std::map<std::string,int> isp_; // isp_[sp_name]: index of species sp_name
@@ -47,13 +44,12 @@ class AtomSet
 
   public:
 
-  AtomSet(const Context& ctxt) : ctxt_(ctxt), nel_(0) {}
+  AtomSet(void) : nel_(0) {}
   ~AtomSet(void);
 
   std::vector<std::vector<Atom *> > atom_list; // atom_list[is][ia]
   std::vector<Species *> species_list;    // species_list[is]
 
-  const Context& context(void) const { return ctxt_; }
   bool addAtom(Atom *a);
   bool delAtom(std::string name);
   bool addSpecies(Species *sp, std::string name);

src/BOSampleStepper.h

@@ -33,7 +33,6 @@ class BOSampleStepper : public SampleStepper
   private:
 
   Wavefunction dwf;
-  Wavefunction* wfv;
   int nitscf_;
   int nite_;
   ChargeDensity cd_;

src/Basis.h

@@ -22,12 +22,7 @@
 #include "D3vector.h"
 #include "UnitCell.h"
 #include <vector>
-
-#ifdef USE_MPI
 #include <mpi.h>
-#else
-  typedef int MPI_Comm;
-#endif
 
 class Basis
 {

src/Bisection.cpp

@@ -97,6 +97,8 @@ Bisection::Bisection(const SlaterDet& sd, const int nlevels[3])
   // largest number of levels
   nlevelsmax_=max(nlevels[0],max(nlevels[1],nlevels[2]));
 
+  assert( nlevelsmax_ >0 );
+
   // real-space grid size for wave functions
   np_[0] = sd.basis().np(0);
   np_[1] = sd.basis().np(1);
@@ -240,8 +242,6 @@ void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)
   for ( int i = 0; i < nmat_; i++ )
     amat_[i]->clear();
 
-  int size=amat_[0]->size();
-
   // allocate matrix for products of projected parts
   DoubleMatrix products(c.context(),c.n(),c.n(),c.nb(),c.nb());
 
@@ -268,6 +268,7 @@ void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)
             // add product to the matrix
             double *coeff_source=products.valptr(0);
             double *coeff_destination=amat_[imat]->valptr(0);
+            const int size=amat_[imat]->size();
             for ( int i=0; i<size; i++ )
               coeff_destination[i]+=coeff_source[i];
           }
@@ -293,6 +294,7 @@ void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)
         // normalize coeffs
         double *coeff=amat_[imat]->valptr(0);
         double fac = 1.0 / (np_[0]*np_[1]*np_[2]);
+        const int size=amat_[imat]->size();
         for ( int i = 0; i < size; i++ )
           coeff[i] *= fac;
         imat++;
@@ -432,9 +434,13 @@ void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)
     double time = (*i).second.real();
     double tmin = time;
     double tmax = time;
-    ctxt_.dmin(1,1,&tmin,1);
-    ctxt_.dmax(1,1,&tmax,1);
-    if ( ctxt_.myproc()==0 )
+    double sbuf = tmin;
+    double rbuf = 0.0;
+    MPI_Reduce(&sbuf,&rbuf,1,MPI_DOUBLE,MPI_MIN,0,MPIdata::comm());
+    sbuf = tmax;
+    rbuf = 0.0;
+    MPI_Reduce(&sbuf,&rbuf,1,MPI_DOUBLE,MPI_MAX,0,MPIdata::comm());
+    if ( MPIdata::onpe0() )
     {
       string s = "name=\"" + (*i).first + "\"";
       cout << "<timing " << left << setw(22) << s

src/Bisection.h

@@ -33,7 +33,7 @@ class Bisection
 {
   private:
 
-    Context ctxt_;
+    const Context& ctxt_;
 
     // bisection levels in each directions
     int nlevels_[3]; // bisection level

src/BisectionCmd.h

@@ -94,51 +94,64 @@ class BisectionCmd : public Cmd
       return 1;
     }
 
-    tm.reset();
-    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
+    if ( ( nLevels[0]+nLevels[1]+nLevels[2] ) == 0)
     {
-      SlaterDet& sd = *wf.sd(0,0);
-      Bisection bisection(sd,nLevels);
-      const int maxsweep = 20;
-      const double tol = 1.e-8;
-      tm.start();
-      bisection.compute_transform(sd,maxsweep,tol);
-      bisection.compute_localization(epsilon);
-      bisection.forward(sd);
-      tm.stop();
-      vector<long int> localization = bisection.localization();
-
       if ( ui->onpe0() )
       {
-        cout << " BisectionCmd: lx=" << nLevels[0]
-             << " ly=" << nLevels[1]
-             << " lz=" << nLevels[2]
-             << " threshold=" << epsilon << endl;
-
-        // print localization vectors and overlaps
-        int sum = 0;
-        const int nst = localization.size();
-        for ( int i = 0; i < nst; i++ )
+        cout << " BisectionCmd: at least one level must be positive" << endl;
+      }
+      return 1;
+    }
+
+    tm.reset();
+    for ( int isp_loc = 0; isp_loc < wf.nsp_loc(); ++isp_loc )
+    {
+      const int ikp_loc = wf.ikp_local(0);
+      if ( ikp_loc >= 0 )
+        {
+        SlaterDet& sd = *wf.sd(isp_loc,ikp_loc);
+        Bisection bisection(sd,nLevels);
+        const int maxsweep = 20;
+        const double tol = 1.e-8;
+        tm.start();
+        bisection.compute_transform(sd,maxsweep,tol);
+        bisection.compute_localization(epsilon);
+        bisection.forward(sd);
+        tm.stop();
+        vector<long int> localization = bisection.localization();
+
+        if ( ui->onpe0() )
         {
-          int count = 0;
-          for ( int j = 0; j < nst; j++ )
+          cout << " BisectionCmd: lx=" << nLevels[0]
+               << " ly=" << nLevels[1]
+               << " lz=" << nLevels[2]
+               << " threshold=" << epsilon << endl;
+
+          // print localization vectors and overlaps
+          int sum = 0;
+          const int nst = localization.size();
+          for ( int i = 0; i < nst; i++ )
           {
-            if ( bisection.overlap(i,j) )
-              count++;
+            int count = 0;
+            for ( int j = 0; j < nst; j++ )
+            {
+              if ( bisection.overlap(i,j) )
+                count++;
+            }
+            cout << "localization[" << i << "]: "
+                 << localization[i] << " "
+                 << bitset<30>(localization[i])
+                 << " size: " << bisection.size(i)
+                 << " overlaps: " << count << endl;
+            sum += count;
           }
-          cout << "localization[" << i << "]: "
-               << localization[i] << " "
-               << bitset<30>(localization[i])
-               << " size: " << bisection.size(i)
-               << " overlaps: " << count << endl;
-          sum += count;
+          cout << " Bisection: total size:    ispin=" << wf.isp_global(isp_loc)
+               << ": " << bisection.total_size() << endl;
+          cout << " Bisection: pair fraction: ispin=" << wf.isp_global(isp_loc)
+               << ": " << bisection.pair_fraction() << endl;
         }
-        cout << " Bisection: total size:    ispin=" << ispin
-             << ": " << bisection.total_size() << endl;
-        cout << " Bisection: pair fraction: ispin=" << ispin
-             << ": " << bisection.pair_fraction() << endl;
-      }
-    }
+      } // if ikp_loc
+    } // isp_loc
     if ( ui->onpe0() )
       cout << " Bisection: time=" << tm.real() << endl;
     return 0;

src/CGCellStepper.cpp

@@ -16,6 +16,7 @@
 //
 ////////////////////////////////////////////////////////////////////////////////
 
+#include "MPIdata.h"
 #include "CGCellStepper.h"
 #include "CGOptimizer.h"
 using namespace std;
@@ -29,7 +30,7 @@ CGCellStepper::CGCellStepper(Sample& s) : CellStepper(s),
   cgopt_.set_alpha_max(0.5);
   cgopt_.set_beta_max(10.0);
 #ifdef DEBUG
-  if ( s.ctxt_.onpe0() )
+  if ( MPIdata::onpe0() )
     cgopt_.set_debug_print();
 #endif
 
@@ -127,7 +128,7 @@ void CGCellStepper::compute_new_cell(double e, const valarray<double>& sigma,
 
   cgopt_.compute_xp(x,e,g,xp);
 
-  if ( s_.ctxt_.onpe0() )
+  if ( MPIdata::onpe0() )
   {
     cout << "  CGCellStepper: alpha = " << cgopt_.alpha() << endl;
   }

src/CGIonicStepper.cpp

@@ -16,6 +16,7 @@
 //
 ////////////////////////////////////////////////////////////////////////////////
 
+#include "MPIdata.h"
 #include "CGIonicStepper.h"
 #include "CGOptimizer.h"
 using namespace std;
@@ -28,7 +29,7 @@ CGIonicStepper::CGIonicStepper(Sample& s) : IonicStepper(s),
   cgopt_.set_alpha_max(50.0);
   cgopt_.set_beta_max(10.0);
 #ifdef DEBUG
-  if ( s.ctxt_.onpe0() )
+  if ( MPIdata::onpe0() )
     cgopt_.set_debug_print();
 #endif
 }
@@ -69,7 +70,7 @@ void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
     largest_disp = max(largest_disp,fabs(xp[i]-x[i]));
   if ( largest_disp > max_disp )
   {
-    if ( s_.ctxt_.onpe0() )
+    if ( MPIdata::onpe0() )
       cout << "  CGIonicStepper: displacement exceeds limit, rescaling" << endl;
     // rescale displacement and reset the CG optimizer
     double fac = max_disp/largest_disp;
@@ -78,7 +79,7 @@ void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
     cgopt_.reset();
   }
 
-  if ( s_.ctxt_.onpe0() )
+  if ( MPIdata::onpe0() )
   {
     cout << "  CGIonicStepper: alpha = " << cgopt_.alpha() << endl;
   }

src/CPSampleStepper.cpp

@@ -16,6 +16,7 @@
 //
 ////////////////////////////////////////////////////////////////////////////////
 
+#include "MPIdata.h"
 #include "CPSampleStepper.h"
 #include "SlaterDet.h"
 #include "MDWavefunctionStepper.h"
@@ -60,16 +61,16 @@ CPSampleStepper::~CPSampleStepper(void)
 {
   delete mdwf_stepper;
   if ( mdionic_stepper != 0 ) delete mdionic_stepper;
+
   for ( TimerMap::iterator i = tmap.begin(); i != tmap.end(); i++ )
   {
-    double time = (*i).second.real();
-    double tmin = time;
-    double tmax = time;
-    s_.ctxt_.dmin(1,1,&tmin,1);
-    s_.ctxt_.dmax(1,1,&tmax,1);
-    if ( s_.ctxt_.myproc()==0 )
+    double time = i->second.real();
+    double tmin, tmax;
+    MPI_Reduce(&time,&tmin,1,MPI_DOUBLE,MPI_MIN,0,MPIdata::comm());
+    MPI_Reduce(&time,&tmax,1,MPI_DOUBLE,MPI_MAX,0,MPIdata::comm());
+    if ( MPIdata::onpe0() && (tmax > 0.0) )
     {
-      string s = "name=\"" + (*i).first + "\"";
+      string s = "name=\"" + i->first + "\"";
       cout << "<timing " << left << setw(22) << s
            << " min=\"" << setprecision(3) << tmin << "\""
            << " max=\"" << setprecision(3) << tmax << "\"/>"
@@ -81,13 +82,13 @@ CPSampleStepper::~CPSampleStepper(void)
 ////////////////////////////////////////////////////////////////////////////////
 void CPSampleStepper::step(int niter)
 {
-  const bool onpe0 = s_.ctxt_.onpe0();
+  const bool onpe0 = MPIdata::onpe0();
 
   // check that there are no fractionally occupied states
   // next line: (3-nspin) = 2 if nspin==1 and 1 if nspin==2
   if ( s_.wf.nel() != (( 3 - s_.wf.nspin() ) * s_.wf.nst()) )
   {
-    if ( s_.ctxt_.onpe0() )
+    if ( onpe0 )
     {
       cout << " CPSampleStepper::step:"
               " fractionally occupied or empty states: cannot run" << endl;
@@ -141,7 +142,7 @@ void CPSampleStepper::step(int niter)
   {
     tm_iter.reset();
     tm_iter.start();
-    if ( s_.ctxt_.mype() == 0 )
+    if ( onpe0 )
       cout << "<iteration count=\"" << iter+1 << "\">\n";
 
     mdwf_stepper->update(dwf);
@@ -255,17 +256,13 @@ void CPSampleStepper::step(int niter)
       ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
     ef_.enthalpy();
 
-    if ( s_.ctxt_.mype() == 0 )
-      cout << "</iteration>" << endl;
-
-    // print iteration time
     tm_iter.stop();
+    // print iteration time
     double time = tm_iter.real();
-    double tmin = time;
-    double tmax = time;
-    s_.ctxt_.dmin(1,1,&tmin,1);
-    s_.ctxt_.dmax(1,1,&tmax,1);
-    if ( s_.ctxt_.myproc()==0 )
+    double tmin, tmax;
+    MPI_Reduce(&time,&tmin,1,MPI_DOUBLE,MPI_MIN,0,MPIdata::comm());
+    MPI_Reduce(&time,&tmax,1,MPI_DOUBLE,MPI_MAX,0,MPIdata::comm());
+    if ( onpe0 )
     {
       string s = "name=\"iteration\"";
       cout << "<timing " << left << setw(22) << s
@@ -273,6 +270,10 @@ void CPSampleStepper::step(int niter)
            << " max=\"" << setprecision(3) << tmax << "\"/>"
            << endl;
     }
+
+    if ( onpe0 )
+      cout << "</iteration>" << endl;
+
     if ( compute_forces )
       s_.constraints.update_constraints(dt);
   } // iter

src/ChargeDensity.h

@@ -38,8 +38,6 @@ class ChargeDensity
 {
   private:
 
-  const Context& ctxt_;
-  MPI_Comm vcomm_;
   const Wavefunction& wf_;
   Basis* vbasis_;
   FourierTransform* vft_;
@@ -63,8 +61,6 @@ class ChargeDensity
   void update_taur(double* taur_up, double* taur_dn) const;
   void update_rhog(void);
 
-  const Context& context(void) const { return ctxt_; }
-  MPI_Comm vcomm(void) const { return vcomm_; }
   Basis* vbasis(void) const { return vbasis_; }
   FourierTransform* vft(void) const { return vft_; }
   FourierTransform* ft(int ikp) const { return ft_[ikp]; }

src/ComputeMLWFCmd.cpp

@@ -20,6 +20,7 @@
 #include<iostream>
 #include "Context.h"
 #include "SlaterDet.h"
+#include "cout0.h"
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -27,10 +28,12 @@ int ComputeMLWFCmd::action(int argc, char **argv)
 {
   Wavefunction& wf = s->wf;
 
+  const bool onpe0 = MPIdata::onpe0();
+
   // Check that only the k=0 point is used
   if ( wf.nkp()>1 || length(wf.kpoint(0)) != 0.0 )
   {
-    if ( ui->onpe0() )
+    if ( onpe0 )
     {
       cout << " ComputeMLWFCmd::action: compute_mlwf can only be used at\n"
            << " the Gamma point (k=0)" << endl;
@@ -38,46 +41,75 @@ int ComputeMLWFCmd::action(int argc, char **argv)
     return 1;
   }
 
-  if ( ui->onpe0() )
+  if ( onpe0 )
     cout << "<mlwfs>" << endl;
 
   D3vector edipole_sum;
-  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
+  for ( int ispin = 0; ispin < wf.nspin(); ++ispin )
   {
-    SlaterDet& sd = *(wf.sd(ispin,0));
-    MLWFTransform* mlwft = new MLWFTransform(sd);
+    const int isp_loc = wf.isp_local(ispin);
+    const int ikp = 0;
+    const int ikp_loc = wf.ikp_local(ikp);
+    ostringstream ostr;
+    int isrc = -1;
+    if ( ( isp_loc >= 0 ) && ( ikp_loc >= 0 ) )
+    {
+      SlaterDet& sd = *(wf.sd(isp_loc,ikp_loc));
+      MLWFTransform* mlwft = new MLWFTransform(sd);
 
-    mlwft->update();
-    mlwft->compute_transform();
-    mlwft->apply_transform(sd);
+      mlwft->update();
+      mlwft->compute_transform();
+      mlwft->apply_transform(sd);
 
-    if ( ui->onpe0() )
-    {
-      cout << " <mlwfset spin=\"" << ispin
-           << "\" size=\"" << sd.nst() << "\">" << endl;
-      for ( int i = 0; i < sd.nst(); i++ )
+      if ( MPIdata::sd_rank() == 0 )
       {
-        D3vector ctr = mlwft->center(i);
-        double sp = mlwft->spread(i);
-        cout.setf(ios::fixed, ios::floatfield);
-        cout.setf(ios::right, ios::adjustfield);
-        cout << "   <mlwf center=\"" << setprecision(6)
-             << setw(12) << ctr.x
-             << setw(12) << ctr.y
-             << setw(12) << ctr.z
-             << " \" spread=\" " << sp << " \"/>"
-             << endl;
-      }
-      D3vector edipole = mlwft->dipole();
-      cout << " <electronic_dipole spin=\"" << ispin << "\"> " << edipole
-           << " </electronic_dipole>" << endl;
-      cout << " </mlwfset>" << endl;
-      edipole_sum += edipole;
+        ostr.str("");
+        isrc = MPIdata::rank();
+        ostr << " <mlwfset spin=\"" << ispin
+             << "\" size=\"" << sd.nst() << "\">" << endl;
+        double total_spread[6];
+        for ( int j = 0; j < 6; j++ )
+           total_spread[j] = 0.0;
+        for ( int i = 0; i < sd.nst(); i++ )
+        {
+          D3vector ctr = mlwft->center(i);
+          double spi[6];
+          for (int j=0; j<3; j++)
+          {
+            spi[j] = mlwft->spread2(i,j);
+            total_spread[j] += spi[j];
+          }
+
+          ostr.setf(ios::fixed, ios::floatfield);
+          ostr.setf(ios::right, ios::adjustfield);
+          ostr << "   <mlwf center=\"" << setprecision(6)
+               << setw(12) << ctr.x
+               << setw(12) << ctr.y
+               << setw(12) << ctr.z
+               << " \" spread=\" "
+               << setw(12) << spi[0]
+               << setw(12) << spi[1]
+               << setw(12) << spi[2] << " \"/>"
+               << endl;