/** \file output.c * Output of forces and energies, Visuals (density and ions for OpenDX). * * Herein all the functions concerning the output of data are gathered: * Initialization of the FileData structure InitOutVisArray() and the files themselves * InitOutputFiles(), * Saving of a calculated state OutputVisSrcFiles() - 1. Psis OutputSrcPsiDensity(), 2. Ions * OutSrcIons() - or retrieving Psis ReadSrcPsiDensity() and Ions ReadSrcIons(), * Preparation (in case of RiemannTensor use) CalculateOutVisDensityPos(), OutVisPosRTransformPosNFRto0(), * Output of visual data (for OpenDx explorer) OutputVis() - 1. OpenDX files CreateDensityOutputGeneral(), * 2. electronic densitiy data OutVisDensity() (uses MPI sending density coefficients OutputOutVisDensity() * and receiving and writing CombineOutVisDensity()), 3. ion data OutVisIons() - * Closing of files CloseOutputFiles(). * * There are some more routines: OutputCurrentDensity() outputs the current density for each magnetic field * direction, OutputVisAllOrbital() outputs all orbitals of a certain minimsation type and * TestReadnWriteSrcDensity() checks whether a certain minimisation group can be written and read again * correctly. * * There are some helpers that open files with certain filenames, making extensive usage of all * the suffixes defined in here: OpenFile(), OpenFileNo(), OpenFileNo2(), OpenFileNoNo() and * OpenFileNoPost(). * * Project: ParallelCarParrinello \author Jan Hamaekers, Frederik Heber \date 2000, 2006 File: output.c $Id: output.c,v 1.51.2.2 2007-04-21 12:55:50 foo Exp $ */ #include #include #include #include //#include #include #include #include"data.h" #include"density.h" #include"errors.h" #include"gramsch.h" #include"helpers.h" #include "init.h" #include"myfft.h" #include"mymath.h" #include"output.h" #include"pdbformat.h" #include"perturbed.h" /* Konvention: Rueckgabe 0 einer Funktion, bedeutet keinen Fehler (entsprechend exitcode 0) */ /* Oeffnet Datei P->mainname+"..." mit what*/ /** Opens a file with FileData::mainname+suffix. * Both the path under which the main programme resides and the path to the config file are * tried subsequently. * \param *P Problem at hand * \param **file file pointer array * \param *suffix suffix after mainname * \param *what access parameter for the file: r, w, rw, r+, w+ ... * \param verbose 1 - print status, 0 - don't print anything * \return 0 - could not open file, 1 - file is open */ int OpenFile(struct Problem *P, FILE** file, const char* suffix, const char* what, int verbose) { char* name; /* Zu erzeugender Dateiname */ name = (char*) Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 1,"OpenFile"); sprintf(name, "%s%s%s", P->Files.default_path, P->Files.mainname, suffix); *file = fopen(name, what); if (*file == 0) { fprintf(stderr,"(%i) Normal access failed: name %s, what %s\n", P->Par.me, name, what); /* hier default file ausprobieren, falls nur gelesen werden soll! */ if(*what == 'r') { name = (char*) Realloc(name,strlen(P->Files.default_path) + strlen(suffix)+1,"OpenFile"); sprintf(name, "%s%s", P->Files.default_path,suffix); *file = fopen(name,what); if (*file != 0) { if (verbose) fprintf(stderr,"(%i) Default file is open: %s\n",P->Par.me,name); Free(name, "OpenFile: name"); return(1); } } fprintf(stderr,"\n(%i)Error: Cannot open neither normal nor default file: %s\n",P->Par.me,name); Free(name, "OpenFile: name"); return(0); } else { if(verbose) fprintf(stderr,"(%i) File is open: %s\n",P->Par.me,name); Free(name, "OpenFile: name"); return(1); } } /** Opens a file with FileData::mainname+suffix+"."+No (2 digits). * Both the path under which the main programme resides and the path to the config file are * tried subsequently. * \param *P Problem at hand * \param **file file pointer array * \param *suffix suffix after mainname * \param No the number with up to two digits * \param *what access parameter for the file: r, w, rw, r+, w+ ... * \param verbose 1 - print status, 0 - don't print anything * \return 0 - could not open file, 1 - file is open */ int OpenFileNo2(struct Problem *P, FILE** file, const char* suffix, int No, const char* what, int verbose) { char* name; /* Zu erzeugender Dateiname */ name = (char*) Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 4,"OpenFile"); sprintf(name, "%s%s%s.%02i", P->Files.default_path, P->Files.mainname, suffix, No); *file = fopen(name, what); if (*file == 0) { /* falls nur gelesen wird, auch default_path ausprobieren */ if (*what == 'r') { name = (char*) Realloc(name,strlen(P->Files.default_path) + strlen(suffix) + 4,"OpenFileNo2"); sprintf(name, "%s%s.%02i", P->Files.default_path, suffix, No); *file = fopen(name, what); if (*file != 0) { if(verbose) fprintf(stderr,"(%i) Default file is open: %s\n", P->Par.me, name); Free(name, "OpenFileNo2: name"); return(1); } } fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name); Free(name, "OpenFileNo2: name"); return(0); } else { if(verbose) fprintf(stderr,"(%i) File is open: %s\n",P->Par.me, name); Free(name, "OpenFileNo2: name"); return(1); } } /* Oeffnet Datei P->Files.mainname+"...".Nr(4stellig) mit what*/ /** Opens a file with FileData::mainname+suffix+"."+No (4 digits). * Both the path under which the main programme resides and the path to the config file are * tried subsequently. * \param *P Problem at hand * \param **file file pointer array * \param *suffix suffix after mainname * \param No the number with up to four digits * \param *what access parameter for the file: r, w, rw, r+, w+ ... * \param verbose 1 - print status, 0 - don't print anything * \return 0 - could not open file, 1 - file is open */ int OpenFileNo(struct Problem *P, FILE** file, const char* suffix, int No, const char* what, int verbose) { char* name; /* Zu erzeugender Dateiname */ name = (char*) Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 6,"OpenFileNo"); sprintf(name, "%s%s%s.%04i", P->Files.default_path, P->Files.mainname, suffix, No); *file = fopen(name, what); if (*file == 0) { /* falls nur gelesen wird, auch default_path ausprobieren */ if (*what == 'r') { name = (char*) Realloc(name,strlen(P->Files.default_path) + strlen(suffix) + 6,"OpenFileNo"); sprintf(name, "%s%s.%04i", P->Files.default_path, suffix, No); *file = fopen(name, what); if (*file != 0) { if(verbose) fprintf(stderr,"(%i) Default file is open: %s\n", P->Par.me, name); Free(name, "OpenFileNo: name"); return(1); } } fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name); Free(name, "OpenFileNo: name"); return(0); } else { if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name); Free(name, "OpenFileNo: name"); return(1); } } /* die nachfolgende Routine wird von seq und par benutzt!!! */ /* Oeffnet Datei P->Files.mainname+"...".No.postfix mit what*/ /** Opens a file with FileData::mainname+suffix+"."+No(4 digits)+postfix. * Only the path under which the main programme resides is tried. * \param *P Problem at hand * \param **file file pointer array * \param *suffix suffix after mainname * \param No the number with up to four digits * \param *postfix post-suffix at the end of filename * \param *what access parameter for the file: r, w, rw, r+, w+ ... * \param verbose 1 - print status, 0 - don't print anything * \return 0 - could not open file, 1 - file is open */ int OpenFileNoPost(struct Problem *P, FILE** file, const char* suffix, int No, const char* postfix, const char* what, int verbose) { char* name; /* Zu erzeugender Dateiname */ name = (char*) Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(postfix) + strlen(suffix) + 6,"OpenFileNoPost"); sprintf(name, "%s%s%s.%04i%s", P->Files.default_path, P->Files.mainname, suffix, No, postfix); *file = fopen(name, what); if (*file == 0) { fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name); Free(name, "OpenFileNoPost: name"); return(0); } else { if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name); Free(name, "OpenFileNoPost: name"); return(1); } } /** Opens a file with FileData::mainname+suffix+"."+No1(4 digits)+"."+No2(4 digits). * Only the path under which the main programme resides is tried. * \param *P Problem at hand * \param **file file pointer array * \param *suffix suffix after mainname * \param No1 first number with up to four digits * \param No2 second number with up to four digits * \param *what access parameter for the file: r, w, rw, r+, w+ ... * \param verbose 1 - print status, 0 - don't print anything * \return 0 - could not open file, 1 - file is open */ int OpenFileNoNo(struct Problem *P, FILE** file, const char* suffix, int No1, int No2, const char* what, int verbose) { /* zuerst suffix; No1: lfd, No2: procId */ char *name; name = (char*) Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 5 + 6 + 1,"OpenFileNoNo"); sprintf(name,"%s%s%s.%04i.%04i", P->Files.default_path, P->Files.mainname, suffix, No1, No2); *file = fopen(name, what); if(*file == 0) { fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name); Free(name, "OpenFileNoNo: name"); return(0); } if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name); Free(name, "OpenFileNoNo: name"); return(1); } /** Transformation of wave function of Riemann level to zeroth before output as Vis. * \param *RT RiemannTensor structure, contains RiemannTensor::LatticeLevel * \param *source source wave function array * \param *dest destination wave function array * \param NF dimensional factor (NDIM, generally 3) */ static void OutVisPosRTransformPosNFRto0(const struct RiemannTensor *RT, fftw_real *source, fftw_real *dest, const int NF) { struct LatticeLevel *Lev = RT->LevR; int es = Lev->NUp0[2]*NF; unsigned int cpyes = sizeof(fftw_real)*es; int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp0[0],Ny=Lev->NUp0[1]; int lx,ly,z,Lx,Ly; for(lx=0; lx < nx; lx++) for(Lx=0; Lx < Nx; Lx++) for(ly=0; ly < ny; ly++) for(Ly=0; Ly < Ny; Ly++) for(z=0; z < nz; z++) { memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))], &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))], cpyes); } } /** prints Norm of each wave function to screen. * \param *out output destination (eg. stdout) * \param *P Problem at hand */ void OutputNorm (FILE *out, struct Problem *P) { struct Lattice *Lat = &P->Lat; struct Psis *Psi = &Lat->Psi; struct LatticeLevel *Lev = P->R.LevS; struct OnePsiElement *OnePsi, *LOnePsi; int i; for (i=0;iMaxPsiOfType+P->Par.Max_me_comm_ST_PsiT;i++) { OnePsi =&Psi->AllPsiStatus[i]; if (OnePsi->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // local one LOnePsi = &Psi->LocalPsiStatus[OnePsi->MyLocalNo]; fprintf(out,"(%i) Norm of Psi %i: %e\n", P->Par.me, OnePsi->MyGlobalNo, GramSchGetNorm2(P,Lev,Lev->LPsi->LocalPsi[OnePsi->MyLocalNo])); } } } /** Output of current Psis state of source level RunStruct::LevS's LevelPsi::LocalPsi to file \ref suffixsrcpsidat. * In case of process 0, the doc file is written in a parsable way with minimisation type RunStruct#CurrentMin, level number * LatticeLevel#LevelNo, and number of grid nodes LatticeLevel#N. * The two (three) different Psis::SpinType's are discerned and in case of SpinUpDown separate data files opened. * * Then for every global wave function of the desired type each coefficient of the reciprocal grid is copied into * Density::DensityCArray[TempDensity], fouriertransformed, copied into Density::DensityArray[TempDensity]. * * As the file output is only handled by process 0, all other coefficient shares are sent within the Psi group to which * the current global wave function belongs to process 0 of the Psi group and from there on to process 0 of all via * MPI. * \param *P Problem at hand * \param type Current minimisation state * \note This serves as a backup file, when the process is terminated and one would like to restart it at the * current calculation lateron, see ReadSrcPsiDensity(). Note also that it is not necessary to specify the * same number of processes on later restart, any number may be used under the condition that the number of * grid nodes match and that there 2 sharing wave functions in case of SpinUpDown. * \sa ReadSrcPsiDensity() - same for Reading the coefficients, TestReadnWriteSrcDensity() - checks both routines against * each other */ void OutputSrcPsiDensity(struct Problem *P, enum PsiTypeTag type) { int i,j,k, Index, zahl, owner; struct Lattice *Lat = &P->Lat; //struct RunStruct *R = &P->R; struct fft_plan_3d *plan = Lat->plan; struct LatticeLevel *LevS = P->R.LevS; struct Psis *Psi = &Lat->Psi; fftw_complex *work = (fftw_complex *)LevS->Dens->DensityArray[TempDensity]; double *destpsiR = (double *)LevS->Dens->DensityArray[TempDensity]; fftw_real *srcpsiR = (fftw_real *)LevS->Dens->DensityCArray[TempDensity]; fftw_complex *srcpsiC = (fftw_complex *)LevS->Dens->DensityCArray[TempDensity]; FILE *SrcPsiData, *SrcPsiDoc; char suffixdat[255], suffixdoc[255]; MPI_Status status; struct OnePsiElement *OnePsiA, *LOnePsiA; fftw_complex *LPsiDatA; int Num = 0, colorNo = 0; int Sent = 0, sent = 0; SpeedMeasure(P,ReadnWriteTime,StartTimeDo); sprintf(suffixdat, ".%.254s.L%i", P->R.MinimisationName[type], LevS->LevelNo); strncpy (suffixdoc, suffixdat, 255); // for the various spin cases, output the doc-file if it's process 0 if (!(P->Par.me_comm_ST)) { switch (Lat->Psi.PsiST) { case SpinDouble: colorNo = 0; strncat (suffixdat, suffixsrcpsidat, 255-strlen(suffixdat)); strncat (suffixdoc, suffixsrcpsidoc, 255-strlen(suffixdoc)); Num = Lat->Psi.GlobalNo[PsiMaxNoDouble]; break; case SpinUp: colorNo = 0; strncat (suffixdat, suffixsrcpsiupdat, 255-strlen(suffixdat)); strncat (suffixdoc, suffixsrcpsiupdoc, 255-strlen(suffixdoc)); Num = Lat->Psi.GlobalNo[PsiMaxNoUp]; break; case SpinDown: colorNo = 1; strncat (suffixdat, suffixsrcpsidowndat, 255-strlen(suffixdat)); strncat (suffixdoc, suffixsrcpsidowndoc, 255-strlen(suffixdoc)); Num = Lat->Psi.GlobalNo[PsiMaxNoDown]; break; } if (!OpenFileNo(P, &SrcPsiData, suffixdat, colorNo, "wb",P->Call.out[ReadOut])) // open SourcePsiData as write binary fprintf(stderr,"(%i) Error opening file with suffix %s for writing!\n",P->Par.me, suffixdat); if (!(P->Par.my_color_comm_ST_Psi)) { // if we are process 0 if (!OpenFile(P, &SrcPsiDoc, suffixdoc, "w",P->Call.out[ReadOut])) // open the (text) doc file fprintf(stderr,"(%i) Error opening file with suffix %s for writing!\n",P->Par.me, suffixdoc); fprintf(SrcPsiDoc, "Mintype\t%i\n", (int)type); fprintf(SrcPsiDoc, "LevelNo\t%i\n", LevS->LevelNo); fprintf(SrcPsiDoc, "GridNodes\t%i\t%i\t%i\n", LevS->N[0], LevS->N[1], LevS->N[2]); fprintf(SrcPsiDoc, "PsiNo\t%i\t%i\n", Num, P->Lat.Psi.GlobalNo[PsiMaxAdd]); fprintf(SrcPsiDoc, "Epsilon\t%lg\t%lg\n", P->R.RelEpsTotalEnergy, P->R.RelEpsKineticEnergy); for (i = 0; i < P->Par.Max_my_color_comm_ST_Psi; i++) { fprintf(SrcPsiDoc, "\t%i", Lat->Psi.RealAllLocalNo[i]); // print number of local Psis } fprintf(SrcPsiDoc, "\n"); fclose(SrcPsiDoc); } } // send/receive around and write share of coefficient array of each wave function MPI_Allreduce(&sent, &Sent, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); // catch all at the starter line fprintf(stderr,"(%i) me (%i/%i) \t Psi (%i/%i)\t PsiT (%i/%i)\n", P->Par.me, P->Par.me_comm_ST, P->Par.Max_me_comm_ST, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi, P->Par.me_comm_ST_PsiT, P->Par.Max_me_comm_ST_PsiT); k = -1; // k is global PsiNo counter for the desired group for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) { // go through all wave functions (plus the extra one for each process) OnePsiA = &Psi->AllPsiStatus[j]; // grab OnePsiA if (OnePsiA->PsiType == type) { // only take desired minimisation group k++; owner = 0; // notes down in process 0 of each psi group the owner of the next coefficient share //fprintf(stderr,"(%i) ST_Psi: OnePsiA %i\tP->Par.me %i\n", P->Par.me,OnePsiA->my_color_comm_ST_Psi,P->Par.my_color_comm_ST_Psi); if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // Belongs to my Psi group? LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo]; LPsiDatA = LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo]; SetArrayToDouble0((double *)srcpsiR, LevS->Dens->TotalSize*2); // zero DensityCArray[TempDensity] for (i=0;iMaxG;i++) { // for each every unique G grid vector Index = LevS->GArray[i].Index; srcpsiC[Index].re = LPsiDatA[i].re; // copy real value srcpsiC[Index].im = LPsiDatA[i].im; // copy imaginary value } for (i=0; iMaxDoubleG; i++) { // also for every doubly appearing G vector (symmetry savings) srcpsiC[LevS->DoubleG[2*i+1]].re = srcpsiC[LevS->DoubleG[2*i]].re; srcpsiC[LevS->DoubleG[2*i+1]].im = -srcpsiC[LevS->DoubleG[2*i]].im; } // do an fft transform from complex to real on these srcPsiR fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNF1, srcpsiC, work); for (i=0; i < LevS->Dens->LocalSizeR; i++) destpsiR[i] = (double)srcpsiR[i]; } else LOnePsiA = NULL; if (P->Par.me_comm_ST == 0) { // if we are process 0 of all, only we may access the file for (i=0; iPar.Max_me_comm_ST_Psi;i++) { // for each share of the coefficient in the PsiGroup if (LOnePsiA == NULL) { // if it's not local, receive coefficients from correct PsiGroup (process 0 within that group) if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, OnePsiA->my_color_comm_ST_Psi, OutputSrcPsiTag, P->Par.comm_ST_PsiT, &status ) != MPI_SUCCESS) Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!"); MPI_Get_count(&status, MPI_DOUBLE, &zahl); if (zahl != LevS->Dens->LocalSizeR) // check number of elements fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, i, zahl, LevS->Dens->LocalSizeR); //else //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, i); } else { // if it's local ... if (i != 0) { // but share of array not for us, receive from owner process within Psi group if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, i, OutputSrcPsiTag, P->Par.comm_ST_Psi, &status ) != MPI_SUCCESS) Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!"); MPI_Get_count(&status, MPI_DOUBLE, &zahl); if (zahl != LevS->Dens->LocalSizeR) // check number of elements fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, i, zahl, LevS->Dens->LocalSizeR); //else //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, i); } // otherwise it was our share already } // store the final share on disc if ((zahl = fwrite(destpsiR, sizeof(double), (size_t)(LevS->Dens->LocalSizeR), SrcPsiData)) != (size_t)(LevS->Dens->LocalSizeR)) { fclose(SrcPsiData); //if (P->Par.me == 0) fprintf(stderr, "(%i)OutputSrcPsiDensity: only %i bytes written instead of expected %i\n", P->Par.me, zahl, LevS->Dens->LocalSizeR); Error(SomeError,"OutputSrcPsiDensity: fwrite Error"); } } } else { // if we are not process 0 of all, we are but a deliverer if (LOnePsiA != NULL) { // send if it's local if (P->Par.me_comm_ST_Psi == 0) { // if we are process 0 in the group, send final share to our process 0 for (owner = 0; owner < P->Par.Max_me_comm_ST_Psi; owner++) { // for all processes of our Psi group if (owner != 0) { // still not our share of coefficients, receive from owner in our Psi group (increasing owner) if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, owner, OutputSrcPsiTag, P->Par.comm_ST_Psi, &status ) != MPI_SUCCESS) Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!"); MPI_Get_count(&status, MPI_DOUBLE, &zahl); if (zahl != LevS->Dens->LocalSizeR) // check number of elements fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, owner, zahl, LevS->Dens->LocalSizeR); //else //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, owner); } else sent++; // only count sent if it was our share if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, OutputSrcPsiTag, P->Par.comm_ST_PsiT) != MPI_SUCCESS) Error(SomeError, "OutputSrcPsiDensity: MPI_Send of loaded coefficients failed!"); //else //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Send to process %i in PsiT group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, 0, k); } } else { sent++; if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, OutputSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS) Error(SomeError, "OutputSrcPsiDensity: MPI_Send of loaded coefficients failed!"); //else //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Send to process %i in Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, 0, k); } } // otherwise we don't have anything to do with this } } } MPI_Allreduce(&sent, &Sent, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); // catch all again at finish fprintf(stderr,"(%i) Out of %i shares %i had to be sent in total, %i from this process alone.\n", P->Par.me, P->Par.Max_me_comm_ST_Psi*Psi->NoOfPsis, Sent, sent); if (!(P->Par.me_comm_ST)) fclose(SrcPsiData); SpeedMeasure(P,ReadnWriteTime,StopTimeDo); } /** Tests whether writing and successant reading of coefficient array is working correctly. * The local wave function array is written to a disc (\sa OutputSrcPsiDensity()), the by \a wavenr * specified coefficient array copied to OldPsiDat and afterwards read again via ReadSrcPsiDensity(). * Comparison per process of each local coefficient shows incorrect read or writes. * \param *P Problem at hand * \param type minimisation type array to test for read and write * \return 1 - successful, 0 - test failed */ int TestReadnWriteSrcDensity(struct Problem *P, enum PsiTypeTag type) { struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; struct Lattice *Lat = &P->Lat; struct Psis *Psi = &Lat->Psi; int i,k; fftw_complex *destpsiC, *srcpsiC; //fprintf(stderr,"(%i)TestReadnWriteSrcDensity\n",P->Par.me); // write whole array of type to disc OutputSrcPsiDensity(P,type); //debug(P,"array written"); // copy specified array to OldPsiDat for (k=0; k < Lat->Psi.LocalNo; k++) // for every local wave function of type, copy coefficients if (Psi->LocalPsiStatus[k].PsiType == type) { // ... yet only for given type srcpsiC = LevS->LPsi->LocalPsi[k]; destpsiC = LevS->LPsi->OldLocalPsi[k - Psi->TypeStartIndex[type]]; for (i=0;iMaxG;i++) { // for each every unique G grid vector destpsiC[i].re = srcpsiC[i].re; // copy real value destpsiC[i].im = srcpsiC[i].im; // copy imaginary value } } //debug(P,"array copied"); // read whole array again if (!ReadSrcPsiDensity(P,type,0,R->LevSNo)) return 0; //debug(P,"array read"); // compare with copied array for (k=0; k < Lat->Psi.LocalNo; k++) // for every local wave function of type, compare coefficients if (Psi->LocalPsiStatus[k].PsiType == type) { // ... yet only for given type srcpsiC = LevS->LPsi->LocalPsi[k]; destpsiC = LevS->LPsi->OldLocalPsi[k - Psi->TypeStartIndex[type]]; for (i=0;iMaxG;i++) // for each every unique G grid vector if ((fabs(destpsiC[i].re - srcpsiC[i].re) >= MYEPSILON) ||(fabs(destpsiC[i].im - srcpsiC[i].im) >= MYEPSILON)) { fprintf(stderr,"(%i)TestReadnWriteSrcDensity: First difference at index %i - %lg+i%lg against loaded %lg+i%lg\n",P->Par.me, i, srcpsiC[i].re, srcpsiC[i].im,destpsiC[i].re,destpsiC[i].im); return 0; } } //debug(P,"array compared"); fprintf(stderr,"(%i)TestReadnWriteSrcDensity: OK!\n",P->Par.me); return 1; } /** Read Psis state from an earlier run. * The doc file is opened, mesh sizes LatticeLevel::N[], global number of Psis read and checked against the known * values in the inital level RunStruct::LevS. * Note, only process 0 handles the files, all read coefficients are transfered to their respective owners via MPI * afterwards. Here, process 0 of a certain Psi group is used as a transferer for the coefficients of the other processes * in this Psi group. He receives them all from process 0 and sends them onward accordingly. The complete set of * coefficients on the real grid for one wave function in the Psi group are transformed into complex wave function * coefficients by the usual fft procedure (see ChangeToLevUp()). * * \param *P Problem at hand * \param type minimisation type to read * \param test whether to just test for presence of files (1) or really read them (0) * \param LevSNo level number to be read * \note This is the counterpart to OutputSrcPsiDensity(). * \return 1 - success, 0 - failure * \note It is not necessary to specify the same number of processes on later restart, any number may be used under * the condition that the number of grid nodes match and that there at least 2 processes sharing wave functions * in case of SpinUpDown. * \sa OutputSrcPsiDensity() - same for writing the coefficients, TestReadnWriteSrcDensity() - checks both routines against * each other */ int ReadSrcPsiDensity(struct Problem *P, enum PsiTypeTag type, int test, int LevSNo) { int i, j, k, Index, owner; struct RunStruct *R = &P->R; struct Lattice *Lat = &P->Lat; struct fft_plan_3d *plan = Lat->plan; struct LatticeLevel *LevS = R->LevS; // keep open for LevelNo read from file struct Psis *Psi = &Lat->Psi; //struct Energy *E = Lat->E; fftw_complex *work; double *destpsiR; fftw_real *srcpsiR; fftw_complex *srcpsiC; FILE *SrcPsiData, *SrcPsiDoc; int N[NDIM], GlobalNo[2]; int LevelNo, readnr=0; int zahl, signal = test ? 1 : 2; // 0 - ok, 1 - test failed, 2 - throw Error char suffixdat[255], suffixdoc[255]; int read_type, Num = 0, colorNo = 0; char spin[20]; double Eps[2]; MPI_Status status; struct OnePsiElement *OnePsiA, *LOnePsiA; int Recv=0, recv=0; SpeedMeasure(P,ReadnWriteTime,StartTimeDo); sprintf(suffixdat, ".%.254s.L%i", P->R.MinimisationName[type], LevSNo); strncpy (suffixdoc, suffixdat, 255); // Depending on Psis::SpinType the source psi doc file is opened and header written switch (Lat->Psi.PsiST) { case SpinDouble: colorNo = 0; strncat (suffixdat, suffixsrcpsidat, 255-strlen(suffixdat)); strncat (suffixdoc, suffixsrcpsidoc, 255-strlen(suffixdoc)); strncpy (spin, "GlobalNoSpinDouble", 20); Num = Lat->Psi.GlobalNo[PsiMaxNoDouble]; break; case SpinUp: colorNo = 0; strncat (suffixdat, suffixsrcpsiupdat, 255-strlen(suffixdat)); strncat (suffixdoc, suffixsrcpsiupdoc, 255-strlen(suffixdoc)); strncpy (spin, "GlobalNoSpinUp", 20); Num = Lat->Psi.GlobalNo[PsiMaxNoUp]; break; case SpinDown: colorNo = 1; strncat (suffixdat, suffixsrcpsidowndat, 255-strlen(suffixdat)); strncat (suffixdoc, suffixsrcpsidowndoc, 255-strlen(suffixdoc)); strncpy (spin, "GlobalNoSpinDown", 20); Num = Lat->Psi.GlobalNo[PsiMaxNoDown]; break; } // open doc file ... if (!(P->Par.me_comm_ST)) { if (!OpenFile(P, &SrcPsiDoc, suffixdoc, "r", test ? 0 : P->Call.out[ReadOut])) { // open doc file debug(P,"ReadSrcPsiDensity: doc file pointer NULL\n"); if (test) { if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); return 0; } if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); Error(SomeError,"ReadSrcPsiDensity: cannot open doc file!"); } // ... and parse critical ... readnr += ParseForParameter(0,SrcPsiDoc,"Mintype",0,1,1,int_type,(int *)&read_type, 1, test ? optional : critical); readnr += ParseForParameter(0,SrcPsiDoc,"LevelNo",0,1,1,int_type,&LevelNo,1, test ? optional : critical); readnr += 3*ParseForParameter(0,SrcPsiDoc,"GridNodes",0,3,1,row_int,&N[0], 1, test ? optional : critical); readnr += 2*ParseForParameter(0,SrcPsiDoc,"PsiNo",0,2,1,row_int,&GlobalNo[0], 1, test ? optional : critical); // and optional items ... if (ParseForParameter(0,SrcPsiDoc,"Epsilon",0,2,1,row_double,&Eps[0],1,optional)) if ((P->Call.ReadSrcFiles == 1) && ((Eps[1] < R->RelEpsKineticEnergy) || (Eps[0] < R->RelEpsTotalEnergy))) { //fprintf(stderr,"(%i) Eps %lg %lg\tRelEps %lg %lg\n", P->Par.me, Eps[0], Eps[1], R->RelEpsTotalEnergy, R->RelEpsKineticEnergy); fprintf(stderr,"(%i) NOTE: Doing minimization after source file parsing due to smaller specified epsilon stop conditions.\n",P->Par.me); P->Call.ReadSrcFiles = 2; // do minimisation even after loading } if (readnr != 7) { // check number of items read debug(P, "ReadSrcPsiDensity: too few doc items in file\n"); fclose(SrcPsiDoc); if (test) { if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); return 0; } fprintf(stderr,"ReadSrcPsiDensity: Only %i items read!\n",readnr); if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); Error(SomeError, "ReadSrcPsiDensity: read error"); } // check if levels match if (LevSNo != LevelNo) { debug(P,"ReadSrcPsiDensity: mismatching levels\n"); fclose(SrcPsiDoc); if (test) { if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); return 0; } if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); Error(SomeError,"ReadSrcPsiDensity: Mismatching levels!"); } else { LevS = &P->Lat.Lev[LevelNo]; } // check if systems in memory and file match if ((read_type != R->CurrentMin) || (N[0] != LevS->N[0]) || (N[1] != LevS->N[1]) || (N[2] != LevS->N[2]) || (GlobalNo[0] != Num) || (GlobalNo[1] != Lat->Psi.GlobalNo[PsiMaxAdd])) { // check read system debug(P,"ReadSrcPsiDensity: srcpsi file does not fit to system\n"); fclose(SrcPsiDoc); if (test) { fprintf(stderr,"(%i) Min\t N(x,y,z)\tPsiNo+AddNo\n file: %s\t %i %i %i\t %i + %i\nsystem: %s\t %d %d %d\t %d + %d\n",P->Par.me, R->MinimisationName[read_type], N[0], N[1], N[2], GlobalNo[0], GlobalNo[1], R->MinimisationName[R->CurrentMin], LevS->N[0] , LevS->N[1], LevS->N[2], Lat->Psi.GlobalNo[PsiMaxNoDouble], Lat->Psi.GlobalNo[PsiMaxAdd]); if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); return 0; } fprintf(stderr,"ReadSrcPsiDensity: Type %i != CurrentMin %i || N[0] %i != %i || N[1] %i != %i || N[2] %i != %i || %s %i + %i != %i + %i\n", read_type, R->CurrentMin, N[0], LevS->N[0], N[1], LevS->N[1], N[2], LevS->N[2], spin, GlobalNo[0], GlobalNo[1], Num, P->Lat.Psi.GlobalNo[PsiMaxAdd]); if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); Error(SomeError,"ReadSrcPsiDensity: srcpsi file does not fit to system"); } signal = 0; // everything went alright, signal ok if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); } else { // others wait for signal from root process if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n"); if (signal == 1) return 0; else if (signal == 2) Error(SomeError, "ReadSrcPsiDensity: Something went utterly wrong, see root process"); else fprintf(stderr,"(%i) ReadSrcPsiDensity: Everything went alright so far\n", P->Par.me); } if (MPI_Bcast(&LevelNo,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS) Error(SomeError,"ReadSrcPsiDensity: Bcast of LevelNo failed\n"); LevS = &P->Lat.Lev[LevelNo]; //if (!test) fprintf(stderr,"(%i) LevelSNo %i\n", P->Par.me, LevS->LevelNo); if (!test) { // set some pointers for work to follow work = (fftw_complex *)LevS->Dens->DensityArray[TempDensity]; destpsiR = (double *)LevS->Dens->DensityArray[TempDensity]; srcpsiR = (fftw_real *)LevS->Dens->DensityCArray[TempDensity]; srcpsiC = (fftw_complex *)LevS->Dens->DensityCArray[TempDensity]; // read share of coefficient array for each wave function and send/receive around owner = 0; MPI_Allreduce (&recv, &Recv, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); //fprintf(stderr,"(%i) me (%i/%i) \t Psi (%i/%i)\t PsiT (%i/%i)\n", P->Par.me, P->Par.me_comm_ST, P->Par.Max_me_comm_ST, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi, P->Par.me_comm_ST_PsiT, P->Par.Max_me_comm_ST_PsiT); k = -1; // k is global PsiNo counter for the desired group for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) { // go through all wave functions (plus the extra one for each process) OnePsiA = &Psi->AllPsiStatus[j]; // grab OnePsiA if (OnePsiA->PsiType == type) { // only take desired minimisation group k++; //fprintf(stderr,"(%i) ST_Psi: OnePsiA %i\tP->Par.me %i\n", P->Par.me,OnePsiA->my_color_comm_ST_Psi,P->Par.my_color_comm_ST_Psi); if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // Belongs to my Psi group? LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo]; else LOnePsiA = NULL; if (P->Par.me_comm_ST == 0) { // if we are process 0 of all, we may access file if (!OpenFileNo(P, &SrcPsiData, suffixdat, colorNo, "r", test ? 0 : P->Call.out[ReadOut])) { Error(SomeError,"ReadSrcPsiDensity: cannot open data file"); } for (i=P->Par.Max_me_comm_ST_Psi-1; i>=0;i--) { // load coefficients fseek( SrcPsiData, (long)((k*N[0]*N[1]*N[2]+i*((long)LevS->Dens->LocalSizeR))*sizeof(double)), SEEK_SET); // seek to beginning of this process' coefficients // readin if ((zahl = fread(destpsiR, sizeof(double), (size_t)(LevS->Dens->LocalSizeR), SrcPsiData)) != (size_t)LevS->Dens->LocalSizeR) { fclose(SrcPsiData); fprintf(stderr, "(%i)ReadSrcPsiDensity: only %i bytes read instead of expected %i\n", P->Par.me, zahl, LevS->Dens->LocalSizeR); Error(SomeError,"ReadSrcPsiDensity: fread Error"); } if (LOnePsiA == NULL) { // if it's not local, send away coefficients to correct PsiGroup (process 0 within that group) if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, OnePsiA->my_color_comm_ST_Psi, ReadSrcPsiTag, P->Par.comm_ST_PsiT) != MPI_SUCCESS) Error(SomeError, "ReadSrcPsiDensity: MPI_Send of loaded coefficients failed!"); //else //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i of loaded coefficients GlobalNo %i, owner %i succeeded!\n", P->Par.me, OnePsiA->my_color_comm_ST_Psi, k, i); } else { // if it's local ... if (i != 0) { // but share of array not for us, send to owner process within Psi group if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, i, ReadSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS) Error(SomeError, "ReadSrcPsiDensity: MPI_Send within Psi group of loaded coefficients failed!"); //else //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i within Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, i, k); } // otherwise it was our share already } } } else { if (LOnePsiA != NULL) { // receive if (P->Par.me_comm_ST_Psi == 0) { // if we are process 0 in the group, receive share from process 0 of all for (owner = P->Par.Max_me_comm_ST_Psi -1; owner >=0; owner--) { // for all processes of our Psi group if (MPI_Recv(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, ReadSrcPsiTag, P->Par.comm_ST_PsiT, &status) != MPI_SUCCESS) Error(SomeError, "ReadSrcPsiDensity: MPI_Recv of loaded coefficients failed!"); MPI_Get_count(&status, MPI_DOUBLE, &zahl); if (zahl != LevS->Dens->LocalSizeR) // check number of elements fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv from process 0 of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, owner, zahl, LevS->Dens->LocalSizeR); //else //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv from process 0 of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, owner); if (owner != 0) { // not our share of coefficients, send to owner in our Psi group if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, owner, ReadSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS) Error(SomeError, "ReadSrcPsiDensity: MPI_Send within Psi group of loaded coefficients failed!"); //else //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i within Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, owner, k); } else recv++; } // otherwise it's our share! } else { // our share within Psi Group not belonging to process 0 of all recv++; if (MPI_Recv(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, ReadSrcPsiTag, P->Par.comm_ST_Psi, &status) != MPI_SUCCESS) Error(SomeError, "ReadSrcPsiDensity: MPI_Recv of loaded coefficients failed!"); MPI_Get_count(&status, MPI_DOUBLE, &zahl); if (zahl != LevS->Dens->LocalSizeR) // check number of elements fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, P->Par.me_comm_ST_Psi, zahl, LevS->Dens->LocalSizeR); //else //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, P->Par.me_comm_ST_Psi); } } // otherwise we don't have anything to do with this } if (LOnePsiA != NULL) { SetArrayToDouble0((double *)srcpsiR, LevS->Dens->TotalSize*2); for (i=0; i < LevS->Dens->LocalSizeR; i++) // copy dest to src srcpsiR[i] = (fftw_real)destpsiR[i]; fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, srcpsiC, work); // fft transform //if (P->Call.out[PsiOut]) //fprintf(stderr,"(%i) LevSNo %i\t LocalPsi %p\n", P->Par.me, LevS->LevelNo, LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo]); for (i=0;iMaxG;i++) { // and copy into wave functions coefficients Index = LevS->GArray[i].Index; LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo][i].re = srcpsiC[Index].re/LevS->MaxN; LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo][i].im = srcpsiC[Index].im/LevS->MaxN; } } if ((P->Par.me_comm_ST == 0) && (SrcPsiData != NULL)) fclose(SrcPsiData); } } MPI_Allreduce (&recv, &Recv, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); fprintf(stderr,"(%i) Out of %i shares %i had to be received in total, %i from this process alone.\n", P->Par.me, P->Par.Max_me_comm_ST_Psi*Psi->NoOfPsis, Recv, recv); SpeedMeasure(P,ReadnWriteTime,StopTimeDo); } return 1; // everything went well till the end } /** Creates the density \ref suffixdensdoc and \ref suffixdensdx files for OpenDx. * Opens \ref suffixdensdoc, fills (pos&data file name, byte order, max mesh points, matrix alignment, steps) * and closes it. * Opens \ref suffixdensdx, then for every FileData::OutVisStep a describing structure for DX is written and * the file closed again. * \param *P Problem at hand * \param me my process number in communicator Psi (0 - do nothing, else - do) */ static void CreateDensityOutputGeneral(struct Problem *P, const int me) { FILE *DensityDoc, *DensityDx; char *posname, *datname; struct LatticeLevel *Lev = &P->Lat.Lev[STANDARTLEVEL]; unsigned int i, MaxPoints, N[NDIM]; double *RB = P->Lat.RealBasis; if (me) return; N[0] = Lev->N[0]*Lev->NUp[0]; N[1] = Lev->N[1]*Lev->NUp[1]; N[2] = Lev->N[2]*Lev->NUp[2]; MaxPoints = (N[0]+1)*(N[1]+1)*(N[2]+1); posname = (char*) Malloc(strlen(P->Files.mainname) + strlen(suffixdenspos) + 1,"OpenFile"); sprintf(posname, "%s%s", P->Files.mainname, suffixdenspos); datname = (char*) Malloc(strlen(P->Files.mainname) + strlen(suffixdensdat) + 1,"OpenFile"); sprintf(datname, "%s%s", P->Files.mainname, suffixdensdat); // write doc file OpenFile(P, &DensityDoc, suffixdensdoc, "w",P->Call.out[ReadOut]); fprintf(DensityDoc,"DensityPositions file = %s.####\n", posname); fprintf(DensityDoc,"DensityData file = %s.####\n", datname); fprintf(DensityDoc,"format = ieee float (Bytes %lu) %s\n",(unsigned long) sizeof(float),msb); fprintf(DensityDoc,"points = %i\n", MaxPoints); fprintf(DensityDoc,"majority = row\n"); fprintf(DensityDoc,"TimeSeries = %i\n",P->Files.OutVisStep+1); fclose(DensityDoc); // write DX file OpenFile(P, &DensityDx, suffixdensdx, "w",P->Call.out[ReadOut]); for (i=0; i < (unsigned int)P->Files.OutVisStep+1; i++) { // for every OutVis step if (i==0) { fprintf(DensityDx,"object \"gridcon\" class gridconnections counts %i %i %i\n\n",(N[0]+1),(N[1]+1),(N[2]+1)); if (P->Files.OutputPosType[i] != active) fprintf(DensityDx, "object \"posdens\" class gridpositions counts %i %i %i\norigin 0.0 0.0 0.0\ndelta %f %f %f\ndelta %f %f %f\ndelta %f %f %f\n\n", (N[0]+1),(N[1]+1),(N[2]+1), (float)(RB[0]/N[0]),(float)(RB[1]/N[0]),(float)RB[2]/N[0], (float)(RB[3]/N[1]),(float)(RB[4]/N[1]),(float)RB[5]/N[1], (float)(RB[6]/N[2]),(float)(RB[7]/N[2]),(float)RB[8]/N[2]); } if (P->Files.OutputPosType[i] == active) { fprintf(DensityDx, "object \"pos.%04u\" class array type float rank 1 shape 3 items %i %s binary\n",i,MaxPoints,msb); fprintf(DensityDx,"data file %s.%04u,0\n",posname,i); } fprintf(DensityDx,"attribute \"dep\" string \"positions\"\n"); fprintf(DensityDx,"# %lu - %lu Bytes\n\n",MaxPoints*i*(unsigned long)sizeof(float)*NDIM,MaxPoints*(i+1)*(unsigned long)sizeof(float)*NDIM-1); fprintf(DensityDx,"object \"dat.%04u\" class array type float rank 0 items %i %s binary\n",i,MaxPoints,msb); fprintf(DensityDx,"data file %s.%04u,0\n",datname,i); fprintf(DensityDx,"attribute \"dep\" string \"positions\"\n"); fprintf(DensityDx,"# %lu - %lu Bytes\n\n",MaxPoints*i*(unsigned long)sizeof(float),MaxPoints*(i+1)*(unsigned long)sizeof(float)-1); fprintf(DensityDx,"object \"obj.%04u\" class field\n",i); if (P->Files.OutputPosType[i] == active) fprintf(DensityDx,"component \"positions\" \"pos.%04i\"\n",i); if (P->Files.OutputPosType[i] != active) fprintf(DensityDx,"component \"positions\" \"posdens\"\n"); fprintf(DensityDx,"component \"connections\" \"gridcon\"\n"); fprintf(DensityDx,"component \"data\" \"dat.%04i\"\n",i); } fprintf(DensityDx,"\nobject \"series\" class series\n"); for (i=0; i < (unsigned int)P->Files.OutVisStep+1; i++) fprintf(DensityDx,"member %i \"obj.%04u\" position %f\n",i,i,(float)i); fprintf(DensityDx,"end\n"); fclose(DensityDx); Free(posname, "CreateDensityOutputGeneral: posname"); Free(datname, "CreateDensityOutputGeneral: datname"); } /** Calculates the OutVis density of the RiemannTensor level. * The usual pattern arises when a density is fftransformed: * -# over all grid vectors up to MaxG * -# over all doubly grid vectors up to MaxDoubleG * -# call to fft_3d_complex_to_real() * * In this case here followed by call to OutVisPosRTransformPosNFRto0() and finally FileData::work * is copied to FileData::PosR. * \param *Lat Lattice structure, containing Lattice::plan and LatticeLevel * \param *F FileData structure, containing FileData::PosC, FileData::PosR, FileData::work, FileData::Totalsize, FileData::LocalSizeR */ static void CalculateOutVisDensityPos(struct Lattice *Lat, struct FileData *F/*, const double FactorC_R, const double FactorR_C*/) { struct fft_plan_3d *plan = Lat->plan; struct RiemannTensor *RT = &Lat->RT; struct LatticeLevel *LevR = RT->LevR; fftw_complex *destC = F->PosC; fftw_real *destR = F->PosR; fftw_complex *work = F->work; fftw_real *workR = (fftw_real*)work; fftw_complex *PosFactor = F->PosFactor; fftw_complex *posfac, *destpos, *destRCS, *destRCD; fftw_complex *coeff; fftw_complex source; int i, Index, pos, n; int NF = NDIM, MaxNUp = F->MaxNUp, TotalSize = F->TotalSize, LocalSizeR = F->LocalSizeR; SetArrayToDouble0((double *)destC, TotalSize*2); for (i=0;i < LevR->MaxG;i++) { Index = LevR->GArray[i].Index; posfac = &PosFactor[MaxNUp*i]; destpos = &destC[MaxNUp*Index*NF]; coeff = &RT->Coeff[i]; for (pos=0; pos < MaxNUp; pos++) { for (n=0; n < NF; n++) { source.re = coeff[n].re; source.im = coeff[n].im; destpos[n+NF*pos].re = source.re*posfac[pos].re-source.im*posfac[pos].im; destpos[n+NF*pos].im = source.re*posfac[pos].im+source.im*posfac[pos].re; } } } for (i=0; i < LevR->MaxDoubleG; i++) { destRCS = &destC[LevR->DoubleG[2*i]*MaxNUp*NF]; destRCD = &destC[LevR->DoubleG[2*i+1]*MaxNUp*NF]; for (pos=0; pos < MaxNUp; pos++) { for (n=0; n < NF; n++) { destRCD[n+NF*pos].re = destRCS[n+NF*pos].re; destRCD[n+NF*pos].im = -destRCS[n+NF*pos].im; } } } fft_3d_complex_to_real(plan, LevR->LevelNo, FFTNFRVecUp0, destC, work); OutVisPosRTransformPosNFRto0(RT, destR, workR, NF); memcpy(destR,workR,sizeof(fftw_real)*LocalSizeR); } /** Prepare Density::DensityArray for output. * Into FileData::work subsequently each node (all z, all y, all x) is written as \f$\log(1+x)\f$, * where x is Density::DensityArray[TotalDensity]. In the end result is send to process 0 (yet not * received here, see CombineOutVisArray()). In case of RiemannTensor use, some more complex calculations * are made: FileData::PosR is used, the coefficient offset'ed to the current node and the log taken there. * \param *P Problem at hand * \param myPE this ranks process in the Psi communcator ParallelSimulationData::me_comm_ST_Psi * \param *srcdens Pointer to DensityArray which is to be displayed */ static void OutputOutVisDensity(struct Problem *P, const int myPE, fftw_real *srcdens) { int N[NDIM], n[NDIM], pos[NDIM]; int destpos = 0; double fac[NDIM], posd[NDIM]; float posf[NDIM+1]; struct Lattice *Lat = &P->Lat; struct LatticeLevel *Lev0 = &Lat->Lev[0]; fftw_real *srcpos = P->Files.PosR; //fftw_real *srcdens = Lev0->Dens->DensityArray[ActualDensity]; //[TotalDensity]; trick to display single density float *dest = (float *)P->Files.work; int Nx = Lev0->Plan0.plan->N[0]; int i; double min, max; N[0] = Lev0->Plan0.plan->local_nx; N[1] = Lev0->Plan0.plan->N[1]; N[2] = Lev0->Plan0.plan->N[2]; max = min = srcdens[0]; for (i=1;iR.Lev0->Dens->LocalSizeR;i++) { if (srcdens[i] < min) min = srcdens[i]; if (srcdens[i] > max) max = srcdens[i]; } if (P->Call.out[PsiOut]) fprintf(stderr,"(%i)OutputOutVisDensity: min %e\tmax %e\n",P->Par.me, min, max); // go through all nodes for (n[0]=0; n[0] < N[0]; n[0]++) { pos[0] = (n[0] == N[0] ? 0 : n[0]); for (n[1]=0; n[1] <= N[1]; n[1]++) { pos[1] = (n[1] == N[1] ? 0 : n[1]); for (n[2]=0; n[2] <= N[2]; n[2]++) { pos[2] = (n[2] == N[2] ? 0 : n[2]); // depending on RiemannTensor use, fill FileData::work switch (Lat->RT.ActualUse) { case inactive: case standby: if ((srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]) > 0.) dest[destpos] = log(1.0+(srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])])); else dest[destpos] = 0.; destpos++; break; case active: posf[0] = srcpos[0+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))]; posf[1] = srcpos[1+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))]; posf[2] = srcpos[2+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))]; fac[0] = ((n[0]+N[0]*myPE)/(double)Nx); fac[1] = (n[1]/(double)N[1]); fac[2] = (n[2]/(double)N[2]); RMat33Vec3(posd, Lat->RealBasis, fac); posf[0] += posd[0]; posf[1] += posd[1]; posf[2] += posd[2]; if ((srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]) > 0.) posf[3] = log(1.0+(srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])])); else posf[3] = 0.; dest[destpos+0] = posf[0]; dest[destpos+1] = posf[1]; dest[destpos+2] = posf[2]; dest[destpos+3] = posf[3]; destpos += 4; break; } } } } if (myPE) MPI_Send(dest, destpos, MPI_FLOAT, 0, OutputDensTag, P->Par.comm_ST_Psi); } /** Combines prepared electronic Psis density and output to file. * If we are process 0, open file suffixdensdat (only when RiemannTensor is used) and suffixdenspos, receive * FileData::work logarithmic coefficients sent by the other processes in OutputOutVisDensity(), go through all * nodes and save the coefficient to file - again depending on RiemannTensor use - followed by FileData::PosTemp * (for y and z nodes), close file(s). * \param *P Problem at hand * \param me this ranks process in the Psi communcator ParallelSimulationData::me_comm_ST_Psi * \param Maxme number of processes in this Psi communcator ParallelSimulationData::Max_me_comm_ST_Psi */ static void CombineOutVisDensity(struct Problem *P, const int me, const int Maxme) { int i,n[NDIM], N[NDIM]; float posf[NDIM+1]; float *source = (float *)P->Files.work; double posd[NDIM], fac[NDIM]; float *Temp = (float *)P->Files.PosTemp; struct Lattice *Lat = &P->Lat; struct LatticeLevel *Lev0 = &Lat->Lev[0]; float step = P->Files.OutVisStep; int No=0, destpos; FILE *DensityData, *DensityPos; int Nx = Lev0->Plan0.plan->N[0]+1; MPI_Status status; if (me) return; // if we are process 0! N[0] = Lev0->Plan0.plan->local_nx; N[1] = Lev0->Plan0.plan->N[1]+1; N[2] = Lev0->Plan0.plan->N[2]+1; // Open respective file depending on RiemannTensor use switch (Lat->RT.ActualUse) { case active: OpenFileNo(P, &DensityPos, suffixdenspos, (int)step, "wb",P->Call.out[ReadOut]); case inactive: case standby: OpenFileNo(P, &DensityData, suffixdensdat, (int)step, "wb",P->Call.out[ReadOut]); break; } // for all processes in the communicator for (i=0; i< Maxme; i++) { if (i) { // if process != 0, receive from this process /* MPI_Probe( i, OutputDensTag, P->Par.comm_ST_Psi, &status );*/ switch (Lat->RT.ActualUse) { case inactive: case standby: MPI_Recv( source, N[0]*N[1]*N[2], MPI_FLOAT, i, OutputDensTag, P->Par.comm_ST_Psi, &status ); break; case active: MPI_Recv( source, N[0]*N[1]*N[2]*4, MPI_FLOAT, i, OutputDensTag, P->Par.comm_ST_Psi, &status ); break; } } destpos = 0; // go through all nodes and save the coefficient to file DensityData, depending on RiemannTensor for (n[0]=0; n[0] < N[0]; n[0]++) { for (n[1]=0; n[1] < N[1]; n[1]++) { for (n[2]=0; n[2] < N[2]; n[2]++) { switch (Lat->RT.ActualUse) { case inactive: case standby: posf[3] = source[destpos]; destpos++; (void)fwrite(&posf[3], sizeof(float), (size_t)(1), DensityData); No++; if (i==0 && n[0] == 0) Temp[(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[3]; break; case active: posf[0] = source[destpos+0]; posf[1] = source[destpos+1]; posf[2] = source[destpos+2]; posf[3] = source[destpos+3]; destpos += 4; (void)fwrite(posf, sizeof(float), (size_t)(NDIM), DensityPos); (void)fwrite(&posf[3], sizeof(float), (size_t)(1), DensityData); No++; if (i==0 && n[0] == 0) { fac[0] = ((n[0]+N[0]*i)/(double)(Nx-1)); fac[1] = (n[1]/(double)(N[1]-1)); fac[2] = (n[2]/(double)(N[2]-1)); RMat33Vec3(posd, Lat->RealBasis, fac); posf[0] -= posd[0]; posf[1] -= posd[1]; posf[2] -= posd[2]; fac[0] = ((Nx-1)/(double)(Nx-1)); fac[1] = (n[1]/(double)(N[1]-1)); fac[2] = (n[2]/(double)(N[2]-1)); RMat33Vec3(posd, Lat->RealBasis, fac); posf[0] += posd[0]; posf[1] += posd[1]; posf[2] += posd[2]; Temp[0+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[0]; Temp[1+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[1]; Temp[2+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[2]; Temp[3+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[3]; } break; } } } } } n[0] = N[0]; for (n[1]=0; n[1] < N[1]; n[1]++) { for (n[2]=0; n[2] < N[2]; n[2]++) { switch (Lat->RT.ActualUse) { case inactive: case standby: (void)fwrite(&Temp[n[2]+N[2]*(n[1])], sizeof(float), (size_t)(1), DensityData); No++; break; case active: (void)fwrite(&Temp[(NDIM+1)*(n[2]+N[2]*(n[1]))], sizeof(float), (size_t)(NDIM), DensityPos); (void)fwrite(&Temp[(NDIM+1)*(n[2]+N[2]*(n[1]))+3], sizeof(float), (size_t)(1), DensityData); No++; break; } } } if (No != Nx*N[1]*N[2]) Error(SomeError,"CombineOutVisDensity: No != points"); switch (Lat->RT.ActualUse) { case active: fclose(DensityPos); case inactive: case standby: fclose(DensityData); break; } } /** Main output electronic Psis density for OpenDX. * If FileData::MeOutVis is set, calls OutputOutVisDensity() followed by CombineOutVisDensity(). * Beforehand CalculateOutVisDensityPos() is called if RiemannTensor is used. * \param *P Problem at hand * \param *src_dens Pointer to DensityArray which is to be displayed */ static void OutVisDensity(struct Problem *P, fftw_real *src_dens) { if (!P->Files.MeOutVis) return; if (P->Lat.RT.ActualUse == active) CalculateOutVisDensityPos(&P->Lat, &P->Files/*, P->Lat.FactorDensityR, P->Lat.FactorDensityC*/); OutputOutVisDensity(P, P->Par.me_comm_ST_Psi, src_dens); /* Achtung hier: P->Files.work (RT->TempC, Dens->DensityCArray[TempDensity]) fuer myPE == 0 nicht veraendern !!! */ CombineOutVisDensity(P, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi); } /** Opening and Initializing of output measurement files. * If this is process 0, opens and writes top line of FileData::ForcesFile, FileData::EnergyFile. * and sets FileData::MeOutVis and FileData::MeOutMes (if output desired) to 1, otherwise 0. * \param *P Problem at hand */ void InitOutputFiles(struct Problem *P) { struct FileData *F = &P->Files; F->ForcesFile = NULL; F->EnergyFile = NULL; F->HamiltonianFile = NULL; F->MinimisationFile = NULL; F->SpreadFile = NULL; // process 0 ? F->MeOutVis = ((P->Par.my_color_comm_ST == 0 && P->Par.my_color_comm_ST_Psi == 0 && F->DoOutVis) ? 1 : 0); F->MeOutCurr = ((P->Par.my_color_comm_ST == 0 && P->Par.my_color_comm_ST_Psi == 0 && F->DoOutCurr) ? 1 : 0); F->MeOutMes = ((P->Par.me == 0 && F->DoOutMes) ? 1 : 0); OpenFile(P, &F->HamiltonianFile, suffixhamiltonianall, "w",P->Call.out[ReadOut]); if (!F->MeOutMes) return; OpenFile(P, &F->ForcesFile, suffixforcesall, "w",P->Call.out[ReadOut]); if (F->ForcesFile == NULL) fprintf(stderr,"Error opening ForcesFile\n"); // write header of forces file fprintf(F->ForcesFile, "%s\t%s\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\n", "Type", "No", "Pos0", "Pos1", "Pos2", "Total0", "Total1", "Total2", "Local0", "Local1", "Local2", "NLocal0", "NLocal1", "NLocal2", "Ewald0", "Ewald1", "Ewald2"); OpenFile(P, &F->EnergyFile, suffixenergyall, "w",P->Call.out[ReadOut]); if (F->EnergyFile == NULL) fprintf(stderr,"Error opening EnergyFile\n"); // write header of energy file if (P->R.DoUnOccupied) { fprintf(F->EnergyFile, "%s\t\t%s\t\t%s\t%s\t\t%s\t%s\t\t%s\t%s\t%s\t\t%s\t\t%s\t%s\t\t%s\t\t%s\t\t%s\n", "Time", "Total", "Total+Gap", "Kinetic", "NonLocal", "GapPsi", "Correlation", "Exchange", "Pseudo", "Hartree", "GapDensity", "-Gauss", "Ewald", "IonKin", "ETotal"); } else { fprintf(F->EnergyFile, "%s\t\t%s\t\t%s\t\t%s\t%s\t%s\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\n", "Time", "Total", "Kinetic", "NonLocal", "Correlation", "Exchange", "Pseudo", "Hartree", "-Gauss", "Ewald", "IonKin", "ETotal"); } OpenFile(P, &F->MinimisationFile, suffixminall, "w",P->Call.out[ReadOut]); if (F->MinimisationFile == NULL) fprintf(stderr,"Error opening MinimsationFile\n"); // write header of minimsation file fprintf(F->MinimisationFile, "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n","Step", "Psi", "PsiStep", "TE", "ATE", "delta", "dEdt0", "ddEdtt0"); } /** Closes all output measurement files. * Closes FileData::ForcesFile and FileData::EnergyFile * \param *P Problem at hand */ void CloseOutputFiles(struct Problem *P) { struct FileData *F = &P->Files; if (!F->MeOutMes) return; if (F->ForcesFile != NULL) fclose(F->ForcesFile); // only they've been opened (thus not pointing to NULL) if (F->EnergyFile != NULL) fclose(F->EnergyFile); if (F->HamiltonianFile != NULL) fclose(F->HamiltonianFile); if (F->MinimisationFile != NULL) fclose(F->MinimisationFile); if (F->SpreadFile != NULL) fclose(F->SpreadFile); } /** Initialization of Problem::FileData structure for visual output. * If this is process 0 (and OutVis desired), allocate memory for FileData::PosTemp, FileData::work, * set the entries of FileData all to their corresponding values from RiemannTensor, * FileData::OutVisStep to zero. * \param *P Problem at hand */ void InitOutVisArray(struct Problem *P) { struct FileData *F = &P->Files; struct Lattice *Lat = &P->Lat; struct RiemannTensor *RT = &Lat->RT; struct LatticeLevel *Lev0 = &Lat->Lev[0]; F->OutputPosType = NULL; if (!F->MeOutVis) return; F->OutVisStep = 0; F->PosTemp = (fftw_complex *) Malloc(sizeof(float)*(Lev0->Plan0.plan->N[1]+1)*(Lev0->Plan0.plan->N[2]+1)* ((Lat->RT.Use != UseRT ? 0 : NDIM)+1), "InitOutVisArray: TempC"); F->work = (fftw_complex *)Lev0->Dens->DensityCArray[TempDensity]; if (Lat->RT.Use != UseRT) return; F->TotalSize = RT->TotalSize[RTAIRT]/NDIM; F->LocalSizeR = RT->LocalSizeR[RTAiGcg]; F->LocalSizeC = RT->LocalSizeC[RTAiGcg]; F->MaxNUp = RT->MaxNUp[RTPFRto0]; F->PosC = RT->DensityC[RTAiGcg]; F->PosR = (fftw_real *)F->PosC; F->work = RT->TempC; F->PosFactor = RT->PosFactor[RTPFRto0]; } static const char suffixionfor[] = ".ions.force"; //!< Suffix for ion forces file static const char suffixionZ[] = ".ions.datZ"; //!< Suffix for ion datZ file static const char suffixionpos[] = ".ions.pos"; //!< Suffix for ion positions file static const char suffixiondx[] = ".ions.dx"; //!< Suffix for ions dx file static const char suffixiondoc[] = ".ions.doc"; //!< Suffix for ions doc file static const char suffixsrciondoc[] = ".srcion.doc"; //!< Suffix for state ions doc file static const char suffixsrciondat[] = ".srcion.data"; //!< Suffix for state ions data file /** Output current ions state. * If this is process0, open file suffixsrciondoc for writing, output Ions::Max_Types and * Ions::Max_IonsOfType of each type - each in a new line - closes it. * Then opens suffixsrciondat for binary writing, outputs Lattice:RealBasis vectors and * position IonType::R and speed IonType::U, closes it. * \param *P Problem at hand * \note This is the ionic counterpart to the elecontric OutputSrcPsiDensity(), storing a so far made * calculation to file. */ static void OutSrcIons(struct Problem *P) { struct Ions *I = &P->Ion; double *U, *pos; double data[2*NDIM]; int is,ia,i; FILE *SrcIonDoc, *SrcIonData; if (!(P->Par.me == 0)) return; // output of ion types and numbers per type OpenFile(P, &SrcIonDoc, suffixsrciondoc, "w",P->Call.out[ReadOut]); fprintf(SrcIonDoc,"%i\n", I->Max_Types); for (is=0; is < I->Max_Types; is++) fprintf(SrcIonDoc,"%i\n", I->I[is].Max_IonsOfType); fclose(SrcIonDoc); OpenFile(P, &SrcIonData, suffixsrciondat, "wb",P->Call.out[ReadOut]); (void)fwrite(P->Lat.RealBasis, sizeof(double), (size_t)(NDIM_NDIM), SrcIonData); for (is=0; is < I->Max_Types; is++) { for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { U = &I->I[is].U[NDIM*ia]; pos = &I->I[is].R[NDIM*ia]; for (i=0;iIon; double *U, *pos; double data[2*NDIM]; int is,ia,i; int Max_Types; int *Max_IonsOfType = NULL; double RealBasis[NDIM_NDIM]; FILE *SrcIonDoc, *SrcIonData; // read the doc file and check if (OpenFile(P, &SrcIonDoc, suffixsrciondoc, "r",P->Call.out[ReadOut])) { if (fscanf(SrcIonDoc,"%i", &Max_Types) != 1) //Error(SomeError, "ReadSrcIons: read error"); return 0; if (Max_Types != I->Max_Types) //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, MaxTypes"); return 0; Max_IonsOfType = (int *) Malloc(Max_Types*sizeof(int), "ReadSrcIons: Max_IonsOfType"); for (is=0; is < Max_Types; is++) { if (fscanf(SrcIonDoc,"%i", &Max_IonsOfType[is]) != 1) //Error(SomeError, "ReadSrcIons: read error"); return 0; if (Max_IonsOfType[is] != I->I[is].Max_IonsOfType) //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, Max_IonsOfType"); return 0; } fclose(SrcIonDoc); // read basis, then positions and speeds of ions if (OpenFile(P, &SrcIonData, suffixsrciondat, "rb",P->Call.out[ReadOut])) { if (fread(RealBasis, sizeof(double), (size_t)(NDIM_NDIM), SrcIonData) != NDIM_NDIM) //Error(SomeError, "ReadSrcIons: read error"); return 0; for (i=0; i < NDIM_NDIM; i++) if (RealBasis[i] != P->Lat.RealBasis[i]) //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, RealBasis"); return 0; for (is=0; is < I->Max_Types; is++) { for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { if (fread(&data, sizeof(double), (size_t)(2*NDIM), SrcIonData) != 2*NDIM) //Error(SomeError, "ReadSrcIons: read error"); return 0; U = &I->I[is].U[NDIM*ia]; pos = &I->I[is].R[NDIM*ia]; for (i=0;iIon; struct FileData *F = &P->Files; int i,is,ia; double *fion, *pos; float data[6]; // holds temporarily twice NDIM values as write buffer int Z; char *datnamef, *datnameZ, *posname; FILE *IonsDataF, *IonsDataZ, *IonsPos, *IonsDoc, *IonsDx; if (!P->Files.MeOutVis && P->Par.me == 0) return; // generate file names datnamef = (char*) malloc(strlen(P->Files.mainname)+strlen(suffixionfor) + 1); sprintf(datnamef, "%s%s", P->Files.mainname, suffixionfor); datnameZ = (char*) malloc(strlen(P->Files.mainname)+strlen(suffixionZ) + 1); sprintf(datnameZ, "%s%s", P->Files.mainname, suffixionZ); posname = (char*) malloc(strlen(P->Files.mainname)+strlen(suffixionpos) + 1); sprintf(posname, "%s%s", P->Files.mainname, suffixionpos); // open, fill and close doc file if (OpenFile(P, &IonsDoc, suffixiondoc, "w",P->Call.out[ReadOut])) { fprintf(IonsDoc,"IonsPos file = %s.####\n", posname); fprintf(IonsDoc,"IonsForce file = %s.####\n", datnamef); fprintf(IonsDoc,"format = ieee float (Bytes %lu) %s = Force(3)\n",(unsigned long) sizeof(float),msb); fprintf(IonsDoc,"IonsZ file = %s.####\n", datnameZ); fprintf(IonsDoc,"format = int (Bytes %lu) %s = Z(1)\n",(unsigned long) sizeof(int),msb); fprintf(IonsDoc,"points = %i\n", I->Max_TotalIons); fprintf(IonsDoc,"majority = row\n"); fprintf(IonsDoc,"TimeSeries = %i\n",F->OutVisStep+1); fclose(IonsDoc); } // open dx file and fill it with each output step, close it if (OpenFile(P, &IonsDx, suffixiondx, "w",P->Call.out[ReadOut])) { for (i=0; i < F->OutVisStep+1; i++) { if (i==0) { /* fprintf(IonsDx,"object \"ioncon\" class array type int rank 1 shape 2 items 0 data follows\nattribute \"element type\" string \"lines\"\nattribute \"ref\" string \"positions\"\n\n"); */ fprintf(IonsDx,"object \"iondatZ\" class array type int rank 0 items %i %s binary\n",I->Max_TotalIons,msb); fprintf(IonsDx,"data file %s,0\n",datnameZ); fprintf(IonsDx,"attribute \"dep\" string \"positions\"\n\n"); } fprintf(IonsDx,"object \"ionpos.%04i\" class array type float rank 1 shape 3 items %i %s binary\n",i,I->Max_TotalIons,msb); fprintf(IonsDx,"data file %s.%04i,0\n\n",posname,i); fprintf(IonsDx,"object \"iondatF.%04i\" class array type float rank 1 shape 3 items %i %s binary\n",i,I->Max_TotalIons,msb); fprintf(IonsDx,"data file %s.%04i,0\n",datnamef,i); fprintf(IonsDx,"attribute \"dep\" string \"positions\"\n\n"); fprintf(IonsDx,"object \"ionobjF.%04i\" class field\n",i); fprintf(IonsDx,"component \"positions\" \"ionpos.%04i\"\n",i); fprintf(IonsDx,"component \"data\" \"iondatF.%04i\"\n",i); /*fprintf(IonsDx,"component \"connections\" \"ioncon\"\n\n");*/ fprintf(IonsDx,"object \"ionobjZ.%04i\" class field\n",i); fprintf(IonsDx,"component \"positions\" \"ionpos.%04i\"\n",i); fprintf(IonsDx,"component \"data\" \"iondatZ\"\n"); /* fprintf(IonsDx,"component \"connections\" \"ioncon\"\n\n");*/ } fprintf(IonsDx,"\nobject \"ionseriesF\" class series\n"); for (i=0; i < F->OutVisStep+1; i++) fprintf(IonsDx,"member %i \"ionobjF.%04i\" position %f\n",i,i,(float)i); fprintf(IonsDx,"\nobject \"ionseriesZ\" class series\n"); for (i=0; i < F->OutVisStep+1; i++) fprintf(IonsDx,"member %i \"ionobjZ.%04i\" position %f\n",i,i,(float)i); fprintf(IonsDx,"end\n"); fclose(IonsDx); } free(datnamef); free(datnameZ); free(posname); // open IonForces, IonZ and IonPosition file, write forces respectively positions for each ion of each type, close them if (OpenFileNo(P, &IonsDataF, suffixionfor, F->OutVisStep, "wb",P->Call.out[ReadOut])) { if (F->OutVisStep == 0) OpenFile(P, &IonsDataZ, suffixionZ, "wb",P->Call.out[ReadOut]); if (OpenFileNo(P, &IonsPos, suffixionpos, F->OutVisStep, "wb",P->Call.out[ReadOut])) { for (is=0; is < I->Max_Types; is++) { for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { fion = &I->I[is].FIon[NDIM*ia]; pos = &I->I[is].R[NDIM*ia]; for (i=0;i<3;i++) { data[i+3] = fion[i]; data[i] = pos[i]; } Z = I->I[is].Z; if (fwrite(&data[0],sizeof(float), (size_t)(3),IonsPos) != 3) Error(FileOpenParams, "Error writing ion positions!"); if (F->OutVisStep == 0) (void)fwrite(&Z,sizeof(int), (size_t)(1),IonsDataZ); if (fwrite(&data[3],sizeof(float), (size_t)(3),IonsDataF) != 3) Error(FileOpenParams, "Error writing ionic forces!"); } } fclose(IonsPos); } if (F->OutVisStep == 0) fclose(IonsDataZ); fclose(IonsDataF); } } /** Output electronic density and ion state files with so far made calculations. * If CallOptions::WriteSrcFiles is set, OutputSrcPsiDensity() and OutSrcIons() are called. * \param *P Problem at hand * \param type which PsiTypeTag should be put to file */ void OutputVisSrcFiles(struct Problem *P, enum PsiTypeTag type) { if (P->Call.WriteSrcFiles) { if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Write srcpsi to disk\n", P->Par.me); OutputSrcPsiDensity(P, type); if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Write srcion to disk\n", P->Par.me); OutSrcIons(P); } // if (!P->Files.MeOutVis) return; // P->Files.OutputPosType = // Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+1),"OutputVis"); // P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse; // CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi); // OutVisDensity(P); // OutVisIons(P); // if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep); // /* P->Files.OutVisStep++; Genau ebend nicht hochzaehlen - wird immer ueberschrieben */ } /** Main output total electronic density and ion data for OpenDX. * Calls subsequently preparing CreateDensityOutputGeneral(), then output of electronic * densities OutVisDensity() and ion data OutVisIons(), increasing finally FileData::OutVisStep. * \param *P Problem at hand * \param *srcdens Pointer to DensityArray which is to be displayed * \note Output is made only RunStruct::OutVisStep steps and if FileData::MeOutVis is set. */ void OutputVis(struct Problem *P, fftw_real *srcdens) { if (!P->Files.MeOutVis) return; P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+1),"OutputVis"); P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse; CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi); OutVisDensity(P, srcdens); OutVisIons(P); if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep); P->Files.OutVisStep++; } /** Output of each current density component for OpenDX. * \param *P Problem at hand * \note Output is made only RunStruct::OutVisStep steps and if FileData::MeOutVis is set. */ void OutputCurrentDensity(struct Problem *P) { int index, i, r; fftw_real *density = P->R.Lev0->Dens->DensityArray[ActualDensity]; fftw_real *CurrentDensity[NDIM*NDIM]; if (!P->Files.MeOutCurr) return; P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+(1)*NDIM),"OutputVis"); for (i=0;i<(1)*NDIM;i++) P->Files.OutputPosType[P->Files.OutVisStep+i] = P->Lat.RT.ActualUse; fprintf(stderr,"(%i) OutVisStep %i, OutputPosType %p\n",P->Par.me, P->Files.OutVisStep, P->Files.OutputPosType); // due to preprocessor values we can't put the following stuff into a loop CurrentDensity[0] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity0]; CurrentDensity[1] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity1]; CurrentDensity[2] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity2]; CurrentDensity[3] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity3]; CurrentDensity[4] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity4]; CurrentDensity[5] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity5]; CurrentDensity[6] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity6]; CurrentDensity[7] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity7]; CurrentDensity[8] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity8]; // output current density, not vector component for (index=0;indexR.Lev0->Dens->TotalSize*2); // reset for (r=0;rR.Lev0->Dens->LocalSizeR;r++) { for (i=0;iPar.me_comm_ST_Psi); OutVisDensity(P, density); OutVisIons(P); if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep); P->Files.OutVisStep++; } } /** Output each orbital in a certain step order for OpenDX. * \param *P Problem at hand * \param offset from which step do we start * \param increment by which increment do we advance step-wise * \param type Only PsiTypeTag orbitals are displayed */ void OutputVisAllOrbital(struct Problem *P, int offset, int increment, enum PsiTypeTag type) { struct Lattice *Lat = &P->Lat; struct Psis *Psi = &Lat->Psi; struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; struct LatticeLevel *Lev0 = R->Lev0; struct Density *Dens0 = Lev0->Dens; struct OnePsiElement *OnePsiA, *LOnePsiA; MPI_Status status; int ElementSize = (sizeof(fftw_complex) / sizeof(double)); fftw_complex *LPsiDatA; int i, p, RecvSource; int Num = Psi->NoOfPsis; if (!P->Files.MeOutVis) return; P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+(Num)),"OutputVis"); P->Files.OutVisStep += offset; P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse; for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // go through all wave functions (plus the extra one for each process) OnePsiA = &Psi->AllPsiStatus[i]; // grab the desired OnePsiA if (OnePsiA->PsiType == type) { // drop if extra one if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local? LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo]; else LOnePsiA = NULL; if (LOnePsiA == NULL) { // if it's not local ... receive it from respective process into TempPsi RecvSource = OnePsiA->my_color_comm_ST_Psi; MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, WannierTag1, P->Par.comm_ST_PsiT, &status ); LPsiDatA=LevS->LPsi->TempPsi; } else { // .. otherwise send it to all other processes (Max_me... - 1) for (p=0;pPar.Max_me_comm_ST_PsiT;p++) if (p != OnePsiA->my_color_comm_ST_Psi) MPI_Send( LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, p, WannierTag1, P->Par.comm_ST_PsiT); LPsiDatA=LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo]; } // LPsiDatA is now set to the coefficients of OnePsi either stored or MPI_Received if (P->Files.MeOutVis) { P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse; // recalculate density for the specific wave function ... CalculateOneDensityR(Lat, LevS, Dens0, LPsiDatA, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0); // ... and output (wherein ActualDensity is used instead of TotalDensity) //OutputVis(P); CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi); OutVisDensity(P, Dens0->DensityArray[ActualDensity]); OutVisIons(P); if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep); P->Files.OutVisStep+=increment; P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse; } } } } /** Read source files containing stored calculations. * Calls ReadSrcPsiDensity() and ReadSrcIons(). * \param *P Problem at hand */ void ReadSrcFiles(struct Problem *P) { if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadSrcPsiDensity\n", P->Par.me); ReadSrcPsiDensity(P, Occupied, 0, P->R.LevSNo); ReadSrcPsiDensity(P, UnOccupied, 0, P->R.LevSNo); if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadSrcIons\n", P->Par.me); ReadSrcIons(P); } /** Plots a cut plane of the real density of one wave function. * \param *P Problem at hand * \param index index of axis (vector orthogonal to plane) * \param node node specifying where to cut at the given axis * \param wavenr global number of wave function * \param *density density array to plot * \sa PlotVectorPlane() - very similar */ void PlotSrcPlane(struct Problem *P, int index, double node, int wavenr, fftw_real *density) { struct RunStruct *R = &P->R; struct Lattice *Lat = &P->Lat; struct LatticeLevel *Lev0 = R->Lev0; const int myPE = P->Par.me_comm_ST_Psi; char filename[255], spin[12]; char *suchpointer; FILE *PlotFile = NULL; time_t seconds; // open file time(&seconds); // get current time switch (Lat->Psi.PsiST) { case SpinDouble: sprintf(&filename[0], ".psi%i_cut%i.csv", wavenr, index); strncat(spin,"SpinDouble",12); break; case SpinUp: sprintf(&filename[0], ".psi%i_cut%i_up.csv", wavenr, index); strncat(spin,"SpinUp",12); break; case SpinDown: sprintf(&filename[0], ".psi%i_cut%i_down.csv", wavenr, index); strncat(spin,"SpinDown",12); break; } if (!myPE) { // only process 0 writes to file OpenFile(P, &PlotFile, filename, "w", P->Call.out[ReadOut]); strcpy(filename, ctime(&seconds)); suchpointer = strchr(filename, '\n'); if (suchpointer != NULL) *suchpointer = '\0'; if (PlotFile != NULL) { fprintf(PlotFile,"# Psi %i, type %s (real density) plot of plane perpendicular to direction e_%i at node %lg, seed %i, config %s, run on %s, #cpus %i", wavenr, spin, index, node, R->Seed, P->Files.default_path, filename, P->Par.Max_me_comm_ST_Psi); fprintf(PlotFile,"\n"); } else { Error(SomeError, "PlotSrcPlane: Opening Plot File"); } } // plot density PlotRealDensity(P, Lev0, PlotFile, index, node, density, density); if (PlotFile != NULL) { // close file fclose(PlotFile); } } /** plots a cut plane of a given 3d real density. * \param *P Problem at hand, contains pointer to Lattice structure * \param *Lev LatticeLevel of the real density * \param PlotFile file pointer (already open and valid) * \param index index of lattice axis * \param n_orth position on lattice axis where to cut * \param *density1 first real density array * \param *density2 second real density array (point to \a *density1 if not needed) */ void PlotRealDensity(struct Problem *P, struct LatticeLevel *Lev, FILE *PlotFile, int index, double n_orth, fftw_real *density1, fftw_real *density2) { struct Lattice *Lat = &P->Lat; int n[NDIM], n0; int N[NDIM]; N[0] = Lev->Plan0.plan->N[0]; N[1] = Lev->Plan0.plan->N[1]; N[2] = Lev->Plan0.plan->N[2]; const int N0 = Lev->Plan0.plan->local_nx; const int myPE = P->Par.me_comm_ST_Psi; double fac[NDIM], x[NDIM]; int i0, i = 0; int PE, zahl; double *buffer; int node = ceil(n_orth/Lat->RealBasisQ[index]*N[index]); MPI_Status status; int sizes[P->Par.Max_me_comm_ST_Psi], c0, c1; switch (index) { case 0: zahl = 4*N[1]*N[2]; break; case 1: zahl = 4*N0*N[2]; break; case 2: zahl = 4*N0*N[1]; break; } buffer = Malloc(sizeof(double)*zahl,"PlotRealDensity: buffer"); c0 = cross(index,0); c1 = cross(index,1); // then for every point on the grid in real space ... for (n0=0;n01 case if (n[index] == node) { // only on the correct plane orthogonal to desired axis and at desired node ... fac[0] = (double)n[0]/(double)N[0]; fac[1] = (double)n[1]/(double)N[1]; fac[2] = (double)n[2]/(double)N[2]; RMat33Vec3(x, Lat->RealBasis, fac); // relative coordinate times basis matrix gives absolute ones i0 = n[2]+N[2]*(n[1]+N[1]*n0); // index to local density array buffer[i++] = x[c0]; // fill buffer buffer[i++] = x[c1]; buffer[i++] = density1[i0]; buffer[i++] = density2[i0]; if (i > zahl) Error(SomeError, "PlotRealDensity: buffer too small!"); } } // exchange sizes of each buffer MPI_Allgather(&i, 1, MPI_INT, sizes, 1, MPI_INT, P->Par.comm_ST_Psi); if (myPE == 0) { for (PE=0; PE < P->Par.Max_me_comm_ST_Psi; PE++) { if (PE != 0) { // receive them if (MPI_Recv(buffer, sizes[PE], MPI_DOUBLE, PE, PlotRealDensityTag, P->Par.comm_ST_Psi, &status) != MPI_SUCCESS) Error(SomeError, "PlotRealDensity: MPI_Recv failure!"); MPI_Get_count(&status, MPI_DOUBLE, &zahl); if (zahl != sizes[PE]) Error(SomeError, "PlotRealDensity: received unexpected amount of elements!"); } //write them: local one (still in buffer) and received ones for (i0 = 0; i0 < sizes[PE];) { fprintf(PlotFile,"%e", buffer[(i0)++]); if ((i0 % 4) == 0) { fprintf(PlotFile,"\n"); } else { fprintf(PlotFile,"\t"); } } } } else { // send them if (MPI_Send(buffer, i, MPI_DOUBLE, 0, PlotRealDensityTag, P->Par.comm_ST_Psi) != MPI_SUCCESS) Error(SomeError, "PlotRealDensity: MPI_Send failure!"); } Free(buffer, "PlotRealDensity: buffer"); }