/*****************************************/
/*  Atomistic simulation of Argon using  */
/*  pair-potential.                      */
/*                                       */
/*           version 1.1 May.27, 1998    */
/*           developed by Ju Li (MIT)    */
/*****************************************/

/*
  Version Notes:
  1.0  TTY interface;
       Woon's pair potential;
       NEV, NTP ensemble;
       Canonical ensemble
       fluctutation formulas;
       Thermal conductivity;
  1.1  Micro-canonical ensemble
       fluctuation formulas;
       Shear viscosity;
*/

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

void main()
{
    int i,j,k,l;
    char sentence[100]={(char)0};
    double tmp, cc, dd, delta2, Treal, Tnow; 
    double p[3], sumkine[3], factorT, factorHT;
    int kb, kHsum=0, kupdate, kw_lp;
    int scale_T, scale_H, scale_H_last, update;
    double rdel, kinewall, HTnow, dispmx2; 
    double sk[3], xfr, xfi, xf, sxij, syij, szij, rx, ry, rz, r2;
    double HT[3][3], HIT[3][3], Hlast[3][3], G[3][3], H1T[3][3], 
	A[3][3], A1[3][3], GI[3][3], GA[3][3], G1[3][3], fH[3][3]; 
    /* temporary averages */
    double kinetum, Htum[3][3], strestum[3][3];
    double k1tum, k2tum, v1k1tum[3][3], v1k2tum[3][3], v2k2tum[3][3][3][3],
	ck1tum[6][6];
    /* overall averages */
    double kinesum, potesum, Hsum[3][3], stressum[3][3], Tnowsum;
    double k1sum, k2sum, v1k1sum[3][3], v1k2sum[3][3], v2k2sum[3][3][3][3],
	ck1sum[6][6];
    double pressure, dE2, dV2, dEdV, bulk_modulus, alpha, Cv, Cp;
    double msd_kstart;
	
    printf ("\n=====================  Reading in instructions  =====================\n\n");
    
    /* Read instructions from standard input, such as "con" file */
    fscanf (stdin, "Whether or not to calculate the thermal conductivity:\n%d\n\n",
	    &calc_conductivity); 
    if (calc_conductivity)
    {
	/* open file handles for heat currents and shear stresses */
	for (i=0;i<3;i++)
	{
	    heat_current[i] = fopen(fn_heat_current[i],"w+");
	    shear_stress[i] = fopen(fn_shear_stress[i],"w+");
	}	    
	printf ("Thermal conductivity and shear viscosity will be calculated for\n");
    }
    
    fscanf (stdin, "Total number of atoms:\n%d\n\n", &np);
    /* Because center of mass is fixed, we "lose" one particle */
    neff = np-1;
    
    fscanf (stdin, "Default geometry of the unit cell (11,12,13,21,22,23,31,32,33), in Angstrom:\n");
    fscanf (stdin, "%lf %lf %lf %lf %lf %lf %lf %lf %lf\n\n",
	    &H[0][0], &H[0][1], &H[0][2],
	    &H[1][0], &H[1][1], &H[1][2],
	    &H[2][0], &H[2][1], &H[2][2]);
    
    fscanf (stdin, "Translation of the unit cell in three dimensions, NC(1,2,3):\n%d %d %d\n\n",
	    &nc[0], &nc[1], &nc[2]);
  
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    H[i][j] = nc[i]*H[i][j]/ULENGTH_IN_ANGSTROM;
    volume = matinv(H, HI);
    
    fscanf (stdin, "Type of structure (CRYSTAL_FCC or LIQUID):\n%s\n\n",  type_of_structure);
    CRYSTAL_FCC = !strcmp(type_of_structure,"CRYSTAL_FCC");
    LIQUID = !strcmp(type_of_structure,"LIQUID");
    if (!CRYSTAL_FCC && !LIQUID)
    {
	fprintf (stdout, "\n **Error: Currently the structure type must either CRYSTAL_FCC or LIQUID\n");
	return;
  }
    else
    {
      printf ("Argon in %s structure with %d atoms in supercell\n", type_of_structure, np);
      printf ("              | %lf %f %f |\n",
	      H[0][0]*ULENGTH_IN_ANGSTROM,
	      H[0][1]*ULENGTH_IN_ANGSTROM,
	      H[0][2]*ULENGTH_IN_ANGSTROM);
      printf ("H(angstrom) = | %f %f %f |\n",
	      H[1][0]*ULENGTH_IN_ANGSTROM,
	      H[1][1]*ULENGTH_IN_ANGSTROM,
	      H[1][2]*ULENGTH_IN_ANGSTROM);
      printf ("              | %f %f %f |\n",
	      H[2][0]*ULENGTH_IN_ANGSTROM,
	      H[2][1]*ULENGTH_IN_ANGSTROM, 
	      H[2][2]*ULENGTH_IN_ANGSTROM);
    }
    printf ("(reduced density = %f = %.1f mol/m^3)\n",  np/volume,
	    np/volume/AVO/ULENGTH_IN_METER/ULENGTH_IN_METER/ULENGTH_IN_METER);
    fscanf (stdin, "Choose 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");
      return;
    }
    else
	printf ("using %s pair potential,", potential_function);
  
    fscanf (stdin, "NEH or NTP mode:\n%s\n\n", sentence);
    NEH = !strcmp(sentence,"NEH");
    NTP = !strcmp(sentence,"NTP");
    if (!NEH && !NTP)
    {
	fprintf (stdout, "\n **Error: no %s mode\n", sentence);
	fprintf (stdout, "\n **Error: Must input a running mode: NEH or NTP\n");
	return;
    }
    else
	printf (" running under %s mode. \n\n", sentence);
  
    fscanf (stdin, "Read from configuration file (if 'NULL' then start anew):\n%s\n\n",
	    fn_config_read);
    read_from_config = strcmp(fn_config_read, "NULL");
    if (read_from_config)
      printf ("The initial configuration will be loaded from file \"%s\"\n", fn_config_read);
    else
	printf ("The initial configuration will be constructed from scratch\n");
    
    fscanf (stdin, "Write on configuration file:\n%s\n\n", fn_config_write);
    printf ("and it will be saved on \"%s\" in the end.\n\n", fn_config_write);
    
    fscanf (stdin, "Real temperature in Kelvin:\n%lf\n\n", &TR);
    T = T_real_to_T_md (TR, &dTdR);
    scaler = BOLZ*T/UENERGY_IN_JOULE; /* desired kT in reduced energy unit */
    
    printf ("The system temperature (real) will be maintained at T = %.2f K\n", TR);
    printf ("or %f in reduced unit, which corresponds to %.2f K\n", BOLZ*TR/UENERGY_IN_JOULE, T);
    printf ("in MD temperature with dT_md/dT_real = %.5f\n\n", dTdR);
    
    fscanf (stdin, "Total number of MD steps:\n%ld\n\n", &maxkb);
    printf ("maxkb (total number of MD steps) = %d\n", maxkb);
    
    fscanf (stdin, "DELTA (integration timestep) in ps:\n%lf\n\n", &delta);
    printf ("delta (integration timestep in ps) = %.4f\n", delta);
    delta /= UTIME_IN_PS;
    delta2 = delta*delta;
    
    fscanf (stdin, "TAUT (particle temperature coupling time constant), in ps:\n%lf\n\n",	&tauT);
    printf ("tauT (particle temperature coupling time constant in ps) = %.2f\n", tauT);
    tauT /= UTIME_IN_PS;
    
    fscanf (stdin, "TAUH (wall temperature coupling time constant), in ps:\n%lf\n\n", &tauH);
    printf ("tauH (wall temperature coupling time constant in ps) = %.2f\n", tauH);
    tauH /= UTIME_IN_PS;
    
    fscanf (stdin, "Wall mass (in Argon mass):\n%lf\n\n",	&wallmass);
    printf ("wallmass (in Argon mass) = %.2f\n\n", wallmass);
    wallmass *= ARGON_MASS/UMASS_IN_KG;
      
    fscanf (stdin, "Step at which stops temperature coupling and starts property averaging\n");
    fscanf (stdin, "(also starts to write J files if calculating the thermal conductivity):\n%ld\n\n",
	    &kstart);
    printf ("At step %d, we will stop temperature coupling and start property averaging\n", kstart);
    if (calc_conductivity)
    {
	printf ("and also begin to save heat current data on \"%s\", \"%s\", \"%s\",\n",
		fn_heat_current[0], fn_heat_current[1], fn_heat_current[2]);
	printf ("and shear stress data on \"%s\", \"%s\", \"%s\".\n",
		fn_shear_stress[0], fn_shear_stress[1], fn_shear_stress[2]);
    }	
    printf ("\n");
    
    fscanf (stdin, "Frequency of running conditions ouput on screen and disk:\n%d\n\n", &kw_lp);
    printf ("Every %d steps, we will print some properties on screen and also \nsave to \"%s\", \"%s\".\n\n",
	    kw_lp, fn_property_output, fn_structure_output);
    
    fscanf (stdin, "Random number generator seed (long integer):\n%ld\n\n", &iseed);
    printf ("Random number generator seed is %ld\n", iseed);
    srandom (iseed);
    
    fscanf (stdin, "Default pressure, in MPa:\n%lf\n\n", &pout[0][0]);
    pout[0][0] /= USTRESS_IN_MPA;
    if (NTP)
    {
	printf ("(NTP mode) The system will be under pressure %.3f MPa\n",
		pout[0][0]*USTRESS_IN_MPA);
	pout[1][1] = pout[0][0];
	pout[2][2] = pout[0][0];
	pout[0][1] = 0.;
	pout[0][2] = 0.;
	pout[1][0] = 0.;
	pout[1][2] = 0.;
	pout[2][0] = 0.;
	pout[2][1] = 0.;
    }
    
    printf ("\n=====================  Reading in complete  =========================\n\n");
    
    /* allocate memory for global variables of the run */
    f    = (double *) malloc (3*np*sizeof(double));
    fs   = (double *) malloc (3*np*sizeof(double));
    s    = (double *) malloc (3*np*sizeof(double));
    v    = (double *) malloc (3*np*sizeof(double));
    s1   = (double *) malloc (3*np*sizeof(double));
    s2   = (double *) malloc (3*np*sizeof(double));
    s3   = (double *) malloc (3*np*sizeof(double));
    s4   = (double *) malloc (3*np*sizeof(double));
    s5   = (double *) malloc (3*np*sizeof(double));
    x    = (double *) malloc (3*np*sizeof(double));
    xold = (double *) malloc (3*np*sizeof(double));
    d    = (double *) malloc (3*np*sizeof(double));
    dold = (double *) malloc (3*np*sizeof(double));
  
    ep = (double *) malloc (np*sizeof(double));
    mass = (double *) malloc (np*sizeof(double));

    idx = (int *) malloc (np*sizeof(int));
    list = (int **) malloc (np*sizeof(int *));
    for (i=0;i<np;i++)
	if ((list[i] = (int *) malloc(MAX_NEIGHBOR*sizeof(int)))==NULL)
	    fprintf (stdout, "\n **Error: Not enough memory for neighbor list\n");
    printf ("Memory allocation complete\n");
    
    if (read_from_config)
    {
	/* read from configuration file: */
	read_config (fn_config_read);
	printf ("Loading initial configuration from \"%s\" complete\n", fn_config_read);
    }
    else
    {
	assign_fcc_lattice();
	printf ("Constructing initial configuration as fcc lattice\n");
	/* initialize relevant quantities for dynamics */
	for (i=0;i<3;i++)
	    for (j=0;j<3;j++)
	    {
		H1[i][j] = 0.;
		H2[i][j] = 0.;
		H3[i][j] = 0.;
		H4[i][j] = 0.;
		H5[i][j] = 0.;
	    }
	for (i=0;i<np;i++) 
	    mass[i] = ARGON_MASS/UMASS_IN_KG;
	for (i=0;i<3*np;i++)
	{
	    s1[i] = 0.;
	    s2[i] = 0.;
	    s3[i] = 0.;
	    s4[i] = 0.;
	    s5[i] = 0.;
      }
    }
    
    /* if (LIQUID && read_from_config) goto clear; */
    if (read_from_config) goto clear;
    
    /*******************************************************/
    /* assign initial velocities to atoms (very important) */
    /*******************************************************/
    printf ("Assigning initial velocities to atoms\n");
    for (i=0;i<3;i++)
    {
      p[i] = 0.;
      sumkine[i] = 0.;
    }
    totalmass = 0.;
    
    /* calculate the total momentum */
    randnorm (v,3*np);
    for (i=0;i<3*np;i++)
    {
	v[i] /= sqrt(mass[i/3]);
	p[i%3] += v[i]*mass[i/3];
	totalmass += mass[i/3];
    }
    totalmass /= 3.;
    
    /* get the average velocity */
    for (i=0;i<3;i++) p[i] /= totalmass;  

    /* scale velocities so that total momentum = 0. */
    for (i=0;i<3*np;i++) v[i] -= p[i%3];
    
    /* calculate the kinetic energies in each direction */
    for (i=0;i<3*np;i++) sumkine[i%3] += v[i]*v[i]*mass[i/3]/2.;
    
    /* scale velocities to desired temperature: the center of mass is
       stationary so we lose three degrees of freedom; phonon density is
       made isotropic by giving equal kinetic energy in each direction;
       and we start with twice the desired kinetic energy since
       initially we are at potential energy minima. */
    
    for (i=0;i<3;i++) sumkine[i] /= 2*0.5*neff*scaler; 
    for (i=0;i<3*np;i++)
    {
	v[i] /= sqrt(sumkine[i%3]);
    }
    
    /* calculate velocities*delta in scaled space */
    for (i=0;i<3*np;i+=3)
    {
	s1[i]   = (HI[0][0]*v[i]+HI[1][0]*v[i+1]+HI[2][0]*v[i+2])*delta;
	s1[i+1] = (HI[0][1]*v[i]+HI[1][1]*v[i+1]+HI[2][1]*v[i+2])*delta;
	s1[i+2] = (HI[0][2]*v[i]+HI[1][2]*v[i+1]+HI[2][2]*v[i+2])*delta;
    }
    
 clear:   /* set zero the property cumulants */
    
    kinetum = 0.;  
    k1tum = 0.;
    k2tum = 0.;
    
    kinesum = 0.;
    potesum = 0.;
    k1sum = 0.;
    k2sum = 0.;
    
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	{
	    Htum[i][j] = 0.; 
	    strestum[i][j] = 0.; 
	    v1k1tum[i][j] = 0.;
	    v1k2tum[i][j] = 0.;
	    
	    Hsum[i][j] = 0.;  
	    stressum[i][j] = 0.;
	    v1k1sum[i][j] = 0.;
	    v1k2sum[i][j] = 0.;
	    for (k=0;k<3;k++)
		for (l=0;l<3;l++)
		{
		    v2k2tum[i][j][k][l] = 0.;
		    v2k2sum[i][j][k][l] = 0.;
		}
	    Hlast[i][j] = 0.;  
	}
    
    for (i=0;i<6;i++)
	for (j=0;j<6;j++)
	{
	    ck1tum[i][j] = 0.;
	    ck1sum[i][j] = 0.; 
	}  
    
    update = 1;
    r2del = pow(volume,2./3.)/NRDFD/4.;
    for (i=0; i<NRDFD; i++) igr[i] = 0;
    kupdate = 0;
    
    /* calculate the real space positions */
    for (i=0;i<3*np;i+=3)
    {
	xold[i]   = H[0][0]*s[i]+H[1][0]*s[i+1]+H[2][0]*s[i+2];
	xold[i+1] = H[0][1]*s[i]+H[1][1]*s[i+1]+H[2][1]*s[i+2];
	xold[i+2] = H[0][2]*s[i]+H[1][2]*s[i+1]+H[2][2]*s[i+2];
    }
    
    for (i=0; i<3*np; i++)
    {
	d[i] = 0.;
	dold[i] = 0.;
    }

    property_output = fopen (fn_property_output, "w+");
    structure_output = fopen (fn_structure_output, "w+");

    printf ("Main loop starts\n\n");
    /********************/
    /* main loop starts */
    /********************/

    for (kb=1;kb<=maxkb;kb++)
    {
	scale_H_last = scale_H && (kb > 1);
	scale_H = NTP && ( !(calc_conductivity && (kb>kstart)) );
	scale_T = (NEH && (kb<=kstart)) || scale_H;
      
	/* predictor */
	for (i=0;i<3*np;i++)
	{
	    s[i]  = s[i]+s1[i]+s2[i]+s3[i]+s4[i]+s5[i];
	    s1[i] = s1[i]+2.*s2[i]+3.*s3[i]+4.*s4[i]+5.*s5[i];
	    s2[i] = s2[i]+3.*s3[i]+6.*s4[i]+10.*s5[i];
	    s3[i] = s3[i]+4.*s4[i]+10.*s5[i];
	    s4[i] = s4[i]+5.*s5[i];
	}

	kinewall = 0.;
	if (scale_H)
	{
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		{
		    H[i][j]  = H[i][j]+H1[i][j]+H2[i][j]+H3[i][j]+H4[i][j]+H5[i][j];
		    H1[i][j] = H1[i][j]+2.*H2[i][j]+3.*H3[i][j]+4.*H4[i][j]+5.*H5[i][j];
		    H2[i][j] = H2[i][j]+3.*H3[i][j]+6.*H4[i][j]+10.*H5[i][j];
		    H3[i][j] = H3[i][j]+4.*H4[i][j]+10.*H5[i][j];
		    H4[i][j] = H4[i][j]+5.*H5[i][j];
		    kinewall += wallmass*H1[i][j]*H1[i][j]/delta2/2.;
		}
	    /* scale the wall by kT/20, 9kT/2 is too big    */
	    /* cell configuration should change adiabaticly */
	    HTnow = kinewall/(scaler/20.)*T;
	    volume = matinv (H,HI);
	}
	
	if (kb > 1)
	{
	    /* check if we need to update neighbor list */
	    tmp = 0.;
	    for (i=0;i<3*np;i+=3)
	    {
		/* calculate maximum displacements of atoms */
		cc = (d[i]-dold[i])*(d[i]-dold[i])
		    +(d[i+1]-dold[i+1])*(d[i+1]-dold[i+1])
		    +(d[i+2]-dold[i+2])*(d[i+2]-dold[i+2]);
		if (cc>tmp) tmp=cc;
	    }
	    update = (4.05*tmp >= (RLIST-CUTOFF)*(RLIST-CUTOFF));
	}
	
	/***********************************************/
	/* update the Verlet list, calculate g(r) and  */
	/* <111> structure factor of the system (which */
	/* should be 1 in a perfect lattice).          */
	/***********************************************/
	
	if (update)
	{
	    kupdate++;
	    sk[0] = 2.*PI*nc[0];
	    sk[1] = 2.*PI*nc[1];
	    sk[2] = 2.*PI*nc[2];
	    xfr = 0.;
	    xfi = 0.;
	    for (i=0; i<np; i++) idx[i] = 0;
	    for (i=0; i<3*np; i+=3)
	    {
		xfr += cos(sk[0]*s[i]+sk[1]*s[i+1]+sk[2]*s[i+2]);
		xfi += sin(sk[0]*s[i]+sk[1]*s[i+1]+sk[2]*s[i+2]);
		for (j=i+3; j<3*np; j+=3)
		{
		    sxij = s[i]   - s[j];
		    syij = s[i+1] - s[j+1];
		    szij = s[i+2] - s[j+2];
		    if (sxij>0.5)  sxij--;
		    else if (sxij<-0.5) sxij++;
		    if(syij>0.5)  syij--;
		    else if (syij<-0.5) syij++;
		    if (szij>0.5)  szij--;
		    else if (szij<-0.5) szij++;
		    rx = H[0][0]*sxij + H[1][0]*syij + H[2][0]*szij;
		    ry = H[0][1]*sxij + H[1][1]*syij + H[2][1]*szij;
		    rz = H[0][2]*sxij + H[1][2]*syij + H[2][2]*szij;
		    r2 = rx*rx+ry*ry+rz*rz;
		    if (r2 < RLIST*RLIST)
		    {
			/* uniformly distribute the list */
			if ((i+j)&1)
			{
			    *(list[i/3]+idx[i/3]) = j;
			    idx[i/3]++;
			}
			else
			{
			    *(list[j/3]+idx[j/3]) = i;
			    idx[j/3]++;
			}
		    }
		    k = floor(r2/r2del);
		    if (k<NRDFD) igr[k]++;
		}
	    }
	    tmp = 0.; k = 0;
	    for (i=0; i<np; i++)
	    {
		if (idx[i]>MAX_NEIGHBOR)
		{
		    fprintf (stdout, "\n **Error: neighbor list out of bound for atom %d\n", i);
		    exit(1);
		}
		tmp += idx[i];
		if (idx[i]>k) k=idx[i];
	    }
	    /* printf ("\nAt step %d, the neighbor list is updated,\n", kb); */
	    /* printf ("average occupancy = %.1f (%3.1f percent), maximum occupancy = %d\n", */
	    /* tmp/np, tmp/np/MAX_NEIGHBOR*100., k); */
	    xf = sqrt(xfr*xfr+xfi*xfi)/np;
	    for (i=0; i<3*np; i++) dold[i] = d[i];
	}
	
	/* interatomic potential function */
	pair_potential();
	Tnow = kine*UENERGY_IN_JOULE/1.5/neff/BOLZ;
      
	/* TEMPORARY fluctuation averages, */
	/* released every KW_LP steps.     */
	kinetum += kine;
	k1tum += 1./kine;
	k2tum += 1./kine/kine;
	for (i=0;i<3;i++)
	    for (j=0;j<3;j++)
	    {
		Htum[i][j] += H[i][j]; 
		strestum[i][j] += stress[i][j]; 
		v1k1tum[i][j] += virial[i][j]/kine*sqrt(volume);
		v1k2tum[i][j] += virial[i][j]/kine/kine*sqrt(volume);
		for (k=0;k<3;k++)
		    for (l=0;l<3;l++)
			v2k2tum[i][j][k][l] += virial[i][j]*virial[k][l]/kine/kine*volume;
	    }
	for (i=0;i<6;i++)
	    for (j=0;j<6;j++)
		ck1tum[i][j] += c[i][j]/kine;
	
	/* OVERALL fluctuation averages,      */
	/* released at end of the simulation. */
	if (kb > kstart)
	{
	    potesum += pote;
	    kinesum += kine;
	    k1sum += 1./kine;
	    k2sum += 1./kine/kine;
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		{
		    Hsum[i][j] += H[i][j];
		    stressum[i][j] += stress[i][j];
		    v1k1sum[i][j] += virial[i][j]/kine*sqrt(volume);
		    v1k2sum[i][j] += virial[i][j]/kine/kine*sqrt(volume);
		    for (k=0;k<3;k++)
			for (l=0;l<3;l++)
			    v2k2sum[i][j][k][l] += virial[i][j]*virial[k][l]/kine/kine*volume;
		}
	    for (i=0;i<6;i++)
		for (j=0;j<6;j++)
		    ck1sum[i][j] += c[i][j]/kine;
	}
	
	if (kb%kw_lp==0)
	{
	    /* screen output and write to data file, also */
	    /* calculate instantaneous elastic constants  */
	    kinetum /= kw_lp;
	    k1tum /= kw_lp;
	    k2tum /= kw_lp;
	    cc = kinetum/(1.5*neff); /* average kT during period */
	    tmp = matinv(Htum,A)/kw_lp/kw_lp/kw_lp; /* average volume */
	    
	    /* Born contribution */
	    for (i=0;i<6;i++)
		for (j=0;j<6;j++)
		    ck1tum[i][j] *= (1.5*neff-1)*cc/kw_lp;
	    
	    /* ideal gas contribution */
	    ck1tum[0][0] += 2*neff*cc/tmp;
	    ck1tum[1][1] += 2*neff*cc/tmp;
	    ck1tum[2][2] += 2*neff*cc/tmp;
	    ck1tum[3][3] += neff*cc/tmp;
	    ck1tum[4][4] += neff*cc/tmp;
	    ck1tum[5][5] += neff*cc/tmp;
	    
	    dE2 = (1.5*neff-1)*(1.5*neff-2)*k2tum - 
		(1.5*neff-1)*(1.5*neff-1)*k1tum*k1tum;
	    
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		{
		    A[i][j] =
			 -(1.5*neff-1)*(1.5*neff-1)*k1tum*v1k1tum[i][j]/kw_lp
			 +(1.5*neff-1)*(1.5*neff-2)*v1k2tum[i][j]/kw_lp; 
		    v1k1tum[i][j] *= (1.5*neff-1)*sqrt(cc)/kw_lp;
		    for (k=0;k<3;k++)
			for (l=0;l<3;l++)		  
			    v2k2tum[i][j][k][l] *= (1.5*neff-1)*(1.5*neff-2)*cc/kw_lp;
		}

	    ck1tum[0][0] += v1k1tum[0][0]*v1k1tum[0][0]+A[0][0]*A[0][0]*cc/dE2-v2k2tum[0][0][0][0]; 
	    ck1tum[0][1] += v1k1tum[0][0]*v1k1tum[1][1]+A[0][0]*A[1][1]*cc/dE2-v2k2tum[0][0][1][1]; 
	    ck1tum[0][2] += v1k1tum[0][0]*v1k1tum[2][2]+A[0][0]*A[2][2]*cc/dE2-v2k2tum[0][0][2][2]; 
	    ck1tum[1][1] += v1k1tum[1][1]*v1k1tum[1][1]+A[1][1]*A[1][1]*cc/dE2-v2k2tum[1][1][1][1]; 
	    ck1tum[1][2] += v1k1tum[1][1]*v1k1tum[2][2]+A[1][1]*A[2][2]*cc/dE2-v2k2tum[1][1][2][2];
	    ck1tum[2][2] += v1k1tum[2][2]*v1k1tum[2][2]+A[2][2]*A[2][2]*cc/dE2-v2k2tum[2][2][2][2]; 
	    
	    ck1tum[0][3] += v1k1tum[0][0]*v1k1tum[1][2]+A[0][0]*A[1][2]*cc/dE2-v2k2tum[0][0][1][2]; 
	    ck1tum[0][4] += v1k1tum[0][0]*v1k1tum[0][2]+A[0][0]*A[0][2]*cc/dE2-v2k2tum[0][0][0][2]; 
	    ck1tum[0][5] += v1k1tum[0][0]*v1k1tum[0][1]+A[0][0]*A[0][1]*cc/dE2-v2k2tum[0][0][0][1];
	    ck1tum[1][3] += v1k1tum[1][1]*v1k1tum[1][2]+A[1][1]*A[1][2]*cc/dE2-v2k2tum[1][1][1][2]; 
	    ck1tum[1][4] += v1k1tum[1][1]*v1k1tum[0][2]+A[1][1]*A[0][2]*cc/dE2-v2k2tum[1][1][0][2]; 
	    ck1tum[1][5] += v1k1tum[1][1]*v1k1tum[0][1]+A[1][1]*A[0][1]*cc/dE2-v2k2tum[1][1][0][1];
	    ck1tum[2][3] += v1k1tum[2][2]*v1k1tum[1][2]+A[2][2]*A[1][2]*cc/dE2-v2k2tum[2][2][1][2]; 
	    ck1tum[2][4] += v1k1tum[2][2]*v1k1tum[0][2]+A[2][2]*A[0][2]*cc/dE2-v2k2tum[2][2][0][2]; 
	    ck1tum[2][5] += v1k1tum[2][2]*v1k1tum[0][1]+A[2][2]*A[0][1]*cc/dE2-v2k2tum[2][2][0][1];
	    
	    ck1tum[3][3] += v1k1tum[1][2]*v1k1tum[1][2]+A[1][2]*A[1][2]*cc/dE2-v2k2tum[1][2][1][2]; 
	    ck1tum[3][4] += v1k1tum[1][2]*v1k1tum[0][2]+A[1][2]*A[0][2]*cc/dE2-v2k2tum[1][2][0][2]; 
	    ck1tum[3][5] += v1k1tum[1][2]*v1k1tum[0][1]+A[1][2]*A[0][1]*cc/dE2-v2k2tum[1][2][0][1]; 
	    ck1tum[4][4] += v1k1tum[0][2]*v1k1tum[0][2]+A[0][2]*A[0][2]*cc/dE2-v2k2tum[0][2][0][2]; 
	    ck1tum[4][5] += v1k1tum[0][2]*v1k1tum[0][1]+A[0][2]*A[0][1]*cc/dE2-v2k2tum[0][2][0][1]; 
	    ck1tum[5][5] += v1k1tum[0][1]*v1k1tum[0][1]+A[0][1]*A[0][1]*cc/dE2-v2k2tum[0][1][0][1]; 
	    
	    for (i=0;i<6;i++)
		for (j=i+1;j<6;j++)
		    ck1tum[j][i] = ck1tum[i][j];
	    
	    pressure = (strestum[0][0]+strestum[1][1]+strestum[2][2])/kw_lp/3.;
	    if (NTP) hamiltonian = pote+kine+pout[0][0]*volume+kinewall;
	    else hamiltonian = pote+kine;
	    
	    fprintf (property_output, "%7d %.3f %.6f %.6f %.5f\n",
		     kb, Tnow, pote/np, hamiltonian/np, pressure);

	    printf ("\n\nKB = %ld\n", kb);
	    printf ("T_md = %f K, T_real = %f K\n", Tnow, T_md_to_T_real(Tnow));
	    printf ("Potential energy = %f = %f eV = %f K = %.1f J/mol\n",
		    pote/np, pote/np*UENERGY_IN_JOULE/EV_IN_JOULE,
		    pote/np*UENERGY_IN_JOULE/BOLZ, pote*UENERGY_IN_JOULE*AVO/np);
	    printf ("Hamiltonian = %f = %f eV = %f K = %.1f J/mol\n",
		    hamiltonian/np, hamiltonian/np*UENERGY_IN_JOULE/EV_IN_JOULE,
		    hamiltonian/np*UENERGY_IN_JOULE/BOLZ, hamiltonian*UENERGY_IN_JOULE*AVO/np);
	    printf ("Pressure = %f = %f MPa\n", pressure, pressure*USTRESS_IN_MPA);
	    printf ("Density = %f = %f 1/A^3 = %f mol/m^3\n\n", np/volume,
		    np/volume/ULENGTH_IN_ANGSTROM/ULENGTH_IN_ANGSTROM/ULENGTH_IN_ANGSTROM, 
		    np/volume/AVO/ULENGTH_IN_METER/ULENGTH_IN_METER/ULENGTH_IN_METER);
	    
	    printf ("       | %.4f %.4f %.4f |",
		    Htum[0][0]*ULENGTH_IN_ANGSTROM/kw_lp,
		    Htum[0][1]*ULENGTH_IN_ANGSTROM/kw_lp,
		    Htum[0][2]*ULENGTH_IN_ANGSTROM/kw_lp);
	    printf ("            | %.4f %.4f %.4f |\n",
		    strestum[0][0]*USTRESS_IN_MPA/kw_lp,
		    strestum[0][1]*USTRESS_IN_MPA/kw_lp,
		    strestum[0][2]*USTRESS_IN_MPA/kw_lp);	  
	    printf ("H(A) = | %.4f %.4f %.4f |",
		    Htum[1][0]*ULENGTH_IN_ANGSTROM/kw_lp,
		    Htum[1][1]*ULENGTH_IN_ANGSTROM/kw_lp,
		    Htum[1][2]*ULENGTH_IN_ANGSTROM/kw_lp);
	    printf ("   t(MPa) = | %.4f %.4f %.4f |\n",
		    strestum[1][0]*USTRESS_IN_MPA/kw_lp,
		    strestum[1][1]*USTRESS_IN_MPA/kw_lp,
		    strestum[1][2]*USTRESS_IN_MPA/kw_lp);
	    printf ("       | %.4f %.4f %.4f |",
		    Htum[2][0]*ULENGTH_IN_ANGSTROM/kw_lp,
		    Htum[2][1]*ULENGTH_IN_ANGSTROM/kw_lp, 
		    Htum[2][2]*ULENGTH_IN_ANGSTROM/kw_lp);
	    printf ("            | %.4f %.4f %.4f |\n\n",
		    strestum[2][0]*USTRESS_IN_MPA/kw_lp,
		    strestum[2][1]*USTRESS_IN_MPA/kw_lp,
		    strestum[2][2]*USTRESS_IN_MPA/kw_lp);
	    
	    printf ("         | %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f |\n",
		    ck1tum[0][0]*USTRESS_IN_MPA, ck1tum[0][1]*USTRESS_IN_MPA, ck1tum[0][2]*USTRESS_IN_MPA,
		    ck1tum[0][3]*USTRESS_IN_MPA, ck1tum[0][4]*USTRESS_IN_MPA, ck1tum[0][5]*USTRESS_IN_MPA);
	    printf ("         | %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f |\n",
		    ck1tum[1][0]*USTRESS_IN_MPA, ck1tum[1][1]*USTRESS_IN_MPA, ck1tum[1][2]*USTRESS_IN_MPA,
		    ck1tum[1][3]*USTRESS_IN_MPA, ck1tum[1][4]*USTRESS_IN_MPA, ck1tum[1][5]*USTRESS_IN_MPA);
	    printf ("C(MPa) = | %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f |\n",
		    ck1tum[2][0]*USTRESS_IN_MPA, ck1tum[2][1]*USTRESS_IN_MPA, ck1tum[2][2]*USTRESS_IN_MPA,
		    ck1tum[2][3]*USTRESS_IN_MPA, ck1tum[2][4]*USTRESS_IN_MPA, ck1tum[2][5]*USTRESS_IN_MPA);
	    printf ("         | %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f |\n",
		    ck1tum[3][0]*USTRESS_IN_MPA, ck1tum[3][1]*USTRESS_IN_MPA, ck1tum[3][2]*USTRESS_IN_MPA,
		    ck1tum[3][3]*USTRESS_IN_MPA, ck1tum[3][4]*USTRESS_IN_MPA, ck1tum[3][5]*USTRESS_IN_MPA);
	    printf ("         | %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f |\n",
		    ck1tum[4][0]*USTRESS_IN_MPA, ck1tum[4][1]*USTRESS_IN_MPA, ck1tum[4][2]*USTRESS_IN_MPA,
		    ck1tum[4][3]*USTRESS_IN_MPA, ck1tum[4][4]*USTRESS_IN_MPA, ck1tum[4][5]*USTRESS_IN_MPA);
	    printf ("         | %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f |\n",
		    ck1tum[5][0]*USTRESS_IN_MPA, ck1tum[5][1]*USTRESS_IN_MPA, ck1tum[5][2]*USTRESS_IN_MPA,
		    ck1tum[5][3]*USTRESS_IN_MPA, ck1tum[5][4]*USTRESS_IN_MPA, ck1tum[5][5]*USTRESS_IN_MPA);
	    
	    fprintf (structure_output, "%f %f %f %f %f %f %f %f %f ",
		     kb*delta, 
		     Htum[0][0]*ULENGTH_IN_ANGSTROM/kw_lp,
		     Htum[1][1]*ULENGTH_IN_ANGSTROM/kw_lp,
		     Htum[2][2]*ULENGTH_IN_ANGSTROM/kw_lp,
		     Htum[0][1]*ULENGTH_IN_ANGSTROM/kw_lp,
		     Htum[0][2]*ULENGTH_IN_ANGSTROM/kw_lp,
		     Htum[2][1]*ULENGTH_IN_ANGSTROM/kw_lp,
		     xf,
		     msd*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
	    fprintf (structure_output, "%f %f %f %f %f %f",
		     strestum[0][0]*USTRESS_IN_MPA/kw_lp,
		     strestum[1][1]*USTRESS_IN_MPA/kw_lp,
		     strestum[2][2]*USTRESS_IN_MPA/kw_lp,
		     strestum[0][1]*USTRESS_IN_MPA/kw_lp,
		     strestum[0][2]*USTRESS_IN_MPA/kw_lp,
		     strestum[2][1]*USTRESS_IN_MPA/kw_lp);
	    fprintf (structure_output,
		     "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
		     ck1tum[0][0]*USTRESS_IN_MPA, ck1tum[0][1]*USTRESS_IN_MPA, ck1tum[0][2]*USTRESS_IN_MPA,
		     ck1tum[0][3]*USTRESS_IN_MPA, ck1tum[0][4]*USTRESS_IN_MPA, ck1tum[0][5]*USTRESS_IN_MPA,
		     ck1tum[1][1]*USTRESS_IN_MPA, ck1tum[1][2]*USTRESS_IN_MPA, ck1tum[1][3]*USTRESS_IN_MPA,
		     ck1tum[1][4]*USTRESS_IN_MPA, ck1tum[1][5]*USTRESS_IN_MPA, ck1tum[2][2]*USTRESS_IN_MPA,
		     ck1tum[2][3]*USTRESS_IN_MPA, ck1tum[2][4]*USTRESS_IN_MPA, ck1tum[2][5]*USTRESS_IN_MPA,
		     ck1tum[3][3]*USTRESS_IN_MPA, ck1tum[3][4]*USTRESS_IN_MPA, ck1tum[3][5]*USTRESS_IN_MPA,
		     ck1tum[4][4]*USTRESS_IN_MPA, ck1tum[4][5]*USTRESS_IN_MPA, ck1tum[5][5]*USTRESS_IN_MPA);
	    
	    kinetum = 0.;
	    k1tum = 0.;
	    k2tum = 0.;
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		{
		    Htum[i][j] = 0.;
		    strestum[i][j] = 0.; 
		    v1k1tum[i][j] = 0.;
		    v1k2tum[i][j] = 0.;
		    for (k=0;k<3;k++)
			for (l=0;l<3;l++)
			    v2k2tum[i][j][k][l] = 0.;
		}
	    for (i=0;i<6;i++)
		for (j=0;j<6;j++)
		    ck1tum[i][j] = 0.;
	} /* finish common properties output */
	
	for (i=0;i<3*np;i+=3)
	{
	    fs[i]   = (HI[0][0]*f[i]+HI[1][0]*f[i+1]+HI[2][0]*f[i+2])/mass[i/3];
	    fs[i+1] = (HI[0][1]*f[i]+HI[1][1]*f[i+1]+HI[2][1]*f[i+2])/mass[i/3];
	    fs[i+2] = (HI[0][2]*f[i]+HI[1][2]*f[i+1]+HI[2][2]*f[i+2])/mass[i/3];
	}
	
	if (scale_H)
	{
	    /* calculate the generalized forces */
	    matran (H, HT);
	    matran (HI, HIT);
	    matmul (HT, H, G);
	    matinv (G, GI);
	    matran (H1, H1T);
	    matmul (H1T, H, A);
	    matmul (HT, H1, A1);
	    matadd (A, A1, G1);
	    matmul (GI, G1, GA);
	    
	    /* generalized forces on sx,sy,sz */
	    for (i=0;i<3*np;i+=3)
	    {
		fs[i]   -= (GA[0][0]*s1[i]+
			    GA[1][0]*s1[i+1]+
			    GA[2][0]*s1[i+2])/delta2;
		fs[i+1] -= (GA[0][1]*s1[i]+
			    GA[1][1]*s1[i+1]+
			    GA[2][1]*s1[i+2])/delta2;
		fs[i+2] -= (GA[0][2]*s1[i]+
			    GA[1][2]*s1[i+1]+
			    GA[2][2]*s1[i+2])/delta2;
	    }
	    
	    matsub (stress, pout, A);
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		    A[i][j] = volume*A[i][j]/wallmass;
	    matmul(A, HIT, fH);
	}
	
	factorT = 1.+delta/2./tauT*((T+1.)/(Tnow+1.)-1.);
	factorHT = 1.+delta/2./tauH*((T+1.)/(HTnow+1.)-1.);
      
	/* corrector */
	cc = delta*delta/2.;
	if (scale_H)
	{
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		{
		    tmp = H2[i][j]-cc*fH[i][j];
		    H[i][j]  -= tmp*F02;
		    H1[i][j] -= tmp*F12;
		    H2[i][j] -= tmp;
		    H3[i][j] -= tmp*F32;
		    H4[i][j] -= tmp*F42;
		    H5[i][j] -= tmp*F52;
		    H1[i][j] *= factorHT;
		}
	    if (kb > kstart/2)
	    {
		/* wait till initial transients are over */
		matadd (H, Hlast, Hlast);
		kHsum++;
	    }
	    volume = matinv (H, HI);
	}
	
	if (scale_H_last && (!scale_H) && (kHsum>0))
	{
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		{
		    H[i][j] = Hlast[i][j]/kHsum;
		    H1[i][j] = 0.;
		    H2[i][j] = 0.;
		    H3[i][j] = 0.;
		    H4[i][j] = 0.;
		    H5[i][j] = 0.;
		}
	    volume = matinv (H, HI);
	}
	
	for (i=0;i<3*np;i++)
	{
	    tmp = s2[i]-cc*fs[i];
	    s[i]  -= tmp*F02;
	    s1[i] -= tmp*F12;
	    s2[i] -= tmp;
	    s3[i] -= tmp*F32;
	    s4[i] -= tmp*F42;
	    s5[i] -= tmp*F52;
	}
	
	msd  = 0.;
	for (i=0;i<3*np;i+=3)
	{
	    x[i]   = H[0][0]*s[i]+H[1][0]*s[i+1]+H[2][0]*s[i+2];
	    x[i+1] = H[0][1]*s[i]+H[1][1]*s[i+1]+H[2][1]*s[i+2];
	    x[i+2] = H[0][2]*s[i]+H[1][2]*s[i+1]+H[2][2]*s[i+2];
	    d[i]   += x[i]-xold[i];
	    d[i+1] += x[i+1]-xold[i+1];
	    d[i+2] += x[i+2]-xold[i+2];
	    msd   += d[i]*d[i]+d[i+1]*d[i+1]+d[i+2]*d[i+2];	  
	}
	msd /= np;
	if (kb==kstart) msd_kstart = msd;
	
	/* periodic boundary conditions */
	for (i=0;i<3*np;i++)
	    if (s[i]<0.) s[i]++;
	    else if (s[i]>1.) s[i]--;
	
	/* recalculate real space positions */
	for (i=0;i<3*np;i+=3)
	{
	    xold[i]   = H[0][0]*s[i]+H[1][0]*s[i+1]+H[2][0]*s[i+2];
	    xold[i+1] = H[0][1]*s[i]+H[1][1]*s[i+1]+H[2][1]*s[i+2];
	    xold[i+2] = H[0][2]*s[i]+H[1][2]*s[i+1]+H[2][2]*s[i+2];
	}
	
	if (scale_T)
	    for (i=0;i<3*np;i++) s1[i] *= factorT;
	
	/* sum up the heat current */
	for (i=0;i<3*np;i++) cj[i%3] += ep[i/3]*v[i];
	if (calc_conductivity && (kb>kstart))
	{
	    /* write into files after kstart */
	    for (i=0;i<3;i++)
	    {
		cj[i] *= sqrt(dTdR*UCONDUCTIVITY_IN_SI*delta/(scaler*T*volume));
		fwrite (&cj[i], sizeof(double), 1, heat_current[i]);
	    }
	    for (i=0;i<3;i++)
		for (j=0;j<3;j++)
		    stress[i][j] *= USTRESS_IN_PASCAL*sqrt(delta*UTIME_IN_SECOND*volume*ULENGTH_IN_METER
					     *ULENGTH_IN_METER*ULENGTH_IN_METER/scaler/UENERGY_IN_JOULE);
	    fwrite (&(stress[1][2]), sizeof(double), 1, shear_stress[0]);
	    fwrite (&(stress[0][2]), sizeof(double), 1, shear_stress[1]);
	    fwrite (&(stress[0][1]), sizeof(double), 1, shear_stress[2]);
	}
    } /* main loop ends; */
    
    /* OVERALL averages */
    tmp = (double) (maxkb-kstart);
    potesum /= tmp;
    kinesum /= tmp;
    k1sum /= tmp;
    k2sum /= tmp;
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	{
	    Hsum[i][j] /= tmp;
	    stressum[i][j] /= tmp;
	    v1k1sum[i][j] /= tmp;
	    v1k2sum[i][j] /= tmp;
	    for (k=0;k<3;k++)
		for (l=0;l<3;l++)		  
		    v2k2sum[i][j][k][l] /= tmp;
	}
    for (i=0;i<6;i++)
	for (j=0;j<6;j++)
	    ck1sum[i][j] /= tmp;
    
    volume = matinv(Hsum,HI);
    cc = kinesum/(1.5*neff); /* average kT during run */
    Treal = T_md_to_T_real(cc*UENERGY_IN_JOULE/BOLZ);
    Tnowsum = T_real_to_T_md (Treal, &dTdR);
    
    dE2 = (1.5*neff-1)*(1.5*neff-2)*k2sum - 
	(1.5*neff-1)*(1.5*neff-1)*k1sum*k1sum;
    dV2 = dEdV = dd = tmp = 0.;
    for (i=0;i<3;i++)
    {
	dV2 += v1k1sum[i][i]/3./sqrt(volume);
	dEdV += (1.5*neff-1)*(1.5*neff-2)*v1k2sum[i][i]/3./sqrt(volume) - 
	    (1.5*neff-1)*(1.5*neff-1)*v1k1sum[i][i]*k1sum/3./sqrt(volume);
	for (k=0;k<3;k++)
	{
	    dd  += v2k2sum[i][i][k][k];
	    tmp += ck1sum[i][k];
	}
    }
    
    dV2 = -neff/volume/volume
	-(1.5*neff-1)*(1.5*neff-1)*dV2*dV2
	+(1.5*neff-1)*(1.5*neff-2)*dd/9./volume
	-2.*(1.5*neff-1)/3./volume*dV2
	-(1.5*neff-1)/9./volume*(tmp-3*dV2);
    
    pressure = (stressum[0][0]+stressum[1][1]+stressum[2][2])/3.;
    bulk_modulus = cc*volume*(dEdV*dEdV/dE2-dV2);
    alpha = dTdR*(pressure-dEdV/dE2)/Tnowsum/bulk_modulus/3.;
    Cv = -scaler/T*(1.5*neff-1)*(1.5*neff-1)*k1sum*k1sum/dE2*dTdR;
    Cp = Cv+9.*volume*Treal*bulk_modulus*alpha*alpha;
    
    /* Born contribution */
    for (i=0;i<6;i++)
	for (j=0;j<6;j++)
	    ck1sum[i][j] *= (1.5*neff-1)*cc;
    
    /* ideal gas contribution */
    ck1sum[0][0] += 2*neff*cc/volume;
    ck1sum[1][1] += 2*neff*cc/volume;
    ck1sum[2][2] += 2*neff*cc/volume;
    ck1sum[3][3] += neff*cc/volume;
    ck1sum[4][4] += neff*cc/volume;
    ck1sum[5][5] += neff*cc/volume;
    
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	{
	    A[i][j] =
		-(1.5*neff-1)*(1.5*neff-1)*k1sum*v1k1sum[i][j]
		+(1.5*neff-1)*(1.5*neff-2)*v1k2sum[i][j];
	    v1k1sum[i][j] *= (1.5*neff-1)*sqrt(cc);
	    for (k=0;k<3;k++)
		for (l=0;l<3;l++)		  
		    v2k2sum[i][j][k][l] *= (1.5*neff-1)*(1.5*neff-2)*cc;
	}
    
    ck1sum[0][0] += v1k1sum[0][0]*v1k1sum[0][0]+A[0][0]*A[0][0]*cc/dE2-v2k2sum[0][0][0][0]; 
    ck1sum[0][1] += v1k1sum[0][0]*v1k1sum[1][1]+A[0][0]*A[1][1]*cc/dE2-v2k2sum[0][0][1][1]; 
    ck1sum[0][2] += v1k1sum[0][0]*v1k1sum[2][2]+A[0][0]*A[2][2]*cc/dE2-v2k2sum[0][0][2][2]; 
    ck1sum[1][1] += v1k1sum[1][1]*v1k1sum[1][1]+A[1][1]*A[1][1]*cc/dE2-v2k2sum[1][1][1][1]; 
    ck1sum[1][2] += v1k1sum[1][1]*v1k1sum[2][2]+A[1][1]*A[2][2]*cc/dE2-v2k2sum[1][1][2][2];
    ck1sum[2][2] += v1k1sum[2][2]*v1k1sum[2][2]+A[2][2]*A[2][2]*cc/dE2-v2k2sum[2][2][2][2]; 
    
    ck1sum[0][3] += v1k1sum[0][0]*v1k1sum[1][2]+A[0][0]*A[1][2]*cc/dE2-v2k2sum[0][0][1][2]; 
    ck1sum[0][4] += v1k1sum[0][0]*v1k1sum[0][2]+A[0][0]*A[0][2]*cc/dE2-v2k2sum[0][0][0][2]; 
    ck1sum[0][5] += v1k1sum[0][0]*v1k1sum[0][1]+A[0][0]*A[0][1]*cc/dE2-v2k2sum[0][0][0][1];
    ck1sum[1][3] += v1k1sum[1][1]*v1k1sum[1][2]+A[1][1]*A[1][2]*cc/dE2-v2k2sum[1][1][1][2]; 
    ck1sum[1][4] += v1k1sum[1][1]*v1k1sum[0][2]+A[1][1]*A[0][2]*cc/dE2-v2k2sum[1][1][0][2]; 
    ck1sum[1][5] += v1k1sum[1][1]*v1k1sum[0][1]+A[1][1]*A[0][1]*cc/dE2-v2k2sum[1][1][0][1];
    ck1sum[2][3] += v1k1sum[2][2]*v1k1sum[1][2]+A[2][2]*A[1][2]*cc/dE2-v2k2sum[2][2][1][2]; 
    ck1sum[2][4] += v1k1sum[2][2]*v1k1sum[0][2]+A[2][2]*A[0][2]*cc/dE2-v2k2sum[2][2][0][2]; 
    ck1sum[2][5] += v1k1sum[2][2]*v1k1sum[0][1]+A[2][2]*A[0][1]*cc/dE2-v2k2sum[2][2][0][1];
    
    ck1sum[3][3] += v1k1sum[1][2]*v1k1sum[1][2]+A[1][2]*A[1][2]*cc/dE2-v2k2sum[1][2][1][2]; 
    ck1sum[3][4] += v1k1sum[1][2]*v1k1sum[0][2]+A[1][2]*A[0][2]*cc/dE2-v2k2sum[1][2][0][2]; 
    ck1sum[3][5] += v1k1sum[1][2]*v1k1sum[0][1]+A[1][2]*A[0][1]*cc/dE2-v2k2sum[1][2][0][1]; 
    ck1sum[4][4] += v1k1sum[0][2]*v1k1sum[0][2]+A[0][2]*A[0][2]*cc/dE2-v2k2sum[0][2][0][2]; 
    ck1sum[4][5] += v1k1sum[0][2]*v1k1sum[0][1]+A[0][2]*A[0][1]*cc/dE2-v2k2sum[0][2][0][1]; 
    ck1sum[5][5] += v1k1sum[0][1]*v1k1sum[0][1]+A[0][1]*A[0][1]*cc/dE2-v2k2sum[0][1][0][1]; 
    
    for (i=0;i<6;i++)
	for (j=i+1;j<6;j++)
	    ck1sum[j][i] = ck1sum[i][j];
    
    for (i=0,totalmass=0.;i<np;i++) totalmass+=mass[i];
    
    /* write in common property output file */
    fprintf (property_output, "\n   This is an MD run of %d Argon particles in %s structure,\n",
	     np, type_of_structure);
    fprintf (property_output, "   using %s pair potential with cutoff radius = %.2f = %.4f A.\n\n",
	     potential_function, CUTOFF, CUTOFF*ULENGTH_IN_ANGSTROM);
    
    fprintf (property_output, "   During the entire run (%ld x %.4f ps)\n", maxkb, delta*UTIME_IN_PS);
    fprintf (property_output, "   %d counts of Verlet list update\n", kupdate);
    fprintf (property_output, "   was observed.\n\n");
    
    fprintf (property_output, "   The running average begins from step %d\n", kstart);
    fprintf (property_output, "                              to   step %d\n\n", maxkb);
    
    fprintf (property_output, "   The average MD temperature is %f = %f K,\n", Tnowsum/T*scaler, Tnowsum);
    fprintf (property_output, "   corresponding to a real temperature of %f K,\n", Treal);
    fprintf (property_output, "   with dT_md / dT_real = %f.\n\n", dTdR);
    
    fprintf (property_output, "   H[1][1] = %8.4f A    stress[1][1] = %9.4f MPa\n",
	     Hsum[0][0]*ULENGTH_IN_ANGSTROM, stressum[0][0]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[1][2] = %8.4f A    stress[1][2] = %9.4f MPa\n",
	     Hsum[0][1]*ULENGTH_IN_ANGSTROM, stressum[0][1]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[1][3] = %8.4f A    stress[1][3] = %9.4f MPa\n",
	     Hsum[0][2]*ULENGTH_IN_ANGSTROM, stressum[0][2]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[2][1] = %8.4f A    stress[2][1] = %9.4f MPa\n",
	     Hsum[1][0]*ULENGTH_IN_ANGSTROM, stressum[1][0]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[2][2] = %8.4f A    stress[2][2] = %9.4f MPa\n",
	     Hsum[1][1]*ULENGTH_IN_ANGSTROM, stressum[1][1]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[2][3] = %8.4f A    stress[2][3] = %9.4f MPa\n",
	     Hsum[1][2]*ULENGTH_IN_ANGSTROM, stressum[1][2]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[3][1] = %8.4f A    stress[3][1] = %9.4f MPa\n",
	     Hsum[2][0]*ULENGTH_IN_ANGSTROM, stressum[2][0]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[3][2] = %8.4f A    stress[3][2] = %9.4f MPa\n",
	     Hsum[2][1]*ULENGTH_IN_ANGSTROM, stressum[2][1]*USTRESS_IN_MPA);
    fprintf (property_output, "   H[3][3] = %8.4f A    stress[3][3] = %9.4f MPa\n\n",
	     Hsum[2][2]*ULENGTH_IN_ANGSTROM, stressum[2][2]*USTRESS_IN_MPA);
    
    fprintf (property_output, "   The averaged density = %f = %.1f mol/m^3\n", np/volume,
	     np/volume/AVO/ULENGTH_IN_METER/ULENGTH_IN_METER/ULENGTH_IN_METER);
    fprintf (property_output, "   The averaged volume per particle = %f = %f A^3\n",
	     volume/np, volume*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM/np);
    fprintf (property_output, "   The averaged pressure is = %f = %f MPa\n\n",
	     pressure, pressure*USTRESS_IN_MPA);
    
    fprintf (property_output,
	     "   The average kinetic energy is %.4f = %.4f eV = %.1f K = %.1f J/mol\n",
	     kinesum/np, kinesum/np*UENERGY_IN_JOULE/EV_IN_JOULE,
	     kinesum/np*UENERGY_IN_JOULE/BOLZ, kinesum*UENERGY_IN_JOULE*AVO/np);
    fprintf (property_output,
	     "   the average potential energy is %.4f = %.4f eV = %.1f K = %.1f J/mol\n",
	     potesum/np, potesum/np*UENERGY_IN_JOULE/EV_IN_JOULE,
	     potesum/np*UENERGY_IN_JOULE/BOLZ, potesum*UENERGY_IN_JOULE*AVO/np);
    potesum += kinesum; 
    fprintf (property_output,
	   "   the average internal energy is %.4f = %.4f eV = %.1f K = %.1f J/mol\n",
	     potesum/np, potesum/np*UENERGY_IN_JOULE/EV_IN_JOULE,
	     potesum/np*UENERGY_IN_JOULE/BOLZ, potesum*UENERGY_IN_JOULE*AVO/np);
    potesum += pressure * volume; 
    fprintf (property_output,
	     "   the average enthalpy is %.4f = %.4f eV = %.1f K = %.1f J/mol\n\n",
	     potesum/np, potesum/np*UENERGY_IN_JOULE/EV_IN_JOULE,
	     potesum/np*UENERGY_IN_JOULE/BOLZ, potesum*UENERGY_IN_JOULE*AVO/np);
    
    fprintf (property_output, "            | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n",
	     ck1sum[0][0]*USTRESS_IN_MPA, ck1sum[0][1]*USTRESS_IN_MPA, ck1sum[0][2]*USTRESS_IN_MPA,
	     ck1sum[0][3]*USTRESS_IN_MPA, ck1sum[0][4]*USTRESS_IN_MPA, ck1sum[0][5]*USTRESS_IN_MPA);
    fprintf (property_output, "            | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n",
	     ck1sum[1][0]*USTRESS_IN_MPA, ck1sum[1][1]*USTRESS_IN_MPA, ck1sum[1][2]*USTRESS_IN_MPA,
	     ck1sum[1][3]*USTRESS_IN_MPA, ck1sum[1][4]*USTRESS_IN_MPA, ck1sum[1][5]*USTRESS_IN_MPA);
    fprintf (property_output, "   C(MPa) = | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n",
	     ck1sum[2][0]*USTRESS_IN_MPA, ck1sum[2][1]*USTRESS_IN_MPA, ck1sum[2][2]*USTRESS_IN_MPA,
	     ck1sum[2][3]*USTRESS_IN_MPA, ck1sum[2][4]*USTRESS_IN_MPA, ck1sum[2][5]*USTRESS_IN_MPA);
    fprintf (property_output, "            | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n",
	     ck1sum[3][0]*USTRESS_IN_MPA, ck1sum[3][1]*USTRESS_IN_MPA, ck1sum[3][2]*USTRESS_IN_MPA,
	     ck1sum[3][3]*USTRESS_IN_MPA, ck1sum[3][4]*USTRESS_IN_MPA, ck1sum[3][5]*USTRESS_IN_MPA);
    fprintf (property_output, "            | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n",
	     ck1sum[4][0]*USTRESS_IN_MPA, ck1sum[4][1]*USTRESS_IN_MPA, ck1sum[4][2]*USTRESS_IN_MPA,
	     ck1sum[4][3]*USTRESS_IN_MPA, ck1sum[4][4]*USTRESS_IN_MPA, ck1sum[4][5]*USTRESS_IN_MPA);
    fprintf (property_output, "            | %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f |\n\n",
	     ck1sum[5][0]*USTRESS_IN_MPA, ck1sum[5][1]*USTRESS_IN_MPA, ck1sum[5][2]*USTRESS_IN_MPA,
	     ck1sum[5][3]*USTRESS_IN_MPA, ck1sum[5][4]*USTRESS_IN_MPA, ck1sum[5][5]*USTRESS_IN_MPA);
    
    fprintf (property_output, "   The finite stress bulk modulus BT = (c11+2*c12+p)/3 = %f MPa\n",
	     bulk_modulus*USTRESS_IN_MPA);
    fprintf (property_output, "   The line thermal expansion coefficient alpha = %f/K\n", alpha);
    fprintf (property_output, "   The heat capacity Cv = %f J/g/K = %f J/mol/K\n",
	     Cv*UENERGY_IN_JOULE/(totalmass*UMASS_IN_KG*1000.),
	     Cv*UENERGY_IN_JOULE/neff*AVO);
    fprintf (property_output, "                     Cp = %f J/g/K = %f J/mol/K\n",  
	     Cp*UENERGY_IN_JOULE/(totalmass*UMASS_IN_KG*1000.),
	     Cp*UENERGY_IN_JOULE/neff*AVO);
    fprintf (property_output, "\n   The self-diffusion coefficient D = %f 10^(-10) m^2/s\n",
	     (msd-msd_kstart)*ULENGTH_IN_METER*ULENGTH_IN_METER/
	     ((maxkb-kstart)*delta*UTIME_IN_SECOND)/6.*1E10);
    /*******************************************************************************************/
    
    fclose (structure_output);
    
    gr_output = fopen(fn_gr_output, "w+");
    for (i=0; i<NRDFD; i++)
    {
	cc = 2*PI*sqrt((i+0.5)*r2del)*r2del;
	tmp = kupdate*np*cc*(np-1)/volume/2.;
	fprintf (gr_output, "%lf %lf\n", 
		 sqrt((i+0.5)*r2del)*ULENGTH_IN_ANGSTROM, igr[i]/tmp); 
    }
    fclose (gr_output);
    
    /* write on configuration file; */
    write_config(fn_config_write);
    
    /* free all allocated arrays */
    free (f);
    free (fs);
    free (s);
    free (v);
    free (s1);
    free (s2);
    free (s3);
    free (s4);
    free (s5);
    free (x);
    free (xold);
    free (d);
    free (dold);
    
    free (ep);
    free (mass);
    
    free (idx);
    for (i=0;i<np;i++) free (list[i]);
    free (list);
    
    if (calc_conductivity)
    {
	for (i=0;i<3;i++)
	{
	    cc = delta*UTIME_IN_PS;
	    /* the last number is the integration stepsize in picoseconds */
	    fwrite (&cc, sizeof(double), 1, heat_current[i]);
	    fclose (heat_current[i]);
	    fwrite (&cc, sizeof(double), 1, shear_stress[i]);
	    fclose (shear_stress[i]);
	}
	/* subroutine that reads in the shear stress disk file */
	/* and calculate the shear viscosity using FFT */
	shear_smile (property_output);
	/* subroutine that reads in the heat current disk file */
	/* and calculate the thermal conductivity using FFT */
	smile (property_output);
    }
    
    fclose (property_output);
    return;
} /* end main() */


