klmcfacet.f90

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

subroutine event_pickup_dimer
!*********************************************************************
!       this subroutine pick up a dimer event and perform it.
!       after one of the dimer atoms is moved to a new site, the 
!       atom is checked that it will relax or not. if it would 
!       relax, then relax it and update another dimer atom; if it 
!       would not relax, let another dimer atom catch up.         
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   dimer atom #2 o -- o #1 (for the case x1==x2)
!
!   dimer atom #1 o  
!                 |    (for the case y1==y2)
!              #2 o 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   the structure of the table read as:
!        dimer_events(1:2,1:dimer_number,1:2)
!                     /        |          \
!                  x,y  the sequency    dimer atom #1 or #2 
!                         of dimer    
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!              Author: Zhiyong Wang 01/26/98
!*********************************************************************
 use algorithm_variables
 use input_variables
 implicit none
  
 integer :: selected_atom,possible_event,case_number,x1,y1,x2,y2,h1,h2, &
            xplus,xminus,yplus,yminus,initial_x1,initial_x2,initial_y1,initial_y2
 real :: random
 interface
  function bond_counting(x, y) result(bond)
   integer::x,y,bond
  end function bond_counting
 end interface
 case_number=0 
 call random_number(random)
 selected_atom = int(random * dimer_number) + 1
 possible_event = mod(int(random*dimer_number*dimer_motion_number)+1,dimer_motion_number)+1
 x1=dimer_events(1, selected_atom,1)
 y1=dimer_events(2, selected_atom,1)
 x2=dimer_events(1, selected_atom,2)
 y2=dimer_events(2, selected_atom,2)
 initial_x1=x1
 initial_x2=x2
 initial_y1=y1
 initial_y2=y2 
 h1=height(x1,y1)
 h2=height(x2,y2)
 call delete_dimer_event(selected_atom)
 type(x1,y1)=0
 location(x1,y1)=0
 type(x2,y2)=0
 location(x2,y2)=0          
 possible_motion: select case (possible_event)

   case(1) 
      height(x1,y1)=height(x1,y1)-1 
      call boundary_condition(x1,y1,xplus,xminus,yplus,yminus)
      if (y1==y2) then
        xminus=x1
        y1=yminus
      end if 
      call facet_dimer(initial_x1,initial_y1,xminus,y1,x2,y2,initial_x1,initial_y1, &
                       initial_x2,initial_y2,1)

   case(2) 
      height(x2,y2)=height(x2,y2)-1    
      call boundary_condition(x2,y2,xplus,xminus,yplus,yminus)
      if (y1==y2) then                               !first  
        xminus=x2                                    ! /|\ /|\       <--- o#1
        y2=yminus                                    !  |   |   
      end if                                         !#2o   o#1 first<--- o#2
 
      call facet_dimer(initial_x2,initial_y2,xminus,y2,x1,y1,initial_x1,initial_y1, &
                       initial_x2,initial_y2,2)
 
   case(3) 
      height(x1,y1)=height(x1,y1)-1   
      call boundary_condition(x1,y1,xplus,xminus,yplus,yminus)
      if (y1==y2) then
        xminus=x1
        y1=yminus
      end if                !     /|\ first <--- o #1
                            !#2    |             | 
                            !o --> o #1          o #2       
      call facet_dimer(initial_x1,initial_y1,xminus,y1,x2,y2,initial_x1,initial_y1, &
                       initial_x2,initial_y2,3)
    
   case(4) 
      height(x2,y2)=height(x2,y2)-1
      call boundary_condition(x2,y2,xplus,xminus,yplus,yminus)
      if (y1==y2) then                             !first  
        xminus=x2                                  !   /|\            o #1
        y2=yminus                                  !    |    #1       |
      end if                                       ! #2 o<--- o  <--- o #2
                                                   !             first
      call facet_dimer(initial_x2,initial_y2,xminus,y2,x1,y1,initial_x1,initial_y1, &
                       initial_x2,initial_y2,4)

   case(5) 
      height(x1,y1)=height(x1,y1)-1 
      call boundary_condition(x1,y1,xplus,xminus,yplus,yminus)
      if (y1==y2) then                              !   o   o
        xplus=x1                                    !   |   | first  o ---> #1
        y1=yplus                                    !  \|/ \|/     
      end if                                        !#2       #1     o ---> #2

      call facet_dimer(initial_x1,initial_y1,xplus,y1,x2,y2,initial_x1,initial_y1, &
                       initial_x2,initial_y2,5)

   case(6) 
      height(x2,y2)=height(x2,y2)-1
      call boundary_condition(x2,y2,xplus,xminus,yplus,yminus)
      if (y1==y2) then                            !first o   o
        xplus=x2                                  !      |   |        o --> #1
        y2=yplus                                  !     \|/ \|/     
      end if                                      !   #2     #1 first o --> #2

      call facet_dimer(initial_x2,initial_y2,xplus,y2,x1,y1,initial_x1,initial_y1, &
                       initial_x2,initial_y2,6)

   case(7) 
      height(x1,y1)=height(x1,y1)-1      
      call boundary_condition(x1,y1,xplus,xminus,yplus,yminus)
      if (y1==y2) then                         !   o --> o
        xplus=x1                               !         | first o --> #1
        y1=yplus                               !        \|/      |     
      end if                                   !#2        #1     o #2

      call facet_dimer(initial_x1,initial_y1,xplus,y1,x2,y2,initial_x1,initial_y1, &
                       initial_x2,initial_y2,7)

   case(8) 
      height(x2,y2)=height(x2,y2)-1
      call boundary_condition(x2,y2,xplus,xminus,yplus,yminus)
      if (y1==y2) then                     !first o <-- o
        xplus=x2                           !      |         o #1
        y2=yplus                           !     \|/        |     
      end if                               !     #2   first o --> #2

      call facet_dimer(initial_x2,initial_y2,xplus,y2,x1,y1,initial_x1,initial_y1, &
                       initial_x2,initial_y2,8)

 end select possible_motion

 if (height(initial_x1,initial_y1)<h1) then
    bond_number=bond_counting(initial_x1,initial_y1)
    if (bond_number==1) then
       call add_ledge_event(initial_x1,initial_y1)
    end if
 end if

 if (height(initial_x2,initial_y2)<h2) then
    bond_number=bond_counting(initial_x2,initial_y2)
    if (bond_number==1) then
       call add_ledge_event(initial_x2,initial_y2)
    end if
 end if

 return

