Index: pcp/src/data.h
===================================================================
--- pcp/src/data.h	(revision 9bfb34dd562c18421160ddc443d8b87f8fda90c2)
+++ pcp/src/data.h	(revision 08e223d407a011d49395a8d992387a0c0cec0e60)
@@ -197,4 +197,5 @@
 #define WannierTag3 172   //!< Message consists of matrix elements A_ij sent/received during evaluating of operator matrices
 #define WannierTag4 173   //!< Message consists of matrix elements B_ij sent/received during evaluating of operator matrices
+#define WannierTag5 174   //!< Message consists of wannier centres and spread sent/received during ChangeWannierCentres()
 #define OtherPsiTag1 180  //!< Message consists of OnePsiElement data from the Spinype#SpinUp group
 #define OtherPsiTag2 181  //!< Message consists of OnePsiElement data from the Spinype#SpinDown group
Index: pcp/src/wannier.c
===================================================================
--- pcp/src/wannier.c	(revision 9bfb34dd562c18421160ddc443d8b87f8fda90c2)
+++ pcp/src/wannier.c	(revision 08e223d407a011d49395a8d992387a0c0cec0e60)
@@ -660,14 +660,16 @@
 /** Changes Wannier Centres according to RunStruct#CommonWannier.
  * \param *P Problem at hand.
- * \param Num number of Psis
- * \param **WannierCentre 2D array (NDIM, \a Num) with wannier centres
- * \param *WannierSpread array with wannier spread per wave function
- */
-void ChangeWannierCentres(struct Problem *P, int Num, double **WannierCentre, double *WannierSpread)
+ */
+void ChangeWannierCentres(struct Problem *P)
 {
   struct RunStruct *R = &P->R;
   struct Lattice *Lat = &P->Lat;
   struct LatticeLevel *LevS = R->LevS;
+  struct Psis *Psi = &Lat->Psi;
+  struct OnePsiElement *OnePsiA;
   int *marker, **group;
+  int Num = Psi->NoOfPsis;
+  double **WannierCentre = NULL;
+  double *WannierSpread = NULL;
   int partner[Num];
   int i,j,l,k;
@@ -676,5 +678,80 @@
   double Spread;
   int *N = LevS->Plan0.plan->N;
-  
+
+  int MPISignal = 0, mpisignal = 0;
+  unsigned int membersize, maxsize;
+  int msgsize;
+  double *buffer = NULL;
+  int position;
+  MPI_Status status;
+  
+  // determine necessary pack/buffer sizes for the wannier stuff
+  MPI_Pack_size(3, MPI_DOUBLE, P->Par.comm_ST_PsiT, &membersize);
+  maxsize = membersize;
+  MPI_Pack_size(1, MPI_DOUBLE, P->Par.comm_ST_PsiT, &membersize);
+  maxsize += membersize;
+  
+  // allocate memory
+  buffer = (double *) Malloc(sizeof(double)*maxsize, "ChangeWannierCentres: *buffer");
+  WannierSpread = (double *) Malloc(sizeof(double)*Num, "ChangeWannierCentres: *WannierSpread");
+  WannierCentre = (double **) Malloc(sizeof(double *)*Num, "ChangeWannierCentres: **WannierCentre");
+  for (j=0;j<Num;j++)
+    WannierCentre[j] = (double *) Malloc(sizeof(double)*(NDIM+1), "ChangeWannierCentres: *WannierCentre[]");
+
+  // exchange wannier centers and spread among processes by going over all wave functions
+  i=-1;
+  for (l=0; l < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; l++) {
+    OnePsiA = &Psi->AllPsiStatus[l];    // grab OnePsiA
+    if (OnePsiA->PsiType == type) {   // drop all but occupied ones
+      i++;  // increase l if it is occupied wave function
+      
+      // send and receive centres and spread (we use the MPI_Pack'ing scheme to preserve correct order of data)
+      if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) {// is this local?
+        position = 0;
+        MPI_Pack(Psi->AddData[OnePsiA->MyLocalNo].WannierCentre, 3, MPI_DOUBLE, buffer, maxsize, &position, P->Par.comm_ST_PsiT);
+        MPI_Pack(&Psi->AddData[OnePsiA->MyLocalNo].WannierSpread, 1, MPI_DOUBLE, buffer, maxsize, &position, P->Par.comm_ST_PsiT);
+        for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++) 
+          if (k != P->Par.my_color_comm_ST_Psi) {  // send to all other processes
+            fprintf(stderr, "(%i) Send goes from %d to %d\n", P->Par.me, k, P->Par.my_color_comm_ST_Psi);
+            if (MPI_Send(buffer, position, MPI_PACKED, k, WannierTag5, P->Par.comm_ST_PsiT) != MPI_SUCCESS) {
+              mpisignal = 1;
+            }
+          } else {  // "send" to us
+            for(j=0;j<NDIM;j++)  // copy to own array as well
+              WannierCentre[i][j] = Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j];
+            WannierSpread[i] = Psi->AddData[OnePsiA->MyLocalNo].WannierSpread;
+          }
+      } else {
+        fprintf(stderr, "(%i) Receive is from %d to %d\n", P->Par.me, OnePsiA->my_color_comm_ST_Psi, P->Par.my_color_comm_ST_Psi);
+        if (MPI_Recv(buffer, maxsize, MPI_PACKED, OnePsiA->my_color_comm_ST_Psi, WannierTag5, P->Par.comm_ST_PsiT, &status) != MPI_SUCCESS)
+          mpisignal = 1;
+        position = 0;
+        MPI_Get_count(&status, MPI_PACKED, &msgsize);
+        MPI_Unpack (buffer, msgsize, &position, WannierCentre[i], 3, MPI_DOUBLE, P->Par.comm_ST_PsiT);
+        MPI_Unpack (buffer, msgsize, &position, &WannierSpread[i], 1, MPI_DOUBLE, P->Par.comm_ST_PsiT);
+      }
+      
+      // check for MPI errors
+      if (MPI_Allreduce(&mpisignal, &MPISignal, 1, MPI_INT, MPI_SUM, P->Par.comm_ST_PsiT) != MPI_SUCCESS) // allreduce signals
+        Error(SomeError,"ChangeWannierCentres: Bcast of signal failed\n");
+      if (MPISignal) { // greater 0 means there was at least one hickup somewhere
+        Error(SomeError,"ChangeWannierCentres: Scattering of Wanniers failed\n");
+      }
+    }
+  }
+  // free the buffer memory
+  Free(buffer, "ChangeWannierCentres: *buffer");
+
+  if ((P->Call.out[NormalOut]) && (1)) { // root prints gathered wannier centres for debugging
+    fprintf(stderr, "(%i) Printing gathered Wannier centres:\n", P->Par.me);
+    for(l=0;l<Num;l++) {
+      fprintf(stderr, "(%i) ", P->Par.me);
+      for(j=0;j<NDIM;j++)
+        fprintf(stderr, "%lg\t", WannierCentre[l][j]);
+      fprintf(stderr, "-\t%lg\n", WannierSpread[l]);
+    }
+  }
+      
+  // change centres
   switch (R->CommonWannier) {
    case 4:
@@ -728,7 +805,7 @@
           Spread += (WannierCentre[l][i] - WannierCentre[k][i])*(WannierCentre[l][i] - WannierCentre[k][i]);
         }
