/******************************************/
/* Conjugate Gradient Relaxation of Argon */
/* Solid at Constant H or Constant Stress */
/*                                        */
/*                   Ver 1.1 Mar.20, 1999 */
/*               developed by Ju Li (MIT) */
/******************************************/

/*
  Version Notes:
  1.0  TTY interface;
       Woon's pair potential;
  1.1  Minor revisions;
*/

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

/** universal constants **/
#ifndef PI
#define PI 3.14159265358979323846
#endif
#define AVO 6.0221367E+23
#define BOLZ 1.380658E-23  /* in J/K */
#define HBAR 1.054589E-34  /* in J.s */
#define EV_IN_J 1.60217733E-19
#define HARTREE_IN_J 4.3597482E-18
#define BOHR_RADIUS_IN_A 0.529177249
/** argon LJ6-12 constants **/
#define SIGMA_IN_M (3.405E-10)
#define EPSILON_IN_J (119.8*BOLZ) 
#define ARGON_MASS_IN_KG (39.948*1E-3/AVO) 
/** fundamental reduced units **/
#define ULENGTH_IN_M SIGMA_IN_M
#define UENERGY_IN_J EPSILON_IN_J
#define UMASS_IN_KG  ARGON_MASS_IN_KG     
/** derived reduced units **/
#define ULENGTH_IN_A (ULENGTH_IN_M*1E10)
#define 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)
/* interaction cutoff */    
#define RCUT (2.5*SIGMA_IN_M/ULENGTH_IN_M)
#define MAX_STRLEN 100
#define BUILTIN_DEFECTS 3
#define MAX_SHADOW_CAPTURE 256
/* reference state */
#define BENCHMARK_NC 4
#define BENCHMARK_LATTICE_CONSTANT (5.3/ULENGTH_IN_A)
#define BENCHMARK_EPS (1E-10)
#define DEFAULT_ATOM_TP "Ar"
#define TP_STRLEN 2
#define MAX_M_EXPANSION (10)
/* numerical recipes declarations */
#define EPS (1.0E-10)
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY (1.0E-20)
#define CGOLD 0.3819660
#define ZEPS (1.0E-10)
#define TOL (2.0E-4)
#define ITMAX 2000
#define SIGN(a,b) ((b)>=0.0?fabs(a):-fabs(a))
#define MOV3(a,b,c,d,e,f) (a)=(d);(b)=(e);(c)=(f)
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d)
#define max(x,y) ((x)>(y)?(x):(y))
#define cube(x) ((x)*(x)*(x))

static int minimize_builtin_defect, minimize_user_defined_defect,
    LJ6_12, WOON, constant_H, constant_S, H0_enforced, write_config, np,
    npfcc, builtin_defect_index, save_screen_output, iseed,
    apply_shadow_positioning, use_default_iatom, iatom, nc[3], 
    linmin_ncom, builtin_defect_excess[] = {0, -1, 1};
static double volume, pote, hamiltonian, shadow_radius,
	perturbation_amplitude, minimization_tolerance, ep_benchmark,
	volumep_benchmark, *f, *s, *mass, *linmin_pcom, *linmin_xicom,
	*f1dim_xt, iatom_s[3], s_benchmark[12*cube(BENCHMARK_NC)],
	f_benchmark[12*cube(BENCHMARK_NC)], H_benchmark[3][3],
	stress_benchmark[3][3], H0[3][3], H[3][3], HI[3][3], eta[3][3],
	M[3][3], sout[3][3], stress[3][3];
static char defect_name[MAX_STRLEN], potential_function[MAX_STRLEN],
    fn_config_read[MAX_STRLEN], fn_config_write[MAX_STRLEN],
    fn_pdb_bef[MAX_STRLEN], fn_pdb_aft[MAX_STRLEN],
    fn_screen_output[MAX_STRLEN], **tp,
    *builtin_defect_names[]={"perfect_crystal", "vacancy", "interstitial"};
FILE *config, *screen_output; 
static double (*linmin_nrfunc)(double []);

static double value (double *s);
static int value_already_called; 
static void force (double *s, double *f);
static double new_H (double H0[3][3], double eta[3][3], double M[3][3],
		     double H[3][3], double HI[3][3]);
static void eta_to_M (double eta[3][3], double M[3][3]);
static double driving_force
(double H0[3][3], double eta[3][3], double c);
static double work_function
(double H0[3][3], double eta[3][3], double sout[3][3]);
static void shadow_positioning (char *fn_pdb);
static void bubble_sort (int nitems, int *index, double *d);
static double pair_potential_energy (int np, double H[3][3], double *s);
static void pair_potential_force
(int np, double H[3][3], double *s, double *f, double stress[3][3]);
static void assign_fcc_lattice
(int nc1, int nc2, int nc3, double *s, char **tp);
static double frandom();
static void mateqv (double A[3][3], double B[3][3]);
static void matran (double A[3][3], double B[3][3]);
static double matnorm (double A[3][3]);
static void matadd (double A[3][3], double B[3][3], double C[3][3]);
static void matadd_mul(double A[3][3], double CC, double B[3][3],
		       double C[3][3]);
static void matsub (double A[3][3], double B[3][3], double C[3][3]);
static double matinv (double A[3][3], double B[3][3]);
static void matmul (double A[3][3], double B[3][3], double C[3][3]);
static void polint (double xa[], double ya[], int n, double x,
		    double *y, double *dy);
static void print (const char *p, ...);
static void cg(double p[], int n, double ftol, int *iter,
	       double *fret, double (*func)(double []),
	       void (*dfunc)(double [], double []));
static void linmin (double p[], double xi[], int n,
                    double *fret, double (*func)(double []));
static double f1dim(double x);
static void mnbrak (double *ax, double *bx, double *cx, double *fa,
                    double *fb, double *fc, double (*func)(double));
static double brent (double ax, double bx, double cx, 
                     double (*f)(double), double tol, double *xmin);
static void nrerror (char error_text[]);

