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

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

#define NELEMS(x)  (sizeof(x) / sizeof(x[0]))


/* needed to define type Result to interface */
static const MD_Member resultMemberList[] = {
  { MD_DOUBLE, 1, "energy" },
  { MD_DOUBLE, 1, "ke" },
  { MD_DOUBLE, 1, "pe" },     /* "pe" through "bound" must be in order */
  { MD_DOUBLE, 1, "bond" },   /* corresponding to ForceEnergy type */
  { MD_DOUBLE, 1, "angle" },
  { MD_DOUBLE, 1, "dihed" },
  { MD_DOUBLE, 1, "impr" },
  { MD_DOUBLE, 1, "elec" },
  { MD_DOUBLE, 1, "vdw" },
  { MD_DOUBLE, 1, "bound" },
  { MD_DOUBLE, 1, "temp" },
  { MD_DVEC, 1, "linmo" },  /* linear momentum */
};


/* needed to define type Param to interface */
static const MD_Member paramMemberList[] = {
  { MD_DOUBLE, 1, "cutoff" },
  { MD_DOUBLE, 1, "switchdist" },
  { MD_CHAR, 4, "switching" },
  { MD_CHAR, 12, "exclude" },
  { MD_DOUBLE, 1, "dielectric" },
  { MD_DOUBLE, 1, "1-4scaling" },    /* for "scaling" member of Param */
  { MD_CHAR, 4, "fulldirect" },
  { MD_CHAR, 4, "sphericalbc" },
  { MD_DVEC, 1, "sphericalbccenter" },
  { MD_DOUBLE, 1, "sphericalbcr1" },
  { MD_DOUBLE, 1, "sphericalbcr2" },
  { MD_DOUBLE, 1, "sphericalbck1" },
  { MD_DOUBLE, 1, "sphericalbck2" },
  { MD_INT32, 1, "sphericalbcexp1" },
  { MD_INT32, 1, "sphericalbcexp2" },
  { MD_CHAR, 4, "cylindricalbc" },
  { MD_CHAR, 4, "cylindricalbcaxis" },
  { MD_DVEC, 1, "cylindricalbccenter" },
  { MD_DOUBLE, 1, "cylindricalbcr1" },
  { MD_DOUBLE, 1, "cylindricalbcr2" },
  { MD_DOUBLE, 1, "cylindricalbcl1" },
  { MD_DOUBLE, 1, "cylindricalbcl2" },
  { MD_DOUBLE, 1, "cylindricalbck1" },
  { MD_DOUBLE, 1, "cylindricalbck2" },
  { MD_INT32, 1, "cylindricalbcexp1" },
  { MD_INT32, 1, "cylindricalbcexp2" },
  { MD_DVEC, 1, "cellbasisvector1" },
  { MD_DVEC, 1, "cellbasisvector2" },
  { MD_DVEC, 1, "cellbasisvector3" },
  { MD_DVEC, 1, "cellorigin" },
/*** mgrid params ***/
  { MD_CHAR, 8, "mgrid" },
  { MD_DOUBLE, 1, "mgridlength" },
  { MD_DOUBLE, 1, "mgridcutoff" },
  { MD_DOUBLE, 1, "mgridspacing" },
  { MD_INT32, 1, "mgridnspacings" },
  { MD_INT32, 1, "mgridnlevels" },
  { MD_CHAR, 12, "mgridapprox" },
  { MD_CHAR, 12, "mgridsplit" },
};


