/*
 * helper function to calculate mean square displacement
 */

void calculateMSD(transformDiffusionInfo *info, double *box,
		  double *x, double *y, double *z, int idx1, int idx2)
{
  double xx, yy, zz;
  double delx, dely, delz;



  /* calculate distance to previous frames coordinates */
  delx = x[idx1] - info->prevx[idx2];
  dely = y[idx1] - info->prevy[idx2];
  delz = z[idx1] - info->prevz[idx2];
 
  /* if the particle moved more than half the box, assume it was imaged and adjust
     the distance of the total movement with respect to the original frame */
  if (box[0] > 0.0) {
    if (delx > box[0]/2.0)
      info->deltax[idx2] -= box[0];
    else if ( delx < -box[0]/2.0)
      info->deltax[idx2] += box[0];
    if (dely > box[1]/2.0)
      info->deltay[idx2] -= box[1];
    else if ( dely < -box[1]/2.0)
      info->deltay[idx2] += box[1];
    if (delz > box[2]/2.0)
      info->deltaz[idx2] -= box[2];
    else if (delz < -box[2]/2.0)
      info->deltaz[idx2] += box[2];
  }

  /* set the current x with reference to the un-imaged trajectory */
  xx = x[idx1] + info->deltax[idx2];
  yy = y[idx1] + info->deltay[idx2];
  zz = z[idx1] + info->deltaz[idx2];

  /* calculate the distance between this "fixed" coordinate and the
     reference (initial) frame */
  delx = xx - info->dx[idx1];
  dely = yy - info->dy[idx1];
  delz = zz - info->dz[idx1];

  /* store the distance for this atom */
  info->distancex[idx2] = delx*delx;
  info->distancey[idx2] = dely*dely;
  info->distancez[idx2] = delz*delz;
  info->distance[idx2] = delx*delx + dely*dely + delz*delz;

  info->dSum1[idx2] += xx*xx + yy*yy + zz*zz;  /* sum r^2 */
  info->dSum2[idx2] += sqrt(xx*xx + yy*yy + zz*zz);	 /* sum r */
  //info->dSum1[idx2] += delx*delx + dely*dely + delz*delz;
  //info->dSum2[idx2] += sqrt(delx*delx + dely*dely + delz*delz);

  /* update the previous coordinate set to match the current coordinates */
  info->prevx[idx2] = x[idx1];
  info->prevy[idx2] = y[idx1];
  info->prevz[idx2] = z[idx1];
}

/** ACTION ROUTINE *************************************************************
 *
 *  transformDiffusion()   --- calculate mean squared displacements vs. time
 *
 *  rewritten 2008 by Hannes Loeffler, Daresbury Laboratory, UK
 *
 ******************************************************************************/

