/** \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 <gsl/gsl_randist.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 Thermostat related values from parameter file.
 * \param *P Problem at hand
 * \param *source parameter file
 */
void InitThermostats(struct Problem *P, FILE *source)
{
  struct FileData *F = &P->Files;
  struct Ions *I = &P->Ion;
  char *thermo = MallocString(12, "IonsInitRead: thermo");
  int j;

  // initialize the naming stuff and all
  I->ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented");
  I->ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames");
  for (j=0;j<MaxThermostats;j++)
    I->ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "IonsInitRead: ThermostatNames[]");
  strcpy(I->ThermostatNames[0],"None");
  I->ThermostatImplemented[0] = 1;
  strcpy(I->ThermostatNames[1],"Woodcock"); 
  I->ThermostatImplemented[1] = 1;
  strcpy(I->ThermostatNames[2],"Gaussian"); 
  I->ThermostatImplemented[2] = 1;
  strcpy(I->ThermostatNames[3],"Langevin"); 
  I->ThermostatImplemented[3] = 1;
  strcpy(I->ThermostatNames[4],"Berendsen"); 
  I->ThermostatImplemented[4] = 1;
  strcpy(I->ThermostatNames[5],"NoseHoover"); 
  I->ThermostatImplemented[5] = 1;
  I->Thermostat = -1;
  if (F->MeOutMes) fprintf(F->TemperatureFile, "# %s\t%s\t%s --- ","Step", "ActualTemp(1)", "ActualTemp(2)");

  // read desired Thermostat from file along with needed additional parameters
  if (ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    if (strcmp(thermo, I->ThermostatNames[0]) == 0) { // None
      if (I->ThermostatImplemented[0] == 1) {
        I->Thermostat = None;
      } else {
        if (P->Par.me == 0)
          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[0]);
        I->Thermostat = None;
      }
    } else if (strcmp(thermo, I->ThermostatNames[1]) == 0) { // Woodcock
      if (I->ThermostatImplemented[1] == 1) {
        I->Thermostat = Woodcock;
        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, int_type, &P->R.ScaleTempStep, 1, critical); // read scaling frequency
      } else {
        if (P->Par.me == 0)
          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[1]);
        I->Thermostat = None;
      }
    } else if (strcmp(thermo, I->ThermostatNames[2]) == 0) { // Gaussian
      if (I->ThermostatImplemented[2] == 1) {
        I->Thermostat = Gaussian;
        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, int_type, &P->R.ScaleTempStep, 1, critical); // read collision rate
      } else {
        if (P->Par.me == 0)
          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[2]);
        I->Thermostat = None;
      }
    } else if (strcmp(thermo, I->ThermostatNames[3]) == 0) { // Langevin
      if (I->ThermostatImplemented[3] == 1) {
        I->Thermostat = Langevin;
        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.TempFrequency, 1, critical); // read gamma
        if (ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 3, 1, double_type, &P->R.alpha, 1, optional)) {
          if (P->Par.me == 0) 
            fprintf(stderr,"(%i) Extended Stochastic Thermostat detected with interpolation coefficient %lg\n", P->Par.me, P->R.alpha);
        } else {
          P->R.alpha = 1.;
        }
      } else {
        if (P->Par.me == 0)
          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[3]);
        I->Thermostat = None;
      }
    } else if (strcmp(thermo, I->ThermostatNames[4]) == 0) { // Berendsen
      if (I->ThermostatImplemented[4] == 1) {
        I->Thermostat = Berendsen;
        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.TempFrequency, 1, critical); // read \tau_T
      } else {
        if (P->Par.me == 0)
          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[4]);
        I->Thermostat = None;
      }
    } else if (strcmp(thermo, I->ThermostatNames[5]) == 0) { // Nose-Hoover
      if (I->ThermostatImplemented[5] == 1) {
        I->Thermostat = NoseHoover;
        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.HooverMass, 1, critical); // read Hoovermass
        P->R.alpha = 0.;
      } else {
        if (P->Par.me == 0)
          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[5]);
        I->Thermostat = None;
      }
    }
    if (I->Thermostat == -1) {
      if (P->Par.me == 0) 
        fprintf(stderr,"(%i) Warning: thermostat name was not understood!\n", P->Par.me);
      I->Thermostat = None;
    }
  } else {
    if ((P->R.MaxOuterStep > 0) && (P->R.TargetTemp != 0)) 
      debug(P, "No thermostat chosen despite finite temperature MD, falling back to None.");
    I->Thermostat = None;
  }
  if (F->MeOutMes) fprintf(F->TemperatureFile, "%s Thermostat\n",I->ThermostatNames[(int)I->Thermostat]); 
  Free(thermo, "InitThermostats: thermo");
}

