/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under 
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing author: Paul Crozier (SNL)
------------------------------------------------------------------------- */

#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_gcmc_molecule.h"
#include "atom.h"
#include "atom_vec.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "group.h"
#include "domain.h"
#include "random_park.h"
#include "force.h"
#include "pair.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;
using namespace MathConst;

/* ---------------------------------------------------------------------- */

FixGCMCMolecule::FixGCMCMolecule(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg)
{
   
  if (narg < 12) error->all(FLERR,"Illegal fix GCMC/molecule command");

  vector_flag = 1;
  size_vector = 6;
  global_freq = 1;
  extvector = 0;
  restart_global = 1;
  time_depend = 1;

  // required args
  
  nevery = atoi(arg[3]);
  nexchanges = atoi(arg[4]);
  nmcmoves = atoi(arg[5]);
  ntype = atoi(arg[6]);
  seed = atoi(arg[7]);
  reservoir_temperature = atof(arg[8]);
  chemical_potential = atof(arg[9]);
  displace = atof(arg[10]);
  nummol = atoi(arg[11]);

  if (ntype <= 0 || ntype > atom->ntypes) 
    error->all(FLERR,"Invalid atom type in fix GCMC command");
  if (nexchanges < 0) error->all(FLERR,"Illegal fix GCMC command");
  if (nmcmoves < 0) error->all(FLERR,"Illegal fix GCMC command");
  if (seed <= 0) error->all(FLERR,"Illegal fix GCMC command");
  if (reservoir_temperature < 0.0) error->all(FLERR,"Illegal fix GCMC command");  
  if (displace < 0.0) error->all(FLERR,"Illegal fix GCMC command"); 

  // compute beta, lambda, sigma, and the zz factor
  
  beta = 1.0/(force->boltz*reservoir_temperature);
  double gas_mass = atom->mass[ntype] + 2.0*atom->mass[ntype+1];  // change atom mass to mass of H2O molecule O ntype=6, H ntype=7 in data-file
  double lambda = sqrt(force->hplanck*force->hplanck/  
		       (2.0*MY_PI*gas_mass*force->mvv2e*
			force->boltz*reservoir_temperature));
  sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e);
  zz = exp(beta*chemical_potential)/(pow(lambda,3));

  
  // set defaults

  molflag = 0;

  // read options from end of input line

  options(narg-12,&arg[12]);

  // random number generator, same for all procs

  random_equal = new RanPark(lmp,seed);

  // random number generator, not the same for all procs

  random_unequal = new RanPark(lmp,seed);
  
  // compute the number of MC cycles that occur nevery timesteps
  
  ncycles = nexchanges + nmcmoves;
  
  // set up reneighboring

  force_reneighbor = 1;
  next_reneighbor = update->ntimestep + 1;

  // zero out counters
  
  nmove_attempts = 0.0;   
  nmove_successes = 0.0;  
  ndel_attempts = 0.0;    
  ndel_successes = 0.0;   
  ninsert_attempts = 0.0; 
  ninsert_successes = 0.0;

  nmax = 0;
  local_gas_list = NULL;
}

/* ---------------------------------------------------------------------- */

FixGCMCMolecule::~FixGCMCMolecule()
{
  delete random_equal;
  delete random_unequal;
  memory->sfree(local_gas_list);
}

/* ---------------------------------------------------------------------- */

int FixGCMCMolecule::setmask()
{
  int mask = 0;
  mask |= PRE_EXCHANGE;
  return mask;
}

/* ---------------------------------------------------------------------- */

void FixGCMCMolecule::init()
{
  // check that no deletable atoms are in atom->firstgroup
  // deleting such an atom would not leave firstgroup atoms first


  int *type = atom->type;
  
  if (atom->firstgroup >= 0) {
    int *mask = atom->mask;
    int nlocal = atom->nlocal;
    int firstgroupbit = group->bitmask[atom->firstgroup];

    int flag = 0;
    for (int i = 0; i < nlocal; i++) {
      if ((type[i] == ntype) && (mask[i] && firstgroupbit)) flag = 1;
      if ((type[i] == ntype+1) && (mask[i] && firstgroupbit)) flag = 1;
    }
    int flagall;
    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);

    if (flagall)
      error->all(FLERR,"Cannot do GCMC on atoms in atom_modify first group");
  }


  if (molflag == 0 && atom->molecule_flag) {
    int *molecule = atom->molecule;
    int *mask = atom->mask;
    int nlocal = atom->nlocal;
    int flag = 0;
    int nmol = 0;
    for (int i = 0; i < nlocal; i++)
      if (type[i] == ntype) {
        if (molecule[i]) flag = 1;
      }
    int flagall;
    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
    //    if (flagall && comm->me == 0)
    //      error->warning(FLERR,"Fix GCMC may delete atom with non-zero molecule ID");
  }    
  //  if (molflag != 0) error->all(FLERR,"Fix GCMC molecule feature does not yet work"); 
  
  if (force->pair->single_enable == 0) 
    error->all(FLERR,"Fix GCMC incompatible with given pair_style");
}

/* ----------------------------------------------------------------------
   attempt particle insertions and deletions
   done before exchange, borders, reneighbor
   so that ghost atoms and neighbor lists will be correct
------------------------------------------------------------------------- */

