1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8      subroutine transition_rate( ng, psi, ek, efermi,
9     .                           ekloc, efermiloc, temp,
10     .                           smooth, wmin, wmax,
11     .                           Sr, Aux, Aux2, numh, listhptr, listh,
12     .                           indxuo, no, nuo, nuotot, xij,
13     .                           maxnh, nbands, kpoint, matrix,
14     .                           intraband )
15C *********************************************************************
16C Finds the matrix element for the dipolar transition between to states
17C Written by DSP. August 1999
18C Restyled for f90 version by JDG. June 2004
19C Modified by DSP. April 2010
20C **************************** INPUT **********************************
21C integer ng                  : first dimension of psi, Aux and Aux2
22C real*8  psi(ng,nuotot,nuo)  : Wavefunctions in current k point
23C real*8  ek(nuotot)          : Eigenvalues
24C real*8  temp                : Electronic temperature
25C real*8  efermi              : Fermi level
26C real*8  ekloc(nuotot)       : Eigenvalues (possibly modified by escissor)
27C real*8  efermiloc           : Fermi level (possibly modified by escissor)
28C real*8  wmin                : minimum transition energy required
29C real*8  wmax                : maximum transition energy required
30C real*8  smooth              : artificial width given to transitions
31C real*8  Sr(maxnh)           : Position operator matrix elements (sparse)
32C real*8  Aux(ng,nuotot,nuo)  : Auxiliary space
33C real*8  Aux2(ng,nuotot,nuo)  : Auxiliary space
34C integer numh(nuo)           : Number of nonzero elements of each row
35C                               of hamiltonian matrix
36C integer listhptr(nuo)       : Pointer to start of row in listh
37C integer listh(maxnh)        : Nonzero hamiltonian-matrix element
38C                               column indexes for each matrix row
39C integer indxuo(no)          : Index of equivalent orbital in unit cell
40C                               Unit cell orbitals must be the first in
41C                               orbital lists, i.e. indxuo.le.nuo, with
42C                               nuo the number of orbitals in unit cell
43C integer no                  : Number of basis orbitals in supercell
44C integer nuo                 : Number of orbitals in the cell (locally)
45C integer nuotot              : Number of orbitals in the cell (globally)
46C integer maxnh               : Maximum dimension of listh
47C real*8  xij(3,maxnh)        : Vectors between orbital centers (sparse)
48C                               (not used if only gamma point)
49C integer nocc                : number of occupied states
50C real*8  kpoint(3)           : Current kpoint
51C real*8  dk(3)               : Vector joining the previous and current
52C                               kpoint
53C character matrix*1          : 'R' or 'P' for position or momentum operator
54C *************************** OUTPUT ****************************
55C real*8  Aux(2,nuotot,nuo)   : matrix elements of the dipolar transition
56C real*8  intraband(2,nuo)    : weight of intraband transitions (if a band
57C                               crosses the Fermi energy
58C *************************** UNITS ***********************************
59C Lengths in atomic units (Bohr).
60C k vectors in reciprocal atomic units.
61C Energies in Rydbergs.
62C *********************************************************************
63
64      use precision
65      use parallel,     only : BlockSize, Node, Nodes
66      use parallelsubs, only : LocalToGlobalOrb
67      use m_fermid,     only : stepf
68      use sys,          only : die
69      use alloc,        only : re_alloc, de_alloc
70#ifdef MPI
71      use mpi_siesta
72#endif
73
74      implicit none
75
76C Passed variables
77      integer nuo, nuotot, maxnh, no,
78     .  listh(maxnh), numh(nuo), listhptr(nuo),
79     .  indxuo(no), nbands, ng
80
81      real(dp)
82     .  xij(3,maxnh), Sr(maxnh),
83     .  psi(ng,nuotot,nuo), kpoint(3), Aux(ng,nuotot,nuo),
84     .  Aux2(ng,nuotot,nuo), ek(nuotot), temp, efermi,
85     .  intraband(ng,nuo), smooth, wmin, wmax, ekloc(nuotot),
86     .  efermiloc, ediff
87
88      character
89     .  matrix*1
90
91C Internal variables
92      integer
93     .  ind, iuo, juo, j, ie, iie, iio, jje, je, jo,
94     .  BNodei, Bnodej, BTest, ig
95      real(dp)
96     .  kxij, skxij, ckxij, pipj1, pipj2,
97     .  f1,  f2, intrb(2), intrb2(2)
98
99#ifdef MPI
100      integer                     :: MPIerror
101      real(dp), pointer           :: AuxLocal2(:,:)
102#endif
103      real(dp), pointer           :: AuxLocal(:,:)
104      real(dp), pointer           :: psibandi(:,:)
105      real(dp), pointer           :: psibandj(:,:)
106      real(dp),              save :: ediffmin = 1.0d-3
107      real(dp),              save :: tiny = 1.0d-9
108
109
110C Start timer
111      call timer('transrate',1)
112
113C Check input matrix
114      if(matrix.ne.'P'.and.matrix.ne.'R')
115     $  call die('transrate: matrix only can take values R or P')
116
117C Initialise matrix elements to zero
118      do iuo = 1,nuo
119        do juo = 1,nuotot
120         do ig=1,ng
121          Aux(ig,juo,iuo) = 0.0_dp
122          Aux2(ig,juo,iuo) = 0.0_dp
123         enddo
124        enddo
125        do ig=1,ng
126           intraband(ig,iuo)=0.0_dp
127        enddo
128      enddo
129
130C Compute matrix elements
131      do iuo = 1,nuo
132        do j = 1,numh(iuo)
133          ind = listhptr(iuo) + j
134          jo = listh(ind)
135          juo = indxuo(jo)
136          if(ng.gt.1) then
137           kxij = kpoint(1) * xij(1,ind) +
138     .            kpoint(2) * xij(2,ind) +
139     .            kpoint(3) * xij(3,ind)
140           ckxij = dcos(kxij)
141           skxij = dsin(kxij)
142           Aux2(1,juo,iuo) = Aux2(1,juo,iuo) + Sr(ind)*ckxij
143           Aux2(2,juo,iuo) = Aux2(2,juo,iuo) + Sr(ind)*skxij
144          else
145           Aux2(1,juo,iuo) = Aux2(1,juo,iuo) + Sr(ind)
146          endif
147        enddo
148      enddo
149
150C Allocate workspace array
151#ifdef MPI
152      nullify( AuxLocal2 )
153      call re_alloc( AuxLocal2, 1, ng, 1, nuotot,
154     &               name    = 'AuxLocal2',
155     &               routine = 'transition_rate' )
156#endif
157      nullify( AuxLocal )
158      call re_alloc( AuxLocal, 1, ng, 1, nuotot,
159     &               name    = 'AuxLocal',
160     &               routine = 'transition_rate' )
161      nullify( psibandi )
162      call re_alloc( psibandi, 1, ng, 1, nuotot,
163     &               name    = 'psibandi',
164     &               routine = 'transition_rate' )
165      nullify( psibandj )
166      call re_alloc( psibandj, 1, ng, 1, nuotot,
167     &               name    = 'psibandj',
168     &               routine = 'transition_rate' )
169
170      BNodei = 0
171      iie = 0
172      do ie = 1,nbands
173        f1 = 2.0d0*stepf((ekloc(ie)-efermiloc)/temp)
174        if (Node.eq.BNodei) then
175          iie = iie + 1
176          do j = 1,nuotot
177            do ig=1,ng
178              psibandi(ig,j) = psi(ig,j,iie)
179            enddo
180          enddo
181        endif
182#ifdef MPI
183        call MPI_Bcast(psibandi(1,1),ng*nuotot,MPI_double_precision,
184     .    BNodei,MPI_Comm_World,MPIerror)
185#endif
186
187C Compute the intra-band contribution to include a Drude-like term
188C for metals
189       if(matrix.eq.'P') then
190             intrb(:)=0.0_dp
191             ediff=ek(ie)-efermi
192             if(dabs(ediff).lt.2.0d0*smooth) then
193              do iuo = 1,nuo
194                call LocalToGlobalOrb(iuo,Node,Nodes,iio)
195                do juo = 1,nuotot
196                 if(ng.eq.2) then
197                   pipj1 = psibandi(1,iio)*psibandi(1,juo) +
198     .                     psibandi(2,iio)*psibandi(2,juo)
199                   pipj2 = psibandi(1,iio)*psibandi(2,juo) -
200     .                     psibandi(2,iio)*psibandi(1,juo)
201
202C The factor 1/2 is necessary because the in subroutine phirphi_opt
203C we have added a factor of 2, to be canceled by the energy denominator
204C in Ry (rather than in Ha). For the calculation of the
205C Drude term we do not have the energy denominator
206                   intrb(1) = intrb(1)
207     .                + 0.5_dp * (pipj1*Aux2(1,juo,iuo)
208     .                - pipj2*Aux2(2,juo,iuo))
209
210                   intrb(2) = intrb(2)
211     .                + 0.5_dp * (pipj1*Aux2(2,juo,iuo)
212     .                + pipj2*Aux2(1,juo,iuo))
213                  else
214                   pipj1 = psibandi(1,iio)*psibandi(1,juo)
215
216                   intrb(1) = intrb(1)
217     .                + 0.5_dp* (pipj1*Aux2(1,juo,iuo))
218                   intrb(2) = 0.0_dp
219                  endif
220                 enddo
221               enddo
222              endif
223#ifdef MPI
224        call MPI_Reduce(intrb,intrb2,2,
225     .    MPI_double_precision,MPI_sum,BNodei,MPI_Comm_World,MPIerror)
226        if (Node.eq.BNodei) then
227          intraband(1:ng,iie) = intrb2(1:ng)
228        endif
229#else
230        if (Node.eq.BNodei) then
231          intraband(1:ng,iie) = intrb(1:ng)
232        endif
233#endif
234        endif
235
236        AuxLocal(1:ng,1:nuotot) = 0.0d0
237
238        BNodej = 0
239        jje = 0
240        do je = 1,nbands
241          if (Node.eq.BNodej) then
242            jje = jje + 1
243          endif
244          if (dabs(ek(ie)-ek(je)).gt.ediffmin) then
245            f2 = 2.0d0*stepf((ekloc(je)-efermiloc)/temp)
246
247            if (f1*(2.0d0-f2).gt.tiny) then
248
249             ediff = ekloc(je) - ekloc(ie)
250             if ((ediff.ge.max(wmin-2.0d0*smooth,0.0d0)).and.
251     .               (ediff.le.wmax+2.0d0*smooth)) then
252
253              if (Node.eq.BNodej) then
254                do j = 1,nuotot
255                  do ig=1,ng
256                    psibandj(ig,j) = psi(ig,j,jje)
257                  enddo
258                enddo
259              endif
260#ifdef MPI
261              call MPI_Bcast(psibandj(1,1),ng*nuotot,
262     $        MPI_double_precision,BNodej,MPI_Comm_World,MPIerror)
263#endif
264
265              do iuo = 1,nuo
266                call LocalToGlobalOrb(iuo,Node,Nodes,iio)
267                do juo = 1,nuotot
268                 if(ng.eq.2) then
269                   pipj1 = psibandi(1,iio)*psibandj(1,juo) +
270     .                     psibandi(2,iio)*psibandj(2,juo)
271                   pipj2 = psibandi(1,iio)*psibandj(2,juo) -
272     .                     psibandi(2,iio)*psibandj(1,juo)
273
274                   AuxLocal(1,je) = AuxLocal(1,je)
275     .                + pipj1*Aux2(1,juo,iuo)
276     .                - pipj2*Aux2(2,juo,iuo)
277
278                   AuxLocal(2,je) = AuxLocal(2,je)
279     .                + pipj1*Aux2(2,juo,iuo)
280     .                + pipj2*Aux2(1,juo,iuo)
281                  else
282                   pipj1 = psibandi(1,iio)*psibandj(1,juo)
283
284                   AuxLocal(1,je) = AuxLocal(1,je)
285     .                + pipj1*Aux2(1,juo,iuo)
286                  endif
287
288                enddo
289              enddo
290              if (matrix.eq.'P') then
291                do ig=1,ng
292                 AuxLocal(ig,je) = AuxLocal(ig,je)/(ek(je)-ek(ie))
293                enddo
294              endif
295             endif
296            endif
297          endif
298          BTest = je/BlockSize
299          if (BTest*BlockSize.eq.je) then
300            BNodej = BNodej + 1
301            if (BNodej .gt. Nodes-1) BNodej = 0
302          endif
303        enddo
304#ifdef MPI
305        call MPI_Reduce(AuxLocal(1,1),AuxLocal2(1,1),ng*nuotot,
306     .    MPI_double_precision,MPI_sum,BNodei,MPI_Comm_World,MPIerror)
307        if (Node.eq.BNodei) then
308          Aux(1:ng,1:nuotot,iie) = AuxLocal2(1:ng,1:nuotot)
309        endif
310#else
311        if (Node.eq.BNodei) then
312          Aux(1:ng,1:nuotot,iie) = AuxLocal(1:ng,1:nuotot)
313        endif
314#endif
315        BTest = ie/BlockSize
316        if (BTest*BlockSize.eq.ie) then
317          BNodei = BNodei + 1
318          if (BNodei .gt. Nodes-1) BNodei = 0
319        endif
320      enddo
321
322C Free workspace array
323      call de_alloc( psibandj, 'psibandj', 'transition_rate')
324      call de_alloc( psibandi, 'psibandi', 'transition_rate')
325      call de_alloc( AuxLocal, 'AuxLocal', 'transition_rate')
326#ifdef MPI
327      call de_alloc( AuxLocal2, 'AuxLocal2', 'transition_rate')
328#endif
329
330C Stop timer
331      call timer('transrate',2)
332
333      end
334