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

#include "z.h"

void print_param (FILE *out)
{
    int i, j, m_op;
    for (m_op=i=0; i<n_category; i++)
    {
        fprintf (out, "[%d]%s\n", category_ops[i], &category_name[i][0]);
        for (j=0; j<category_ops[i]; j++)
            fprintf (out, "%16.14f ", bound[1][m_op+j]);
        fprintf (out, "\n");
        for (j=0; j<category_ops[i]; j++)
            fprintf (out, "%16.14f ", param[m_op+j]);
        fprintf (out, "\n");
        for (j=0; j<category_ops[i]; j++)
            fprintf (out, "%16.14f ", bound[0][m_op+j]);
        fprintf (out, "\n");
        m_op += category_ops[i];
    }
    fprintf (out, "\n");
    return;
} /* end print_param() */


#define rho(mass,volume,a) ((mass)/(volume)/CUBE(a*1E-8))
void print_property (FILE *out)
{
    int i;
    double B, dE, DX[3], DS[3];
    fprintf (out, "At 0K in equilibrium, the current model predicts\n\n");
    for (i=0; i<MAX_XTAL_STRUCTURE; i++)
    {
        fprintf (out, "%s lattice constant = %f [%f] A (%.2f%%)\n",
                 my[i].name, my[i].a, xtal_target[i][XTAL_EQ_A],
                 100*rerr(my[i].a, xtal_target[i][XTAL_EQ_A]));
        fprintf (out, " -> density = %f [%f] g/cm^3 (%.2f%%)\n",
                 rho(my[i].mass, my[i].volume, my[i].a),
                 rho(my[i].mass, my[i].volume, xtal_target[i][XTAL_EQ_A]),
                 100*rerr(rho(my[i].mass, my[i].volume, my[i].a),
                          rho(my[i].mass, my[i].volume,
                              xtal_target[i][XTAL_EQ_A])));
        fprintf (out, "cohesive energy = %f [%f] eV/unit cell (%.2f%%);\n",
                 -my[i].pote, NaCl_cohesive_energy_target +
                 xtal_target[NaCl][XTAL_EQ_POTE]-xtal_target[i][XTAL_EQ_POTE],
                 100*rerr(-my[i].pote,  NaCl_cohesive_energy_target+
                          xtal_target[NaCl][XTAL_EQ_POTE] -
                          xtal_target[i][XTAL_EQ_POTE] ) );
        fprintf (out, "\n");
    }

    fprintf (out, "ZnS/NaCl separation = %f [%f] eV (%.2f%%);\n",
             my[ZnS].pote - my[NaCl].pote,
             xtal_target[ZnS][XTAL_EQ_POTE] - xtal_target[NaCl][XTAL_EQ_POTE],
             100*rerr(my[ZnS].pote - my[NaCl].pote,
                      xtal_target[ZnS][XTAL_EQ_POTE] -
                      xtal_target[NaCl][XTAL_EQ_POTE]));
    fprintf (out, "CsCl/ZnS separation = %f [%f] eV (%.2f%%);\n\n",
             my[CsCl].pote - my[ZnS].pote,
             xtal_target[CsCl][XTAL_EQ_POTE] - xtal_target[ZnS][XTAL_EQ_POTE],
             100*rerr(my[CsCl].pote - my[ZnS].pote,
                      xtal_target[CsCl][XTAL_EQ_POTE] -
                      xtal_target[ZnS][XTAL_EQ_POTE]));

    fprintf (out, "Bulk modulus = %f [%f] GPa (%.2f%%);\n",
             bulk_modulus, bulk_modulus_target, 100*rerr
             (bulk_modulus, bulk_modulus_target));
    fprintf (out, "C11 = %f [%f] GPa (%.2f%%);\n",
             c11, c11_target, 100*rerr(c11, c11_target));
    fprintf (out, "C12 = %f [%f] GPa (%.2f%%);\n",
             c12, c12_target, 100*rerr(c12, c12_target));
    fprintf (out, "C44 = %f [%f] GPa (%.2f%%);\n",
             c44, c44_target, 100*rerr(c44, c44_target));
    fprintf (out, "\n");
    fprintf (out, "Schottky energy = %f [%f] eV (%.2f%%);\n\n",
             schottky_energy, schottky_energy_target,
             100*rerr(schottky_energy, schottky_energy_target));

    this_structure = &my[PureZr];
    fprintf (out, "hcp Zr lattice constant = %f A;\n", this_structure->a);
    fprintf (out, "Zr atomic volume = %f [%f] A^3 (%.2f%%);\n",
             CUBE(this_structure->a) * this_structure->volume /
             this_structure->n, PureZr_atomic_volume_target, 100*rerr
             ( CUBE(this_structure->a) * this_structure->volume /
               this_structure->n, PureZr_atomic_volume_target));
    B = Bulk_Modulus (this_structure);
    fprintf (out, "Zr bulk modulus = %f [%f] GPa (%.2f%%);\n",
             B, PureZr_bulk_modulus_target, 100*
             rerr(B, PureZr_bulk_modulus_target));
    /* printf ("C11 - C12 = %f GPa;\n", C11_C12 (this_structure)); */
    /* printf ("C44 = %f GPa;\n", C44_0 (this_structure)); */
    fprintf (out, "hcp Zr cohesive energy = %f [%f] eV/atom (%.2f%%);\n\n",
             -this_structure->pote / this_structure->n,
             PureZr_cohesive_energy_target, 100*
             rerr(-this_structure->pote / this_structure->n,
                  PureZr_cohesive_energy_target));

    this_structure = &my[PureC];
    fprintf (out, "fcc C lattice constant = %f A;\n", this_structure->a);
    fprintf (out, "C atomic volume = %f [%f] A^3 (%.2f%%);\n",
             CUBE(this_structure->a) * this_structure->volume /
             this_structure->n, PureC_atomic_volume_target, 100*rerr
             ( CUBE(this_structure->a) * this_structure->volume /
               this_structure->n, PureC_atomic_volume_target));
    B = Bulk_Modulus (this_structure);
    fprintf (out, "C bulk modulus = %f [%f] GPa (%.2f%%);\n",
             B, PureC_bulk_modulus_target, 100*
             rerr(B, PureC_bulk_modulus_target));
    fprintf (out, "fcc C cohesive energy = %f [%f] eV/atom (%.2f%%);\n\n",
             -this_structure->pote / this_structure->n,
             PureC_cohesive_energy_target, 100*
             rerr(-this_structure->pote / this_structure->n,
                  PureC_cohesive_energy_target));

    fprintf (out, "ZrC heat of formation = %f [%f] eV/pair (%.2f%%);\n\n",
             my[PureZr].pote / my[PureZr].n + my[PureC].pote / my[PureC].n
             - my[NaCl].pote, heat_of_formation_target, 100*rerr
             (my[PureZr].pote / my[PureZr].n + my[PureC].pote / my[PureC].n
              - my[NaCl].pote, heat_of_formation_target));

    fprintf (out, "In configuration \"%s\" (lattice constant [A] = %f),\n",
             my[ZrC_CDIS].name, my[ZrC_CDIS].a);
    fprintf (out, "the atomic coordinates and force constants [N/m] are\n");
    for (i=0; i<my[ZrC_CDIS].n; i++)
        fprintf (out, "%2s: s= %.2f %.2f %.2f,\tfx= %8.5f [%8.5f]"
                 " (%7.2f%%);\n", ATOM_SYMBOL(Z[my[ZrC_CDIS].t[i]]),
                 my[ZrC_CDIS].s[3*i], my[ZrC_CDIS].s[3*i+1],
                 my[ZrC_CDIS].s[3*i+2], my[ZrC_CDIS].f[3*i] * UENERGY_IN_J /
                 SQUARE(ULENGTH_IN_M) / DIS_DEL / my[ZrC_CDIS].a,
                 ZrC_CDIS_fx_target[i], 100*rerr
                 (my[ZrC_CDIS].f[3*i] * UENERGY_IN_J /
                  SQUARE(ULENGTH_IN_M) / DIS_DEL / my[ZrC_CDIS].a,
                  ZrC_CDIS_fx_target[i]) );
    fprintf (out, "\n");

    fprintf (out, "In configuration \"%s\" (lattice constant [A] = %f),\n",
             my[ZrC_ZrDIS].name, my[ZrC_ZrDIS].a);
    fprintf (out, "the atomic coordinates and force constants [N/m] are\n");
    for (i=0; i<my[ZrC_ZrDIS].n; i++)
        fprintf (out, "%2s: s= %.2f %.2f %.2f,\tfx= %8.5f [%8.5f]"
                 " (%7.2f%%);\n", ATOM_SYMBOL(Z[my[ZrC_ZrDIS].t[i]]),
                 my[ZrC_ZrDIS].s[3*i], my[ZrC_ZrDIS].s[3*i+1],
                 my[ZrC_ZrDIS].s[3*i+2], my[ZrC_ZrDIS].f[3*i] * UENERGY_IN_J /
                 SQUARE(ULENGTH_IN_M) / DIS_DEL / my[ZrC_ZrDIS].a,
                 ZrC_ZrDIS_fx_target[i], 100*rerr
                 (my[ZrC_ZrDIS].f[3*i] * UENERGY_IN_J /
                  SQUARE(ULENGTH_IN_M) / DIS_DEL / my[ZrC_ZrDIS].a,
                  ZrC_ZrDIS_fx_target[i]) );
    fprintf (out, "\n");

    this_structure = &my[ZrC_TEST];
    this_structure->a = my[ZrC_ZrDIS].a;
    for (dE=i=0; i<this_structure->n; i++)
    {
        V3SUB ( &this_structure->s[3*i], &my[ZrC_ZrDIS].s[3*i], DS );
        V3mM3 ( DS, this_structure->H, DX );
        dE -= V3DOT ( &my[ZrC_ZrDIS].f[3*i], DX );
    }
    printf ("%e\n", dE * this_structure->a);
    printf ("%e\n", potential_energy(this_structure, param) -
            potential_energy(&my[ZrC_ZrDIS], param));
    
    return;
} /* end print_property() */


