/** \file ions.c
 * Ionic Force, speed, coordinate and energy calculation.
 * Herein are all the routines for the molecular dynamics calculations such as
 * reading of configuration files IonsInitRead() and
 * free'ing RemoveIonsRead(),
 * summing up forces CalculateIonForce(),
 * correcting for temperature CorrectVelocities(),
 * or for fixation center of balance CorrectForces(),
 * outputting them to file OutputIonForce(),
 * calculating kinetic CalculateEnergyIonsU(), ewald CalculateEwald() energies,
 * moving ion UpdateIonsR() (position) UpdateIonsU() (speed),
 * scaling to match some temperature ScaleTemp()
 * are gathered.
 * Finally, needed for structure optimization GetOuterStop() evaluates the stop
 * condition of a minimal mean force acting on each ion.
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers
 \date 2000

  File: ions.c
  $Id: ions.c,v 1.34 2007/02/05 16:15:27 foo Exp $
*/

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include "mymath.h"


// use double precision fft when we have it
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef HAVE_DFFTW_H
#include "dfftw.h"
#else
#include "fftw.h"
#endif

#include "data.h"
#include "errors.h"
#include "helpers.h"
#include "mymath.h"
#include "ions.h"
#include "init.h"

/** Readin of the ion section of parameter file.
 * Among others the following paramaters are read from file:
 * Ions::MaxTypes, IonType::Max_IonsOfType, IonType::Z, IonType::R,
 * IonType:IMT, ...
 * \param P Problem at hand (containing Ions and Lattice structures)
 * \param *source file pointer
 */
