/*  _______________________________________________________________________
 *
 *                        RDPARM/PTRAJ: APR 2000
 *  _______________________________________________________________________
 *
 *  $Header: /thr/gamow/cvsroot/amber7/src/ptraj/energy.c,v 1.1 2001/09/18 20:20:57 cheatham Exp $
 *
 *  Revision: $Revision: 1.1 $
 *  Date: $Date: 2001/09/18 20:20:57 $
 *  Last checked in by $Author: cheatham $
 *
 *  This source is now archived under CVS at Scripps in the amber7 tree.
 *
 *  NOTE: this is a beta, pre-release version, is under constant development,
 *  probably contains some bugs and is under constant revision; therefore
 *  take care.  Please report bugs (and fixes!) to cheatham@chpc.utah.edu or 
 *  cheatham@cgl.ucsf.edu
 *
 *  Do not distribute this code without explicit permission.
 *  Do not incorporate into other codes without explicit permission.
 *
 *  Current contact information:
 *
 *     Thomas Cheatham, III
 *     2000 East, 30 South, Skaggs Hall 201
 *     Department of Medicinal Chemistry
 *     University of Utah
 *     Salt Lake City, UT 84112-5820
 *     cheatham@chpc.utah.edu
 *     (801) 587-9653
 *     FAX: (801) 585-5366
 *
 *  Other contributors:
 *
 *    David Case (Scripps)
 *    Michael Crowley (Scripps)
 *    Jed Pitera (UCSF)
 *    Vickie Tsui (Scripps)
 *  _______________________________________________________________________
 */



#include <stdio.h>
#include <string.h>
#include <math.h>

#define ENERGY_MODULE
#include "ptraj.h"


   void
initializeAtomInfo(atomInfo *atom)
{
  atom->energy = 0.0;
  atom->force = 0.0;
  atom->bonds = 0;
  atom->angles = 0;
  atom->dihedrals = 0;
}



   rtfInfo *
loadEnergyInfoFromPrmtop(char *filename)
{
  FILE *fp;
  char *buffer;
  Parm *oldparm, *newparm;
  rtfInfo *rtf;

  int i, ai, aj, ak, al;

  /*
   *  create a link to the current "global" parm
   */
  oldparm = parm;

  /*
   *  open up the file and set up a new parm
   */
  if ( openFile(&fp, filename, "r") == 0 )
    error("loadEnergyInfoFromPrmtop()", "Cannot open prmtop file %s", filename);
  newparm = (Parm *) safe_malloc(sizeof(Parm));
  newparm->fp = fp;
  newparm->filename = (char *) safe_malloc(sizeof(char) * (strlen(filename)+1));
  strcpy(newparm->filename, filename);
  parm = newparm;
  readParm();
  parm = oldparm;

  /*
   *  allocate rtfInfo *
   */
  rtf = (rtfInfo *) safe_malloc(sizeof(rtfInfo));
  INITIALIZE_rtfInfo(rtf);

  /*
   *  setup atom info
   */

  rtf->atoms = parm->NTOTAT;
  rtf->atom = (atomInfo *) safe_malloc(sizeof(atomInfo) * rtf->atoms);
  for (i=0; i<rtf->atoms; i++)
    initializeAtomInfo(&rtf->atom[i]);

  /*
   *  process bonds
   *    for each bond:
   *       copy over the parameters
   *       update the atom list which points back to the potential bonds
   *       zero out energy/force
   */

  rtf->bonds = newparm->NBONH + newparm->MBONA;
  rtf->bond = (Bond *) safe_malloc(sizeof(Bond) * rtf->bonds);
  for (i = 0; i < rtf->bonds; i++) {

    rtf->bond[i].rk = newparm->bond[i].rk;
    rtf->bond[i].req = newparm->bond[i].req;
    ai = newparm->bond[i].atom[0];
    rtf->bond[i].atom[0] = ai;
    aj = newparm->bond[i].atom[1];
    rtf->bond[i].atom[1] = aj;
    rtf->bond[i].scale = newparm->bond[i].scale;
    rtf->bond[i].active = 1;

    if ( rtf->atom[ai].bonds == MAX_BONDS_PER_ATOM ) {
      fprintf(stderr, "WARNING in loadEnergyInfoFromPrmtop().  Too many bonds per\n");
      fprintf(stderr, "atom were found for atom %i.  Max is %i\n.  Increase\n",
	      ai, MAX_BONDS_PER_ATOM);
      fprintf(stderr, "MAX_BONDS_PER_ATOM in energy.h and recompile\n");
    }
    if ( rtf->atom[aj].bonds == MAX_BONDS_PER_ATOM ) {
      fprintf(stderr, "WARNING in loadEnergyInfoFromPrmtop().  Too many bonds per\n");
      fprintf(stderr, "atom were found for atom %i.  Max is %i\n.  Increase\n",
	      aj, MAX_BONDS_PER_ATOM);
      fprintf(stderr, "MAX_BONDS_PER_ATOM in energy.h and recompile\n");
    }
	
    rtf->atom[ai].bondP[ rtf->atom[ai].bonds++ ] = i;
    rtf->atom[aj].bondP[ rtf->atom[aj].bonds++ ] = i;

    rtf->Ebond[i] = 0.0;

  }  

  return(rtf);


}

   rtfInfo *
loadEnergyInfoFromPSF()
{
  return(NULL);
}



   double
calculateBondEnergy(rtfInfo *rtf, double *x, double *y, double *z) 
{
  double rx, ry, rz, r2, r, delr;
  double totalE;

  int i, ai, aj;

  totalE = 0;
  for (i=0; i < rtf->bonds; i++) {

    if (rtf->bond[i].active) {

      ai = rtf->bond[i].atom[0];
      aj = rtf->bond[i].atom[1];

      rx = x[ai] - x[aj];
      ry = y[ai] - y[aj];
      rz = z[ai] - z[aj];

      r2 = rx*rx + ry*ry + rz*rz;
      r = sqrt(r2);
      delr = r - rtf->bond[i].req;

      rtf->Ebond[i] = rtf->bond[i].rk * delr * delr;
      totalE += rtf->Ebond[i];

      /*
      df = rtf->bond[i].rk * delr * 2.0 / r;
      fxi =  rx * df;
      fxj = -rx * df;
      fyi =  ry * df;
      fyj = -ry * df;
      fzi =  rz * df;
      fzj = -rz * df;
      */
    }
  }
  
  return (totalE);
}