static double xtal_total_energy (int i, double density)
{
    density *= CUBE(LATTICE_CONSTANT/4.711804);
    switch (i)
    {
        case NaCl:
            return (0.419379 * density * density -
                    5.498075 * density - 1534.978880);
        case ZnS:
            return (0.512786 * density * density -
                    5.369512 * density - 1537.816979);
        case CsCl:
            return (0.334708 * density * density -
                    4.585749 * density - 1534.799776);
        default:
            pe("total energy of \"%s\" is unavailable.\n", my[i].name);
    }
    return(0);
} /* end xtal_total_energy() */


#define MIN_POTE           -1553.1
/* #define MAX_POTE           -1549.7 */
#define MAX_POTE           -1549.5
#define MIN_DENSITY         4.
#define MAX_DENSITY         8.
#define DENSITY_MESH        100

void plot_cohesive_energy_curves (char filename[])
{
    int i, j;
    double density, density_del = (MAX_DENSITY-MIN_DENSITY)/DENSITY_MESH;
    double curve_pote_shift;
    TermString buf;
    char mycolor[MAX_STRUCTURE] = {'g', 'b', 'r', 'c', 'm'};
    FILE *out;
    out = fopen (filename, "w+");
    buf[0] = (char)0;
    curve_pote_shift = my[NaCl].pote - xtal_target[NaCl][XTAL_EQ_POTE];
    fprintf (out, "clf;\n\n");
    for (i=0; i<MAX_XTAL_STRUCTURE; i++)
    {
        fprintf (out, "%s = [\n", my[i].name);
        for (j=0; j<=DENSITY_MESH; j++)
        {
            density = MIN_DENSITY + j * density_del;
            my[i].a = pow (my[i].mass / density / my[i].volume, 1./3) * 1E8;
            fprintf (out, "%lf %lf %lf\n", density,
                     xtal_total_energy(i, density) + curve_pote_shift,
                     potential_energy(&my[i], param));
        }
        fprintf (out, "];\n");
        fprintf (out, "plot(%s(:,1), %s(:,2), '%c');\n",
                 my[i].name, my[i].name, mycolor[i]);
        sprintf (buf+strlen(buf), "'%s LDA',", my[i].name);
        if (i==0) fprintf (out, "hold on;\n");
        fprintf (out, "plot(%s(:,1), %s(:,3), '%c--');\n\n",
                 my[i].name, my[i].name, mycolor[i]);
        /* sprintf (buf+strlen(buf), "'%s fit',", my[i].name); */
        sprintf (buf+strlen(buf), "'%s empirical',", my[i].name);
    }
    fprintf (out, "axis([%lf %lf %lf %lf]);\n",
             MIN_DENSITY, MAX_DENSITY,
             MIN_POTE+curve_pote_shift, MAX_POTE+curve_pote_shift);
    buf[strlen(buf)-1] = (char)0;
    fprintf (out, "legend(%s, 3);\n", buf);
    fprintf (out, "xlabel('Crystal density [g/cm^3]');\n");
    fprintf (out, "ylabel('Cohesive energy [eV/pair]');\n");
    fprintf (out, "title('LDA (shifted) {\\it vs} empirical potential "
             " for ZrC crystals');\n");
    fclose (out);
    return;
} /* end plot_cohesive_energy_curves() */


