/********************************************/
/* Program To Fit Finnis-Sinclair Potential */
/*                                          */
/*           V2.0  May 28 2000  Ju Li (MIT) */
/********************************************/

#include "z.h"

double potential_energy (Structure *S, double param[MAX_OP])
{
    int i,j,k,l,nn[MAX_ATOM],ni[MAX_ATOM][MAX_FIRST],swapped;
    double *ds, *dx, RCUT, AA, BB, CC, DD, ff, pote, *DS, *DX,
        *DX_NEXT, H[3][3], K, U, W, ALPHA_ZrC, BETA_ZrC, ALPHA_CZr,
        BETA_CZr, A[MAX_TYPE][MAX_TYPE], B[MAX_TYPE][MAX_TYPE], LAMBDA,
        C[MAX_TYPE][MAX_TYPE], D[MAX_TYPE][MAX_TYPE], sc, cc, cz, ss,
        tt, m[4], xx[3], wk[3], wx[3], RC[MAX_TYPE][MAX_TYPE], h,
        neigh[MAX_ATOM][MAX_FIRST*4], GAMMA;

    A[0][0] = 2.9296875;
    B[0][0] = 2.58787395638939;
    C[0][0] = 1.32308467741935;
    D[0][0] = 4.3672464262563;
    RC[0][0] = 7.0;
    K = 0.1;
    RC[0][1] = RC[1][0] = W = 3.5;
    U = 0.1;
    RC[1][1] = 0.;

    i = 0;
    A[0][1] = A[1][0] = param[i++];
    B[0][1] = B[1][0] = param[i++];
    C[0][1] = C[1][0] = param[i++];
    D[0][1] = D[1][0] = param[i++];
    LAMBDA  = 0;
    ALPHA_ZrC = ALPHA_CZr = param[i++];
    BETA_ZrC  = BETA_CZr  = param[i++];
    GAMMA = 2.;
    
    IZERO (S->n, nn);
    M3MULTIPLY (S->a, S->H, H);
    for (DS=S->ds,DX=S->dx,i=0; i<S->n; i++)
        for (j=i; j<S->n; j++,DS+=3*MAX_IMAGE,DX=DX_NEXT)
        {
            DX_NEXT = DX + 4*MAX_IMAGE;
            RCUT = RC[S->t[i]][S->t[j]];
            swapped = (S->t[i] == C_ATOM) && (S->t[j] == Zr_ATOM);
            for (ds=DS,dx=DX; dx<DX_NEXT; ds+=3,dx+=4)
            {
                V3M3LENGTH (ds, H, dx);
                if (dx[3] >= RCUT) break;
                if (S->t[i] == S->t[j]) continue;
                if (dx[3] >= W) continue;
                if ((nn[i]==MAX_FIRST) || (nn[j]==MAX_FIRST))
                    pe("potential_energy: MAX_FIRST=%d reached.\n",MAX_FIRST);
                V3DIV ( dx, dx[3], &(neigh[i][4*nn[i]]) );
                neigh[i][4*nn[i]+3] = dx[3];
                ni[i][nn[i]++] = j;
                V3DIV ( dx, -dx[3], &(neigh[j][4*nn[j]]) );
                neigh[j][4*nn[j]+3] = dx[3];
                ni[j][nn[j]++] = i;
                if (swapped) V3NEG(dx,dx);
            }
        }

    pote = 0;
    VZERO (S->n, S->X);
    for (DX=S->dx,i=0; i<S->n; i++)
        for (j=i; j<S->n; j++,DX=DX_NEXT)
        {
            DX_NEXT = DX + 4*MAX_IMAGE;
            RCUT = RC[S->t[i]][S->t[j]];
            swapped = (S->t[i] == C_ATOM) && (S->t[j] == Zr_ATOM);
            if (swapped) SWAP(i,j,k);
            AA = A[S->t[i]][S->t[j]];
            BB = B[S->t[i]][S->t[j]];
            CC = C[S->t[i]][S->t[j]];
            DD = D[S->t[i]][S->t[j]];
            for (dx=DX; dx<DX_NEXT; dx+=4)
            {
                if (dx[3] >= RCUT) break;
                if ( S->t[i] == S->t[j] )
                    if ( S->t[i] == Zr_ATOM )
                    {
                        pote += 2 * exp( AA * (BB - dx[3]) +
                                         K / (dx[3] - RCUT) );
                        cc = exp( CC * (DD - dx[3]) + K / (dx[3] - RCUT) );
                        S->X[i] += cc;
                        S->X[j] += cc;
                    }
                    else {}
                else
                {
                    if (dx[3] >= W) continue;
                    pote += 2 * exp( AA * (BB - dx[3]) +
                                     U / (dx[3] - W) );
                    for (sc=k=0; k<nn[i]; k++)
                    { /* screening at Zr */
                        tt = 1 + V3DOT(dx, &(neigh[i][4*k])) / dx[3];
                        if (tt <= 0) continue;
                        if ((tt > 2-TINY)&&(neigh[i][4*k+3]==dx[3])) continue;
                        sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                            exp ( GAMMA * ( U / (neigh[i][4*k+3] - W) -
                                            CC * neigh[i][4*k+3] ) );
                    }
                    for (k=0; k<nn[j]; k++)
                    { /* screening at C */
                        tt = 1 - V3DOT(dx, &(neigh[j][4*k])) / dx[3];
                        if (tt <= 0) continue;
                        if ((tt > 2-TINY)&&(neigh[j][4*k+3]==dx[3])) continue;
                        sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                            exp ( GAMMA * ( U / (neigh[j][4*k+3] - W) -
                                            CC * neigh[j][4*k+3] ) );
                    }
                    sc = pow ( sc , 1. / GAMMA ) /
                        exp( U / (dx[3] - W) - CC * dx[3] );
                    h = exp( CC * (DD - dx[3]) + U / (dx[3] - W) - sc );
                    S->X[i] += h * h;
                    S->X[j] += h * h;
                    pote -= LAMBDA * h;
                }
            }
            if (swapped) SWAP(i,j,k);
        }
    for (i=0; i<S->n; i++) pote -= (S->X[i]=sqrt(S->X[i]));
    if (!S->calculate_force) return (pote);

    VZERO (3*S->n, S->f);
    for (DX=S->dx,i=0; i<S->n; i++)
        for (j=i; j<S->n; j++,DX=DX_NEXT)
        {
            DX_NEXT = DX + 4*MAX_IMAGE;
            RCUT = RC[S->t[i]][S->t[j]];
            swapped = (S->t[i] == C_ATOM) && (S->t[j] == Zr_ATOM);
            if (swapped) SWAP(i,j,k);
            AA = A[S->t[i]][S->t[j]];
            BB = B[S->t[i]][S->t[j]];
            CC = C[S->t[i]][S->t[j]];
            DD = D[S->t[i]][S->t[j]];
            for (dx=DX; dx<DX_NEXT; dx+=4)
            {
                if (dx[3] >= RCUT) break;
                if ( S->t[i] == S->t[j] )
                    if ( S->t[i] == Zr_ATOM )
                    {
                        ff = ( 2 * ( AA + K / (dx[3] - RCUT) / (dx[3] - RCUT) )
                               * exp( AA * (BB - dx[3]) + K / (dx[3] - RCUT) )
                               - ( CC + K / (dx[3] - RCUT) / (dx[3] - RCUT) )
                               * exp( CC * (DD - dx[3]) + K / (dx[3] - RCUT) )
                               * ( 0.5 / S->X[i] + 0.5 / S->X[j] ) ) / dx[3];
                        V3ADDmuL (ff, dx, &(S->f[3*j]));
                        V3SUBmuL (&S->f[3*i], ff, dx);
                    }
                    else {}
                else
                {
                    if (dx[3] >= W) continue;
                    for (sc=k=0; k<nn[i]; k++)
                    { /* screening at Zr */
                        tt = 1 + V3DOT(dx, &(neigh[i][4*k])) / dx[3];
                        if (tt <= 0) continue;
                        if ((tt > 2-TINY)&&(neigh[i][4*k+3]==dx[3])) continue;
                        sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                            exp ( GAMMA * ( U / (neigh[i][4*k+3] - W) -
                                            CC * neigh[i][4*k+3] ) );
                    }
                    for (k=0; k<nn[j]; k++)
                    { /* screening at C */
                        tt = 1 - V3DOT(dx, &(neigh[j][4*k])) / dx[3];
                        if (tt <= 0) continue;
                        if ((tt > 2-TINY)&&(neigh[j][4*k+3]==dx[3])) continue;
                        sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                            exp ( GAMMA * ( U / (neigh[j][4*k+3] - W) -
                                            CC * neigh[j][4*k+3] ) );
                    }
                    sc = pow ( sc , 1. / GAMMA ) /
                        exp( U / (dx[3] - W) - CC * dx[3] );
                    h = exp( CC * (DD - dx[3]) + U / (dx[3] - W) - sc );
                    cc = h * h / S->X[i] + h * h / S->X[j] + LAMBDA * h;
                    ff = ( 2 * ( AA + U / (dx[3] - W) / (dx[3] - W) ) *
                           exp( AA * (BB - dx[3]) + U / (dx[3] - W) ) -
                           ( CC + U / (dx[3] - W) / (dx[3] - W) )
                           * cc * ( 1 + sc ) ) / dx[3];
                    V3ADDmuL (ff, dx, &(S->f[3*j]));
                    V3SUBmuL (&S->f[3*i], ff, dx);
                    if (sc < TINY) continue;
                    V3DIV ( dx, dx[3], xx );
                    for (k=0; k<nn[i]; k++)
                    {
                        tt = 1 + V3DOT(xx, &(neigh[i][4*k]));
                        if (tt <= 0) continue;
                        if ((tt > 2-TINY)&&(neigh[i][4*k+3]==dx[3])) continue;
                        ss = cc / pow(sc, GAMMA-1) / GAMMA *
                            pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                            exp ( GAMMA * ( U / (neigh[i][4*k+3] - W) -
                                            CC * neigh[i][4*k+3] -
                                            U / (dx[3] - W) +
                                            CC * dx[3] ) );
                        ff = ss * GAMMA * ( CC + U /
                                            (neigh[i][4*k+3] - W) /
                                            (neigh[i][4*k+3] - W) );
                        V3ADDmuL (ff, &(neigh[i][4*k]), &(S->f[3*ni[i][k]]) );
                        V3SUBmuL (&S->f[3*i], ff, &(neigh[i][4*k]));
                        ss *= BETA_ZrC / tt;
                        V3CROSS ( xx, &(neigh[i][4*k]), m );
                        V3CROSS ( m, &(neigh[i][4*k]), wk );
                        V3CROSS ( xx, m, wx );
                        tt = ss / neigh[i][4*k+3];
                        V3ADDmuL (tt, wk, &(S->f[3*ni[i][k]]) );
                        V3SUBmuL (&S->f[3*i], tt, wk);
                        tt = ss / dx[3];
                        V3ADDmuL (tt, wx, &(S->f[3*j]) );
                        V3SUBmuL (&S->f[3*i], tt, wx);
                    }
                    for (k=0; k<nn[j]; k++)
                    {
                        tt = 1 - V3DOT(xx, &(neigh[j][4*k]));
                        if (tt <= 0) continue;
                        if ((tt > 2-TINY)&&(neigh[j][4*k+3]==dx[3])) continue;
                        ss = cc / pow(sc, GAMMA-1) / GAMMA *
                            pow ( tt /  ALPHA_CZr, BETA_CZr ) *
                            exp ( GAMMA * ( U / (neigh[j][4*k+3] - W) -
                                            CC * neigh[j][4*k+3] -
                                            U / (dx[3] - W) +
                                            CC * dx[3] ) );
                        ff = ss * GAMMA * ( CC + U /
                                            (neigh[j][4*k+3] - W) /
                                            (neigh[j][4*k+3] - W) );
                        V3ADDmuL (ff, &(neigh[j][4*k]), &(S->f[3*ni[j][k]]) );
                        V3SUBmuL (&(S->f[3*j]), ff, &(neigh[j][4*k]));
                        ss *= BETA_CZr / tt;
                        V3CROSS ( &(neigh[j][4*k]), xx, m );
                        V3CROSS ( m, &(neigh[j][4*k]), wk );
                        V3CROSS ( m, xx, wx );
                        tt = ss / neigh[j][4*k+3];
                        V3ADDmuL (tt, wk, &(S->f[3*ni[j][k]]) );
                        V3SUBmuL (&(S->f[3*j]), tt, wk);
                        tt = ss / dx[3];
                        V3ADDmuL (tt, wx, &(S->f[3*i]) );
                        V3SUBmuL (&(S->f[3*j]), tt, wx);
                    }
                }
            }
            if (swapped) SWAP(i,j,k);
        }
    return (pote);
} /* end potential_energy() */


