/**********/
/* diag.h */
/**********/

#ifndef _DIAG_H
#define _DIAG_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 */

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

/** reduced units based on 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 UFREQ_IN_HZ (1./UTIME_IN_SECOND/2./PI)
#define UFREQ_IN_THZ (1e-12*UFREQ_IN_HZ)   
#define EV_IN_JOULE 1.60217733E-19
#define UFREQ_IN_MEV (HBAR*2*PI*UFREQ_IN_HZ/EV_IN_JOULE*1000.)

/* Woon's ab initio pair potential */
#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) */

/* potential cutoff in reduced length unit */
#define CUTOFF (2.5*3.39/ULENGTH_IN_ANGSTROM)
/* maximum possible reduced density of the solid */
#define MAXDENSITY 2.0
#define MAXNEIGHBOR ((int)(4./3.*PI*CUTOFF*CUTOFF*CUTOFF*MAXDENSITY))
#define MAX_STRING_LENGTH 150

/** common varaibles **/
int np; /* number of particles */
double H0[3][3], M[3][3], H[3][3], HI[3][3], volume;

/** particle variables **/
double *s;   /* positions in rescaled supercell [0,1) */
double *mass; /* particle masses in reduced unit */

/** control variables **/
int LJ6_12, WOON, random_k_sampling;
int save_dm, save_eig, save_dos, save_ldos;
char config_name[MAX_STRING_LENGTH];
char potential_function[MAX_STRING_LENGTH]; 
long iseed; /* random number generator seed */

/** real-space dynamical matrices **/
double *dm;
int num_dm, *dm_index;
double delta_s[3*MAXNEIGHBOR];
/** compressed 1D storage **/

/* sample k-points */
int num_kpts;
double *kpts;

/* dense diagonalization zhpev_ storages */
double *ap, *w, *z, *work, *rwork;

/** LDOS atoms **/
#define MAX_LDOS_ATOM_BRACKETS 10
char ldos_atom_list[MAX_STRING_LENGTH];
int ldos_atom_head[MAX_LDOS_ATOM_BRACKETS];
int ldos_atom_tail[MAX_LDOS_ATOM_BRACKETS];
int num_ldos_atom_brackets, num_ldos_atoms;
double minfreq, maxfreq, delfreq;
int num_bins;
double *dos, *ldos;

/** input/output files **/
char fn_config_read[MAX_STRING_LENGTH];
char fn_dm[MAX_STRING_LENGTH];
char fn_kpt[MAX_STRING_LENGTH];
char fn_eig[MAX_STRING_LENGTH];
char fn_dos[MAX_STRING_LENGTH];
char fn_ldos[MAX_STRING_LENGTH];
FILE *filehandle;

/***************************/
/** function declarations **/
/***************************/
void force_constants
(double H[3][3], double *s, double *force_constant);
int images_within_radius
(double H[3][3], double HI[3][3], double *s,
 int j, int i, double delta_s[3*MAXNEIGHBOR]);
int parse_number_string(char *sentence, int head[], int tail[]);
double matinv(double A[3][3], double B[3][3]);
double frandom();

#endif