/***********/
/* argon.h */
/***********/

#ifndef _ARGON_H
#define _ARGON_H

/** universal constants **/
#ifndef PI
#define PI 3.14159265358979323846
#endif
#define AVO   6.0221367E+23
#define BOLZ  1.380658E-23    /* in Joule */
#define HBAR  1.054589E-34    /* in Joule.s */

/** argon constants **/
#define ARGON_MASS  (39.948*1E-3/AVO) /* in kg */
#define DEBYE_TEMP  92 /* Debye temperature in Kelvin */

/** argon LJ6-12 constants **/
#define SIGMA     3.405E-10     /* in meter */
#define EPSILON   (119.8*BOLZ)  /* in Joule */

/** reduced units benchmarked to argon LJ6-12 **/
#define UENERGY_IN_JOULE EPSILON
#define ULENGTH_IN_METER SIGMA
#define UMASS_IN_KG ARGON_MASS     
/* from now on are all derived reduced units --> */
#define ULENGTH_IN_ANGSTROM (SIGMA*1E10)
#define UTIME_IN_SECOND sqrt(UMASS_IN_KG*ULENGTH_IN_METER*ULENGTH_IN_METER/UENERGY_IN_JOULE)
#define UTIME_IN_PS (UTIME_IN_SECOND/1.E-12) 
#define USTRESS_IN_PASCAL (UENERGY_IN_JOULE/ULENGTH_IN_METER/ULENGTH_IN_METER/ULENGTH_IN_METER)
#define USTRESS_IN_MPA (USTRESS_IN_PASCAL/1.E6)
#define USTRESS_IN_BAR (USTRESS_IN_PASCAL/1.E5)
#define UCONDUCTIVITY_IN_SI (UENERGY_IN_JOULE/ULENGTH_IN_METER/UTIME_IN_SECOND) /* in W/m/K */

/* Woon's ab initio pair potential */
#define EV_IN_JOULE 1.60217733E-19
#define HARTREE_IN_JOULE 4.3597482E-18
#define BOHR_RADIUS_IN_ANGSTROM 0.529177249
#define WOON_EPSILON (0.4228*HARTREE_IN_JOULE/1000./UENERGY_IN_JOULE)
#define WOON_RE (7.1714*BOHR_RADIUS_IN_ANGSTROM/ULENGTH_IN_ANGSTROM)
#define WOON_Q  0.6060
#define WOON_KE (0.6627*HARTREE_IN_JOULE/1000./UENERGY_IN_JOULE/BOHR_RADIUS_IN_ANGSTROM/BOHR_RADIUS_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM)
#define WOON_K1 sqrt((WOON_Q+1)*WOON_KE/WOON_Q/WOON_EPSILON)
#define WOON_K2 sqrt(WOON_Q*WOON_KE/(WOON_Q+1)/WOON_EPSILON)
/* D. E. Woon, Chem. Phys. Lett. 204, 29 (1993) */

/** predictor-corrector integration parameters **/
#define F02  (3./16.)
#define F12  (251./360.)
#define F32  (11./18.)
#define F42  (1./6.)
#define F52  (1./60.)

/** common varaibles **/
int np, neff; /* number of particles */
int nc[3];   /* number of unit cells in three dimensions */
double TR;   /* real temperature in Kelvin */
double T, dTdR; /* MD temperature in Kelvin and its derivative */
double H[3][3]; /* shape of the box: H[0][*] is the first edge */
double HI[3][3], volume;  /* current supercell volume in reduced units */
double H1[3][3], H2[3][3], H3[3][3], H4[3][3], H5[3][3]; 

char type_of_structure[12];
int CRYSTAL_FCC, LIQUID;
char potential_function[10]; 
int LJ6_12, WOON;
/* potential cutoff in reduced length unit */
#define CUTOFF (2.5*3.39/ULENGTH_IN_ANGSTROM)
#define RLIST (1.15*CUTOFF) /* in reduced length unit */
#define MAX_NEIGHBOR  256
int *idx, **list;  /* neighbour list */
#define NRDFD 1000 /* g(r) mesh */
int igr[NRDFD];
double r2del;

/** particle variables **/
double *f, *fs;   /* forces in real and rescaled space */
double *s;        /* positions in rescaled space [0,1) */
double *v, *s1;   /* velocities in real and rescaled space */
double *s2, *s3, *s4, *s5;   /* higher order derivatives */
double *mass, totalmass; /* particle masses in reduced unit */
double *x, *xold;    /* particle positions in real space */
double *d, *dold; /* particle displacements in real space */
double *ep; /* individual particle energies */

/** property cumulants **/
double pote, kine, hamiltonian;  /* potential and kinetic energy */
double pout[3][3], virial[3][3], stress[3][3]; /* external and instantaneous system stress */
double cj[3];  /* heat current cumulant */
double c[6][6]; /* elastic contants in Voigt notation */
double msd; /* mean squared displacement */

/** control variables **/
long maxkb; /* total number of MD steps */
long kstart; /* step at which to do property averaging */
long iseed; /* random number generator seed */
double scaler;  /* kT in reduced unit */
double delta;   /* integration time constant */
double tauT, tauH; /* particle and wall temperature coupling constants */
double wallmass; /* mass of the wall in reduced unit */
/* from here on are Boolean variables --> */
int calc_conductivity;
int NEH, NTP; 
int read_from_config, write_lp; 

/** thermal conductivity calculation parameters **/
#define NPART   256

/** input/output files **/
char fn_config_read[100], fn_config_write[100];
/* filenames of input and output configurations */
#define fn_property_output "temp.out"
/* common properties output file */
FILE *property_output;
#define fn_structure_output "struc.out"
/* structural properties output file */
FILE *structure_output;
#define fn_gr_output "gr.out"
/* pair correlation function output file */
FILE *gr_output;

/* from here on concerns thermal conductivity calculations --> */
static char *fn_heat_current[3]={"jx.out","jy.out","jz.out"};
/* heat current files */
FILE *heat_current[3];
#define fn_freq_output "freq.out"
/* heat current power spectrum output file */
FILE *freq_output;
#define fn_corr_output "corr.out"
/* heat current correlation function output file */
FILE *corr_output;

/* from here on concerns shear viscosity calculations --> */
static char *fn_shear_stress[3]={"syz.out","sxz.out","sxy.out"};
/* shear stress files */
FILE *shear_stress[3];
#define fn_shear_freq_output "shear_freq.out"
/* shear stress power spectrum output file */
FILE *shear_freq_output;
#define fn_shear_corr_output "shear_corr.out"
/* shear stress correlation function output file */
FILE *shear_corr_output;

/***************************/
/** function declarations **/
/***************************/
double frand();
void randnorm(double *x, int n);
void matran (double A[3][3], double B[3][3]);
void matadd (double A[3][3], double B[3][3], double C[3][3]);
void matsub (double A[3][3], double B[3][3], double C[3][3]);
double matinv (double A[3][3], double B[3][3]);
void matmul (double A[3][3], double B[3][3], double C[3][3]);
void assign_fcc_lattice();
void pair_potential();
void shear_smile(FILE *output);
void smile(FILE *output);
void realft(double data[], unsigned long n, int isign);
void four1(double data[], unsigned long nn, int isign);
double T_real_to_T_md (double TR, double *dTdR);
double T_md_to_T_real (double T);
void write_config (char filename[]);
void read_config (char filename[]);

#endif