自学教程：C++ ERROR_SPRINTF函数代码示例

Particle *get_mol_com_particle(Particle *calling_p) {    int mol_id;    int i;    Particle *p;    mol_id=calling_p->p.mol_id;    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) {            char *errtxt = runtime_error(128 + 3*TCL_INTEGER_SPACE);            ERROR_SPRINTF(errtxt,"Particle does not exist in put_mol_force_on_parts! id=%i/n",topology[mol_id].part.e[i]);            return NULL;        }#endif        if (ifParticleIsVirtual(p)) {            return p;        }    }#ifdef VIRTUAL_SITES_DEBUG    char *errtxt = runtime_error(128 + 3*TCL_INTEGER_SPACE);    ERROR_SPRINTF(errtxt,"No com found in get_mol_com_particleParticle does not exist in put_mol_force_on_parts! pnr=%i/n",calling_p->p.identity);    return NULL;#endif    return calling_p;}

int updatePartCfg(int bonds_flag){  int j;  if(partCfg)    return 1;  partCfg = (Particle*)malloc(n_part*sizeof(Particle));  if (bonds_flag != WITH_BONDS)    mpi_get_particles(partCfg, NULL);  else    mpi_get_particles(partCfg,&partCfg_bl);  for(j=0; j<n_part; j++)    unfold_position(partCfg[j].r.p,partCfg[j].l.i);  partCfgSorted = 0;#ifdef VIRTUAL_SITES  if (!sortPartCfg()) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{094 could not sort partCfg} ");    return 0;  }  if (!updatePartCfg(bonds_flag)) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{094 could not update positions of virtual sites in partcfg } ");    return 0;  }#endifreturn 1;}

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) {     char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);     ERROR_SPRINTF(errtxt,"Particle does not have a mol id! pnr=%i/n",		   calling_p->p.identity);     return NULL;   }   for (i=0;i<topology[mol_id].part.n;i++){      p=local_particles[topology[mol_id].part.e[i]];      if (p==NULL){         char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);         ERROR_SPRINTF(errtxt,"Particle does not exist in put_mol_force_on_parts! id=%i/n",topology[mol_id].part.e[i]);         return NULL;      }      if (ifParticleIsVirtual(p)) {          return p;       }   }   char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);   ERROR_SPRINTF(errtxt,"No com found in get_mol_com_particleParticle does not exist in put_mol_force_on_parts! pnr=%i/n",calling_p->p.identity);   return NULL;   return calling_p;}

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) {    char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);    ERROR_SPRINTF(errtxt, "{analyze stress_profile: getintersection: intersection is not between the two given particles - %e is not between %e and %e and box size is %e, given is %d/n ",value,0.0,p2r[given],box_size[given],given);    return 0;   } else if (given == get) {    *answer =  drem_down(value + pos1[given],box_size[given]);;  } else if (0==p2r[given]) {    char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);    ERROR_SPRINTF(errtxt, "{analyze stress_profile: getintersection: intersection is a line, not a point - value is %g same as %g and %g/n",value,0.0,p2r[given]);    return 0;     } else {    *answer =  drem_down(pos1[get]+p2r[get]/p2r[given]*value,box_size[get]);  }  return 1;}

