/* ----------------------------------------------------------------------
   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: Aidan Thompson (SNL)
------------------------------------------------------------------------- */

#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "tad.h"
#include "universe.h"
#include "update.h"
#include "atom.h"
#include "domain.h"
#include "region.h"
#include "comm.h"
#include "velocity.h"
#include "integrate.h"
#include "min.h"
#include "neighbor.h"
#include "modify.h"
#include "neb.h"
#include "compute.h"
#include "fix.h"
#include "fix_event_tad.h"
#include "fix_store_state.h"
#include "force.h"  
#include "pair.h"
#include "random_park.h"
#include "random_mars.h"
#include "output.h"
#include "dump.h"
#include "finish.h"
#include "timer.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

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

TAD::TAD(LAMMPS *lmp) : Pointers(lmp) {}

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

TAD::~TAD()
{
  memory->sfree(fix_event_list);
  if (neb_logfilename != NULL) delete [] neb_logfilename;
  delete [] min_style;
  delete [] min_style_neb;
}

/* ----------------------------------------------------------------------
   perform TAD simulation on root proc
   other procs only used for NEB calcs
------------------------------------------------------------------------- */

void TAD::command(int narg, char **arg)
{

  fix_event_list = NULL;
  n_event_list = 0;
  nmax_event_list = 0;
  nmin_event_list = 10;

  // error checks

  if (domain->box_exist == 0) 
    error->all("Tad command before simulation box is defined");
  if (universe->nworlds == 1) 
    error->all("Cannot use TAD with a single replica for NEB");
  if (universe->nworlds != universe->nprocs)
    error->all("Can only use TAD with 1-processor replicas for NEB");
  if (atom->sortfreq > 0)
    error->all("Cannot use TAD with atom_modify sort enabled for NEB");
  if (atom->map_style == 0) 
    error->all("Cannot use TAD unless atom map exists for NEB");

  if (narg < 7) error->universe_all("Illegal tad command");

  nsteps = atoi(arg[0]);
  t_event = atoi(arg[1]);
  templo = atof(arg[2]);
  temphi = atof(arg[3]);
  delta_conf = atof(arg[4]);
  tmax = atof(arg[5]);
  
  char id_compute[strlen(arg[6])+1];
  strcpy(id_compute,arg[6]);
  
  options(narg-7,&arg[7]);

  // total # of timesteps must be multiple of t_event

  if (t_event <= 0) error->universe_all("Invalid t_event in tad command");
  if (nsteps % t_event) 
    error->universe_all("TAD nsteps must be multiple of t_event");

  if (delta_conf <= 0.0 || delta_conf >= 1.0) 
    error->universe_all("Invalid delta_conf in tad command");

  if (tmax <= 0.0) 
    error->universe_all("Invalid tmax in tad command");

  // deltconf = (ln(1/delta))/freq_min (timestep units)

  deltconf = -log(delta_conf)*tmax/update->dt;

  // local storage

  int me_universe = universe->me;
  int nprocs_universe = universe->nprocs;

  MPI_Comm_rank(world,&me);
  MPI_Comm_size(world,&nprocs);

  delta_beta = (1.0/templo - 1.0/temphi) / force->boltz;
  ratio_beta = templo/temphi; 

  // create FixEventTAD object to store last event

  int narg2 = 3;
  char **args = new char*[narg2];
  args[0] = (char *) "tad_event";
  args[1] = (char *) "all";
  args[2] = (char *) "EVENT/TAD";
  modify->add_fix(narg2,args);
  fix_event = (FixEventTAD *) modify->fix[modify->nfix-1];
  delete [] args;

  // create FixStoreState object to store revert state

  narg2 = 13;
  args = new char*[narg2];
  args[0] = (char *) "tad_revert";
  args[1] = (char *) "all";
  args[2] = (char *) "store/state";
  args[3] = (char *) "0";
  args[4] = (char *) "x";
  args[5] = (char *) "y";
  args[6] = (char *) "z";
  args[7] = (char *) "ix";
  args[8] = (char *) "iy";
  args[9] = (char *) "iz";
  args[10] = (char *) "vx";
  args[11] = (char *) "vy";
  args[12] = (char *) "vz";
  modify->add_fix(narg2,args);
  fix_revert = (FixStoreState *) modify->fix[modify->nfix-1];
  delete [] args;

  // create Finish for timing output

  finish = new Finish(lmp);

  // assign FixEventTAD to event-detection compute
  // necessary so it will know atom coords at last event

  int icompute = modify->find_compute(id_compute);
  if (icompute < 0) error->all("Could not find compute ID for TAD");
  compute_event = modify->compute[icompute];
  compute_event->reset_extra_compute_fix("tad_event");

  // reset reneighboring criteria since will perform minimizations

  neigh_every = neighbor->every;
  neigh_delay = neighbor->delay;
  neigh_dist_check = neighbor->dist_check;

  if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
    if (me_universe == 0) 
      error->warning("Resetting reneighboring criteria during TAD");
  }

  neighbor->every = 1;
  neighbor->delay = 0;
  neighbor->dist_check = 1;

  // initialize TAD as if one long dynamics run

  update->whichflag = 1;
  update->nsteps = nsteps;
  update->beginstep = update->firststep = update->ntimestep;
  update->endstep = update->laststep = update->firststep + nsteps;
  update->restrict_output = 1;
  if (update->laststep < 0 || update->laststep > MAXBIGINT)
    error->all("Too many timesteps");

  lmp->init();

  // set minimize style for quench

  narg2 = 1;
  args = new char*[narg2];
  args[0] = min_style;

  update->create_minimize(narg2,args);

  delete [] args;

  // init minimizer settings and minimizer itself

  update->etol = etol;
  update->ftol = ftol;
  update->max_eval = maxeval;

  update->minimize->init();

  // perform TAD simulation

  if (me_universe == 0 && universe->uscreen) 
    fprintf(universe->uscreen,"Setting up TAD ...\n");

  if (me_universe == 0) {
    if (universe->uscreen) 
      fprintf(universe->uscreen,
	      "Step CPU N M Status Barrier Margin t_lo delt_lo\n"
	      );
    if (universe->ulogfile) 
      fprintf(universe->ulogfile,
	      "Step CPU N M Status Barrier Margin t_lo delt_lo\n"
	      );
  }

  ulogfile_lammps = universe->ulogfile;
  uscreen_lammps = universe->uscreen;
  ulogfile_neb = NULL;
  uscreen_neb = NULL;
  if (me_universe == 0 && neb_logfilename) 
    ulogfile_neb = fopen(neb_logfilename,"w");

  // store hot state and quenched event, only on replica 0

  // need this line if quench() does only setup_minimal()
  // update->minimize->setup();

  // This should work with if uncommented, but does not
  // if (universe->iworld == 0) {
  fix_event->store_state();
  quench();
  timer->barrier_start(TIME_LOOP);
  time_start = timer->array[TIME_LOOP];
  fix_event->store_event_tad(update->ntimestep);
  log_event(0);
  fix_event->restore_state();

  // do full init/setup
  
  update->whichflag = 1;
  lmp->init();
  update->integrate->setup();
  //   }
    
  // main loop: look for events until out of time
  // (1) dynamics, store state, quench, check event, restore state
  // (2) if event, perform NEB, record in fix_event_list
  // (3) if confident, pick earliest event 

  nbuild = ndanger = 0;
  time_neb = time_dynamics = time_quench = time_comm = time_output = 0.0;

  timer->barrier_start(TIME_LOOP);
  time_start = timer->array[TIME_LOOP];
  
  int confident_flag, event_flag;

  if (universe->iworld == 0) {
    while (update->ntimestep < update->endstep) {

      // initialize list of possible events
      
      initialize_event_list();
      confident_flag = 0;

      while (update->ntimestep < update->endstep) {
	event_flag = 0;
	while (update->ntimestep < update->endstep) {

	  dynamics();


	  fix_event->store_state();


	  quench();


	  event_flag = check_event();
	  MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld);

	  if (event_flag) break;

	  // restore hot state

	  fix_event->restore_state();

	  // store hot state in revert

	  fix_revert->end_of_step();
	}
	if (!event_flag) break;

	add_event();

	perform_neb(n_event_list-1);
	compute_tlo(n_event_list-1);
	confident_flag = check_confidence();
	MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld);
	if (confident_flag) break;
	if (universe->iworld == 0) revert();
      }
      if (!confident_flag) break;
    
      perform_event(event_first);

      // need to sync timestep with TAD 

      MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld);
    
      int restart_flag = 0;
      if (output->restart_every && universe->iworld == 0)
	if (fix_event->event_number % output->restart_every == 0)
	  restart_flag = 1;

      // full init/setup since are starting after event

      update->whichflag = 1;
      lmp->init();
      update->integrate->setup();

    // write restart file of hot coords

      if (restart_flag) {
	timer->barrier_start(TIME_LOOP);
	output->write_restart(update->ntimestep);
	timer->barrier_stop(TIME_LOOP);
	time_output += timer->array[TIME_LOOP];
      }
    }

  } else {

    while (update->ntimestep < update->endstep) {
      confident_flag = 0;
      while (update->ntimestep < update->endstep) {
	event_flag = 0;
	while (update->ntimestep < update->endstep) {
	  update->ntimestep += t_event;
	  MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld);

	  if (event_flag) break;
	}
	if (!event_flag) break;
	perform_neb(-1);
	MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld);
	if (confident_flag) break;
      }
      if (!confident_flag) break;
    
      // need to sync timestep with TAD 

      MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld);
    }
  }

  // set total timers and counters so Finish() will process them

  timer->array[TIME_LOOP] = time_start;
  timer->barrier_stop(TIME_LOOP);

  timer->array[TIME_PAIR] = time_neb;
  timer->array[TIME_BOND] = time_dynamics;
  timer->array[TIME_KSPACE] = time_quench;
  timer->array[TIME_COMM] = time_comm;
  timer->array[TIME_OUTPUT] = time_output;

  neighbor->ncalls = nbuild;
  neighbor->ndanger = ndanger;

  if (me_universe == 0) {
    if (universe->uscreen) 
      fprintf(universe->uscreen,
	      "Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT 
	      " atoms\n",
	      timer->array[TIME_LOOP],nprocs_universe,nsteps,atom->natoms);
    if (universe->ulogfile) 
      fprintf(universe->ulogfile,
	      "Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT 
	      " atoms\n",
              timer->array[TIME_LOOP],nprocs_universe,nsteps,atom->natoms);
  }

  if (me_universe == 0) fclose(ulogfile_neb);
  
  finish->end(3);
  
  update->whichflag = 0;
  update->firststep = update->laststep = 0;
  update->beginstep = update->endstep = 0;
  update->restrict_output = 0;

  // reset reneighboring criteria

  neighbor->every = neigh_every;
  neighbor->delay = neigh_delay;
  neighbor->dist_check = neigh_dist_check;

  delete finish;
  modify->delete_fix("tad_event");
  modify->delete_fix("tad_revert");
  delete_event_list();

  compute_event->reset_extra_compute_fix(NULL);
}

