/*
 * test_mdio.c
 *
 * Summary:  Test program for MDIO library.
 *
 * Author:  David Hardy
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mdio.h"
#include "pdata.h"

#define  INITLEN  5
#define  TBUFLEN  5

int test_param(Param *p, Pbuf *pb);
int test_topo(Topo *t, Tbuf *tbuf);

int main(int argc, char **argv)
{
  Param p;
  Pbuf pbuf;
  Topo t;
  Tbuf tbuf;
  MD_Atom atom[TBUFLEN];
  MD_Bond bond[TBUFLEN];
  MD_Angle angle[TBUFLEN];
  MD_Dihed dihed[TBUFLEN];
  MD_Impr impr[TBUFLEN];
  MD_Excl excl[TBUFLEN];
  const char *paramfile;
  const char *topofile;

/* get file names from command line */
  if (argc != 3) {
    fprintf(stderr, "syntax:  %s paramfile topofile\n", argv[0]);
    exit(1);
  }
  paramfile = argv[1];
  topofile = argv[2];

/* init Pbuf */
  pbuf.atomp = (MD_Atom_Param *) malloc(INITLEN * sizeof(MD_Atom_Param));
  pbuf.atomp_len = INITLEN;
  pbuf.atomp_num = 0;

  pbuf.bondp = (MD_Bond_Param *) malloc(INITLEN * sizeof(MD_Bond_Param));
  pbuf.bondp_len = INITLEN;
  pbuf.bondp_num = 0;

  pbuf.anglep = (MD_Angle_Param *) malloc(INITLEN * sizeof(MD_Angle_Param));
  pbuf.anglep_len = INITLEN;
  pbuf.anglep_num = 0;

  pbuf.dihedp = (MD_Dihed_Param *) malloc(INITLEN * sizeof(MD_Dihed_Param));
  pbuf.dihedp_len = INITLEN;
  pbuf.dihedp_num = 0;

  pbuf.imprp = (MD_Impr_Param *) malloc(INITLEN * sizeof(MD_Impr_Param));
  pbuf.imprp_len = INITLEN;
  pbuf.imprp_num = 0;

  pbuf.nbfixp = (MD_Nbfix_Param *) malloc(INITLEN * sizeof(MD_Nbfix_Param));
  pbuf.nbfixp_len = INITLEN;
  pbuf.nbfixp_num = 0;

/* init Tbuf */
  tbuf.atom = atom;
  tbuf.atom_len = TBUFLEN;
  tbuf.atom_num = 0;
  tbuf.atom_max = 0;

  tbuf.bond = bond;
  tbuf.bond_len = TBUFLEN;
  tbuf.bond_num = 0;
  tbuf.bond_max = 0;

  tbuf.angle = angle;
  tbuf.angle_len = TBUFLEN;
  tbuf.angle_num = 0;
  tbuf.angle_max = 0;

  tbuf.dihed = dihed;
  tbuf.dihed_len = TBUFLEN;
  tbuf.dihed_num = 0;
  tbuf.dihed_max = 0;

  tbuf.impr = impr;
  tbuf.impr_len = TBUFLEN;
  tbuf.impr_num = 0;
  tbuf.impr_max = 0;

  tbuf.excl = excl;
  tbuf.excl_len = TBUFLEN;
  tbuf.excl_num = 0;
  tbuf.excl_max = 0;

/* test param reading */
  param_init(&p, &pbuf);
  param_open(&p, paramfile, PARAM_XPLOR);
  test_param(&p, &pbuf);
  param_close(&p);

/* test topo reading */
  topo_init(&t, &tbuf, &p);
  topo_open(&t, topofile, TOPO_XPLOR);
  test_topo(&t, &tbuf);
  topo_close(&t);

/* clean up */
  topo_destroy(&t);
  param_destroy(&p);

  free(pbuf.atomp);
  free(pbuf.bondp);
  free(pbuf.anglep);
  free(pbuf.dihedp);
  free(pbuf.imprp);
  free(pbuf.nbfixp);

  return 0;
}