int main()
{
    int i, j, k, l;
    double tmp, cc, dd, ee;
    char buf[MAX_STRLEN] = {0};
    double pote_before_relaxation, volume_before_relaxation; 
    double cm_shift[3] = {0};
    
    printf ("Reading in instructions...\n");
    /* Read instructions from input, such as "ron" file */
    fscanf (stdin, "To relax (0) a built-in defect "
	    "(1) a user-defined configuration:\n%d\n\n", &i);
    minimize_builtin_defect = (i==0);
    minimize_user_defined_defect = (i==1);
    if ((!minimize_builtin_defect) && (!minimize_user_defined_defect))
    {
	printf("Error: (To relax) select options 0 or 1.\n");
	exit(1);
    }
    
    fscanf (stdin, "Name of the defect configuration (->*):\n%s\n\n",
	    defect_name);
	
    fscanf (stdin, "User-defined configuration file that could be read in "
	    "(default=*):\n%s\n\n", fn_config_read);
    if (!strcasecmp(fn_config_read, "default"))
	sprintf(fn_config_read, "%s", defect_name);
    
    fscanf (stdin, "Pair potential to be used (LJ6_12 or WOON):\n%s\n\n",
	    potential_function);
    LJ6_12 = !strcasecmp(potential_function, "LJ6_12");
    WOON   = !strcasecmp(potential_function, "WOON");
    if ((!LJ6_12) && (!WOON))
    {
	printf ("Error: Choose between LJ6_12 or WOON pair potential\n");
	exit(1);
    }

    fscanf (stdin, "constant_H or constant_S mode:\n%s\n\n", buf);
    constant_H = !strcasecmp(buf, "constant_H");
    constant_S = !strcasecmp(buf, "constant_S");
    if ((!constant_H) && (!constant_S))
    {
	printf ("Error: do not have \"%s\" mode, choose between "
		"constant_H and constant_S mode.\n", buf);
	exit(1);
    }

    fscanf (stdin, "Default reference unit cell (11,12,13,22,23,33) [A]:\n");
    fscanf (stdin, "%lf %lf %lf %lf %lf %lf\n\n", &H0[0][0], &H0[0][1],
	    &H0[0][2], &H0[1][1], &H0[1][2], &H0[2][2]);
    /* symmetric cell does not lose any freedom */
    H0[1][0] = H0[0][1]; H0[2][0] = H0[0][2]; H0[2][1] = H0[1][2];
    
    fscanf (stdin, "Default multiplications of the unit cell:\n");
    fscanf (stdin, "%d %d %d\n\n", &nc[0], &nc[1], &nc[2]);
    npfcc = 4 * nc[0] * nc[1] * nc[2];
    
    /* convert to reduced units */
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	    H0[i][j] *= nc[i] / ULENGTH_IN_A;

    fscanf (stdin, "In constant_H mode, use the above?\n%s\n\n", buf);
    H0_enforced = constant_H &&
	((buf[0]=='1') || (buf[0]=='y') || (buf[0]=='Y'));
    
    fscanf (stdin, "Default external stress (11,12,13,22,23,33) [MPa]:\n");
    fscanf (stdin, "%lf %lf %lf %lf %lf %lf\n\n", &sout[0][0], &sout[0][1],
	    &sout[0][2], &sout[1][1],  &sout[1][2], &sout[2][2]);
    sout[1][0] = sout[0][1];
    sout[2][0] = sout[0][2];
    sout[2][1] = sout[1][2];
    /* convert to reduced units */
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	    sout[i][j] /= USTRESS_IN_MPA;
    
    fscanf (stdin, "Index of the built-in defect that could be relaxed:\n");
    fscanf (stdin, "(0. perfect_crystal  1. vacancy  2. interstitial)\n%d\n\n",
	    &builtin_defect_index);
    if ( (builtin_defect_index<0) || (builtin_defect_index>=BUILTIN_DEFECTS) )
    {
	printf ("Error: select built-in defect index from 0 to %d.\n",
		BUILTIN_DEFECTS-1);
	exit(1);
    }
    
    fscanf (stdin, "Amplitude of random position perturbations "
	    "(in lattice constant):\n%lf\n\n", &perturbation_amplitude);

    fscanf (stdin, "Random number generator seed (0->time seed):\n");
    fscanf (stdin, "%d\n", &iseed); 
    if (iseed == 0) iseed = time(NULL);
    srandom ((unsigned)iseed);
    
    fscanf (stdin, "Minimization tolerance:\n%lf\n\n",
	    &minimization_tolerance);

    fscanf (stdin, "Save to .pdb files the particle positions around "
	    "iatom (-2=don't,\n-1=use_default) within R[A], before (default "
	    "= *.u) and after\nrelaxation (default = *.r):\n");
    fscanf (stdin, "%d %lf %s %s\n\n",
	    &iatom, &shadow_radius, fn_pdb_bef, fn_pdb_aft);
    
    apply_shadow_positioning = (iatom != -2); 
    use_default_iatom = (iatom == -1);
    shadow_radius /= ULENGTH_IN_A;
    if (!strcasecmp (fn_pdb_bef, "default"))
	sprintf (fn_pdb_bef, "%s.u", defect_name);
    if (!strcasecmp (fn_pdb_aft, "default"))
	sprintf (fn_pdb_aft, "%s.r", defect_name);
    
    fscanf (stdin, "Save screen output to file (default=*.out, "
	    "NULL=don't):\n%s\n\n", fn_screen_output);
    if ( save_screen_output=strcasecmp(fn_screen_output, "NULL") )
	if ( !strcasecmp(fn_screen_output, "default") )
	    sprintf (fn_screen_output, "%s.out", defect_name);
    
    fscanf (stdin, "Save optimized configuration on file "
	    "(default=*, NULL=don't):\n%s\n", fn_config_write);
    if ( write_config=strcasecmp(fn_config_write, "NULL") )
	if ( !strcasecmp(fn_config_write, "default") )
	    sprintf (fn_config_write, "%s", defect_name);
    
    printf ("Instructions read complete.\n\n");
    
    if (save_screen_output) screen_output=fopen(fn_screen_output, "w+");
    
    print("=======================  Mission Statement  "
	  "=======================\n");
    print("Free energy will be minimized for defect \"%s\"\n", defect_name);
    print("using %s pair potential with RCUT = %f = %f A\n",
	  potential_function, RCUT, RCUT*ULENGTH_IN_A);
    print("under \"%s\" mode.\n\n",
	  constant_H?"constant_H":"constant_S");
    
    if (constant_S)
    {
	print("            | %.5f %.5f %.5f |\n",
	      sout[0][0]*USTRESS_IN_MPA,
	      sout[0][1]*USTRESS_IN_MPA,
	      sout[0][2]*USTRESS_IN_MPA);	  
	print("sout[MPa] = | %.5f %.5f %.5f |\n",
	      sout[1][0]*USTRESS_IN_MPA,
	      sout[1][1]*USTRESS_IN_MPA,
	      sout[1][2]*USTRESS_IN_MPA);
	print("            | %.5f %.5f %.5f |\n\n",
	      sout[2][0]*USTRESS_IN_MPA,
	      sout[2][1]*USTRESS_IN_MPA,
	      sout[2][2]*USTRESS_IN_MPA);
    }
    
    /* benchmark system */
    print("The defect will be benchmarked against perfect crystal\n");
    print("(BENCHMARK_NP = %d) at P = 0:\n", 4*cube(BENCHMARK_NC));
    assign_fcc_lattice (BENCHMARK_NC, BENCHMARK_NC, BENCHMARK_NC,
			s_benchmark, NULL);
    
    H_benchmark[0][1] = H_benchmark[1][0] = 0.;
    H_benchmark[0][2] = H_benchmark[2][0] = 0.;
    H_benchmark[1][2] = H_benchmark[2][1] = 0.;
    cc = BENCHMARK_LATTICE_CONSTANT * 0.5;
    dd = BENCHMARK_LATTICE_CONSTANT * 2.0;
    for (i=0; ;i++)
    { /* binary search for zero pressure */
	tmp = (cc + dd) / 2.;
	H_benchmark[0][0] = BENCHMARK_NC * tmp;
	H_benchmark[1][1] = BENCHMARK_NC * tmp;
	H_benchmark[2][2] = BENCHMARK_NC * tmp;
	pair_potential_force (4*cube(BENCHMARK_NC), H_benchmark,
			      s_benchmark, f_benchmark, stress);
	if (fabs(stress[0][0])/USTRESS_IN_MPA<BENCHMARK_EPS) break;
	else if (stress[0][0]>0.0) cc = tmp;
	else dd = tmp;
	if (i >= 100)
	{
	    print("Error: main(): binary search for zero pressure fails.\n");
	    exit(1);
	}
    }

    ep_benchmark = pair_potential_energy
	(4*cube(BENCHMARK_NC),H_benchmark,s_benchmark) /4/cube(BENCHMARK_NC);
    volumep_benchmark = cube(tmp) / 4;
    print("benchmarked lattice constant = %f = %f A,\n",
	  tmp, tmp*ULENGTH_IN_A);
    print("benchmarked volume per particle = %f = %f A^3,\n", 
	  cube(tmp)/4, cube(tmp/ULENGTH_IN_A)/4);
    print("benchmarked energy per particle = %f = %f eV,\n",
	  ep_benchmark, ep_benchmark*UENERGY_IN_J/EV_IN_J);
    print("benchmark done.\n\n");

    /* determine the maximum number of particles */
    if (minimize_builtin_defect)
	np = (builtin_defect_excess[builtin_defect_index]>0)?
	    npfcc+builtin_defect_excess[builtin_defect_index]:npfcc;
    else
    {  /* user-defined configuration */
	if ((config=fopen(fn_config_read, "r"))==NULL)
	{
	    printf("Error: cannot open file \"%s\"\n.",
		   fn_config_read);
	    exit(1);
	}
	fscanf (config, "Number of particles = %d\n\n", &np);
    }
    
    print("Allocating memories for maximum of %d particles.\n",np);
    s = (double *) malloc((3*np+6)*sizeof(double));
    f = (double *) malloc((3*np+6)*sizeof(double));
    /* keep mass[] for isotope defects, but don't do anything with it */
    mass = (double *) malloc(np*sizeof(double));
    /* name space for atom type */
    tp   = (char **) malloc (np*sizeof(char *));
    tp[0] = (char *) malloc (np*(TP_STRLEN+1)*sizeof(char));
    for (i=1; i<np; i++)
	tp[i] = tp[i-1] + TP_STRLEN + 1;
    
    /* assign reduced coordinates */
    if (minimize_builtin_defect)
    {
	print("The initial configuration is created from scratch.\n");
	assign_fcc_lattice (nc[0], nc[1], nc[2], s, tp);
	np = npfcc; 
    }
    else 
    {
	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",  &nc[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",  &nc[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",  &nc[2]);
	for (i=0; i<3; i++)
	    for (j=0; j<3; j++)
		H[i][j] /= ULENGTH_IN_A; 
	volume = matinv(H, HI);
	fscanf (config, "TP Mass(amu)     sx        sy        sz"
		"     sx1(1/ns)  sy1(1/ns)  sz1(1/ns)\n");
	for (i=0; i<3*np; i+=3)
	    fscanf (config, "%s %lf %lf %lf %lf %lf %lf %lf\n",
		    tp[i/3], &mass[i/3], &s[i], &s[i+1], &s[i+2],
		    &tmp, &tmp, &tmp);
	print("Loading initial configuration from \"%s\" complete.\n",
	      fn_config_read);
	print("Found %d particles.\n\n", np);
    }
    
    /* H matrix and accessories */
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	{
	    if (minimize_builtin_defect || H0_enforced)
		H[i][j] = H0[i][j];
	    else H0[i][j] = H[i][j];
	    M[i][j] = 1.;
	    eta[i][j] = 0.;
	}
    volume = new_H (H0, eta, M, H, HI);
    
    print("               | %.5f %.5f %.5f |\n",
	  H[0][0]*ULENGTH_IN_A,
	  H[0][1]*ULENGTH_IN_A,
	  H[0][2]*ULENGTH_IN_A);	  
    print("Initial H[A] = | %.5f %.5f %.5f |\n",
	  H[1][0]*ULENGTH_IN_A,
	  H[1][1]*ULENGTH_IN_A,
	  H[1][2]*ULENGTH_IN_A);
    print("               | %.5f %.5f %.5f |\n\n",
	  H[2][0]*ULENGTH_IN_A,
	  H[2][1]*ULENGTH_IN_A,
	  H[2][2]*ULENGTH_IN_A);

    print("Total supercell volume = %f = %f A^3.\n\n", volume,
	  volume*cube(ULENGTH_IN_A));

    /* creating defect configuration */
    if (minimize_builtin_defect)
    {
	switch(builtin_defect_index)
	{
	case 0: /* perfect crystal */
	    iatom = 0;
	    iatom_s[0] = s[3*iatom];
	    iatom_s[1] = s[3*iatom+1];
	    iatom_s[2] = s[3*iatom+2];
	    break;
	case 1: /* vacancy */
	    iatom = 0;
	    iatom_s[0] = s[3*iatom];
	    iatom_s[1] = s[3*iatom+1];
	    iatom_s[2] = s[3*iatom+2];
	    s[0] = s[3*np-3];
	    s[1] = s[3*np-2];
	    s[2] = s[3*np-1];
	    np--;
	    break;
	case 2: /* interstitial */
	    iatom = 0;
	    iatom_s[0] = s[3*iatom];
	    iatom_s[1] = s[3*iatom+1];
	    iatom_s[2] = s[3*iatom+2];
	    s[3*np]   = s[0] + 0.25/nc[0];
	    s[3*np+1] = s[1] + 0.25/nc[1];
	    s[3*np+2] = s[2] + 0.25/nc[2];
	    sprintf(tp[np++], "%2s", DEFAULT_ATOM_TP);
	    break;
	default:
	    print("Error: builtin_defect[%d] does not exist.\n",
		  builtin_defect_index);
	    exit(1);
	}
    }
    
    print("Defect \"%s\" generated successfully on the computer.\n\n",
	  defect_name);

    if (apply_shadow_positioning)
    {
	print("Particles within distance %.3f A of original particle "
	      "%d would\n", shadow_radius*ULENGTH_IN_A, iatom);
	print("be recorded in \"%s\" before minimization and\n", fn_pdb_bef);
	print("\"%s\" after minimization, in .pdb format\n", fn_pdb_aft);
	print("that can be viewed by Rasmol (freeware).\n\n");
    }
    
    print("******************** Before Minimization ********************\n");
    pote_before_relaxation = pair_potential_energy (np, H, s); 
    print("The total energy = %f eV, energy/atom = %f eV,\n",
	  pote_before_relaxation*UENERGY_IN_J/EV_IN_J,
	  pote_before_relaxation/np*UENERGY_IN_J/EV_IN_J);
    print("so the defect formation energy (unrelaxed)\n");
    print("= %f - %d x %f = %f eV.\n",
	  pote_before_relaxation*UENERGY_IN_J/EV_IN_J, np,
	  ep_benchmark*UENERGY_IN_J/EV_IN_J,
	  (pote_before_relaxation-ep_benchmark*np)*UENERGY_IN_J/EV_IN_J);
    print("(free energy reference is perfect crystal at zero stress)\n\n");
    
    volume_before_relaxation = volume;
    print("The total volume = %.3f A^3, volume/atom = %f A^3,\n", 
	  volume*ULENGTH_IN_A*ULENGTH_IN_A*ULENGTH_IN_A,
	  volume/np*ULENGTH_IN_A*ULENGTH_IN_A*ULENGTH_IN_A);
    print("so the excess volume (unrelaxed)\n");
    print("= %f - %d x %f = %f A^3.\n", volume*cube(ULENGTH_IN_A),
	  np, volumep_benchmark*cube(ULENGTH_IN_A),
	  (volume-np*volumep_benchmark)*cube(ULENGTH_IN_A));
    print("(volume reference is perfect crystal at zero stress)\n\n");
    
    pair_potential_force(np, H, s, f, stress);
    print("                      | %.5f %.5f %.5f |\n",
	  stress[0][0]*USTRESS_IN_MPA,
	  stress[0][1]*USTRESS_IN_MPA,
	  stress[0][2]*USTRESS_IN_MPA);	  
    print("Initial stress[MPa] = | %.5f %.5f %.5f |\n",
	  stress[1][0]*USTRESS_IN_MPA,
	  stress[1][1]*USTRESS_IN_MPA,
	  stress[1][2]*USTRESS_IN_MPA);
    print("                      | %.5f %.5f %.5f |\n\n",
	  stress[2][0]*USTRESS_IN_MPA,
	  stress[2][1]*USTRESS_IN_MPA,
	  stress[2][2]*USTRESS_IN_MPA);
    
    if (apply_shadow_positioning) shadow_positioning(fn_pdb_bef);
    print("*************************************************************\n\n");
    
    print("Random number seed = %d, using which we give each particle\n",
	  iseed);
    cc = perturbation_amplitude*cbrt(volume/np*4)/2.; 
    print("a random kick of (%f, %f) A in three directions.\n\n",
	  -cc*ULENGTH_IN_A, cc*ULENGTH_IN_A);
    for (i=0; i<3*np; i++)
    {
	cc = perturbation_amplitude*cbrt(1./np*4)*(frandom()-0.5);
	s[i] += cc;
	cm_shift[i%3] += cc;
    }
    for (i=0; i<3*np; i++)
    {
	s[i] -= cm_shift[i%3]/np;
	/* make sure that s[i] is between 0 and 1 */
	s[i] -= floor(s[i]);
    }
    
    print("Minimization tolerance = %.10f.\n\n", minimization_tolerance);
    
    if (write_config)
	print("The relaxed configuration would be saved on \"%s\".\n",
	      fn_config_write);

    if (save_screen_output)
	print("The screen output would be saved on \"%s\".\n",
	      fn_screen_output);
    
    print("==================== End of Mission Statement "
	  "=====================\n\n");

    print("Minimization starts...\n\n");
    
    value_already_called = 0;
	
    if (constant_H)
	cg (s, 3*np, minimization_tolerance, &i, &cc, value, force);
    else
    { /* Lagrangian strains are the additional degrees of freedom */ 
	s[3*np]   = eta[0][0];
	s[3*np+1] = eta[1][1];
	s[3*np+2] = eta[2][2];
	s[3*np+3] = eta[0][1];
	s[3*np+4] = eta[0][2];
	s[3*np+5] = eta[1][2];
	cg (s, 3*np+6, minimization_tolerance, &i, &cc, value, force);
	eta[0][0] = s[3*np];
	eta[1][1] = s[3*np+1];
	eta[2][2] = s[3*np+2];
	eta[0][1] = s[3*np+3];
	eta[0][2] = s[3*np+4];
	eta[1][2] = s[3*np+5];
	eta[1][0] = eta[0][1];
	eta[2][0] = eta[0][2];
	eta[2][1] = eta[1][2];
	volume = new_H (H0, eta, M, H, HI);
    }
    
    print("\nMinimization converged after %d iterations.\n\n",i);
    
    print("******************** After Minimization ********************\n");
    print("The total energy = %f eV, energy/atom = %f eV,\n",
	  cc*UENERGY_IN_J/EV_IN_J, cc*UENERGY_IN_J/EV_IN_J/np);
    print("so the defect formation energy (relaxed)\n");
    print("= %f - %d x %f = %f eV.\n", cc*UENERGY_IN_J/EV_IN_J, np,
	  ep_benchmark*UENERGY_IN_J/EV_IN_J,
	  (cc-ep_benchmark*np)*UENERGY_IN_J/EV_IN_J);
    print("(free energy reference is perfect crystal at zero stress)\n");
    print("and it has increased by %f eV during the relaxation.\n\n",
	  (cc-pote_before_relaxation)*UENERGY_IN_J/EV_IN_J);

    print("             | %.5f %.5f %.5f |\n",
	  H[0][0]*ULENGTH_IN_A,
	  H[0][1]*ULENGTH_IN_A,
	  H[0][2]*ULENGTH_IN_A);	  
    print("Final H[A] = | %.5f %.5f %.5f |\n",
	  H[1][0]*ULENGTH_IN_A,
	  H[1][1]*ULENGTH_IN_A,
	  H[1][2]*ULENGTH_IN_A);
    print("             | %.5f %.5f %.5f |\n\n",
	  H[2][0]*ULENGTH_IN_A,
	  H[2][1]*ULENGTH_IN_A,
	  H[2][2]*ULENGTH_IN_A);
    
    print("The total volume = %.3f A^3, volume/atom = %f A^3,\n", 
	  volume*cube(ULENGTH_IN_A), volume/np*cube(ULENGTH_IN_A));
    print("so the excess volume (relaxed)\n");
    print("= %f - %d x %f = %f A^3.\n", volume*cube(ULENGTH_IN_A),
	  np, volumep_benchmark*cube(ULENGTH_IN_A),
	  (volume-np*volumep_benchmark)*cube(ULENGTH_IN_A));
    print("(volume reference is perfect crystal at zero stress)\n");     
    print("and it has increased by %f A^3 during the run.\n\n",
	  (volume-volume_before_relaxation)*cube(ULENGTH_IN_A));
    
    pair_potential_force(np, H, s, f, stress);
    print("                    | %.5f %.5f %.5f |\n",
	  stress[0][0]*USTRESS_IN_MPA,
	  stress[0][1]*USTRESS_IN_MPA,
	  stress[0][2]*USTRESS_IN_MPA);	  
    print("Final stress[MPa] = | %.5f %.5f %.5f |\n",
	  stress[1][0]*USTRESS_IN_MPA,
	  stress[1][1]*USTRESS_IN_MPA,
	  stress[1][2]*USTRESS_IN_MPA);
    print("                    | %.5f %.5f %.5f |\n\n",
	  stress[2][0]*USTRESS_IN_MPA,
	  stress[2][1]*USTRESS_IN_MPA,
	  stress[2][2]*USTRESS_IN_MPA);
    
    if (apply_shadow_positioning) shadow_positioning(fn_pdb_aft);
    print("************************************************************\n\n");

    if (write_config)
    { /* save system configuration in readable ASCII format */
	print("Saving relaxed configuration on file \"%s\".\n",
	      fn_config_write);
	config = fopen (fn_config_write, "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], s[i], s[i+1], s[i+2], 0., 0., 0.);
	fclose (config);
    }
    
    /* free all allocated memories */
    free (tp[0]);
    free (tp);
    free (f);
    free (s);
    free (mass);
    
    if (save_screen_output)
    {
	print("screen output saved on file \"%s\".\n", fn_screen_output);
	close(screen_output);
    }

    printf("Mission terminates successfully.\n\n");
    return(0);
} /* end main() */


double value (double *s)
{
    double free_energy;
    if (constant_H) free_energy = pair_potential_energy (np, H, s);
    else if (constant_S)
    {
	eta[0][0] = s[3*np];
	eta[1][1] = s[3*np+1];
	eta[2][2] = s[3*np+2];
	eta[0][1] = s[3*np+3];
	eta[0][2] = s[3*np+4];
	eta[1][2] = s[3*np+5];
	eta[1][0] = eta[0][1];
	eta[2][0] = eta[0][2];
	eta[2][1] = eta[1][2];
	volume = new_H (H0, eta, M, H, HI);
    	free_energy = pair_potential_energy (np, H, s)
	    - work_function (H0, eta, sout);
    }
    print("value = %f\n", free_energy*UENERGY_IN_J/EV_IN_J);
    value_already_called = 1;
    return (free_energy);
} /* end value() */


void force (double *s, double *f)
{
    int i;
    double fsx, fsy, fsz;
    double MI[3][3], tmp[3][3], tmpp[3][3];
    /* should call value() before calling force() */
    if (!value_already_called) value(s);
    else value_already_called = 0;
    /* now that H, M, volume etc. should be ready */
    print("\nforce evaluation...\n");
    pair_potential_force (np, H, s, f, stress);
    for (i=0;i<3*np;i+=3)
    {
	fsx = H[0][0]*f[i]+H[1][0]*f[i+1]+H[2][0]*f[i+2];
	fsy = H[0][1]*f[i]+H[1][1]*f[i+1]+H[2][1]*f[i+2];
	fsz = H[0][2]*f[i]+H[1][2]*f[i+1]+H[2][2]*f[i+2];
	f[i]   = fsx;
	f[i+1] = fsy;
	f[i+2] = fsz;
    }
    if (constant_S)
    {
	matinv (M, MI);
	matadd (sout, stress, tmp);
	matmul (MI, tmp, tmpp);
	matmul (tmpp, MI, tmp);
	print("stress discrepancy = %f MPa\n",
	      matnorm(tmp)*USTRESS_IN_MPA);
	f[3*np]   = tmp[0][0] * volume;
	f[3*np+1] = tmp[1][1] * volume;
	f[3*np+2] = tmp[2][2] * volume;
	f[3*np+3] = tmp[0][1] * volume * 2.;
	f[3*np+4] = tmp[0][2] * volume * 2.;
	f[3*np+5] = tmp[1][2] * volume * 2.;
	print("f = %f %f %f %f %f %f\n", f[3*np], f[3*np+1], f[3*np+2],
	      f[3*np+3], f[3*np+4], f[3*np+5]);
    }
    return;
} /* end force() */


/* update M, H, HI, volume from reference state H0 and Lagrangian strain eta */
double new_H (double H0[3][3], double eta[3][3], double M[3][3],
	      double H[3][3], double HI[3][3])
{
    double cc;
    eta_to_M (eta, M);
    matmul (H0, M, H);
    cc = matinv (H, HI);
    return (cc);
} /* end new_H() */


/* calculate affine transformation matrix M from ETA */
void eta_to_M (double eta[3][3], double M[3][3])
{
    int i, j;
    double tmp[3][3], tmpp[3][3];
    const double coeff[] =
    {1, 1, -1./2, 1./2, -5./8,
     7./8, -21./16, 33./16, -429./128, 715./128,
     -2431./256, 4199./256, -29393./1024, 52003./1024, -185725./2048,
     334305./2048, -9694845./32768, 17678835./32768, -64822395./65536,
     119409675./65536, -883631595./262144, 1641030105./262144};
    /* Taylor expansion coefficients of M = sqrt(1+2*ETA) */
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	    tmp[i][j] = (i==j)?1:0;
    mateqv (tmp, M);
    for (i=1; i<=MAX_M_EXPANSION; i++)
    {
	matmul (tmp, eta, tmpp);
	mateqv (tmpp, tmp);
	matadd_mul (M, coeff[i], tmp, M);
    }
    return;
} /* end eta_to_M() */


/* how much work is done by changing c*eta (c from 0 to 1) */
double driving_force (double H0[3][3], double eta[3][3], double c)
{
    int i, j;
    double etac[3][3], Mc[3][3], MIc[3][3], volumec; 
    double tmp[3][3], tmpp[3][3], cc;
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	    etac[i][j] = c * eta[i][j];
    eta_to_M (etac, Mc);
    matmul (H0, Mc, tmp);
    volumec = matinv (tmp, tmpp);
    matinv (Mc, MIc);
    matmul (MIc, sout, tmpp);
    matmul (tmpp, MIc, tmp);
    for (cc=0.,i=0; i<3; i++)
	for (j=0; j<3; j++)
	    cc += volumec * tmp[i][j] * eta[i][j];
    return (cc);
} /* end driving_force() */


/* calculate the work done by external stress in symmetric */
/* deformation space by straight path Romberg integration: */
#define JMAX 25
#define KMAX 8
#define REPS (1.0E-10)
double work_function
(double H0[3][3], double eta[3][3], double sout[3][3])
{
    int i, j, k; 
    double  h[JMAX+2], s[JMAX+2], del, sum, error, cc, t;
    /* if the strain is too large, prevent collapse */
    /* if (matnorm(eta)>0.5) return(-matnorm(eta)*matnorm(eta)); */
    if (matnorm(eta) > 0.5) return (0.);
    h[1] = 1.;
    s[1] = 0.5 * (driving_force(H0, eta, 0.)+
		  driving_force(H0, eta, 1.));
    for (j=2, i=1; j<=JMAX; j++, i<<=1)
    { /* integration by trapzoidal rule with progressively decreasing h */
	h[j] = 0.25 * h[j-1];
	del = 1. / i;
	t = 0.5*del;
	sum = 0.;
	for (k=1; k<=i; k++)
	{
	    sum += driving_force(H0, eta, t);
	    t += del;
	}
	s[j] = 0.5*(s[j-1]+sum/i);
	if (j >= KMAX)
	{
	    polint(&h[j-KMAX],&s[j-KMAX],KMAX,0.,&cc,&error);
	    /* print("%f \n", fabs(error)/fabs(cc)); */
	    if (fabs(error) <= REPS*fabs(cc)) return(cc);
	}
    }
    nrerror("work_function(): can not converge in Romberg integration");
    return(0.);
} /* end work_function() */
#undef REPS
#undef KMAX
#undef JMAX


/* print out particle coordinates around */
/* iatom_s[] and save to a .pdb file.    */
#define SHIFT_TINY (0.0002)
void shadow_positioning (char *fn_pdb)
{
    int i, j, k, ip;
    double dsx, dsy, dsz, rx, ry, rz, r2;
    double dx[MAX_SHADOW_CAPTURE], dy[MAX_SHADOW_CAPTURE],
	dz[MAX_SHADOW_CAPTURE], d[MAX_SHADOW_CAPTURE];
    int index[MAX_SHADOW_CAPTURE], particle_index[MAX_SHADOW_CAPTURE];
    FILE *pdb;

    for (ip=0, i=0; i<3*np; i+=3)
    {
	dsx = s[i]   - iatom_s[0];
	dsy = s[i+1] - iatom_s[1];
	dsz = s[i+2] - iatom_s[2];
	if (dsx>0.5) dsx--;
	else if (dsx<-0.5) dsx++;
	if (dsy>0.5) dsy--;
	else if (dsy<-0.5) dsy++;
	if (dsz>0.5) dsz--;
	else if (dsz<-0.5) dsz++;
	rx = H[0][0]*dsx + H[1][0]*dsy + H[2][0]*dsz;
	ry = H[0][1]*dsx + H[1][1]*dsy + H[2][1]*dsz;
	rz = H[0][2]*dsx + H[1][2]*dsy + H[2][2]*dsz;
	r2 = rx*rx+ry*ry+rz*rz;
	if (r2<shadow_radius*shadow_radius)
	{
	    if (ip>=MAX_SHADOW_CAPTURE)
	    {
		print("shadow_positioning(): out of bound, "
		      "change MAX_SHADOW_CAPTURE\n");
		exit(1);
	    }
	    dx[ip] = rx * ULENGTH_IN_A;
	    dy[ip] = ry * ULENGTH_IN_A;
	    dz[ip] = rz * ULENGTH_IN_A;
	    d[ip] = sqrt(r2) * ULENGTH_IN_A;
	    particle_index[ip] = i/3;
	    ip++;
	}
    }
    /* sort the radii in ascending order */
    bubble_sort (ip, index, d);
    print("-------- shadow positioning around tag particle %d --------\n",
	  iatom);
    print(" neigh  index     D[A]     Dx[A]     Dy[A]     Dz[A]\n");
    print("----------------------------------------------------------\n");
    pdb = fopen (fn_pdb, "w+");
    fprintf (pdb, "HEADER    PDB file generated by Relax <Ju Li>");
    fprintf (pdb, "HEADER    R = %f A around atom %d (%2s).\n",
	     shadow_radius*ULENGTH_IN_A, iatom, tp[iatom]);
    for (i=0; i<ip; i++)
    {
	k = index[i]; 
	fprintf (pdb, "ATOM  %5d %2s          %2s    "
		 "%8.3f%8.3f%8.3f\n", i+1, tp[particle_index[k]],
		 " 1", dx[k], dy[k], dz[k]);
	j = (d[index[0]]<=SHIFT_TINY)? i : i+1;
	print("%5d  %5d  %8.5f  %8.5f  %8.5f  %8.5f\n",
	      j, particle_index[k], d[k], dx[k], dy[k], dz[k]);
    }
    print("----------------------------------------------------------\n");
    print("(above is saved on \"%s\" in .pdb format)\n", fn_pdb);
    fclose(pdb);
    return;
} /* end shadow_positioning() */
#undef SHIFT_TINY


/* simplest realization of ascending sorting */
void bubble_sort (int nitems, int *index, double *d)
{
    int i,j,k;
    for (i=0; i<nitems; i++) index[i] = i;
    for (i=1; i<nitems; i++)
    	for (j=i; j>0; j--)
	{
	    if (d[index[j]]<d[index[j-1]])
	    {
		k = index[j-1];
		index[j-1] = index[j];
		index[j] = k;
	    }
	    else break;
	}
    return;
} /* bubble_sort */


/* 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 */
double pair_potential_energy (int np, double H[3][3], double *s)
{
    int i, j;
    double sxij, syij, szij, rx, ry, rz, r, r2, r6, r12, pote, a, b;
    double res = LJ6_12?LJ_RESIDUE:WOON_RESIDUE;
    double ref = LJ6_12?LJ_RESIDUE_FORCE:WOON_RESIDUE_FORCE;
    pote = 0.;
    for (i=0; i<3*np; i+=3)
	for (j=i+3; j<3*np; j+=3)
	{
	    sxij = s[i]   - s[j];
	    syij = s[i+1] - s[j+1];
	    szij = s[i+2] - s[j+2];
	    if (sxij>0.5) sxij--;
	    else if (sxij<-0.5) sxij++;
	    if (syij>0.5) syij--;
	    else if (syij<-0.5) syij++;
	    if (szij>0.5) szij--;
	    else if (szij<-0.5) szij++;
	    rx = H[0][0]*sxij + H[1][0]*syij + H[2][0]*szij;
	    ry = H[0][1]*sxij + H[1][1]*syij + H[2][1]*szij;
	    rz = H[0][2]*sxij + H[1][2]*syij + H[2][2]*szij;
	    r2 = rx*rx+ry*ry+rz*rz;
 	    if ( r2 < RCUT*RCUT )
	    {
		r = sqrt(r2);
		if (LJ6_12)
		{
		    r2  = 1./r2;
		    r6  = r2*r2*r2;
		    r12 = r6*r6;
		    pote += 4*(r12-r6) - res + ref*(r-RCUT);
		}
		else
		{
		    a = WOON_EPSILON * WOON_Q *
			exp(-WOON_K1 * (r-WOON_RE));
		    b = WOON_EPSILON * (WOON_Q+1) *
			exp(-WOON_K2 * (r-WOON_RE)); 
		    pote += a - b - res + ref*(r-RCUT);
		}
	    }
	}
    return (pote);
} /* end pair_potential_energy() */


void pair_potential_force
(int np, double H[3][3], double *s, double *f, double stress[3][3])
{
    int i, j, k;
    double volume, HI[3][3], sxij, syij, szij, rx, ry, rz;
    double r, r2, r6, r12, vij, ff, ffx, ffy, ffz, cc, a, b;
    double ref = LJ6_12?LJ_RESIDUE_FORCE:WOON_RESIDUE_FORCE;

    volume = matinv (H, HI);
    for (i=0; i<3*np; i++) f[i] = 0.;
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	    stress[i][j] = 0.;
    
    for (i=0; i<3*np; i+=3)
	for (j=i+3; j<3*np; j+=3)
	{
	    sxij = s[i]   - s[j];
	    syij = s[i+1] - s[j+1];
	    szij = s[i+2] - s[j+2];
	    if (sxij>0.5) sxij--;
	    else if (sxij<-0.5) sxij++;
	    if (syij>0.5) syij--;
	    else if (syij<-0.5) syij++;
	    if (szij>0.5) szij--;
	    else if (szij<-0.5) szij++;

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

	    if ( r2 < RCUT*RCUT )
	    {
		r = sqrt(r2);
		if (LJ6_12)
		{
		    r2  = 1./r2;
		    r6  = r2*r2*r2;
		    r12 = r6*r6;
		    ff  = (48.*r12-24.*r6)*r2 - ref/r;
		}
		else
		{
		    a = WOON_EPSILON * WOON_Q *
			exp(-WOON_K1 * (r-WOON_RE));
		    b = WOON_EPSILON * (WOON_Q+1) *
			exp(-WOON_K2 * (r-WOON_RE)); 
		    ff = (WOON_K1*a - WOON_K2*b - ref) / r;
		}
		ffx = ff*rx;
		ffy = ff*ry;
		ffz = ff*rz;
		f[i]   += ffx;
		f[i+1] += ffy;
		f[i+2] += ffz;
		f[j]   -= ffx;
		f[j+1] -= ffy;
		f[j+2] -= ffz;
		/* calculate the stress tensor */
		stress[0][0] += ffx*rx;
		stress[1][1] += ffy*ry;
		stress[2][2] += ffz*rz;
		stress[0][1] += (ffx*ry+ffy*rx)/2.;
		stress[0][2] += (ffx*rz+ffz*rx)/2.;
		stress[1][2] += (ffy*rz+ffz*ry)/2.;
	    }
	}
    stress[0][0] /= volume;
    stress[1][1] /= volume;
    stress[2][2] /= volume;
    stress[0][1] /= volume;
    stress[0][2] /= volume;
    stress[1][2] /= volume;
    stress[1][0] = stress[0][1];
    stress[2][0] = stress[0][2];
    stress[2][1] = stress[1][2];

    return;
} /* end pair_potential_force() */


/* assign fcc reduced coordinates according to nc[] */
void assign_fcc_lattice (int nc1, int nc2, int nc3, double *s, char **tp)
{
    int i, j, k, l, ip;
    double x[12];
    x[0] = 0.;  x[1] = 0.;  x[2] = 0.;
    x[3] = 0.5; x[4] = 0.5; x[5] = 0.;
    x[6] = 0.;  x[7] = 0.5; x[8] = 0.5;
    x[9] = 0.5; x[10] = 0.; x[11] = 0.5;
    for (ip=i=0; i<nc1; i++)
	for (j=0; j<nc2; j++)
	    for (k=0; k<nc3; k++,ip+=12)
		for (l=0; l<12; l+=3)
		{
		    s[l+ip]   = (x[l]  +i+0.25) / nc1;
		    s[l+1+ip] = (x[l+1]+j+0.25) / nc2;
		    s[l+2+ip] = (x[l+2]+k+0.25) / nc3;
		    if (tp != NULL)
			strncpy (tp[(ip+l)/3], DEFAULT_ATOM_TP, TP_STRLEN);
		}
    return;
} /* end assign_fcc_lattice() */


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


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


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


double matnorm(double A[3][3])
{
    int i,j;
    double result;
    result = 0.;
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    result += A[i][j]*A[i][j];
    result = sqrt(result);
    return(result);
} /* end matnorm() */


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


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


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


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


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


void cg(double p[], int n, double ftol, int *iter, double *fret,
	double (*func)(double []), void (*dfunc)(double [], double []))
{
    int j,its;
    double gg,gam,fp,dgg;
    double *g,*h,*xi;
    /* initialize cg() */
    g = (double *)malloc(n*sizeof(double));
    /* cg() allocates forces */
    h = (double *)malloc(n*sizeof(double));
    /* cg() allocates search vectors */
    xi = (double *)malloc(n*sizeof(double));
    /* cg() allocates search directions */
    /* initialize linmin() */
    f1dim_xt = (double *)malloc(n*sizeof(double));
    fp=(*func)(p);  /* value */
    (*dfunc)(p,xi);  /* force */
    for (j=0;j<n;j++)
    {
	g[j] = -xi[j];
	xi[j]=h[j]=g[j];
    }
    for (its=1;its<=ITMAX;its++)
    {
	*iter=its;
	linmin(p,xi,n,fret,func);
	if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) 
	    goto fine;
	fp=(*func)(p);
	(*dfunc)(p,xi);
	dgg=gg=0.0;
	for (j=0;j<n;j++)
	{
	    gg += g[j]*g[j];
	    dgg += (xi[j]+g[j])*xi[j];
	}
	if (gg == 0.0) goto fine;
	gam=dgg/gg;
	for (j=0;j<n;j++)
	{
	    g[j] = -xi[j];
	    xi[j]=h[j]=g[j]+gam*h[j];
	}
    }
error:
    nrerror("Too many iterations in cg()");
    return;
fine:
    free(f1dim_xt);
    free(xi);
    free(h);
    free(g);
    return;
} /* end cg() */


void linmin (double p[], double xi[], int n,
	     double *fret, double (*func)(double []))
{
    int j;
    double xx,xmin,fx,fb,fa,bx,ax;
    linmin_ncom=n;
    linmin_nrfunc=func;
    linmin_pcom = p;
    linmin_xicom = xi;
    ax=0.0;
    xx=1.0;
    mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
    *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
    for (j=0;j<n;j++)
    {
	xi[j] *= xmin;
	p[j] += xi[j];
    }
}  /* end linmin() */


double f1dim(double x)
{
    int j;
    for (j=0;j<linmin_ncom;j++)
	f1dim_xt[j]=linmin_pcom[j]+x*linmin_xicom[j];
    return ((*linmin_nrfunc)(f1dim_xt));
} /* end f1dim() */


void mnbrak (double *ax, double *bx, double *cx,
	     double *fa, double *fb, double *fc,
	     double (*func)(double))
{
    double ulim,u,r,q,fu,dum;
    *fa=(*func)(*ax);
    *fb=(*func)(*bx);
    if (*fb > *fa)
    {
	SHFT(dum,*ax,*bx,dum);
	SHFT(dum,*fb,*fa,dum);
    }
    *cx=(*bx)+GOLD*(*bx-*ax);
    *fc=(*func)(*cx);
    while (*fb > *fc)
    {
	r=(*bx-*ax)*(*fb-*fc);
	q=(*bx-*cx)*(*fb-*fa);
	u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
	    (2.0*SIGN(max(fabs(q-r),TINY),q-r));
	ulim=(*bx)+GLIMIT*(*cx-*bx);
	if ((*bx-u)*(u-*cx) > 0.0)
	{
	    fu=(*func)(u);
	    if (fu < *fc)
	    {
		*ax=(*bx);
		*bx=u;
		*fa=(*fb);
		*fb=fu;
		return;
	    }
	    else if (fu > *fb)
	    {
		*cx=u;
		*fc=fu;
		return;
	    }
	    u=(*cx)+GOLD*(*cx-*bx);
	    fu=(*func)(u);
	}
	else if ((*cx-u)*(u-ulim) > 0.0)
	{
	    fu=(*func)(u);
	    if (fu < *fc)
	    {
		SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx));
		SHFT(*fb,*fc,fu,(*func)(u));
	    }
	}
	else if ((u-ulim)*(ulim-*cx) >= 0.0)
	{
	    u=ulim;
	    fu=(*func)(u);
	}
	else
	{
	    u=(*cx)+GOLD*(*cx-*bx);
	    fu=(*func)(u);
	}
	SHFT(*ax,*bx,*cx,u);
	SHFT(*fa,*fb,*fc,fu);
    }
} /* end mnbrak() */


