klmcfacet.f90

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

          location(x,yminus)=0  
        end if
      end if

      if (height(xplus, y)== height(x,y)) then        
        if(type(xplus,y)==1) then                      
          call delete_adatom_event(location(xplus,y)) 
          call add_ledge_event(xplus,y)                
        else                                          
          if (type(xplus,y)==2) then
            x1=dimer_events(1,location(xplus,y),1)
            y1=dimer_events(2,location(xplus,y),1)
            x2=dimer_events(1,location(xplus,y),2)
            y2=dimer_events(2,location(xplus,y),2)
            call delete_dimer_event(location(xplus,y))
            if((x1==xplus) .and. (y1==y)) then  
              call add_ledge_event(x2,y2)
            else
              call add_ledge_event(x1,y1)
            end if    
          else if(type(xplus,y)==3) then 
            call delete_ledge_event(location(xplus,y))
          end if
          type(xplus,y)=0
          location(xplus,y)=0 
        end if
      end if

      if (height(x,yplus)==height(x,y)) then                                                       
        if(type(x,yplus)==1) then                          
          call delete_adatom_event(location(x,yplus))       
          call add_ledge_event(x,yplus)                     
        else
          if(type(x,yplus)==2) then
            x1=dimer_events(1,location(x,yplus),1)
            y1=dimer_events(2,location(x,yplus),1)
            x2=dimer_events(1,location(x,yplus),2)
            y2=dimer_events(2,location(x,yplus),2)
            call delete_dimer_event(location(x,yplus))
            if((x1==x) .and. (y1==yplus)) then  
              call add_ledge_event(x2,y2)
            else
              call add_ledge_event(x1,y1)
            end if         
          else if(type(x,yplus)==3) then
            call delete_ledge_event(location(x,yplus))
          end if
          type(x,yplus)=0
          location(x,yplus)=0
        end if
      end if

  end select update
  
  return

end subroutine local_update



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


subroutine add_adatom_event(x,y)
!*********************************************************************
!    this subroutine adds the possible events of a single adatom to 
!    the table of adatom_events.  
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    --------------> y axis
!    |       1
!    |       o                 x ----- the position of the atom
!    |   2 o x o 4             1,2,3,4----possible moving direction
!    |       o                 o ----- possible ending position
!    |       3
!   \|/
!    x axis
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!                    Author: Zhiyong Wang 01/26/98
!*********************************************************************
  use algorithm_variables
  implicit none

  integer :: x,y

  adatom_number = adatom_number + 1
  type(x,y)=1
  location(x,y)=adatom_number   
  adatom_events(1,adatom_number)=x
  adatom_events(2,adatom_number)=y    
  return

end subroutine add_adatom_event


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


subroutine delete_adatom_event(n)
!*********************************************************************
!    this subroutine removes the possible events of a single adatom  
!    from the table of adatom_events
!    
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    --------------> y axis
!    |       1
!    |       o                 x ----- the position of the atom
!    |   2 o x o 4             1,2,3,4----possible moving direction
!    |       o                 o ----- possible ending position
!    |       3
!   \|/
!    x axis
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!               Author: Zhiyong Wang 01/26/98
!*********************************************************************  
  use algorithm_variables
  implicit none

  integer :: n 
  if (n == adatom_number) then
    adatom_events(1, n) = 0
    adatom_events(2, n) = 0  
  else if (n < adatom_number) then   
    adatom_events(1, n) = adatom_events(1,adatom_number)
    adatom_events(2, n) = adatom_events(2,adatom_number)            
    location(adatom_events(1,n), adatom_events(2,n))=n
  end if
  
  adatom_number = adatom_number - 1
  
  return

end subroutine delete_adatom_event


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