static double this_pote (double lattice_constant)
{
    this_structure->a = lattice_constant;
    return (potential_energy(this_structure, param));
} /* end this_pote() */


#define error_inc(x,target,weight) \
  ((weight)*pow(fabs(rerr(x,target)),exponent))

/* the error function */
double fe (double guess[MAX_OP])
{
    int i, j;
    double a, b, c, pote, error, B;
    memcpy (param, guess, MAX_OP*sizeof(double));
    error = 0.;
    for (i=0; i<MAX_XTAL_STRUCTURE; i++)
    {
        this_structure = &my[i];
        pote = onemin (this_pote, lattice_constant_range[i][0],
                     lattice_constant_range[i][1], EPS, &a);
        this_structure->pote = pote;
        this_structure->a = a;
        error += error_inc(a, xtal_target[i][XTAL_EQ_A],
                           xtal_weight[i][XTAL_EQ_A]);
        if (i == NaCl)
        {
            c44 = C44_0 (this_structure);
            bulk_modulus = Bulk_Modulus (this_structure);
            c = C11_C12 (this_structure);
            c11 = (3*bulk_modulus + 2*c) / 3;
            c12 = c11 - c;
            error +=
                error_inc(bulk_modulus, bulk_modulus_target,
                          bulk_modulus_weight) +
                error_inc(c11,c11_target,c11_weight) +
                error_inc(c12,c12_target,c12_weight) +
                error_inc(c44,c44_target,c44_weight);
            error += error_inc( -pote, NaCl_cohesive_energy_target,
                                xtal_weight[NaCl][XTAL_EQ_POTE] );
        }
        else error += error_inc( my[i].pote - my[i-1].pote,
                                 xtal_target[i][XTAL_EQ_POTE] -
                                 xtal_target[i-1][XTAL_EQ_POTE],
                                 xtal_weight[i][XTAL_EQ_POTE] );
    }
    my[ZrVAC].a = my[NaCl].a;
    my[ZrVAC].pote = potential_energy (&my[ZrVAC], param);
    my[CVAC].a  = my[NaCl].a;
    my[CVAC].pote = potential_energy (&my[CVAC], param);
    schottky_energy = my[ZrVAC].pote + my[CVAC].pote -
        (my[ZrVAC].n + my[CVAC].n) * my[NaCl].pote / my[NaCl].n;
    error += error_inc( schottky_energy, schottky_energy_target,
                        schottky_energy_weight );
    /* pure Zr */
    this_structure = &my[PureZr];
    this_structure->pote = onemin (this_pote, 1, 10, EPS, &a);
    this_structure->a = a;
    error += error_inc( CUBE(my[PureZr].a) * my[PureZr].volume / my[PureZr].n,
                        PureZr_atomic_volume_target,
                        PureZr_atomic_volume_weight );
    B = Bulk_Modulus (this_structure);
    error += error_inc( B, PureZr_bulk_modulus_target,
                        PureZr_bulk_modulus_weight );
    error += error_inc( -this_structure->pote / this_structure->n,
                        PureZr_cohesive_energy_target,
                        PureZr_cohesive_energy_weight );
    /* pure C */
    /* this_structure = &my[PureC]; */
    /* this_structure->pote = onemin (this_pote, 1, 6, EPS, &a); */
    /* this_structure->a = a; */
    /* error += error_inc( CUBE(my[PureC].a)/4, PureC_atomic_volume_target, */
    /* PureC_atomic_volume_weight ); */
    /* B = Bulk_Modulus (this_structure); */
    /* error += error_inc( B, PureC_bulk_modulus_target, */
    /* PureC_bulk_modulus_weight ); */
    /* error += error_inc( -this_structure->pote, PureC_cohesive_energy_target, */
    /* PureC_cohesive_energy_weight ); */
    my[PureC].n = 1;
    my[PureC].pote = -7.43;
    /* ZrC heat of formation */
    error += error_inc ( my[PureZr].pote / my[PureZr].n +
                         my[PureC].pote / my[PureC].n - my[NaCl].pote,
                         heat_of_formation_target, heat_of_formation_weight );
    /* ZrC atomic force constants */
    my[ZrC_CDIS].a = my[NaCl].a;
    potential_energy (&my[ZrC_CDIS], param);
    for (i=0; i<my[ZrC_CDIS].n; i++)
        error += error_inc ( my[ZrC_CDIS].f[3*i] * UENERGY_IN_J /
                             SQUARE(ULENGTH_IN_M) / DIS_DEL / my[ZrC_CDIS].a,
                             ZrC_CDIS_fx_target[i], ZrC_CDIS_fx_weight[i] );
    my[ZrC_ZrDIS].a = my[NaCl].a;
    potential_energy (&my[ZrC_ZrDIS], param);
    for (i=0; i<my[ZrC_ZrDIS].n; i++)
        error += error_inc ( my[ZrC_ZrDIS].f[3*i] * UENERGY_IN_J /
                             SQUARE(ULENGTH_IN_M) / DIS_DEL / my[ZrC_ZrDIS].a,
                             ZrC_ZrDIS_fx_target[i], ZrC_ZrDIS_fx_weight[i] );
    error = pow(error, 1./exponent);
    if ( (kb % print_freq == 0) || just_print_it )
    {
        print_param (stdout);
        print_property (stdout);
        printf ("step %d:  error = %lf\n\n", kb, error);
    }
    kb++;
    return (error);
} /* end fe() */