void pair_potential()
{
    /* interaction potential */
    int i, j, k, l;
    double residue, tmp;
    double sxij, syij, szij, rx, ry, rz;
    double r, r2, r6, r12, vij, ff, ffx, ffy, ffz, cc, a, b;
    
    if (LJ6_12)
	residue = 2.*(pow(CUTOFF,-12.)-pow(CUTOFF,-6.)); 
    else if (WOON)
	residue = WOON_EPSILON/2 *
	    (WOON_Q*exp(-WOON_K1*(CUTOFF-WOON_RE))-
	     (WOON_Q+1)*exp(-WOON_K2*(CUTOFF-WOON_RE))); 
    
    /* initialize the properties cumulants */
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	{
	    /* interaction potential contribution to the stress */
	    virial[i][j] = 0.;
	    stress[i][j] = 0.;
	}
        
    for (i=0; i<6; i++)
	for (j=0; j<6; j++)
	    c[i][j] = 0.;
    
    pote = 0.;
    kine = 0.;
    for (i=0; i<3; i++)	cj[i] = 0;
    for (i=0; i<3*np; i++) f[i] = 0.;
    for (i=0; i<3*np; i+=3)
    {
	v[i]   = (H[0][0]*s1[i]+H[1][0]*s1[i+1]+H[2][0]*s1[i+2])/delta;
	v[i+1] = (H[0][1]*s1[i]+H[1][1]*s1[i+1]+H[2][1]*s1[i+2])/delta;
	v[i+2] = (H[0][2]*s1[i]+H[1][2]*s1[i+1]+H[2][2]*s1[i+2])/delta;
	cc = (v[i]*v[i]+v[i+1]*v[i+1]+v[i+2]*v[i+2])*mass[i/3]/2.;
	kine += cc;
	ep[i/3] = cc;
	for (k=0; k<3; k++)
	    for (l=0; l<3; l++)
		stress[k][l] += mass[i/3]*v[i+k]*v[i+l];
    }
    
    for (i=0; i<3*np; i+=3)
    {
	/* loops over atoms*/
	for (k=0; k<idx[i/3]; k++)
	{
	    /* inner loop: run through all the neigbors of i */
	    j = *(list[i/3]+k);
	    
	    sxij = s[i]   - s[j];
	    syij = s[i+1] - s[j+1];
	    szij = s[i+2] - s[j+2];

	    if (sxij>0.5) sxij--;
	    else if (sxij<-0.5) sxij++;
	    if (syij>0.5) syij--;
	    else if (syij<-0.5) syij++;
	    if (szij>0.5) szij--;
	    else if (szij<-0.5) szij++;

	    rx = H[0][0]*sxij + H[1][0]*syij + H[2][0]*szij;
	    ry = H[0][1]*sxij + H[1][1]*syij + H[2][1]*szij;
	    rz = H[0][2]*sxij + H[1][2]*syij + H[2][2]*szij;
	    r2 = rx*rx+ry*ry+rz*rz;

	    if (r2<CUTOFF*CUTOFF)
	    {
		if (LJ6_12)
		{
		    r2  = 1./r2;
		    r6  = r2*r2*r2;
		    r12 = r6*r6;
		    vij = 2.*(r12-r6)-residue;
		    ff  = (48.*r12-24.*r6)*r2;
		    cc = (672.*r12-192.*r6)*r2*r2;
		}
		else if (WOON)
		{
		    r = sqrt(r2);
		    a = WOON_EPSILON * WOON_Q * exp(-WOON_K1 * (r-WOON_RE));
		    b = WOON_EPSILON * (WOON_Q+1) * exp(-WOON_K2 * (r-WOON_RE)); 
		    vij = (a-b)/2.-residue;
		    ff = (WOON_K1*a-WOON_K2*b)/r;
		    cc = (WOON_K1*WOON_K1*a-WOON_K2*WOON_K2*b+ff)/r2;
		}
		
		ep[i/3] += vij;
		ep[j/3] += vij;
		pote  += vij+vij;
		
		ffx = ff*rx;
		ffy = ff*ry;
		ffz = ff*rz;
		f[i]   += ffx;
		f[i+1] += ffy;
		f[i+2] += ffz;
		f[j]   -= ffx;
		f[j+1] -= ffy;
		f[j+2] -= ffz;
		
		/* calculate the virial tensor */
		virial[0][0] += ffx*rx;
		virial[1][1] += ffy*ry;
		virial[2][2] += ffz*rz;
		virial[0][1] += (ffx*ry+ffy*rx)/2.;
		virial[0][2] += (ffx*rz+ffz*rx)/2.;
		virial[1][2] += (ffy*rz+ffz*ry)/2.;

		/* accumulate the heat current */
		tmp = (ffx*(v[i]+v[j])+
		       ffy*(v[i+1]+v[j+1])+
		       ffz*(v[i+2]+v[j+2]))/2.;
		cj[0] += tmp * rx;
		cj[1] += tmp * ry;
		cj[2] += tmp * rz;

		/* accumulate the elastic constants */
		/* 21 independent components */
		
		/* Born contribution */
		c[0][0] += cc*rx*rx*rx*rx;
		c[0][1] += cc*rx*rx*ry*ry;
		c[0][2] += cc*rx*rx*rz*rz;
		c[1][1] += cc*ry*ry*ry*ry;
		c[1][2] += cc*ry*ry*rz*rz;
		c[2][2] += cc*rz*rz*rz*rz;

		c[0][3] += cc*rx*rx*ry*rz;
		c[0][4] += cc*rx*rx*rx*rz;
		c[0][5] += cc*rx*rx*rx*ry;
		c[1][3] += cc*ry*ry*ry*rz;
		c[1][4] += cc*ry*ry*rx*rz;
		c[1][5] += cc*ry*ry*rx*ry;
		c[2][3] += cc*rz*rz*ry*rz;
		c[2][4] += cc*rz*rz*rx*rz;
		c[2][5] += cc*rz*rz*rx*ry;

		c[3][3] += cc*ry*rz*ry*rz;
		c[3][4] += cc*ry*rz*rx*rz;
		c[3][5] += cc*ry*rz*rx*ry;
		c[4][4] += cc*rx*rz*rx*rz;
		c[4][5] += cc*rx*rz*rx*ry;
		c[5][5] += cc*rx*ry*rx*ry;
	    }
	} /* end inner loop */
    } /* end outer loop */

    /*******************************************************/
    /* calculate the elastic constants without the virial  */
    /* fluctuation part, will be added in the main program */
    /*******************************************************/
    for (i=0; i<6; i++)
	for (j=i; j<6; j++)
	{
	    c[i][j] /= volume;
	    c[j][i] = c[i][j];
	}
    
    stress[0][0] = (stress[0][0]+virial[0][0])/volume;
    stress[1][1] = (stress[1][1]+virial[1][1])/volume;
    stress[2][2] = (stress[2][2]+virial[2][2])/volume;
    stress[0][1] = (stress[0][1]+virial[0][1])/volume;
    stress[0][2] = (stress[0][2]+virial[0][2])/volume;
    stress[1][2] = (stress[1][2]+virial[1][2])/volume;
    stress[1][0] = stress[0][1];
    stress[2][0] = stress[0][2];
    stress[2][1] = stress[1][2];
    
    virial[0][0] /= volume;
    virial[1][1] /= volume;
    virial[2][2] /= volume;
    virial[0][1] /= volume;
    virial[0][2] /= volume;
    virial[1][2] /= volume;
    virial[1][0] = virial[0][1];
    virial[2][0] = virial[0][2];
    virial[2][1] = virial[1][2];
    
    return;
} /* end pair-potential() */