/* ----------------------------------------------------------------------
   single short dynamics run
------------------------------------------------------------------------- */

void TAD::dynamics()
{
  update->whichflag = 1;
  update->nsteps = t_event;

  lmp->init();
  update->integrate->setup();
  // this may be needed if don't do full init
  //modify->addstep_compute_all(update->ntimestep);
  int ncalls = neighbor->ncalls;

  timer->barrier_start(TIME_LOOP);
  update->integrate->run(t_event);
  timer->barrier_stop(TIME_LOOP);
  time_dynamics += timer->array[TIME_LOOP];

  nbuild += neighbor->ncalls - ncalls;
  ndanger += neighbor->ndanger;

  update->integrate->cleanup();
  finish->end(0);
}

/* ----------------------------------------------------------------------
   quench minimization
------------------------------------------------------------------------- */

void TAD::quench()
{
  bigint ntimestep_hold = update->ntimestep;
  bigint endstep_hold = update->endstep;

  // need to change whichflag so that minimize->setup() calling 
  // modify->setup() will call fix->min_setup()

  update->whichflag = 2;
  update->nsteps = maxiter;
  update->endstep = update->laststep = update->firststep + maxiter;
  if (update->laststep < 0 || update->laststep > MAXBIGINT)
    error->all("Too many iterations");

  // full init works

  lmp->init();
  update->minimize->setup();

  // partial init does not work

  //modify->addstep_compute_all(update->ntimestep);
  //update->minimize->setup_minimal(1);

  int ncalls = neighbor->ncalls;

  timer->barrier_start(TIME_LOOP);
  update->minimize->run(maxiter);
  timer->barrier_stop(TIME_LOOP);
  time_quench += timer->array[TIME_LOOP];

  if (neighbor->ncalls == ncalls) quench_reneighbor = 0;
  else quench_reneighbor = 1;

  update->minimize->cleanup();
  finish->end(1);

  // reset timestep as if quench did not occur
  // clear timestep storage from computes, since now invalid

  update->ntimestep = ntimestep_hold;
  update->endstep = update->laststep = endstep_hold;
  for (int i = 0; i < modify->ncompute; i++)
    if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
}

