/** \file grad.c * Gradient Calculation. * * The most important structure is Gradient with its Gradient::GradientArray. * It contains the various stages of the calculated conjugate direction, beginning * from the Psis::LocalPsiStatus. Subsequently, CalculateCGGradient(), CalculatePreConGrad(), * CalculateConDir() and CalculateLineSearch() are called which use earlier stages to * calculate the latter ones: ActualGradient, PreConDirGradient, ConDirGradient and OldConDirGradient, * do the one-dimensional minimisation and finally update the coefficients of the wave function, * whilst these ought to remain orthogonal. * * Initialization of the Gradient structure and arrays is done in InitGradients(), deallocation * in RemoveGradients(). * * Small functions calculate the scalar product between two gradient directions GradSP(), find * the minimal angle in the line search CalculateDeltaI() or evaluate complex expressions * CalculateConDirHConDir(), down to basic calculation of the gradient of the wave function * CalculateGradientNoRT(). CalculateHamiltonian() sets up the Hamiltonian for the Kohn-Sham * orbitals and calculates the eigenvalues. * * Other functions update the coefficients after the position of the ion has changed - * CalcAlphaForUpdateWaveAfterIonMove() and UpdateWaveAfterIonMove(). * * Project: ParallelCarParrinello \author Jan Hamaekers \date 2000 File: grad.c $Id: grad.c,v 1.90 2007-03-29 13:38:04 foo Exp $ */ #include #include #include #include #include #include #include #include #include #include #include #include #include "mymath.h" #include "data.h" #include "errors.h" #include "energy.h" #include "helpers.h" #include "myfft.h" #include "density.h" #include "gramsch.h" #include "perturbed.h" #include "pseudo.h" #include "grad.h" #include "run.h" #include "wannier.h" /** Initialization of Gradient Calculation. * Gradient::GradientArray[GraSchGradient] is set onto "extra" Psis::LocalPsi (needed * for orthogonal conjugate gradient directions!), * for the remaining GradientTypes memory is allocated and zerod. * \param *P Problem at hand * \param *Grad Gradient structure */ void InitGradients(struct Problem *P, struct Gradient *Grad) { struct Lattice *Lat = &P->Lat; struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; struct LatticeLevel *ILevS = R->InitLevS; struct Psis *Psi = &Lat->Psi; int i; Grad->GradientArray[GraSchGradient] = LevS->LPsi->LocalPsi[Psi->LocalNo]; for (i=1; iGradientArray[i] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILevS->MaxG, "InitGradients"); SetArrayToDouble0((double*)Grad->GradientArray[i], 2*ILevS->MaxG); } } /** Deallocates memory used in Gradient::GradientArray. * Gradient::GradientArray[GraSchGradient] is set to null. * \param *P Problem at hand * \param *Grad Gradient structure */ void RemoveGradients(struct Problem *P, struct Gradient *Grad) { int i; P=P; Grad->GradientArray[GraSchGradient] = NULL; for (i=1; iGradientArray[i], "RemoveGradients: Grad->GradientArray[i]"); } /** Calculates the scalar product of two wave functions with equal symmetry. * \param *P Problem at hand * \param *Lev LatticeLevel structure * \param *LPsiDatA first array of complex wave functions coefficients * \param *LPsiDatB second array of complex wave functions coefficients * \return complex scalar product of the two wave functions. * \warning MPI_Allreduce has to be done still, remember coefficients are shared among processes! */ double GradSP(struct Problem *P, struct LatticeLevel *Lev, const fftw_complex *LPsiDatA, const fftw_complex *LPsiDatB) { double LocalSP=0.0; int i = 0; if (Lev->GArray[0].GSq == 0.0) { LocalSP += LPsiDatA[0].re * LPsiDatB[0].re; i++; } for (; i < Lev->MaxG; i++) { LocalSP += 2.*(LPsiDatA[i].re * LPsiDatB[i].re + LPsiDatA[i].im * LPsiDatB[i].im); // factor 2 because of the symmetry arising from gamma point! } return(LocalSP); } /** Calculates the scalar product of two wave functions with unequal symmetry. * \param *P Problem at hand * \param *Lev LatticeLevel structure * \param *LPsiDatA first array of complex wave functions coefficients * \param *LPsiDatB second array of complex wave functions coefficients * \return complex scalar product of the two wave functions. * \warning MPI_Allreduce has to be done still, remember coefficients are shared among processes! */ double GradImSP(struct Problem *P, struct LatticeLevel *Lev, fftw_complex *LPsiDatA, fftw_complex *LPsiDatB) { double LocalSP=0.0; int i = 0; if (Lev->GArray[0].GSq == 0.0) { i++; } for (; i < Lev->MaxG; i++) { LocalSP += 2.*(LPsiDatA[i].re * LPsiDatB[i].im - LPsiDatA[i].im * LPsiDatB[i].re); // factor 2 because of the symmetry arising from gamma point! } return(LocalSP); } /** Calculation of gradient of wave function(without RiemannTensor). * We want to calculate the variational derivative \f$\frac{\delta}{\delta \psi^*_j} E[{\psi^*_j}]\f$, * which we need for the minimisation of the energy functional. In the plane wave basis this evaluates to * \f[ * \langle \chi_G| -\frac{1}{2}\nabla^2 + \underbrace{V^H + V^{ps,loc} + V^{XC}}_{V^{loc}} + V^{ps,nl} | \psi_i \rangle * \qquad (section 4.3) * \f] * Therefore, DensityType#HGDensity is back-transformed, exchange correlation potential added up (CalculateXCPotentialNoRT()), * and vector multiplied with the wave function Psi (per node R), transformation into reciprocal grid space and we are back * in the plane wave basis set as desired. CalculateAddNLPot() adds the non-local potential. Finally, the kinetic term * is added up (simple in reciprocal space). The energy eigenvalue \a *Lambda is computed by summation over all reciprocal grid * vectors and in the meantime also the gradient direction (negative derivative) stored in \a *grad while the old values * are moved to \a *oldgrad. At last, \a *Lambda is summed up among all processes. * \f[ * \varphi_i^{(m)} (G) = \frac{1}{2} |G|^2 c_{i,G} + \underbrace{V^H (G) + V^{ps,loc} (G)}_{+V^{HG} (G)} + V^{ps,nl} (G) | \psi_i (G) \rangle * \qquad (\textnormal{end of section 4.5.2}) * \f] * \param *P Problem at hand, contains Lattice and RunStruct * \param *source source wave function coefficients * \param PsiFactor occupation number, see Psis::PsiFactor * \param *Lambda eigenvalue \f$\lambda_i = \langle \psi_i |h|\psi_i \rangle\f$, see OnePsiElementAddData::Lambda * \param ***fnl non-local form factors over ion type, ion of type and reciprocal grid vector, see PseudoPot::fnl * \param *grad return array for (negative) gradient coefficients \f$\varphi_i^{(m)} (G)\f$ * \param *oldgrad return array for (negative) gradient coefficients \f$\varphi_i^{(m-1)} (G)\f$ * \param *Hc_grad return array for complex coefficients of \f$H | \psi_i (G) \rangle\f$, see GradientArray#HcGradient * \sa CalculateConDirHConDir() - similar calculation for another explanation. */ void CalculateGradientNoRT(struct Problem *P, fftw_complex *source, const double PsiFactor, double *Lambda, fftw_complex ***fnl, fftw_complex *grad, fftw_complex *oldgrad, fftw_complex *Hc_grad) { struct Lattice *Lat = &P->Lat; struct RunStruct *R = &P->R; struct LatticeLevel *Lev0 = R->Lev0; struct LatticeLevel *LevS = R->LevS; struct Density *Dens0 = Lev0->Dens; struct fft_plan_3d *plan = Lat->plan; int Index, g; fftw_complex *work = Dens0->DensityCArray[TempDensity]; fftw_real *HGcR = Dens0->DensityArray[HGcDensity]; fftw_complex *HGcRC = (fftw_complex*)HGcR; fftw_complex *HGC = Dens0->DensityCArray[HGDensity]; fftw_real *HGCR = (fftw_real *)HGC; fftw_complex *PsiC = Dens0->DensityCArray[ActualPsiDensity]; fftw_real *PsiCR = (fftw_real *)PsiC; int nx,ny,nz,iS,i0,s; const int Nx = LevS->Plan0.plan->local_nx; const int Ny = LevS->Plan0.plan->N[1]; const int Nz = LevS->Plan0.plan->N[2]; const int NUpx = LevS->NUp[0]; const int NUpy = LevS->NUp[1]; const int NUpz = LevS->NUp[2]; double lambda; const double HGcRCFactor = 1./LevS->MaxN; SpeedMeasure(P, LocTime, StartTimeDo); // back-transform HGDensity //if (isnan(HGC[0].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGC[%i]_%i.re = NaN!\n", 0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); } fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF1, HGC, work); // evaluate exchange potential with this density, add up onto HGCR //if (isnan(HGCR[0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGcR[:%i:]_%i = NaN!\n", 0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); } CalculateXCPotentialNoRT(P, HGCR); // add V^{xc} or \epsilon^{xc} upon V^H + V^{ps} for (nx=0;nxActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); } //if (isnan(HGCR[i0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGCR[%i]_%i = NaN!\n", i0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); } } fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGcRC, work); SpeedMeasure(P, LocTime, StopTimeDo); /* NonLocalPP */ //SpeedMeasure(P, NonLocTime, StartTimeDo); CalculateAddNLPot(P, Hc_grad, fnl, PsiFactor); // also resets Hc_grad beforehand //SetArrayToDouble0((double *)Hc_grad,2*LevS->MaxG); SpeedMeasure(P, NonLocTime, StopTimeDo); /* Lambda und Grad */ for (g=0;gMaxG;g++) { Index = LevS->GArray[g].Index; /* FIXME - factoren */ Hc_grad[g].re += PsiFactor*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[g].GSq*source[g].re); Hc_grad[g].im += PsiFactor*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[g].GSq*source[g].im); //if (isnan(source[g].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): source_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); } } if (grad != NULL) { lambda = 0.0; s = 0; if (LevS->GArray[0].GSq == 0.0) { lambda += Hc_grad[0].re*source[0].re + Hc_grad[0].im*source[0].im; oldgrad[0].re = grad[0].re; // store away old preconditioned steepest decent gradient oldgrad[0].im = grad[0].im; grad[0].re = -Hc_grad[0].re; grad[0].im = -Hc_grad[0].im; s++; } for (g=s;gMaxG;g++) { lambda += 2.*(Hc_grad[g].re*source[g].re + Hc_grad[g].im*source[g].im); oldgrad[g].re = grad[g].re; // store away old preconditioned steepest decent gradient oldgrad[g].im = grad[g].im; grad[g].re = -Hc_grad[g].re; grad[g].im = -Hc_grad[g].im; //if (isnan(Hc_grad[g].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): Hc_grad_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); } } MPI_Allreduce ( &lambda, Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi); //fprintf(stderr,"(%i) Lambda = %e\n",P->Par.me, *Lambda); } } /** Calculates the Hamiltonian matrix in Kohn-Sham basis set and outputs to file. * In a similiar manner to TestGramSch() all necessary coefficients are retrieved from the * other processes and \f$\langle Psi_i | H | Psi_j \rangle\f$ = \f$\lambda_{i,j}\f$ calculated. * \param *P Problem at hand * \sa CalculateGradientNoRT() - same calculation (of \f$\lambda_i\f$), yet only for diagonal elements */ void CalculateHamiltonian(struct Problem *P) { struct Lattice *Lat = &P->Lat; struct FileData *F = &P->Files; struct RunStruct *R = &P->R; struct Psis *Psi = &Lat->Psi; struct LatticeLevel *Lev0 = R->Lev0; struct LatticeLevel *LevS = R->LevS; //struct LevelPsi *LPsi = LevS->LPsi; struct Density *Dens0 = Lev0->Dens; int Index; struct fft_plan_3d *plan = Lat->plan; fftw_complex *work = Dens0->DensityCArray[TempDensity]; fftw_real *HGcR = Dens0->DensityArray[HGcDensity]; fftw_complex *HGcRC = (fftw_complex*)HGcR; fftw_complex *HGC = Dens0->DensityCArray[HGDensity]; fftw_real *HGCR = (fftw_real *)HGC; fftw_complex *PsiC = Dens0->DensityCArray[ActualPsiDensity]; fftw_real *PsiCR = (fftw_real *)PsiC; fftw_complex ***fnl; fftw_complex *Hc_grad = P->Grad.GradientArray[HcGradient]; int nx,ny,nz,iS,i0,s; const int Nx = LevS->Plan0.plan->local_nx; const int Ny = LevS->Plan0.plan->N[1]; const int Nz = LevS->Plan0.plan->N[2]; const int NUpx = LevS->NUp[0]; const int NUpy = LevS->NUp[1]; const int NUpz = LevS->NUp[2]; double lambda, Lambda; const double HGcRCFactor = 1./LevS->MaxN; double pot_xc, Pot_xc; int PsiFactorA = 1; int PsiFactorB = 1; int NoPsis = 0; int NoOfPsis = 0; if (P->Call.out[LeaderOut]) fprintf(stderr, "(%i)Setting up Hamiltonian for Level %d...",P->Par.me, LevS->LevelNo); if (P->Call.out[StepLeaderOut]) fprintf(stderr, "\n(%i) H = \n",P->Par.me); ApplyTotalHamiltonian(P,LevS->LPsi->LocalPsi[0],Hc_grad, P->PP.fnl[0], 1., 1); // update HGDensity-> // go through all pairs of orbitals int i,k,l,m,j=0,RecvSource; MPI_Status status; struct OnePsiElement *OnePsiB, *OnePsiA, *LOnePsiA, *LOnePsiB; int ElementSize = (sizeof(fftw_complex) / sizeof(double)); fftw_complex *LPsiDatA=NULL, *LPsiDatB; l = -1; switch (Psi->PsiST) { case SpinDouble: NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoDouble]; NoOfPsis += Psi->GlobalNo[PsiMaxAdd]; break; case SpinUp: NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoUp]; NoOfPsis += Psi->GlobalNo[PsiMaxAdd]; break; case SpinDown: NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoDown]; NoOfPsis += Psi->GlobalNo[PsiMaxAdd]; break; }; gsl_matrix *H = gsl_matrix_calloc(NoOfPsis,NoOfPsis); for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // go through all wave functions //fprintf(stderr,"(%i) GlobalNo: %d\tLocalNo: %d\n", P->Par.me,Psi->AllPsiStatus[i].MyGlobalNo,Psi->AllPsiStatus[i].MyLocalNo); OnePsiA = &Psi->AllPsiStatus[i]; // grab OnePsiA if (OnePsiA->PsiType <= UnOccupied) { // drop the extra and perturbed ones l++; 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, HamiltonianTag, P->Par.comm_ST_PsiT, &status ); LPsiDatA=LevS->LPsi->TempPsi; } else { // .. otherwise send it to all other processes (Max_me... - 1) for (k=0;kPar.Max_me_comm_ST_PsiT;k++) if (k != OnePsiA->my_color_comm_ST_Psi) MPI_Send( LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, k, HamiltonianTag, 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 (LOnePsiA != NULL) { // if fnl (!) are locally available! // Trick is: fft wave function by moving it to higher level and there take only each ...th element into account PsiFactorA = 1; //Psi->LocalPsiStatus[OnePsiA->MyLocalNo].PsiFactor; CalculateOneDensityR(Lat, LevS, Dens0, LPsiDatA, Dens0->DensityArray[ActualDensity], R->FactorDensityR*PsiFactorA, 1); // note: factor is not used when storing result in DensityCArray[ActualPsiDensity] in CalculateOneDensityR()! for (nx=0;nxLevelNo, FFTNF1, HGcRC, work); // V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle /* NonLocalPP */ fnl = P->PP.fnl[OnePsiA->MyLocalNo]; // fnl only locally available! CalculateAddNLPot(P, Hc_grad, fnl, PsiFactorA); // sets Hc_grad to zero beforehand // V^{ps,nl} (G)| \Psi_b (G) + [V^{HG} (G) + V^{XC} (G)] | \Psi_b (G) \rangle /* Lambda und Grad */ s = 0; if (LevS->GArray[0].GSq == 0.0) { Index = LevS->GArray[0].Index; /* FIXME - factoren */ Hc_grad[0].re += PsiFactorA*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[0].GSq*LPsiDatA[0].re); Hc_grad[0].im += PsiFactorA*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[0].GSq*LPsiDatA[0].im); s++; } for (;sMaxG;s++) { Index = LevS->GArray[s].Index; /* FIXME - factoren */ Hc_grad[s].re += PsiFactorA*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[s].GSq*LPsiDatA[s].re); Hc_grad[s].im += PsiFactorA*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[s].GSq*LPsiDatA[s].im); } // 1/2 |G|^2 + V^{ps,nl} (G) + V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle } else { SetArrayToDouble0((double *)Hc_grad,2*LevS->MaxG); } m = -1; // former border of following loop: Psi->GlobalNo[Psi->PsiST+1]+Psi->GlobalNo[Psi->PsiST+4]+P->Par.Max_me_comm_ST_PsiT for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) { // go through all wave functions OnePsiB = &Psi->AllPsiStatus[j]; // grab OnePsiA if (OnePsiB->PsiType <= UnOccupied) { // drop the extra ones m++; if (OnePsiB->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local? LOnePsiB = &Psi->LocalPsiStatus[OnePsiB->MyLocalNo]; else LOnePsiB = NULL; if (LOnePsiB == NULL) { // if it's not local ... receive it from respective process into TempPsi RecvSource = OnePsiB->my_color_comm_ST_Psi; MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, HamiltonianTag, P->Par.comm_ST_PsiT, &status ); LPsiDatB=LevS->LPsi->TempPsi; } else { // .. otherwise send it to all other processes (Max_me... - 1) for (k=0;kPar.Max_me_comm_ST_PsiT;k++) if (k != OnePsiB->my_color_comm_ST_Psi) MPI_Send( LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, k, HamiltonianTag, P->Par.comm_ST_PsiT); LPsiDatB=LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo]; } // LPsiDatB is now set to the coefficients of OnePsi either stored or MPI_Received PsiFactorB = 1; //Psi->LocalPsiStatus[OnePsiB->MyLocalNo].PsiFactor; lambda = 0; //PsiFactorB * GradSP(P, LevS, LPsiDatB, Hc_grad); Lambda = 0; s=0; if (LevS->GArray[0].GSq == 0.0) { lambda += GSL_REAL(gsl_complex_mul(gsl_complex_conjugate(convertComplex(LPsiDatB[s])), convertComplex(Hc_grad[s]))); s++; } for (; s < LevS->MaxG; s++) { lambda += GSL_REAL(gsl_complex_mul(gsl_complex_conjugate(convertComplex(LPsiDatB[s])), convertComplex(Hc_grad[s]))); lambda += GSL_REAL(gsl_complex_mul(convertComplex(LPsiDatB[s]), gsl_complex_conjugate(convertComplex(Hc_grad[s])))); // c(-G) = c^\ast (G) due to gamma point symmetry ! } // due to the symmetry the resulting lambda is both real and symmetric in l,m (i,j) ! lambda *= PsiFactorB; // \langle \Psi_a | 1/2 |G|^2 + V^{ps,nl} (G) + V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle MPI_Allreduce ( &lambda, &Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST); // gather from all(!) (coeffs- and Psis-sharing) groups //if (P->Par.me == 0) { // and Output if (P->Call.out[StepLeaderOut]) fprintf(stderr,"%2.5lg\t",Lambda); if (P->Par.me == 0) fprintf(F->HamiltonianFile, "%lg\t", Lambda); //} if (m < NoOfPsis && l < NoOfPsis) { //fprintf(stderr,"Index (%i/%i,%i/%i): %e \n",l,NoOfPsis,m,NoOfPsis, Lambda); gsl_matrix_set(H,l,m,Lambda); } else fprintf(stderr,"Index (%i,%i) out of range %i\n",l,m,NoOfPsis); } if ((P->Par.me == 0) && (j == NoOfPsis-1)) fprintf(F->HamiltonianFile, "\n"); if ((P->Call.out[StepLeaderOut]) && (j == NoOfPsis-1)) fprintf(stderr,"\n"); } } } if (P->Par.me == 0) { fprintf(F->HamiltonianFile, "\n"); fflush(F->HamiltonianFile); } if ((P->Call.out[LeaderOut]) && (!P->Call.out[StepLeaderOut])) { fprintf(stderr,"finished.\n"); } // storing H matrix for latter use in perturbed calculation for (i=0;iNoOfPsis;i++) for (j=0;jNoOfPsis;j++) Psi->lambda[i][j] = gsl_matrix_get(H,i,j); // retrieve EW from H gsl_vector *eval = gsl_vector_alloc(NoOfPsis); gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc(NoOfPsis); gsl_eigen_symm(H, eval, w); gsl_eigen_symm_free(w); gsl_matrix_free(H); // print eigenvalues if(P->Call.out[LeaderOut]) { fprintf(stderr,"(%i) Unsorted Eigenvalues:", P->Par.me); for (i=0;iCall.out[LeaderOut]) { fprintf(stderr,"(%i) All Eigenvalues:", P->Par.me); for (i=0;iPar.me); for (i=0;iAddData[i].Lambda); fprintf(stderr,"\n"); } // print difference between highest occupied and lowest unoccupied if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0) { Lat->Energy[UnOccupied].bandgap = gsl_vector_get(eval, Psi->NoOfPsis) - gsl_vector_get(eval, Psi->NoOfPsis-1); //Lat->Energy[UnOccupied].homolumo = gsl_vector_get(eval, Psi->NoOfPsis/2 + 1) - gsl_vector_get(eval, (Psi->NoOfPsis-1)/2); if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) Band gap: \\epsilon_{%d+%d}(%d+1) - \\epsilon_{%d+%d}(%d) = %lg Ht\n", P->Par.me, NoPsis,P->Lat.Psi.GlobalNo[PsiMaxAdd],NoPsis, NoPsis, P->Lat.Psi.GlobalNo[PsiMaxAdd], NoPsis, Lat->Energy[UnOccupied].bandgap); //fprintf(stderr,"(%i) HOMO-LUMO gap (HOMO is: %i, LUMO is %i): %lg Ht\n", P->Par.me, Psi->NoOfPsis/2 + 1, (Psi->NoOfPsis-1)/2, Lat->Energy[UnOccupied].homolumo); } // now, calculate \int v_xc (r) n(r) d^3r SetArrayToDouble0((double *)HGCR,2*Dens0->TotalSize); CalculateXCPotentialNoRT(P, HGCR); /* FIXME - kann man halbieren */ pot_xc = 0; Pot_xc = 0; for (i=0;iLocalSizeR;i++) pot_xc += HGCR[i]*Dens0->DensityArray[TotalDensity][i]; MPI_Allreduce ( &pot_xc, &Pot_xc, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi); if (P->Call.out[StepLeaderOut]) { fprintf(stderr,"(%i) Exchange-Correlation potential energy: \\int v_xc[n] (r) n(r) d^3r = %lg\n", P->Par.me, HGcRCFactor*Pot_xc*Lat->Volume); } gsl_vector_free(eval); } /** Calculation of direction of steepest descent. * The Gradient has been evaluated in CalculateGradient(). Now the SD-direction * \f$\zeta_i = - (H-\lambda_i)\psi_i\f$ is stored in * Gradient::GradientArray[GradientTypes#GraSchGradient]. * Afterwards, the "extra" local Psi (which is the above gradient) is orthogonalized * under exclusion of the to be minimialized Psi to ascertain Orthonormality of this * gradient with respect to all other states/orbits. * \param *P Problem at hand, contains Lattice and RunStruct * \param *source source wave function coefficients, \f$\psi_i\f$ * \param *grad gradient coefficients, see CalculateGradient(), on return contains conjugate gradient coefficients, \f$\varphi_i^{(m)}\f$ * \param *Lambda eigenvalue \f$\lambda_i\f$ = \f$\langle \psi_i^{(m)}|H|\psi_i^{(m)} \rangle\f$ * \warning coefficients in \a *grad are overwritten */ static void CalculateCGGradient(struct Problem *P, fftw_complex *source, fftw_complex *grad, double Lambda) { struct Lattice *Lat = &P->Lat; struct RunStruct *R = &P->R; struct Psis *Psi = &Lat->Psi; struct LatticeLevel *LevS = R->LevS; fftw_complex *GradOrtho = P->Grad.GradientArray[GraSchGradient]; int g; int ElementSize = (sizeof(fftw_complex) / sizeof(double)); for (g=0;gMaxG;g++) { GradOrtho[g].re = grad[g].re+Lambda*source[g].re; GradOrtho[g].im = grad[g].im+Lambda*source[g].im; } // include the ExtraPsi (which is the GraSchGradient!) SetGramSchExtraPsi(P, Psi, NotOrthogonal); // exclude the minimised Psi SetGramSchActualPsi(P, Psi, NotUsedToOrtho); SpeedMeasure(P, GramSchTime, StartTimeDo); // makes conjugate gradient orthogonal to all other orbits //fprintf(stderr,"CalculateCGGradient: GramSch() for extra orbital\n"); GramSch(P, LevS, Psi, Orthogonalize); SpeedMeasure(P, GramSchTime, StopTimeDo); SetGramSchActualPsi(P, Psi, IsOrthonormal); memcpy(grad, GradOrtho, ElementSize*LevS->MaxG*sizeof(double)); } /** Preconditioning for the gradient calculation. * A suitable matrix for preconditioning is [TPA89] * \f[ * M_{G,G'} = \delta_{G,G'} \frac{27+18x+12x^2+8x^3}{27+18x+12x^2+8x^3+16x^4} \qquad (5.3) * \f] * where \f$x= \frac{\langle \xi_G|-\frac{1}{2}\nabla^2|\xi_G \rangle} * {\langle \psi_i^{(m)}|-\frac{1}{2}\nabla^2|\psi_i^{(m)} \rangle}\f$. * (This simple form is due to the diagonal dominance of the kinetic term in Hamiltonian, and * x is much simpler to evaluate in reciprocal space!) * With this one the higher plane wave states are being transformed such that they * nearly equal the error vector and converge at an equal rate. * So that the new steepest descent direction is \f$M \tilde{\zeta}_i\f$. * Afterwards, the new direction is orthonormalized via GramSch() and the "extra" * local Psi in a likewise manner to CalculateCGGradient(). * * \param *P Problem at hand, contains LatticeLevel and RunStruct * \param *grad array of steepest descent coefficients, \f$\varphi_i^{(m)}\f$ * \param *PCgrad return array for preconditioned steepest descent coefficients, \f$\tilde{\varphi}_i^{(m)}\f$ * \param T kinetic eigenvalue \f$\langle \psi_i | - \frac{1}{2} \nabla^2 | \psi_i \rangle\f$ */ static void CalculatePreConGrad(struct Problem *P, fftw_complex *grad, fftw_complex *PCgrad, double T) { struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; int g; double x, K, dK; struct Psis *Psi = &P->Lat.Psi; fftw_complex *PCOrtho = P->Grad.GradientArray[GraSchGradient]; int ElementSize = (sizeof(fftw_complex) / sizeof(double)); if (fabs(T) < MYEPSILON) T = 1; for (g=0;gMaxG;g++) { x = LevS->GArray[g].GSq; x /= T; K = (((8*x+12)*x +18)*x+27); dK = ((((16*x+8)*x +12)*x+18)*x+27); K /= dK; c_re(PCOrtho[g]) = K*c_re(grad[g]); c_im(PCOrtho[g]) = K*c_im(grad[g]); } SetGramSchExtraPsi(P, Psi, NotOrthogonal); SpeedMeasure(P, GramSchTime, StartTimeDo); // preconditioned direction is orthogonalized //fprintf(stderr,"CalculatePreConGrad: GramSch() for extra orbital\n"); GramSch(P, LevS, Psi, Orthogonalize); SpeedMeasure(P, GramSchTime, StopTimeDo); memcpy(PCgrad, PCOrtho, ElementSize*LevS->MaxG*sizeof(double)); } /** Calculation of conjugate direction. * The new conjugate direction is according to the algorithm of Fletcher-Reeves: * \f[ * \varphi^{(m)}_i = * \tilde{\eta}^{(m)}_i + \frac{\langle \tilde{\eta}^{(m)}_i | \tilde{\varphi}_i^{(m)} \rangle}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} \varphi_i^{(m-1)} * \f] * Instead we use the algorithm of Polak-Ribiere: * \f[ * \varphi^{(m)}_i = * \tilde{\eta}^{(m)}_i + \frac{\langle \tilde{\eta}^{(m)}_i | \bigl (\tilde{\varphi}_i^{(m)} \rangle - \tilde{\varphi}_i^{(m-1)} \rangle \bigr)}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} \varphi_i^{(m-1)} * \f] * with an additional automatic restart: * \f[ * \varphi^{(m)}_i = \tilde{\eta}^{(m)}_i \quad \textnormal{if}\ \frac{\langle \tilde{\eta}^{(m)}_i | \bigl (\tilde{\varphi}_i^{(m)} \rangle - \tilde{\varphi}_i^{(m-1)} \rangle \bigr)}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} <= 0 * \f] * That's why ConDirGradient and OldConDirGradient are stored, the scalar product - summed up among all processes - * is saved in OnePsiAddData::Gamma. The conjugate direction is again orthonormalized by GramSch() and stored * in \a *ConDir. * \param *P Problem at hand, containing LatticeLevel and RunStruct * \param step minimisation step m * \param *GammaDiv scalar product \f$\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle\f$, on calling expects (m-1), on return contains (m) * \param *grad array for coefficients of steepest descent direction \f$\tilde{\zeta}_i^{(m)}\f$ * \param *oldgrad array for coefficients of steepest descent direction of \a step-1 \f$\tilde{\zeta}_i^{(m-1)}\f$ * \param *PCgrad array for coefficients of preconditioned steepest descent direction \f$\tilde{\eta}_i^{(m)}\f$ * \param *ConDir return array for coefficients of conjugate array of \a step, \f$\tilde{\varphi}_i^{(m)}\f$ = \f$\frac{\tilde{\varphi}_i^{(m)}} {||\tilde{\varphi}_i^{(m)}||}\f$ * \param *ConDir_old array for coefficients of conjugate array of \a step-1, \f$\tilde{\varphi}_i^{(m-1)}\f$ */ static void CalculateConDir(struct Problem *P, int step, double *GammaDiv, fftw_complex *grad, fftw_complex *oldgrad, fftw_complex *PCgrad, fftw_complex *ConDir, fftw_complex *ConDir_old) { struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; struct Psis *Psi = &P->Lat.Psi; fftw_complex *Ortho = P->Grad.GradientArray[GraSchGradient]; int ElementSize = (sizeof(fftw_complex) / sizeof(double)); int g; double S[3], Gamma, GammaDivOld = *GammaDiv; double LocalSP[2], PsiSP[2]; LocalSP[0] = GradSP(P, LevS, PCgrad, grad); LocalSP[1] = GradSP(P, LevS, PCgrad, oldgrad); MPI_Allreduce ( LocalSP, PsiSP, 2, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi); *GammaDiv = S[0] = PsiSP[0]; S[1] = GammaDivOld; S[2] = PsiSP[1]; if (step) { // only in later steps is the scalar product used, but always condir stored in oldcondir and Ortho (working gradient) if (fabs(S[1]) < MYEPSILON) fprintf(stderr,"CalculateConDir: S[1] = %lg\n",S[1]); Gamma = (S[0]-S[2])/S[1]; // the final CG beta coefficient (Polak-Ribiere) if (fabs(S[1]) < MYEPSILON) { if (fabs((S[0]-S[2])) < MYEPSILON) Gamma = 1.0; else Gamma = 0.0; } if (Gamma < 0) { // Polak-Ribiere with automatic restart: gamma = max {0, gamma} Gamma = 0.; if (P->Call.out[StepLeaderOut]) fprintf(stderr,"(%i) Polak-Ribiere: Restarting iteration by setting gamma = 0.\n", P->Par.me); } for (g=0; g < LevS->MaxG; g++) { c_re(ConDir[g]) = c_re(PCgrad[g]) + Gamma*c_re(ConDir_old[g]); c_im(ConDir[g]) = c_im(PCgrad[g]) + Gamma*c_im(ConDir_old[g]); c_re(ConDir_old[g]) = c_re(ConDir[g]); c_im(ConDir_old[g]) = c_im(ConDir[g]); c_re(Ortho[g]) = c_re(ConDir[g]); c_im(Ortho[g]) = c_im(ConDir[g]); } } else { Gamma = 0.0; for (g=0; g < LevS->MaxG; g++) { c_re(ConDir[g]) = c_re(PCgrad[g]); c_im(ConDir[g]) = c_im(PCgrad[g]); c_re(ConDir_old[g]) = c_re(ConDir[g]); c_im(ConDir_old[g]) = c_im(ConDir[g]); c_re(Ortho[g]) = c_re(ConDir[g]); c_im(Ortho[g]) = c_im(ConDir[g]); } } // orthonormalize SetGramSchExtraPsi(P, Psi, NotOrthogonal); SpeedMeasure(P, GramSchTime, StartTimeDo); //fprintf(stderr,"CalculateConDir: GramSch() for extra orbital\n"); GramSch(P, LevS, Psi, Orthonormalize); SpeedMeasure(P, GramSchTime, StopTimeDo); memcpy(ConDir, Ortho, ElementSize*LevS->MaxG*sizeof(double)); } /** Calculation of \f$\langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\widetilde{\widetilde{\varphi}}_i^{(m)} \rangle\f$. * Evaluates most of the parts of the second derivative of the energy functional \f$\frac{\delta^2 E}{\delta \Theta^2}|_{\Theta=0}\f$, * most importantly the energy expection value of the conjugate direction. * \f[ * \langle \widetilde{\widetilde{\varphi}}_i^{(m)} (G) | -\frac{1}{2} \nabla^2 + \underbrace{V^H (G) + V^{ps,loc} (G)}_{=V^{HG}(G)} + V^{ps,nl} | \widetilde{\widetilde{\varphi}} (G) \rangle * \f] * The conjugate direction coefficients \f$\varphi_i^{(m)} (G)\f$ are fftransformed to \f$\varphi_i^{(m)} (R)\f$, * times the DoubleDensityTypes#HGDensity \f$V^{HG} (R) \f$we have the first part, * then \f$\rho (R) = 2\frac{f_i}{V} \Re (\varphi_i^{(m)} (R) \cdot \psi_i (R))\f$is evaluated, * HGDensity is fftransformed to \f$V^{HG} | \varphi_i^{(m)} (G) \rangle\f$, * CalculateAddNLPot() adds non-local potential \f$V^{ps,nl} (G)\f$ to it, * kinetic term is added in the reciprocal base as well \f$-\frac{1}{2} |G|^2\f$, * and at last scalar product between it and the conjugate direction, and done. * * Meanwhile the fftransformed ConDir density \f$\rho (R)\f$ is used for the ddEddt0s \f$A_H\f$, \f$A_{XC}\f$. * In case of RunStruct::DoCalcCGGauss in the following expression additionally \f$\hat{n}^{Gauss}\f$ * is evaluted, otherwise left out: * \f[ * A^H = \frac{4\pi}{V} \sum_G \frac{|\frac{\rho(G)}{V^2 N} + \hat{n}^{Gauss}|}{|G|^2} * \qquad (\textnormal{section 5.1, line search}) * \f] * \param *P Problem at hand * \param *ConDir conjugate direction \f$\tilde{\tilde{\varphi}}_i^{(m)}\f$ * \param PsiFactor occupation number, see Psis::PsiFactor * \param *ConDirHConDir first return value: energy expectation value for conjugate direction * \param *HartreeddEddt0 second return value: second derivative of Hartree energy \f$A_H\f$ * \param *XCddEddt0 third return value: second derivative of exchange correlation energy \f$A_{XC}\f$ * \warning The MPI_Allreduce for the scalar product in the end has not been done! * \sa CalculateLineSearch() - used there */ void CalculateConDirHConDir(struct Problem *P, fftw_complex *ConDir, const double PsiFactor, double *ConDirHConDir, double *HartreeddEddt0, double *XCddEddt0) { struct Lattice *Lat = &P->Lat; struct PseudoPot *PP = &P->PP; struct RunStruct *R = &P->R; struct LatticeLevel *Lev0 = R->Lev0; struct LatticeLevel *LevS = R->LevS; struct Density *Dens0 = Lev0->Dens; int Index, g, i, pos, s=0; struct fft_plan_3d *plan = Lat->plan; fftw_complex *work = Dens0->DensityCArray[TempDensity]; fftw_real *PsiCD = (fftw_real *)work; fftw_complex *PsiCDC = work; fftw_real *PsiR = (fftw_real *)Dens0->DensityCArray[ActualPsiDensity]; fftw_complex *tempdestRC = (fftw_complex *)Dens0->DensityArray[TempDensity]; fftw_real *HGCDR = Dens0->DensityArray[HGcDensity]; fftw_complex *HGCDRC = (fftw_complex*)HGCDR; fftw_complex *HGC = Dens0->DensityCArray[HGDensity]; fftw_real *HGCR = (fftw_real *)HGC; fftw_complex *posfac, *destpos, *destRCS, *destRCD; const double HartreeFactorR = 2.*R->FactorDensityR*PsiFactor; double HartreeFactorC = 4.*PI*Lat->Volume; double HartreeFactorCGauss = 1./Lev0->MaxN; //R->FactorDensityC*R->HGcFactor/Lev0->MaxN; // same expression as first and second factor cancel each other int nx,ny,nz,iS,i0; const int Nx = LevS->Plan0.plan->local_nx; const int Ny = LevS->Plan0.plan->N[1]; const int Nz = LevS->Plan0.plan->N[2]; const int NUpx = LevS->NUp[0]; const int NUpy = LevS->NUp[1]; const int NUpz = LevS->NUp[2]; const double HGcRCFactor = 1./LevS->MaxN; fftw_complex *HConDir = P->Grad.GradientArray[TempGradient]; int DoCalcCGGauss = R->DoCalcCGGauss; fftw_complex rp, rhog; int it; struct Ions *I = &P->Ion; /* ConDir H ConDir */ // (G)->(R): retrieve density from wave functions into PsiCD in the usual manner (yet within the same level!) (see CalculateOneDensityR()) SetArrayToDouble0((double *)tempdestRC, Dens0->TotalSize*2); for (i=0;iMaxG;i++) { Index = LevS->GArray[i].Index; posfac = &LevS->PosFactorUp[LevS->MaxNUp*i]; destpos = &tempdestRC[LevS->MaxNUp*Index]; for (pos=0; pos < LevS->MaxNUp; pos++) { destpos[pos].re = ConDir[i].re*posfac[pos].re-ConDir[i].im*posfac[pos].im; destpos[pos].im = ConDir[i].re*posfac[pos].im+ConDir[i].im*posfac[pos].re; } } for (i=0; iMaxDoubleG; i++) { destRCS = &tempdestRC[LevS->DoubleG[2*i]*LevS->MaxNUp]; destRCD = &tempdestRC[LevS->DoubleG[2*i+1]*LevS->MaxNUp]; for (pos=0; pos < LevS->MaxNUp; pos++) { destRCD[pos].re = destRCS[pos].re; destRCD[pos].im = -destRCS[pos].im; } } fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNFUp, tempdestRC, work); DensityRTransformPos(LevS,(fftw_real*)tempdestRC, PsiCD); // HGCR contains real(!) HGDensity from CalculateGradientNoRT() ! for (nx=0;nxLocalSizeR; i++) PsiCD[i] *= HartreeFactorR*PsiR[i]; // is rho(r) = PsiFactor times \psi_i is in PsiR times 1/V times \varphi_i^{(m)} *XCddEddt0 = CalculateXCddEddt0NoRT(P, PsiCD); // calcs A^{XC} fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF1, PsiCDC, tempdestRC); s = 0; *HartreeddEddt0 = 0.0; if (Lev0->GArray[0].GSq == 0.0) s++; if (DoCalcCGGauss) { // this part is wrong, n^Gauss is not part of HartreeddEddt0 (only appears in 2*( - ...)) for (g=s; g < Lev0->MaxG; g++) { Index = Lev0->GArray[g].Index; rp.re = 0.0; rp.im = 0.0; for (it = 0; it < I->Max_Types; it++) { c_re(rp) += (c_re(I->I[it].SFactor[g])*PP->FacGauss[it][g]); c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]); } c_re(rhog) = c_re(PsiCDC[Index])*HartreeFactorCGauss+c_re(rp); c_im(rhog) = c_im(PsiCDC[Index])*HartreeFactorCGauss+c_im(rp); if (Lev0->GArray[g].GSq < MYEPSILON) fprintf(stderr,"CalculateConDirHConDir: Lev0->GArray[g].GSq = %lg",Lev0->GArray[g].GSq); *HartreeddEddt0 += (c_re(rhog)*c_re(rhog)+c_im(rhog)*c_im(rhog))/(Lev0->GArray[g].GSq); } } else { HartreeFactorC *= HartreeFactorCGauss*HartreeFactorCGauss; for (g=s; g < Lev0->MaxG; g++) { Index = Lev0->GArray[g].Index; if (Lev0->GArray[g].GSq < MYEPSILON) fprintf(stderr,"CalculateConDirHConDir: Lev0->GArray[g].GSq = %lg",Lev0->GArray[g].GSq); *HartreeddEddt0 += (c_re(PsiCDC[Index])*c_re(PsiCDC[Index])+c_im(PsiCDC[Index])*c_im(PsiCDC[Index]))/(Lev0->GArray[g].GSq); } } *HartreeddEddt0 *= HartreeFactorC; /* fertig mit ddEddt0 */ fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGCDRC, work); CalculateCDfnl(P, ConDir, PP->CDfnl); // calculate needed non-local form factors CalculateAddNLPot(P, HConDir, PP->CDfnl, PsiFactor); // add non-local potential, sets to zero beforehand //SetArrayToDouble0((double *)HConDir,2*LevS->MaxG); // now we have V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle for (g=0;gMaxG;g++) { // factor by one over volume and add kinetic part Index = LevS->GArray[g].Index; /* FIXME - factoren */ HConDir[g].re += PsiFactor*(HGCDRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[g].GSq*ConDir[g].re); HConDir[g].im += PsiFactor*(HGCDRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[g].GSq*ConDir[g].im); } // now, T^{kin} + V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle *ConDirHConDir = GradSP(P, LevS, ConDir, HConDir); // finally, \langle \varphi | V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle which is returned } /** Find the minimal \f$\Theta_{min}\f$ for the line search. * We find it as * \f[ * \widetilde{\Theta}_{min} = \frac{1}{2} \tan^{-1} \Bigl ( -\frac{ \frac{\delta E}{\delta \Theta}_{|\Theta =0} } * { \frac{1}{2} \frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0} } \Bigr ) * \f] * Thus, with given \f$E(0)\f$, \f$\frac{\delta E}{\delta \Theta}_{|\Theta =0}\f$ and \f$\frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0}\f$ * \f$\Theta_{min}\f$is calculated and the unknowns A, B and C (see CalculateLineSearch()) obtained. So that our * approximating formula \f$\widetilde{E}(\Theta) = A + B cos(2\Theta) + CSin(2\Theta)\f$ is fully known and we * can evaluate the energy functional and its derivative at \f$\Theta_{min} \f$. * \param E0 energy eigenvalue at \f$\Theta = 0\f$: \f$A= E(0) - B\f$ * \param dEdt0 first derivative of energy eigenvalue at \f$\Theta = 0\f$: \f$C = \frac{1}{2} \frac{\delta E}{\delta \Theta}|_{\Theta=0}\f$ * \param ddEddt0 first derivative of energy eigenvalue at \f$\Theta = 0\f$: \f$B = -\frac{1}{4} \frac{\delta^2 E}{\delta \Theta^2}|_{\Theta=0}\f$ * \param *E returned energy functional at \f$\Theta_{min}\f$ * \param *dE returned first derivative energy functional at \f$\Theta_{min}\f$ * \param *ddE returned second derivative energy functional at \f$\Theta_{min}\f$ * \param *dcos cosine of the found minimal angle * \param *dsin sine of the found minimal angle * \return \f$\Theta_{min}\f$ * \sa CalculateLineSearch() - used there */ double CalculateDeltaI(double E0, double dEdt0, double ddEddt0, double *E, double *dE, double *ddE, double *dcos, double *dsin) { double delta, dcos2, dsin2, A, B, Eavg; if (fabs(ddEddt0) < MYEPSILON) fprintf(stderr,"CalculateDeltaI: ddEddt0 = %lg",ddEddt0); delta = atan(-2.*dEdt0/ddEddt0)/2.0; //.5* *dcos = cos(delta); *dsin = sin(delta); dcos2 = cos(delta*2.0); //2.* dsin2 = sin(delta*2.0); //2.* A = -ddEddt0/4.; B = dEdt0/2.; //fprintf(stderr,"CalculateDeltaI: A = %lg, B = %lg\n", A, B); Eavg = E0-A; *E = Eavg + A*dcos2 + B*dsin2; *dE = -2.*A*dsin2+2.*B*dcos2; *ddE = -4.*A*dcos2-4.*B*dsin2; if (*ddE < 0) { // if it's a maximum (indicated by positive second derivative at found theta) if (delta < 0) { // go the minimum - the two are always separated by PI/2 by the ansatz delta += PI/2.0; ///2. } else { delta -= PI/2.0; ///2. } *dcos = cos(delta); *dsin = sin(delta); dcos2 = cos(delta*2.0); //2.* dsin2 = sin(delta*2.0); //2.* *E = Eavg + A*dcos2 + B*dsin2; *dE = -2.*A*dsin2+2.*B*dcos2; *ddE = -4.*A*dcos2-4.*B*dsin2; } return(delta); } /** Line search: One-dimensional minimisation along conjugate direction. * We want to minimize the energy functional along the calculated conjugate direction. * States \f$\psi_i^{(m)}(\Theta) = \widetilde{\psi}_i^{(m)} \cos(\Theta) + \widetilde{\widetilde{\varphi}}_i^{(m)} * \sin(\Theta)\f$ are normalized and orthogonal for every \f$\Theta \in {\cal R}\f$. Thus we * want to find the minimal \f$\Theta_{min}\f$. Intenting to evaluate the energy functional as * seldom as possible, we approximate it by: * \f$\widetilde{E}(\Theta) = A + B cos(2\Theta) + CSin(2\Theta)\f$ * where the functions shall match the energy functional and its first and second derivative at * \f$\Theta = 0\f$, fixing the unknowns A, B and C.\n * \f$\frac{\delta E}{\delta \Theta}_{|\Theta =0} = * 2 \Re \bigl ( \langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\psi_i^{(m)} \rangle \bigr )\f$ * and * \f$\frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0} = * 2\bigl( \langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\widetilde{\widetilde{\varphi}}_i^{(m)} \rangle * - \langle \psi_i^{(m)}|H|\psi_i^{(m)} \rangle * + \underbrace{\int_\Omega \nu_H[\rho](r)\rho(r)dr}_{A_H} * + \underbrace{\int_\Omega \rho(r)^2 \frac{\delta V^{xc}}{\delta n}(r) dr}_{A_{xc}}\f$. \n * \f$A_H\f$ and \f$A_{xc}\f$ can be calculated from the Hartree respectively the exchange * correlation energy.\n * Basically, this boils down to using CalculateConDirHConDir() (we already have OnePsiAddData::Lambda), * in evaluation of the energy functional and its derivatives as shown in the above formulas * and calling with these CalculateDeltaI() which finds the above \f$\Theta_{min} \f$. \n * The coefficients of the current wave function are recalculated, thus we have the new wave function, * the total energy and its derivatives printed to stderr.\n * * \param *P Problem at hand * \param *source current minimised wave function * \param *ConDir conjugate direction, \f$\tilde{\tilde{\varphi}}_i^{(m)}\f$ * \param *Hc_grad complex coefficients of \f$H | \psi_i (G) \rangle\f$, see GradientArray#HcGradient * \param x if NULL do approximative line search by second order taylor expansion, otherwise use given * parameter between -M_PI...M_PI. Is simply put through to CalculateLineSearch() * \param PsiFactor occupation number */ static void CalculateLineSearch(struct Problem *P, fftw_complex *source, fftw_complex *ConDir, fftw_complex *Hc_grad, const double PsiFactor, const double *x) { struct Lattice *Lat = &P->Lat; struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; struct FileData *F = &P->Files; //struct LatticeLevel *Lev0 = R->Lev0; //struct Density *Dens = Lev0->Dens; //struct Psis *Psi = &P->Lat.Psi; struct Energy *En = Lat->E; //int ElementSize = (sizeof(fftw_complex) / sizeof(double)); double dEdt0, ddEddt0, ConDirHConDir, Eavg, A, B;//, sourceHsource; double E0, E, dE, ddE, delta, dcos, dsin, dcos2, dsin2; double EI, dEI, ddEI, deltaI, dcosI, dsinI; double HartreeddEddt0, XCddEddt0; double d[4],D[4], Diff; //double factorDR = R->FactorDensityR; int g, i; //int ActNum; // ========= dE / dt | 0 ============ //fprintf(stderr,"(%i) |ConDir| = %lg, |Hc_grad| = %lg\n", P->Par.me, sqrt(GradSP(P, LevS, ConDir, ConDir)), sqrt(GradSP(P, LevS, Hc_grad, Hc_grad))); d[0] = 2.*GradSP(P, LevS, ConDir, Hc_grad); // ========== ddE / ddt | 0 ========= CalculateConDirHConDir(P, ConDir, PsiFactor, &d[1], &d[2], &d[3]); MPI_Allreduce ( &d, &D, 4, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi); dEdt0 = D[0]; for (i=MAXOLD-1; i > 0; i--) En->dEdt0[i] = En->dEdt0[i-1]; En->dEdt0[0] = dEdt0; ConDirHConDir = D[1]; HartreeddEddt0 = D[2]; XCddEddt0 = D[3]; ddEddt0 = 0.0; switch (R->CurrentMin) { case Occupied: ddEddt0 = HartreeddEddt0+XCddEddt0; // these terms drop out in unoccupied variation due to lacking dependence of n^{tot} case UnOccupied: ddEddt0 += 2.*(ConDirHConDir - Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda); break; default: ddEddt0 = 0.0; break; } //if (R->CurrentMin <= UnOccupied) //fprintf(stderr,"ConDirHConDir = %lg\tLambda = %lg\t Hartree = %lg\t XC = %lg\n",ConDirHConDir,Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda,HartreeddEddt0,XCddEddt0); for (i=MAXOLD-1; i > 0; i--) En->ddEddt0[i] = En->ddEddt0[i-1]; En->ddEddt0[0] = ddEddt0; E0 = En->TotalEnergy[0]; if (x == NULL) { // delta //if (isnan(E0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): E0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); } //if (isnan(dEdt0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): dEdt0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); } //if (isnan(ddEddt0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): ddEddt0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); } deltaI = CalculateDeltaI(E0, dEdt0, ddEddt0, &EI, &dEI, &ddEI, &dcosI, &dsinI); } else { deltaI = *x; dcosI = cos(deltaI); dsinI = sin(deltaI); dcos2 = cos(deltaI*2.0); //2.* dsin2 = sin(deltaI*2.0); //2.* A = -ddEddt0/4.; B = dEdt0/2.; //fprintf(stderr,"CalculateDeltaI: A = %lg, B = %lg\n", A, B); Eavg = E0-A; EI = Eavg + A*dcos2 + B*dsin2; dEI = -2.*A*dsin2+2.*B*dcos2; ddEI = -4.*A*dcos2-4.*B*dsin2; } delta = deltaI; E = EI; dE = dEI; ddE = ddEI; dcos = dcosI; dsin = dsinI; if (R->ScanPotential >= R->ScanAtStep) { delta = (R->ScanPotential-(R->ScanAtStep-1))*R->ScanDelta; dcos = cos(delta); dsin = sin(delta); //fprintf(stderr,"(%i) Scaning Potential form with delta %lg\n", P->Par.me, delta); } // shift energy delta values for (i=MAXOLD-1; i > 0; i--) { En->delta[i] = En->delta[i-1]; En->ATE[i] = En->ATE[i-1]; } // store new one En->delta[0] = delta; En->ATE[0] = E; if (R->MinStep != 0) Diff = fabs(En->TotalEnergy[1] - E0)/(En->TotalEnergy[1] - E0)*fabs((E0 - En->ATE[1])/E0); else Diff = 0.; // rotate local Psi and ConDir according to angle Theta for (g = 0; g < LevS->MaxG; g++) { // Here all coefficients are updated for the new found wave function //if (isnan(ConDir[g].re)) { fprintf(stderr,"WARNGING: CalculateLineSearch(): ConDir_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); } c_re(source[g]) = c_re(source[g])*dcos + c_re(ConDir[g])*dsin; c_im(source[g]) = c_im(source[g])*dcos + c_im(ConDir[g])*dsin; } // print to stderr if (P->Call.out[StepLeaderOut] && x == NULL) { if (1) switch (R->CurrentMin) { case UnOccupied: fprintf(stderr, "(%i,%i,%i)S(%i,%i,%i):\tTGE: %e\tATGE: %e\t Diff: %e\t --- d: %e\tdEdt0: %e\tddEddt0: %e\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, Lat->Energy[UnOccupied].TotalEnergy[0], E, Diff, delta, dEdt0, ddEddt0); break; default: fprintf(stderr, "(%i,%i,%i)S(%i,%i,%i):\tTE: %e\tATE: %e\t Diff: %e\t --- d: %e\tdEdt0: %e\tddEddt0: %e\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, Lat->Energy[Occupied].TotalEnergy[0], E, Diff, delta, dEdt0, ddEddt0); break; } } if ((P->Par.me == 0) && (R->ScanFlag || (R->ScanPotential < R->ScanAtStep))) { fprintf(F->MinimisationFile, "%i\t%i\t%i\t%e\t%e\t%e\t%e\t%e\n",R->MinStep, R->ActualLocalPsiNo, R->PsiStep, E0, E, delta, dEdt0, ddEddt0); fflush(F->MinimisationFile); } } /** Find the new minimised wave function. * Note, the gradient wave the function is expected in GradientTypes#ActualGradient. Then, * subsequently calls (performing thereby a conjugate gradient minimisation algorithm * of Fletchers&Reeves): * - CalculateCGGradient()\n * find direction of steepest descent: GradientTypes#ActualGradient * - CalculatePreConGrad()\n * precondition steepest descent direction: GradientTypes#PreConGradient * - CalculateConDir()\n * calculate conjugate direction: GradientTypes#ConDirGradient and GradientTypes#OldConDirGradient * - CalculateLineSearch()\n * perform the minimisation in this one-dimensional subspace: LevelPsi::LocalPsi[RunStruct::ActualLocalPsiNo] * \note CalculateGradient() is not re-done. * \param *P Problem at hand * \param x if NULL do approximative line search by second order taylor expansion, otherwise use given paramter between -M_PI...M_PI */ void CalculateNewWave(struct Problem *P, double *x) { struct Lattice *Lat = &P->Lat; struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; const int i = R->ActualLocalPsiNo; const int ConDirStep = R->PsiStep; if (((!R->ScanPotential) || (R->ScanPotential < R->ScanAtStep)) && (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent != 1)) { // only evaluate CG once //if (isnan(P->Grad.GradientArray[ActualGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): ActualGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); } //fprintf(stderr,"(%i) Evaluating ConDirGradient from SD direction...\n", P->Par.me); CalculateCGGradient(P, LevS->LPsi->LocalPsi[i], P->Grad.GradientArray[ActualGradient], Lat->Psi.AddData[i].Lambda); //fprintf(stderr,"CalculateCGGradient: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[ActualGradient], P->Grad.GradientArray[ActualGradient])); CalculatePreConGrad(P, P->Grad.GradientArray[ActualGradient], P->Grad.GradientArray[PreConGradient], Lat->Psi.AddData[i].T); //if (isnan(P->Grad.GradientArray[PreConGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): PreConGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); } //fprintf(stderr,"CalculatePreConGrad: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[PreConGradient], P->Grad.GradientArray[PreConGradient])); CalculateConDir(P, ConDirStep, &Lat->Psi.AddData[i].Gamma, P->Grad.GradientArray[ActualGradient], P->Grad.GradientArray[OldActualGradient], P->Grad.GradientArray[PreConGradient], P->Grad.GradientArray[ConDirGradient], P->Grad.GradientArray[OldConDirGradient]); //if (isnan(P->Grad.GradientArray[ConDirGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): ConDirGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); } //fprintf(stderr,"CalculateConDir: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[ConDirGradient], P->Grad.GradientArray[ConDirGradient])); } //else //TestForOrth(P,R->LevS,P->Grad.GradientArray[GraSchGradient]); //fprintf(stderr,"(%i) DoBrent[%i] = %i\n", P->Par.me, R->ActualLocalPsiNo, Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent); if (R->ScanPotential) R->ScanPotential++; // increment by one CalculateLineSearch(P, LevS->LPsi->LocalPsi[i], P->Grad.GradientArray[ConDirGradient], P->Grad.GradientArray[HcGradient], Lat->Psi.LocalPsiStatus[i].PsiFactor, x); //if (isnan(LevS->LPsi->LocalPsi[i][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): {\\wildetilde \\varphi}_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); } SetGramSchExtraPsi(P, &Lat->Psi, NotUsedToOrtho); //TestForOrth(P,R->LevS,P->Grad.GradientArray[GraSchGradient]); } /** Calculates the parameter alpha for UpdateWaveAfterIonMove() * Takes ratios of differences of old and older positions with old and current, * respectively with old and older. Alpha is then the ratio between the two: * \f[ * \alpha = \frac{v_x[old] \cdot v_x[now] + v_y[old] \cdot v_y[now] + v_z[old] \cdot v_z[now]}{v_x[old]^2 + v_y[old]^2 + v_z[old]^2} * \f] * Something like a rate of change in the ion velocity. * \param *P Problem at hand * \return alpha */ static double CalcAlphaForUpdateWaveAfterIonMove(struct Problem *P) { struct Ions *I = &P->Ion; struct RunStruct *Run = &P->R; double alpha=0; double Sum1=0; double Sum2=0; int is,ia; double *R,*R_old,*R_old_old; Run->AlphaStep++; if (Run->AlphaStep <= 1) return(0.0); for (is=0; is < I->Max_Types; is++) { for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { //if (I->I[is].IMT[ia] == MoveIon) { // Even FixedIon moves, only not by other's forces R = &I->I[is].R[NDIM*ia]; R_old = &I->I[is].R_old[NDIM*ia]; R_old_old = &I->I[is].R_old_old[NDIM*ia]; Sum1 += (R_old[0]-R_old_old[0])*(R[0]-R_old[0]); Sum1 += (R_old[1]-R_old_old[1])*(R[1]-R_old[1]); Sum1 += (R_old[2]-R_old_old[2])*(R[2]-R_old[2]); Sum2 += (R_old[0]-R_old_old[0])*(R_old[0]-R_old_old[0]); Sum2 += (R_old[1]-R_old_old[1])*(R_old[1]-R_old_old[1]); Sum2 += (R_old[2]-R_old_old[2])*(R_old[2]-R_old_old[2]); //} } } if (fabs(Sum2) > MYEPSILON) { alpha = Sum1/Sum2; } else { fprintf(stderr,"(%i) Sum2 <= MYEPSILON\n", P->Par.me); alpha = 0.0; } //if (isnan(alpha)) { fprintf(stderr,"(%i) CalcAlphaForUpdateWaveAfterIonMove: alpha is NaN!\n", P->Par.me); Error(SomeError, "NaN-Fehler!");} return (alpha); } /** Updates Wave functions after UpdateIonsR. * LocalPsi coefficients are shifted by the difference between LevelPsi::LocalPsi and LevelPsi::OldLocalPsi * multiplied by alpha, new value stored in old one. * Afterwards all local wave functions are set to NotOrthogonal and then orthonormalized via * GramSch(). * * \param *P Problem at hand * \sa CalcAlphaForUpdateWaveAfterIonMove() - alpha is calculated here */ double UpdateWaveAfterIonMove(struct Problem *P) #if 0 static double CalculateDeltaII(double E0, double dEdt0, double E1, double deltastep, double *E, double *dE, double *ddE, double *dcos, double *dsin) { double delta, dcos2, dsin2, A, B, Eavg; if (fabs(1-cos(2*deltastep) < MYEPSILON) fprintf(stderr,"CalculateDeltaII: (1-cos(2*deltastep) = %lg",(1-cos(2*deltastep)); A = (E0-E1+0.5*dEdt0*sin(2*deltastep))/(1-cos(2*deltastep)); B = 0.5*dEdt0; Eavg = E0-A; if (fabs(A) < MYEPSILON) fprintf(stderr,"CalculateDeltaII: A = %lg",A); delta = 0.5*atan(B/A); *dcos = cos(delta); *dsin = sin(delta); dcos2 = cos(2.*delta); dsin2 = sin(2.*delta); *E = Eavg + A*dcos2 + B*dsin2; *dE = -2.*A*dsin2+2.*B*dcos2; *ddE = -4.*A*dcos2-4.*B*dsin2; if (*ddE < 0) { if (delta < 0) { delta += PI/2.; } else { delta -= PI/2.; } *dcos = cos(delta); *dsin = sin(delta); dcos2 = cos(2.*delta); dsin2 = sin(2.*delta); *E = Eavg + A*dcos2 + B*dsin2; *dE = -2.*A*dsin2+2.*B*dcos2; *ddE = -4.*A*dcos2-4.*B*dsin2; } return(delta); } #endif { struct RunStruct *R = &P->R; struct LatticeLevel *LevS = R->LevS; struct Lattice *Lat = &P->Lat; struct Psis *Psi = &Lat->Psi; double alpha=CalcAlphaForUpdateWaveAfterIonMove(P); int i,g,p,ResetNo=0; fftw_complex *source, *source_old; struct OnePsiElement *OnePsi = NULL; //fprintf(stderr, "%lg\n", alpha); if (P->R.AlphaStep > 1 && alpha != 0.0) { for (i=Psi->TypeStartIndex[Occupied]; i < Psi->TypeStartIndex[Perturbed_P0]; i++) { source = LevS->LPsi->LocalPsi[i]; source_old = LevS->LPsi->OldLocalPsi[i]; for (g=0; g < LevS->MaxG; g++) { c_re(source[g]) = c_re(source[g])+alpha*(c_re(source[g])-c_re(source_old[g])); c_im(source[g]) = c_im(source[g])+alpha*(c_im(source[g])-c_im(source_old[g])); c_re(source_old[g])=c_re(source[g]); c_im(source_old[g])=c_im(source[g]); } } } else { for (i=Psi->TypeStartIndex[Occupied]; i < Psi->TypeStartIndex[Perturbed_P0]; i++) { source = LevS->LPsi->LocalPsi[i]; source_old = LevS->LPsi->OldLocalPsi[i]; for (g=0; g < LevS->MaxG; g++) { c_re(source_old[g])=c_re(source[g]); c_im(source_old[g])=c_im(source[g]); } } } //fprintf(stderr, "(%i) alpha = %lg\n", P->Par.me, alpha); if (alpha != 0.0) { for (p=0; p< Psi->LocalNo; p++) { Psi->LocalPsiStatus[p].PsiGramSchStatus = (Psi->LocalPsiStatus[p].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho; //if (R->CurrentMin > UnOccupied) //fprintf(stderr,"(%i) Setting L-Status of %i to %i\n",P->Par.me,p, Psi->LocalPsiStatus[p].PsiGramSchStatus); } ResetNo=0; for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) { for (p=ResetNo; p < ResetNo+Psi->AllLocalNo[i]-1; p++) { Psi->AllPsiStatus[p].PsiGramSchStatus = (Psi->AllPsiStatus[p].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho; //if (R->CurrentMin > UnOccupied) //fprintf(stderr,"(%i) Setting A-Status of %i to %i\n",P->Par.me,p, Psi->AllPsiStatus[p].PsiGramSchStatus); } ResetNo += Psi->AllLocalNo[i]; OnePsi = &Psi->AllPsiStatus[ResetNo-1]; // extra wave function is set to NotUsed //if (R->CurrentMin > UnOccupied) fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo-1, NotUsedToOrtho); OnePsi->PsiGramSchStatus = (int)NotUsedToOrtho; } SpeedMeasure(P, InitGramSchTime, StartTimeDo); //fprintf(stderr,"(%i) UpdateWaveAfterIonMove: ResetGramSch() for %i orbitals\n",P->Par.me, p); GramSch(P, LevS, Psi, Orthonormalize); SpeedMeasure(P, InitGramSchTime, StartTimeDo); } return(alpha); }