/** Parses for a given keyword the ion coordinates and velocites.
 * \param *P Problem at hand, contains pointer to Ion and IonType structure
 * \param *source
 * \param *keyword keyword string
 * \param repetition which repeated version of the keyword shall be read by ReadParameters()
 * \param relative whether atom coordinates are relative to unit cell size or absolute
 * \param Bohr conversion factor to atomiclength (e.g. 1.0 if coordinates in bohr radius, 0.529... if given in Angstrom)
 * \param sigma variance of maxwell-boltzmann distribution
 * \param *R where to put the read coordinates
 * \param *U where to put the read velocities
 * \param *IMT whether its MoveType#FixedIon or not
 * \return if 1 - success, 0 - failure
 */
int ParseIonsCoordinatesAndVelocities(struct Problem *P, FILE *source, char *keyword, int repetition, int relative, double sigma, double *R, double *U, int *IMT) 
{
  //struct RunStruct *Run = &P->R;
  struct Ions *I = &P->Ion;
  struct Lattice *L = &P->Lat;
  int k;
  double R_trafo[3];
  gsl_rng * r;
  const gsl_rng_type * T;
  double Bohr = (I->IsAngstroem) ? ANGSTROEMINBOHRRADIUS : 1.;
  
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
  if (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 1, 1, double_type, &R_trafo[0], repetition, optional) &&
    ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 2, 1, double_type, &R_trafo[1], repetition, optional) &&
    ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 3, 1, double_type, &R_trafo[2], repetition, optional) &&
    ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 4, 1, int_type, IMT, repetition, optional)) {
    // if read successful, then try also parsing velocities
    if ((ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 5, 1, double_type, &U[0], repetition, optional)) &&
      (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 6, 1, double_type, &U[1], repetition, optional)) && 
      (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 7, 1, double_type, &U[2], repetition, optional))) {
      //nothing
    } else {  // if no velocities specified, set to maxwell-boltzmann distributed
      if ((I->Thermostat != None)) {
        if (P->Par.me == 0) fprintf(stderr, "(%i) Missing velocities of ion %s: Maxwell-Boltzmann with %lg\n", P->Par.me, keyword, sigma);
        for (k=0;k<NDIM;k++)
          U[k] = gsl_ran_gaussian (r, sigma);
      } else {
        for (k=0;k<NDIM;k++)
          U[k] = 0.;
      }
    }
    // change if coordinates were relative
    if (relative) {
      //fprintf(stderr,"(%i)Ion coordinates are relative %i ... \n",P->Par.me, relative);
      RMat33Vec3(R, L->RealBasis, R_trafo); // multiply with real basis
    } else {
      for (k=0;k<NDIM;k++) // otherweise just copy to vector
        R[k] = R_trafo[k];
    }
    // check if given MoveType is valid
    if (*IMT < 0 || *IMT >= MaxIonMoveType) {
      fprintf(stderr,"Bad Ion MoveType set to MoveIon for %s\n", keyword);
      *IMT = 0;
    }
    if (!relative) {
      //fprintf(stderr,"converting with %lg\n", Bohr);
      SM(R, Bohr, NDIM); // convert angstroem to bohrradius
      SM(U, Bohr, NDIM); // convert angstroem to bohrradius
    }
    // finally periodicly correct coordinates to remain in unit cell
    RMat33Vec3(R_trafo, L->InvBasis, R);
    for (k=0; k <NDIM; k++) { // periodic correction until within unit cell
      while (R_trafo[k] < 0) 
        R_trafo[k] += 1.0;
      while (R_trafo[k] >= 1.0) 
        R_trafo[k] -= 1.0;
    }
    RMat33Vec3(R, L->RealBasis, R_trafo);
    // print coordinates if desired
    if ((P->Call.out[ReadOut])) {
      fprintf(stderr,"(%i) coordinates of Ion %s: (x,y,z) = (",P->Par.me, keyword);
      for(k=0;k<NDIM;k++) {
        fprintf(stderr,"%lg ",R[k]);
        if (k != NDIM-1) fprintf(stderr,", ");
        else fprintf(stderr,")\n");
      }
    }
    gsl_rng_free (r);
    return 1;
  } else {
    for (k=0;k<NDIM;k++) {
      U[k] = R[k] = 0.;
    }
    *IMT = 0; 
    gsl_rng_free (r);
    return 0;
  }
}