void write_config (char filename[])
{
    /* save system configuration in readable ASCII format */
    int i; 
    FILE *config; 
    config = fopen(filename, "w+");
    fprintf (config, "Number of particles = %d\n\n", np);

    fprintf (config, "H(1,1) = %f Angstrom\n",  H[0][0]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(1,2) = %f Angstrom\n",  H[0][1]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(1,3) = %f Angstrom\n\n", H[0][2]*ULENGTH_IN_ANGSTROM);
    
    fprintf (config, "H(2,1) = %f Angstrom\n",  H[1][0]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(2,2) = %f Angstrom\n",  H[1][1]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(2,3) = %f Angstrom\n\n", H[1][2]*ULENGTH_IN_ANGSTROM);
    
    fprintf (config, "H(3,1) = %f Angstrom\n",  H[2][0]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(3,2) = %f Angstrom\n",  H[2][1]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(3,3) = %f Angstrom\n\n", H[2][2]*ULENGTH_IN_ANGSTROM);

    fprintf (config,
	     "Mass(amu)     sx        sy        sz     sx1(1/ps)  sy1(1/ps)  sz1(1/ps)\n");
    for (i=0;i<3*np;i+=3)
	fprintf (config, "%f  %f  %f  %f  %.6f  %.6f  %.6f\n",
		 mass[i/3]*UMASS_IN_KG*1000*AVO,
		 s[i], s[i+1], s[i+2], s1[i]/delta/UTIME_IN_PS,
		 s1[i+1]/delta/UTIME_IN_PS, s1[i+2]/delta/UTIME_IN_PS);
    fclose (config);
    return;
}  /* end write_config() */


