/*****************************/
/* Lattice Dynamics of Argon */
/*                           */
/*      ver 1.2 Mar.22, 1999 */
/*  developed by Ju Li (MIT) */
/*****************************/

/*
  Version Notes:
  1. TTY interface.
  2. Woon's pair potential.
  3. Upon request, save the real
     space dynamical matrices of 
     atom(s) on a .dm file so if
     it is a perfect crystal, it
     can be used by dm.m to plot
     out the phonon dispersion
     curves.
  4. Calculate and diagonalize
     the dynamical matrix for
     arbitrary configuration and
     save the eigenmodes on a 
     .eig file, which can be
     handled by collect.m.
  5. User-defined supercell k
     sampling or random sampling.
  6. Upon request, evaluate DOS
     and LDOS(s).   
*/

#include <stdio.h>
#include <math.h>
#include <string.h>

/** universal constants **/
#ifndef PI
#define PI 3.14159265358979323846
#endif
#define AVO 6.0221367E+23
#define BOLZ 1.380658E-23  /* in J/K */
#define HBAR 1.054589E-34  /* in J.s */
#define EV_IN_J 1.60217733E-19
#define HARTREE_IN_J 4.3597482E-18
#define BOHR_RADIUS_IN_A 0.529177249
/** argon LJ6-12 constants **/
#define SIGMA_IN_M (3.405E-10)
#define EPSILON_IN_J (119.8*BOLZ) 
#define ARGON_MASS_IN_KG (39.948*1E-3/AVO) 
/** fundamental reduced units **/
#define ULENGTH_IN_M SIGMA_IN_M
#define UENERGY_IN_J EPSILON_IN_J
#define UMASS_IN_KG  ARGON_MASS_IN_KG     
/** derived reduced units **/
#define ULENGTH_IN_A (ULENGTH_IN_M*1E10)
#define UTIME_IN_SECOND (ULENGTH_IN_M*sqrt(UMASS_IN_KG/UENERGY_IN_J))
#define UTIME_IN_PS (UTIME_IN_SECOND*1E12)
#define UFREQ_IN_HZ (1./UTIME_IN_SECOND/2./PI)
#define UFREQ_IN_THZ (1E-12*UFREQ_IN_HZ)   
#define UFREQ_IN_MEV (HBAR*2*PI*UFREQ_IN_HZ/EV_IN_J*1000.)
/* interaction cutoff */    
#define RCUT (2.5*SIGMA_IN_M/ULENGTH_IN_M)
/* maximum possible reduced density of the solid */
#define MAXDENSITY 2.0
#define MAXNEIGHBOR ((int)(4./3.*PI*RCUT*RCUT*RCUT*MAXDENSITY))
#define MAX_STRLEN 150
#define MAX_LDOS_ATOM_BRACKETS 10
#define DEFAULT_ATOM_TP "Ar"
#define TP_STRLEN 2

static int np, LJ6_12, WOON, random_k_sampling, save_dm, 
    save_eig, save_dos, save_ldos, num_dm, num_kpts, iseed,
    num_bins, num_ldos_atom_brackets, num_ldos_atoms, *dm_index,
    nc[3], ldos_atom_head[MAX_LDOS_ATOM_BRACKETS],
    ldos_atom_tail[MAX_LDOS_ATOM_BRACKETS];
static double volume, minfreq, maxfreq, delfreq, *s, *mass, *dm,
    *dos, *ldos, *kpts, *ap, *w, *z, *work, *rwork,
    delta_s[3*MAXNEIGHBOR], H0[3][3], M[3][3], H[3][3], HI[3][3];
static char config_name[MAX_STRLEN], potential_function[MAX_STRLEN],
    ldos_atom_list[MAX_STRLEN], fn_config_read[MAX_STRLEN],
    fn_dm[MAX_STRLEN], fn_kpt[MAX_STRLEN], fn_eig[MAX_STRLEN],
    fn_dos[MAX_STRLEN], fn_ldos[MAX_STRLEN], **tp;
static FILE *filehandle;

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

