Reverted to original definition of qTF; put the deallocation of elect…

…ron eigenvectors after e-ch. imp interaction calculation in the regression test code.
nakib · Sep 22, 2023 · d557069 · d557069
1 parent c75a8bd
commit d557069
Show file tree

Hide file tree

Showing 3 changed files with 23 additions and 25 deletions.
diff --git a/src/bz_sums.f90 b/src/bz_sums.f90
@@ -62,14 +62,15 @@ subroutine calculate_qTF(crys, el)
           end do
        end do
 
-!!$       !Free-electron gas Thomas-Fermi model
-!!$       ! qTF**2 = spindeg*e^2*beta/nptq/vol_pcell/perm0/epsilon0*Sum_{BZ}f0_{k}(1-f0_{k})
-!!$       crys%qTF = sqrt(1.0e9_r64*crys%qTF*el%spindeg*beta*qe**2/product(el%wvmesh)&
-!!$            /crys%volume/perm0/crys%epsilon0) !nm^-1
-
-       !Eq. 22 of PRB 107, 125207 (2023) with the typo in the power of left hand side fixed
+       !Free-electron gas Thomas-Fermi model
+       ! qTF**2 = spindeg*e^2*beta/nptq/vol_pcell/perm0/epsilon0*Sum_{BZ}f0_{k}(1-f0_{k})
        crys%qTF = sqrt(1.0e9_r64*crys%qTF*el%spindeg*beta*qe**2/product(el%wvmesh)&
-            /crys%volume/perm0) !nm^-1
+            /crys%volume/perm0/crys%epsilon0) !nm^-1
+
+       !Compare the above to Eq. 22 of PRB 107, 125207 (2023) with the typo
+       !in the power of left hand side fixed.
+       !crys%qTF = sqrt(1.0e9_r64*crys%qTF*el%spindeg*beta*qe**2/product(el%wvmesh)&
+       !     /crys%volume/perm0) !nm^-1
 
        if(this_image() == 1) then
           write(*, "(A, 1E16.8, A)") ' Thomas-Fermi screening wave vector = ', crys%qTF, ' 1/nm'

diff --git a/src/interactions.f90 b/src/interactions.f90
@@ -129,11 +129,11 @@ real(r64) function gchimp2(el, crys, qcrys, evec_k, evec_kp)
     qcart = matmul(crys%reclattvecs, qcrys)
     qmag = twonorm(qcart)
 
-    !if(qmag == 0.0_r64) then
-    !   eps_3x3 = crys%epsilon0*eye(3_i64)
-    !else
-    eps_3x3 = (crys%epsilon0 + (crys%qTF/qmag)**2)*eye(3_i64)
-    !end if
+    !Note the difference compared to Eq. 23 of PRB 107, 125207 (2023).
+    !This is because in my expression for qTF^2 (module bz_sums > calculate_qTF),
+    !the denominator contains a factor of epsilon0.
+    !As such, (q_TF_effective)^2 in the above paper = epsilon0 x my qTF^2.
+    eps_3x3 = crys%epsilon0*(1.0_r64 + (crys%qTF/qmag)**2)*eye(3_i64)
 
     prefac = 1.0e-3_r64/crys%volume/perm0**2*&
          (el%chimp_conc_n*(qe*el%Zn)**2 + el%chimp_conc_p*(qe*el%Zp)**2)
@@ -146,17 +146,14 @@ real(r64) function gchimp2(el, crys, qcrys, evec_k, evec_kp)
     do ik1 = -3, 3
        do ik2 = -3, 3
           do ik3 = -3, 3
-             G = (  ik1*crys%reclattvecs(:, 1) &
-                  + ik2*crys%reclattvecs(:, 2) &
-                  + ik3*crys%reclattvecs(:, 3)  )
-
-             Gplusq = G + qcart
+             Gplusq = (  ik1*crys%reclattvecs(:, 1) &
+                       + ik2*crys%reclattvecs(:, 2) &
+                       + ik3*crys%reclattvecs(:, 3)  ) + qcart
 
              denom = abs(dot_product(Gplusq, matmul(eps_3x3, Gplusq)))**2
 
              !Only want G /= -q in the sum over G
              if(denom > 0.0_r64) &
-             !if(all(G + qcart /= 0.0_r64)) &
                   Gsum = Gsum + 1.0_r64/denom
           end do
        end do
@@ -2482,7 +2479,7 @@ subroutine calculate_echimp_interaction_ibzk(crys, el, num)
 
                 !Calculate matrix element
                 g2 = gchimp2(el, crys, q_indvec/dble(el%wvmesh), &
-                     el%evecs_irred(ik, m, :), el%evecs(ikp, n, :))  
+                     el%evecs_irred(ik, m, :), el%evecs(ikp, n, :))
 
                 !Evaulate delta function
                 if(num%tetrahedra) then

diff --git a/test/bte_regression.f90 b/test/bte_regression.f90
@@ -273,12 +273,6 @@ program bte_regression
      call t_event%end_timer('IBZ e-ph transition probabilities')
   end if
 
-  if(num%onlyebte .or. num%drag .or. num%phe .or. num%drag &
-       .or. num%plot_along_path) then
-     !Deallocate Wannier quantities
-     call wann%deallocate_wannier(num)
-  end if
-
   if(num%onlyebte .or. num%drag .or. num%phe) then
      !After this point the electron eigenvectors are not needed
      call el%deallocate_eigenvecs
@@ -295,6 +289,12 @@ program bte_regression
      end if
   end if
 
+  if(num%onlyebte .or. num%drag .or. num%phe .or. num%drag &
+       .or. num%plot_along_path) then
+     !Deallocate Wannier quantities
+     call wann%deallocate_wannier(num)
+  end if
+
   if(num%onlyphbte .or. num%drag) then
      if(.not. num%read_V) then
         call t_event%start_timer('IBZ q ph-ph interactions')