/********************************************/
/*  Conjugate gradient relaxation of Argon  */
/*  solid at constant H or constant stress. */
/*                                          */
/*               Version 1.0 May.29, 1998   */
/*               Developed by Ju Li (MIT)   */
/********************************************/


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

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

int main()
{
    int i,j,k,l;
    double tmp, cc, dd, ee;
    char sentence[MAX_STRING_LENGTH]={0};
    double pote_before_relaxation, volume_before_relaxation; 
    double cm_shift[3]={0.,0.,0.};
    
    printf ("Reading in instructions...\n");
    /* Read instructions from input, such as "ron" file */
    fscanf (stdin,
	    "To relax  0. a provided defect  1. a user-defined configuration:\n%d\n\n",
	    &i);
    minimize_provided_defect = (i==0);
    minimize_user_defined_defect = (i==1);
    if ((!minimize_provided_defect) && (!minimize_user_defined_defect))
    {
	fprintf (stdout,"\n **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 = !strcmp(potential_function,"LJ6_12");
    WOON = !strcmp(potential_function,"WOON");
    if ((!LJ6_12) && (!WOON))
    {
	fprintf (stdout, "\n **Error: Choose between LJ6_12 or WOON pair potential\n");
	exit(1);
    }

    fscanf (stdin, "const_H or const_s mode:\n%s\n\n", sentence);
    const_H = !strcmp(sentence,"const_H");
    const_s = !strcmp(sentence,"const_s");
    if ((!const_H) && (!const_s))
    {
	fprintf (stdout,
		 "\n **Error: does not have \"%s\" mode, choose NEH or NTP.\n",
		 sentence);
	exit(1);
    }

    fscanf (stdin,
	    "Number of atoms in corresponding perfect crystal(fcc) state:\n%d\n\n",
	    &npfcc);
    
    fscanf (stdin, "Default supercell (11,22,33,12,13,23), in Angstrom:\n");
    fscanf (stdin, "%lf %lf %lf %lf %lf %lf\n\n",
	    &H0[0][0], &H0[1][1], &H0[2][2],
	    &H0[0][1], &H0[0][2], &H0[1][2]);
    /* convert to reduced units */
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    if (j>=i) H0[i][j] /= ULENGTH_IN_ANGSTROM;
	    else H0[j][i] = H0[i][j];

    fscanf (stdin,
	    "In const_s mode, should the above be enforced initially (yes/no)?\n%s\n\n",
	    sentence);
    H0_enforced = !strcasecmp (sentence,"yes");
    
    fscanf (stdin, "Default external stress (11,22,33,12,13,23), in MPa:\n");
    fscanf (stdin, "%lf %lf %lf %lf %lf %lf\n\n",
	    &sout[0][0], &sout[1][1], &sout[2][2],
	    &sout[0][1], &sout[0][2], &sout[1][2]);
    /* convert to reduced units */
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    if (j>=i) sout[i][j] /= USTRESS_IN_MPA;
	    else sout[j][i] = sout[i][j];
    
    fscanf (stdin, "Index of the provided defect that could be relaxed:\n");
    fscanf (stdin,
	    "(0. perfect_crystal  1. vacancy  2. interstitial)\n%d\n\n",
	    &provided_defect_index);
    i = 2;
    if ((provided_defect_index<0) || (provided_defect_index>i))
    {
	fprintf (stdout,"\n **Error: select provided defect index from 0 to %d.\n",i);
	exit(1);
    }
    
    fscanf (stdin,
	    "Amplitude of random position perturbations (in lattice constant):\n%lf\n\n",
	    &perturbation_amplitude);

    fscanf (stdin, "Random number generator seed (long integer):\n%ld\n\n", &iseed);
    srandom (iseed);
    
    fscanf (stdin, "Minimization tolerance:\n%lf\n\n", &minimization_tolerance);

    fscanf (stdin,
	    "Save to .xyz files the position of atoms around iatom (-2=don't, -1=use_default)\n");
    fscanf (stdin,
	    "within R[A] before (default = *.bef) and after minimization (default = *.aft):\n");
    fscanf (stdin, "%d %lf %s %s\n\n",
	    &iatom, &shadow_radius, fn_xyz_bef, fn_xyz_aft);
    apply_shadow_positioning = (iatom!=-2); 
    use_default_iatom = (iatom==-1);
    shadow_radius /= ULENGTH_IN_ANGSTROM;
    if (!strcasecmp (fn_xyz_bef,"default")) sprintf(fn_xyz_bef, "%s.bef", defect_name);
    if (!strcasecmp (fn_xyz_aft,"default")) sprintf(fn_xyz_aft, "%s.aft", 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\n",
	  potential_function, CUTOFF, CUTOFF*ULENGTH_IN_ANGSTROM);
    print("Angstrom under \"%s\" mode.\n\n", const_H?"const_H":"const_s");
    
    if (const_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 */
    print("The defect will be benchmarked against perfect crystal\n");
    print("(BENCHMARK_NP=%d) at P=0:\n", BENCHMARK_NP);
    assign_fcc_lattice (BENCHMARK_NC, s_benchmark);
    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 ((int)BENCHMARK_NP, 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 ((int)BENCHMARK_NP, H_benchmark, s_benchmark)/BENCHMARK_NP;
    volumep_benchmark =
	H_benchmark[0][0]*H_benchmark[0][0]*H_benchmark[0][0]/BENCHMARK_NP;
    print("benchmarked lattice constant = %f = %f A,\n", tmp, tmp*ULENGTH_IN_ANGSTROM);
    print("benchmarked volume per particle = %f = %f A^3,\n", tmp*tmp*tmp/4,
	  tmp*tmp*tmp/4*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
    print("benchmarked energy per particle = %f = %f eV,\n", ep_benchmark,
	  ep_benchmark*UENERGY_IN_JOULE/EV_IN_JOULE);
    print("benchmark finished.\n\n");

    /* deduce the maximum number of particles */
    if (minimize_provided_defect)
    {
	nc = rint(pow(npfcc/4.,1./3.)); 
	if (4*nc*nc*nc!=npfcc)
	{
	    print("assign_fcc_lattice(): original npfcc=%d, truncated to %d.\n",
		  npfcc, 4*nc*nc*nc);
	    npfcc = 4*nc*nc*nc;
	}
	np = (provided_defect_excess[provided_defect_index]>0)?
	    npfcc+provided_defect_excess[provided_defect_index]:npfcc;
    }
    else if (minimize_user_defined_defect)
    {
	if ((config_read=fopen(fn_config_read, "r"))==NULL)
	{
	    fprintf(stdout, "\n **Error: cannot locate file \"%s\"\n.", fn_config_read);
	    exit(1);
	}
	fscanf(config_read, "Number of particles = %d\n\n", &np);
	close(config_read);
    }
    else exit(1);
    print("Allocating memories for maximum of %d particles...\n",np);
    
    /* allocate memories for the main program */
    s = (double *) malloc((3*np+6)*sizeof(double));
    f = (double *) malloc((3*np+6)*sizeof(double));
    /* keep mass[] to preserve isotope defects, but don't do anything with it */
    mass = (double *) malloc(np*sizeof(double));
    print("Memory allocation complete.\n\n");

    /* assign reduced coordinates */
    if (minimize_provided_defect)
    {
	print("The initial configuration will be constructed from scratch...\n");
	/* assign fcc structure according to nc=i */
	assign_fcc_lattice (nc,s);
	print("fcc lattice construction complete, npfcc = %d.\n\n", npfcc);
	np = npfcc; 
    }
    else if (minimize_user_defined_defect)
    {
	read_from_config(fn_config_read);
	print("Loading initial configuration from \"%s\" complete,\n", fn_config_read);
	print("Actual number of particles = %d.\n\n", np);
    }
    else exit(1);

    /* H matrix and accessories */
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	{
	    if (minimize_provided_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_ANGSTROM,
	  H[0][1]*ULENGTH_IN_ANGSTROM,
	  H[0][2]*ULENGTH_IN_ANGSTROM);	  
    print("Initial H[A] = | %.5f %.5f %.5f |\n",
	  H[1][0]*ULENGTH_IN_ANGSTROM,
	  H[1][1]*ULENGTH_IN_ANGSTROM,
	  H[1][2]*ULENGTH_IN_ANGSTROM);
    print("               | %.5f %.5f %.5f |\n\n",
	  H[2][0]*ULENGTH_IN_ANGSTROM,
	  H[2][1]*ULENGTH_IN_ANGSTROM,
	  H[2][2]*ULENGTH_IN_ANGSTROM);

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

    /* constructing defect configuration */
    if (minimize_provided_defect)
    {
	switch(provided_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;
	    s[3*np+1] = s[1]+0.25/nc;
	    s[3*np+2] = s[2]+0.25/nc;
	    np++;
	    break;
	   
	default:
	    print("\n **Error: provided_defect[%d] does not exist.\n",provided_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_ANGSTROM, iatom);
	print("be recorded in \"%s\" before minimization and\n", fn_xyz_bef);
	print("\"%s\" after minimization, in .xyz format\n", fn_xyz_aft);
	print("that can be rendered by Rasmol (free software on web).\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_JOULE/EV_IN_JOULE,
	  pote_before_relaxation/np*UENERGY_IN_JOULE/EV_IN_JOULE);
    print("so the defect formation energy (unrelaxed)\n");
    print("= %f - %d x %f = %f eV.\n",
	  pote_before_relaxation*UENERGY_IN_JOULE/EV_IN_JOULE, np,
	  ep_benchmark*UENERGY_IN_JOULE/EV_IN_JOULE,
	  (pote_before_relaxation-ep_benchmark*np)*UENERGY_IN_JOULE/EV_IN_JOULE);
    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_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM,
	  volume/np*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
    print("so the excess volume (unrelaxed)\n");
    print("= %f - %d x %f = %f A^3.\n",
	  volume*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM,
	  np, volumep_benchmark*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM,
	  (volume-np*volumep_benchmark)*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
    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_xyz_bef);
    print("*************************************************************\n\n");
    
    print("Random number seed = %ld, 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_ANGSTROM, cc*ULENGTH_IN_ANGSTROM);
    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 (const_H)
	cg (s, 3*np, minimization_tolerance, &i, &cc, value, force);
    else if (const_s)
    {
	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);
    }
    else
    {
	print("Minimization mode does not exist.\n");
	exit(1);
    }
    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_JOULE/EV_IN_JOULE, cc*UENERGY_IN_JOULE/EV_IN_JOULE/np);
    print("so the defect formation energy (relaxed)\n");
    print("= %f - %d x %f = %f eV.\n",
	  cc*UENERGY_IN_JOULE/EV_IN_JOULE, np,
	  ep_benchmark*UENERGY_IN_JOULE/EV_IN_JOULE,
	  (cc-ep_benchmark*np)*UENERGY_IN_JOULE/EV_IN_JOULE);
    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_JOULE/EV_IN_JOULE);

    print("             | %.5f %.5f %.5f |\n",
	  H[0][0]*ULENGTH_IN_ANGSTROM,
	  H[0][1]*ULENGTH_IN_ANGSTROM,
	  H[0][2]*ULENGTH_IN_ANGSTROM);	  
    print("Final H[A] = | %.5f %.5f %.5f |\n",
	  H[1][0]*ULENGTH_IN_ANGSTROM,
	  H[1][1]*ULENGTH_IN_ANGSTROM,
	  H[1][2]*ULENGTH_IN_ANGSTROM);
    print("             | %.5f %.5f %.5f |\n\n",
	  H[2][0]*ULENGTH_IN_ANGSTROM,
	  H[2][1]*ULENGTH_IN_ANGSTROM,
	  H[2][2]*ULENGTH_IN_ANGSTROM);
    
    print("The total volume = %.3f A^3, volume/atom = %f A^3,\n", 
	  volume*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM,
	  volume/np*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
    print("so the excess volume (relaxed)\n");
    print("= %f - %d x %f = %f A^3.\n",
	  volume*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM,
	  np, volumep_benchmark*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM,
	  (volume-np*volumep_benchmark)*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
    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)
	  *ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM*ULENGTH_IN_ANGSTROM);
    
    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_xyz_aft);
    print("************************************************************\n\n");

    /* write on configuration file */
    if (write_config)
    {
	print("Saving relaxed configuration on file \"%s\".\n", fn_config_write);
	write_on_config (fn_config_write);
    }
    
    /* free all allocated memories */
    free (f);
    free (s);
    free (mass);
    print("Memories freed.\n");
    
    if (save_screen_output)
    {
	print("Saving screen output 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 (const_H)
	free_energy = pair_potential_energy (np, H, s);
    else if (const_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_JOULE/EV_IN_JOULE);
    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 (const_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;
}


/* 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,k;
    double tmp[3][3], tmpp[3][3];
    
    for (i=0; i<3; i++)
	for (k=0; k<3; k++)
	{
	    M[i][k] = eta[i][k];
	    if (i==k) M[i][k]++;
	}

    /* 2nd order */
    matmul (eta,eta,tmp);
    mateqv (tmp,tmpp);
    matadd_mul (M,-0.5,tmp,M);

    /* 3rd order */
    matmul (tmpp,eta,tmp);
    mateqv (tmp,tmpp);
    matadd_mul (M,+0.5,tmp,M);
    
    /* 4th order */
    matmul (tmpp,eta,tmp);
    mateqv (tmp,tmpp);
    matadd_mul (M,-5./8.,tmp,M);

    /* 5th order */
    matmul (tmpp,eta,tmp);
    mateqv (tmp,tmpp);
    matadd_mul (M,+7./8.,tmp,M);

    return;
} /* end eta_to_M() */


/* see 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, k;
    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 (k=0; k<3; k++)
	    etac[i][k] = c*eta[i][k];
    eta_to_M (etac, Mc);
    matmul (H0, Mc, tmp);
    volumec = matinv (tmp, tmpp);
    matinv (Mc, MIc);
    matmul (MIc, sout, tmpp);
    matmul (tmpp, MIc, tmp);
    cc = 0.;
    for (i=0; i<3; i++)
	for (k=0; k<3; k++)
	    cc += tmp[i][k]*eta[i][k]*volumec;
    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.));
    /* integration by trapzoidal rule with progressively decreasing h */
    for (j=2, i=1; j<=JMAX; j++, i<<=1)
    {
	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[3] and save to a file
   in .xyz format */
#define SHIFT_TINY 0.0002
void shadow_positioning (char *fn_xyz)
{
    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 *xyz; /* .xyz files */

    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("** error: shadow_positioning(): out of bound, change MAX_SHADOW_CAPTURE\n");
		exit(1);
	    }
	    dx[ip] = rx*ULENGTH_IN_ANGSTROM;
	    dy[ip] = ry*ULENGTH_IN_ANGSTROM;
	    dz[ip] = rz*ULENGTH_IN_ANGSTROM;
	    d[ip] = sqrt(r2)*ULENGTH_IN_ANGSTROM;
	    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");
    xyz = fopen(fn_xyz, "w+");
    fprintf(xyz, "   %d\n", ip); 
    fprintf(xyz, "   Defect %s\n", defect_name); 
    for (i=0; i<ip; i++)
    {
	k = index[i]; 
	fprintf(xyz, "   Ar %f %f %f\n", 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 .xyz format)\n", fn_xyz);
    fclose(xyz);
    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 */

double pair_potential_energy (int np, double H[3][3], double *s)
{
    int i, j;
    double residue, sxij, syij, szij, rx, ry, rz, r, r2, r6, r12, pote, a, b;
    
    if (LJ6_12)
	residue = 4*(pow(CUTOFF,-12.)-pow(CUTOFF,-6.)); 
    else if (WOON)
	residue = WOON_EPSILON*
	    (WOON_Q*exp(-WOON_K1*(CUTOFF-WOON_RE))-
	     (WOON_Q+1)*exp(-WOON_K2*(CUTOFF-WOON_RE))); 
    
    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<CUTOFF*CUTOFF)
	    {
		if (LJ6_12)
		{
		    r2  = 1./r2;
		    r6  = r2*r2*r2;
		    r12 = r6*r6;
		    pote += 4*(r12-r6)-residue;
		}
		else if (WOON)
		{
		    r = sqrt(r2);
		    a = WOON_EPSILON * WOON_Q * exp(-WOON_K1 * (r-WOON_RE));
		    b = WOON_EPSILON * (WOON_Q+1) * exp(-WOON_K2 * (r-WOON_RE)); 
		    pote += a-b-residue;
		}
	    }
	} 
    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;

    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<CUTOFF*CUTOFF)
	    {
		if (LJ6_12)
		{
		    r2  = 1./r2;
		    r6  = r2*r2*r2;
		    r12 = r6*r6;
		    ff  = (48.*r12-24.*r6)*r2;
		}
		else if (WOON)
		{
		    r = sqrt(r2);
		    a = WOON_EPSILON * WOON_Q * exp(-WOON_K1 * (r-WOON_RE));
		    b = WOON_EPSILON * (WOON_Q+1) * exp(-WOON_K2 * (r-WOON_RE)); 
		    ff = (WOON_K1*a-WOON_K2*b)/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() */


/* save system configuration in readable ASCII format */
void write_on_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 Angstrom\n",  H[0][0]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(1,2) = %f Angstrom\n",  H[0][1]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(1,3) = %f Angstrom\n\n", H[0][2]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(2,1) = %f Angstrom\n",  H[1][0]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(2,2) = %f Angstrom\n",  H[1][1]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(2,3) = %f Angstrom\n\n", H[1][2]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(3,1) = %f Angstrom\n",  H[2][0]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(3,2) = %f Angstrom\n",  H[2][1]*ULENGTH_IN_ANGSTROM);
    fprintf (config, "H(3,3) = %f Angstrom\n\n", H[2][2]*ULENGTH_IN_ANGSTROM);
    fprintf (config,
	     "Mass(amu)     sx        sy        sz     sx1(1/ps)  sy1(1/ps)  sz1(1/ps)\n");
    for (i=0;i<3*np;i+=3)
	fprintf (config, "%f  %f  %f  %f  %.6f  %.6f  %.6f\n",
		 mass[i/3], s[i], s[i+1], s[i+2], 0., 0., 0.);
    fclose (config);
    return;
}  /* end write_config() */


/* load system coordinate configuration from file */
void read_from_config(char filename[])
{
    int i,j;
    double dummy; 
    FILE *config; 
    if ((config=fopen(filename, "r"))==NULL)
    {
      fprintf(stdout, "\n **Error: cannot locate file \"%s\"\n.", filename);
      exit(1);
    }
    fscanf (config, "Number of particles = %d\n\n", &np);
    fscanf (config, "H(1,1) = %lf Angstrom\n",  &H[0][0]);
    fscanf (config, "H(1,2) = %lf Angstrom\n",  &H[0][1]);
    fscanf (config, "H(1,3) = %lf Angstrom\n\n", &H[0][2]);
    fscanf (config, "H(2,1) = %lf Angstrom\n",  &H[1][0]);
    fscanf (config, "H(2,2) = %lf Angstrom\n",  &H[1][1]);
    fscanf (config, "H(2,3) = %lf Angstrom\n\n", &H[1][2]);
    fscanf (config, "H(3,1) = %lf Angstrom\n",  &H[2][0]);
    fscanf (config, "H(3,2) = %lf Angstrom\n",  &H[2][1]);
    fscanf (config, "H(3,3) = %lf Angstrom\n\n", &H[2][2]);
    for (i=0;i<3;i++)
	for (j=0;j<3;j++)
	    H[i][j] /= ULENGTH_IN_ANGSTROM; 
    volume = matinv(H, HI);
    fscanf(config,
	   "Mass(amu)     sx        sy        sz     sx1(1/ps)  sy1(1/ps)  sz1(1/ps)\n");
    for (i=0;i<3*np;i+=3)
	fscanf (config, "%lf  %lf  %lf  %lf  %lf  %lf  %lf\n",
		&mass[i/3], &s[i], &s[i+1], &s[i+2], &dummy,  &dummy, &dummy);
    /* keep mass[] to preserve isotope defects, but don't do anything with it */
    fclose (config);
    return;
} /* end read_config() */


/* Assign fcc reduced coordinates according to NC; s[] must be long enough */
void assign_fcc_lattice(int nc, double *s)
{
    static int i,j,k,l,ip;
    static 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;
    ip = 0;
    for (i=0; i<nc; i++)
	for (j=0; j<nc; j++)
	    for (k=0; k<nc; k++, ip+=12)
		for (l=0; l<12; l+=3)
		{
		    s[l+ip]   = (x[l]+i+1./8.)/nc;
		    s[l+1+ip] = (x[l+1]+j+1./8.)/nc;
		    s[l+2+ip] = (x[l+2]+k+1./8.)/nc;
		}
    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(FMAX(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_STRING_LENGTH];
    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(" **error: 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() */