/*
 * Copyright (C) 2004-2005 by David J. Hardy.  All rights reserved.
 *
 * run.c
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
#include "mdapi/mdengine.h"
#include "force/force.h"
#include "deven/engine.h"
#include "debug/debug.h"


/* prototypes for internal functions */
static int reductions(Engine *);


int32 deven_run(MD_Front *front, int32 numsteps, int32 flags)
{
  Engine *eng = MD_engine_data(front);
  Step *step = &(eng->step);
  StepParam *step_param = &(eng->step_param);
  StepSystem *step_system = &(eng->step_system);
  int32 n;

  /* update internal data and force library */
  if (flags & MD_UPDATE) {

    /* call update routine */
    if (deven_update(front)) return MD_FAIL;
  }

  /*** take advantage of fact that step object points to step param ***/
  /*** (this is a kludge!) ***/
  step_param->timestep = eng->timestep;

  /* see if callbacks need to be processed before stepping */
  if (MD_ready_callback(front)) {

    /* make sure that forces are up-to-date */
    if (step_compute(step, step_system, 0)) {
      return MD_FAIL;
    }

    /* reductions might be needed by callback */
    if (reductions(eng)) return MD_FAIL;

    /* execute callbacks */
    if (MD_exec_callback(front) || MD_wait_callback(front)) return MD_FAIL;
  }

  /* call step library for numsteps */
  for (n = 0;  n < numsteps;  n++) {

    /* take single steps so we can monitor callbacks */
    if (step_compute(step, step_system, 1)) {
      return MD_FAIL;
    }

    /* increment step counter with interface */
    MD_incrstep(front);

    /* see if callbacks need to be processed for this step */
    if (MD_ready_callback(front)) {

      /* reductions might be needed by callback */
      if (reductions(eng)) return MD_FAIL;

      /* execute callbacks */
      if (MD_exec_callback(front) || MD_wait_callback(front)) return MD_FAIL;
    }
  }

  /* compute reductions before returning */
  if (reductions(eng)) return MD_FAIL;

  return 0;
}


int reductions(Engine *e)
{
  Result *result = &(e->result);
  Step *step = &(e->step);
  StepSystem *step_system = &(e->step_system);

  /* adjust energies from mgrid evaluation */
  if (e->ismgrid) {
    result->fe.elec += e->mgrid_system.u_elec;
    result->fe.pe += e->mgrid_system.u_elec;
  }

  /* call step library for kinetic energy and temperature */
  if (step_find_reductions(step, step_system)) {
    return MD_FAIL;
  }

  result->ke = step_system->kinetic_energy;
  result->energy = result->ke + result->fe.pe;
  result->temp = step_system->temperature;
  return 0;
}