void read_config (char filename[])
{
    /* load in system configuration from file */
    int i,j;
    FILE *config; 
    if ((config=fopen(filename, "r"))==NULL)
    {
      fprintf(stdout, "\n **Error: cannot locate file \"%s\"\n.", filename);
      exit(1);
    }

    fscanf (config, "Number of particles = %d\n\n", &np);

    fscanf (config, "H(1,1) = %lf Angstrom\n",  &H[0][0]);
    fscanf (config, "H(1,2) = %lf Angstrom\n",  &H[0][1]);
    fscanf (config, "H(1,3) = %lf Angstrom\n\n", &H[0][2]);
    
    fscanf (config, "H(2,1) = %lf Angstrom\n",  &H[1][0]);
    fscanf (config, "H(2,2) = %lf Angstrom\n",  &H[1][1]);
    fscanf (config, "H(2,3) = %lf Angstrom\n\n", &H[1][2]);
    
    fscanf (config, "H(3,1) = %lf Angstrom\n",  &H[2][0]);
    fscanf (config, "H(3,2) = %lf Angstrom\n",  &H[2][1]);
    fscanf (config, "H(3,3) = %lf Angstrom\n\n", &H[2][2]);
    
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	{
	    H[i][j] /= ULENGTH_IN_ANGSTROM; 
	    H1[i][j] = 0.;
	    H2[i][j] = 0.;
	    H3[i][j] = 0.;
	    H4[i][j] = 0.;
	    H5[i][j] = 0.;
	}
    volume = matinv(H, HI);
      
    fscanf (config, "Mass(amu)     sx        sy        sz     sx1(1/ps)  sy1(1/ps)  sz1(1/ps)\n");
    for (i=0;i<3*np;i+=3)
    {
	fscanf (config, "%lf  %lf  %lf  %lf  %lf  %lf  %lf\n",
		&mass[i/3], &s[i], &s[i+1], &s[i+2],
		 &s1[i],  &s1[i+1], &s1[i+2]);
	mass[i/3] /= (UMASS_IN_KG*1000*AVO);
	/* printf ("%d %f %f %f %f\n", i/3, mass[i/3], s[i], s[i+1], s[i+2]); */
    }
    fclose (config);

    for (i=0;i<3*np;i++)
    {
      s1[i] *= delta*UTIME_IN_PS;
      s2[i] = 0.;
      s3[i] = 0.;
      s4[i] = 0.;
      s5[i] = 0.;
    }
    return;
} /* end read_config() */

