« back to all changes in this revision
Viewing changes to Src/m_ts_contour_eq.f90
- browse files at revision 694
- compare with another revision
- download diff
- download tarball
- view history from revision 694
- Committer: Nick Papior
- Date: 2018-05-04 19:17:31 UTC
- mfrom: (560.1.344 4.1)
- Revision ID: nickpapior@gmail.com-20180504191731-jeidvxb26hxal6d2
Merged r903, pivoting fixes, DM in tbtrans and custom energy-points
- files modified:
- Docs/tbtrans.tex
added
removed
Src/m_ts_contour_eq.f90
442
442
! and not in "random" order
443
443
idx = get_c_io_index(mu%Eq_seg(i))
444
444
445
if ( leqi(Eq_c(idx)%c_io%part,'circle') ) then
445
if ( leqi(Eq_c(idx)%c_io%part,'user') ) then
446
447
call contour_file(Eq_c(idx),mu,lift)
448
449
else if ( leqi(Eq_c(idx)%c_io%part,'circle') ) then
446
450
447
451
call contour_Circle(Eq_c(idx),mu,R,cR)
448
452
513
517
514
518
! the radius can be calculated using two triangles in the circle
515
519
! there is no need to use the cosine-relations
516
R = .5_dp * cR / cos(alpha) ** 2
520
R = 0.5_dp * cR / cos(alpha) ** 2
517
521
518
522
! the real-axis center
519
523
cR = a + R
1127
1131
call TanhSinh_Exact(c%c_io%N,ce,cw,a,b, p=tmp)
1128
1132
1129
1133
end if
1130
1131
1134
1135
case ( CC_USER )
1136
1137
! Read the file information
1138
call contour_file(c,mu,Eta)
1139
1132
1140
case default
1133
1141
write(*,*) 'Method for contour ',trim(c%c_io%name), &
1134
1142
' could not be deciphered: ', c%c_io%method
1135
1143
call die('Could not determine the line-integral')
1136
1144
end select
1137
1138
! I know this is "bad" practice, however, zero is a well-defined quantity.
1139
set_c = sum(abs(c%c(:))) == 0._dp
1140
1145
1141
1146
! get the index in the ID array (same index in w-array)
1142
1147
call ID2idx(c,mu%ID,idx)
1143
1144
do i = 1 , c%c_io%N
1145
if ( set_c ) then
1148
1149
if ( method(c%c_io) == CC_USER ) then
1150
1151
do i = 1 , c%c_io%N
1152
c%w(idx,i) = c%w(idx,i) * nf((real(c%c(i), dp) - mu%mu) / mu%kT)
1153
end do
1154
1155
else
1156
1157
! I know this is "bad" practice, however, zero is a well-defined quantity.
1158
set_c = sum(abs(c%c(:))) == 0._dp
1159
1160
do i = 1 , c%c_io%N
1161
if ( set_c ) then
1146
1162
c%c(i) = dcmplx(ce(i),Eta)
1147
else
1163
else
1148
1164
if ( abs(c%c(i) - dcmplx(ce(i),Eta)) > 1.e-10_dp ) then
1149
call die('contour_line: Error on contour match')
1165
call die('contour_line: Error on contour match')
1150
1166
end if
1151
end if
1152
1153
!ztmp = (c%c(i) - mu%mu) / mu%kT
1154
c%w(idx,i) = cw(i) * nf((ce(i) - mu%mu) / mu%kT)
1155
1156
end do
1157
1167
end if
1168
1169
c%w(idx,i) = cw(i) * nf((ce(i) - mu%mu) / mu%kT)
1170
1171
end do
1172
1173
end if
1174
1158
1175
deallocate(ce,cw)
1159
1176
1160
1177
end subroutine contour_line
1256
1273
1257
1274
case default
1258
1275
1259
! we revert so that we can actually use the line-integral
1276
! we revert so that we can actually use the line-integral
1277
! The tail and line are equivalent in the sense that the
1278
! fermi functions are applied to the weights
1260
1279
c%c_io%part = 'line'
1261
1280
1262
1281
call contour_line(c,mu,Eta)
1438
1457
1439
1458
end subroutine contour_continued_fraction
1440
1459
1460
! This routine will read the contour points from a file
1461
subroutine contour_file(c,mu,Eta)
1462
use m_io, only: io_assign, io_close
1463
use fdf, only: fdf_convfac
1464
1465
type(ts_cw), intent(inout) :: c
1466
type(ts_mu), intent(in) :: mu
1467
! The lifting into the complex plane
1468
real(dp), intent(in) :: Eta
1469
1470
integer :: iu, iostat, ne, idx
1471
logical :: exist
1472
complex(dp) :: E , W
1473
real(dp) :: rE, iE, rW, iW, conv
1474
character(len=512) :: file, line
1475
character(len=16) :: unit
1476
1477
! The contour type contains the file name in:
1478
! c%c_io%cN (weirdly enough)
1479
file = c%c_io%cN
1480
call ID2idx(c,mu%ID,idx)
1481
1482
call io_assign(iu)
1483
1484
! Open the contour file
1485
inquire(file=trim(file), exist=exist)
1486
if ( .not. exist ) then
1487
call die('The file: '//trim(file)//' could not be found to read contour points!')
1488
end if
1489
1490
! Open the file
1491
open(iu, file=trim(file), form='formatted', status='old')
1492
1493
ne = 0
1494
! The default unit is eV.
1495
! On every line an optional unit-specificer may be used to specify the
1496
! subsequent lines units (until another unit is specified)
1497
conv = fdf_convfac('eV', 'Ry')
1498
do
1499
! Now we have the line
1500
read(iu, '(a)', iostat=iostat) line
1501
if ( iostat /= 0 ) exit
1502
if ( len_trim(line) == 0 ) cycle
1503
line = trim(adjustl(line))
1504
if ( line(1:1) == '#' ) cycle
1505
1506
! We have a line with energy and weight
1507
ne = ne + 1
1508
! There are three optional ways of reading this
1509
! 1. ReE ImE, ReW ImW [unit]
1510
read(line, *, iostat=iostat) rE, iE, rW, iW, unit
1511
if ( iostat == 0 ) then
1512
conv = fdf_convfac(unit, 'Ry')
1513
else
1514
read(line, *, iostat=iostat) rE, iE, rW, iW
1515
end if
1516
if ( iostat == 0 ) then
1517
c%c(ne) = dcmplx(rE,iE) * conv
1518
c%w(idx,ne) = dcmplx(rW,iW) * conv
1519
cycle
1520
end if
1521
1522
! 2. ReE ImE, ReW [unit]
1523
iW = 0._dp
1524
read(line, *, iostat=iostat) rE, iE, rW, unit
1525
if ( iostat == 0 ) then
1526
conv = fdf_convfac(unit, 'Ry')
1527
else
1528
read(line, *, iostat=iostat) rE, iE, rW
1529
end if
1530
if ( iostat == 0 ) then
1531
c%c(ne) = dcmplx(rE,iE) * conv
1532
c%w(idx,ne) = dcmplx(rW,iW) * conv
1533
cycle
1534
end if
1535
1536
! 3. ReE , ReW [unit]
1537
iE = Eta
1538
iW = 0._dp
1539
read(line, *, iostat=iostat) rE, rW, unit
1540
if ( iostat == 0 ) then
1541
conv = fdf_convfac(unit, 'Ry')
1542
else
1543
read(line, *, iostat=iostat) rE, rW
1544
end if
1545
if ( iostat == 0 ) then
1546
c%c(ne) = dcmplx(rE * conv,iE)
1547
c%w(idx,ne) = dcmplx(rW,iW) * conv
1548
cycle
1549
end if
1550
1551
call die('Contour file: '//trim(file)//' is not formatted correctly. &
1552
&Please read the documentation!')
1553
1554
end do
1555
1556
call io_close(iu)
1557
1558
if ( c%c_io%N /= ne ) then
1559
call die('Error in reading the contour points from file: '//trim(file))
1560
end if
1561
1562
end subroutine contour_file
1563
1441
1564
function Eq_E(id,step) result(c)
1442
1565
integer, intent(in) :: id
1443
1566
integer, intent(in), optional :: step
Loggerhead is a web-based interface for Breezy
Version: 2.0.1