/** 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,d,me, count;
  int relative; // whether read-in coordinate are relative (1) or absolute  
  double Bohr = ANGSTROEMINBOHRRADIUS; /* Angstroem */
  double sigma, R[3], U[3];
  struct IonConstrained *ptr = NULL;
     
  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, 1, critical);
  ParseForParameter(P->Call.out[ReadOut],source,"IsAngstroem", 0, 1, 1, int_type, &I->IsAngstroem, 1, 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, 1, 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, 1, optional))
    relative = 0;
  if (!ParseForParameter(P->Call.out[ReadOut],source,"StructOpt", 0, 1, 1, int_type, &I->StructOpt, 1, optional))
    I->StructOpt = 0;
  
  InitThermostats(P, source);
  
  /* 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, 1, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, int_type, &I->I[i].Z, 1, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &I->I[i].rgauss, 1, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &I->I[i].l_max, 1, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 5, 1, int_type, &I->I[i].l_loc, 1, critical);
		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 6, 1, double_type, &I->I[i].IonMass, 1, critical);
    if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 7, 1, string_type, &I->I[i].Name[0], 1, 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], 1, 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].moment = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: moment");
    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].moment[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: moment");
      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].IonMass *= Units2Electronmass;  // in config given in amu not in "atomic units" (electron mass, m_e = 1)
                                            // proportional factor between electron and proton mass
    I->I[i].IonFac = Units2Electronmass/I->I[i].IonMass;
    I->TotalMass += I->I[i].Max_IonsOfType*I->I[i].IonMass;
    sigma = sqrt(P->R.TargetTemp/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].Constraints = (struct IonConstrained **) Malloc(sizeof(struct IonConstrained *)*I->I[i].Max_IonsOfType, "IonsInitRead: **Constraints");
    for (it=0;it<I->I[i].Max_IonsOfType;it++)
      I->I[i].Constraints[it] = NULL;
    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].FMagnetic = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FMagnetic");
    I->I[i].FConstraint = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FConstraint");
    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");
    // now parse ion coordination
    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);
      for (d=0; d <NDIM; d++) // transfor to old and old_old coordinates as well
        I->I[i].FMagnetic[d+j*NDIM] = 0.;
      // first ones are starting positions
      if ((count = ParseIonsCoordinatesAndVelocities(P, source, name, 1, relative, sigma, &I->I[i].R[NDIM*j], &I->I[i].U[NDIM*j], &IMT))) {
        for (d=0; d <NDIM; d++) // transfor to old and old_old coordinates as well
          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->I[i].IMT[j] = IMT;
        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) R[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,I->I[i].R[0+NDIM*j],I->I[i].R[1+NDIM*j],I->I[i].R[2+NDIM*j]);
        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) U[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,I->I[i].U[0+NDIM*j],I->I[i].U[1+NDIM*j],I->I[i].U[2+NDIM*j]);
        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) IMT[%d,%d] = %i\n", P->Par.me, i,j,I->I[i].IMT[j]);
      } else {
        Error(SomeError, "Could not find or parse coordinates of an ion");
      }
      // following ones are constrained motions
      while (ParseIonsCoordinatesAndVelocities(P, source, name, ++count, relative, sigma, R, U, &IMT)) {
        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) Constrained read %d ...\n", P->Par.me, count-1);
        ptr = AddConstraintItem(&I->I[i], j);  // allocate memory
        for (d=0; d< NDIM; d++) { // fill the constraint item
          ptr->R[d] = R[d];
          ptr->U[d] = U[d];
        }
        ptr->IMT = IMT;
        ptr->step = count-1;
        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) R[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,ptr->R[0],ptr->R[1],ptr->R[2]);
        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) U[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,ptr->U[0],ptr->U[1],ptr->U[2]);
        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) IMT[%d,%d] = %i\n", P->Par.me, i,j,ptr->IMT);
      }
      //if (P->Call.out[ReadOut]) 
      fprintf(stderr,"(%i) A total of %d additional constrained moves for ion (%d,%d) was parsed.\n", P->Par.me, count-2, i,j);

      I->Max_TotalIons++;
    }
  }
  Free(free_name, "IonsInitRead: 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:");
}

/** Allocates memory for a IonConstrained item and adds to end of list.
 * \param *I IonType structure
 * \param ion which ion of the type
 * \return pointer to created item
 */