/* ----------------------------------------------------------------------
   check for an event
   return 0 if no event
   return 1 if event
------------------------------------------------------------------------- */

int TAD::check_event()
{
  int flag;

  flag = 0;
  if (compute_event->compute_scalar() > 0.0) flag = 1;

  return flag;
}

/* ----------------------------------------------------------------------
   universe proc 0 prints event info
------------------------------------------------------------------------- */

void TAD::log_event(int ievent)
{
  timer->array[TIME_LOOP] = time_start;
  if (universe->me == 0) {
    double tfrac = 0.0;
    if (universe->uscreen)
      fprintf(universe->uscreen,
	      BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
	      fix_event->event_timestep,
	      timer->elapsed(TIME_LOOP),
	      fix_event->event_number,ievent,
	      "E ",
	      fix_event->ebarrier,tfrac,
	      fix_event->tlo,deltfirst);
    if (universe->ulogfile)
      fprintf(universe->ulogfile,
	      BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
	      fix_event->event_timestep,
	      timer->elapsed(TIME_LOOP),
	      fix_event->event_number,ievent,
	      "E ",
	      fix_event->ebarrier,tfrac,
	      fix_event->tlo,deltfirst);
  }

  // dump snapshot of quenched coords
  // must reneighbor and compute forces before dumping
  // addstep_compute_all insures eng/virial are calculated if needed
  
  if (output->ndump && universe->iworld == 0) {
    timer->barrier_start(TIME_LOOP);
    modify->addstep_compute_all(update->ntimestep);
    update->integrate->setup_minimal(1);
    output->write_dump(update->ntimestep);
    timer->barrier_stop(TIME_LOOP);
    time_output += timer->array[TIME_LOOP];
  }

}

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

