Particle *get_mol_com_particle(Particle *calling_p){
   int mol_id;
   int i;
   Particle *p;

   mol_id=calling_p->p.mol_id;

   if (mol_id < 0) {
     ostringstream msg;
     msg <<"Particle does not have a mol id! pnr= " << calling_p->p.identity << "\n";
     runtimeError(msg);
     return NULL;
   }
   for (i=0;i<topology[mol_id].part.n;i++){
      p=local_particles[topology[mol_id].part.e[i]];

      if (p==NULL){
          ostringstream msg;
          msg <<"Particle does not exist in put_mol_force_on_parts! id= " << topology[mol_id].part.e[i] << "\n";
          runtimeError(msg);
         return NULL;
      }

      if (ifParticleIsVirtual(p)) {
          return p;
       }
   }

   ostringstream msg;
   msg <<"No com found in get_mol_com_particleParticle does not exist in put_mol_force_on_parts! pnr= " << calling_p->p.identity << "\n";
   runtimeError(msg);
   return NULL;

   return calling_p;
}

void lb_GPU_sanity_checks()
{
  if(this_node == 0){
    if (lbpar_gpu.agrid < 0.0) {
        ostringstream msg;
        msg <<"Lattice Boltzmann agrid not set";
        runtimeError(msg);
    }
    if (lbpar_gpu.tau < 0.0) {
        ostringstream msg;
        msg <<"Lattice Boltzmann time step not set";
        runtimeError(msg);
    }
    for(int i=0;i<LB_COMPONENTS;i++){
      if (lbpar_gpu.rho[0] < 0.0) {
          ostringstream msg;
          msg <<"Lattice Boltzmann fluid density not set";
          runtimeError(msg);
      }
      if (lbpar_gpu.viscosity[0] < 0.0) {
          ostringstream msg;
          msg <<"Lattice Boltzmann fluid viscosity not set";
          runtimeError(msg);
      }
    }
  }
}

int getintersection(double pos1[3], double pos2[3],int given, int get, double value, double *answer, double box_size[3])
{
  /*pos1 and pos2 are two particle positions.                                                  */
  /*given and get are integers from 0 to 2. 0 = x direction. 1 = y direction. 2 = z direction  */
  /*there is a point on the line between the two particles p1 and p2 such that r[given]=value  */
  /*this procedure returns the value of r[get] at that point                                   */

  double p2r[3];
  int i;

  for (i=0;i<3;i++) {
    p2r[i] = drem_down((pos2[i]-pos1[i])+box_size[i]/2.0,box_size[i])-box_size[i]/2.0;
  }
  value = drem_down((value-pos1[given])+box_size[given]/2.0,box_size[given])-box_size[given]/2.0;
  //PTENSOR_TRACE(fprintf(stderr,"%d: getintersection: p1 is %f %f %f p2 is %f %f %f p2r is %f %f %f newvalue is %f\n",this_node,pos1[0],pos1[1],pos1[2],pos2[0],pos2[1],pos2[2],p2r[0],p2r[1],p2r[2],value););
  
  if ((value)*(p2r[given]) < -0.0001) {
      ostringstream msg;
      msg <<"analyze stress_profile: getintersection: intersection is not between the two given particles - " << value << " is not between " << 0.0 << " and " << p2r[given] << " and box size is " << box_size[given] << ", given is " << given << "\n";
      runtimeError(msg);
    return 0; 
  } else if (given == get) {
    *answer =  drem_down(value + pos1[given],box_size[given]);;
  } else if (0==p2r[given]) {
      ostringstream msg;
      msg <<"analyze stress_profile: getintersection: intersection is a line, not a point - value is " << value << " same as " << 0.0 << " and " << p2r[given] << "\n";
      runtimeError(msg);
    return 0;   
  } else {
    *answer =  drem_down(pos1[get]+p2r[get]/p2r[given]*value,box_size[get]);
  }
  return 1;
}