int main()
{
    int i, j, k, l, m, ip;
    double cc, dd, ee;
    char *p;
    
    /* read instructions from input, such as "don" file */
    printf ("Reading in instructions...\n");
    
    fscanf (stdin, "Name of configuration (->*, "
	    "default=Data/fcc_unit_cell):\n%s\n\n",
	    config_name);
    if (!strcasecmp(config_name,"default"))
	sprintf(config_name, "Data/fcc_unit_cell");
    
    fscanf (stdin,
            "Pair potential to be used (LJ6_12 or WOON):\n%s\n\n",
            potential_function);
    LJ6_12 = !strcmp(potential_function, "LJ6_12");
    WOON = !strcmp(potential_function, "WOON");
    if ( (!LJ6_12) && (!WOON) )
    {
        fprintf (stdout, "Error: choose between LJ6_12 and WOON "
		 "pair potential\n");
        exit(1);
    }

    fscanf (stdin, "Load configuration from file (default=*):\n%s\n\n",
	    fn_config_read);
    if (!strcasecmp(fn_config_read,"default"))
	sprintf(fn_config_read, "%s", config_name);
    
    fscanf (stdin, "Apply affine transformation "
	    "(11,22,33,21=12,31=13,32=23):\n");
    fscanf (stdin, "%lf %lf %lf %lf %lf %lf\n\n",
	    &M[0][0], &M[1][1], &M[2][2],
	    &M[0][1], &M[0][2], &M[1][2]);
    /* symmetrize */
    M[1][0] = M[0][1];
    M[2][0] = M[0][2];
    M[2][1] = M[1][2];
    
    fscanf (stdin, "Save real-space dynamical matrices on file\n");
    fscanf (stdin, "(NULL=don't, default=*.dm):\n%s\n\n", fn_dm);
    if (save_dm=strcasecmp(fn_dm, "NULL"))
	if (!strcasecmp(fn_dm, "default"))
	    sprintf(fn_dm, "%s.dm", config_name);
    
    fscanf (stdin, "Supercell k-point sampling file (if numeric then\n");
    fscanf (stdin,
	    "apply random sampling in BZ, default=Kpts/fcckpt.10):\n%s\n\n",
	    fn_kpt);
    if (random_k_sampling=((*fn_kpt>='0')&&(*fn_kpt<='9')))
    {	
	/* user inputs a integer, do random k sampling */
	for (p=fn_kpt;(*p>='0')&&(*p<='9');p++); *p = 0;
	sscanf(fn_kpt, "%d", &num_kpts);
    }
    else if (!strcasecmp(fn_kpt, "default"))
	sprintf(fn_kpt, "Kpts/fcckpt.10");
	  
    fscanf (stdin,
	    "Save normal modes on file (NULL=don't, default=*.eig):\n%s\n\n",
	    fn_eig);
    if (save_eig=strcasecmp(fn_eig, "NULL"))
	if (!strcasecmp(fn_eig, "default"))
	    sprintf(fn_eig, "%s.eig", config_name);
    
    fscanf (stdin,
	    "List atoms (NULL=don't, default=all atoms):\n%s\n\n",
	    ldos_atom_list);

    fscanf (stdin,
	    "Save DOS on file (NULL=don't, default=*.out):\n%s\n\n",
	    fn_dos);
    if (save_dos=strcasecmp(fn_dos, "NULL"))
	if (!strcasecmp(fn_dos, "default"))
	    sprintf(fn_dos, "%s.out", config_name);
    
    fscanf (stdin,
	    "Save LDOS on file (NULL=don't, default=*_ldos.out):\n%s\n\n",
	    fn_ldos);
    if (save_ldos=strcasecmp(fn_ldos, "NULL"))
	if (!strcasecmp(fn_ldos, "default"))
	    sprintf(fn_ldos, "%s_ldos.out", config_name);
    
    fscanf (stdin, "Minimum frequency (in THz):\n%lf\n\n", &minfreq);
    fscanf (stdin, "Maximum frequency (in THz):\n%lf\n\n", &maxfreq);
    fscanf (stdin, "How many bins to represent DOS and LDOS:\n%d\n\n",
	    &num_bins);
    delfreq = (maxfreq-minfreq)/num_bins;
    
    fscanf (stdin, "Random number generator seed (0->time seed):\n");
    fscanf (stdin, "%d\n", &iseed); 
    if (iseed == 0) iseed = time(NULL);
    srandom ((unsigned)iseed);
        
    printf("Instructions read complete.\n\n");
    
    printf ("Vibrational eigenmodes will be calculated for \"%s\"\n",
	    config_name);
    printf ("using %s potential with Rcut = %f = %f Angstrom.\n\n",
	  potential_function, RCUT, RCUT*ULENGTH_IN_A);
    
    printf ("The configuration is being loaded from \"%s\".\n",
	    fn_config_read); 
    if ((filehandle=fopen(fn_config_read, "r"))==NULL)
    {
	fprintf(stdout, "Error: cannot open file \"%s\"\n.",
		fn_config_read);
	exit(1);
    }
    fscanf (filehandle, "Number of particles = %d\n\n", &np);
    printf ("Allocate simple arrays for %d particles.\n\n", np);
    s = (double *) malloc(3*np*sizeof(double));
    mass = (double *) malloc(np*sizeof(double));
    /* name space for atom type */
    tp   = (char **) malloc (np*sizeof(char *));
    tp[0] = (char *) malloc (np*(TP_STRLEN+1)*sizeof(char));
    for (i=1; i<np; i++)
	tp[i] = tp[i-1] + TP_STRLEN + 1;
    
    if (save_dos)
    {
	dos = (double *) malloc (num_bins*sizeof(double));
	for (i=0; i<num_bins; i++) dos[i] = 0.;
    }
    
    /* read in supercell geometry from configuration file */ 
    fscanf (filehandle, "H(1,1) = %lf A\n", &H0[0][0]);
    fscanf (filehandle, "H(1,2) = %lf A\n", &H0[0][1]);
    fscanf (filehandle, "H(1,3) = %lf A\n", &H0[0][2]);
    fscanf (filehandle, "nc(1)  = %d\n\n",  &nc[0]);
    fscanf (filehandle, "H(2,1) = %lf A\n", &H0[1][0]);
    fscanf (filehandle, "H(2,2) = %lf A\n", &H0[1][1]);
    fscanf (filehandle, "H(2,3) = %lf A\n", &H0[1][2]);
    fscanf (filehandle, "nc(2)  = %d\n\n",  &nc[1]);
    fscanf (filehandle, "H(3,1) = %lf A\n", &H0[2][0]);
    fscanf (filehandle, "H(3,2) = %lf A\n", &H0[2][1]);
    fscanf (filehandle, "H(3,3) = %lf A\n", &H0[2][2]);
    fscanf (filehandle, "nc(3)  = %d\n\n",  &nc[2]);
    /* convert to reduced units */
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	    H0[i][j] /= ULENGTH_IN_A;
    
    /* apply affine transformation */
    H[0][0] = H0[0][0]*M[0][0] + H0[0][1]*M[1][0] + H0[0][2]*M[2][0];
    H[0][1] = H0[0][0]*M[0][1] + H0[0][1]*M[1][1] + H0[0][2]*M[2][1];
    H[0][2] = H0[0][0]*M[0][2] + H0[0][1]*M[1][2] + H0[0][2]*M[2][2];
    H[1][0] = H0[1][0]*M[0][0] + H0[1][1]*M[1][0] + H0[1][2]*M[2][0];
    H[1][1] = H0[1][0]*M[0][1] + H0[1][1]*M[1][1] + H0[1][2]*M[2][1];
    H[1][2] = H0[1][0]*M[0][2] + H0[1][1]*M[1][2] + H0[1][2]*M[2][2];
    H[2][0] = H0[2][0]*M[0][0] + H0[2][1]*M[1][0] + H0[2][2]*M[2][0];
    H[2][1] = H0[2][0]*M[0][1] + H0[2][1]*M[1][1] + H0[2][2]*M[2][1];
    H[2][2] = H0[2][0]*M[0][2] + H0[2][1]*M[1][2] + H0[2][2]*M[2][2];
    volume = fabs(matinv(H, HI));

    printf ("After transformation the supercell geometry is\n");
    printf ("        | %.5f %.5f %.5f |\n",
	  H[0][0]*ULENGTH_IN_A,
	  H[0][1]*ULENGTH_IN_A,
	  H[0][2]*ULENGTH_IN_A);	  
    printf (" H[A] = | %.5f %.5f %.5f |\n",
	  H[1][0]*ULENGTH_IN_A,
	  H[1][1]*ULENGTH_IN_A,
	  H[1][2]*ULENGTH_IN_A);
    printf ("        | %.5f %.5f %.5f |\n",
	  H[2][0]*ULENGTH_IN_A,
	  H[2][1]*ULENGTH_IN_A,
	  H[2][2]*ULENGTH_IN_A);
    
    printf ("Supercell volume = %f = %f A^3.\n\n", volume,
	    volume*ULENGTH_IN_A*ULENGTH_IN_A*ULENGTH_IN_A);
    if (np/volume > 0.85*MAXDENSITY)
    {
	printf("Current density %f is above 85%% of MAXDENSITY = %f.\n",
	       np/volume, MAXDENSITY);
	printf("Overflow possible. Please increase MAXDENSITY.\n");
	exit(1);
    }
    
    /* read in reduced coordinates from configuration file */ 
    fscanf(filehandle, "TP Mass(amu)     sx        sy        sz"
	   "     sx1(1/ns)  sy1(1/ns)  sz1(1/ns)\n");
    for (i=0; i<3*np; i+=3)
    {
	fscanf (filehandle, "%2s %lf %lf %lf %lf %lf %lf %lf\n",
		tp[i/3], &mass[i/3], &s[i], &s[i+1], &s[i+2],
		&cc,  &cc, &cc);
	/* convert from AMU to reduced mass unit */
	mass[i/3] = mass[i/3]/AVO*1E-3/UMASS_IN_KG;
    }
    fclose (filehandle);
    
    printf ("Configuration \"%s\" generated successfully on the computer.\n\n",
	    config_name);
    
    /* count the total number of interactions */
    printf ("Counting all possible (image) interactions within Rcut = %fA.\n",
	    RCUT*ULENGTH_IN_A);
    for (num_dm=0,k=0,i=0; i<np; i++)
	for (j=i; j<np; j++)
	{
	    l = images_within_radius (H, HI, s, j, i, delta_s);
	    num_dm += l;
	    if (j==i) k+=l;
	}
    
    printf("There are %d (image) interactions in the entire system,\n",
	   num_dm);
    printf("and on the average %d neighbors per atom.\n\n",
	   (2*num_dm-k)/np);
    
    printf("Allocating dynamical matrices for all these interactions.\n");
    dm_index = (int *)malloc(2*num_dm*sizeof(int));
    /* save the particle index j and i*/
    dm = (double *)malloc((3+6)*num_dm*sizeof(double));
    /* the first 3 are delta_s between j and i; the next 6 are (xx,
       xy, xz, yy, yz, zz) dynamical matrix components, since in pair
       potentials they are symmetric */
    if ((dm==NULL)||(dm_index==NULL))
    {
	printf("Not enough memory (%d bytes).\n",
	       (3+6)*num_dm*sizeof(double)+2*num_dm*sizeof(int));
	exit(1);
    }
    
    /* evaluating the dynamical matrices */
    for (ip=0,i=0; i<np; i++)
	for (j=i; j<np; j++)
	{
	    k = images_within_radius (H, HI, s, j, i, delta_s);
	    for (l=0; l<k; l++)
	    {
		/* save the particle indices */
		dm_index[2*ip] = j;
		dm_index[2*ip+1] = i;
		/* save the relative delta_s */
		dm[9*ip]   = delta_s[3*l];
		dm[9*ip+1] = delta_s[3*l+1];
		dm[9*ip+2] = delta_s[3*l+2];
		/* calculate the force constant matrix */
		force_constants (H, delta_s+3*l, dm+9*ip+3);
		dm[9*ip+3] *= UFREQ_IN_THZ*UFREQ_IN_THZ/sqrt(mass[i]*mass[j]);
		dm[9*ip+4] *= UFREQ_IN_THZ*UFREQ_IN_THZ/sqrt(mass[i]*mass[j]);
		dm[9*ip+5] *= UFREQ_IN_THZ*UFREQ_IN_THZ/sqrt(mass[i]*mass[j]);
		dm[9*ip+6] *= UFREQ_IN_THZ*UFREQ_IN_THZ/sqrt(mass[i]*mass[j]);
		dm[9*ip+7] *= UFREQ_IN_THZ*UFREQ_IN_THZ/sqrt(mass[i]*mass[j]);
		dm[9*ip+8] *= UFREQ_IN_THZ*UFREQ_IN_THZ/sqrt(mass[i]*mass[j]);
		ip++;
	    }
	}

    /* save the real-space dynamical matrices on .dm file */
    if (save_dm)
    {
	printf ("Saving real space dynamical matrices (A and THz^2)\n");
	printf ("on \"%s\".\n\n", fn_dm);
	filehandle = fopen(fn_dm, "w+");
	/* number of records */
	fprintf (filehandle, "%d\n", num_dm);
	for (ip=0; ip<num_dm; ip++)
	{
	    fprintf (filehandle, "%d\n", dm_index[2*ip]); /* j */
	    fprintf (filehandle, "%d\n", dm_index[2*ip+1]); /* i */
	    /* ji delta_s */
	    fprintf (filehandle, "%f\n", dm[9*ip]); 
	    fprintf (filehandle, "%f\n", dm[9*ip+1]);
	    fprintf (filehandle, "%f\n", dm[9*ip+2]);
	    /* dynamical matrix elements in THz^2 */
	    fprintf (filehandle, "%f\n", dm[9*ip+3]); 
	    fprintf (filehandle, "%f\n", dm[9*ip+4]);
	    fprintf (filehandle, "%f\n", dm[9*ip+5]);
	    fprintf (filehandle, "%f\n", dm[9*ip+6]);
	    fprintf (filehandle, "%f\n", dm[9*ip+7]);
	    fprintf (filehandle, "%f\n", dm[9*ip+8]);
	}
	fclose(filehandle);
    }

    /* recognize which atomic components should be recorded */
    printf ("Parsing input of LDOS indices \"%s\":\n", ldos_atom_list);
    if (!strcasecmp(ldos_atom_list, "NULL"))
	num_ldos_atom_brackets = 0;
    else if (!strcasecmp(ldos_atom_list, "default"))
    {
	/* all atoms inside supercell */
	num_ldos_atom_brackets = 1;
	ldos_atom_head[0] = 0;
	ldos_atom_tail[0] = np-1;
    }
    else
    {
	/* no numerical bracket is found */
	if (!(num_ldos_atom_brackets=parse_number_string
	      (ldos_atom_list, ldos_atom_head, ldos_atom_tail)))
	{
	    printf ("syntax error, no recognizable integer brackets.\n");
	    exit(1);
	}
	else if (num_ldos_atom_brackets>MAX_LDOS_ATOM_BRACKETS)
	{
	    printf ("number of brackets more than MAX_LDOS_ATOM_BRACKETS.\n");
	    exit(1);
	}
    }

    /* find out how many atoms we need to keep track of */
    printf ("Number of brackets recognized = %d. They are\n",
	    num_ldos_atom_brackets);
    for (num_ldos_atoms=0,i=0; i<num_ldos_atom_brackets; i++)
    {
	if (ldos_atom_tail[i]<ldos_atom_head[i])
	{
	    k = ldos_atom_tail[i];
	    ldos_atom_tail[i] = ldos_atom_head[i];
	    ldos_atom_head[i] = k;
	}
	printf ("atoms %d to %d \n", ldos_atom_head[i], ldos_atom_tail[i]);
	num_ldos_atoms += ldos_atom_tail[i]-ldos_atom_head[i]+1;
    }
    if (save_ldos)
    {
	ldos = (double *)
	    malloc (3*num_ldos_atoms*num_bins*sizeof(double));
	for (i=0; i<3*num_ldos_atoms*num_bins; i++) ldos[i] = 0.;
    }	
    printf ("\n");
    
    /* get sampling k-points in BZ */
    if (random_k_sampling)
    {
	printf("Randomly generating %d k-points in supercell BZ.\n\n",
	       num_kpts);
	kpts = (double *) malloc(4*num_kpts*sizeof(double));
	for (i=0; i<num_kpts; i++)
	{
	    /* coefficients of reciprocal lattice vectors */ 
	    kpts[4*i]   = frandom();
	    kpts[4*i+1] = frandom();
	    kpts[4*i+2] = frandom();
	    /* sampling weight */
	    kpts[4*i+3] = 1./num_kpts;
	}
    }
    else
    {
	printf("Reading supercell sampling k-points from \"%s\":\n",
	       fn_kpt);
	if ((filehandle=fopen(fn_kpt, "r"))==NULL)
	{
	    printf ("sampling k-point file \"%s\" does not exist.\n",
		    fn_kpt);
	    exit(1);
	}
	fscanf(filehandle, "%d\n", &num_kpts);
	kpts = (double *) malloc(4*num_kpts*sizeof(double));
	printf("there are %d k-points stored in \"%s\"\n",
	       num_kpts, fn_kpt);
	printf("(in reciprocal lattice vector coefficients).\n\n");
	for (i=0; i<num_kpts; i++)
	    fscanf (filehandle, "%lf %lf %lf %lf\n",
		    kpts+4*i, kpts+4*i+1, kpts+4*i+2, kpts+4*i+3);
	fclose(filehandle);
    }	
    
    printf ("Allocate memory for k-space dynamical matrix "
	    "and diagonalization.\n");
    ap = (double *) malloc(3*np*(3*np+1)*sizeof(double));
    w = (double *) malloc(3*np*sizeof(double));
    z = (double *) malloc(2*3*np*3*np*sizeof(double));
    work = (double *) malloc(4*3*np*sizeof(double));
    rwork = (double *) malloc(6*3*np*sizeof(double));

    /* save overall information on .eig file */
    if (save_eig)
    {
	printf ("Eigenmodes will be saved on \"%s\".\n\n", fn_eig);
	filehandle = fopen(fn_eig, "w+");
	fprintf (filehandle, "Number of atoms = %d\n\n", np);
	fprintf (filehandle, "Number of k-points = %d\n\n", num_kpts);
	fprintf (filehandle,
		 "Vibrational normal modes will be recorded for %d atoms:\n",
		 num_ldos_atoms);
	for (i=0; i<num_ldos_atom_brackets; i++)
	    for (j=ldos_atom_head[i]; j<=ldos_atom_tail[i]; j++) 
		fprintf (filehandle, "atom %d\n", j);
	fprintf (filehandle, "\n");
    }
    
    printf ("Calculation starts (diagonalize D(k) with "
	    "LAPACK routine zhpev_):\n\n");
    
    /* go over all supercell k-points */
    for (k=0; k<num_kpts; k++)
    {
	/* save k-point information on .eig file */
	if (save_eig)
	{	
	    fprintf (filehandle, "k-point number %d, weight = %f\n",
		     k, kpts[4*k+3]);
	    fprintf (filehandle, "(%f %f %f) in reciprocal lattice vectors\n",
		     kpts[4*k], kpts[4*k+1], kpts[4*k+2]);
	    fprintf (filehandle, "freq(THz)  Re(u_x)  Im(u_x)  Re(u_y)  "
		     "Im(u_y)  Re(u_z)  Im(u_z) ->\n");
	}
	if (k%(int)ceil(10000./np/np/np)==0)
	    printf ("Doing k-point %d.\n", k);
	/* assemble the k-space dynamical matrix */
	for (i=0; i<3*np*(3*np+1); i++) ap[i]=0.;
	for (ip=0; ip<num_dm; ip++)
	{
	    /* j is always greater than or equal to i */
	    j = dm_index[2*ip];
	    i = dm_index[2*ip+1];
	    
	    /* (ix, ix) */
	    m = 9*i*i+9*i;
	    ap[m] += dm[9*ip+3];
	    /* (ix, iy) */
	    ap[m+6*i+2] += dm[9*ip+4];
	    /* (ix, iz) */
	    ap[m+12*i+6] += dm[9*ip+5];
	    /* (iy, iy) */
	    ap[m+6*i+4] += dm[9*ip+6];
	    /* (iy, iz) */
	    ap[m+12*i+8] += dm[9*ip+7];
	    /* (iz, iz) */
	    ap[m+12*i+10] += dm[9*ip+8];

	    /* phase */
	    ee = 2*PI*(kpts[4*k]*dm[9*ip]+
		       kpts[4*k+1]*dm[9*ip+1]+
		       kpts[4*k+2]*dm[9*ip+2]);
	    cc = cos(ee);
	    dd = sin(ee);
	    
	    /* (ix, jx) */
	    m = 9*j*j+3*j+6*i;
	    ap[m] -= dm[9*ip+3]*cc;
	    ap[m+1] -= dm[9*ip+3]*dd;
	    /* (ix, jy) */
	    ap[m+6*j+2] -= dm[9*ip+4]*cc;
	    ap[m+6*j+3] -= dm[9*ip+4]*dd;
	    /* (ix, jz) */
	    ap[m+12*j+6] -= dm[9*ip+5]*cc;
	    ap[m+12*j+7] -= dm[9*ip+5]*dd;
	    /* (iy, jy) */
	    ap[m+6*j+4] -= dm[9*ip+6]*cc;
	    ap[m+6*j+5] -= dm[9*ip+6]*dd;
	    /* (iy, jz) */
	    ap[m+12*j+8] -= dm[9*ip+7]*cc;
	    ap[m+12*j+9] -= dm[9*ip+7]*dd;
	    /* (iz, jz) */
	    ap[m+12*j+10] -= dm[9*ip+8]*cc;	    
	    ap[m+12*j+11] -= dm[9*ip+8]*dd;	    
	    
	    if (j!=i)
	    {
		/* (alpha, beta) is symmetric for pair interaction */
		/* (iy, jx) */
		ap[m+2] -= dm[9*ip+4]*cc;
		ap[m+3] -= dm[9*ip+4]*dd;	    
		/* (iz, jx) */
		ap[m+4] -= dm[9*ip+5]*cc;
		ap[m+5] -= dm[9*ip+5]*dd;	    
		/* (iz, jy) */
		ap[m+6*j+6] -= dm[9*ip+7]*cc;
		ap[m+6*j+7] -= dm[9*ip+7]*dd;
		
		/* do the same for j */
		m = 9*j*j+9*j;
		ap[m] += dm[9*ip+3];
		/* (jx, jy) */
		ap[m+6*j+2] += dm[9*ip+4];
		/* (jx, jz) */
		ap[m+12*j+6] += dm[9*ip+5];
		/* (jy, jy) */
		ap[m+6*j+4] += dm[9*ip+6];
		/* (jy, jz) */
		ap[m+12*j+8] += dm[9*ip+7];
		/* (jz, jz) */
		ap[m+12*j+10] += dm[9*ip+8];
	    }
	}	    

	ip = 3*np;
	i = 0;
	zhpev_ ("V", "U", &ip, ap, w, z, &ip, work, rwork, &i);
	if (save_eig) printf ("Save eigenmodes in \"%s\".\n\n", fn_eig);

	/* the eigenvalues are already assorted in ascending order */
	for (l=0; l<3*np; l++)
	{
	    /* frequency in THz, if imaginary than save as negative */
	    w[l] = (w[l]>=0.)?sqrt(w[l]):-sqrt(-w[l]);
	    if (save_eig) fprintf (filehandle, "%f  ", w[l]);
	    /* frequency bin index */
	    ip = floor((w[l]-minfreq)/delfreq);
	    if ( (ip<0) || (ip>=num_bins) )
	    {
		printf ("mode w=%f THz out of freq. bound\n",
			w[l]);
		ip = -1;
	    }
	    if (save_dos && (ip!=-1)) dos[ip] += kpts[4*k+3];
	    for (m=0,i=0; i<num_ldos_atom_brackets; i++)
		for (j=ldos_atom_head[i]; j<=ldos_atom_tail[i]; j++)
		{
		    if (save_eig)
		    {
			fprintf (filehandle,
				 "%8.5f %8.5f %8.5f %8.5f %8.5f %8.5f ",
				 z[2*(l*3*np+3*j)],   z[2*(l*3*np+3*j)+1],
				 z[2*(l*3*np+3*j)+2], z[2*(l*3*np+3*j)+3],
				 z[2*(l*3*np+3*j)+4], z[2*(l*3*np+3*j)+5]);
			fprintf (filehandle, "\n");
		    }
		    if (save_ldos && (ip!=-1))
		    {
			ldos[m*3*num_bins+ip] +=
			    (z[2*(l*3*np+3*j)]*z[2*(l*3*np+3*j)]+
			     z[2*(l*3*np+3*j)+1]*z[2*(l*3*np+3*j)+1])*
			    kpts[4*k+3];
			ldos[(m*3+1)*num_bins+ip] +=
			    (z[2*(l*3*np+3*j)+2]*z[2*(l*3*np+3*j)+2]+
			    z[2*(l*3*np+3*j)+3]*z[2*(l*3*np+3*j)+3])*
			    kpts[4*k+3];
			ldos[(m*3+2)*num_bins+ip] +=
			    (z[2*(l*3*np+3*j)+4]*z[2*(l*3*np+3*j)+4]+
			    z[2*(l*3*np+3*j)+5]*z[2*(l*3*np+3*j)+5])*
			    kpts[4*k+3];
		    }
		    m++;
		}
	}
	if (save_eig) fprintf (filehandle, "\n");
    }
    if (save_eig) fclose(filehandle);
    printf ("\n");
    
    if (save_dos)
    {
	printf ("Saving DOS on \"%s\":\n", fn_dos);
	printf ("total of %d eigenvalues, on the average %d per bin.\n\n",
		3*np*num_kpts, 3*np*num_kpts/num_bins);
	filehandle = fopen(fn_dos, "w+");
	for (i=0; i<num_bins; i++)
	    fprintf (filehandle, "%9.6f %9.6f\n",
		     minfreq+(i+0.5)*delfreq, dos[i]/3./np);
	fclose(filehandle);
    }
    
    if (save_ldos)
    {
	printf ("Saving LDOS on \"%s\":\n", fn_ldos);
	printf ("total of %d eigenvalues, on the average %d per bin.\n\n",
		3*np*num_kpts, 3*np*num_kpts/num_bins);
	filehandle = fopen(fn_ldos, "w+");
	fprintf (filehandle, "%9d ", num_ldos_atoms);
	for (i=0; i<num_ldos_atom_brackets; i++)
	    for (j=ldos_atom_head[i]; j<=ldos_atom_tail[i]; j++)
		for (l=0; l<3; l++)
		    fprintf(filehandle, "%9.6f ", j+(double)l/10.);
	fprintf (filehandle, "\n");
	for (k=0; k<num_bins; k++)
	{
	    fprintf (filehandle, "%9.6f ", minfreq+(k+0.5)*delfreq);
	    for (m=0; m<3*num_ldos_atoms; m++)
		fprintf(filehandle, "%9.6f ", ldos[m*num_bins+k]);
	    fprintf (filehandle, "\n");
	}
	fclose(filehandle);
    }
    
    /* free all allocated memories */
    printf ("Freeing all allocated memories.\n");
    free (rwork);
    free (work);
    free (w);
    free (z);
    free (ap);
    free (kpts);
    free (dm_index);
    free (dm);
    free (mass);
    free (s);
    if (save_dos) free(dos);
    if (save_ldos) free(ldos);
    
    printf("Mission terminates successfully.\n\n");
    return 0;
} /* end main() */

