Index: pcp/src/grad.c
===================================================================
--- pcp/src/grad.c	(revision 7853f74d2acdb0b0a9aaedb67635b8dd47d82bb4)
+++ pcp/src/grad.c	(revision b76b9719af2de1d1c075afe06b09e2200400a72c)
@@ -189,7 +189,8 @@
   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!"); }
+  //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;nx<Nx;nx++)  
@@ -199,5 +200,6 @@
       	iS = nz+Nz*(ny+Ny*nx);
       	HGcR[iS] = HGCR[i0]*PsiCR[i0]; /* Matrix Vector Mult */
-        //if (isnan(HGCR[i0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGcR[%i]_%i = NaN!\n", i0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
+        //if (isnan(PsiCR[i0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): PsiCR[%i]_%i = NaN!\n", i0, R->ActualLocalPsiNo); 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);
@@ -467,6 +469,8 @@
   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
@@ -1044,4 +1048,5 @@
   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], 
@@ -1049,5 +1054,4 @@
   		      Lat->Psi.AddData[i].Lambda);
     //fprintf(stderr,"CalculateCGGradient: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[ActualGradient], P->Grad.GradientArray[ActualGradient]));
-    //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!"); }
     CalculatePreConGrad(P, P->Grad.GradientArray[ActualGradient], 
   		      P->Grad.GradientArray[PreConGradient], 
@@ -1199,9 +1203,10 @@
     }
   }
-  //fprintf(stderr, "%lg\n", alpha);
+  //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,"Setting L-Status of %i to %i\n",p, Psi->LocalPsiStatus[p].PsiGramSchStatus);
+      //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;
@@ -1209,5 +1214,6 @@
       for (p=ResetNo; p < ResetNo+Psi->AllLocalNo[i]-1; p++) {
 				Psi->AllPsiStatus[p].PsiGramSchStatus = (Psi->LocalPsiStatus[p].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho;
-        //if (R->CurrentMin > UnOccupied) fprintf(stderr,"Setting A-Status of %i to %i\n",p, Psi->AllPsiStatus[p].PsiGramSchStatus);
+        //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];
@@ -1217,5 +1223,5 @@
     }
     SpeedMeasure(P, InitGramSchTime, StartTimeDo);
-    //fprintf(stderr,"UpdateWaveAfterIonMove: ResetGramSch() for %i orbital\n",p);
+    //fprintf(stderr,"(%i) UpdateWaveAfterIonMove: ResetGramSch() for %i orbitals\n",P->Par.me, p);
     GramSch(P, LevS, Psi, Orthonormalize);
     SpeedMeasure(P, InitGramSchTime, StartTimeDo);