/* default values for params */
static const Param defaultParam = {
  /* cutoff = */ 0.0,
  /* switchDist = */ 0.0,
  /* switching = */ "off",
  /* exclude = */ "",
  /* dielectric = */ 1.0,
  /* 1-4scaling = */ 1.0,
  /* fullDirect = */ "no",
  /* sphericalBC = */ "off",
  /* sphericalBCCenter = */ { 0.0, 0.0, 0.0 },
  /* sphericalBCr1 = */ 0.0,
  /* sphericalBCr2 = */ 0.0,
  /* sphericalBCk1 = */ 0.0,
  /* sphericalBCk2 = */ 0.0,
  /* sphericalBCexp1 = */ 2,
  /* sphericalBCexp2 = */ 2,
  /* cylindricalBC = */ "off",
  /* cylindricalBCAxis = */ "",
  /* cylindricalBCCenter = */ { 0.0, 0.0, 0.0 },
  /* cylindricalBCr1 = */ 0.0,
  /* cylindricalBCr2 = */ 0.0,
  /* cylindricalBCl1 = */ 0.0,
  /* cylindricalBCl2 = */ 0.0,
  /* cylindricalBCk1 = */ 0.0,
  /* cylindricalBCk2 = */ 0.0,
  /* cylindricalBCexp1 = */ 2,
  /* cylindricalBCexp2 = */ 2,
  /* cellBasisVector1 = */ { 0.0, 0.0, 0.0 },
  /* cellBasisVector2 = */ { 0.0, 0.0, 0.0 },
  /* cellBasisVector3 = */ { 0.0, 0.0, 0.0 },
  /* cellOrigin = */ { 0.0, 0.0, 0.0 },
/*** mgrid params ***/
  /* mgrid = */ "off",
  /* mgridLength = */ 0.0,
  /* mgridCutoff = */ 0.0,
  /* mgridSpacing = */ 0.0,
  /* mgridNspacings = */ 0,
  /* mgridNlevels = */ 0,
  /* mgridApprox = */ "cubic",
  /* mgridSplit = */ "taylor2",
};


