Changeset 08e223


Ignore:
Timestamp:
May 30, 2008, 9:04:37 AM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
6df7db
Parents:
9bfb34
Message:

ChangeWannierCentres can now be called from outside ComputeMLWF()

Wannier centres and spread are exchanged via MPI to allow for external calling (necessary if read src densities shall be perturbed with different CommonWannier values!)

Location:
pcp/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • pcp/src/data.h

    r9bfb34 r08e223  
    197197#define WannierTag3 172   //!< Message consists of matrix elements A_ij sent/received during evaluating of operator matrices
    198198#define WannierTag4 173   //!< Message consists of matrix elements B_ij sent/received during evaluating of operator matrices
     199#define WannierTag5 174   //!< Message consists of wannier centres and spread sent/received during ChangeWannierCentres()
    199200#define OtherPsiTag1 180  //!< Message consists of OnePsiElement data from the Spinype#SpinUp group
    200201#define OtherPsiTag2 181  //!< Message consists of OnePsiElement data from the Spinype#SpinDown group

  • pcp/src/wannier.c

    r9bfb34 r08e223  
    660660/** Changes Wannier Centres according to RunStruct#CommonWannier.
    661661 * \param *P Problem at hand.
    662  * \param Num number of Psis
    663  * \param **WannierCentre 2D array (NDIM, \a Num) with wannier centres
    664  * \param *WannierSpread array with wannier spread per wave function
    665  */
    666 void ChangeWannierCentres(struct Problem *P, int Num, double **WannierCentre, double *WannierSpread)
     662 */
     663void ChangeWannierCentres(struct Problem *P)
    667664{
    668665  struct RunStruct *R = &P->R;
    669666  struct Lattice *Lat = &P->Lat;
    670667  struct LatticeLevel *LevS = R->LevS;
     668  struct Psis *Psi = &Lat->Psi;
     669  struct OnePsiElement *OnePsiA;
    671670  int *marker, **group;
     671  int Num = Psi->NoOfPsis;
     672  double **WannierCentre = NULL;
     673  double *WannierSpread = NULL;
    672674  int partner[Num];
    673675  int i,j,l,k;
    … …  
    676678  double Spread;
    677679  int *N = LevS->Plan0.plan->N;
    678  
     680
     681  int MPISignal = 0, mpisignal = 0;
     682  unsigned int membersize, maxsize;
     683  int msgsize;
     684  double *buffer = NULL;
     685  int position;
     686  MPI_Status status;
     687 
     688  // determine necessary pack/buffer sizes for the wannier stuff
     689  MPI_Pack_size(3, MPI_DOUBLE, P->Par.comm_ST_PsiT, &membersize);
     690  maxsize = membersize;
     691  MPI_Pack_size(1, MPI_DOUBLE, P->Par.comm_ST_PsiT, &membersize);
     692  maxsize += membersize;
     693 
     694  // allocate memory
     695  buffer = (double *) Malloc(sizeof(double)*maxsize, "ChangeWannierCentres: *buffer");
     696  WannierSpread = (double *) Malloc(sizeof(double)*Num, "ChangeWannierCentres: *WannierSpread");
     697  WannierCentre = (double **) Malloc(sizeof(double *)*Num, "ChangeWannierCentres: **WannierCentre");
     698  for (j=0;j<Num;j++)
     699    WannierCentre[j] = (double *) Malloc(sizeof(double)*(NDIM+1), "ChangeWannierCentres: *WannierCentre[]");
     700
     701  // exchange wannier centers and spread among processes by going over all wave functions
     702  i=-1;
     703  for (l=0; l < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; l++) {
     704    OnePsiA = &Psi->AllPsiStatus[l];    // grab OnePsiA
     705    if (OnePsiA->PsiType == type) {   // drop all but occupied ones
     706      i++;  // increase l if it is occupied wave function
     707     
     708      // send and receive centres and spread (we use the MPI_Pack'ing scheme to preserve correct order of data)
     709      if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) {// is this local?
     710        position = 0;
     711        MPI_Pack(Psi->AddData[OnePsiA->MyLocalNo].WannierCentre, 3, MPI_DOUBLE, buffer, maxsize, &position, P->Par.comm_ST_PsiT);
     712        MPI_Pack(&Psi->AddData[OnePsiA->MyLocalNo].WannierSpread, 1, MPI_DOUBLE, buffer, maxsize, &position, P->Par.comm_ST_PsiT);
     713        for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++)
     714          if (k != P->Par.my_color_comm_ST_Psi) {  // send to all other processes
     715            fprintf(stderr, "(%i) Send goes from %d to %d\n", P->Par.me, k, P->Par.my_color_comm_ST_Psi);
     716            if (MPI_Send(buffer, position, MPI_PACKED, k, WannierTag5, P->Par.comm_ST_PsiT) != MPI_SUCCESS) {
     717              mpisignal = 1;
     718            }
     719          } else {  // "send" to us
     720            for(j=0;j<NDIM;j++)  // copy to own array as well
     721              WannierCentre[i][j] = Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j];
     722            WannierSpread[i] = Psi->AddData[OnePsiA->MyLocalNo].WannierSpread;
     723          }
     724      } else {
     725        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);
     726        if (MPI_Recv(buffer, maxsize, MPI_PACKED, OnePsiA->my_color_comm_ST_Psi, WannierTag5, P->Par.comm_ST_PsiT, &status) != MPI_SUCCESS)
     727          mpisignal = 1;
     728        position = 0;
     729        MPI_Get_count(&status, MPI_PACKED, &msgsize);
     730        MPI_Unpack (buffer, msgsize, &position, WannierCentre[i], 3, MPI_DOUBLE, P->Par.comm_ST_PsiT);
     731        MPI_Unpack (buffer, msgsize, &position, &WannierSpread[i], 1, MPI_DOUBLE, P->Par.comm_ST_PsiT);
     732      }
     733     
     734      // check for MPI errors
     735      if (MPI_Allreduce(&mpisignal, &MPISignal, 1, MPI_INT, MPI_SUM, P->Par.comm_ST_PsiT) != MPI_SUCCESS) // allreduce signals
     736        Error(SomeError,"ChangeWannierCentres: Bcast of signal failed\n");
     737      if (MPISignal) { // greater 0 means there was at least one hickup somewhere
     738        Error(SomeError,"ChangeWannierCentres: Scattering of Wanniers failed\n");
     739      }
     740    }
     741  }
     742  // free the buffer memory
     743  Free(buffer, "ChangeWannierCentres: *buffer");
     744
     745  if ((P->Call.out[NormalOut]) && (1)) { // root prints gathered wannier centres for debugging
     746    fprintf(stderr, "(%i) Printing gathered Wannier centres:\n", P->Par.me);
     747    for(l=0;l<Num;l++) {
     748      fprintf(stderr, "(%i) ", P->Par.me);
     749      for(j=0;j<NDIM;j++)
     750        fprintf(stderr, "%lg\t", WannierCentre[l][j]);
     751      fprintf(stderr, "-\t%lg\n", WannierSpread[l]);
     752    }
     753  }
     754     
     755  // change centres
    679756  switch (R->CommonWannier) {
    680757   case 4:
    … …  
    728805          Spread += (WannierCentre[l][i] - WannierCentre[k][i])*(WannierCentre[l][i] - WannierCentre[k][i]);
    729806        }
    730         Spread = sqrt(Spread);  // distance in Spread
     807        //Spread = sqrt(Spread);  // distance in Spread
    731808        //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]);
    732         if (Spread < 1.5*(WannierSpread[l]+WannierSpread[k])) {// if distance smaller than sum of spread
     809        if (Spread < 5*(WannierSpread[l]+WannierSpread[k])*(WannierSpread[l]+WannierSpread[k])) {// if distance smaller than sum of spread times two plus safety
    733810          group[l][j++] = k;  // add k to group of l
    734811          partner[l]++;
    … …  
    821898   break;
    822899  }
     900
     901  if ((P->Call.out[NormalOut]) && (1)) { // root prints gathered wannier centres for debugging
     902    fprintf(stderr, "(%i) Printing changed Wannier centres:\n", P->Par.me);
     903    for(l=0;l<Num;l++) {
     904      fprintf(stderr, "(%i) ", P->Par.me);
     905      for(j=0;j<NDIM;j++)
     906        fprintf(stderr, "%lg\t", WannierCentre[l][j]);
     907      fprintf(stderr, "-\t%lg\n", WannierSpread[l]);
     908    }
     909  }
     910     
     911  // store info into local arrays again
     912  i=-1;
     913  for (l=0; l < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; l++) {  // go through all wave functions
     914    OnePsiA = &Psi->AllPsiStatus[l];    // grab OnePsiA
     915    if (OnePsiA->PsiType == type) {   // drop all but occupied ones
     916      i++;  // increase l if it is occupied wave function
     917      if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) {// is this local
     918        for(j=0;j<NDIM;j++)
     919          Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j] = WannierCentre[i][j];
     920      }
     921    }
     922  }
     923
     924  // show locally stored info for debugging
     925  if (P->Call.out[NormalOut]) {
     926    i=-1;
     927    for (l=0; l < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; l++) {  // go through all wave functions
     928      OnePsiA = &Psi->AllPsiStatus[l];    // grab OnePsiA
     929      if (OnePsiA->PsiType == type) {   // drop all but occupied ones
     930        i++;  // increase l if it is occupied wave function
     931        if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // is this local?
     932           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);
     933      }
     934    }
     935  }
     936
     937  // free memory
     938  Free(WannierSpread, "ChangeWannierCentres: *WannierSpread");
     939  for (j=0;j<Num;j++)
     940    Free(WannierCentre[j], "ChangeWannierCentres: *WannierCentre[]");
     941  Free(WannierCentre, "ChangeWannierCentres: **WannierCentre");
    823942}
    824943
    … …  
    11751294  ComputeWannierCentresfromVarianceMatrices(P, DiagData.A, &spread, &spreadSQ, WannierCentre, WannierSpread);
    11761295
    1177   // join Wannier orbital to groups with common centres under certain conditions
    1178   //debug (P,"Changing Wannier Centres according to CommonWannier");
    1179   ChangeWannierCentres(P, Num, WannierCentre, WannierSpread);
    1180  
    11811296  // write to file
    11821297  //debug (P,"Writing spread file");
    … …  
    18301945        for (j=0;j<NDIM;j++)
    18311946          Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j] = WannierCentre[i][j];
     1947        Psi->AddData[OnePsiA->MyLocalNo].WannierSpread = WannierSpread[i];
    18321948      }
    18331949      if (P->Par.me_comm_ST == 0)
    … …  
    19342050          //return 0;
    19352051        // and store 'em (for all who have this Psi local)
    1936         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]);
     2052        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]);
    19372053        for (j=0;j<NDIM;j++) Psi->AddData[OnePsiA->MyLocalNo].WannierCentre[j] = WannierCentre[j];
    19382054        Psi->AddData[OnePsiA->MyLocalNo].WannierSpread = WannierCentre[NDIM];

  • pcp/src/wannier.h

    r9bfb34 r08e223  
    3434};
    3535
     36
    3637void PrintGSLMatrix(struct Problem *P, gsl_matrix *U, int Num, const char *msg);
    3738void ComputeMLWF(struct Problem *P);
    3839void WriteWannierFile(struct Problem *P, double spread, double old_spread, double **WannierCentre, double *WannierSpread);
    3940int ParseWannierFile(struct Problem *P);
     41void ChangeWannierCentres(struct Problem *P);
    4042void SerialDiagonalization(struct Problem *P, struct DiagonalizationData *DiagData);
    4143void ParallelDiagonalization(struct Problem *P, struct DiagonalizationData *DiagData);
Note: See TracChangeset for help on using the changeset viewer.