/** This function takes a given grid supplied by the user and    determines the correct orientation in which to do the fourier    transform.  In this regard one of the grid dimensions must be 0    and the other two must be integer multiples of two and equal to    each other.  The dimension that is 0 will be assigned an internal    reference called zdir and will be used to calculate the height    function used in the fft*/void map_to_2dgrid() {  int i;  STAT_TRACE(fprintf(stderr,"%d,executing map_to_2dgrid /n",this_node));  /* Reset values of mapping */  xdir = -1;  ydir = -1;  zdir = -1;  /* Find the grid normal */  for ( i = 0 ; i < 3 ; i++) {    if ( mode_grid_3d[i] == 0 ) {      if (zdir != -1 ) { /* grid normal must be unique */ 	char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);	ERROR_SPRINTF(errtxt, "{029 fft_modes_init: grid dimensions are <%d,%d,%d>, but one and only one must be = 0} ",		mode_grid_3d[0],mode_grid_3d[1],mode_grid_3d[2]);	return;      } else {	zdir = i;      }    }     else if ( mode_grid_3d[i] < 0 ) {      char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);      ERROR_SPRINTF(errtxt, "{030 fft_modes_init: grid dimensions are <%d,%d,%d>, but all must be >= 0} ",	      mode_grid_3d[0],mode_grid_3d[1],mode_grid_3d[2]);      return;    }    else {      if (  xdir == -1 ) {xdir = i;}      else {ydir = i;}          }      }  /* Check that grid normal was found */  if ( zdir == -1 ) {    char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);    ERROR_SPRINTF(errtxt, "{031 fft_modes_init: grid dimensions are <%d,%d,%d>, but one and only one must be = 0} ",	    mode_grid_3d[0],mode_grid_3d[1],mode_grid_3d[2]);    return;  }  STAT_TRACE(fprintf(stderr,		     "%d,map_to_2dgrid found the following mapping: xdir = %d, ydir = %d, zdir = %d /n",		     this_node, xdir, ydir, zdir));  /* Now that we know the grid normal check that the other two dimensions are equal and multiples of 2 */  if ( mode_grid_3d[xdir] != mode_grid_3d[ydir] ) {    char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);    ERROR_SPRINTF(errtxt, "{032 fft_modes_init: grid dimensions are <%d,%d,%d>, but two must be equal and the other 0} ",	    mode_grid_3d[xdir],mode_grid_3d[ydir],mode_grid_3d[zdir]);    return;  }  if ( (mode_grid_3d[xdir]/2.0 - floor(mode_grid_3d[xdir]/2.0) > MODES2D_NUM_TOL)        || (mode_grid_3d[ydir]/2.0 - floor(mode_grid_3d[ydir]/2.0) > MODES2D_NUM_TOL) ) {    char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);    ERROR_SPRINTF(errtxt, "{033 fft_modes_init: grid dimensions are <%d,%d,%d>. All non zero values must be integer multiples of 2} ",	    mode_grid_3d[xdir],mode_grid_3d[ydir],mode_grid_3d[zdir]);    return;  }}

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){         char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);         ERROR_SPRINTF(errtxt,"Particle does not exist in put_mol_force_on_parts! id=%i/n",topology[mol_id].part.e[i]);         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){         char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);         ERROR_SPRINTF(errtxt,"Particle does not exist in put_mol_force_on_parts! id=%i/n",topology[mol_id].part.e[i]);         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){      char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);      ERROR_SPRINTF(errtxt,"There is more than one COM input_mol_force_on_parts! mol_id=%i/n",mol_id);      return;   }#endif}

/* Get a complete list of the orientations of every lipid assuming a   bilayer structure.  Requires grid*/int get_lipid_orients(IntList* l_orient) {  int i,gi,gj, atom;  double zreflocal,zref;    double dir[3];  double refdir[3] = {0,0,1};  double grid_size[2];  double* height_grid;  if ( xdir + ydir + zdir == -3 || mode_grid_3d[xdir] <= 0 || mode_grid_3d[ydir] <= 0 ) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{036 cannot calculate lipid orientations with uninitialized grid} ");    return ES_ERROR;  }  /* Allocate memory for height grid arrays and initialize these arrays */  height_grid = (double*) malloc((mode_grid_3d[xdir])*sizeof(double)*mode_grid_3d[ydir]);  /* Calculate physical size of grid mesh */  grid_size[xdir] = box_l[xdir]/(double)mode_grid_3d[xdir];  grid_size[ydir] = box_l[ydir]/(double)mode_grid_3d[ydir];  /* Update particles */  updatePartCfg(WITHOUT_BONDS);  //Make sure particles are sorted  if (!sortPartCfg()) {    fprintf(stderr,"%d,could not sort partCfg /n",this_node);    return -1;  }    if ( !calc_fluctuations(height_grid, 1) ) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{034 calculation of height grid failed } ");    return -1;  }  zref = calc_zref( zdir );  for ( i = 0 ; i < n_molecules ; i++) {    atom = topology[i].part.e[0];    gi = floor( partCfg[atom].r.p[xdir]/grid_size[xdir] );    gj = floor( partCfg[atom].r.p[ydir]/grid_size[ydir] );    zreflocal = height_grid[gj+gi*mode_grid_3d[xdir]] + zref;    l_orient->e[i] = lipid_orientation(atom,partCfg,zreflocal,dir,refdir);  }  free(height_grid);  return 1;}

