← Back to branch summary

~ubuntu-branches/debian/stretch/abinit/stretch

« back to all changes in this revision

Viewing changes to src/17ddb/eliashberg_1d.F90

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2007-09-14 13:05:00 UTC
  • Revision ID: james.westby@ubuntu.com-20070914130500-1kzh2mrgo6ir4b6i
Tags: upstream-5.3.4.dfsg
Import upstream version 5.3.4.dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/17ddb/eliashberg_1d.F90
 
1
!{\src2tex{textfont=tt}}
 
2
!!****f* ABINIT/eliashberg_1d
 
3
!!
 
4
!! NAME
 
5
!! eliashberg_1d
 
6
!!
 
7
!! FUNCTION
 
8
!!  Solve the Eliashberg equations in the isotropic case
 
9
!!   First the linearized case, which allows the estimation of Tc
 
10
!!   then the full case which gives the gap as a function of temperature.
 
11
!!
 
12
!! COPYRIGHT
 
13
!! Copyright (C) 2004-2007 ABINIT group (MVer)
 
14
!! This file is distributed under the terms of the
 
15
!! GNU General Public Licence, see ~abinit/COPYING
 
16
!! or http://www.gnu.org/copyleft/gpl.txt .
 
17
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
 
18
!!
 
19
!! INPUTS
 
20
!!  a2f_1d = alpha^2F function averaged over the FS (only energy dependence)
 
21
!!  elph_ds = datastructure with phonon matrix elements
 
22
!!  gkk2 = gkk2 matrix elements on full FS grid for each phonon mode
 
23
!!  n0 = DOS at the Fermi level calculated from the FSkpt integration weights
 
24
!!  natom = number of atoms
 
25
!!  nsym = number of symops
 
26
!!  phon_ds = datastructure with interatomic force constants
 
27
!!
 
28
!! OUTPUT
 
29
!!
 
30
!! SIDE EFFECTS
 
31
!!
 
32
!! PARENTS
 
33
!!      elphon
 
34
!!
 
35
!! CHILDREN
 
36
!!      eli_diag_m_1d,eli_lambda_1d,eli_z_1d
 
37
!!
 
38
!! NOTES
 
39
!!  na2f = number of frequency points for alpha^2F function
 
40
!!
 
41
!! SOURCE
 
42
 
 
43
#if defined HAVE_CONFIG_H
 
44
#include "config.h"
 
45
#endif
 
46
 
 
47
subroutine eliashberg_1d(a2f_1d,elph_ds,mustar,n0,natom,nsym,phon_ds)
 
48
 
 
49
 use defs_basis
 
50
 use defs_datatypes
 
51
 use defs_elphon
 
52
 
 
53
!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
 
54
#ifdef HAVE_FORTRAN_INTERFACES
 
55
 use interfaces_17ddb, except_this_one => eliashberg_1d
 
56
#endif
 
57
!End of the abilint section
 
58
 
 
59
 implicit none
 
60
 
 
61
!Arguments ------------------------------------
 
62
  ! needed for phonon interpolation
 
63
!scalars
 
64
 integer,intent(in) :: natom,nsym
 
65
 real(dp),intent(in) :: mustar,n0
 
66
 type(elph_type),intent(in) :: elph_ds
 
67
 type(phon_type),intent(in) :: phon_ds
 
68
!arrays
 
69
 real(dp),intent(in) :: a2f_1d(elph_ds%na2f)
 
70
 
 
71
!Local variables-------------------------------
 
72
  ! for diagonalization of gammma matrix
 
73
  ! output variables for gtdyn9+phfrq3
 
74
!scalars
 
75
 integer :: eivec=1,i1,i2,iatom,ib1,ib2,ibranch,idir,ier,ii,iiter,imatsu,indx
 
76
 integer :: iost,ip,ipoint,iqpt,iqptfull,iseg,jbranch,jmatsu,k1,kbranch,kdir
 
77
 integer :: maxiter,maxmappl,mu,nmatsu,nu,qtor,unit_del,unit_lam,unit_z
 
78
 real(dp) :: dist,gaussfactor,gaussprefactor,gaussval,maxeigval,omega_cutoff
 
