source: pcp/src/run.c@ ccd028

/** \file run.c
 * Initialization of levels and calculation super-functions.
 * 
 * Most important functions herein are CalculateForce() and CalculateMD(), which calls various
 * functions in order to perfom the Molecular Dynamics simulation. MinimiseOccupied() and MinimiseUnOccupied()
 * call various functions to perform the actual minimisation for the occupied and unoccupied wave functions.
 * CalculateMinimumStop() evaluates the stop condition for desired precision or step count (or external signals).
 * 
 * Minor functions are ChangeToLevUp() (which takes the calculation to the next finer level),
 * UpdateActualPsiNo() (next Psi is minimized and an additional orthonormalization takes place) and UpdateToNewWaves() 
 * (which reinitializes density calculations after the wave functions have changed due to the ionic motion).
 * OccupyByFermi() re-occupies orbitals according to Fermi distribution if calculated with additional orbitals.
 * InitRun() and InitRunLevel() prepare the RunStruct with starting values. UpdateIon_PRCG() implements a CG
 * algorithm to minimize the ionic forces and thus optimize the structure. 
 * 
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers
 \date 2000

  File: run.c
  $Id: run.c,v 1.101.2.2 2007-04-21 13:01:13 foo Exp $
*/

#include <signal.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_randist.h>
#include "mpi.h"
#include "data.h"
#include "errors.h"
#include "helpers.h"
#include "init.h"
#include "opt.h"
#include "myfft.h"
#include "gramsch.h"
#include "output.h"
#include "energy.h"
#include "density.h"
#include "ions.h"
#include "run.h"
#include "riemann.h"
#include "mymath.h"
#include "pcp.h"
#include "perturbed.h"
#include "wannier.h"


/** Initialization of the (initial) zero and simulation levels in RunStruct structure.
 * RunStruct::InitLevS is set onto the STANDARTLEVEL in Lattice::Lev[], RunStruct::InitLev0 on
 * level 0, RunStruct::LevS onto Lattice::MaxLevel-1 (maximum level) and RunStruct::Lev0 onto
 * Lattice::MaxLevel-2 (one below).
 * In case of RiemannTensor use an additional Riemann level is intertwined.
 * \param *P Problem at hand
 */
void InitRunLevel(struct Problem *P)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct RiemannTensor *RT = &Lat->RT;
  int d,i;

  switch (Lat->RT.Use) {
  case UseNotRT:
    R->InitLevSNo = STANDARTLEVEL;
    R->InitLev0No = 0;
    R->InitLevS = &P->Lat.Lev[R->InitLevSNo];
    R->InitLev0 = &P->Lat.Lev[R->InitLev0No];
    R->LevSNo = Lat->MaxLevel-1;
    R->Lev0No = Lat->MaxLevel-2;
    R->LevS = &P->Lat.Lev[R->LevSNo];
    R->Lev0 = &P->Lat.Lev[R->Lev0No];
    break;
  case UseRT:
    R->InitLevSNo = STANDARTLEVEL;
    R->InitLev0No = 0;
    R->InitLevS = &P->Lat.Lev[R->InitLevSNo];
    R->InitLev0 = &P->Lat.Lev[R->InitLev0No];

    /*    R->LevSNo = Lat->MaxLevel-1;
          R->Lev0No = Lat->MaxLevel-2;*/
    R->LevSNo = Lat->MaxLevel-2;
    R->Lev0No = Lat->MaxLevel-3;

    R->LevRNo = P->Lat.RT.RiemannLevel;
    R->LevRSNo = STANDARTLEVEL;
    R->LevR0No = 0;
    R->LevS = &P->Lat.Lev[R->LevSNo];
    R->Lev0 = &P->Lat.Lev[R->Lev0No];
    R->LevR = &P->Lat.Lev[R->LevRNo];
    R->LevRS = &P->Lat.Lev[R->LevRSNo];
    R->LevR0 = &P->Lat.Lev[R->LevR0No];
    for (d=0; d<NDIM; d++) {
      RT->NUpLevRS[d] = 1;
      for (i=R->LevRNo-1; i >=  R->LevRSNo; i--) 
                                RT->NUpLevRS[d] *= Lat->LevelSizes[i];
      RT->NUpLevR0[d] = 1;
      for (i=R->LevRNo-1; i >=  R->LevR0No; i--) 
                                RT->NUpLevR0[d] *= Lat->LevelSizes[i];
    }
    break;
  }
}


/** Initialization of RunStruct structure.
 * Most of the actual entries in the RunStruct are set to their starter no-nonsense
 * values (init if LatticeLevel is not STANDARTLEVEL otherwise normal max): FactorDensity,
 * all Steps, XCEnergyFactor and HGcFactor, current and archived energie values are zeroed.
 * \param *P problem at hand
 */
void InitRun(struct Problem *P)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct Psis *Psi = &Lat->Psi;
  int i,j;

#ifndef SHORTSPEED
  R->MaxMinStepFactor = Psi->AllMaxLocalNo;
#else
  R->MaxMinStepFactor = SHORTSPEED;
#endif
  if (R->LevSNo == STANDARTLEVEL) {
    R->ActualMaxMinStep = R->MaxMinStep*R->MaxPsiStep*R->MaxMinStepFactor; 
    R->ActualRelEpsTotalEnergy = R->RelEpsTotalEnergy;
    R->ActualRelEpsKineticEnergy = R->RelEpsKineticEnergy;
    R->ActualMaxMinStopStep = R->MaxMinStopStep;
    R->ActualMaxMinGapStopStep = R->MaxMinGapStopStep;
  } else {
    R->ActualMaxMinStep = R->MaxInitMinStep*R->MaxPsiStep*R->MaxMinStepFactor;
    R->ActualRelEpsTotalEnergy = R->InitRelEpsTotalEnergy;
    R->ActualRelEpsKineticEnergy = R->InitRelEpsKineticEnergy;
    R->ActualMaxMinStopStep = R->InitMaxMinStopStep;
    R->ActualMaxMinGapStopStep = R->InitMaxMinGapStopStep;
  }
  
  R->FactorDensityR = 1./Lat->Volume;
  R->FactorDensityC = Lat->Volume;
  
  R->OldActualLocalPsiNo = R->ActualLocalPsiNo = 0;
        R->UseForcesFile = 0;
  R->UseOldPsi = 1;
  R->MinStep = 0;
  R->PsiStep = 0;
  R->AlphaStep = 0;
  R->DoCalcCGGauss = 0;
  R->CurrentMin = Occupied;

  R->MinStopStep = 0;
  
  R->ScanPotential = 0; // in order to deactivate, simply set this to 0
  R->ScanAtStep = 6;    // must not be set to same as ScanPotential (then gradient is never calculated)
  R->ScanDelta = 0.01;  // step size on advance
  R->ScanFlag = 0;      // flag telling that we are scanning
  
  //R->DoBrent = 0;   // InitRun() occurs after ReadParameters(), thus this deactivates DoBrent line search
  
  /*  R->MaxOuterStep = 1;
      R->MeanForceEps = 0.0;*/
  
  R->NewRStep = 1;
  /* Factor */
  R->XCEnergyFactor = 1.0/R->FactorDensityR;
  R->HGcFactor = 1.0/Lat->Volume;

  /* Sollte auch geaendert werden */
  /*Grad->GradientArray[GraSchGradient] = LevS->LPsi->LocalPsi[Psi->LocalNo];*/

  for (j=Occupied;j<Extra;j++)
    for (i=0; i < RUNMAXOLD; i++) {
      R->TE[j][i] = 0;
      R->KE[j][i] = 0;
    }
    
  R->MinimisationName = (char **) Malloc((perturbations+3)*(sizeof(char *)), "InitRun: *MinimisationName");
  for (j=Occupied;j<=Extra;j++)
    R->MinimisationName[j] = (char *) MallocString(6*(sizeof(char)), "InitRun: MinimisationName[]");
  strncpy(R->MinimisationName[0],"Occ",6); 
  strncpy(R->MinimisationName[1],"UnOcc",6); 
  strncpy(R->MinimisationName[2],"P0",6); 
  strncpy(R->MinimisationName[3],"P1",6); 
  strncpy(R->MinimisationName[4],"P2",6); 
  strncpy(R->MinimisationName[5],"RxP0",6); 
  strncpy(R->MinimisationName[6],"RxP1",6); 
  strncpy(R->MinimisationName[7],"RxP2",6); 
  strncpy(R->MinimisationName[8],"Extra",6); 
}

/** Re-occupy orbitals according to Fermi (bottom-up energy-wise).
 * All OnePsiElementAddData#PsiFactor's are set to zero. \a electrons is set to Psi#Use-dependent
 * Psis#GlobalNo.
 * Then we go through OnePsiElementAddData#Lambda, find biggest, put one or two electrons into
 * its PsiFactor, withdraw from \a electrons. Go on as long as there are \a electrons left.
 * \param *P Problem at hand
 */
void OccupyByFermi(struct Problem *P) {
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  int i, index, electrons = 0;
  double lambda, electronsperorbit;
  
  for (i=0; i< Psi->LocalNo; i++) {// set all PsiFactors to zero
    Psi->LocalPsiStatus[i].PsiFactor = 0.0;
    Psi->LocalPsiStatus[i].PsiType = UnOccupied;
    //Psi->LocalPsiStatus[i].PsiGramSchStatus = (R->DoSeparated) ? NotUsedToOrtho : NotOrthogonal;
  }

  electronsperorbit = (Psi->Use == UseSpinUpDown) ? 1 : 2;
  switch (Psi->PsiST) { // how many electrons may we re-distribute
    case SpinDouble:
      electrons = Psi->GlobalNo[PsiMaxNoDouble];
      break;
    case SpinUp:
      electrons = Psi->GlobalNo[PsiMaxNoUp];
      break;
    case SpinDown:
      electrons = Psi->GlobalNo[PsiMaxNoDown];
      break;
  } 
  while (electrons > 0) {
    lambda = 0.0;
    index = 0;
    for (i=0; i< Psi->LocalNo; i++) // seek biggest unoccupied one
      if ((lambda < Psi->AddData[i].Lambda) && (Psi->LocalPsiStatus[i].PsiFactor == 0.0)) { 
        index = i;
        lambda = Psi->AddData[i].Lambda;
      }
    Psi->LocalPsiStatus[index].PsiFactor = electronsperorbit; // occupy state
    Psi->LocalPsiStatus[index].PsiType = Occupied;    
    electrons--;   // one electron less
  }
  for (i=0; i< Psi->LocalNo; i++) // set all PsiFactors to zero
    if (Psi->LocalPsiStatus[i].PsiType == UnOccupied) Psi->LocalPsiStatus[i].PsiFactor = 1.0;
    
  SpeedMeasure(P, DensityTime, StartTimeDo);
  UpdateDensityCalculation(P);
  SpeedMeasure(P, DensityTime, StopTimeDo);
  InitPsiEnergyCalculation(P,Occupied);  // goes through all orbitals calculating kinetic and non-local
  CalculateDensityEnergy(P, 0);
  EnergyAllReduce(P);
//  SetCurrentMinState(P,UnOccupied);
//  InitPsiEnergyCalculation(P,UnOccupied); /* STANDARTLEVEL */
//  CalculateGapEnergy(P); /* STANDARTLEVEL */
//  EnergyAllReduce(P);
//  SetCurrentMinState(P,Occupied);
}

