#include "MD.h"
#include "common.c"
#include "Relax/potential.c"
#include "Relax/elastic_constants.c"
#include "Relax/reference.c"
#include "Relax/work.c"
#include "auxiliary.c"
#include "interatomic.c"
#include "engine.c"

#ifndef _TEST
int main (int argc, char *argv[])
{
    int i,j;
    double cc,dd,ee;
    TermString buf;
    M3 eta;

    for (buf[0]=(char)0,i=0; i<argc; i++)
        strcat(strcat(buf,argv[i]),(i==argc-1)?"":" ");
    printf ("Program \"%s\" starts (pid = %d).\n",
            buf, getpid());

    fscanf (stdin, "SAVE SCREEN OUTPUT TO FILE (DEFAULT=md, "
            "NO=DON'T):\n%s\n", buf);
    if (strcasecmp(buf, "NO"))
    {
        if (!strcasecmp(buf, "default")) strcpy (buf, "md");
        fp_screen = wopen (buf);
        ftie (stdout, fp_screen);
    }
    else ft = stdout;
    
    fscanf (stdin, "POTENTIAL FUNCTION TO BE USED (NOC):\n%s\n",
            potential_name);
    convert_to_uppercase (potential_name);
    NOC = !strcasecmp(potential_name, "NOC");
    if (!NOC) pe ("Please choose NOC potential.\n");
    Fprintf (ft, "Potential = %s, RCUT_ZrZr = %g [A], RCUT_ZrC = %g [A]\n",
             potential_name, RCUT_ZrZr*ULENGTH_IN_A, RCUT_ZrC*ULENGTH_IN_A);

    fscanf (stdin, "INTEGRATION ALGORITHM:  g (TIME/EVALUATION) [fs]:  "
            "halt [g]:\n%s %lf %d\n", buf, &g, &halt_in_g);
    convert_to_uppercase (buf);
    for (integrator_id=0; integrator_id<NEWTONIAN_MAX; integrator_id++)
        if (!strcasecmp(buf,newtonian[integrator_id].name)) break;
    if (integrator_id == NEWTONIAN_MAX)
        pe ("Integration scheme \"%s\" does not exist.\n", buf);
    Fprintf (ft, "Integrator = %s,  K = %d.\n",
             buf, newtonian[integrator_id].K);
    g /= UTIME_IN_FS;
    Fprintf (ft, "g = %g = %g [fs],  full step = %g = %g [fs].\n",
             g, g*UTIME_IN_FS, newtonian[integrator_id].K*g,
             newtonian[integrator_id].K*g*UTIME_IN_FS);
    if ( halt_in_g % newtonian[integrator_id].K != 0 )
    {
        Fprintf (ft, "** halt [g] = %d is not integer multiple"
                 " of K = %d, truncated ->\n", halt_in_g,
                 newtonian[integrator_id].K);
        halt_in_g = (halt_in_g / newtonian[integrator_id].K + 1) *
            newtonian[integrator_id].K;
    }
    Fprintf (ft, "halt = %d [g] = %d full steps.\n",
             halt_in_g, halt_in_g / newtonian[integrator_id].K);

    fscanf (stdin, "MD TIME [g]:  REPORT FREQ [g]:  "
            "FINAL AVERAGE STARTS AT [g]:\n%d %d %d\n",
            &md_in_g, &report_freq_in_g, &final_avg_start_in_g);
    md_in_g = md_in_g / halt_in_g * halt_in_g;
    Fprintf (ft, "MD runs %d [g] = %d [halt] = %g = %g [ps].\n",
             md_in_g, md_in_g / halt_in_g, md_in_g*g,
             md_in_g*g*UTIME_IN_PS);
    report_freq_in_g = report_freq_in_g / halt_in_g * halt_in_g;
    Fprintf (ft, "report freq = %d [g] = %d [halt] = %g = %g [ps].\n",
             report_freq_in_g, report_freq_in_g / halt_in_g,
             report_freq_in_g*g, report_freq_in_g*g*UTIME_IN_PS);
    Fprintf (ft, "MD evolution -> \"%s\".\n", fn_draw);
    fp_draw = wopen (fn_draw);
    final_avg_start_in_g = final_avg_start_in_g / halt_in_g * halt_in_g;
    Fprintf (ft, "Final average begins at %d [g] = %d [halt] = "
             "%g = %g [ps].\n", final_avg_start_in_g,
             final_avg_start_in_g / halt_in_g, final_avg_start_in_g*g,
             final_avg_start_in_g*g*UTIME_IN_PS);
    if (final_avg_start_in_g >= md_in_g)
        pe ("final_avg_start = %d [g]  >=  MD TIME = %d [g],\n"
            "there is no room left for property averaging.\n",
            final_avg_start_in_g, md_in_g);

    fscanf (stdin, "SAVE FINAL CONFIGURATION ON FILE (default=md.cfg):\n");
    fscanf (stdin, "%s\n", fn_final_config);
    if (!strcasecmp(fn_final_config, "default"))
        strcpy (fn_final_config, "md.cfg");
    Fprintf (ft, "The final configuration is saved on \"%s\".\n",
             fn_final_config);
    fscanf (stdin, "SAVE TEMPORARY CONFIGURATION "
            "(.*=NON-OVERWRITING,NO=DON'T):  WITH FREQ [g]:\n%s %s\n",
            fn_tmp_config, buf);
    if (!strcasecmp(fn_tmp_config, "default"))
        strcpy (fn_tmp_config, fn_final_config);
    Fprintf (ft, "And temporary configuration -> \"%s\",\n",
             fn_tmp_config);
    tmp_config_serial_mode = strstr(fn_tmp_config, ".*");
    if (tmp_config_serial_mode)
    {
        *(++tmp_config_serial_mode) = EOS;
        ++tmp_config_serial_mode;
        Fprintf (ft, "with serial mode turned on, and saved\n");
    }
    else Fprintf (ft, "with serial mode turned off, and saved\n");
    tmp_config_freq_in_g = (!strcasecmp(buf, "default")) ?
        DEF_TMP_CONFIG_FREQ * report_freq_in_g :
        atoi(buf) / halt_in_g * halt_in_g;
    if (tmp_config_freq_in_g == 0) tmp_config_freq_in_g = halt_in_g;
    Fprintf (ft, "every %d [g] = %d [halt] = %g = %g [ps].\n",
             tmp_config_freq_in_g, tmp_config_freq_in_g / halt_in_g,
             tmp_config_freq_in_g*g, tmp_config_freq_in_g*g*UTIME_IN_PS);
    Fcr (ft);
    fscanf (stdin, "\n");

    if ( tmp_config_serial_mode && (md_in_g / tmp_config_freq_in_g >= 64) )
    {
        Fprintf (ft, "** Warning: %d copies of %s*%s will be saved,\n"
                 "PERHAPS TOO MUCH?\n\n", md_in_g / tmp_config_freq_in_g,
                 fn_tmp_config, tmp_config_serial_mode);
        waitfor (SUFFICIENT_TIME_TO_NOTICE_IN_SECONDS);
    }

    rebuild_reference (ft); Fcr(ft);

    fscanf (stdin, "LOAD CONFIGURATION FROM FILE (NO->FRESHSTART):\n");
    fscanf (stdin, "%s\n", fn_input_config);
    fscanf (stdin, "IF FRESHSTART, THEN NC1:  NC2:  NC3:  \n"
            "LATTICE CONSTANT [A] (default->P=0):\n%d %d %d %s\n",
            &nc[0], &nc[1], &nc[2], buf);
    cc = (!strcasecmp(buf, "default")) ?
        lattice_constant_reference * ULENGTH_IN_A : atof(buf);
    if (!strcasecmp(fn_input_config, "NO"))
        Config_rebuild_Xtal
            (Config_Aapp_to_Alib, nc[0],nc[1],nc[2],
             Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1., cc);
    else Config_Load (fn_input_config, ft, Config_Aapp_to_Alib);
    Config_fold_into_PBC (Config_Aapp_to_Alib);
    M3EQV( H, H0 );
    rebind_ct (Config_Aapp_to_Alib, "Zr C", ct, &tp, ft); Fcr(ft);
    Neighborlist_Recreate_Form (Config_Aapp_to_Alib, ct, N);
    fscanf (stdin, "NTH or NTS mode:\n%s\n", buf);
    NTH = !strcasecmp(buf, "NTH");
    NTS = !strcasecmp(buf, "NTS");
    if ( (!NTH) && (!NTS) )
        pe ("main: must select between NTH or NTS mode.\n");
    fscanf (stdin, "IF NTH, THEN APPLY STRAIN SCHEDULE "
            "(11,12,13,22,23,33):\n");
    if (NTH)
    {
        Fprintf (ft, "Starting Mode = NTH.\n");
        if ( grab_schedules
             (stdin, 0, md_in_g, DIMSYM, eta_schedule_tokens,
              eta_schedules, ft) == SCHEDULES_ALL_CONSTANTS )
        {
            Fcr (ft);
            get_scheduled_symmat (eta_schedules, 0, eta);
            pure_deform (H0, eta, H);
            S3fPR (ft, "We impose Lagrangian strain = %M.\n ", eta);
            S3fPRmul (ft, "new H = %M [A],\nwhich remains fixed.",
                      H, ULENGTH_IN_A);
            N->min_shrinkage = 1;
        }
        else N->min_shrinkage = DEF_MIN_SHRINKAGE;
        Fcr (ft);
        freedom = DIMENSION * np;
        REALLOC (main, f, freedom, double);
    }
    fscanf_skipline (stdin);
    fscanf (stdin, "IF NTS, THEN APPLY STRESS [MPa] SCHEDULE:  "
            "VOLUMETRIC RESPONSE CONSTRAINT?\n");
    if (NTS)
    {
        Fprintf (ft, "Starting Mode = NTS.\n");
        grab_schedules
            (stdin, 0, md_in_g, DIMSYM, sout_schedule_tokens,
             sout_schedules, ft);
        Fcr (ft);
        for (i=0; i<DIMSYM; i++)
            schedulelist_scale_value (sout_schedules[i], MPA_IN_USTRESS);
        N->min_shrinkage = DEF_MIN_SHRINKAGE;
        fscanf (stdin, "%s", buf);
        VOLUMETRIC =
            (!strcasecmp(buf, "VOLUMETRIC")) ||
            (!strcasecmp(buf, "yes"));
        Fprintf (ft, "Volumetric constraint = %s.\n", BWORD(VOLUMETRIC));
        Fcr (ft);
        i = VOLUMETRIC ? 1 : DIMSYM;
        freedom = DIMENSION * np + i;
        REALLOC (main, s,  freedom, double);
        VZERO ( i, &(s[DIMENSION*np]) );
        REALLOC (main, s1, freedom, double);
        VZERO ( i, &(s1[DIMENSION*np]) );
        REALLOC (main, f, freedom, double);
    }
    fscanf_skipline (stdin);
    fscanf (stdin, "T [K] SCHEDULE:  T/STRESS COUPLING TIME CONSTANT [g] "
            "(0->CHECK CONSERVATION):\n");
    grab_schedules (stdin, 0, md_in_g, 1, &T_schedule_token, &T_schedule, ft);
    Fcr (ft);
    fscanf (stdin, "%s\n", buf);
    coupling_time_constant = (!strcasecmp(buf, "default")) ?
        DEF_COUPLING_TIME_CONSTANT*g : atof(buf)*g;
    if (coupling_time_constant > 0)
        Fprintf (ft, "T/stress coupling time constant = %g [g] = %g "
                 "= %g [ps].\n", coupling_time_constant / g,
                 coupling_time_constant, coupling_time_constant*UTIME_IN_PS);
    else
    {
        Fprintf (ft, "Disengage T coupling%s to\n"
                 "verify conservative Hamiltonian dynamics.\n",
                 NTH ? "" : " and wall kinetic energy protection");
        waitfor (SUFFICIENT_TIME_TO_NOTICE_IN_SECONDS);
    }
    if (NTS)
    {
        if (coupling_time_constant > 0)
            wall = B_reference * M3VOLUME(H0) *
                SQUARE(coupling_time_constant);
        else wall = B_reference * M3VOLUME(H0) *
                 SQUARE(DEF_COUPLING_TIME_CONSTANT*g);
        Fprintf (ft, "NTS wall constant W = %e.\n", wall);
        /* sout,T conditions reach stasis before final_avg_start? */
        preliminary.NTS_calm_in_g =
            schedules_calm (final_avg_start_in_g, md_in_g,
                            DIMSYM, sout_schedules, 1, &T_schedule, 0);
        preliminary.NTS_calm_and_start_average_in_g =
            rint( (1-NTS_CALM_AND_START_AVERAGE) * preliminary.NTS_calm_in_g
                  +  NTS_CALM_AND_START_AVERAGE  * final_avg_start_in_g )
            / halt_in_g * halt_in_g;
        preliminary.NTS_calm_and_stop_average_in_g =
            rint( (1-NTS_CALM_AND_STOP_AVERAGE) * preliminary.NTS_calm_in_g
                  +  NTS_CALM_AND_STOP_AVERAGE  * final_avg_start_in_g )
            / halt_in_g * halt_in_g;
        preliminary.do_preliminary =
            ( final_avg_start_in_g >=
              preliminary.NTS_calm_and_stop_average_in_g ) &&
            ( preliminary.NTS_calm_and_stop_average_in_g >
              preliminary.NTS_calm_and_start_average_in_g );
        if ( preliminary.do_preliminary )
        {
            Fprintf (ft, "\nYou have chosen to de-activate P.R. at %d [g]\n"
                     "with the H averaged between %d and %d [g], and\n"
                     "finally to decouple T at %d g, right?\n",
                     preliminary.NTS_calm_and_stop_average_in_g,
                     preliminary.NTS_calm_and_start_average_in_g,
                     preliminary.NTS_calm_and_stop_average_in_g,
                     final_avg_start_in_g);
            waitfor (SUFFICIENT_TIME_TO_NOTICE_IN_SECONDS);
        }
    }
    Fcr (ft);
    fscanf (stdin, "\n");

    fscanf (stdin, "REQUEST CALCULATING MSD:  G(R):  "
            "THERMAL CONDUCTIVITY:  SHEAR VISCOSITY:\n");
    fscanf (stdin, "%s", buf);
    N->track_dx = !strcasecmp(buf, "yes");
    fscanf (stdin, "%s", buf);
    CALC_GR = !strcasecmp(buf, "yes");
    Fprintf (ft, "MSD calculation = %s,  g(r) calculation = %s.\n",
             BWORD(N->track_dx), BWORD(CALC_GR));
    if (CALC_GR)
    {
        Gr_Recreate_Form (Config_Aapp_to_Alib, ct, GR);
        i = Gr_Reset (Config_Aapp_to_Alib, ct, GR, ft); Fcr(ft);
        GR_freq_in_g = ( (md_in_g - final_avg_start_in_g)
                         / i / halt_in_g + 1 ) * halt_in_g;
        Fprintf (ft, "Accumulate g(r) instance every %d [g] = %d [halt] = %g\n"
                 "= %g [ps] after final average starts at %d [g] = %g [ps].\n",
                 GR_freq_in_g, GR_freq_in_g / halt_in_g,
                 GR_freq_in_g*g, GR_freq_in_g*g*UTIME_IN_PS,
                 final_avg_start_in_g, final_avg_start_in_g*g*UTIME_IN_PS);
        Fcr (ft);
    }
    fscanf (stdin, "%s", buf);
    CALC_CONDUCTIVITY = !strcasecmp(buf, "yes");
    if ( CALC_CONDUCTIVITY )
    {
        REALLOC (main, v,  DIMENSION * np, double);
        REALLOC (main, ep, DIMENSION * np, double);
        spectral_openflux (halt_in_g*g*UTIME_IN_S, SPECTRAL_TYPE_HEAT);
    }
    fscanf (stdin, "%s\n", buf);
    CALC_VISCOSITY = !strcasecmp(buf, "yes");
    if ( CALC_VISCOSITY )
        spectral_openflux (halt_in_g*g*UTIME_IN_S, SPECTRAL_TYPE_STRESS);
    Fprintf (ft, "Thermal conductivity = %s,  shear viscosity = %s.\n",
             BWORD(CALC_CONDUCTIVITY), BWORD(CALC_VISCOSITY));
    Fcr (ft);

    fscanf (stdin, "RANDOM NUMBER SEED (0->TIME SEED, "
            "default->FIXED SEED):\n%s\n", buf);
    i = (!strcasecmp(buf, "default")) ? DEFAULT_RANDOM_SEED : atoi(buf);
    if (i == 0) i = time(NULL);
    Randomize (i);
    Fprintf (ft, "Random number seed = %d.\n", i);
    Fcr(ft);

    Fprintf (ft, "Instructions read complete.\n");
    Fcr(ft);

    N->max_atom_displacement = DEF_MAX_ATOM_DISPLACEMENT;
    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, ft, ct, &tp, N);
    Fcr(ft);
    
    if (!strcasecmp(fn_input_config, "NO"))
    {
        Fprintf (ft, "As this is a freshstart, we assign thermal\n"
                 "velocities with %g x kinetic energy overdrive.\n",
                 FRESHSTART_OVERDRIVE);
        Fcr(ft);
        cc = FRESHSTART_OVERDRIVE * schedulelist_value(T_schedule,0);
        motion = Config_motion_initialize (cc, Config_Aapp_to_Alib);
    }
    else
    {
        Fprintf (ft, "Configuration \"%s\" has the following motion "
                 "properties:\n", fn_input_config);
        motion = Config_analyze_motion (Config_Aapp_to_Alib, ft);
        Fcr(ft);
    }

    actual_T.md = motion->T_in_K;
    actual_T.real = Tmd_to_Treal (actual_T.md, &actual_T.dTmd__dTreal);
    if (actual_T.real < JUMPSTART_T_IN_K)
    {
        desired_T.real = JUMPSTART_T_IN_K;
        desired_T.md = Treal_to_Tmd (desired_T.real, &desired_T.dTmd__dTreal);
        Fprintf (ft, "T = %g [K] is too small to start, Sir.\n"
                 "We give it a starting temperature of %g [K].\n",
                 actual_T.real, desired_T.real);
        Fcr(ft);
        motion = Config_motion_initialize (desired_T.md, Config_Aapp_to_Alib);
        waitfor (SUFFICIENT_TIME_TO_NOTICE_IN_SECONDS);
    }

    if ( V3NEZERO( motion->s1_avg ) )
    {
        motion = Config_motion_zero_drift (motion, Config_Aapp_to_Alib);
        Fprintf (ft, "Captain Picard: average velocity zapped!\n");
        Fcr(ft);
    }

    t_in_g = -1;
    potential_realloc(np);
    newtonian_integrator_prep
        (integrator_id, g, freedom, s, s1, f, &engine);

    taskwatch ("MD", TASKWATCH_INIT);
    Fprintf (ft, "MD simulation starts on %s.\n", date_string());
    Fcr(ft);

    t_in_g = 0;
    for (i=0; i < md_in_g / halt_in_g; i++)
    {
        if ( NTS && preliminary.do_preliminary )
        {
            if ( t_in_g < preliminary.NTS_calm_and_stop_average_in_g )
            {
                preliminary.preliminary_done = FALSE;
                if ( t_in_g >= preliminary.NTS_calm_and_start_average_in_g )
                    avg_add ( Avg_what(Geometry, preliminary.geometry),
                              M3e(H) );
            }
            else if ( ! preliminary.preliminary_done )
            { /* one shot */
                avg_done ( Avg_what(Geometry, preliminary.geometry) );
                M3EQV ( preliminary.geometry.H, H );
                freedom = DIMENSION * np;
                preliminary.preliminary_done = TRUE;
            }
        }
        newtonian_integrator_firststep
            (integrator_id, g, freedom, s, s1, f, &engine);
        for (j=1; j<halt_in_g/newtonian[integrator_id].K; j++)
            newtonian_integrator_nextstep
                (integrator_id, g, freedom, s, s1, f, &engine);
        newtonian_integrator_halt
            (integrator_id, g, freedom, s, s1, f, &engine);
        if (t_in_g % tmp_config_freq_in_g == 0)
            if (tmp_config_serial_mode)
            {
                sprintf (buf, "%s%d%s", fn_tmp_config, t_in_g,
                         tmp_config_serial_mode);
                save_config (buf);
            }
            else save_config (fn_tmp_config);
        if ( N->track_dx && (t_in_g >= final_avg_start_in_g) &&
             (!final_avg_start.snapshot_taken) )
        {
            final_avg_start.t_in_g = t_in_g;
            final_avg_start.MSD = MSD();
            final_avg_start.snapshot_taken = TRUE;
        }
        if ( CALC_GR && (t_in_g >= final_avg_start_in_g) &&
             (t_in_g % GR_freq_in_g == 0) )
            Gr_Accumulate (Config_Aapp_to_Alib, ct, tp, GR);
    }
    fclose (fp_draw);
    save_config (fn_final_config);

    avg_done ( Avg_what(Property, final) );
    Fprintf (ft, "The final average was from %d to %d [g] = %g [ps]:\n",
             final_avg_start_in_g, md_in_g,
             (md_in_g - final_avg_start_in_g) * g * UTIME_IN_PS );
    S3fPRmul (ft, "H = %M [A];\n", final.H, ULENGTH_IN_A);
    Fprintf (ft, "volume = %g [A^3] = %g [A^3/atom];\n",
             M3VOLUME(final.H) * UVOLUME_IN_A3,
             M3VOLUME(final.H) * UVOLUME_IN_A3 / np);
    Fprintf (ft, "density = %g [mol/m^3] = %g [g/cm^3];\n",
             np/M3VOLUME(final.H)*UNUMBERDENSITY_IN_MOL__M3,
             motion->mobilemass/M3VOLUME(final.H)*UMASSDENSITY_IN_G__CM3);
    Fprintf (ft, "kine = %.8g [eV] = %.8g [eV/atom] = %g [kJ/mol];\n",
             final.kine*UENERGY_IN_EV, final.kine/np*UENERGY_IN_EV,
             final.kine/np*UENERGY_IN_KJ__MOL);
    cc = 2 * final.kine / DIMENSION / (motion->nmobile-1) / BOLZ_IN_UENERGY__K;
    dd = Tmd_to_Treal(cc, &ee);
    Fprintf (ft, "Tmd = %g [K], Treal = %g [K], dTmd__dTreal = %g;\n",
             cc, dd, ee);
    Fprintf (ft, "pote = %.8g [eV] = %.8g [eV/atom] = %g [kJ/mol];\n",
             final.pote*UENERGY_IN_EV, final.pote/np*UENERGY_IN_EV,
             final.pote/np*UENERGY_IN_KJ__MOL);
    cc = final.pote + final.kine;
    Fprintf (ft, "tote = %.8g [eV] = %.8g [eV/atom] = %g [kJ/mol];\n",
             cc*UENERGY_IN_EV, cc/np*UENERGY_IN_EV, cc/np*UENERGY_IN_KJ__MOL);
    S3fPRmul (ft, "internal stress = %M [MPa];\n",
              final.stress, USTRESS_IN_MPA);
    M3TRACELESS ( final.stress, eta, cc );
    dd = M32NORM (eta) / SQRT2;
    cc /= 3;
    Fprintf (ft, "pressure = %g [MPa], von Mises shear stress "
             "= %g [MPa];\n", cc*USTRESS_IN_MPA, dd*USTRESS_IN_MPA);
    if (cc != 0)
    {
        dd = final.pote + final.kine + cc * M3VOLUME(final.H);
        Fprintf (ft, "enthalpy = %.8g [eV] = %.8g [eV/atom] = "
                 "%g [kJ/mol];\n", dd*UENERGY_IN_EV,
                 dd/np*UENERGY_IN_EV, dd/np*UENERGY_IN_KJ__MOL);
    }
    Fcr(ft);

    if ( N->track_dx )
    {
        cc = MSD();
        Fprintf (ft, "In this simulation of %g ps, each atom has on average"
                 "drifted\n%g A -> self-diffusion coefficient = %e [m^2/s]."
                 "\n\n", t_in_g*g*UTIME_IN_PS, sqrt(cc)*ULENGTH_IN_A,
                 (cc - final_avg_start.MSD) /
                 (t_in_g - final_avg_start.t_in_g) / g *
                 UAREA_IN_M2 / UTIME_IN_S / 6.);
    }
    if (CALC_GR)
    {
        Gr_Save (Config_Aapp_to_Alib, ct, GR, ft);
        Fcr (ft);
    }
    if (CALC_CONDUCTIVITY)
    {
        spectral_closeflux (SPECTRAL_TYPE_HEAT);
        Fcr (ft);
        smile (SPECTRAL_TYPE_HEAT, ft);
        Fcr (ft);
    }
    if (CALC_VISCOSITY)
    {
        spectral_closeflux (SPECTRAL_TYPE_STRESS);
        Fcr (ft);
        smile (SPECTRAL_TYPE_STRESS, ft);
        Fcr (ft);
    }

    Fcr (ft);
    Neighborlist_print_statistics (np, N, ft);
    Fcr (ft);

    Fprintf (ft, "MD simulation ends on %s.\n", date_string());
    Fcr (ft);
    cc = taskwatch ("MD", TASKWATCH_CHECK);
    Fprintf (ft, "This simulation took %s, which converts to\n"
             "%g step/s, or %e particle x step / s.\n\n",
             earthtime(cc), md_in_g/cc, md_in_g*np/cc);

    fbrk ();
    fclose (fp_screen);
    return (0);
} /* end main() */
#endif /* ifndef _TEST */


