File: nonlinopt.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 (254 lines) | stat: -rw-r--r-- 8,487 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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254

! Copyright (C) 2002-2010 S. Sharma, J. K. Dewhurst 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.

!BOP
! !ROUTINE: nonlinopt
! !INTERFACE:
subroutine nonlinopt
! !USES:
use modmain
use modomp
! !DESCRIPTION:
!   Calculates susceptibility tensor for non-linear optical second-harmonic
!   generation (SHG). The terms (ztm) are numbered according to Eqs. (49)-(51)
!   of the article {\it Physica Scripta} {\bf T109}, 128 (2004). Other good
!   references are {\it Phys. Rev. B} {\bf 48}, 11705 (1993) and
!   {\it Phys. Rev. B} {\bf 53}, 10751 (1996).
!
! !REVISION HISTORY:
!   Rewrote earlier version, June 2010 (Sharma)
!EOP
!BOC
implicit none
! local variables
integer ik,jk,ist,jst,kst
integer iw,a,b,c,l,nthd
! smallest eigenvalue difference allowed in denominator
real(8), parameter :: etol=1.d-4
real(8) eji,eki,ekj,t1
complex(8) pii(3),dji(3),vji(3),vik(3),vkj(3)
complex(8) eta,ztm(3,3),z1
character(256) fname
! allocatable arrays
real(8), allocatable :: w(:)
complex(8), allocatable :: pmat(:,:,:)
complex(8), allocatable :: chiw(:,:),chi2w(:,:)
! initialise universal variables
call init0
call init1
! read Fermi energy from file
call readfermi
! get the eigenvalues and occupancies from file
do ik=1,nkpt
  call getevalsv(filext,ik,vkl(:,ik),evalsv(:,ik))
  call getoccsv(filext,ik,vkl(:,ik),occsv(:,ik))
end do
! generate energy grid (starting from zero)
allocate(w(nwplot))
t1=wplot(2)/dble(nwplot)
do iw=1,nwplot
  w(iw)=t1*dble(iw-1)
end do
! allocate response function arrays
allocate(chiw(nwplot,3))
allocate(chi2w(nwplot,2))
! i divided by the complex relaxation time
eta=cmplx(0.d0,swidth,8)
! begin loop over components
do l=1,noptcomp
  a=optcomp(1,l)
  b=optcomp(2,l)
  c=optcomp(3,l)
  chiw(:,:)=0.d0
  chi2w(:,:)=0.d0