void smile(FILE *output)
{
    /* calculate the heat current correlation function averaged over */
    /* three orthogonal directions, and write results in OUTPUT */
    int i, j, k, l; 
    unsigned long size, n, ncorr, nfreq;
    double *data, *freq, *corr;
    double delta, total, first_dip;
    
    heat_current[0] = fopen (fn_heat_current[0], "r");
    fseek (heat_current[0], 0L, SEEK_END);
    size = ftell(heat_current[0]);
    
    /* can only use a chunk that is integer power of 2 */
    for (n=2; (size/sizeof(double)) > n; n<<=1); n>>=1;
    printf ("\nUsable heat current length = %ld,", n);
    if ((data=(double *)malloc(n*sizeof(double))) == NULL) 
    {
	fprintf (stdout, "\n **Error: not enough memory to load in %s\n.",
		 fn_heat_current[0]);
	exit(1);
    }

    /* the last double precision number is delta in ps */
    fseek (heat_current[0], size-sizeof(double), SEEK_SET);
    fread (&delta, sizeof(double), 1, heat_current[0]);
    fclose(heat_current[0]);
    printf ("  delta = %.4f ps\n\n", delta);
    
    ncorr = n/NPART;
    corr = (double *) malloc(ncorr*sizeof(double));
    for (i=0; i<ncorr; i++) corr[i] = 0.;
    nfreq = n/NPART;
    freq = (double *) malloc(nfreq*sizeof(double));
    for (i=0; i<nfreq; i++) freq[i] = 0.;
    /* conform to Fortran convention */
    data--; corr--; freq--;
    
    for (j=0; j<3; j++)
    {
	printf ("Loading in %s...\n", fn_heat_current[j]);
	heat_current[j] = fopen (fn_heat_current[j], "r");
	fseek (heat_current[j], size-(n+1)*sizeof(double), SEEK_SET);
	fread (data+1, sizeof(double), n, heat_current[j]);
	fclose (heat_current[j]);
	printf ("%s loaded in, FFT the heat current...\n", fn_heat_current[j]);

	realft (data, n, 1);
	printf ("FFT complete, getting the power spectrum...\n");
	    
	/* get rid of net drift, pure fluctuation */
	data[1] = 0.;
	data[2] = data[2]*data[2]/n;
	for (i=3; i<=n; i+=2)
	{
	    data[i] = (data[i]*data[i]+data[i+1]*data[i+1])/n;
	    data[i+1] = 0.;
	}
	
	/* divide by 2 to be the single sided integral */
	for (i=1; i<=nfreq; i++) freq[i] += data[2*i-1]/3./2.;

	printf ("inverse FFT the power spectrum...\n");
	realft (data, n, -1);
	printf ("inverse FFT complete, accumulate to counters.\n\n");
	for (i=1; i<=ncorr; i++) corr[i] += data[i]/3.;
	
    }  /* end processing principle directions */
    
    /* look for point where the correlation function first turns negative */
    for (i=1; (corr[i]>0.) && (i<=ncorr); i++);
    for (total=0.,l=1; l<=i; l++) total += corr[l];
    
    fprintf (output, "\n   From a piece of %.1f ps heat current data,\n",
	     n*delta);
    fprintf (output, "   the first dip in correlation function occurs at %.3f ps\n",
	     i*delta);
    fprintf (output, "   and by then the thermal conductivity is %.3f mW/m/K.\n\n",
	     1000*total);
    
    /* save the correlation function */
    corr_output = fopen (fn_corr_output, "w+");
    for (i=1; i<=ncorr; i++)
	fprintf (corr_output, "%.4f %E\n", (i-1)*delta, corr[i]);
    fclose (corr_output);
    
    /* save the power spectrum */
    freq_output = fopen (fn_freq_output, "w+");
    for (i=1; i<=nfreq; i++) fprintf (freq_output, "%d %f\n", i-1, freq[i]);
    fclose (freq_output);

    free(++data);
    free(++corr);
    free(++freq);
    
    return; 
} /* end smile() */

