klmcfacet.f90

用fortran和C写的薄膜生长模拟程序

    height(1:int(0.2*width), 1:length)=1
    height(int(0.2*width)+1:int(0.4*width), 1:length)=2
    height(int(0.4*width)+1:int(0.6*width), 1:length)=3 
    height(int(0.6*width)+1:int(0.8*width), 1:length)=4 
    height(int(0.8*width)+1:width, 1:length)=5
!    do i=1, int(0.*width*length)
!      number=number+1
!      call random_number(random)   
!      x = int(random*width) + 1     
!     y = mod(int(random*length*width),length) + 1
!     height(x,y)=height(x,y)-1
!      if (height(x,y)<0) then
!         height(x,y)=0
!         number=number-1
!      end if 
!    end do   
    initial_atom_number= (int(0.4*width)-int(0.2*width))*length*15   !  &
               !         +(int(0.6*width)-int(0.4*width))*length*3     &
               !         - number
  else    !!!!!!!! this feature finally completed on 8/12/98
    open(9,file=in_structure)
    read(9,*)w,l,initial_atom_number, initial_total 
    read_file: do
       if((width==w) .and.(length==l)) exit
       close(9,status="keep")
       print *, "The size of the structure is not right,"
       print *, "please enter a file name with the stucture of size ",width,"*",length   
       read *, in_structure
       open(9,file=in_structure)
       read(9,*)w,l,initial_atom_number
    end do read_file 
    do i=1, width
      do j=1, length
        read(9,*) height(i,j), type(i,j), location(i,j)
        if (type(i,j)==1) then
           adatom_number = adatom_number + 1
           adatom_events(1,location(i,j))=i
           adatom_events(2,location(i,j))=j
        else if ( type(i,j)==3 ) then
           ledge_number = Ledge_number + 1
           ledge_events(1,location(i,j))=i
           ledge_events(2,location(i,j))=j
        end if
      end do
    end do
    read(9,*)dimer_number
    if (dimer_number/=0) then
       do i=1, dimer_number
          read(9,*) dimer_events(1,i,1),dimer_events(2,i,1),dimer_events(1,i,2),dimer_events(2,i,2)   
       end do
    end if

  end if 

  remain_number=initial_atom_number
  
  return

end subroutine initialization


!       ====== * ====== * ====== * ====== * ====== * ======


subroutine klmc_algorithm
!*********************************************************************
!  the following subroutine is the key algorithm of klmc simulation
! ------------------------------------------------------------------
!               Author: Zhiyong Wang 01/26/98
!*********************************************************************

use algorithm_variables
use input_variables
implicit none

real :: random
 
if (facet=="110") then
  do while (total_number <int(monolayer*atoms_monolayer)) 
    steps=steps+1 
    call rate_calculation
    call random_number(random)
    if (random <= deposition_rate/total_rate) then
       call deposition
       time = time + time_step       
    else  
       call random_number(random)
       do while (random==0.0)
          call random_number(random)
       end do
       if ((random > 0.0) .and. (random <= pa)) then
          call event_pickup_adatom 
       else if ((random > pa) .and. (random <= (pa + pd))) then
          call event_pickup_dimer
       else if (random>pa+pd) then
          call event_pickup_ledge
       end if 
       time=time+time_step
    end if
  end do

else 
  do while (total_number < int(monolayer*width*length)) 
    steps=steps+1 
    call rate_calculation
    call random_number(random)
    if (random <= deposition_rate/total_rate) then
       call deposition
       time = time + time_step       
    else  
       call random_number(random)
       do while (random==0.0)
         call random_number(random)
       end do
       if ((random > 0.0) .and. (random <= pa)) then
          call event_pickup_adatom 
       else if ((random > pa) .and. (random <= (pa + pd))) then
          call event_pickup_dimer
       else if (random>pa+pd) then
          call event_pickup_ledge
       end if  
       time=time+time_step
    end if
  end do

end if  

return  

end subroutine klmc_algorithm


!       ====== * ====== * ====== * ====== * ====== * ======


subroutine deposition
!*********************************************************************
!       the following subroutine deposits an atom on surface 
!       and relaxes it to a stable site
! ------------------------------------------------------------------
!                  Author: Zhiyong Wang 01/26/98
!*********************************************************************

