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 9module m_ts_options 10 11 use precision, only : dp 12 13 use m_mixing, only: tMixer 14 15 use m_ts_electype, only : Elec 16 use m_ts_chem_pot, only : ts_mu 17 use m_ts_tdir 18 19 implicit none 20 21 public 22 save 23 24 ! ###### SIESTA-options ###### 25 ! The following options override the siesta settings upon 26 ! entering the transiesta SCF 27 28 ! The tolerance 29 real(dp) :: ts_Dtol ! = tolerance for density matrix 30 real(dp) :: ts_Htol ! = tolerance for Hamiltonian 31 logical :: ts_converge_dQ = .true. ! whether we should converge the charge 32 real(dp) :: ts_dQtol ! = tolerance for charge in the device region (after SCF) 33 integer :: ts_hist_keep = 0 34 35 ! ###### end SIESTA-options ##### 36 37 ! Whether we should stop before transiesta begins... 38 logical :: TS_siesta_stop = .false. 39 40 !< Whether TranSiesta is allowed to start without the 0-bias calculation. 41 logical :: TS_start_bias = .false. 42 43 ! Controls to save the TSHS file 44 logical :: TS_HS_save = .true. 45 logical :: TS_DE_save = .false. 46 ! whether we will use the bias-contour 47 logical :: IsVolt = .false. 48 ! maximum difference between chemical potentials 49 real(dp) :: Volt = 0._dp 50 ! The temperature for transiesta calculations 51 real(dp) :: ts_kT 52 ! Electrodes and different chemical potentials 53 integer :: N_Elec = 0 54 type(Elec), allocatable, target :: Elecs(:) 55 integer :: N_mu = 0 56 type(ts_mu), allocatable, target :: mus(:) 57 58 integer :: TS_scf_mode = 0 59 60 ! Flag to control whether we should update the forces (i.e. calculate energy-density matrix) 61 logical :: Calc_Forces = .true. 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 ! File name for reading in the grid for the Hartree potential 69 character(len=150) :: Hartree_fname = ' ' 70 71 ! A quantity describing the accuracy of the coordinates of the 72 ! electrodes. 73 ! * Should only be edited by experienced users * 74 real(dp) :: Elecs_xa_EPS = 1.e-4_dp 75 76 ! The user can request to analyze the system, returning information about the 77 ! tri-diagonalization partition and the contour 78 logical :: TS_Analyze = .false. 79 80 ! List of private formats for printing information 81 character(len=*), parameter, private :: f1 ='(''ts: '',a,t53,''='',tr4,l1)' 82 character(len=*), parameter, private :: f10='(''ts: '',a,t53,''='',tr4,a)' 83 character(len=*), parameter, private :: f11='(''ts: '',a)' 84 character(len=*), parameter, private :: f5 ='(''ts: '',a,t53,''='',i5,a)' 85 character(len=*), parameter, private :: f20='(''ts: '',a,t53,''='',i0,'' -- '',i0)' 86 character(len=*), parameter, private :: f6 ='(''ts: '',a,t53,''='',f10.4,tr1,a)' 87 character(len=*), parameter, private :: f7 ='(''ts: '',a,t53,''='',f12.6,tr1,a)' 88 character(len=*), parameter, private :: f8 ='(''ts: '',a,t53,''='',f10.4)' 89 character(len=*), parameter, private :: f9 ='(''ts: '',a,t53,''='',tr1,e9.3)' 90 character(len=*), parameter, private :: f15='(''ts: '',a,t53,''='',2(tr1,i0,'' x''),'' '',i0)' 91 92 ! set copy of SCF mixing 93 type(tMixer), pointer :: ts_scf_mixs(:) => null() 94 95contains 96 97 98 ! We have separated the options routines to only deal 99 ! with their respective parts 100 101 ! > Read generic options for transiesta which has 102 ! nothing to do with the electrodes or the chemical potentials 103 subroutine read_ts_generic( cell ) 104 105 use fdf, only : fdf_get, leqi 106 use intrinsic_missing, only : VNORM 107 108 use siesta_options, only : dDtol, dHtol 109 110 use m_ts_global_vars, only: TSmode, onlyS 111 use m_ts_method, only: TS_FULL, TS_BTD, TS_MUMPS, ts_method 112 113 use m_ts_weight, only : read_ts_weight 114 use ts_dq_m, only : ts_dq_read 115 116#ifdef SIESTA__MUMPS 117 use m_ts_mumps_init, only : read_ts_mumps 118#endif 119 120 use m_mixing, only: mixers_init 121 use m_mixing_scf, only: scf_mixs 122 123 ! Input variables 124 real(dp), intent(in) :: cell(3,3) 125 126 ! Local variables 127 character(len=200) :: chars 128 129 ! This has to be the first routine to be read 130 if ( N_mu /= 0 ) call die('read_ts_generic: error in programming') 131 if ( N_Elec /= 0 ) call die('read_ts_generic: error in programming') 132 133 ! Read in general values that should be used in the electrode generation 134 ! I.e. these FDF-parameters are used for diagon runs with transiesta 135#ifdef TRANSIESTA 136 TS_HS_save = fdf_get('TS.HS.Save',.true.) 137 TS_DE_save = fdf_get('TS.DE.Save',.true.) 138#else 139 TS_HS_save = fdf_get('TS.HS.Save',.false.) 140 TS_DE_save = fdf_get('TS.DE.Save',.false.) 141#endif 142 onlyS = fdf_get('TS.onlyS',.false.) 143 onlyS = fdf_get('TS.S.Save',onlyS) 144 145 ! Immediately return from transiesta when this occurs 146 ! no settings from the intrinsic transiesta routines 147 ! are needed. 148 if ( onlyS .or. .not. TSmode ) return 149 150 ! When running TSmode we FORCE TS.HS.Save and TS.DE.Save 151 TS_HS_save = .true. 152 TS_DE_save = .true. 153 154 ! Force the run of a biased TranSiesta run 155 ! from a pristine siesta calculation. 156 TS_start_bias = fdf_get('TS.Voltage.FromSiesta', .false.) 157 158 ! Read in the transiesta SCF mixing options 159 call mixers_init('TS.SCF', ts_scf_mixs ) 160 if ( .not. associated(ts_scf_mixs) ) then 161 ts_scf_mixs => scf_mixs 162 end if 163 164 ! Read in the mixing for the transiesta cycles 165 ts_Dtol = fdf_get('TS.SCF.DM.Tolerance',dDTol) 166 ts_Htol = fdf_get('TS.SCF.H.Tolerance',dHTol) 167 ts_hist_keep = fdf_get('TS.SCF.Mix.History.Keep',0) 168 169 ! Stop after siesta has converged 170 TS_siesta_stop = fdf_get('TS.SIESTA.Only',.false.) 171 172 ! Reading the Transiesta solution method 173 chars = fdf_get('TS.SolutionMethod','BTD') 174 if ( leqi(chars,'full') ) then 175 ts_method = TS_FULL 176 else if ( leqi(chars,'BTD') .or. leqi(chars,'tri') ) then 177 ts_method = TS_BTD 178#ifdef SIESTA__MUMPS 179 else if ( leqi(chars,'mumps') ) then 180 ts_method = TS_MUMPS 181#endif 182 else 183 call die('Unrecognized TranSiesta solution method: '//trim(chars)) 184 end if 185 186 ! currently this does not work 187 chars = fdf_get('SCF.Initialize','diagon') 188 chars = fdf_get('TS.SCF.Initialize',chars) 189 if ( leqi(chars,'diagon') ) then 190 TS_scf_mode = 0 191 chars = 'none' 192 else if ( leqi(chars,'transiesta') ) then 193 TS_scf_mode = 1 194 chars = 'init' 195 end if 196 197 chars = fdf_get('TS.BTD.Optimize','speed') 198 if ( leqi(chars,'speed') ) then 199 BTD_method = 0 200 else if ( leqi(chars,'memory') ) then 201 BTD_method = 1 202 else 203 call die('Could not determine flag TS.BTD.Optimize, please & 204 &see manual.') 205 end if 206 207 ! Determine whether the user wishes to only do an analyzation 208 TS_Analyze = fdf_get('TS.Analyze',.false.) 209 210 ! Read charge-correction methods 211 call ts_dq_read( ) 212 213#ifdef SIESTA__MUMPS 214 call read_ts_mumps( ) 215#endif 216 217 call read_ts_weight( ) 218 219 ! whether to calculate the forces or not (default calculate everything) 220 Calc_Forces = fdf_get('TS.Forces',.true.) 221 222 end subroutine read_ts_generic 223 224 225 ! > Reads the chemical potentials as well as the applied 226 ! Bias. 227 ! The bias is an intricate part of the chemical potential why it 228 ! is read in here. 229 subroutine read_ts_chem_pot( ) 230 231 use fdf, only : fdf_get 232 use units, only: eV, Kelvin 233 234 use siesta_options, only : kT => Temp 235 236 use m_ts_global_vars, only: TSmode, onlyS 237 238 use m_ts_chem_pot, only : fdf_nmu, fdffake_mu, fdf_mu, name 239 240 implicit none 241 242! ******************* 243! * LOCAL variables * 244! ******************* 245 logical :: err 246 integer :: i, j 247 real(dp) :: rtmp 248 249 if ( onlyS .or. .not. TSmode ) return 250 251 ts_kT = fdf_get('TS.ElectronicTemperature',kT,'Ry') 252 if ( ts_kT / Kelvin < 10._dp ) then 253 call die('TranSiesta electronic temperature *must* & 254 &be larger than 10 kT') 255 end if 256 257 ! The sign can not be chosen from this (several mu, where to define it) 258 Volt = fdf_get('TS.Voltage',0._dp,'Ry') 259 ! Voltage situation is above 0.01 mV 260 IsVolt = abs(Volt) > 0.00001_dp * eV 261 262 ! Read in the chemical potentials 263 N_mu = fdf_nmu('TS',ts_kT,mus) 264 err = .true. 265 if ( N_mu < 1 ) then 266 err = .false. 267 N_mu = fdffake_mu(mus,ts_kT,Volt) 268 end if 269 do i = 1 , N_mu 270 ! Default things that could be of importance 271 if ( err .and. .not. fdf_mu('TS',mus(i),ts_kT,Volt) ) then 272 call die('Could not find chemical potential: ' & 273 //trim(name(mus(i)))) 274 end if 275 276 ! We do not allow the electronic temperature 277 ! to be below 10 kT 278 if ( mus(i)%kT / Kelvin < 10._dp ) then 279 call die('TranSiesta electronic temperature *must* & 280 &be larger than 10 kT') 281 end if 282 283 end do 284 285 ! We consider 10 Kelvin to be the minimum allowed 286 ! temperature difference of the leads. 287 rtmp = 10._dp * Kelvin 288 do j = 1 , N_mu - 1 289 do i = j + 1 , N_mu 290 if ( abs(mus(i)%kT - mus(j)%kT) > rtmp ) then 291 ! If there exists a temperature gradient 292 ! we are in non-equilibrium, hence we need the 293 ! bias-setup no matter V! 294 IsVolt = .true. 295 exit 296 end if 297 end do 298 end do 299 300 end subroutine read_ts_chem_pot 301 302 303 ! Reads all information regarding the electrodes, nothing more. 304 subroutine read_ts_elec( cell, na_u, xa, lasto) 305 306 use fdf, only : fdf_get, fdf_obsolete, fdf_deprecated, leqi 307 use parallel, only : IONode 308 use intrinsic_missing, only : IDX_SPC_PROJ, EYE, VNORM 309 use intrinsic_missing, only : VEC_PROJ_SCA 310 311 use m_os, only : file_exist 312 313 use files, only: slabel 314 use units, only: eV, Ang 315 316 use m_ts_global_vars, only: TSmode, onlyS 317 318 use m_ts_chem_pot, only : copy, chem_pot_add_Elec 319 320 use m_ts_electype, only : fdf_nElec, fdf_Elec 321 use m_ts_electype, only : Name, TotUsedOrbs, TotUsedAtoms 322 use m_ts_electype, only : init_Elec_sim 323 324 use m_ts_method, only : ts_init_electrodes, a_isBuffer 325 use m_cite, only : add_citation 326 327 implicit none 328 329! ******************* 330! * INPUT variables * 331! ******************* 332 real(dp), intent(in) :: cell(3,3) 333 integer, intent(in) :: na_u, lasto(0:na_u) 334 real(dp), intent(in) :: xa(3,na_u) 335 336! ******************* 337! * LOCAL variables * 338! ******************* 339 integer :: i, j 340 real(dp) :: rtmp 341 logical :: err, bool 342 character(len=200) :: chars, c 343 type(ts_mu) :: tmp_mu 344 345 if ( onlyS .or. .not. TSmode ) return 346 347 if ( N_mu == 0 ) call die('read_ts_elecs: error in programming') 348 349 ! the title of the green's functions are now non-generic 350 call fdf_obsolete('TS.GFTitle') 351 352 ! The regular options for describing the electrodes can not be 353 ! used anymore... 354 call fdf_obsolete('TS.HSFileLeft') 355 call fdf_obsolete('TS.GFFileLeft') 356 call fdf_obsolete('TS.NumUsedAtomsLeft') 357 call fdf_obsolete('TS.ReplicateA1Left') 358 call fdf_obsolete('TS.ReplicateA2Left') 359 call fdf_obsolete('TS.HSFileRight') 360 call fdf_obsolete('TS.GFFileRight') 361 call fdf_obsolete('TS.NumUsedAtomsRight') 362 call fdf_obsolete('TS.ReplicateA1Right') 363 call fdf_obsolete('TS.ReplicateA2Right') 364 365 ! notice that this does not have the same meaning... 366 call fdf_deprecated('TS.UpdateDMCROnly','TS.Elecs.DM.Update') 367 call fdf_deprecated('TS.UseBulk','TS.Elecs.Bulk') 368 369 ! whether or not the electrodes should be re-instantiated 370 call fdf_deprecated('TS.CalcGF','TS.Elecs.GF.ReUse') 371 call fdf_deprecated('TS.ReUseGF','TS.Elecs.GF.ReUse') 372 373 ! To determine the same coordinate nature of the electrodes 374 Elecs_xa_EPS = fdf_get('TS.Elecs.Coord.Eps',0.001_dp*Ang, 'Bohr') 375 376 ! detect how many electrodes we have 377 N_Elec = fdf_nElec('TS',Elecs) 378 if ( N_Elec < 1 ) then 379 ! We initialize to 2 electrodes (Left/Right) 380 N_Elec = 2 381 allocate(Elecs(N_Elec)) 382 Elecs(1)%name = 'Left' 383 Elecs(1)%ID = 1 384 Elecs(2)%name = 'Right' 385 Elecs(2)%ID = 2 386 ! if they do-not exist, the user will be told 387 if ( IONode ) then 388 c = '(''transiesta: *** '',a,/)' 389 write(*,c)'No electrode names were found, default Left/Right are expected' 390 end if 391 end if 392 393 ! If only one electrode you are not allowed to move the Fermi-level 394 ! of the electrode. That should be done by other means (i.e. use NetCharge) 395 if ( N_Elec == 1 ) then 396 ! Notice that below the chemical potential gets corrected 397 ! EVEN if the user supplied a bias. 398 if ( IsVolt .and. IONode ) then 399 c = '(''transiesta: *** '',a)' 400 write(*,c) 'Single electrode calculations does not allow shifting the chemical potential.' 401 write(*,c) 'You should do that by changing the states filled in the system.' 402 write(*,c) 'Consult the manual on how to do this.' 403 call die('Please set the chemical potential to zero for your one electrode') 404 end if 405 end if 406 407 ! Setup default parameters for the electrodes 408 ! first electrode is the "left" 409 ! last electrode is the "right" 410 ! the remaining electrodes have their chemical potential at 0 411 ! We should probably warn if +2 electrodes are used and t_dir is the 412 ! same for all electrodes... Then the user needs to know what (s)he is doing... 413 ! Accuracy required for self-energy convergence 414 Elecs(:)%accu = fdf_get('TS.Elecs.Accuracy',1.e-13_dp*eV,'Ry') 415 Elecs(:)%Eta = fdf_get('TS.Elecs.Eta',0.001_dp*eV,'Ry') 416 Elecs(:)%Bulk = fdf_get('TS.Elecs.Bulk',.true.) ! default everything to bulk electrodes 417 if ( Elecs(1)%Bulk ) then 418 ! Default is cross-terms if we use bulk electrodes 419 chars = 'cross-terms' 420 else 421 ! For non-bulk systems, we default to all 422 chars = 'all' 423 end if 424 c = fdf_get('TS.Elecs.DM.Update',trim(chars)) 425 if ( leqi(c,'none') ) then 426 Elecs(:)%DM_update = 0 427 else if ( leqi(c,'cross-terms') .or. & 428 leqi(c,'cross-term') ) then 429 Elecs(:)%DM_update = 1 430 else if ( leqi(c,'all') ) then 431 Elecs(:)%DM_update = 2 432 else 433 call die('TS.Elecs.DM.Update [cross-terms,none,all]: & 434 &unrecognized option: '//trim(c)) 435 end if 436 437 ! Whether we should always set the DM to bulk 438 ! values (by reading in from electrode DM) 439 if ( TS_scf_mode == 1 .and. .not. IsVolt ) then 440 chars = 'bulk' 441 else 442 chars = 'diagon' 443 end if 444 chars = fdf_get('TS.Elecs.DM.Init',trim(chars)) 445 if ( leqi(chars,'diagon') ) then 446 Elecs(:)%DM_init = 0 447 else if ( leqi(chars,'bulk') ) then 448 Elecs(:)%DM_init = 1 449 else if ( leqi(chars,'force-bulk') ) then 450 Elecs(:)%DM_init = 2 451 else 452 call die('TS.Elecs.DM.Init unknown value [diagon,bulk,force-bulk]') 453 end if 454 if ( IsVolt ) then 455 if ( Elecs(1)%DM_init == 1 .and. IONode) then 456 if ( IONode ) then 457 c = '(''transiesta: *** '',a,/)' 458 write(*,c)'Will default to not read in electrode DM, only applicable for V = 0 calculations' 459 end if 460 end if 461 Elecs(:)%DM_init = 0 462 end if 463 464 ! Whether we should try and re-use the surface Green function 465 ! files 466 Elecs(:)%ReUseGF = fdf_get('TS.Elecs.GF.ReUse',.true.) 467 468 ! whether all calculations should be performed 469 ! "out-of-core" i.e. whether the GF files should be created or not 470 Elecs(:)%out_of_core = fdf_get('TS.Elecs.Out-of-core',.true.) 471 472 do i = 1 , N_Elec 473 474 ! If we only have 2 electrodes we take them 475 ! as though the atomic indices are the first and last 476 ! respectively. 477 if ( N_Elec == 2 ) then 478 if ( i == 1 ) then 479 err = fdf_Elec('TS',slabel,Elecs(i),N_mu,mus,idx_a= 1) 480 else if ( i == 2 ) then 481 err = fdf_Elec('TS',slabel,Elecs(i),N_mu,mus,idx_a=-1) 482 end if 483 else 484 ! Default things that could be of importance 485 err = fdf_Elec('TS',slabel,Elecs(i),N_mu,mus) 486 end if 487 if ( .not. err ) then 488 call die('Could not find electrode: '//trim(name(Elecs(i)))) 489 end if 490 491 if ( Elecs(i)%idx_a < 0 ) & 492 Elecs(i)%idx_a = na_u + Elecs(i)%idx_a + 1 493 if ( Elecs(i)%idx_a < 1 .or. & 494 na_u < Elecs(i)%idx_a ) then 495 print *,Elecs(i)%idx_a,na_u 496 call die("Electrode position does not exist") 497 end if 498 if ( N_Elec == 2 ) then 499 ! Correct for buffer atoms, first electrode steps "up" 500 ! second electrode steps "down" 501 if ( i == 1 ) then 502 j = Elecs(i)%idx_a 503 do while ( a_isBuffer(j) ) 504 j = j + 1 505 end do 506 Elecs(i)%idx_a = j 507 else 508 j = Elecs(i)%idx_a + TotUsedAtoms(Elecs(i)) - 1 509 do while ( a_isBuffer(j) ) 510 j = j - 1 511 end do 512 Elecs(i)%idx_a = j - TotUsedAtoms(Elecs(i)) + 1 513 end if 514 end if 515 ! set the placement in orbitals 516 Elecs(i)%idx_o = lasto(Elecs(i)%idx_a-1)+1 517 518 ! Initialize the electrode quantities for the 519 ! stored values 520 call init_Elec_sim(Elecs(i), cell, na_u, xa) 521 522 end do 523 524 ! If many electrodes, no transport direction can be specified 525 ! Hence we use this as an error-check (also for N_Elec == 1) 526 if ( any(Elecs(:)%t_dir > 3) ) then 527 ts_tidx = - N_Elec 528 529 ! We add the real-space self-energy article 530 if ( IONode ) then 531 call add_citation("10.1103/PhysRevB.100.195417") 532 end if 533 534 else 535 select case ( N_Elec ) 536 case ( 1 ) 537 ! The easy case 538 ! We simple need to figure out if the electrode 539 ! has its transport direction aligned with the 540 ! lattice vectors 541 542 i = Elecs(1)%pvt(Elecs(1)%t_dir) 543 544 ! For a single transport direction to be true, 545 ! both the projections _has_ to be 1, exactly! 546 rtmp = VEC_PROJ_SCA(cell(:,i), Elecs(1)%cell(:,Elecs(1)%t_dir)) 547 rtmp = rtmp / VNORM(Elecs(1)%cell(:,Elecs(1)%t_dir)) 548 bool = abs(abs(rtmp) - 1._dp) < 1.e-5_dp 549 550 if ( bool ) then 551 552 ! The transport direction for the electrodes are the same... 553 ! And fully encompassed! We have a single transport 554 ! direction. 555 ts_tidx = i 556 557 else 558 559 ! In case we have a skewed transport direction 560 ! we have some restrictions... 561 ts_tidx = - N_Elec 562 563 end if 564 565 case ( 2 ) 566 567 ! Retrieve the indices of the unit-cell directions 568 ! according to the electrode transport directions. 569 ! We have already calculated the pivoting table for 570 ! the electrodes 571 i = Elecs(1)%pvt(Elecs(1)%t_dir) 572 j = Elecs(2)%pvt(Elecs(2)%t_dir) 573 bool = i == j 574 575 ! For a single transport direction to be true, 576 ! both the projections _has_ to be 1, exactly! 577 rtmp = VEC_PROJ_SCA(cell(:,i), Elecs(1)%cell(:,Elecs(1)%t_dir)) 578 rtmp = rtmp / VNORM(Elecs(1)%cell(:,Elecs(1)%t_dir)) 579 bool = bool .and. abs(abs(rtmp) - 1._dp) < 1.e-5_dp 580 rtmp = VEC_PROJ_SCA(cell(:,j), Elecs(2)%cell(:,Elecs(2)%t_dir)) 581 rtmp = rtmp / VNORM(Elecs(2)%cell(:,Elecs(2)%t_dir)) 582 bool = bool .and. abs(abs(rtmp) - 1._dp) < 1.e-5_dp 583 584 if ( bool ) then 585 586 ! The transport direction for the electrodes are the same... 587 ! And fully encompassed! We have a single transport 588 ! direction. 589 ts_tidx = i 590 591 else 592 593 ! In case we have a skewed transport direction 594 ! we have some restrictions... 595 ts_tidx = - N_Elec 596 597 end if 598 599 case default 600 601 ! N_Elec > 2 602 ! Here we always have these settings 603 ts_tidx = - N_Elec 604 605 end select 606 end if 607 608 ! The user can selectively decide how the bias 609 ! is applied. 610 ! For N-terminal calculations we advice the user 611 ! to use a Poisson solution they add. 612 if ( ts_tidx > 0 ) then 613 ! We have a single unified semi-inifinite direction 614 chars = fdf_get('TS.Poisson','ramp') 615 else 616 chars = fdf_get('TS.Poisson','elec-box') 617 end if 618 619#ifdef NCDF_4 620 if ( file_exist(chars, Bcast = .true.) ) then 621 622 Hartree_fname = trim(chars) 623 ts_tidx = 0 624 625 else 626#endif 627 Hartree_fname = ' ' 628 if ( leqi(chars,'ramp') .or. leqi(chars, 'ramp-cell') ) then 629 if ( ts_tidx <= 0 ) then 630 call die('TS.Poisson cannot be ramp for & 631 &anything but 2-electrodes with aligned transport direction.') 632 end if 633 else if ( leqi(chars,'elec-box') ) then 634 ts_tidx = - N_Elec 635 else 636#ifdef NCDF_4 637 call die('Error in specifying how the Hartree potential & 638 &should be placed. [ramp|elec-box|NetCDF-file]') 639#else 640 call die('Error in specifying how the Hartree potential & 641 &should be placed. [ramp|elec-box]') 642#endif 643 end if 644#ifdef NCDF_4 645 end if 646#endif 647 648 ! In case we are doing equilibrium, fix the chemical potential to the first 649 if ( .not. IsVolt ) then 650 651 ! force it to be zero... can be necessary if considering single electrode 652 ! calculations (assures V == 0) 653 Volt = 0._dp 654 655 ! copy over electrode... 656 call copy(mus(1),tmp_mu) 657 658 ! Deallocate all strings 659 do i = 1 , N_mu 660 deallocate(mus(i)%Eq_seg) 661 end do 662 deallocate(mus) 663 664 ! create the first chemical potential again 665 N_mu = 1 666 allocate(mus(1)) 667 call copy(tmp_mu,mus(1)) 668 deallocate(tmp_mu%Eq_seg) 669 670 ! Firmly assure the chemical potential to be zero 671 mus(1)%mu = 0._dp 672 mus(1)%ID = 1 673 674 ! Assign all electrodes to the same chemical potential 675 do i = 1 , N_Elec 676 Elecs(i)%mu => mus(1) 677 end do 678 679 end if 680 681 ! Populate the electrodes in the chemical potential type 682 do i = 1 , N_Elec 683 err = .true. 684 do j = 1 , N_mu 685 if ( associated(Elecs(i)%mu,target=mus(j)) ) then 686 call chem_pot_add_Elec(mus(j),i) 687 err = .false. 688 exit 689 end if 690 end do 691 if ( err ) then 692 call die('We could not attribute a chemical potential & 693 &to electrode: '//trim(Elecs(i)%name)) 694 end if 695 end do 696 697 if ( na_u <= sum(TotUsedAtoms(Elecs)) ) then 698 write(*,'(a)') 'Please stop this madness. What where you thinking?' 699 write(*,*) na_u, sum(TotUsedAtoms(Elecs)) 700 call die('Electrodes occupy the entire device!!!') 701 end if 702 703 ! Initialize the electrode regions 704 call ts_init_electrodes(na_u,lasto,N_Elec,Elecs) 705 706 end subroutine read_ts_elec 707 708 709 ! This routine reads options for transiesta after 710 ! having read in the electrodes 711 ! It allows one to do things between reading options 712 subroutine read_ts_after_Elec( cell, nspin, na_u, xa, lasto, & 713 ts_kscell, ts_kdispl) 714 715 use fdf, only : fdf_get, leqi 716 use intrinsic_missing, only : VNORM, IDX_SPC_PROJ, EYE 717 use atomlist, only : qa 718 use siesta_options, only : charnet 719 720 use m_ts_global_vars, only: TSmode, onlyS 721 722 use m_ts_method, only: TS_BTD_A_PROPAGATION, TS_BTD_A_COLUMN 723 use m_ts_method, only: ts_A_method, a_isBuffer 724 725 use m_ts_electype, only: check_Elec_sim 726 727 use m_ts_contour, only: read_contour_options 728 use m_ts_contour_eq, only: N_Eq, Eq_c 729 730 use m_ts_weight, only : read_ts_weight 731 use ts_dq_m, only : ts_dq_read 732 733 use m_ts_hartree, only: read_ts_hartree_options 734 735#ifdef SIESTA__MUMPS 736 use m_ts_mumps_init, only : read_ts_mumps 737#endif 738 739 ! Input variables 740 real(dp), intent(in) :: cell(3,3) 741 integer, intent(in) :: nspin, na_u 742 real(dp), intent(in) :: xa(3,na_u) 743 integer, intent(in) :: lasto(0:na_u) 744 integer, intent(in) :: ts_kscell(3,3) 745 real(dp), intent(in) :: ts_kdispl(3) 746 747 ! Local variables 748 character(len=50) :: chars 749 integer :: i 750 real(dp) :: dev_q 751 logical :: Gamma3(3) 752 753 if ( onlyS .or. .not. TSmode ) return 754 755 if ( N_Elec == 0 ) call die('read_ts_after_Elecs: error in programming') 756 757 ! Read spectral calculation method for BTD method 758 if ( N_Elec > 3 ) then 759 chars = fdf_get('TS.BTD.Spectral','column') 760 else 761 chars = fdf_get('TS.BTD.Spectral','propagation') 762 end if 763 if ( leqi(chars,'column') ) then 764 ts_A_method = TS_BTD_A_COLUMN 765 else if ( leqi(chars,'propagation') ) then 766 ts_A_method = TS_BTD_A_PROPAGATION 767 else 768 call die('TS.BTD.Spectral option is not column or propagation. & 769 &Please correct input.') 770 end if 771 772 ! Read in options again, at this point we have 773 ! the correct ts_tidx 774 call read_ts_hartree_options() 775 776 ! read in contour options 777 call read_contour_options( N_Elec, Elecs, N_mu, mus, ts_kT, IsVolt, Volt ) 778 779 ! Quickly update the forces update if continued fraction 780 ! is used, one need an extra Green function evaluation to 781 ! get the forces correct (iR and -R) 782 do i = 1 , N_Eq 783 if ( leqi(Eq_c(i)%c_io%part,'cont-frac') ) then 784 ! the forces are not updated, dispite the user requests 785 calc_forces = .false. 786 end if 787 end do 788 789 ! Check for Gamma in each direction 790 do i = 1 , 3 791 if ( ts_kdispl(i) /= 0._dp ) then 792 Gamma3(i) = .false. 793 else if ( sum(ts_kscell(i,:)) > 1 ) then 794 ! Note it is the off-diagonal for this direction 795 Gamma3(i) = .false. 796 else 797 Gamma3(i) = .true. 798 end if 799 end do 800 801 do i = 1 , N_Elec 802 ! Initialize the electrode quantities for the 803 ! stored values 804 call check_Elec_sim(Elecs(i), nspin, cell, na_u, xa, & 805 Elecs_xa_EPS, lasto, Gamma3, ts_kscell, ts_kdispl) 806 end do 807 808 ! Now we know which atoms are buffer's and electrodes 809 ! This allows us to decide the tolerance for convergence of 810 ! the charges. 811 ! By default we allow a difference up to q(device) / 1000 812 dev_q = - charnet 813 do i = 1, na_u 814 if ( .not. a_isBuffer(i) ) then 815 ! count this atom 816 dev_q = dev_q + qa(i) 817 end if 818 end do 819 ts_dQtol = fdf_get('TS.SCF.dQ.Tolerance',dev_q / 1000._dp) 820 ts_converge_dQ = fdf_get('TS.SCF.dQ.Converge', .true.) 821 822 end subroutine read_ts_after_Elec 823 824 825 subroutine print_ts_options( cell ) 826 827 use fdf, only: fdf_get, leqi 828 use parallel, only: IOnode 829 830 use intrinsic_missing, only : IDX_SPC_PROJ, EYE 831 use units, only: eV, Kelvin 832 833 use m_mixing, only: mixers_print 834 use m_mixing_scf, only: scf_mixs 835 836 use m_ts_electype, only: print_settings 837 838 use m_ts_global_vars, only: TSmode, onlyS 839 840 use m_ts_contour, only: print_contour_options 841 842 use m_ts_method, only: TS_FULL, TS_BTD, TS_MUMPS, ts_method, na_Buf 843 use m_ts_method, only: TS_BTD_A_COLUMN, TS_BTD_A_PROPAGATION 844 use m_ts_method, only: ts_A_method 845 846 use ts_dq_m, only: TS_DQ_METHOD, TS_DQ_METHOD_BUFFER, TS_DQ_METHOD_FERMI 847 use ts_dq_m, only: TS_DQ_FACTOR, TS_DQ_FERMI_TOLERANCE 848 use ts_dq_m, only: TS_DQ_FERMI_MAX, TS_DQ_FERMI_ETA 849 850 use m_ts_weight, only: TS_W_METHOD, TS_W_CORRELATED 851 use m_ts_weight, only: TS_W_ORB_ORB, TS_W_TR_ATOM_ATOM, TS_W_SUM_ATOM_ATOM 852 use m_ts_weight, only: TS_W_TR_ATOM_ORB, TS_W_SUM_ATOM_ORB, TS_W_MEAN 853 use m_ts_weight, only: TS_W_K_METHOD 854 use m_ts_weight, only: TS_W_K_CORRELATED, TS_W_K_UNCORRELATED 855 856#ifdef SIESTA__MUMPS 857 use m_ts_mumps_init, only: MUMPS_mem, MUMPS_ordering, MUMPS_block 858#endif 859 860 use m_ts_hartree, only: TS_HA_PLANES, TS_HA_frac, TS_HA_offset 861 862 implicit none 863 864 ! The unit-cell 865 real(dp), intent(in) :: cell(3,3) 866 867! ******************* 868! * LOCAL variables * 869! ******************* 870 character(len=200) :: chars 871 logical :: ltmp 872 real(dp) :: tmp33(3,3) 873 integer :: i, tdir 874 875 if ( .not. IONode ) return 876 877 write(*,*) 878 write(*,f11) repeat('*', 62) 879 880 if ( onlyS .or. .not. TSmode ) then 881 882 write(*,f1) 'Save H and S matrices', TS_HS_save 883 write(*,f1) 'Save DM and EDM matrices', TS_DE_save 884 write(*,f1) 'Only save the overlap matrix S', onlyS 885 886 write(*,f11) repeat('*', 62) 887 write(*,*) 888 889 return 890 end if 891 892 write(*,f6) 'Voltage', Volt/eV,'Volts' 893 write(*,f1) 'Save H and S matrices', TS_HS_save 894 if ( TS_Analyze ) then 895 write(*,f11)'Will analyze bandwidth of LCAO sparse matrix and quit' 896 end if 897 write(*,f7) 'Electronic temperature',ts_kT/Kelvin,'K' 898 if ( ts_tidx < 1 ) then 899 write(*,f11) 'Transport individually selected for electrodes' 900 else 901 write(chars,'(a,i0)') 'A',ts_tidx 902 write(*,f10) 'Transport along unit-cell vector',trim(chars) 903 ! Calculate Cartesian transport direction 904 call eye(3,tmp33) 905 tdir = IDX_SPC_PROJ(tmp33,cell(:,ts_tidx),mag=.true.) 906 select case ( tdir ) 907 case ( 1 ) 908 write(*,f10) 'Transport along Cartesian vector','X' 909 case ( 2 ) 910 write(*,f10) 'Transport along Cartesian vector','Y' 911 case ( 3 ) 912 write(*,f10) 'Transport along Cartesian vector','Z' 913 end select 914 end if 915 chars = ' ' 916 if ( TS_HA_PLANES(1, 1) ) chars = trim(chars) // '-A' 917 if ( TS_HA_PLANES(2, 1) ) chars = trim(chars) // '+A' 918 if ( TS_HA_PLANES(1, 2) ) chars = trim(chars) // '-B' 919 if ( TS_HA_PLANES(2, 2) ) chars = trim(chars) // '+B' 920 if ( TS_HA_PLANES(1, 3) ) chars = trim(chars) // '-C' 921 if ( TS_HA_PLANES(2, 3) ) chars = trim(chars) // '+C' 922 write(*,f10) 'Fixing Hartree potential at cell boundary', trim(chars) 923 write(*,f8) 'Fix Hartree potential fraction', TS_HA_frac 924 write(*,f7) 'Hartree potential offset', TS_HA_offset/eV, 'eV' 925 926 if ( ts_method == TS_FULL ) then 927 write(*,f10)'Solution method', 'Full inverse' 928 else if ( ts_method == TS_BTD ) then 929 write(*,f10)'Solution method', 'BTD' 930 chars = fdf_get('TS.BTD.Pivot','atom+'//trim(Elecs(1)%name)) 931 write(*,f10)'BTD pivoting method', trim(chars) 932 if ( BTD_method == 0 ) then 933 chars = 'speed' 934 else if ( BTD_method == 1 ) then 935 chars = 'memory' 936 end if 937 write(*,f10)'BTD creation algorithm', trim(chars) 938 if ( IsVolt ) then 939 select case ( ts_A_method ) 940 case ( TS_BTD_A_PROPAGATION ) 941 write(*,f10)'BTD spectral function algorithm','propagation' 942 case ( TS_BTD_A_COLUMN ) 943 write(*,f10)'BTD spectral function algorithm','column' 944 case default 945 call die('Error in setup BTD. A calc') 946 end select 947 end if 948#ifdef SIESTA__MUMPS 949 else if ( ts_method == TS_MUMPS ) then 950 write(*,f10)'Solution method', 'MUMPS' 951 write(*,f5)'MUMPS extra memory', MUMPS_mem,'%' 952 write(*,f5)'MUMPS blocking factor', MUMPS_block,'' 953 select case ( MUMPS_ordering ) 954 case ( 7 ) 955 write(*,f10)'MUMPS ordering', 'auto' 956 case ( 6 ) 957 write(*,f10)'MUMPS ordering', 'QAMD' 958 case ( 5 ) 959 write(*,f10)'MUMPS ordering', 'METIS' 960 case ( 4 ) 961 write(*,f10)'MUMPS ordering', 'PORD' 962 case ( 3 ) 963 write(*,f10)'MUMPS ordering', 'SCOTCH' 964 case ( 2 ) 965 write(*,f10)'MUMPS ordering', 'AMF' 966 case ( 0 ) 967 write(*,f10)'MUMPS ordering', 'AMD' 968 end select 969#endif 970 end if 971 write(*,f9) 'SCF DM tolerance',ts_Dtol 972 write(*,f7) 'SCF Hamiltonian tolerance',ts_Htol/eV, 'eV' 973 write(*,f1) 'SCF converge charge',ts_converge_dQ 974 write(*,f9) 'SCF charge tolerance',ts_dQtol 975 976 select case ( TS_scf_mode ) 977 case ( 0 ) 978 write(*,f10) 'Initialize DM by','diagon' 979 case ( 1 ) 980 write(*,f10) 'Initialize DM by','transiesta' 981 end select 982 if ( IsVolt ) then 983 if ( len_trim(Hartree_fname) > 0 ) then 984 write(*,f10) 'User supplied Hartree potential', & 985 trim(Hartree_fname) 986 else 987 if ( ts_tidx > 0 ) then 988 write(*,f11) 'Hartree potential ramp across entire cell' 989 else 990 write(*,f11) 'Hartree potential will be placed in electrode box' 991 end if 992 end if 993 994 chars = 'Non-equilibrium contour weight method' 995 select case ( TS_W_METHOD ) 996 case ( TS_W_ORB_ORB ) 997 write(*,f10) trim(chars),'orb-orb' 998 case ( TS_W_CORRELATED + TS_W_TR_ATOM_ATOM ) 999 write(*,f10) trim(chars),'Correlated Tr[atom]-Tr[atom]' 1000 case ( TS_W_TR_ATOM_ATOM ) 1001 write(*,f10) trim(chars),'Uncorrelated Tr[atom]-Tr[atom]' 1002 case ( TS_W_CORRELATED + TS_W_TR_ATOM_ORB ) 1003 write(*,f10) trim(chars),'Correlated Tr[atom]-orb' 1004 case ( TS_W_TR_ATOM_ORB ) 1005 write(*,f10) trim(chars),'Uncorrelated Tr[atom]-orb' 1006 case ( TS_W_CORRELATED + TS_W_SUM_ATOM_ATOM ) 1007 write(*,f10) trim(chars),'Correlated Sum[atom]-Sum[atom]' 1008 case ( TS_W_SUM_ATOM_ATOM ) 1009 write(*,f10) trim(chars),'Uncorrelated Sum[atom]-Sum[atom]' 1010 case ( TS_W_CORRELATED + TS_W_SUM_ATOM_ORB ) 1011 write(*,f10) trim(chars),'Correlated Sum[atom]-orb' 1012 case ( TS_W_SUM_ATOM_ORB ) 1013 write(*,f10) trim(chars),'Uncorrelated Sum[atom]-orb' 1014 case ( TS_W_MEAN ) 1015 write(*,f10) trim(chars),'Algebraic mean' 1016 case default 1017 ! This is an easy place for cathing mistakes 1018 call die('Error in code, weighting method unrecognized.') 1019 end select 1020 chars = 'Non-equilibrium contour weight k-method' 1021 select case ( TS_W_K_METHOD ) 1022 case ( TS_W_K_CORRELATED ) 1023 write(*,f10) trim(chars),'Correlated k-points' 1024 case ( TS_W_K_UNCORRELATED ) 1025 write(*,f10) trim(chars),'Uncorrelated k-points' 1026 end select 1027 end if 1028 if ( .not. Calc_Forces ) then 1029 write(*,f11) '*** TranSiesta will NOT update forces ***' 1030 end if 1031 1032 if ( TS_DQ_METHOD == 0 ) then 1033 write(*,f11)'Will not correct charge fluctuations' 1034 else if ( TS_DQ_METHOD == TS_DQ_METHOD_BUFFER ) then ! Correct in buffer 1035 if ( 0 < na_Buf ) then 1036 write(*,f10)'Charge correction','buffer' 1037 else 1038 call die('Charge correction can not happen in buffer as no buffer & 1039 &atoms exist.') 1040 end if 1041 write(*,f8)'Charge correction factor',TS_DQ_FACTOR 1042 else if ( TS_DQ_METHOD == TS_DQ_METHOD_FERMI ) then ! Correct fermi-lever 1043 write(*,f10)'Charge correction','Fermi-level' 1044 write(*,f8)'Charge correction dQ tolerance',TS_DQ_FERMI_TOLERANCE 1045 write(*,f7)'Fermi-level extrapolation eta value ',TS_DQ_FERMI_ETA/eV, 'eV' 1046 write(*,f8)'Charge correction factor',TS_DQ_FACTOR 1047 write(*,f7)'Max change in Fermi-level allowed', & 1048 TS_DQ_FERMI_MAX / eV,'eV' 1049 end if 1050 1051 ! Print mixing options 1052 if ( associated(ts_scf_mixs, target=scf_mixs) ) then 1053 write(*,f11)'TS.SCF mixing options same as SCF' 1054 else 1055 call mixers_print('TS.SCF', ts_scf_mixs) 1056 end if 1057 1058 write(*,f11)' >> Electrodes << ' 1059 ltmp = ts_tidx < 1 .and. IsVolt 1060 do i = 1 , size(Elecs) 1061 call print_settings(Elecs(i), 'ts', & 1062 box = ltmp) 1063 end do 1064 1065 ! Print the contour information 1066 call print_contour_options( 'TS' , IsVolt ) 1067 1068 write(*,f11) repeat('*', 62) 1069 write(*,*) 1070 1071 end subroutine print_ts_options 1072 1073 1074 subroutine print_ts_warnings( Gamma, cell, na_u, xa, Nmove ) 1075 1076 use parallel, only: IONode, Nodes 1077 use intrinsic_missing, only : VNORM, VEC_PROJ_SCA 1078 1079 use m_os, only: file_exist 1080 1081 use units, only: Kelvin, eV, Ang 1082 use siesta_options, only: FixSpin 1083 1084 use m_ts_global_vars, only: TSmode, onlyS 1085 use m_ts_chem_pot, only : Name, Eq_segs 1086 use m_ts_electype, only : TotUsedAtoms, Name, Elec_frac 1087 1088 use m_ts_method, only : a_isElec, a_isBuffer 1089 use m_ts_method, only : ts_A_method, TS_BTD_A_COLUMN 1090 1091 use m_ts_contour_eq, only: N_Eq_E 1092 use m_ts_contour_neq, only: contour_neq_warnings 1093 1094 use m_ts_hartree, only: TS_HA_frac 1095 1096 ! Input variables 1097 logical, intent(in) :: Gamma 1098 real(dp), intent(in) :: cell(3,3) 1099 integer, intent(in) :: na_u 1100 real(dp), intent(in) :: xa(3,na_u) 1101 integer, intent(in) :: Nmove 1102 1103 ! Local variables 1104 integer :: i, j, iEl, idx, idx1, idx2, itmp3(3) 1105 real(dp) :: rtmp, tmp3(3), tmp33(3,3), bdir(2) 1106 real(dp) :: p(3) 1107 logical :: err, warn, ltmp 1108 1109 if ( .not. IONode ) return 1110 1111 warn = .false. 1112 err = .false. 1113 1114 write(*,'(3a)') repeat('*',24),' Begin: TS CHECKS AND WARNINGS ',repeat('*',24) 1115 if ( FixSpin ) then 1116 if ( TSmode ) then 1117 write(*,'(a)') 'Fixed spin for transiesta calculations is not implemented!' 1118 call die('Fixing spin is not possible in TranSiesta') 1119 end if 1120 if ( TS_HS_save ) then 1121 write(*,'(a)') 'Fixed spin aligns the Fermi-levels in the output TSHS to spin-UP!' 1122 end if 1123 if ( TS_DE_save ) then 1124 write(*,'(a)') 'Fixed spin aligns the Fermi-levels in the output TSDE to spin-UP!' 1125 end if 1126 end if 1127 1128 if ( TS_HA_frac /= 1._dp ) then 1129 write(*,'(a)') 'Fraction of Hartree potential is NOT 1.' 1130 warn = .true. 1131 end if 1132 if ( TS_HA_frac < 0._dp .or. 1._dp < TS_HA_frac ) then 1133 write(*,'(a)') 'Fraction of Hartree potential is below 0.' 1134 write(*,'(a)') ' MUST be in range [0;1]' 1135 call die('Vha fraction erronously set.') 1136 end if 1137 1138 ! Return if not a transiesta calculation 1139 if ( onlyS .or. .not. TSmode ) then 1140 write(*,'(3a,/)') repeat('*',24),' End: TS CHECKS AND WARNINGS ',repeat('*',26) 1141 return 1142 end if 1143 1144 if ( ts_tidx < 1 .and. len_trim(Hartree_fname) == 0 .and. IsVolt ) then 1145 write(*,'(a)') 'Hartree potiental correction is the box solution & 1146 &which is not advised. Please supply your own Poisson solution.' 1147 end if 1148 1149 if ( ts_A_method == TS_BTD_A_COLUMN ) then 1150 write(*,'(a)') 'Memory usage can be reduced by setting:' 1151 write(*,'(a)') ' TS.BTD.Spectral propagation' 1152 end if 1153 1154 1155 ! Check that all chemical potentials are really different 1156 ! their energy difference has to be below 0.1 meV and the 1157 ! temperature difference has to be below 10 K. 1158 rtmp = 10._dp * Kelvin 1159 do i = 1 , N_mu - 1 1160 do j = i + 1 , N_mu 1161 if ( abs(mus(i)%mu - mus(j)%mu) < 0.0001_dp * eV .and. & 1162 abs(mus(i)%kT - mus(j)%kT) <= rtmp ) then 1163 write(*,'(a)') 'Two chemical potentials: '//trim(name(mus(i)))//' and ' & 1164 //trim(name(mus(j)))//' are the same, in bias calculations this & 1165 &is not allowed.' 1166 err = .true. 1167 end if 1168 end do 1169 end do 1170 1171 ! Check that all chemical potentials are in use 1172 if ( any(mus(:)%N_El == 0) ) then 1173 write(*,'(a)') 'A/Some chemical potential(s) have not been assigned any electrodes. & 1174 &All chemical potentials *MUST* be assigned an electrode' 1175 err = .true. 1176 end if 1177 1178 ! check that all have at least 2 contour points on the equilibrium contour 1179 ! the 3rd is the fictive pole segment 1180 if ( .not. all(Eq_segs(mus(:)) > 0) ) then 1181 write(*,'(a)') 'All chemical potentials does not have at least & 1182 &1 equilibrium contours' 1183 err = .true. 1184 end if 1185 if ( .not. all(Eq_segs(mus(:)) /= 2) ) then 1186 write(*,'(a)') 'No chemical potential can have only two & 1187 &equilibrium contours' 1188 write(*,'(a)')'Either of these:' 1189 write(*,'(a)')' 1. continued fraction' 1190 write(*,'(a)')' or' 1191 write(*,'(a)')' 1. Circle contour' 1192 write(*,'(a)')' 2. Tail contour' 1193 write(*,'(a)')' 3. Residuals (poles)' 1194 err = .true. 1195 end if 1196 1197 ! we need to check that they indeed do not overlap 1198 do i = 1 , N_Elec 1199 idx1 = Elecs(i)%idx_a 1200 idx2 = idx1 + TotUsedAtoms(Elecs(i)) - 1 1201 ! we need to check every electrode, 1202 ! specifically because if one of the electrodes is fully located 1203 ! inside the other and we check the "small" one 1204 ltmp = .false. 1205 do j = 1 , N_Elec 1206 if ( i == j ) cycle 1207 idx = Elecs(j)%idx_a 1208 if ( (idx <= idx1 .and. & 1209 idx1 < idx + TotUsedAtoms(Elecs(j))) ) then 1210 ltmp = .true. 1211 end if 1212 if ( (idx <= idx2 .and. & 1213 idx2 < idx + TotUsedAtoms(Elecs(j))) ) then 1214 ltmp = .true. 1215 end if 1216 if ( ltmp ) then 1217 write(*,'(a)') 'Electrode: '//trim(Elecs(i)%name) 1218 write(*,'(a,i0,a,i0)') 'Positions: ',idx1,' -- ',idx2 1219 idx1 = Elecs(j)%idx_a 1220 idx2 = idx1 + TotUsedAtoms(Elecs(j)) - 1 1221 write(*,'(a)') 'Electrode: '//trim(Elecs(j)%name) 1222 write(*,'(a,i0,a,i0)') 'Positions: ',idx1,' -- ',idx2 1223 write(*,'(a)') 'Overlapping electrodes is not physical, please correct.' 1224 err = .true. 1225 end if 1226 end do 1227 1228 ! Warn if using non-bulk electrodes and one does not update everything! 1229 if ( .not. Elecs(i)%Bulk ) then 1230 select case ( Elecs(i)%DM_update ) 1231 case ( 0 ) 1232 write(*,'(a)') 'Electrode: '//trim(Elecs(i)%name) 1233 write(*,'(a)') ' is using a non-bulk electrode Hamiltonian and does not' 1234 write(*,'(a)') ' update the electrode region or the cross-term DM!' 1235 write(*,'(a)') ' Please consider setting "DM-update all"' 1236 case ( 1 ) 1237 write(*,'(a)') 'Electrode: '//trim(Elecs(i)%name) 1238 write(*,'(a)') ' is using a non-bulk electrode Hamiltonian and does not' 1239 write(*,'(a)') ' update the electrode region DM!' 1240 write(*,'(a)') ' Please consider setting "DM-update all"' 1241 end select 1242 end if 1243 1244 end do 1245 1246 ! CHECK THIS (we could allow it by only checking the difference...) 1247 if ( maxval(mus(:)%mu) - minval(mus(:)%mu) - abs(Volt) > 1.e-4_dp * eV ) then 1248 write(*,'(a)') 'Chemical potentials [eV]:' 1249 do i = 1 , N_Elec 1250 write(*,'(a,f10.5,a)') trim(Elecs(i)%name)//' at ',Elecs(i)%mu%mu/eV,' eV' 1251 end do 1252 write(*,'(a)') 'The difference must satisfy: "max(ChemPots)-min(ChemPots) - abs(Volt) < 1e-4 eV"' 1253 write(*,'(a,f10.5,a)') 'max(ChemPots) at ', maxval(mus(:)%mu)/eV,' eV' 1254 write(*,'(a,f10.5,a)') 'min(ChemPots) at ', minval(mus(:)%mu)/eV,' eV' 1255 write(*,'(a,f10.5,a)') '|V| at ', abs(Volt)/eV,' eV' 1256 write(*,'(a)') 'Chemical potentials are not consistent with the bias applied.' 1257 err = .true. 1258 end if 1259 1260 ! Check that the bias does not introduce a gating 1261 if ( any(abs(mus(:)%mu) - abs(Volt) > 1.e-9_dp ) ) then 1262 write(*,'(a)') 'Chemical potentials must lie in the range [-V;V] with the maximum & 1263 &difference being V' 1264 write(*,'(a)') 'Chemical potentials must not introduce consistent Ef shift to the system.' 1265 err = .true. 1266 end if 1267 1268 1269 ! Check that we can actually start directly in transiesta 1270 if ( TS_scf_mode == 1 ) then ! TS-start 1271 if ( .not. all(Elecs(:)%DM_update >= 1) ) then 1272 write(*,'(a)')'WARNING: Responsibility is now on your side' 1273 write(*,'(a)')'WARNING: Requesting immediate start, yet we & 1274 &do not update cross-terms.' 1275 warn = .true. 1276 end if 1277 end if 1278 1279 ! Calculate the number of optimal contour points 1280 i = mod(N_Eq_E(), Nodes) ! get remaining part of equilibrium contour 1281 if ( i /= 0 ) then 1282 i = Nodes - i 1283 write(*,'(a)')'Without loosing performance you can increase & 1284 &the equilibrium integration precision.' 1285 write(*,'(a,i0,a)')'You can add ',i,' more energy points in the & 1286 &equilibrium contours, for FREE!' 1287 if ( i/N_mu > 0 ) then 1288 write(*,'(a,i0,a)')'This is ',i/N_mu, & 1289 ' more energy points per chemical potential.' 1290 end if 1291 end if 1292 1293 call contour_nEq_warnings() 1294 1295 if ( .not. Calc_Forces ) then 1296 write(*,f11) '*** TranSiesta will NOT update forces ***' 1297 write(*,f11) '*** ALL FORCES AFTER TRANSIESTA HAS RUN ARE WRONG ***' 1298 if ( Nmove > 0 ) then 1299 write(*,'(a)')'Relaxation with TranSiesta *REQUIRES* an update of & 1300 &the energy density matrix. Will continue at your request.' 1301 err = .true. 1302 end if 1303 end if 1304 1305 if ( Nmove > 0 .and. .not. all(Elecs(:)%DM_update > 0) ) then 1306 write(*,'(a)') 'TranSiesta relaxation is only allowed if you also & 1307 &update, at least, the cross terms, please set: & 1308 &TS.Elecs.DM.Update [cross-terms|all]' 1309 err = .true. 1310 end if 1311 1312 ! A transiesta calculation requires that all atoms 1313 ! are within the unit-cell. 1314 ! However, for N_Elec > 2 calculations this 1315 ! is not a requirement as the potential is placed 1316 ! irrespective of the actual box. 1317 1318 ltmp = .false. 1319 call reclat(cell,tmp33,0) 1320 1321 do i = 1 , na_u 1322 1323 ! If we use a boxed or user-supplied potential 1324 ! profile, we need not check the coordinates. 1325 ! It is the users responsibility for user-defined 1326 ! potential grids. 1327 ! The boxed one takes care of grid-points in the periodic picture 1328 ! And if the grid-plane does not cut the actual grid, it will 1329 ! error out. 1330 if ( ts_tidx < 1 ) cycle 1331 1332 ! Buffer atoms can be where-ever 1333 if ( a_isBuffer(i) ) cycle 1334 1335 ! Check the index 1336 do j = 1 , 3 1337 itmp3(j) = floor( dot_product(xa(:,i),tmp33(:,j)) ) 1338 end do 1339 1340 ! Only check the transport direction 1341 ! Note that we have projected onto the unit-cell 1342 ! vector, hence tidx and not tdir 1343 if ( itmp3(ts_tidx) /= 0 ) then 1344 write(*,'(i0,''('',i0,''='',i0,'')'',tr1)',advance='no')& 1345 i,ts_tidx,itmp3(ts_tidx) 1346 ltmp = .true. 1347 end if 1348 1349 end do 1350 if ( ltmp ) then 1351 write(*,'(/a)')'atom(cell-dir=neighbour-cell)' 1352 write(*,'(a)') '*** Device atomic coordinates are not inside unit-cell.' 1353 write(*,'(a)') '*** This is a requirement for bias calculations' 1354 write(*,'(a)') ' as the Poisson equation cannot be correctly handled' 1355 write(*,'(a)') ' due to inconsistencies with the grid and atomic coordinates' 1356 if ( IsVolt ) then 1357 write(*,'(a)') 'Please move device atoms inside the device region in & 1358 &the transport direction.' 1359 err = .true. 1360 else 1361 write(*,'(a)') '*** Will continue, but will die when running V /= 0 ***' 1362 warn = .true. 1363 end if 1364 end if 1365 1366 if ( ts_tidx < 1 ) then 1367 1368 write(*,'(a)') '*** TranSiesta semi-infinite directions are individual ***' 1369 write(*,'(a)') '*** It is heavily advised to have any electrodes with no & 1370 &periodicity' 1371 write(*,'(a)') ' in the transverse directions be located as far from any & 1372 &cell-boundaries' 1373 write(*,'(a)') ' as possible. This has to do with the electrostatic potential & 1374 &correction. ***' 1375 if ( IsVolt ) then 1376 write(*,'(a)') '*** Please ensure electrode unit-cells are as confined as possible.' 1377 write(*,'(a)') ' I.e. do not add superfluous vacuum if not needed in the & 1378 &electrode calculation.' 1379 write(*,'(a)') ' The initial guess for the potential profile is heavily influenced' 1380 write(*,'(a)') ' by the electrode unit-cell sizes! ***' 1381 end if 1382 1383 else 1384 1385 ! The transport direction is well-defined 1386 ! and the Hartree potential is fixed at the bottom of the 1387 ! unit-cell of the A[ts_tidx] direction. 1388 ! We will let the user know if any atoms co-incide with 1389 ! the plane as that might hurt convergence a little. 1390 1391 if ( N_Elec <= 2 ) then 1392 1393 ! Get the electrode fraction of the position 1394 if ( N_Elec == 1 ) then 1395 1396 call Elec_frac(Elecs(1),cell,na_u,xa,ts_tidx, fmin = bdir(1)) 1397 call Elec_frac(Elecs(1),cell,na_u,xa,ts_tidx, fmax = bdir(2)) 1398 iEl = 1 1399 if ( bdir(1) < bdir(2) ) then 1400 i = 1 1401 else 1402 i = 2 1403 end if 1404 1405 else 1406 1407 ! Get the electrode fraction of the position 1408 call Elec_frac(Elecs(1),cell,na_u,xa,ts_tidx, fmin = bdir(1)) 1409 call Elec_frac(Elecs(2),cell,na_u,xa,ts_tidx, fmin = bdir(2)) 1410 1411 ! Determine the electrode closest to the 1412 ! lower boundary 1413 if ( bdir(1) < bdir(2) ) then 1414 ! The first electrode is closest 1415 i = 1 1416 iEl = 1 1417 call Elec_frac(Elecs(2),cell,na_u,xa,ts_tidx, fmax = bdir(2)) 1418 else 1419 ! The second electrode is closest 1420 i = 2 1421 iEl = 2 1422 call Elec_frac(Elecs(1),cell,na_u,xa,ts_tidx, fmax = bdir(1)) 1423 end if 1424 1425 end if 1426 1427 ! Get the fraction 1428 tmp3 = cell(:,ts_tidx) / VNORM(cell(:,ts_tidx)) 1429 1430 ! Check lower limit 1431 if ( bdir(i) < 0._dp ) then 1432 1433 ! Tell the user to shift the entire structure 1434 ! by the offset to origo + 1/2 a bond-length 1435 1436 ! Origo offset: 1437 p = - cell(:,ts_tidx) * bdir(i) 1438 p = p + tmp3 * Elecs(iEl)%dINF_layer * 0.5_dp 1439 ! Bond-length 1440 write(*,'(a,f9.5,a)') 'Electrode: '//trim(Elecs(iEl)%name)//' lies & 1441 &outside the unit-cell.' 1442 write(*,'(a)')'Please shift the entire structure using the & 1443 &following recipe:' 1444 write(*,'(a)') 'If you already have AtomicCoordinatesFormat, add these' 1445 write(*,'(tr1,a)') 'AtomicCoordinatesFormat Ang' 1446 write(*,'(tr1,a)') '%block AtomicCoordinatesOrigin' 1447 write(*,'(tr1,3(tr2,f12.4))') p / Ang 1448 write(*,'(tr1,a)') '%endblock AtomicCoordinatesOrigin' 1449 err = .true. 1450 1451 end if 1452 1453 ! Check upper limit 1454 if ( i == 1 ) then 1455 i = 2 1456 else 1457 i = 1 1458 end if 1459 if ( N_Elec == 2 ) iEl = i 1460 if ( 1._dp < bdir(i) ) then 1461 1462 bdir(i) = bdir(i) - 1._dp 1463 1464 ! Tell the user to shift the entire structure 1465 ! by the offset to origo + 1/2 a bond-length 1466 1467 ! Origo offset: 1468 p = - cell(:,ts_tidx) * bdir(i) 1469 p = p - tmp3 * Elecs(iEl)%dINF_layer * 0.5_dp 1470 ! Bond-length 1471 write(*,'(a,f9.5,a)') 'Electrode: '//trim(Elecs(iEl)%name)//' lies & 1472 &outside the unit-cell.' 1473 write(*,'(a)')'Please shift the entire structure using the & 1474 &following recipe:' 1475 write(*,'(a)') 'If you already have AtomicCoordinatesFormat, add these' 1476 write(*,'(tr1,a)') 'AtomicCoordinatesFormat Ang' 1477 write(*,'(tr1,a)') '%block AtomicCoordinatesOrigin' 1478 write(*,'(tr1,3(tr2,f12.4))') p / Ang 1479 write(*,'(tr1,a)') '%endblock AtomicCoordinatesOrigin' 1480 err = .true. 1481 1482 end if 1483 1484 else 1485 1486 call die("ts_options: ts_tidx < 0 with N_Elec > 2") 1487 1488 end if 1489 1490 end if 1491 1492 ! If the user has requested to initialize using transiesta 1493 ! and the user does not utilize the bulk DM, they should be 1494 ! warned 1495 if ( TS_scf_mode == 1 .and. any(Elecs(:)%DM_init == 0) ) then 1496 write(*,'(a)') 'You are not initializing the electrode DM/EDM. & 1497 &This may result in very wrong electrostatic potentials close to & 1498 &the electrode/device boundary region.' 1499 if ( IsVolt ) then 1500 write(*,'(a)') ' This warning is only applicable for V == 0 calculations!' 1501 write(*,'(a)') ' I give this warning because it is not clear how your V = 0 calcualtion was done.' 1502 end if 1503 warn = .true. 1504 end if 1505 1506 ! warn the user about suspicous work regarding the electrodes 1507 do i = 1 , N_Elec 1508 1509 idx1 = Elecs(i)%idx_a 1510 idx2 = idx1 + TotUsedAtoms(Elecs(i)) - 1 1511 1512 if ( .not. Elecs(i)%Bulk ) then 1513 write(*,'(a)') 'Electrode '//trim(Elecs(i)%name)//' will & 1514 ¬ use bulk Hamiltonian.' 1515 warn = .true. 1516 end if 1517 1518 if ( Elecs(i)%DM_update == 0 ) then 1519 write(*,'(a)') 'Electrode '//trim(Elecs(i)%name)//' will & 1520 ¬ update cross-terms or local region.' 1521 warn = .true. 1522 end if 1523 1524 if ( .not. Elecs(i)%Bulk .and. Elecs(i)%DM_update /= 2 ) then 1525 write(*,'(3a)') 'Electrode ',trim(Elecs(i)%name), & 1526 ' has non-bulk Hamiltonian and does not update all' 1527 warn = .true. 1528 end if 1529 1530 if ( .not. Elecs(i)%kcell_check ) then 1531 write(*,'(a)') 'Electrode '//trim(Elecs(i)%name)//' will & 1532 ¬ check the k-grid sampling vs. system k-grid & 1533 &sampling. Ensure appropriate sampling.' 1534 end if 1535 if ( Elecs(i)%V_frac_CT /= 0._dp ) then 1536 write(*,'(a)') 'Electrode '//trim(Elecs(i)%name)//' will & 1537 &shift coupling Hamiltonian with a shift in energy & 1538 &corresponding to the applied bias. Be careful here.' 1539 warn = .true. 1540 end if 1541 if ( Elecs(i)%delta_Ef /= 0._dp ) then 1542 write(*,'(a)') 'Electrode '//trim(Elecs(i)%name)//' will & 1543 &shift the electronic structure manually. Be careful here.' 1544 warn = .true. 1545 end if 1546 1547 if ( Elecs(i)%repeat .and. Elecs(i)%bloch%size() > 1 ) then 1548 write(*,'(a)') 'Electrode '//trim(Elecs(i)%name)//' is & 1549 &using Bloch unfolding using the repeat scheme! & 1550 &Please use the tiling scheme (it is orders of magnitudes faster!).' 1551 warn = .true. 1552 end if 1553 1554 ! if any buffer atoms exist, we should suggest to the user 1555 ! to use TS.Elec.<elec> [DM-update cross-terms|all] 1556 ! in case any buffer atoms are too close 1557 ltmp = .false. 1558 do j = 1 , na_u 1559 ! skip non-buffer atoms 1560 if ( .not. a_isBuffer(j) ) cycle 1561 do idx = idx1 , idx2 1562 ! Proximity of 4 Ang enables this check 1563 ltmp = VNORM(xa(:,idx)-xa(:,j)) < 4._dp * Ang 1564 if ( ltmp ) exit 1565 end do 1566 if ( ltmp ) exit 1567 end do 1568 if ( ltmp .and. Elecs(i)%DM_update == 0 ) then 1569 ! some buffer atoms are close to this electrode 1570 ! Advice to use dm_update 1571 write(*,'(a,/,a)') 'Electrode '//trim(Elecs(i)%name)//' is & 1572 &likely terminated by buffer atoms. It is HIGHLY recommended to add this:', & 1573 ' TS.Elec.'//trim(Elecs(i)%name)//' DM-update [cross-terms|all]' 1574 warn = .true. 1575 end if 1576 1577 ! In case DM_bulk is requested we assert that the file exists 1578 ltmp = file_exist(Elecs(i)%DEfile) 1579 if ( Elecs(i)%DM_init > 0 .and. .not. ltmp ) then 1580 write(*,'(a,/,a)') 'Electrode '//trim(Elecs(i)%name)//' TSDE & 1581 &file cannot be located in: '//trim(Elecs(i)%DEfile)//'.', & 1582 ' Please add TS.DE.Save T to the electrode calculation or & 1583 &specify the exact file position using ''TSDE-file'' in the& 1584 & TS.Elec block.' 1585 err = .true. 1586 end if 1587 1588 end do 1589 1590 if ( N_Elec /= 2 .and. any(Elecs(:)%DM_update == 0) ) then 1591 write(*,'(a,/,a)') 'Consider updating more elements when doing & 1592 &N-electrode calculations. The charge conservation typically & 1593 &increases.',' TS.Elecs.DM.Update [cross-terms|all]' 1594 warn = .true. 1595 end if 1596 1597 ! Check that the pivoting table is unique 1598 do iEl = 1, N_Elec 1599 if ( sum(Elecs(iEl)%pvt) /= 6 .or. count(Elecs(iEl)%pvt==2) /= 1 ) then 1600 write(*,'(a,/,a)') 'The pivoting table for the electrode unit-cell, & 1601 &onto the simulation unit-cell is not unique: '//trim(Elecs(iEl)%name), & 1602 ' Please check your electrode and device cell parameters!' 1603 write(*,'(a)') ' Combining this with electric fields or dipole-corrections is NOT advised!' 1604 warn = .true. 1605 end if 1606 end do 1607 1608 write(*,'(3a,/)') repeat('*',24),' End: TS CHECKS AND WARNINGS ',repeat('*',26) 1609 1610 if ( warn ) then 1611 ! Print BIG warning sign 1612 1613 write(*,'(tr18,a)') repeat('*',40) 1614 write(*,'(tr19,a)') 'TRANSIESTA REPORTED IMPORTANT WARNINGS' 1615 write(*,'(tr18,a)') repeat('*',40) 1616 1617 end if 1618 1619 if ( err ) then 1620 write(*,'(/tr18,a)') repeat('*',30) 1621 write(*,'(tr19,a)') 'TRANSIESTA REPORTED ERRORS' 1622 write(*,'(tr18,a)') repeat('*',30) 1623 1624 call die('One or more errors have occured doing TranSiesta & 1625 &initialization, check the output') 1626 end if 1627 1628 end subroutine print_ts_warnings 1629 1630 1631 subroutine print_ts_blocks( na_u, xa ) 1632 1633 use parallel, only : IONode 1634 use files, only : slabel 1635 1636 use m_ts_global_vars, only: TSmode, onlyS 1637 use m_ts_chem_pot, only: print_mus_block 1638 use m_ts_contour, only: print_contour_block, io_contour 1639 1640 1641 ! Input variables 1642 integer, intent(in) :: na_u 1643 real(dp), intent(in) :: xa(3,na_u) 1644 1645 ! Local variables 1646 integer :: i 1647 1648 if ( .not. IONode ) return 1649 1650 if ( onlyS .or. .not. TSmode ) return 1651 1652 write(*,'(/,a,/)') '>>> TranSiesta block information for FDF-file START <<<' 1653 1654 call print_mus_block( 'TS' , N_mu , mus) 1655 1656 call print_contour_block( 'TS' , IsVolt ) 1657 1658 write(*,'(/,a,/)') '>>> TranSiesta block information for FDF-file END <<<' 1659 1660 ! write out the contour 1661 call io_contour(IsVolt, mus, slabel) 1662 1663 end subroutine print_ts_blocks 1664 1665 subroutine val_swap(v1,v2) 1666 real(dp), intent(inout) :: v1, v2 1667 real(dp) :: tmp 1668 tmp = v1 1669 v1 = v2 1670 v2 = tmp 1671 end subroutine val_swap 1672 1673end module m_ts_options 1674