int transformDiffusion(actionInformation *action, double *x, double *y, double *z,
		       double *box, int mode)
{
  stackType **argumentStackPointer;
  ptrajState *state;
  transformDiffusionInfo *info;

  int i, j, currentAtom, nInside;

  char *mask;

  double delx, dely, delz;
  double xx, yy, zz;
  double X[1], Y[1], Z[1], XX[2], YY[2], ZZ[2];
  double avgx, avgy, avgz, average;
  double total_mass, time, dist2, minDist, avMSD;
  double ucell[9], recip[9];

  /*
   *  USAGE:
   *
   *     diffusion [mask <mask>] [time <time per frame>] [out <file>]
   *               ([mask2 <mask>] [lower <distance>] [upper <distance>]
   *                [nwout <file>]) [avout <file>] [distances] [com]
   *
   *
   *  action argument usage:
   *
   *  mask:
   *    atoms for which the diffusion is calculated
   *  iarg1: 
   *    0 -- default, only print out averages
   *    1 -- print out un-imaged distances r(t)-r_0 for each of the
   *         active atoms
   *  carg1:
   *    the transformDiffusionInfo structure 
   *
   */



  state = (ptrajState *) action->state;

  if (mode == PTRAJ_SETUP) {
    argumentStackPointer = (stackType **) action->carg1;
    action->carg1 = NULL;

    info = safe_malloc(sizeof(transformDiffusionInfo));
    INITIALIZE_transformDiffusionInfo(info);

    /* check for mask */
    mask = argumentStackKeyToString(argumentStackPointer, "mask", NULL);

    if (mask == NULL) {
      fprintf(stderr, "WARNING in ptraj(), DIFFUSION: no mask specified. Ignoring command!\n");
      safe_free(info);
      return -1;
    } else {
      info->mask = processAtomMask(mask, state);
    }

    info->outputFilename =
      argumentStackKeyToString(argumentStackPointer, "out", "diffusion.dat");

    if (info->outputFilename == NULL) {
      fprintf(stderr, "WARNING in ptraj(), DIFFUSION: no root filename given. Ignoring command!\n");
      safe_free(info);
      return -1;
    }

    info->outputAverDist =
      argumentStackKeyToString(argumentStackPointer, "avout", NULL);

    info->timePerFrame = argumentStackKeyToDouble(argumentStackPointer,
							   "time", 1.0L);
    if (info->timePerFrame < 0) {
      fprintf(stderr, "WARNING in ptraj(), DIFFUSION: diffusion time per frame incorrectly specified. Ignoring command!\n");
      safe_free(info);
      return -1;
    }

    action->iarg1 = argumentStackContains(argumentStackPointer, "distances");
    action->iarg2 = argumentStackContains(argumentStackPointer, "com");

    mask = argumentStackKeyToString(argumentStackPointer, "mask2", NULL);

    if (mask) {
      info->mask2 = processAtomMask(mask, state);

      info->lowerCutoff = argumentStackKeyToDouble(argumentStackPointer,
							    "lower", 0.01L);
      info->upperCutoff = argumentStackKeyToDouble(argumentStackPointer,
							    "upper", 3.5L);
      info->lowerCutoff *= info->lowerCutoff;
      info->upperCutoff *= info->upperCutoff;

      info->outputNumWat =
	argumentStackKeyToString(argumentStackPointer, "nwout", "nw.dat");

      action->iarg2 = 2;
    }

    if (action->iarg2)
      action->iarg1 = 0;	/* no output of individual diffusion */

    action->carg1 = (void *) info;
    return 0;
  }


  info = (transformDiffusionInfo *) action->carg1;

  if (mode == PTRAJ_STATUS) {
    fprintf(stdout, "  DIFFUSION\n");

    fprintf(stdout, "      The atoms in mask1 are ");
    printAtomMask(info->mask, state);
    if (action->iarg2 == 1) {
      fprintf(stdout, "      Center of mass diffusion of atoms in mask1 will be computed\n");
    } else if (action->iarg2 == 2) {
      fprintf(stdout,
	      "      Atoms in mask2 in the range %7.3f to %7.3f Angstrom will be used\n",
	      sqrt(info->lowerCutoff), sqrt(info->upperCutoff));
      fprintf(stdout, "      The atoms in mask2 are ");
      printAtomMask(info->mask2, state);
    }

    if (action->iarg1 == 0) {
      fprintf(stdout, "      Only the average");
    } else {
      fprintf(stdout, "      The average and individual");
    }

    fprintf(stdout, " results will be dumped to %s\n", 
	    info->outputFilename);

    if (action->iarg2 == 2) {
      fprintf(stdout, "      The number of atoms in the shell will be dumped to %s\n",
	      info->outputNumWat);
    }

    if (info->outputAverDist) {
      fprintf(stdout, "      <dr^2> will be dumped to %s\n",
	      info->outputAverDist);
    }

    fprintf(stdout, "      The time step between frames is %7.3f ps.\n",
	    info->timePerFrame);
    fprintf(stdout, 
	    "\n      To calculate diffusion constants, calculate the slope of the lines(s)\n      and multiply by 10.0/6.0; this will give units of 1x10**-5 cm**2/s\n");

    return 0;

  } else if (mode == PTRAJ_PRINT) {
    if (info->outputAverDist) {
      for (i = 0; i < state->atoms; i++) {
 	if (info->dSum1[i] > 0.0) {
	  avMSD = (info->dSum1[i] - info->dSum2[i] * info->dSum2[i] /
		   info->elapsedFrames) / info->elapsedFrames;
	  fprintf(info->outputad, "%7i %9.3f\n", i, avMSD);
	}
      }
    }
  } else if (mode == PTRAJ_CLEANUP) {
    safe_free(info->dx);
    safe_free(info->dy);
    safe_free(info->dz);
    safe_free(info->prevx);
    safe_free(info->prevy);
    safe_free(info->prevz);
    safe_free(info->deltax);
    safe_free(info->deltay);
    safe_free(info->deltaz);
    safe_free(info->distance);
    safe_free(info->dSum1);
    safe_free(info->dSum2);
    safe_free(info->mask);
    safe_free(info->mask2);
    safe_free(info->outputFilename);
    safe_fclose(info->output);

    if (action->iarg2 == 2) {
      safe_free(info->nInside);
      safe_fclose(info->outputnw);
    }

    if (info->outputAverDist) {
      safe_fclose(info->outputad);
    }

    INITIALIZE_transformDiffusionInfo(info);
    safe_free(info);
  }


  if (mode != PTRAJ_ACTION) return 0;

  /* update local state information */
  for (i = 0; i < 6; i++)
    state->box[i] = box[i];

  info->elapsedFrames++;

  /* load up initial frame if necessary */
  if (info->dx == NULL) {
    info->dSum1 = safe_malloc(sizeof(double) * state->atoms);
    info->dSum2 = safe_malloc(sizeof(double) * state->atoms);

    for (i = 0; i < state->atoms; i++) {
      info->dSum1[i] = 0.0;
      info->dSum2[i] = 0.0;
    }

    if (action->iarg2 == 0) {  /* all */
      info->dx = safe_malloc(sizeof(double) * state->atoms);
      info->dy = safe_malloc(sizeof(double) * state->atoms);
      info->dz = safe_malloc(sizeof(double) * state->atoms);

      for (i = 0; i < state->atoms; i++) {
	if (info->mask[i])
	  info->activeAtoms++;
      }

      info->prevx = 
	safe_malloc(sizeof(double) * info->activeAtoms);
      info->prevy = 
	safe_malloc(sizeof(double) * info->activeAtoms);
      info->prevz = 
	safe_malloc(sizeof(double) * info->activeAtoms);

      info->distancex = 
	safe_malloc(sizeof(double) * info->activeAtoms);
      info->distancey = 
	safe_malloc(sizeof(double) * info->activeAtoms);
      info->distancez = 
	safe_malloc(sizeof(double) * info->activeAtoms);
      info->distance = 
	safe_malloc(sizeof(double) * info->activeAtoms);

      info->deltax = 
	safe_malloc(sizeof(double) * info->activeAtoms);
      info->deltay = 
	safe_malloc(sizeof(double) * info->activeAtoms);
      info->deltaz = 
	safe_malloc(sizeof(double) * info->activeAtoms);

      for (currentAtom = 0, i = 0; i < state->atoms; i++) {
	info->dx[i] = x[i];
	info->dy[i] = y[i];
	info->dz[i] = z[i];

	if (info->mask[i]) {
	  info->distancex[currentAtom] = 0.0;
	  info->distancey[currentAtom] = 0.0;
	  info->distancez[currentAtom] = 0.0;
	  info->distance[currentAtom] = 0.0;

	  info->deltax[currentAtom] = 0.0;
	  info->deltay[currentAtom] = 0.0;
	  info->deltaz[currentAtom] = 0.0;

	  info->prevx[currentAtom] = x[i];
	  info->prevy[currentAtom] = y[i];
	  info->prevz[currentAtom] = z[i];

	  currentAtom++;
	}
      }
    } else if (action->iarg2 == 1) { 	/* center of mass */
      info->dx = safe_malloc(sizeof(double));
      info->dy = safe_malloc(sizeof(double));
      info->dz = safe_malloc(sizeof(double));

      info->prevx = safe_malloc(sizeof(double));
      info->prevy = safe_malloc(sizeof(double));
      info->prevz = safe_malloc(sizeof(double));

      X[0] = Y[0] = Z[0] = total_mass = 0.0L;
      for (i = 0; i < state->atoms; i++) {
	if (info->mask == NULL || info->mask[i] == 1) {
	  X[0] += state->masses[i] * x[i];
	  Y[0] += state->masses[i] * y[i];
	  Z[0] += state->masses[i] * z[i];
	  total_mass += state->masses[i];
	}
      }

      info->prevx[0] = info->dx[0] = X[0] / total_mass;
      info->prevy[0] = info->dy[0] = Y[0] / total_mass;
      info->prevz[0] = info->dz[0] = Z[0] / total_mass;

      info->distancex = safe_malloc(sizeof(double));
      info->distancey = safe_malloc(sizeof(double));
      info->distancez = safe_malloc(sizeof(double));
      info->distance = safe_malloc(sizeof(double));

      info->deltax = safe_malloc(sizeof(double));
      info->deltay = safe_malloc(sizeof(double));
      info->deltaz = safe_malloc(sizeof(double));

      info->distancex[0] = 0.0;
      info->distancey[0] = 0.0;
      info->distancez[0] = 0.0;
      info->distance[0] = 0.0;

      info->deltax[0] = 0.0;
      info->deltay[0] = 0.0;
      info->deltaz[0] = 0.0;
    } else if (action->iarg2 == 2) {  /* region based */
      info->nInside = safe_malloc(sizeof(double) * state->atoms);

      info->dx = safe_malloc(sizeof(double) * state->atoms);
      info->dy = safe_malloc(sizeof(double) * state->atoms);
      info->dz = safe_malloc(sizeof(double) * state->atoms);

      info->prevx = 
	safe_malloc(sizeof(double) * state->atoms);
      info->prevy = 
	safe_malloc(sizeof(double) * state->atoms);
      info->prevz = 
	safe_malloc(sizeof(double) * state->atoms);

      info->distancex = 
	safe_malloc(sizeof(double) * state->atoms);
      info->distancey = 
	safe_malloc(sizeof(double) * state->atoms);
      info->distancez = 
	safe_malloc(sizeof(double) * state->atoms);
      info->distance = 
	safe_malloc(sizeof(double) * state->atoms);

      info->deltax = 
	safe_malloc(sizeof(double) * state->atoms);
      info->deltay = 
	safe_malloc(sizeof(double) * state->atoms);
      info->deltaz = 
	safe_malloc(sizeof(double) * state->atoms);

      for (i = 0; i < state->atoms; i++) {
	info->dx[i] = x[i];
	info->dy[i] = y[i];
	info->dz[i] = z[i];

	info->prevx[i] = x[i];
	info->prevy[i] = y[i];
	info->prevz[i] = z[i];

	info->distancex[i] = 0.0;
	info->distancey[i] = 0.0;
	info->distancez[i] = 0.0;
	info->distance[i] = 0.0;

	info->deltax[i] = 0.0;
	info->deltay[i] = 0.0;
	info->deltaz[i] = 0.0;
      }
    }

    info->elapsedFrames = 0;

    if (!openFile(&info->output, info->outputFilename, "w") ) {
      fprintf(stderr,
	      "WARNING in ptraj(), DIFFUSION: Cannot open diffusion output file %s\n",
	      info->outputFilename);
    }

    if (action->iarg2 == 2) {
      if (!openFile(&info->outputnw, info->outputNumWat, "w") ) {
	fprintf(stderr,
		"WARNING in ptraj(), DIFFUSION: Cannot open diffusion output file %s.\n",
		info->outputNumWat);
      }
    }

    if (info->outputAverDist) {
      if (!openFile(&info->outputad, info->outputAverDist, "w") ) {
	fprintf(stderr,
		"WARNING in ptraj(), DIFFUSION: Cannot open diffusion output file %s.\n",
		info->outputAverDist);
      }
    }

    return 1;
  }


  if (action->iarg2 == 0) {	/* all */
    for (currentAtom = 0, i = 0; i < state->atoms; i++) {
      if (info->mask[i] == 1) {
	calculateMSD(info, state->box, x, y, z, i, currentAtom);
	currentAtom++;
      }
    }
  } else if (action->iarg2 == 1) {  /* center of mass */
    X[0] = Y[0] = Z[0] = total_mass = 0.0L;
    for (i = 0; i < state->atoms; i++) {
      if (info->mask == NULL || info->mask[i] == 1) {
	X[0] += state->masses[i] * x[i];
	Y[0] += state->masses[i] * y[i];
	Z[0] += state->masses[i] * z[i];
	total_mass += state->masses[i];
      }
    }

    X[0] /= total_mass;
    Y[0] /= total_mass;
    Z[0] /= total_mass;

    calculateMSD(info, state->box, X, Y, Z, 0, 0);
  } else if (action->iarg2 == 2) {  /* region based */
    /* check which atoms from mask1 are inside or outside the region defined
       by lower and upper */
    for (i = 0; i < state->atoms; i++) {
      info->nInside[i] = 0;

      if (info->mask[i] == 1) {
	XX[0] =  x[i];
	YY[0] =  y[i];
	ZZ[0] =  z[i];

	minDist = info->upperCutoff;

	for (j = 0; j < state->atoms; j++) {
	  if (info->mask2[j] == 1) {
	    XX[1] =  x[j];
	    YY[1] =  y[j];
	    ZZ[1] =  z[j];

	    dist2 = calculateDistance2(0, 1, XX, YY, ZZ, box, ucell, recip, 0.0);

	    /* find minimum distance */
	    if (dist2 < minDist)
	      minDist = dist2;
	  }
	}

	/* treat only solvent molecules in the selected region */
	if (minDist > info->lowerCutoff && minDist < info->upperCutoff) {
	  info->nInside[i] = 1;
	  calculateMSD(info, state->box, x, y, z, i, i);
	}
      }
    }
  }

  time = info->elapsedFrames * info->timePerFrame;

  /* accumulate averages */
  average = 0.0L;
  avgx = 0.0L;
  avgy = 0.0L;
  avgz = 0.0L;

  if (action->iarg2 == 0) {	/* all */
    for (i = 0; i < info->activeAtoms; i++) {
      average += info->distance[i];
      avgx += info->distancex[i];
      avgy += info->distancey[i];
      avgz += info->distancez[i];
    }

    average /= (double) info->activeAtoms;
    avgx /= (double) info->activeAtoms;
    avgy /= (double) info->activeAtoms;
    avgz /= (double) info->activeAtoms;
  } else if (action->iarg2 == 1) {  /* center of mass */
    average = info->distance[0];
    avgx = info->distancex[0];
    avgy = info->distancey[0];
    avgz = info->distancez[0];
  } else if (action->iarg2 == 2) {  /* region based */
    for (nInside = 0, i = 0; i < state->atoms; i++) {
      if (info->nInside[i]) {
	average += info->distance[i];
	avgx += info->distancex[i];
	avgy += info->distancey[i];
	avgz += info->distancez[i];

	nInside += info->nInside[i];
      }
    }

    if (nInside == 0)
      error("diffusion()", "No atoms of mask 1 left for processing.  Aborting!\n");

    average /= (double) nInside;
    avgx /= (double) nInside;
    avgy /= (double) nInside;
    avgz /= (double) nInside;

    fprintf(info->outputnw, "%9.3f %7i\n", time, nInside);
  }

  /* dump output */
  fprintf(info->output, "%10.3f %10.3f %10.3f %10.3f %10.3f",
	  time, avgx, avgy, avgz, average);

  /* dump un-imaged distances if requested */
  if (action->iarg1) {
    for (i = 0; i < info->activeAtoms; i++) {
      fprintf(info->output, " %9.3f %9.3f %9.3f %9.3f",
	      info->distancex[i],
	      info->distancey[i],
	      info->distancez[i],
	      info->distance[i]);
    }
  }

  fprintf(info->output, "\n");

  return 1;
}