void put_mol_force_on_parts(Particle *p_com){
   int i,j,mol_id;
   Particle *p;
   double force[3],M;
#ifdef VIRTUAL_SITES_DEBUG
   int count=0;
#endif
  mol_id=p_com->p.mol_id;
   for (i=0;i<3;i++){
      force[i]=p_com->f.f[i];
      p_com->f.f[i]=0.0;
   }
#ifdef MASS
   M=0;
   for (i=0;i<topology[mol_id].part.n;i++){
      p=local_particles[topology[mol_id].part.e[i]];
      #ifdef VIRTUAL_SITES_DEBUG
      if (p==NULL){
          ostringstream msg;
          msg <<"Particle does not exist in put_mol_force_on_parts! id= " << topology[mol_id].part.e[i] << "\n";
          runtimeError(msg);
         return;
      }
      #endif
       if (ifParticleIsVirtual(p)) continue;
      M+=PMASS(*p);
   }
#else
   M=topology[mol_id].part.n-1;
#endif
   for (i=0;i<topology[mol_id].part.n;i++){
      p=local_particles[topology[mol_id].part.e[i]];
      #ifdef VIRTUAL_SITES_DEBUG
      if (p==NULL){
          ostringstream msg;
          msg <<"Particle does not exist in put_mol_force_on_parts! id= " << topology[mol_id].part.e[i] << "\n";
          runtimeError(msg);
         return;
      }
      #endif
      if (!ifParticleIsVirtual(p)) {
         for (j=0;j<3;j++){
            p->f.f[j]+=PMASS(*p)*force[j]/M;
         }
#ifdef VIRTUAL_SITES_DEBUG
         count++;
#endif
      }
   }
#ifdef VIRTUAL_SITES_DEBUG
   if (count!=topology[mol_id].part.n-1){
       ostringstream msg;
       msg <<"There is more than one COM input_mol_force_on_parts! mol_id= " << mol_id << "\n";
       runtimeError(msg);
      return;
   }
#endif
}

/* new version for new topology structure */
int analyze_fold_molecules(float *coord, double shift[3])
{
  int m,p,i, tmp;
  int mol_size, ind;
  double cm_tmp, com[3];

  #ifdef LEES_EDWARDS
  if(lees_edwards_offset != 0.0 ){
    fprintf(stderr, "Error: Folding molecules not supported under Lees-Edwards.\n");
    exit(8);
  }
  #endif
  
  /* check molecule information */
  if ( n_molecules < 0 ) return ES_ERROR;

  if (!sortPartCfg()) {
      ostringstream msg;
      msg <<"analyze_fold_molecules: could not sort particle config, particle ids not consecutive?";
      runtimeError(msg);
    return ES_ERROR;
  }

  /* loop molecules */
  for(m=0; m<n_molecules; m++) {
    mol_size = topology[m].part.n;
    if(mol_size > 0) {
      /* calc center of mass */
      calc_mol_center_of_mass(topology[m],com);
      /* fold coordinates */
      for(i=0; i<3; i++) {
	if ( PERIODIC(i) ) { 
	  tmp = (int)floor((com[i]+shift[i])*box_l_i[i]);
	  cm_tmp =0.0;
	  for(p=0; p<mol_size; p++) {
	    ind        = 3*topology[m].part.e[p] + i;
	    coord[ind] -= tmp*box_l[i];
	    coord[ind] += shift[i];
	    cm_tmp     += coord[ind];
	  }
	  cm_tmp /= (double)mol_size;
      if(cm_tmp < -10e-6 || cm_tmp > box_l[i]+10e-6) {
          ostringstream msg;
          msg <<"analyze_fold_molecules: chain center of mass is out of range (coord " << i << ": " << cm_tmp << " not in box_l " << box_l[i] << ")";
          runtimeError(msg);
	    return ES_ERROR;
	  }
	}
      }
    }
  }
  return ES_OK;
}

/*!
  \internal
*/
void QSInterpreter::parseError()
{
#ifdef QSDEBUGGER
  int sourceId = d->interpreter->debuggerEngine()->sourceId();
  QString msg = d->interpreter->errorMessages().first();
  QString script = d->interpreter->nameOfSourceId(sourceId);
  int line = d->interpreter->errorLines().first();
  runtimeError(msg, script, line);
#else
  runtimeError(d->interpreter->errorMessages().first(), QString::null,
               d->interpreter->errorLines().first());
#endif
}

