/***************************************************/
/* Conjugate gradient relaxation of Zr/C atomistic */
/* configuration at constant H or stress.          */
/*                                                 */
/*           Ver 2.0,  Apr 25 2000, Ju Li (MIT).   */
/***************************************************/

#include "relax.h"

#include "potential.c"
#include "elastic_constants.c"
#include "reference.c"
#include "work.c"

Aapp_Define_Config;
Chemtab ct[1]={0};
Tp *tp=NULL;
Neighborlist N[1]={0};
M3 sout, stress, H0;
double pote, volume, workdone, *f=NULL;

bool minimize_builtin_defect, minimize_user_defined_defect;
bool constant_H, constant_S;
int builtin_defect_index;
double perturbation_amplitude;
unsigned int random_seed;
double minimization_tolerance;
int report_frequency;
bool save_screen_output;
bool write_config;

TermString defect_name;
TermString fname_read;
TermString fname_out;
TermString fname_write;

const char *builtin_name [BUILTIN_MAX] =
{
    "perfect_crystal",
    "WC",
    "C_vacancy",
    "Zr_vacancy",
    "surface_001",
    "ledge_001_a",
    "ledge_001_b",
    "ledge_001_c",
    "ledge_001_d",
    "ZrC_995",
    "ZrC_99",
    "ZrC_98",
    "ZrC_97",
    "ZrC_96",
    "ZrC_95",
    "ZrC_90",
    "ZrC_80",
    "ZrC_70",
    "ZrC_60",
    "ZrC_50",
    "ZrC_40",
    "ZrC_30",
    "ZrC_20",
    "ZrC_10",
    "ZrC_0",
};

static double pote_before_relaxation, volume_before_relaxation;
static double atom_force_norm;

void report (FILE *ft)
{
    M3 eta;
    Fprintf (ft, "** Status Report at %s Potential Evaluation **\n\n",
             word_for_order(CG_state.number_of_evaluations+1));
    Neighborlist_print_statistics (np, N, ft); Fcr(ft);
    Fprintf (ft, "********************* "
             "After Minimization *********************\n");
    Fprintf (ft, "The total energy = %f eV, energy/atom = %f eV,\n",
             pote*UENERGY_IN_EV, pote*UENERGY_IN_EV/np);
    Fprintf (ft, "so the defect formation energy (relaxed)\n");
    Fprintf (ft, "= %f - %d x %f = %f eV.\n", pote*UENERGY_IN_EV, np,
             ep_reference*UENERGY_IN_EV, (pote-ep_reference*np)*UENERGY_IN_EV);
    Fprintf (ft, "(energy reference is perfect crystal at zero stress)\n");
    Fprintf (ft, "and it has increased by %f eV during the relaxation.\n",
             (pote-pote_before_relaxation)*UENERGY_IN_EV);
    if (constant_S)
    {
        M3MULTIPLY (ULENGTH_IN_A, H, eta);
        S3fPR (ft, "Current H = %M A.\n", eta);
        Fprintf (ft, "The total volume = %.3f A^3, volume/atom = %f A^3,\n",
                 volume*UVOLUME_IN_A3, volume/np*UVOLUME_IN_A3);
        Fprintf (ft, "so the excess volume (relaxed)\n");
        Fprintf (ft, "= %f - %d x %f = %f A^3.\n", volume*UVOLUME_IN_A3,
                 np, volumep_reference*UVOLUME_IN_A3,
                 (volume-np*volumep_reference)*UVOLUME_IN_A3);
        Fprintf (ft, "(reference is perfect crystal at zero stress)\n");
        Fprintf (ft, "and it has increased by %f A^3 during the run,\n",
                 (volume-volume_before_relaxation)*UVOLUME_IN_A3);
        Lagrangian_strain(H0, H, eta);
        S3fPR (ft, "with Lagrangian strain = %M.\n", eta);
        Fprintf(ft, "work done = %f eV, G = %f eV.\n",
                workdone*UENERGY_IN_EV, (pote-workdone)*UENERGY_IN_EV);
    }
    Fprintf (ft, "atom residual |f| = %e eV/A,\n",
             atom_force_norm*UENERGY_IN_EV/ULENGTH_IN_A);
    M3MULTIPLY (USTRESS_IN_MPA, stress, eta);
    S3fPR (ft, "Current stress = %M MPa.\n", eta);
    Fprintf (ft, "**********************"
             "****************************************\n\n");
    if (write_config)
    {
        Fprintf (ft, "Saving relaxed configuration on file \"%s\".\n",
                 fname_write);
        Config_save (Config_Aapp_to_Alib, TRUE, fname_write);
    }
    return;
} /* end report() */