/*
 * helper routine to dynamically allocate memory by a fixed amount
 */

void *increaseMem(void *mem, size_t tsize, int need,
		  int *cur_size, int inc_size)
{
  void *ptr;


  if (need >= *cur_size) {
    if (*cur_size == 0) {
      if (inc_size) {
	if ( (ptr = (void *) malloc(tsize * inc_size)) == NULL)
	  error("increaseMem", "Error in allocating of %i bytes",
		tsize * inc_size);
	*cur_size = inc_size;
      } else {
	ptr = NULL;
      }
    } else {
      if ( (ptr = (void *)
	    realloc(mem, tsize * (*cur_size + inc_size))) == NULL ) {
	error("increaseMem", "Error in reallocating of %i bytes",
	      tsize * inc_size);
      } else {
	*cur_size += inc_size;
      }
    }

    mem = ptr;
    return mem;
  } else {
    return mem;
  }
}


/** ACTION ROUTINE ************************************************************
 *
 *  transformMRT()   --- perform mean residence time calculations
 *
 *  Supplementary routines:
 *    calculateDistance2 --- calculate a distance**2 and do imaging (ptraj.c)
 *    increaseMem --- dynamically increase memory by a fixed amount
 *
 *  written by Hannes H. Loeffler, 2005-2007: Laboratory of Molecular Design, ICMB,
 *                                 the University of Tokyo
 *                                 2008: Daresbury Laboratory, UK
 *
 *****************************************************************************/