-        Spread = sqrt(Spread);  // distance in Spread
+        //Spread = sqrt(Spread);  // distance in Spread
         //fprintf(stderr,"(%i) %i to %i: distance %e, SpreadSum = %e + %e = %e \n", P->Par.me, l, k, Spread, WannierSpread[l], WannierSpread[k], WannierSpread[l]+WannierSpread[k]); 
-        if (Spread < 1.5*(WannierSpread[l]+WannierSpread[k])) {// if distance smaller than sum of spread
+        if (Spread < 5*(WannierSpread[l]+WannierSpread[k])*(WannierSpread[l]+WannierSpread[k])) {// if distance smaller than sum of spread times two plus safety
           group[l][j++] = k;  // add k to group of l
           partner[l]++;
@@ -821,4 +898,46 @@
    break;
   }
+
+  if ((P->Call.out[NormalOut]) && (1)) { // root prints gathered wannier centres for debugging
+    fprintf(stderr, "(%i) Printing changed Wannier centres:\n", P->Par.me);
+    for(l=0;l<Num;l++) {
+      fprintf(stderr, "(%i) ", P->Par.me);
+      for(j=0;j<NDIM;j++)
+        fprintf(stderr, "%lg\t", WannierCentre[l][j]);
+      fprintf(stderr, "-\t%lg\n", WannierSpread[l]);
+    }
+  }
+      
+  // store info into local arrays again
+  i=-1;
+  for (l=0; l < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; l++) {  // go through all wave functions
+    OnePsiA = &Psi->AllPsiStatus[l];    // grab OnePsiA
+    if (OnePsiA->PsiType == type) {   // drop all but occupied ones
+      i++;  // increase l if it is occupied wave function
+      if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) {// is this local
+        for(j=0;j<NDIM;j++)
+          Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j] = WannierCentre[i][j];
+      }
+    }
+  }
+
+  // show locally stored info for debugging
+  if (P->Call.out[NormalOut]) {
+    i=-1;
+    for (l=0; l < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; l++) {  // go through all wave functions
+      OnePsiA = &Psi->AllPsiStatus[l];    // grab OnePsiA
+      if (OnePsiA->PsiType == type) {   // drop all but occupied ones
+        i++;  // increase l if it is occupied wave function
+        if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // is this local?
+           fprintf(stderr,"(%i) Psi %i, L %i: (x,y,z) = (%lg, %lg, %lg), Spread %lg\n",P->Par.me,i, R->LevSNo, Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[0], Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[1], Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[2], Psi->AddData[OnePsiA->MyLocalNo].WannierSpread);
+      }
+    }
+  }
+
+  // free memory
+  Free(WannierSpread, "ChangeWannierCentres: *WannierSpread");
+  for (j=0;j<Num;j++)
+    Free(WannierCentre[j], "ChangeWannierCentres: *WannierCentre[]");
+  Free(WannierCentre, "ChangeWannierCentres: **WannierCentre");
 }
 
