/*******************************************************************************
*                                                                              *
*                                   Viewmol                                    *
*                                                                              *
*                               T H E R M O . C                                *
*                                                                              *
*                 Copyright (c) Joerg-R. Hill, October 2003                    *
*                                                                              *
********************************************************************************
*
* $Id: thermo.c,v 1.3 2003/11/07 11:16:44 jrh Exp $
* $Log: thermo.c,v $
* Revision 1.3  2003/11/07 11:16:44  jrh
* Release 2.4
*
* Revision 1.2  2000/12/10 15:37:13  jrh
* Release 2.3
*
* Revision 1.1  1999/05/24 01:29:53  jrh
* Initial revision
*
*/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<Xm/Xm.h>
#include "viewmol.h"

#define AVOGADRO 6.022045e23

extern void balanceDialog(Widget, caddr_t, XmAnyCallbackStruct *);
extern int checkReaction(void);
extern char *getStringResource(Widget, char *);
extern void *getmem(size_t, size_t);
extern void fremem(void **);

void thermo(struct MOLECULE *);
int  reaction(int, int *);
double getMass(struct MOLECULE *);
double getVolume(struct MOLECULE *);
double getDensity(struct MOLECULE *mol);
int  isLinear(struct MOLECULE *);
double getSymmetryNumber(struct MOLECULE *);
void getSumFormula(struct MOLECULE *, char *, int);
void printCoefficients(char *, double *, char *, int, int);

extern Widget topShell;
extern struct WINDOW windows[];
extern struct MOLECULE *molecules;
extern double temp, pressure;
extern int nmolecule, ne, debug;

void thermo(struct MOLECULE *mol)
{
  const double N=6.022169e23, k=1.380622e-23;
  const double h=6.626196e-34, c=2.997925e10;
  double R=k*N, pi=4.0*atan(1.0), RT=R*temp;
  double sigma, hc, u=0.0, r, s;
  register int i;

  /* Translation */
  mol->htrl=1.5*RT;
  if (mol->unitcell)
    mol->pv=AVOGADRO*pressure*1e-30*getVolume(mol)*101325.0;
  else
    mol->pv=RT;
  mol->ctrl=2.5*R;
  if (temp > 0.0)
/*  mol->strl=R*(2.5+1.5*log(2.0*pi*k/(h*h))+log(k)+1.5*log(getMass(mol))+2.5*log(temp)-log(pressure));
    mol->strl=R*(2.5+1.5*log(2.0*pi*k*1.0e7/N)+log(k/(1.0e20*h*h*h))+1.5*log(getMass(mol))+2.5*log(temp)-log(pressure)); */
    mol->strl=R*(2.5*log(temp)+1.5*log(getMass(mol))-log(pressure))-2.31482*4.184;
  else
    mol->strl=0.0;

  /* Rotation */
  mol->hrot=0.0;
  mol->srot=0.0;
  mol->crot=0.0;
  mol->cvib=0.0;
  mol->hvib=0.0;
  mol->svib=0.0;
  if (mol->na > 1)
  {
    sigma=getSymmetryNumber(mol);
    if (isLinear(mol))
    {
      mol->hrot=RT;
      mol->crot=R;
      for (i=0; i<3; i++)
      {
        if (mol->rotConstants[i] != 0.0)
        {
          if (temp > 0.0)
            mol->srot=R*log(k*temp/(sigma*h*mol->rotConstants[i]))+R;
          else
            mol->srot=R;
          break;
        }
      }
    }
    else
    {
      mol->hrot=1.5*RT;
      mol->crot=1.5*R;
      if (temp > 0.0)
        mol->srot=0.5*R*log(pi*k*k*k*temp*temp*temp/(sigma*sigma*h*h*h
                 *c*c*c*mol->rotConstants[0]*mol->rotConstants[1]
                 *mol->rotConstants[2]))+1.5*R;
      else
        mol->srot=1.5*R;
    }

    /* Vibration */
    hc=h*c;
    for (i=0; i<mol->nmodes; i++)
    {
      if (mol->normal_modes[i].wavenumber < 10.0) continue;
      if (temp > 0.0)
      {
        u=hc*mol->normal_modes[i].wavenumber/(k*temp);
        r=exp(-u);
      }
      else
        r=0.0;
      s=1.0/(1.0-r);
      mol->cvib+=u*u*r*s*s;
      mol->hvib+=hc*mol->normal_modes[i].wavenumber*(r*s+0.5);
      mol->svib+=u*r*s-log(1.0/s);
    }
    mol->cvib*=R;
    mol->hvib*=N;
    mol->svib*=R;
  }
  if (debug)
  {
    printf("              C/[J/mol]    H/[J/mol]  S[J/(mol K)]\n");
    printf("Translation %12.6f %12.6f %12.6f\n", mol->ctrl, mol->htrl, mol->strl);
    printf("Rotation    %12.6f %12.6f %12.6f\n", mol->crot, mol->hrot, mol->srot);
    printf("Vibration   %12.6f %12.6f %12.6f\n", mol->cvib, mol->hvib, mol->svib);
  }
}