/* to fool my configuration viewer */
#define MAGNIFICATION  3
void wrap_up (int sig)
{
    int i, j;
    FILE *rep;
    double r0, e0, error;
    double x0[3]={0,0,0}, H0[3][3]={{1,0,0},{0,1,0},{0,0,1}};
    if (sig != SIG_SUCCESS) printf ("Interrupt received, wrapping up ...\n\n");
    else printf ("Optimization terminates successfully ...\n\n");
    printf ("---------------------------------------------------------\n");
    just_print_it = 1;
    error = fe (param);
    rep = fopen(REP_FILENAME, "w");
    fprintf (rep, "Using binary Morse potential model\n");
    fprintf (rep, "V(r) = D0 * ( exp(2 * BETA * (R0-r)) - "
             "2 * exp(BETA * (R0-r)) ),\n");
    fprintf (rep, "and cutoff scheme W(r) = V(r) - "
             "V(RC) - V'(RC) * (r - RC),\n");
    fprintf (rep, "we arrive at optimized parameter set "
             "(upper_bound  fitted  lower_bound):\n");
    print_param (rep);
    print_property (rep);
    fprintf (rep, "step %d:  error = %lf\n\n", kb-1, error);
    for (i=0; i<MAX_STRUCTURE; i++)
    {
        fprintf (rep, "Saving atomic configuration in \"%s\".\n",
                 my[i].name);
        printf ("Saving atomic configuration in \"%s\".\n", my[i].name);
        if (i==PureZr)
            write_config (my+i, x0, my[i].H, MAGNIFICATION, my[i].name);
        else write_config (my+i, x0, H0, MAGNIFICATION, my[i].name);
    }
    fprintf (rep, "Plotting cohesive energy curves in \"%s\".\n", "coh.m");
    printf ("Plotting cohesive energy curves in \"%s\".\n", "coh.m");
    plot_cohesive_energy_curves ("coh.m");
    printf ("The above summary is saved in \"%s\".\n", REP_FILENAME);
    printf ("---------------------------------------------------------\n");
    fclose(rep);
    exit(0);
} /* end wrap_up() */