void IonsInitRead(struct Problem *P, FILE *source) {
	struct Ions *I = &P->Ion;
  //struct RunStruct *Run = &P->R;
  int i,it,j,IMT,BorAng, d,me, count;
  int k;
  struct Lattice *L = &P->Lat;
  int relative; // whether read-in coordinate are relative (1) or absolute  
  double Bohr = ANGSTROEMINBOHRRADIUS; /* Angstroem */
  double sigma, R[3], U[3];

  I->Max_TotalIons = 0;
  I->TotalZval = 0;
  MPI_Comm_rank(MPI_COMM_WORLD, &me);
  ParseForParameter(P->Call.out[ReadOut],source,"RCut", 0, 1, 1, double_type, &I->R_cut, critical);
  ParseForParameter(P->Call.out[ReadOut],source,"IsAngstroem", 0, 1, 1, int_type, &BorAng, critical);
  if (!I->IsAngstroem) {
  	Bohr = 1.0;
	} else {  // adapt RealBasis (is in Angstroem as well)
    SM(P->Lat.RealBasis, Bohr, NDIM_NDIM); 
    RMatReci3(P->Lat.InvBasis,  P->Lat.RealBasis);
  }
  ParseForParameter(P->Call.out[ReadOut],source,"MaxTypes", 0, 1, 1, int_type, &I->Max_Types, critical);
  I->I = (struct IonType *) Malloc(sizeof(struct IonType)*I->Max_Types, "IonsInitRead: IonType");
  if (!ParseForParameter(P->Call.out[ReadOut],source,"RelativeCoord", 0, 1, 1, int_type, &relative, optional))
    relative = 0;
  if (!ParseForParameter(P->Call.out[ReadOut],source,"StructOpt", 0, 1, 1, int_type, &I->StructOpt, optional))
    I->StructOpt = 0;

  /* Ions Data */
  I->RLatticeVec = NULL;
  I->TotalMass = 0;
  char *free_name, *name;
  name = free_name = Malloc(255*sizeof(char),"IonsInitRead: Name");
  for (i=0; i < I->Max_Types; i++) {
  	sprintf(name,"Ion_Type%i",i+1);
    I->I[i].corecorr = NotCoreCorrected;
    I->I[i].Name = MallocString(255, "IonsInitRead: Name");
    I->I[i].Symbol = MallocString(255, "IonsInitRead: Symbol");
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, int_type, &I->I[i].Max_IonsOfType, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, int_type, &I->I[i].Z, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &I->I[i].rgauss, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &I->I[i].l_max, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 5, 1, int_type, &I->I[i].l_loc, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 6, 1, double_type, &I->I[i].IonMass, critical);
    if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 7, 1, string_type, &I->I[i].Name[0], optional))
      sprintf(&I->I[i].Name[0],"type%i",i);
    if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 8, 1, string_type, &I->I[i].Symbol[0], optional))
      sprintf(&I->I[i].Symbol[0],"t%i",i);
    if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) Element name: %s, symbol: %s\n",P->Par.me, I->I[i].Name, I->I[i].Symbol);
    I->I[i].sigma = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma");
    I->I[i].sigma_rezi = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi");
    I->I[i].sigma_PAS = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_PAS");
    I->I[i].sigma_rezi_PAS = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi_PAS");
    for (it=0;it<I->I[i].Max_IonsOfType;it++) {
      I->I[i].sigma[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma");
      I->I[i].sigma_rezi[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma_rezi");
      I->I[i].sigma_PAS[it] = (double *) Malloc(sizeof(double) * NDIM, "IonsInitRead: sigma_PAS");
      I->I[i].sigma_rezi_PAS[it] = (double *) Malloc(sizeof(double) * NDIM, "IonsInitRead: sigma_rezi_PAS");
    }
    //I->I[i].IonFac = I->I[i].IonMass;
    I->I[i].IonFac = 1/I->I[i].IonMass;
    I->TotalMass += I->I[i].Max_IonsOfType*I->I[i].IonMass;
    I->I[i].ZFactor = -I->I[i].Z*4.*PI;
    I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
    I->I[i].R = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: R");
    I->I[i].R_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old");
    I->I[i].R_old_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old_old");
    I->I[i].FIon = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon");
    I->I[i].FIon_old = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon_old");
    SetArrayToDouble0(I->I[i].FIon_old, NDIM*I->I[i].Max_IonsOfType);
    I->I[i].SearchDir = Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: SearchDir");
    I->I[i].GammaA = Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: GammaA");
    SetArrayToDouble0(I->I[i].GammaA, I->I[i].Max_IonsOfType);
    I->I[i].U = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: U");
    SetArrayToDouble0(I->I[i].U, NDIM*I->I[i].Max_IonsOfType);
    I->I[i].SFactor = NULL;
    I->I[i].IMT = (enum IonMoveType *) Malloc(sizeof(enum IonMoveType) * I->I[i].Max_IonsOfType, "IonsInitRead: IMT");
    I->I[i].FIonL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonL");
    I->I[i].FIonNL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonNL");
    I->I[i].FEwald = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FEwald");
    I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
    for (j=0; j < I->I[i].Max_IonsOfType; j++) {
      I->I[i].alpha[j] = 2.;
    	sprintf(name,"Ion_Type%i_%i",i+1,j+1);
			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, double_type, &R[0], critical);
			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, double_type, &R[1], critical);
			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &R[2], critical);
			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &IMT, critical);
      // change if coordinates were relative
      if (relative) {
        //fprintf(stderr,"(%i)Ion coordinates are relative %i ... \n",P->Par.me, relative);
        RMat33Vec3(&I->I[i].R[0+NDIM*j], L->RealBasis, R); // multiply with real basis
      } else {
        for (k=0;k<NDIM;k++) // and copy to old vector
          I->I[i].R[k+NDIM*j] = R[k];
      }
      if (IMT < 0 || IMT >= MaxIonMoveType) {
				fprintf(stderr,"Bad Ion MoveType set to MoveIon for Ion (%i,%i)\n",i,j);
				IMT = 0;
      }
        
      if ((P->Call.out[ReadOut])) {  // rotate coordinates
        fprintf(stderr,"(%i) coordinates of Ion %i: (x,y,z) = (",P->Par.me,j);
        for(k=0;k<NDIM;k++) {
          fprintf(stderr,"%lg ",I->I[i].R[k+NDIM*j]);
          if (k != NDIM-1) fprintf(stderr,", ");
          else fprintf(stderr,")\n");
        }
      }

      I->I[i].IMT[j] = (enum IonMoveType)IMT;
      SM(&I->I[i].R[NDIM*j], 1./Bohr, NDIM);
      RMat33Vec3(R, L->InvBasis, &I->I[i].R[NDIM*j]);
      for (d=0; d <NDIM; d++) {
				while (R[d] < 0) 
				  R[d] += 1.0;
				while (R[d] >= 1.0) 
				  R[d] -= 1.0;
      }
      RMat33Vec3(&I->I[i].R[NDIM*j], L->RealBasis, R);
      for (d=0; d <NDIM; d++) {
				I->I[i].R_old_old[d+NDIM*j] = I->I[i].R_old[d+NDIM*j] = I->I[i].R[d+NDIM*j];
      }
      I->Max_TotalIons++;
    }
  }
  Free(free_name);
  I->Max_Max_IonsOfType = 0;
  for (i=0; i < I->Max_Types; i++) 
    I->Max_Max_IonsOfType = MAX(I->Max_Max_IonsOfType, I->I[i].Max_IonsOfType);
  I->FTemp = (double *) Malloc(sizeof(double) * NDIM * I->Max_Max_IonsOfType, "IonsInitRead:");
}