int transformMRT(actionInformation *action, double *x, double *y, double *z,
		 double *box, int mode)
{
  stackType **argumentStackPointer;
  ptrajState *state;
  transformMRTInfo *info;

  FILE *outFile, *acfFile;

  int i, j, jj, k, kk, l, m, memSize, memSize2;
  int t, B, intLen, longestInterval, inside;
  int newWin, nWat, nWin, theta;
  int *atomMask, *nWater, *memCap;

  char *mask, *buffer;

  double dist2, minDist, total_mass0, total_mass1;
  double X[2], Y[2], Z[2];
  double ucell[9], recip[9];

  /*
   *  USAGE
   *
   *  mrt out <filename> [autocorr <filename> tcorr <time> toffset <time>]
   *    [solvent <mask>] [lower <dist>] [upper <dist>] [time <t>] [tstar <t>]
   *    [noimage] (siteatoms <mask> | [onemol] <sitemask> |
   *                                  <sitemask1> ... <sitemaskN>)
   *
   *  action argument usage:
   *
   *  iarg1: disable imaging?: 0 = image, 1 = noimage
   *  darg1: time interval between frames
   *  carg1: a transformMRTInfo structure (stored in actionInformation for
   *         repeated calls)
   */


  state = action->state;

  if (mode == PTRAJ_SETUP) {
    argumentStackPointer = (stackType **) action->carg1;
    action->carg1 = NULL;

    info = (transformMRTInfo *) safe_malloc(sizeof(transformMRTInfo));
    INITIALIZE_transformMRTInfo(info);
    info->allocMem += sizeof(transformMRTInfo);

    /*
     * pick arguments from the stack
     */

    /* imaging? */
    action->iarg1 = argumentStackContains(argumentStackPointer, "noimage");

    /* output file */
    info->filename =
      argumentStackKeyToString(argumentStackPointer, "out", NULL);
    info->allocMem += sizeof(info->filename);
    
    if (info->filename == NULL) {
      fprintf(stderr, "WARNING in ptraj(), MRT: no output file specified. Ignoring command!\n");
      safe_free(info);
      return -1;
    }

    /* check for solvent mask */
    mask = argumentStackKeyToString(argumentStackPointer, "solvent", NULL);

    if (mask == NULL)
      info->solventMask = state->solventMask;  /* default solvent */
    else
      info->solventMask = processAtomMask(mask, state);

    if (info->solventMask == NULL) {
      fprintf(stderr, "WARNING in ptraj(), MRT: no solvent mask specified. Ignoring command!\n");
      safe_free(info);
      return -1;
    }

    info->allocMem += sizeof(info->solventMask);

    /* time resolution (dt) in trajectory */
    action->darg1 = argumentStackKeyToDouble(argumentStackPointer, "time", 1.0L);

    /* time water can be inside/outside without counting it as left/entered */
    info->nStar = (int) floor(
      argumentStackKeyToDouble(argumentStackPointer, "tstar", 0.0L) /
      action->darg1 + 0.5);

    /* lower and upper distance limits */
    info->lowerCutoff =		/* default prevents binning of 0 distances */
      argumentStackKeyToDouble(argumentStackPointer, "lower", 0.01L);

    info->upperCutoff =		/* upper limit approx. RDF minimum */
      argumentStackKeyToDouble(argumentStackPointer, "upper", 3.5L);

    /* store squared values */
    info->lowerCutoff *= info->lowerCutoff;
    info->upperCutoff *= info->upperCutoff;

    /* calculate autocorrelation function too? */
    info->autoCorr =
      argumentStackKeyToString(argumentStackPointer, "autocorr", NULL);

    /* autocorrelation parameters */
    info->wSize = (int) floor(
      argumentStackKeyToDouble(argumentStackPointer, "tcorr", 400.0L) /
      action->darg1 + 0.5);

    info->nOffset = (int) floor(
      argumentStackKeyToDouble(argumentStackPointer, "toffset", 10.0L) /
      action->darg1 + 0.5);

    if (info->nOffset < 1 || info->nOffset > info->wSize) {
      fprintf(stderr, "WARNING in ptraj(), MRT: toffset must be in the range from %8.3f to %8.3f. Ignoring command!\n", action->darg1, info->wSize * action->darg1);
      safe_free(info);
      return -1;
    }

    if ( (int) fmod(info->wSize, info->nOffset) != 0) {
      fprintf(stderr, "WARNING in ptraj(), MRT: tcorr must be multiple of toffset. Ignoring command!\n");
      safe_free(info);
      return -1;
    }

    nWin = info->wSize / info->nOffset;
    info->idxMaxWin = nWin - 1;

    /* treat the siteatoms as a single entity? */
    action->iarg2 = argumentStackContains(argumentStackPointer, "onemol");

    /* process solvent mask */
    mask = argumentStackKeyToString(argumentStackPointer, "siteatoms", NULL);
    memSize = memSize2 = 0;

    if (mask != NULL) {
      atomMask = processAtomMask(mask, state);

      if (atomMask == NULL) {
	fprintf(stderr, "WARNING in ptraj(), MRT: no solvent mask specified. Ignoring command!\n");
	safe_free(info);
	return -1;
      }

      for (i = 0, j = 0; j < state->atoms; j++) {
	if (atomMask[j] > 0) {
	  info->soluteMask = (int **)
	    increaseMem(info->soluteMask, sizeof(int), i, &memSize, MRT_INCSIZE);

	  /* store only indexes of "active" atoms */
	  info->soluteMask[i] = (int *) safe_malloc(sizeof(int) * 2);
	  memSize2 += 2 * sizeof(int);

	  info->soluteMask[i][0] = j;
	  info->soluteMask[i][1] = -1;
	  i++;
	}
      }
    } else {	     /* try to read rest of command line as individual masks */
      for (i = 0;
	   (buffer = getArgumentString(argumentStackPointer, NULL)) != NULL;
	   i++) {
	info->soluteMask = (int **)
	  increaseMem(info->soluteMask, sizeof(int), i, &memSize, MRT_INCSIZE);
	memCap = (int *)
	  increaseMem(memCap, sizeof(int), i, &memSize2, MRT_INCSIZE);

	atomMask = processAtomMask(buffer, state);

	if (atomMask == NULL) {
	  fprintf(stderr, "WARNING in ptraj(), MRT: error in solute mask. Ignoring command!\n");
	  safe_free(info->soluteMask);
	  safe_free(info);
	  return -1;
	}

	j = memCap[i] = 0;

	for (k = 0; k < state->atoms; k++) {
	  if (atomMask[k] > 0) {
	    info->soluteMask[i] = (int *)
	      increaseMem(info->soluteMask[i],
			  sizeof(int), j, &memCap[i], MRT_INCSIZE);
 
	    /* store only indexes of "active" atoms */
	    info->soluteMask[i][j] = k;
	    j++;
	  }
	}

	info->soluteMask[i] = (int *)
	  increaseMem(info->soluteMask[i],
		      sizeof(int), j, &memCap[i], MRT_INCSIZE);

	info->soluteMask[i][j] = -1;
	info->allocMem += memCap[i];
      }

      if (i == 1 && action->iarg2) {
	info->solventMask2 = (int *) safe_malloc(sizeof(int) *
						 state->solventMolecules);
	info->allocMem += sizeof(info->solventMask2);
      } else {
	action->iarg2 = 0;
      }
    }  

    if ( (info->numberSoluteSites = i ) < 1) {
      fprintf(stderr, "WARNING in ptraj(), MRT: no solute mask(s) specified. Ignoring command!\n");
      safe_free(info);
      return -1;
    }

    info->allocMem += memSize + memSize2;

    /*
     * allocate memory
     */

    /*** setup memory for inside/outside counts */
    /* FIXME: allocate memory only for inside waters, elminate dependence
       on state->solventMolecules */
    info->in_ = (int *)
      safe_malloc(sizeof(int) * info->numberSoluteSites * state->solventMolecules);
    info->out_ = (int *)
      safe_malloc(sizeof(int) * info->numberSoluteSites * state->solventMolecules);
    info->allocMem += 2 * sizeof(int) * info->numberSoluteSites * state->solventMolecules;

    /* allocate memory for pointers to each row */
    info->inCount = (int **)
      safe_malloc(sizeof(int) * info->numberSoluteSites);
    info->outCount = (int **)
      safe_malloc(sizeof(int) * info->numberSoluteSites);

    /* point each row pointer to its row */
    for (i = 0; i < info->numberSoluteSites; i++) {
      info->inCount[i] = info->in_ + (i * state->solventMolecules);
      info->outCount[i] = info->out_ + (i * state->solventMolecules);
    }
    /*** end setup */


    info->tCounter = (int *)
      safe_malloc(sizeof(int) * info->numberSoluteSites);

    for (i = 0; i < info->numberSoluteSites; i++)
      info->tCounter[i] = -1;

    info->memCap1 = (int *)
      safe_malloc(sizeof(int) * info->numberSoluteSites);

    /* preallocate memory for first dimension */
    info->timeInside = (int **)
      safe_malloc(sizeof(int) * info->numberSoluteSites);
    info->allocMem += 5 * sizeof(int) * info->numberSoluteSites;


    if (info->autoCorr) {
      /* preallocate memory for the first two dimensions */
      info->nw_ = (int **)
	safe_malloc(sizeof(int) * info->numberSoluteSites * nWin );


      info->mem2_ = (int *)
	safe_malloc(sizeof(int) * info->numberSoluteSites * nWin);

      info->nwc_ = (int *)
	safe_malloc(sizeof(int) * info->numberSoluteSites * nWin);
      info->allocMem += 3 * sizeof(int) * info->numberSoluteSites * nWin;

      info->nVisits_ = (int *)
	safe_malloc(sizeof(int) * info->numberSoluteSites * info->wSize);

      info->acf_ = (double *)
	safe_malloc(sizeof(double) * info->numberSoluteSites * info->wSize);
      info->allocMem += 2 * sizeof(int) * info->numberSoluteSites * info->wSize;

      info->nWater0 = (int ***)
	safe_malloc(sizeof(int) * info->numberSoluteSites);

      info->nwCounter0 = (int **)
	safe_malloc(sizeof(int) * info->numberSoluteSites);

      info->memCap2 = (int **)
	safe_malloc(sizeof(int) * info->numberSoluteSites);

      info->nVisits = (int **)
	safe_malloc(sizeof(int) * info->numberSoluteSites);
      info->allocMem += 4 * sizeof(int) * info->numberSoluteSites;

      info->acf = (double **)
	safe_malloc(sizeof(double) * info->numberSoluteSites);
      info->allocMem += sizeof(double) * info->numberSoluteSites;

      for (i = 0; i < info->numberSoluteSites; i++) {
	info->memCap2[i] = info->mem2_ + i * nWin;
	info->nwCounter0[i] = info->nwc_ + i * nWin;
	info->nWater0[i] = info->nw_ + i * nWin;
	info->nVisits[i] = info->nVisits_ + i * info->wSize;
	info->acf[i] = info->acf_ + i * info->wSize;
      }
    }

    /* store info for later calls to this function */
    action->carg1 = (void *) info;

    return 0;
  }

  info = (transformMRTInfo *) action->carg1;

  if (mode == PTRAJ_STATUS) {
    fprintf(stdout,
            "  MRT: mean residence time calculation\n");
    fprintf(stdout,
            "      Solvent molecules considered in the range %6.2f to %6.2f Angstrom\n",
            sqrt(info->lowerCutoff), sqrt(info->upperCutoff) );
    fprintf(stdout, "      Solvent atom selection is "),
    printAtomMask(info->solventMask, state);

    if (action->iarg2)
      fprintf(stdout, "      Solute mask ist treated as single molecule\n");

    fprintf(stdout, "      Solute atom selection is (%i site%s)\n",
	    info->numberSoluteSites,
	    info->numberSoluteSites > 1 ? "s" : "");

    for (i = 0; i < info->numberSoluteSites; i++) {
      fprintf(stdout, "        %5i ", i+1);

      k = 0;
      atomMask = (int *) safe_malloc(sizeof(int) * state->atoms);

      /* reconstruct mask from index list */
      for (j = 0; j < state->atoms; j++) {
	if (info->soluteMask[i][k] == j) {
	  atomMask[j] = 1;
	  k++;
	} else {
	  atomMask[j] = 0;
	}
      }

      printAtomMask(atomMask, state);

      safe_free(atomMask);
    }

    fprintf(stdout, "      The time step between frames is %7.3f ps\n", action->darg1);
    fprintf(stdout, "      Survival function: output to %s\n", info->filename);
    fprintf(stdout, "        Water molecules may leave/enter for %7.3f ps\n",
	    info->nStar * action->darg1);

    if (info->autoCorr) {
      fprintf(stdout, "      Autocorrelation function: output to %s\n",
	      info->autoCorr);
      fprintf(stdout, "        Window correlation length is %8.3f ps\n",
	      info->wSize * action->darg1);
      fprintf(stdout, "        Time offset between windows is %8.3f ps\n",
	      info->nOffset * action->darg1);
    }

    if (action->iarg1)
      fprintf(stdout, "      Imaging is disabled\n");

    fprintf(stdout, "      Memory allocated so far: %8.2f MB\n",
            (double) info->allocMem / TO_MB);

  } else if (mode == PTRAJ_PRINT) {

    if (openFile(&outFile, info->filename, "w") == 0) {
      fprintf(stderr, "WARNING in ptraj(), MRT: couldn't open output file %s.\n",
              info->filename);
      return -1;
    }

    fprintf(stdout, "PTRAJ MRT: dumping correlation data to output file%s\n",
	    info->autoCorr ? "s" : "");

    /* flush counters */
    for (i = 0; i < info->numberSoluteSites; i++) {
      for (j = 0; j < state->solventMolecules; j++) {
        if (info->inCount[i][j] > 0) {
	  info->timeInside[i] = (int *)
	    increaseMem(info->timeInside[i],
			sizeof(int), ++info->tCounter[i], &info->memCap1[i],
			MRT_INCSIZE);

          info->timeInside[i][info->tCounter[i]] =
            info->inCount[i][j];
        }
      }
    }

    if (prnlev > 5) {
      for (i = 0; i < info->numberSoluteSites; i++) {
        fprintf(stdout, "site #%i (len=%i):", i+1, info->tCounter[i]);

        for (j = 0; j < info->tCounter[i]; j++)
          fprintf(stdout, " %i", info->timeInside[i][j]);

        fprintf(stdout, "\n");
      }
    }

    /* determine longest interval */
    longestInterval = 0;

    for (i = 0; i < info->numberSoluteSites; i++)
      for (j = 0; j < info->tCounter[i]; j++)
        if (info->timeInside[i][j] > longestInterval)
          longestInterval = info->timeInside[i][j];


    /*
     * calculate mean residence time  for each site;
     * see RW Impey, PA Madden and IR McDonald, JPhysChem 87(1983), 5071
     * for algorithm see eqs. 2, 3 and Appendix in
     * AE Garcia and L Stiller, JComputChem 14(1993), 1396
     */

    fprintf(outFile, "# survival function\n");

    for (t = 0; t < longestInterval; t++) {  /* time loop */
      fprintf(outFile, "%9.3f", t * action->darg1);

      for (i = 0; i < info->numberSoluteSites; i++) {
        B = 0;

        for (j = 0; j < info->tCounter[i]; j++) {
          intLen = info->timeInside[i][j];

          if (intLen > t && intLen > info->nStar)
            B += intLen - t;    /* (intLen - t) equals number of substrings */
        }

        fprintf(outFile, " %7.4f",
                (double) B / ( (double) state->maxFrames - (double) t));
      }

      fprintf(outFile, "\n");
    }

    safe_fclose(outFile);


    if (info->autoCorr) {
      /*
       * write autocorrelation function
       */

      if (openFile(&acfFile, info->autoCorr, "w") == 0) {
	fprintf(stderr, "WARNING in ptraj(), MRT: couldn't open output file %s.\n",
		info->autoCorr);
	return -1;
      }

      fprintf(acfFile, "# autocorrelation function\n");

      for (t = 0; t < info->wSize; t++) {  /* time loop */
	fprintf(acfFile, "%9.3f", t * action->darg1);

	for (i = 0; i < info->numberSoluteSites; i++) {
	  /* print average value */
	  fprintf(acfFile, " %7.4f", info->acf[i][t] /
		  (double) info->nVisits[i][t]);
	}

	fprintf(acfFile, "\n");
      }

      safe_fclose(acfFile);
    }

    fprintf(stdout, "PTRAJ MRT: %8.2f MB of memory allocated\n",
            (double) info->allocMem / TO_MB);

    memSize = info->memCap0;
    for (i = 0; i < info->numberSoluteSites; i++) {
      memSize += info->memCap1[i];

      if (info->autoCorr) {
	for (j = 0; j < info->idxMaxWin + 1; j++)
	  memSize += info->memCap2[i][j];
      }
    }

    fprintf(stdout, "PTRAJ MRT: %8.2f MB of memory allocated dynamically\n",
            (double) memSize / TO_MB);

    fprintf(stdout, "PTRAJ MRT: %8.2f MB total memory allocated\n",
            (double) (info->allocMem + memSize) / TO_MB);

  } else if (mode == PTRAJ_CLEANUP) {
    safe_free(info->memCap1);

    safe_free(info->in_);
    safe_free(info->out_);
    safe_free(info->inCount);
    safe_free(info->outCount);
    safe_free(info->tCounter);

    for (i = 0; i < info->numberSoluteSites; i++) {
      safe_free(info->soluteMask[i]);
      safe_free(info->timeInside[i]);
    }

    safe_free(info->soluteMask);

    if (action->iarg2)
      safe_free(info->solventMask2);

    safe_free(info->timeInside);

    if (info->autoCorr) {
      safe_free(info->mem2_);
      safe_free(info->nwc_);
      safe_free(info->nw_);
      safe_free(info->nVisits_);
      safe_free(info->acf_);
      safe_free(info->memCap2);
      safe_free(info->nWater);
      safe_free(info->nWater0);
      safe_free(info->nwCounter0);
      safe_free(info->nVisits);
      safe_free(info->acf);
    }

    safe_free(info->filename);
    safe_free(info);
  }

  if (mode != PTRAJ_ACTION) return 0;

  /* update local state information */
  for (i = 0; i < 6; i++)
    state->box[i] = box[i];

  /* set up information for imaging non-orthorhombic cells if necessary */
  if (action->iarg1 == 0 &&
      (box[3] != 90.0 || box[4] != 90.0 || box[5] != 90.0))
    boxToRecip(box, ucell, recip);


  /*
   * calculate the mean residence time based on the survival function
   * (RW Impey, PA Madden and IR McDonald, JPhysChem 87(1983), 5071)
   * and the autocorrelation function (optional)
   */

  if (info->autoCorr) {
    t = info->timeStep;
    newWin = ( (int) fmod(t, info->nOffset) == 0 );

    /* manage window counters */
    if (newWin) {
      if (++info->ntor > info->idxMaxWin) /* reuse windows */
	info->ntor = 0;		/* equals (int) fmod(t, info->idxMaxWin + 1) */

      if (info->ncWin <= info->idxMaxWin)  /* current # of parallel windows */
	info->ncWin++;
    }
  }


  /* iterate over solute sites */
  for (i = 0; i < info->numberSoluteSites; i++) {  /* i loop */

    /* reset counters for acf calculation */
    if (info->autoCorr) {
      nWat = -1;

      if (newWin)
	info->nwCounter0[i][info->ntor] = -1;
    }

    if (!action->iarg2) {
      X[0] = Y[0] = Z[0] = total_mass0 = 0.0L;

      /* calculate center-of-mass of active atoms */
      for (j = 0; (k = info->soluteMask[i][j]) >= 0; j++) {
	X[0] += state->masses[k] * x[k];
	Y[0] += state->masses[k] * y[k];
	Z[0] += state->masses[k] * z[k];
	total_mass0 += state->masses[k];
      }

      X[0] = X[0] / total_mass0;
      Y[0] = Y[0] / total_mass0;
      Z[0] = Z[0] / total_mass0;

      atomMask = info->solventMask;
    } else { /* sitemask is single molecule */

      /* filter water molecules according to selected region */
      for (k = 0; k < state->solventMolecules; k++) {
	info->solventMask2[k] = 0;
	X[1] = Y[1] = Z[1] = total_mass1 = 0.0L;

	/* calculate center-of-mass */
	for (kk = state->solventMoleculeStart[k];
	     kk < state->solventMoleculeStop[k];
	     kk++) {
	  if (info->solventMask[kk] > 0) {
	    X[1] += state->masses[kk] * x[kk];
	    Y[1] += state->masses[kk] * y[kk];
	    Z[1] += state->masses[kk] * z[kk];
	    total_mass1 += state->masses[kk];
	  }
	}

	X[1] = X[1] / total_mass1;
	Y[1] = Y[1] / total_mass1;
	Z[1] = Z[1] / total_mass1;

	minDist = info->upperCutoff;

	for (j = 0; (l = info->soluteMask[i][j]) >= 0; j++) {
	  X[0] = x[l];
	  Y[0] = y[l];
	  Z[0] = z[l];

	  /* calculate squared distance (imaged if necessary) */
	  dist2 = calculateDistance2(0, 1, X, Y, Z, box, ucell, recip, 0.0);

	  /* find minimum distance */
	  if (dist2 < minDist)
	    minDist = dist2;

          /* the following code is faster but works only for nearest waters */
/*           if (dist2 > info->lowerCutoff && dist2 < info->upperCutoff) { */
/*             info->solventMask2[k] = 1; */
/*             break; */
/*           } */
	}

	/* store only solvent molecules in the selected region */
	if (minDist > info->lowerCutoff && minDist < info->upperCutoff) {
	  info->solventMask2[k] = 1;
	}
      }

      atomMask = info->solventMask2;
    }


    /* iterate over solvent molecules */
    for (j = 0; j < state->solventMolecules; j++) {  /* j loop */

      if (!action->iarg2) {	/* distance based criterion */
        X[1] = Y[1] = Z[1] = total_mass1 = 0.0L;

        /* calculate center-of-mass */
        for (jj = state->solventMoleculeStart[j];
             jj < state->solventMoleculeStop[j];
             jj++) {
          if (atomMask[jj] > 0) {
            X[1] += state->masses[jj] * x[jj];
            Y[1] += state->masses[jj] * y[jj];
            Z[1] += state->masses[jj] * z[jj];
            total_mass1 += state->masses[jj];
          }
        }

        X[1] = X[1] / total_mass1;
        Y[1] = Y[1] / total_mass1;
        Z[1] = Z[1] / total_mass1;

	/* calculate squared distance (imaged if necessary) */
	dist2 = calculateDistance2(0, 1, X, Y, Z, box, ucell, recip, 0.0);
	inside = (dist2 > info->lowerCutoff && dist2 < info->upperCutoff);
      } else {			/* predefined region */
	inside = atomMask[j];
      }

      /* this code works for t* = 0 */
      /* CHECK: t* > 0 */
      if (inside) {  /* inside */
        if (info->outCount[i][j] <= info->nStar) {
          info->inCount[i][j] += info->outCount[i][j] + 1;
        } else {
	  info->timeInside[i] = (int *)
	    increaseMem(info->timeInside[i],
			sizeof(int), ++info->tCounter[i], &info->memCap1[i],
			MRT_INCSIZE);

	  /* store the time a water spent inside */
          info->timeInside[i][info->tCounter[i]] = info->inCount[i][j];
          info->inCount[i][j] = 1;  /* reset in-counter; we are still inside */
        }

        info->outCount[i][j] = 0; /* reset out-counter */

	if (info->autoCorr) {
	  k = info->ntor;

	  /* new time origin every nOffset steps */
	  if (newWin) {
	    info->nWater0[i][k] = (int *)
	      increaseMem(info->nWater0[i][k], sizeof(int),
			  ++info->nwCounter0[i][k], &info->memCap2[i][k],
			  MRT_INCSIZE2);

	    info->nWater0[i][k][info->nwCounter0[i][k]] = j;
	  }

          /* current time step and site */
	  info->nWater = (int *)
	    increaseMem(info->nWater, sizeof(int), ++nWat,
			&info->memCap0, MRT_INCSIZE2);

          info->nWater[nWat] = j;
	}
      } else {                  /* outside */
	if (info->inCount[i][j] > 0)  /* ignore leading zeros */
	  info->outCount[i][j]++;  /* store the time a water spent outside */
      }
    }

    if (info->autoCorr) {
      /* 
       * compute auto-correlation function on-the-fly
       * use a sliding window approach to improve statistics
       */

      /* iterate over all parallel windows */
      for (k = 0; k < info->ncWin; k++) {
	theta = 0;

	/* compare water indexes of current time step to each time origin */
        for (l = 0; l <= info->nwCounter0[i][k]; l++)
          for (m = 0; m <= nWat; m++)
            if (info->nWater0[i][k][l] == info->nWater[m])
              theta++;

	l = (int) fmod(t - k * info->nOffset, info->wSize);

	info->acf[i][l] += (double) theta;
	info->nVisits[i][l]++;
      }
    }
  }

  info->timeStep++;

  return 1;
}