double brent
(double ax, double bx, double cx, 
 double (*f)(double), double tol, double *xmin)
{
    int iter;
    double a,b,d,etemp,fu,fv,fw,fx,p,q,
	r,tol1,tol2,u,v,w,x,xm;
    double e=0.0;
    
    a=(ax < cx ? ax : cx);
    b=(ax > cx ? ax : cx);
    x=w=v=bx;
    fw=fv=fx=(*f)(x);
    for (iter=1;iter<=ITMAX;iter++)
    {
	xm=0.5*(a+b);
	tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
	if (fabs(x-xm) <= (tol2-0.5*(b-a)))
	{
	    *xmin=x;
	    return fx;
	}
	if (fabs(e) > tol1)
	{
	    r=(x-w)*(fx-fv);
	    q=(x-v)*(fx-fw);
	    p=(x-v)*q-(x-w)*r;
	    q=2.0*(q-r);
	    if (q > 0.0) p = -p;
	    q=fabs(q);
	    etemp=e;
	    e=d;
	    if (fabs(p) >= fabs(0.5*q*etemp)
		|| p <= q*(a-x) || p >= q*(b-x))
		d=CGOLD*(e=(x >= xm ? a-x : b-x));
	    else
	    {
		d=p/q;
		u=x+d;
		if (u-a < tol2 || b-u < tol2)
		    d=SIGN(tol1,xm-x);
	    }
	}
	else
	{
	    d=CGOLD*(e=(x >= xm ? a-x : b-x));
	}
	u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
	fu=(*f)(u);
	if (fu <= fx)
	{
	    if (u >= x) a=x; else b=x;
	    SHFT(v,w,x,u);
	    SHFT(fv,fw,fx,fu);
	}
	else
	{
	    if (u < x) a=u; else b=u;
	    if (fu <= fw || w == x)
	    {
		v=w;
		w=u;
		fv=fw;
		fw=fu;
	    }
	    else if (fu <= fv || v == x || v == w)
	    {
		v=u;
		fv=fu;
	    }
	}
    }
    nrerror("Too many iterations in brent");
    *xmin=x;
    return fx;
}  /* end brent() */


