Implementation of RPA dielectric screening #170
Conversation
| ! NOTE: Now the wannier deallocation is happening after solve_BTE | ||
| if(num%need_Wannier) then | ||
| !Deallocate Wannier quantities | ||
| call wann%deallocate_wannier(num) | ||
| end if | ||
|
|
There was a problem hiding this comment.
Just leaving a note here to mention that this is a temporary solution. The code basically ends after the BTE solution, so there is no need to deallocate anything manually at that stage. Earlier deallocation of those quantities was needed to free up memory. But I'll handle this later.
|
|
||
| prefac = 1.0e9_r64*qe/(perm0*crys%epsilon0) ! ev.nm | ||
|
|
||
| !Transfer wave vector in Cartesian coordinates |
| qcart = matmul(crys%reclattvecs, qcrys) | ||
|
|
||
| !Use a safe range for the G vector sums | ||
| !Assuming G = G' (neglecting local-field effects) |
There was a problem hiding this comment.
We are not assuming G = G', we are ignoring G /= G' terms.
| do ik1 = -3, 3 | ||
| do ik2 = -3, 3 | ||
| do ik3 = -3, 3 | ||
| G_plusq = ( ik1*crys%reclattvecs(:, 1) & |
There was a problem hiding this comment.
What is the meaning of the underscore? We are calculating G + q, so Gplusq is a perfectly good name.
| diel_qw = 1 - prefac*X0_qw/twonorm(G_plusq)**2 | ||
| W_qw = W_qw + 1.0_r64/diel_qw/twonorm(G_plusq)**2 |
There was a problem hiding this comment.
Here you are calculating the same 2-norm twice. Better define a variable and evaluate the 2-norm once.
| diel_qw = 1 - prefac*X0_qw/twonorm(G_plusq)**2 | ||
| W_qw = W_qw + 1.0_r64/diel_qw/twonorm(G_plusq)**2 |
There was a problem hiding this comment.
There should be comments above each of these lines explaining what exactly is being done here physically. Also, please check the style guide regarding using spaces to separate conceptually distinct statements.
| !! linear interpolation from 1d array evaluated on a fine, continuous mesh. | ||
| !! to a sample mesh |
There was a problem hiding this comment.
Here document what the arguments and the output are.
| @@ -1694,4 +1694,35 @@ subroutine Hilbert_transform(fx, Hfx) | |||
| Hfx(k + 1) = -(fx(k + 2) - fx(k) + term2 + term3)/pi | |||
| end do | |||
| end subroutine Hilbert_transform | |||
|
|
|||
| pure function interpolator_1d(samp, cont, res_cont) result(res_samp) | |||
There was a problem hiding this comment.
What does res_cont mean? By looking at res_samp from its context, it seems like the prefix res is supposed to mean result. But why would an argument of a function contain a result?
| @@ -1694,4 +1694,35 @@ subroutine Hilbert_transform(fx, Hfx) | |||
| Hfx(k + 1) = -(fx(k + 2) - fx(k) + term2 + term3)/pi | |||
| end do | |||
| end subroutine Hilbert_transform | |||
|
|
|||
| pure function interpolator_1d(samp, cont, res_cont) result(res_samp) | |||
There was a problem hiding this comment.
These names are all quite cryptic to me.
| do isamp = 1, nsamp | ||
| w = samp(isamp) | ||
| ileft = minloc(abs(cont - w), dim = 1) | ||
| if(ileft == ncont) then | ||
| res_samp(isamp) = res_cont(ileft) | ||
| else | ||
| iright = ileft + 1 | ||
| res_samp(isamp) = ((cont(iright) - w)*res_cont(ileft) + & | ||
| (w - cont(ileft))*res_cont(iright))/dcont | ||
| end if | ||
| end do |
There was a problem hiding this comment.
Does this 1-d interpolator work generally? Can you add some documentation to briefly explain how this works?
| call test_array(itest)%assert([9, 19, 29]*1.0_r64, interpolator_1d([2, 4, 6]*1.0_r64, & | ||
| [1, 3, 5, 7]*1.0_r64, [4, 14, 24, 34]*1.0_r64)) |
There was a problem hiding this comment.
It would help to explain a bit what this test it about. For example, what linear function are you testing?
Added a new function to compute RPA screened Coulomb interaction, used to compute the electron electron scattering rates. Also there is a new function added in misc.f90, which is a 1d interpolator (along with its test).
A small note: In some places I have added "?" in the comments which should be looked at in the future.