double ep_reference;
double volumep_reference;
double lattice_constant_reference;
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() */

#define REFERENCE_NC  5
void rebuild_reference (FILE *info)
{
    double ratio;
    double delta = pow(EPS,1./4.);
    Fprintf (info, "**** 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.711804);
    rebind_ct (Config_Aapp_to_Alib, "Zr C", ct, &tp, info); Fcr(info);
    Neighborlist_Recreate_Form (Config_Aapp_to_Alib, ct, N);
    N->min_shrinkage = 0.8;
    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, info, ct, &tp, N); Fcr(info);
    REALLOC (rebuild_reference, 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 (info, "stress = %M MPa.\n", stress);
        fprintf (info, "trace = %e MPa.\n", M3TR(stress));
        pe ("rebuild_reference: stress at minimum is nonzero.\n");
    }
    volumep_reference = M3VOLUME(H) / np;
    lattice_constant_reference = cbrt(M3VOLUME(H)) / REFERENCE_NC;
    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;
    if (info)
    {
        fprintf (info, "reference lattice constant = %f = %f A,\n",
                 lattice_constant_reference,
                 lattice_constant_reference * ULENGTH_IN_A);
        fprintf (info, "reference volume per particle = %f = %f A^3,\n",
                 volumep_reference,
                 volumep_reference * UVOLUME_IN_A3);
        fprintf (info, "reference energy per particle = %f = %f eV,\n",
                 ep_reference, ep_reference * UENERGY_IN_EV);
        fprintf (info, "C11 = %f = %f MPa;\n",
                 C11_reference, C11_reference * USTRESS_IN_MPA);
        fprintf (info, "C12 = %f = %f MPa;\n",
                 C12_reference, C12_reference * USTRESS_IN_MPA);
        fprintf (info, "C44 = %f = %f MPa;\n",
                 C44_reference, C44_reference * USTRESS_IN_MPA);
        fprintf (info, "B = %f = %f MPa;\n",
                 B_reference, B_reference * USTRESS_IN_MPA);
    }
    return;
} /* end rebuild_reference() */

#ifdef _reference_TEST
Aapp_Define_Config;
Chemtab ct[1]={0};
Tp *tp=NULL;
Neighborlist N[1]={0};
double pote, volume, *f=NULL;
M3 stress;
int main (int argc, char *argv[])
{
    rebuild_reference (stdout);
    return (0);
}
#endif /* _reference_TEST */