void Lattice::interpolate_linear(double* pos, double* value) {    int left_halo_index[3];    double d[3];    if (this->halo_size <= 0) {        char* c = runtime_error(128);        ERROR_SPRINTF(c, "Error in interpolate_linear: halo size is 0");        return;    }    for (int dim = 0; dim<3; dim++) {        left_halo_index[dim]=(int) floor((pos[dim]-this->local_offset[dim])/this->agrid[dim]) + this->halo_size;        d[dim]=((pos[dim]-this->local_offset[dim])/this->agrid[dim] - floor((pos[dim]-this->local_offset[dim])/this->agrid[dim]));        if (left_halo_index[dim] < 0 || left_halo_index[dim] >= this->halo_grid[dim]) {            char* c = runtime_error(128);            ERROR_SPRINTF(c, "Error in interpolate_linear: Particle out of range");            return;        }    }    double w[8];    index_t index[8];    w[0] = (1-d[0])*(1-d[1])*(1-d[2]);    index[0]=get_linear_index(   left_halo_index[0], left_halo_index[1], left_halo_index[2], this->halo_grid);    w[1] = ( +d[0])*(1-d[1])*(1-d[2]);    index[1]=get_linear_index(   left_halo_index[0]+1, left_halo_index[1], left_halo_index[2], this->halo_grid);    w[2] = (1-d[0])*( +d[1])*(1-d[2]);    index[2]=get_linear_index(   left_halo_index[0], left_halo_index[1]+1, left_halo_index[2], this->halo_grid);    w[3] = ( +d[0])*( +d[1])*(1-d[2]);    index[3]=get_linear_index(   left_halo_index[0]+1, left_halo_index[1]+1, left_halo_index[2], this->halo_grid);    w[4] = (1-d[0])*(1-d[1])*( +d[2]);    index[4]=get_linear_index(   left_halo_index[0], left_halo_index[1], left_halo_index[2]+1, this->halo_grid);    w[5] = ( +d[0])*(1-d[1])*( +d[2]);    index[5]=get_linear_index(   left_halo_index[0]+1, left_halo_index[1], left_halo_index[2]+1, this->halo_grid);    w[6] = (1-d[0])*( +d[1])*( +d[2]);    index[6]=get_linear_index(   left_halo_index[0], left_halo_index[1]+1, left_halo_index[2]+1, this->halo_grid);    w[7] = ( +d[0])*( +d[1])*( +d[2]);    index[7]=get_linear_index(   left_halo_index[0]+1, left_halo_index[1]+1, left_halo_index[2]+1, this->halo_grid);    for (unsigned int i = 0; i<this->dim; i++) {        value[i] = 0;    }    double* local_value;    for (unsigned int i=0; i<8; i++) {        get_data_for_linear_index(index[i], (void**) &local_value);        for (unsigned int j = 0; j<this->dim; j++) {            value[j]+=w[i]*local_value[j];        }    }}

