///////////////////////////////////////////////////////////////////////////////
//
// gOOO.cpp
//
////////////////////////////////////////////////////////////////////////////////
//
// Compute distribution of OOO angles in a water sample
//
// input: xyz file including cell size information
// first line: natoms
// second line: iteration_number ax ay az  bx by bz cx cy cz (Angstrom)
// next lines: atom_name x y z (Angstrom)
// [.. repeat all above lines for each iteration]
//
// use: gOOO rcut[Angstrom] da < file.xyz > gOOO.dat
//   rcut: cutoff distance [Angstrom]
//   da: bin size in angle histogram [degrees]
// output:
// for each time step:
//   for each atom O1 at each time step:
//     for each pair of neighbors O2 O3:
//       print the angle O2-O1-O3
//
// Dependencies: D3vector.h, UnitCell.C from Qbox

#include<iostream>
#include<vector>
#include<algorithm>
#include<cassert>
#include "D3vector.h"
#include "UnitCell.h"
#include <string>
using namespace std;

int main(int argc, char** argv)
{
  UnitCell cell;
  string nm;
  vector<string> name;
  vector<D3vector> tau;
  unsigned long int na,nconfig=0,nbonds=0;
  vector<int> index;

  if ( argc != 3 )
  {
    cerr << " use: gOOO rcut[Angstrom] da < file.xyz > gOOO.dat" << endl;
    return 1;
  }
  const double rcut = atof(argv[1]);
  const double da = atof(argv[2]);
  const int nbin = (int) (180.0/da);
  vector<int> count(nbin);

  D3vector a0,a1,a2;
  cin >> na;

  tau.resize(na);
  name.resize(na);

  double sumcoord = 0.0;

  while ( !cin.eof() )
  {
    nconfig++;

    // read iter number and cell parameters on same line as na
    int iter;
    cin >> iter >> a0 >> a1 >> a2;
    while( !cin.eof() && cin.get() != '\n');
    cell.set(a0,a1,a2);

    // gather indices of O atoms in index vector
    for ( int i = 0; i < na; i++ )
    {
      cin >> nm >> tau[i];
      while( !cin.eof() && cin.get() != '\n');
      if ( nconfig == 1 )
      {
        name[i] = nm;
        // name[i] now contains the name of the species
        // check if nm starts with 'O'
        if ( nm[0] == 'O' ) index.push_back(i);
      }
    }

    // the number of oxygen atoms is index.size()
    // the position of atom i is tau[i]
    // the name of atom i is name[i]

    // build the neighbor list of atom i
    for ( int i = 0; i < index.size(); i++ )
    {
      vector<int> neighbor_list;
      for ( int j = 0; j < index.size(); j++ )
      {
        if ( i != j )
        {
          D3vector r12 = tau[index[j]] - tau[index[i]];
          cell.fold_in_ws(r12);
          if ( length(r12) < rcut )
          neighbor_list.push_back(j);
        }
      }
      //cout << "neighbor_list[" << i << "] ";
      //for ( int nj = 0; nj < neighbor_list.size(); nj++ )
      //  cout << " " << neighbor_list[nj];
      //cout << endl;

      sumcoord += neighbor_list.size();

      // process all pairs that are neighbors of tau[index[i]]
      for ( int nj = 0; nj < neighbor_list.size(); nj++ )
      for ( int nk = nj+1; nk < neighbor_list.size(); nk++ )
      {
        int j = neighbor_list[nj];
        int k = neighbor_list[nk];
        D3vector r12 = tau[index[j]] - tau[index[i]];
        cell.fold_in_ws(r12);
        D3vector r32 = tau[index[k]] - tau[index[i]];
        cell.fold_in_ws(r32);
        const double sp = normalized(r12) * normalized(r32);
        const double c = max(-1.0,min(1.0,sp));
        const double a = (180.0/M_PI)*acos(c);
        int ia = (int) (a/da+0.5);
        // angle is in bin ia, i.e. angle in [ia*da-0.5*da, ia*da+0.5*da]
        if ( ia < nbin ) count[ia]++;
      }
    }

    cin >> na;
  }

  // compute average coordination number
  double avgcoord = sumcoord/(nconfig*index.size());

  // normalization
  double sum = 0;
  for ( int i = 1; i < nbin; i++ )
    sum += da*count[i];

  // output distribution
  cout << "# gOOO rcut=" << rcut << " avgcoord=" << avgcoord << endl;
  for ( int i = 1; i < nbin; i++ )
  {
    const double a = i*da;
    cout << a << " " << count[i]/sum << endl;
  }
  return 0;
}