use algorithm_variables 
use input_variables
implicit none

 integer :: x,y,y1,y2,selected_pos,x1,x2,x3,x4
 real :: random

 interface
  function bond_counting(x, y) result(bond)
   integer::x,y,bond
  end function bond_counting
 end interface

 call random_number(random)
 x=int(1.6)
 if (facet=="110") then
     selected_pos=int(random*atoms_monolayer)+1
     if (mod(selected_pos,(2*facet110_div_2-1))==0) then
        x=int(selected_pos/(2*facet110_div_2-1))+facet110_div_2-1
        y=width-int(selected_pos/(2*facet110_div_2-1))
     else if (mod(selected_pos,(2*facet110_div_2-1)) <= facet110_div_2) then
        x=int(selected_pos/(2*facet110_div_2-1))+mod(selected_pos,(2*facet110_div_2-1))
        y=facet110_div_1-int(selected_pos/(2*facet110_div_2-1))+mod(selected_pos,(2*facet110_div_2-1))-1
     else
        x=int(selected_pos/(2*facet110_div_2-1))+1+mod(selected_pos,(2*facet110_div_2-1))-facet110_div_2
        y=facet110_div_1-int(selected_pos/(2*facet110_div_2-1))+mod(selected_pos,(2*facet110_div_2-1))-facet110_div_2-1
     end if
 else  
     x = int(random*width) + 1 
     y = mod(int(random*length*width),length) + 1
 end if
   
 height(x, y) = height(x, y) + 1
 if (facet=="110") then
   if (layer=="yes") then 
     x1=facet110_div_1+height(x,y)
     x4=facet110_div_2+width-height(x,y)+1
     if (mod(height(x,y),2)==0) then
        x2=1-(facet110_div_1-height(x,y))
        x3=facet110_div_1-height(x,y)+1
     else
        x3=facet110_div_1-height(x,y)
        x2=1-(facet110_div_1-height(x,y)+1)
     end if 
   else
     x1=facet110_div_1+1
     x4=facet110_div_2+width
     x2=1-facet110_div_1
     x3=facet110_div_1-1
   end if

   if (y<-x+x1 .or. y<x+x2 .or. y>x+x3 .or. y>-x+x4) then
       height(x,y)=height(x,y)-1
       partial_flux=partial_flux+1
   else
       total_number=total_number+1
       remain_number=remain_number+1
       bond_number=bond_counting(x,y) 
       if (bond_number == 0) then
          call relax(x,y)
       else
          call local_update(x,y)
       end if
   end if 

 else if (facet=="111") then
   if (layer=="yes") then 
     if (mod(height(x,y),2)==0) then
       x3=int(height(x,y)/2)
       x1=int(height(x,y)/2)+1
       x2=length-int(height(x,y)/2)+1
       x4=width-int(height(x,y)/2)
     else
       x3=int(height(x,y)/2)+1
       x1=int(height(x,y)/2)+1
       x2=length-int(height(x,y)/2)
       x4=width-int(height(x,y)/2)
     end if
   else
     x1=1
     x2=length
     x3=1
     x4=width
   end if

   if (y<x1 .or. y>x2 .or. x<x3 .or. x>x4) then
       height(x,y)=height(x,y)-1
       partial_flux=partial_flux+1
   else
       total_number=total_number+1
       remain_number=remain_number+1
       bond_number=bond_counting(x,y) 
       if (bond_number == 0) then
          call relax(x,y)
       else
          call local_update(x,y)
       end if
   end if

 else
    total_number=total_number+1
    remain_number=remain_number+1
    bond_number=bond_counting(x,y) 
    if (bond_number == 0) then
       call relax(x,y)
    else
       call local_update(x,y)
    end if 

 end if

 return

end subroutine deposition


!       ====== * ====== * ====== * ====== * ====== * ======


subroutine relax(x,y)
!*********************************************************************
!     this subroutine relaxes an uper layer atom to a more stable 
!     lower layer site
! ------------------------------------------------------------------
!   ----------------> y axis
!   |        /(x,y) for even plane 
!   |  1 o ,/                              
!   |       x                              
!   |  4 o    o 2      
!   |
!   |        1 o    o4
!  \|/           x 
!              /'   o2
!   x axis    /(x,y) for odd plane  
!
! , & ' moving direction for atom jump form upper layer to lower layer
!
! 1,2,4 case number increment, increse the number if the site is empty 
!
! x --- upper layer atom 
! o --- lower layer atom
! ------------------------------------------------------------------
!                Author: Zhiyong Wang 01/26/98       
!*********************************************************************  
  
  use algorithm_variables         
  implicit none                    
                                   
  integer :: case_number,x,y,xplus,xminus,yplus,yminus               
  real :: random
                                  
  case_number=0                   
  call boundary_condition(x,y,xplus,xminus,yplus,yminus)
         
  if ( mod(height(x,y),2)== 1 ) then 
    xplus=xminus
    yminus=yplus
  end if
                             
  if (height(x,yminus) < height(x,y)-1) case_number = case_number+1 
  if (height(xplus,y) < height(x,y)-1) case_number = case_number+2
  if (height(xplus,yminus)<height(x,y)-1) case_number = case_number+4 
  
  if (case_number==0) then
    call add_adatom_event(x,y) 
  else
    height(x,y) =height(x,y) - 1 
    relaxation : select case (case_number) 
         case(1)  
             call facet_relax(x,yminus)
         case(2)
	     call facet_relax(xplus,y)
         case(3,7)
             call random_number(random)
             if (random<0.5) then
                call facet_relax(x,yminus)
             else
	        call facet_relax(xplus,y)
             end if
         case(4)
	     call facet_relax(xplus,yminus)
         case(5)
             call random_number(random)
             if (random<0.5) then
                call facet_relax(xplus,yminus)
             else
	        call facet_relax(x,yminus)
             end if 
         case(6)
             call random_number(random)
             if (random<0.5) then
	        call facet_relax(xplus,yminus)
             else
	        call facet_relax(xplus,y)
             end if
    end select relaxation
  end if

 return

 end subroutine relax



subroutine facet_relax(x,y)
!********************************************************************
!                For Facet Growth Rate
!-------------------------------------------------------------------     
!                Author: Zhiyong Wang
!********************************************************************
 use input_variables
 use algorithm_variables 
 implicit none

 integer::x,y,x1,x2,x3,x4

  if (facet=="110") then
    height(x,y)=height(x,y)+1
    if (layer=="yes") then 
      x1=facet110_div_1+height(x,y)