int observable_calc_blocked_com_force(observable* self) {  double* A = self->last_value;  unsigned int i;  unsigned int block;  unsigned int n_blocks;  unsigned int blocksize;  unsigned int id;  IntList* ids;  if (!sortPartCfg()) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{094 could not sort partCfg} ");    return -1;  }  ids=(IntList*) self->container;  n_blocks=self->n/3;   blocksize=ids->n/n_blocks;  for ( block = 0; block < n_blocks; block++ ) {    for ( i = 0; i < blocksize; i++ ) {      id = ids->e[block*blocksize+i];      if (ids->e[i] >= n_part)        return 1;      A[3*block+0] +=  partCfg[id].f.f[0]/time_step/time_step*2;      A[3*block+1] +=  partCfg[id].f.f[1]/time_step/time_step*2;      A[3*block+2] +=  partCfg[id].f.f[2]/time_step/time_step*2;    }  }  return 0;}

int observable_calc_blocked_com_position(observable* self) {  double* A = self->last_value;  unsigned int i;  unsigned int block;  unsigned int n_blocks;  unsigned int blocksize;  unsigned int id;  double total_mass = 0;  IntList* ids;  if (!sortPartCfg()) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{094 could not sort partCfg} ");    return -1;  }  ids=(IntList*) self->container;  n_blocks=self->n/3;   blocksize=ids->n/n_blocks;  for ( block = 0; block < n_blocks; block++ ) {    total_mass = 0;    for ( i = 0; i < blocksize; i++ ) {      id = ids->e[block*blocksize+i];      if (ids->e[i] >= n_part)        return 1;      A[3*block+0] +=  PMASS(partCfg[id])*partCfg[id].r.p[0];      A[3*block+1] +=  PMASS(partCfg[id])*partCfg[id].r.p[1];      A[3*block+2] +=  PMASS(partCfg[id])*partCfg[id].r.p[2];      total_mass += PMASS(partCfg[ids->e[i]]);    }    A[3*block+0] /=  total_mass;    A[3*block+1] /=  total_mass;    A[3*block+2] /=  total_mass;  }  return 0;}

void cells_on_max_cut_change(int shrink){  double old_max_range = max_range;  calc_maximal_cutoff();  if (max_cut > 0.0) {    if (skin >= 0.0)      max_range = max_cut + skin;    else      /* if the skin is not yet set, assume zero. */      max_range = max_cut;  }  else    /* if no interactions yet, we also don't need a skin */    max_range = 0.0;  /* no need to do something if     1. the range didn't change numerically (<= necessary for the start case,     when max_range and old_max_range == 0.0)     2. it shrank, and we shouldn't shrink (NpT) */  if ((fabs(max_range - old_max_range) <= ROUND_ERROR_PREC * max_range) ||      (!shrink && (max_range < old_max_range)))    return;    cells_re_init(CELL_STRUCTURE_CURRENT);  for (int i = 0; i < 3; i++)    if (local_box_l[i] < max_range) {      char *errtext = runtime_error(128 + TCL_INTEGER_SPACE);      ERROR_SPRINTF(errtext,"{013 box_l in direction %d is still too small} ", i);    }}

int MMM1D_sanity_checks(){  char *errtxt;  if (PERIODIC(0) || PERIODIC(1) || !PERIODIC(2)) {    errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt, "{022 MMM1D requires periodicity 0 0 1} ");    return 1;  }  if (cell_structure.type != CELL_STRUCTURE_NSQUARE) {    errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt, "{023 MMM1D requires n-square cellsystem} ");    return 1;  }  return 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 ) {	char *errtxt = runtime_error(128 + 3*TCL_INTEGER_SPACE);	ERROR_SPRINTF(errtxt, "{ 094 can't calculate molforces no such molecule as %d }",mol);	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;	}      }    }  }}

/** 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)) {    ERROR_SPRINTF(runtime_error(128),         "{ Error in lb_calc_local_rho in %s %d: CPU LB not switched on. } ", __FILE__, __LINE__);    *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];}

/** 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)) {    ERROR_SPRINTF(runtime_error(128),         "{ Error in lb_calc_local_j in %s %d: CPU LB not switched on. } ", __FILE__, __LINE__);    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];}

void correct_pos_shake(){   int    repeat_,  cnt=0;   int repeat=1;   while (repeat!=0 && cnt<SHAKE_MAX_ITERATIONS)   {     init_correction_vector();     repeat_ = 0;     compute_pos_corr_vec(&repeat_);     ghost_communicator(&cell_structure.collect_ghost_force_comm);     app_pos_correction();     /**Ghost Positions Update*/     ghost_communicator(&cell_structure.update_ghost_pos_comm);     if(this_node==0)         MPI_Reduce(&repeat_, &repeat, 1, MPI_INT, MPI_SUM, 0, comm_cart);     else         MPI_Reduce(&repeat_, NULL, 1, MPI_INT, MPI_SUM, 0, comm_cart);     MPI_Bcast(&repeat, 1, MPI_INT, 0, comm_cart);       cnt++;   }// while(repeat) loop   if (cnt >= SHAKE_MAX_ITERATIONS) {     char *errtxt = runtime_error(100 + ES_INTEGER_SPACE);     ERROR_SPRINTF(errtxt, "{053 RATTLE failed to converge after %d iterations} ", cnt);   }   check_resort_particles();}

