klmcfacet.f90

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

! $Id: klmcfacet.f90,v 1.2 1998/01/28 23:25:21 zywang Exp zywang $:
!         ***********************************************
!             KLMC SIMULATION OF COPPER DEPOSITION ON
!                COPPER (100) SURFACE (Version 1.1)
!         ***********************************************                   
!                               by
!
!           Zhiyong Wang, Youhong Li and James B. Adams
!                     Arizona State University
!
!                        Acknowledgements
!
!                       David F. Richards
!                             NCSA                                     
!                     last update: 05/23/99
!         ***********************************************   
! @ This is the test version of KLMC--AWL. All Rights Reserved.    


module input_variables
!**********************************************************************
!  this module defines the input and basic variables which only are 
!  calculated once 
! ------------------------------------------------------------------
!            Author: Zhiyong Wang  03/01/98
!**********************************************************************

  implicit none

!**********************************************************************
! adatom_energy_barrier--------activation energy for the diffusion of 
!                              an isolated adatom (eV)
! adatom_hop_rate--------------the hopping rate of isolated adatoms
! adatom_motion_number---------possible motions of an isolated adatom
! adatom_prefactor-------------prefactor for the diffusion of isolated 
!                              adatoms (Hz) 
! deposition_flux--------------deposition flux (monolayer/sec)
! deposition_rate--------------deposition rate (atoms/sec)
! dimer_energy_barrier---------activation energy for the diffusion of 
!                              dimer (eV)
! dimer_hop_rate---------------the hopping rate of dimers
! dimer_motion_number----------possible motions of a dimer
! dimer_prefactor--------------prefactor for the diffusion of the 
!                              dimers (Hz)
! initial_atom_number----------the atom number of a initial structure
!                              above a flat surface   
! length,  width---------------simulation dimension 
! ledge_motion_number----------possible motions of a ledge atom
! simulation_time--------------simulation time (sec)
! sub_ad_barrier---------------activation energy for the diffusion of 
!                              an isolated adatom on substrate(eV)
! sub_ad_prefactor-------------prefactor for the diffusion of isolated 
!                              adatoms on substrate (Hz) 
! sub_ad_rate------------------the hopping rate of isolated adatoms on 
!                              substrate
! sub_di_barrier---------------activation energy for the diffusion of 
!                              a dimer on substrate(eV)
! sub_di_prefactor-------------prefactor for the diffusion of dimer 
!                              on substrate(Hz) 
! sub_di_rate------------------the hopping rate of dimers on substrate
! sub_material-----------------substrate material
! sub_up_barrier---------------the energy barrier for the diffusion of 
!                              a ledge adatom climbing up from substrate
!                              to an upper layer   
! temperature------------------temperature(k)
! up_barrier-------------------the energy barrier for the diffusion of 
!                              a ledge adatom climbing up to an upper
!                              layer 
!**********************************************************************
! date_time--------------------the date & time of running the simulation
! big_ben----------------------the parameter used by intrinsic subroutine
!                              date_and_time 
!**********************************************************************

  integer :: length,width,adatom_motion_number=4,dimer_motion_number=8,&
          &  ledge_motion_number=2, deposition_rate, date_time(8),     &
          &  initial_atom_number

  real :: deposition_flux, adatom_energy_barrier, temperature,        &
       & dimer_energy_barrier, adatom_prefactor, dimer_prefactor,     &
       & adatom_hop_rate,dimer_hop_rate,sub_ad_rate, sub_di_rate,     &
       & sub_ad_prefactor,sub_di_prefactor,sub_ad_barrier,up_barrier, &
       & sub_di_barrier,sub_up_barrier,facet110_barrier_adatom,       &
       & facet110_barrier_dimer,facet110_barrier_ledge,ledge_hop_rate,&  
       & facet111_barrier_adatom,facet111_barrier_dimer, monolayer,   &
       & ledge_prefactor,facet111_barrier_ledge,ledge_energy_barrier,&
       & facet110_rate_adatom,facet110_rate_dimer,facet110_rate_ledge,&
       & facet111_rate_adatom,facet111_rate_dimer,facet111_rate_ledge,&
       & facet110_back_barrier,facet111_back_barrier, initial_total,  &
       & facet110_self_barrier,facet111_self_barrier  
    
  real, parameter :: kb=8.617e-5    !*** boltzmann constant in eV/K ***

  character (20) output_file, output_file_xyz, in_structure, restart, &
               & restart_file,facet,facet_plot,layer,output_file_pdb
  character (2)  material
  character (5)  sub_material
  character (10) big_ben(3)
  namelist /input/ deposition_flux,adatom_energy_barrier,sub_material, &
            & dimer_energy_barrier, temperature, adatom_prefactor,     &
            & dimer_prefactor,length, width,ledge_prefactor,facet_plot,&
            & adatom_motion_number, dimer_motion_number, material,     &
            & ledge_motion_number, output_file, output_file_xyz,       &
            & in_structure, sub_ad_prefactor, sub_di_prefactor,        &
            & sub_ad_barrier,sub_di_barrier,sub_up_barrier,up_barrier, &
            & restart, restart_file, facet110_barrier_adatom,          &
            & facet110_barrier_dimer, facet110_barrier_ledge,layer,    &
            & facet111_barrier_adatom,facet111_barrier_dimer,facet,    &
            & facet111_barrier_ledge,monolayer,ledge_energy_barrier,   &
            & facet110_back_barrier,facet111_back_barrier,             &
            & facet110_self_barrier,facet111_self_barrier,output_file_pdb             