void TAD::options(int narg, char **arg)
{
  if (narg < 0) error->all("Illegal tad command");

  // set defaults
  
  etol = 0.1;
  ftol = 0.1;
  maxiter = 40;
  maxeval = 50;

  etol_neb = 0.01;
  ftol_neb = 0.01;
  n1steps_neb = 100;
  n2steps_neb = 100;
  nevery_neb = 10;

  min_style = new char[3];
  strcpy(min_style,"cg");
  min_style_neb = new char[9];
  strcpy(min_style_neb,"quickmin");
  neb_logfilename = NULL;

  int iarg = 0;
  while (iarg < narg) {
    if (strcmp(arg[iarg],"min") == 0) {
      if (iarg+5 > narg) error->all("Illegal tad command");
      etol = atof(arg[iarg+1]);
      ftol = atof(arg[iarg+2]);
      maxiter = atoi(arg[iarg+3]);
      maxeval = atoi(arg[iarg+4]);
      if (maxiter < 0 || maxeval < 0 || 
	  etol < 0.0 || ftol < 0.0 ) 
	error->all("Illegal tad command");
      iarg += 5;

    } else if (strcmp(arg[iarg],"neb") == 0) {
      if (iarg+6 > narg) error->all("Illegal tad command");
      etol_neb = atof(arg[iarg+1]);
      ftol_neb = atof(arg[iarg+2]);
      n1steps_neb = atoi(arg[iarg+3]);
      n2steps_neb = atoi(arg[iarg+4]);
      nevery_neb = atoi(arg[iarg+5]);
      if (etol_neb < 0.0 || ftol_neb < 0.0 || 
	  n1steps_neb < 0 || n2steps_neb < 0 ||
	  nevery_neb < 0) error->all("Illegal tad command");
      iarg += 6;
      
    } else if (strcmp(arg[iarg],"min_style") == 0) {
      if (iarg+2 > narg) error->all("Illegal tad command");
      int n = strlen(arg[iarg+1]) + 1;
      delete [] min_style;
      min_style = new char[n];
      strcpy(min_style,arg[iarg+1]);
      iarg += 2;

    } else if (strcmp(arg[iarg],"neb_style") == 0) {
      if (iarg+2 > narg) error->all("Illegal tad command");
      int n = strlen(arg[iarg+1]) + 1;
      delete [] min_style_neb;
      min_style_neb = new char[n];
      strcpy(min_style_neb,arg[iarg+1]);
      iarg += 2;

    } else if (strcmp(arg[iarg],"neb_log") == 0) {
      delete [] neb_logfilename;
      if (iarg+2 > narg) error->all("Illegal tad command");
      if (strcmp(arg[iarg+1],"none") == 0) neb_logfilename = NULL;
      else {
	int n = strlen(arg[iarg+1]) + 1;
	neb_logfilename = new char[n];
	strcpy(neb_logfilename,arg[iarg+1]);
      }
      iarg += 2;
    } else error->all("Illegal tad command");
  }
}