/* Woon's pair potential, Chem. Phys. Lett. 204, 29 (1993) */
#define WOON_Q 0.6060
#define WOON_EPSILON (0.4228*HARTREE_IN_J/1000./UENERGY_IN_J)
#define WOON_RE (7.1714*BOHR_RADIUS_IN_A/ULENGTH_IN_A)
#define WOON_KE (0.6627*HARTREE_IN_J/1000./UENERGY_IN_J/ \
		 pow(BOHR_RADIUS_IN_A/ULENGTH_IN_A, 2.))
#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)
#define WOON_RESIDUE (WOON_EPSILON*(WOON_Q*exp(-WOON_K1*(RCUT-WOON_RE))- \
				    (WOON_Q+1)*exp(-WOON_K2*(RCUT-WOON_RE))))
#define WOON_RESIDUE_FORCE \
    (WOON_EPSILON*(WOON_Q*WOON_K1*exp(-WOON_K1*(RCUT-WOON_RE))- \
		   (WOON_Q+1)*WOON_K2*exp(-WOON_K2*(RCUT-WOON_RE))))
#define LJ_RESIDUE (4*(pow(RCUT,-12.)-pow(RCUT,-6.)))
#define LJ_RESIDUE_FORCE ((48*pow(RCUT,-12.)-24*pow(RCUT,-6.))/RCUT)
/* Evaluate the force constants (xx, xy, xz, yy, yz, zz) */
/* of a pair interaction. */
void force_constants
(double H[3][3], double *s, double *force_constant)
{
    double r1, r2, a, b; 
    double ir2, ir4, ir6, ir12;
    double first, second; /* V'/r and V''/r */
    double rx, ry, rz;
    double res = LJ6_12?LJ_RESIDUE:WOON_RESIDUE;
    double ref = LJ6_12?LJ_RESIDUE_FORCE:WOON_RESIDUE_FORCE;
    
    rx = H[0][0]*s[0] + H[1][0]*s[1] + H[2][0]*s[2];
    ry = H[0][1]*s[0] + H[1][1]*s[1] + H[2][1]*s[2];
    rz = H[0][2]*s[0] + H[1][2]*s[1] + H[2][2]*s[2];
    r2 = rx*rx+ry*ry+rz*rz;
    r1 = sqrt(r2);
    
    if (LJ6_12)
    {
	ir2 = 1./r2;
	ir4 = ir2 * ir2;
	ir6 = ir4 * ir2;
	ir12 = ir6 * ir6;
	first = (-48.*ir12+24.*ir6)*ir2 + ref/r1;
	second = (624.*ir12-168.*ir6)*ir4;
    }
    else if(WOON)
    {
	a = WOON_EPSILON * WOON_Q * exp(-WOON_K1 * (r1-WOON_RE));
	b = WOON_EPSILON * (WOON_Q+1) * exp(-WOON_K2 * (r1-WOON_RE)); 
	first = (-WOON_K1*a+WOON_K2*b+ ref)/r1;
	second = (WOON_K1*WOON_K1*a-WOON_K2*WOON_K2*b)/r2;
    }
    force_constant[0] = (second-first/r2)*rx*rx+first;
    force_constant[1] = (second-first/r2)*rx*ry;
    force_constant[2] = (second-first/r2)*rx*rz;
    force_constant[3] = (second-first/r2)*ry*ry+first;
    force_constant[4] = (second-first/r2)*ry*rz;
    force_constant[5] = (second-first/r2)*rz*rz+first;
    return;
} /* end force_constants() */