static void constructor(AbstractActorInstance *pBase) 
{
  ActorInstance_art_Source *thisActor = (ActorInstance_art_Source*)pBase;

  if (thisActor->filename == NULL) {
    runtimeError(pBase,"Parameter not set: fileName");
  } else {
    thisActor->file = fopen(thisActor->filename, "r");
    if (thisActor->file == NULL) {
      runtimeError(pBase,"Cannot open file for output: %s: %s", 
                   thisActor->filename, strerror(errno));
    }
  }
}

/** 
    Calculate the center of mass, total mass, velocity, total force, and trap force on all trapped molecules 
*/
void calc_mol_info () {

  /* list of trapped molecules on this node */
  IntList local_trapped_mols;

  /* check to see if all the topology information has been synced to the various slave nodes */
  if ( !topo_part_info_synced ) {
      ostringstream msg;
      msg << "can't calculate molforces: must execute analyse set topo_part_sync first";
      runtimeError(msg);
    return;
  }

  init_intlist(&local_trapped_mols);

  /* Find out which trapped molecules are on this node */
  get_local_trapped_mols(&local_trapped_mols);

  /* Calculate the center of mass, mass, velocity, force of whatever fraction of each trapped molecule is on this node*/
  calc_local_mol_info(&local_trapped_mols);

  /* Communicate all this molecular information between nodes.
     It is all sent to the master node which combines it, calculates the trap forces,
     and sends the information back */
  if (this_node == 0) { 
    mpi_comm_mol_info(&local_trapped_mols);
  } else {
    mpi_comm_mol_info_slave(&local_trapped_mols);
  }

  realloc_intlist(&local_trapped_mols,0);
}

/* A list of trapped molecules present on this node is created (local_trapped_mols)*/
void get_local_trapped_mols (IntList *local_trapped_mols)
{
  int c, i, mol, j, fixed;

  for (c = 0; c < local_cells.n; c++) {
    for(i = 0; i < local_cells.cell[c]->n; i++) {
      mol = local_cells.cell[c]->part[i].p.mol_id;
      if ( mol >= n_molecules ) {
      ostringstream msg;
      msg <<"can't calculate molforces no such molecule as " << mol ;
      runtimeError(msg);
	return;
      }

      /* Check to see if this molecule is fixed */
      fixed =0;
      for(j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
	if (topology[mol].trap_flag & COORD_FIXED(j)) fixed = 1;
	if (topology[mol].noforce_flag & COORD_FIXED(j)) fixed = 1;
#endif
      }  
      if (fixed) {
	/* if this molecule isn't already in local_trapped_mols then add it in */
	if (!intlist_contains(local_trapped_mols,mol)) {
	  realloc_intlist(local_trapped_mols, local_trapped_mols->max + 1);
	  local_trapped_mols->e[local_trapped_mols->max-1] = mol;
	  local_trapped_mols->n = local_trapped_mols->max;
	}
      }
    }
  }
}

//---------------------------------------------------
bool skInterpreter::insertArrayValue(skStackFrame& frame,skRValue& robject, skExprNode * array_index,const skString& attr,const skRValue& value)
//---------------------------------------------------
{
    bool found=false;
    if (robject.type()==skRValue::T_Object) {
        skRValue rExpr;
        SAVE_VARIABLE(rExpr);
        rExpr=evaluate(frame,array_index);
        found=robject.obj()->setValueAt(rExpr,attr,value);
        if (found==false)
            runtimeError(frame,skSTR("Setting array member failed\n"));
        RELEASE_VARIABLE(rExpr);
    } else
        runtimeError(frame,skSTR("Cannot set an array member on a non-object\n"));
    return found;
}

TypePtr execSubscripts(TypePtr typePtr)
{
	//----------------------------------------
	// Loop to execute bracketed subscripts...
	while (codeToken == TKN_LBRACKET)
	{
		do
		{
			getCodeToken();
			execExpression();
			int32_t subscriptValue = tos->integer;
			pop();
			//-------------------------
			// Range check the index...
			if ((subscriptValue < 0) || (subscriptValue >= typePtr->info.array.elementCount))
				runtimeError(ABL_ERR_RUNTIME_VALUE_OUT_OF_RANGE);
			tos->address += (subscriptValue * typePtr->info.array.elementTypePtr->size);
			if (codeToken == TKN_COMMA)
				typePtr = typePtr->info.array.elementTypePtr;
		} while (codeToken == TKN_COMMA);
		getCodeToken();
		if (codeToken == TKN_LBRACKET)
			typePtr = typePtr->info.array.elementTypePtr;
	}
	return (typePtr->info.array.elementTypePtr);
}

