Worked on GaAs model RPA.

fpm.toml

@@ -50,10 +50,10 @@ name="superconda"
 source-dir="app"
 main="superconda.f90"
 
-#[[executable]]
-#name = "screening_comparison"
-#source-dir="test"
-#main = "screening_comparison.f90"
+[[executable]]
+name = "screening_comparison"
+source-dir="src"
+main = "screening_comparison.f90"
 
 [[test]]
 name = "test_misc"
@@ -71,9 +71,9 @@ source-dir="test"
 main = "test_periodictable.f90"
 
 [[test]]
-name = "test_vector_allreps"
+name = "test_vector"
 source-dir="test"
-main = "test_vector_allreps.f90"
+main = "test_vector.f90"
 
 [[test]]
 name = "check_interactions_symmetries"
@@ -84,8 +84,3 @@ main = "check_interactions_symmetries.f90"
 name = "bte_regression"
 source-dir="test"
 main = "bte_regression.f90"
-
-#[[test]]
-#name = "screening_comparison"
-#source-dir="test"
-#main = "screening_comparison.f90"

test/screening_comparison.f90 → src/screening_comparison.f90

@@ -67,12 +67,12 @@ program screening_comparison
 
   !Create wave vector mesh
   numq = 400
-  call linspace(qmags, 0.0_r64, 2.0_r64*kF, numq)
+  call linspace(qmags, 0.0_r64, 3.0_r64*kF, numq)
   call write2file_rank1_real("RPA_test_qmags", qmags)
 
   !Create bosonic energy mesh
   numomega = 400
-  call linspace(Omegas, 0.0_r64, 2.0_r64*eF, numomega)
+  call linspace(Omegas, 0.0_r64, 3.0_r64*eF, numomega)
   call write2file_rank1_real("RPA_test_Omegas", Omegas)
 
   !Calculate analytic Im RPA dielectric function
@@ -212,7 +212,8 @@ subroutine calculate_Imeps(qmags, ens, chempot, m_eff, eF, kF, beta, Imeps)
     end do
     Imeps(1, :) = 0.0_r64
     Imeps = (m_eff/me/bohr2nm/kF/eF/beta)*Imeps
-
+
+!!$    !This gives the same result as the expression above:
 !!$    !Locals
 !!$    integer :: iOmega, iq
 !!$    real(r64) :: E1(size(qmags), size(ens)), E2(size(qmags), size(ens)), Eq(size(qmags))
@@ -222,8 +223,8 @@ subroutine calculate_Imeps(qmags, ens, chempot, m_eff, eF, kF, beta, Imeps)
 !!$    Eq = energy_parabolic(qmags, m_eff)
 !!$
 !!$    do iq = 1, size(qmags)
-!!$       E1(iq, :) = (ens - Eq)**2/4.0/Eq
-!!$       E2(iq, :) = (ens + Eq)**2/4.0/Eq
+!!$       E1(iq, :) = (ens(:) - Eq(iq))**2/4.0/Eq(iq)
+!!$       E2(iq, :) = (ens(:) + Eq(iq))**2/4.0/Eq(iq)
 !!$    end do    
 !!$
 !!$    do iOmega = 1, size(ens)
@@ -267,8 +268,8 @@ subroutine calculate_Reeps(qmags, ens, chempot, m_eff, eF, eplasmon, &
     !Here we need an extra factor of ms/me.
 
     !Magic numbers?
-    ngrid = 10000
-    ymax = 20.0_r64
+    ngrid = 1000
+    ymax = 10.0!20.0_r64
 
     allocate(y(ngrid), I0(ngrid), Reeps(size(qmags), size(ens)))
 
@@ -301,19 +302,15 @@ subroutine calculate_Reeps(qmags, ens, chempot, m_eff, eF, eplasmon, &
           Reeps(iq, iOmega) = aux0*(aux1 - aux2)
        end do
     end do
-
-    !Reeps = epsinf + (0.25_r64/pi/kF/bohr2nm*m_eff/me)*Reeps
 
-    !DBG
-    Reeps = 1.0_r64 + (0.25_r64/pi/kF/bohr2nm*m_eff/me)*Reeps
-    
+    !       valence band contribution + electronic contribution
+    Reeps = epsinf + (0.25_r64/pi/kF/bohr2nm*m_eff/me)*Reeps
+
     !Omega -> 0 limit
     !Reeps(2:size(qmags), 1) = epsinf*&
     !     (1.0_r64 + ks_squared/qmags(2:size(qmags))**2)
 
     !q -> limit
-    !Reeps(1, :) = epsinf*(1.0_r64 - (eplasmon/ens)**2)
-    !DBG
     Reeps(1, :) = 1.0_r64 - (eplasmon/ens)**2
 
 !!$    !Check where the plasmon mode is at the Gamma point

test/test_vector_allreps.f90 → test/test_vector.f90

@@ -1,4 +1,4 @@
-program test_vector_allreps
+program test_vector
 
   use precision, only : r64, i64
   use testify_m, only : testify
@@ -196,4 +196,4 @@ program test_vector_allreps
   call tests_all%report
 
   if(tests_all%get_status() .eqv. .false.) error stop -1
-end program test_vector_allreps
+end program test_vector