/* objective function of CG minimization */
static double free_energy (double *s_new, double *f)
{
    register int i;
    double equalizer;
    V3 fs;
    M3 U, M, HI, MI, SUM, FSUM;
    s = s_new;
    if (constant_S)
    {
        equalizer = pow((double)np, 2./3.);
        U[0][0] = s[DIMENSION*np] / equalizer;
        U[1][1] = s[DIMENSION*np+1] / equalizer;
        U[2][2] = s[DIMENSION*np+2] / equalizer;
        U[1][0] = U[0][1] = s[DIMENSION*np+3] / equalizer;
        U[2][0] = U[0][2] = s[DIMENSION*np+4] / equalizer;
        U[2][1] = U[1][2] = s[DIMENSION*np+5] / equalizer;
        M3ADDDIAG (U, 1, M);
        M3MUL (H0, M, H);
    }
    M3InV (H, HI, volume);
    Neighborlist_Maintenance (Config_Aapp_to_Alib, stdout, ct, &tp, N);
    pote = potential (np,H,s,tp,N,f,stress);
    VLENGTH(DIMENSION*np, f, i, atom_force_norm);
    for (i=0; i<np; i++)
    {
        M3mV3 (H, &f[DIMENSION*i], fs);
        V3NEG (fs, &f[DIMENSION*i]);
    }
    if (constant_S)
    {
        M3ADD (sout, stress, SUM);
        M3INV (M, MI, workdone);
        M3MUL (SUM, MI, FSUM);
        M3Symmetrize (FSUM);
        f[DIMENSION*np]   = -FSUM[0][0] * volume / equalizer;
        f[DIMENSION*np+1] = -FSUM[1][1] * volume / equalizer;
        f[DIMENSION*np+2] = -FSUM[2][2] * volume / equalizer;
        f[DIMENSION*np+3] = -FSUM[0][1] * volume * 2. / equalizer;
        f[DIMENSION*np+4] = -FSUM[0][2] * volume * 2. / equalizer;
        f[DIMENSION*np+5] = -FSUM[1][2] * volume * 2. / equalizer;
        workdone = work (H0, sout, U);
    }
    fprintf (stdout, "k = %d, |f| = %e eV/A, pote = %f eV,\n",
             CG_state.number_of_evaluations,
             atom_force_norm*UENERGY_IN_EV/ULENGTH_IN_A, pote*UENERGY_IN_EV);
    if (constant_S)
        fprintf(stdout, "workdone=%f eV, G=%f eV, "
                "stress discrepancy=%e MPa.\n",
                workdone*UENERGY_IN_EV, (pote-workdone)*UENERGY_IN_EV,
                M32NORM(SUM)*USTRESS_IN_MPA);
    if ( (report_frequency > 0) &&
         (CG_state.number_of_evaluations % report_frequency ==
          report_frequency-1) )
    {
        fcr(stdout);
        report(ft);
        Fcr(ft);
    }
    if (constant_S) return (pote-workdone);
    else return (pote);
} /* end free_energy() */


/* Create random carbon vacancies to make ZrC_X */
static double Random_C_vacancies (double X, Alib_Declare_Config)
{
    int i,j,k,Zr_count,C_count,new_C_count,kill_count,*C_idx;
    MALLOC( Random_C_vacancies, C_idx, *np, int );
    for (C_count=i=0; i<*np; i++)
        if (!strcmp(SYMBOL(i)," C")) C_idx[C_count++] = i;
    Zr_count = (*np) - C_count;
    new_C_count = floor(Zr_count * X);
    kill_count = C_count - new_C_count;
    if (kill_count < 0)
        pe ("Random_C_vacancies: kill_count = %d\n", kill_count);
    for (i=0; i<kill_count; i++)
    {
        j = fran (i, C_count);
        SWAP (C_idx[i], C_idx[j], k);
    }
    Config_rid_of_blacklist (kill_count, C_idx, Config_Alib_to_Alib);
    free (C_idx);
    return (DOUBLE(new_C_count)/Zr_count);
} /* end Random_C_vacancies() */