void shear_smile(FILE *output)
{
    /* calculate the shear stress correlation function averaged */
    /* over yz, xz, xy components, and write results in OUTPUT  */
    int i, j, k, l; 
    unsigned long size, n, ncorr, nfreq;
    double *data, *freq, *corr;
    double delta, total, first_dip;
    
    shear_stress[0] = fopen (fn_shear_stress[0], "r");
    fseek (shear_stress[0], 0L, SEEK_END);
    size = ftell(shear_stress[0]);
    
    /* can only use a chunk that is integer power of 2 */
    for (n=2; (size/sizeof(double)) > n; n<<=1); n>>=1;
    printf ("\nUsable shear stress length = %ld,", n);
    if ((data=(double *)malloc(n*sizeof(double))) == NULL) 
    {
	fprintf (stdout, "\n **Error: not enough memory to load in %s\n.",
		 fn_shear_stress[0]);
	exit(1);
    }

    /* the last double precision number is delta in ps */
    fseek (shear_stress[0], size-sizeof(double), SEEK_SET);
    fread (&delta, sizeof(double), 1, shear_stress[0]);
    fclose(shear_stress[0]);
    printf ("  delta = %.4f ps\n\n", delta);
    
    ncorr = n/NPART;
    corr = (double *) malloc(ncorr*sizeof(double));
    for (i=0; i<ncorr; i++) corr[i] = 0.;
    nfreq = n/NPART;
    freq = (double *) malloc(nfreq*sizeof(double));
    for (i=0; i<nfreq; i++) freq[i] = 0.;
    /* conform to Fortran convention */
    data--; corr--; freq--;
    
    for (j=0; j<3; j++)
    {
	printf ("Loading in %s...\n", fn_shear_stress[j]);
	shear_stress[j] = fopen (fn_shear_stress[j], "r");
	fseek (shear_stress[j], size-(n+1)*sizeof(double), SEEK_SET);
	fread (data+1, sizeof(double), n, shear_stress[j]);
	fclose (shear_stress[j]);
	printf ("%s loaded in, FFT the shear stress...\n", fn_shear_stress[j]);

	realft (data, n, 1);
	printf ("FFT complete, getting the power spectrum...\n");
	    
	/* get rid of net drift, pure fluctuation */
	data[1] = 0.;
	data[2] = data[2]*data[2]/n;
	for (i=3; i<=n; i+=2)
	{
	    data[i] = (data[i]*data[i]+data[i+1]*data[i+1])/n;
	    data[i+1] = 0.;
	}
	
	/* divide by 2 to be the single sided integral */
	for (i=1; i<=nfreq; i++) freq[i] += data[2*i-1]/3./2.;

	printf ("inverse FFT the power spectrum...\n");
	realft (data, n, -1);
	printf ("inverse FFT complete, accumulate to counters.\n\n");
	for (i=1; i<=ncorr; i++) corr[i] += data[i]/3.;
	
    }  /* end processing principle directions */
    
    /* look for point where the correlation function first turns negative */
    for (i=1; (corr[i]>0.) && (i<=ncorr); i++);
    for (total=0.,l=1; l<=i; l++) total += corr[l];
    
    fprintf (output, "\n   From a piece of %.1f ps shear stress data,\n",
	     n*delta);
    fprintf (output, "   the first dip in correlation function occurs at %.3f ps\n",
	     i*delta);
    fprintf (output, "   and by then the shear viscosity is %.3f uPa x second.\n\n",
	     1000000*total);
    
    /* save the correlation function */
    shear_corr_output = fopen (fn_shear_corr_output, "w+");
    for (i=1; i<=ncorr; i++)
	fprintf (shear_corr_output, "%.4f %E\n", (i-1)*delta, corr[i]);
    fclose (shear_corr_output);
    
    /* save the power spectrum */
    shear_freq_output = fopen (fn_shear_freq_output, "w+");
    for (i=1; i<=nfreq; i++) fprintf (shear_freq_output, "%d %f\n", i-1, freq[i]);
    fclose (shear_freq_output);

    free(++data);
    free(++corr);
    free(++freq);
    
    return; 
} /* end shear_smile() */