void construct_neighbor_list (Structure *S)
{
    double *ds, dx[3], *dr2, *DS;
    int i, j, m, n, maxnc[3], nc[3], *idx;
    M3InV (S->H, S->HI, S->volume);
    M3rows_to_cover_sphere (S->H, cbrt(MAX_IMAGE*S->volume*3./4/PI), maxnc);
    V3aDd (maxnc, 4);
    m = total_replica(maxnc) * S->n;
    MALLOC ( construct_neighbor_list,  ds, 3*m,  double );
    MALLOC ( construct_neighbor_list, dr2,   m,  double );
    MALLOC ( construct_neighbor_list, idx,   m,  int );
    S->mass = 0;
    for (i=0; i<S->n; i++) S->mass += ATOM_MASS_IN_G(Z[S->t[i]]);
    for (DS=S->ds,i=0; i<S->n; i++)
        for (j=i; j<S->n; j++)
        {
            for (m=0,nc[0]=-maxnc[0]; nc[0]<=maxnc[0]; nc[0]++)
                for (nc[1]=-maxnc[1]; nc[1]<=maxnc[1]; nc[1]++)
                    for (nc[2]=-maxnc[2]; nc[2]<=maxnc[2]; nc[2]++)
                    {
                        if (i==j)
                            if (nc[0] < 0) continue;
                            else if (nc[0] == 0)
                                if (nc[1] < 0) continue;
                                else if (nc[1] == 0)
                                    if (nc[2] <= 0) continue;
                        V3ADDSUB (nc, &S->s[3*j], &S->s[3*i], &ds[3*m]);
                        V3mM3 (&ds[3*m], S->H, dx);
                        dr2[m] = V3LENGTH2 (dx);
                        m++;
                    }
            qsort_numerical_recipes ( m, dr2, idx, USE_NEW_IDX );
            if (m < MAX_IMAGE)
            {
                for (n=0; n<m; n++,DS+=3)
                    V3EQV( &ds[3*idx[n]], DS );
                V3EYE( HUGE, DS );
                DS += 3 * (MAX_IMAGE - m);
            }
            else for (n=0; n<MAX_IMAGE; n++,DS+=3)
                V3EQV( &ds[3*idx[n]], DS );
        }
    free (idx);
    free (dr2);
    free (ds);
    return;
} /* end construct_neighbor_list() */