void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
    fprintf(stderr,"Numerical Recipes run-time error...\n");
    fprintf(stderr,"%s\n",error_text);
    fprintf(stderr,"...now exiting to system...\n");
    exit(1);
}  /* end nrerror() */


void polint(double xa[], double ya[], int n, double x, double *y, double *dy)
{
	int i,m,ns=1;
	double den,dif,dift,ho,hp,w;
	double *c,*d;
	dif=fabs(x-xa[1]);
	c = (double *) malloc((n+1)*(sizeof(double)));
	d = (double *) malloc((n+1)*(sizeof(double)));
	for (i=1;i<=n;i++) {
		if ( (dift=fabs(x-xa[i])) < dif) {
			ns=i;
			dif=dift;
		}
		c[i]=ya[i];
		d[i]=ya[i];
	}
	*y=ya[ns--];
	for (m=1;m<n;m++) {
		for (i=1;i<=n-m;i++) {
			ho=xa[i]-x;
			hp=xa[i+m]-x;
			w=c[i+1]-d[i];
			if ((den=ho-hp) == 0.0)
			{
			    printf("Error in routine polint().\n");
			    exit(1);
			}
			den=w/den;
			d[i]=hp*den;
			c[i]=ho*den;
		}
		*y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
	}
	free(d);
	free(c);
}  /* end polint() */