void save_config (char *fn)
{
    int i;
    FILE *fp;
    if (!(Fvalid(fn) && strncasecmp(fn,"NO",2) && strncasecmp(fn,"NULL",4)))
        return;
    fp = wopen(fn);
    fprintf (fp, "Number of particles = %d\n", np);
    fprintf (fp, "H0(1,1) = %.15g A\n", H[0][0]);
    fprintf (fp, "H0(1,2) = %.15g A\n", H[0][1]);
    fprintf (fp, "H0(1,3) = %.15g A\n", H[0][2]);
    fprintf (fp, "H0(2,1) = %.15g A\n", H[1][0]);
    fprintf (fp, "H0(2,2) = %.15g A\n", H[1][1]);
    fprintf (fp, "H0(2,3) = %.15g A\n", H[1][2]);
    fprintf (fp, "H0(3,1) = %.15g A\n", H[2][0]);
    fprintf (fp, "H0(3,2) = %.15g A\n", H[2][1]);
    fprintf (fp, "H0(3,3) = %.15g A\n", H[2][2]);
    fprintf (fp, ".NO_VELOCITY.\n");
    fprintf (fp, "entry_count = %d\n", 3+3);
    fprintf (fp, "auxiliary[0] = Fx [eV/A]\n");
    fprintf (fp, "auxiliary[1] = Fy [eV/A]\n");
    fprintf (fp, "auxiliary[2] = Fz [eV/A]\n");
    for (i=0; i<np; i++)
    {
        if ( (i==0) || ( (i>0) && (tp[i]!=tp[i-1]) ) )
        {
            fprintf (fp, "%f\n", mass[i]*UMASS_IN_AMU);
            fprintf (fp, "%s\n", SYM(i));
        }
        fprintf (fp, "%.5g %.5g %.5g %.5g %.5g %.5g\n",
                 s[3*i], s[3*i+1], s[3*i+2],
                 f[3*i], f[3*i+1], f[3*i+2]);
    }
    fclose(fp);
    Fprintf (ft, "Configuration saved on file \"%s\".\n\n", fn);
    return;
} /* end save_config() */