struct IonConstrained * AddConstraintItem(struct IonType *I, int ion)
{
  struct IonConstrained **ptr = &(I->Constraints[ion]);
  while (*ptr != NULL) { // advance to end of list
    ptr = &((*ptr)->next);
  }
  // allocated memory for structure and items within
  *ptr = (struct IonConstrained *) Malloc(sizeof(struct IonConstrained), "CreateConstraintItem: IonConstrained");
  (*ptr)->R = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *R");
  (*ptr)->U = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *U");
  (*ptr)->next = NULL;
  return *ptr;
}

/** Removes first of item of IonConstrained list for given \a ion.
 * Frees memory of items within and then
 * \param *I IonType structure
 * \param ion which ion of the type
 * \return 1 - success, 0 - failure (list is already empty, start pointer was NULL)
 */ 
int RemoveConstraintItem(struct IonType *I, int ion)
{
  struct IonConstrained **ptr = &(I->Constraints[ion]); 
  struct IonConstrained *next_ptr;
  if (*ptr != NULL) {
    next_ptr = (*ptr)->next;
    Free((*ptr)->R, "RemoveConstraintItem: R");
    Free((*ptr)->U, "RemoveConstraintItem: U");
    Free((*ptr), "RemoveConstraintItem: IonConstrained structure");
    *ptr = next_ptr; // set start of list to next item (automatically is null if ptr was last item)
    return 1;
  } else {
    return 0;
  }
}

/** Calculated Ewald force.
 * The basic idea of the ewald method is to add a countercharge - here of gaussian form -
 * which shields the long-range charges making them effectively short-ranged, while the
 * countercharges themselves can be analytically (here by (hidden) Fouriertransform: it's
 * added to the electron density) integrated and withdrawn. 
 * \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 if != 0 table mirror cells to take into (local!) account of ewald summation 
 */