/** Calculated Ewald force.
 * \f[
 *  E_{K-K} = \frac{1}{2} \sum_L \sum_{K=\{(I_s,I_a,J_s,J_a)|R_{I_s,I_a} - R_{J_s,J_a} - L\neq 0\}} \frac{Z_{I_s} Z_{J_s}} {|R_{I_s,I_a} - R_{J_s,J_a} - L|} 
 *    \textnormal{erfc} \Bigl ( \frac{|R_{I_s,I_a} - R_{J_s,J_a} - L|} {\sqrt{r_{I_s}^{Gauss}+r_{J_s}^{Gauss}}}\Bigr )
 *    \qquad (4.10)
 * \f]
 * Calculates the core-to-core force between the ions of all super cells.
 * In order for this series to converge there must be a certain summation
 * applied, which is the ewald summation and which is nothing else but to move
 * in a circular motion from the current cell to the outside up to Ions::R_cut.
 * \param *P Problem at hand
 * \param first additional calculation beforehand if != 0
 */
void CalculateEwald(struct Problem *P, int first) 
{
  int r,is,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec;
  int MaxPar=P->Par.procs, ParMe=P->Par.me, MaxLocalVec, StartVec;
  double rcut2,r2,erre2,addesr,addpre,arg,fac,gkl,esr,fxx,repand;
  double R1[3];
  double R2[3];
  double R3[3];
  double t[3];
  struct Lattice *L = &P->Lat;
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  if (I->R_cut != 0.0) {
    if (first) {
      SpeedMeasure(P, InitSimTime, StopTimeDo);
      rcut2 = I->R_cut*I->R_cut;
      ActualVec =0;
      MaxCell = (int)5.*I->R_cut/pow(L->Volume,1./3.);
      if (MaxCell < 2) MaxCell = 2;
      for (i = -MaxCell; i <= MaxCell; i++) {
        for (j = -MaxCell; j <= MaxCell; j++) {
  				for (k = -MaxCell; k <= MaxCell; k++) {
  				  r2 = 0.0;
  				  for (ir=0; ir <3; ir++) {
  				    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
  				    r2 = t[ir]*t[ir];
  				  }
  				  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
  				  if ((r2 <= rcut2) || Is_Neighb_Cell) {
  				    ActualVec++;
  				  }
  				}
        }
      }
      MaxVec = ActualVec;
      I->MaxVec = ActualVec;
      MaxLocalVec = MaxVec / MaxPar;
      StartVec = ParMe*MaxLocalVec;
      r = MaxVec % MaxPar;
      if (ParMe >= r) {
        StartVec += r;
      } else {
        StartVec += ParMe;
        MaxLocalVec++;
      }
      I->MaxLocalVec = MaxLocalVec;
      LocalActualVec = ActualVec = 0;
      I->RLatticeVec = (double *)
        Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald:");
      for (i = -MaxCell; i <= MaxCell; i++) {
        for (j = -MaxCell; j <= MaxCell; j++) {
  				for (k = -MaxCell; k <= MaxCell; k++) {
  				  r2 = 0.0;
  				  for (ir=0; ir <3; ir++) {
  				    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
  				    r2 = t[ir]*t[ir];
  				  }
  				  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
  				  if ((r2 <= rcut2) || Is_Neighb_Cell) {
  				    if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) {
  				      I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
  				      I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
  				      I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
  				      LocalActualVec++;
  				    }
  				    ActualVec++;
  				  }
  				}
        }
      }
      SpeedMeasure(P, InitSimTime, StartTimeDo);
    }
    SpeedMeasure(P, EwaldTime, StartTimeDo);
    esr = 0.0;
    for (is1=0;is1 < I->Max_Types; is1++) 
      for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
        for (i=0; i< NDIM; i++)
  				I->FTemp[i+NDIM*(ia1)] = 0;
    for (is1=0;is1 < I->Max_Types; is1++) {
      for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
        for (i=0;i<3;i++) {
  				R1[i]=I->I[is1].R[i+NDIM*ia1];
        }
        for (ir2=0;ir2<I->MaxLocalVec;ir2++) {
  				for (is2=0;is2 < I->Max_Types; is2++) {
  				  gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
  				  for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
  				    for (i=0;i<3;i++) {
  				      R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2];
  				    }
  				    erre2=0.0;
  				    for (i=0;i<3;i++) {
  				      R3[i] = R1[i]-R2[i];
  				      erre2+=R3[i]*R3[i];
  				    }
  				    if (erre2 > MYEPSILON) {
  				      arg=sqrt(erre2);
  				      fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;
  				      
  				      arg *= gkl;
  				      addesr = derf(arg);
  				      addesr = (1.0-addesr)*fac;
  				      esr += addesr;
  				      addpre=exp(-arg*arg)*gkl;
  				      addpre=PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;
  				      repand=(addesr+addpre)/erre2;
  				      for (i=0;i<3;i++) {
  								fxx=repand*R3[i];
  								/*fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);*/
  								I->FTemp[i+NDIM*(ia1)] += fxx;
  								I->FTemp[i+NDIM*(ia2)] -= fxx;
  				      }
  				    }
  				  }
  				}
        }
      }
    }
  } else {
    esr = 0.0;
    for (is1=0;is1 < I->Max_Types; is1++) 
      for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
        for (i=0; i< NDIM; i++)
          I->FTemp[i+NDIM*(ia1)] = 0.;
  }
  SpeedMeasure(P, EwaldTime, StopTimeDo);
  for (is=0;is < I->Max_Types; is++) {
    MPI_Allreduce (I->FTemp , I->I[is].FEwald, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);
  }
  MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
}