int main (int argc, char *argv[])
{
    int i, j, n_op, ndiv;
    double guess[MAX_OP], total_weight, del;

    my[NaCl].n       = 2;
    my[ZnS].n        = 2;
    my[CsCl].n       = 2;
    my[ZrVAC].n      = 7;
    my[CVAC].n       = 7;
    my[PureZr].n     = 2;
    my[PureC].n      = 1;
    my[ZrC_CDIS].n   = 8;
    my[ZrC_ZrDIS].n  = 8;
    my[ZrC_TEST].n   = 8;

    for (i=0; i<MAX_STRUCTURE; i++)
    {
        if (my[i].n > MAX_ATOM)
            pe ("Structure %d exceeds MAX_ATOM = %d\n", i, MAX_ATOM);
        MALLOC (main, my[i].t,    my[i].n, int);
        MALLOC (main, my[i].s,  3*my[i].n, double);
        MALLOC (main, my[i].ds, 3*my[i].n*(my[i].n+1)/2*MAX_IMAGE, double);
        MALLOC (main, my[i].dx, 4*my[i].n*(my[i].n+1)/2*MAX_IMAGE, double);
        MALLOC (main, my[i].X,    my[i].n, double);
        MALLOC (main, my[i].f,  3*my[i].n, double);
    }

    my[i=NaCl].name = "NaCl";
    M3ASSIGN ( 0,.5,.5, .5,0,.5, .5,.5,0, my[i].H );
    my[i].t[j=0] = Zr_ATOM;
    V3ZERO ( &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( -.5,.5,.5, &my[i].s[3*j] );

    my[i=ZnS].name = "ZnS";
    M3ASSIGN ( 0,.5,.5, .5,0,.5, .5,.5,0, my[i].H );
    my[i].t[j=0] = Zr_ATOM;
    V3ZERO ( &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( .25,.25,.25, &my[i].s[3*j] );
    
    my[i=CsCl].name = "CsCl";
    M3IDENTITY ( my[i].H );
    my[i].t[j=0] = Zr_ATOM;
    V3ZERO ( &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( .5,.5,.5, &my[i].s[3*j] );

    my[i=ZrVAC].name = "ZrVAC";
    M3IDENTITY ( my[i].H );
    my[i].t[j=0] = Zr_ATOM;
    V3ASSIGN ( .5,.5,0, &my[i].s[3*j] );
    my[i].t[++j] = Zr_ATOM;
    V3ASSIGN ( 0,.5,.5, &my[i].s[3*j] );
    my[i].t[++j] = Zr_ATOM;
    V3ASSIGN ( .5,0,.5, &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( .5,0,0,  &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( 1,.5,0,  &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( .5,.5,.5, &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( 1,0,.5,   &my[i].s[3*j] );

    my[i=CVAC].name = "CVAC";
    M3EQV (my[ZrVAC].H, my[i].H);
    for (j=0; j<my[i].n; j++)
        my[i].t[j] = Zr_ATOM + C_ATOM - my[ZrVAC].t[j];
    memcpy (my[i].s, my[ZrVAC].s, 3*my[i].n*sizeof(double));

    my[i=PureZr].name = "PureZr";
    M3ASSIGN ( 1,0,0, -.5,SQRT3/2,0, 0,0,1.6298, my[i].H );
    my[i].t[j=0] = Zr_ATOM;
    V3ZERO ( &my[i].s[3*j] );
    my[i].t[++j] = Zr_ATOM;
    V3ASSIGN ( 2*THIRD,THIRD,HALF, &my[i].s[3*j] );

    my[i=PureC].name = "PureC";
    M3ASSIGN ( 0,.5,.5, .5,0,.5, .5,.5,0, my[i].H );
    my[i].t[j=0] = C_ATOM;
    V3ZERO ( &my[i].s[3*j] );

    my[i=ZrC_CDIS].name = "ZrC_CDIS";
    my[i].calculate_force = TRUE;
    M3IDENTITY ( my[i].H );
    my[i].t[j=0] = Zr_ATOM;
    V3ZERO ( &my[i].s[3*j] );
    my[i].t[++j] = Zr_ATOM;
    V3ASSIGN ( 0,.5,.5,  &my[i].s[3*j] );
    my[i].t[++j] = Zr_ATOM;
    V3ASSIGN ( .5,0,.5,  &my[i].s[3*j] );
    my[i].t[++j] = Zr_ATOM;
    V3ASSIGN ( .5,.5,0,  &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( .5,.5,.5, &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( .5-DIS_DEL,0,0, &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( 0,.5,0,   &my[i].s[3*j] );
    my[i].t[++j] = C_ATOM;
    V3ASSIGN ( 0,0,.5,   &my[i].s[3*j] );

    my[i=ZrC_ZrDIS].name = "ZrC_ZrDIS";
    my[i].calculate_force = TRUE;
    M3EQV (my[ZrC_CDIS].H, my[i].H);
    for (j=0; j<my[i].n; j++)
        my[i].t[j] = Zr_ATOM + C_ATOM - my[ZrC_CDIS].t[j];
    memcpy (my[i].s, my[ZrC_CDIS].s, 3*my[i].n*sizeof(double));

    my[ZrC_TEST].name = "ZrC_TEST";
    my[ZrC_TEST].calculate_force = TRUE;
    M3EQV ( my[i].H , my[ZrC_TEST].H );
    TimeRandomize();
    for (j=0; j<my[ZrC_TEST].n; j++)
    {
        my[ZrC_TEST].t[j] = my[i].t[j];
        my[ZrC_TEST].s[3*j]   = my[i].s[3*j]   + TEST_DEL * FRANDOM();
        my[ZrC_TEST].s[3*j+1] = my[i].s[3*j+1] + TEST_DEL * FRANDOM();
        my[ZrC_TEST].s[3*j+2] = my[i].s[3*j+2] + TEST_DEL * FRANDOM();
    }

    for (i=0; i<MAX_STRUCTURE; i++) construct_neighbor_list(my+i);

    scanf ("[%d] Fitting Parameters (Upper_bound  Guess  Lower_bound):\n",
           &n_category);
    if (n_category > MAX_CATEGORY)
        pe("number of categories (%d) > MAX_CATEGORY (%d)\n",
           n_category, MAX_CATEGORY);
    for (n_op=i=0; i<n_category; i++)
    {
        scanf ("[%d]", &category_ops[i]);
        for (j=0; (category_name[i][j]=(char)getc(stdin))!='\n'; j++);
        category_name[i][j]=(char)0;
        if (n_op+category_ops[i] > MAX_OP)
            pe("number of parameters > MAX_OP (%d)\n", MAX_OP);
        for (j=0; j<category_ops[i]; j++)
            scanf ("%lf", &bound[1][n_op+j]); scanf ("\n");
        for (j=0; j<category_ops[i]; j++)
            scanf ("%lf", &guess[n_op+j]); scanf ("\n");
        for (j=0; j<category_ops[i]; j++)
            scanf ("%lf", &bound[0][n_op+j]); scanf ("\n");
        n_op += category_ops[i];
    }
    
    scanf ("\nCrystal equilibrium lattice constant [A] range, weight:\n");
    scanf ("%lf %lf %lf\n", &lattice_constant_range[NaCl][0],
           &lattice_constant_range[NaCl][1], &xtal_weight[NaCl][XTAL_EQ_A]);
    scanf ("%lf %lf %lf\n", &lattice_constant_range[ZnS][0],
           &lattice_constant_range[ZnS][1], &xtal_weight[ZnS][XTAL_EQ_A]);
    scanf ("%lf %lf %lf\n", &lattice_constant_range[CsCl][0],
           &lattice_constant_range[CsCl][1], &xtal_weight[CsCl][XTAL_EQ_A]);

    scanf ("\nNaCl cohesive energy weight:\n%lf\n",
           &xtal_weight[NaCl][XTAL_EQ_POTE]);
    scanf ("\nZnS/NaCl equilibrium energy difference weight:\n%lf\n",
           &xtal_weight[ZnS][XTAL_EQ_POTE]);
    scanf ("\nCsCl/ZnS equilibrium energy difference weight:\n%lf\n",
           &xtal_weight[CsCl][XTAL_EQ_POTE]);
    total_weight = 0.;
    for (i=0; i<MAX_XTAL_STRUCTURE; i++)
        for (j=0; j<MAX_XTAL_PROPERTY; j++)
            total_weight += xtal_weight[i][j];

    scanf ("\nBulk modulus, C11, C12, C44 weights:\n");
    scanf ("%lf %lf %lf %lf\n", &bulk_modulus_weight, &c11_weight,
           &c12_weight, &c44_weight);
    total_weight += bulk_modulus_weight + c11_weight + c12_weight + c44_weight;

    scanf ("\nSchottky energy weight:\n");
    scanf ("%lf\n", &schottky_energy_weight);
    total_weight += schottky_energy_weight;

    scanf ("\nPure Zr atomic volume, bulk modulus "
           "and cohesive energy weights:\n");
    scanf ("%lf %lf %lf\n", &PureZr_atomic_volume_weight,
           &PureZr_bulk_modulus_weight, &PureZr_cohesive_energy_weight);
    total_weight += PureZr_atomic_volume_weight + PureZr_bulk_modulus_weight +
        PureZr_cohesive_energy_weight;

    scanf ("\nPure C atomic volume, bulk modulus "
           "and cohesive energy weights:\n");
    scanf ("%lf %lf %lf\n", &PureC_atomic_volume_weight,
           &PureC_bulk_modulus_weight, &PureC_cohesive_energy_weight);
    total_weight += PureC_atomic_volume_weight + PureC_bulk_modulus_weight +
        PureC_cohesive_energy_weight;

    scanf ("\nZrC heat of formation weight:\n");
    scanf ("%lf\n", &heat_of_formation_weight);
    total_weight += heat_of_formation_weight;

    scanf ("\nPerturbed ZrC_CDIS atomic forces weights:\n");
    for (i=0; i<my[ZrC_CDIS].n; i++)
    {
        scanf ("%lf", &ZrC_CDIS_fx_weight[i]);
        total_weight += ZrC_CDIS_fx_weight[i];
    }

    scanf ("\nPerturbed ZrC_ZrDIS atomic forces weights:\n");
    for (i=0; i<my[ZrC_ZrDIS].n; i++)
    {
        scanf ("%lf", &ZrC_ZrDIS_fx_weight[i]);
        total_weight += ZrC_ZrDIS_fx_weight[i];
    }

    /* normalize all weights */
    for (i=0; i<MAX_XTAL_STRUCTURE; i++)
        for (j=0; j<MAX_XTAL_PROPERTY; j++)
            xtal_weight[i][j] /= total_weight;
    bulk_modulus_weight /= total_weight;
    c11_weight /= total_weight;
    c12_weight /= total_weight;
    c44_weight /= total_weight;
    schottky_energy_weight /= total_weight;
    PureZr_atomic_volume_weight /= total_weight;
    PureZr_bulk_modulus_weight /= total_weight;
    PureZr_cohesive_energy_weight /= total_weight;
    PureC_atomic_volume_weight /= total_weight;
    PureC_bulk_modulus_weight /= total_weight;
    PureC_cohesive_energy_weight /= total_weight;
    heat_of_formation_weight /= total_weight;
    for (i=0; i<my[ZrC_CDIS].n; i++) ZrC_CDIS_fx_weight[i] /= total_weight;
    for (i=0; i<my[ZrC_ZrDIS].n; i++) ZrC_ZrDIS_fx_weight[i] /= total_weight;

    scanf ("\nRelative error module exponent:\n%lf\n", &exponent);
    scanf ("\nMina del, print-out frequency, maximum steps:\n");
    scanf ("%lf %d %d\n", &del, &print_freq, &maxkb);

    if (maxkb > 0)
    {
        just_print_it = 0;
        signal (SIGINT, wrap_up);
        ndiv = 100;
        kb = 0;
        FORTRAN_SYMBOL(minaf77)
            (fe, &n_op, &ndiv, &del, bound[0], bound[1], guess, param);
    }
    else memcpy (param, guess, MAX_OP*sizeof(double));

    wrap_up (SIG_SUCCESS);
    return (0);
} /* end main() */