void FixGCMCMolecule::pre_exchange()
{
  double **x = atom->x;
  int *molecule = atom->molecule;
  double dimz;

  // just return if should not be called on this timestep
  
  if (next_reneighbor != update->ntimestep) return;

  if (domain->triclinic == 0) {
    xlo = domain->boxlo[0];
    xhi = domain->boxhi[0];
    ylo = domain->boxlo[1];
    yhi = domain->boxhi[1];
    zlo = domain->boxlo[2];
    zhi = domain->boxhi[2];
    sublo = domain->sublo;
    subhi = domain->subhi;
  } else {
    xlo = domain->boxlo_bound[0];
    xhi = domain->boxhi_bound[0];
    ylo = domain->boxlo_bound[1];
    yhi = domain->boxhi_bound[1];
    zlo = domain->boxlo_bound[2];
    zhi = domain->boxhi_bound[2];
    sublo = domain->sublo_lamda;
    subhi = domain->subhi_lamda;
  }
    volume = domain->xprd * domain->yprd * domain->zprd;

  // grow local_gas_list array if necessary

  if (atom->nlocal > nmax) {
    memory->sfree(local_gas_list);
    nmax = atom->nmax;
    local_gas_list = (int *) memory->smalloc(nmax*sizeof(int),
					     "GCMCMolecule:local_gas_list");
  }
 
  int *type = atom->type;

  ngas_local = 0;
  for (int i = 0; i < atom->nlocal; i++) {  
    if (type[i] == ntype) {
      local_gas_list[ngas_local] = i;
      ngas_local = ngas_local + 1;
    }      
  }

  MPI_Allreduce(&ngas_local,&ngas,1,MPI_INT,MPI_SUM,world);
  MPI_Scan(&ngas_local,&ngas_before,1,MPI_INT,MPI_SUM,world);
  ngas_before -= ngas_local;

  // perform ncycles MC cycles
  
  for (int i = 0; i < ncycles; i++) {
    int random_int_fraction = 
      static_cast<int>(random_equal->uniform()*ncycles) + 1;	
    if (random_int_fraction <= nmcmoves) {
      attempt_move();
    } else {
     if (random_equal->uniform() < 0.5) attempt_deletion();
      else attempt_insertion();
    }
  }
  
  next_reneighbor = update->ntimestep + nevery;
} 
/*---------------------------------------------------------------------
     Select hydrogens connected to chosen oxygen
-----------------------------------------------------------------------*/
void FixGCMCMolecule::connect(int mol, int o, int *h1, int *h2)
{
  int    iflag;
  iflag=1;
  for (int j=0; j < atom->nlocal; j++) {
    if (atom->molecule[j] == mol) {
      if (j != o) {
        if (iflag == 1) {
          *h1 = j;
          iflag = 0;
        } else {
          *h2 = j;
        }
      }
    }
  }
  /**  int *oxygen;               // there is missing something like this, however with this setup the simulation just freezes
  int *hydrogen1;
  int *hydrogen2;
  MPI_Allreduce(&o,&oxygen,1,MPI_INT,MPI_MAX,world);
  MPI_Allreduce(&h1,&hydrogen1,1,MPI_INT,MPI_MAX,world);
  MPI_Allreduce(&h2,&hydrogen2,1,MPI_INT,MPI_MAX,world);*/

}


/* ----------------------------------------------------------------------
   choose particle randomly across all procs and attempt displacement
------------------------------------------------------------------------- */

void FixGCMCMolecule::attempt_move()
{  
  int i,iwhichglobal,iwhichlocal,j,kH1,kH2,iflag,kH1_local,kH2_local;
  double rx,ry,rz;
  double coordO[3],coordH1[3],coordH2[3];
  double **x = atom->x;
  
  int *molecule = atom->molecule;
 
  nmove_attempts += 1.0;
  
  int success = 0;

  iwhichglobal = static_cast<int> (ngas*random_equal->uniform());
  
  if ((iwhichglobal >= ngas_before) && 
      (iwhichglobal < ngas_before + ngas_local)) {
    iwhichlocal = iwhichglobal - ngas_before;
    i = local_gas_list[iwhichlocal];   // Select random O (water molecule)
    
    connect(molecule[i],i,&kH1, &kH2);  // Look for connecting hydrogens 

    if (molecule[kH1] == molecule[kH2] && molecule[kH1] == molecule[i]) { 

      double energy_before = energy(i,kH1,kH2,ntype,x[i]) + energy(kH1,i,kH2,ntype+1,x[kH1]) + energy(kH2,i,kH1,ntype+1,x[kH2]); //Add up old energies of O and H
      
      double rsq = 1.1;
      
      while (rsq > 1.0) {
	rx = 2*random_unequal->uniform() - 1.0;
	ry = 2*random_unequal->uniform() - 1.0;
	if (domain->dimension == 3) rz = 2*random_unequal->uniform() - 1.0;
	else rz = 0.0;
	rsq = rx*rx + ry*ry + rz*rz;
      }
   
      coordO[0] = x[i][0] + displace*rx;   // calculate new coordinates of displaced water
      coordO[1] = x[i][1] + displace*ry;
      coordO[2] = x[i][2] + displace*rz;
      coordH1[0] = x[kH1][0] + displace*rx;
      coordH1[1] = x[kH1][1] + displace*ry;
      coordH1[2] = x[kH1][2] + displace*rz;
      coordH2[0] = x[kH2][0] + displace*rx;
      coordH2[1] = x[kH2][1] + displace*ry;
      coordH2[2] = x[kH2][2] + displace*rz;
      
      double energy_after = energy(i,kH1,kH2,ntype,coordO) + energy(kH1,i,kH2,ntype+1,coordH1) + energy(kH2,i,kH1,ntype+1,coordH2); // Add up new energies   
      
      if (random_unequal->uniform() < exp(-beta*(energy_after - energy_before))) {  

	x[i][0] = coordO[0];    // Accept new coordinates
	x[i][1] = coordO[1]; 
	x[i][2] = coordO[2];
	x[kH1][0] = coordH1[0]; 
	x[kH1][1] = coordH1[1]; 
	x[kH1][2] = coordH1[2];
	x[kH2][0] = coordH2[0]; 
	x[kH2][1] = coordH2[1]; 
	x[kH2][2] = coordH2[2];
	
	success = 1;

      }
    }
  } 
  
  int success_all = 0;
  MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
  
  if (success_all) {
    nmove_successes += 1.0;  
    comm->borders(); 
  }
}