void CalculateEwald(struct Problem *P, int first) 
{
  int r,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec; //,is
  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 (first) {
    SpeedMeasure(P, InitSimTime, StopTimeDo);
    rcut2 = I->R_cut*I->R_cut;
    ActualVec = 0;  // counts number of cells taken into account
    MaxCell = (int)(5.*I->R_cut/pow(L->Volume,1./3.));    // the number of cells in each direction is prop. to axis length over cutoff
    if (MaxCell < 2) MaxCell = 2;     // take at least 2
    for (i = -MaxCell; i <= MaxCell; i++) {
      for (j = -MaxCell; j <= MaxCell; j++) {
				for (k = -MaxCell; k <= MaxCell; k++) {
				  r2 = 0.;
				  for (ir=0; ir <NDIM; ir++) {  // t is the offset vector pointing to distant cell (only diagonal entries)
				    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 squared distance to other cell 
				  }
				  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));   // whether it's directly adjacent
				  if ((r2 <= rcut2) || Is_Neighb_Cell) {    // if it's either adjacent or within cutoff, count it in
				    ActualVec++;
				  }
				}
      }
    }
    MaxVec = ActualVec;   // store number of cells
    I->MaxVec = ActualVec;
    MaxLocalVec = MaxVec / MaxPar;  // share among processes
    StartVec = ParMe*MaxLocalVec; // for this process start at its patch
    r = MaxVec % MaxPar;  // rest not fitting...
    if (ParMe >= r) {   // ones who don't get extra cells have to shift their start vector
      StartVec += r;
    } else {    // others get extra cell and have to shift also by a bit
      StartVec += ParMe;
      MaxLocalVec++;
    }
    I->MaxLocalVec = MaxLocalVec;
    LocalActualVec = ActualVec = 0;
    // now go through the same loop again and note down every offset vector for cells in this process' patch
    if (MaxLocalVec != 0) {
      I->RLatticeVec = (double *)
        Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald: I->RLatticeVec");
    } else {
      fprintf(stderr,"(%i) Warning in Ewald summation: Got MaxLocalVec == 0\n", P->Par.me);
    }
    for (i = -MaxCell; i <= MaxCell; i++) {
      for (j = -MaxCell; j <= MaxCell; j++) {
				for (k = -MaxCell; k <= MaxCell; k++) {
				  r2 = 0.;
				  for (ir=0; ir <NDIM; ir++) {  // create offset vector from diagonal entries
				    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 adjacent or within cutoff ...
				    if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) { // ... and belongs to patch, store 'em
				      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);

  // first set everything to zero
  esr = 0.0;
  //for (is1=0;is1 < I->Max_Types; is1++) 
  // then take each ion of each type once 
  for (is1=0;is1 < I->Max_Types; is1++) {
    // reset FTemp for IonType
    for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
      for (i=0; i< NDIM; i++)
        I->FTemp[i+NDIM*(ia1)] = 0.;
    // then sum for each on of the IonType
    for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
      for (i=0;i<NDIM;i++) {
				R1[i]=I->I[is1].R[i+NDIM*ia1];    // R1 is local coordinate of current first ion
      }
      for (ir2=0;ir2<I->MaxLocalVec;ir2++) {  // for every cell of local patch (each with the same atoms)
				for (is2=0;is2 < I->Max_Types; is2++) { // and for each ion of each type a second time
          // gkl = 1./(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
				  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]; // R2 is coordinate of current second ion in the distant cell!
				    }

				    erre2=0.0;
				    for (i=0;i<NDIM;i++) {
				      R3[i] = R1[i]-R2[i];    // R3 = R_is1,ia1 - R_is2,ia2 - L
				      erre2+=R3[i]*R3[i];     // erre2 = |R_is1,ia1 - R_is2,ia2 - L|^2
				    }
				    if (erre2 > MYEPSILON) {  // appears as one over ... in fac, thus check if zero
              arg = sqrt(erre2); // arg = |R_is1,ia1 - R_is2,ia2 - L|
				      fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;    // fac = -1/2 * Z_is1 * Z_is2 / arg
				      
				      arg *= gkl;             // arg = |R_is1,ia1 - R_is2,ia2 - L|/(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
				      addesr = derf(arg);  
              addesr = (1.-addesr)*fac; // complementary error function: addesr = 1/2*erfc(arg)/|R_is1,ia1 - R_is2,ia2 - L|
				      esr += addesr;
              addpre=exp(-arg*arg)*gkl;
              addpre = PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;  // addpre = exp(-|R_is1,ia1 - R_is2,ia2 - L|^2/(r_is1,gauss^2 + r_is2,gauss^2))/(sqrt(pi)*sqrt(r_is1,gauss^2 + r_is2,gauss^2))
              repand = (addesr+addpre)/erre2;
//              if ((fabs(repand) > MYEPSILON)) {
//                fprintf(stderr,"(%i) Ewald[is%d,ia%d/is%d,ia%d,(%d)]: Z1 %lg, Z2 %lg\tpre %lg\t esr %lg\t fac %lg\t arg %lg\tR3 (%lg,%lg,%lg)\n", P->Par.me,  is1, ia1, is2, ia2, ir2, PP->zval[is1], PP->zval[is2],addpre,addesr, fac, arg, R3[0], R3[1], R3[2]);
//              }
				      for (i=0;i<NDIM;i++) {
								fxx=repand*R3[i];
//                if ((fabs(repand) > MYEPSILON)) {
//								  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*2.0;      // actio = reactio?
								//I->FTemp[i+NDIM*(ia2)] -= fxx;
				      }
				    }
				  }
				}
      }
    }
    MPI_Allreduce (I->FTemp , I->I[is1].FEwald, NDIM*I->I[is1].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);  
  }
  SpeedMeasure(P, EwaldTime, StopTimeDo);
  //for (is=0;is < I->Max_Types; is++) {
  //}
  MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