/* save structure expanded in parallelepiped H0[][] (in a) */
/* with respect to x0[] (in a), in readable text format.   */
void write_config (Structure *S, double x0[3], double H0[3][3],
                   double magnification, char filename[])
{
    int i, j, n, *t, maxnc[3], nc[3];
    double *s, s0[3], ds[3], dx[3], H0I[3][3];
    FILE *config;
    M3rows_to_cover_sphere (S->H, M3maxrowradius(H0), maxnc);
    V3aDd (maxnc, 4);
    i = total_replica(maxnc) * S->n;
    MALLOC ( write_config, t,   i, int );
    MALLOC ( write_config, s, 3*i, double );
    V3mM3 (x0, S->HI, s0);
    V3TRIM (s0, s0);
    M3INV (H0, H0I, ds[0]);
    for (i=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]++)
                for (n=0; n<S->n; n++)
                {
                    V3ADDSUB (nc, &S->s[3*n], s0, ds);
                    V3mM3 (ds, S->H, dx);
                    V3mM3 (dx, H0I,  ds);
                    if (V3XOU(ds, -TINY, 1+TINY)) continue;
                    t[i] = S->t[n];
                    /* to fit in the box */
                    for (j=0; j<3; j++) ds[j] = .5 + (ds[j]-.5) * .99999;
                    V3EQV (ds, s+3*i);
                    i++;
                }
    config = wopen(filename);
    fprintf (config, "Number of particles = %d\n\n", i);
    magnification *= S->a;
    fprintf (config, "H(1,1) = %f A\n",   H0[0][0]*magnification);
    fprintf (config, "H(1,2) = %f A\n",   H0[0][1]*magnification);
    fprintf (config, "H(1,3) = %f A\n\n", H0[0][2]*magnification);
    fprintf (config, "H(2,1) = %f A\n",   H0[1][0]*magnification);
    fprintf (config, "H(2,2) = %f A\n",   H0[1][1]*magnification);
    fprintf (config, "H(2,3) = %f A\n\n", H0[1][2]*magnification);
    fprintf (config, "H(3,1) = %f A\n",   H0[2][0]*magnification);
    fprintf (config, "H(3,2) = %f A\n",   H0[2][1]*magnification);
    fprintf (config, "H(3,3) = %f A\n\n", H0[2][2]*magnification);
    fprintf (config, "TP Mass(amu)     sx        sy        sz ");
    fprintf (config, "    sx1(1/ns)  sy1(1/ns)  sz1(1/ns)\n");
    for (j=0; j<3*i; j+=3)
        fprintf (config, "%2s %.6f  %.6f  %.6f  %.6f  %.6f  %.6f  %.6f\n",
                 ATOM_SYMBOL(Z[t[j/3]]), ATOM_MASS_IN_AMU(Z[t[j/3]]),
                 s[j], s[j+1], s[j+2], 0., 0., 0.);
    fclose (config);
    free(s);
    free(t);
    return;
}  /* end write_config() */