// Update the pos of the given virtual particle as defined by the real // particles in the same moleculevoid update_mol_pos_particle(Particle *p){ // First obtain the real particle responsible for this virtual particle: // Find the 1st real particle in the topology for the virtual particle's mol_id Particle *p_real = vs_relative_get_real_particle(p); // Check, if a real particle was found if (!p_real) {   char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);   ERROR_SPRINTF(errtxt,"virtual_sites_relative.cpp - update_mol_pos_particle(): No real particle associated with virtual site./n");   return; }  // Calculate the quaternion defining the orientation of the vecotr connectinhg // the virtual site and the real particle // This is obtained, by multiplying the quaternion representing the director // of the real particle with the quaternion of the virtual particle, which  // specifies the relative orientation. double q[4]; multiply_quaternions(p_real->r.quat,p->r.quat,q); // Calculate the director resulting from the quaternions double director[3]; convert_quat_to_quatu(q,director); // normalize double l =sqrt(sqrlen(director)); // Division comes in the loop below  // Calculate the new position of the virtual sites from // position of real particle + director  int i; double new_pos[3]; double tmp; for (i=0;i<3;i++) {  new_pos[i] =p_real->r.p[i] +director[i]/l*p->p.vs_relative_distance;  // Handle the case that one of the particles had gone over the periodic  // boundary and its coordinate has been folded  #ifdef PARTIAL_PERIODIC   if (PERIODIC(i))   #endif  {    tmp =p->r.p[i] -new_pos[i];    //printf("%f/n",tmp);    if (tmp > box_l[i]/2.) {     //printf("greater than box_l/2 %f/n",tmp);     p->r.p[i] =new_pos[i] + box_l[i];    }    else if (tmp < -box_l[i]/2.) {     //printf("smaller than box_l/2 %f/n",tmp);     p->r.p[i] =new_pos[i] - box_l[i];    }    else p->r.p[i] =new_pos[i];   }   #ifdef PARTIAL_PERIODIC    else p->r.p[i] =new_pos[i];   #endif//   fold_coordinate(p->r.p,p->l.i,i); }}

/**     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 ) {    char *errtxt = runtime_error(128 + 3*TCL_INTEGER_SPACE);    ERROR_SPRINTF(errtxt, "{ 093 can't calculate molforces: must execute analyse set topo_part_sync first }");    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);}

int observable_calc_flux_density_profile(observable* self) {  double* A = self->last_value;  int binx, biny, binz;  double ppos[3];  double x, y, z;  int img[3];  double bin_volume;  IntList* ids;#ifdef LEES_EDWARDS  double v_le[3];#endif  if (!sortPartCfg()) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{094 could not sort partCfg} ");    return -1;  }  profile_data* pdata;  pdata=(profile_data*) self->container;  ids=pdata->id_list;  double xbinsize=(pdata->maxx - pdata->minx)/pdata->xbins;  double ybinsize=(pdata->maxy - pdata->miny)/pdata->ybins;  double zbinsize=(pdata->maxz - pdata->minz)/pdata->zbins;  double v[3];  double v_x, v_y, v_z;      for (int i = 0; i< self->n; i++ ) {    A[i]=0;  }  for (int i = 0; i<ids->n; i++ ) {    if (ids->e[i] >= n_part)      return 1;/* We use folded coordinates here */    v[0]=partCfg[ids->e[i]].m.v[0]*time_step;    v[1]=partCfg[ids->e[i]].m.v[1]*time_step;    v[2]=partCfg[ids->e[i]].m.v[2]*time_step;    memcpy(ppos, partCfg[ids->e[i]].r.p, 3*sizeof(double));    memcpy(img, partCfg[ids->e[i]].l.i, 3*sizeof(int));    fold_position(ppos, img);    x=ppos[0];    y=ppos[1];    z=ppos[2];    binx  =(int)floor((x-pdata->minx)/xbinsize);    biny  =(int)floor((y-pdata->miny)/ybinsize);    binz  =(int)floor((z-pdata->minz)/zbinsize);    if (binx>=0 && binx < pdata->xbins && biny>=0 && biny < pdata->ybins && binz>=0 && binz < pdata->zbins) {      bin_volume=xbinsize*ybinsize*zbinsize;      v_x=v[0];      v_y=v[1];      v_z=v[2];      A[3*(binx*pdata->ybins*pdata->zbins + biny*pdata->zbins + binz) + 0] += v_x/bin_volume;      A[3*(binx*pdata->ybins*pdata->zbins + biny*pdata->zbins + binz) + 1] += v_y/bin_volume;      A[3*(binx*pdata->ybins*pdata->zbins + biny*pdata->zbins + binz) + 2] += v_z/bin_volume;    }   }  return 0;}

void Lattice::interpolate(double* pos, double* value) {    if (this->interpolation_type == INTERPOLATION_LINEAR) {        interpolate_linear(pos, value);    } else {        char* c = runtime_error(128);        ERROR_SPRINTF(c, "Unknown interpolation type");    }}