int change_particle_bond(int part, int *bond, int _delete)
{
  int pnode;
  if (!particle_node)
    build_particle_node();

  if (part < 0 || part > max_seen_particle)
    return ES_ERROR;
  pnode = particle_node[part];

  if (pnode == -1)
    return ES_ERROR;
  if(_delete != 0 || bond == NULL)
    _delete = 1;

  if (bond != NULL) {
    if (bond[0] < 0 || bond[0] >= n_bonded_ia) {
        ostringstream msg;
        msg <<"invalid/unknown bonded interaction type " << bond[0];
        runtimeError(msg);
      return ES_ERROR;
    }
  }
  return mpi_send_bond(pnode, part, bond, _delete);
}

void execRepeatStatement(void)
{
	PSTR loopStartLocation = codeSegmentPtr;
	int32_t iterations = 0;
	do
	{
		getCodeToken();
		if(codeToken != TKN_UNTIL)
			do
			{
				execStatement();
				if(ExitWithReturn)
					return;
			}
			while(codeToken != TKN_UNTIL);
		//---------------------------
		// Check for infinite loop...
		iterations++;
		if(iterations == MaxLoopIterations)
			runtimeError(ABL_ERR_RUNTIME_INFINITE_LOOP);
		//-------------------------------
		// Eval the boolean expression...
		getCodeToken();
		execExpression();
		if(tos->integer == 0)
			codeSegmentPtr = loopStartLocation;
		//--------------------------
		// Grab the boolean value...
		pop();
	}
	while(codeSegmentPtr == loopStartLocation);
}

/** Calculate the local fluid momentum.
 * The calculation is implemented explicitly for the special case of D3Q19.
 * @param index The local lattice site (Input).
 * @param j local fluid speed
 */
inline void lb_calc_local_j(index_t index, double *j) {

#ifndef D3Q19
#error Only D3Q19 is implemened!
#endif
  if (!(lattice_switch & LATTICE_LB)) {
      ostringstream msg;
      msg <<"Error in lb_calc_local_j in " << __FILE__ << __LINE__ << ": CPU LB not switched on.";
      runtimeError(msg);
    j[0]=j[1]=j[2]=0;
    return;
  }

  j[0] =   lbfluid[0][1][index]  - lbfluid[0][2][index]
         + lbfluid[0][7][index]  - lbfluid[0][8][index]  
         + lbfluid[0][9][index]  - lbfluid[0][10][index] 
         + lbfluid[0][11][index] - lbfluid[0][12][index] 
         + lbfluid[0][13][index] - lbfluid[0][14][index];
  j[1] =   lbfluid[0][3][index]  - lbfluid[0][4][index]
         + lbfluid[0][7][index]  - lbfluid[0][8][index]  
         - lbfluid[0][9][index]  + lbfluid[0][10][index]
         + lbfluid[0][15][index] - lbfluid[0][16][index] 
         + lbfluid[0][17][index] - lbfluid[0][18][index]; 
  j[2] =   lbfluid[0][5][index]  - lbfluid[0][6][index]  
         + lbfluid[0][11][index] - lbfluid[0][12][index] 
         - lbfluid[0][13][index] + lbfluid[0][14][index]
         + lbfluid[0][15][index] - lbfluid[0][16][index] 
         - lbfluid[0][17][index] + lbfluid[0][18][index];

}