/** Use next local Psi: Update RunStruct::ActualLocalPsiNo.
 * Increases OnePsiElementAddData::Step, RunStruct::MinStep and RunStruct::PsiStep.
 * RunStruct::OldActualLocalPsiNo is set to current one and this distributed
 * (UpdateGramSchOldActualPsiNo()) among process.
 * Afterwards RunStruct::ActualLocalPsiNo is increased (modulo Psis::LocalNo of
 * this process) and again distributed (UpdateGramSchActualPsiNo()).
 * Due to change in the GramSchmidt-Status, GramSch() is called for Orthonormalization.
 * \param *P Problem at hand#
 * \param IncType skip types PsiTypeTag#UnOccupied or PsiTypeTag#Occupied we only want next(thus we can handily advance only through either type)
 */
void UpdateActualPsiNo(struct Problem *P, enum PsiTypeTag IncType)
{
  struct RunStruct *R = &P->R;
  if (R->CurrentMin != IncType) {
    SetCurrentMinState(P,IncType);
    R->PsiStep = R->MaxPsiStep; // force step to next Psi
  }
  P->Lat.Psi.AddData[R->ActualLocalPsiNo].Step++;
  R->MinStep++;
  R->PsiStep++;
  if (R->OldActualLocalPsiNo != R->ActualLocalPsiNo) {
    R->OldActualLocalPsiNo = R->ActualLocalPsiNo;
    UpdateGramSchOldActualPsiNo(P, &P->Lat.Psi);
  }
  if (R->PsiStep >= R->MaxPsiStep) {
    R->PsiStep=0;
    do { 
      R->ActualLocalPsiNo++;
      R->ActualLocalPsiNo %= P->Lat.Psi.LocalNo;
    } while (P->Lat.Psi.AllPsiStatus[R->ActualLocalPsiNo].PsiType != IncType);
    UpdateGramSchActualPsiNo(P, &P->Lat.Psi);
    //fprintf(stderr,"(%i) ActualLocalNo: %d\n", P->Par.me, R->ActualLocalPsiNo);
  }
  if ((R->UseAddGramSch == 1 && (R->OldActualLocalPsiNo != R->ActualLocalPsiNo || P->Lat.Psi.NoOfPsis == 1)) || R->UseAddGramSch == 2) {
    if (P->Lat.Psi.LocalPsiStatus[R->OldActualLocalPsiNo].PsiGramSchStatus != NotUsedToOrtho) // don't reset by accident last psi of former minimisation run 
      SetGramSchOldActualPsi(P, &P->Lat.Psi, NotOrthogonal); 
    SpeedMeasure(P, GramSchTime, StartTimeDo);
    //OrthogonalizePsis(P);
    if (R->CurrentMin <= UnOccupied)
      GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
    else
      GramSch(P, R->LevS, &P->Lat.Psi, Orthogonalize); //Orthogonalize
    SpeedMeasure(P, GramSchTime, StopTimeDo);  
  }
}

/** Resets all OnePsiElement#DoBrent.\
 * \param *P Problem at hand
 * \param *Psi pointer to wave functions
 */
void ResetBrent(struct Problem *P, struct Psis *Psi) {
  int i;
  for (i=0; i< Psi->LocalNo; i++) {// set all PsiFactors to zero
    //fprintf(stderr,"(%i) DoBrent[%i] = %i\n", P->Par.me, i, Psi->LocalPsiStatus[i].DoBrent);
    if (Psi->LocalPsiStatus[i].PsiType == Occupied) Psi->LocalPsiStatus[i].DoBrent = 4;
  }
}

/** Sets current minimisation state.
 * Stores given \a state in RunStruct#CurrentMin and sets pointer Lattice#E accordingly.
 * \param *P Problem at hand
 * \param state given PsiTypeTag state
 */
void SetCurrentMinState(struct Problem *P, enum PsiTypeTag state) {
  P->R.CurrentMin = state;
  P->R.TotalEnergy = &(P->R.TE[state][0]);
  P->R.KineticEnergy = &(P->R.KE[state][0]);
  P->R.ActualRelTotalEnergy = &(P->R.ActualRelTE[state][0]);
  P->R.ActualRelKineticEnergy = &(P->R.ActualRelKE[state][0]);
  P->Lat.E = &(P->Lat.Energy[state]);  
}
/*{  
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  P->Lat.Psi.AddData[R->ActualLocalPsiNo].Step++;
  R->MinStep++;
  R->PsiStep++;
  if (R->OldActualLocalPsiNo != R->ActualLocalPsiNo) {  // remember old actual local number
    R->OldActualLocalPsiNo = R->ActualLocalPsiNo;
    UpdateGramSchOldActualPsiNo(P, &P->Lat.Psi);
  }
  if (R->PsiStep >= R->MaxPsiStep) {  // done enough minimisation steps for this orbital?
    R->PsiStep=0;
    do {  // step on as long as we are still on a SkipType orbital
      R->ActualLocalPsiNo++;
      R->ActualLocalPsiNo %= P->Lat.Psi.LocalNo;
    } while ((P->Lat.Psi.LocalPsiStatus[R->ActualLocalPsiNo].PsiType == SkipType));
    UpdateGramSchActualPsiNo(P, &P->Lat.Psi);
    if (R->UseAddGramSch >= 1) {
      SetGramSchOldActualPsi(P,Psi,NotOrthogonal);
      // setze von OldActual bis bla auf nicht orthogonal
      GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
    }
  } else if (R->UseAddGramSch == 2) {
    SetGramSchOldActualPsi(P, &P->Lat.Psi, NotOrthogonal); 
    //if (SkipType == UnOccupied)
      //ResetGramSch(P,Psi);
    //fprintf(stderr,"UpdateActualPsiNo: GramSch() for %i\n",R->OldActualLocalPsiNo);
    GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
  }
}*/

/** Upgrades the calculation to the next finer level.
 * If we are below the initial level, 
 * ChangePsiAndDensToLevUp() prepares density and Psi coefficients.
 * Then the level change is made as RunStruct::LevSNo and RunStruct::Lev0No are decreased.
 * The RunStruct::OldActualLocalPsi is set to current one and both are distributed
 * (UpdateGramSchActualPsiNo(), UpdateGramSchOldActualPsiNo()).
 * The PseudoPot'entials adopt the level up by calling ChangePseudoToLevUp().
 * Now we are prepared to reset Energy::PsiEnergy and local and total density energy and
 * recalculate them: InitPsiEnergyCalculation(), CalculateDensityEnergy() and CalculateIonsEnergy().
 * Results are gathered EnergyAllReduce() and the output made EnergyOutput().
 * Finally, the stop condition are reset for the new level (depending if it's the STANDARTLEVEL or
 * not).
 * \param *P Problem at hand
 * \param *Stop is set to zero if we are below or equal to init level (see CalculateForce())
 * \sa UpdateToNewWaves() very similar in the procedure, only the update of the Psis and density
 *     (ChangePsiAndDensToLevUp()) is already made there.
 * \bug Missing TotalEnergy shifting for other PsiTypeTag's!
 */
static void ChangeToLevUp(struct Problem *P, int *Stop)
{
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct Energy *E = Lat->E;
  struct RiemannTensor *RT = &Lat->RT;
  int i;
  if (R->LevSNo <= R->InitLevSNo) {
    if (P->Call.out[LeaderOut] && (P->Par.me == 0))
      fprintf(stderr, "(%i) ChangeLevUp: LevSNo(%i) <= InitLevSNo(%i)\n", P->Par.me, R->LevSNo, R->InitLevSNo);
    *Stop = 1;
    return;
  }
  if (P->Call.out[LeaderOut] && (P->Par.me == 0))
    fprintf(stderr, "(0) ChangeLevUp: LevSNo(%i) InitLevSNo(%i)\n", R->LevSNo, R->InitLevSNo); 
  *Stop = 0;
  P->Speed.LevUpSteps++;
  SpeedMeasure(P, SimTime, StopTimeDo);
  SpeedMeasure(P, InitSimTime, StartTimeDo);
  SpeedMeasure(P, InitDensityTime, StartTimeDo);
  ChangePsiAndDensToLevUp(P);
  SpeedMeasure(P, InitDensityTime, StopTimeDo);
  R->LevSNo--;
  R->Lev0No--;
  if (RT->ActualUse == standby && R->LevSNo == STANDARTLEVEL) {
    P->Lat.RT.ActualUse = active;
    CalculateRiemannTensorData(P);
    Error(SomeError, "Calculate RT: Not further implemented");
  }
  R->LevS = &P->Lat.Lev[R->LevSNo];
  R->Lev0 = &P->Lat.Lev[R->Lev0No];
  R->OldActualLocalPsiNo = R->ActualLocalPsiNo;
  UpdateGramSchActualPsiNo(P, &P->Lat.Psi);
  UpdateGramSchOldActualPsiNo(P, &P->Lat.Psi);
  ResetBrent(P, &P->Lat.Psi);
  R->PsiStep=0;
  R->MinStep=0;
  P->Grad.GradientArray[GraSchGradient] = R->LevS->LPsi->LocalPsi[Psi->LocalNo];
  ChangePseudoToLevUp(P);
  for (i=0; i<MAXALLPSIENERGY; i++) 
    SetArrayToDouble0(E->PsiEnergy[i], Psi->LocalNo);
  SetArrayToDouble0(E->AllLocalDensityEnergy, MAXALLDENSITYENERGY);
  SetArrayToDouble0(E->AllTotalDensityEnergy, MAXALLDENSITYENERGY);
  for (i=MAXOLD-1; i > 0; i--) {
    E->TotalEnergy[i] = E->TotalEnergy[i-1];
    Lat->Energy[UnOccupied].TotalEnergy[i] = Lat->Energy[UnOccupied].TotalEnergy[i-1];
  }
  InitPsiEnergyCalculation(P,Occupied);
  CalculateDensityEnergy(P,1); 
  CalculateIonsEnergy(P);
  EnergyAllReduce(P);
/*  SetCurrentMinState(P,UnOccupied);
  InitPsiEnergyCalculation(P,UnOccupied);
  CalculateGapEnergy(P);
  EnergyAllReduce(P);
  SetCurrentMinState(P,Occupied);*/
  EnergyOutput(P,0);
  if (R->LevSNo == STANDARTLEVEL) {
    R->ActualMaxMinStep = R->MaxMinStep*R->MaxPsiStep*R->MaxMinStepFactor; 
    R->ActualRelEpsTotalEnergy = R->RelEpsTotalEnergy;
    R->ActualRelEpsKineticEnergy = R->RelEpsKineticEnergy;
    R->ActualMaxMinStopStep = R->MaxMinStopStep;
    R->ActualMaxMinGapStopStep = R->MaxMinGapStopStep;
  } else {
    R->ActualMaxMinStep = R->MaxInitMinStep*R->MaxPsiStep*R->MaxMinStepFactor;
    R->ActualRelEpsTotalEnergy = R->InitRelEpsTotalEnergy;
    R->ActualRelEpsKineticEnergy = R->InitRelEpsKineticEnergy;
    R->ActualMaxMinStopStep = R->InitMaxMinStopStep; 
    R->ActualMaxMinGapStopStep = R->InitMaxMinGapStopStep; 
  }
  R->MinStopStep = 0;
  SpeedMeasure(P, InitSimTime, StopTimeDo);
  SpeedMeasure(P, SimTime, StartTimeDo);
  if (P->Call.out[LeaderOut] && (P->Par.me == 0))
    fprintf(stderr, "(0) ChangeLevUp: LevSNo(%i) InitLevSNo(%i) Done\n", R->LevSNo, R->InitLevSNo); 
}