/** Sum up all forces acting on Ions.
 * IonType::FIon = IonType::FEwald  + IonType::FIonL + IonType::FIonNL for all
 * dimensions #NDIM and each Ion of IonType
 * \param *P Problem at hand
 */
void CalculateIonForce(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  int i,j,d;
  for (i=0; i < I->Max_Types; i++) 
    for (j=0; j < I->I[i].Max_IonsOfType; j++)
      for (d=0; d<NDIM;d++)
				I->I[i].FIon[d+j*NDIM] = I->I[i].FEwald[d+j*NDIM] + I->I[i].FIonL[d+j*NDIM] + I->I[i].FIonNL[d+j*NDIM];
}

/** Write Forces to FileData::ForcesFile.
 * goes through all Ions per IonType and each dimension of #NDIM and prints in one line:
 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL, 
 * non-local IonType::FIonNL, ewald force IonType::FEwald
 * \param *P Problem at hand
 */
void OutputIonForce(struct Problem *P)
{
  struct RunStruct *R = &P->R;
  struct FileData *F = &P->Files;
  struct Ions *I = &P->Ion;
  FILE *fout = NULL;
  int i,j;
  if (F->MeOutMes != 1) return;
  fout = F->ForcesFile;
  for (i=0; i < I->Max_Types; i++) 
    for (j=0; j < I->I[i].Max_IonsOfType; j++) 
      fprintf(fout, "%i\t%i\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", i,j, 
	      I->I[i].R[0+j*NDIM], I->I[i].R[1+j*NDIM], I->I[i].R[2+j*NDIM],
	      I->I[i].FIon[0+j*NDIM], I->I[i].FIon[1+j*NDIM], I->I[i].FIon[2+j*NDIM],
	      I->I[i].FIonL[0+j*NDIM], I->I[i].FIonL[1+j*NDIM], I->I[i].FIonL[2+j*NDIM],
	      I->I[i].FIonNL[0+j*NDIM], I->I[i].FIonNL[1+j*NDIM], I->I[i].FIonNL[2+j*NDIM],
	      I->I[i].FEwald[0+j*NDIM], I->I[i].FEwald[1+j*NDIM], I->I[i].FEwald[2+j*NDIM]);
  fflush(fout);
  fprintf(fout, "MeanForce:\t%e\n", R->MeanForce[0]);
} 

