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
9! This code segment has been fully created by:
10! Nick Papior, 2014
11
12module m_tbt_options
13
14  use precision, only : dp
15
16  use units, only: Ang
17
18  use m_ts_tdir, only: ts_tidx
19
20  use m_ts_electype
21  use m_ts_chem_pot
22
23  use dictionary
24
25  implicit none
26
27  ! Common flags for parameters
28  public
29  save
30
31  ! The standard name_prefix
32#ifdef TBT_PHONON
33  character(len=*), parameter :: name_prefix = 'PHT'
34#else
35  character(len=*), parameter :: name_prefix = 'TBT'
36#endif
37
38  ! The temperature
39  real(dp) :: kT
40
41  ! Electrodes and different chemical potentials
42  integer :: N_Elec = 0
43  type(Elec), allocatable, target :: Elecs(:)
44  integer :: N_mu = 0
45  type(ts_mu), allocatable, target :: mus(:)
46
47  ! Whether we should stop right after having created
48  ! the Green's function files
49  logical :: stop_after_GS = .false.
50
51  ! Dictionary to contain the data saving methods
52  ! Each key corresponds to some data calculation
53  ! algorithm.
54  ! To check whether data should be calculated do:
55  ! if ( 'DOS-Gf' .in. save_DATA ) then
56  !   calculate DOS of Gf
57  ! end fi
58  type(dictionary_t) :: save_DATA
59
60  ! Number of eigenchannels to calculate
61  integer :: N_eigen = 0
62
63  ! If the energy-contour is not perfectly divisable by the number of nodes then adjust
64  integer :: BTD_method = 0 ! Optimization method for determining the best tri-diagonal matrix split
65  ! 0  == We optimize for speed
66  ! 1  == We optimize for memory
67
68  ! A quantity describing the accuracy of the coordinates of the
69  ! electrodes.
70  ! * Should only be edited by experienced users *
71  real(dp) :: Elecs_xa_EPS = 1.e-4_dp
72
73  ! Every 5% of the calculation progress it will print an estimation
74  real :: percent_tracker = 5.
75
76#ifdef NCDF_4
77  ! Save file names for data files
78  character(len=250) :: cdf_fname = ' '
79  character(len=250) :: cdf_fname_sigma = ' '
80  character(len=250) :: cdf_fname_proj = ' '
81#endif
82
83
84  ! List of private formats for printing information
85  character(len=*), parameter, private :: f1 ='(''tbt: '',a,t53,''='',tr4,l1)'
86  character(len=*), parameter, private :: f10='(''tbt: '',a,t53,''='',tr4,a)'
87  character(len=*), parameter, private :: f11='(''tbt: '',a)'
88  character(len=*), parameter, private :: f12='(''tbt: '',a,t53,''='',tr2,i0)'
89  character(len=*), parameter, private :: f5 ='(''tbt: '',a,t53,''='',i5,a)'
90  character(len=*), parameter, private :: f20='(''tbt: '',a,t53,''='',i0,'' -- '',i0)'
91  character(len=*), parameter, private :: f6 ='(''tbt: '',a,t53,''='',f10.4,tr1,a)'
92  character(len=*), parameter, private :: f7 ='(''tbt: '',a,t53,''='',f12.6,tr1,a)'
93  character(len=*), parameter, private :: f8 ='(''tbt: '',a,t53,''='',f10.4)'
94  character(len=*), parameter, private :: f9 ='(''tbt: '',a,t53,''='',tr1,e9.3)'
95  character(len=*), parameter, private :: f15='(''tbt: '',a,t53,''='',2(tr1,i0,'' x''),'' '',i0)'
96
97
98contains
99
100  subroutine read_tbt_generic(na_u, lasto)
101
102    use m_region, only: rgn_delete
103    use fdf, only: fdf_defined
104    use m_ts_method, only: ts_init_regions
105    ! This array is never used used, so we delete it
106    use m_ts_method, only: r_pvt
107
108    use m_tbt_diag, only: init_diag
109
110    ! The number of atoms
111    integer, intent(in) :: na_u
112    ! A summated list of last orbitals on atoms.
113    integer, intent(in) :: lasto(0:na_u)
114
115    ! Initialize the buffer regions
116    if ( fdf_defined('TBT.Atoms.Buffer') ) then
117       call ts_init_regions('TBT',na_u,lasto)
118    else
119       call ts_init_regions('TS',na_u,lasto)
120    end if
121
122    call rgn_delete(r_pvt)
123
124    ! Initialize the diagonalization method.
125    call init_diag( )
126
127  end subroutine read_tbt_generic
128
129
130  ! > Reads the chemical potentials as well as the applied
131  ! Bias.
132  ! The bias is an intricate part of the chemical potential why it
133  ! is read in here.
134  subroutine read_tbt_chem_pot( )
135
136    use fdf, only : fdf_get
137    use units, only: eV, Kelvin
138
139    use m_ts_chem_pot, only : fdf_nmu, fdffake_mu, fdf_mu, name
140
141    use m_tbt_hs, only: Volt, IsVolt
142
143    implicit none
144
145    ! *******************
146    ! * LOCAL variables *
147    ! *******************
148    logical :: err
149    integer :: i
150
151    ! Read in the temperature
152    kT = fdf_get('ElectronicTemperature',1.9e-3_dp,'Ry')
153    kT = fdf_get('TS.ElectronicTemperature',kT,'Ry')
154    kT = fdf_get('TBT.ElectronicTemperature',kT,'Ry')
155
156    ! Read in the chemical potentials
157    N_mu = fdf_nmu('TBT',kT,mus)
158    if ( N_mu < 1 ) then
159       N_mu = fdf_nmu('TS',kT,mus)
160    end if
161    err = .true.
162    if ( N_mu < 1 ) then
163       err = .false.
164       if ( IsVolt ) then
165          ! There is a bias: default
166          ! to two chemical potentials with the
167          ! applied bias
168          N_mu = fdffake_mu(mus,kT,Volt)
169
170       else
171          ! There is no bias. Simply create
172          ! one chemical potential.
173          ! This will make the electrodes
174          ! default to the one chemical potential.
175          N_mu = 1
176          allocate(mus(1))
177          mus(1)%kT = kT
178          mus(1)%ckT = ' '
179          mus(1)%name = 'Fermi-level'
180          mus(1)%mu = 0._dp
181          mus(1)%cmu = '0. eV'
182          mus(1)%ID = 1
183          ! These are not used, but we populate
184          ! them anyway
185          mus(1)%N_poles = 1
186          allocate(mus(1)%Eq_seg(1))
187          mus(1)%Eq_seg(1) = '*NONE'
188
189       end if
190
191    end if
192    do i = 1 , N_mu
193       ! Default things that could be of importance
194       if ( fdf_mu('TBT',mus(i),kT,Volt) ) then
195          ! success
196       else if ( fdf_mu('TS',mus(i),kT,Volt) ) then
197          ! success
198       else if ( err ) then
199          ! only error out if it couldn't be found and forced
200          ! created
201          call die('Could not find chemical potential: ' &
202               //trim(name(mus(i))))
203       end if
204    end do
205
206#ifdef TBT_PHONON
207    ! Phonon transport cannot define different chemical potentials
208    ! Furthermore, they should be zero
209    do i = 1 , N_mu
210       if ( abs(mus(i)%mu) > 1.e-10 * eV ) then
211          call die('Phonon transport does not define chemical &
212               &potentials. I.e. you cannot lift the frequency spectra.')
213       end if
214    end do
215#endif
216
217  end subroutine read_tbt_chem_pot
218
219
220  ! Reads all information regarding the electrodes, nothing more.
221  subroutine read_tbt_elec( cell, na_u, xa, lasto)
222
223    use fdf, only : fdf_get, fdf_obsolete, fdf_deprecated, leqi
224    use parallel, only : IONode
225    use intrinsic_missing, only : IDX_SPC_PROJ, EYE
226    use intrinsic_missing, only : VEC_PROJ_SCA, VNORM
227
228    use m_os, only : file_exist
229
230    use files, only: slabel
231    use units, only: eV
232
233    use m_tbt_hs, only: spin_idx
234
235    use m_ts_chem_pot, only : copy, chem_pot_add_Elec
236
237    use m_ts_electype, only : fdf_nElec, fdf_Elec
238    use m_ts_electype, only : Name, TotUsedOrbs, TotUsedAtoms
239    use m_ts_electype, only : init_Elec_sim
240
241    use m_ts_method, only : ts_init_electrodes, a_isBuffer
242
243    implicit none
244
245    ! *******************
246    ! * INPUT variables *
247    ! *******************
248    real(dp), intent(in) :: cell(3,3)
249    integer,  intent(in) :: na_u, lasto(0:na_u)
250    real(dp), intent(in) :: xa(3,na_u)
251
252    ! *******************
253    ! * LOCAL variables *
254    ! *******************
255    integer :: i, j
256    real(dp) :: rtmp
257    logical :: err, bool
258
259    if ( N_mu == 0 ) call die('read_tbt_elecs: error in programming')
260
261    ! To determine the same coordinate nature of the electrodes
262    Elecs_xa_EPS = fdf_get('TS.Elecs.Coord.Eps',0.001_dp*Ang, 'Bohr')
263    Elecs_xa_EPS = fdf_get('TBT.Elecs.Coord.Eps',Elecs_xa_EPS,'Bohr')
264
265    ! detect how many electrodes we have
266    N_Elec = fdf_nElec('TBT', Elecs)
267    if ( N_Elec < 1 ) N_Elec = fdf_nElec('TS', Elecs)
268    if ( N_Elec < 1 ) then
269       ! We initialize to 2 electrodes (Left/Right)
270       N_Elec = 2
271       allocate(Elecs(N_Elec))
272       Elecs(1)%name = 'Left'
273       Elecs(1)%ID = 1
274       Elecs(2)%name = 'Right'
275       Elecs(2)%ID = 2
276       ! if they do-not exist, the user will be told
277       if ( IONode ) then
278          write(*,'(/,''tbt: *** '',a)') 'No electrode names were found, &
279               &default Left/Right are expected'
280       end if
281    end if
282
283    ! Setup default parameters for the electrodes
284    ! first electrode is the "left"
285    ! last electrode is the "right"
286    ! the remaining electrodes have their chemical potential at 0
287    ! Currently the transport direction for all electrodes is the default
288    ! We should probably warn if +2 electrodes are used and t_dir is the
289    ! same for all electrodes... Then the user needs to know what (s)he is doing...
290    Elecs(:)%Bulk = fdf_get('TS.Elecs.Bulk',.true.) ! default everything to bulk electrodes
291    Elecs(:)%Bulk = fdf_get('TBT.Elecs.Bulk',Elecs(1)%Bulk)
292
293    rtmp = fdf_get('TS.Elecs.Eta',0.001_dp*eV,'Ry')
294    rtmp = fdf_get('TBT.Elecs.Eta',rtmp,'Ry')
295#ifdef TBT_PHONON
296    ! eta value needs to be squared as it is phonon spectrum
297    if ( rtmp > 0._dp ) rtmp = rtmp ** 2
298#endif
299    Elecs(:)%Eta = rtmp
300
301    rtmp = fdf_get('TS.Elecs.Accuracy',1.e-13_dp*eV,'Ry')
302    rtmp = fdf_get('TBT.Elecs.Accuracy',rtmp,'Ry')
303    Elecs(:)%accu = rtmp
304
305    ! whether all calculations should be performed
306    ! "out-of-core" i.e. whether the GF files should be created or not
307    ! In tbtrans this is now defaulted to in-core
308    Elecs(:)%out_of_core = fdf_get('TBT.Elecs.Out-of-core',.false.)
309
310    ! Whether we should try and re-use the surface Green function
311    ! files
312    Elecs(:)%ReUseGF = fdf_get('TS.Elecs.GF.ReUse',.true.)
313    Elecs(:)%ReUseGF = fdf_get('TBT.Elecs.GF.ReUse',Elecs(1)%ReUseGF)
314
315    ! Will stop after creating the GF files.
316    stop_after_GS = fdf_get('TBT.Elecs.GF.Only',.false.)
317
318    do i = 1 , N_Elec
319
320       ! If we only have 2 electrodes we take them
321       ! as though the atomic indices are the first and last
322       ! respectively.
323       if ( N_Elec == 2 ) then
324          if ( i == 1 ) then
325             err = fdf_Elec('TBT',slabel,Elecs(i),N_mu,mus,idx_a= 1, &
326                  name_prefix = name_prefix)
327             if ( .not. err ) &
328                  err = fdf_Elec('TS',slabel,Elecs(i),N_mu,mus,idx_a= 1, &
329                  name_prefix = name_prefix)
330          else
331             err = fdf_Elec('TBT',slabel,Elecs(i),N_mu,mus,idx_a=-1, &
332                  name_prefix = name_prefix)
333             if ( .not. err ) &
334                  err = fdf_Elec('TS',slabel,Elecs(i),N_mu,mus,idx_a=-1, &
335                  name_prefix = name_prefix)
336          end if
337       else
338          ! Default things that could be of importance
339          err = fdf_Elec('TBT',slabel,Elecs(i),N_mu,mus)
340          if ( .not. err ) &
341               err = fdf_Elec('TS',slabel,Elecs(i),N_mu,mus, &
342               name_prefix = name_prefix)
343       end if
344       if ( .not. err ) then
345          call die('Could not find electrode: '//trim(name(Elecs(i))))
346       end if
347       if ( Elecs(i)%idx_a < 0 ) &
348            Elecs(i)%idx_a = na_u + Elecs(i)%idx_a + 1
349       if ( Elecs(i)%idx_a < 1 .or. &
350            na_u < Elecs(i)%idx_a ) then
351          print *,Elecs(i)%idx_a,na_u
352          call die("Electrode position does not exist")
353       end if
354       if ( N_Elec == 2 ) then
355          ! Correct for buffer atoms, first electrode steps "up"
356          ! second electrode steps "down"
357          if ( i == 1 ) then
358             j = Elecs(i)%idx_a
359             do while ( a_isBuffer(j) )
360                j = j + 1
361             end do
362             Elecs(i)%idx_a = j
363          else
364             j = Elecs(i)%idx_a + TotUsedAtoms(Elecs(i)) - 1
365             do while ( a_isBuffer(j) )
366                j = j - 1
367             end do
368             Elecs(i)%idx_a = j - TotUsedAtoms(Elecs(i)) + 1
369          end if
370       end if
371       ! set the placement in orbitals
372       Elecs(i)%idx_o = lasto(Elecs(i)%idx_a-1)+1
373
374       ! Initialize electrode parameters
375       call init_Elec_sim(Elecs(i),cell,na_u,xa)
376
377    end do
378
379    ! Initialize the electrode regions
380    call ts_init_electrodes(na_u,lasto,N_Elec,Elecs)
381
382    ! If many electrodes, no transport direction can be specified
383    ! Hence we use this as an error-check (also for N_Elec == 1)
384    if ( any(Elecs(:)%t_dir > 3) ) then
385      ts_tidx = - N_Elec
386    else
387
388      if ( N_Elec /= 2 ) then
389        ! Signals no specific unit-cell direction of transport
390        ts_tidx = - N_Elec
391      else
392
393        ! Retrieve the indices of the unit-cell directions
394        ! according to the electrode transport directions.
395        ! We have already calculated the pivoting table for
396        ! the electrodes
397        i = Elecs(1)%pvt(Elecs(1)%t_dir)
398        j = Elecs(2)%pvt(Elecs(2)%t_dir)
399
400        bool = i == j
401
402        ! For a single transport direction to be true,
403        ! both the projections _has_ to be 1, exactly!
404        rtmp = VEC_PROJ_SCA(cell(:,i), Elecs(1)%cell(:,Elecs(1)%t_dir))
405        rtmp = rtmp / VNORM(Elecs(1)%cell(:,Elecs(1)%t_dir))
406        bool = bool .and. abs(abs(rtmp) - 1._dp) < 1.e-5_dp
407        rtmp = VEC_PROJ_SCA(cell(:,j), Elecs(2)%cell(:,Elecs(2)%t_dir))
408        rtmp = rtmp / VNORM(Elecs(2)%cell(:,Elecs(2)%t_dir))
409        bool = bool .and. abs(abs(rtmp) - 1._dp) < 1.e-5_dp
410
411        if ( bool ) then
412
413          ! The transport direction for the electrodes are the same...
414          ! And fully encompassed! We have a single transport
415          ! direction.
416          ts_tidx = i
417
418        else
419
420          ! In case we have a skewed transport direction
421          ! we have some restrictions...
422          ts_tidx = -N_Elec
423
424        end if
425      end if
426    end if
427
428    ! Populate the electrodes in the chemical potential type
429    do i = 1 , N_Elec
430       err = .true.
431       do j = 1 , N_mu
432          if ( associated(Elecs(i)%mu,target=mus(j)) ) then
433             call chem_pot_add_Elec(mus(j),i)
434             err = .false.
435             exit
436          end if
437       end do
438       if ( err ) then
439          call die('We could not attribute a chemical potential &
440               &to electrode: '//trim(Elecs(i)%name))
441       end if
442    end do
443
444    ! check that all electrodes and chemical potentials are paired in
445    ! some way.
446    if ( any(mus(:)%N_El == 0) ) then
447       call die('A/Some chemical potential(s) has/have not been assigned any electrodes. &
448            &All chemical potentials *MUST* be assigned an electrode')
449    end if
450
451    if ( na_u <= sum(TotUsedAtoms(Elecs)) ) then
452       write(*,'(a)') 'Please stop this madness. What where you thinking?'
453       call die('Electrodes occupy the entire device!!!')
454    end if
455
456  end subroutine read_tbt_elec
457
458
459  subroutine read_tbt_after_Elec(nspin, cell, na_u, lasto, xa, no_u, kscell, kdispl)
460
461    use fdf, only: fdf_get, leqi
462    use parallel, only: IONode
463
464    use m_ts_method, only: ts_method, TS_BTD
465    use m_ts_method, only: ts_A_method, TS_BTD_A_COLUMN, TS_BTD_A_PROPAGATION
466
467    use m_tbt_contour, only: read_contour_options
468
469    use m_tbt_sigma_save, only: init_Sigma_options
470    use m_tbt_dH, only: init_dH_options
471    use m_tbt_dSE, only: init_dSE_options
472
473    ! *******************
474    ! * INPUT variables *
475    ! *******************
476    integer, intent(in) :: nspin
477    real(dp), intent(in) :: cell(3,3)
478    integer,  intent(in) :: na_u, lasto(0:na_u)
479    real(dp), intent(in) :: xa(3,na_u)
480    integer, intent(in) :: no_u
481    integer,  intent(in) :: kscell(3,3)
482    real(dp), intent(in) :: kdispl(3)
483
484    integer :: i
485    logical :: ltmp, Gamma3(3), only_T_Gf
486    character(len=150) :: chars
487
488    ! we must have read the electrodes first
489    if ( N_Elec == 0 ) call die('read_tbt_options: Error in programming')
490
491    percent_tracker = fdf_get('TBT.Progress',5.)
492    percent_tracker = max(0., percent_tracker)
493
494    ! Reading the Transiesta solution method
495    chars = fdf_get('TBT.SolutionMethod','BTD')
496    if ( leqi(chars,'BTD').or.leqi(chars,'tri') ) then
497      ts_method = TS_BTD
498    else
499      call die('Unrecognized TBtrans solution method: '//trim(chars))
500    end if
501
502    chars = fdf_get('TS.BTD.Optimize','speed')
503    chars = fdf_get('TBT.BTD.Optimize',trim(chars))
504    if ( leqi(chars,'speed') ) then
505      BTD_method = 0
506    else if ( leqi(chars,'memory') ) then
507      BTD_method = 1
508    else
509      call die('Could not determine flag TBT.BTD.Optimize, please &
510          &see manual.')
511    end if
512
513    ! Read spectral calculation method for BTD method
514    if ( N_Elec > 3 ) then
515      chars = fdf_get('TS.BTD.Spectral','column')
516    else
517      chars = fdf_get('TS.BTD.Spectral','propagation')
518    end if
519    chars = fdf_get('TBT.BTD.Spectral',trim(chars))
520    if ( leqi(chars,'column') ) then
521      ts_A_method = TS_BTD_A_COLUMN
522    else if ( leqi(chars,'propagation') ) then
523      ts_A_method = TS_BTD_A_PROPAGATION
524    else
525      call die('TBT.BTD.Spectral option is not column or propagation. &
526          &Please correct input.')
527    end if
528
529    ! initial
530    only_T_Gf = .true.
531
532    ! Whether we should assert and calculate
533    ! all transmission amplitudes
534    ltmp = fdf_get('TBT.T.Elecs.All',N_Elec == 1)
535    ltmp = fdf_get('TBT.T.All',N_Elec == 1)
536    if ( ltmp ) then
537      save_DATA = save_DATA // ('T-all'.kv.1)
538    end if
539
540    N_eigen = fdf_get('TBT.T.Eig',0)
541    if ( N_eigen > 0 ) then
542      save_DATA = save_DATA // ('T-eig'.kv.N_eigen)
543      only_T_Gf = .false.
544    end if
545
546    ! Should we calculate DOS of electrode bulk Green function
547    ltmp = fdf_get('TBT.DOS.Elecs', .false. )
548    if ( ltmp ) then
549      save_DATA = save_DATA // ('DOS-Elecs'.kv.1)
550    end if
551
552    ltmp = fdf_get('TBT.T.Bulk', N_Elec == 1)
553    if ( ltmp ) then
554      ! when calculating the DOS for the electrode
555      ! we also get the bulk transmission.
556      save_DATA = save_DATA // ('DOS-Elecs'.kv.1)
557    end if
558
559    ! Should we calculate DOS of Green function
560    ltmp = fdf_get('TBT.DOS.Gf', N_Elec == 1)
561    if ( ltmp ) then
562      save_DATA = save_DATA // ('DOS-Gf'.kv.1)
563      only_T_Gf = .false.
564    end if
565
566    ! Should we calculate DOS of spectral function
567    ltmp = fdf_get('TBT.DOS.A', N_Elec == 1)
568    if ( ltmp ) then
569      save_DATA = save_DATA // ('DOS-A'.kv.1)
570      only_T_Gf = .false.
571    end if
572
573    ! Should we calculate DOS of all spectral functions
574    ltmp = fdf_get('TBT.DOS.A.All', N_Elec == 1)
575    if ( ltmp ) then
576      save_DATA = save_DATA // ('DOS-A'.kv.1)
577      save_DATA = save_DATA // ('DOS-A-all'.kv.1)
578      only_T_Gf = .false.
579    end if
580
581    ! Should we calculate orbital current
582    ltmp = fdf_get('TBT.Current.Orb', .false. )
583    if ( ltmp ) then
584      save_DATA = save_DATA // ('DOS-A'.kv.1)
585      save_DATA = save_DATA // ('orb-current'.kv.1)
586      only_T_Gf = .false.
587    end if
588
589    ! Options for density-matrix calculations
590    ltmp = fdf_get('TBT.DM.Gf', .false.)
591    if ( ltmp ) then
592      save_DATA = save_DATA // ('DOS-Gf'.kv.1)
593      save_DATA = save_DATA // ('DM-Gf'.kv.1)
594    end if
595
596    ltmp = fdf_get('TBT.DM.A', .false.)
597    if ( ltmp ) then
598      save_DATA = save_DATA // ('DOS-A'.kv.1)
599      save_DATA = save_DATA // ('DM-A'.kv.1)
600    end if
601
602
603    ! Options for COOP and COHP curves.
604    ! These are orbital (energy) populations that can be used to determine the
605    ! bonding nature of the material.
606    ltmp = fdf_get('TBT.COOP.Gf', .false.)
607    if ( ltmp ) then
608      save_DATA = save_DATA // ('DOS-Gf'.kv.1)
609      save_DATA = save_DATA // ('COOP-Gf'.kv.1)
610    end if
611
612    ltmp = fdf_get('TBT.COOP.A', .false. )
613    if ( ltmp ) then
614      save_DATA = save_DATA // ('DOS-A'.kv.1)
615      save_DATA = save_DATA // ('COOP-A'.kv.1)
616    end if
617
618    ltmp = fdf_get('TBT.COHP.Gf', .false. )
619    if ( ltmp ) then
620      save_DATA = save_DATA // ('DOS-Gf'.kv.1)
621      save_DATA = save_DATA // ('COHP-Gf'.kv.1)
622    end if
623
624    ltmp = fdf_get('TBT.COHP.A', .false. )
625    if ( ltmp ) then
626      save_DATA = save_DATA // ('DOS-A'.kv.1)
627      save_DATA = save_DATA // ('COHP-A'.kv.1)
628    end if
629
630    ltmp = fdf_get('TBT.T.Out',N_Elec == 1)
631    if ( ltmp ) then
632      save_DATA = save_DATA // ('T-sum-out'.kv.1)
633    end if
634
635    ! We cannot calculate the transmission for more than 3
636    ! electrodes if using the diagonal
637    if ( only_T_Gf .and. N_Elec > 3 ) then
638      only_T_Gf = .false.
639    end if
640
641    if ( only_T_Gf ) then
642      ! TODO, consider changing this to .true.
643      only_T_Gf = fdf_get('TBT.T.Gf',.false.)
644    end if
645    if ( only_T_Gf ) then
646      save_DATA = save_DATA // ('T-Gf'.kv.1)
647      ts_A_method = TS_BTD_A_COLUMN
648
649      ! We get sum-out for free, so save it
650      save_DATA = save_DATA // ('T-sum-out'.kv.1)
651
652      if ( N_Elec > 2 ) then
653        ! When having more than one electrode
654        ! we *must* calculate all transmissions
655        ! to be able to get the actual transmissions
656        ! using the diagonal Green function
657        save_DATA = save_DATA // ('T-all'.kv.1)
658      end if
659
660    end if
661
662    call init_Sigma_options( save_DATA )
663
664    if ( IONode ) write(*,*) ! new-line
665    ! The init_dH_options also checks for the
666    ! consecutive sparse patterns.
667    call init_dH_options( no_u )
668    call init_dSE_options( )
669
670
671    ! read in contour options
672    call read_contour_options( N_Elec, Elecs, N_mu, mus )
673
674    ! Check for Gamma in each direction
675    do i = 1 , 3
676      if ( kdispl(i) /= 0._dp ) then
677        Gamma3(i) = .false.
678      else if ( sum(kscell(i,:)) > 1 ) then
679        ! Note it is the off-diagonal for this unit-cell
680        ! direction
681        Gamma3(i) = .false.
682      else
683        Gamma3(i) = .true.
684      end if
685    end do
686
687    do i = 1 , N_Elec
688      ! Initialize the electrode quantities for the stored values
689      call check_Elec_sim(Elecs(i), nspin, cell, na_u, xa, &
690          Elecs_xa_EPS, lasto, Gamma3)
691    end do
692
693  end subroutine read_tbt_after_Elec
694
695  subroutine print_tbt_options(nspin)
696
697    use units, only: Kelvin, eV
698    use parallel, only: IONode
699    use files, only: slabel
700
701    use m_ts_method, only: ts_A_method, TS_BTD_A_COLUMN, TS_BTD_A_PROPAGATION
702
703    use m_tbt_contour, only: print_contour_tbt_options, io_contour_tbt
704    use m_tbt_contour, only: print_contour_tbt_block
705    use m_tbt_save, only: print_save_options
706    use m_tbt_diag, only: print_diag
707    use m_tbt_dH, only: print_dH_options
708    use m_tbt_dSE, only: print_dSE_options
709    use m_tbt_sigma_save, only: print_Sigma_options
710    use m_tbt_proj, only: print_proj_options
711    use m_tbt_hs, only: Volt, IsVolt, spin_idx
712
713    integer, intent(in) :: nspin
714
715    integer :: i
716
717    if ( .not. IONode ) return
718
719    write(*,*)
720    write(*,f11) repeat('*', 62)
721
722    write(*,f6) 'Electronic temperature (reference)',kT/Kelvin,'K'
723    if ( IsVolt ) then
724       write(*,f6) 'Voltage', Volt/eV,'Volts'
725    else
726       write(*,f11) 'No applied bias'
727    end if
728    write(*,f1) 'Calculate transmission only using diag(Gf)',('T-Gf'.in.save_DATA)
729    write(*,f1) 'Saving bulk transmission for electrodes',('DOS-Elecs'.in.save_DATA)
730    write(*,f1) 'Saving DOS from bulk electrodes',('DOS-Elecs'.in.save_DATA)
731    write(*,f1) 'Saving DOS from Green function',('DOS-Gf'.in.save_DATA)
732    if ( 'DOS-A-all' .in. save_DATA ) then
733       write(*,f1) 'Saving DOS from all spectral functions',.true.
734    else
735       write(*,f1) 'Saving DOS from spectral functions',('DOS-A' .in. save_DATA)
736    end if
737    write(*,f1) 'Saving bond currents (orb-orb)',('orb-current'.in.save_DATA)
738
739    ! DM
740    write(*,f1) 'Saving DM from Green function',('DM-Gf'.in.save_DATA)
741    write(*,f1) 'Saving DM from spectral functions',('DM-A'.in.save_DATA)
742
743    ! COOP/COHP curves
744    write(*,f1) 'Saving COOP from Green function',('COOP-Gf'.in.save_DATA)
745    write(*,f1) 'Saving COOP from spectral functions',('COOP-A'.in.save_DATA)
746    write(*,f1) 'Saving COHP from Green function',('COHP-Gf'.in.save_DATA)
747    write(*,f1) 'Saving COHP from spectral functions',('COHP-A'.in.save_DATA)
748
749    write(*,f12) 'Calc. # transmission eigenvalues',N_eigen
750    write(*,f1) 'Calc. T between all electrodes',('T-all'.in.save_DATA)
751    write(*,f1) 'Calc. total T out of electrodes',('T-sum-out'.in.save_DATA)
752    if ( nspin > 1 ) then
753       if ( spin_idx == 0 ) then
754          write(*,f11) 'Calculate spin UP and DOWN'
755       else if ( spin_idx == 1 ) then
756          write(*,f10) 'Calculate spin ','UP'
757       else if ( spin_idx == 2 ) then
758          write(*,f10) 'Calculate spin ','DOWN'
759       else
760          call die('Error in spin_idx')
761       end if
762    else
763       write(*,f11) 'Non-polarized Hamiltonian'
764    end if
765
766
767    ! Algorithm choices...
768    select case ( BTD_method )
769    case ( 0 )
770       write(*,f10)'BTD creation algorithm', 'speed'
771    case ( 1 )
772       write(*,f10)'BTD creation algorithm', 'memory'
773    end select
774    select case ( ts_A_method )
775    case ( TS_BTD_A_PROPAGATION )
776       write(*,f10)'BTD spectral function algorithm','propagation'
777    case ( TS_BTD_A_COLUMN )
778       write(*,f10)'BTD spectral function algorithm','column'
779    case default
780       call die('Error in setup BTD. A calc')
781    end select
782
783
784    call print_diag()
785    call print_Sigma_options( save_DATA )
786    call print_dH_options( )
787    call print_dSE_options( )
788
789    call print_save_options()
790    call print_proj_options( save_DATA )
791
792    write(*,f11)'          >> Electrodes << '
793    do i = 1 , size(Elecs)
794       call print_settings(Elecs(i),'tbt')
795    end do
796
797    call print_contour_tbt_options( 'TBT' )
798
799    write(*,f11) repeat('*', 62)
800    write(*,*)
801
802    call io_contour_tbt(slabel)
803
804    write(*,f11) repeat('<', 62)
805
806    call print_contour_tbt_block( 'TBT' )
807
808    write(*,f11) repeat('<', 62)
809
810  end subroutine print_tbt_options
811
812  subroutine print_tbt_warnings( Gamma )
813
814    use fdf, only: fdf_get, fdf_defined
815    use units, only: eV
816    use parallel, only: IONode
817
818    use m_tbt_save, only: save_parallel
819
820    use m_tbt_dH, only: print_dH_warnings
821    use m_tbt_hs, only: Volt
822
823    ! Whether the user requests a Gamma calculation
824    logical, intent(in) :: Gamma
825
826    integer :: i
827    logical :: ltmp, rem_DOS_Elecs, rem_T_Gf, has
828
829    ! Removal of keys
830    rem_DOS_Elecs = 'DOS-Elecs' .in. save_DATA
831    if ( rem_DOS_Elecs ) then
832      ! Currently the TBTGF files does not contain the DOS
833      if ( all(Elecs(:)%out_of_core) ) &
834          call delete(save_DATA,'DOS-Elecs')
835      if ( 'Sigma-only' .in. save_DATA ) &
836          call delete(save_DATA,'DOS-Elecs')
837    end if
838
839    rem_T_Gf = 'T-Gf' .in. save_DATA
840    if ( N_Elec > 3 ) then
841      call delete(save_DATA,'T-Gf')
842    end if
843
844    if ( .not. IONode ) return
845
846    write(*,*) ! new-line
847    write(*,'(3a)') repeat('*',24),' Begin: TBT CHECKS AND WARNINGS ',repeat('*',24)
848
849    if ( N_eigen < 0 ) then
850       call die('Number of transmission eigenvalues MUST be &
851            &zero or positive.')
852    end if
853
854    if ( .not. Gamma ) then
855       ltmp = .not. fdf_get('SpinSpiral',.false.)
856       ltmp = fdf_get('TBT.Symmetry.TimeReversal',ltmp)
857       if ( ('orb-current' .in.save_DATA) ) then
858          if ( IONode .and. ltmp ) then
859             write(*,'(a,/,a)') 'WARNING: k-averaging orbital currents with &
860                  &time-reversal symmetry will not reproduce','the correct &
861                  &orbital current. Set TBT.Symmetry.TimeReversal F'
862          end if
863       end if
864
865       ! Also for COOP analysis
866       has = 'COOP-Gf' .in. save_DATA
867       has = has .or. ('COOP-A' .in. save_DATA)
868       has = has .or. ('COHP-Gf' .in. save_DATA)
869       has = has .or. ('COHP-A' .in. save_DATA)
870       if ( IONode .and. ltmp .and. has ) then
871          write(*,'(a,/,a)') 'WARNING: k-averaging COOP/COHP with &
872               &time-reversal symmetry will not reproduce','the correct &
873               &populations. Set TBT.Symmetry.TimeReversal F'
874       end if
875    end if
876
877    ! CHECK THIS (we could allow it by only checking the difference...)
878    if (  maxval(mus(:)%mu) - minval(mus(:)%mu) - abs(Volt) > 1.e-8_dp * eV ) then
879       if ( IONode ) then
880          write(*,'(a)') 'Chemical potentials [eV]:'
881          do i = 1 , N_Elec
882             write(*,'(a,f10.5,a)') trim(Name(Elecs(i)))//' at ',Elecs(i)%mu%mu/eV,' eV'
883          end do
884          write(*,'(a)') 'The difference must satisfy: "max(ChemPots)-min(ChemPots) - abs(Volt) < 1e-8 eV"'
885          write(*,'(a,f10.5,a)') 'max(ChemPots) at ', maxval(mus(:)%mu)/eV,' eV'
886          write(*,'(a,f10.5,a)') 'min(ChemPots) at ', minval(mus(:)%mu)/eV,' eV'
887          write(*,'(a,f10.5,a)') '|V| at ', abs(Volt)/eV,' eV'
888       end if
889       call die('Chemical potentials are not consistent with the bias applied.')
890    end if
891
892
893    ! Check that the bias does not introduce a gating
894    if ( any(abs(mus(:)%mu) - abs(Volt) > 1.e-9_dp ) ) then
895       write(*,'(a)') 'Chemical potentials must lie in the range [-V;V] with the maximum &
896            &difference being V'
897       call die('Chemical potentials must not introduce consistent Ef shift to the system.')
898    end if
899
900    if ( rem_DOS_Elecs ) then
901       if ( any(Elecs(:)%out_of_core) ) then
902          write(*,'(a)')' Disabling electrode DOS calculation, only &
903               &enabled for in-core self-energy calculations.'
904       end if
905       if ( 'Sigma-only' .in. save_DATA ) then
906          write(*,'(a)')' Disabling electrode DOS calculation, only &
907               &enabled when calculating transmission (not TBT.SelfEnergy.Only).'
908       end if
909    end if
910
911    if ( rem_T_Gf .and. N_Elec > 3 ) then
912       write(*,'(a)')' ** Disabling transport calculation using diagonal, &
913            &not possible with N_elec > 3.'
914    end if
915
916    do i = 1, N_Elec
917      if ( Elecs(i)%repeat .and. Elecs(i)%bloch%size() > 1 ) then
918        write(*,'(a)') 'Electrode '//trim(Elecs(i)%name)//' is &
919            &using Bloch unfolding using the repeat scheme! &
920            &Please use the tiling scheme (it is orders of magnitudes faster!).'
921      end if
922    end do
923
924#ifdef MPI
925#ifdef NCDF_PARALLEL
926    if ( .not. save_parallel ) then
927       write(*,'(a)') ' ** Speed up the execution by utilizing parallel I/O'
928       write(*,'(a)') '  > TBT.CDF.MPI true'
929    end if
930#endif
931#endif
932
933    call print_dH_warnings( save_DATA )
934
935    write(*,'(3a,/)') repeat('*',24),' End: TBT CHECKS AND WARNINGS ',repeat('*',26)
936
937  end subroutine print_tbt_warnings
938
939end module m_tbt_options
940