int test_topo(Topo *t, Tbuf *tbuf)
{
  int retval;
  int natom = 0, nbond = 0, nangle = 0, ndihed = 0, nimpr = 0, nexcl = 0;
  int i;

  printf("\n\nList topology from Tbuf:\n\n");
  retval = topo_read(t);
  while (retval != 0) {
    switch (retval) {
    case MD_ATOM:
      printf("atom_num = %d   atom_max = %d\n", 
             tbuf->atom_num, tbuf->atom_max);
      for (i = 0;  i < tbuf->atom_num;  i++) {
        printf("%d  m=%f  q=%f  type=%d  atomname=%s\n", 
               natom++, tbuf->atom[i].m, tbuf->atom[i].q,
               tbuf->atom[i].param, tbuf->atom[i].name);
      }
      break;
    case MD_BOND:
      printf("bond_num = %d   bond_max = %d\n", 
             tbuf->bond_num, tbuf->bond_max);
      for (i = 0;  i < tbuf->bond_num;  i++) {
        printf("%d  atom1=%d  atom2=%d  type=%d\n", 
               nbond++, tbuf->bond[i].atom[0], tbuf->bond[i].atom[1],
               tbuf->bond[i].param);
      }
      break;

    case MD_ANGLE:
      printf("angle_num = %d   angle_max = %d\n", 
             tbuf->angle_num, tbuf->angle_max);
      for (i = 0;  i < tbuf->angle_num;  i++) {
        printf("%d  atom1=%d  atom2=%d  atom3=%d  type=%d\n", 
               nangle++, tbuf->angle[i].atom[0], tbuf->angle[i].atom[1],
               tbuf->angle[i].atom[2], tbuf->angle[i].param);
      }
      break;

    case MD_DIHED:
      printf("dihed_num = %d   dihed_max = %d\n", 
             tbuf->dihed_num, tbuf->dihed_max);
      for (i = 0;  i < tbuf->dihed_num;  i++) {
        printf("%d  atom1=%d  atom2=%d  atom3=%d  atom4=%d  type=%d\n", 
               ndihed++, tbuf->dihed[i].atom[0], tbuf->dihed[i].atom[1],
               tbuf->dihed[i].atom[2], tbuf->dihed[i].atom[3],
               tbuf->dihed[i].param);
      }
      break;

    case MD_IMPR:
      printf("impr_num = %d   impr_max = %d\n", 
             tbuf->impr_num, tbuf->impr_max);
      for (i = 0;  i < tbuf->impr_num;  i++) {
        printf("%d  atom1=%d  atom2=%d  atom3=%d  atom4=%d  type=%d\n", 
               nimpr++, tbuf->impr[i].atom[0], tbuf->impr[i].atom[1],
               tbuf->impr[i].atom[2], tbuf->impr[i].atom[3], 
               tbuf->impr[i].param);
      }
      break;

    case MD_EXCL:
      printf("excl_num = %d   excl_max = %d\n", 
             tbuf->excl_num, tbuf->excl_max);
      for (i = 0;  i < tbuf->excl_num;  i++) {
        printf("%d  atom1=%d  atom2=%d\n", 
               nexcl++, tbuf->excl[i].atom[0], tbuf->excl[i].atom[1]);
      }
      break;

    default:
      fprintf(stderr, "test_topo: an error occurred\n");
      exit(1);
    }

    retval = topo_read(t);
  }

  return 0;
}