/* ----------------------------------------------------------------------
   perform NEB calculation
------------------------------------------------------------------------- */

void TAD::perform_neb(int ievent)
{

  double **x = atom->x;
  int nlocal = atom->nlocal;

  double *buf_final;
  memory->create(buf_final,3*nlocal,"tad:buffinal");

  // set system to quenched state of event ievent

  if (universe->iworld == 0) {

    fix_event_list[ievent]->restore_event();

    int ii = 0;
    for (int i = 0; i < nlocal; i++) {
      buf_final[ii++] = x[i][0]; 
      buf_final[ii++] = x[i][1]; 
      buf_final[ii++] = x[i][2]; 
    }
  }

  MPI_Bcast(buf_final,3*nlocal,MPI_DOUBLE,universe->root_proc[0],
	    universe->uworld);

  double *buf_init;
  memory->create(buf_init,3*nlocal,"tad:bufinit");
    
  // set system to quenched state of fix_event

  if (universe->iworld == 0) {

    fix_event->restore_event();

    int ii = 0;
    for (int i = 0; i < nlocal; i++) {
      buf_init[ii++] = x[i][0]; 
      buf_init[ii++] = x[i][1]; 
      buf_init[ii++] = x[i][2]; 
    }
  }

  MPI_Bcast(buf_init,3*nlocal,MPI_DOUBLE,universe->root_proc[0],
	    universe->uworld);

  // create FixNEB object to support NEB

  int narg2 = 4;
  char **args = new char*[narg2];
  args[0] = (char *) "neb";
  args[1] = (char *) "all";
  args[2] = (char *) "neb";
  char str[128];
  args[3] = str; 
  double kspring = 1.0;
  sprintf(args[3],"%f",kspring);
  modify->add_fix(narg2,args);
  fix_neb = (Fix *) modify->fix[modify->nfix-1];
  delete [] args;

  // switch minimize style to quickmin for NEB

  narg2 = 1;
  args = new char*[narg2];
  args[0] = min_style_neb;

  update->create_minimize(narg2,args);

  delete [] args;

  // create NEB object

  neb = new NEB(lmp,etol_neb,ftol_neb,n1steps_neb, 
		n2steps_neb,nevery_neb,buf_init,buf_final);

  // free up temporary arrays

  memory->destroy(buf_init);
  memory->destroy(buf_final);

  // Run NEB

  int beginstep_hold = update->beginstep;
  int endstep_hold = update->endstep;
  int ntimestep_hold = update->ntimestep;
  int nsteps_hold = update->nsteps;

  if (universe->me == 0) {
    universe->ulogfile = ulogfile_neb;
    universe->uscreen = uscreen_neb;
  }

  // Had to bypass timer interface
  // because timer->array is reset
  // inside neb->run()

//    timer->barrier_start(TIME_LOOP);
//    neb->run();
//    timer->barrier_stop(TIME_LOOP);
//    time_neb += timer->array[TIME_LOOP];

  MPI_Barrier(world);
  double time_tmp = MPI_Wtime();
  neb->run();
  MPI_Barrier(world);
  time_neb += MPI_Wtime() - time_tmp;

  if (universe->me == 0) {
    universe->ulogfile = ulogfile_lammps;
    universe->uscreen = uscreen_lammps;
  }

  // Extract barrier energy from NEB

  if (universe->iworld == 0)
    fix_event_list[ievent]->ebarrier = neb->ebf;

  update->beginstep = update->firststep = beginstep_hold;
  update->endstep = update->laststep = endstep_hold;
  update->ntimestep = ntimestep_hold;
  update->nsteps = nsteps_hold;

  // switch minimize style back for quench

  narg2 = 1;
  args = new char*[narg2];
  args[0] = min_style;

  update->create_minimize(narg2,args);

  update->etol = etol;
  update->ftol = ftol;

  delete [] args;

  // Clean up 

  modify->delete_fix("neb");
  delete neb;
}

