/******************************************/
/* Molecular Dynamics Simulation of Argon */
/* using Pair Potential.                  */
/*                                        */
/*               Version 2.1 Feb.23, 1999 */
/*                            Ju Li (MIT) */
/******************************************/

/*
  Version Notes:
  1.0  Q/A interface;
       Woon's pair potential;
       NTH, NTS ensembles;
       Thermal conductivity;
  1.1  Micro-canonical ensemble
       fluctuation formulas;
       Shear viscosity;
  2.0  O(N) implementation:
       Bin mesh and list;
       Particle anchor list;
       X Window visualization.
  2.1  Temperature, stress or
       unit cell H schedules;
       Volumetric constraints;
       Immobile atoms;
       Binary interaction and
       arbitrary masses compound;
       NTS runaway protection;
*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include "utils.h"
/* spectral method */
#include "smile.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 UVOLUME_IN_M3 pow(ULENGTH_IN_M,3.)
#define USTRESS_IN_PA (UENERGY_IN_J/UVOLUME_IN_M3)
#define USTRESS_IN_MPA (USTRESS_IN_PA/1.E6)
#define USTRESS_IN_BAR (USTRESS_IN_PA/1.E5)
#define UCONDUCTIVITY_IN_WMK (UENERGY_IN_J/UTIME_IN_SECOND/ULENGTH_IN_M)
#define DEFAULT_UNIT_BASE (6)
/* 5th-order predictor-corrector algorithm constants */
#define F02 (3./16.)
#define F12 (251./360.)
#define F32 (11./18.)
#define F42 (1./6.)
#define F52 (1./60.)
/* interaction cutoff (binary interactions will be scaled to it) */
#define RCUT (2.5*SIGMA_IN_M/ULENGTH_IN_M)
/* max(fabs(Lagrangian strain eigenvalues)) as re-mesh trigger */
#define DEFAULT_MAX_STRAIN 0.02
/* maximum particle drift from anchoring s[] in pure Ar */
#define DEFAULT_MAX_DRIFT 0.08
/* g(r) 1/sqrt(N) statistical accuracy */
#define GRACCU 0.001
#define GRMESH 200
#define GRMAXBINS 125
#define DEFAULT_GRCUT_IN_A (15.)
#define GRFILE "gr.m"
/* coefficient of 1/sqrt(N) fluctuation tolerance */
#define MAX_FLUCTUATION (4.)
#define MAX_TAG 256
#define TINY (1.E-8)
#define MAX_STRLEN 128
#define DEFAULT_TMP_CONFIG "tmp."
#define DEFAULT_TMP_CONFIG_FREQ 10 /* multiplied by output freq. */
#define TIME_TO_NOTICE (0.)
#define DEFAULT_ATOM_TP "Ar"
#define TP_STRLEN 2
#define JUMPSTART_KINE (0.1) /* in Kelvin */
#define ELASTIC_CONSTANT_IN_MPA (2000.)
/* maximum number of different kinds of atoms */
#define MAX_KIND 128 
/* just a PSEUDONYM to make Rasmol render nice pictures */
#define DEFAULT_NAR_TP "He"
#define DEFAULT_NAR_NUMBER rint(np*0) /* number of NAR's */
#define DEFAULT_NAR_MASS_RATIO (9999.0) /* mass ratio */
#define DEFAULT_NN_L_RATIO (1.0) /* interaction length ratio */
#define DEFAULT_NN_E_RATIO (1.0) /* interaction energy ratio */
#define DEFAULT_AN_L_RATIO sqrt(DEFAULT_NN_L_RATIO)
#define DEFAULT_AN_E_RATIO ((DEFAULT_NN_E_RATIO+1)/2.)
/* mass greater than this is a token for immobile atom */
#define IMMOBILE_MASS (9999.)
#define MAX_VOIDS 128
#define NTS_MAX_STRAIN2 (3*square(0.5))
#define DEFAULT_MESH_FILE "me.out"
#define MESH 100
/* cubic ijk -> linear bin index */
#define nbin(i,j,k) (((i+(maxbin[0]<<8))%maxbin[0]*maxbin[1]+ \
		      (j+(maxbin[1]<<8))%maxbin[1])*maxbin[2]+ \
		     (k+(maxbin[2]<<8))%maxbin[2])
/* i as pseudo-randomized "owner" of any i\neq j pair interaction */
#define own_it(i,j) ((((i+j)%2)&&(i>j))||((!((i+j)%2))&&(i<j)))
#define min(x,y) ((x)<(y)?(x):(y))
#define max(x,y) ((x)>(y)?(x):(y))
/* symmetric 3 x 3 -> 0-5 */
#define od(i,j) ((5-min(i,j))*min(i,j)/2+max(i,j))
/* binary interaction potential length and energy scaler */
#define L(i,j) (l_ratio[tpd[i]+tpd[j]])
#define E(i,j) (e_ratio[tpd[i]+tpd[j]])
#define square(i) ((i)*(i))

static int LJ6_12, WOON, CALC, NTH, NTS, HSCH, SSCH, TSCH, VOLUMETRIC,
    HCALM, SCALM, TCALM, USE_CONFIG_H, maxkb, kstart, kw_lp, iseed, np,
    np3, neff, kw_gr, nw_repartition, nw_anchor, anchortum, kb, nkind, 
    ntpd, nk, nar, scale_H, scale_T, Havg_start, nmb, *mb, *tpd, *mybin, 
    *binbindx, *binbinlist, *bindx, *binlist, *idx, *list, nc[4], 
    maxbin[4], kind[MAX_KIND][2];
static double MAX_STRAIN, MAX_DRIFT, RLIST, GRCUT, TR, delta, tau,
    wallmass, volume, volumeold, hydro, shear, T, dTdR, totalmass,
    r2del, kinewall, HTnow, msd, msd_kstart, structure_factor, pote,
    kine, Tnow, dE2, pressure, hamiltonian, factorT, factorHT, Treal,
    dV2, dEdV, bulk_modulus, alpha, Cv, Cp, *sa, *s, *s1, *s2, *s3,
    *s4, *s5, *v, *f, *fs, *d, *mass, *ep, l_ratio[3], e_ratio[3],
    maxl_ratio, minl_ratio[2], r[4], x[12], cj[3], c[6][6], H[3][3],
    H1[3][3], H2[3][3], H3[3][3], H4[3][3], H5[3][3], HI[3][3],
    HIOLD[3][3], A[3][3], B[3][3], G[3][3], GI[3][3], GA[3][3],
    G1[3][3], HT[3][3], H1T[3][3], HIT[3][3], Havg[3][3], sout[3][3],
    virial[3][3], stress[3][3], fH[3][3], HIBEGIN[3][3];
static int nvoids, tmp_config_freq;
static double voids[MAX_VOIDS][4], df[4], af[5], density[MESH], vx[MESH],
    temp[MESH], taum[MESH][3][3];
static char tmp_config[MAX_STRLEN];
/* current and long-term averages */
static double kinetum, k1tum, k2tum, Htum[3][3], strestum[3][3], 
    v1k1tum[3][3], v1k2tum[3][3], v2k2tum[3][3][3][3], ck1tum[6][6];
static double kinesum, k1sum, k2sum, potesum, Hsum[3][3], stressum[3][3],
    v1k1sum[3][3], v1k2sum[3][3], v2k2sum[3][3][3][3], ck1sum[6][6];
static char buf[512], fn_config_read[MAX_STRLEN], **tp,
    fn_config_write[MAX_STRLEN], potfun_name[MAX_STRLEN],
    fn_property[] = "temp.out", fn_structure[] = "struc.out",
    *Hschname[6] = {"H11","H12","H13","H22","H23","H33"},
    *Sschname[6] = {"S11","S12","S13","S22","S23","S33"};
static FILE *config,*property,*structure,*heat_current[3],*shear_stress[3];

static void write_config (char filename[]);
static void save_tmp_config ();
static void save_mesh (char filename[]);
static void ds (double s1[], double s2[], double d[]);
static double r2 (int i, int j, double si[], double sj[]);
static void pair_potential();
static double T_real_to_T_md (double TR, double *dTdR);
static double T_md_to_T_real (double T);
static void divide_supercell_into_bins();
static void particles_move_in();
static void transfer_particle (int i, int o, int n);
static void re_anchor (int i);
static void initialize_all_lists();
static void print_list_statistics (FILE *out);
static void accumulate_and_save_gr ();
static void print_atom_statistics (FILE *out, char *sp);
static int TFE (double sz1, double sz2);

/* save system configuration in readable text format */
void write_config (char filename[])
{
    int i;
    FILE *config; 
    config = fopen (filename, "w+");
    fprintf (config, "Number of particles = %d\n\n", np);
    fprintf (config, "H(1,1) = %f A\n", H[0][0]*ULENGTH_IN_A);
    fprintf (config, "H(1,2) = %f A\n", H[0][1]*ULENGTH_IN_A);
    fprintf (config, "H(1,3) = %f A\n", H[0][2]*ULENGTH_IN_A);
    fprintf (config, "nc(1)  = %d\n\n", nc[0]);
    fprintf (config, "H(2,1) = %f A\n", H[1][0]*ULENGTH_IN_A);
    fprintf (config, "H(2,2) = %f A\n", H[1][1]*ULENGTH_IN_A);
    fprintf (config, "H(2,3) = %f A\n", H[1][2]*ULENGTH_IN_A);
    fprintf (config, "nc(2)  = %d\n\n", nc[1]);
    fprintf (config, "H(3,1) = %f A\n", H[2][0]*ULENGTH_IN_A);
    fprintf (config, "H(3,2) = %f A\n", H[2][1]*ULENGTH_IN_A);
    fprintf (config, "H(3,3) = %f A\n", H[2][2]*ULENGTH_IN_A);
    fprintf (config, "nc(3)  = %d\n\n", nc[2]);
    fprintf (config, "TP Mass(amu)     sx        sy        sz ");
    fprintf (config, "    sx1(1/ns)  sy1(1/ns)  sz1(1/ns)\n");
    for (i=0; i<3*np; i+=3)
	fprintf (config, "%s %f  %f  %f  %f  %.6f  %.6f  %.6f\n",
		 tp[i/3], mass[i/3]*UMASS_IN_KG*1000*AVO,
		 s[i], s[i+1], s[i+2],
		 s1[i]  /delta/UTIME_IN_PS*1000,
		 s1[i+1]/delta/UTIME_IN_PS*1000,
		 s1[i+2]/delta/UTIME_IN_PS*1000);
    fclose (config);
    return;
}  /* end write_config() */


void save_tmp_config ()
{
    char buf[MAX_STRLEN];
    static int tmp_config_counter = 1;
    if (strcasecmp(tmp_config, "NULL"))
    {
	sprintf (buf, "%s%d", tmp_config, tmp_config_counter);
	write_config (buf);
	printf ("saved on \"%s\" at step %d, t = %.2f ps.\n\n",
		buf, kb, kb*delta*UTIME_IN_PS);
	tmp_config_counter++;
    }
    return;
} /* end save_tmp_config() */


/* meshed properties output */
void save_mesh (char filename[])
{
    int i;
    FILE *out;
    double tot = volume / MESH * (kb-kstart+1);
    if ( strcasecmp(filename, "NULL") )
    {
	out = fopen (filename, "w");
	fprintf (out, "%d %f %d %d %f %f\n", kb, kb*delta*UTIME_IN_PS,
		 kstart, kb-kstart+1, H[2][2]*nc[3]/nc[2], H[2][2]*df[2]);
	for (i=0; i<MESH; i++)
	    fprintf( out, "%f %f %f %f %f %f\n",
		     (i+0.5)/MESH*H[2][2],
		     density[i] / tot, 
		     density[i]>0? vx[i] / density[i] : 0,
		     density[i]>0? UENERGY_IN_J / BOLZ *
		     (temp[i]/density[i]-ARGON_MASS_IN_KG/
		      UMASS_IN_KG*square(vx[i]/density[i]))/3. : 0,
		     density[i]>0? taum[i][0][2]/tot : 0,
		     density[i]>0? (taum[i][0][0]+taum[i][1][1]+
				    taum[i][2][2])/3./tot : 0 );
    }
    fclose(out);
    return;
} /* end save_mesh() */


/* ds() and r2() both take first argument as origin */
void ds (double s1[], double s2[], double d[])
{
    d[0] = s2[0] - s1[0];
    d[1] = s2[1] - s1[1];
    d[2] = s2[2] - s1[2];
    /* unique j image in H-parallelepiped */
    while (d[0] > 0.5)  d[0]--;
    while (d[0] < -0.5) d[0]++; 
    while (d[1] > 0.5)  d[1]--;
    while (d[1] < -0.5) d[1]++;
    while (d[2] > 0.5)  d[2]--;
    while (d[2] < -0.5) d[2]++; 
    return;
} /* end ds() */