/** Updates after the wave functions have changed (e.g.\ Ion moved).
 * Old and current RunStruct::ActualLocalPsiNo are zeroed and distributed among the processes.
 * InitDensityCalculation() is called, afterwards the pseudo potentials update to the new
 * wave functions UpdatePseudoToNewWaves().
 * Energy::AllLocalDensityEnergy, Energy::AllTotalDensityEnergy, Energy::AllTotalIonsEnergy and
 * Energy::PsiEnergy[i] are set to zero.
 * We are set to recalculate all of the following energies: Psis InitPsiEnergyCalculation(), density
 * CalculateDensityEnergy(), ionic CalculateIonsEnergy() and ewald CalculateEwald().
 * Results are gathered from all processes EnergyAllReduce() and EnergyOutput() is called.
 * Finally, the various conditons in the RunStruct for stopping the calculation are set: number of
 * minimisation steps, relative total or kinetic energy change or how often stop condition was
 * evaluated.
 * \param *P Problem at hand
 */
static void UpdateToNewWaves(struct Problem *P)
{
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct Energy *E = Lat->E;
  int i;//,type;
  R->OldActualLocalPsiNo = R->ActualLocalPsiNo = 0;
  //if (isnan((double)R->LevS->LPsi->LocalPsi[R->OldActualLocalPsiNo][0].re)) { fprintf(stderr,"(%i) WARNING in UpdateGramSchActualPsiNo(): LPsi->LocalPsi[%i]_[%i] = NaN!\n", P->Par.me, R->OldActualLocalPsiNo, 0); Error(SomeError, "NaN-Fehler!"); }
  UpdateGramSchActualPsiNo(P, &P->Lat.Psi);
  UpdateGramSchOldActualPsiNo(P, &P->Lat.Psi);
  R->PsiStep=0;
  R->MinStep=0;
  SpeedMeasure(P, InitDensityTime, StartTimeDo);
  //if (isnan((double)R->LevS->LPsi->LocalPsi[R->OldActualLocalPsiNo][0].re)) { fprintf(stderr,"(%i) WARNING in Update../InitDensityCalculation(): LPsi->LocalPsi[%i]_[%i] = NaN!\n", P->Par.me, R->OldActualLocalPsiNo, 0); Error(SomeError, "NaN-Fehler!"); }
  InitDensityCalculation(P);
  SpeedMeasure(P, InitDensityTime, StopTimeDo);
  UpdatePseudoToNewWaves(P);
  for (i=0; i<MAXALLPSIENERGY; i++) 
    SetArrayToDouble0(E->PsiEnergy[i], Psi->LocalNo);
  SetArrayToDouble0(E->AllLocalDensityEnergy, MAXALLDENSITYENERGY);
  SetArrayToDouble0(E->AllTotalDensityEnergy, MAXALLDENSITYENERGY);
  SetArrayToDouble0(E->AllTotalIonsEnergy, MAXALLIONSENERGY);
  InitPsiEnergyCalculation(P,Occupied);
  CalculateDensityEnergy(P,1); 
  CalculateIonsEnergy(P);
  CalculateEwald(P, 0);
  EnergyAllReduce(P);
/*  if (R->DoUnOccupied) {
          SetCurrentMinState(P,UnOccupied);
          InitPsiEnergyCalculation(P,UnOccupied); 
          CalculateGapEnergy(P);        
          EnergyAllReduce(P);
  }
  if (R->DoPerturbation)
    for(type=Perturbed_P0;type <=Perturbed_RxP2;type++) {
      SetCurrentMinState(P,type);
      InitPerturbedEnergyCalculation(P,1);
      EnergyAllReduce(P);
    }
  SetCurrentMinState(P,Occupied);*/
  E->TotalEnergyOuter[0] = E->TotalEnergy[0];
  EnergyOutput(P,0);
  R->ActualMaxMinStep = R->MaxMinStep*R->MaxPsiStep*R->MaxMinStepFactor; 
  R->ActualRelEpsTotalEnergy = R->RelEpsTotalEnergy;
  R->ActualRelEpsKineticEnergy = R->RelEpsKineticEnergy;
  R->ActualMaxMinStopStep = R->MaxMinStopStep;
  R->ActualMaxMinGapStopStep = R->MaxMinGapStopStep;
  R->MinStopStep = 0;
}

/** Evaluates the stop condition and returns 0 or 1 for occupied states.
 * Stop is set when:
 * - SuperStop at best possible point (e.g.\ LevelChange): RunStruct::PsiStep == 0 && SuperStop == 1
 * - RunStruct::PsiStep && RunStruct::MinStopStep modulo RunStruct::ActualMaxMinStopStep == 0 
 *    - To many minimisation steps: RunStruct::MinStep > RunStruct::ActualMaxMinStopStep
 *    - below relative rate of change:
 *        - Remember old values: Shift all RunStruct::TotalEnergy and RunStruct::KineticEnergy by
 *          one and transfer current one from Energy::TotalEnergy and Energy::AllTotalPsiEnergy[KineticEnergy].
 *        - if more than one minimisation step was made, calculate the relative changes of total
 *          energy and kinetic energy and store them in RunStruct::ActualRelTotalEnergy and
 *          RunStruct::ActualRelKineticEnergy and check them against the sought for minimum
 *          values RunStruct::ActualRelEpsTotalEnergy and RunStruct::ActualRelEpsKineticEnergy.
 * - if RunStruct::PsiStep is zero (default), increase RunStruct::MinStopStep
 * \param *P Problem at hand
 * \param SuperStop 1 - external signal: ceasing calculation, 0 - no signal
 * \return Stop: 1 - stop, 0 - continue
 */
int CalculateMinimumStop(struct Problem *P, int SuperStop)
{
  int Stop = 0, i;
  struct RunStruct *R = &P->R;
  struct Energy *E = P->Lat.E;
  if (R->PsiStep == 0 && SuperStop) Stop = 1;
  if (R->PsiStep == 0 && ((R->MinStopStep % R->ActualMaxMinStopStep == 0 && R->CurrentMin != UnOccupied) || (R->MinStopStep % R->ActualMaxMinGapStopStep == 0 && R->CurrentMin == UnOccupied))) {
    if (R->MinStep >= R->ActualMaxMinStep) Stop = 1;
    for (i=RUNMAXOLD-1; i > 0; i--) {
      R->TotalEnergy[i] = R->TotalEnergy[i-1];
      R->KineticEnergy[i] = R->KineticEnergy[i-1];
    }
    R->TotalEnergy[0] = E->TotalEnergy[0];
    R->KineticEnergy[0] = E->AllTotalPsiEnergy[KineticEnergy];
    if (R->MinStopStep) {
      //if (R->TotalEnergy[1] < MYEPSILON) fprintf(stderr,"CalculateMinimumStop: R->TotalEnergy[1] = %lg\n",R->TotalEnergy[1]);
      R->ActualRelTotalEnergy[0] = fabs((R->TotalEnergy[0]-R->TotalEnergy[1])/R->TotalEnergy[1]);
      //if (R->KineticEnergy[1] < MYEPSILON) fprintf(stderr,"CalculateMinimumStop: R->KineticEnergy[1] = %lg\n",R->KineticEnergy[1]);
      //if (R->CurrentMin < Perturbed_P0) 
      R->ActualRelKineticEnergy[0] = fabs((R->KineticEnergy[0]-R->KineticEnergy[1])/R->KineticEnergy[1]);
      //else R->ActualRelKineticEnergy[0] = 0.;
      if (P->Call.out[LeaderOut] && (P->Par.me == 0))
      switch (R->CurrentMin) {
        default:
                fprintf(stderr, "ARelTE: %e\tARelKE: %e\n", R->ActualRelTotalEnergy[0], R->ActualRelKineticEnergy[0]); 
          break;
        case UnOccupied:
          fprintf(stderr, "ARelTGE: %e\tARelKGE: %e\n", R->ActualRelTotalEnergy[0], R->ActualRelKineticEnergy[0]); 
          break;
      }
      //fprintf(stderr, "(%i) Comparing TE: %lg < %lg\tKE: %lg < %lg?\n", P->Par.me, R->ActualRelTotalEnergy[0], R->ActualRelEpsTotalEnergy, R->ActualRelKineticEnergy[0], R->ActualRelEpsKineticEnergy);
      if ((R->ActualRelTotalEnergy[0] < R->ActualRelEpsTotalEnergy) &&
          (R->ActualRelKineticEnergy[0]  < R->ActualRelEpsKineticEnergy))
        Stop = 1;
    }
  }
  if (R->PsiStep == 0)
    R->MinStopStep++;
  if (P->Call.WriteSrcFiles == 2)
    OutputVisSrcFiles(P, R->CurrentMin);
  return(Stop);
}

/** Evaluates the stop condition and returns 0 or 1 for gap energies.
 * Stop is set when:
 * - SuperStop at best possible point (e.g.\ LevelChange): RunStruct::PsiStep == 0 && SuperStop == 1
 * - RunStruct::PsiStep && RunStruct::MinStopStep modulo RunStruct::ActualMaxMinStopStep == 0 
 *    - To many minimisation steps: RunStruct::MinStep > RunStruct::ActualMaxMinStopStep
 *    - below relative rate of change:
 *        - Remember old values: Shift all RunStruct::TotalEnergy and RunStruct::KineticEnergy by
 *          one and transfer current one from Energy::TotalEnergy and Energy::AllTotalPsiEnergy[KineticEnergy].
 *        - if more than one minimisation step was made, calculate the relative changes of total
 *          energy and kinetic energy and store them in RunStruct::ActualRelTotalEnergy and
 *          RunStruct::ActualRelKineticEnergy and check them against the sought for minimum
 *          values RunStruct::ActualRelEpsTotalEnergy and RunStruct::ActualRelEpsKineticEnergy.
 * - if RunStruct::PsiStep is zero (default), increase RunStruct::MinStopStep
 * \param *P Problem at hand
 * \param SuperStop 1 - external signal: ceasing calculation, 0 - no signal
 * \return Stop: 1 - stop, 0 - continue
 * \sa CalculateMinimumStop() - same procedure for occupied states
 *//*
static double CalculateGapStop(struct Problem *P, int SuperStop)
{
  int Stop = 0, i;
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = P->Lat.E;
  if (R->PsiStep == 0 && SuperStop) Stop = 1;
  if (R->PsiStep == 0 && (R->MinStopStep % R->ActualMaxMinGapStopStep) == 0) {
    if (R->MinStep >= R->ActualMaxMinStep) Stop = 1;
    for (i=RUNMAXOLD-1; i > 0; i--) {
      R->TotalGapEnergy[i] = R->TotalGapEnergy[i-1];
      R->KineticGapEnergy[i] = R->KineticGapEnergy[i-1];
    }
    R->TotalGapEnergy[0] = Lat->Energy[UnOccupied].TotalEnergy[0];
    R->KineticGapEnergy[0] = E->AllTotalPsiEnergy[GapPsiEnergy];
    if (R->MinStopStep) {
      if (R->TotalGapEnergy[1] < MYEPSILON) fprintf(stderr,"CalculateMinimumStop: R->TotalGapEnergy[1] = %lg\n",R->TotalGapEnergy[1]);
      R->ActualRelTotalGapEnergy[0] = fabs((R->TotalGapEnergy[0]-R->TotalGapEnergy[1])/R->TotalGapEnergy[1]);
      if (R->KineticGapEnergy[1] < MYEPSILON) fprintf(stderr,"CalculateMinimumStop: R->KineticGapEnergy[1] = %lg\n",R->KineticGapEnergy[1]);
      R->ActualRelKineticGapEnergy[0] = fabs((R->KineticGapEnergy[0]-R->KineticGapEnergy[1])/R->KineticGapEnergy[1]);
      if (P->Call.out[LeaderOut] && (P->Par.me == 0))
        fprintf(stderr, "(%i) -------------------------> ARelTGE: %e\tARelKGE: %e\n", P->Par.me, R->ActualRelTotalGapEnergy[0], R->ActualRelKineticGapEnergy[0]); 
      if ((R->ActualRelTotalGapEnergy[0] < R->ActualRelEpsTotalGapEnergy) &&
          (R->ActualRelKineticGapEnergy[0]  < R->ActualRelEpsKineticGapEnergy))
        Stop = 1;
    }
  }
  if (R->PsiStep == 0)
    R->MinStopStep++;

  return(Stop);
}*/