/** Frees memory IonType.
 * All pointers from IonType and the pointer on it are free'd
 * \param *I Ions structure to be free'd
 * \sa IonsInitRead where memory is allocated
 */
void RemoveIonsRead(struct Ions *I)
{
  int i, it;
  for (i=0; i < I->Max_Types; i++) {
    //fprintf(stderr,"Name\n");
    Free(I->I[i].Name);
    //fprintf(stderr,"Abbreviation\n");
    Free(I->I[i].Symbol); 
    for (it=0;it<I->I[i].Max_IonsOfType;it++) {
    //fprintf(stderr,"sigma[it]\n");
      Free(I->I[i].sigma[it]);
    //fprintf(stderr,"sigma_rezi[it]\n");
      Free(I->I[i].sigma_rezi[it]);
    //fprintf(stderr,"sigma_PAS[it]\n");
      Free(I->I[i].sigma_PAS[it]);
    //fprintf(stderr,"sigma_rezi_PAS[it]\n");
      Free(I->I[i].sigma_rezi_PAS[it]);
    }
    //fprintf(stderr,"sigma\n");
    Free(I->I[i].sigma);
    //fprintf(stderr,"sigma_rezi\n");
    Free(I->I[i].sigma_rezi);
    //fprintf(stderr,"sigma_PAS\n");
    Free(I->I[i].sigma_PAS);
    //fprintf(stderr,"sigma_rezi_PAS\n");
    Free(I->I[i].sigma_rezi_PAS);
    //fprintf(stderr,"R\n");
    Free(I->I[i].R);
    //fprintf(stderr,"R_old\n");
    Free(I->I[i].R_old);
    //fprintf(stderr,"R_old_old\n");
    Free(I->I[i].R_old_old);
    //fprintf(stderr,"FIon\n");
    Free(I->I[i].FIon);
    //fprintf(stderr,"FIon_old\n");
    Free(I->I[i].FIon_old);
    //fprintf(stderr,"SearchDir\n");
    Free(I->I[i].SearchDir);
    //fprintf(stderr,"GammaA\n");
    Free(I->I[i].GammaA);
    //fprintf(stderr,"FIonL\n");
    Free(I->I[i].FIonL);
    //fprintf(stderr,"FIonNL\n");
    Free(I->I[i].FIonNL);
    //fprintf(stderr,"U\n");
    Free(I->I[i].U);
    //fprintf(stderr,"SFactor\n");
    Free(I->I[i].SFactor);
    //fprintf(stderr,"IMT\n");
    Free(I->I[i].IMT);
    //fprintf(stderr,"FEwald\n");
    Free(I->I[i].FEwald);
    Free(I->I[i].alpha);
  }
  //fprintf(stderr,"RLatticeVec\n");
  if (I->R_cut == 0.0)
    Free(I->RLatticeVec);
  //fprintf(stderr,"FTemp\n");
  Free(I->FTemp);
  //fprintf(stderr,"I\n");
  Free(I->I);
}

/** Shifts center of gravity of ion forces IonType::FIon.
 * First sums up all forces of movable ions to a "force temperature",
 * then reduces each force by this temp, so that all in all
 * the net force is 0.
 * \param *P Problem at hand
 * \sa CorrectVelocity
 * \note why is FTemp not divided by number of ions: is probably correct, but this
 *       function was not really meaningful from the beginning.
 */
