klmcfacet.f90

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

      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
       remain_number=remain_number-1
    else
       call local_update(x,y)
    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
       remain_number=remain_number-1
    else
       call local_update(x,y)
    end if
  else
    height(x,y) = height(x,y) + 1
    call local_update(x,y)
  end if

 return

 end subroutine facet_relax



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


subroutine rate_calculation
!********************************************************************
!        the following subroutine calculate the total_rate 
!        and time_step of all the surface events
!  ---------------------------------------------------------------
!               Author: Zhiyong Wang 01/26/98
!********************************************************************

 use input_variables
 use algorithm_variables 
 implicit none

 adatom_rate = adatom_number * adatom_hop_rate
 dimer_rate = dimer_number * dimer_hop_rate
 ledge_rate = ledge_number * ledge_hop_rate
 diffusion_rate = adatom_rate+dimer_rate*0.5+ledge_rate
 if (diffusion_rate /= 0) then
    pa = adatom_rate/diffusion_rate
    pd = 0.5*dimer_rate/diffusion_rate
    pl = ledge_rate/diffusion_rate
 end if
 total_rate = deposition_rate + diffusion_rate
 time_step=1.0/total_rate

 return

end subroutine rate_calculation 


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

function bond_counting(x, y) result(bond)
!********************************************************************
!   this subroutine counts the inplane bond number of a given atom
! -----------------------------------------------------------------
!               Author: Zhiyong Wang 01/26/98  
!********************************************************************

  use algorithm_variables
  implicit none
 
  integer :: x,y,bond,xplus,xminus,yplus,yminus
  
  bond = 0
  call boundary_condition(x,y,xplus,xminus,yplus,yminus)

  if (height(x,yminus) >= height(x,y)) bond = bond + 1
  if (height(xminus,y) >= height(x,y)) bond = bond + 1
  if (height(x,yplus) >= height(x,y)) bond = bond + 1
  if (height(xplus,y) >= height(x,y)) bond = bond + 1

end function bond_counting



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


subroutine boundary_condition(x,y,xp,xm,yp,ym)
!********************************************************************
!    this subroutine set perodic boundary conditions
!    xm----xminus, xp----xplus, ym----yminus, yp----yplus
! -----------------------------------------------------------------
!              Author: Zhiyong Wang 01/26/98 
!********************************************************************
  use input_variables

  implicit none
 
  integer,intent(in) :: x,y
  integer,intent(out):: xp,xm,yp,ym 
   
  xp = x+1
  xm = x-1
  yp = y+1
  ym = y-1
  if (facet=="none") then
     xp = mod(xp,width+1)+int(xp/(width+1))
     xm = mod(xm,width)+width*int((width-xm)/width) 
     yp = mod(yp,length+1)+int(yp/(length+1))
     ym = mod(ym,length)+length*int((length-ym)/length)
  end if

  return

end subroutine boundary_condition


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


subroutine local_update(x, y)
!********************************************************************
!  this subroutine updates the local mobile event after an atom is 
!  deposited or diffused
! -----------------------------------------------------------------
!  ------------------> y axis
!  |                    
!  |        1                 1,2,4,8---bond_case_number increment
!  |        o                 o---local neighbors for a deposited atom 
!  |    2 o x o 8             x---deposited atom (x,y) 
!  |        o 
!  |        4
! \|/
!  x axis
! ----------------------------------------------------------------
!                 Author: Zhiyong Wang 01/26/98
!********************************************************************
  use algorithm_variables
  implicit none

  integer :: bond_case_number,x,y,xplus,xminus,yplus,yminus,x1,x2,y1,y2

  bond_case_number = 0
  call boundary_condition(x,y,xplus,xminus,yplus,yminus)
  if (height(xminus,y)>=height(x,y)) bond_case_number=bond_case_number+1
  if (height(x,yminus)>=height(x,y)) bond_case_number=bond_case_number+2 
  if (height(xplus,y)>=height(x,y)) bond_case_number=bond_case_number+4
  if (height(x,yplus)>=height(x,y)) bond_case_number=bond_case_number+8

  update: select case(bond_case_number)
          
