/*******************************/
/*  Lattice Dynamics of Argon  */
/*                             */
/*      Ver 1.1 June.13, 1998  */
/*    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>
#include "diag.h"

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=fcc_unit_cell):\n%s\n\n",
	    config_name);
    if (!strcasecmp(config_name,"default"))
	sprintf(config_name, "fcc_unit_cell");
    
    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 (long integer):\n%ld\n\n",
	    &iseed);
    srandom (iseed);
    
    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, "\n **Error: Choose between LJ6_12 or WOON pair potential\n");
        exit(1);
    }
    printf("Instructions read complete.\n\n");
    /* Instructions read complete */
    
    printf ("Vibrational eigenmodes will be calculated for \"%s\"\n", config_name);
    printf ("using %s potential with Rcut = %f = %f Angstrom.\n\n",
	  potential_function, CUTOFF, CUTOFF*ULENGTH_IN_ANGSTROM);
    
    printf ("The configuration is being loaded from \"%s\".\n", fn_config_read); 
    if ((filehandle=fopen(fn_config_read, "r"))==NULL)
    {
	fprintf(stdout, "\n **Error: cannot locate 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));
    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 Angstrom\n",   &H0[0][0]);
    fscanf (filehandle, "H(1,2) = %lf Angstrom\n",   &H0[0][1]);
    fscanf (filehandle, "H(1,3) = %lf Angstrom\n\n", &H0[0][2]);
    fscanf (filehandle, "H(2,1) = %lf Angstrom\n",   &H0[1][0]);
    fscanf (filehandle, "H(2,2) = %lf Angstrom\n",   &H0[1][1]);
    fscanf (filehandle, "H(2,3) = %lf Angstrom\n\n", &H0[1][2]);
    fscanf (filehandle, "H(3,1) = %lf Angstrom\n",   &H0[2][0]);
    fscanf (filehandle, "H(3,2) = %lf Angstrom\n",   &H0[2][1]);
    fscanf (filehandle, "H(3,3) = %lf Angstrom\n\n", &H0[2][2]);
    /* convert to reduced units */
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    H0[i][j] /= ULENGTH_IN_ANGSTROM;

    /* 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_ANGSTROM,
	  H[0][1]*ULENGTH_IN_ANGSTROM,
	  H[0][2]*ULENGTH_IN_ANGSTROM);	  
    printf (" H[A] = | %.5f %.5f %.5f |\n",
	  H[1][0]*ULENGTH_IN_ANGSTROM,
	  H[1][1]*ULENGTH_IN_ANGSTROM,
	  H[1][2]*ULENGTH_IN_ANGSTROM);
    printf ("        | %.5f %.5f %.5f |\n",
	  H[2][0]*ULENGTH_IN_ANGSTROM,
	  H[2][1]*ULENGTH_IN_ANGSTROM,
	  H[2][2]*ULENGTH_IN_ANGSTROM);
    
    printf ("Supercell volume = %f = %f A^3.\n\n", volume,
	    volume*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
    if (np/volume > 0.85*MAXDENSITY)
    {
	printf("Current density %f is above 85%% of MAXDENSITY %f in \"diag.h\".\n",
	       np/volume, MAXDENSITY);
	printf ("Overflow feared. Please increase MAXDENSITY in \"diag.h\".\n");
	exit(1);
    }
    
    /* read in reduced coordinates from configuration file */ 
    fscanf(filehandle,
	   "Mass(amu)     sx        sy        sz     sx1(1/ps)  sy1(1/ps)  sz1(1/ps)\n");
    for (i=0;i<3*np;i+=3)
    {
	fscanf (filehandle, "%lf  %lf  %lf  %lf  %lf  %lf  %lf\n",
		&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",
	    CUTOFF*ULENGTH_IN_ANGSTROM);
    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 (Angstrom 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 (save_dos) 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)
		    {
			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() */


/* 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;
    
    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;
    
    if (LJ6_12)
    {
	ir2 = 1./r2;
	ir4 = ir2 * ir2;
	ir6 = ir4 * ir2;
	ir12 = ir6 * ir6;
	first = (-48.*ir12+24.*ir6)*ir2;
	second = (624.*ir12-168.*ir6)*ir4;
    }
    else if(WOON)
    {
	r1 = sqrt(r2);
	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)/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 CUTOFF 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 CUTOFF, but that thickness is just 1/|G_i|, the
       corresponding reciprocal vector, if G^T * H = I.  */
    nx = ceil(CUTOFF*sqrt(HI[0][0]*HI[0][0]+
			  HI[1][0]*HI[1][0]+HI[2][0]*HI[2][0]));
    ny = ceil(CUTOFF*sqrt(HI[0][1]*HI[0][1]+
			  HI[1][1]*HI[1][1]+HI[2][1]*HI[2][1]));
    nz = ceil(CUTOFF*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<CUTOFF*CUTOFF)
		    {
			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_STRING_LENGTH], *p, *q, *r;
    int number_of_brackets = 0;
    p = sentence;
    while(1)
    {
	if ((*p==0)||(p-sentence>MAX_STRING_LENGTH)) 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() */