//  fprintf(stderr,"\n");
}

/** Sum up all forces acting on Ions.
 * IonType::FIon = IonType::FEwald  + IonType::FIonL + IonType::FIonNL + IonType::FMagnetic 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] + I->I[i].FMagnetic[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\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].FMagnetic[0+j*NDIM], I->I[i].FMagnetic[1+j*NDIM], I->I[i].FMagnetic[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]);
} 

/** Reads Forces from CallOptions::ForcesFile.
 * goes through all Ions per IonType and each dimension of #NDIM and parses 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 ParseIonForce(struct Problem *P)
{
	//struct RunStruct *Run = &P->R;
  struct Ions *I = &P->Ion;
  FILE *finput = NULL;
	char line[1024];
  int i,j;
  double i1, j1;
  double R[NDIM];
  
  if (P->Call.ForcesFile == NULL) return;
  finput = fopen(P->Call.ForcesFile, "r");
	//fprintf(stderr, "File pointer says %p\n", finput);
	if (finput == NULL)
		Error(InitReading, "ForcesFile does not exist.");
	fscanf(finput, "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n", line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line);	// skip header line
  for (i=0; i < I->Max_Types; i++) 
    for (j=0; j < I->I[i].Max_IonsOfType; j++)
			if (feof(finput)) {
				Error(InitReading, "ForcesFile does not contain enough lines.");
			} else {
				//fscanf(finput, "%s\n", line);
				//if (strlen(line) < 100) {
					//fprintf(stderr, "i %d, j %d, length of line %ld, line %s.\n", i,j, strlen(line), line);
					//Error(InitReading, "ForcesFile does not contain enough lines or line is faulty.");
				//} else
					//fprintf(stderr, "Parsed line: '%s'\n", line);
	      	fscanf(finput, "%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n", &i1,&j1, 
		      &R[0], &R[1], &R[2],
		      &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].FMagnetic[0+j*NDIM], &I->I[i].FMagnetic[1+j*NDIM], &I->I[i].FMagnetic[2+j*NDIM],
		      &I->I[i].FEwald[0+j*NDIM], &I->I[i].FEwald[1+j*NDIM], &I->I[i].FEwald[2+j*NDIM]);
				if ((i != (int)i1) || (j != (int)j1))
					Error(InitReading, "Line does not match to desired ion.");
		}
	fclose(finput);
  //fprintf(stderr, "MeanForce:\t%e\n", Run->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++) {
    Free(I->I[i].Name, "RemoveIonsRead: I->I[i].Name");
    Free(I->I[i].Symbol, "RemoveIonsRead: I->I[i].Symbol"); 
    for (it=0;it<I->I[i].Max_IonsOfType;it++) {
      Free(I->I[i].moment[it], "RemoveIonsRead: I->I[i].moment[it]");
      Free(I->I[i].sigma[it], "RemoveIonsRead: I->I[i].sigma[it]");
      Free(I->I[i].sigma_rezi[it], "RemoveIonsRead: I->I[i].sigma_rezi[it]");
      Free(I->I[i].sigma_PAS[it], "RemoveIonsRead: I->I[i].sigma_PAS[it]");
      Free(I->I[i].sigma_rezi_PAS[it], "RemoveIonsRead: I->I[i].sigma_rezi_PAS[it]");
      while (RemoveConstraintItem(&I->I[i], it)); // empty the constrained item list
    }
    Free(I->I[i].sigma, "RemoveIonsRead: I->I[i].sigma");
    Free(I->I[i].sigma_rezi, "RemoveIonsRead: I->I[i].sigma_rezi");
    Free(I->I[i].sigma_PAS, "RemoveIonsRead: I->I[i].sigma_PAS");
    Free(I->I[i].sigma_rezi_PAS, "RemoveIonsRead: I->I[i].sigma_rezi_PAS");
    Free(I->I[i].R, "RemoveIonsRead: I->I[i].R");
    Free(I->I[i].R_old, "RemoveIonsRead: I->I[i].R_old");
    Free(I->I[i].R_old_old, "RemoveIonsRead: I->I[i].R_old_old");
    Free(I->I[i].Constraints, "RemoveIonsRead: I->I[i].Constraints");
    Free(I->I[i].FIon, "RemoveIonsRead: I->I[i].FIon");
    Free(I->I[i].FIon_old, "RemoveIonsRead: I->I[i].FIon_old");
    Free(I->I[i].SearchDir, "RemoveIonsRead: I->I[i].SearchDir");
    Free(I->I[i].GammaA, "RemoveIonsRead: I->I[i].GammaA");
    Free(I->I[i].FIonL, "RemoveIonsRead: I->I[i].FIonL");
    Free(I->I[i].FIonNL, "RemoveIonsRead: I->I[i].FIonNL");
    Free(I->I[i].FMagnetic, "RemoveIonsRead: I->I[i].FMagnetic");
    Free(I->I[i].U, "RemoveIonsRead: I->I[i].U");
    Free(I->I[i].SFactor, "RemoveIonsRead: I->I[i].SFactor");
    Free(I->I[i].IMT, "RemoveIonsRead: I->I[i].IMT");
    Free(I->I[i].FEwald, "RemoveIonsRead: I->I[i].FEwald");
    Free(I->I[i].alpha, "RemoveIonsRead: I->I[i].alpha");
  }
  if (I->R_cut == 0.0)
    Free(I->RLatticeVec, "RemoveIonsRead: I->RLatticeVec");
  Free(I->FTemp, "RemoveIonsRead: I->FTemp");
  Free(I->I, "RemoveIonsRead: 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
 * \param actual (1 - don't at current RunStruct#OuterStep, 0 - do)
 */