/* ----------------------------------------------------------------------
   check if confidence criterion for tstop is satisfied
   return 0 if not satisfied
   return 1 if satisfied
------------------------------------------------------------------------- */

int TAD::check_confidence()
{
  int flag;

  // update stopping time

  deltstop = deltconf*pow(deltfirst/deltconf, ratio_beta);

  flag = 0;
  if (deltstop < update->ntimestep - fix_event->event_timestep) flag = 1;

  return flag;
}

/* ----------------------------------------------------------------------
   reflect back in to starting state
------------------------------------------------------------------------- */

void TAD::revert()
{
  double **x = atom->x;
  double **v = atom->v;
  int *image = atom->image;
  int nlocal = atom->nlocal;

  double **array_atom = fix_revert->array_atom;

  for (int i = 0; i < nlocal; i++) {
    x[i][0] = array_atom[i][0];
    x[i][1] = array_atom[i][1];
    x[i][2] = array_atom[i][2];
    image[i] = ((int(array_atom[i][5]) + 512 & 1023) << 20) |
      ((int(array_atom[i][4]) + 512 & 1023) << 10) |
      (int(array_atom[i][3]) + 512 & 1023);
    v[i][0] = -array_atom[i][6];
    v[i][1] = -array_atom[i][7];
    v[i][2] = -array_atom[i][8];
  }
}

/* ----------------------------------------------------------------------
   Initialize list of possible events
------------------------------------------------------------------------- */

void TAD::initialize_event_list() {

  // First delete old events, if any

  delete_event_list();

  // Create new list

  n_event_list = 0;
  grow_event_list(nmin_event_list);
}

/* ----------------------------------------------------------------------
   Delete list of possible events
------------------------------------------------------------------------- */

void TAD::delete_event_list() {

  for (int i = 0; i < n_event_list; i++) {
    char str[128];
    sprintf(str,"tad_event_%d",i);
    modify->delete_fix(str);
  }
  memory->sfree(fix_event_list);
  fix_event_list = NULL;
  n_event_list = 0;
  nmax_event_list = 0;

}

