int ComputeAtomsPerCell(int ***CellMap, 
			int NRows,int NCols, 
			int NParts){
//  Uncomment Lines for a parallel version
//  int MY_PE;
  int max = 0;
//  MPI_Comm_rank(MPI_COMM_WORLD,&MY_PE);
//  if(MY_PE == 0){
    double *rCellSizes = (double*)malloc(sizeof(double)*NRows*NCols);
    double ransum = 0.0;
    
    for(int i=0;i<NRows;i++){
      for(int j=0;j<NCols;j++){
	double ran = drand48();
	rCellSizes[i*NRows + j] = ran;
	ransum += ran;
      }
    }
    int molsum = 0;
    for(int i=0;i<NRows;i++){
      for(int j=0;j<NCols;j++){
	double ran = rCellSizes[i*NRows+j]/ransum;
	CellMap[i][j][3] = int(NParts*ran);
	molsum += CellMap[i][j][3];
	if(CellMap[i][j][3] > max) max = CellMap[i][j][3];
      }
    }
    free(rCellSizes);
    
    while(molsum < NParts){
      for(int i=0;i<NRows;i++){
	for(int j=0;j<NCols;j++){
	  CellMap[i][j][3]++;
	  molsum++;
	  if(CellMap[i][j][3] > max) max = CellMap[i][j][3];
	  if(molsum >= NParts) {
	    break;
	  }
	}
      }
    }
  }
 // MPI_Bcast(&max,1,MPI_INT,0,MPI_COMM_WORLD);
 // for(int i=0;i<NRows;i++){
 //   for(int j=0;j<NCols;j++){
 //     MPI_Bcast(&CellMap[i][j][3],1,MPI_INT,0,MPI_COMM_WORLD);
 //   }
 // }
  return max;
}