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