void OutputIonCoordinates(struct Problem *P, int actual)
{
  //struct RunStruct *R = &P->R;
  struct Ions *I = &P->Ion;
  char filename[255];
  FILE *output;
  int is, ia, nr = 0;
  double Bohr = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
/*  if (P->Par.me == 0 && P->Call.out[ReadOut]) { 
    fprintf(stderr, "(%i) ======= Differential 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_old[NDIM*ia+0],I->I[is].R[NDIM*ia+1]-I->I[is].R_old[NDIM*ia+1],I->I[is].R[NDIM*ia+2]-I->I[is].R_old[NDIM*ia+2]);
    fprintf(stderr, "(%i) =========================================\n",P->Par.me);
  }*/
  if (actual)
    sprintf(filename, "%s.MD", P->Call.MainParameterFile); 
  else
    sprintf(filename, "%s.opt", P->Call.MainParameterFile);
  if (((actual) && (P->R.OuterStep <= 0)) || ((!actual) && (P->R.StructOptStep == 1))) { // on first step write complete config file 
    WriteParameters(P, filename);
  } else { // afterwards simply add the current ion coordinates as constrained steps
    if ((P->Par.me == 0) && (output = fopen(filename,"a"))) {
      fprintf(output,"# ====== MD step %d ========= \n", (actual) ? P->R.OuterStep : P->R.StructOptStep);
      for (is=0; is < I->Max_Types; is++)
        for (ia=0;ia<I->I[is].Max_IonsOfType;ia++) {
          fprintf(output,"Ion_Type%i_%i\t%2.9f\t%2.9f\t%2.9f\t%i\t%2.9f\t%2.9f\t%2.9f\t# Number in molecule %i\n", is+1, ia+1, I->I[is].R[0+NDIM*ia]*Bohr, I->I[is].R[1+NDIM*ia]*Bohr, I->I[is].R[2+NDIM*ia]*Bohr, I->I[is].IMT[ia], I->I[is].U[0+NDIM*ia]*Bohr, I->I[is].U[1+NDIM*ia]*Bohr, I->I[is].U[2+NDIM*ia]*Bohr, ++nr);
        }
      fflush(output);
    }
  }
}

/** 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) // even FixedIon moves, only not by other's forces
				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) // even FixedIon moves, only not by other's forces
				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]);
}