!     case(0)                                           !       |
!       call add_adatom_event(x,y)                      !    -- x --
                                                        !       |
     case(1)
      if(type(x,y)==3) call delete_ledge_event(location(x,y)) 
      if (height(xminus, y)== height(x,y)) then                        
                                                       !     1 o #1 of dimer
        if(type(xminus,y)==1) then                     !       |
          call delete_adatom_event(location(xminus,y)) !    -- x --
          call add_dimer_event(xminus,y,x,y)           !       |
        else 
          if (type(xminus,y)==2) then
            x1=dimer_events(1,location(xminus,y),1)
            y1=dimer_events(2,location(xminus,y),1)
            x2=dimer_events(1,location(xminus,y),2)
            y2=dimer_events(2,location(xminus,y),2)
            call delete_dimer_event(location(xminus,y))
            if((x1==xminus) .and. (y1==y)) then  
              call add_ledge_event(x2,y2)
            else             
              call add_ledge_event(x1,y1)
            end if    
          else if(type(xminus,y)==3) then 
            call delete_ledge_event(location(xminus,y))
          end if
          call add_ledge_event(x,y)
          type(xminus,y)=0
          location(xminus,y)=0        
        end if

      else

        call add_ledge_event(x,y)

      end if  
    
     case(2)
      if(type(x,y)==3) call delete_ledge_event(location(x,y)) 
      if (height(x,yminus)==height(x,y)) then       
                                                
        if(type(x,yminus)==1) then                            !        |
          call delete_adatom_event(location(x,yminus))        ! 2 o -- x --
          call add_dimer_event(x,y,x,yminus)                  !        |
        else
          if(type(x,yminus)==2) then
            x1=dimer_events(1,location(x,yminus),1)
            y1=dimer_events(2,location(x,yminus),1)
            x2=dimer_events(1,location(x,yminus),2)
            y2=dimer_events(2,location(x,yminus),2)
            call delete_dimer_event(location(x,yminus))
            if((x1==x) .and. (y1==yminus)) then  
              call add_ledge_event(x2,y2)
            else
              call add_ledge_event(x1,y1)
            end if     
          else if(type(x,yminus)==3) then
            call delete_ledge_event(location(x,yminus))
          end if
          call add_ledge_event(x,y) 
          type(x,yminus)=0
          location(x,yminus)=0  
        end if

      else

        call add_ledge_event(x,y)

      end if  

     case(4)               
      if(type(x,y)==3) call delete_ledge_event(location(x,y))
 
      if (height(xplus, y)== height(x,y)) then                             
        if(type(xplus,y)==1) then                      !       |
          call delete_adatom_event(location(xplus,y))  !    -- x --
          call add_dimer_event(x,y,xplus,y)            !       |
        else                                           !       o 4
          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
          call add_ledge_event(x,y)
          type(xplus,y)=0
          location(xplus,y)=0 
        end if

      else

        call add_ledge_event(x,y)

      end if  

     case(8)
      if(type(x,y)==3) call delete_ledge_event(location(x,y)) 
      if (height(x,yplus)==height(x,y)) then       
                                                
        if(type(x,yplus)==1) then                      !    |    8
          call delete_adatom_event(location(x,yplus))  ! -- x -- o #1 of dimer
          call add_dimer_event(x,yplus,x,y)            !    |
        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
          call add_ledge_event(x,y) 
          type(x,yplus)=0
          location(x,yplus)=0
        end if

      else

        call add_ledge_event(x,y)

      end if  

     case(3,5,6,7,9,10,11,12,13,14,15)

      if (height(xminus, y)== height(x,y)) then        !       1 o 
        if(type(xminus,y)==1) then                     !         |
          call delete_adatom_event(location(xminus,y)) !  2 o -- x -- o 8
          call add_ledge_event(xminus,y)               !         |
        else                                           !         o 4
          if (type(xminus,y)==2) then
            x1=dimer_events(1,location(xminus,y),1)
            y1=dimer_events(2,location(xminus,y),1)
            x2=dimer_events(1,location(xminus,y),2)
            y2=dimer_events(2,location(xminus,y),2)
            call delete_dimer_event(location(xminus,y))
            if((x1==xminus) .and. (y1==y)) then  
              call add_ledge_event(x2,y2)
            else       
              call add_ledge_event(x1,y1)
            end if    
          else if(type(xminus,y)==3) then 
            call delete_ledge_event(location(xminus,y))
          end if
          type(xminus,y)=0
          location(xminus,y)=0 
        end if
      end if

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