/* ----------------------------------------------------------------------
   attempt particle deletion
------------------------------------------------------------------------- */

void FixGCMCMolecule::attempt_deletion()
{ 
  int j,kH1,kH2,iflag,k,molID_old,atomO,atomH1,atomH2,check_mol,tagO,tagH1,tagH2,aflag[3][3],l,oxy,H1,H2;
  int angle_atom1,angle_atom2,angle_atom3,bond_atomH1,bond_atomH2,jwhichlocal,lastO,lastH1,lastH2;
  int bond_atomO[2],max,sortflag,val,second_oxy;
  int *molecule=atom->molecule;
  double **x=atom->x, vec[3],ab1,ab2,alpha,period,diff1[3],diff2[3];
  int *type=atom->type;
  int **bond_type = atom->bond_type;
  int **angle_type = atom->angle_type;

  
  ndel_attempts += 1.0;
  
  if (ngas == 1) return;  // Do not delete molecules if only one is left

  int i,iwhichglobal,iwhichlocal,ind,indH1,indH2;
  AtomVec *avec = atom->avec;

  // choose particle randomly across all procs and delete it
  // keep ngas, ngas_local, ngas_before, and local_gas_list current
  // after each deletion
    
  int success = 0;
  iwhichglobal = static_cast<int> (ngas*random_equal->uniform());
  
  if ((iwhichglobal >= ngas_before) && 
      (iwhichglobal < ngas_before + ngas_local)) {
    
    iwhichlocal = iwhichglobal - ngas_before;
    i = local_gas_list[iwhichlocal];   // Select random molecule to be deleted    

    connect(molecule[i],i,&kH1,&kH2);

    molID_old=molecule[i];        // Save molID of selected molecule for later use

    max = 1;
    for (int j=0; j < ngas;j++) {
      if (atom->tag[local_gas_list[j]] > max) {          // Look for oxygen with the highest tagID or molID in the gas list
	atomO = local_gas_list[j];
	max = atom->tag[local_gas_list[j]];
	jwhichlocal = j;
      }
    }
    connect(molecule[atomO],atomO,&atomH1,&atomH2);    
    
    
    if (molecule[kH1] == molecule[kH2] && molecule[kH1] == molecule[i]) {
      
      double deletion_energy = energy(i,kH1,kH2,ntype,atom->x[i]) + energy(kH1,i,kH2,ntype+1,atom->x[kH1]) + energy(kH2,i,kH1,ntype+1,atom->x[kH2]); // Add up energy of O and H     
      if (random_unequal->uniform()  < 
	  ngas*exp(beta*deletion_energy)/(zz*volume)) { 

// if i equals atom with the highest tagID skip this part
	
	if (i != atomO) {             
	  tagO = atom->tag[i];     // Save atomID of selected molecule
	  tagH1 = atom->tag[kH1];
	  tagH2 = atom->tag[kH2];
	  
	  angle_atom1 = atom->angle_atom1[i][0];  // Save connectivities
	  angle_atom2 = atom->angle_atom2[i][0];
	  angle_atom3 = atom->angle_atom3[i][0];
	  bond_atomO[0] = atom->bond_atom[i][0];
	  bond_atomO[1] = atom->bond_atom[i][1];
	  bond_atomH1 = atom->tag[i]; //tagO;
	  bond_atomH2 = atom->tag[i]; //tagO;
	  
	  avec->copy(atomO,i,1);    // copy attributes of molecule with highest tagID onto selected molecule, but keep only their coordinates at the end
	  
	  avec->copy(atomH1,kH1,1);
	  
	  avec->copy(atomH2,kH2,1);
	  
	  check_mol = molecule[atomO];    
	  
	  if (molID_old < molecule[atomO]) { // Keep smaller molID 
	    molecule[i]=molID_old;     
	    molecule[kH1]=molID_old;
	    molecule[kH2]=molID_old;
	  }
	  
	  atom->tag[i] = tagO;          // Recover old tagID of selected molecule
	  atom->tag[kH1] = tagH1;
	  atom->tag[kH2] = tagH2;
	  
	  atom->angle_atom1[i][0] = angle_atom1;  // Recover old connectivities of selected molecule
	  atom->angle_atom2[i][0] = angle_atom2;
	  atom->angle_atom3[i][0] = angle_atom3;
	  atom->bond_atom[i][0] = bond_atomO[0];
	  atom->bond_atom[i][1] = bond_atomO[1];
	  atom->bond_atom[kH1][0] = bond_atomH1;
	  atom->bond_atom[kH2][0] = bond_atomH2;
	  
	  
	  for (int j = 0; j < atom->nlocal; j++) {  // if molID of atomO is not largest
	    if (molecule[j] > check_mol) {
	      molecule[j] = molecule[j]-1;
            }
	  }  
	  
	}	  
	
	lastO = local_gas_list[ngas_local-1];  // Select last molecule in the gaslist      
	
	connect(molecule[lastO],lastO,&lastH1,&lastH2);
	
       	molID_old = molecule[lastO];
	
	avec->copy(lastO,atomO,1);    // copy attributes of last molecule on molecule with highest tagIDs to eliminate it 
	
	avec->copy(lastH1,atomH1,1);
	
	avec->copy(lastH2,atomH2,1);	
	
	
       	// if all the atoms of last correspond to the last three atoms in the atom list, we are done
        // -----------------------------------------------------------------------------------------------------------------------
        // otherwise copy the attributes of the atoms not belonging to the last molecule

	sortflag=0;
	second_oxy = 0;

	for (int j = 0; j < 3; j++) {
	  for (int l = 0; l < 3; l++) {
	    aflag[j][l] = 0;
	  }
	}
	
	if (atom->nlocal-3 == lastO) aflag[0][0] = -1;    // flag the atoms of last with -1 if they coincide with the last three atoms in the atoms list 
	if (atom->nlocal-2 == lastO) aflag[0][1] = -1;   
	if (atom->nlocal-1 == lastO) aflag[0][2] = -1;   
	if (atom->nlocal-3 == lastH1) aflag[1][0] = -1;   
	if (atom->nlocal-2 == lastH1) aflag[1][1] = -1;
	if (atom->nlocal-1 == lastH1) aflag[1][2] = -1;
	if (atom->nlocal-3 == lastH2) aflag[2][0] = -1;   
	if (atom->nlocal-2 == lastH2) aflag[2][1] = -1;
	if (atom->nlocal-1 == lastH2) aflag[2][2] = -1;


	if (molecule[atom->nlocal-2] != molID_old ) {           
	  if (aflag[0][0] == 0 && aflag[0][2] == 0) {  // if lastO is not n-3 nor n-1, copy attributes of n-2 onto it
	    avec->copy(atom->nlocal-2,lastO,1);
	    aflag[0][0] = -1;                    // You should not overwrite lastO in the following steps
	    aflag[0][1] = -1;
	    aflag[0][2] = -1;          
	    
	  } else if (aflag[1][0] == 0 && aflag[1][2] == 0) { // if lastH1 is not n-3 nor n-1, copy atributes of n-2 onto it
	    avec->copy(atom->nlocal-2,lastH1,1);
	    if (type[atom->nlocal-2] == ntype) {     // if an hydrogen has been over written by an oxygen account for the new atomID in the gaslist
	      local_gas_list[ngas_local-2] = lastH1;
	      sortflag = 1;                         // if necessary sort new gaslist
	      second_oxy = 1;
	    }
	    aflag[1][0] = -1;
	    aflag[1][1] = -1;
	    aflag[1][2] = -1;
	    
	  } else if (aflag[2][0] == 0 && aflag[2][2] == 0)  {  // if lastH2 is not n-3 nor n-1, copy atributes of n-2 onto it
	    avec->copy(atom->nlocal-2,lastH2,1);
	    if (type[atom->nlocal-2] == ntype) {      // if an hydrogen has been over written by an oxygen account for the new atomID in the gaslist 
	      local_gas_list[ngas_local-2] = lastH2;
	      sortflag=1;                               // if necessary sort new gaslist 
	      second_oxy = 1;
	    }
	    aflag[2][0] = -1;
	    aflag[2][1] = -1;
	    aflag[2][2] = -1;
	    
	  }
	}
	
	if (molecule[atom->nlocal-3] != molID_old ) {    // The same as above with n-3 not part of last
	  if (aflag[0][1] == 0 && aflag[0][2] == 0) {
	    avec->copy(atom->nlocal-3,lastO,1);
	    aflag[0][0] = -1;
	    aflag[0][1] = -1;
	    aflag[0][2] = -1;
	  } else if (aflag[1][1] == 0 && aflag[1][2] == 0) {
	    avec->copy(atom->nlocal-3,lastH1,1); 
	    if (type[atom->nlocal-3] == ntype) {
	      if (second_oxy == 1) {
		local_gas_list[ngas_local-3] = lastH1;    // if it happens a second time that a hydrogen has been overwritten by an oxygen
	      } else {
		local_gas_list[ngas_local-2] = lastH1;
	      }
	      sortflag=1;
	    }
	    aflag[1][0] = -1;
	    aflag[1][1] = -1;
	    aflag[1][2] = -1;
	    
	  } else if (aflag[2][1] == 0 && aflag[2][2] == 0) {
	    avec->copy(atom->nlocal-3,lastH2,1);
	    if (type[atom->nlocal-3] == ntype) {
	      if (second_oxy == 1) {
		local_gas_list[ngas_local-3] = lastH2;
	      } else {
		local_gas_list[ngas_local-2] = lastH2;
	      }
	      sortflag =1;
	    }
	    aflag[2][0] = -1;
	    aflag[2][1] = -1;
	    aflag[2][2] = -1;
	    
	  }
	}    
	
	if (molecule[atom->nlocal-1] != molID_old ) {    // The same as above with n-1 not part of last
	  if (aflag[0][0] == 0 && aflag[0][1] == 0) {
	    avec->copy(atom->nlocal-1,lastO,1); 
	    aflag[0][0] = -1;
	    aflag[0][1] = -1;
	    aflag[0][2] = -1;
	    
	  } else if (aflag[1][0] == 0 && aflag[1][1] == 0) {
	    avec->copy(atom->nlocal-1,lastH1,1);
	    aflag[1][0] = -1;
	    aflag[1][1] = -1;
	    aflag[1][2] = -1;
	    
	  } else if (aflag[2][0] == 0 && aflag[2][1] == 0) {
	    avec->copy(atom->nlocal-1,lastH2,1);
	    aflag[2][0] = -1;
	    aflag[2][1] = -1;
	    aflag[2][2] = -1;
	    
	  }
	}
	
	// If it happens that atoms with the same tagID are among the last three atoms, special treatment it necessary

	if (atom->tag[atom->nlocal-3] == atom->tag[atom->nlocal-2]) {
	  if (type[atom->nlocal-3] == ntype) {    // if the coinciding atoms are lastO and atomO
	    if (sortflag == 1) second_oxy = 1;
	    sortflag=1;
	    if (molecule[atom->nlocal-1] != molID_old) {  // if n-1 does not belong to last, copy lastO onto lastH2
	      avec->copy(atom->nlocal-3,lastH2,1);
              if (second_oxy == 1) {
                local_gas_list[ngas_local-3] = lastH2;  // Correct id in gaslist
              } else {
                local_gas_list[ngas_local-2] = lastH2;
              }
              aflag[2][2] = -1;
	    } else {                                     // if n-1 belongs to last
	      if (aflag[1][2] == 0 && aflag[2][2] == 0) {  // if not lastH1 nor lastH2 are among the last three atoms
		if (atom->nlocal-1 == atomH1) {          // if n-1 is atomH1 copy n-3 onto atomH2 to eliminate molecule with highest tagID completely
		  avec->copy(atom->nlocal-3,atomH2,1);  
                  if (second_oxy == 1) {                 
                    local_gas_list[ngas_local-3] = atomH2;
                  } else {
                    local_gas_list[ngas_local-2] = atomH2;
                  }
                  aflag[2][2] = -1;
		} else {                               // if n-1 is atomH2 copy n-3 onto atomH1
		  avec->copy(atom->nlocal-3,atomH1,1);
                  if (second_oxy == 1) {
		    local_gas_list[ngas_local-3] = atomH1;
                  } else {
                    local_gas_list[ngas_local-2] = atomH1;
                  }
                  aflag[1][2] = -1;
		}
	      } else if (aflag[1][2] == 0) { // if lastH2 is n-1
		avec->copy(atom->nlocal-3,lastH1,1);
                if (second_oxy == 1) {
                  local_gas_list[ngas_local-3] = lastH1;
                } else {
                  local_gas_list[ngas_local-2] = lastH1;
                }
                aflag[1][2] = -1;
	      } else if (aflag[2][2] == 0) {  // if lastH1 is n-1
		avec->copy(atom->nlocal-3,lastH2,1);
                if (second_oxy == 1) {
                  local_gas_list[ngas_local-3] = lastH2;
                } else {
                  local_gas_list[ngas_local-2] = lastH2;
                } 
                aflag[2][2] = -1;
	      }
	    }
	  } else if (type[atom->nlocal-3] == ntype +1) { // if n-3 and n-2  atoms are atomH1 and lastH1 (or atomH2 and lastH2)
            if (type[atom->nlocal-1] == ntype) {    // if n-1 is lastO
	      if (aflag[1][2] == 0) {
		avec->copy(atom->nlocal-3,lastH1,1);  //if lastH1 is not among last atoms, copy lastH2 onto lastH1
		aflag[1][2] = -1;
	      } else if (aflag[2][2] == 0) {          // if lastH2 is not among last atoms, copy lastH1 onto lastH2
		avec->copy(atom->nlocal-3,lastH2,1);
		aflag[2][2] = -1;
	      }
	    } else if (type[atom->nlocal-1] == ntype + 1) { // if n-1 is a hydrogen
	      if (aflag[0][2] == 0) {
		avec->copy(atom->nlocal-3,lastO,1);   // if n-1 is part of last
		aflag[0][2] = -1;
	      } else if (aflag[1][2] == 0) {         // if n-1 is not part of last
		avec->copy(atom->nlocal-3,lastH1,1);
		aflag[1][2] = -1;
	      
	      }
	    }
	  } 
	}
	if (atom->tag[atom->nlocal-1] == atom->tag[atom->nlocal-2]) {  // the same as above for n-1 equals n-2
	  if (type[atom->nlocal-2] == ntype) {
	    if (sortflag == 1) second_oxy = 1; 
	    sortflag = 1;
	    if (molecule[atom->nlocal-3] != molID_old) {
	      avec->copy(atom->nlocal-2,lastH2,1);
              if (second_oxy == 1) {
                local_gas_list[ngas_local-3] = lastH2;
              } else {
                local_gas_list[ngas_local-2] = lastH2;
              }
              aflag[2][0] = -1;
	    } else {
	      if (aflag[1][0] == 0 && aflag[2][0] == 0) {
		if (atom->nlocal-3 == atomH1) {
		  avec->copy(atom->nlocal-2,atomH2,1);
                  if (second_oxy == 1) {
                    local_gas_list[ngas_local-3] = atomH2;
                  } else {
                    local_gas_list[ngas_local-2] = atomH2;
                  }
                  aflag[2][0] = -1;
		} else {
		  avec->copy(atom->nlocal-2,atomH1,1);
                  if (second_oxy == 1) {
                    local_gas_list[ngas_local-3] = atomH1;
                  } else {
                    local_gas_list[ngas_local-2] = atomH1;
                  }
                  aflag[1][0] = -1;
		}
	      } else if (aflag[1][0] == 0) {
		avec->copy(atom->nlocal-2,lastH1,1);
                if (second_oxy == 1) {
                  local_gas_list[ngas_local-3] = lastH1;
                } else {
                  local_gas_list[ngas_local-2] = lastH1;
                }
                aflag[1][0] = -1;
	      } else if (aflag[2][0] == 0) {
		avec->copy(atom->nlocal-2,lastH2,1);
                if (second_oxy == 1) {
		  local_gas_list[ngas_local-3] = lastH2;
		} else {
		  local_gas_list[ngas_local-2] = lastH2;
                }
                aflag[2][0] = -1;
	      }
	    }
	  } else if (type[atom->nlocal-2] == ntype +1) {
            if (type[atom->nlocal-3] == ntype) {
	      if (aflag[1][0] == 0) {
		avec->copy(atom->nlocal-2,lastH1,1);
		aflag[1][0] = -1;
	      }  else if (aflag[2][0] == 0) {
		avec->copy(atom->nlocal-2,lastH2,1);
		aflag[2][0] = -1;
	      }
	    } else if (type[atom->nlocal-3] == ntype +1) {
	      if (aflag[0][0] == 0) {
		avec->copy(atom->nlocal-2,lastO,1);
		aflag[0][0] = -1;
	      } else if (aflag[1][0] == 0) {
		avec->copy(atom->nlocal-3,lastH1,1);
		aflag[1][0] = -1;
		
	      } 
	    }
	  }    
	}
	if (atom->tag[atom->nlocal-3] == atom->tag[atom->nlocal-1]) {  // the same as above for n-3 equals n-1                                                                                          
	  if (type[atom->nlocal-3] == ntype) {
	    if (sortflag == 1) second_oxy = 1; 
	    sortflag=1;
	    if (molecule[atom->nlocal-2] != molID_old) {
	      avec->copy(atom->nlocal-3,lastH2,1);
              if (second_oxy == 1) {
                local_gas_list[ngas_local-3] = lastH2;
              } else {
                local_gas_list[ngas_local-2] = lastH2;
              }
              aflag[2][1] = -1;
	    } else {
	      if (aflag[1][1] == 0 && aflag[2][1] == 0) {
		if (atom->nlocal-2 == atomH1) {
		  avec->copy(atom->nlocal-3,atomH2,1);
                  if (second_oxy == 1) {
                    local_gas_list[ngas_local-3] = atomH2;
                  } else {
                    local_gas_list[ngas_local-2] = atomH2;
                  }
                  aflag[2][1] = -1;
		} else {
		  avec->copy(atom->nlocal-3,atomH1,1);
                  if (second_oxy == 1) {
                    local_gas_list[ngas_local-3] = atomH1;
                  } else {
                    local_gas_list[ngas_local-2] = atomH1;
                  }
                  aflag[1][1] = -1;
		}
	      } else if (aflag[1][1] == 0) {
		avec->copy(atom->nlocal-3,lastH1,1);
                if (second_oxy == 1) {
		  local_gas_list[ngas_local-3] = lastH1;
                } else {
                  local_gas_list[ngas_local-2] = lastH1;
                }
                aflag[1][1] = -1;
	      } else if (aflag[2][1] == 0) {
		avec->copy(atom->nlocal-3,lastH2,1);
                if (second_oxy == 1) {
                  local_gas_list[ngas_local-3] = lastH2;
                } else {
                  local_gas_list[ngas_local-2] = lastH2;
                }
                aflag[2][1] = -1;
	      }
	    }
	  } else if (type[atom->nlocal-3] == ntype +1) {
            if (type[atom->nlocal-2] == ntype) {
	      if (aflag[1][1] == 0) {
	        avec->copy(atom->nlocal-3,lastH1,1);
                aflag[1][1] = -1;
	      } else if (aflag[2][1] == 0) {
	        avec->copy(atom->nlocal-3,lastH2,1);
		  aflag[2][1] = -1;
              }
            } else if (type[atom->nlocal-2] == ntype +1) {
	      if (aflag[0][1] == 0) {
		avec->copy(atom->nlocal-3,lastO,1);
		aflag[0][1] = -1;
	      } else if (aflag[1][1] == 0) {
		avec->copy(atom->nlocal-3,lastH1,1);
		aflag[1][1] = -1;
		
	      }
	      
	    }
	  }
	}     
	
	
	atom->nbonds = atom->nbonds-2;
	atom->nangles = atom->nangles-1;
	
	if (sortflag == 1) {   // Sort new gaslist if hydrogen has been replaced by oxygen
	  if (second_oxy == 1) {
	    j=4;
	    while (j <= ngas_local) {
	      if (local_gas_list[ngas_local-j+1] < local_gas_list[ngas_local-j]) {
		val = local_gas_list[ngas_local-j];
		local_gas_list[ngas_local-j] = local_gas_list[ngas_local-j+1];
		local_gas_list[ngas_local-j+1] = val;
		j=j+1;
	      } else {
		j=ngas_local+1;
	      }
	    }
	  }
	  j=3;
	  while (j <= ngas_local) {
	    if (local_gas_list[ngas_local-j+1] < local_gas_list[ngas_local-j]) {
	      val = local_gas_list[ngas_local-j];
	      local_gas_list[ngas_local-j] = local_gas_list[ngas_local-j+1];
	      local_gas_list[ngas_local-j+1] = val;
	      j=j+1;
	    } else {
	      j=ngas_local+1;
	    }
	  }
	}
	
	// Skip last molecule of the gaslist

	ngas_local = ngas_local - 1;
	
        // Skip last three atoms of atoms list	
	
	atom->nlocal = atom->nlocal -3;
	
	success = 1; 
	
      }
    }
  }  
  
  int success_all = 0;
  MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
  
  if (success_all) {
    ngas = ngas - 1;
    ndel_successes += 1.0;
    atom->natoms = atom->natoms - 3;
    if (iwhichglobal < ngas_before) ngas_before = ngas_before - 1;      
    comm->borders();
    if (atom->tag_enable) {
      atom->tag_extend();
      if (atom->map_style) {
	atom->map_init();
	atom->map_set();
      }
      
    }
  }
  
}
  
  
/* ----------------------------------------------------------------------
   attempt particle insertion
------------------------------------------------------------------------- */

