API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 20th, 2017
72
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 28.12 KB | None | 0 0
raw download clone embed print report
  1. ! Chuong trinh se tinh n(t), p(t) và u(t)
  2. ! tham so dau vao: e(k), f0   (cac tham so nay se dc nhap vao, hoac viet mot truong chinh tinh toan rieng.
  3. ! Cac gia tri ban dau: n(t=0)=0 ; p(t=0) = 10^-6; u(t=0) =1.; t0=0;
  4.  
  5. Module global      ! Khai bao cac hang so se su dung trong chuong trinh
  6.  implicit none
  7.  save
  8.  real, parameter :: pi= 3.141593
  9.  real, parameter :: hbar = 658.d-3
  10.  complex, parameter :: iu = (0.,1.) , ui = (1.,0.)
  11.  real(8) :: t0
  12.  integer, parameter, public:: estep = 30  ! Luoi' k (cac gia tri cua k)
  13.  complex, dimension(2,2) :: sigma_z, sigma_y  ! khai bao 2 ma tran sigma
  14.  
  15. ! Pho gia tri nang luong theo k.
  16.  real(8) ::    h
  17. real(8),  dimension(estep) ::E
  18.   real(8)::  P_n                                                       ! bien do bom pump_p va pump_n
  19.   complex(8) :: P_p
  20.  real, parameter :: Lc = 10000.0                                           ! don vi nm
  21.  real, parameter :: E0 = 1.616d3 , Ec0 = 0.1 + 1.616d3                                       ! don vi meV
  22.  real, parameter :: hbar_omega = 4.3                                            ! don vi meV
  23.  real(8), parameter :: m_e = 0.067* 5.7d-3 , m_h = 0.45 *5.7d-3
  24.  real(8), parameter :: tau_c = 5.0    , tau_x= 10.0                                                 ! don vi ps
  25.  real(8), parameter :: S = 10.d8                                                 !nm^2
  26. End module global
  27.  
  28. !============================================================================================================
  29.  
  30.  
  31. Module Runge                                                ! Runge-kutta gan dung bac 4.
  32.  use global
  33.  implicit none
  34.  contains                                                   ! Lenh nay dung de bao lai cac subroutine con cua subroutine Runge
  35.    subroutine rk4(t, h, e, p, n,n0, u,A_in, B_in, C_in, D_in, E_in, F_in,G_in, H_in,I_in, J_in, K_in, L_in,&
  36.                               A_out, B_out, C_out, D_out, E_out, F_out, G_out, H_out, I_out, J_out, K_out, L_out, &
  37.                                t_out, p_out, n_out,n0_out, u_out)
  38.                      
  39.  
  40.     implicit none
  41.     real(8), intent(in), dimension(estep) ::e
  42.     complex(8), intent(in) :: p  , n0                              ! Vi p là ham chi phu thuoc t, ko phu thuoc k.
  43.     complex(8), intent(out):: p_out , n0_out
  44.  
  45.     complex(8), intent(in) , dimension(estep) :: n             ! mang n gom estep phan tu, moi phan tu ung voi mot gia tri cua k.
  46.     complex(8), intent(out), dimension(estep) :: n_out
  47.  
  48.     real(8), intent(in) :: t ,  h  
  49.     real(8), intent(out) :: t_out                        ! t la thoi gian, h la buoc nhay cua t
  50.    
  51.     complex(8), intent (in), dimension(2,2,estep):: u
  52.     complex(8), intent (out), dimension(2,2,estep):: u_out  !Ung voi moi gia tri cua k, ham U la mot matran 2x2
  53.  
  54.     complex(8), intent(in),dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in, E_in, F_in,&
  55.                                                                   G_in, H_in, I_in, J_in, K_in, L_in
  56.     complex(8), intent(out),dimension(2,2,2,estep, estep, estep) ::A_out, B_out, C_out, D_out, E_out, F_out, G_out, H_out, &
  57.                                                  I_out, J_out, K_out, L_out
  58.  
  59.  
  60.     ! Khai bao cac he so tam thoi trong subroutine:
  61.     real(8) :: t1, t3
  62.     complex(8) :: kp1, kp2, kp3, kp4, p1, p2, p3, kn01, kn02, kn03, kn04, n01,n02,n03
  63.     complex(8), dimension(estep) :: kn1, kn2, kn3, kn4, n1, n2, n3, x
  64.     complex(8), dimension(2,2,estep):: ku1, ku2, ku3, ku4, u1, u2, u3
  65.     complex(8), dimension(2,2,2, estep, estep, estep) :: A_temp1, A_temp2, A_temp3, A_temp4, B_temp1, B_temp2, B_temp3, B_temp4,  &
  66.                                                          C_temp1, C_temp2, C_temp3, C_temp4, D_temp1, D_temp2, D_temp3, D_temp4  
  67.  
  68.  
  69.     complex(8), dimension(2,2,2, estep, estep, estep) ::  E_temp1, E_temp2, E_temp3, E_temp4, F_temp1, F_temp2, F_temp3, F_temp4, &
  70.                                                    G_temp1, G_temp2, G_temp3, G_temp4, H_temp1, H_temp2, H_temp3, H_temp4, &
  71.                             I_temp1, I_temp2, I_temp3, I_temp4, J_temp1, J_temp2, J_temp3, J_temp4, &
  72.                             K_temp1, K_temp2, K_temp3, K_temp4, L_temp1, L_temp2, L_temp3, L_temp4
  73.    
  74.     complex(8), dimension(2,2,2,estep,estep,estep):: A_1, A_2, A_3, B_1, B_2, B_3, C_1, C_2, C_3, D_1, D_2, D_3, &
  75.                                                      E_1, E_2, E_3, F_1, F_2, F_3, G_1, G_2, G_3, H_1, H_2, H_3, &
  76.                                                      I_1, I_2, I_3, J_1, J_2, J_3, K_1, K_2, K_3, L_1, L_2, L_3                                                      
  77.            x=(0.0,0.0)
  78.  
  79.     call u_deriv(t,e,p,n0,u,ku1)
  80.  
  81.     Call A_deriv(t,p,n,u,A_temp1)
  82.     call B_deriv(t,p,n,u,B_temp1)
  83.     call C_deriv(t,p,n,u,C_temp1)
  84.     call D_deriv(t,p,n,u,D_temp1)
  85.  
  86.     call p_deriv(t,p,n,u, A_in, B_in, C_in, D_in,kp1)
  87.         call n0_deriv(t,p, x, kn01)
  88.  
  89.     call E_deriv(t,p,n,u,E_temp1)
  90.     call F_deriv(t,p,n,u,F_temp1)
  91.     call G_deriv(t,p,n,u,G_temp1)
  92.     call H_deriv(t,p,n,u,H_temp1)
  93.     call I_deriv(t,p,n,u,I_temp1)
  94.     call J_deriv(t,p,n,u,J_temp1)
  95.     call K_deriv(t,p,n,u,K_temp1)
  96.     call L_deriv(t,p,n,u,L_temp1)
  97.  
  98.     call n_deriv(t,p,n,u, E_1, F_1, G_1, H_1, I_1, J_1, K_1, L_1, kn1)
  99.  
  100.       t1  = t+ h/2.0
  101.       u1 = u + h*ku1/2.0
  102.       n1 = n + h*kn1/2.0
  103.       p1 = p + h*kp1/2.0
  104.           n01 =n0 + h*kn01/2.0
  105.       A_1 = A_in + h*A_temp1/2.0
  106.       B_1 = B_in + h*B_temp1/2.0
  107.       C_1 = C_in + h*C_temp1/2.0
  108.       D_1 = D_in + h*D_temp1/2.0
  109.       E_1 = E_in + h*E_temp1/2.0
  110.       F_1 = F_in + h*F_temp1/2.0
  111.       G_1 = G_in + h*G_temp1/2.0
  112.       H_1 = H_in + h*H_temp1/2.0
  113.       I_1 = I_in + h*I_temp1/2.0
  114.       J_1 = J_in + h*J_temp1/2.0
  115.       K_1 = K_in + h*K_temp1/2.0
  116.       L_1 = L_in + h*L_temp1/2.0
  117.  
  118.  
  119.  
  120.     call u_deriv(t1, e, p1, n01, u1, ku2)
  121.  
  122.     Call A_deriv(t1,p1,n1,u1,A_temp2)
  123.     call B_deriv(t1,p1,n1,u1,B_temp2)
  124.     call C_deriv(t1,p1,n1,u1,C_temp2)
  125.     call D_deriv(t1,p1,n1,u1,D_temp2)
  126.  
  127.     call p_deriv(t1,p1,n1,u1, A_1, B_1, C_1, D_1,kp2)
  128.         call n0_deriv(t1,p1, kn1, kn02)
  129.  
  130.     call E_deriv(t1,p1,n1,u1,E_temp2)
  131.     call F_deriv(t1,p1,n1,u1,F_temp2)
  132.     call G_deriv(t1,p1,n1,u1,G_temp2)
  133.     call H_deriv(t1,p1,n1,u1,H_temp2)
  134.     call I_deriv(t1,p1,n1,u1,I_temp2)
  135.     call J_deriv(t1,p1,n1,u1,J_temp2)
  136.     call K_deriv(t1,p1,n1,u1,K_temp2)
  137.     call L_deriv(t1,p1,n1,u1,L_temp2)
  138.  
  139.     call n_deriv(t1,p1,n1,u1, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kn2)
  140.      
  141.       u2 = u + h*ku2/2.0
  142.       n2 = n + h*kn2/2.0
  143.       p2 = p + h*kp2/2.0
  144.           n02= n0+ h*kn02/2.0
  145.  
  146.       A_2 = A_in + h*A_temp2/2.0
  147.       B_2 = B_in + h*B_temp2/2.0
  148.       C_2 = C_in + h*C_temp2/2.0
  149.       D_2 = D_in + h*D_temp2/2.0
  150.       E_2 = E_in + h*E_temp2/2.0
  151.       F_2 = F_in + h*F_temp2/2.0
  152.       G_2 = G_in + h*G_temp2/2.0
  153.       H_2 = H_in + h*H_temp2/2.0
  154.       I_2 = I_in + h*I_temp2/2.0
  155.       J_2 = J_in + h*J_temp2/2.0
  156.       K_2 = K_in + h*K_temp2/2.0
  157.       L_2 = L_in + h*L_temp2/2.0     
  158.      
  159.  
  160.  
  161.     call u_deriv(t1, e, p2, n02, u2, ku3)
  162.  
  163.     Call A_deriv(t1,p2,n2,u2,A_temp3)
  164.     call B_deriv(t1,p2,n2,u2,B_temp3)
  165.     call C_deriv(t1,p2,n2,u2,C_temp3)
  166.     call D_deriv(t1,p2,n2,u2,D_temp3)
  167.  
  168.     call p_deriv(t1, p2, n2, u2, A_2, B_2, C_2, D_2, kp3)
  169.         call n0_deriv(t1, p2, kn2, kn03)
  170.     call E_deriv(t1,p2,n2,u2,E_temp3)
  171.     call F_deriv(t1,p2,n2,u2,F_temp3)
  172.     call G_deriv(t1,p2,n2,u2,G_temp3)
  173.     call H_deriv(t1,p2,n2,u2,H_temp3)
  174.     call I_deriv(t1,p2,n2,u2,I_temp3)
  175.     call J_deriv(t1,p2,n2,u2,J_temp3)
  176.     call K_deriv(t1,p2,n2,u2,K_temp3)
  177.     call L_deriv(t1,p2,n2,u2,L_temp3)
  178.  
  179.     call n_deriv(t1, p2, n2, u2, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kn3)
  180.  
  181.      t3 = t + h
  182.      u3 = u + h*ku3
  183.      n3 = n + h*kn3
  184.      p3 = p + h*kp3
  185.          n03 =n0 + h*kn03
  186.      A_3 = A_in + h*A_temp3
  187.      B_3 = B_in + h*B_temp3
  188.      C_3 = C_in + h*C_temp3
  189.      D_3 = D_in + h*D_temp3
  190.      E_3 = E_in + h*E_temp3
  191.      F_3 = F_in + h*F_temp3
  192.      G_3 = G_in + h*G_temp3
  193.      H_3 = H_in + h*H_temp3
  194.      I_3 = I_in + h*I_temp3
  195.      J_3 = J_in + h*J_temp3
  196.      K_3 = K_in + h*K_temp3
  197.      L_3 = L_in + h*L_temp3
  198.  
  199.      
  200.  
  201.     call u_deriv(t3, e, p3, n03, u3, ku4)
  202.  
  203.     Call A_deriv(t3,p3,n3,u3,A_temp4)
  204.     call B_deriv(t3,p3,n3,u3,B_temp4)
  205.     call C_deriv(t3,p3,n3,u3,C_temp4)
  206.     call D_deriv(t3,p3,n3,u3,D_temp4)
  207.  
  208.     call p_deriv(t3, p3, n3, u3, A_3, B_3, C_3, D_3, kp4)
  209.         call n0_deriv(t3,p3, kn3, kn04)
  210.     call E_deriv(t3,p3,n3,u3,E_temp4)
  211.         call F_deriv(t3,p3,n3,u3,F_temp4)
  212.     call G_deriv(t3,p3,n3,u3,G_temp4)
  213.     call H_deriv(t3,p3,n3,u3,H_temp4)
  214.     call I_deriv(t3,p3,n3,u3,I_temp4)
  215.     call J_deriv(t3,p3,n3,u3,J_temp4)
  216.     call K_deriv(t3,p3,n3,u3,K_temp4)
  217.     call L_deriv(t3,p3,n3,u3,L_temp4)
  218.  
  219.     call n_deriv(t3, p3, n3, u3, E_3, F_3, G_3, H_3, I_3, J_3, K_3, L_3, kn4)
  220.  
  221.      t_out = t3
  222.      u_out = u + h*(ku1 +ku4 + 2.0*(ku2 + ku3))/6.0
  223.      n_out = n + h*(kn1 +kn4 + 2.0*(kn2 + kn3))/6.0
  224.      p_out = p + h*(kp1 +kp4 + 2.0*(kp2 + kp3))/6.0
  225.         n0_out = n0 + h*(kn01 +kn04 +2.0*(kn02+kn03))/6.0
  226.      A_out = A_in + h*(A_temp1 + A_temp4 +2.0*(A_temp2 + A_temp3))/6.0
  227.      B_out = B_in + h*(B_temp1 + B_temp4 +2.0*(B_temp2 + B_temp3))/6.0
  228.      C_out = C_in + h*(C_temp1 + C_temp4 +2.0*(C_temp2 + C_temp3))/6.0
  229.      D_out = D_in + h*(D_temp1 + D_temp4 +2.0*(D_temp2 + D_temp3))/6.0
  230.      E_out = E_in + h*(E_temp1 + E_temp4 +2.0*(E_temp2 + E_temp3))/6.0
  231.      F_out = F_in + h*(F_temp1 + F_temp4 +2.0*(F_temp2 + F_temp3))/6.0
  232.      G_out = G_in + h*(G_temp1 + G_temp4 +2.0*(G_temp2 + G_temp3))/6.0
  233.      H_out = H_in + h*(H_temp1 + H_temp4 +2.0*(H_temp2 + H_temp3))/6.0
  234.      I_out = I_in + h*(I_temp1 + I_temp4 +2.0*(I_temp2 + I_temp3))/6.0
  235.      J_out = J_in + h*(J_temp1 + J_temp4 +2.0*(J_temp2 + J_temp3))/6.0
  236.      K_out = K_in + h*(K_temp1 + K_temp4 +2.0*(K_temp2 + K_temp3))/6.0
  237.      L_out = L_in + h*(L_temp1 + L_temp4 +2.0*(L_temp2 + L_temp3))/6.0
  238.  
  239.     end subroutine rk4
  240.  
  241.     subroutine n0_deriv(t, p, kn, kn0)
  242.      implicit none
  243.      real(8), intent(in) :: t
  244.      complex(8), intent(in), dimension(estep):: kn
  245.      complex(8), intent(in)::p
  246.      complex(8), intent(out):: kn0
  247.      real(8), external ::pump_n
  248.  
  249.      kn0 =   -p*conjg(p)/tau_c
  250.  
  251.     end subroutine n0_deriv
  252.    
  253.    
  254.  
  255.    Subroutine u_deriv(t,e,p,n0,u,ku)
  256.  
  257.     implicit none
  258.  
  259.     real(8), intent(in) :: t
  260.     complex(8), intent(in) ::p ,n0
  261.     complex(8), intent(in), dimension(2,2,estep) :: u
  262.     complex(8), intent(out), dimension(2,2,estep) ::ku
  263.     real(8), intent(in), dimension(estep) ::e
  264.     real(8), external :: f0
  265.     integer :: i,a,b
  266.  
  267.     complex(8), dimension(2,2, estep) :: H                         !H la hamiltonian
  268.     complex(8), dimension(2,2) ::H_temp , u_temp , ku_temp    
  269.                              
  270.       Do i =1 , estep                                             ! Thiet lap H, moi phan tu cua H theo k la mot matran 2x2
  271.       H(:,:,i) = (e(i) - 1614. + f0(0,0)*p*conjg(p))* sigma_z + iu*f0(0,0)*p*p*sigma_y               !tam thoi su dung n(1) thay cho n(0)
  272.      end do
  273.  
  274.          Do i=1 , estep
  275.           Do a=1, 2
  276.            Do b=1, 2
  277.       ku(a,b,i) = -iu* (H(a,1,i)*U(1,b,i) + H(a,2,i)*U(2,b,i))/hbar                 !Day la pt A.11 , Lenh matmul() la nhan ma tran
  278.       end do
  279.          end do
  280.         end do
  281.  
  282.     end subroutine u_deriv
  283.  
  284.  
  285.     Subroutine A_deriv(t, p, n, u, A_temp)
  286.      implicit none
  287.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: A_temp
  288.       real(8), intent(in) :: t
  289.       complex(8), intent(in):: p
  290.       complex(8), dimension(estep), intent(in) :: n
  291.       complex(8), intent(in), dimension(2,2,estep) :: u
  292.       integer :: b,d,m,k,l,q
  293.  
  294.        Do b=1, 2
  295.         Do d=1, 2
  296.          Do m=1, 2
  297.           Do k=1, estep
  298.            Do q=1, estep
  299.             do l=1, estep
  300.             A_temp(b,d,m,k,l,q) = conjg(u(b,1,k))*u(d,1,l)*u(m,1,q)*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p)
  301.             end do
  302.            end do
  303.           end do
  304.          end do
  305.         end do
  306.        end do
  307.      End subroutine A_deriv
  308.  
  309.  
  310.  
  311.      Subroutine B_deriv(t, p, n, u, B_temp)
  312.       implicit none
  313.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: B_temp
  314.       real(8), intent(in) :: t
  315.       complex(8), intent(in):: p
  316.       complex(8), dimension(estep), intent(in) :: n
  317.       complex(8), intent(in), dimension(2,2,estep) :: u
  318.       integer :: b,d,m,k,l,q
  319.        
  320.        Do b=1, 2
  321.         Do d=1, 2
  322.          Do m=1, 2
  323.           Do k=1, estep
  324.            Do q=1, estep
  325.             do l=1, estep
  326.             B_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l) + n(l)*n(q)+n(l)*n(k) - n(k)*n(q))*conjg(p)
  327.             end do
  328.            end do
  329.           end do
  330.          end do
  331.         end do
  332.        end do
  333.      End subroutine B_deriv
  334.  
  335.  
  336.     Subroutine C_deriv(t, p, n, u, C_temp)
  337.      implicit none
  338.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: C_temp
  339.       real(8), intent(in) :: t
  340.       complex(8), intent(in):: p
  341.       complex(8), dimension(estep), intent(in) :: n
  342.       complex(8), intent(in), dimension(2,2,estep) :: u
  343.       integer :: b,d,m,k,l,q
  344.  
  345.        Do b=1, 2
  346.         Do d=1, 2
  347.          Do m=1, 2
  348.           Do k=1, estep
  349.            Do q=1, estep
  350.             do l=1, estep
  351.             C_temp(b,d,m,k,l,q) = conjg(u(b,1,k))*u(d,1,l)*u(m,1,q)*(n(k)+n(k)*n(l)+n(k)*n(q)-n(l)*n(q))*p
  352.              end do
  353.            end do
  354.           end do
  355.          end do
  356.         end do
  357.        end do
  358.      End subroutine C_deriv
  359.  
  360.  
  361.  
  362.     Subroutine D_deriv(t, p, n, u, D_temp)
  363.      implicit none
  364.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: D_temp
  365.       real(8), intent(in) :: t
  366.       complex(8), intent(in):: p
  367.       complex(8), dimension(estep), intent(in) :: n
  368.       complex(8), intent(in), dimension(2,2,estep) :: u
  369.       integer :: b,d,m,k,l,q
  370.  
  371.        Do b=1, 2
  372.         Do d=1, 2
  373.          Do m=1, 2
  374.           Do k=1, estep
  375.            Do q=1, estep
  376.             do l=1, estep
  377.             D_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p
  378.              end do
  379.            end do
  380.           end do
  381.          end do
  382.         end do
  383.        end do
  384.      End subroutine D_deriv
  385.  
  386.  
  387.  
  388.     Subroutine E_deriv(t, p, n, u, E_temp)
  389.      implicit none
  390.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: E_temp
  391.       real(8), intent(in) :: t
  392.       complex(8), intent(in):: p
  393.       complex(8), dimension(estep), intent(in) :: n
  394.       complex(8), intent(in), dimension(2,2,estep) :: u
  395.       integer :: b,d,m,k,l,q
  396.  
  397.        Do b=1, 2
  398.         Do d=1, 2
  399.          Do m=1, 2
  400.           Do k=1, estep
  401.            Do q=1, estep
  402.             do l=1, estep
  403.             E_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*p
  404.             end do
  405.            end do
  406.           end do
  407.          end do
  408.         end do
  409.        end do
  410.      End subroutine E_deriv
  411.  
  412.  
  413. !thuc ra F_deriv chi la conjg cua E_deriv, co the luoc bo cac giai doan tinh cac F,H,J,L nay de lam chuong trinh nhe hon.
  414.     Subroutine F_deriv(t, p, n, u, F_temp)
  415.      implicit none
  416.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: F_temp
  417.       real(8), intent(in) :: t
  418.       complex(8), intent(in):: p
  419.       complex(8), dimension(estep), intent(in) :: n
  420.       complex(8), intent(in), dimension(2,2,estep) :: u
  421.       integer :: b,d,m,k,l,q
  422.  
  423.        Do b=1, 2
  424.         Do d=1, 2
  425.          Do m=1, 2
  426.           Do k=1, estep
  427.            Do q=1, estep
  428.             do l=1, estep
  429.             F_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*p)
  430.              end do
  431.            end do
  432.           end do
  433.          end do
  434.         end do
  435.        end do
  436.      End subroutine F_deriv
  437.  
  438.  
  439.  
  440.  
  441.     Subroutine G_deriv(t, p, n, u, G_temp)
  442.      implicit none
  443.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: G_temp
  444.       real(8), intent(in) :: t
  445.       complex(8), intent(in):: p
  446.       complex(8), dimension(estep), intent(in) :: n
  447.       complex(8), intent(in), dimension(2,2,estep) :: u
  448.       integer :: b,d,m,k,l,q
  449.  
  450.        Do b=1, 2
  451.         Do d=1, 2
  452.          Do m=1, 2
  453.           Do k=1, estep
  454.            Do q=1, estep
  455.              do l=1, estep
  456.             G_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*conjg(p)
  457.             end do
  458.            end do
  459.           end do
  460.          end do
  461.         end do
  462.        end do
  463.      End subroutine G_deriv
  464.  
  465.  
  466.  
  467.  
  468.     Subroutine H_deriv(t, p, n, u, H_temp)
  469.      implicit none
  470.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: H_temp
  471.       real(8), intent(in) :: t
  472.       complex(8), intent(in):: p
  473.       complex(8), dimension(estep), intent(in) :: n
  474.       complex(8), intent(in), dimension(2,2,estep) :: u
  475.       integer :: b,d,m,k,l,q
  476.  
  477.        Do b=1, 2
  478.         Do d=1, 2
  479.          Do m=1, 2
  480.           Do k=1, estep
  481.            Do q=1, estep
  482.             do l=1, estep
  483.             H_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*conjg(p))
  484.             end do
  485.            end do
  486.           end do
  487.          end do
  488.         end do
  489.        end do
  490.      End subroutine H_deriv
  491.  
  492.  
  493.  
  494.     Subroutine I_deriv(t, p, n, u, I_temp)
  495.      implicit none
  496.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: I_temp
  497.       real(8), intent(in) :: t
  498.       complex(8), intent(in):: p
  499.       complex(8), dimension(estep), intent(in) :: n
  500.       complex(8), intent(in), dimension(2,2,estep) :: u
  501.       integer :: b,d,m,k,l,q
  502.  
  503.        Do b=1, 2
  504.         Do d=1, 2
  505.          Do m=1, 2
  506.           Do k=1, estep
  507.            Do q=1, estep
  508.             Do l=1, estep
  509.             I_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p)
  510.                      end do
  511.            end do
  512.           end do
  513.          end do
  514.         end do
  515.        end do
  516.      End subroutine I_deriv
  517.  
  518.  
  519.  
  520.     Subroutine J_deriv(t, p, n, u, J_temp)
  521.      implicit none
  522.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: J_temp
  523.       real(8), intent(in) :: t
  524.       complex(8), intent(in):: p
  525.       complex(8), dimension(estep), intent(in) :: n
  526.       complex(8), intent(in), dimension(2,2,estep) :: u
  527.       integer :: b,d,m,k,l,q
  528.  
  529.        Do b=1, 2
  530.         Do d=1, 2
  531.          Do m=1, 2
  532.           Do k=1, estep
  533.            Do q=1, estep
  534.             do l=1 ,estep
  535.             J_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p))
  536.              end do
  537.            end do
  538.           end do
  539.          end do
  540.         end do
  541.        end do
  542.      End subroutine J_deriv
  543.  
  544.  
  545.  
  546.     Subroutine K_deriv(t, p, n, u, K_temp)
  547.      implicit none
  548.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: K_temp
  549.       real(8), intent(in) :: t
  550.       complex(8), intent(in):: p
  551.       complex(8), dimension(estep), intent(in) :: n
  552.       complex(8), intent(in), dimension(2,2,estep) :: u
  553.       integer :: b,d,m,k,l,q
  554.  
  555.        Do b=1, 2
  556.         Do d=1, 2
  557.          Do m=1, 2
  558.           Do k=1, estep
  559.            Do q=1, estep
  560.             do l= 1, estep
  561.             K_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p
  562.             end do
  563.            end do
  564.           end do
  565.          end do
  566.         end do
  567.        end do
  568.      End subroutine K_deriv
  569.  
  570.  
  571.  
  572.     Subroutine L_deriv(t, p, n, u, L_temp)
  573.      implicit none
  574.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: L_temp
  575.       real(8), intent(in) :: t
  576.       complex(8), intent(in):: p
  577.       complex(8), dimension(estep), intent(in) :: n
  578.       complex(8), intent(in), dimension(2,2,estep) :: u
  579.       integer :: b,d,m,k,l,q
  580.  
  581.        Do b=1, 2
  582.         Do d=1, 2
  583.          Do m=1, 2
  584.           Do k=1, estep
  585.            Do q=1, estep
  586.             do l=1, estep
  587.             L_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p)
  588.              end do
  589.            end do
  590.           end do
  591.          end do
  592.         end do
  593.        end do
  594.      End subroutine L_deriv
  595.  
  596.  
  597.  
  598. ! Trong khai bao duoi day, su dung E_in, F_in...de tranh nham lan voi cac chi so dem i,k,m...
  599.  
  600.      Subroutine n_deriv(t, p, n, u, E_in, F_in, G_in, H_in, I_in, J_in, K_in, L_in, kn)
  601.       implicit none
  602.       complex(8), intent(in), dimension(2,2,2,estep,estep,estep) :: E_in, F_in, G_in, H_in, I_in, J_in, K_in, L_in
  603.       real(8), intent(in) :: t
  604.       complex(8), intent(in):: p
  605.       complex(8), dimension(estep), intent(in) :: n
  606.       complex(8), intent(in), dimension(2,2,estep) :: u
  607.       complex(8), intent(out), dimension(estep):: kn
  608.       integer :: a,b,c,d,i,m,k,l,q
  609.       complex(8), dimension(estep) :: summ
  610.       real(8), external :: f0 , integral    , pump_n                                    ! khai bao ham f0
  611.  
  612.       kn=(0.0,0.0)
  613.       summ = (0.0,0.0)
  614.              
  615.         do k=1, estep
  616.          do q=1, estep
  617.           do l=0, estep-1
  618.            do a=1,2
  619.             do b=1,2
  620.              do c=1,2
  621.               do d=1,2
  622.                do i=1,2
  623.                 do m=1,2
  624.                  
  625.                  if (l > abs(k-q) .and.  l < abs(k+q) .and. (4*k*k*q*q-(k*k+q*q-l*l)**2) > 0) then
  626.                  
  627.                   summ(k) = -(f0(k,q)*f0(k,q)*0.09*0.09/(2.0*pi*pi*pi))*                                                                    &
  628.                             ( (sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  629.                                 conjg(u(1,a,k))*u(1,c,l+1)*u(1,i,q)*E_in(b,d,m,k,l+1,q)*conjg(p)                              &
  630.                                        +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  631.                                 u(1,a,k)*conjg(u(1,c,l+1))*conjg(u(1,i,q))*F_in(b,d,m,k,l+1,q)*p                      &
  632.                                 +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  633.                                 conjg(u(1,a,k))*u(1,c,l+1)*conjg(u(1,i,q))*G_in(b,d,m,k,l+1,q) *p                     &
  634.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  635.                                 u(1,a,k)*conjg(u(1,c,l+1))*u(1,i,q)*H_in(b,d,m,k,l+1,q)*conjg(p)                      &
  636.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  637.                                 conjg(u(2,a,k))*u(2,c,l+1)*u(2,i,q)*I_in(b,d,m,k,l+1,q)*conjg(p)                          &
  638.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  639.                                 u(2,a,k)*conjg(u(2,c,l+1))*conjg(u(2,i,q))*J_in(b,d,m,k,l+1,q) *p                     &
  640.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  641.                                 conjg(u(2,a,k))*u(2,c,l+1)*conjg(u(2,i,q))*K_in(b,d,m,k,l+1,q) *p                     &
  642.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  643.                                 u(2,a,k)*conjg(u(2,c,l+1))*u(2,i,q)*L_in(b,d,m,k,l+1,q)*conjg(p)  )
  644.                  kn(k) = kn(k) + summ(k)
  645.                   end if
  646.                 end do
  647.                end do
  648.               end do
  649.              end do
  650.             end do
  651.            end do
  652.           end do
  653.          end do
  654.         end do
  655.  
  656.          Do k=1, estep
  657.          kn(k) =  kn(k)  - n(k)/tau_x + pump_n(t,k)
  658.          end do
  659.  
  660.      end subroutine n_deriv
  661.  
  662.  
  663.  
  664.  
  665.      subroutine p_deriv(t, p, n, u, A_in, B_in, C_in, D_in, kp)
  666.  
  667.       implicit none
  668.       complex(8), intent(in), dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in
  669.       real(8), intent(in) :: t
  670.       complex(8), intent(in):: p
  671.       complex(8), dimension(estep), intent(in) :: n
  672.       complex(8), intent(in), dimension(2,2,estep) :: u
  673.       complex(8), intent(out):: kp
  674.       integer :: a,b,c,d,i,m,k,l,q
  675.       complex(8):: sump
  676.       real(8), external :: f0, integral
  677.      complex(8), external :: pump_p
  678.  
  679.  
  680.       kp=(0.0,0.0)
  681.       sump = (0.0,0.0)
  682.  
  683.         do k=1, estep
  684.          do q=1, estep
  685.           do l=0, estep-1
  686.            do a=1,2
  687.             do b=1,2
  688.              do c=1,2
  689.               do d=1,2
  690.                do i=1,2
  691.                 do m=1,2
  692.                  
  693.                  if (l > abs(k-q) .and.  l < abs(k+q).and.(4*k*k*q*q-(k*k+q*q-l*l)**2) > 0) then
  694.                         sump =     (f0(k,q)*f0(k,q)*0.09*0.09/(2.0*pi**3))*&    
  695.                        (-(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*3./30.)*integral(k,q,l)* &
  696.                            u(1,a,k)*conjg(u(1,c,l+1))*conjg(u(1,i,q))*A_in(b,d,m,k,l+1,q)                        &
  697.                       - (sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*3./30.)*integral(k,q,l)* &
  698.                            conjg(u(1,a,k))*u(1,c,l+1)*conjg(u(1,i,q))*B_in(b,d,m,k,l+1,q)                        &
  699.                      - (sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*3./30.)*integral(k,q,l)* &
  700.                            u(2,a,k)*conjg(u(2,c,l+1))*conjg(u(2,i,q))*C_in(b,d,m,k,l+1,q)                        &
  701.                        -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*3./30.)*integral(k,q,l)* &
  702.                            conjg(u(2,a,k))*u(2,c,l+1)*conjg(u(2,i,q))*D_in(b,d,m,k,l+1,q) )
  703.                   kp = kp + sump
  704.  
  705.                   end if
  706.                 end do
  707.                end do
  708.               end do
  709.              end do
  710.             end do
  711.            end do
  712.           end do
  713.          end do
  714.         end do
  715.     kp=  conjg(kp)  - p/tau_c + pump_p(t)
  716.  
  717.    
  718.     end subroutine p_deriv
  719.   End module runge
  720.  
  721.  
  722. !===================================================================================================
  723.  
  724. !Bat dau chuong trinh chinh.
  725. !Chuong trinh nay se doc cac tham so dau vao: e(k) tu file. Doc cac gia tri t, h, f0 tu nguoi dung
  726. ! Sau do se xuat ra file các gia tri cua p(t), n(k,t)
  727.  
  728.   Program process
  729.     use global
  730.     use runge
  731.     implicit none
  732.  
  733.     complex(8) :: p , n0                              
  734.     complex(8):: p_out, n0_out
  735.  
  736.     complex(8) , dimension(estep) :: n            
  737.     complex(8), dimension(estep) :: n_out
  738.  
  739.     real(8) :: t  , t_max
  740.     real(8) :: t_out  , n_tot                      
  741.    
  742.     complex(8),  dimension(2,2,estep):: u
  743.     complex(8),  dimension(2,2,estep):: u_out  
  744.  
  745.     complex(8), dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in, E_in, F_in,&
  746.                                                                   G_in, H_in, I_in, J_in, K_in, L_in
  747.     complex(8), dimension(2,2,2,estep, estep, estep) ::A_out, B_out, C_out, D_out, E_out, F_out, G_out, H_out,&
  748.                                                  I_out, J_out, K_out, L_out
  749.         real(8), external :: f0, pump_n
  750.     integer :: i
  751.  
  752.  sigma_z(1,1) = (1.0,0.0)
  753.  sigma_z(1,2) = (0.0,0.0)
  754.  sigma_z(2,1) = (0.0,0.0)
  755.  sigma_z(2,2) =(-1.0,0.0)
  756.  sigma_y(1,1) = (0.0,0.0)
  757.  sigma_y(1,2) =-iu
  758.  sigma_y(2,1) = iu
  759.  sigma_y(2,2) = (0.0,0.0)
  760.  
  761.  
  762. Print *,"Nhap t0"
  763. read *, t
  764.  
  765. Print *, "Nhap h, voi h la buoc cua t"
  766. read * , h
  767.  
  768. Print *, "Nhap t_max"
  769. read *, t_max
  770.  
  771. Print *, "Nhap Pump_p"
  772.  read *, p_p
  773. Print *, "nhap pump_n"
  774.  read *, p_n
  775.  
  776. Print *, "Nhap e(k)"
  777.   open (unit =2, file="list_E.txt")
  778.   do i=1, estep
  779.     read (2,*) e(i)
  780.  
  781.     print *, e(i)
  782.   end do
  783.  
  784. p= (1.d-30,1.d-30)
  785. n=(0.0,0.0)
  786. n0 =(0.0,0.0)
  787. u=(1.0,0.0)
  788. A_in=(0.0,0.0)
  789. B_in=(0.0,0.0)
  790.  C_in=(0.0,0.0)
  791. D_in=(0.0,0.0)
  792. E_in=(0.0,0.0)
  793. F_in=(0.0,0.0)
  794. G_in=(0.0,0.0)
  795. H_in=(0.0,0.0)
  796. I_in=(0.0,0.0)
  797. J_in=(0.0,0.0)
  798. K_in=(0.0,0.0)
  799. L_in=(0.0,0.0)
  800.  
  801.  
  802.  
  803. open (unit=1, file="result1.txt")
  804. open (unit=3, file="result2.txt")
  805. open (unit=4, file="result3.txt")
  806.   write (1,*)    real(n_tot) , real(n_out(1)), real(n_out(2))
  807.   write (3,*)   t_out , real(p_out*conjg(p_out))
  808.   write (4,*)   t_out, p_out , u_out(:,:,1)
  809.  
  810. do while(t <= t_max)
  811.  
  812.  call rk4(t, h, e ,p, n, n0, u,A_in, B_in, C_in, D_in, E_in, F_in, G_in, H_in,I_in, J_in, K_in, L_in,&
  813.                               A_out, B_out, C_out, D_out, E_out, F_out,G_out, H_out, I_out, J_out, K_out, L_out, &
  814.                                t_out, p_out, n_out,n0_out, u_out)
  815.   print *, t_out
  816.   n_tot = real(p_out*conjg(p_out)) + real( Sum( n_out(1:(estep-1))))
  817.   write (1,*)   real(n_tot) , real(n_out(1)), real(n_out(9))
  818.   write (3,*)   t_out , real(p_out*conjg(p_out))
  819.   write (4,*)   t_out,  p_out
  820.   t = t_out
  821.   p = p_out
  822.   n0 =n0_out
  823.   n = n_out
  824.   u = u_out
  825.   A_in = A_out
  826.   B_in = B_out
  827.   C_in = C_out
  828.   D_in = D_out
  829.   E_in = E_out
  830.   F_in = F_out
  831.   G_in = G_out
  832.   H_in = H_out
  833.   I_in = I_out
  834.   J_in = J_out
  835.   K_in = K_out
  836.   L_in = L_out
  837.  
  838. end do
  839.  
  840.  
  841. end program process
  842.  
  843.  
  844.  
  845. !=================================================================================  
  846.  function f0(a,b)
  847.   use global
  848.   implicit none
  849.    real(8) :: f0, k, q, l, Ec_k, Ex_k, Ec_q, Ex_q, Ec_l, Ex_l, Ec_0, Ex_0
  850.    integer :: a,b
  851.  
  852.  k = real(a)*0.3/30.
  853.  
  854.       q = real(b)*0.3/30.
  855.  
  856.       l=k-q
  857.                                                          
  858.  
  859.    Ec_k = sqrt(Ec0**2+ (hbar*k*k/sqrt(13.5))**2)
  860.  
  861.    Ex_k = E0 + hbar*(k**2)/(2.0*(m_e + m_h))
  862.  
  863.    Ec_q = sqrt(Ec0**2+ (hbar*q*q/sqrt(13.5))**2)
  864.  
  865.    Ex_q = E0 + hbar*(q**2)/(2.0*(m_e + m_h))
  866.  
  867.    Ec_l = sqrt(Ec0**2+ (hbar*l*l/sqrt(13.5))**2)
  868.  
  869.    Ex_l = E0 + hbar*(l**2)/(2.0*(m_e + m_h))
  870.  
  871.    Ec_0 = sqrt(Ec0**2)
  872.  
  873.    Ex_0 = E0
  874.  
  875.    f0 = 6.0*10.0*100.0 * 1.0/ ( sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_k-Ex_k)))**2)*  &
  876.                              sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_q-Ex_q)))**2)* &
  877.                              sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_l-Ex_l)))**2)* &
  878.                              sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_0-Ex_0)))**2)*S )
  879.  
  880. End function f0
  881.  
  882. !=======================================================================
  883. ! ham nay rut gon bieu thuc tinh k, p, l
  884.  
  885. Function integral(a,b,d)  
  886.  use global
  887.  implicit none
  888.  real(8) :: integral, k, q, l
  889.  integer :: a,b,d
  890.  
  891.       k = real(a)*0.3/30.
  892.  
  893.       q = real(b)*0.3/30.
  894.  
  895.       l = real(d)*0.3/30.
  896.  
  897.     integral = q*l/(sqrt(4.0*k*k*q*q-(k*k+q*q-l*l)**2))
  898.  
  899. end function integral
  900.  
  901.  
  902.  
  903.  
  904. ! ======================================================================
  905. ! Ham bom cua p(t)
  906.  
  907. Function pump_p(t)
  908.  use global
  909.  implicit none
  910.  
  911.  real, parameter ::  t_p= 3.0    ! gia tri tam thoi de test
  912.  complex(8), parameter::delta_P = (1.d3,0.)
  913.  real(8)::     t
  914.  complex(8) :: pump_p
  915.  
  916.  pump_p = P_p*exp(-2.0*log(2.0)*(t**2)/t_p**2)
  917.  
  918. end function pump_p
  919.  
  920. !======================================================================
  921. ! Ham bom cua n(t)
  922.  
  923. function pump_n(t, a)
  924.  
  925.  use global
  926.  implicit none
  927.  
  928.  real, parameter:: P_0 =5. , t_p =3.0 , delta_E =0.25    ! gia tri tam thoi
  929.  real(8) :: pump_n , t ,k
  930.  integer :: a
  931. if (t  <= 40.) then
  932.  pump_n = P_n *exp(-((E(a)-1.61579d3)**2)/delta_E)
  933. else
  934.  pump_n = 0.
  935.  
  936. end if
  937.  
  938.  end function pump_n
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!