Index: pcp/src/run.c
===================================================================
--- pcp/src/run.c	(revision 02ba605529127c06a6a320265561553a54e81414)
+++ pcp/src/run.c	(revision 4871821ec548a5d804253c2387c34663c0a2872d)
@@ -711,4 +711,5 @@
  * \param *Stop flag to determine if epsilon stop conditions have met
  * \param *SuperStop flag to determinte whether external signal's required end of calculations
+ * \bug ResetGramSch() not allowed after reading orthonormal values from file
  */
 static void MinimiseOccupied(struct Problem *P, int *Stop, int *SuperStop)
@@ -747,13 +748,14 @@
       InitPsisValue(P, Psi->TypeStartIndex[Occupied], Psi->TypeStartIndex[Occupied+1]);  // initialize perturbed array for this run
       ResetGramSchTagType(P, Psi, Occupied, NotOrthogonal); // loaded values are orthonormal
-      SpeedMeasure(P, InitGramSchTime, StartTimeDo);  
-      GramSch(P, R->LevS, Psi, Orthonormalize);
-      SpeedMeasure(P, InitGramSchTime, StopTimeDo);  
     } else {
       SpeedMeasure(P, InitSimTime, StartTimeDo);  
-      if(P->Call.out[MinOut]) fprintf(stderr,"(%i) Reading from file...\n", P->Par.me);
+      if(P->Call.out[MinOut]) fprintf(stderr, "(%i) Re-initializing %s psi array from source file of recent calculation\n", P->Par.me, R->MinimisationName[R->CurrentMin]);
       ReadSrcPsiDensity(P, Occupied, 0, R->LevSNo);
-      ResetGramSchTagType(P, Psi, Occupied, IsOrthonormal); // loaded values are orthonormal
+      //ResetGramSchTagType(P, Psi, Occupied, IsOrthonormal); // loaded values are orthonormal
+      // note: this did not work and is currently not clear why not (as TestGramSch says: OK, but minimisation goes awry without the following GramSch)
     }
+    SpeedMeasure(P, InitGramSchTime, StartTimeDo);  
+    GramSch(P, R->LevS, Psi, Orthonormalize);
+    SpeedMeasure(P, InitGramSchTime, StopTimeDo);  
     SpeedMeasure(P, InitDensityTime, StartTimeDo);
     InitDensityCalculation(P);
@@ -770,5 +772,5 @@
     R->LevS->Step++;
     EnergyOutput(P,0);
-  }
+  } 
   if (P->Call.ReadSrcFiles != 1) {  // otherwise minimise oneself
     if(P->Call.out[LeaderOut]) fprintf(stderr,"(%i)Beginning minimisation of type %s ...\n", P->Par.me, R->MinimisationName[Occupied]);
@@ -783,8 +785,10 @@
       if (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent != 1) {
         //SetArrayToDouble0((double *)LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo],LevS->MaxG*2);
-        memcpy(LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo], LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], ElementSize*LevS->MaxG*sizeof(double)); // restore old Psi from OldPsi
-        //fprintf(stderr,"(%i) Psi %i at %p stored in OldPsi at %p: Old[0] %lg+i%lg\n", P->Par.me, R->ActualLocalPsiNo, LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo], LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo][0].re, LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo][0].im);
-        f_m = P->Lat.E->TotalEnergy[0]; // grab first value
-        m = 0.;
+        if (R->DoBrent == 1) {
+          memcpy(LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo], LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], ElementSize*LevS->MaxG*sizeof(double)); // restore old Psi from OldPsi
+          //fprintf(stderr,"(%i) Psi %i at %p stored in OldPsi at %p: Old[0] %lg+i%lg\n", P->Par.me, R->ActualLocalPsiNo, LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo], LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo][0].re, LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo][0].im);
+          f_m = P->Lat.E->TotalEnergy[0]; // grab first value
+          m = 0.;
+        }
         CalculateNewWave(P,NULL);
         if ((R->DoBrent == 1) && (fabs(Lat->E->delta[0]) < M_PI/4.)) 
@@ -798,5 +802,5 @@
             // if we stepped on to a new Psi, which is already down at DoBrent=1 unlike the last one,
             // then an up-to-date gradient is missing for the following Brent line search
-            if(P->Call.out[LeaderOut]) fprintf(stderr,"(%i) We stepped on to a new Psi, which is already in the Brent regime ...re-calc delta\n", P->Par.me);
+            if(P->Call.out[MinOut]) fprintf(stderr,"(%i) We stepped on to a new Psi, which is already in the Brent regime ...re-calc delta\n", P->Par.me);
             memcpy(LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo], LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], ElementSize*LevS->MaxG*sizeof(double)); // restore old Psi from OldPsi
             //fprintf(stderr,"(%i) Psi %i at %p stored in OldPsi at %p: Old[0] %lg+i%lg\n", P->Par.me, R->ActualLocalPsiNo, LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo], LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo][0].re, LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo][0].im);
