main.cpp

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

#include "mpi.h"
#include <iostream>
#include <cstdlib> //for the random number generator
#include <ctime>
#include <fstream>
#include <vector>
#include "printvector.h"
#include "material.h"
#include "lattice.h"

using namespace std;
using namespace MPI;

//-----------------------------------------------------------------------------------------------//
int main(int argc,char **argv)
{
  Init(argc,argv);
  int myid = COMM_WORLD.Get_rank();
  int numprocs = COMM_WORLD.Get_size();
  

  //initialize the matrix
  vector<int> temp(6);
  int x=10,y=10,z=10;//100x100x100 lattice
  int temperature = 273;
  int t_start = 273;
  int t_finish = 300;
  int settle = 1000;
  int i=0,j=0,k=0,proc=-1,node=0;
  double E_tot=0;
  int count=0;
  Lattice lat(x,y,z,numprocs,myid,temp);

  ofstream data ("energy.dat",ios::out);//open a file for output
  if (!data)
    cout <<"WARNING!!!  couldn't open the energy.dat file"<<endl;
  ofstream T ("temperature.dat",ios::out);//open a file for output
  if (!T)
    cout <<"WARNING!!!  couldn't open the energy.dat file"<<endl;
  int test = lat.Init();//initialize the cube
  if(test == 0)
    {
      lat.exchange();    //exchange our data with our neighboors
      lat.inner_comps(); //we should do our inner computation here
      lat.wait();        //wait for our exchange to occur if needed
      lat.outer_comps(); //now do the comps on the outer walls
      lat.set_temperature(temperature);//set the current lattice temperature

      for(temperature = t_start; temperature < t_finish; temperature++)
	{
	  for(int ii = 0; ii< settle; ii++)
	    {
	      lat.get_W();//calculate W vector
	      lat.rand_flip(i,j,k,proc,node);//pick random number/find and flip the node
	      
	      //	      if(myid==proc)
	      //		cout <<"we flipped node " <<node<<" at ijk = " <<i <<" "<<j<<" "<<k<<" on process "<<proc<<endl;
	      lat.adjust_E_crystal(i,j,k,proc,node);//adjust the crystal energy since we just flipped the node
	      proc = -1;	      
	      double temp = lat.get_E_crystal();
	      E_tot += temp;
 	      if(myid==0)
 		{
 		  cout <<"iteration "<<ii<<" \r";
 		  data <<temp<<endl;
		}
	      COMM_WORLD.Barrier();
	    }
	  lat.set_temperature(temperature);//adjust the temperature
	  E_tot /= (settle-1);
 	  if(myid==0)
 	    T << E_tot<<endl;
	  E_tot = 0;
	}

    }
  Finalize();
  return 0;
}