end module input_variables


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


module algorithm_variables
!**********************************************************************
!           this module defines the algorithm variables
! ------------------------------------------------------------------
!              Author: Zhiyong Wang 01/26/98
!**********************************************************************
  
  implicit none

!**********************************************************************
! adatom_number----------------the number of all the isolated adatoms
! adatom_events----------------the event table of isolated adatoms
! adatom_rate------------------the rate sum of all the isolated adatoms
! bond_number------------------the inplane bond number of an atom
! diffusion_rate---------------the rate sum of all the possible 
!                              diffusion events
! dimer_events-----------------the event table of dimers
! dimer_number-----------------the number of all the dimers
! dimer_rate-------------------the rate sum of all the dimers
! height-----------------------surface height (zposition)
! ledge_number-----------------the number of all the ledge atoms
! ledge_events-----------------the event table of ledge atoms
! ledge_rate-------------------the rate sum of all the ledge atoms
! location---------------------the location of a given mobile atom in 
!                              event table 
! pa---------------------------probability for a mobile atom to be in 
!                              adatom group
! pb---------------------------probability for a mobile atom to be in
!                              dimer group
! pl---------------------------probability for a mobile atom to be in
!                              ledge atom group
! steps------------------------the number of simulation steps
! time-------------------------the actual time 
! total_rate-------------------the sum of all the rates of surface 
!                              events
! time_step--------------------the inverse of total_rate
! type-------------------------the type of the atoms(could be 0(immobile)
!                              1(single adatom),2(dimer),3(ledge adatom))
! xpos, ypos-------------------the position of an atom
! xpos1,ypos1------------------the position of dimer atom#1
! xpos2,ypos2---- -------------the position of dimer atom#2
! xminus,yminus,xplus,yplus----the first nearest neighbor position of a
!                              given atom           
!**********************************************************************

  integer :: adatom_number,dimer_number,bond_number,ledge_number,      &
             total_number,steps,facet110_div_1, facet110_div_2,        &
             atoms_monolayer,remain_number, partial_flux     

  real :: time,time_step,total_rate,diffusion_rate,dimer_rate,         &
          adatom_rate,ledge_rate,pa,pd,pl

  integer, allocatable :: height(:,:), type(:,:), location(:,:)

  integer, allocatable :: adatom_events(:,:),dimer_events(:,:,:),      &
                          ledge_events(:,:)

end module algorithm_variables


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


module analysis_variables
!**********************************************************************
!    this module defines the analysis variables, will be added soon
!  -----------------------------------------------------------------
!              Author: Zhiyong Wang 01/26/98
!**********************************************************************
  
implicit none

!**********************************************************************
! island_number----------------the number of islands on surface
! island_size_distribution-----the final distribution of island size
! roughness--------------------surface roughness  
! step_number------------------the number of steps on surface
!**********************************************************************

  integer :: island_number, step_number, island_size_distribution     

  real :: roughness
 
end module analysis_variables


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


program klmc_main
!*********************************************************************
!    kinetic lattice monte carlo code for cu thin film growth on 
!    cu(100) surface
!  ----------------------------------------------------------------
!            Author: Zhiyong Wang 02/05/98
!*********************************************************************
 
use input_variables
use algorithm_variables
use analysis_variables

implicit none
  call date_and_time(big_ben(1),big_ben(2),big_ben(3),date_time)
  call random_seed()
  call data_input
  call basic_rate
  call initialization 
  call klmc_algorithm 
  call data_output
  call data_output_pdb