void realft (double data[], unsigned long n, int isign)
{
        void four1(double data[], unsigned long nn, int isign);
        unsigned long i,i1,i2,i3,i4,np3;
        double c1=0.5,c2,h1r,h1i,h2r,h2i;
        double wr,wi,wpr,wpi,wtemp,theta;

        theta=3.141592653589793/(n>>1);
        if (isign == 1) {
                c2 = -0.5;
                four1(data,n>>1,1);
        } else {
                c2=0.5;
                theta = -theta;
        }
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0+wpr;
        wi=wpi;
        np3=n+3;
        for (i=2;i<=(n>>2);i++) {
                i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
                h1r=c1*(data[i1]+data[i3]);
                h1i=c1*(data[i2]-data[i4]);
                h2r = -c2*(data[i2]+data[i4]);
                h2i=c2*(data[i1]-data[i3]);
                data[i1]=h1r+wr*h2r-wi*h2i;
                data[i2]=h1i+wr*h2i+wi*h2r;
                data[i3]=h1r-wr*h2r+wi*h2i;
                data[i4] = -h1i+wr*h2i+wi*h2r;
                wr=(wtemp=wr)*wpr-wi*wpi+wr;
                wi=wi*wpr+wtemp*wpi+wi;
        }
        if (isign == 1) {
                data[1] = (h1r=data[1])+data[2];
                data[2] = h1r-data[2];
        } else {
                data[1]=c1*((h1r=data[1])+data[2]);
                data[2]=c1*(h1r-data[2]);
                four1(data,n>>1,-1);
        }

	/* Added by Ju Li: factor of 2/n is */
	/* multiplied for inverse realft(); */
	if (isign == -1)
	    for (i=1; i<=n; i++)
		data[i] *= 2./n;
	
	return;
} /* end realft() */