int mindist3(int part_id, double r_catch, int *ids) {
  Particle *partCfgMD;
  double dx,dy,dz;
  int i, me, caught=0;

  partCfgMD = (Particle*)malloc(n_part*sizeof(Particle));
  mpi_get_particles(partCfgMD, NULL);
  me = -1; /* Since 'mpi_get_particles' returns the particles unsorted, it's most likely that 'partCfgMD[i].p.identity != i'
	      --> prevent that! */
  for(i=0; i<n_part; i++) if (partCfgMD[i].p.identity == part_id) me = i; 
  if (me == -1) {
      ostringstream msg;
      msg <<"failed to find desired particle " << part_id;
      runtimeError(msg);
    return 0;
  }
  for (i=0; i<n_part; i++) {
    if (i != me) {
      dx = partCfgMD[me].r.p[0] - partCfgMD[i].r.p[0];   dx -= dround(dx/box_l[0])*box_l[0];
      dy = partCfgMD[me].r.p[1] - partCfgMD[i].r.p[1];   dy -= dround(dy/box_l[1])*box_l[1];
      dz = partCfgMD[me].r.p[2] - partCfgMD[i].r.p[2];   dz -= dround(dz/box_l[2])*box_l[2];
      if (sqrt(SQR(dx)+SQR(dy)+SQR(dz)) < r_catch) ids[caught++]=partCfgMD[i].p.identity;
    }
  }
  free(partCfgMD); 
  return (caught);
}

void freeLocal(SymTableNodePtr idPtr)
{
	//---------------------------------------
	// Frees data allocated on local stack...
	TypePtr typePtr = (TypePtr)(idPtr->typePtr);
	StackItemPtr itemPtr = nullptr;
	if(((typePtr->form == FRM_ARRAY) /* || (typePtr->form == FRM_RECORD)*/) &&
			(idPtr->defn.key != DFN_REFPARAM))
	{
		switch(idPtr->defn.info.data.varType)
		{
			case VAR_TYPE_NORMAL:
				itemPtr = stackFrameBasePtr + idPtr->defn.info.data.offset;
				break;
//			case VAR_TYPE_ETERNAL:
//				itemPtr = stack + idPtr->defn.info.data.offset;
//				break;
//			case VAR_TYPE_STATIC:
//				itemPtr = StaticDataPtr + idPtr->defn.info.data.offset;
//				break;
		}
		if(!itemPtr)
			runtimeError(0);
		else
			ABLStackFreeCallback(itemPtr->address);
	}
}

/** Calculate the local fluid density.
 * The calculation is implemented explicitly for the special case of D3Q19.
 * @param index the local lattice site (Input).
 * @param rho   local fluid density
 */
inline void lb_calc_local_rho(index_t index, double *rho) {

#ifndef D3Q19
#error Only D3Q19 is implemened!
#endif

  // unit conversion: mass density
  if (!(lattice_switch & LATTICE_LB)) {
      ostringstream msg;
      msg <<"Error in lb_calc_local_rho in " << __FILE__ << __LINE__ << ": CPU LB not switched on.";
      runtimeError(msg);
    *rho =0;
    return;
  }

  double avg_rho = lbpar.rho[0]*lbpar.agrid*lbpar.agrid*lbpar.agrid;

  *rho =   avg_rho
         + lbfluid[0][0][index]
         + lbfluid[0][1][index]  + lbfluid[0][2][index]
         + lbfluid[0][3][index]  + lbfluid[0][4][index]
         + lbfluid[0][5][index]  + lbfluid[0][6][index] 
         + lbfluid[0][7][index]  + lbfluid[0][8][index]  
	       + lbfluid[0][9][index]  + lbfluid[0][10][index]
         + lbfluid[0][11][index] + lbfluid[0][12][index] 
	       + lbfluid[0][13][index] + lbfluid[0][14][index] 
         + lbfluid[0][15][index] + lbfluid[0][16][index] 
	       + lbfluid[0][17][index] + lbfluid[0][18][index];

}

void pushAddress(Address address)
{
	StackItemPtr valuePtr = ++tos;
	if(valuePtr >= &stack[MAXSIZE_STACK])
		runtimeError(ABL_ERR_RUNTIME_STACK_OVERFLOW);
	valuePtr->address = address;
}

void pushBoolean(bool value)
{
	StackItemPtr valuePtr = ++tos;
	if(valuePtr >= &stack[MAXSIZE_STACK])
		runtimeError(ABL_ERR_RUNTIME_STACK_OVERFLOW);
	valuePtr->integer = (value ? 1 : 0);
}

void pushReal(float value)
{
	StackItemPtr valuePtr = ++tos;
	if(valuePtr >= &stack[MAXSIZE_STACK])
		runtimeError(ABL_ERR_RUNTIME_STACK_OVERFLOW);
	valuePtr->real = value;
}