#define StepTolerance 1e-4

static void CalculateEnergy(struct Problem *P) {
  SpeedMeasure(P, DensityTime, StartTimeDo);
  UpdateDensityCalculation(P);
  SpeedMeasure(P, DensityTime, StopTimeDo);
  UpdatePsiEnergyCalculation(P);
  CalculateDensityEnergy(P, 0);
  //CalculateIonsEnergy(P);
  EnergyAllReduce(P);
}

/** Energy functional depending on one parameter \a x (for a certain Psi in a certain conjugate direction).
 * \param x parameter for the which the function must be minimized
 * \param *params additional params
 * \return total energy if Psi is changed according to the given parameter
 */
static double fn1 (double x, void * params) {
  struct Problem *P = (struct Problem *)(params);
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *LevS = R->LevS;
  int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  int i=R->ActualLocalPsiNo;
  double ret;
  
  //fprintf(stderr,"(%i) Evaluating fnl at %lg ...\n",P->Par.me, x);
  //TestForOrth(P,R->LevS,P->Grad.GradientArray[GraSchGradient]);
  CalculateNewWave(P, &x); // also stores Psi to oldPsi
  //TestGramSch(P,R->LevS,&P->Lat.Psi,Occupied);
  //fprintf(stderr,"(%i) Testing for orthogonality of %i against ...\n",P->Par.me, R->ActualLocalPsiNo);
  //TestForOrth(P, LevS, LevS->LPsi->LocalPsi[R->ActualLocalPsiNo]);
  //UpdateActualPsiNo(P, Occupied);
  //UpdateEnergyArray(P);
  CalculateEnergy(P);
  ret = Lat->E->TotalEnergy[0];
  memcpy(LevS->LPsi->LocalPsi[i], LevS->LPsi->OldLocalPsi[i], ElementSize*LevS->MaxG*sizeof(double)); // restore old Psi from OldPsi
  //fprintf(stderr,"(%i) Psi %i at %p retrieved from 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);
  CalculateEnergy(P);
  //fprintf(stderr,"(%i) fnl(%lg) = %lg\n", P->Par.me, x, ret);
  return ret;
}

#ifdef HAVE_INLINE
inline void flip(double *a, double *b) {
#else
void flip(double *a, double *b) {
#endif
  double tmp = *a;
  *a = *b;
  *b = tmp;
}


/** Minimise PsiType#Occupied orbitals.
 * It is checked whether CallOptions#ReadSrcFiles is set and thus coefficients for the level have to be
 * read from file and afterwards initialized.
 * 
 * Then follows the main loop, until a stop condition is met:
 * -# CalculateNewWave()\n
 *    Over a conjugate gradient method the next (minimal) wave function is sought for.
 * -# UpdateActualPsiNo()\n
 *    Switch local Psi to next one.
 * -# UpdateEnergyArray()\n
 *    Shift archived energy values to make space for next one.
 * -# UpdateDensityCalculation(), SpeedMeasure()'d in DensityTime\n
 *    Calculate TotalLocalDensity of LocalPsis and gather results as TotalDensity.
 * -# UpdatePsiEnergyCalculation()\n
 *    Calculate kinetic and non-local energy contributons from the Psis.
 * -# CalculateDensityEnergy()\n
 *    Calculate remaining energy contributions from the Density and adds \f$V_xc\f$ onto DensityTypes#HGDensity.
 * -# CalculateIonsEnergy()\n
 *    Calculate the Gauss self energy of the Ions.
 * -# EnergyAllReduce()\n
 *    Gather PsiEnergy results from all processes and sum up together with all other contributions to TotalEnergy.
 * -# CheckCPULIM()\n
 *    Check if external signal has been received (e.g. end of time slit on cluster), break operation at next possible moment.
 * -# CalculateMinimumStop()\n
 *    Evaluates stop condition if desired precision or steps or ... have been reached. Otherwise go to
 *    CalculateNewWave().
 * 
 * Before return orthonormality is tested.
 * \param *P Problem at hand
 * \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)
{
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  //struct FileData *F = &P->Files;
//  int i;
//  double norm;
  //double dEdt0,ddEddt0,HartreeddEddt0,XCddEddt0, d[4], D[4],ConDirHConDir;
  struct LatticeLevel *LevS = R->LevS;
  int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  int iter = 0, status, max_iter=10;
  const gsl_min_fminimizer_type *T;
  gsl_min_fminimizer *s;
  double m, a, b;
  double f_m = 0., f_a, f_b;
  double dcos, dsin;
  int g;
  fftw_complex *ConDir = P->Grad.GradientArray[ConDirGradient];
  fftw_complex *source = NULL, *oldsource = NULL; 
  gsl_function F;
  F.function = &fn1;
  F.params = (void *) P;
  T = gsl_min_fminimizer_brent;
  s = gsl_min_fminimizer_alloc (T);
  int DoBrent, StartLocalPsiNo;
 
  ResetBrent(P,Psi);
  *Stop = 0;
  if (P->Call.ReadSrcFiles != DoNotParse) {
    if (!ReadSrcPsiDensity(P,Occupied,1, R->LevSNo)) {  // if file for level exists and desired, read from file
      P->Call.ReadSrcFiles = DoNotParse; // -r was bogus, remove it, have to start anew
      if(P->Call.out[MinOut]) fprintf(stderr,"(%i) Re-initializing, files are missing/corrupted...\n", P->Par.me);
      InitPsisValue(P, Psi->TypeStartIndex[Occupied], Psi->TypeStartIndex[Occupied+1]);  // initialize perturbed array for this run
      ResetGramSchTagType(P, Psi, Occupied, NotOrthogonal); // loaded values are orthonormal
    } else {
      SpeedMeasure(P, InitSimTime, StartTimeDo);  
      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
      // 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);
    SpeedMeasure(P, InitDensityTime, StopTimeDo);
    InitPsiEnergyCalculation(P, Occupied);  // go through all orbitals calculating kinetic and non-local
    StartLocalPsiNo = R->ActualLocalPsiNo;
    do { // otherwise OnePsiElementAddData#Lambda is calculated only for current Psi not for all
      CalculateDensityEnergy(P, 0);
      UpdateActualPsiNo(P, Occupied);
    } while (R->ActualLocalPsiNo != StartLocalPsiNo);
    CalculateIonsEnergy(P);
    EnergyAllReduce(P);
    SpeedMeasure(P, InitSimTime, StopTimeDo);  
    R->LevS->Step++;
    EnergyOutput(P,0);
  } 
  if ((P->Call.ReadSrcFiles != DoReadAllSrcDensities) && (P->Call.ReadSrcFiles != DoReadOccupiedSrcDensities)) {  // otherwise minimise oneself
    if(P->Call.out[LeaderOut]) fprintf(stderr,"(%i)Beginning minimisation of type %s ...\n", P->Par.me, R->MinimisationName[Occupied]);
    while (*Stop != 1) {  // loop testing condition over all Psis
      // in the following loop, we have two cases:
      // 1) still far away and just guessing: Use the normal CalculateNewWave() to improve Psi
      // 2) closer (DoBrent=-1): use brent line search instead
      // and due to these two cases, we also have two ifs inside each in order to catch stepping from one case
      // to the other - due to decreasing DoBrent and/or stepping to the next Psi (which may not yet be DoBrent==1)

      // case 1)
      if (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent != 1) {
        //SetArrayToDouble0((double *)LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo],LevS->MaxG*2);
        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.)) 
          Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent--;
        if (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent != 1) {
          UpdateActualPsiNo(P, Occupied);
          UpdateEnergyArray(P);
          CalculateEnergy(P); // just to get a sensible delta
          if ((R->ActualLocalPsiNo != R->OldActualLocalPsiNo) && (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent == 1)) {
             R->OldActualLocalPsiNo = R->ActualLocalPsiNo;
            // 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[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);
            f_m = P->Lat.E->TotalEnergy[0]; // grab first value
            m = 0.;
            DoBrent = Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent;
            Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent = 2;
            CalculateNewWave(P,NULL);
            Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent = DoBrent;
          }
          //fprintf(stderr,"(%i) fnl(%lg) = %lg\n", P->Par.me, m, f_m);
        } 
      }
      
      // case 2)
      if (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent == 1) {
        R->PsiStep=R->MaxPsiStep;  // no more fresh gradients from this point for current ActualLocalPsiNo
        a = b = 0.5*fabs(Lat->E->delta[0]);
        // we have a meaningful first minimum guess from above CalculateNewWave() resp. from end of this if of last step: Lat->E->delta[0]
        source = LevS->LPsi->LocalPsi[R->ActualLocalPsiNo];
        oldsource = LevS->LPsi->OldLocalPsi[R->ActualLocalPsiNo];
        //SetArrayToDouble0((double *)source,LevS->MaxG*2);
        do { 
          a -= fabs(Lat->E->delta[0]) == 0 ? 0.1 : fabs(Lat->E->delta[0]);
          if (a < -M_PI/2.) a = -M_PI/2.;// for this to work we need the pre-estimation which leads us into a nice regime (without gradient being the better _initial_ guess for a Psi)
          dcos = cos(a);
          dsin = sin(a);
          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(oldsource[g])*dcos + c_re(ConDir[g])*dsin;
            c_im(source[g]) = c_im(oldsource[g])*dcos + c_im(ConDir[g])*dsin;
          }
          CalculateEnergy(P);
          f_a = P->Lat.E->TotalEnergy[0]; // grab second value at left border
          //fprintf(stderr,"(%i) fnl(%lg) = %lg, Check ConDir[0] = %lg+i%lg, source[0] = %lg+i%lg, oldsource[0] = %lg+i%lg, TotDens[0] = %lg\n", P->Par.me, a, f_a, ConDir[0].re, ConDir[0].im, source[0].re, source[0].im, oldsource[0].re, oldsource[0].im, R->Lev0->Dens->DensityArray[TotalDensity][0]);
        } while (f_a < f_m);

        //SetArrayToDouble0((double *)source,LevS->MaxG*2);
        do {
          b += fabs(Lat->E->delta[0]) == 0 ? 0.1 : fabs(Lat->E->delta[0]);
          if (b > M_PI/2.) b = M_PI/2.;
          dcos = cos(b);
          dsin = sin(b);
          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(oldsource[g])*dcos + c_re(ConDir[g])*dsin;
            c_im(source[g]) = c_im(oldsource[g])*dcos + c_im(ConDir[g])*dsin;
          }
          CalculateEnergy(P);
          f_b = P->Lat.E->TotalEnergy[0]; // grab second value at left border
          //fprintf(stderr,"(%i) fnl(%lg) = %lg\n", P->Par.me, b, f_b);
        } while (f_b < f_m);
        
        memcpy(source, oldsource, ElementSize*LevS->MaxG*sizeof(double)); // restore old Psi from OldPsi
        //fprintf(stderr,"(%i) Psi %i at %p retrieved from 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);
        CalculateEnergy(P);

        if(P->Call.out[ValueOut]) fprintf(stderr,"(%i) Preparing brent with f(a) (%lg,%lg)\t f(b) (%lg,%lg)\t f(m) (%lg,%lg) ...\n", P->Par.me,a,f_a,b,f_b,m,f_m);
        iter=0; 
        gsl_min_fminimizer_set_with_values (s, &F, m, f_m, a, f_a, b, f_b);
        if(P->Call.out[ValueOut]) fprintf (stderr,"(%i) using %s method\n",P->Par.me, gsl_min_fminimizer_name (s));
        if(P->Call.out[ValueOut]) fprintf (stderr,"(%i) %5s [%9s, %9s] %9s %9s\n",P->Par.me, "iter", "lower", "upper", "min", "err(est)");
        if(P->Call.out[ValueOut]) fprintf (stderr,"(%i) %5d [%.7f, %.7f] %.7f %.7f\n",P->Par.me, iter, a, b, m, b - a);
        do {
          iter++;
          status = gsl_min_fminimizer_iterate (s);
       
          m = gsl_min_fminimizer_x_minimum (s);
          a = gsl_min_fminimizer_x_lower (s);
          b = gsl_min_fminimizer_x_upper (s);
       
          status = gsl_min_test_interval (a, b, 0.001, 0.0);
       
          if (status == GSL_SUCCESS)
            if(P->Call.out[ValueOut]) fprintf (stderr,"(%i) Converged:\n",P->Par.me);
       
          if(P->Call.out[ValueOut]) fprintf (stderr,"(%i) %5d [%.7f, %.7f] %.7f %.7f\n",P->Par.me,
            iter, a, b, m, b - a);
        } while (status == GSL_CONTINUE && iter < max_iter);
        CalculateNewWave(P,&m);
        TestGramSch(P,LevS,Psi,Occupied);
        UpdateActualPsiNo(P, Occupied); // step on due setting to MaxPsiStep further above
        UpdateEnergyArray(P);
        CalculateEnergy(P);
        //fprintf(stderr,"(%i) Final value for Psi %i: %lg\n", P->Par.me, R->ActualLocalPsiNo, P->Lat.E->TotalEnergy[0]);
        R->MinStopStep =  R->ActualMaxMinStopStep; // check stop condition every time
        if (*SuperStop != 1)
          *SuperStop = CheckCPULIM(P);
        *Stop = CalculateMinimumStop(P, *SuperStop);
        R->OldActualLocalPsiNo = R->ActualLocalPsiNo;
        if (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent == 1) { // new wave function means new gradient!
          DoBrent = Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent;
          Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent = 2;
          //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.;
          CalculateNewWave(P,NULL);
          Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent = DoBrent;
        }
      }

      if (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent != 1) { // otherwise the following checks eliminiate stop=1 from above
        if (*SuperStop != 1)
          *SuperStop = CheckCPULIM(P);
        *Stop = CalculateMinimumStop(P, *SuperStop);
      }
      /*EnergyOutput(P, Stop);*/
      P->Speed.Steps++;
      R->LevS->Step++;
      /*ControlNativeDensity(P);*/
      //fprintf(stderr,"(%i) Stop %i\n",P->Par.me, Stop);
    }
    if (*SuperStop == 1) OutputVisSrcFiles(P, Occupied); // is now done after localization (ComputeMLWF())
  }
  TestGramSch(P,R->LevS,Psi, Occupied);
}