#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void four1(double data[], unsigned long nn, int isign)
{
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        double tempr,tempi;

        n=nn << 1;
        j=1;
        for (i=1;i<n;i+=2) {
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        mmax=2;
        while (n > mmax) {
                istep=mmax << 1;
                theta=isign*(6.28318530717959/mmax);
                wtemp=sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=sin(theta);
                wr=1.0;
                wi=0.0;
                for (m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*data[j]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
} /* end four1() */
#undef SWAP

void assign_fcc_lattice()
{
    int i,j,k,l,ip;
    double x[12];
    
    if (4*nc[0]*nc[1]*nc[2] != np)
    {
	fprintf (stdout,
		 "\n **Error: number of particles not equal to %d, as defined by nc(1,2,3).\n",
		 4*nc[0]*nc[1]*nc[2]);
	exit(1);
    }
    
    x[0] = 0.;
    x[1] = 0.;
    x[2] = 0.;

    x[3] = 0.5;
    x[4] = 0.5;
    x[5] = 0.;

    x[6] = 0.;
    x[7] = 0.5;
    x[8] = 0.5;

    x[9] = 0.5;
    x[10] = 0.;
    x[11] = 0.5;
      
    ip = 0;
    for (i=0; i<nc[0]; i++)
	for (j=0; j<nc[1]; j++)
	    for (k=0; k<nc[2]; k++, ip+=12)
		for (l=0; l<12; l+=3)
		{
		    s[l+ip]   = (x[l]+i+1./8.)/nc[0];
		    s[l+1+ip] = (x[l+1]+j+1./8.)/nc[1];
		    s[l+2+ip] = (x[l+2]+k+1./8.)/nc[2];
		}
    
    return;
} /* end assign_fcc_lattice() */


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


void randnorm(double *x, int n)
{
    /* generate n normal-distributed random numbers in array x */
    int ip = 0;
    double s1, s2;
    while (ip < n) 
    {
	s1 = sqrt(-2.0*log(frandom()));
	s2 = 12*PI*frandom();
	x[ip++] = s1 * cos(s2);
	if (ip<n) x[ip++] = s1 * sin(s2);
    }
    return;
} /* end randnorm() */


double T_real_to_T_md (double TR, double *dTdR)
{
    /* convert real temperature to MD temperature: */
    /* as a first correction to QM effects */
    /* double cc = exp(DEBYE_TEMP/TR)-1; */
    /* *dTdR = exp(DEBYE_TEMP/TR)/cc/cc*DEBYE_TEMP*DEBYE_TEMP/TR/TR; */
    /* return ((1/cc+0.5)*DEBYE_TEMP); */
    *dTdR = 1.0;
    return (TR);
}

double T_md_to_T_real (double T)
{
    /* return (DEBYE_TEMP/log(1/(T/DEBYE_TEMP-0.5)+1.)); */
    return (T);
}

void matran(double A[3][3], double B[3][3])
{
    B[0][0] = A[0][0];
    B[0][1] = A[1][0];
    B[0][2] = A[2][0];
    B[1][0] = A[0][1];
    B[1][1] = A[1][1];
    B[1][2] = A[2][1];
    B[2][0] = A[0][2];
    B[2][1] = A[1][2];
    B[2][2] = A[2][2];
    return;
} /* end matran() */

void matadd(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0] + B[0][0];
    C[0][1] = A[0][1] + B[0][1];
    C[0][2] = A[0][2] + B[0][2];
    C[1][0] = A[1][0] + B[1][0];
    C[1][1] = A[1][1] + B[1][1];       
    C[1][2] = A[1][2] + B[1][2];       
    C[2][0] = A[2][0] + B[2][0];           
    C[2][1] = A[2][1] + B[2][1];      
    C[2][2] = A[2][2] + B[2][2];
    return;
} /* end matadd() */

void matsub(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0] - B[0][0];
    C[0][1] = A[0][1] - B[0][1];
    C[0][2] = A[0][2] - B[0][2]; 
    C[1][0] = A[1][0] - B[1][0];         
    C[1][1] = A[1][1] - B[1][1];       
    C[1][2] = A[1][2] - B[1][2];       
    C[2][0] = A[2][0] - B[2][0];           
    C[2][1] = A[2][1] - B[2][1];      
    C[2][2] = A[2][2] - B[2][2];
    return;
} /* end matsub() */

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() */

void matmul(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0]*B[0][0] + A[0][1]*B[1][0] + A[0][2]*B[2][0];
    C[0][1] = A[0][0]*B[0][1] + A[0][1]*B[1][1] + A[0][2]*B[2][1];
    C[0][2] = A[0][0]*B[0][2] + A[0][1]*B[1][2] + A[0][2]*B[2][2];
    C[1][0] = A[1][0]*B[0][0] + A[1][1]*B[1][0] + A[1][2]*B[2][0];
    C[1][1] = A[1][0]*B[0][1] + A[1][1]*B[1][1] + A[1][2]*B[2][1];
    C[1][2] = A[1][0]*B[0][2] + A[1][1]*B[1][2] + A[1][2]*B[2][2];
    C[2][0] = A[2][0]*B[0][0] + A[2][1]*B[1][0] + A[2][2]*B[2][0];
    C[2][1] = A[2][0]*B[0][1] + A[2][1]*B[1][1] + A[2][2]*B[2][1];
    C[2][2] = A[2][0]*B[0][2] + A[2][1]*B[1][2] + A[2][2]*B[2][2];
    return;
} /* end matmul() */