/*********************************************/
/* Force constant calculation for Zr/C cell. */
/*                                           */
/*       Jun. 1, 2000 Ju Li <liju99@mit.edu> */
/*********************************************/

#include "../relax.h"
#include <Phonon.h>

#include "../potential.c"
#include "../elastic_constants.c"

#define DEF_FNAME  "fconst"

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

static double pote_before_relaxation;
static double atom_force_norm;
static double ep_reference;
static double volumep_reference;
static double lattice_constant_reference;
static double B_reference, C11_C12_reference, C44_reference,
    C11_reference, C12_reference;

double hydrostatic (double ratio)
{
    double value;
    M3 HH;
    M3MULTIPLY (ratio, H, HH);
    value = potential(np,HH,s,tp,N,f,stress);
    value = M3TR(stress);
    value *= value;
    return(value);
} /* end hydrostatic() */

/* objective function of CG minimization */
static double free_energy (double *s_new, double *f)
{
    register int i;
    V3 fs;
    M3 HI;
    s = s_new;
    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]);
    }
    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);
    return (pote);
} /* end free_energy() */


#define REFERENCE_NC   8
#define REFERENCE_NPA  2

int main (int argc, char *argv[])
{
    int i, j, k, l, m, n, w, *idx;
    double cc, ratio, delta=DELTA, *old_f, *force, *magnitude;
    TermString buf, fname;
    M3 eta, HI;
    V3 old_s, ds;
    Phonon_force_constant P[1]={0};

    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");

    for (i=0; i<3; i++)
        for (j=i; j<3; j++)
        {
            fscanf (stdin, "%s = %lf\n", buf, &eta[i][j]);
            eta[j][i] = eta[i][j];
        }
    fscanf (stdin, "\n");

    fscanf (stdin, "Force constant file (default=\"fconst\"):\n%s\n\n",
            fname);
    if (!strcasecmp(fname, "default")) strcpy (fname, DEF_FNAME);

    Fprintf (stdout, "**** Building crystal reference under %s potential\n",
             potential_name);
    Config_rebuild_Xtal
        (Config_Aapp_to_Alib, REFERENCE_NC,REFERENCE_NC,REFERENCE_NC,
         Xtal_abstracts+XTAL_nacl, "Zr", -1., "C", -1., 4.7);
    rebind_ct (Config_Aapp_to_Alib, "Zr C", ct, &tp, stdout); cr();
    Neighborlist_Recreate_Form (Config_Aapp_to_Alib, ct, N);
    N->min_shrinkage = 0.9;
    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;
    Neighborlist_Recreate (Config_Aapp_to_Alib, stdout, ct, &tp, N); cr();
    REALLOC (main, f, DIMENSION*np, double);
    onemin (hydrostatic, N->min_shrinkage, 2-N->min_shrinkage,
            SQUARE(EPS), &ratio);
    M3MultiplY (ratio, H);
    pote = potential(np,H,s,tp,N,f,stress);
    ep_reference = pote / np;
    if (NOTSMALL(M3TR(stress)))
    {
        M3MultiplY (USTRESS_IN_MPA, stress);
        S3fPR (stdout, "stress = %M MPa.\n", stress);
        fprintf (stdout, "trace = %e MPa.\n", M3TR(stress));
        pe ("main: stress at minimum is nonzero.\n");
    }
    volume = M3VOLUME(H);
    volumep_reference = volume / np;
    lattice_constant_reference = cbrt(8 * volumep_reference);
    B_reference = Calculate_Unrelaxed_Bulk_Modulus (delta);
    C11_C12_reference = Calculate_Unrelaxed_C11_C12 (delta);
    C44_reference = Calculate_Unrelaxed_C44 (delta);
    C11_reference = (3*B_reference + 2*C11_C12_reference) / 3;
    C12_reference = C11_reference - C11_C12_reference;
    fprintf (stdout, "reference lattice constant = %f = %f A,\n",
             lattice_constant_reference,
             lattice_constant_reference * ULENGTH_IN_A);
    fprintf (stdout, "reference volume per particle = %f = %f A^3,\n",
             volumep_reference,
             volumep_reference * UVOLUME_IN_A3);
    fprintf (stdout, "reference energy per particle = %f = %f eV,\n",
             ep_reference, ep_reference * UENERGY_IN_EV);
    fprintf (stdout, "C11 = %f = %f MPa;\n",
             C11_reference, C11_reference * USTRESS_IN_MPA);
    fprintf (stdout, "C12 = %f = %f MPa;\n",
             C12_reference, C12_reference * USTRESS_IN_MPA);
    fprintf (stdout, "C44 = %f = %f MPa;\n",
             C44_reference, C44_reference * USTRESS_IN_MPA);
    fprintf (stdout, "B = %f = %f MPa;\n",
             B_reference, B_reference * USTRESS_IN_MPA);
    cr();

    if (M3NEZERO(eta))
    {
        pure_deform (H, eta, H);
        S3fPR (stdout, "We impose Lagrangian strain = %M.\n ", eta);
        S3fPR (stdout, "new H = %M A.\n ", H);
        N->min_shrinkage = 1;
        N->max_atom_displacement = DEF_MAX_ATOM_DISPLACEMENT;
        Neighborlist_Recreate (Config_Aapp_to_Alib, stdout, ct, &tp, N); cr();
        Fprintf (stdout, "********************* "
                 "Before Minimization **********************\n");
        pote_before_relaxation = potential(np,H,s,tp,N,f,stress);
        Fprintf (stdout, "The total energy = %f eV, energy/atom = %f eV,\n",
                 pote_before_relaxation*UENERGY_IN_EV,
                 pote_before_relaxation/np*UENERGY_IN_EV);
        Fprintf (stdout, "so the defect formation energy (unrelaxed)\n");
        Fprintf (stdout, "= %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 (stdout,
                 "(energy reference is perfect crystal at zero stress)\n");
        Fprintf (stdout, "**********************"
                 "******************************************\n\n");
        zxcgr (free_energy, DIMENSION*np, s); cr();
        Fprintf (stdout, "********************* "
                 "After Minimization ***********************\n");
        Fprintf (stdout, "The total energy = %f eV, energy/atom = %f eV,\n",
                 pote*UENERGY_IN_EV, pote*UENERGY_IN_EV/np);
        Fprintf (stdout, "so the defect formation energy (relaxed)\n");
        Fprintf (stdout, "= %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 (stdout,
                 "and it has increased by %f eV during the relaxation.\n",
                 (pote-pote_before_relaxation)*UENERGY_IN_EV);
        Fprintf (stdout, "atom residual |f| = %e eV/A,\n",
                 atom_force_norm*UENERGY_IN_EV/ULENGTH_IN_A);
        M3MULTIPLY (-USTRESS_IN_MPA, stress, eta);
        S3fPR (stdout, "stress = %M MPa.\n", eta);
        Fprintf (stdout, "**********************"
                 "******************************************\n\n");
    }

    PHONON_force_constant_CALC
        (Config_Aapp_to_Alib, ct, tp, N, REFERENCE_NC, &potential, P);
    P->fconst_fn_basename = fname;
    phonon_force_constant_SAVE (P); cr();

    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() */