/*
 * data.c
 *
 * Contains init, destroy, and validate routines for deven 'Data'
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mdutil.h"
#include "data.h"
#include "force.h"


/*
 * Exclusion strings
 */
static const char *EXCL_STRING[] = {
  "none",
  "1-2",
  "1-3",
  "1-4",
  "scaled1-4"
};


/*
 * Electrostatic evaluation strings
 */
const char *ELEC_EVAL_STRING[] = {
  "Cutoff",
  "MG",
  "Mgrid",
  "DPMTA",
  "FullDirect",
};


/*
 * MG switch type strings
 */
const char *SWITCH_TYPE_STRING[] = {
  "dummy",
  "C1",
  "C2",
  "C3",
  "dummy"
};


/*
 * Define init for Data.
 *
 * Returns 0 on success, MD_FAIL on failure.
 */
MD_Errcode deven_data_init(Data *data)
{
  MD_Sim *sim = data->sim;

  if ((data->force = (Force *) calloc(1, sizeof(Force))) == NULL) {
    MD_ERR(sim, MD_engine_name(sim), MD_ERR_MEMALLOC,
           "deven_data_init call to calloc");
    return MD_FAIL;
  }

  /* initialize reduction labels */
  data->nreducs = NREDUCS;
  data->p_reduc = &(data->reduc[0]);
  data->p_reduclabel = &(data->reduclabel[0]);
  strncpy(data->reduclabel[PE_BOND],  "bond",          sizeof(MD_String)-1);
  strncpy(data->reduclabel[PE_ANGLE], "angle",         sizeof(MD_String)-1);
  strncpy(data->reduclabel[PE_DIHED], "dihedral",      sizeof(MD_String)-1);
  strncpy(data->reduclabel[PE_IMPR],  "improper",      sizeof(MD_String)-1);
  strncpy(data->reduclabel[PE_ELEC],  "electrostatic", sizeof(MD_String)-1);
  strncpy(data->reduclabel[PE_VDW],   "van der Waals", sizeof(MD_String)-1);
  strncpy(data->reduclabel[PE_BC],    "boundary",      sizeof(MD_String)-1);
  strncpy(data->reduclabel[PE],       "potential",     sizeof(MD_String)-1);
  strncpy(data->reduclabel[KE],       "kinetic",       sizeof(MD_String)-1);
  strncpy(data->reduclabel[TE],       "energy",        sizeof(MD_String)-1);
  strncpy(data->reduclabel[SHADOW2],  "shadow2",       sizeof(MD_String)-1);
  strncpy(data->reduclabel[SHADOW4],  "shadow4",       sizeof(MD_String)-1);
  strncpy(data->reduclabel[SHADOW6],  "shadow6",       sizeof(MD_String)-1);
  strncpy(data->reduclabel[SHADOW8],  "shadow8",       sizeof(MD_String)-1);

  /* initialize simulation parameters */
  data->simparam_len = 1;

  data->timestep = DEFAULT_TIMESTEP;
  data->p_timestep = &(data->timestep);

  data->dielectric = DEFAULT_DIELECTRIC;
  data->p_dielectric = &(data->dielectric);

  data->scaling14 = DEFAULT_SCALING14;
  data->p_scaling14 = &(data->scaling14);

  data->switchdist = DEFAULT_SWITCHDIST;
  data->p_switchdist = &(data->switchdist);

  data->cutoff = DEFAULT_CUTOFF;
  data->p_cutoff = &(data->cutoff);

  strncpy(data->exclude, DEFAULT_EXCLUDE, sizeof(MD_String)-1);
  data->p_exclude = &(data->exclude);

  data->is_switching = DEFAULT_IS_SWITCHING;
  data->p_is_switching = &(data->is_switching);

  data->stepnum = DEFAULT_STEPNUM;
  data->p_stepnum = &(data->stepnum);

  data->firststepnum = DEFAULT_FIRSTSTEPNUM;
  data->p_firststepnum = &(data->firststepnum);

  strncpy(data->elec_eval, DEFAULT_ELEC_EVAL, sizeof(MD_String)-1);
  data->p_elec_eval = &(data->elec_eval);

  data->is_shadow = DEFAULT_IS_SHADOW;
  data->p_is_shadow = &(data->is_shadow);

  data->mg_divisions = DEFAULT_MG_DIVISIONS;
  data->p_mg_divisions = &(data->mg_divisions);

  data->mg_gridsize = DEFAULT_MG_GRIDSIZE;
  data->p_mg_gridsize = &(data->mg_gridsize);

  data->mg_cutoff = DEFAULT_MG_CUTOFF;
  data->p_mg_cutoff = &(data->mg_cutoff);

  data->mg_levels = DEFAULT_MG_LEVELS;
  data->p_mg_levels = &(data->mg_levels);

  strncpy(data->mg_switchtype_string, DEFAULT_MG_SWITCHTYPE_STRING,
      sizeof(MD_String)-1);
  data->p_mg_switchtype_string = &(data->mg_switchtype_string);

  data->dpmta_levels = DEFAULT_DPMTA_LEVELS;
  data->p_dpmta_levels = &(data->dpmta_levels);

  data->dpmta_mpterms = DEFAULT_DPMTA_MPTERMS;
  data->p_dpmta_mpterms = &(data->dpmta_mpterms);

  data->dpmta_is_fft = DEFAULT_DPMTA_IS_FFT;
  data->p_dpmta_is_fft = &(data->dpmta_is_fft);

  data->dpmta_blocking_fft = DEFAULT_DPMTA_BLOCKING_FFT;
  data->p_dpmta_blocking_fft = &(data->dpmta_blocking_fft);

  data->dpmta_theta = DEFAULT_DPMTA_THETA;
  data->p_dpmta_theta = &(data->dpmta_theta);

  data->is_sphericalbc = DEFAULT_IS_SPHERICALBC;
  data->p_is_sphericalbc = &(data->is_sphericalbc);

  data->sphericalbc_center.x = DEFAULT_SPHERICALBC_CENTER_X;
  data->sphericalbc_center.y = DEFAULT_SPHERICALBC_CENTER_Y;
  data->sphericalbc_center.z = DEFAULT_SPHERICALBC_CENTER_Z;
  data->p_sphericalbc_center = &(data->sphericalbc_center);

  data->sphericalbc_r1 = DEFAULT_SPHERICALBC_R1;
  data->p_sphericalbc_r1 = &(data->sphericalbc_r1);

  data->sphericalbc_k1 = DEFAULT_SPHERICALBC_K1;
  data->p_sphericalbc_k1 = &(data->sphericalbc_k1);

  data->sphericalbc_exp1 = DEFAULT_SPHERICALBC_EXP1;
  data->p_sphericalbc_exp1 = &(data->sphericalbc_exp1);

  data->sphericalbc_r2 = DEFAULT_SPHERICALBC_R2;
  data->p_sphericalbc_r2 = &(data->sphericalbc_r2);

  data->sphericalbc_k2 = DEFAULT_SPHERICALBC_K2;
  data->p_sphericalbc_k2 = &(data->sphericalbc_k2);

  data->sphericalbc_exp2 = DEFAULT_SPHERICALBC_EXP2;
  data->p_sphericalbc_exp2 = &(data->sphericalbc_exp2);

  return 0;
}