79
 real(dp) :: phnow,qphnrm=one,res,tc,total_weight,weight
 
80
 character(len=fnlen) :: fname
 
81
!arrays
 
82
 real(dp) :: displ(2,3*natom,3*natom)
 
83
 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
 
84
 real(dp) :: eigval(elph_ds%nbranch),eigvec(2*elph_ds%nbranch*elph_ds%nbranch)
 
85
 real(dp) :: gam_now(2,elph_ds%nbranch,elph_ds%nbranch)
 
86
 real(dp) :: gam_now2(2,elph_ds%nbranch,elph_ds%nbranch)
 
87
 real(dp) :: gammafact(elph_ds%nbranch),pheigval(3*natom)
 
88
 real(dp) :: pheigvec(2*3*natom*3*natom),phfrq_tmp(3*natom),qpt(3),redkpt(3)
 
89
 real(dp),allocatable :: delta_1d(:),lambda_1d(:),matrx(:,:),mm_1d(:,:),z_1d(:)
 
90
 real(dp),allocatable :: zhpev1(:,:),zhpev2(:)
 
91
 
 
92
! *********************************************************************
 
93
 
 
94
#ifdef VMS
 
95
!DEC$ ATTRIBUTES ALIAS:'ZHPEV' :: zhpev
 
96
#endif
 
97
 
 
98
 write (*,*) 'eliashberg_1d : enter '
 
99
 
 
100
! maximum number of iterations to find T_c
 
101
 maxiter=30
 
102
 
 
103
! Fix nmatsu. Should add test in iiter loop to check if
 
104
!  omega_cutoff is respected
 
105
 nmatsu = 50
 
106
! write (*,*) ' eliashberg_1d : nmatsu = ', nmatsu
 
107
 
 
108
 allocate (lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu))
 
109
 allocate (delta_1d(-nmatsu:nmatsu),mm_1d(-nmatsu:nmatsu,-nmatsu:nmatsu))
 
110
 
 
111
 unit_lam=108
 
112
 fname=trim(elph_ds%elph_base_name) // "_LAM"
 
113
 open (UNIT=unit_lam,FILE=fname,STATUS='REPLACE')
 
114
 unit_z=109
 
115
 fname=trim(elph_ds%elph_base_name) // "_Z"
 
116
 open (UNIT=unit_z,FILE=fname,STATUS='REPLACE')
 
117
 unit_del=110
 
118
 fname=trim(elph_ds%elph_base_name) // "_DEL"
 
119
 open (UNIT=unit_del,FILE=fname,STATUS='REPLACE')
 
120
 
 
121
!
 
122
!  1) use linearized Eliashberg equation to find Tc
 
123
! $ \sum_j \mathbf{M}_{ij} \Delta_j = \zeta \cdot \Delta_i $ $i,j = 1 .. n_{\mathrm{Matsubara}}$
 
124
! $\zeta = 1$ gives T$_c$ $\beta = \frac{1}{\mathrm{T}}$ $\omega_i = (2 i + 1) \pi \mathrm{T}$
 
125
! $ \mathbf{M}_{ij} = \frac{\pi}{\beta} \frac{\lambda (\omega_i - \omega_j)}{Z (\omega_i)}$
 
126
! $ Z (\omega_i) = 1 + \frac{\pi}{\beta \omega_i} \sum_j \lambda(\omega_i - \omega_j) \mathrm{sgn}(\omega_j)$
 
127
!
 
128
 
 
129
! initial guess for T$_c$ in Hartree (1Ha =3.067e5 K)
 
130
 tc = 0.0001
 
131
!
 
132
!   big iterative loop
 
133
!
 
134
 do iiter=1,maxiter
 
135
 
 
136
  omega_cutoff = (two*nmatsu+one) * pi * tc
 
137
!  write (*,*) 'eliashberg_1d : omega_cutoff = ', omega_cutoff
 
138
!  write (*,*) 'eliashberg_1d : tc = ', tc
 
139
 
 
140
!
 
141
! calculate array of lambda values
 
142
!
 
143
  call eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc)
 
144
  write (unit_lam,'(a)') '#'
 