#ifdef _hc_TEST
void energy_moment (V3 em)
{
    register int i;
    V3 x;
    V3ZERO (em);
    for (i=np; i--;)
    {
        V3mM3 (&(s[3*i]), H, x);
        V3ADDmuL (ep[i], x, em);
    }
    return;
} /* end energy_moment() */
int main (int argc, char *argv[])
{
    register int i;
    double cc, *ss, *ss1, *ff, g = 5e-5;
    V3 em_1, em1, tmp, hc0;
    M3 HI;
    ft = stdout;
    Config_rebuild_Xtal
        (Config_Aapp_to_Alib, REFERENCE_NC,REFERENCE_NC,REFERENCE_NC,
         Xtal_abstracts+XTAL_NACL, "Zr", -1., "C", -1., 4.711804);
    Config_Opencell (Config_Aapp_to_Alib, 2.);
    Config_perturb_position (Config_Aapp_to_Alib, 0.05);
    rebind_ct (Config_Aapp_to_Alib, "Zr C", ct, &tp, ft); Fcr(ft);
    Neighborlist_Recreate_Form (Config_Aapp_to_Alib, ct, N);
    N->max_atom_displacement = DEF_MAX_ATOM_DISPLACEMENT;
    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, ft, ct, &tp, N);
    Fcr(ft);
    motion = Config_motion_initialize (2000., Config_Aapp_to_Alib);
    REALLOC (main, f,  DIMENSION * np, double);
    REALLOC (main, v,  DIMENSION * np, double);
    REALLOC (main, ep, DIMENSION * np, double);
    NOC = TRUE;
    potential_realloc(np);
    potential_current();
    V3EQV (hc, hc0);
    VCLONE (DIMENSION * np, s,  ss);
    VCLONE (DIMENSION * np, s1, ss1);
    VCLONE (DIMENSION * np, f,  ff);
    M3InV (H, HI, cc);
    for (i=np; i--;)
    {
        V3mM3 ( &(ff[DIMENSION*i]), HI, tmp );
        V3DiV (tmp, mass[i]);
        V3ADDMUL ( &(ss[3*i]),  -g/2, &(ss1[3*i]), &(s[3*i])  );
        V3ADDMUL ( &(ss1[3*i]), -g/2, tmp,         &(s1[3*i]) );
    }
    potential_current();
    energy_moment (em_1);
    for (i=np; i--;)
    {
        V3mM3 ( &(ff[DIMENSION*i]), HI, tmp );
        V3DiV (tmp, mass[i]);
        V3ADDMUL ( &(ss[3*i]),  g/2, &(ss1[3*i]), &(s[3*i])  );
        V3ADDMUL ( &(ss1[3*i]), g/2, tmp,         &(s1[3*i]) );
    }
    potential_current();
    energy_moment (em1);
    for (i=0; i<DIMENSION; i++)
        printf ("%e %e\n", (em1[i]-em_1[i])/g, hc0[i]);
    save_config ("hc.cfg");
    return (0);
}
#endif /* _hc_TEST */