/** ACTION ROUTINE ************************************************************
 *
 *  transformImgDist()   --- calculate distances to all neighbour images and
 *                           reported those below a user-specified threshold
 *
 *  Supplementary routines:
 *
 *    
 *  written by Hannes H. Loeffler, Laboratory of Molecular Design, ICMB,
 *                                 the University of Tokyo
 *
 *****************************************************************************/

int transformImgDist(actionInformation *action, double *x, double *y, double *z, 
		     double *box, int mode)
{
  stackType **argumentStackPointer;
  ptrajState *state;
  transformImgDistInfo *info;

  int i, ii, j, jj, ix, iy, iz;

  char *mask;

  double xim, yim, zim;
  double dx, dy, dz, dist;

  /*
   *  USAGE
   *
   *  imgdist mask <mask> threshold <distance>
   *
   *  action argument usage:
   *
   *  carg1: a transformImgDist structure (stored in actionInformation for
   *         repeated calls)
   */


  state = action->state;

  if (mode == PTRAJ_SETUP) {
    argumentStackPointer = (stackType **) action->carg1;
    action->carg1 = NULL;

    info = (transformImgDistInfo *) safe_malloc(sizeof(transformImgDistInfo));
    INITIALIZE_transformImgDistInfo(info);


    /*
     * pick arguments from the stack
     */

    mask = argumentStackKeyToString(argumentStackPointer, "mask", NULL);

    if (mask == NULL) {
      fprintf(stderr, "WARNING in ptraj(), DISTIMG: no mask specified. Ignoring command!\n");
      safe_free(info);
      return -1;
    } else {
      info->mask = processAtomMask(mask, state);
      info->activeAtoms = (int *) safe_malloc(sizeof(int) * state->atoms);

      j = 0;
      for (i = 0; i < state->atoms; i++) {
	if (info->mask[i] > 0) {
	  info->activeAtoms[j] = i;
	  j++;
	}
      }

      info->numberActiveAtoms = j;
    }

    /* threshold */
    info->threshold =
      argumentStackKeyToDouble(argumentStackPointer, "threshold", 5.0L);
    info->threshold *= info->threshold;

    /* output file */
    info->filename =
      argumentStackKeyToString(argumentStackPointer, "out", NULL);

    if (info->filename == NULL) {
      fprintf(stderr, "WARNING in ptraj(), IMGDIST: no output file specified. Ignoring command!\n");
      safe_free(info);
      return -1;
    }

    /* store info for later calls to this function */
    action->carg1 = (void *) info;

    return 0;
  }

  info = (transformImgDistInfo *) action->carg1;

  if (mode == PTRAJ_STATUS) {
    fprintf(stdout, "  IMGDIST: distances to images calculated for mask ");
    printAtomMask(info->mask, state);
    fprintf(stdout,
            "      Distances shorter than %6.2f Angstrom written to %s\n",
	    sqrt(info->threshold), info->filename);

    if (openFile(&info->outFile, info->filename, "w") == 0) {
      fprintf(stderr, "WARNING in ptraj(), IMGDIST: couldn't open output file %s.\n",
              info->filename);
      safe_free(info);
      return -1;
    }
  } else if (mode == PTRAJ_PRINT) {
    fprintf(stdout, "PTRAJ IMGDIST: smallest distance is %6.2f\n", sqrt(info->min) );
  } else if (mode == PTRAJ_CLEANUP) {
    safe_free(info);
  }

  if (mode != PTRAJ_ACTION) return 0;

  /* update local state information */
  for (i = 0; i < 6; i++)
    state->box[i] = box[i];

  /*
   * find distances to images smaller than the threshold
   * assume orthorombic cell
   */

  for (ix = -1; ix <= 0; ix++) {  /* we only need images in one direction */
    for (iy = -1; iy <= 0; iy++) {
      for (iz = -1; iz <= 0; iz++) {
	if ( !(ix == 0 && iy == 0 && iz == 0) ) {  /* skip self-image */
	  for (i = 0; i < info->numberActiveAtoms; i++) {
	    ii = info->activeAtoms[i];

	    /* create image */
	    xim = x[ii] + state->box[0]*ix;
	    yim = y[ii] + state->box[1]*iy;
	    zim = z[ii] + state->box[2]*iz;

	    for (j = i + 1; j < info->numberActiveAtoms; j++) {
	      jj = info->activeAtoms[j];

	      dx = xim - x[jj];
	      dy = yim - y[jj];
	      dz = zim - z[jj];

	      dist = dx*dx + dy*dy + dz*dz;

	      /* write distances shorter than the threshold to the output file */
	      if (dist < info->threshold)
		fprintf(info->outFile, "%i: %f\n",
			info->timeStep, sqrt(dist) );

	      /* remember shortest distance */
	      if (dist < info->min)
		info->min = dist;
	    }
	  }
	}
      }
    }
  }

  info->timeStep++;

  return 1;
}