int reaction(int which, int *nreactions)
{
  struct MOLECULE *mol;
  double *matrix, *c, *denom, x;
  int nelem, nmol, found=FALSE, notInt, ret, type=0, *molIndex;
  char *elements;
  register int i, j, k;

  for (i=0; i<nmolecule; i++)
  {
    if (molecules[i].reaction != 0)
    {
      type=molecules[i].reaction;
      break;
    }
  }
  for (j=i+1; j<nmolecule; j++)
  {
    if (molecules[j].reaction != 0 && abs(type) != abs(molecules[j].reaction))
      return(INCONSISTENT_TYPE);
  }
  /* Construct stoichiometric matrix */
  matrix=(double *)getmem((size_t)(nmolecule*ne), sizeof(double));
  elements=(char *)getmem((size_t)ne, 4*sizeof(char));
  molIndex=(int *)getmem((size_t)nmolecule, sizeof(int));
  nelem=0;
  nmol=0;
  for (i=0; i<nmolecule; i++)
  {
    mol=&molecules[i];
    molIndex[i]=i;
    if (molecules[i].reaction == REACTANT)
    {
      mol->stoichioNumber=(-1);
      for (j=0; j<mol->na; j++)
      {
        found=FALSE;
        for (k=0; k<nelem; k++)
        {
          if (!strcmp(&elements[4*k], mol->atoms[j].element->symbol))
          {
            found=TRUE;
            matrix[nmol*ne+k]+=(double)(mol->reaction);
          }
        }
        if (!found)
        {
          strcpy(&elements[4*nelem], mol->atoms[j].element->symbol);
          matrix[nmol*ne+nelem]=(double)(mol->reaction);
          nelem++;
        }
      }
      nmol++;
    }
    else
      mol->stoichioNumber=0;
  }
  for (i=0; i<nmolecule; i++)
  {
    if (molecules[i].reaction == PRODUCT || molecules[i].reaction == ALLREACTIONS)
    {
      mol=&molecules[i];
      mol->stoichioNumber=1;
      for (j=0; j<mol->na; j++)
      {
        found=FALSE;
        for (k=0; k<nelem; k++)
        {
          if (!strcmp(&elements[4*k], mol->atoms[j].element->symbol))
          {
            found=TRUE;
            if (mol->reaction == PRODUCT)
              matrix[nmol*ne+k]+=(double)(mol->reaction);
            else
              matrix[nmol*ne+k]+=1.0;
          }
        }
        if (!found)
        {
          strcpy(&elements[4*nelem], mol->atoms[j].element->symbol);
          if (mol->reaction == PRODUCT)
            matrix[nmol*ne+nelem]=(double)(mol->reaction);
          else
            matrix[nmol*ne+nelem]=1.0;
          nelem++;
        }
      }
      nmol++;
    }
  }
  if (type != ALLREACTIONS)
  {
    for (i=0; i<nelem; i++)
    {
      found=FALSE;
      k=FALSE;
      for (j=0; j<nmol; j++)
      {
        if (matrix[j*ne+i] < 0) found=TRUE;
        if (matrix[j*ne+i] > 0) k=TRUE;
      }
      if (!found || !k)
      {
        fremem((void *)&elements);
        fremem((void *)&matrix);
          return(MISSING_ATOMS);
      }
    }
  }
  if (debug) printCoefficients("Start coefficients", matrix, elements, nelem,
                               FALSE);
  ret=REACTION_OK;
  for (i=0; i<nelem; i++)
  {
    while (matrix[i*ne+i] == 0.0 && nelem > i)
    {
      /* diagonal element is zero, try to find another row */
      /* where this is not the case, then swap the two rows */
      found=FALSE;
      for (j=i+1; j<nelem; j++)
      {
        if (matrix[i*ne+j] != 0.0)
        {
          found=TRUE;
          for (k=0; k<nmol; k++)
          {
            x=matrix[k*ne+j];
            matrix[k*ne+j]=matrix[k*ne+i];
            matrix[k*ne+i]=x;
          }
        }
        break;
      }
      if (!found)
      {
        /* there is no other row where the diagonal element would be non-zero */
        /* check whether the entire row is zero */
        found=FALSE;
        for (j=i+1; j<nmol; j++)
        {
          if (matrix[j*ne+i] != 0.0)
          {
            found=TRUE;
            break;
          }
        }
        if (!found)
        {
          /* this row is completely zero, eliminate it */
          for (j=i+1; j<nelem; j++)
          {
            for (k=0; k<nmol; k++)
            {
              x=matrix[k*ne+j];
              matrix[k*ne+j]=matrix[k*ne+i];
              matrix[k*ne+i]=x;
            }
          }
          nelem--;
        }
        else
        {
          /* there is no other row we can swap with and the current row is also */
          /* not completely zero, try to find a column we can swap with */
          found=FALSE;
          for (j=i+1; j<nmol; j++)
              {
            if (matrix[j*ne+i] != 0.0)
            {
              found=TRUE;
              k=molIndex[i];
              molIndex[i]=molIndex[j];
              molIndex[j]=k;
              for (k=0; k<nelem; k++)
              {
                x=matrix[j*ne+k];
                matrix[j*ne+k]=matrix[i*ne+k];
                matrix[i*ne+k]=x;
              }
              break;
            }
          }
          if (!found)
            ret=CANT_BALANCE;
          break;
        }
      }
    }
    if (ret || !found) break;
    x=1.0/matrix[i*ne+i];
    for (j=0; j<nmol; j++)
      matrix[j*ne+i]*=x;
    for (j=0; j<nelem; j++)
    {
      if (j == i) continue;
      for (k=nmol-1; k>=i; k--)
        matrix[k*ne+j]-=matrix[k*ne+i]*matrix[i*ne+j];
    }
/*  if (debug) printCoefficients("Iteration", matrix, elements, nelem,
                                     FALSE); */
  }
  c=(double *)getmem((size_t)nmol, sizeof(double));
  if (type == ALLREACTIONS)
  {
    *nreactions=nmol-nelem;
    *nreactions=*nreactions > 0 ? *nreactions : 1;
    for (j=0; j<nelem; j++)
      c[j]=-matrix[(which+nelem)*ne+j];
    for (j=nelem; j<nmol; j++)
    {
      if (which+nelem == j)
        c[j]=1.0;
      else
        c[j]=0.0;
    }
  }
  else
  {
    *nreactions=1;
    for (i=0; i<nmol; i++)
    {
      c[i]=0.0;
      for (j=nelem; j<nmol; j++)
      {
        c[i]+=matrix[j*ne+i];
      }
      if (c[i] == 0.0)
          c[i]=1.0;
      else
        c[i]=(-c[i]);
    }
  }
  notInt=TRUE;
  while (notInt)
  {
    notInt=FALSE;
    for (i=0; i<nmol; i++)
    {
      j=(int)rint(c[i]);
      if (fabs(c[i]-(double)j) > 1.0e-6)
      {
        notInt=TRUE;
        break;
      }
    }
    if (notInt)
    {
      denom=(double *)getmem((size_t)nmol, sizeof(double));
      for (i=0; i<nmol; i++)
      {
        denom[i]=0.0;
        j=1;
        while (TRUE)
        {
          x=c[i]*(double)j;
          k=(int)rint(x);
          if (fabs(x-(double)k) < 1.0e-6)
          {
            denom[i]=(double)j;
            break;
          }
          j++;
        }
        if (denom[i] > 1.0) break;
      }
      for (j=0; j<nmol; j++)
        c[j]*=denom[i];
      fremem((void **)&denom);
    }
  }
  /* Reaction(s) have been balanced, assign stoichiometric numbers */
  j=0;
  if (type == ALLREACTIONS)
  {
    for (i=0; i<nmolecule; i++)
    {
      if (molecules[i].reaction == ALLREACTIONS)
        molecules[i].stoichioNumber=(int)rint(c[molIndex[j++]]);
    }
  }
  else
  {
    /* Assign stoichiometric numbers for reactants first, for products second, since */
    /* reactants are in the stoichiometric matrix first, followed by products */
    for (i=0; i<nmolecule; i++)
    {
      if (molecules[i].reaction == REACTANT)
        molecules[i].stoichioNumber=(int)rint(-c[molIndex[j++]]);
    }
    for (i=0; i<nmolecule; i++)
    {
      if (molecules[i].reaction == PRODUCT)
        molecules[i].stoichioNumber=(int)rint(c[molIndex[j++]]);
    }
    /* Check whether balancing screwed up the original, user desired, assignment of */
    /* reactants and products, if it did generate error */
    for (i=0; i<nmolecule; i++)
    {
      if ((molecules[i].reaction == REACTANT && molecules[i].stoichioNumber > 0) ||
          (molecules[i].reaction == PRODUCT  && molecules[i].stoichioNumber < 0))
      {
        ret=CANT_BALANCE;
        break;
      }
    }
  }
  if (debug) printCoefficients("Resulting coefficients", matrix, elements,
                               nelem, TRUE);
  fremem((void **)&c);
  fremem((void **)&molIndex);
  fremem((void *)&elements);
  fremem((void *)&matrix);
/*if (ret == CANT_BALANCE && type != ALLREACTIONS) */
  if (type != ALLREACTIONS && !checkReaction())
  {
    balanceDialog((Widget)0, (caddr_t)0, (XmAnyCallbackStruct *)0);
    ret=REACTION_OK;
  }
  return(ret);
}