void FixGCMCMolecule::attempt_insertion()
{  
  int flag,success,ind,indH1,indH2,iflag,k;
  double coordO[3],lamdaO[3],coordH1[3],lamdaH1[3],coordH2[3],lamdaH2[3],vec[3],ab1,ab2,alpha,delta;
  double *newcoordO, *newcoordH1, *newcoordH2; 
  ninsert_attempts += 1.0;
  
  // choose random position for new atom within box

  coordO[0] = xlo + random_equal->uniform() * (xhi-xlo);
  coordO[1] = ylo + random_equal->uniform() * (yhi-ylo);
  coordO[2] = zlo + random_equal->uniform() * (zhi-zlo);

// Define coordinates of H such that angle is around 110 degrees and bond length is 1 A

  coordH1[0] = coordO[0] + 0.5736;
  coordH1[1] = coordO[1];
  coordH1[2] = coordO[2] + 0.81915;
 
  coordH2[0] = coordO[0] + 0.5736;
  coordH2[1] = coordO[1];
  coordH2[2] = coordO[2] - 0.81915;

  // check if new atoms are in my sub-box or above it if I'm highest proc
  // if so, add to my list via create_atom()
  // initialize info about the atoms
  // set group mask to "all" plus fix group
  
  if (domain->triclinic) {
    domain->x2lamda(coordO,lamdaO);
    domain->x2lamda(coordH1,lamdaH1);
    domain->x2lamda(coordH2,lamdaH2);
    newcoordO = lamdaO;
    newcoordH1 = lamdaH1;
    newcoordH2 = lamdaH2;
  } else {
    newcoordO = coordO;
    newcoordH1 = coordH1;
    newcoordH2 = coordH2;
  }


  flag = 0;
 
  if (newcoordO[0] >= sublo[0] && newcoordO[0] < subhi[0] &&
    newcoordO[1] >= sublo[1] && newcoordO[1] < subhi[1] &&
    newcoordO[2] >= sublo[2]  && newcoordO[2] < subhi[2] &&     
    newcoordH1[0] >= sublo[0] && newcoordH1[0] < subhi[0] &&
    newcoordH1[1] >= sublo[1] && newcoordH1[1] < subhi[1] &&
    newcoordH1[2] >= sublo[2]  && newcoordH1[2] < subhi[2] &&   
    newcoordH2[0] >= sublo[0] && newcoordH2[0] < subhi[0] &&
    newcoordH2[1] >= sublo[1] && newcoordH2[1] < subhi[1] &&
      newcoordH2[2] >= sublo[2] && newcoordH2[2] < subhi[2] ) flag = 1; 
    
  success = 0;

  if (flag) {
    int nall = atom->nlocal + atom->nghost;
    double insertion_energy = energy(nall,nall + 1,nall + 2,ntype,coordO) + energy(nall+1,nall, nall +2,ntype+1,coordH1) + energy(nall+2,nall, nall+1,ntype+1,coordH2); // Add up energy of inserted molecule

    if (random_unequal->uniform() < 
	zz*volume*exp(-beta*insertion_energy)/(ngas+1)) {
  
      atom->avec->create_atom(ntype,coordO);
      atom->avec->create_atom(ntype+1,coordH1);  // create new water molecule typeID of hydrogen must be ntype+1
      atom->avec->create_atom(ntype+1,coordH2);

      int m = atom->nlocal - 3;     // define new tagID of oxygen
   
      atom->type[m] = ntype;
      atom->molecule[m] = ngas + 1 + nummol; // allocate new moleculeID, thus it was important to eliminate the atom with the highest molID from the system, nummol gives the number of molecules other than water 
      atom->mask[m] = 1 | groupbit;
      atom->v[m][0] = random_unequal->gaussian()*sigma;
      atom->v[m][1] = random_unequal->gaussian()*sigma;
      atom->v[m][2] = random_unequal->gaussian()*sigma;

      atom->type[m+1] = ntype +1;
      atom->molecule[m+1] = ngas + 1 + nummol;  
      atom->mask[m+1] = 1 | groupbit;
      atom->v[m+1][0] = random_unequal->gaussian()*sigma;
      atom->v[m+1][1] = random_unequal->gaussian()*sigma;
      atom->v[m+1][2] = random_unequal->gaussian()*sigma;

      atom->type[m+2] = ntype +1;
      atom->molecule[m+2] = ngas + 1 +nummol; 
      atom->mask[m+2] = 1 | groupbit;
      atom->v[m+2][0] = random_unequal->gaussian()*sigma;
      atom->v[m+2][1] = random_unequal->gaussian()*sigma;
      atom->v[m+2][2] = random_unequal->gaussian()*sigma;
      
      k = local_gas_list[0]; // Take attributes from the first molecule in the gaslist


     connect(atom->molecule[k],k,&indH1,&indH2);
        
     atom->q[m]=atom->q[k];     // Define charges
     atom->q[m+1]=atom->q[k]/-2.0;
     atom->q[m+2]=atom->q[k]/-2.0;

     atom->bond_type[m][0]=atom->bond_type[k][0];  // Define bond types of oxygen
     atom->bond_type[m][1]=atom->bond_type[k][1];
     atom->bond_atom[m][0]=m+2;                    // Define connectivity
     atom->bond_atom[m][1]=m+3;
     atom->num_bond[m]=atom->num_bond[k];         // Define number of bonds per oxygen

     atom->bond_type[m+1][0]=atom->bond_type[indH1][0];  // bond types of hydrogen
     atom->bond_type[m+2][0]=atom->bond_type[indH2][0];
     atom->bond_atom[m+1][0]=m+1;       // Connectivity
     atom->bond_atom[m+2][0]=m+1;
     atom->num_bond[m+1]=atom->num_bond[indH1];  // Number of bonds per hydrogen
     atom->num_bond[m+2]=atom->num_bond[indH2];

     atom->nbonds = atom->nbonds+2;             

     atom->angle_type[m][0]=atom->angle_type[k][0]; // define angle types
     atom->angle_atom1[m][0]=m+2;  // Connectivities
     atom->angle_atom2[m][0]=m+1;
     atom->angle_atom3[m][0]=m+3;
     atom->num_angle[m]=atom->num_angle[k]; // Number of angles per oxygen

     atom->nangles=atom->nangles+1;
    
     int nfix = modify->nfix;
     Fix **fix = modify->fix;
     for (int j = 0; j < nfix; j++) {
       if (fix[j]->create_attribute) fix[j]->set_arrays(m);
       if (fix[j]->create_attribute) fix[j]->set_arrays(m+1);
       if (fix[j]->create_attribute) fix[j]->set_arrays(m+2);
     }
     if (atom->nlocal > nmax) { // Size of maxmol would be sufficient
       nmax = atom->nmax;
       local_gas_list = (int *) 
	 memory->srealloc(local_gas_list,nmax*sizeof(int),
			  "GCMC:local_gas_list");
     }
     ngas_local = ngas_local + 1;      
     local_gas_list[ngas_local-1] = atom->nlocal-3;
     
     
     success = 1; 
     
    }
  } 
  int success_all = 0;
  MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
  
  if (success_all) {
    ngas=ngas + 1; // Increase water molecules by 1
    count_mol=count_mol + 1;
    ninsert_successes += 1.0;
    MPI_Scan(&ngas_local,&ngas_before,1,MPI_INT,MPI_SUM,world);
    ngas_before = ngas_before - ngas_local;
    comm->borders();
    if (atom->tag_enable) {
      atom->natoms=atom->natoms + 3; // Increase number of atoms by 3
      atom->tag_extend();
      if (atom->map_style) {
        atom->map_init();
        atom->map_set();
      }
    }   
  } 
}