/* "tee" screen output to a file */
void print (const char *format, ...)
{
    va_list ap, bp;
    char *p, *q, buffer[MAX_STRLEN];
    va_start (ap, format); va_start (bp, format);
    p = (char *) format;
    while ((q=strchr(p, '%'))!=NULL)
    {
	for (q++; (*q!='d') && (*q!='f') && (*q!='s'); q++);
	strcpy (buffer, p);
	buffer[q-p+1] = (char) 0;
	switch (*q)
	{
	case 'd':
	    printf ((const char *) buffer, va_arg(ap, int));
	    if (save_screen_output)
		fprintf (screen_output, (const char *) buffer,
			 va_arg(bp, int));
	    break;
	case 'f':
	    printf ((const char *) buffer, va_arg(ap, double));
	    if (save_screen_output)
		fprintf (screen_output, (const char *) buffer,
			 va_arg(bp, double));
	    break;
	case 's':
	    printf ((const char *) buffer, va_arg(ap, char *));
	    if (save_screen_output)
		fprintf (screen_output, buffer, va_arg(bp, char *));
	    break;
	default:
	    printf ("print(): \"%s\": no such format\n", buffer);
	    exit(1);
	}
	p = q+1;
    } 
    printf ("%s", p);
    if (save_screen_output) fprintf (screen_output, "%s", p);
    va_end (ap);
    va_end (bp);
    return;
} /* end print() */