# ifndef _WATER_H_
# define _WATER_H_

# define MAX_HOH 1200
typedef vector9 Water;
typedef double *Water_ptr;
typedef Water  Water_array[MAX_HOH];

# define HTOD(hoh) {prt_hoh(hoh,"hoh");}

/*       --                                                            */
# define OX(hoh) (hoh)
# define O (hoh) (hoh)
# define H1(hoh) ((hoh)+3)
# define H2(hoh) ((hoh)+6)
/*       --                                                             */

# ifdef OFF
# define HOHCPY(hoh1,hoh2) bcopy((char *) hoh1, (char *) hoh2, sizeof (Water))
# endif
# define HOHCPY(hoh1,hoh2) {    \
   VECT3CPY(OX(hoh1),OX(hoh2)); \
   VECT3CPY(H1(hoh1),H1(hoh2)); \
   VECT3CPY(H2(hoh1),H2(hoh2)); \
  }

# define HOHSUBTRACT(nhoh,ohoh,v)\
{\
  OX(nhoh)[0]=OX(ohoh)[0] - v[0];\
  OX(nhoh)[1]=OX(ohoh)[1] - v[1];\
  OX(nhoh)[2]=OX(ohoh)[2] - v[2];\
  H1(nhoh)[0]=H1(ohoh)[0] - v[0];\
  H1(nhoh)[1]=H1(ohoh)[1] - v[1];\
  H1(nhoh)[2]=H1(ohoh)[2] - v[2];\
  H2(nhoh)[0]=H2(ohoh)[0] - v[0];\
  H2(nhoh)[1]=H2(ohoh)[1] - v[1];\
  H2(nhoh)[2]=H2(ohoh)[2] - v[2];\
}

/**
***  ---> Minimum-image Convention Macros <---
**/
 
# define CROUND(x) ((x) > 0 ? (int) ( (x) + 0.5 ) : - ( (int) ( 0.5 - (x) ) ) ) 
# define CEILFLOOR(x) ((x) > 0 ? floor(x) : ceil(x) ) 

# define MK_NBOX(nbox,rij,CELL) { \
    nbox[0]=CROUND(rij[0]/CELL[0]); \
    nbox[1]=CROUND(rij[1]/CELL[1]); \
    nbox[2]=CROUND(rij[2]/CELL[2]); \
  } 

# ifdef OFF
# define MINIM(rij,nbox,cell) (rij - ((double) nbox)*cell) 
# endif
# define MINIM(rij,nbox,cell) (nbox ? (rij - ((double) nbox)*cell) : rij)
# define MINIM3(mij,rij,nbox,CELL)  \
  { \
    mij[0] = MINIM(rij[0],nbox[0],CELL[0]); \
    mij[1] = MINIM(rij[1],nbox[1],CELL[1]); \
    mij[2] = MINIM(rij[2],nbox[2],CELL[2]); \
  }

/*       ---------                                            */
# define MINIM_DSQ(rij,nbox,CELL,temp) ( \
/*       ---------                                            */\
   SQR(MINIM(rij[0],nbox[0],CELL[0]),temp[0]) + \
   SQR(MINIM(rij[1],nbox[1],CELL[1]),temp[1]) + \
   SQR(MINIM(rij[2],nbox[2],CELL[2]),temp[2])   \
  )

/*       --------                                            */
# define HOHMINIM(hoh,CELL) { \
/*       --------                                            */\
    int nbox[3];              \
    Water mhoh;               \
    MK_NBOX(nbox,hoh,CELL);   \
    MINIM3(OX(mhoh),OX(hoh),nbox,CELL); \
    MINIM3(H1(mhoh),H1(hoh),nbox,CELL); \
    MINIM3(H2(mhoh),H2(hoh),nbox,CELL); \
    HOHCPY(mhoh,hoh);                   \
  }  