int32 deven_init(MD_Front *f, int32 flags)
{
  /* pointer to internal data */
  Engine *eng = NULL;

  /* access permissions for engine data arrays */
  const int32 aw = MD_RESIZE | MD_WRITE | MD_NOTIFY;
  const int32 arwc = MD_RESIZE | MD_READ | MD_WRITE | MD_CBREAD | MD_NOTIFY;
  const int32 arc = MD_ERESIZE | MD_READ | MD_CBREAD;
  const int32 ae = MD_ERESIZE;

  /* attributes static engine data arrays */
  /* for result and param, change 0 to type number after definition */
  MD_Attrib attr_result = { 0, 1, 1, MD_READ | MD_CBREAD };
  MD_Attrib attr_param = { 0, 1, 1, MD_READ | MD_WRITE | MD_NOTIFY };
  MD_Attrib attr_timestep = { MD_DOUBLE, 1, 1, MD_READ | MD_WRITE };
  MD_Attrib attr_init_temp
    = { MD_DOUBLE, 1, 1, MD_READ | MD_WRITE | MD_NOTIFY };
  MD_Attrib attr_com_motion = { MD_CHAR, 0, 4, MD_READ | MD_WRITE | MD_SETLEN };
  MD_Attrib attr_seed = { MD_INT32, 1, 1, MD_READ | MD_WRITE };
  MD_Attrib attr_ndegfreedom = { MD_INT32, 1, 1, MD_READ };

  /* need clock time to init seed for random number generator */
  struct timeval tv;

  /* allocate internal data structure */
  if ((eng = (Engine *) calloc(1, sizeof(Engine))) == NULL) {
    return MD_error(f, MD_ERR_MEMALLOC);
  }

  /* setup engine and internal data with interface */
  if (MD_setup_engine(f, eng)) return MD_FAIL;

  /* define Result and Param types with interface */
  attr_result.type = MD_new_type(f, "Result", resultMemberList,
      NELEMS(resultMemberList), sizeof(Result));
  if (attr_result.type == MD_FAIL) return MD_FAIL;

  attr_param.type = MD_new_type(f, "Param", paramMemberList,
      NELEMS(paramMemberList), sizeof(Param));
  if (attr_param.type == MD_FAIL) return MD_FAIL;

  /* define errors with interface */
  eng->err_param = MD_new_error(f, "Invalid engine parameters", 0);
  if (eng->err_param == MD_FAIL) return MD_FAIL;

  eng->err_unsupported = MD_new_error(f, "Unsupported engine feature", 0);
  if (eng->err_unsupported == MD_FAIL) return MD_FAIL;

  eng->err_force_init = MD_new_error(f, "Force initialization failed", -1);
  if (eng->err_force_init == MD_FAIL) return MD_FAIL;

  eng->err_force = MD_new_error(f, "Force computation failed", -1);
  if (eng->err_force == MD_FAIL) return MD_FAIL;

  /* set engine params to default values */
  eng->param = defaultParam;

  /* set default values for additional sim params */
  eng->timestep = 1.0;
  eng->init_temp = 0.0;
  strcpy(eng->com_motion, "no");
  attr_com_motion.len = 3;  /* strlen("no") + nil-terminator */

  /* choose random seed using clock */
  if (gettimeofday(&tv, NULL)) return MD_FAIL;
  eng->seed = tv.tv_sec;

  /* establish engine data arrays with interface */
  if (/* force field params and topology arrays allow write and resize */
      (eng->atomprm = MD_engdata(f, "atomprm", MD_ATOMPRM, aw)) == NULL
      || (eng->bondprm = MD_engdata(f, "bondprm", MD_BONDPRM, aw)) == NULL
      || (eng->angleprm = MD_engdata(f, "angleprm", MD_ANGLEPRM, aw)) == NULL
      || (eng->dihedprm = MD_engdata(f, "dihedprm", MD_TORSPRM, aw)) == NULL
      || (eng->imprprm = MD_engdata(f, "imprprm", MD_TORSPRM, aw)) == NULL
      || (eng->nbfixprm = MD_engdata(f, "nbfixprm", MD_NBFIXPRM, aw)) == NULL
      || (eng->atom = MD_engdata(f, "atom", MD_ATOM, aw)) == NULL
      || (eng->bond = MD_engdata(f, "bond", MD_BOND, aw)) == NULL
      || (eng->angle = MD_engdata(f, "angle", MD_ANGLE, aw)) == NULL
      || (eng->dihed = MD_engdata(f, "dihed", MD_TORS, aw)) == NULL
      || (eng->impr = MD_engdata(f, "impr", MD_TORS, aw)) == NULL
      || (eng->excl = MD_engdata(f, "excl", MD_EXCL, aw)) == NULL
      /* pos and vel allow read, write, and callback read */
      || (eng->pos = MD_engdata(f, "pos", MD_DVEC, arwc)) == NULL
      || (eng->vel = MD_engdata(f, "vel", MD_DVEC, arwc)) == NULL
      /* engine resized force array allows read and callback read */
      || (eng->force_engdata = MD_engdata(f, "force", MD_DVEC, arc)) == NULL
      /* result and param point to pre-allocated memory */
      || (eng->result_engdata = MD_engdata_buffer(f, "result", attr_result,
          &(eng->result))) == NULL
      || (eng->param_engdata = MD_engdata_buffer(f, "param", attr_param,
          &(eng->param))) == NULL
      /* additional sim params */
      || (eng->timestep_engdata = MD_engdata_buffer(f, "timestep",
          attr_timestep, &(eng->timestep))) == NULL
      || (eng->init_temp_engdata = MD_engdata_buffer(f, "temperature",
          attr_init_temp, &(eng->init_temp))) == NULL
      || (eng->com_motion_engdata = MD_engdata_buffer(f, "commotion",
          attr_com_motion, &(eng->com_motion))) == NULL
      || (eng->seed_engdata = MD_engdata_buffer(f, "seed", attr_seed,
          &(eng->seed))) == NULL
      /* additional results */
      || (eng->numdegreesfreedom = MD_engdata_buffer(f, "numdegreesfreedom",
          attr_ndegfreedom, &(eng->ndegfreedom))) == NULL
    /*** data arrays for mgrid - permit only engine resize ***/
      || (eng->wrap = MD_engdata(f, "wrap", MD_DVEC, ae)) == NULL
      || (eng->poswrap = MD_engdata(f, "poswrap", MD_DVEC, ae)) == NULL
      || (eng->f_elec = MD_engdata(f, "f_elec", MD_DVEC, ae)) == NULL
      || (eng->charge = MD_engdata(f, "charge", MD_DOUBLE, ae)) == NULL
    /*** for debugging electrostatics ***/
   /* || (eng->force_elec = MD_engdata(f, "force_elec", MD_DVEC, ae)) == NULL */
      ) {
    return MD_FAIL;
  }

  /* setup run routine with interface */
  if (MD_setup_run(f, deven_run)) return MD_FAIL;

  /* init force library */
  if (force_init(&(eng->force))) return MD_FAIL;

  /* init mgrid library */
  if (mgrid_init(&(eng->mgrid))) return MD_FAIL;

  /*** step library is initialized later ***/

  printf("# deven started successfully\n");

  return 0;
}


void deven_done(MD_Front *f)
{
  Engine *eng = MD_engine_data(f);

  /* done with mgrid library */
  mgrid_done(&(eng->mgrid));

  /* done with force library */
  force_done(&(eng->force));

  /* free memory allocated during init */
  free(eng);

  printf("# deven stopped successfully\n");
}