void CorrectForces(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  double *FIon;
  double FTemp[NDIM];
  int is, ia, d;
  for (d=0; d<NDIM; d++)
    FTemp[d] = 0;
  for (is=0; is < I->Max_Types; is++) {		// calculate force temperature for each type ...
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {	// .. and each ion of this type ...
      FIon = &I->I[is].FIon[NDIM*ia];
      if (I->I[is].IMT[ia] == MoveIon) {		// .. if it's movable
				for (d=0; d<NDIM; d++) 
			 	  FTemp[d] += FIon[d];
      }
    }
  }
  for (is=0; is < I->Max_Types; is++) {		// and then reduce each by this value
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
      FIon = &I->I[is].FIon[NDIM*ia];
      if (I->I[is].IMT[ia] == MoveIon) {
				for (d=0; d<NDIM; d++) 
			 	 FIon[d] -= FTemp[d];						// (?) why not -= FTemp[d]/I->Max_TotalIons
      }
    }
  }
}


/** Shifts center of gravity of ion velocities (rather momentums).
 * First sums up ion speed U times their IonMass (summed up for each
 * dimension), then reduces the velocity by temp per Ions::TotalMass (so here
 * total number of Ions is included)
 * \param *P Problem at hand
 * \sa CorrectForces
 */
void CorrectVelocity(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  double *U;
  double UTemp[NDIM];
  int is, ia, d;
  for (d=0; d<NDIM; d++)
    UTemp[d] = 0;
  for (is=0; is < I->Max_Types; is++) {		// all types ...
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {	// ... all ions per type ...
      U = &I->I[is].U[NDIM*ia];
      if (I->I[is].IMT[ia] == MoveIon) {	// ... if it's movable
				for (d=0; d<NDIM; d++) 
				  UTemp[d] += I->I[is].IonMass*U[d];
      }
    }
  }
  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) {
				for (d=0; d<NDIM; d++) 
				  U[d] -= UTemp[d]/I->TotalMass;		// shift by mean velocity over mass and number of ions
      }
    }
  }
}

/** Moves ions according to calculated force.
 * Goes through each type of IonType and each ion therein, takes mass and
 * Newton to move each ion to new coordinates IonType::R, while remembering the last
 * two coordinates and the last force of the ion. Coordinates R are being
 * transformed to inverted base, shifted by +-1.0 and back-transformed to
 * real base.
 * \param *P Problem at hand
 * \sa CalculateIonForce
 * \sa UpdateIonsU
 * \warning U is not changed here, only used to move the ion.
 */
void UpdateIonsR(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  int is, ia, d;
  double *R, *R_old, *R_old_old, *FIon, *FIon_old, *U, *FIonL, *FIonNL, *FEwald;
  double IonMass, a;
  double sR[NDIM], sRold[NDIM], sRoldold[NDIM];
  const int delta_t = P->R.delta_t;
  for (is=0; is < I->Max_Types; is++) {	// go through all types ...
    IonMass = I->I[is].IonMass;
    if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsR: IonMass = %lg",IonMass);
    a = delta_t*0.5/IonMass;		// F/m = a and thus: s = 0.5 * F/m * t^2 + v * t =: t * (F * a + v)
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
      FIon = &I->I[is].FIon[NDIM*ia];
      FIonL = &I->I[is].FIonL[NDIM*ia];
      FIonNL = &I->I[is].FIonNL[NDIM*ia];
      FEwald = &I->I[is].FEwald[NDIM*ia];
      FIon_old = &I->I[is].FIon_old[NDIM*ia];
      U = &I->I[is].U[NDIM*ia];
      R = &I->I[is].R[NDIM*ia];
      R_old = &I->I[is].R_old[NDIM*ia];
      R_old_old = &I->I[is].R_old_old[NDIM*ia];
      if (I->I[is].IMT[ia] == MoveIon) {
				for (d=0; d<NDIM; d++) {
				  R_old_old[d] = R_old[d];				// shift old values
				  R_old[d] = R[d];
				  R[d] += delta_t*(U[d]+a*FIon[d]);			// F = m * a and s = 0.5 * a * t^2
				  FIon_old[d] = FIon[d];					// store old force
				  FIon[d] = 0;										// zero all as a sign that's been moved
				  FIonL[d] = 0;
				  FIonNL[d] = 0;
				  FEwald[d] = 0;
				}
				RMat33Vec3(sR, P->Lat.InvBasis, R);
				RMat33Vec3(sRold, P->Lat.InvBasis, R_old);
				RMat33Vec3(sRoldold, P->Lat.InvBasis, R_old_old);
				for (d=0; d <NDIM; d++) {
				  while (sR[d] < 0) {
				    sR[d] += 1.0;
				    sRold[d] += 1.0;
				    sRoldold[d] += 1.0;
				  }
				  while (sR[d] >= 1.0) {
				    sR[d] -= 1.0;
				    sRold[d] -= 1.0;
				    sRoldold[d] -= 1.0;
				  }
				}
				RMat33Vec3(R, P->Lat.RealBasis, sR);
				RMat33Vec3(R_old, P->Lat.RealBasis, sRold);
				RMat33Vec3(R_old_old, P->Lat.RealBasis, sRoldold);
      }
    }
  }
}