end subroutine event_pickup_dimer



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



subroutine facet_dimer(initial_x,initial_y,x,y,stay_x,stay_y,initial_x1, &
                       initial_y1,initial_x2,initial_y2,case)
!*************************************************************************
! this subroutine just make the code easier to write when add the facet 
! feature
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!                Author: Zhiyong Wang
!*************************************************************************
 use algorithm_variables
 use input_variables
 implicit none

 integer:: x,y,stay_x,stay_y,initial_x1,initial_x2,initial_y1,initial_y2,  &
           case_number,xminus,xplus,yminus,yplus,initial_x,initial_y,case, &
           x1,x2,x3,x4
 real:: random

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

 case_number=0

 if (facet=="110") then
    height(x,y)=height(x,y)+1
    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
       call random_number(random)
       if (random<=facet110_rate_dimer) then 
          remain_number=remain_number-1
          call add_adatom_event(stay_x, stay_y)
       else
          height(initial_x, initial_y)=height(initial_x, initial_y)+1
          call add_dimer_event(initial_x1, initial_y1, initial_x2, initial_y2)
       end if
    else
                                            ! /|\ /|\ first<--- o#1
       bond_number=bond_counting(x,y)       !  |   |     
                                            !#2o   o #1    <--- o#2  
       if(bond_number==0) then                                                         
         call relax_dimer(x,y,stay_x,stay_y,initial_x1,initial_y1,initial_x2,initial_y2,case)
       else                          ! here case is possible event number                                                  
         call local_update(x,y)
         call boundary_condition(stay_x,stay_y,xplus,xminus,yplus,yminus)
         if ( mod(height(stay_x,stay_y),2)== 1 ) then 
           xplus=xminus
           yminus=yplus
         end if
         if (height(stay_x,yminus) < height(stay_x,stay_y)-1) case_number = case_number+1 
         if (height(xplus,stay_y) < height(stay_x,stay_y)-1) case_number = case_number+2
         if (height(xplus,yminus)<height(stay_x,stay_y)-1) case_number = case_number+4
         if (case_number==0) then        
            call dimer_update(case,initial_x1,initial_y1,initial_x2,initial_y2) 
         else                          ! here case is possible event number
            call relax(stay_x,stay_y)
         end if           
       end if
    end if
       
 else if (facet=="111") then
    height(x,y)=height(x,y)+1
    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
       call random_number(random)
       if (random<=facet111_rate_dimer) then 
          remain_number=remain_number-1
          call add_adatom_event(stay_x,stay_y)
       else
          height(initial_x,initial_y)=height(initial_x,initial_y)+1
          call add_dimer_event(initial_x1, initial_y1, initial_x2, initial_y2)
       end if
    else
                                            ! /|\ /|\ first<--- o#1
       bond_number=bond_counting(x,y)       !  |   |     
                                           !#2o   o #1    <--- o#2  
       if(bond_number==0) then                                                         
         call relax_dimer(x,y,stay_x,stay_y,initial_x1,initial_y1,initial_x2,initial_y2,case)
       else                          ! here case is possible event number                                                  
         call local_update(x,y)
         call boundary_condition(stay_x,stay_y,xplus,xminus,yplus,yminus)
         if ( mod(height(stay_x,stay_y),2)== 1 ) then 
            xplus=xminus
            yminus=yplus
         end if
         if (height(stay_x,yminus) < height(stay_x,stay_y)-1) case_number = case_number+1 
         if (height(xplus,stay_y) < height(stay_x,stay_y)-1) case_number = case_number+2
         if (height(xplus,yminus)<height(stay_x,stay_y)-1) case_number = case_number+4
         if (case_number==0) then        
            call dimer_update(case,initial_x1,initial_y1,initial_x2,initial_y2) 
         else                          ! here case is possible event number
            call relax(stay_x,stay_y)
         end if           
       end if
    end if 
 else
    height(x,y)=height(x,y)+1            ! /|\ /|\ first<--- o#1
    bond_number=bond_counting(x,y)       !  |   |     
                                         !#2o   o #1    <--- o#2  
    if(bond_number==0) then                                                         
      call relax_dimer(x,y,stay_x,stay_y,initial_x1,initial_y1,initial_x2,initial_y2,case)
    else                          ! here case is possible event number                                                  
      call local_update(x,y)
      call boundary_condition(stay_x,stay_y,xplus,xminus,yplus,yminus)
      if ( mod(height(stay_x,stay_y),2)== 1 ) then 
         xplus=xminus
         yminus=yplus
      end if
      if (height(stay_x,yminus) < height(stay_x,stay_y)-1) case_number = case_number+1 
      if (height(xplus,stay_y) < height(stay_x,stay_y)-1) case_number = case_number+2
      if (height(xplus,yminus)<height(stay_x,stay_y)-1) case_number = case_number+4
      if (case_number==0) then        
         call dimer_update(case,initial_x1,initial_y1,initial_x2,initial_y2) 
      else                          ! here case is possible event number
         call relax(stay_x,stay_y)
      end if  
    end if
 end if    

 return

 end subroutine facet_dimer


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

  
subroutine relax_dimer(mx1,my1,mx2,my2,x1,y1,x2,y2,n)
!*********************************************************************
!     this subroutine relaxes one or both  of the uper layer dimer 
!     atoms to a more stable lower layer site after one of them was 
!     diffused, and updates the neighbors
! -------------------------------------------------------------------
!    mx1----the new xposition after the fist-motion atom was moved
!    my1----the new yposition after the fist-motion atom was moved
!    mx2----the xposition of one of the dimer atom which did not move 
!    my2----the yposition of one of the dimer atom which did not move
!    n------the nth dimer event 
! -------------------------------------------------------------------
!               Author: Zhiyong Wang 01/26/98
!*********************************************************************  
  
  use algorithm_variables         
  implicit none                    
                                   
  integer :: case_number,mx1,mx2,my1,my2,x1,y1,x2,y2,n,xplus,xminus,yplus,yminus                
                                         
  case_number=0                   
  call boundary_condition(mx1,my1,xplus,xminus,yplus,yminus)
         
  if ( mod(height(mx1,my1),2)== 1 ) then 
    xplus=xminus
    yminus=yplus
  end if
                              
  if (height(mx1,yminus) < height(mx1,my1)-1) case_number = case_number+1 
  if