!  call data_output_xyz(1)

!  call data_analysis

end program klmc_main


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


subroutine data_input
!*********************************************************************
!              the following subroutine gets the input data
!  ---------------------------------------------------------------
!                Author: Zhiyong Wang 01/26/98
!*********************************************************************

use input_variables 
implicit none
open(3,file='sample.in',status='old')
 read(3,nml=input)
 close(3,status='keep')
return

end subroutine data_input


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


subroutine basic_rate
!*********************************************************************
!  this subroutine calculate the basic rate by using input data
! -----------------------------------------------------------------
!             Author: Zhiyong Wang 01/26/98
!*********************************************************************
 use input_variables
 implicit none

 real:: r

 deposition_rate = width * length * deposition_flux
 adatom_hop_rate=adatom_prefactor*exp(-adatom_energy_barrier/(kb*temperature))
 ledge_hop_rate= ledge_prefactor*exp(-ledge_energy_barrier/(kb*temperature))
 dimer_hop_rate = dimer_prefactor*exp(-dimer_energy_barrier/(kb*temperature))
 if (facet=="110") then
    r=exp(-facet110_self_barrier/(kb*temperature))/     &
      (exp(-facet110_back_barrier/(kb*temperature))+    &
       exp(-facet110_self_barrier/(kb*temperature)))    
    facet110_rate_adatom=exp(-facet110_barrier_adatom/(kb*temperature))/      &
                         (exp(-adatom_energy_barrier/(kb*temperature))+       &
                          exp(-facet110_barrier_adatom/(kb*temperature)))*r   
    facet110_rate_dimer=exp(-facet110_barrier_dimer/(kb*temperature))/        &
                        (exp(-dimer_energy_barrier/(kb*temperature))+         &
                         exp(-facet110_barrier_dimer/(kb*temperature)))*r     
    facet110_rate_ledge=exp(-facet110_barrier_ledge/(kb*temperature))/        &
                        (exp(-ledge_energy_barrier/(kb*temperature))+         &
                         exp(-facet110_barrier_ledge/(kb*temperature)))*r     

 else if (facet=="111") then

    r=exp(-facet111_self_barrier/(kb*temperature))/     &
      (exp(-facet111_back_barrier/(kb*temperature))+    &
       exp(-facet111_self_barrier/(kb*temperature)))    
    facet111_rate_adatom=exp(-facet111_barrier_adatom/(kb*temperature))/      &
                         (exp(-adatom_energy_barrier/(kb*temperature))+       &
                          exp(-facet111_barrier_adatom/(kb*temperature)))*r   
    facet111_rate_dimer=exp(-facet111_barrier_dimer/(kb*temperature))/        &
                        (exp(-dimer_energy_barrier/(kb*temperature))+         &
                         exp(-facet111_barrier_dimer/(kb*temperature)))*r     
    facet111_rate_ledge=exp(-facet111_barrier_ledge/(kb*temperature))/        &
                        (exp(-ledge_energy_barrier/(kb*temperature))+         &
                         exp(-facet111_barrier_ledge/(kb*temperature)))*r      
 end if

 return

end subroutine basic_rate


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


subroutine initialization
!*********************************************************************
! the following subroutine sets the initial simulation size and sets 
! initial height, type,location, table elements and time equal to 0
! ------------------------------------------------------------------
!             Author: Zhiyong Wang 02/05/98
!*********************************************************************

use input_variables
use algorithm_variables
implicit none

  integer :: i, j, w, l, x, y, number
  real :: random
  number=0  
  total_number=0
  partial_flux=0
  adatom_number=0
  dimer_number=0
  ledge_number=0 

  if (facet=="110") then
     atoms_monolayer=(width*length)+(width-1)*(length-1)
     if (width<=length) then
        facet110_div_1=width
        facet110_div_2=length
     else
        facet110_div_1=length
        facet110_div_2=width
     end if 
     length=(width+length)-1
     width=length
  end if
          
  allocate(height(0:width+1, 0:length+1))!!!!new for facet
  allocate(type(1:width, 1:length))
  allocate(location(1:width, 1:length))
  allocate(adatom_events(1:2,1:length*width/2))
  allocate(dimer_events(1:2, 1:length*width/6, 1:2))
  allocate(ledge_events(1:2, 1:length*width/6))

  height=0
  type=0
  location=0
  adatom_events = 0
  dimer_events = 0
  ledge_events = 0
  time = 0.0 
  steps=0

  if(in_structure=="none") then 
    initial_atom_number=0
  else if(in_structure=="epitaxy") then