/*
 * Define destroy for Data.
 */
void deven_data_destroy(Data *data)
{
  MD_Int k;

  if (data->force) deven_force_destroy(data->force);
  free(data->force);
  free(data->atomprm);
  free(data->bondprm);
  free(data->angleprm);
  free(data->dihedprm);
  free(data->imprprm);
  free(data->nbfixprm);
  free(data->atom);
  free(data->bond);
  free(data->angle);
  free(data->dihed);
  free(data->impr);
  free(data->excl);
  for (k = 0;  k < NUMARR;  k++) {
    free(data->posarr[k]);
    free(data->velarr[k]);
  }
  free(data->f);
}


/*
 * validate the simulation parameters
 */
MD_Errcode deven_data_validate_simparam(Data *data)
{
  MD_Sim *sim = data->sim;
  MD_Int k;
  char buf[BUFSIZE];

  if (data->dielectric < 1.0) {
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam,
           "deven_data_validate_simparam: "
           "dielectric constant should be >= 1");
    return MD_FAIL;
  }
  if (data->scaling14 > 1.0 || data->scaling14 < 0.0) {
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam,
           "deven_data_validate_simparam: "
           "1-4 scaling should be between 0 and 1 inclusive");
    return MD_FAIL;
  }
  if (data->cutoff <= 0.0) {
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam,
           "deven_data_validate_simparam: "
           "cutoff should be > 0");
    return MD_FAIL;
  }
  if (data->switchdist > data->cutoff) {
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam,
           "deven_data_validate_simparam: "
           "switchdist should be <= cutoff");
    return MD_FAIL;
  }

  /* set exclusion policy */
  data->excl_policy = -1;
  for (k = EXCL_NONE;  k <= SCALED_14;  k++) {
    if (mdutil_strcasecmp(data->exclude, EXCL_STRING[k]) == 0) {
      data->excl_policy = k;
    }
  }
  if (data->excl_policy < 0) {
    /* issue an error, no valid string matched */
    sprintf(buf, "deven_data_validate_simparam: "
            "unrecognized exclusion policy: \"%s\"", data->exclude);
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam, buf);
    return MD_FAIL;
  }

  /* set electrostatic evaluation */
  data->elec_eval_policy = -1;
  for (k = ELEC_EVAL_CUTOFF;  k <= ELEC_EVAL_FULLDIRECT;  k++) {
    if (mdutil_strcasecmp(data->elec_eval, ELEC_EVAL_STRING[k]) == 0) {
      data->elec_eval_policy = k;
    }
  }
  if (data->elec_eval_policy < 0) {
    /* issue an error, no valid string matched */
    sprintf(buf, "deven_data_validate_simparam: "
            "unrecognized electrostatic evaluation policy: \"%s\"",
            data->elec_eval);
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam, buf);
    return MD_FAIL;
  }

  if (data->elec_eval_policy == ELEC_EVAL_MG
      || data->elec_eval_policy == ELEC_EVAL_MGRID
      || data->elec_eval_policy == ELEC_EVAL_DPMTA) {
    /* issue an error, unsupported feature in this version */
    sprintf(buf, "deven_data_validate_simparam: "
            "this version does not support "
            "electrostatic evaluation policy: \"%s\"",
            data->elec_eval);
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam, buf);
    return MD_FAIL;
  }

  /* set MG switch type string */
  data->mg_switchtype = -1;
  for (k = MG_SWITCH_FIRST + 1;  k < MG_SWITCH_LAST;  k++) {
    if (mdutil_strcasecmp(data->mg_switchtype_string, 
          SWITCH_TYPE_STRING[k]) == 0) {
      data->mg_switchtype = k;
    }
  }
  if (data->mg_switchtype < 0) {
    /* issue an error, no valid string matched */
    sprintf(buf, "deven_data_validate_simparam: "
            "unrecognized MG switch type string: \"%s\"",
            data->mg_switchtype_string);
    MD_ERR(sim, MD_engine_name(sim), data->err_simparam, buf);
    return MD_FAIL;
  }

  return 0;
}