/*****************************************************************************
 *                                                                           *
 *               This script reads in a CPMD 3.4.1 output file               *
 *               and calculates the angle    between specified               *
 *               atoms.  It does not check  whether  there are               *
 *               enough atoms in the system   (e.g. requesting               *
 *               angle between atoms 1, 2 and n+1 in an n-atom               *
 *                 system WILL give SOME number), so beware.                 *
 *                                                                           *
 *                           Veronika Brazdova 2002                          *
 *                           vb@chemie.hu-berlin.de                          *
 *                                                                           *
 *****************************************************************************/

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

#define PI 3.141592654
#define BUFFER 150

/***********************************main**************************************/

/// Return codes for function main:
//  0 successful run
//  1 no command line arguments or file not found
//  2 two or three atoms of the same number

int main(int argc, char *argv[])
{
  int i, j, sequence[3], at[3], save_at;
  double atom[3][3], save;
  double vector[2][3], amplitude[2]; 
  double scalar, cosine;
  int atom_no[3], switched;
  double angle;
  char   line[BUFFER];
  FILE   *fp;

  // argv[0]=name-of-the-executable
  // argv[1]=input file
  // argv[2]=number of the first atom
  // argv[3]=number of the second atom
  // argv[4]=number of the third atom
  
  /// Handle the command line arguments
  if (argc < 5)
  {
    fprintf(stderr,\
      "angle CPMD_3.4.1_output no_1st_atom no_2nd_atom no_3rd_atom\n");
   return 1;
  }
  
  for (i = 0; i < 3; i++)
  {
    atom_no[i] = atoi(argv[i+2]);
    sequence[i] = atom_no[i];
  }

  if (atom_no[0] == atom_no[1]\
    || atom_no[0] == atom_no[2] || atom_no[1] == atom_no[2])
  {
    fprintf(stderr,"angle: two or three atoms of the same number\n");
    return 2;
  }

  switched = 1;
  for (i = 0; i < 3; i++)
    at[i] = i;

  /// It is necessary to sort the atoms in order of their numbers
  //  in the CPMD output (into the array sequence), but it is also 
  //  necessary to keep the original order (in the array atom_no), 
  //  so as to obtain the correct angle, as well as to have array
  //  indices 0 -- 2 switched so that the elements atom[at[0]][i]
  //  point at coordinates of the first atom which was 
  //  _requested_from_the_command_line_, the elements atom[at[1]][i]
  //  at coordinates of the second atom, etc.

  while (switched)
  {
    switched = 0;
    for (i = 0; i < 3; i++)
    {
      if ((sequence[i+1] - sequence[i]) < 0)
      {
        save = sequence[i];
        sequence[i] = sequence[i+1];
        sequence[i+1] = save;
        save_at = at[i];
        at[i] = at[i+1];
        at[i+1] = save_at;
        switched = 1;
      }
    }
  }

  ///Print out information about the input parameters

  printf("%s \n","# Angle: changes in angles during a CPMD run.");
  printf("%s %i %i %i\n","# Atoms ",atom_no[0],atom_no[1],atom_no[2]);
  printf("%s %s\n","# CPMD output file: ",argv[1]);
      
  if (!(fp = fopen(argv[1],"r")))
  {
    fprintf(stderr,"angle: cannot open file %s\n",argv[1]);
    return 1;
  } 

  i = 0;

  while (fgets(line,BUFFER,fp))
  {
    if (strstr(line,"   ATOM          COORDINATES        "))
    {
      for (j = 0; j < sequence[0]; j++)
        fgets(line,BUFFER,fp);
      sscanf(line," %*i %*s %lf %lf %lf %*s",&atom[at[0]][0],\
        &atom[at[0]][1],&atom[at[0]][2]);
      for (j = sequence[0]; j < sequence[1]; j++)
        fgets(line,BUFFER,fp);
      sscanf(line," %*i %*s %lf %lf %lf %*s",&atom[at[1]][0],\
        &atom[at[1]][1],&atom[at[1]][2]);
      for (j = sequence[1]; j < sequence[2]; j++)
        fgets(line,BUFFER,fp);
      sscanf(line," %*i %*s %lf %lf %lf %*s",&atom[at[2]][0],\
        &atom[at[2]][1],&atom[at[2]][2]);
      scalar = 0; amplitude[0] = 0; amplitude[1] = 0;
      for (j = 0; j < 3; j++)
      {
        vector[0][j] = atom[0][j] - atom[1][j];
        vector[1][j] = atom[2][j] - atom[1][j];
        scalar = scalar + vector[0][j] * vector[1][j];
        amplitude[0] = amplitude[0] + vector[0][j] * vector[0][j];
        amplitude[1] = amplitude[1] + vector[1][j] * vector[1][j];
      }
      amplitude[0] = sqrt(amplitude[0]);
      amplitude[1] = sqrt(amplitude[1]); 
      cosine = scalar / (amplitude[0] * amplitude[1]);
      angle = acos(cosine) * 180 / PI; 
      printf("%i %lf \n",i,angle); 
      i = i + 1;
    }
  }
  return 0;
} /// main