subroutine event_pickup_adatom
!*********************************************************************
!    this subroutine pick up an single adatom event and perform it
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    --------------> y axis
!    |       1
!    |       o              x -------the position of the selected atom
!    |   2 o x o 4          1,2,3,4--possible moving direction
!    |       o              o -------possible ending position
!    |       3                       (possible_event)
!   \|/
!    x axis
!  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!              Author: Zhiyong Wang 01/26/98
!*********************************************************************
  use algorithm_variables
  use input_variables
  implicit none
  
  integer :: selected_atom, possible_event,x,y,xplus,xminus,yplus,yminus, &
             initial_x,initial_y,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)
  selected_atom = int(random*adatom_number)+1
  possible_event=mod(int(random*adatom_number*adatom_motion_number)+1,adatom_motion_number)+1
  x = adatom_events(1,selected_atom)
  y = adatom_events(2,selected_atom)
  height(x,y) = height(x,y)-1
  initial_x=x
  initial_y=y
  call delete_adatom_event(selected_atom)
  type(x,y)=0         !*** change type to immoble ***
  location(x,y)=0
  call boundary_condition(x,y,xplus,xminus,yplus,yminus)

  possible_motion: select case (possible_event)

     case(1)  
        x=xminus
     case(2)
        y=yminus
     case(3)
        x=xplus
     case(4)
        y=yplus

  end select possible_motion

  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_adatom) then 
          remain_number=remain_number-1
       else
          height(initial_x, initial_y)=height(initial_x, initial_y)+1
          call add_adatom_event(initial_x, initial_y)
       end if
    else
       bond_number = bond_counting(x,y)          
                                                 !     o 1
       if(bond_number==0) then                   !    /|\
          call relax(x,y)                        !     |   this figure is the 
       else                                      !     x   example for case 1
          call local_update(x,y)
       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 
       if (random<=facet111_rate_adatom) then
          remain_number=remain_number-1
       else
          height(initial_x,initial_y)=height(initial_x,initial_y)+1
          call add_adatom_event(initial_x, initial_y)
       end if
    else
       bond_number = bond_counting(x,y)          
                                                 !     o 1
       if(bond_number==0) then                   !    /|\
          call relax(x,y)                        !     |   this figure is the 
       else                                      !     x   example for case 1
          call local_update(x,y)
       end if 
    end if
  
  else
    height(x,y) = height(x,y) + 1
    bond_number = bond_counting(x,y)          
                                              !     o 1
    if(bond_number==0) then                   !    /|\
       call relax(x,y)                        !     |   this figure is the 
    else                                      !     x   example for case 1
       call local_update(x,y)
    end if 

  end if    

 return

end subroutine event_pickup_adatom
 

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


subroutine add_dimer_event(x1,y1,x2,y2)
!*********************************************************************
!    this subroutine adds the possible events of a dimer to 
!    the table of dimer_events
!  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   dimer atom #2 o -- o #1 (for the case x1==x2)
!
!   dimer atom #1 o  
!                 |    (for the case y1==y2)
!              #2 o 
!  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
!               Author: Zhiyong Wang 01/26/98
!*********************************************************************
use algorithm_variables
implicit none

  integer :: x1,x2,y1,y2 

  dimer_number = dimer_number+1

  type(x1,y1)=2                             ! * motion #8 *
  type(x2,y2)=2                             !                o #1
  location(x1,y1)=dimer_number              !           then |
  location(x2,y2)=dimer_number              !               \|/
  dimer_events(1,dimer_number,1)=x1         ! move first <-- o #2
  dimer_events(2,dimer_number,1)=y1
  dimer_events(1,dimer_number,2)=x2
  dimer_events(2,dimer_number,2)=y2  

  return                      
                                               
end subroutine add_dimer_event     


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


subroutine delete_dimer_event(n)
!*********************************************************************
!    this subroutine delets the possible events of a dimer from 
!    the table of dimer_events
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   dimer atom #2 o -- o #1 (for the case x1==x2)
!
!   dimer atom #1 o  
!                 |    (for the case y1==y2)
!              #2 o 
!  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!               Author: Zhiyong Wang 01/26/98
!*********************************************************************
  use algorithm_variables
  use input_variables
  implicit none

  integer :: n    
  if (n == dimer_number) then
    dimer_events(1, dimer_number, 1) = 0
    dimer_events(2, dimer_number, 1) = 0
    dimer_events(1, dimer_number, 2) = 0
    dimer_events(2, dimer_number, 2) = 0
  else if (n<dimer_number) then
    dimer_events(1, n, 1) = dimer_events(1,dimer_number,1)
    dimer_events(2, n, 1) = dimer_events(2,dimer_number,1)
    dimer_events(1, n, 2) = dimer_events(1,dimer_number,2)
    dimer_events(2, n, 2) = dimer_events(2,dimer_number,2)
    location(dimer_events(1,dimer_number,1), dimer_events(2,dimer_number,1)) = n
    location(dimer_events(1,dimer_number,2), dimer_events(2,dimer_number,2)) = n
  end if
 
  dimer_number = dimer_number - 1

  return

end subroutine delete_dimer_event


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