/** Minimisation of the PsiTagType#UnOccupied orbitals in the field of the occupied ones.
 * Assuming RunStruct#ActualLocalPsiNo is currenlty still an occupied wave function, we stop onward to the first
 * unoccupied and reset RunStruct#OldActualLocalPsiNo. Then it is checked whether CallOptions#ReadSrcFiles is set
 * and thus coefficients for the level have to be read from file and afterwards initialized.
 * 
 * Then follows the main loop, until a stop condition is met:
 * -# CalculateNewWave()\n
 *    Over a conjugate gradient method the next (minimal) wave function is sought for.
 * -# UpdateActualPsiNo()\n
 *    Switch local Psi to next one.
 * -# UpdateEnergyArray()\n
 *    Shift archived energy values to make space for next one.
 * -# UpdateDensityCalculation(), SpeedMeasure()'d in DensityTime\n
 *    Calculate TotalLocalDensity of LocalPsis and gather results as TotalDensity.
 * -# UpdatePsiEnergyCalculation()\n
 *    Calculate kinetic and non-local energy contributons from the Psis.
 * -# CalculateGapEnergy()\n
 *    Calculate Gap energies (Hartreepotential, Pseudo) and the gradient.
 * -# EnergyAllReduce()\n
 *    Gather PsiEnergy results from all processes and sum up together with all other contributions to TotalEnergy.
 * -# CheckCPULIM()\n
 *    Check if external signal has been received (e.g. end of time slit on cluster), break operation at next possible moment.
 * -# CalculateMinimumStop()\n
 *    Evaluates stop condition if desired precision or steps or ... have been reached. Otherwise go to
 *    CalculateNewWave().
 * 
 * Afterwards, the coefficients are written to file by OutputVisSrcFiles() if desired. Orthonormality is tested, we step
 * back to the occupied wave functions and the densities are re-initialized. 
 * \param *P Problem at hand
 * \param *Stop flag to determine if epsilon stop conditions have met
 * \param *SuperStop flag to determinte whether external signal's required end of calculations
 */
static void MinimiseUnoccupied (struct Problem *P, int *Stop, int *SuperStop) {
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  int StartLocalPsiNo;

  *Stop = 0;
  R->PsiStep = R->MaxPsiStep; // in case it's zero from CalculateForce()
  UpdateActualPsiNo(P, UnOccupied); // step on to next unoccupied one
  R->OldActualLocalPsiNo = R->ActualLocalPsiNo; // reset, otherwise OldActualLocalPsiNo still points to occupied wave function
  UpdateGramSchOldActualPsiNo(P,Psi);
  if ((P->Call.ReadSrcFiles == DoReadAllSrcDensities) && ReadSrcPsiDensity(P,UnOccupied,1, R->LevSNo)) {
    SpeedMeasure(P, InitSimTime, StartTimeDo);  
    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, UnOccupied, 0, R->LevSNo);
    if (P->Call.ReadSrcFiles < DoReadAndMinimise) {
      ResetGramSchTagType(P, Psi, UnOccupied, IsOrthonormal); // loaded values are orthonormal
      SpeedMeasure(P, DensityTime, StartTimeDo);
      InitDensityCalculation(P);
      SpeedMeasure(P, DensityTime, StopTimeDo);
      InitPsiEnergyCalculation(P,UnOccupied);  // go through all orbitals calculating kinetic and non-local
      //CalculateDensityEnergy(P, 0);
      StartLocalPsiNo = R->ActualLocalPsiNo;
      do { // otherwise OnePsiElementAddData#Lambda is calculated only for current Psi not for all
        CalculateGapEnergy(P);
        UpdateActualPsiNo(P, Occupied);
      } while (R->ActualLocalPsiNo != StartLocalPsiNo);
      EnergyAllReduce(P);
    }
    SpeedMeasure(P, InitSimTime, StopTimeDo);  
  }
  if (P->Call.ReadSrcFiles != DoReadAllSrcDensities) {
    SpeedMeasure(P, InitSimTime, StartTimeDo);  
    ResetGramSchTagType(P, Psi, UnOccupied, NotOrthogonal);
    SpeedMeasure(P, GramSchTime, StartTimeDo);  
    GramSch(P, R->LevS, Psi, Orthonormalize);
    SpeedMeasure(P, GramSchTime, StopTimeDo);  
    SpeedMeasure(P, InitDensityTime, StartTimeDo);
    InitDensityCalculation(P);
    SpeedMeasure(P, InitDensityTime, StopTimeDo);
    InitPsiEnergyCalculation(P,UnOccupied);  // go through all orbitals calculating kinetic and non-local
    //CalculateDensityEnergy(P, 0);
    CalculateGapEnergy(P);
    EnergyAllReduce(P);
    SpeedMeasure(P, InitSimTime, StopTimeDo);  
    R->LevS->Step++;
    EnergyOutput(P,0);
    if(P->Call.out[LeaderOut]) fprintf(stderr,"(%i)Beginning minimisation of type %s ...\n", P->Par.me, R->MinimisationName[R->CurrentMin]);
    while (*Stop != 1) {
      CalculateNewWave(P,NULL);
      UpdateActualPsiNo(P, UnOccupied);
      /* New */
      UpdateEnergyArray(P);
      SpeedMeasure(P, DensityTime, StartTimeDo);
      UpdateDensityCalculation(P);
      SpeedMeasure(P, DensityTime, StopTimeDo);
      UpdatePsiEnergyCalculation(P);
      //CalculateDensityEnergy(P, 0);
      CalculateGapEnergy(P);  //calculates XC, HGDensity, afterwards gradient, where V_xc is added upon HGDensity
      EnergyAllReduce(P);
      if (*SuperStop != 1)
        *SuperStop = CheckCPULIM(P);
      *Stop = CalculateMinimumStop(P, *SuperStop);
      /*EnergyOutput(P, Stop);*/
      P->Speed.Steps++;
      R->LevS->Step++;
      /*ControlNativeDensity(P);*/
    }
    OutputVisSrcFiles(P, UnOccupied);
//    if (!TestReadnWriteSrcDensity(P,UnOccupied))
//      Error(SomeError,"TestReadnWriteSrcDensity failed!");
  }
  TestGramSch(P,R->LevS,Psi, UnOccupied);
  ResetGramSchTagType(P, Psi, UnOccupied, NotUsedToOrtho);
  // re-calculate Occupied values (in preparation for perturbed ones)
  UpdateActualPsiNo(P, Occupied); // step on to next occupied one
  SpeedMeasure(P, DensityTime, StartTimeDo);
  InitDensityCalculation(P);  // re-init densities to occupied standard
  SpeedMeasure(P, DensityTime, StopTimeDo);
//  InitPsiEnergyCalculation(P,Occupied);
//  CalculateDensityEnergy(P, 0);
//  EnergyAllReduce(P);
}


/** Calculate the forces.
 * From RunStruct::LevSNo downto RunStruct::InitLevSNo the following routines are called in a loop:
 * -# In case of RunStruct#DoSeparated another loop begins for the unoccupied states with some reinitalization
 *    before and afterwards. The loop however is much the same as the one above.
 * -# ChangeToLevUp()\n
 *    Repeat the loop or when the above stop is reached, the level is changed and the loop repeated.
 * 
 * Afterwards comes the actual force and energy calculation by calling:
 * -# ControlNativeDensity()\n
 *    Checks if the density still reproduces particle number.
 * -# CalculateIonLocalForce(), SpeedMeasure()'d in LocFTime\n
 *    Calculale local part of force acting on Ions.
 * -# CalculateIonNonLocalForce(), SpeedMeasure()'d in NonLocFTime\n
 *    Calculale local part of force acting on Ions.
 * -# CalculateEwald()\n
 *    Calculate Ewald force acting on Ions.
 * -# CalculateIonForce()\n
 *    Sum up those three contributions.
 * -# CorrectForces()\n
 *    Shifts center of gravity of all forces for each Ion, so that the cell itself remains at rest.
 * -# GetOuterStop()
 *    Calculates a mean force per Ion.
 * \param *P Problem at hand
 * \return 1 - cpulim received, break operation, 0 - continue as normal
 */
