API tools faq
paste
Advertisement
Guest User

Untitled

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