//
// h2ofill.C: try to fill a cell with n molecules by monte-carlo
//
// use: h2ofill seed a b c nmol rcut
//
// Algorithm: try to put an extra molecule at a random position and 
// reject it if it touches any other molecule within rcut
//
#include<iostream>
#include<cstdlib>
#include "D3vector.h"
using namespace std;

D3vector pbc(D3vector r, D3vector cell);

int main(int argc, char **argv)
{
  const double dOH = 0.97 / 0.529177;
  const double tHOH = 104.0 * 3.1415926 / 180.0;
  const double ninety = 90.0 * 3.1415926 / 180.0;
 
  if ( argc != 7 )
  {
    cerr << " use: h2ofill seed a b c nmol rcut" << endl;
    return 1;
  }
  int seed=atoi(argv[1]);
  double a=atof(argv[2]);
  double b=atof(argv[3]);
  double c=atof(argv[4]);
  D3vector cell(a,b,c);
  int ntot=atoi(argv[5]);
  double rc=atof(argv[6]);
  double rc2 = rc * rc;

  srand48(seed);

  D3vector *tauO = new D3vector[ntot];
  D3vector *tauH1 = new D3vector[ntot];
  D3vector *tauH2 = new D3vector[ntot];

  cout << "set cell " << a << " 0 0  0 " << b << " 0  0 0 " << c << endl;
  int n = 0;
  while ( n < ntot )
  {
    double x = a * ( drand48() - 0.5 );
    double y = b * ( drand48() - 0.5 );
    double z = c * ( drand48() - 0.5 );
    D3vector trialO(x,y,z);

    x = drand48() - 0.5;
    y = drand48() - 0.5;
    z = drand48() - 0.5;
    D3vector eOH1 = normalized(D3vector(x,y,z));
    D3vector trialH1 = trialO + dOH * eOH1;

    x = drand48() - 0.5;
    y = drand48() - 0.5;
    z = drand48() - 0.5;
    D3vector eOH2 = normalized(D3vector(x,y,z));
    eOH2 = normalized(eOH2 - (eOH2*eOH1) * eOH1);
    // H1-O-H2 is ninety degrees
    eOH2 = rotate(eOH2, (tHOH-ninety)*(eOH1^eOH2));
    D3vector trialH2 = trialO + dOH * eOH2;

    int accepted = 1;
    for ( int i = 0; i < n && accepted; i++ )
    {
      accepted = ( norm(pbc(trialO-tauO[i],cell)) > rc2 &&
                   norm(pbc(trialH1-tauH1[i],cell)) > rc2 &&
                   norm(pbc(trialH1-tauH2[i],cell)) > rc2 &&
                   norm(pbc(trialH2-tauH1[i],cell)) > rc2 &&
                   norm(pbc(trialH2-tauH2[i],cell)) > rc2 );
    } 
      
    if ( accepted )
    {

      tauO[n] = trialO;
      tauH1[n] = trialH1;
      tauH2[n] = trialH2;

      n++;

      cout << "atom O" << n << " oxygen " << trialO << endl;
      cout << "atom H" << 2*n-1 << " hydrogen " << trialH1 << endl;
      cout << "atom H" << 2*n << " hydrogen " << trialH2 << endl;
    }

  }

  return 0;
}

D3vector pbc(D3vector r, D3vector cell)
{
  double xt = fmod( r.x, cell.x );
  double yt = fmod( r.y, cell.y );
  double zt = fmod( r.z, cell.z );

  double x = xt - cell.x * ( (int) (2*xt/cell.x) );
  double y = yt - cell.y * ( (int) (2*yt/cell.y) );
  double z = zt - cell.z * ( (int) (2*zt/cell.z) );

  return D3vector(x,y,z);
}