/** Resets the forces array.
 * \param *P Problem at hand
 */
void ResetForces(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  double *FIon, *FIonL, *FIonNL, *FEwald;
  int is, ia, d;
  for (is=0; is < I->Max_Types; is++) { // go through all types ...
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
      FIon = &I->I[is].FIon[NDIM*ia];
      FIonL = &I->I[is].FIonL[NDIM*ia];
      FIonNL = &I->I[is].FIonNL[NDIM*ia];
      FEwald = &I->I[is].FEwald[NDIM*ia];
      for (d=0; d<NDIM; d++) {
        FIon[d] = 0.;                   // zero all as a sign that's been moved
        FIonL[d] = 0.;
        FIonNL[d] = 0.;
        FEwald[d] = 0.;
      }
    }
  }
};

/** Changes ion's velocity according to acting force.
 * IonType::U is changed by the weighted force of actual step and last one
 * according to Newton.
 * \param *P Problem at hand
 * \sa UpdateIonsR
 * \sa CalculateIonForce
 * \warning R is not changed here.
 */
void UpdateIonsU(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  int is, ia, d;
  double *FIon, *FIon_old, *U;
  double IonMass, a;
  const int delta_t = P->R.delta_t;
  for (is=0; is < I->Max_Types; is++) {
    IonMass = I->I[is].IonMass;
    if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsU: IonMass = %lg",IonMass);
    a = delta_t*0.5/IonMass;				// (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
      FIon = &I->I[is].FIon[NDIM*ia];
      FIon_old = &I->I[is].FIon_old[NDIM*ia];
      U = &I->I[is].U[NDIM*ia];
      if (I->I[is].IMT[ia] == MoveIon)
				for (d=0; d<NDIM; d++) {
			  		U[d] += a * (FIon[d]+FIon_old[d]);
				}
    }
  }
}

/** CG process to optimise structure.
 * \param *P Problem at hand
 */
void UpdateIons(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  int is, ia, d;
  double *R, *R_old, *R_old_old, *FIon, *S, *GammaA;
  double IonFac, GammaB; //, GammaT;
  double FNorm, StepNorm;
  for (is=0; is < I->Max_Types; is++) {
    IonFac = I->I[is].IonFac;
    GammaA = I->I[is].GammaA;
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
      FIon = &I->I[is].FIon[NDIM*ia];
      S = &I->I[is].SearchDir[NDIM*ia];
      GammaB = GammaA[ia];
      GammaA[ia] = RSP3(FIon,S);
      FNorm = sqrt(RSP3(FIon,FIon));
      StepNorm = log(1.+IonFac*FNorm); // Fix Hack
      if (FNorm != 0)
        StepNorm /= sqrt(RSP3(FIon,FIon));
      else
        StepNorm = 0;
      //if (GammaB != 0.0) {
      //        GammaT = GammaA[ia]/GammaB;
      //        for (d=0; d>NDIM; d++)
      //          S[d] = FIon[d]+GammaT*S[d];
      //      } else {
        for (d=0; d<NDIM; d++)
          S[d] = StepNorm*FIon[d];
        //      }
      R = &I->I[is].R[NDIM*ia];
      R_old = &I->I[is].R_old[NDIM*ia];
      R_old_old = &I->I[is].R_old_old[NDIM*ia];
      if (I->I[is].IMT[ia] == MoveIon)
        for (d=0; d<NDIM; d++) {
          R_old_old[d] = R_old[d];
          R_old[d] = R[d];
          R[d] += S[d]; // FixHack*IonFac;
        }
    }
  }
}

