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! This code segment has been improved or fully created by: 9! Nick Papior Andersen, 2012, nickpapior@gmail.com 10! 11module m_ts_voltage 12 13 ! Module for containing various ways of implementing the voltage in a 14 ! TranSIESTA run. 15 ! As a standard method TranSIESTA applies the voltage ramp across the 16 ! entire unit cell. 17 ! As a future progress the voltage ramp might be situated in the contact region 18 ! only. This makes more physical sense as the voltage drop does not occur 19 ! in the electrode regions. 20 ! 21 ! Created and copyrighted by: Nick Papior Andersen, 2012 22 ! The use of this program is allowed for not-for-profit research only. 23 ! Copy or disemination of all or part of this package is not 24 ! permitted without prior and explicit authorization by the author. 25 26 use precision, only : dp 27 use m_ts_tdir 28 29 implicit none 30 31 private 32 33 ! The idea is to have sub routines in this module to do 34 ! various voltage layouts 35 public :: ts_voltage 36 public :: ts_init_voltage 37 38 ! The corresponding bias for either the left electrode 39 real(dp), save :: V_low = 0._dp 40 real(dp), save :: V_high = 0._dp 41 42contains 43 44 subroutine ts_init_voltage(cell, na_u, xa, nmesh) 45 use parallel, only : IONode 46 use m_ts_electype, only : TotUsedAtoms, Elec 47 use m_ts_options, only : Elecs, N_Elec 48 use m_ts_options, only : Volt, Hartree_fname 49 use units, only : eV, ang 50 51 use m_geom_box, only : voxel_in_box 52 53 ! *********************** 54 ! * INPUT variables * 55 ! *********************** 56 real(dp), intent(in) :: cell(3,3) 57 integer, intent(in) :: na_u 58 real(dp), intent(in) :: xa(3,na_u) 59 integer, intent(in) :: nmesh(3) ! total number of mesh-points. 60 61 logical :: bool 62 real(dp) :: tmp, ll(3), zero(3) 63 integer :: iElL, iElR, iEl, ia, iia 64 65 if ( IONode ) then 66 write(*,*) 67 end if 68 69 if ( ts_tidx < 1 ) then 70 71 if ( IONode .and. len_trim(Hartree_fname) == 0 ) then 72 write(*,'(a)')'ts-voltage: Lifted locally on each electrode' 73 write(*,'(a)')'ts-voltage: WARNING ****************************************' 74 write(*,'(a)')'ts-voltage: WARNING Prefer to use a custom Poisson solution!' 75 write(*,'(a)')'ts-voltage: WARNING ****************************************' 76 77 else if ( IONode .and. len_trim(Hartree_fname) > 0 ) then 78 write(*,'(3a)')'ts-voltage: ', & 79 'User supplied Poisson solution in file ', & 80 trim(Hartree_fname) 81#ifdef NCDF_4 82 call ts_ncdf_voltage_assert(Hartree_fname,cell,nmesh) 83#else 84 call die('NetCDF4 has not been compiled in, this functionality & 85 &is not supported.') 86#endif 87 end if 88 89 ! Find the lowest and highest chemical potential 90 V_low = Elecs(1)%mu%mu 91 V_high = Elecs(1)%mu%mu 92 do iEl = 2 , size(Elecs) 93 V_low = min(Elecs(iEl)%mu%mu,V_low) 94 V_high = max(Elecs(iEl)%mu%mu,V_high) 95 end do 96 if ( Volt < 0._dp ) then 97 ! with a negative bias, we have to reverse 98 ! the high-low 99 ! Remember that in this case V_high and V-low are 100 ! only used for the external Hartree potential. 101 tmp = V_high 102 V_high = V_low 103 V_low = tmp 104 end if 105 106 ! Check that all electrode atoms are residing in the boxes 107 ! defined by the electrodes 108 bool = .false. 109 zero = 0._dp 110 do iEl = 1 , N_Elec 111 do ia = 1 , TotUsedAtoms(Elecs(iEl)) 112 iia = Elecs(iEl)%idx_a - 1 + ia 113 ll = xa(:,iia) 114 if (.not.voxel_in_box(Elecs(iEl)%box,ll,zero)) then 115 if ( IONode ) & 116 write(*,'(3a,i0)')'Electrode ', & 117 trim(Elecs(iEl)%name),& 118 ' does not reside within the defined Hartree box & 119 &for the potential. Please ensure unit-cell & 120 &vectors have the _correct_ direction: ',iia 121 bool = .true. 122 end if 123 end do 124 end do 125 126 if ( bool .and. len_trim(Hartree_fname) == 0 ) then 127 call die('ts-voltage: Check output, an electrode cannot be & 128 &correctly applied a bias.') 129 end if 130 131 return 132 133 end if 134 135 ! set the left chemical potential 136 call get_elec_indices(Elecs, cell, na_u, xa, iElL, iElR) 137 138 ! this will be the applied bias in the "lower" 139 ! electrode in the unit-cell 140 V_high = Elecs(iElL)%mu%mu 141 V_low = Elecs(iElR)%mu%mu 142 143 ! Print out the coordinates of the ramp placement 144 call print_ts_voltage() 145 146 contains 147 148 subroutine get_elec_indices(Elecs, cell, na_u, xa, iElL, iElR) 149 use m_ts_electype, only : Elec_frac 150 151 ! *********************** 152 ! * INPUT variables * 153 ! *********************** 154 type(Elec), intent(in) :: Elecs(2) 155 real(dp), intent(in) :: cell(3,3) 156 integer, intent(in) :: na_u 157 real(dp), intent(in) :: xa(3,na_u) 158 159 ! *********************** 160 ! * OUTPUT variables * 161 ! *********************** 162 integer, intent(out) :: iElL, iElR 163 164 ! *********************** 165 ! * LOCAL variables * 166 ! *********************** 167 integer :: i 168 real(dp) :: r1, r2 169 170 call Elec_frac(Elecs(1), cell, na_u, xa, ts_tidx, fmin = r1) 171 call Elec_frac(Elecs(2), cell, na_u, xa, ts_tidx, fmin = r2) 172 173 if ( r1 < r2 ) then 174 iElL = 1 175 iElR = 2 176 else 177 iElL = 2 178 iElR = 1 179 end if 180 181 end subroutine get_elec_indices 182 183 end subroutine ts_init_voltage 184 185 subroutine ts_voltage(cell, nmesh, nmeshl, Vscf) 186 use precision, only : grid_p 187 use m_ts_options, only : Hartree_fname 188 ! *********************** 189 ! * INPUT variables * 190 ! *********************** 191 real(dp), intent(in) :: cell(3,3) 192 integer, intent(in) :: nmesh(3), nmeshl(3) 193 ! *********************** 194 ! * OUTPUT variables * 195 ! *********************** 196 real(grid_p), intent(inout) :: Vscf(:) 197 198 call timer('ts_volt',1) 199 200 ! Voltage drop in between the electrodes 201 ! The indices for the full cell is set 202 ! correctly to not have two routines doing the 203 ! same 204 if ( ts_tidx > 0 ) then 205 call ts_ramp(cell, nmesh, nmeshl, Vscf) 206#ifdef NCDF_4 207 else if ( len_trim(Hartree_fname) > 0 ) then 208 call ts_ncdf_Voltage(cell, Hartree_fname, 'V', nmesh, nmeshl, Vscf) 209#endif 210 else 211 call ts_elec_only(nmesh, nmeshl, Vscf) 212 end if 213 214 call timer('ts_volt',2) 215 216 end subroutine ts_voltage 217 218 ! Apply a ramp to the electrostatic potential. 219 subroutine ts_ramp(cell, nmesh, nmeshl, Vscf) 220 use precision, only : grid_p 221 use m_ts_options, only : Volt 222 use m_mesh_node, only : offset_i 223 ! *********************** 224 ! * INPUT variables * 225 ! *********************** 226 real(dp), intent(in) :: cell(3,3) 227 integer, intent(in) :: nmesh(3) 228 integer, intent(in) :: nmeshl(3) 229 ! *********************** 230 ! * OUTPUT variables * 231 ! *********************** 232 real(grid_p), intent(inout) :: Vscf(:) 233 234 ! *********************** 235 ! * LOCAL variables * 236 ! *********************** 237 integer :: i1, i2, i3, idT, imesh 238 real(dp) :: dF, dV 239 240 ! field in [0;end]: v = e*x = f*index 241 dF = - (V_high - V_low) / real(nmesh(ts_tidx) - 1,dp) 242 243 ! Find quantities in mesh coordinates 244 if ( product(nmeshl) /= size(Vscf) ) & 245 call die('ramp_elec: Vscf size not correct') 246 247 ! Add the electric field potential to the input potential 248 imesh = 0 249 if ( ts_tidx == 1 ) then 250 251 do i3 = 1 , nmeshl(3) 252 do i2 = 1 , nmeshl(2) 253 do i1 = offset_i(1), offset_i(1)+nmeshl(1)-1 254 255 dV = V_high + dF*i1 256 257 imesh = imesh + 1 258 Vscf(imesh) = Vscf(imesh) + dV 259 260 end do 261 end do 262 end do 263 264 else if ( ts_tidx == 2 ) then 265 266 do i3 = 1,nmeshl(3) 267 do i2 = offset_i(2),offset_i(2)+nmeshl(2)-1 268 269 dV = V_high + dF*i2 270 271 do i1 = 1 , nmeshl(1) 272 imesh = imesh + 1 273 Vscf(imesh) = Vscf(imesh) + dV 274 end do 275 276 end do 277 end do 278 279 else 280 281 do i3 = offset_i(3),offset_i(3)+nmeshl(3)-1 282 283 dV = V_high + dF*i3 284 285 do i2 = 1 , nmeshl(2) 286 do i1 = 1 , nmeshl(1) 287 imesh = imesh + 1 288 Vscf(imesh) = Vscf(imesh) + dV 289 end do 290 end do 291 292 end do 293 294 end if 295 296 end subroutine ts_ramp 297 298 299 subroutine ts_elec_only(nmesh, nmeshl, Vscf) 300 use precision, only : grid_p 301#ifdef MPI 302 use mpi_siesta 303#endif 304 305 use m_ts_electype 306 use m_ts_options, only : N_Elec, Elecs 307 use m_mesh_node, only : dL 308 use m_mesh_node, only : dMesh, offset_r, offset_i, mesh_correct_idx 309 use m_geom_box, only : voxel_in_box 310 311 ! *********************** 312 ! * INPUT variables * 313 ! *********************** 314 integer, intent(in) :: nmesh(3), nmeshl(3) 315 ! *********************** 316 ! * OUTPUT variables * 317 ! *********************** 318 real(grid_p), intent(inout) :: Vscf(:) 319 320 ! *********************** 321 ! * LOCAL variables * 322 ! *********************** 323 integer :: i1, i2, i3, iEl, imesh 324 integer :: imin(3), imax(3), idx(3) 325 real(dp) :: ll(3) 326#ifdef TRANSIESTA_BOX 327 integer, allocatable :: n_V(:) 328 329 allocate(n_V(N_Elec)) 330 n_V = 0 331#endif 332 333 ! We do a loop in the local grid 334 do iEl = 1 , N_Elec 335 336 call Elec_box2grididx(Elecs(iEl),nmesh,dL,imin,imax) 337 338 ! Now we have the minimum index for the box encompassing 339 ! the electrode 340 341 ! Loop the indices, and figure out whether 342 ! each of them lies in the local grid 343!$OMP parallel do default(shared), private(i1,i2,i3,ll,idx,imesh) 344 do i3 = imin(3) , imax(3) 345 do i2 = imin(2) , imax(2) 346 do i1 = imin(1) , imax(1) 347 348 ! Transform the index to the unit-cell index 349 idx(1) = i1 350 idx(2) = i2 351 idx(3) = i3 352 call mesh_correct_idx(nmesh,idx) 353 idx = idx - offset_i - 1 354 355 if ( idx(1) >= 0 .and. idx(2) >= 0 .and. idx(3) >= 0 .and. & 356 idx(1) < nmeshl(1) .and. idx(2) < nmeshl(2) .and. & 357 idx(3) < nmeshl(3) ) then 358 359 ll(:) = offset_r(:) + & 360 idx(1)*dL(:,1) + idx(2)*dL(:,2) + idx(3)*dL(:,3) 361 362 if ( voxel_in_box(Elecs(iEl)%box, ll, dMesh)) then 363 imesh = idx(1) + (idx(2) + idx(3)*nmeshl(2))*nmeshl(1) + 1 ! to correct for -1 364#ifdef TRANSIESTA_BOX 365 n_V(iEl) = n_V(iEl) + 1 366#endif 367 Vscf(imesh) = Vscf(imesh) + Elecs(iEl)%mu%mu 368 end if 369 end if 370 371 end do 372 end do 373 end do 374!$OMP end parallel do 375 376 end do 377 378#ifdef TRANSIESTA_BOX 379#ifdef MPI 380 call MPI_AllReduce(MPI_In_Place,n_V,N_Elec,MPI_Integer,MPI_Sum, & 381 MPI_Comm_World,i1) 382#endif 383 print '(10(tr1,i8))',n_V,size(vscf) 384 385 if ( any(n_V == 0) ) then 386 iEl = minloc(n_V,1) 387 write(*,'(3a)') 'ts-voltage: Elec-Box ',trim(Elecs(iEl)%name), & 388 ' is not within the device grid box.' 389 ! call die('ts-voltage: Elec-Box has an electrode outside & 390 ! &the grid. Check TS output') 391 end if 392 393 deallocate(n_V) 394#endif 395 396 end subroutine ts_elec_only 397 398 399 ! Print out the voltage direction dependent on the cell parameters. 400 subroutine print_ts_voltage() 401 use parallel, only : IONode 402 use units, only : eV 403 404 if ( .not. IONode ) return 405 406 ! Print the ramp coordinates 407 write(*,'(a,2(f6.3,tr1,a),a)') & 408 'ts-voltage: Ramp ', V_high/eV, 'eV to ', V_low/eV, 'eV ', & 409 'placed in cell' 410 411 end subroutine print_ts_voltage 412 413#ifdef NCDF_4 414 ! Read in a potential file from a NetCDF file. 415 ! Thus the potential landscape can be fully customized by the user. 416 ! We note that the potential landscape need only be calculated 417 ! for one V, direct interpolation is possible as 418 ! the solution to the Poisson equation is linearly dependent on the BC 419 subroutine ts_ncdf_voltage(cell, fname, V_name, nmesh, nmeshl, V) 420 use precision, only: grid_p 421#ifdef MPI 422 use mpi_siesta, only : MPI_Comm_World, MPI_Bcast, MPI_Grid_Real 423#endif 424 use m_ncdf_io, only : cdf_r_grid 425 use netcdf_ncdf 426 427#ifdef TRANSIESTA_VOLTAGE_DEBUG 428 use iogrid_netcdf, only: write_grid_netcdf 429#endif 430 431 real(dp), intent(in) :: cell(3,3) 432 character(len=*), intent(in) :: fname, V_name 433 ! global and local number of mesh-divisions 434 integer, intent(in) :: nmesh(3), nmeshl(3) 435 real(grid_p), intent(inout) :: V(:) 436 437 type(hNCDF) :: ncdf 438 real(grid_p) :: Vmm(2), fact 439 real(grid_p), allocatable :: tmpV(:) 440#ifdef MPI 441 integer :: MPIerror 442#endif 443 444 ! All nodes should allocate an auxilliary grid 445 allocate(tmpV(product(nmeshl))) 446 if ( size(V) /= size(tmpV) ) & 447 call die('ncdf_voltage: Vscf size not correct') 448 449 ! Open the file 450 call ncdf_open(ncdf,fname, mode=NF90_NOWRITE) 451 452 ! Read in the grid (should be in Ry) 453 call cdf_r_grid(ncdf,trim(V_name),nmeshl,tmpV) 454 455 ! retrieve the max min from the file (should be in Ry) 456 call ncdf_get_var(ncdf,trim(V_name)//'min',Vmm(1)) 457 call ncdf_get_var(ncdf,trim(V_name)//'max',Vmm(2)) 458#ifdef MPI 459 call MPI_Bcast(Vmm,2,MPI_Grid_Real,0,MPI_Comm_World,MPIerror) 460#endif 461 462 call ncdf_close(ncdf) 463 464 ! Correct the limits so that we align to the current potential 465 fact = ( V_high - V_low ) / ( Vmm(2) - Vmm(1) ) 466 ! Align the bottom potentials so that the range becomes correct 467 Vmm(1) = V_low - Vmm(1) * fact 468 469#ifdef TRANSIESTA_VOLTAGE_DEBUG 470 tmpV(:) = tmpV(:) * fact + Vmm(1) 471 call write_grid_netcdf( cell, nmesh, 1, product(nmeshl), tmpV, & 472 "TransiestaHartreePotential") 473 call bye('transiesta debug for Hartree potential') 474#endif 475 476 V(:) = Vmm(1) + V(:) + tmpV(:) * fact 477 478 deallocate(tmpV) 479 480 end subroutine ts_ncdf_voltage 481 482 483 subroutine ts_ncdf_voltage_assert(fname, cell, nmesh) 484 485 use parallel, only : IONode 486 487 use netcdf_ncdf 488 use dictionary 489 490 character(len=*), intent(in) :: fname 491 real(dp), intent(in) :: cell(3,3) 492 integer, intent(in) :: nmesh(3) 493 494 type(hNCDF) :: ncdf 495 type(dictionary_t) :: dic 496 character(len=3), parameter :: XYZ = 'ABC' 497 integer :: lnmesh(3), i 498 logical :: found 499 500 if ( .not. IONode ) return 501 502 ! Ensure that it will fail if not found 503 lnmesh = -1 504 call ncdf_open(ncdf, trim(fname), mode=NF90_NOWRITE) 505 506 ! Get the grid-size, here we just fetch all dimensions 507 ! and find them afterwards 508 call ncdf_inq(ncdf,dict_dim=dic) 509 510 ! We allow the names to be n1/x/a 511 if ( 'a' .in. dic ) then 512 call assign(lnmesh(1),dic,'a') 513 else if ( 'n1' .in. dic ) then 514 call assign(lnmesh(1),dic,'n1') 515 else if ( 'x' .in. dic ) then 516 call assign(lnmesh(1),dic,'x') 517 end if 518 519 ! We allow the names to be n2/y/b 520 if ( 'b' .in. dic ) then 521 call assign(lnmesh(2),dic,'b') 522 else if ( 'n2' .in. dic ) then 523 call assign(lnmesh(2),dic,'n2') 524 else if ( 'y' .in. dic ) then 525 call assign(lnmesh(2),dic,'y') 526 end if 527 528 ! We allow the names to be n3/z/c 529 if ( 'c' .in. dic ) then 530 call assign(lnmesh(3),dic,'c') 531 else if ( 'n3' .in. dic ) then 532 call assign(lnmesh(3),dic,'n3') 533 else if ( 'z' .in. dic ) then 534 call assign(lnmesh(3),dic,'z') 535 end if 536 537 ! Clean up 538 call ncdf_close(ncdf) 539 call delete(dic) 540 541 ! Check variables 542 if ( any(lnmesh /= nmesh) ) then 543 544 write(*,*) 545 write(*,'(a)') 'TS.Poisson file cannot be used!' 546 write(*,'(a)') 'Please carefully read the below error message:' 547 548 write(*,'(/,a)') 'TranSiesta internal grid dimensions are:' 549 do i = 1, 3 550 write(*,'(a,i0)') ' '//trim(XYZ(i:i))//': ', nmesh(i) 551 end do 552 553 write(*,'(/,3a)') 'File: ', trim(fname), ' grid dimensions are:' 554 found = .true. 555 do i = 1, 3 556 write(*,'(a,i0)') ' '//trim(XYZ(i:i))//': ', lnmesh(i) 557 found = found .and. lnmesh(i) > 0 558 end do 559 560 write(*,*) 561 if ( found ) then 562 write(*,'(a)') 'The grid dimensions *MUST* be equivalent.' 563 else 564 write(*,'(3a)') 'The lattice dimensions in ', trim(fname), ' *MUST* be named:' 565 write(*,'(a)') ' 1st lattice vector: a|n1|x' 566 write(*,'(a)') ' 2nd lattice vector: b|n2|y' 567 write(*,'(a)') ' 3rd lattice vector: c|n3|z' 568 end if 569 write(*,*) 570 571 call die('Incorrect grid in NetCDF file for user supplied & 572 &Poisson solution.') 573 574 else 575 write(*,'(a)') 'ts-voltage: User Poisson solution asserted...' 576 end if 577 578 end subroutine ts_ncdf_voltage_assert 579 580#endif 581 582end module m_ts_voltage 583