/* Calculate unrelaxed elastic constants of configuration: */
/* np,H,s,tp,N, f(scratched),stress(scratched) must exist, */
/* and also the corresponding volume and pote.             */

double Calculate_Unrelaxed_Bulk_Modulus (double delta)
{
    M3 HH;
    double pote_P, JP[3][3], pote_M, JM[3][3];
    M3diagonal (1+delta/3, 1+delta/3, 1+delta/3, JP);
    M3diagonal (1-delta/3, 1-delta/3, 1-delta/3, JM);
    M3MUL (H, JP, HH);
    pote_P = potential(np,HH,s,tp,N, f,stress);
    M3MUL (H, JM, HH);
    pote_M = potential(np,HH,s,tp,N, f,stress);
    return ((pote_P + pote_M - 2 * pote ) / delta / delta / M3VOLUME(H));
} /* end Calculate_Unrelaxed_Bulk_Modulus() */


double Calculate_Unrelaxed_C11_C12 (double delta)
{
    M3 HH;
    double pote_P, JP[3][3], pote_M, JM[3][3];
    M3diagonal (1+delta, 1+delta, 1/(1+delta)/(1+delta), JP);
    M3diagonal (1-delta, 1-delta, 1/(1-delta)/(1-delta), JM);
    M3MUL (H, JP, HH);
    pote_P = potential(np,HH,s,tp,N,f,stress);
    M3MUL (H, JM, HH);
    pote_M = potential(np,HH,s,tp,N,f,stress);
    return ((pote_P + pote_M - 2 * pote) / 6 / delta / delta / M3VOLUME(H));
} /* end Calculate_Unrelaxed_C11_C12() */


double Calculate_Unrelaxed_C44 (double delta)
{
    M3 HH;
    double pote_P, JP[3][3], pote_M, JM[3][3];
    M3ASSIGN(1,delta,delta,delta,1,delta,delta,delta,1,JP);
    M3ASSIGN(1,-delta,-delta,-delta,1,-delta,-delta,-delta,1,JM);
    M3MUL (H, JP, HH);
    pote_P = potential(np,HH,s,tp,N,f,stress);
    M3MUL (H, JM, HH);
    pote_M = potential(np,HH,s,tp,N,f,stress);
    return ((pote_P + pote_M - 2 * pote) / 12 / delta / delta / M3VOLUME(H));
} /* end Calculate_Unrelaxed_C44() */