int isLinear(struct MOLECULE *mol)
{
  if (mol->rotConstants[0] == 0.0 || mol->rotConstants[1] == 0.0 ||
      mol->rotConstants[2] == 0.0)
    return(TRUE);
  else
    return(FALSE);
}

double getSymmetryNumber(struct MOLECULE *mol)
{
  if (!strcmp(mol->pgroup, "c1") || !strcmp(mol->pgroup, "cs") ||
      !strcmp(mol->pgroup, "ci") || !strcmp(mol->pgroup, "c*v"))
    return(1.0);
  if (!strcmp(mol->pgroup, "d*h")) return(2.0);
  if (!strcmp(mol->pgroup, "td")) return(12.0);
  if (!strcmp(mol->pgroup, "oh")) return(24.0);
  if (!strcmp(mol->pgroup, "ih")) return(60.0);
  if (mol->pgroup[0] == 'c') return(atof(&mol->pgroup[1]));
  if (mol->pgroup[0] == 'd') return(2*atof(&mol->pgroup[1]));
  if (mol->pgroup[0] == 's') return(atof(&mol->pgroup[1])/2.0);
  return(1);
}

double getMass(struct MOLECULE *mol)
{
  double mass=0.0;
  register int i;

  for (i=0; i<mol->na; i++)
  {
    if (mol->atoms[i].flags & ORIGINAL)
      mass+=mol->atoms[i].mass;
  }

  return(mass);
}

