File: epsinv.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (109 lines) | stat: -rw-r--r-- 3,157 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

! Copyright (C) 2010 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine epsinv
use modmain
use modmpi
use modomp
implicit none
! local variables
integer iq,ik,ig,iw
integer n,info,nthd
! allocatable arrays
integer, allocatable :: ipiv(:)
integer(8), allocatable :: lock(:)
real(8), allocatable :: vgqc(:,:),gqc(:),gclgq(:),jlgqr(:,:,:)
complex(8), allocatable :: ylmgq(:,:),sfacgq(:,:)
complex(8), allocatable :: epsi(:,:,:),work(:)
! allocate local arrays
allocate(vgqc(3,ngrf),gqc(ngrf),gclgq(ngrf))
allocate(jlgqr(njcmax,nspecies,ngrf))
allocate(ylmgq(lmmaxo,ngrf),sfacgq(ngrf,natmtot))
allocate(epsi(ngrf,ngrf,nwrf))
! initialise the OpenMP locks
allocate(lock(nwrf))
do iw=1,nwrf
  call omp_init_lock(lock(iw))
end do
if (mp_mpi) write(*,*)
! synchronise MPI processes
call mpi_barrier(mpicom,ierror)
! loop over q-points
do iq=1,nqpt
  if (mp_mpi) write(*,'("Info(epsinv): ",I6," of ",I6," q-points")') iq,nqpt
! generate the G+q-vectors and related quantities
  call gengqrf(vqc(:,iq),vgqc,gqc,jlgqr,ylmgq,sfacgq)
! generate the regularised Coulomb Green's function in G+q-space
  call gengclgq(.true.,iq,ngrf,gqc,gclgq)
! use the symmetric form
  gclgq(:)=sqrt(gclgq(:))
! zero the response function (stored in epsi)
  epsi(:,:,:)=0.d0
  call omp_hold(nkptnr/np_mpi,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ik=1,nkptnr
! distribute among MPI processes
    if (mod(ik-1,np_mpi).ne.lp_mpi) cycle
! compute v^1/2 chi0 v^1/2
    call genvchi0(.false.,ik,lock,0.d0,vql(:,iq),gclgq,jlgqr,ylmgq,sfacgq, &
     ngrf,epsi)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
! add epsi from each process and redistribute
  if (np_mpi.gt.1) then
    n=nwrf*ngrf*ngrf
    call mpi_allreduce(mpi_in_place,epsi,n,mpi_double_complex,mpi_sum,mpicom, &
     ierror)
  end if
! negate and add delta(G,G')
  epsi(:,:,:)=-epsi(:,:,:)
  do ig=1,ngrf
    epsi(ig,ig,:)=epsi(ig,ig,:)+1.d0
  end do
!-------------------------------------!
!     invert epsilon over G-space     !
!-------------------------------------!
  call omp_hold(nwrf,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(ipiv,work,info) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do iw=1,nwrf
    allocate(ipiv(ngrf),work(ngrf))
    call zgetrf(ngrf,ngrf,epsi(:,:,iw),ngrf,ipiv,info)
    if (info.eq.0) call zgetri(ngrf,epsi(:,:,iw),ngrf,ipiv,work,ngrf,info)
    if (info.ne.0) then
      write(*,*)
      write(*,'("Error(epsinv): unable to invert epsilon")')
      write(*,'(" for q-point ",I6)') iq
      write(*,'(" and frequency ",I6)') iw
      write(*,*)
      stop
    end if
    deallocate(ipiv,work)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
! write inverse RPA epsilon to EPSINV.OUT
  if (mp_mpi) call putepsinv(iq,epsi)
! end loop over q-points
end do
! destroy the OpenMP locks
do iw=1,nwrf
  call omp_destroy_lock(lock(iw))
end do
deallocate(lock)
deallocate(vgqc,gqc,gclgq,jlgqr)
deallocate(ylmgq,sfacgq,epsi)
! synchronise MPI processes
call mpi_barrier(mpicom,ierror)
return
end subroutine