/***** 
This file is part of the Babel Program
Copyright (C) 1992-96 W. Patrick Walters and Matthew T. Stahl 
All Rights Reserved 
All Rights Reserved 
All Rights Reserved 
All Rights Reserved 

For more information please contact :

babel@mercury.aichem.arizona.edu
--------------------------------------------------------------------------------

FILE : wrgau.c
AUTHOR(S) : Pat Walters, Matthew Stahl
DATE : 1-93
PURPOSE : Routines to write a Gaussian input file

******/


#include "bbltyp.h"

#define DELIMSTR " \n\t,"


int 
write_gaussian(FILE *file1, ums_type *mol)
{
  int i=0;
  char type_name[5];
  double a,d;
  int result;
  FILE *file2;
  char buffer[BUFF_SIZE];
  char babel_dir[BUFF_SIZE];

  if (Atoms > 0)
  {
    initialize_internal(&mol);
    cartint(mol);  
    cartgeom(mol);
  }
  
  if ((file2 = fopen("gauss.hdr","r")) != NULL)
  {
    while (fgets(buffer,sizeof(buffer),file2))
      fprintf(file1,"%s",buffer);

    if (file2)
      fclose(file2);
  }
  else 
    if (getenv("BABEL_DIR"))
    {
      strcpy(babel_dir,getenv("BABEL_DIR"));
      strcat(babel_dir,"/gauss.hdr");

      if ((file2 = fopen(babel_dir,"r")) != NULL)
      {
	while (fgets(buffer,sizeof(buffer),file2))
	  fprintf(file1,"%s",buffer);

	if (file2)
	  fclose(file2);
      }
    }
    else
    {
      fprintf(file1,"%cmem=3700000 words\n",'\045');
      fprintf(file1,"#%s\n\n",OutputKeywords);
      fprintf(file1,"This Gaussian input file generated by Babel %s\n\n",BABEL_VERSION);
      fprintf(file1,"Ch Mu\n");
    }

  for(i=1; i <= Atoms; i++)
  {
    strcpy(type_name,Type(i));
    result = get_output_type(i,"XYZ",Type(i),type_name,all_caps);
    clean_atom_type(type_name);
    switch(i)
    {
    case 1 :
      fprintf(file1,"  %s \n",type_name);
      break;
    case 2 :
      fprintf(file1,"  %s    %d r%d \n",
	      type_name,
	      mol->internal[i].na,i);
      break;
    case 3 :
      fprintf(file1,"  %s    %d r%d    %d a%d\n",
	      type_name,
	      mol->internal[i].na,i,
	      mol->internal[i].nb,i);
      break;
    default :
      fprintf(file1,"  %s    %d r%d    %d a%d    %d d%d\n",
	      type_name,
	      mol->internal[i].na,i,
	      mol->internal[i].nb,i,
	      mol->internal[i].nc,i); 
    }
  }
  
  fprintf(file1,"Variables:\n");

  for(i=2; i <= Atoms; i++)
  {
    a = fix_g90_angle(mol->internal[i].w);
    d = fix_g90_angle(mol->internal[i].t);
    switch(i)
    {
    case 2 :
      fprintf(file1,"r2= %6.4f\n",mol->internal[2].r);
      break;
    case 3 :
      fprintf(file1,"r3= %6.4f\na3= %6.2f\n",
	      mol->internal[3].r, a);
      break;
    default :
	fprintf(file1,"r%d= %6.4f\na%d= %6.2f\nd%d= %6.2f\n",
		i,mol->internal[i].r,
		i,a,
		i,d);
      break;
    }
  }
  fprintf(file1,"\n");

  return(TRUE);
}

double 
fix_g90_angle(double angle)
{
  if (angle < 0) 
    return(360 + angle);
  return(angle);
}


int
  write_gaussian_template(FILE *file1,ums_type *mol)
{
  int i=0;
  FILE *file2;
  char buffer[BUFF_SIZE];

  if (NOTEQ("gauout",InputInfo.code))
  {
    show_warning("Template Z-matrix may only be used if input file is gaussian output type");
    return(FALSE);
  }
  
  if (Atoms > 0)
    initialize_internal(&mol);

  file2 = open_read(InfileName);
  if (!file2)
  {
    show_warning("Unable to open file");
    return(FALSE);
  }
  while (fgets(buffer,sizeof(buffer),file2))
  {
    if (strstr(buffer,"Z-matrix") || strstr(buffer,"Z-Matrix"))
    {
      if (strstr(buffer,"Symbolic"))
      {
	for (i = 0;i < 2;i++)
	  fgets(buffer,sizeof(buffer),file2);
	break;
      }
      else
	if (strstr(buffer,"checkpoint"))
	{
	  for (i = 0;i < 3;i++)
	    fgets(buffer,sizeof(buffer),file2);
	  break;
	}
      else
      {
	show_warning("Unable to locate start point of template z-matrix");
        if (file2)
          fclose(file2);
        return(FALSE);
      }
    }
  }
  
  if (check_for_eof(file2))
  {
    show_warning("Unable to find template z-matrix");
    if (file2)
      fclose(file2);
    return(FALSE);
  }


  NA(1) = NB(1) = NC(1) = 0;

  fgets(buffer,sizeof(buffer),file2);  
  NA(2) = atoi(gettoken(buffer,DELIMSTR,2));
  NB(2) = NC(2) = 0;
  
  fgets(buffer,sizeof(buffer),file2);
  NA(3) = atoi(gettoken(buffer,DELIMSTR,2));
  NB(3) = atoi(gettoken(buffer,DELIMSTR,4));
  NC(3) = 0;
  

  for (i = 4;i <= Atoms;i++)
  {
    fgets(buffer,sizeof(buffer),file2);
    NA(i) = atoi(gettoken(buffer,DELIMSTR,2));
    NB(i) = atoi(gettoken(buffer,DELIMSTR,4));
    NC(i) = atoi(gettoken(buffer,DELIMSTR,6));
  }

  if (file2)
    fclose(file2);
  
  tmplt_cartgeom(mol);
  write_gaussian(file1,mol);

  return(TRUE);
}


void 
  tmplt_cartgeom(ums_type *mol)
{  
  int i;

  R(1) = 0.0;
  W(1) = 0.0;
  T(1) = 0.0;
  
  R(2) = distance(Point(NA(2)),Point(2));
  W(2) = 0.0;
  T(2) = 0.0;

  R(3) = distance(Point(mol->internal[3].na),Point(3));
  W(3) = bond_angle(Point(3),Point(NA(3)),Point(NB(3)));
  T(3) = 0.0;    

  for (i = 4;i <= Atoms;i++)
  {
    R(i) = distance(Point(i),Point(NA(i)));
    W(i) = bond_angle(Point(i),Point(NA(i)),Point(NB(i)));
    T(i) = torsion(Point(i),Point(NA(i)),Point(NB(i)),Point(NC(i)));
    T(i) = fix_g90_angle(T(i));
  }

}