int CalculateForce(struct Problem *P)
{
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct LatticeLevel *LevS = R->LevS;
  struct FileData *F = &P->Files;
  struct Ions *I = &P->Ion;
  int Stop=0, SuperStop = 0, OuterStop = 0;
  //int i, j;
  SpeedMeasure(P, SimTime, StartTimeDo);
  if ((F->DoOutVis == 2) || (P->Call.ForcesFile == NULL)) {  // if we want to draw those pretty density pictures, we have to solve the ground state in any case
    while ((R->LevSNo > R->InitLevSNo) || (!Stop && R->LevSNo == R->InitLevSNo)) {
      // occupied
      R->PsiStep = R->MaxPsiStep; // reset in-Psi-minimisation-counter, so that we really advance to the next wave function
      R->OldActualLocalPsiNo = R->ActualLocalPsiNo; // reset OldActualLocalPsiNo, as it might still point to a perturbed wave function from last level
      UpdateGramSchOldActualPsiNo(P,Psi);
      ControlNativeDensity(P);
      MinimiseOccupied(P, &Stop, &SuperStop);
      if (!I->StructOpt) {
        if (((P->Call.ReadSrcFiles != DoReadAllSrcDensities) || (P->Call.ReadSrcFiles != DoReadOccupiedSrcDensities)) || (!ParseWannierFile(P))) { // only localize and store if they have just been minimised (hence don't come solely from file), otherwise read stored values from file (if possible)
          SpeedMeasure(P, WannierTime, StartTimeDo);
          ComputeMLWF(P);   // localization of orbitals
          SpeedMeasure(P, WannierTime, StopTimeDo);
          OutputVisSrcFiles(P, Occupied); // rewrite now localized orbitals
    //      if (!TestReadnWriteSrcDensity(P,Occupied))
    //        Error(SomeError,"TestReadnWriteSrcDensity failed!");
        }
      
//      // plot psi cuts
//      for (i=0; i < Psi->MaxPsiOfType; i++)  // go through all wave functions (here without the extra ones for each process)
//        if ((Psi->AllPsiStatus[i].PsiType == Occupied) && (Psi->AllPsiStatus[i].my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi))
//          for (j=0;j<NDIM;j++) {
//              //fprintf(stderr,"(%i) Plotting Psi %i/%i cut axis %i at coordinate %lg \n",P->Par.me, i, Psi->AllPsiStatus[i].MyGlobalNo, j, Lat->Psi.AddData[Psi->AllPsiStatus[i].MyLocalNo].WannierCentre[j]);
//              CalculateOneDensityR(Lat, R->LevS, R->Lev0->Dens, R->LevS->LPsi->LocalPsi[Psi->AllPsiStatus[i].MyLocalNo], R->Lev0->Dens->DensityArray[ActualDensity], R->FactorDensityR, 0);
//              PlotSrcPlane(P, j, Lat->Psi.AddData[Psi->AllPsiStatus[i].MyLocalNo].WannierCentre[j], Psi->AllPsiStatus[i].MyGlobalNo, R->Lev0->Dens->DensityArray[ActualDensity]);
//            }
          
        // unoccupied calc 
        if (R->DoUnOccupied) {
          MinimiseUnoccupied(P, &Stop, &SuperStop);
        }
        // hamiltonian
        CalculateHamiltonian(P);  // lambda_{kl} needed (and for bandgap after UnOccupied)

        //TestSawtooth(P, 0);
        //TestSawtooth(P, 1);
        //TestSawtooth(P, 2);
                
        // perturbed calc
        if ((R->DoPerturbation)) { // && R->LevSNo <= R->InitLevSNo) {
          AllocCurrentDensity(R->Lev0->Dens);// lock current density arrays
          MinimisePerturbed(P, &Stop, &SuperStop); // herein InitDensityCalculation() is called, thus no need to call it beforehand

          SpeedMeasure(P, CurrDensTime, StartTimeDo);
          if (SuperStop != 1) {
            if ((R->DoFullCurrent == 1) || ((R->DoFullCurrent == 2) && (CheckOrbitalOverlap(P) == 1))) { //test to check whether orbitals have mutual overlap and thus \\DeltaJ_{xc} must not be dropped
              R->DoFullCurrent = 1; // set to 1 if it was 2 but Check...() yielded necessity
              //debug(P,"Filling with Delta j ...");
              FillDeltaCurrentDensity(P);
            }
          }
          SpeedMeasure(P, CurrDensTime, StopTimeDo);
          TestCurrent(P,0);
          TestCurrent(P,1);
          TestCurrent(P,2);
          if (F->DoOutCurr && R->Lev0->LevelNo == 0)  // only output in uppermost level)
            OutputCurrentDensity(P);
          if (R->VectorPlane != -1)
            PlotVectorPlane(P,R->VectorPlane,R->VectorCut);
          CalculateMagneticSusceptibility(P);
          debug(P,"Normal calculation of shielding over R-space");
          CalculateMagneticMoment(P);
          CalculateChemicalShieldingByReciprocalCurrentDensity(P);
          SpeedMeasure(P, CurrDensTime, StopTimeDo);
          DisAllocCurrentDensity(R->Lev0->Dens);    // unlock current density arrays
        }  // end of if perturbation
        InitDensityCalculation(P);  // all unperturbed(!) wave functions've "changed" from ComputeMLWF(), thus reinit density
      } else  // end of if StructOpt or MaxOuterStep
                OutputVisSrcFiles(P, Occupied); // in structopt or MD write for every level

      if ((!I->StructOpt) && (!R->MaxOuterStep))  // print intermediate levels energy results if we don't do MD or StructOpt
            EnergyOutput(P, 1);
      // next level
      ChangeToLevUp(P, &Stop);
      //if (isnan(LevS->LPsi->LocalPsi[R->ActualLocalPsiNo][0].re)) { fprintf(stderr,"(%i) WARNING in ChangeToLevUp(): LPsi->LocalPsi[%i]_[%i] = NaN!\n", P->Par.me, R->ActualLocalPsiNo, 0); Error(SomeError, "NaN-Fehler!"); }      
      LevS = R->LevS; // re-set pointer that's changed from LevUp
    }
          SpeedMeasure(P, SimTime, StopTimeDo);
          ControlNativeDensity(P);
          TestGramSch(P,LevS,Psi, Occupied);
          // necessary for correct ionic forces ...
                SpeedMeasure(P, LocFTime, StartTimeDo);
                CalculateIonLocalForce(P);
                SpeedMeasure(P, LocFTime, StopTimeDo);
                SpeedMeasure(P, NonLocFTime, StartTimeDo);
                CalculateIonNonLocalForce(P);
                SpeedMeasure(P, NonLocFTime, StopTimeDo);
                CalculateEwald(P, 1);
                CalculateIonForce(P);
  }
  if (P->Call.ForcesFile != NULL) { // if we parse forces from file, values are written over (in case of DoOutVis)
    fprintf(stderr, "Parsing Forces from file.\n");
    ParseIonForce(P); 
    //CalculateIonForce(P);
  }
  CorrectForces(P);
  // ... on output of densities
  if (F->DoOutOrbitals) { // output of each orbital
    debug(P,"OutputVisAllOrbital");
    OutputVisAllOrbital(P,0,1,Occupied);
  }
  //OutputNorm(P);
  //fprintf(stderr,"(%i) DoubleG: %p, CArray[22]: %p, OldLocalPsi: %p\n", P->Par.me, R->LevS->DoubleG, R->Lev0->Dens->DensityCArray[22], R->LevS->LPsi->OldLocalPsi);
  //OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
  /*TestGramSch(P, R->LevS, &P->Lat.Psi); */
  GetOuterStop(P);
  //fprintf(stderr,"(%i) DoubleG: %p, CArray[22]: %p, OldLocalPsi: %p\n", P->Par.me, R->LevS->DoubleG, R->Lev0->Dens->DensityCArray[22], R->LevS->LPsi->OldLocalPsi); 
  if (SuperStop) OuterStop = 1;
  return OuterStop;
}

/** Checks whether the given positions \a *v have changed wrt stored in IonData structure.
 * \param *P Problem at hand
 * \param *v gsl_vector storing new positions
 */
int CheckForChangedPositions(struct Problem *P, const gsl_vector *v)
{
  struct Ions *I = &P->Ion;
  int is,ia,k, index=0;
  int diff = 0;
  double *R_ion;
  for (is=0; is < I->Max_Types; is++) // for all elements
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
      R_ion = &I->I[is].R[NDIM*ia];
      for (k=0;k<NDIM;k++) { // for all dimensions
        if (fabs(R_ion[k] - gsl_vector_get (v, index++)) > MYEPSILON)
          diff++;
      }      
    }
  return diff;  
}

/** Wrapper for CalculateForce() for simplex minimisation of total energy.
 * \param *v vector with degrees of freedom
 * \param *params additional arguments, here Problem at hand
 */
double StructOpt_func(const gsl_vector *v, void *params)
{
  struct Problem *P = (struct Problem *)params;
  struct RunStruct *R = &P->R;
  struct Ions *I = &P->Ion;
  struct Energy *E = P->Lat.E;
  int i;
  double *R_ion, *R_old, *R_old_old;//, *FIon;
  //double norm = 0.;
  int is,ia,k,index = 0;
  int OuterStop;
  double diff = 0., tmp;
  debug (P, "StructOpt_func");
  if (CheckForChangedPositions(P,v)) {
    // update ion positions from vector coordinates
    for (is=0; is < I->Max_Types; is++) // for all elements
      for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
        R_ion = &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];
        tmp = 0.;
        for (k=0;k<NDIM;k++) { // for all dimensions
          R_old_old[k] = R_old[k];
          R_old[k] = R_ion[k];
          tmp += (R_ion[k]-gsl_vector_get (v, index))*(R_ion[k]-gsl_vector_get (v, index));
          R_ion[k] = gsl_vector_get (v, index++);        
        }
        diff += sqrt(tmp);
      }  
    if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Summed Difference to former position %lg\n", P->Par.me, diff);
    // recalculate ionic forces (do electronic minimisation)
    R->OuterStep++;
    R->NewRStep++;
    UpdateWaveAfterIonMove(P);
    for (i=MAXOLD-1; i > 0; i--)    // store away old energies
      E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
    UpdateToNewWaves(P);
    E->TotalEnergyOuter[0] = E->TotalEnergy[0];
    OuterStop = CalculateForce(P);
    //UpdateIonsU(P);
    //CorrectVelocity(P);
    //CalculateEnergyIonsU(P);
  /*  if ((P->R.ScaleTempStep > 0) && ((R->OuterStep-1) % P->R.ScaleTempStep == 0))
      ScaleTemp(P);*/
    if ((R->OuterStep-1) % P->R.OutSrcStep == 0) 
      OutputVisSrcFiles(P, Occupied);
    if ((R->OuterStep-1) % P->R.OutVisStep == 0) {      
  /*      // recalculate density for the specific wave function ...
        CalculateOneDensityR(Lat, LevS, Dens0, PsiDat, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0);
        // ... and output (wherein ActualDensity is used instead of TotalDensity)
        OutputVis(P);
        OutputIonForce(P);
        EnergyOutput(P, 1);*/
    }
  }
  if (P->Par.me == 0) fprintf(stderr,"(%i) TE %e\n",P->Par.me, E->TotalEnergy[0]);
  return E->TotalEnergy[0];
}

