1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8#:include "common.fypp" 9 10!> Interface to libPoisson routines 11module poisson_init 12 use dftbp_accuracy, only : dp 13 use dftbp_constants, only : pi 14 use dftbp_commontypes, only : TOrbitals 15 use dftbp_globalenv, only : stdOut 16 use dftbp_message 17#:if WITH_MPI 18 use libmpifx_module 19#:endif 20 use poisson 21 use libnegf_vars, only : TTransPar 22 use system_calls, only: create_directory 23 implicit none 24 private 25 26 public :: poiss_init 27 public :: poiss_updcharges, poiss_getshift 28 public :: poiss_destroy 29 public :: poiss_updcoords 30 public :: poiss_savepotential 31 public :: TPoissonInfo 32 public :: TPoissonStructure 33 34 35 !> Geometry of the atoms for the Poisson solver 36 type TPoissonStructure 37 38 !> number of atoms in central cell 39 integer :: nAtom 40 41 !> number of species 42 integer :: nSpecies 43 44 !> type of the atoms (nAtom) 45 integer, pointer :: specie0(:) 46 47 !> atom START pos for squared H/S 48 integer, pointer :: iatomstart(:) 49 50 !> coordinates in central cell 51 real(dp), pointer :: x0(:,:) 52 53 !> total number of electrons 54 real(dp) :: nel 55 56 !> lattice vectors 57 real(dp) :: latVecs(3,3) 58 59 !> electron temperature 60 real(dp) :: tempElec 61 62 !> tells whether the system is periodic 63 logical :: isperiodic 64 65 end type TPoissonStructure 66 67 68 !> Information for the Poisson solver 69 type TPoissonInfo 70 71 !> verbosity level of the library 72 integer :: verbose 73 74 !> solve or not Poisson 75 logical :: defined = .false. 76 77 !> Poisson box 78 real(dp) :: poissBox(3) 79 80 !> Minimal grid spacing 81 real(dp) :: poissGrid(3) 82 83 !> Has a containing box been identified for the structure? (.false. for periodic systems!) 84 logical :: foundBox 85 86 !> Maximum radius of atom charge density 87 real(dp) :: maxRadAtomDens 88 89 !> Required solution accuracy 90 real(dp) :: poissAcc 91 92 !> use bulk potential as boundary condition 93 logical :: bulkBC 94 95 !> read bulk potential from file 96 logical :: readBulkPot 97 98 !> activates local boundary conditions mode (C|S) 99 character(1) :: localBCType 100 101 !> buffer spacing of local boundary condition 102 real(dp) :: bufferLocBC 103 104 !> forced boundary conditions in each direction 105 integer :: overrideBC(6) 106 107 !> forced boundary conditions on bulk potential 108 integer :: overrBulkBC(6) 109 110 !> save the potential on a file 111 logical :: savePotential = .false. 112 113 !> Recompute poisson after D.M. 114 logical :: solveTwice = .false. 115 116 !> maximum number of poisson iter 117 integer :: maxPoissIter 118 119 !> Planar or cylindrical gate 120 character(1) :: gateType 121 122 !> gate direction 123 integer :: gatedir 124 125 !> Gate potential 126 real(dp) :: gatePot 127 128 !> Gate length along transport direction 129 real(dp) :: gateLength_l 130 131 !> Gate length in transverse direction 132 real(dp) :: gateLength_t 133 134 !> Insulator length 135 real(dp) :: insLength 136 137 !> Radius of insulator 138 real(dp) :: insRad 139 140 !> Radius of gate 141 real(dp) :: gateRad 142 143 !> Insulator relative dielect. 144 real(dp) :: eps_r 145 146 !> Buffer layer between dielectric and vacuum 147 real(dp) :: dr_eps 148 149 !> Box buffer inside the contact region 150 real(dp) :: bufferBox 151 152 !> Use new numerical renormalization volume (preserves total charge) 153 logical :: numericNorm 154 155 !> Check whether density cut off fits into the PLs 156 logical :: cutoffcheck 157 158 !> Number of Nodes for parallel P. 159 integer :: maxNumNodes 160 161 !> scratch folder name 162 character(:), allocatable :: scratch 163 164 end type TPoissonInfo 165 166 167contains 168 169 !> Initialise gDFTB environment and variables 170#:if WITH_MPI 171 subroutine poiss_init(structure, orb, hubbU, poissoninfo, transpar, mpicomm, initinfo) 172#:else 173 subroutine poiss_init(structure, orb, hubbU, poissoninfo, transpar, initinfo) 174#:endif 175 !> initialisation choices for poisson solver 176 Type(TPoissonStructure), intent(in) :: structure 177 178 !> data structure with atomic orbital information 179 type(TOrbitals), intent(in) :: orb 180 181 !> Hubbard Us stored as (orbital, atom) 182 real(dp), intent(in) :: hubbU(:,:) 183 184 !> Solver settings 185 Type(TPoissonInfo), intent(in) :: poissoninfo 186 187 !> Transport parameters 188 Type(TTransPar), intent(in) :: transpar 189#:if WITH_MPI 190 !> MPI details 191 Type(mpifx_comm), intent(in) :: mpicomm 192#:endif 193 !> Success of initialisation 194 logical, intent(out) :: initinfo 195 196 ! local variables 197 integer :: i, iErr 198 199 iErr = 0 200 initinfo = .true. 201 202 #:if WITH_MPI 203 call poiss_mpi_init(mpicomm) 204 call poiss_mpi_split(min(poissoninfo%maxNumNodes, mpicomm%size)) 205 call mpi_barrier(mpicomm, iErr) 206 !#:else 207 !call error("The Poisson solver currently requires MPI parallelism to be enabled") 208 #:endif 209 210 write(stdOut,*) 211 write(stdOut,*) 'Poisson Initialisation:' 212 write(stdOut,'(a,i0,a)') ' Poisson parallelized on ',numprocs,' node(s)' 213 write(stdOut,*) 214 215 ! notify solver of standard out unit 216 call set_stdout(stdOut) 217 218 ! Directory for temporary files 219 call set_scratch(poissoninfo%scratch) 220 221 if (id0) then 222 ! only use a scratch folder on the master node 223 call create_directory(trim(scratchfolder),iErr) 224 end if 225 226 if (active_id) then 227 ! processors over which the right hand side of the Poisson equation is parallelised 228 229 call init_structure(structure%nAtom, structure%nSpecies, structure%specie0, structure%x0,& 230 & structure%latVecs, structure%isperiodic) 231 232 call init_skdata(orb%nShell, orb%angShell, hubbU, iErr) 233 234 if (iErr.ne.0) then 235 call error("I am sorry... cannot proceed. orbital shells should be in the order s,p,d") 236 endif 237 238 call init_charges() 239 240 ! Initialise renormalization factors for grid projection 241 242 if (iErr.ne.0) then 243 call poiss_destroy() 244 initinfo = .false. 245 return 246 endif 247 248 !init default values for parameters 249 call init_defaults() 250 251 call set_verbose(poissonInfo%verbose) 252 253 call set_temperature(0.0_dp) 254 255 !----------------------------------------------------------------------------- 256 ! GP: verify if the calculation is on an open system (cluster=false) or not 257 ! 258 ! note: in previous versions only ncont was checked but this is not 259 ! sufficient as it does not take into account a contact calculation 260 ! with poisson solver, where the transport block is defined 261 !----------------------------------------------------------------------------- 262 if (transpar%defined .and. transpar%taskUpload) then 263 !init the number of contacts as in dftb+ 264 call set_ncont(transpar%ncont) 265 else 266 call set_ncont(0) 267 endif 268 call set_cluster(.false.) 269 if (.not.transpar%defined) then 270 call set_cluster(.true.) 271 endif 272 if (transpar%defined .and. (.not. transpar%taskUpload)) then 273 call set_cluster(.true.) 274 endif 275 276 !-----------------------------------------------------------------------------+ 277 ! TRANSPORT PARAMETER NEEDED FOR POISSON (contact partitioning) 278 !-----------------------------------------------------------------------------+ 279 call set_mol_indeces(transpar%idxdevice(1:2), structure%natom) 280 call set_cont_indeces(transpar%contacts(1:ncont)%idxrange(1), 1) 281 call set_cont_indeces(transpar%contacts(1:ncont)%idxrange(2), 2) 282 call set_contdir(transpar%contacts(1:ncont)%dir) 283 call set_fermi(transpar%contacts(1:ncont)%eFermi(1)) 284 call set_potentials(transpar%contacts(1:ncont)%potential) 285 call set_builtin() 286 287 call set_dopoisson(poissoninfo%defined) 288 call set_poissonbox(poissoninfo%poissBox) 289 call set_poissongrid(poissoninfo%poissGrid) 290 call set_accuracy(poissoninfo%poissAcc) 291 292 FoundBox = poissoninfo%foundBox 293 InitPot = poissoninfo%bulkBC 294 ReadBulk = poissoninfo%readBulkPot 295 MaxPoissIter = poissoninfo%maxPoissIter 296 select case (poissoninfo%localBCType) 297 case('G') 298 localBC = 0 299 case('C') 300 localBC = 1 301 case('S') 302 localBC = 2 303 end select 304 deltaR_max = poissoninfo%maxRadAtomDens 305 306 ! if deltaR_max > 0 is a radius cutoff, if < 0 a tolerance 307 if (deltaR_max < 0.0_dp) then 308 write(stdOut,*) "Atomic density tolerance: ", -deltaR_max 309 deltaR_max = getAtomDensityCutoff(-deltaR_max, uhubb) 310 end if 311 312 write(stdOut,*) "Atomic density cutoff: ", deltaR_max, "a.u." 313 314 if (ncont /= 0 .and. poissoninfo%cutoffcheck) then 315 call checkDensityCutoff(deltaR_max, transpar%contacts(:)%length) 316 end if 317 318 dR_cont = poissoninfo%bufferLocBC 319 bufferBox = poissoninfo%bufferBox 320 SavePot = poissoninfo%savePotential 321 overrideBC = poissoninfo%overrideBC 322 overrBulkBC = poissoninfo%overrBulkBC 323 !-----------------------------------------------------------------------------+ 324 ! Gate settings 325 DoGate=.false. 326 DoCilGate=.false. 327 if (poissoninfo%gateType.eq.'C') then 328 DoCilGate = .true. 329 end if 330 if (poissoninfo%gateType.eq.'P') then 331 DoGate = .true. 332 end if 333 334 Gate = poissoninfo%gatePot 335 GateLength_l = poissoninfo%gateLength_l 336 GateLength_t = poissoninfo%gateLength_t 337 338 ! Planar gate must be along y 339 GateDir = poissoninfo%gatedir 340 341 OxLength = poissoninfo%insLength 342 Rmin_Gate = poissoninfo%gateRad 343 Rmin_Ins = poissoninfo%insRad 344 eps_r = poissoninfo%eps_r 345 dr_eps = poissoninfo%dr_eps 346 347 ! Use fixed analytical renormalization (approximate) or numerical 348 fixed_renorm = .not.(poissoninfo%numericNorm) 349 350 ! Performs parameters checks 351 call check_biasdir() 352 call check_poisson_box() 353 if (any(overrideBC.ne.0)) then 354 period = .false. 355 end if 356 call check_parameters() 357 call check_localbc() 358 call write_parameters() 359 call check_contacts() 360 361 !-----------------------------------------------------------------------------+ 362 ! Parameters from old PAR.in not yet parsed: 363 ! 364 ! PoissPlane 365 ! DoTip,tip_atom,base_atom1,base_atom2 366 !-----------------------------------------------------------------------------+ 367 368 write(stdOut,'(79(">"))') 369 370 endif 371 372 end subroutine poiss_init 373 374 375 !> Release gDFTB varibles in poisson library 376 subroutine poiss_destroy() 377 378 if (active_id) then 379 write(stdOut,'(A)') 380 write(stdOut,'(A)') 'Release Poisson Memory:' 381 call poiss_freepoisson() 382 endif 383 384 end subroutine poiss_destroy 385 386 387 !> Interface subroutine to call Poisson 388 subroutine poiss_getshift(V_L_atm,grad_V) 389 390 !> potential at atom sites 391 real(dp), intent(inout) :: V_L_atm(:,:) 392 393 !> Gradient of potential 394 real(dp), intent(out), optional :: grad_V(:,:) 395 396 real(dp) :: fakegrad(3,1) 397 398 integer :: ndim, ierr, array_size 399 integer :: PoissFlag 400 401 ! The subroutine gives the atom shifts. 402 403 ! Poiss Flag is a control flag: 404 ! 0 - potential in SCC 405 ! 1 - atomic shift component of gradient 406 407 ! This information is needed by all nodes 408 PoissFlag=0 409 if(present(grad_V)) then 410 PoissFlag=1 411 end if 412 413 if (active_id) then 414 415 select case(PoissFlag) 416 case(0) 417 if (verbose.gt.30) then 418 write(stdOut,*) 419 write(stdOut,'(80("="))') 420 write(stdOut,*) ' SOLVING POISSON EQUATION ' 421 write(stdOut,'(80("="))') 422 end if 423 call init_PoissBox(iErr) 424 if (iErr /= 0) then 425 call error("Failure during initialisation of the Poisson box") 426 end if 427 call mudpack_drv(PoissFlag,V_L_atm,fakegrad) 428 case(1) 429 call mudpack_drv(PoissFlag,V_L_atm,grad_V) 430 end select 431 432 if (verbose.gt.30) then 433 write(stdOut,'(80("*"))') 434 end if 435 end if 436 437 #:if WITH_MPI 438 call mpifx_barrier(global_comm, ierr) 439 select case(PoissFlag) 440 ! Note: V_L_atm and grad_V are allocated for all processes in dftb+.F90 441 case(0) 442 call mpifx_bcast(global_comm, V_L_atm) 443 case(1) 444 call mpifx_bcast(global_comm, V_L_atm) 445 call mpifx_bcast(global_comm, grad_V) 446 end select 447 call mpifx_barrier(global_comm) 448 #:endif 449 450 end subroutine poiss_getshift 451 452 453 !> Interface subroutine to overload Mulliken charges stored in libPoisson 454 subroutine poiss_updcharges(q,q0) 455 456 !> populations 457 real(dp), intent(in) :: q(:,:) 458 459 !> reference charges 460 real(dp), intent(in) :: q0(:,:) 461 462 integer :: nsh, l, i, o, orb 463 real(dp) :: Qtmp 464 465 if (active_id) then 466 467 if (size(q, dim=2).ne.natoms) then 468 call error('ERROR in udpcharges size of q') 469 end if 470 471 do i = 1, natoms 472 nsh = lmax(izp(i))+1 473 orb=0 474 do l = 1, nsh 475 Qtmp = 0.0_dp 476 do o= 1,2*l-1 477 orb = orb + 1 478 ! - is for negative electron charge density 479 Qtmp = Qtmp + q0(orb,i)-q(orb,i) 480 enddo 481 dQmat(l,i) = Qtmp 482 enddo 483 enddo 484 485 endif 486 487 end subroutine poiss_updcharges 488 489 490 !> Calculates the atom density cutoff from the density tolerance. 491 function getAtomDensityCutoff(denstol, uhubb) result(res) 492 493 !> Density tolerance. 494 real(dp), intent(in) :: denstol 495 496 !> List of atomic Hubbard U values. 497 real(dp), intent(in) :: uhubb(:,:) 498 499 !> Maximal atom density cutoff. 500 real(dp) :: res 501 502 integer :: typ, nsh 503 real(dp) :: tau 504 505 res = 0.0_dp 506 do typ = 1, size(uhubb,2) 507 nsh = lmax(typ)+1 508 tau = 3.2_dp * minval(uhubb(1:nsh,typ)) 509 res = max(res, -log(8.0_dp * pi / tau**3 * denstol) / tau) 510 end do 511 512 end function getAtomDensityCutoff 513 514 515 !> Checks whether density cutoff fits into the PLs and stop if not. 516 subroutine checkDensityCutoff(rr, pllens) 517 518 !> Density cutoff. 519 real(dp), intent(in) :: rr 520 521 !> Lengths of the principal layers in the contacts. 522 real(dp), intent(in) :: pllens(:) 523 524 integer :: ii 525 526 ! GP In Poisson both the contact layers are used, which is the reason why we 527 ! have a factor 2 in front of pllens 528 do ii = 1, size(pllens) 529 if (rr > 2.0_dp * pllens(ii) + 1e-12_dp) then 530 write(stdOut,"(A,I0,A)") "!!! ERROR: Atomic density cutoff incompatible with the principle& 531 & layer width in contact ", ii, "." 532 write(stdOut,"(A,G10.3,A,G10.3,A)") " (", rr, ">", pllens(ii), ")" 533 call error("Either enlarge PL width in the contact or increase AtomDensityCutoff or& 534 & AtomDensityTolerance.") 535 end if 536 end do 537 538 end subroutine checkDensityCutoff 539 540 541end module poisson_init 542