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

#include "z.h"

double Bulk_Modulus (Structure *S)
{
    Structure D;
    double pote_P, JP[3][3]={{1+DELTA/3,0,0},
                             {0,1+DELTA/3,0},
                             {0,0,1+DELTA/3}};
    double pote_M, JM[3][3]={{1-DELTA/3,0,0},
                             {0,1-DELTA/3,0},
                             {0,0,1-DELTA/3}};
    D = *S;
    M3MUL (S->H, JP, D.H);
    pote_P = potential_energy (&D, param);
    M3MUL (S->H, JM, D.H);
    pote_M = potential_energy (&D, param);
    return ((pote_P + pote_M - 2 * S->pote) / DELTA / DELTA /
            (S->volume * CUBE(S->a)) * USTRESS_IN_GPA);
} /* end Bulk_Modulus() */


double C11_C12 (Structure *S)
{
    Structure D;
    double pote_P, JP[3][3]={{1+DELTA,0,0},
                             {0,1+DELTA,0},
                             {0,0,1/(1+DELTA)/(1+DELTA)}};
    double pote_M, JM[3][3]={{1-DELTA,0,0},
                             {0,1-DELTA,0},
                             {0,0,1/(1-DELTA)/(1-DELTA)}};
    D = *S;
    M3MUL (S->H, JP, D.H);
    pote_P = potential_energy (&D, param);
    M3MUL (S->H, JM, D.H);
    pote_M = potential_energy (&D, param);
    return ((pote_P + pote_M - 2 * S->pote) / 6 / DELTA / DELTA /
            (S->volume * CUBE(S->a)) * USTRESS_IN_GPA);
} /* end C11_C12() */


double C44_0 (Structure *S)
{
    Structure D;
    double pote_P, JP[3][3]={{1,DELTA,DELTA},
                             {DELTA,1,DELTA},
                             {DELTA,DELTA,1}};
    double pote_M, JM[3][3]={{1,-DELTA,-DELTA},
                             {-DELTA,1,-DELTA},
                             {-DELTA,-DELTA,1}};
    D = *S;
    M3MUL (S->H, JP, D.H);
    pote_P = potential_energy (&D, param);
    M3MUL (S->H, JM, D.H);
    pote_M = potential_energy (&D, param);
    return ((pote_P + pote_M - 2 * S->pote) / 12 / DELTA / DELTA /
            (S->volume * CUBE(S->a)) * USTRESS_IN_GPA);
} /* end C44_0() */