double getVolume(struct MOLECULE *mol)
{
  double volume=0.0, x, y, z;
  int n=1;

  if (mol->unitcell)
  {
    x=mol->unitcell->corners[n-1].x;
    y=mol->unitcell->corners[n-1].y;
    z=mol->unitcell->corners[n-1].z;
    volume=(mol->unitcell->corners[n].x-x)*(mol->unitcell->corners[n+1].y-y)*(mol->unitcell->corners[n+2].z-z)
          +(mol->unitcell->corners[n].y-y)*(mol->unitcell->corners[n+1].z-z)*(mol->unitcell->corners[n+2].x-x)
          +(mol->unitcell->corners[n].z-z)*(mol->unitcell->corners[n+1].x-x)*(mol->unitcell->corners[n+2].y-y)
          -(mol->unitcell->corners[n].z-z)*(mol->unitcell->corners[n+1].y-y)*(mol->unitcell->corners[n+2].x-x)
          -(mol->unitcell->corners[n].x-x)*(mol->unitcell->corners[n+1].z-z)*(mol->unitcell->corners[n+2].y-y)
          -(mol->unitcell->corners[n].y-y)*(mol->unitcell->corners[n+1].x-x)*(mol->unitcell->corners[n+2].z-z);
  }
  return(volume);
}