/* calculate all images of particle j to particle i */
/* (excluding self-interaction when j = i)          */
int images_within_radius
(double H[3][3], double HI[3][3], double *s,
 int j, int i, double delta_s[3*MAXNEIGHBOR])
{
    double sxji, syji, szji, rx, ry, rz;
    int nx, ny, nz, ix, iy, iz, neighbor;
    
    /* the RCUT sphere has to be contained in the parallelepiped
       formed by (nx*H(1,:), ny*H(2,:), nz*H(3,:)), which requires the
       thickness opposite to any of the parallelepiped surfaces to be
       greater than RCUT, but that thickness is just 1/|G_i|, the
       corresponding reciprocal vector, if G^T * H = I.  */
    nx = ceil(RCUT*sqrt(HI[0][0]*HI[0][0]+
			  HI[1][0]*HI[1][0]+HI[2][0]*HI[2][0]));
    ny = ceil(RCUT*sqrt(HI[0][1]*HI[0][1]+
			  HI[1][1]*HI[1][1]+HI[2][1]*HI[2][1]));
    nz = ceil(RCUT*sqrt(HI[0][2]*HI[0][2]+
			  HI[1][2]*HI[1][2]+HI[2][2]*HI[2][2]));
    neighbor = 0;
    for (ix=-nx; ix<=nx; ix++)
	for (iy=-ny; iy<=ny; iy++)
	    for (iz=-nz; iz<=nz; iz++)
		if (!((j==i)&&(ix==0)&&(iy==0)&&(iz==0)))
		{
		    sxji = s[3*j]   - s[3*i]   + ix;
		    syji = s[3*j+1] - s[3*i+1] + iy;
		    szji = s[3*j+2] - s[3*i+2] + iz;
		    rx = H[0][0]*sxji + H[1][0]*syji + H[2][0]*szji;
		    ry = H[0][1]*sxji + H[1][1]*syji + H[2][1]*szji;
		    rz = H[0][2]*sxji + H[1][2]*syji + H[2][2]*szji;
		    if (rx*rx+ry*ry+rz*rz<RCUT*RCUT)
		    {
			delta_s[3*neighbor]   = sxji;
			delta_s[3*neighbor+1] = syji;
			delta_s[3*neighbor+2] = szji;
			neighbor++;
		    }
		}
    return (neighbor);
} /* end images_within_radius() */


