#ifndef ions_h
#define ions_h

/** \file ions.h
 * Header file for \ref ions.c
 * 
 * Contains declarations of the functions implemented in \ref ions.c,
 * enumerations such CoreCorrType - whether the config file has values for the
 * core electron density - or IonMoveType - is the ion allowed to move or not -
 * and structures for IonType, containg all Ions of a atomic type, which are all
 * attached to Ions.
 * 
  Project: ParallelCarParrinello
  Jan Hamaekers
  2000

  File: ions.h
  $Id: ions.h,v 1.19 2007/02/09 09:13:48 foo Exp $
*/


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

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

enum CoreCorrType { NotCoreCorrected, CoreCorrected};	 //!< Pseudopotential: Whether there is a correction factor for the potential ansatz or not
#define MaxIonMoveType 2									//!< Number of different IonMoveType
enum IonMoveType { MoveIon, FixedIon}; 		//!< Whether the Ion is movable or not (influences kinetic energy among others)

/** Structure containing data for a single ion type.
 * Contains Mass, atomic number, number of types, coordinates (current and the
 * last two old positions), forces (and last one), velocity, some Pseudopotential
 * values
 */
struct IonType { /* A single Ion */
  int Z;							//!< Atomic number of ion
  double IonMass;			//!< Mass of ion
  int Max_IonsOfType;	//!< Maximum number of ions of this type (needed for memory allocation)
  double *R;					//!< coordinate (x,y,z) array for all ions of this type R[i][j] (i=x,y,z;j=1,..)
  double *R_old;			//!< (?) Remembers last position (two-dim array, see IonType::R)
  double *R_old_old;	//!< (?) Remembers position before the last one (two-dim array, see IonType::R)
  double *FIon;				//!< (?) Overall force acting on the ion (two-dim array, see IonType::R)
  double *FIonL;			//!< (?) Linear force acting on the ion (two-dim array, see IonType::R)
  double *FIonNL;			//!< (?) Nonlinear force acting on the ion (two-dim array, see IonType::R)
  double *U; 					//!< Velocity of the ion (two-dim array, see IonType::R)
  double *FIon_old;		//!< (?) Remembers prior force acting on the ion (two-dim array, see IonType::R)
  double *FEwald;			//!< (?) Ewald force (two-dim array, see IonType::R)
  double IonFac;      //!< mass-like factor in CG structure optimization to determinte step width
  double *GammaA;     //!< 
  double *SearchDir;  //!< 
  double rgauss;			//!< Gauss radius
  double *alpha;      //!< parameter for CSDGT gauge, see ShiftGaugeOrigin()
  int l_max;					//!< (?) Pseudopotential specific value
  int l_loc;					//!< (?) Pseudopotential specific value
  double ZFactor; 		//!< (?) CalcVIFactor
  fftw_complex *SFactor;      //!< structure factor, \f$S_{I_s} (G) \f$, thereby the centre of a potential is shifted away from the origin
  enum CoreCorrType corecorr;	//!< (?) Core correction due to pseudo potential ansatz
  enum IonMoveType *IMT;		//!< Ion is moving or fixed
  char *Name;         //!< name of ion's element
  char *Symbol; //!< short form of ion's element
  double **sigma;     //!< shielding tensor
  double **sigma_rezi;      //!< shielding tensor for reciprocal calculation
  double **sigma_PAS;       //!< Principal Axis System shielding tensor
  double **sigma_rezi_PAS;  //!< Principal Axis System  shielding tensor for reciprocal calculation
  double chi[NDIM*NDIM];  //!< magnetic susceptibility
  double chi_PAS[NDIM];  //!< Principal Axis System  magnetic susceptibility
};

/** Containing structure for all ion data.
 * Contains max types, counts, cut values, temperature and total mass among others
 */
struct Ions { /* All Ions */
  int Max_Types;			//!< Number of types overall (needed for memory allocation)
  int Max_TotalIons;	//!< Maximum of of total ions (needed for memory allocation)
  int Max_Max_IonsOfType; //!< Maximum of ions per type (needed for memory allocation)
  int TotalZval;
  struct IonType *I;	//!< Structure containing data for each type
  /* Ewald */
  double R_cut;				//!< Radial cut value
  int MaxVec;
	int MaxLocalVec;
  double *RLatticeVec;
  double *FTemp;			//!< 
  double EKin;				//!< Kinetic energy of electrons
  double ActualTemp;	//!< Actual temperature
  double TotalMass;		//!< Total mass of all ions in the structure
  int StructOpt;      //!< whether structure optimization (1) (with CG) or MD (0) shall be done
};


/* Functions */
void IonsInitRead(struct Problem *P, FILE *source);
void CalculateEwald(struct Problem *P, int first);
void RemoveIonsRead(struct Ions *I);
/*void CalculateIonLocalForce(struct Problem *P);*/
void CalculateIonForce(struct Problem *P);
void OutputIonForce(struct Problem *P);
void OutputIonCoordinates(struct Problem *P);
void UpdateIons(struct Problem *P);
void UpdateIonsR(struct Problem *P);
void UpdateIonsU(struct Problem *P);
void CalculateEnergyIonsU(struct Problem *P);
void ScaleTemp(struct Problem *P);
void GetOuterStop(struct Problem *P);
void CorrectForces(struct Problem *P);
void CorrectVelocity(struct Problem *P);
#endif