int main (int argc, char *argv[])
{
    int i, nc[3], swapped=FALSE;
    double cc, dd;
    char *taglist=NULL;
    TermString buf;
    FILE *fp;
    V3 ss, xx, cm;
    M3 HI;
    
    fscanf (stdin, "Potential function to be used (NOC):\n%s\n\n",
            potential_name);
    NOC = !strcasecmp(potential_name, "NOC");
    if (!NOC) pe ("Please choose NOC potential\n");

    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))
        pe ("(To relax) select options 0 or 1.\n");
    
    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", fname_read);
    if (!strcasecmp(fname_read, "default"))
        sprintf(fname_read, "%s", defect_name);

    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))
        pe ("do not have \"%s\" mode,\n"
            "choose between constant_H and constant_S mode.\n", buf);

    fscanf (stdin, "Default multiplications of the unit cell:\n");
    fscanf (stdin, "%d %d %d\n\n", &nc[0], &nc[1], &nc[2]);

    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];
    M3DividE (sout, USTRESS_IN_MPA);  /* convert to reduced units */

    fscanf(stdin, "Index of the built-in defect that could be relaxed:\n");
    for (i=0; i<6; i++) fscanf_skipline(stdin);
    fscanf(stdin, "%d\n\n", &builtin_defect_index);
    if (minimize_builtin_defect)
    {
        if (OUW(builtin_defect_index, BUILTIN_MAX))
            pe ("select built-in defect index from 0 to %d.\n", BUILTIN_MAX-1);
        sprintf(defect_name, "%s", builtin_name[builtin_defect_index]);
    }

    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", &random_seed);
    if (random_seed == 0) random_seed = time((time_t *)NULL);
    Randomize(random_seed);

    fscanf (stdin, "Minimization tolerance (default=EPS):\n%s\n\n", buf);
    if (!strcasecmp(buf,"default")) minimization_tolerance = EPS;
    else sscanf(buf, "%lf", &minimization_tolerance);

    fscanf (stdin, "Maximum number of evaluations & report frequency "
            "(default,inf):\n%s", buf);
    if (!strcasecmp(buf,"default"))
        CG_option.max_potential_evaluations=DEFAULT_MAX_POTENTIAL_EVALUATIONS;
    else if (!strcasecmp(buf,"inf"))
        CG_option.max_potential_evaluations = HUGE_INT;
    else sscanf(buf, "%d", &CG_option.max_potential_evaluations);
    fscanf (stdin, "%s\n\n", buf);
    if (!strcasecmp(buf,"default")) report_frequency=DEFAULT_REPORT_FREQUENCY;
    else if (!strcasecmp(buf,"INF")) report_frequency = HUGE_INT;
    else sscanf(buf, "%d", &report_frequency);

    fscanf (stdin, "Save screen output to file (default=*.out, "
            "NULL=don't):\n%s\n\n", fname_out);
    if ( save_screen_output=strcasecmp(fname_out, "NULL") )
    {
        if ( !strcasecmp(fname_out, "default") )
            sprintf (fname_out, "%s.out", defect_name);
        fp = wopen (fname_out);
        ftie (stdout, fp);
    }
    else ft = stdout;

    fscanf (stdin, "Save optimized configuration on file "
            "(default=*, NULL=don't):\n%s\n", fname_write);
    if ( write_config=strcasecmp(fname_write, "NULL") )
        if ( !strcasecmp(fname_write, "default") )
            sprintf (fname_write, "%s", defect_name);

    rebuild_reference (ft); Fcr(ft);

    Fprintf (ft, "Free energy will be minimized for defect \"%s\"\n",
             defect_name);
    Fprintf (ft, "using %s potential with RCUT_ZrZr=%g, RCUT_ZrC=%g\n",
             potential_name, RCUT_ZrZr, RCUT_ZrC);
    Fprintf (ft, "under \"%s\" condition.\n\n",
             constant_H ? "constant_H" : "constant_S");

    if (constant_S)
    {
        M3MULTIPLY (USTRESS_IN_MPA, sout, stress);
        S3fPR (ft, "sout = %M MPa.\n ", stress);
    }

    if (minimize_builtin_defect)
    {
        Fprintf (ft, "To construct \"%s\",\n"
                 "we first build %d x %d x %d crystal.\n\n",
                 builtin_name[builtin_defect_index], nc[0],nc[1],nc[2]);
        switch(builtin_defect_index)
        {
            case BUILTIN_XTAL: /* perfect crystal */
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                break;
            case BUILTIN_WC:  /* perfect crystal */
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_WC, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A / SQRT2, -1.);
                break;
            case BUILTIN_C_VACANCY: /* carbon vacancy */
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                for (i=0; i<np; i++) if (!strcmp(SYM(i)," C")) break;
                Fprintf (ft, "**** Then atom %d (1-%d), a %s of mass = %g"
                         " amu,\n**** at s = (%g %g %g), is REMOVED.\n\n",
                         i+1, np, SYM(i), mass[i]*UMASS_IN_AMU,
                         s[DIMENSION*i], s[DIMENSION*i+1], s[DIMENSION*i+2]);
                Config_DECIMATE (i,Config_Aapp_to_Alib);
                break;
            case BUILTIN_Zr_VACANCY:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                for (i=0; i<np; i++) if (!strcmp(SYM(i),"Zr")) break;
                Fprintf (ft, "**** Then atom %d (1-%d), a %s of mass = %g"
                         " amu,\n**** at s = (%g %g %g), is REMOVED.\n\n",
                         i+1, np, SYM(i), mass[i]*UMASS_IN_AMU,
                         s[DIMENSION*i], s[DIMENSION*i+1], s[DIMENSION*i+2]);
                Config_DECIMATE (i,Config_Aapp_to_Alib);
                break;
            case BUILTIN_LEDGE_001_C:
            case BUILTIN_LEDGE_001_D:
                if (builtin_defect_index == BUILTIN_LEDGE_001_C)
                    builtin_defect_index = BUILTIN_LEDGE_001_A;
                if (builtin_defect_index == BUILTIN_LEDGE_001_D)
                    builtin_defect_index = BUILTIN_LEDGE_001_B;
                swapped = TRUE;
            case BUILTIN_LEDGE_001_A:
            case BUILTIN_LEDGE_001_B:
            case BUILTIN_SURFACE_001:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                if ( (builtin_defect_index == BUILTIN_LEDGE_001_A) ||
                     (builtin_defect_index == BUILTIN_LEDGE_001_B) )
                {
                    dd = DOUBLE(nc[0]/2) / nc[0];
                    if (builtin_defect_index == BUILTIN_LEDGE_001_A)
                        cc = (nc[1]/2 + 0.5) / nc[1];
                    else cc = (nc[1]/2 + 1.) / nc[1];
                    CALLOC ( main, taglist, np, char );
                    for (i=0; i<np; i++)
                        if ( (s[3*i+2] < 0.5 / nc[2]) &&
                             ( (s[3*i] > dd ) ||
                               ( (s[3*i] < 0.5 / nc[0] ) &&
                                 (s[3*i+1] > cc ) ) ) ) taglist[i] = 1;
                    if (builtin_defect_index == BUILTIN_LEDGE_001_B)
                    {
                        cc = (nc[2]/2+.5) / nc[2];
                        dd = (nc[2]/2+1.) / nc[2];
                        for (i=0; i<np; i++)
                            if ( (s[3*i+2] > cc) && (s[3*i+2] < dd) &&
                                 (!strcmp(SYM(i)," C")) )
                            {
                                taglist[i] = 1;
                                break;
                            }
                    }
                    Config_rid_of_taglist (taglist, Config_Aapp_to_Alib);
                    Free(taglist);
                    if (swapped)
                    {
                        rebind_ct (Config_Aapp_to_Alib, "Zr C", ct, &tp, NULL);
                        for (i=0; i<np; i++) tp[i] = Zr_ATOM + C_ATOM - tp[i];
                        rematch_ct(ct, &tp, Config_Aapp_to_Alib);
                        swapped = FALSE;
                    }
                }
                for (i=0; i<np; i++)
                    s[3*i+2] = .5+(s[3*i+2]-.5)*H[2][2]/(H[2][2]+20.);
                H[2][2] += 20.;
                break;
            case BUILTIN_ZrC_995:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.995, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_99:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.99, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_98:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.98, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_97:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.97, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_96:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.96, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_95:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.95, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_90:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.90, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_80:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.80, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_70:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.70, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_60:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.60, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_50:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.50, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_40:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.40, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_30:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.30, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_20:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.20, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_10:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.10, Config_Aapp_to_Alib) );
                break;
            case BUILTIN_ZrC_0:
                Config_rebuild_Xtal
                    (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
                     Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1.,
                     lattice_constant_reference * ULENGTH_IN_A);
                Fprintf ( ft, "**** ZrC_X = %g\n\n",
                          Random_C_vacancies (0.0, Config_Aapp_to_Alib) );
                break;
            default:
                pe ("defect \"%s\" is not implemented yet.\n",
                    builtin_name[builtin_defect_index]);
        }                    
    }
    else Config_Load (fname_read, ft, Config_Aapp_to_Alib);
    rebind_ct (Config_Aapp_to_Alib, "Zr C", ct, &tp, ft); Fcr(ft);
    Config_fold_into_PBC (Config_Aapp_to_Alib);
    REALLOC (main, f, DIMENSION*np+DIMENSION*(DIMENSION+1)/2, double);
    REALLOC (main, s, DIMENSION*np+DIMENSION*(DIMENSION+1)/2, double);
    Neighborlist_Recreate_Form (Config_Aapp_to_Alib, ct, N);
    if (constant_S) N->min_shrinkage = DEF_MIN_SHRINKAGE;
    N->rcut[0] = RCUT_ZrZr;
    N->rcut[1] = N->rcut[2] = RCUT_ZrC;
    N->rcut[3] = RCUT_CC;
    N->neighbormax[0] = N->neighbormax[1] = DEF_NEIGHBORMAX;
    N->max_atom_displacement = DEF_MAX_ATOM_DISPLACEMENT;
    Neighborlist_Recreate (Config_Aapp_to_Alib, ft, ct, &tp, N); Fcr(ft);
    
    Fprintf (ft, "********************* "
             "Before Minimization *********************\n");
    pote_before_relaxation = potential(np,H,s,tp,N,f,stress);
    Fprintf (ft, "The total energy = %f eV, energy/atom = %f eV,\n",
             pote_before_relaxation*UENERGY_IN_EV,
             pote_before_relaxation/np*UENERGY_IN_EV);
    Fprintf (ft, "so the defect formation energy (unrelaxed)\n");
    Fprintf (ft, "= %f - %d x %f = %f eV.\n",
             pote_before_relaxation*UENERGY_IN_EV,
             np, ep_reference*UENERGY_IN_EV,
             (pote_before_relaxation-ep_reference*np)*UENERGY_IN_EV);
    Fprintf (ft, "(energy reference is perfect crystal at zero stress)\n");
    M3InV (H, HI, volume_before_relaxation);
    Fprintf (ft, "The total volume = %.3f A^3, volume/atom = %f A^3,\n",
             volume_before_relaxation*UVOLUME_IN_A3,
             volume_before_relaxation/np*UVOLUME_IN_A3);
    Fprintf (ft, "so the excess volume (unrelaxed)\n");
    Fprintf (ft, "= %f - %d x %f = %f A^3.\n",
             volume_before_relaxation*UVOLUME_IN_A3,
             np, volumep_reference*UVOLUME_IN_A3,
             (volume_before_relaxation-np*volumep_reference)*UVOLUME_IN_A3);
    Fprintf (ft, "(volume reference is perfect crystal at zero stress)\n");
    M3MultiplY (USTRESS_IN_MPA, stress);
    S3fPR (ft, "Initial stress = %M MPa.\n", stress);
    Fprintf (ft, "**********************"
             "****************************************\n\n");
    save_config ("start.cfg");
    
    Fprintf (ft, "Random seed = %d, using which we give each particle\n",
             random_seed);
    cc = perturbation_amplitude * lattice_constant_reference;
    Fprintf (ft, "a random kick of (%g, %g) A in all %d directions.\n\n",
             -cc/2*ULENGTH_IN_A, cc/2*ULENGTH_IN_A, DIMENSION);
    V3ZERO(cm);
    for (dd=i=0; i<np; i++)
    {
        xx[0] = FRANDOM() * cc;
        xx[1] = FRANDOM() * cc;
        xx[2] = FRANDOM() * cc;
        V3mM3 (xx, HI, ss);
        V3AdD (ss, &s[DIMENSION*i]);
        dd += mass[i];
        V3ADDmuL (mass[i], ss, cm); 
    }
    V3DiV (cm, dd);
    for (i=0; i<np; i++)
    {
        V3SuB (&s[DIMENSION*i], cm);
        V3TRIM (&s[DIMENSION*i], &s[DIMENSION*i]);
    }

    Fprintf (ft, "Minimization tolerance = %e.\n\n", minimization_tolerance);
    CG_option.gradient_square_sum_tolerance = minimization_tolerance;
    
    if (save_screen_output)
        Fprintf (ft, "The screen output would be saved on \"%s\".\n",
                 fname_out);
    if (write_config)
        Fprintf (ft, "The relaxed configuration would be saved on \"%s\".\n",
                 fname_write);
    
    Fprintf (ft, "Minimization starts...\n\n");

    if (constant_S)
    {
        VZERO (DIMENSION*(DIMENSION+1)/2, s+DIMENSION*np);
        M3EQV(H,H0);
        zxcgr (free_energy, DIMENSION*np+DIMENSION*(DIMENSION+1)/2, s);
    }
    else zxcgr (free_energy, DIMENSION*np, s);
    
    cr();
    if (CG_state.error) Fprintf (ft, "%s\n", CG_state.error_mesg);
    else Fprintf (ft, "Minimization converged after %d evaluations.\n\n",
                  CG_state.number_of_evaluations);
    report(ft); Fcr(ft);

    fbrk();
    fclose(fp);
    if (save_screen_output)
        printf("The screen output has been saved on \"%s\".\n\n", fname_out);

    printf("Performing last-minute neighborlist consistency check,\n"
           "you may Ctrl-C to interrupt...\n");
    Neighborlist_Check (Config_Aapp_to_Alib, stdout, ct, &tp, N);

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