int test_param(Param *p, Pbuf *pb)
{
  MD_Atom_Param *at;
  MD_Bond_Param *bn;
  MD_Angle_Param *an;
  MD_Dihed_Param *di;
  MD_Impr_Param *im;
  MD_Nbfix_Param *nb;
  int i, j, k, retval;
  int kat, kbn, kan, kdi, kim, knb;
  char *s, stype[64], a1[64], a2[64], a3[64], a4[64];

  /* read the file */
  retval = param_read(p);
  while (retval > 0) {
    switch (retval) {
    case MD_ATOM_PARAM:
      pb->atomp = realloc(pb->atomp, 
          (pb->atomp_len + INITLEN) * sizeof(MD_Atom_Param));
      pb->atomp_len += INITLEN;
      break;
    case MD_BOND_PARAM:
      pb->bondp = realloc(pb->bondp, 
          (pb->bondp_len + INITLEN) * sizeof(MD_Bond_Param));
      pb->bondp_len += INITLEN;
      break;
    case MD_ANGLE_PARAM:
      pb->anglep = realloc(pb->anglep, 
          (pb->anglep_len + INITLEN) * sizeof(MD_Angle_Param));
      pb->anglep_len += INITLEN;
      break;
    case MD_DIHED_PARAM:
      pb->dihedp = realloc(pb->dihedp, 
          (pb->dihedp_len + INITLEN) * sizeof(MD_Dihed_Param));
      pb->dihedp_len += INITLEN;
      break;
    case MD_IMPR_PARAM:
      pb->imprp = realloc(pb->imprp, 
          (pb->imprp_len + INITLEN) * sizeof(MD_Impr_Param));
      pb->imprp_len += INITLEN;
      break;
    case MD_NBFIX_PARAM:
      pb->nbfixp = realloc(pb->nbfixp, 
          (pb->nbfixp_len + INITLEN) * sizeof(MD_Nbfix_Param));
      pb->nbfixp_len += INITLEN;
      break;
    default:
      fprintf(stderr, "unknown return value %d from param_read()\n", retval);
      exit(1);
    }
    retval = param_read(p);
  }
  if (retval < 0) {
    fprintf(stderr, "an error occurred in param_read()\n");
    exit(1);
  }

  /* print out the parameters from Pbuf */
  printf("\n\nList parameters from Pbuf:\n\n");
  printf("atomp_len = %d     atomp_num = %d\n",
         pb->atomp_len, pb->atomp_num);
  at = pb->atomp;
  for (i = 0;  i < pb->atomp_num;  i++) {
    printf("%3d  emin=%f  rmin=%f  emin14=%f  rmin14=%f  atomtype=%s\n", i,
           at[i].emin, at[i].rmin, at[i].emin14, at[i].rmin14, at[i].type);
  }
  printf("\n");
  printf("bondp_len = %d     bondp_num = %d\n",
         pb->bondp_len, pb->bondp_num);
  bn = pb->bondp;
  for (i = 0;  i < pb->bondp_num;  i++) {
    printf("%3d  k=%f  r0=%f\n", i, bn[i].k, bn[i].r0);
  }
  printf("\n");
  printf("anglep_len = %d     anglep_num = %d\n",
         pb->anglep_len, pb->anglep_num);
  an = pb->anglep;
  for (i = 0;  i < pb->anglep_num;  i++) {
    printf("%3d  k_theta=%f  theta0=%f  k_ub=%f  r_ub=%f\n", i,
           an[i].k_theta, an[i].theta0, an[i].k_ub, an[i].r_ub);
  }
  printf("\n");
  printf("dihedp_len = %d     dihedp_num = %d\n",
         pb->dihedp_len, pb->dihedp_num);
  di = pb->dihedp;
  for (i = 0;  i < pb->dihedp_num;  i++) {
    printf("%3d  mult=%d\n", i, di[i].mult);
    for (j = 0;  j < di[i].mult;  j++) {
      printf("     k=%f  phi=%f  n=%d\n", 
             di[i].k[j], di[i].phi[j], di[i].n[j]);
    }
  }
  printf("\n");
  printf("imprp_len = %d     imprp_num = %d\n",
         pb->imprp_len, pb->imprp_num);
  im = pb->imprp;
  for (i = 0;  i < pb->imprp_num;  i++) {
    printf("%3d  mult=%d\n", i, im[i].mult);
    for (j = 0;  j < im[i].mult;  j++) {
      printf("     k=%f  phi=%f  n=%d\n", 
             im[i].k[j], im[i].phi[j], im[i].n[j]);
    }
  }
  printf("\n");
  printf("nbfixp_len = %d     nbfixp_num = %d\n",
         pb->nbfixp_len, pb->nbfixp_num);
  nb = pb->nbfixp;
  for (i = 0;  i < pb->nbfixp_num;  i++) {
    printf("%3d  emin=%f  rmin=%f  emin14=%f  rmin14=%f  type1=%d  type2=%d\n",
           i, nb[i].emin, nb[i].rmin, nb[i].emin14, nb[i].rmin14,
           nb[i].param[0], nb[i].param[1]);
  }
  printf("\n");

  /* print out the parameters based on keys from Pdata */
/*     normally an MDIO user shouldn't access Pdata    */
/*     directly, but we do here for debugging purposes */
  printf("\n\nContents of parameter database:\n\n");
  kat = 0;
  kbn = 0;
  kan = 0;
  kdi = 0;
  kim = 0;
  knb = 0;
  for (i = 0;  i < p->pdata->key_num;  i++) {
    s = p->pdata->key[i];
    sscanf(s, "%s", stype);
    if (strcmp(stype, "atom") == 0) {
      sscanf(s, "%*s %s", a1);
      k = pdata_find_atom(p->pdata, a1);
      if (k < 0) {
        fprintf(stderr, "not in database: atom %s\n", a1);
        exit(1);
      }
      printf("%3d index: %2d %s %f %f %f %f %s\n", kat++, k, s, 
             at[k].emin, at[k].rmin, at[k].emin14, at[k].rmin14,
             at[k].type);
    }
    else if (strcmp(stype, "bond") == 0) {
      sscanf(s, "%*s %s %s", a1, a2);
      k = pdata_find_bond(p->pdata, a1, a2);
      if (k < 0) {
        fprintf(stderr, "not in database: bond %s %s\n", a1, a2);
        exit(1);
      }
      printf("%3d index: %2d %s %f %f\n", kbn++, k, s, 
             bn[k].k, bn[k].r0);
    }
    else if (strcmp(stype, "angle") == 0) {
      sscanf(s, "%*s %s %s %s", a1, a2, a3);
      k = pdata_find_angle(p->pdata, a1, a2, a3);
      if (k < 0) {
        fprintf(stderr, "not in database: angle %s %s %s\n", a1, a2, a3);
        exit(1);
      }
      printf("%3d index: %2d %s %f %f %f %f\n", kan++, k, s, 
             an[k].k_theta, an[k].theta0, an[k].k_ub, an[k].r_ub);
    }
    else if (strcmp(stype, "dihed") == 0) {
      sscanf(s, "%*s %s %s %s %s", a1, a2, a3, a4);
      k = pdata_find_dihed(p->pdata, a1, a2, a3, a4);
      if (k < 0) {
        fprintf(stderr,"not in database: dihed %s %s %s %s\n", a1, a2, a3, a4);
        exit(1);
      }
      printf("%3d index: %2d %s mult=%d\n", kdi++, k, s, di[k].mult);
      for (j = 0;  j < di[k].mult;  j++) {
        printf("           %f %f %d\n", di[k].k[j], di[k].phi[j], di[k].n[j]);
      }
    }
    else if (strcmp(stype, "impr") == 0) {
      sscanf(s, "%*s %s %s %s %s", a1, a2, a3, a4);
      k = pdata_find_impr(p->pdata, a1, a2, a3, a4);
      if (k < 0) {
        fprintf(stderr,"not in database: impr %s %s %s %s\n", a1, a2, a3, a4);
        exit(1);
      }
      printf("%3d index: %2d %s mult=%d\n", kdi++, k, s, im[k].mult);
      for (j = 0;  j < im[k].mult;  j++) {
        printf("           %f %f %d\n", im[k].k[j], im[k].phi[j], im[k].n[j]);
      }
    }
    else if (strcmp(stype, "nbfix") == 0) {
      sscanf(s, "%*s %s %s", a1, a2);
      k = pdata_find_nbfix(p->pdata, a1, a2);
      if (k < 0) {
        fprintf(stderr, "not in database: nbfix %s %s\n", a1, a2);
        exit(1);
      }
      printf("%3d index: %2d %s %f %f %f %f %d %d\n", knb++, k, s, 
             nb[k].emin, nb[k].rmin, nb[k].emin14, nb[k].rmin14, 
             nb[k].param[0], nb[k].param[1]);
    }
    else {
      fprintf(stderr, "unknown param type: _%s_ in string _%s_\n", stype, s);
      exit(1);
    }
  }

  return 0;
}