/*       -------                                            */
# define FOLDHOH(hoh,mhoh,CELL)                                  \
/*       -------                                            */   \
{                                                                \
  int nbox[3];                                                   \
  MK_NBOX(nbox,OX(hoh),CELL);                                    \
  MINIM3(OX(mhoh),OX(hoh),nbox,CELL);                            \
  MK_NBOX(nbox,H1(hoh),CELL);                                    \
  MINIM3(H1(mhoh),H1(hoh),nbox,CELL);                            \
  MK_NBOX(nbox,H2(hoh),CELL);                                    \
  MINIM3(H2(mhoh),H2(hoh),nbox,CELL);                            \
}  

/**    For water 
***    R(OH) 0.9572 /,
***    theta TORAD(104.52 ),
**/
# define MAG_HOHP .58588228 /* A R(OH)*cos(theta) */

/***
 *** Using these definitions:
 ***    1. The dipole moment p is defined from the midpoint of the
 ***         hydrogens to oxygen.
 ***    2. The radius vector from the origin to the oxygen. 
 ***    3. Therefore,
 ***       nmu[0] = radial projection of dipole moment 
 ***              > 0 when hydrogens face inwards.
 ***
 **/


/*       ------                                           */
# define CALC_P(hohi,p) \
/*       ------                                          */\
{                                                         \
  p[0]=hohi[0] - (H1(hohi)[0]+H2(hohi)[0])/2.0;           \
  p[1]=hohi[1] - (H1(hohi)[1]+H2(hohi)[1])/2.0;           \
  p[2]=hohi[2] - (H1(hohi)[2]+H2(hohi)[2])/2.0;           \
}

/*       -------                                           */
# define NCALC_P(hohi,p) \
/*       -------                                          */\
{                                                         \
  double magp;\
  p[0]=hohi[0] - (H1(hohi)[0]+H2(hohi)[0])/2.0;           \
  p[1]=hohi[1] - (H1(hohi)[1]+H2(hohi)[1])/2.0;           \
  p[2]=hohi[2] - (H1(hohi)[2]+H2(hohi)[2])/2.0;           \
  magp = sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);\
  p[0] /= magp;\
  p[1] /= magp;\
  p[2] /= magp;\
}

/*       --------                                  */
# define CALC_RPP(mhohi,p,rpp,orig)            \
/*       --------                                  */\
{                                                 \
  Vector dir;                                   \
  double magh;                                    \
  dir[0]=OX(mhohi)[0]-orig[0];                         \
  dir[1]=OX(mhohi)[1]-orig[1];                         \
  dir[2]=OX(hohi);\
  magh = hypot(dir[0],dir[1]);                    \
  rpp = p[0]*dir[0] + p[1]*dir[1];                \
/*  if (magh == 0)  */                                \
/*    rpp = 0;      */                                \
/*  else            */                                \
/*    rpp /= magh*MAG_HOHP; */                         \
  rpp /= magh;\
}

/*       ---------                                  */
# define CALC_COST(mhohi,p,pr)            \
/*       ---------                                  */\
{                                                 \
  double magh = hypot(OX(mhohi)[0],OX(mhohi)[1]);\
  /* double magp = hypot(p[0],p[1]); */\
  /* double magp = MAG_HOHP; */\
  double magp = sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);\
  if (magh == 0 || magp == 0)\
     pr = 0;\
  else {\
    pr = p[0]*OX(mhohi)[0] + p[1]*OX(mhohi)[1];\
    pr /= magh*magp;\
  }\
}

/*       -------                                  */
# define CALC_PR(mhohi,p,pr)            \
/*       -------                                  */\
{                                                 \
  double magh = hypot(OX(mhohi)[0],OX(mhohi)[1]);                    \
  pr = p[0]*OX(mhohi)[0] + p[1]*OX(mhohi)[1];                \
  if (magh == 0)                                  \
     pr = 0;                                      \
  else                                            \
     pr /= magh;                         \
}


# endif