int parse_number_string(char *sentence, int head[], int tail[])
{
    char buffer[MAX_STRLEN], *p, *q, *r;
    int number_of_brackets = 0;
    p = sentence;
    while(1)
    {
	if ((*p==0)||(p-sentence>MAX_STRLEN)) break;
	if ((*p>='0')&&(*p<='9'))
	{
	    /* read in a pair */
	    strcpy(buffer, p);
	    for (q=buffer;(*q>='0')&&(*q<='9');q++);*q=0;
	    /* read in head */
	    sscanf(buffer, "%d", head+number_of_brackets);
	    for (q++;(!((*q>='0')&&(*q<='9')))&&(*q!=0);q++);
	    /* if there is no tail */
	    if (*q==0)
	    {
		tail[number_of_brackets]=head[number_of_brackets];
		return (number_of_brackets+1);
	    }
	    else
	    {
		/* read in tail */
		for (r=q;(*q>='0')&&(*q<='9');q++);*q=0;
		sscanf(r,"%d",tail+number_of_brackets);
		number_of_brackets++;
		p+=q-buffer;
	    }
	}
	p++;
    }
    return number_of_brackets;
} /* end parse_number_string() */


double matinv(double A[3][3], double B[3][3])
{
    double D11,D22,D33,D12,D21,D13,D31,D32,D23,volume;
    D11=A[1][1]*A[2][2]-A[1][2]*A[2][1];
    D22=A[2][2]*A[0][0]-A[2][0]*A[0][2];
    D33=A[0][0]*A[1][1]-A[0][1]*A[1][0];
    D12=A[1][2]*A[2][0]-A[1][0]*A[2][2];
    D23=A[2][0]*A[0][1]-A[2][1]*A[0][0];
    D31=A[0][1]*A[1][2]-A[0][2]*A[1][1];
    D13=A[1][0]*A[2][1]-A[2][0]*A[1][1];
    D21=A[2][1]*A[0][2]-A[0][1]*A[2][2];
    D32=A[0][2]*A[1][0]-A[1][2]*A[0][0];
    volume=A[0][0]*D11+A[0][1]*D12+A[0][2]*D13;
    B[0][0]=D11/volume;
    B[1][1]=D22/volume;
    B[2][2]=D33/volume;
    B[0][1]=D21/volume;
    B[1][2]=D32/volume;
    B[2][0]=D13/volume;
    B[1][0]=D12/volume;
    B[2][1]=D23/volume;
    B[0][2]=D31/volume;
    return (volume);
} /* end matinv() */


/* returns pseudo-random number on (0,1) */
double frandom()
{
    return ((random()+0.5)/pow((double)2, (double)31));
} /* end frandom() */