/** Wrapper for CalculateForce() for simplex minimisation of ionic forces.
 * \param *v vector with degrees of freedom
 * \param *params additional arguments, here Problem at hand
 */
double StructOpt_f(const gsl_vector *v, void *params)
{
  struct Problem *P = (struct Problem *)params;
  struct RunStruct *R = &P->R;
  struct Ions *I = &P->Ion;
  struct Energy *E = P->Lat.E;
  int i;
  double *R_ion, *R_old, *R_old_old;//, *FIon;
  //double norm = 0.;
  int is,ia,k,index = 0;
  int OuterStop;
  double diff = 0., tmp;
  //debug (P, "StructOpt_f");
  if (CheckForChangedPositions(P,v)) {
    // update ion positions from vector coordinates
    for (is=0; is < I->Max_Types; is++) // for all elements
      for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
        R_ion = &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];
        tmp = 0.;
        for (k=0;k<NDIM;k++) { // for all dimensions
          R_old_old[k] = R_old[k];
          R_old[k] = R_ion[k];
          tmp += (R_ion[k]-gsl_vector_get (v, index))*(R_ion[k]-gsl_vector_get (v, index));
          R_ion[k] = gsl_vector_get (v, index++);        
        }
        diff += sqrt(tmp);
      }  
    if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Summed Difference to former position %lg\n", P->Par.me, diff);
    // recalculate ionic forces (do electronic minimisation)
    //R->OuterStep++;
    R->NewRStep++;
    UpdateWaveAfterIonMove(P);
    for (i=MAXOLD-1; i > 0; i--)    // store away old energies
      E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
    UpdateToNewWaves(P);
    E->TotalEnergyOuter[0] = E->TotalEnergy[0];
    OuterStop = CalculateForce(P);
    //UpdateIonsU(P);
    //CorrectVelocity(P);
    //CalculateEnergyIonsU(P);
  /*  if ((P->R.ScaleTempStep > 0) && ((R->OuterStep-1) % P->R.ScaleTempStep == 0))
      ScaleTemp(P);*/
    if ((R->OuterStep-1) % P->R.OutSrcStep == 0) 
      OutputVisSrcFiles(P, Occupied);
    /*if ((R->OuterStep-1) % P->R.OutVisStep == 0) {      
        // recalculate density for the specific wave function ...
        CalculateOneDensityR(Lat, LevS, Dens0, PsiDat, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0);
        // ... and output (wherein ActualDensity is used instead of TotalDensity)
        OutputVis(P);
        OutputIonForce(P);
        EnergyOutput(P, 1);
    }*/
  }
  GetOuterStop(P);
  //if (P->Call.out[LeaderOut] && (P->Par.me == 0)) fprintf(stderr,"(%i) Absolute Force summed over all Ions %e\n",P->Par.me, norm);
  return R->MeanForce[0];
  //if (P->Call.out[LeaderOut] && (P->Par.me == 0)) fprintf(stderr,"(%i) Struct_optf returning: %lg\n",P->Par.me,E->TotalEnergy[0]); 
  //return E->TotalEnergy[0];
}

void StructOpt_df(const gsl_vector *v, void *params, gsl_vector *df) 
{
  struct Problem *P = (struct Problem *)params;
  struct Ions *I = &P->Ion;
  double *FIon;
  int is,ia,k, index=0;
  //debug (P, "StructOpt_df");
  // look through coordinate vector if positions have changed sind last StructOpt_f call
  if (CheckForChangedPositions(P,v)) {// if so, recalc to update forces
    debug (P, "Calling StructOpt_f to update");
    StructOpt_f(v, params);
  }
  for (is=0; is < I->Max_Types; is++) // for all elements
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
      FIon = &I->I[is].FIon[NDIM*ia];
      for (k=0;k<NDIM;k++) { // for all dimensions
        gsl_vector_set (df, index++, FIon[k]);
      }      
    } 
  if (P->Call.out[LeaderOut] && (P->Par.me == 0)) {
    fprintf(stderr,"(%i) Struct_Optdf returning",P->Par.me); 
    gsl_vector_fprintf(stderr, df, "%lg");
  } 
}

void StructOpt_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df)
{
  *f = StructOpt_f(x, params);
  StructOpt_df(x, params, df);
}


/** CG implementation for the structure optimization.
 * We follow the example from the GSL manual.
 * \param *P Problem at hand
 */
void UpdateIon_PRCG(struct Problem *P)
{
  //struct RunStruct *Run = &P->R;
  struct Ions *I = &P->Ion;
  size_t np = NDIM*I->Max_TotalIons;  // d.o.f = number of ions times number of dimensions
  int is, ia, k, index;
  double *R;
 
  const gsl_multimin_fdfminimizer_type *T;
  gsl_multimin_fdfminimizer *s;
  gsl_vector *x;
  gsl_multimin_function_fdf minex_func;

  size_t iter = 0;
  int status;

  /* Starting point */
  x = gsl_vector_alloc (np);
  //fprintf(stderr,"(%i) d.o.f. = %i\n", P->Par.me, (int)np);

  index=0;
  for (is=0; is < I->Max_Types; is++) // for all elements
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
      R = &I->I[is].R[NDIM*ia];
      for (k=0;k<NDIM;k++) // for all dimensions
        gsl_vector_set (x, index++, R[k]);      
    }

  /* Initialize method and iterate */
  minex_func.f = &StructOpt_f;
  minex_func.df = &StructOpt_df;
  minex_func.fdf = &StructOpt_fdf;
  minex_func.n = np;
  minex_func.params = (void *)P;
  
  T = gsl_multimin_fdfminimizer_conjugate_pr;
  s = gsl_multimin_fdfminimizer_alloc (T, np);

  gsl_multimin_fdfminimizer_set (s, &minex_func, x, 0.1, 0.001);

  fprintf(stderr,"(%i) Commencing Structure optimization with PRCG: dof %d\n", P->Par.me,(int)np);
  do {
    iter++;
    status = gsl_multimin_fdfminimizer_iterate(s);

    if (status)
      break;
    
    status = gsl_multimin_test_gradient (s->gradient, 1e-2);
       
    if (status == GSL_SUCCESS)
      if (P->Par.me == 0) fprintf (stderr,"(%i) converged to minimum at\n", P->Par.me);

    if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) Commencing '%s' step %i ... \n",P->Par.me, gsl_multimin_fdfminimizer_name(s), P->R.StructOptStep);
    if ((P->Call.out[NormalOut]) && (P->Par.me == 0)) fprintf (stderr, "(%i) %5d %10.5f\n", P->Par.me, (int)iter, s->f);
    //gsl_vector_fprintf(stderr, s->dx, "%lg");
    OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
    OutputIonCoordinates(P, 0);
    P->R.StructOptStep++;
  } while ((status == GSL_CONTINUE) && (P->R.StructOptStep < P->R.MaxStructOptStep));

  gsl_vector_free(x);
  gsl_multimin_fdfminimizer_free (s);
}

/** Simplex implementation for the structure optimization.
 * We follow the example from the GSL manual.
 * \param *P Problem at hand
 */
void UpdateIon_Simplex(struct Problem *P)
{
  struct RunStruct *Run = &P->R;
  struct Ions *I = &P->Ion;
  size_t np = NDIM*I->Max_TotalIons;  // d.o.f = number of ions times number of dimensions
  int is, ia, k, index;
  double *R;
 
  const gsl_multimin_fminimizer_type *T;
  gsl_multimin_fminimizer *s;
  gsl_vector *x, *ss;
  gsl_multimin_function minex_func;

  size_t iter = 0;
  int status;
  double size;

  ss = gsl_vector_alloc (np);
  gsl_vector_set_all(ss, .2);
  /* Starting point */
  x = gsl_vector_alloc (np);
  //fprintf(stderr,"(%i) d.o.f. = %i\n", P->Par.me, (int)np);

  index=0;
  for (is=0; is < I->Max_Types; is++) // for all elements
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
      R = &I->I[is].R[NDIM*ia];
      for (k=0;k<NDIM;k++) // for all dimensions
        gsl_vector_set (x, index++, R[k]);      
    }

  /* Initialize method and iterate */
  minex_func.f = &StructOpt_f;
  minex_func.n = np;
  minex_func.params = (void *)P;
  
  T = gsl_multimin_fminimizer_nmsimplex;
  s = gsl_multimin_fminimizer_alloc (T, np);

  gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

  fprintf(stderr,"(%i) Commencing Structure optimization with NM simplex: dof %d\n", P->Par.me, (int)np);
  do {
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);

    if (status)
      break;
    
    size = gsl_multimin_fminimizer_size (s);
    status = gsl_multimin_test_size (size, 1e-4);
       
    if (status == GSL_SUCCESS)
      if (P->Par.me == 0) fprintf (stderr,"(%i) converged to minimum at\n", P->Par.me);

    if (P->Call.out[MinOut]) fprintf(stderr,"(%i) Commencing '%s' step %i ... \n",P->Par.me, gsl_multimin_fminimizer_name(s), P->R.StructOptStep);
    if ((P->Call.out[MinOut]) && (P->Par.me == 0)) fprintf (stderr, "(%i) %5d %10.5f %10.5f\n", P->Par.me, (int)iter, s->fval, size);
    OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
    OutputIonCoordinates(P, 0);
    P->R.StructOptStep++;
  } while ((status == GSL_CONTINUE) && (Run->OuterStep < Run->MaxOuterStep));

  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);
}

/** Implementation of various thermostats.
 * All these thermostats apply an additional force which has the following forms:
 * -# Woodcock
 *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
 * -# Gaussian
 *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
 * -# Langevin
 *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$  
 * -# Berendsen
 *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
 * -# Nose-Hoover
 *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
 * These Thermostats either simply rescale the velocities, thus Thermostats() should be called after UpdateIonsU(), and/or
 * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
 * belatedly and the constraint force set.
 * \param *P Problem at hand
 * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
 * \sa InitThermostat()
 */
