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;
  }
#endif
return 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];

}