///////////////////////////////////////////////////////////////////////////////
//
// hbonds.cpp
//
////////////////////////////////////////////////////////////////////////////////
//
// Compute average number of hydrogen bonds in a water sample
// Luzar-Chandler criterion
// A. Luzar and D. Chandler
// J. Chem. Phys. 98, p. 8160 (1993).
// dOO < 3.5 Angstron and H-O-O angle < 30 deg
//
// 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: hbonds < file.xyz
// output:
// average number of hydrogen bonds per molecule
//
// Dependencies: D3vector.h, UnitCell.C from Qbox

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

int main(int argc, char** argv)
{
  const double rcut = 3.5;

  UnitCell cell;
  string nm;
  vector<string> name;
  vector<D3vector> tau;
  int na,nconfig=0,npairs=0;
  vector<int> O_index, H_index;
  vector<int> nhbonds(8);

  if ( argc != 1 )
  {
    cerr << " use: hbonds < file.xyz" << endl;
    return 1;
  }

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

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

  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 O_index vector
    // indices of H atoms in H_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' ) O_index.push_back(i);
        if ( nm[0] == 'H' ) H_index.push_back(i);
      }
    }

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

    const double s32 = 0.5*sqrt(3.0);
    // scan all oxygen pairs and select pairs with distance < rcut
    for ( int i = 0; i < O_index.size(); i++ )
    {
      int nhb=0; // number of H bonds for atom i
      for ( int j = 0; j < O_index.size(); j++ )
      {
        if ( i!=j )
        {
          D3vector rOO = tau[O_index[j]] - tau[O_index[i]];
          cell.fold_in_ws(rOO);
          if ( length(rOO) < rcut )
          {
            npairs++;
            // check for H atom within rcut of first oxygen
            for ( int k = 0; k < H_index.size(); k++ )
            {
              D3vector rOiH = tau[H_index[k]] - tau[O_index[i]];
              cell.fold_in_ws(rOiH);
              if ( length(rOiH) < rcut )
              {
                D3vector rOjH = tau[H_index[k]] - tau[O_index[j]];
                cell.fold_in_ws(rOjH);
                if ( length(rOjH) < rcut )
                {
                  D3vector eOiH = normalized(rOiH);
                  D3vector eOjH = normalized(rOjH);
                  D3vector eOO = normalized(rOO);
                  double c1 = eOiH * eOO;
                  double c2 = -eOjH * eOO;
                  if ( (c1 > s32) && (c2 > s32) ) nhb++;
                }
              }
            }
          }
        }
      }
      // atom O i has nhb hydrogen bonds
      assert(nhb<10);
      nhbonds[nhb]++;
    }
    cin >> na;
  }

  int nO = O_index.size();
  //double avgcoord = ((double)npairs)/(nconfig*nO);

  int nhbsum = 0;
  int nhbisum = 0;
  for ( int i = 0; i < nhbonds.size(); i++ )
  {
    nhbsum += nhbonds[i];
    nhbisum += i * nhbonds[i];
  }
  double avgnhbonds = ((double)nhbisum)/(nconfig*nO);

  cout << "#   Nhb   n[0]   n[1]   n[2]   n[3]   n[4]   n[5]   n[6]   n[7]"
       << endl;
  cout << setprecision(4) << setw(7) << fixed << right << avgnhbonds;
  for ( int i = 0; i < nhbonds.size(); i++ )
    cout << setprecision(4) << setw(7) << fixed << right
         << ((double)nhbonds[i])/nhbsum;
  cout << endl;

  return 0;
}