/* ----------------------------------------------------------------------
   add event
------------------------------------------------------------------------- */

void TAD::add_event()
{

  // create FixEventTAD object to store possible event

  int narg = 3;
  char **args = new char*[narg];

  char str[128];
  sprintf(str,"tad_event_%d",n_event_list);

  args[0] = str;
  args[1] = (char *) "all";
  args[2] = (char *) "EVENT/TAD";
  modify->add_fix(narg,args);

  if (n_event_list == nmax_event_list) 
    grow_event_list(nmax_event_list+nmin_event_list);
  n_event_list += 1;
  int ievent = n_event_list-1;
  fix_event_list[ievent] = (FixEventTAD *) modify->fix[modify->nfix-1];

  // store quenched state for new event

  fix_event_list[ievent]->store_event_tad(update->ntimestep);

  // store hot state for new event

  fix_event->restore_state();
  fix_event_list[ievent]->store_state();

  // string clean-up

  delete [] args;

}

/* ----------------------------------------------------------------------
   compute cold time for event ievent
------------------------------------------------------------------------- */

void TAD::compute_tlo(int ievent)
{
  double deltlo,delthi,ebarrier;

  ebarrier = fix_event_list[ievent]->ebarrier;
  delthi = fix_event_list[ievent]->event_timestep 
    - fix_event->event_timestep;
  deltlo = delthi*exp(ebarrier*delta_beta);
  fix_event_list[ievent]->tlo = fix_event->tlo + deltlo;

  // update first event

  char* statstr = "D ";

  if (ievent == 0) {
    deltfirst = deltlo;
    event_first = ievent;
    statstr = "DF";
  } else if (deltlo < deltfirst) {
    deltfirst = deltlo;
    event_first = ievent;
    statstr = "DF";
  }

  // first-replica output about each event
  
  timer->array[TIME_LOOP] = time_start;
  if (universe->me == 0) {
    double tfrac = 0.0;
    if (ievent > 0) tfrac = delthi/deltstop;

    if (universe->uscreen) 
      fprintf(universe->uscreen,
	      BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
	      fix_event_list[ievent]->event_timestep,
	      timer->elapsed(TIME_LOOP),
	      fix_event->event_number,
	      ievent,statstr,ebarrier,tfrac,
	      fix_event->tlo,deltlo);

    if (universe->ulogfile) 
      fprintf(universe->ulogfile,
	      BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
	      fix_event_list[ievent]->event_timestep,
	      timer->elapsed(TIME_LOOP),
	      fix_event->event_number,
	      ievent,statstr,ebarrier,tfrac,
	      fix_event->tlo,deltlo);
  }

}

/* ----------------------------------------------------------------------
   perform event
------------------------------------------------------------------------- */

void TAD::perform_event(int ievent)
{
  // reset timestep to that of event

  update->ntimestep = fix_event_list[ievent]->event_timestep;

  // Copy event to current event
  // Should really use copy constructor for this
  fix_event->tlo = fix_event_list[ievent]->tlo;
  fix_event->ebarrier = fix_event_list[ievent]->ebarrier;
  fix_event->event_number++;
  fix_event->event_timestep = update->ntimestep;
  fix_event_list[ievent]->restore_event();
  fix_event->store_event_tad(fix_event_list[ievent]->event_timestep);

  // output stats and dump for quenched state

  log_event(ievent);

  // load and store hot state

  fix_event_list[ievent]->restore_state();
  fix_event->store_state();
}

/* ----------------------------------------------------------------------
   Allocate list of pointers to events
------------------------------------------------------------------------- */

void TAD::grow_event_list(int nmax) {
  if (nmax_event_list > nmax) return;
  fix_event_list = (FixEventTAD **) 
    memory->srealloc(fix_event_list,nmax*sizeof(FixEventTAD *),"tad:eventlist");
  nmax_event_list = nmax;
}