lattice.cpp

用monte carlo方法计算化学反应的代码Oxygen Ion Ordering

      E_crystal += calc_inner_wall(lat,i,i,starty,endy,startz,endz,nbrs,calc_E_count);
    }

  //right wall sans corners
  if(sxy[0][1]==x)
    {
      vector<int> nbrs(6,0);
      i=X;
      get_start_and_end(starty,endy,sxy,Y,y,1);
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[right] = 1; nbrs[front]=1;//bottom, left,right,front
      E_crystal += calc_inner_wall(lat,i,i,starty,endy,startz,endz,nbrs,calc_E_count);
    }

  //do the top wall
  if(sxy[2][1]==z)
    {
      vector<int> nbrs(6,0);
      k=Z;
      get_start_and_end(startx,endx,sxy,X,x,0);
      get_start_and_end(starty,endy,sxy,Y,y,1);
      nbrs[top] = 1; nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;//top, bottom, left,front
      E_crystal += calc_inner_wall(lat,startx,endx,starty,endy,k,k,nbrs,calc_E_count);
    }
  
  //do the bottom wall
  if(sxy[2][0]==1)
    {
      vector<int> nbrs(6,0);
      k=1;
      get_start_and_end(startx,endx,sxy,X,x,0);
      get_start_and_end(starty,endy,sxy,Y,y,1);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;// bottom, left,front
      E_crystal += calc_inner_wall(lat,startx,endx,starty,endy,k,k,nbrs,calc_E_count);
    }
  
  //now we have the vertical edges left
  //calc the vertical 000 corner
  if((sxy[0][0]==1)&(sxy[1][0]==1))//x_start =  1 and y_start = 1
    {
      vector<int> nbrs(6,0);
      i = 1;
      j = 1;
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;// bottom, left,front
      E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
    }
  if((sxy[0][1]==x)&(sxy[1][0]==1))//x,0,0 corner        x_start =  X and y_start = 1
    {
      i = X;
      j = 1;
      vector<int> nbrs(6,0);
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[right]=1; nbrs[front]=1;// bottom, left,right,front
      E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
    }
  
  if((sxy[0][1]==x)&(sxy[1][1]==y))//x,y,0 corner   x_end =  X and y_end = Y
    {
      i=X;
      j=Y;
      vector<int> nbrs(6,0);
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[right]=1; nbrs[front]=1; nbrs[back] = 1;// bottom, left,right,front,back
      E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
    }
  
  
  if((sxy[0][0]==1)&(sxy[1][1]==y))//0,Y,0 corner        x_start =  1 and y_start = Y
    {
      i=1;
      j=Y;
      vector<int> nbrs(6,0);
      get_start_and_end(startz,endz,sxy,Z,z,2);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1; nbrs[back] = 1;// bottom, left,front,back
      E_crystal += calc_inner_wall(lat,i,i,j,j,startz,endz,nbrs,calc_E_count);
    }


    //now we go through and do the horizontal edges
    //go down x first
  if((sxy[2][0]==1)&(sxy[1][0]==1))//  z_start =  1 and y_start = 1
    {
      j=1;
      k=1;
      vector<int> nbrs(6,0);
      get_start_and_end(startx,endx,sxy,X,x,0);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;// bottom, left,front
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
    }
  if((sxy[2][1]==z)&(sxy[1][0]==1))//  z_end =  Z and y_start = 1
    {
      j=1;
      k=Z;
      vector<int> nbrs(6,0);
      get_start_and_end(startx,endx,sxy,X,x,0);
      nbrs[top]=1;nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;//top, bottom, left,front
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
    }
  if((sxy[2][0]==1)&(sxy[1][1]==y))//  z_start =  1 and y_end = Y
    {
      j=Y;
      k=1;
      vector<int> nbrs(6,0);
      get_start_and_end(startx,endx,sxy,X,x,0);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;nbrs[back]=1;// bottom, left,front,back
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
    }
  if((sxy[2][1]==z)&(sxy[1][1]==y))//  z_end = Z and y_end = Y
    {
      j=Y;
      k=Z;
      vector<int> nbrs(6,0);
      get_start_and_end(startx,endx,sxy,X,x,0);
      nbrs[top]=1;nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;nbrs[back]=1;//top, bottom, left,front,back
      E_crystal += calc_inner_wall(lat,startx,endx,j,j,k,k,nbrs,calc_E_count);
    }
  
  //now go down y
  if((sxy[2][0]==1)&(sxy[0][0]==1))//  z_start =  1 and x_start = 1
    {
      i=1;
      k=1;
      vector<int> nbrs(6,0);
      get_start_and_end(starty,endy,sxy,Y,y,1);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;//top, bottom, left,
      E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
    }
  if((sxy[2][1]==z)&(sxy[0][0]==1))//  z_end =  Z and x_start = 1
    {
      i=1;
      k=Z;
      vector<int> nbrs(6,0);
      get_start_and_end(starty,endy,sxy,Y,y,1);
      nbrs[top]=1; nbrs[bottom]=1; nbrs[left]=1; nbrs[front]=1;//top, bottom, left, front
      E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
    }
  if((sxy[2][0]==1)&(sxy[0][1]==x))//  z_start =  1 and x_end = X
    {
      i=X;
      k=1;
      vector<int> nbrs(6,0);
      get_start_and_end(starty,endy,sxy,Y,y,1);
      nbrs[bottom]=1; nbrs[left]=1; nbrs[right]=1;nbrs[front]=1;//bottom, left, right,front
      E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
    }
  if((sxy[2][1]==z)&(sxy[0][1]==x))//  z_end = Z and x_end = X
    {
      i=X;
      k=Z;
      vector<int> nbrs(6,0);
      get_start_and_end(starty,endy,sxy,Y,y,1);
      nbrs[top]=1;nbrs[bottom]=1; nbrs[left]=1; nbrs[right]=1;nbrs[front]=1;//top,bottom, left, right,front
      E_crystal += calc_inner_wall(lat,i,i,starty,endy,k,k,nbrs,calc_E_count);
    }
  

  //now the 8 corners are left
    if((sxy[0][0]==1)&(sxy[1][0]==1)&(sxy[2][0]==1)) //0,0,0 corner
      {
        i=1; j=1; k=1;
        E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
        E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
        E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
	calc_E_count+=3;
      }
  if((sxy[0][0]==1)&(sxy[1][1]==y)&(sxy[2][0]==1))// 0,Y,0 corner
    {
      i=1; j=Y; k=1;
      E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
      E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
      E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
      E_crystal += lat[i][j][k].calc_E(5,lat);//calc E for back node
      calc_E_count+=4;
    }
  if((sxy[0][1]==x)&(sxy[1][0]==1)&(sxy[2][0]==1))// X,0,0 corner
    {
      i=X; j=1; k=1;
      E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
      E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
      E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
      E_crystal += lat[i][j][k].calc_E(3,lat);//calc E for right node
      calc_E_count+=4;
    }
  if((sxy[0][0]==1)&(sxy[1][0]==1)&(sxy[2][1]==z))// 0,0,Z corner
    {
      i=1; j=1; k=Z;
      E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
      E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
      E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
      E_crystal += lat[i][j][k].calc_E(0,lat);//calc E for top node
      calc_E_count+=4;
    }
  if((sxy[0][1]==x)&(sxy[1][1]==y)&(sxy[2][1]==z))// X,Y,Z corner
    {
      i=X; j=Y; k=Z;
      E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
      E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
      E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
      E_crystal += lat[i][j][k].calc_E(0,lat);//calc E for top node
      E_crystal += lat[i][j][k].calc_E(5,lat);//calc E for back node
      E_crystal += lat[i][j][k].calc_E(3,lat);//calc E for right node
      calc_E_count+=6;
    }
  if((sxy[0][0]==1)&(sxy[1][1]==y)&(sxy[2][1]==z))// 0,Y,Z corner
    {
      i=1; j=Y; k=Z;
      E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
      E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
      E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
      E_crystal += lat[i][j][k].calc_E(0,lat);//calc E for top node
      E_crystal += lat[i][j][k].calc_E(5,lat);//calc E for back node
      calc_E_count+=5;
    }
  if((sxy[0][1]==x)&(sxy[1][0]==1)&(sxy[2][1]==z))// X,0,Z corner
    {
      i=X; j=1; k=Z;
      E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
      E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
      E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
      E_crystal += lat[i][j][k].calc_E(0,lat);//calc E for top node
      E_crystal += lat[i][j][k].calc_E(3,lat);//calc E for right node
      calc_E_count+=5;
    }
  if((sxy[0][1]==x)&(sxy[1][1]==y)&(sxy[2][0]==1))// X,Y,0 corner
    {
      i=X; j=Y; k=1;
      E_crystal += lat[i][j][k].calc_E(1,lat);//calc E for bottom node
      E_crystal += lat[i][j][k].calc_E(2,lat);//calc E for left node
      E_crystal += lat[i][j][k].calc_E(4,lat);//calc E for front node
      E_crystal += lat[i][j][k].calc_E(5,lat);//calc E for back node
      E_crystal += lat[i][j][k].calc_E(3,lat);//calc E for right node
      calc_E_count+=5;
    }


}
//-------------------------------------------------------------------------------//
double Lattice::get_E_crystal()
{
  //we need to do an allreduce with the sum operator
  double temp=0;
  MPI_Allreduce (&E_crystal,&temp, 1, 
  		 MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 //  if(myid==0)
//     cout <<"Energy at process 0 = "<<E_crystal<<" total energy = "<<temp<<endl;
  //we need to multiply this number by the w vector to weight it properly
  cout <<"w_length = "<<w_length<<" total crystal energy " << temp<<endl;
  //  temp = temp * w_length;
  return temp;
}

//-------------------------------decompose the lattice---------------------------
void decompose(const int &num_procs, vector<int> &D)
{
  int dims[3]={0,0,0};

  MPI_Dims_create(num_procs,3,dims);
  for(int i=0;i<3;i++)
    D[i] = dims[i];
};
//--------------------------------get cartesian coords------------------------
void get_cart(const vector<int> &dims, vector<vector<int> > &coords)
{
  vector<int>temp[3];
 
  //start going down x, then y, then z,
  for (int z = 0; z < dims[2]; z++)
    for(int y = 0; y < dims[1]; y++)
      for(int x = 0;x < dims[0]; x++)
	{
	  temp[0].push_back(x);
	  temp[1].push_back(y);
	  temp[2].push_back(z);
	}

  for(int x = 0;x<coords.size();x++)
    {
      coords[x][0] = temp[0][x];
      coords[x][1] = temp[1][x];
      coords[x][2] = temp[2][x];
    } 
}
//----------------------------------find_neighboors--------------------------
void find_neighboors(vector< vector<int> > &nbrs, const int &myid,
		     const vector< vector <int> > & coords, const vector<int> &dims)
{
  nbrs[myid][left  ] = myid - 1;
  nbrs[myid][right ] = myid + 1;

  nbrs[myid][top   ] = myid + dims[0]*dims[1];
  nbrs[myid][bottom] = myid - dims[0]*dims[1];

  nbrs[myid][front ] = myid - dims[0];
  nbrs[myid][back  ] = myid + dims[0];

  //check boundaries
  if(coords[myid][0] == dims[0] - 1)//we are on the right
    nbrs[myid][right] = -1;
  if(coords[myid][0] == 0) //we are on the left
    nbrs[myid][left]  = -1;

  if(coords[myid][1] == 0) //we are on the front
    nbrs[myid][front]  = -1;
  if(coords[myid][1] == dims[1] - 1)  //we are on the back
    nbrs[myid][back] = -1;

  if(coords[myid][2] == dims[2] - 1)  //we are on top
      nbrs[myid][top] = -1;
  if(coords[myid][2] == 0)//we are on the bottom
    nbrs[myid][bottom] = -1;
}
//--------------------------------------------------
void MPE_DECOMP1D(const int &n, const int &numprocs, 
		  const int &myid, int &s, int &e)
{
  int nlocal = n/numprocs;
  int deficit = n%numprocs;
  
  s = myid * nlocal + 1;
  
  if (myid < deficit)
    {
      s+=myid;
      nlocal++;
    }
  else
    s+=deficit;
  e = s + nlocal - 1 ;
  if ((e>n) | (myid==(numprocs-1)))
      e = n;
}
//----------------------------------------------------
void My_send(const int &destination,const int &size,
	     const int &x_len,const int &y_len, const int &z_len,
	     const int &constant_index,
	     vector< vector< vector <Material> > > &lat,
	     MPI_Request &req, vector<int> &node_send)
{
  int count = 0;
  vector<int> temp(6,0);
  node_send.resize(6*size);


  if(destination >= 0)
    {
      //      cout <<"destination " <<destination<<" size "<<6*size<<endl;
      count = 0;

      if(x_len==0)
	for(int j = 1; j<= y_len;j++)
	  for(int k = 1; k<= z_len; k++)
	    {
	      //	      send[count]=lat[constant_index][j][k];
	      lat[constant_index][j][k].get_nodes(temp);
	      for(int ii=0; ii<temp.size(); ii++)
		node_send[6*count+ii] = temp[ii];
	      count++;
	    }
      else if(y_len==0)
	for(int i=1; i<= x_len; i++)
	  for(int k = 1; k<= z_len; k++)
	    {
	      //	      send[count]=lat[i][constant_index][k];
	      lat[i][constant_index][k].get_nodes(temp);
	      for(int ii=0; ii<temp.size(); ii++)
		node_send[6*count+ii] = temp[ii];
	      count++;
	    }
      else
	for(int i=1; i<= x_len; i++)
	  for(int j = 1; j<= y_len;j++)