/* ---------------------------------------------------------------------
    Sort molecules in the water list
------------------------------------------------------------------------*/
void FixGCMCMolecule::sort(int *list,int n)
{
  int val,i,j,temp;
  for (int i=0; i< n; j++) {
      for (int j = n; j< i+2; i--) {
          if (list[j-1] > list[j]) {
	    temp = list[j-1];
            list[j-1]=list[j];
	    list[j] = temp;
          }
       }  
   }
}
/* ----------------------------------------------------------------------
   compute particle's interaction energy with the rest of the system
------------------------------------------------------------------------- */

double FixGCMCMolecule::energy(int i, int k, int l, int wtype,double *coord) 
{
  double delx,dely,delz,rsq;

  double **x = atom->x;
  int *type = atom->type;
  int nall = atom->nlocal + atom->nghost;
  pair = force->pair;
  cutsq = force->pair->cutsq;

  double fpair = 0.0;
  double factor_coul = 1.0;
  double factor_lj = 1.0;
  
  double total_energy = 0.0;
  for (int j = 0; j < nall; j++) {

    if (i == j) continue;
    if (k == j) continue;
    if (l == j) continue;

    delx = coord[0] - x[j][0];
    dely = coord[1] - x[j][1];
    delz = coord[2] - x[j][2];
    rsq = delx*delx + dely*dely + delz*delz;
    int jtype = type[j];

    if (rsq < cutsq[wtype][jtype])
      total_energy += 
	pair->single(i,j,wtype,jtype,rsq,factor_coul,factor_lj,fpair);
  }

  return total_energy;
}