/* 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];  /* check molecule information */  if ( n_molecules < 0 ) return (TCL_ERROR);  if (!sortPartCfg()) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt, "{059 analyze_fold_molecules: could not sort particle config, particle ids not consecutive?} ");    return TCL_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) {	    char *errtxt = runtime_error(128 + TCL_INTEGER_SPACE + 2*TCL_DOUBLE_SPACE);	    ERROR_SPRINTF(errtxt,"{060 analyze_fold_molecules: chain center of mass is out of range (coord %d: %.14f not in box_l %.14f)} ",		    i,cm_tmp,box_l[i]);	    return (TCL_ERROR);	  }	}      }    }  }  return (TCL_OK);}

int observable_stress_tensor(observable* self) {  if (!sortPartCfg()) {    char *errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt,"{094 could not sort partCfg} ");    return -1;  }  observable_compute_stress_tensor(1,self->last_value,self->n);  return 0;}

void integrator_sanity_checks(){  char *errtext;  if ( time_step < 0.0 ) {    errtext = runtime_error(128);    ERROR_SPRINTF(errtext, "{010 time_step not set} ");  }}

inline void check_particle_force(Particle *part) {  for (int i=0; i< 3; i++) {    if (isnan(part->f.f[i])) {      char *errtext = runtime_error(128);      ERROR_SPRINTF(errtext,"{999 force on particle %d was NAN.} ",		    part->p.identity);    }  }#ifdef ROTATION  for (int i=0; i< 3; i++) {    if (isnan(part->f.torque[i])) {      char *errtext = runtime_error(128);      ERROR_SPRINTF(errtext,"{999 force on particle %d was NAN.} ",		    part->p.identity);    }  }#endif}

/* but p_com is a real particle */void calc_mol_pos(Particle *p_com,double r_com[3]){   int i,j,mol_id;   double M=0;   double vec12[3];   Particle *p;#ifdef VIRTUAL_SITES_DEBUG   int count=0;#endif   for (i=0;i<3;i++){      r_com[i]=0.0;   }   mol_id=p_com->p.mol_id;   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){         char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);         ERROR_SPRINTF(errtxt,"Particle does not exist in calc_mol_pos! id=%i/n",topology[mol_id].part.e[i]);         return;      }      #endif      if (ifParticleIsVirtual(p)) continue;      get_mi_vector(vec12,p->r.p, p_com->r.p);      for (j=0;j<3;j++){          r_com[j] += PMASS(*p)*vec12[j];      }      M+=PMASS(*p);#ifdef VIRTUAL_SITES_DEBUG      count++;#endif   }   for (j=0;j<3;j++){      r_com[j] /= M;      r_com[j] += p_com->r.p[j];   }#ifdef VIRTUAL_SITES_DEBUG   if (count!=topology[mol_id].part.n-1){      char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);      ERROR_SPRINTF(errtxt,"There is more than one COM in calc_mol_pos! mol_id=%i/n",mol_id);      return;   }#endif}

inline void lb_local_fields_get_boundary_flag(index_t index, int *boundary) {    if (!(lattice_switch & LATTICE_LB)) {    ERROR_SPRINTF(runtime_error(128),         "{ Error in lb_local_fields_get_boundary_flag in %s %d: CPU LB not switched on. } ", __FILE__, __LINE__);    *boundary = 0;    return;  }  *boundary = lbfields[index].boundary;}

/**Compute positional corrections*/void compute_pos_corr_vec(int *repeat_){  Bonded_ia_parameters *ia_params;  int i, j, k, c, np, cnt=-1;  Cell *cell;  Particle *p, *p1, *p2;  double r_ij_t[3], r_ij[3], r_ij_dot, G, pos_corr, r_ij2;  for (c = 0; c < local_cells.n; c++) {    cell = local_cells.cell[c];    p  = cell->part;    np = cell->n;    for(i = 0; i < np; i++) {      p1 = &p[i];      k=0;      while(k<p1->bl.n) {	ia_params = &bonded_ia_params[p1->bl.e[k++]];	if( ia_params->type == BONDED_IA_RIGID_BOND ) {	  cnt++;	  p2 = local_particles[p1->bl.e[k++]];	  if (!p2) {	    char *errtxt = runtime_error(128 + 2*ES_INTEGER_SPACE);	    ERROR_SPRINTF(errtxt,"{051 rigid bond broken between particles %d and %d (particles not stored on the same node)} ",		    p1->p.identity, p1->bl.e[k-1]);	    return;	  }	  get_mi_vector(r_ij  , p1->r.p    , p2->r.p    );	  r_ij2 = sqrlen(r_ij);	  if(fabs(1.0 - r_ij2/ia_params->p.rigid_bond.d2) > ia_params->p.rigid_bond.p_tol) {            get_mi_vector(r_ij_t, p1->r.p_old, p2->r.p_old);	    r_ij_dot = scalar(r_ij_t, r_ij);	    G = 0.50*(ia_params->p.rigid_bond.d2 - r_ij2 )/r_ij_dot;#ifdef MASS	    G /= (PMASS(*p1)+PMASS(*p2));#else	    G /= 2;#endif	    for (j=0;j<3;j++) {	      pos_corr = G*r_ij_t[j];	      p1->f.f[j] += pos_corr*PMASS(*p2);	      p2->f.f[j] -= pos_corr*PMASS(*p1);	    }	    /*Increase the 'repeat' flag by one */	      *repeat_ = *repeat_ + 1;	  }	}	else	  /* skip bond partners of nonrigid bond */          k+=ia_params->num;      } //while loop    } //for i loop  } //for c loop}

/** Velocity correction vectors are computed*/void compute_vel_corr_vec(int *repeat_){  Bonded_ia_parameters *ia_params;  int i, j, k, c, np;  Cell *cell;  Particle *p, *p1, *p2;  double v_ij[3], r_ij[3], K, vel_corr;  for (c = 0; c < local_cells.n; c++)    {      cell = local_cells.cell[c];      p  = cell->part;      np = cell->n;      for(i = 0; i < np; i++) {        p1 = &p[i];        k=0;	while(k<p1->bl.n)	  {	    ia_params = &bonded_ia_params[p1->bl.e[k++]];	    if( ia_params->type == BONDED_IA_RIGID_BOND )	      {		p2 = local_particles[p1->bl.e[k++]];		if (!p2) {		  char *errtxt = runtime_error(128 + 2*ES_INTEGER_SPACE);		  ERROR_SPRINTF(errtxt,"{054 rigid bond broken between particles %d and %d (particles not stored on the same node)} ",			  p1->p.identity, p1->bl.e[k-1]);		  return;		}		vecsub(p1->m.v, p2->m.v, v_ij);                get_mi_vector(r_ij, p1->r.p, p2->r.p);		if(fabs(scalar(v_ij, r_ij)) > ia_params->p.rigid_bond.v_tol)	        {		  K = scalar(v_ij, r_ij)/ia_params->p.rigid_bond.d2;#ifdef MASS		  K /= (PMASS(*p1) + PMASS(*p2));#else		  K /= 2.0;#endif		for (j=0;j<3;j++)		  {		    vel_corr = K*r_ij[j];		    p1->f.f[j] -= vel_corr*PMASS(*p2);		    p2->f.f[j] += vel_corr*PMASS(*p1);		  }		  *repeat_ = *repeat_ + 1 ;		}	      }	    else	      k += ia_params->num;	  } //while loop      } //for i loop    } //for c loop}

double get_mol_dist(Particle *p1,Particle *p2){   Particle *p1_com,*p2_com;   double dist[3],dist2;   p1_com=get_mol_com_particle(p1);   p2_com=get_mol_com_particle(p2);   #ifdef VIRTUAL_SITES_DEBUG   if (p1_com==NULL){      char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);      ERROR_SPRINTF(errtxt,"COM Particle not found for particle in get_mol_dist id=%i/n",p1->p.identity);      dist[0]=dist[1]=dist[2]=0.0;   }   if (p2_com==NULL){      char *errtxt = runtime_error(128 + 3*ES_INTEGER_SPACE);      ERROR_SPRINTF(errtxt,"COM Particle not found for particle in get_mol_dist id=%i/n",p2->p.identity);      dist[0]=dist[1]=dist[2]=0.0;   }   #endif   get_mi_vector(dist,p1_com->r.p, p2_com->r.p);   dist2=SQR(dist[0])+SQR(dist[1])+SQR(dist[2]);   return sqrt(dist2);}

int mdlc_sanity_checks(){#ifdef PARTIAL_PERIODIC    char *errtxt;  if (!PERIODIC(0) || !PERIODIC(1) || !PERIODIC(2)) {    errtxt = runtime_error(128);    ERROR_SPRINTF(errtxt, "{006 mdlc requires periodicity 1 1 1} ");    return 1;  }#endif    return 0;}