145
  write (unit_lam,'(a)') '# ABINIT package : lambda file'
 
146
  write (unit_lam,'(a)') '#'
 
147
  write (unit_lam,'(a,I10,a)') '# lambda_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
 
148
  write (unit_lam,'(a,E16.6,a,E16.6)') '#  from ', -omega_cutoff, ' to ', omega_cutoff
 
149
  write (unit_lam,'(a)') '#  lambda_1d is the frequency dependent coupling constant '
 
150
  write (unit_lam,'(a)') '#  in the Eliashberg equations '
 
151
  write (unit_lam,'(a)') '#'
 
152
  do imatsu=-nmatsu,nmatsu
 
153
   write (unit_lam,*) imatsu,lambda_1d(imatsu)
 
154
  end do
 
155
  write (unit_lam,*)
 
156
 
 
157
!
 
158
! calculate array of z values
 
159
!
 
160
  call eli_z_1d (lambda_1d,nmatsu,tc,z_1d)
 
161
  write (unit_z,'(a)') '#'
 
162
  write (unit_z,'(a)') '# ABINIT package : Z file'
 
163
  write (unit_z,'(a)') '#'
 
164
  write (unit_z,'(a,I10,a)') '# z_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
 
165
  write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff
 
166
  write (unit_z,'(a)') '# z_1d is the renormalization factor in the Eliashberg equations'
 
167
  write (unit_z,'(a)') '#'
 
168
  do imatsu=-nmatsu,nmatsu
 
169
   write (unit_z,*) imatsu,z_1d(imatsu)
 
170
  end do
 
171
 
 
172
!!
 
173
!! apply M matrix until a maximal eigenvalue is found.
 
174
!!
 
175
!  call eli_m_iter_1d (delta_1d,lambda_1d,maxeigval,nmatsu,tc,z_1d)
 
176
 
 
177
!
 
178
!   diagonalize M brute forcefully
 
179
!
 
180
  call eli_diag_m_1d(delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d)
 
181
 
 
182
  write (unit_del,'(a)') '#'
 
183
  write (unit_del,'(a)') '# eliashberg_1d : delta_1d = '
 
184
  write (unit_del,'(a)') '#'
 
185
  write (unit_del,'(a,i6,a)') '# delta_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
 
186
  write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff
 
187
  write (unit_z,'(a)') '# delta_1d is the gap function in the Eliashberg equations'
 
188
  write (unit_z,'(a)') '#'
 
189
  do imatsu=-nmatsu,nmatsu
 
190
   write (unit_del,*) imatsu,delta_1d(imatsu)
 
191
  end do
 
192
  write (unit_del,*)
 
193
 
 
194
!  write (111,*) tc, maxeigval
 
195
!
 
196
!  if eigenvalue is < 1 increase T
 
197
!  else if eigenvalue is > 1 decrease T
 
198
!  if eigenvalue ~= 1 stop
 
199
!
 
200
  if (abs(maxeigval-one) < tol8) then
 
201
    write (*,*) 'Eliashberg Tc found = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)'
 
202
    exit
 
203
  else if (maxeigval > 0.001_dp) then
 
204
    tc = tc * maxeigval
 
205
  else
 
206
    write (*,*) 'maxeigval is very small'
 
207
    tc = tc * 1000.0_dp
 
208
  end if
 
209
 
 
210
 
 
211
 end do
 
212
!end iiter do
 
213
 
 
214
 if (abs(maxeigval-one) > tol8) then
 
215
   write (*,*) 'eliashberg_1d : Tc not converged. ', maxeigval, ' /= 1'
 
216
   write (*,*) 'Eliashberg Tc nonetheless = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)'
 
217
 end if
 
218
 
 
219
 deallocate (lambda_1d,z_1d,delta_1d,mm_1d)
 
220
 
 
221
 close (UNIT=unit_z)
 
222
 close (UNIT=unit_lam)
 
223
 close (UNIT=unit_del)
 
224
 
 
225
 write (*,*) ' eliashberg_1d : end '
 
226
 
 
227
 
 
228
end subroutine eliashberg_1d
 
229
!!***