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

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

#define PI 3.141592654
#define BUFFER 150
#define BOHR2ANG 0.5291772

/// Return codes for function main:
//  0 successful run
//  1 no command line arguments or file not found

int main(int argc, char *argv[])
{
  int i, j, save_val, atom_no[2];
  double atom1[3], atom2[3], temp[3], distance, scale_factor;
  char   line[BUFFER];
  FILE   *fp;

  save_val = 0;

  // 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] = scale the distance on/off
  
  /// Handle the command line arguments
  if (argc < 4)
   {
     fprintf(stderr,\
       "distance cpmd_output no_1st_atom no_2nd_atom [scale]\n");
     return 1;
   }

  /// Scale the bondlength? Default is no.
  scale_factor = 1;
  if (argc > 4)
    scale_factor = -1;
  
  atom_no[0] = atoi(argv[2]);
  atom_no[1] = atoi(argv[3]);

  if ((atom_no[1] - atom_no[0]) < 0)
  {
    save_val = atom_no[0];
    atom_no[0] = atom_no[1];
    atom_no[1] = save_val;
  }

  /// Print out information about the input parameters
  printf("%s \n","# Distance: changes in bondlengths during a CPMD run.");
  printf("%s %i %s %i\n","# Atoms ",atom_no[0]," and ",atom_no[1]);
  printf("%s %s\n","# CPMD output file: ",argv[1]);
   
      
  if (!(fp = fopen(argv[1],"r")))
  {
    fprintf(stderr,"distance: cannot open file %s\n",argv[1]);
    return 1;
  } 

  i = 0;
  while (fgets(line,BUFFER,fp))
  {
    if (strstr(line,"   ATOM          COORDINATES          ") != NULL)
    {
      for (j = 0; j < atom_no[0]; j++)
      fgets(line,BUFFER,fp);
      sscanf(line," %*i %*s %lf %lf %lf %*s",&atom1[0],&atom1[1],&atom1[2]);
      for (j = atom_no[0]; j<atom_no[1]; j++)
      fgets(line,BUFFER,fp);
      sscanf(line," %*i %*s %lf %lf %lf %*s",&atom2[0],&atom2[1],&atom2[2]);
      temp[0] = atom1[0] - atom2[0];
      temp[1] = atom1[1] - atom2[1]; 
      temp[2] = atom1[2] - atom2[2];
      distance = sqrt(temp[0]*temp[0]+temp[1]*temp[1]+temp[2]*temp[2])*BOHR2ANG;
      if (scale_factor < 0)
        scale_factor = distance;
      printf("%4i %9.6lf \n",i,distance/scale_factor);
      i = i + 1;
    }
  }
  return 0;
} /// main