/** Print coordinates of all ions to stdout.
 * \param *P Problem at hand
 */
void OutputIonCoordinates(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  int is, ia;
  if (P->Par.me == 0) { 
    fprintf(stderr, "(%i) ======= Updated Ion Coordinates =========\n",P->Par.me);
    for (is=0; is < I->Max_Types; is++)
      for (ia=0; ia < I->I[is].Max_IonsOfType; ia++)
        //fprintf(stderr, "(%i) R[%i/%i][%i/%i] = (%e,%e,%e)\n", P->Par.me, is, I->Max_Types, ia, I->I[is].Max_IonsOfType, I->I[is].R[NDIM*ia+0],I->I[is].R[NDIM*ia+1],I->I[is].R[NDIM*ia+2]);  
        fprintf(stderr, "Ion_Type%i_%i\t%.6f\t%.6f\t%.6f\t0\t# Atom from StructOpt\n", is+1, ia+1, I->I[is].R[NDIM*ia+0],I->I[is].R[NDIM*ia+1],I->I[is].R[NDIM*ia+2]);
    fprintf(stderr, "(%i) =========================================\n",P->Par.me);
  }
}

/** Calculates kinetic energy (Ions::EKin) of movable Ions.
 * Does 0.5 / IonType::IonMass * IonTye::U^2 for each Ion, 
 * also retrieves actual temperatur by 2/3 from just
 * calculated Ions::EKin.
 * \param *P Problem at hand
 */
void CalculateEnergyIonsU(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  int is, ia, d;
  double *U;
  double IonMass, a, ekin = 0;
  for (is=0; is < I->Max_Types; is++) {
    IonMass = I->I[is].IonMass;
    a = 0.5*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)
				for (d=0; d<NDIM; d++) {
			  		ekin += a * U[d]*U[d];
				}
    }
  }
  I->EKin = ekin;
  I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
}

/**	Scales ion velocities to match temperature.
 * In order to match Ions::ActualTemp with desired RunStruct::TargetTemp
 * each velocity in each dimension (for all types, all ions) is scaled
 * with the ratio of these two. Ions::EKin is recalculated.
 * \param *P Problem at hand
 */
void ScaleTemp(struct Problem *P)
{
  struct Ions *I = &P->Ion;
  int is, ia, d;
  double *U;
  double IonMass, a, ekin = 0;
  if (I->ActualTemp < MYEPSILON) fprintf(stderr,"ScaleTemp: I->ActualTemp = %lg",I->ActualTemp);
  double ScaleTempFactor = sqrt(P->R.TargetTemp/I->ActualTemp);
  for (is=0; is < I->Max_Types; is++) {
    IonMass = I->I[is].IonMass;
    a = 0.5*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)
				for (d=0; d<NDIM; d++) {
				  U[d] *= ScaleTempFactor;
				  ekin += a * U[d]*U[d];
				}
    }
  }
  I->EKin = ekin;
  I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
}

/** Calculates mean force vector as stop criteria in structure optimization.
 * Calculates a mean force vector per ion as the usual euclidian norm over total
 * number of ions and dimensions, being the sum of each ion force (for all type,
 * all ions, all dims) squared.
 * The mean force is archived in RunStruct::MeanForce and printed to screen.
 * \param *P Problem at hand
 */
void GetOuterStop(struct Problem *P)
{
  struct RunStruct *R = &P->R;
  struct Ions *I = &P->Ion;
  int is, ia, IonNo=0, i, d;
  double MeanForce = 0.0;
  for (is=0; is < I->Max_Types; is++) 
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
      IonNo++;
      for (d=0; d <NDIM; d++)
				MeanForce += I->I[is].FIon[d+NDIM*ia]*I->I[is].FIon[d+NDIM*ia];
    }
  for (i=MAXOLD-1; i > 0; i--) 
    R->MeanForce[i] = R->MeanForce[i-1];
  MeanForce = sqrt(MeanForce/(IonNo*NDIM));
  R->MeanForce[0] = MeanForce;
  if (P->Call.out[LeaderOut] && (P->Par.me == 0))
    fprintf(stderr,"MeanForce = %e\n", R->MeanForce[0]);
}