#define SHOULD_HALT \
  ((stage==0)&&(t_in_g>=0)&&(t_in_g%halt_in_g==0))

#define SHOULD_REPORT \
  ((stage==0)&&(t_in_g>=0)&&(t_in_g%report_freq_in_g==0))

#define SHOULD_FINAL_AVG \
  ((stage==0)&&(t_in_g>=final_avg_start_in_g)&&(t_in_g%halt_in_g==0))

#define SHOULD_CALC_CONDUCTIVITY (CALC_CONDUCTIVITY && \
  SHOULD_HALT && SHOULD_FINAL_AVG)

#define SHOULD_CALC_VISCOSITY (CALC_VISCOSITY && \
  SHOULD_HALT && SHOULD_FINAL_AVG)

void engine (int stage, int freedom, double *s, double *f)
{
    register int i;
    double cc, volume, kine_wall;
    V3 tmp;
    M3 U, K, KI, HI, UD, UDKI, GDGI, sout;

    if (t_in_g < 0)
    {
        VZERO (freedom, f);
        return;
    }
    if (freedom == DIMENSION * np)
    {
        if (NTH)
        {
            get_scheduled_symmat (eta_schedules, t_in_g, U);
            pure_deform (H0, U, H);
        }
        Neighborlist_Maintenance (Config_Aapp_to_Alib, ft, ct, &tp, N);
        if ( SHOULD_HALT )
            if (SHOULD_CALC_CONDUCTIVITY) potential_current();
            else if (SHOULD_CALC_VISCOSITY) potential_stress();
            else if (SHOULD_REPORT) potential_stress();
            else if (SHOULD_FINAL_AVG) potential_stress();
            else potential_force();
        else potential_force();
        M3InV (H, HI, volume);
        for (i=np; i--;)
            if (MOBILE_ATOM(i))
                V3MM3DIV ( &(f[DIMENSION*i]), HI, mass[i], tmp);
    }
    else
    {
        if (freedom == DIMENSION * np + 1)
        {
            M3IdentitY ( s[DIMENSION*np],  U  );
            M3IdentitY ( s1[DIMENSION*np], UD );
        }
        else
        {
            SYMMAT_DOWNLOAD ( &(s[DIMENSION*np]),  U  );
            SYMMAT_DOWNLOAD ( &(s1[DIMENSION*np]), UD );
        }
        M3ADDDIAG (U, 1, K);
        M3MUL (H0, K, H);
        Neighborlist_Maintenance (Config_Aapp_to_Alib, ft, ct, &tp, N);
        if ( SHOULD_HALT )
            if (SHOULD_CALC_CONDUCTIVITY) potential_current();
            else potential_stress();
        else potential_stress();
        M3INV (K, KI, cc);
        M3InV (H, HI, volume);
        M3MUL (UD, KI, UDKI);
        M3Symmetrize2 (UDKI);
        M3MUL2 (H, UDKI, HI, GDGI, sout);
        for (i=np; i--;)
            if (MOBILE_ATOM(i))
            {
                V3MM3DIV ( &(f[DIMENSION*i]), HI, mass[i], tmp );
                V3SUBmM3 ( &(f[DIMENSION*i]), &(s1[DIMENSION*i]), GDGI );
            }
        get_scheduled_symmat (sout_schedules, t_in_g, sout);
        M3ADD (stress, sout, UDKI);
        M3MUL (UDKI, KI, GDGI);
        M3Symmetrize (GDGI);
        cc = volume / wall;
        if (freedom == DIMENSION * np + 1)
            f[DIMENSION*np] = cc * M3TR(GDGI) / 3;
        else
            SYMMAT_MULUPLOAD ( cc, GDGI, &(f[DIMENSION*np]) );
    }
    if ( SHOULD_HALT )
    {
        desired_T.real = schedulelist_value (T_schedule, t_in_g);
        desired_T.md = Treal_to_Tmd (desired_T.real, &desired_T.dTmd__dTreal);
        actual_T.md = motion->T_in_K;
        actual_T.real = Tmd_to_Treal (actual_T.md, &actual_T.dTmd__dTreal);
        if (SHOULD_CALC_CONDUCTIVITY)
        {
            /* Fprintf (ft, "%e\n", hc[0]* */
            /* sqrt( UKAPPA_IN_W__M_K * desired_T.dTmd__dTreal * */
            /* halt_in_g * g / desired_T.real / volume / */
            /* (BOLZ_IN_UENERGY__K * desired_T.real) )); */
            spectral_writeflux
                (hc, sqrt( UKAPPA_IN_W__M_K * desired_T.dTmd__dTreal *
                           halt_in_g * g / desired_T.real / volume /
                           (BOLZ_IN_UENERGY__K * desired_T.real) ),
                 SPECTRAL_TYPE_HEAT);
        }
        if (SHOULD_CALC_VISCOSITY)
        {
            V3ASSIGN ( stress[1][2], stress[0][2], stress[0][1], hc );
            spectral_writeflux
                (hc, sqrt( UVISCOSITY_IN_PA_S * halt_in_g * g * volume /
                           (BOLZ_IN_UENERGY__K * desired_T.real) ),
                 SPECTRAL_TYPE_STRESS);
        }
        if (freedom > DIMENSION * np)
        {
            kine_wall = wall * M32NORM2 ( UD ) / 2;
            cc = sqrt ( (kine_wall / NTS_MAX_WALL_KINE_IN_KT /
                         BOLZ_IN_UENERGY__K + JUMPSTART_T_IN_K) /
                        (desired_T.md + JUMPSTART_T_IN_K) );
            if ( (cc > 1) && (coupling_time_constant > 0) )
            {
                M3DividE (UD, cc);
                if (freedom == DIMENSION * np + 1)
                    s1[DIMENSION*np] = UD[0][0];
                else SYMMAT_UPLOAD( UD, &(s1[DIMENSION*np]) );
            }
        }
        if ( ( (t_in_g < final_avg_start_in_g) ||
               (freedom > DIMENSION * np) ) &&
             (coupling_time_constant > 0) )
        {
            cc = 1. - halt_in_g * g / coupling_time_constant *
                ( (actual_T.md  + JUMPSTART_T_IN_K) /
                  (desired_T.md + JUMPSTART_T_IN_K) - 1. );
            if (cc < 0.5) cc = 0.5;
            else if (cc > 2) cc = 2;
            motion = Config_motion_scale_T (cc, motion, Config_Aapp_to_Alib);
        }
    }
    if ( SHOULD_REPORT )
    {
        if (N->track_dx) DRAW ( MSD, MSD()*UAREA_IN_A2 );
        M3DRAWMUL ( H, H, ULENGTH_IN_A );
        DRAW ( T, actual_T.real );
        DRAW ( Tmd, actual_T.md );
        DRAW ( pote, pote * UENERGY_IN_EV );
        SYMMAT_DRAWMUL ( t, stress, USTRESS_IN_MPA );
        fflush (fp_draw);
        Fprintf (ft, "time = %d [g] = %d [halt] = %g = %g [ps]:\n",
                 t_in_g, t_in_g / halt_in_g, t_in_g * g,
                 t_in_g * g * UTIME_IN_PS);
        S3fPRmul (ft, "H = %M [A];\n", H, ULENGTH_IN_A);
        Fprintf (ft, "Treal = %g [K], Tmd = %g [K];\n",
                 actual_T.real, actual_T.md);
        Fprintf (ft, "pote = %.8g [eV] = %.8g [eV/atom];\n",
                 pote * UENERGY_IN_EV, pote / np * UENERGY_IN_EV);
        cc = pote + motion->kine;
        Fprintf (ft, "tote = %.8g [eV] = %.8g [eV/atom];\n",
                 cc * UENERGY_IN_EV, cc / np * UENERGY_IN_EV);
        if (freedom > DIMENSION * np)
        {
            cc += kine_wall - work(H0, sout, U);
            Fprintf (ft, "Hamiltonian = %.8g [eV] = %.8g [eV/atom];\n",
                     cc * UENERGY_IN_EV, cc / np * UENERGY_IN_EV);
        }
        S3fPRmul (ft, "stress = %M [MPa].", stress, USTRESS_IN_MPA);
        Fcr(ft);
    }
    if ( SHOULD_FINAL_AVG )
        avg_add ( Avg_what(Property, final), M3e(H), motion->kine,
                  pote, M3e(stress) );
    if ((t_in_g >= 0) && (stage < newtonian[integrator_id].K)) t_in_g++;
    /* Neighborlist_Check (Config_Aapp_to_Alib, ft, ct, &tp, N); */
    return;
} /* end engine() */