///////////////////////////////////////////////////////////////////////////////
//
// lsi.cpp
//
////////////////////////////////////////////////////////////////////////////////
// Compute the local structure index defined by Shiratani and Sasai
// E. Shiratani and M. Sasai, J. Chem. Phys. 104, 7671 (1996)
// E. Shiratani and M. Sasai, J. Chem. Phys. 108, 3264 (1998)
//
// 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: lsi < file.xyz > lsi.dat
//
// example: lsi < h2orun.xyz
//
// output:
// for each atom at each time step: distance to the first six O neighbors
// and LSI value.
//
// 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;
  int na,nconfig=0;
  vector<int> index;

  if ( argc != 1 )
  {
    cerr << " use: lsi < file.xyz > lsi.dat" << 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);

    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 atoms of species spname is index.size()
    // the position of atom i is tau[i]

    vector<double> r; // distances from atom i
    for ( int i = 0; i < index.size(); i++ )
    {
      r.clear();
      // loop on j != i
      for ( int j = 0; j < index.size(); j++ )
      {
        if ( i != j )
        {
          D3vector t = tau[index[i]] - tau[index[j]];
          cell.fold_in_ws(t);
          r.push_back(length(t));
        }
      }
      cout.width(4);
      cout << name[index[i]] << " ";
      // r contains the distances from atom i
      sort(r.begin(),r.end());
      cout << r[0] << " " << r[1] << " " << r[2]
           << " " << r[3] << " " << r[4] << " " << r[5];
      // compute the LSI
      // compute number of atoms within 3.7 Angstrom
      int n = 0;
      for ( int i = 0; i < r.size(); i++ )
        if ( r[i] < 3.7 ) n++;

      // check that there is at least one atom outside of the sphere
      // to compute delta[i]
      assert(r.size()>n);

      // compute delta[i] for each atom
      vector<double> d(n);
      for ( int i = 0; i < n; i++ )
        d[i] = r[i+1] - r[i];

      // average d
      double davg = (r[n]-r[0])/n;

      double lsi = 0.0;
      for ( int i = 0; i < n; i++ )
        lsi += ( d[i] - davg ) * ( d[i] - davg );
      lsi /= n;
      cout << " lsi: " << lsi << " Å^2" << endl;
    }
    cin >> na;
  }
  return 0;
}