/* ----------------------------------------------------------------------
   parse optional parameters at end of input line 
------------------------------------------------------------------------- */

void FixGCMCMolecule::options(int narg, char **arg)
{
  if (narg < 0) error->all(FLERR,"Illegal fix GCMC command");

  int iarg = 0;
  while (iarg < narg) {
    if (strcmp(arg[iarg],"molecule") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix GCMC command");
      if (strcmp(arg[iarg+1],"no") == 0) molflag = 0;
      else if (strcmp(arg[iarg+1],"yes") == 0) molflag = 1;
      else error->all(FLERR,"Illegal fix evaporate command");
      iarg += 2;
    } else error->all(FLERR,"Illegal fix GCMC command");
  }
}

/* ----------------------------------------------------------------------
  return acceptance ratios
------------------------------------------------------------------------- */

double FixGCMCMolecule::compute_vector(int n)
{
  if (n == 0) return nmove_attempts;
  if (n == 1) return nmove_successes;
  if (n == 2) return ndel_attempts;
  if (n == 3) return ndel_successes;
  if (n == 4) return ninsert_attempts;
  if (n == 5) return ninsert_successes;
  return 0.0;
}

/* ----------------------------------------------------------------------
   memory usage of local atom-based arrays
------------------------------------------------------------------------- */

double FixGCMCMolecule::memory_usage()
{
  double bytes = nmax * sizeof(int);
  return bytes;
}

/* ----------------------------------------------------------------------
   pack entire state of Fix into one write 
------------------------------------------------------------------------- */

void FixGCMCMolecule::write_restart(FILE *fp)
{
  int n = 0;
  double list[4];
  list[n++] = random_equal->state();
  list[n++] = random_unequal->state();
  list[n++] = next_reneighbor;

  if (comm->me == 0) {
    int size = n * sizeof(double);
    fwrite(&size,sizeof(int),1,fp);
    fwrite(list,sizeof(double),n,fp);
  }
}

/* ----------------------------------------------------------------------
   use state info from restart file to restart the Fix  
------------------------------------------------------------------------- */

void FixGCMCMolecule::restart(char *buf)
{
  int n = 0;
  double *list = (double *) buf;

  seed = static_cast<int> (list[n++]);
  random_equal->reset(seed);
  
  seed = static_cast<int> (list[n++]);
  random_unequal->reset(seed);
  
  next_reneighbor = static_cast<int> (list[n++]);
}