double getDensity(struct MOLECULE *mol)
{
  double volume, mass, density=0.0;

  if ((volume=getVolume(mol)) != 0.0)
  {
    mass=getMass(mol);
    density=mass/(AVOGADRO*volume*1e-24);
  }
  return(density);
}

void getSumFormula(struct MOLECULE *mol, char *formula, int sub)
{
  struct ELEMENT **found;
  int *count, nfound, new;
  char num[10], order[MAXLENLINE], *word;
  register int i, j;

  found=(struct ELEMENT **)getmem((size_t)mol->na, sizeof(struct ELEMENT *));
  count=(int *)getmem((size_t)mol->na, sizeof(int));
  nfound=0;
  for (i=0; i<mol->na; i++)
  {
    new=TRUE;
    for (j=0; j<nfound; j++)
    {
      if (mol->atoms[i].element == found[j])
      {
        count[j]++;
        new=FALSE;
        break;
      }
    }
    if (new && strncmp(mol->atoms[i].element->symbol, "Uc", 2) && strncmp(mol->atoms[i].element->symbol, "Mp", 2))
    {
      found[nfound]=mol->atoms[i].element;
      count[nfound]++;
      nfound++;
    }
  }
  strncpy(order, getStringResource(topShell, "elementSortOrder"), MAXLENLINE-1);
  word=strtok(order, ",");
  *formula='\0';
  while (word != NULL)
  {
    for (i=0; i<nfound; i++)
    {
      if (!strcmp(found[i]->symbol, word))
      {
        strcat(formula, found[i]->symbol);
        if (count[i] > 1)
        {
	  if (sub)
	  {
            if (count[i] > 9)
              sprintf(num, "_{%d}", count[i]);
            else
              sprintf(num, "_%d", count[i]);
	  }
	  else
	    sprintf(num, "%d", count[i]);
          strcat(formula, num);
        }
        count[i]=0;
      }
    }
    word=strtok(NULL, ",");
  }
  for (i=0; i<nfound; i++)
  {
    if (count[i] > 0)
    {
      strcat(formula, found[i]->symbol);
      if (count[i] > 1)
      {
	if (sub)
	{
          if (count[i] > 9)
            sprintf(num, "_{%d}", count[i]);
          else
            sprintf(num, "_%d", count[i]);
	}
	else
          sprintf(num, "%d", count[i]);
        strcat(formula, num);
      }
    }
  }
  fremem((void **)&found);
  fremem((void **)&count);
}

void printCoefficients(char *label, double *matrix, char *elements, int nelem,
                       int final)
{
  char formula[MAXLENLINE];
  register int i, j, k;

  printf("%s\n\t", label);
  for (i=0; i<nelem; i++)
    printf("  %s  ", &elements[4*i]);
  printf("\n");
  k=0;
  for (i=0; i<nmolecule; i++)
  {
    if (molecules[i].reaction == REACTANT || molecules[i].reaction == ALLREACTIONS)
    {
      getSumFormula(&molecules[i], formula, TRUE);
      printf("%s\t", formula);
      for (j=0; j<nelem; j++)
        printf("%5.2f ", matrix[k*ne+j]);
      if (final) printf(" %d", molecules[i].stoichioNumber);
      printf("\n");
      k++;
    }
  }
  for (i=0; i<nmolecule; i++)
  {
    if (molecules[i].reaction == PRODUCT)
    {
      getSumFormula(&molecules[i], formula, TRUE);
      printf("%s\t", formula);
      for (j=0; j<nelem; j++)
        printf("%5.2f ", matrix[k*ne+j]);
      if (final) printf(" %d", molecules[i].stoichioNumber);
      printf("\n");
      k++;
    }
  }
}