! parallel loop over non-reduced k-points
  call omp_hold(nkptnr,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(pmat,jk,ist,jst,kst) &
!$OMP PRIVATE(eji,eki,ekj,t1,z1,pii) &
!$OMP PRIVATE(dji,vji,vik,vkj,ztm,iw) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ik=1,nkptnr
    allocate(pmat(nstsv,nstsv,3))
!$OMP CRITICAL(nonlinopt_1)
    write(*,'("Info(nonlinopt): ",I6," of ",I6," k-points")') ik,nkptnr
!$OMP END CRITICAL(nonlinopt_1)
! equivalent reduced k-point
    jk=ivkik(ivk(1,ik),ivk(2,ik),ivk(3,ik))
! read momentum matrix elements from file
    call getpmat(.false.,vkl(:,ik),pmat)
! scissor correct the matrix elements
    do ist=1,nstsv
      if (evalsv(ist,jk).lt.efermi) then
        do jst=1,nstsv
          if (evalsv(jst,jk).gt.efermi) then
            eji=evalsv(jst,jk)-evalsv(ist,jk)
            t1=(eji+scissor)/eji
            pmat(ist,jst,:)=t1*pmat(ist,jst,:)
          end if
        end do
      end if
    end do
    z1=(wkptnr*occmax/omega)*zi
! loop over valence states
    do ist=1,nstsv
      if (evalsv(ist,jk).lt.efermi) then
        pii(:)=pmat(ist,ist,:)
! loop over conduction states
        do jst=1,nstsv
          if (evalsv(jst,jk).gt.efermi) then
            eji=evalsv(jst,jk)-evalsv(ist,jk)+scissor
            dji(:)=pmat(jst,jst,:)-pii(:)
            vji(:)=pmat(jst,ist,:)
! loop over intermediate states
            do kst=1,nstsv
              if ((kst.ne.ist).and.(kst.ne.jst)) then
                eki=evalsv(kst,jk)-evalsv(ist,jk)
                ekj=evalsv(kst,jk)-evalsv(jst,jk)
                if (evalsv(kst,jk).gt.efermi) then
                  eki=eki+scissor
                else
                  ekj=ekj-scissor
                end if
                vik(:)=pmat(ist,kst,:)
                vkj(:)=pmat(kst,jst,:)
! interband terms
                t1=-eji*eki*(-ekj)*(eki+ekj)
                if (abs(t1).gt.etol) then
                  t1=1.d0/t1
                else
                  t1=0.d0
                end if
                ztm(1,1)=z1*conjg(vji(a))*(conjg(vkj(b))*conjg(vik(c)) &
                 +conjg(vik(b))*conjg(vkj(c)))*t1
                t1=eji*(-eki)*ekj*(-eki-eji)
                if (abs(t1).gt.etol) then
                  t1=1.d0/t1
                else
                  t1=0.d0
                end if
                ztm(1,2)=0.5d0*z1*vkj(c)*(vji(a)*vik(b)+vik(a)*vji(b))*t1
                t1=eji*(-eki)*ekj*(ekj-eji)
                if (abs(t1).gt.etol) then
                  t1=1.d0/t1
                else
                  t1=0.d0
                end if
                ztm(1,3)=0.5d0*z1*vik(b)*(vkj(c)*vji(a)+vji(c)*vkj(a))*t1
! intraband terms
                t1=(-eki)*ekj*eji**3
                if (abs(t1).gt.etol) then
                  t1=1.d0/t1
                else
                  t1=0.d0
                end if
                ztm(2,1)=0.5d0*z1*(eki*vik(b)*(vkj(c)*vji(a)+vji(c)*vkj(a)) &
                 +ekj*vkj(c)*(vji(a)*vik(b)+vik(a)*vji(b)))*t1
                t1=((-eji)*eki*(-ekj)*eji**2)/(-ekj-eki)
                if (abs(t1).gt.etol) then
                  t1=1.d0/t1
                else
                  t1=0.d0
                end if
                ztm(2,3)=z1*conjg(vji(a))*(conjg(vkj(b))*conjg(vik(c)) &
                 +conjg(vik(b))*conjg(vkj(c)))*t1
! modulation terms
                t1=ekj*(-eki)*eji**3
                if (abs(t1).gt.etol) then
                  t1=1.d0/t1
                else
                  t1=0.d0
                end if
                ztm(3,1)=0.25d0*z1*(-eki)*vkj(a)*(vji(b)*vik(c)+vik(b)*vji(c))*t1
                ztm(3,2)=0.25d0*z1*ekj*vik(a)*(vkj(b)*vji(c)+vji(b)*vkj(c))*t1
!$OMP CRITICAL(nonlinopt_2)
                do iw=1,nwplot
! 2w interband
                  chi2w(iw,1)=chi2w(iw,1)+ztm(1,1)/(eji-2.d0*w(iw)+eta)
! 2w intraband
                  chi2w(iw,2)=chi2w(iw,2)+ztm(2,3)/(eji-2.d0*w(iw)+eta)
! w interband
                  chiw(iw,1)=chiw(iw,1)-(ztm(1,2)-ztm(1,3))/(eji-w(iw)+eta)
! w intraband
                  chiw(iw,2)=chiw(iw,2)+ztm(2,1)/(eji-w(iw)+eta)
! w modulation
                  chiw(iw,3)=chiw(iw,3)+(ztm(3,1)-ztm(3,2))/(eji-w(iw)+eta)
                end do
!$OMP END CRITICAL(nonlinopt_2)
              end if
! end loop over kst
            end do
            ztm(2,2)=4.d0*z1*conjg(vji(a))*(dji(b)*vji(c)+vji(b)*dji(c))/eji**4
            ztm(3,3)=0.25d0*z1*vji(a)*(vji(b)*dji(c)+dji(b)*vji(c))/eji**4
!$OMP CRITICAL(nonlinopt_2)
            do iw=1,nwplot
! 2w intraband
              chi2w(iw,2)=chi2w(iw,2)+ztm(2,2)/(eji-2.d0*w(iw)+eta)
! w modulation
              chiw(iw,3)=chiw(iw,3)+ztm(3,3)/(eji-w(iw)+eta)
            end do
!$OMP END CRITICAL(nonlinopt_2)
! end loop over jst
          end if
        end do
! end loop over ist
      end if
    end do
    deallocate(pmat)
! end loop over k-points
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
! write to files
  write(fname,'("CHI_INTER2w_",3I1,".OUT")') a,b,c
  open(51,file=trim(fname),form='FORMATTED')
  write(fname,'("CHI_INTRA2w_",3I1,".OUT")') a,b,c
  open(52,file=trim(fname),form='FORMATTED')
  write(fname,'("CHI_INTERw_",3I1,".OUT")') a,b,c
  open(53,file=trim(fname),form='FORMATTED')
  write(fname,'("CHI_INTRAw_",3I1,".OUT")') a,b,c
  open(54,file=trim(fname),form='FORMATTED')
  write(fname,'("CHI_",3I1,".OUT")') a,b,c
  open(55,file=trim(fname),form='FORMATTED')
  do iw=1,nwplot
    write(51,'(2G18.10)') w(iw),dble(chi2w(iw,1))
    write(52,'(2G18.10)') w(iw),dble(chi2w(iw,2))
    write(53,'(2G18.10)') w(iw),dble(chiw(iw,1))
    write(54,'(2G18.10)') w(iw),dble(chiw(iw,2))
    t1=dble(chi2w(iw,1)+chi2w(iw,2)+chiw(iw,1)+chiw(iw,2)+chiw(iw,3))
    write(55,'(2G18.10)') w(iw),t1
  end do
  write(51,'("     ")')
  write(52,'("     ")')
  write(53,'("     ")')
  write(54,'("     ")')
  write(55,'("     ")')
  do iw=1,nwplot
    write(51,'(2G18.10)') w(iw),aimag(chi2w(iw,1))
    write(52,'(2G18.10)') w(iw),aimag(chi2w(iw,2))
    write(53,'(2G18.10)') w(iw),aimag(chiw(iw,1))
    write(54,'(2G18.10)') w(iw),aimag(chiw(iw,2))
    t1=aimag(chi2w(iw,1)+chi2w(iw,2)+chiw(iw,1)+chiw(iw,2)+chiw(iw,3))
    write(55,'(2G18.10)') w(iw),t1
  end do
  close(51); close(52); close(53); close(54); close(55)
! end loop over components
end do
write(*,*)
write(*,'("Info(nonlinopt):")')
write(*,'(" susceptibility tensor written to CHI_abc.OUT")')
write(*,'(" interband contributions written to CHI_INTERx_abc.OUT")')
write(*,'(" intraband contributions written to CHI_INTRAx_abc.OUT")')
write(*,'(" for components")')
do l=1,noptcomp
  write(*,'("  a = ",I1,", b = ",I1,", c = ",I1)') optcomp(1:3,l)
end do
deallocate(w,chiw,chi2w)
return
end subroutine
!EOC