/* (r[0],r[1],r[2]) is x_{ji}, r[3] is |x_{ji}|^2 */
/* r[] should be a cached global for efficiency.  */
double r2 (int i, int j, double si[], double sj[])
{
    double ds1, ds2, ds3;
    ds1 = sj[3*j]   - si[3*i];
    ds2 = sj[3*j+1] - si[3*i+1];
    ds3 = sj[3*j+2] - si[3*i+2];
    if (ds1 > 0.5) ds1--;
    else if (ds1 < -0.5) ds1++; 
    if (ds2 > 0.5) ds2--;
    else if (ds2 < -0.5) ds2++; 
    if (ds3 > 0.5) ds3--;
    else if (ds3 < -0.5) ds3++;
    /* centered around i */
    r[0] = ds1*H[0][0] + ds2*H[1][0] + ds3*H[2][0];
    r[1] = ds1*H[0][1] + ds2*H[1][1] + ds3*H[2][1];
    r[2] = ds1*H[0][2] + ds2*H[1][2] + ds3*H[2][2];
    return (r[3] = r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
} /* end r2() */


/* 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)
/* Generic pair-wise interatomic potential */
void pair_potential()
{
    int i, j, k, l;
    double r1, ir2, ir6, ir12, a, b, vij;
    double ff, ffx, ffy, ffz, cc, dd;
    double res = LJ6_12?LJ_RESIDUE:WOON_RESIDUE;
    double ref = LJ6_12?LJ_RESIDUE_FORCE:WOON_RESIDUE_FORCE;
    
    /* initialize property cumulants */
    pote = kine = cj[0] = cj[1] = cj[2] = 0;
    for (i=0; i<np3; i++) f[i] = 0;
    /* virial[][]: potential contribution to stress */
    virial[0][0] = stress[0][0] = 0;
    virial[1][1] = stress[1][1] = 0;
    virial[2][2] = stress[2][2] = 0;
    virial[0][1] = stress[0][1] = 0;
    virial[0][2] = stress[0][2] = 0;
    virial[1][2] = stress[1][2] = 0;
    /* elastic constants */
    for (i=0; i<6; i++)
	for (j=i; j<6; j++) c[i][j] = 0.; 
    
    /* real-space velocities and some particle properties */
    for (i=0; i<np3; i+=3)
    {
	v[i]   = ( s1[i]*H[0][0] + s1[i+1]*H[1][0] +
		   s1[i+2]*H[2][0] ) / delta;
	v[i+1] = ( s1[i]*H[0][1] + s1[i+1]*H[1][1] +
		   s1[i+2]*H[2][1] ) / delta;
	v[i+2] = ( s1[i]*H[0][2] + s1[i+1]*H[1][2] +
		   s1[i+2]*H[2][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=k; l<3; l++)
		stress[k][l] += mass[i/3] * v[i+k] * v[i+l];
	if ( mb[i/3] && (!tpd[i/3]) && (kb>=kstart) )
	{ /* begin to accumulate meshed properties */
	    j = floor(s[i+2]*MESH);
	    while (j > MESH ) j -= MESH;
	    while (j < 0)     j += MESH;
	    density[j] += 1;    /* total occurrences */
	    vx[j] += v[i];      /* "total" x-speed   */
	    temp[j] += 2 * cc;  /* "total" m*v*v     */
	    for (k=0; k<3; k++)
		for (l=k; l<3; l++)
		    taum[j][k][l] += mass[i/3] * v[i+k] * v[i+l];
	}
    }
    
    for (i=0; i<np; i++) /* expects neighbor-list support */
	for (k=idx[2*i]; k<idx[2*i+1]; k++)
	    if (r2(i,j=list[k],s,s) < square(L(i,list[k])*RCUT))
	    {
		/* scale to RCUT frame */
		r1 = sqrt(r[3]) / L(i,j);
		if ( LJ6_12 )
		{
		    ir2 = 1. / r1 / r1;
		    ir6 = ir2 * ir2 * ir2;
		    ir12 = ir6 * ir6;
		    vij = 4*(ir12-ir6) - res + ref*(r1-RCUT);
		    ff = (48*ir12-24*ir6)*ir2 - ref/r1;
		    cc = ((624*ir12-168*ir6)*ir2 + ff) * ir2;
		}
		else
		{
		    a = WOON_EPSILON * WOON_Q *
			exp(-WOON_K1 * (r1 - WOON_RE));
		    b = WOON_EPSILON * (WOON_Q+1) *
			exp(-WOON_K2 * (r1 - WOON_RE));
 		    vij = a - b - res + ref*(r1-RCUT);
		    ff = (WOON_K1*a - WOON_K2*b - ref) / r1;
		    cc = (WOON_K1*WOON_K1*a -
			  WOON_K2*WOON_K2*b + ff) * ir2;
		}
		
		vij *= E(i,j);
		ff *= E(i,j) / L(i,j) / L(i,j);
		cc *= E(i,j) / L(i,j) / L(i,j) / L(i,j) / L(i,j);
		    
		pote  += vij;
		ep[i] += vij / 2.;
		ep[j] += vij / 2.;
		
		ffx = ff * r[0];
		ffy = ff * r[1];
		ffz = ff * r[2];
		
		f[3*j]   += ffx;
		f[3*j+1] += ffy;
		f[3*j+2] += ffz;
		f[3*i]   -= ffx;
		f[3*i+1] -= ffy;
		f[3*i+2] -= ffz;
		
		/* calculate the Virial tensor */
		virial[0][0] += ffx * r[0];
		virial[1][1] += ffy * r[1];
		virial[2][2] += ffz * r[2];
		virial[0][1] += ffx * r[1];
		virial[0][2] += ffx * r[2];
		virial[1][2] += ffy * r[2];

		/* meshed internal fluid stress tensor */
		dd = s[3*j+2] - s[3*i+2];
		while (dd < -0.5) dd++;
		while (dd >  0.5) dd--;
		l = floor((s[3*i+2]+dd/2)*MESH);
		while (l > MESH ) l -= MESH;
		while (l < 0)     l += MESH;
		if ( mb[i] && (!tpd[i]) &&
		     mb[j] && (!tpd[j]) && (kb>=kstart) )
		{
		    taum[l][0][0] += ffx * r[0];
		    taum[l][1][1] += ffy * r[1];
		    taum[l][2][2] += ffz * r[2];
		    taum[l][0][1] += ffx * r[1];
		    taum[l][0][2] += ffx * r[2];
		    taum[l][1][2] += ffy * r[2];
		}
						
		/* accumulate the heat currents */
		dd = ( ffx * (v[3*i]   + v[3*j]  ) +
		       ffy * (v[3*i+1] + v[3*j+1]) +
		       ffz * (v[3*i+2] + v[3*j+2]) ) / 2.;
		cj[0] += dd * r[0];
		cj[1] += dd * r[1];
		cj[2] += dd * r[2];
		
		/* accumulate the elastic constants */
		/* 21 independent components */
		
		/* Born contribution */
		c[0][0] += cc*r[0]*r[0]*r[0]*r[0];
		c[0][1] += cc*r[0]*r[0]*r[1]*r[1];
		c[0][2] += cc*r[0]*r[0]*r[2]*r[2];
		c[1][1] += cc*r[1]*r[1]*r[1]*r[1];
		c[1][2] += cc*r[1]*r[1]*r[2]*r[2];
		c[2][2] += cc*r[2]*r[2]*r[2]*r[2];

		c[0][3] += cc*r[0]*r[0]*r[1]*r[2];
		c[0][4] += cc*r[0]*r[0]*r[0]*r[2];
		c[0][5] += cc*r[0]*r[0]*r[0]*r[1];
		c[1][3] += cc*r[1]*r[1]*r[1]*r[2];
		c[1][4] += cc*r[1]*r[1]*r[0]*r[2];
		c[1][5] += cc*r[1]*r[1]*r[0]*r[1];
		c[2][3] += cc*r[2]*r[2]*r[1]*r[2];
		c[2][4] += cc*r[2]*r[2]*r[0]*r[2];
		c[2][5] += cc*r[2]*r[2]*r[0]*r[1];
		
		c[3][3] += cc*r[1]*r[2]*r[1]*r[2];
		c[3][4] += cc*r[1]*r[2]*r[0]*r[2];
		c[3][5] += cc*r[1]*r[2]*r[0]*r[1];
		c[4][4] += cc*r[0]*r[2]*r[0]*r[2];
		c[4][5] += cc*r[0]*r[2]*r[0]*r[1];
		c[5][5] += cc*r[0]*r[1]*r[0]*r[1];
	    }

    /* micro-convection part of the heat current */
    for (i=0; i<np3; i++) cj[i%3] += ep[i/3] * v[i];
    
    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];
    
    /* evaluate 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];
	}
    return;
} /* end pair_potential() */


/* real-world QM temperature to MD simulation temperature */
double T_real_to_T_md (double TR, double *dTdR)
{
    *dTdR = 1.0;
    return (TR);
} /* end T_real_to_T_md() */


/* inverse transform */
double T_md_to_T_real (double T)
{
    return (T);
} /* end T_md_to_T_real() */


/* Divide the supercell parallelepiped into bins according */
/* to RLIST and make linklist in binbindx[], binbinlist[]. */
void divide_supercell_into_bins()
{
    static int first_entry = 1;
    int i, j, k, di, dj, dk, l, m, n;
    double dd[3];
    if (first_entry) first_entry = 0;
    else
    {  /* free prior allocations */
	free (binbindx);
	free (binbinlist);
    }
    /* this is the starting and reference point */
    volume = volumeold = matinv (H, HIOLD);
    /* state of the ensuing strain session */
    for (i=0; i<3; i++)
    {  /* current thickness of H in the ith direction */
	dd[i] = 1./sqrt(HIOLD[0][i]*HIOLD[0][i]+
			HIOLD[1][i]*HIOLD[1][i]+
			HIOLD[2][i]*HIOLD[2][i]);
	/* multiplied by maximum shrink factor during the session */
	dd[i] *= sqrt(1-2*MAX_STRAIN);
	if (dd[i] <= 2*maxl_ratio*RCUT)
	{
	    printf ("divide_supercell_into_bins: possible double "
		    "interaction images:\n");
	    printf ("current thickness[%d] * sqrt(1-2*MAX_STRAIN)\n"
		    "= %.5f <= 2*maxl_ratio*RCUT = %.5f.\n", i, dd[i],
		    2*maxl_ratio*RCUT);
	    exit(1);
	}
	/* necessary thickness of the bins right now for */
	/* its thickness to be ALWAYS > maxl_ratio*RLIST */
	maxbin[i] = floor(dd[i] / RLIST / maxl_ratio);
	if (maxbin[i] < 1) maxbin[i] = 1;
    }
    maxbin[3] = maxbin[0] * maxbin[1] * maxbin[2];
    /* maximal (full) connectivity in 3D */
    binbinlist = (int *) malloc(maxbin[3]*27*sizeof(int));
    binbindx = (int *) malloc((2*maxbin[3]+1)*sizeof(int));
    binbindx[2*maxbin[3]] = maxbin[3]*27;
    /* actual linklist for bins that can interact */
    for (i=0; i<maxbin[0]; i++)
	for (j=0; j<maxbin[1]; j++)
	    for (k=0; k<maxbin[2]; k++)
	    {
		l = nbin(i,j,k);
		binbindx[2*l] = 27 * l;
		binbindx[2*l+1] = binbindx[2*l];
		for (di=-1; di<=1; di++)
		    for (dj=-1; dj<=1; dj++)
			for (dk=-1; dk<=1; dk++)
			{
			    m = nbin(i+di,j+dj,k+dk);
			    /* make sure it is a new record */
			    for (n=binbindx[2*l]; n<binbindx[2*l+1]; n++)
				if (m==binbinlist[n]) break;
			    if (n==binbindx[2*l+1])
				binbinlist[binbindx[2*l+1]++] = m;
			}
	    }
    return;
} /* end divide_supercell_into_bins() */


/* Put particles in bins and construct bin-particle list */
/* bindx[],binlist[] and partic-partic list idx[],list[] */
void particles_move_in ()
{
    static int first_entry = 1;
    int i, j, l;
    double maxp;
    if (first_entry) first_entry = 0;
    else 
    {  /* free prior allocations */
	free (bindx); free (binlist);
	free (idx); free (list);
    }
    /* "most" possible number of atoms appearing in one bin */
    maxp = (double) np / maxbin[3];
    maxp = ceil(maxp * (1+MAX_FLUCTUATION/sqrt(maxp)));
    binlist = (int *) malloc(maxbin[3]*maxp*sizeof(int));
    bindx = (int *) malloc((2*maxbin[3]+1)*sizeof(int));
    bindx[2*maxbin[3]] = maxbin[3]*maxp;
    if ( (binlist == NULL) || (bindx == NULL) )
    {
	printf ("particles_move_in: not enough memory "
		"(%d bytes) for binlist.\n",
		(int)(maxbin[3]*maxp*sizeof(int)));
	exit(1);
    }
    for (i=0; i<maxbin[3]; i++)
    {  /* uniform allocation of initial space */
	bindx[2*i]   = i * maxp;
	bindx[2*i+1] = bindx[2*i];
    }
    /* max. particle remembers nearby anchors in sphere maxl_ratio*RLIST */
    l = ceil(pow(maxl_ratio*RLIST/sqrt(1-2*MAX_STRAIN),3.) *
	     np/volume*4*PI/3/2);
    l = ceil(l * (1+MAX_FLUCTUATION/sqrt((double)l)));
    list = (int *) malloc(np*l*sizeof(int));
    idx = (int *) malloc((2*np+1)*sizeof(int));
    idx[2*np] = np*l;
    if ( (list == NULL) || (idx == NULL) )
    {
	printf ("particles_move_in: not enough memory "
		"(%d bytes) for particle list.\n",
		(int)(np*l*sizeof(int)));
	exit(1);
    }
    for (i=0; i<np; i++)
    {
	idx[2*i]   = i * l;
	idx[2*i+1] = idx[2*i];
    }
    for (i=0; i<np; i++)
    {  /* fill in particles */
	while (s[3*i] < 0)   s[3*i]++;
	while (s[3*i] > 1)   s[3*i]--;
	while (s[3*i+1] < 0) s[3*i+1]++;
	while (s[3*i+1] > 1) s[3*i+1]--;
	while (s[3*i+2] < 0) s[3*i+2]++;
	while (s[3*i+2] > 1) s[3*i+2]--;
	/* the bin that the current atom (anchor) belongs to */
	l = mybin[i] = nbin((int)floor(s[3*i]   * maxbin[0]),
			    (int)floor(s[3*i+1] * maxbin[1]),
			    (int)floor(s[3*i+2] * maxbin[2]));
	safe_append (bindx, binlist, i, l, 0, maxbin[3]);
	/* only privileged functions can modify sa[], */
	/* re_anchor() piecemeal, here wholesale. */
	sa[3*i]   = s[3*i];
	sa[3*i+1] = s[3*i+1];
	sa[3*i+2] = s[3*i+2]; 
    }
    /* generate particle-particle (anchor) list */
    for (i=0; i<np; i++) re_anchor(i);
    return;
}  /* particles_move_in() */


/* Transfer particle from one bin to another */
void transfer_particle (int i, int o, int n)
{
    int m;
    /* delete this particle from bin o */
    for (m=bindx[2*o];
	 (m<bindx[2*o+1]) && (binlist[m]!=i);
	 m++);
    if (m == bindx[2*o+1])
    {
	printf ("transfer_particle: deleting particle %d \n", i);
	printf ("which does not exist in bin No.%d.\n", o);
	exit(1);
    }
    for (m++; m<bindx[2*o+1]; m++)
	binlist[m-1] = binlist[m];
    bindx[2*o+1]--;
    /* insert this particle at the end of the new bin */
    safe_append (bindx, binlist, i, n, 0, maxbin[3]);
    mybin[i] = n;
    return;
} /* end transfer_particle() */


/* particle refreshes its own and new neighbors' anchor lists */
void re_anchor (int i)
{
    int j, k, l, m, n, tag[MAX_TAG], tag_top=0;
    /* set anchors */
    sa[3*i]   = s[3*i];
    sa[3*i+1] = s[3*i+1];
    sa[3*i+2] = s[3*i+2];
    /* the bin that this atom (anchor) belongs to */
    m = nbin((int)floor(sa[3*i]   * maxbin[0]),
	     (int)floor(sa[3*i+1] * maxbin[1]),
	     (int)floor(sa[3*i+2] * maxbin[2]));
    if (m != mybin[i]) transfer_particle (i, mybin[i], m);
    for (j=binbindx[2*mybin[i]]; j<binbindx[2*mybin[i]+1]; j++)
    { /* neighboring bins */
	l = binbinlist[j];
	for (k=bindx[2*l]; k<bindx[2*l+1]; k++)
	{ /* particles in neighboring bins */
	    m = binlist[k];
	    if ( (i != m) &&
		 (r2(i,m,sa,sa) <= square(L(i,m)*RLIST)) )
	    { /* this pair interaction "belongs" to either i or m */
		if ( own_it(i,m) )
		    if (tag_top >= MAX_TAG)
		    { /* overflow detection */
			printf ("re_anchor: MAX_TAG = %d exceeded.\n",
				MAX_TAG);
			exit(1);
		    }
		    else tag[tag_top++] = m;
		else
		{  /* see if i is already in m's list */
		    for (n=idx[2*m];
			 (n<idx[2*m+1]) && (list[n]!=i); n++);
		    if (n == idx[2*m+1]) /* if not */
			safe_append (idx, list, i, m, 0, np);
		}
	    }
	}
    }
    if (idx[2*i]+tag_top <= idx[2*i+2])
    { /* there is plenty of space */
	memcpy(list+idx[2*i], tag, tag_top*sizeof(int));
	idx[2*i+1] = idx[2*i] + tag_top;
    }
    else
    {
	memcpy(list+idx[2*i], tag, (idx[2*i+2]-idx[2*i])*sizeof(int));
	idx[2*i+1] = idx[2*i+2];
	for (n=idx[2*i+2]-idx[2*i]; n<tag_top; n++)
	    safe_append (idx, list, i, tag[n], 0, np);
    }
    return;
} /* end re_anchor() */


/* Bin-bin, bin-particle, particle-bin lists, */
/* partic-partic list; HIOLD[][], sa[] anchors */
void initialize_all_lists()
{
    int i; double dd[3];
    divide_supercell_into_bins();
    /* current bin thicknesses on three faces */
    for (i=0; i<3; i++)
	dd[i] = 1./sqrt(HIOLD[0][i]*HIOLD[0][i]+
			HIOLD[1][i]*HIOLD[1][i]+
			HIOLD[2][i]*HIOLD[2][i]) / maxbin[i];
    printf ("The current (step = %d)\n", kb);
    printf ("    | %f %f %f |\n", H[0][0], H[0][1], H[0][2]);
    printf ("H = | %f %f %f | => reference state\n",
	    H[1][0], H[1][1], H[1][2]);
    printf ("    | %f %f %f |\n", H[2][0], H[2][1], H[2][2]);
    printf ("is divided into %d x %d x %d = %d bins, of thicknesses\n",
	    maxbin[0], maxbin[1], maxbin[2], maxbin[3]);
    printf ("(%.3f %.3f %.3f) * sqrt(1-2*%.3f) >= RLIST = %.3f.\n\n",
	    dd[0], dd[1], dd[2], MAX_STRAIN, RLIST);
    particles_move_in();
    return;
}  /* end initialize_all_lists() */


void print_list_statistics (FILE *out)
{
    int i, j, l, m, n;
    /* collect allocation and occupation statistics */
    for (i=j=m=0,n=2L<<20; i<maxbin[3]; i++)
    {
	if (m < bindx[2*i+1]-bindx[2*i])
	    m = bindx[2*i+1]-bindx[2*i];
	if (n > bindx[2*i+1]-bindx[2*i])
	    n = bindx[2*i+1]-bindx[2*i];
	j += pow(bindx[2*i+1]-bindx[2*i], 2.);
    }
    fprintf (out, "Bin-particle list (%.3f Mbytes allocated, %.1f per "
	     "bin):\n", (double)bindx[2*maxbin[3]]*sizeof(int)/(2L<<20),
	     (double)bindx[2*maxbin[3]]/maxbin[3]);
    fprintf (out, "max = %d, min = %d, avg occup. = %.1f = %.1f%%, "
	     "fluc = %.2f.\n", m, n, (double)np/maxbin[3],
	     100.*np/bindx[2*maxbin[3]],
	     sqrt((double)j/maxbin[3]-pow((double)np/maxbin[3],2.)));
    for (i=j=l=m=0,n=2L<<20; i<np; i++)
    {
	l += idx[2*i+1]-idx[2*i];
	j += pow(idx[2*i+1]-idx[2*i], 2.);
	if (m < idx[2*i+1]-idx[2*i])
	    m = idx[2*i+1]-idx[2*i];
	if (n > idx[2*i+1]-idx[2*i])
	    n = idx[2*i+1]-idx[2*i];
    }
    fprintf (out, "Particle-particle list (%.3f Mbytes allocated, "
	     "%.1f per particle):\n",
	     (double)idx[2*np]*sizeof(int)/(2L<<20),
	     (double)idx[2*np]/np);
    fprintf (out, "max = %d, min = %d, avg occup. = %.1f = %.1f%%, "
	     "fluc = %.2f.\n\n", m, n, (double)l/np, 100.*l/idx[2*np],
	     sqrt((double)j/np-pow((double)l/np,2.)));
    return;
} /* end print_list_statistics() */


void accumulate_and_save_gr (char fn_gr[])
{
    static int nw_gr = 0;
    static long igr[GRMESH+1][3] = {0};
    static double TRsum = 0, densitysum = 0;
    int i, j, k, l, m, n, w, gb[3], over, overbin[GRMAXBINS];
    double cc, dd, ee, ff, gg;
    FILE *gr;
    /* number of necessary bins to cover GRCUT */
    for (i=0; i<3; i++)
	gb[i] = ceil(GRCUT*sqrt(HI[0][i]*HI[0][i]+
				HI[1][i]*HI[1][i]+
				HI[2][i]*HI[2][i])*maxbin[i]);
    for (l=0; l<maxbin[3]; l++)
    {  /* for each bin */
	over = 0;
	/* go over its surrounding bins */
	for (i=-gb[0]; i<=gb[0]; i++) 
	    for (j=-gb[1]; j<=gb[1]; j++)
		for (k=-gb[2]; k<=gb[2]; k++)
		{  /* nearby bin index */
		    m = nbin(l/(maxbin[1]*maxbin[2])+i,
			     l%(maxbin[1]*maxbin[2])/maxbin[2]+j,
			     l%maxbin[2]+k);
		    /* if went over before, just skip */
		    for (w=0; w<over; w++) if (overbin[w] == m) break;
		    if (w == over)
		    {  /* we discover a new bin */
			if (over >= GRMAXBINS)
			{  /* overflow detection */
			    printf ("accumulate_and_save_gr: "
				    "GRMAXBINS = %d exceeded.\n", GRMAXBINS);
			    exit(1);
			}
			for (w=bindx[2*l]; w<bindx[2*l+1]; w++)
			    for (n=bindx[2*m]; n<bindx[2*m+1]; n++)
				if ( own_it(binlist[w], binlist[n]) )
				    if ( r2(binlist[w], binlist[n], s, s)
					 < GRCUT*GRCUT )
					++igr[(int)floor(r[3]/r2del)]
					    [tpd[binlist[w]]+tpd[binlist[n]]];
			/* remember this bin */
			overbin[over++] = m;
		    }
		}
    }
    nw_gr++;  /* number of times igr[][] has been accumulated */
    TRsum += T_md_to_T_real(kine*UENERGY_IN_J/1.5/neff/BOLZ);
    densitysum += np/AVO/volume/UVOLUME_IN_M3;
    gr = fopen (fn_gr, "w+");
    fprintf (gr, "clf; nw_gr = %d; GR = [\n", nw_gr);
    for (i=0; i<GRMESH; i++)
    {  /* shell's exact volume */
	cc = 4*PI/3*(pow((i+1)*r2del,1.5)-pow(i*r2del,1.5));
	/* cc at front to prevent int multiplication overflow */
	dd = cc*nw_gr*(np-nar)*(np-nar-1)/volume/2.;
	ee = cc*nw_gr*nar*(np-nar)/volume;
	ff = cc*nw_gr*nar*(nar-1)/volume/2.;
	gg = dd + ee + ff;
	if (fabs(dd) > TINY) dd = 1./dd;
	if (fabs(ee) > TINY) ee = 1./ee;
	if (fabs(ff) > TINY) ff = 1./ff;
	if (fabs(gg) > TINY) gg = 1./gg;
	fprintf (gr, "%f %f",
		 cbrt((pow((i+1)*r2del,1.5)+pow(i*r2del,1.5))/2.)*
		 ULENGTH_IN_A, igr[i][0]*dd);
	if (ntpd > 0)
	    fprintf (gr, " %f %f %f", igr[i][1]*ee, igr[i][2]*ff,
		     (igr[i][0]+igr[i][1]+igr[i][2])*gg);
	fprintf (gr, "\n");
    }
    fprintf (gr, "]; plot(GR(:,1),GR(:,2));\n");
    fprintf (gr, "v = axis; axis([0 %f v(3) v(4)]);\n",
	     GRCUT*ULENGTH_IN_A);
    fprintf (gr, "xlabel('r [Angstrom]'); ylabel('g(r)');\n");
    fprintf (gr, "title('RDF of Argon at T = %f K, \\rho = %f "
	     "mol/m^3');\n", TRsum/nw_gr, densitysum/nw_gr);
    if (nar > 0)
    {
	fprintf (gr, "hold on; plot(GR(:,1),GR(:,3),'r');\n");
	fprintf (gr, "plot(GR(:,1),GR(:,4),'g');\n");
	fprintf (gr, "plot(GR(:,1),GR(:,5),'b');\n");
    }
    fclose (gr);
    return;
} /* end accumulate_and_save_gr() */


void print_atom_statistics (FILE *out, char *sp)
{
    int j;
    fprintf(out,"%sThere are %d atoms other than \"%2s of %.4f amu\".\n",
	    sp, nk, DEFAULT_ATOM_TP, ARGON_MASS_IN_KG*1000*AVO); 
    fprintf(out,"%saltogether %d different kinds will coexist:\n",
	    sp, nkind);
    fprintf(out,"%s--------------------------------------\n", sp);
    fprintf(out,"%s   Type   Mass(amu)    Occurrences\n", sp);
    for (j=0; j<nkind; j++)
	fprintf(out,"%s    %2s   %8.3f  %7d(%5.2f%%)\n", sp,
		tp[kind[j][0]], mass[kind[j][0]]*UMASS_IN_KG*1000*AVO,
		kind[j][1], 100.*kind[j][1]/np);
    fprintf(out,"%s--------------------------------------\n", sp);
    if (ntpd > 0)
    {  /* even the interactions are different */
	fprintf(out,"%samong which %d kinds: ", sp, ntpd);
	for (j=0; j<nkind; j++)
	    if (tpd[kind[j][0]])
		fprintf(out,"%2s(%.3f), ", tp[kind[j][0]],
			mass[kind[j][0]]*UMASS_IN_KG*1000*AVO);
	fprintf(out,"are \"Not %2s\",\n", DEFAULT_ATOM_TP);
	fprintf(out,"%sand therefore interact with scaled parameters,\n", sp);
	fprintf(out,"%s--------------------------------------\n", sp);
	fprintf(out,"%s         sigma-ratio   epsilon-ratio\n", sp);
	fprintf(out,"%s%2s-%2s     %f       %f\n", sp,
		DEFAULT_ATOM_TP, DEFAULT_ATOM_TP,
		l_ratio[0], e_ratio[0]);
	fprintf(out,"%s%2s-N%2s    %f       %f\n", sp,
		DEFAULT_ATOM_TP, DEFAULT_ATOM_TP,
		l_ratio[1], e_ratio[1]);
	fprintf(out,"%sN%2s-N%2s   %f       %f\n", sp,
		DEFAULT_ATOM_TP, DEFAULT_ATOM_TP,
		l_ratio[2], e_ratio[2]);
	fprintf(out,"%s--------------------------------------\n\n", sp);
	waitfor (TIME_TO_NOTICE);
	/* the bins have to accommodate the worst possibility */
	maxl_ratio = max(l_ratio[0], l_ratio[1]);
	maxl_ratio = max(maxl_ratio, l_ratio[2]);
	fprintf(out,"%sBecause of it, the system bin size has to increase\n"
		"%sfrom pure-%2s by maxl_ratio = %f; and\n", sp, sp,
		DEFAULT_ATOM_TP, maxl_ratio);
	/* particles of either species drift somewhat differently */
	minl_ratio[0] = min(l_ratio[0], l_ratio[1]);
	minl_ratio[1] = min(l_ratio[2], l_ratio[1]);
	fprintf(out,"%sMAX_DRIFT's of %2s, N%2s respectively should decrease "
		"by\n", sp, DEFAULT_ATOM_TP, DEFAULT_ATOM_TP);
	fprintf(out,"%sminl_ratio[0] = %f, minl_ratio[1] = %f.\n", sp,
		minl_ratio[0], minl_ratio[1]);
    }
    if (nmb != np)
    {
	fprintf(out,"\n%sAnd, we have found %d atoms with nominal "
		"masses\n", sp, np-nmb);
	fprintf(out,"%s> %.3f AMU, and they will be treated as "
		"immobile.\n", sp, IMMOBILE_MASS);
    }
    fprintf(out,"\n");
    waitfor (TIME_TO_NOTICE);
    return;
} /* end print_atom_statistics() */


/* Approximate linear regression of vx and quadratic */
/* regression of T for monoatomic fluids in sz1-sz2. */
#define ORDER_V 1
#define ORDER_T 2
#define MAX_ORDER (max(ORDER_V,ORDER_T))
int TFE (double sz1, double sz2)
{ /* approximate thermodynamic field estimator */
    int i, j, k;
    double cc, vxavg, kei, sz[2*MAX_ORDER+1], vx_sz[ORDER_V+1], 
	ke_sz[ORDER_T+1], S[ORDER_T+1][ORDER_T+1],
	SI[ORDER_T+1][ORDER_T+1];
    for (i=0; i<2*ORDER_V+1; i++) sz[i] = 0.;
    for (i=0; i<ORDER_V+1; i++) vx_sz[i] = 0.;
    /* estimate vx */
    for (i=0; i<np3; i+=3)
	if ( mb[i/3] && (!tpd[i/3]) &&  
	     (s[i+2] >= min(sz1,sz2)) &&
	     (s[i+2] <= max(sz1,sz2)) )
	    for (cc=1,k=0; k<2*ORDER_V+1; k++,cc*=s[i+2])
	    { /* sz[k]=\sum (sz_i)^k; vx_sz[k]=\sum vx_i*(sz_i)^k */
		sz[k] += cc;
		if (k<ORDER_V+1) vx_sz[k] += v[i] * cc;
	    }
    if (sz[0] > 0)
    { /* there are atoms in the region */
	sz[1] /= sz[0];     /* avg. s   */
	sz[2] /= sz[0];     /* avg. s*s */
	vx_sz[0] /= sz[0];  /* avg. v   */
	vx_sz[1] /= sz[0];  /* avg. v*s */
	af[0] = ( vx_sz[0]*sz[2] - sz[1]*vx_sz[1] ) /
	    (sz[2] - sz[1]*sz[1]);
	af[1] = ( vx_sz[1] - sz[1]*vx_sz[0] ) /
	    (sz[2] - sz[1]*sz[1]);
	for (i=0; i<2*ORDER_T+1; i++) sz[i] = 0.;
	for (i=0; i<ORDER_T+1; i++) ke_sz[i] = 0.;
	/* then estimate the kinetic energy based on this */
	for (i=0; i<np3; i+=3)
	    if ( mb[i/3] && (!tpd[i/3]) &&
		 (s[i+2] >= min(sz1,sz2)) &&
		 (s[i+2] <= max(sz1,sz2)) )
	    { /* stream velocity */
		vxavg = af[0] + af[1] * s[i+2];
		kei = mass[i/3] / 3 * UENERGY_IN_J / BOLZ *
		    ( square(v[i+1]) + square(v[i+2]) + 
		      square(v[i]-vxavg) );
		for (cc=1,k=0; k<2*ORDER_T+1; k++,cc*=s[i+2])
		{ /* ke_sz[k] = \sum kei * (sz_i)^k */
		    sz[k] += cc;
		    if (k<ORDER_T+1) ke_sz[k] += kei * cc;
		}
		for (j=0; j<ORDER_T+1; j++)
		    for (k=0; k<ORDER_T+1; k++)
			S[j][k] = s[j+k];
		matinv (S, SI); /* only for 3 x 3 matrix */
		for (j=0; j<ORDER_T+1; j++)
		    for (af[2+j]=k=0; k<ORDER_T+1; k++)
			af[2+j] += SI[j][k] * ke_sz[k];
	    }
    }
    return ((int)sz[0]);
} /* end TFE() */


int main (int argc, char *argv[])
{
    int i, j, k, l, m, n;
    double cc, dd, ee, ff, tmp, *calm;
    struct schedule_list *sch;
    
    for (buf[0]=(char)0,i=0; i<argc; i++)
	strcat(strcat(buf,argv[i]),(i==argc-1)?"":" ");
    printf ("\nProgram \"%s\" starts (pid = %d).\n",
	    buf, getpid());
    watch ("mission", 0); /* entire code clock */
    
    /* Read instructions from "con" file */
    /* preprocess only control variables */
    printf ("\nReading in instructions ...\n\n");
    fscanf (stdin, "Total number of MD steps:\n");
    fscanf (stdin, "%d\n\n", &maxkb);
    
    fscanf (stdin, "Starts property averaging at step:\n");
    fscanf (stdin, "%d\n\n", &kstart);
    if (kstart >= maxkb)
    {
	printf ("error: kstart >= MAXKB, "
		"leaving no room for averaging.\n");
	exit(1);
    }
    else if (kstart >= 2*maxkb/3)
    {
	printf ("** Warning: avg. period small compared with MAXKB.\n\n");
	waitfor (TIME_TO_NOTICE);
    }
    
    fscanf (stdin,
	    "Load configuration from file ('NULL'->start anew):\n");
    fscanf (stdin, "%s\n\n", fn_config_read);
    
    fscanf (stdin, "NTH or NTS mode, VOLUMETRIC or free:\n");
    fscanf (stdin, "%s", buf);
    NTH = !strcasecmp(buf, "NTH");
    NTS = !strcasecmp(buf, "NTS");
    if ( (!NTH) && (!NTS) )
    {
	printf("Error: \"%s\" run mode is invalid,\n", buf);
	printf("must select between NTH or NTS mode.\n");
	return(1);
    }
    /* VOLUMETRIC in NTS means keeping orthogonality; */
    /*      in NTH means keeping the volume constant  */
    fscanf (stdin, "%s\n\n", buf);
    VOLUMETRIC = !strcasecmp(buf, "VOLUMETRIC");
    
    fscanf (stdin, "Default schedule of the reference unit cell [A]:\n");
    /* HSCH: if there is an active schedule for H */
    /* HCALM: step at which activities in H cease */
    for (HSCH=HCALM=i=0; i<3; i++)
	for (j=i; j<3; j++)
	{
	    for (k=0; (buf[k]=getc(stdin))!='\n'; k++);
	    buf[k] = '\0';
	    sch = add_schedule (Hschname[od(i,j)], buf, 0, maxkb);
	    /* if none of them intends to change, then no change! */
	    HSCH = HSCH || sch->nontrivial;
	    calm = calm_period (Hschname[od(i,j)], kstart, maxkb);
	    HCALM = (calm==NULL)? maxkb : max(HCALM,calm[0]);
	}
    /* If NTS mode, then H schedules are automatically abandoned */
    if (NTS) HSCH = 0;
    fscanf (stdin, "\n");


    fscanf (stdin, "nc(1,2,3), number of wall layers among nc(3):\n");
    fscanf (stdin, "%d %d %d %d\n\n", &nc[0], &nc[1], &nc[2], &nc[3]);

    fscanf (stdin, "OPC vx [reduced], T [K], sz_OPC, sz_TFE_bottom:\n");
    fscanf (stdin, "%lf %lf %lf %lf\n\n", &df[0], &df[1], &df[2], &df[3]);
    
    fscanf (stdin, "Default multiplications of the unit cell "
	    "(also structure factor):\n");
    /* Irrespective of whether the config is loaded, nc[] will be */
    /* used as Miller's index for a structure factor calculation. */
    fscanf (stdin, "%d %d %d\n\n", &nc[0], &nc[1], &nc[2]);
    
    fscanf (stdin, "In NTH, use the above H schedule "
	    "instead of the configuration file:\n");
    fscanf (stdin, "%s\n\n", buf);
    /* If it is NULL, you don't have any other choice.. */
    if ( strcasecmp(fn_config_read, "NULL") )
    { /* if NTS, then we SHOULD start with config; ignore this question */
	USE_CONFIG_H = NTS ||
	    ((buf[0]!='1') && (buf[0]!='y') && (buf[0]!='Y'));
	if ( USE_CONFIG_H )
	{ /* abort default H schedule */
	    HSCH = 0;
	    HCALM = 0;
	    printf ("Do NOT override H from \"%s\" --\n"
		    "default H schedule abandoned.\n\n",
		    fn_config_read);
	}
	else printf ("OVERRIDE H from \"%s\" --\n"
		     "will use default H schedule (%s).\n\n",
		     fn_config_read,
		     HSCH?(HCALM<=kstart?"nontrivial but calm":
			   "nontrivial and not calm"):"trivial");
	/* let user see the notice */
	waitfor (TIME_TO_NOTICE);
    }
    
    fscanf (stdin, "Default external stress [MPa] schedule:\n");
    for (SSCH=SCALM=i=0; i<3; i++)
	for (j=i; j<3; j++)
	{
	    for (k=0; (buf[k]=getc(stdin))!='\n'; k++);
	    buf[k] = '\0';
	    sch = add_schedule (Sschname[od(i,j)], buf, 0, maxkb);
	    SSCH = SSCH || sch->nontrivial;
	    calm = calm_period (Sschname[od(i,j)], kstart, maxkb);
	    SCALM = (calm==NULL)? maxkb : max(SCALM,calm[0]);	    
	}
    /* of course, only in NTS mode can we apply a stress schedule! */
    if (NTH) SSCH = 0;
    fscanf (stdin, "\n");
    
    fscanf (stdin, "Save configuration on file:\n");
    fscanf (stdin, "%s\n\n", fn_config_write);
    
    fscanf (stdin, "Pair potential to be used ('LJ6_12' or 'WOON'):\n");
    fscanf (stdin, "%s\n\n", potfun_name);
    LJ6_12 = !strcasecmp (potfun_name, "LJ6_12");
    WOON   = !strcasecmp (potfun_name, "WOON");
    if ( (!LJ6_12) && (!WOON) )
    {
	printf("Error: select between LJ6_12 or WOON potential.\n");
	exit(1);
    }
    
    fscanf (stdin, "Request calculating thermal conductivity "
	    "and shear viscosity?\n");
    fscanf (stdin, "%s\n\n", buf);
    CALC = (buf[0]=='1') || (buf[0]=='y') || (buf[0]=='Y');
    
    fscanf (stdin, "MAX_STRAIN for strain-triggered "
	    "bin re-partition:\n");
    fscanf (stdin, "%s\n\n", buf);
    MAX_STRAIN =  strcasecmp(buf, "default")? fabs(atof(buf)) :
	(NTH&&(!HSCH))? 0.001 : DEFAULT_MAX_STRAIN;
    /* even in the dullest case, leave some margin for timestepping */
    
    fscanf (stdin, "MAX_DRIFT for particle re-anchoring:\n");
    fscanf (stdin, "%s\n\n", buf);
    MAX_DRIFT = strcasecmp(buf, "default")?fabs(atof(buf)):
	DEFAULT_MAX_DRIFT;
    /* To make sure that all pairs < L(i,j)RCUT will be found IF their */
    /* anchors < L(i,j)RLIST at a prior moment in this strain session */
    /* and they are currently < min(L(i,j))MAX_DRIFT*RCUT <= */
    /* L(i,j)MAX_DRIFT*RCUT from their anchors. */
    RLIST = sqrt((1+2*MAX_STRAIN)/(1-2*MAX_STRAIN)) *
	(1+2*MAX_DRIFT)*RCUT;
    
    fscanf (stdin, "GRCUT for g(r) calculation [A]:\n");
    fscanf (stdin, "%s\n\n", buf);
    GRCUT = strcasecmp(buf, "default")? fabs(atof(buf)):
	DEFAULT_GRCUT_IN_A;
    GRCUT /= ULENGTH_IN_A;
    
    fscanf (stdin, "Temperature [K] schedule:\n");
    for (k=0; (buf[k]=getc(stdin))!='\n'; k++);
    buf[k] = '\0';
    sch = add_schedule ("REALT", buf, 0, maxkb);
    TSCH = sch->nontrivial;
    calm = calm_period ("REALT", kstart, maxkb);
    TCALM = (calm==NULL)? maxkb : calm[0];
    fscanf (stdin, "\n");
    
    fscanf (stdin, "Integration timestep size [fs]:\n");
    fscanf (stdin, "%lf\n\n", &delta);
    
    fscanf (stdin, "Temperature/stress coupling time constant [ps]:\n");
    fscanf (stdin, "%s\n\n", buf);
    /* debugging purpose: turn off coupling to check conservation */
    tau = ( strncasecmp(buf,"inf",3) && strcasecmp(buf,"huge") )?
	fabs(atof(buf) ): -1;
    
    fscanf (stdin, "Frequency of output on screen and disk:\n");
    fscanf (stdin, "%d\n\n", &kw_lp);
    
    fscanf (stdin, "Random number generator seed (0->time seed):\n");
    fscanf (stdin, "%d\n\n", &iseed); 
    if (iseed == 0) iseed = time(NULL);
    srandom ((unsigned)iseed);
    
    fscanf (stdin, "Create voids at (sx sy sz radius[A]):\n");
    for (nvoids=0; ; nvoids++)
    {
	while ( ((buf[0]=getc(stdin)) != '\n') &&
		(buf[0] != (char)EOF) &&
		(buf[0] != 'n') &&
		(buf[0] != 'N') &&
		(buf[0] != '(') );
	if (buf[0] == '(')
	{
	    if (nvoids >= MAX_VOIDS)
	    {
		printf ("Error: MAX_VOIDS = %d exceeded.\n", MAX_VOIDS);
		exit(1);
	    }
	    fscanf (stdin, "%lf %lf %lf %lf)",
		    &voids[nvoids][0], &voids[nvoids][1],
		    &voids[nvoids][2], &voids[nvoids][3]);
	    voids[nvoids][3] /= ULENGTH_IN_A;
	}
	else
	{
	    if ((buf[0] == 'n')||(buf[0] == 'N'))
		while ( ((buf[0]=getc(stdin) != '\n')) &&
			(buf[0] != (char)EOF) );
	    break;
	}
    }
    fscanf (stdin, "\n");
    
    fscanf (stdin, "Save on temporary config (NULL=don't) "
	    "with frequency:\n%s %s\n", tmp_config, buf);
    if (!strcasecmp(tmp_config, "default"))
	strcpy(tmp_config, DEFAULT_TMP_CONFIG);
    tmp_config_freq = (!strcasecmp(buf, "default"))?
	DEFAULT_TMP_CONFIG_FREQ*kw_lp: atoi(buf);
    
    printf ("Instructions read complete.\n\n");
    
    /* first determine the actual number of particles */
    if (strcasecmp(fn_config_read, "NULL"))
    {
	if ((config=fopen(fn_config_read,"r")) == NULL)
	{
	    printf("Error: cannot open config file \"%s\".\n",
		   fn_config_read);
	    exit(1);
	}
	printf ("The initial configuration is loaded from \"%s\".\n",
		fn_config_read);
	fscanf (config, "Number of particles = %d\n\n", &np);
    }
    else
    {
	printf ("Initial config. is structurally perfect "
		"fcc crystal.\n");
	np = DEFAULT_UNIT_BASE*nc[0]*nc[1]*nc[2];
    }
    
    printf ("Number of particles in this simulation = %d.\n\n", np);
    np3 = 3 * np;
    
    /* allocate memory for individual particle properties */
    sa   = (double *) malloc (np3*sizeof(double));
    s    = (double *) malloc (np3*sizeof(double));
    s1   = (double *) malloc (np3*sizeof(double));
    s2   = (double *) malloc (np3*sizeof(double));
    s3   = (double *) malloc (np3*sizeof(double));
    s4   = (double *) malloc (np3*sizeof(double));
    s5   = (double *) malloc (np3*sizeof(double));
    v    = (double *) malloc (np3*sizeof(double));
    f    = (double *) malloc (np3*sizeof(double));
    fs   = (double *) malloc (np3*sizeof(double));
    d    = (double *) malloc (np3*sizeof(double));
    mass = (double *) malloc (np*sizeof(double));
    ep   = (double *) malloc (np*sizeof(double));
    mybin = (int *) malloc (np*sizeof(int));
    /* 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;
    /* 0 for Ar; 1 for NAr, or "not Ar" */
    tpd = (int *) malloc (np*sizeof(int));
    /* 0 for immobile; 1 for mobile */
    mb  = (int *) malloc (np*sizeof(int));
    
    if ( strcasecmp(fn_config_read, "NULL") )
    {  /* if loading the configuration from a file */
	fscanf (config, "H(1,1) = %lf A\n", &H[0][0]);
	fscanf (config, "H(1,2) = %lf A\n", &H[0][1]);
	fscanf (config, "H(1,3) = %lf A\n", &H[0][2]);
	fscanf (config, "nc(1)  = %d\n\n",  &maxbin[0]);
	fscanf (config, "H(2,1) = %lf A\n", &H[1][0]);
	fscanf (config, "H(2,2) = %lf A\n", &H[1][1]);
	fscanf (config, "H(2,3) = %lf A\n", &H[1][2]);
	fscanf (config, "nc(2)  = %d\n\n",  &maxbin[1]);
	fscanf (config, "H(3,1) = %lf A\n", &H[2][0]);
	fscanf (config, "H(3,2) = %lf A\n", &H[2][1]);
	fscanf (config, "H(3,3) = %lf A\n", &H[2][2]);
	fscanf (config, "nc(3)  = %d\n\n",  &maxbin[2]);
	if (USE_CONFIG_H)
	{
	    nc[0] = maxbin[0];
	    nc[1] = maxbin[1];
	    nc[2] = maxbin[2];
	}
	fscanf (config, "TP Mass(amu)     sx        sy        sz"
		"     sx1(1/ns)  sy1(1/ns)  sz1(1/ns)\n");
	for (i=0; i<np3; i+=3)
	{
	    fscanf (config, "%s %lf %lf %lf %lf %lf %lf %lf\n",
		    tp[i/3], &mass[i/3], &s[i], &s[i+1], &s[i+2],
		    &s1[i],  &s1[i+1], &s1[i+2]);
	    /* mass was in amu -> reduced */
	    mass[i/3] /= (UMASS_IN_KG*1000*AVO);
	    /* saved s1[] are in 1/ns -> dimensionless */
	    s1[i]   *= delta / 1000 / 1000;
	    s1[i+1] *= delta / 1000 / 1000;
	    s1[i+2] *= delta / 1000 / 1000;
	    if (strncmp(tp[i/3], DEFAULT_ATOM_TP, TP_STRLEN))
		tpd[i/3] = 1;
	    else tpd[i/3] = 0;
	}
	fclose (config);
    }
    else
    { /* assign perfect fcc ABC layering if start anew */
	x[0] = 0.;    x[1] = 0.;   x[2] = 0.;
	x[3] = 0.5;   x[4] = 0.5;  x[5] = 0.;
	x[6] = 1./6;  x[7] = 0.5;  x[8] = 1./3;
	x[9] = 2./3;  x[10] = 0.;  x[11] = 1./3;
	x[12] = 5./6; x[13] = 0.5; x[14] = 2./3;
	x[15] = 1./3; x[16] = 0.;  x[17] = 2./3;
	for (i=0; i<np; i++) f[i] = frandom();
	quicksort (np, f, mybin);
	/* mybin[] serving as tmp variable */
	for (i=0; i<np; i++)
	{ /* randomly choose certain particles to be NAR */
	    if (mybin[i] < DEFAULT_NAR_NUMBER) tpd[i] = 1;
	    else tpd[i] = 0;
	    /* if one wants to randomize particle indices for */
	    /* load-balancing purpose, comment the next line: */
	    mybin[i] = i;
	}
	for (m=i=0; i<nc[0]; i++)
	    for (j=0; j<nc[1]; j++)
		for (k=0; k<nc[2]; k++,m+=3*DEFAULT_UNIT_BASE)
		    for (l=0; l<3*DEFAULT_UNIT_BASE; l+=3)
		    {
			n = mybin[(m+l)/3];
			s[3*n]   = (x[l]  +i+0.25)/nc[0];
			s[3*n+1] = (x[l+1]+j+0.25)/nc[1];
			s[3*n+2] = (x[l+2]+k+0.25)/nc[2];
			/* make these layers immobile walls */
			if (k < nc[3]) tpd[n] = 1;
			mass[n]  = ARGON_MASS_IN_KG / UMASS_IN_KG;
			if (tpd[n]) mass[n] *= DEFAULT_NAR_MASS_RATIO;
			strncpy (tp[n], tpd[n]?
				 DEFAULT_NAR_TP:DEFAULT_ATOM_TP, TP_STRLEN);
			*(tp[n]+TP_STRLEN) = (char)0;
		    }
    }
    
    /* formally introduce H[][] */
    for (i=0; i<3; i++)
	for (j=2; j>=0; j--)
	{
	    if ( USE_CONFIG_H && (j>=i) )
	    { /* delete default schedule input from "con" */
		cancel_schedule (Hschname[od(i,j)]);
		/* use those loaded from configuration file and */
		/* symmetrize within the "reference unit cell". */
		sprintf (buf, "%f", (H[i][j]/nc[i]+H[j][i]/nc[j])/2);
		add_schedule (Hschname[od(i,j)], buf, 0, maxkb);
	    }
	    /* H[][] was in A -> reduced */
	    H[i][j] = nc[i] * schedule(Hschname[od(i,j)], 0) /
		ULENGTH_IN_A;
	    H1[i][j] = H2[i][j] = H3[i][j] = H4[i][j] = H5[i][j] = 0.;
	}
    volume = matinv (H, HI);
    mateqv (HI, HIBEGIN);
    
    if (HSCH) printf("There is a nontrivial H schedule, with\n");
    printf ("            | %f %f %f |\n", H[0][0]*ULENGTH_IN_A,
	    H[0][1]*ULENGTH_IN_A, H[0][2]*ULENGTH_IN_A);
    printf ("Initial H = | %f %f %f | A, nc = (%d %d %d)\n",
	    H[1][0]*ULENGTH_IN_A, H[1][1]*ULENGTH_IN_A,
	    H[1][2]*ULENGTH_IN_A, nc[0], nc[1], nc[2]);
    printf ("            | %f %f %f |\n", H[2][0]*ULENGTH_IN_A,
	    H[2][1]*ULENGTH_IN_A, H[2][2]*ULENGTH_IN_A);
    /* sanity check */
    if (nc[0]*nc[1]*nc[2] <= 0)
    {
	printf ("Error: illegal nc[] partition of H into "
		"reference \"unit cells\".");
	exit(1);
    }
    if ( HSCH && (HCALM <= kstart) )
	printf("but which stops changing after KSTART;\n");
    
    /* create voids in the configuration */
    for (i=k=0; i<np; i++)
    {
	for (j=0; j<nvoids; j++)
	{ /* centered around void j */
	    ds (&s[3*i], &voids[j][0], f);
	    r[0] = f[0]*H[0][0] + f[1]*H[1][0] + f[2]*H[2][0];
	    r[1] = f[0]*H[0][1] + f[1]*H[1][1] + f[2]*H[2][1];
	    r[2] = f[0]*H[0][2] + f[1]*H[1][2] + f[2]*H[2][2];
	    r[3] = r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
	    if (r[3] < voids[j][3]*voids[j][3]-TINY) break;
	}
	if (j==nvoids)
	{ /* eligible to survive */
	    s[3*k]   = s[3*i];
	    s[3*k+1] = s[3*i+1];
	    s[3*k+2] = s[3*i+2];
	    mass[k]  = mass[i];
	    strcpy (tp[k], tp[i]);
	    tpd[k]   = tpd[i];
	    k++;
	}
    }
    if ( np > k )
    {
	printf ("\nTo create %d voids, we throw away %d particles,\n",
		nvoids, np-k);
	printf ("and now has only %d particles.\n\n", np=k);
	waitfor (TIME_TO_NOTICE);
    }
    
    /* binary interaction parameters */
    l_ratio[0] = 1.; /* 0+0: Ar-Ar */
    /* 1+0/0+1: Ar-NAr */
    l_ratio[1] = DEFAULT_AN_L_RATIO;
    /* 1/1: NAr-NAr */
    l_ratio[2] = DEFAULT_NN_L_RATIO; 
    e_ratio[0] = 1.; 
    e_ratio[1] = DEFAULT_AN_E_RATIO;
    e_ratio[2] = DEFAULT_NN_E_RATIO; 
    
    /* atom kind statistics */
    for (nkind=i=0; i<np; i++)
    { /* different TP or mass constitutes a "kind" */
	for (j=0; j<nkind; j++)
	    if ( (!strcmp(tp[i], tp[kind[j][0]])) &&
		 (mass[i] == mass[kind[j][0]]) )
	    {
		kind[j][1]++;
		break;
	    }
	if ( j == nkind ) /* found a new kind of atom */
	    if ( nkind < MAX_KIND )
	    {
		kind[nkind][0]   = i;
		kind[nkind++][1] = 1;
	    }
	    else
	    {
		printf ("Error: MAX_KIND = %d exceeded.\n", MAX_KIND);
		exit(1);
	    }
    }
    
    for (nk=ntpd=nar=j=0; j<nkind; j++)
    { /* how many aliens are there? */
	if ( tpd[kind[j][0]] ||
	     (mass[kind[j][0]] != ARGON_MASS_IN_KG/UMASS_IN_KG) )
	    nk += kind[j][1];
	if (tpd[kind[j][0]])
	{ /* also, are they different interactions? */
	    ntpd++;
	    nar += kind[j][1];
	}
    }
    
    /* are there immobile atoms? */
    for (nmb=i=0; i<np; i++)
	if (mass[i]*UMASS_IN_KG*1000.*AVO >= IMMOBILE_MASS-TINY)
	    s1[3*i] = s1[3*i+1] = s1[3*i+2] = mb[i] = 0;
	else {nmb++; mb[i] = 1;}
    /* the center of mass is fixed: we "lose" one particle */
    neff = nmb - 1;
    
    maxl_ratio = minl_ratio[0] = minl_ratio[1] = 1.;
    /* these could be modified by print_atom_statistics() */
    if (nkind > 1) print_atom_statistics(stdout, "");
    
    printf ("Initial reduced density = %f = %.1f mol/m^3.\n\n",
	    np/volume, np/volume/AVO/UVOLUME_IN_M3);
    waitfor (TIME_TO_NOTICE);
    
    if (NTH)
    {
	printf ("Simulation will be carried out in NTH mode\n");
	if (VOLUMETRIC) printf ("keeping the initial volume invariant;\n");
    }
    else
    {
	printf ("Simulation will be carried out in NTS mode\n");
	for (i=0; i<3; i++)
	    for (j=i; j<3; j++)
	    {
		cc = schedule(Sschname[od(i,j)],0) / USTRESS_IN_MPA;
		sout[i][j] = sout[j][i] = cc;
	    }
	if (VOLUMETRIC)
	{
	    printf ("with hydrostatic constraints;\n");
	    trace_decompose (sout, sout, NULL);
	}
	if (SSCH)
	    printf("There is a nontrivial stress schedule with initial\n");
	else
	    printf("The loading stay constant during the simulation, with\n");
	printf ("                  | %f %f %f |\n",
		sout[0][0]*USTRESS_IN_MPA, sout[0][1]*USTRESS_IN_MPA,
		sout[0][2]*USTRESS_IN_MPA);
	printf ("external stress = | %f %f %f | MPa,\n",
		sout[1][0]*USTRESS_IN_MPA, sout[1][1]*USTRESS_IN_MPA,
		sout[1][2]*USTRESS_IN_MPA);
	printf ("                  | %f %f %f |\n",
		sout[2][0]*USTRESS_IN_MPA, sout[2][1]*USTRESS_IN_MPA, 
		sout[2][2]*USTRESS_IN_MPA);
	if ( SSCH && (SCALM <= kstart) )
	    printf("but which stops changing after KSTART;\n");
	/* tensor invariants */
	hydro = trace_decompose (sout, NULL, sout);
	shear = sqrt( sout[0][0]*sout[0][0]/2 + sout[1][1]*sout[1][1]/2 +
		      sout[2][2]*sout[2][2]/2 + sout[0][1]*sout[0][1]   +
		      sout[0][2]*sout[0][2]   + sout[1][2]*sout[1][2] );
	printf ("mean normal stress = %f MPa,\n", hydro*USTRESS_IN_MPA);
	printf ("von Mises shear stress = %f MPa,\n", shear*USTRESS_IN_MPA);
    }
    printf ("using the %s pair-potential.\n\n", potfun_name);
    waitfor (TIME_TO_NOTICE);
    
    if (TSCH) printf("We are planning for a nontrivial T schedule with\n");
    TR = schedule ("REALT", 0);
    printf ("Initial real temperature = %.2f K\n", TR);
    T = T_real_to_T_md (TR, &dTdR);
    if (TSCH) printf ("and final real temperature = %.2f K;\n",
		      schedule ("REALT", maxkb));
    if ( TSCH && (TCALM <= kstart) )
	printf("but the schedule stops changing after KSTART.\n");
    printf("\n");
    waitfor (TIME_TO_NOTICE);
    
    /* wall vibration time: 1/tau = sqrt(C*L/W) */
    wallmass = (tau>0)? ELASTIC_CONSTANT_IN_MPA /
	USTRESS_IN_MPA * cbrt(volume) * square(tau/UTIME_IN_PS) :
	1000. * ARGON_MASS_IN_KG / UMASS_IN_KG;
    
    printf ("MAXKB = %d\t KSTART = %d\n", maxkb, kstart);
    printf ("KW_LP = %d\t ISEED  = %d\n", kw_lp, iseed);
    printf ("delta = %.2f ns  tau    = %.2f ps\n", delta, tau);
    printf ("wallmass = %.2f Argon mass\n\n",
	    wallmass * UMASS_IN_KG / ARGON_MASS_IN_KG);
    waitfor (TIME_TO_NOTICE);
    
    /* time constants -> reduced */
    delta /= (UTIME_IN_PS*1000);
    tau   /= UTIME_IN_PS;
    
    /* system diagnostics */
    printf ("Across the entire simulation,\n");
    if (tau < 0) /* T coupling is turned off */
	if (NTH)
	    if (!HSCH)
		printf ("internal energy is conserved because the\n"
			"temperature coupling is turned off and\n"
			"there is no active H schedule.\n\n");
	    else
		printf ("internal energy may not be conserved although\n"
			"the temperature coupling is turned off because\n"
			"there is an active H schedule.\n\n");
	else  /* must be NTS mode */
	    if (!SSCH)
		if (VOLUMETRIC)
		    printf ("the hydrostatic enthalpy is conserved because\n"
			    "the temperature coupling is turned off, there\n"
			    "is NO active stress schedule + VOLUMETRIC.\n\n");
		else
		    printf ("generalized enthalpy is conserved (piecewise)\n"
			    "because the temperature coupling is turned off\n"
			    "and there is no active stress schedule.\n\n"); 
	    else
		if (VOLUMETRIC)
		    printf ("the hydrostatic enthalpy may not be conserved\n"
			    "although the temperature coupling is turned off\n"
			    "because there is an active stress schedule.\n\n");
		else
		    printf ("the generalized enthalpy may not be conserved\n"
			    "although the temperature coupling is turned off\n"
			    "because there is an active stress schedule.\n\n");
    else
	printf ("no dynamical quantities are conserved because\n"
		"the temperature coupling is turned on.\n\n");
    waitfor (TIME_TO_NOTICE);
    
    printf ("After KSTART = %d,\n", kstart);
    if (NTH)
	if ( (TCALM <= kstart) && (HCALM <= kstart) )
	    printf ("the conserved quantity is the internal energy\n"
		    "because there are no further T or H activities\n"
		    "and the system automatically switch to NEH mode.\n\n");
	else
	    printf ("there are no conserved dynamical quantities\n"
		    "because there are still activities in %s.\n\n",
		    (TCALM<=kstart)?"H":(HCALM<=kstart)?"T":"T and H");
    else /* must be NTS mode */
	if ( (TCALM <= kstart) && (SCALM <= kstart) )
	    printf ("the conserved quantity is the internal energy\n"
		    "as the system switches to NEH mode, since there\n"
		    "are no further T or S activities.\n\n");
	else
	    printf ("there are no conserved dynamical quantities\n"
		    "because there are still activities in %s.\n\n",
		    (TCALM<=kstart)?"S":(SCALM<=kstart)?"T":"T and S");
    waitfor (TIME_TO_NOTICE);
    
    if ( (!(NTH && (TCALM <= kstart) && (HCALM <= kstart))) &&
	 (!(NTS && (TCALM <= kstart) && (SCALM <= kstart))) )
    {
	printf("************************ WARNING! ************************\n");
	printf("Because the system cannot switch to NEH mode after KSTART,\n");
	printf("micro-canonical fluctuation formulas does not give the\n");
	printf("correct results for second-order thermodynamic derivatives\n");
	printf("that appear in the final printouts;\n");
	if ( CALC )
	{
	    CALC = 0;
	    printf("\nFurthermore, the request for thermal conductivity\n"
		   "calculation is DENIED because it requires NEH modes.\n");
	}
	printf("************************ WARNING! ************************\n");
	printf ("\n");
    }
    else if ( CALC )
	printf("Thermal conductivity and shear viscosity will be\n"
	       "correctly calculated.\n\n");
    waitfor (TIME_TO_NOTICE);    
    
    printf ("Property averages   => \"%s\";\n", fn_property);
    property = fopen (fn_property, "w+");
    printf ("Structure evolution => \"%s\";\n", fn_structure);
    structure = fopen (fn_structure, "w+");
    printf ("Temporary configur. => \"%s*\" (every %d steps);\n",
	    tmp_config, tmp_config_freq);
    printf ("Radial distribution => \"%s\";\n", GRFILE);
    if ( CALC )
    {  /* designated filenames in "smile.h" */
	printf ("Heat current => \"%s\", \"%s\", \"%s\";\n",
		fn_flux[SMILE_HEAT][0], fn_flux[SMILE_HEAT][1],
		fn_flux[SMILE_HEAT][2]);
	printf ("Shear stress => \"%s\", \"%s\", \"%s\";\n",
		fn_flux[SMILE_STRESS][0], fn_flux[SMILE_STRESS][1],
		fn_flux[SMILE_STRESS][2]);
	cc = delta * UTIME_IN_PS;
	for (i=0; i<3; i++)
	{ /* file handles are private variables here */
	    heat_current[i] = fopen (fn_flux[SMILE_HEAT][i],"w+");
	    /* the first number is the timestep size in picoseconds */
	    fwrite (&cc, sizeof(double), 1, heat_current[i]);
	    shear_stress[i] = fopen (fn_flux[SMILE_STRESS][i],"w+");
	    fwrite (&cc, sizeof(double), 1, shear_stress[i]);
	}
    }
    printf ("Final configuration => \"%s\".\n\n", fn_config_write);
    waitfor (TIME_TO_NOTICE);
    
    /* calculate the total mass and approximate kine in the system */
    for (cc=totalmass=0.,i=0; i<np; i++)
	if (mb[i])
	{
	    totalmass += mass[i];
	    if (strcasecmp(fn_config_read, "NULL"))
		cc += (s1[3*i]*s1[3*i] + s1[3*i+1]*s1[3*i+1] +
		       s1[3*i+2]*s1[3*i+2]) * pow(volume,2./3) /
		    delta / delta * mass[i] * UENERGY_IN_J / BOLZ / 3 / np;
	}
    
    if ( cc <= JUMPSTART_KINE )
    {
	/* minimum kinetic energy is needed to jump start MD */
	if (strcasecmp(fn_config_read, "NULL"))
	    printf ("Kinetic energy in \"%s\"\n= %.3f K < "
		    "JUMPSTART_KINE = %.3f K  -->\n",
		    fn_config_read, cc, JUMPSTART_KINE);
	printf ("Assign random velocities to particles ...\n");
	x[0] = x[1] = x[2] = 0;
	randnorm (np3, v);
	for (i=0; i<np3; i++)
	    if (mb[i/3])
	    {  /* assign Gaussian with <v^2> = kT/m */
		v[i]   /= sqrt(mass[i/3]);
		x[i%3] += mass[i/3] * v[i];
	    }
	for (cc=0,i=0; i<np3; i++)
	    if (mb[i/3])
	    {  /* make total momentum 0 */
		v[i] -= x[i%3] / totalmass;
		cc   += v[i] * v[i] * mass[i/3] / 2.;
	    }
	/* scale velocities to desired temperature; center of */
	/* mass is fixed so we lose three degrees of freedom; */
	/* phonon density is isotropic, and we start with twice */
	/* the desired kinetic energy since at potential ground */
	for (i=0; i<np3; i++)
	    if (mb[i/3])
		v[i] *= sqrt(3*neff*BOLZ*T/UENERGY_IN_J/cc);

	/* if the user wants to have a flow, give the fluid a push */
	if ( (df[2] >= 0. ) && (df[2] <= 1. ) && 
	     (i=TFE(df[3], df[2])) )
	{
	    printf ("By doing linear regression on the %d %s atoms "
		    "inside\n", i, DEFAULT_ATOM_TP);
	    printf ("sz = (%.5f, %.5f)  =>  vx(%.5f) = %f,\n",
		    df[3], df[2], df[2], cc=af[0]+af[1]*df[2]);
	    printf ("vx(%.5f=%.5f A) = 0 (wall at %.5f=%.5f A).\n\n",
		    -af[0]/af[1], -af[0]/af[1]*H[2][2]*ULENGTH_IN_M*1E10,
		    (double)nc[3]/nc[2], (double)nc[3]/nc[2]*
		    ULENGTH_IN_M*1E10);
	    /* apply corrections to particle velocities */
	    for (i=0; i<np3; i+=3)
		if ( mb[i/3] && (!tpd[i/3]) )
		    v[i] += (df[0]-cc) + (s[i+2]-df[2]) *
			(df[0]/(df[2]-nc[3]/nc[2]) - af[1]);
	    waitfor (TIME_TO_NOTICE);  
	}
		
	for (i=0; i<np3; i+=3)
	    if (mb[i/3])
	    {  /* convert to s1[], assuming H is constant */
		s1[i]   = ( v[i]*HI[0][0] + v[i+1]*HI[1][0] +
			    v[i+2]*HI[2][0] ) * delta;
		s1[i+1] = ( v[i]*HI[0][1] + v[i+1]*HI[1][1] +
			    v[i+2]*HI[2][1] ) * delta;
		s1[i+2] = ( v[i]*HI[0][2] + v[i+1]*HI[1][2] +
			    v[i+2]*HI[2][2] ) * delta;
	    }
    }
    
    for (i=0; i<np3; i++)
	s2[i] = s3[i] = s4[i] = s5[i] = 0.;

    printf ("RCUT = %.3f = %.3f A,  MAX_DRIFT = %.3f RCUT,\n",
	    RCUT, RCUT*ULENGTH_IN_A, MAX_DRIFT);
    printf ("Max |strain eigenvalue| for repartitioning = %.3f,\n",
	    MAX_STRAIN);
    printf ("-> possible shrinkage ratio because of that = %f;\n",
	    sqrt(1-2*MAX_STRAIN)/sqrt(1+2*MAX_STRAIN));
    printf ("=> RLIST = %.3f = %.3f A = %.3f RCUT.\n\n",
	    RLIST, RLIST*ULENGTH_IN_A, RLIST/RCUT);
    waitfor (TIME_TO_NOTICE);
    
    /* bin-bin, bin-particle, particle-bin, particle-particle */
    initialize_all_lists();
    nw_repartition = nw_anchor = 0;
    print_list_statistics (stdout);
    waitfor (TIME_TO_NOTICE);

    i = min(maxbin[0],maxbin[1]);
    i = min(i,maxbin[2]);
    /* maxl_ratio*RLIST is guaranteed to < bin thicknesses at any time */
    if (2*GRCUT > i*maxl_ratio*RLIST) GRCUT = i*maxl_ratio*RLIST/2.;
    /* g(r): how much data we get in a mesh in one instance */
    cc = np/2.*(np-1)/volume*4*PI*GRCUT*GRCUT*GRCUT/3/GRMESH;
    /* so how many instances are needed to achieve GRACCU */
    dd = ceil(1./GRACCU/GRACCU/cc);
    /* thus the sampling rate should be */
    kw_gr = ceil((maxkb-kstart)/dd);
    /* the meshes are regular in \delta r2 */
    r2del = GRCUT*GRCUT/GRMESH;
    
    printf ("GRCUT  = %.3f A\tGRMESH = %d\n",
	    GRCUT*ULENGTH_IN_A, GRMESH); 
    printf ("GRACCU = %.2f%%\t\tkw_gr  = %d\n",
	    100*GRACCU, kw_gr); 
    printf ("nc = [%d %d %d] (structure factor)\n\n",
	    nc[0], nc[1], nc[2]);
    waitfor (TIME_TO_NOTICE);
    
    /* anchor displacements from original positions */
    for (i=0; i<np3; i++) d[i] = 0.;
    
    /* clear the property cumulants */
    kinetum = k1tum = k2tum = 0.;
    kinesum = k1sum = k2sum = potesum = 0.;
    anchortum = 0;
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	{
	    Htum[i][j] = strestum[i][j] = 0.;
	    v1k1tum[i][j] = v1k2tum[i][j] = 0.;
	    Hsum[i][j] = stressum[i][j] = 0.;
	    v1k1sum[i][j] = 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.;
		}
	    Havg[i][j] = 0.;  
	}
    /* wait till transients are over */
    Havg_start = 0.5*max(SCALM,TCALM) + 0.5*kstart;
    for (i=0; i<6; i++)
	for (j=0; j<6; j++)
	{ 
	    ck1tum[i][j] = 0.;
	    ck1sum[i][j] = 0.; 
	}
    
    printf ("Simulation starts ...\n");
    for (kb=0; kb<maxkb; kb++)
    {
	/* NTS vibrates H; HSCH stretches H */
	scale_H = ( NTS && (!((SCALM<=kstart)&&
			      (TCALM<=kstart)&&
			      (kb>=kstart))) )  ||
	    ( HSCH && (!((HCALM<=kstart)&&(kb>=kstart))) );
	/* scale_H -> scale_T, because heat is generated */
	scale_T = (kb<kstart) || scale_H ||
	    ( TSCH && (!((HCALM<=kstart)&&(kb>=kstart))) );
	/* letting tau=inf effectively eliminates scale_T, */
	/* but that may be dangerous to do.                */
	
	/* external conditions */
	T = T_real_to_T_md (TR=schedule("REALT",kb), &dTdR);
	if ( NTH && scale_H )
	{
	    for (i=0; i<3; i++)
		for (j=0; j<3; j++) 
		    H[i][j] = schedule (Hschname[od(i,j)],kb) *
			nc[i] / ULENGTH_IN_A;
	    if (VOLUMETRIC)
	    { /* keep the old volume */
		cc = matinv (H, HI);
		for (i=0; i<3; i++)
		    for (j=0; j<3; j++)
		    { /* to prevent numerical instability, */
			/* volume shouldn't be written */
			H[i][j] = H[i][j] * cbrt(volume/cc);
			HI[i][j] = HI[i][j] / cbrt(volume/cc);
		    }
	    }
	    else volume = matinv (H, HI);
	}
	else if ( NTS && scale_H )
	{
	    for (i=0; i<3; i++)
		for (j=0; j<3; j++)
		    sout[i][j] = schedule (Sschname[od(i,j)],kb) /
			USTRESS_IN_MPA;
	    if (VOLUMETRIC)
		hydro = trace_decompose (sout, sout, NULL);
	    else hydro = trace_decompose (sout, NULL, sout);
	}
	
	for (msd=i=0; i<np; i++)
	{ /* check individual particle drifts from anchor */
	    if (mb[i])
		if (r2(i,i,sa,s) > square(minl_ratio[tpd[i]]*MAX_DRIFT*RCUT))
		{  /* d[] track displacements of anchors */
		    d[3*i]   += r[0];
		    d[3*i+1] += r[1];
		    d[3*i+2] += r[2];
		    re_anchor(i);
		    nw_anchor++;
		    anchortum++;
		    msd += d[3*i]*d[3*i]+d[3*i+1]*d[3*i+1]+d[3*i+2]*d[3*i+2];
		}
		else msd += (d[3*i]+r[0])*(d[3*i]+r[0])+
			 (d[3*i+1]+r[1])*(d[3*i+1]+r[1])+
			 (d[3*i+2]+r[2])*(d[3*i+2]+r[2]);
	    else /* clean out possible garbage */
		if (idx[2*i+1]==idx[2*i+2]) re_anchor(i);
	}
	msd /= neff;
  	if (kb==kstart) msd_kstart = msd;

	/* catastrophic NTS strain protection */
	if ( NTS )
	{ /* Lagrangian strain with respect to starting H */
	    matmul (HIBEGIN, H, A); 
	    matran (A, B);
	    matmul (A, B, G);
	    for (cc=i=0; i<3; i++)
		for (j=0; j<3; j++)
		    cc += square(G[i][j]-((i==j)?1.:0.)) / 4.;
	    if (cc >= NTS_MAX_STRAIN2)
	    {
		printf ("\nAt step %d, t = %.2f ps\n", kb,
			kb * delta * UTIME_IN_PS);
		printf ("    | %f %f %f |\n", H[0][0]*ULENGTH_IN_A,
			H[0][1]*ULENGTH_IN_A, H[0][2]*ULENGTH_IN_A);
		printf ("H = | %f %f %f | A,\n",
			H[1][0]*ULENGTH_IN_A, H[1][1]*ULENGTH_IN_A,
			H[1][2]*ULENGTH_IN_A);
		printf ("    | %f %f %f |\n", H[2][0]*ULENGTH_IN_A,
			H[2][1]*ULENGTH_IN_A, H[2][2]*ULENGTH_IN_A);
		printf ("Lagrangian strain norm2 >= NTS_MAX_STRAIN2 "
			"= %f, exit.\n", NTS_MAX_STRAIN2);
		save_tmp_config();
		exit(1);
	    }
	}
	
	/* Lagrangian strain with respect to old H */
	matmul (HIOLD, H, A); 
	matran (A, B);
	/* G will be re-used for getting hamiltonian */
	matmul (A, B, G);
	diag3 (G, x, GA);
	/* check if principle strains exceed bound */
	for (i=0; i<3; i++)
	{  /* epsilon = (J^T*J - 1)/2 */
	    x[i] = (x[i]-1) / 2.;
	    if (fabs(x[i]) >= MAX_STRAIN)
	    {
		printf ("|strain| in (%.4f %.4f %.4f) direction "
			"is now greater than\n",
			GA[i][0], GA[i][1], GA[i][2]);
		printf ("MAX_STRAIN = %f, repartition the system...\n",
			MAX_STRAIN);
		for (j=0; j<np; j++)
		{ /* since sa[] are going to be changed wholesale */
		    r2(j, j, sa, s);
		    d[3*j]   += r[0];
		    d[3*j+1] += r[1];
		    d[3*j+2] += r[2];
		}
		initialize_all_lists();
		nw_repartition++;
		/* at remesh incident, backup system configuration */
		save_tmp_config();
		break;
	    }
	}
	
	if ( (kb>=kstart) && ((kb-kstart)%kw_gr==0) )
	    accumulate_and_save_gr (GRFILE);
	
	/* Gear 5th-order predictor */
	for (i=0; i<np3; i++)
	    if (mb[i/3])
	    {
		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];
	    }
	
	if ( NTS && scale_H )
	{
	    for (kinewall=0,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]
			/ delta / delta / 2.;
		}
	    volume = matinv (H, HI);
	    /* H should change quasi-statically as if immersed in */
	    /* viscous oil -> 9kT/2 may be too big, let it be kT/4 */
	    HTnow = kinewall * UENERGY_IN_J / (0.25 * BOLZ);
	    factorHT = 1.-delta/2./tau*((HTnow+1.)/(T+1.)-1);
	    if (factorHT < 0.5) factorHT = 0.5;
	    else if (factorHT > 2) factorHT = 2;
	}
	
	/* real-space forces, stresses, elastic constants */
	pair_potential();
	Tnow = kine * UENERGY_IN_J / neff / (1.5 * BOLZ);
	factorT  = 1.-delta/2./tau*((Tnow+1.)/(T+1.)-1);
	if (factorT < 0.5) factorT = 0.5;
	else if (factorT > 2) factorT = 2;
	
	/* temporary fluctuation averages, */
	/* released every KW_LP steps.     */
	if (kb >= 1)
	{
	    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 the end of 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>0) && (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;
	    /* average kT during period */
	    cc = kinetum / (1.5*neff);
	    /* average volume */
	    dd = matinv(Htum,A) / kw_lp / kw_lp / kw_lp;
	    
	    /* 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/dd;
	    ck1tum[1][1] += 2*neff*cc/dd;
	    ck1tum[2][2] += 2*neff*cc/dd;
	    ck1tum[3][3] += neff*cc/dd;
	    ck1tum[4][4] += neff*cc/dd;
	    ck1tum[5][5] += neff*cc/dd;
	    
	    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 is system's response to external stress */
	    pressure = (strestum[0][0]+strestum[1][1]+
			strestum[2][2])/kw_lp/3.;
	    /* Parrinello & Rahman, J. Appl. Phys. 52, 7182 (1981) */
	    hamiltonian = pote + kine;
	    if ( scale_H )
		hamiltonian += kinewall - hydro*volume -
		    volumeold*(sout[0][0]*(G[0][0]-1)/2 +
			       sout[1][1]*(G[1][1]-1)/2 +
			       sout[2][2]*(G[2][2]-1)/2 +
			       sout[0][1]*G[0][1] + sout[0][2]*G[0][2] +
			       sout[1][2]*G[1][2]);
	    
	    fprintf (property, "%7d %.3f %.6f %.6f %.5f\n",
		     kb, Tnow, pote/np, hamiltonian/np, pressure);
	    fflush (property);

	    /* save meshed properties */
	    if (kb >= kstart) save_mesh (DEFAULT_MESH_FILE);
	    
	    /* screen output */
	    printf ("\nKB = %d, T_md = %f K, T_real = %f K\n",
		    kb, Tnow, T_md_to_T_real(Tnow));
	    printf ("Potential energy = %.5f\\%.5f eV"
		    "\\%.4f K\\%.1f J/mol\n",
		    pote/np, pote/np*UENERGY_IN_J/EV_IN_J,
		    pote/np*UENERGY_IN_J/BOLZ,
		    pote/np*UENERGY_IN_J*AVO);
	    printf ("Hamiltonian = %f\\%f eV\\%f K\\%.1f J/mol\n",
		    hamiltonian/np,
		    hamiltonian/np*UENERGY_IN_J/EV_IN_J,
		    hamiltonian/np*UENERGY_IN_J/BOLZ,
		    hamiltonian/np*UENERGY_IN_J*AVO);
	    printf ("System 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/UVOLUME_IN_M3*1E-30, 
		    np/volume/AVO/UVOLUME_IN_M3);
	    printf ("       | %.3f %.3f %.3f |",
		    Htum[0][0]*ULENGTH_IN_A/kw_lp,
		    Htum[0][1]*ULENGTH_IN_A/kw_lp,
		    Htum[0][2]*ULENGTH_IN_A/kw_lp);
	    printf ("            | %.3f %.3f %.3f |\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) = | %.3f %.3f %.3f |",
		    Htum[1][0]*ULENGTH_IN_A/kw_lp,
		    Htum[1][1]*ULENGTH_IN_A/kw_lp,
		    Htum[1][2]*ULENGTH_IN_A/kw_lp);
	    printf ("   t(MPa) = | %.3f %.3f %.3f |\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 ("       | %.3f %.3f %.3f |",
		    Htum[2][0]*ULENGTH_IN_A/kw_lp,
		    Htum[2][1]*ULENGTH_IN_A/kw_lp, 
		    Htum[2][2]*ULENGTH_IN_A/kw_lp);
	    printf ("            | %.3f %.3f %.3f |\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 ("         | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\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 ("         | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\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) = | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\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 ("         | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\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 ("         | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\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 ("         | %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f |\n\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);

	    print_list_statistics (stdout);
	    printf ("Re-anchoring rate = %f /atom/step.\n\n",
		    (double)anchortum/kw_lp/np);

	    /* <111> static structure factor of the system which */
	    /* should be 1 for perfect fcc lattice (differs with */
	    /* S(Q) in classical liquids by NP). If one needs    */
	    /* sensible result here for a loaded configuration   */
	    /* then input nc[] accordingly in "con".             */
	    for (cc=dd=i=0; i<np3; i+=3)
	    {
		ee = 2*PI*(nc[0]*s[i]+nc[1]*s[i+1]+nc[2]*s[i+2]);
		cc += cos(ee); dd += sin(ee);
	    }
	    structure_factor = (cc*cc+dd*dd)/np/np;
	    
	    /* structure file output */
	    fprintf (structure, "%f %f %f %f %f %f %f ",
		     kb*delta*UTIME_IN_PS,
		     T_md_to_T_real(kinetum*UENERGY_IN_J/neff/1.5/BOLZ),
		     pote/np*UENERGY_IN_J*AVO,
		     (pote+kine)/np*UENERGY_IN_J*AVO,
		     hamiltonian/np*UENERGY_IN_J*AVO,
		     structure_factor,
		     msd*ULENGTH_IN_A*ULENGTH_IN_A);
	    fprintf (structure, "%f %f %f %f %f %f %f ",
		     Htum[0][0]*ULENGTH_IN_A/kw_lp,
		     Htum[1][1]*ULENGTH_IN_A/kw_lp,
		     Htum[2][2]*ULENGTH_IN_A/kw_lp,
		     Htum[1][2]*ULENGTH_IN_A/kw_lp,
		     Htum[0][2]*ULENGTH_IN_A/kw_lp,
		     Htum[0][1]*ULENGTH_IN_A/kw_lp,
		     matinv(Htum,A)*pow(ULENGTH_IN_A/kw_lp,3.));
	    cc = trace_decompose (strestum, A, NULL);
	    dd = sqrt( A[0][0]*A[0][0]/2 + A[1][1]*A[1][1]/2 +
		       A[2][2]*A[2][2]/2 + A[0][1]*A[0][1]   +
		       A[0][2]*A[0][2]   + A[1][2]*A[1][2] );
	    fprintf (structure, "%f %f %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[1][2]*USTRESS_IN_MPA/kw_lp,
		     strestum[0][2]*USTRESS_IN_MPA/kw_lp,
		     strestum[0][1]*USTRESS_IN_MPA/kw_lp,
		     cc*USTRESS_IN_MPA/kw_lp,
		     dd*USTRESS_IN_MPA/kw_lp);
	    fprintf (structure, "%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);
	    fflush (structure);
	    
	    kinetum = k1tum = k2tum = 0.;
	    anchortum = 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 runtime outputs */
	
	/* once for a while, backup system configuration */
	if ( (kb > 0) && (kb % tmp_config_freq == 0) )
	    save_tmp_config();
	
	for (i=0; i<np3; i+=3)
	{
	    fs[i]   = ( f[i]*HI[0][0] + f[i+1]*HI[1][0] +
			f[i+2]*HI[2][0] ) / mass[i/3];
	    fs[i+1] = ( f[i]*HI[0][1] + f[i+1]*HI[1][1] +
			f[i+2]*HI[2][1] ) / mass[i/3];
	    fs[i+2] = ( f[i]*HI[0][2] + f[i+1]*HI[1][2] +
			f[i+2]*HI[2][2] ) /mass[i/3];
	}
	
	if ( NTS && scale_H )
	{
	    if (VOLUMETRIC) trace_decompose (stress, stress, NULL);
	    matran (H, HT);
	    matmul (H, HT, G); /* as my H is transpose of P&R's */
	    matinv (G, GI);
	    matran (H1, H1T);
	    matmul (H1, HT, A);
	    matmul (H, H1T, B);
	    matadd (A, B, G1);
	    matmul (GI, G1, GA); /* G^{-1}\dot{G}*delta in P&R */
	    /* generalized forces on s[] */
	    for (i=0; i<np3; i+=3)
	    {
		fs[i]   -= ( GA[0][0]*s1[i] + GA[0][1]*s1[i+1] +
			     GA[0][2]*s1[i+2] ) / delta / delta;
		fs[i+1] -= ( GA[1][0]*s1[i] + GA[1][1]*s1[i+1] +
			     GA[1][2]*s1[i+2] ) /delta / delta;
		fs[i+2] -= ( GA[2][0]*s1[i] + GA[2][1]*s1[i+1] +
			     GA[2][2]*s1[i+2] ) / delta / delta;
	    }
	    /* generalized forces on H[][]: hydrostatic term */
	    for (i=0; i<3; i++)
		for (j=0; j<3; j++)
		{
		    A[i][j] = stress[i][j];
		    if (i==j) A[i][j] += hydro;
		    A[i][j] *= volume / wallmass;
		}
	    matran (HI, HIT);
	    matmul (HIT, A, fH);
	    /* shear term */
	    matmul (HIOLD, H, G);
	    matmul (sout, G, G1);
	    matran (HIOLD, B);
	    matmul (B, G1, G);
	    for (i=0; i<3; i++)
		for (j=0; j<3; j++)
		{
		    fH[i][j] += G[i][j] * volumeold / wallmass;
		    tmp = H2[i][j] - fH[i][j] * delta * delta / 2.;
		    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;
		    /* prevent explosion by improper initial condition */
		    if (tau > 0) H1[i][j] *= factorHT;
		}
	    volume = matinv (H, HI);
	    if ( (kb >= Havg_start) && (kb < kstart) )
		matadd (H, Havg, Havg);
	}
	
	if ( NTS && (SCALM<=kstart) && (TCALM<=kstart) &&
	     (kb==kstart) && (kstart > Havg_start) )
	{
	    for (i=0; i<3; i++)
		for (j=0; j<3; j++)
		{ /* at this step scale_H is false */
		    H[i][j] = Havg[i][j] / (kstart-Havg_start);
		    H1[i][j] = 0.;
		    H2[i][j] = 0.;
		    H3[i][j] = 0.;
		    H4[i][j] = 0.;
		    H5[i][j] = 0.;
		}
	    volume = matinv (H, HI);
	}
	
	/* Gear 5th-order corrector */
	for (i=0; i<np3; i++)
	    if (mb[i/3])
	    {
		tmp = s2[i] - fs[i] * delta * delta / 2.;
		s[i]  -= tmp * F02;
		s1[i] -= tmp * F12;
		s2[i] -= tmp;
		s3[i] -= tmp * F32;
		s4[i] -= tmp * F42;
		s5[i] -= tmp * F52;
		/* periodic boundary conditions */
		if (s[i]<0.) s[i]++;
		else if (s[i]>1.) s[i]--;
		if ( scale_T && (tau>0) ) s1[i] *= factorT;
	    }
	
	if ( CALC && (kb>=kstart) )
	{  /* write to current files */
	    for (i=0; i<3; i++)
	    {   /* numerically massage the Green-Kubo formula */
		/* such that k = \sum_i<cj[0]cj[i]> in W/M/K. */
		cj[i] *= sqrt(UCONDUCTIVITY_IN_WMK*dTdR*
			      delta/(BOLZ*T/UENERGY_IN_J*T*volume));
		fwrite (&cj[i], sizeof(double), 1, heat_current[i]);
		for (j=0; j<3; j++)
		    stress[i][j] *=
			sqrt(USTRESS_IN_PA*UTIME_IN_SECOND*
			     delta*volume/(BOLZ*T/UENERGY_IN_J));
	    }
	    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]);
	    if ( kb % kw_lp == 0 )
	    {
		fflush(heat_current[0]);
		fflush(heat_current[1]);
		fflush(heat_current[2]);
		fflush(shear_stress[0]);
		fflush(shear_stress[1]);
		fflush(shear_stress[2]);
	    }
	}
    } /* main loop ends; */
    
    /* harvest 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 */
    Tnow = cc * UENERGY_IN_J / BOLZ;
    Treal = T_md_to_T_real (Tnow);
    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)/Tnow/bulk_modulus/3.;
    Cv = -(BOLZ*T/UENERGY_IN_J)/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];
    
    /* write in common property output file */
    fprintf (property, "\n   This is an MD simulation of "
	     "%d particles using\n", np);
    fprintf (property, "   %s pair potential with RCUT = %.3f "
	     "= %.4f A.\n\n", potfun_name, RCUT, RCUT*ULENGTH_IN_A);
    
    if (nkind>1) print_atom_statistics(property, "   ");
    
    fprintf (property, "   During the entire simulation (%d x %.4f ps),\n",
	     maxkb, delta*UTIME_IN_PS);
    fprintf (property, "   there are %d strain-triggered repartitions"
	     " (MAX_STRAIN = %.3f),\n", nw_repartition, MAX_STRAIN);
    fprintf (property,
	     "   and %d re-anchors, or %f /atom/step.\n\n",
	     nw_anchor, (double)nw_anchor/np/maxkb);
    
    fprintf (property, "   The final average starts from step "
	     "%d to %d.\n\n", kstart, maxkb);
    
    fprintf (property, "   The average MD temperature is %f = %f K,\n",
	     Tnow/T*(BOLZ*T/UENERGY_IN_J), Tnow);
    fprintf (property,
	     "   corresponding to a real temperature of %f K,\n", Treal);
    fprintf (property, "   with dT_md / dT_real = %f.\n\n", dTdR);
    
    fprintf (property,
	     "   H[1][1] = %8.4f A    stress[1][1] = %9.4f MPa\n",
	     Hsum[0][0]*ULENGTH_IN_A, stressum[0][0]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[1][2] = %8.4f A    stress[1][2] = %9.4f MPa\n",
	     Hsum[0][1]*ULENGTH_IN_A, stressum[0][1]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[1][3] = %8.4f A    stress[1][3] = %9.4f MPa\n",
	     Hsum[0][2]*ULENGTH_IN_A, stressum[0][2]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[2][1] = %8.4f A    stress[2][1] = %9.4f MPa\n",
	     Hsum[1][0]*ULENGTH_IN_A, stressum[1][0]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[2][2] = %8.4f A    stress[2][2] = %9.4f MPa\n",
	     Hsum[1][1]*ULENGTH_IN_A, stressum[1][1]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[2][3] = %8.4f A    stress[2][3] = %9.4f MPa\n",
	     Hsum[1][2]*ULENGTH_IN_A, stressum[1][2]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[3][1] = %8.4f A    stress[3][1] = %9.4f MPa\n",
	     Hsum[2][0]*ULENGTH_IN_A, stressum[2][0]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[3][2] = %8.4f A    stress[3][2] = %9.4f MPa\n",
	     Hsum[2][1]*ULENGTH_IN_A, stressum[2][1]*USTRESS_IN_MPA);
    fprintf (property,
	     "   H[3][3] = %8.4f A    stress[3][3] = %9.4f MPa\n\n",
	     Hsum[2][2]*ULENGTH_IN_A, stressum[2][2]*USTRESS_IN_MPA);
    
    fprintf (property, "   Avg. density = %f\\%f 1/A^3\\%f mol/m^3\n",
	     np/volume, np/volume/UVOLUME_IN_M3*1E-30,
	     np/volume/AVO/UVOLUME_IN_M3);
    fprintf (property,
	     "   Avg. pressure = %f\\%f MPa\n\n",
	     pressure, pressure*USTRESS_IN_MPA);
    
    fprintf (property, "   Avg. kinetic energy = "
	     "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n",
	     kinesum/np, kinesum/np*UENERGY_IN_J/EV_IN_J,
	     kinesum/np*UENERGY_IN_J/BOLZ, kinesum*UENERGY_IN_J*AVO/np);
    fprintf (property, "   Avg. potential energy = "
	     "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n",
	     potesum/np, potesum/np*UENERGY_IN_J/EV_IN_J,
	     potesum/np*UENERGY_IN_J/BOLZ, potesum*UENERGY_IN_J*AVO/np);
    potesum += kinesum; 
    fprintf (property, "   Avg. internal energy = "
	     "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n",
	     potesum/np, potesum/np*UENERGY_IN_J/EV_IN_J,
	     potesum/np*UENERGY_IN_J/BOLZ, potesum*UENERGY_IN_J*AVO/np);
    potesum += pressure * volume; 
    fprintf (property, "   Avg. enthalpy = "
	     "%.4f\\%.4f eV\\%.1f K\\%.1f J/mol\n\n",
	     potesum/np, potesum/np*UENERGY_IN_J/EV_IN_J,
	     potesum/np*UENERGY_IN_J/BOLZ, potesum*UENERGY_IN_J*AVO/np);
    
    fprintf (property,
	     "            | %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,
	     "            | %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,
	     "   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,
	     "            | %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,
	     "            | %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,
	     "            | %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, "   The finite stress bulk modulus "
	     "BT = (c11+2*c12+p)/3 = %f MPa\n", bulk_modulus*USTRESS_IN_MPA);
    fprintf (property, 
	     "   The line thermal expansion coefficient alpha = %f/K\n",
	     alpha);
    fprintf (property,
	     "   The heat capacity Cv = %f J/g/K = %f J/mol/K\n",
	     Cv*UENERGY_IN_J/(totalmass*UMASS_IN_KG*1000.),
	     Cv*UENERGY_IN_J/neff*AVO);
    fprintf (property,
	     "                     Cp = %f J/g/K = %f J/mol/K\n",  
	     Cp*UENERGY_IN_J/(totalmass*UMASS_IN_KG*1000.),
	     Cp*UENERGY_IN_J/neff*AVO);
    fprintf (property,
	     "\n   The self-diffusion coefficient D = %f 10^(-10) m^2/s\n",
	     (msd-msd_kstart)*ULENGTH_IN_M*ULENGTH_IN_M /
	     ((maxkb-kstart-1)*delta*UTIME_IN_SECOND)/6.*1E10);
    
    fclose (structure);
    
    /* save final configuration on file */
    write_config (fn_config_write);
	    
    /* free all allocated arrays */
    free (idx); free (list);
    free (bindx); free (binlist);
    free (binbindx); free (binbinlist);

    free (mb);
    free (tpd);
    free (tp[0]);
    free (tp);
    free (mybin);
    free (ep);
    free (mass);
    free (d);
    free (fs);
    free (f);
    free (v);
    free (s5);
    free (s4);
    free (s3);
    free (s2);
    free (s1);
    free (s);

    cancel_all_schedules();
    
    if ( CALC )
    {
	for (i=0; i<3; i++)
	{
	    fclose (heat_current[i]);
	    fclose (shear_stress[i]);
	}
	/* re-load heat currents and shear stresses from files and */
	/* calculate the average scalar thermal conductivity and */
	/* shear viscosity using FFT and the spectral method. */
	smile (SMILE_HEAT, property);
	smile (SMILE_STRESS, property);
    }
    
    cc = watch ("mission", 3);
    fprintf (property, "\n   This simulation took %s --"
	     "\n   or %f s/step, %f ms/step/particle.\n\n",
	     earth_time(cc), cc/maxkb, 1E3*cc/np/maxkb);
    fclose (property);
    return(0);
} /* end main() */