@@ -1175,8 +1294,4 @@
   ComputeWannierCentresfromVarianceMatrices(P, DiagData.A, &spread, &spreadSQ, WannierCentre, WannierSpread);
 
-  // join Wannier orbital to groups with common centres under certain conditions
-  //debug (P,"Changing Wannier Centres according to CommonWannier");
-  ChangeWannierCentres(P, Num, WannierCentre, WannierSpread); 
-  
   // write to file
   //debug (P,"Writing spread file");
@@ -1830,4 +1945,5 @@
         for (j=0;j<NDIM;j++)
           Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j] = WannierCentre[i][j];
+        Psi->AddData[OnePsiA->MyLocalNo].WannierSpread = WannierSpread[i];
       }
       if (P->Par.me_comm_ST == 0) 
@@ -1934,5 +2050,5 @@
           //return 0;
         // and store 'em (for all who have this Psi local)
-        if ((P->Par.me == 0) && P->Call.out[ValueOut]) fprintf(stderr,"(%i) Psi %i, L %i: (x,y,z) = (%lg, %lg, %lg), Spread %lg\n",P->Par.me,i, R->LevSNo, WannierCentre[0], WannierCentre[1], WannierCentre[2], WannierCentre[NDIM]);
+        if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) Psi %i, L %i: (x,y,z) = (%lg, %lg, %lg), Spread %lg\n",P->Par.me,i, R->LevSNo, WannierCentre[0], WannierCentre[1], WannierCentre[2], WannierCentre[NDIM]);
         for (j=0;j<NDIM;j++) Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j] = WannierCentre[j];
         Psi->AddData[OnePsiA->MyLocalNo].WannierSpread = WannierCentre[NDIM];
Index: pcp/src/wannier.h
===================================================================
--- pcp/src/wannier.h	(revision 9bfb34dd562c18421160ddc443d8b87f8fda90c2)
+++ pcp/src/wannier.h	(revision 08e223d407a011d49395a8d992387a0c0cec0e60)
@@ -34,8 +34,10 @@
 };
 
+
 void PrintGSLMatrix(struct Problem *P, gsl_matrix *U, int Num, const char *msg);
 void ComputeMLWF(struct Problem *P);
 void WriteWannierFile(struct Problem *P, double spread, double old_spread, double **WannierCentre, double *WannierSpread);
 int ParseWannierFile(struct Problem *P);
+void ChangeWannierCentres(struct Problem *P);
 void SerialDiagonalization(struct Problem *P, struct DiagonalizationData *DiagData);
 void ParallelDiagonalization(struct Problem *P, struct DiagonalizationData *DiagData);