void Thermostats(struct Problem *P, enum thermostats i)
{
  struct FileData *Files = &P->Files;
  struct Ions *I = &P->Ion;
  int is, ia, d;
  double *U;
  double a, ekin = 0.;
  double E = 0., F = 0.;
  double delta_alpha = 0.;
  const int delta_t = P->R.delta_t;
  double ScaleTempFactor;
  double sigma;
  gsl_rng * r;
  const gsl_rng_type * T;
  
  // calculate current temperature
  CalculateEnergyIonsU(P); // Temperature now in I->ActualTemp
  ScaleTempFactor = P->R.TargetTemp/I->ActualTemp;
  //if ((P->Par.me == 0) && (I->ActualTemp < MYEPSILON)) fprintf(stderr,"Thermostat: (1) I->ActualTemp = %lg",I->ActualTemp);
  if (Files->MeOutMes) fprintf(Files->TemperatureFile, "%d\t%lg",P->R.OuterStep, I->ActualTemp);
    
  // differentating between the various thermostats
  switch(i) {
     case None:
      debug(P, "Applying no thermostat...");
      break;
     case Woodcock:
      if ((P->R.ScaleTempStep > 0) && ((P->R.OuterStep-1) % P->R.ScaleTempStep == 0)) {
        debug(P, "Applying Woodcock thermostat...");
        for (is=0; is < I->Max_Types; is++) {
          a = 0.5*I->I[is].IonMass;
          for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
            U = &I->I[is].U[NDIM*ia];
            if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
              for (d=0; d<NDIM; d++) {
                U[d] *= sqrt(ScaleTempFactor);
                ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
              }
          }
        }
      }
      break;
     case Gaussian:
      debug(P, "Applying Gaussian thermostat...");
      for (is=0; is < I->Max_Types; is++) { // sum up constraint constant
        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
          U = &I->I[is].U[NDIM*ia];
          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
            for (d=0; d<NDIM; d++) {
              F += U[d] * I->I[is].FIon[d+NDIM*ia]; 
              E += U[d]*U[d]*I->I[is].IonMass;
            }
        }
      }
      if (P->Call.out[ValueOut]) fprintf(stderr, "(%i) Gaussian Least Constraint constant is %lg\n", P->Par.me, F/E);
      for (is=0; is < I->Max_Types; is++) { // apply constraint constant on constraint force and on velocities
        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
          U = &I->I[is].U[NDIM*ia];
          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
            for (d=0; d<NDIM; d++) {
              I->I[is].FConstraint[d+NDIM*ia] = (F/E) * (U[d]*I->I[is].IonMass);
              U[d] += delta_t/I->I[is].IonMass * (I->I[is].FConstraint[d+NDIM*ia]);
              ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
            }
        }
      }
      break;
     case Langevin:
      debug(P, "Applying Langevin thermostat...");
      // init random number generator
      gsl_rng_env_setup();
      T = gsl_rng_default;
      r = gsl_rng_alloc (T);
      // Go through each ion
      for (is=0; is < I->Max_Types; is++) {
        sigma  = sqrt(P->R.TargetTemp/I->I[is].IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
          U = &I->I[is].U[NDIM*ia];
          // throw a dice to determine whether it gets hit by a heat bath particle
          if (((((rand()/(double)RAND_MAX))*P->R.TempFrequency) < 1.)) {  // (I->I[is].IMT[ia] == MoveIon) &&  even FixedIon moves, only not by other's forces
            if (P->Par.me == 0) fprintf(stderr,"(%i) Particle %i,%i was hit (sigma %lg): %lg -> ", P->Par.me, is, ia, sigma, sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]));
            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
            for (d=0; d<NDIM; d++) {
              U[d] = gsl_ran_gaussian (r, sigma);
            }
            if (P->Par.me == 0) fprintf(stderr,"%lg\n", sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]));
          }
          for (d=0; d<NDIM; d++)
            ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
        }
      }
      break;
     case Berendsen:
      debug(P, "Applying Berendsen-VanGunsteren thermostat...");
      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];
          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
            for (d=0; d<NDIM; d++) {
              U[d] *= sqrt(1+(P->R.delta_t/P->R.TempFrequency)*(ScaleTempFactor-1));
              ekin += 0.5*I->I[is].IonMass * U[d]*U[d];
            }
        }
      }
      break;
     case NoseHoover:
      debug(P, "Applying Nose-Hoover thermostat...");
      // dynamically evolve alpha (the additional degree of freedom)
      delta_alpha = 0.;
      for (is=0; is < I->Max_Types; is++) { // sum up constraint constant
        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
          U = &I->I[is].U[NDIM*ia];
          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
            for (d=0; d<NDIM; d++) {
              delta_alpha += U[d]*U[d]*I->I[is].IonMass;
            }
        }
      }
      delta_alpha = (delta_alpha - (3.*I->Max_TotalIons+1.) * P->R.TargetTemp)/(P->R.HooverMass*Units2Electronmass);
      P->R.alpha += delta_alpha*delta_t;
      if (P->Par.me == 0) fprintf(stderr,"(%i) alpha = %lg * %i = %lg\n", P->Par.me, delta_alpha, delta_t, P->R.alpha);
      // apply updated alpha as additional force
      for (is=0; is < I->Max_Types; is++) { // apply constraint constant on constraint force and on velocities
        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
          U = &I->I[is].U[NDIM*ia];
          if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
            for (d=0; d<NDIM; d++) {
              I->I[is].FConstraint[d+NDIM*ia] = - P->R.alpha * (U[d] * I->I[is].IonMass);
              U[d] += delta_t/I->I[is].IonMass * (I->I[is].FConstraint[d+NDIM*ia]);
              ekin += (0.5*I->I[is].IonMass) * U[d]*U[d];
            }
        }
      }
      break;
  }    
  I->EKin = ekin;
  I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);    
  //if ((P->Par.me == 0) && (I->ActualTemp < MYEPSILON)) fprintf(stderr,"Thermostat: (2) I->ActualTemp = %lg",I->ActualTemp);
  if (Files->MeOutMes) { fprintf(Files->TemperatureFile, "\t%lg\n", I->ActualTemp); fflush(Files->TemperatureFile); }
}

/** Does the Molecular Dynamics Calculations.
 * All of the following is SpeedMeasure()'d in SimTime.
 * Initialization by calling:
 * -# CorrectVelocity()\n
 *    Shifts center of gravity of Ions momenta, so that the cell itself remains at rest.
 * -# CalculateEnergyIonsU(), SpeedMeasure()'d in TimeTypes#InitSimTime\n
 *    Calculates kinetic energy of "movable" Ions.
 * -# CalculateForce()\n
 *    Does the minimisation, calculates densities, then energies and finally the forces.
 * -# OutputVisSrcFiles()\n
 *    If desired, so-far made calculations are stored to file for later restarting.
 * -# OutputIonForce()\n
 *    Write ion forces to file.
 * -# EnergyOutput()\n
 *    Write calculated energies to screen or file.
 * 
 * The simulation phase begins: 
 * -# UpdateIonsR()\n
 *    Move Ions according to the calculated force.
 * -# UpdateWaveAfterIonMove()\n
 *    Update wave functions by averaging LocalPsi coefficients after the Ions have been shifted.
 * -# UpdateToNewWaves()\n
 *    Update after wave functions have changed.
 * -# CalculateForce()\n
 *    Does the minimisation, calculates densities, then energies and finally the forces.
 * -# UpdateIonsU()\n
 *    Change ion's velocities according to the calculated acting force.
 * -# CorrectVelocity()\n
 *    Shifts center of gravity of Ions momenta, so that the cell itself remains at rest.
 * -# CalculateEnergyIonsU()\n
 *    Calculates kinetic energy of "movable" Ions.
 * -# ScaleTemp()\n
 *    The temperature is scaled, so the systems energy remains constant (they must not gain momenta out of nothing)
 * -# OutputVisSrcFiles()\n
 *    If desired, so-far made calculations are stored to file for later restarting.
 * -# OutputVis()\n
 *    Visulization data for OpenDX is written at certain steps if desired.
 * -# OutputIonForce()\n
 *    Write ion forces to file.
 * -# EnergyOutput()\n
 *    Write calculated energies to screen or file.
 * 
 * After the ground state is found:
 * -# CalculateUnOccupied()\n
 *    Energies of unoccupied orbitals - that have been left out completely so far -
 *    are calculated.
 * -# TestGramSch()\n
 *    Test if orbitals are still orthogonal.
 * -# CalculateHamiltonian()\n
 *    Construct Hamiltonian and calculate Eigenvalues.
 * -# ComputeMLWF()\n
 *    Localize orbital wave functions.
 * 
 * \param *P Problem at hand
 */
void CalculateMD(struct Problem *P)
{
  struct RunStruct *R = &P->R;
  struct Ions *I = &P->Ion;
  struct Energy *E = P->Lat.E;
  int OuterStop = 0;
  int i;
  
  SpeedMeasure(P, SimTime, StartTimeDo);
  // initial calculations (bring density on BO surface and output start energies, coordinates, densities, ...)
  SpeedMeasure(P, InitSimTime, StartTimeDo);
  R->OuterStep = 0;
  CorrectVelocity(P);
  CalculateEnergyIonsU(P);
  OuterStop = CalculateForce(P);
  //R->OuterStep++;
  P->Speed.InitSteps++;
  SpeedMeasure(P, InitSimTime, StopTimeDo);

  OutputIonCoordinates(P, 1);
  OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
  OutputIonForce(P);
  EnergyOutput(P, 1);
  
  // if desired perform beforehand a structure relaxation/optimization
  if (I->StructOpt) {
    debug(P,"Commencing minimisation on ionic structure ...");
    R->StructOptStep = 0;
    //UpdateIon_PRCG(P);
    //UpdateIon_Simplex(P);
    while ((R->MeanForce[0] > 1e-4) && (R->StructOptStep < R->MaxStructOptStep)) {
      R->StructOptStep++; 
      OutputIonCoordinates(P, 1);
      UpdateIons(P); 
      UpdateWaveAfterIonMove(P);
      for (i=MAXOLD-1; i > 0; i--)    // store away old energies
        E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
      UpdateToNewWaves(P);
      E->TotalEnergyOuter[0] = E->TotalEnergy[0]; 
      OuterStop = CalculateForce(P);
      CalculateEnergyIonsU(P);
      if ((R->StructOptStep-1) % P->R.OutSrcStep == 0) 
        OutputVisSrcFiles(P, Occupied);
      if ((R->StructOptStep-1) % P->R.OutVisStep == 0) {      
        OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
        OutputIonForce(P);
        EnergyOutput(P, 1);
      }
      if (P->Par.me == 0) fprintf(stderr,"(%i) Mean force is %lg\n", P->Par.me, R->MeanForce[0]); 
    }
    OutputIonCoordinates(P, 1);
  }
  if (I->StructOpt && !OuterStop) {
    I->StructOpt = 0;
    OuterStop = CalculateForce(P);
  }
  
  // and now begin with the molecular dynamics simulation
  debug(P,"Commencing MD simulation ...");
  while (!OuterStop && R->OuterStep < R->MaxOuterStep) {
    R->OuterStep++;
    if (P->Par.me == 0) {
      if (R->OuterStep > 1) fprintf(stderr,"\b\b\b\b\b\b\b\b\b\b\b\b");
      fprintf(stderr,"Time: %f fs\r", R->t*Atomictime2Femtoseconds);
      fflush(stderr);
    }
    OuterStop = CalculateForce(P);
    P->R.t += P->R.delta_t;   // increase current time by delta_t
    R->NewRStep++;

    UpdateIonsU(P);
    CorrectVelocity(P);
    Thermostats(P, I->Thermostat);
    CalculateEnergyIonsU(P);

    UpdateIonsR(P);
    OutputIonCoordinates(P, 1);

    UpdateWaveAfterIonMove(P);
    for (i=MAXOLD-1; i > 0; i--)    // store away old energies
      E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
    UpdateToNewWaves(P);
    E->TotalEnergyOuter[0] = E->TotalEnergy[0]; 
    //if ((P->R.ScaleTempStep > 0) && ((R->OuterStep-1) % P->R.ScaleTempStep == 0))
    //  ScaleTemp(P);
    if ((R->OuterStep-1) % P->R.OutSrcStep == 0) 
      OutputVisSrcFiles(P, Occupied);
    if ((R->OuterStep-1) % P->R.OutVisStep == 0) {      
      OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
      OutputIonForce(P);
      EnergyOutput(P, 1);
    }
    ResetForces(P);
  }
  SpeedMeasure(P, SimTime, StopTimeDo);
  CloseOutputFiles(P);
}