Line data Source code
1 : !!****m* ABINIT/m_psolver
2 : !! NAME
3 : !! m_psolver
4 : !!
5 : !! FUNCTION
6 : !! Poisson solver
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR,TRangel).
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_psolver
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_abicore
27 : use m_errors
28 : use m_abi2big
29 : use m_cgtools
30 : use m_xmpi
31 :
32 : use defs_abitypes, only : mpi_type
33 : use m_geometry, only : metric
34 : use m_drivexc, only : mkdenpos
35 :
36 : implicit none
37 :
38 : private
39 : !!***
40 :
41 : public :: psolver_rhohxc
42 : public :: psolver_hartree
43 : public :: psolver_kernel
44 : !!***
45 :
46 : contains
47 : !!***
48 :
49 : !!****f* ABINIT/psolver_rhohxc
50 : !! NAME
51 : !! psolver_rhohxc
52 : !!
53 : !! FUNCTION
54 : !! Given rho(r), compute Hartree potential considering the system as
55 : !! an isolated one. This potential is obtained from the convolution
56 : !! of 1/r and rho(r), treated in Fourier space. This method is a wrapper around
57 : !! Psolver() developped for BigDFT.
58 : !! It can compute the xc energy and potential if required. This computation is
59 : !! built on the drivexc() routine of ABINIT but access it directly from real
60 : !! space. The present routine is a real space counter part to rhotoxc().
61 : !!
62 : !! INPUTS
63 : !! dtset <type(dataset_type)>=all input variables in this dataset
64 : !! mpi_enreg=MPI-parallelisation information.
65 : !! rhor(nfft,nspden)=electron density in real space in electrons/bohr**3
66 : !!
67 : !! OUTPUT
68 : !! enhartr=returned Hartree energy (hartree).
69 : !! enxc=returned exchange and correlation energy (hartree).
70 : !! envxc=returned energy of the Vxc potential (hartree).
71 : !! vhartr(nfft)=Hartree potential.
72 : !! vxc(nfft,nspden)=xc potential
73 : !! vxcavg=<Vxc>=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r].
74 : !!
75 : !! NOTE
76 : !! In psolver, with nspden == 2, rhor(:,1) = density up and
77 : !! rhor(:,2) = density down.
78 : !! But in ABINIT (dtset%usewvl != 1) rhor(:,1) = total density and
79 : !! rhor(:,2) = density up .
80 : !! In ABINIT (dtset%usewvl != 1), the same convention is used as in psolver.
81 : !!
82 : !! SOURCE
83 :
84 0 : subroutine psolver_rhohxc(enhartr, enxc, envxc, icoulomb, ixc, &
85 0 : & mpi_enreg, nfft, ngfft, nhat,nhatdim,&
86 0 : & nscforder, nspden, n3xccc, rhor, rprimd,&
87 0 : & usexcnhat,usepaw,usewvl,vhartr, vxc, vxcavg, wvl,wvl_den,wvl_e,&
88 0 : & xccc3d,xclevel,xc_denpos)
89 :
90 : #if defined HAVE_BIGDFT
91 : use BigDFT_API, only : XC_potential,ELECTRONIC_DENSITY,coulomb_operator
92 : use poisson_solver, only : H_potential
93 : #endif
94 :
95 : !Arguments ------------------------------------
96 : !scalars
97 : integer, intent(in) :: nhatdim,nspden,n3xccc
98 : integer, intent(in) :: nfft, icoulomb, ixc, nscforder, usewvl
99 : integer,intent(in) :: usexcnhat,usepaw,xclevel
100 : real(dp),intent(in) :: rprimd(3,3)
101 : real(dp), intent(in) :: xc_denpos
102 : real(dp), intent(out) :: enxc, envxc, enhartr, vxcavg
103 : type(mpi_type), intent(in) :: mpi_enreg
104 : type(wvl_internal_type), intent(in) :: wvl
105 : type(wvl_denspot_type), intent(inout) :: wvl_den
106 : type(wvl_energy_terms), intent(inout) :: wvl_e
107 : !arrays
108 : integer, intent(in) :: ngfft(18)
109 : real(dp),intent(in) :: xccc3d(n3xccc)
110 : real(dp),intent(in) :: nhat(nfft,nspden*nhatdim)
111 : real(dp),intent(inout) :: rhor(nfft, nspden)
112 : real(dp),intent(out) :: vhartr(nfft)
113 : real(dp),intent(out) :: vxc(nfft, nspden)
114 :
115 : !Local variables-------------------------------
116 : #if defined HAVE_BIGDFT
117 : ! n_c and \hat{n} can be added/rested inside bigdft by passing
118 : ! them as pointers (rhocore and rhohat):
119 : logical, parameter :: add_n_c_here=.true. !Add n_c here or inside bigdft
120 : logical, parameter :: rest_hat_n_here=.true. !Rest \hat{n} here or inside bigdft
121 : !scalars
122 : integer :: me,nproc,comm
123 : integer :: ifft,ispin
124 : integer :: iwarn, opt_mkdenpos
125 : integer :: nfftot,ngrad
126 : integer :: n1i,n2i,n3d,n3i
127 : real(dp) :: tmpDown, tmpUp, tmpPot,ucvol,ucvol_local
128 : logical :: sumpion,test_nhat,use_psolver=.false.
129 : character(len=500) :: message
130 : character(len = 1) :: datacode, bndcode
131 : !arrays
132 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
133 : real(dp) :: hgrid(3)
134 : real(dp) :: vxcmean(1)
135 : real(dp), pointer :: rhocore(:,:,:,:),rhohat(:,:,:,:)
136 : real(dp), pointer :: pot_ion(:,:,:,:),rhonow(:,:)
137 : real(dp), dimension(6) :: xcstr
138 : type(coulomb_operator) :: kernel
139 : #endif
140 :
141 : ! *********************************************************************
142 :
143 : DBG_ENTER("COLL")
144 :
145 : #if defined HAVE_BIGDFT
146 :
147 : nfftot=PRODUCT(ngfft(1:3))
148 : comm=mpi_enreg%comm_fft
149 : if(usewvl==1) comm=mpi_enreg%comm_wvl
150 : me=xmpi_comm_rank(comm)
151 : nproc=xmpi_comm_size(comm)
152 :
153 : if(n3xccc>0) then
154 : if(nfft .ne. n3xccc)then
155 : write(message,'(a,a,a,2(i0,1x))')&
156 : & 'nfft and n3xccc should be equal,',ch10,&
157 : & 'however, nfft and n3xccc=',nfft,n3xccc
158 : ABI_BUG(message)
159 : end if
160 : end if
161 : if(nspden==4) then
162 : ABI_ERROR('nspden==4 not coded yet')
163 : end if
164 :
165 : if (ixc==0) then
166 : vxcavg=zero
167 : test_nhat=.false.
168 :
169 : ! No xc at all is applied (usually for testing)
170 : ABI_WARNING('Note that no xc is applied (ixc=0).')
171 :
172 : else if (ixc/=20) then
173 :
174 : ! ngrad=1 is for LDAs or LSDs, ngrad=2 is for GGAs
175 : ngrad=1;if(xclevel==2)ngrad=2
176 : ! ixc 31 to 35 are for mgga test purpose only (fake functionals based on LDA but need the gradients too)
177 : if(ixc>=31 .and. ixc<=35)ngrad=2
178 : ! Test: has a compensation density to be added/substracted (PAW) ?
179 : ! test_nhat=((nhatdim==1).and.(usexcnhat==0.or.(ngrad==2.and.nhatgrdim==1)))
180 : test_nhat=((nhatdim==1).and.(usexcnhat==0))
181 : end if
182 :
183 :
184 : !Compute different geometric tensor, as well as ucvol, from rprimd
185 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
186 :
187 : if (icoulomb == 0) then
188 : ! The kernel is built with 'P'eriodic boundary counditions.
189 : bndcode = 'P'
190 : else if (icoulomb == 1) then
191 : ! The kernel is built with 'F'ree boundary counditions.
192 : bndcode = 'F'
193 : else if (icoulomb == 2) then
194 : ! The kernel is built with 'S'urface boundary counditions.
195 : bndcode = 'S'
196 : end if
197 :
198 : !This makes the tests fail?
199 : !For NC and n_c=0, call psolver, which uses less memory:
200 : !if(usepaw==0 .and. n3xccc==0) use_psolver=.true.
201 :
202 : if(nspden > 2)then
203 : write(message, '(a,a,a,i0)' )&
204 : & 'Only non-spin-polarised or collinear spin is allowed,',ch10,&
205 : & 'while the argument nspden = ', nspden
206 : ABI_ERROR(message)
207 : end if
208 :
209 : !We do the computation.
210 : write(message, "(A,A,A,3I6)") "psolver_rhohxc(): compute potentials (Vhartree and Vxc)...", ch10, &
211 : & " | dimension:", ngfft(1:3)
212 : call wrtout(std_out, message,'COLL')
213 :
214 : if(usewvl==1) then
215 : hgrid=(/wvl_den%denspot%dpbox%hgrids(1),wvl_den%denspot%dpbox%hgrids(2),wvl_den%denspot%dpbox%hgrids(3)/)
216 : else
217 : hgrid=(/ rprimd(1,1) / ngfft(1), rprimd(2,2) / ngfft(2), rprimd(3,3) / ngfft(3) /)
218 : end if
219 :
220 : if (usewvl == 0) then
221 : ! We get the kernel.
222 : call psolver_kernel( hgrid, 2, icoulomb, me, kernel, comm, ngfft, nproc, nscforder)
223 : elseif(usewvl==1) then
224 : ! In this case, the kernel is already computed.
225 : ! We just shallow copy it.
226 : kernel = wvl_den%denspot%pkernel
227 : end if
228 :
229 : if(usewvl==1) then
230 : if(wvl_den%denspot%rhov_is .ne. ELECTRONIC_DENSITY) then
231 : message= "psolver_rhohxc: rhov should contain the electronic density"
232 : ABI_ERROR(message)
233 : end if
234 : end if
235 :
236 : if(usewvl==1) then
237 : n1i=wvl%Glr%d%n1i; n2i=wvl%Glr%d%n2i; n3i=wvl%Glr%d%n3i
238 : n3d=wvl_den%denspot%dpbox%n3d
239 : else
240 : n1i=ngfft(1); n2i=ngfft(2) ; n3i=ngfft(3)
241 : n3d=ngfft(13)
242 : end if
243 :
244 : if (usewvl == 0) then
245 : ! ucvol_local=product(hgrid)*half**3*real(n1i*n2i*n3i,dp)
246 : ! write(*,*)'hgrid, n1i,n2i,n3i',hgrid,ngfft(1:3)
247 : ! write(*,*)'ucvol_local',ucvol_local
248 : ucvol_local = ucvol
249 : ! write(*,*)'ucvol_local',ucvol_local
250 : else
251 : ! We need to tune the volume when wavelets are used because, not
252 : ! all FFT points are used.
253 : ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
254 : ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(product(wvl_den%denspot%dpbox%ndims), dp)
255 : end if
256 :
257 : !Core density array
258 : if(n3xccc==0 .or. add_n_c_here) then
259 : nullify(rhocore)
260 : ! Pending, next line should follow the same logic that the rest
261 : if(usewvl==1 .and. usepaw==0) rhocore=> wvl_den%denspot%rho_C
262 : else
263 : if(usepaw==1) then
264 : ABI_MALLOC(rhocore,(n1i,n2i,n3d,1)) !not spin dependent
265 : call wvl_rhov_abi2big(1,xccc3d,rhocore)
266 :
267 : ! Make rhocore positive to avoid numerical instabilities in V_xc
268 : iwarn=0 ; opt_mkdenpos=0
269 : call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhocore, tol20 )
270 : end if
271 : end if
272 :
273 : !write(*,*)'psolver_rhohxc, erase me, set rhocore=0'
274 : !if( associated(wvl_den%denspot%rho_C))wvl_den%denspot%rho_C=zero
275 : !if(associated(rhocore))rhocore=zero
276 :
277 : !Rhohat array:
278 : if(test_nhat .and. .not. rest_hat_n_here) then
279 : ! rhohat => nhat !do not know how to point 4 index to 2 index
280 : ! here we have to copy since convention for spin changes.
281 : ABI_MALLOC(rhohat,(n1i,n2i,n3d,nspden))
282 : call wvl_rhov_abi2big(1,nhat,rhohat)
283 : else
284 : nullify(rhohat)
285 : end if
286 :
287 : !Data are always distributed when using the wavelets, even if nproc = 1.
288 : !The size given is the complete size of the box, not the distributed size
289 : !stored in ngfft.
290 : if (nproc > 1 .or. usewvl > 0) then
291 : datacode = 'D'
292 : else
293 : datacode = 'G'
294 : end if
295 :
296 : !If usewvl=1, vpsp(or v_ext) will be summed to vhartree
297 : if(usewvl==1) then
298 : pot_ion=>wvl_den%denspot%V_ext
299 : sumpion=.false.
300 : ! Note:
301 : ! if sumpion==.true.
302 : ! call wvl_newvtr in setvtr and rhotov
303 : ! if sumpion==.false.
304 : ! modify setvtr and rhotov to not use wvl_newvtr and follow the normal ABINIT flow.
305 : else
306 : ! This is not allowed
307 : ! pot_ion=>vxc !this is a dummy variable here
308 : sumpion=.false.
309 : end if
310 :
311 :
312 : !To make this work, make sure that xc_init has been called
313 : !in gstate.
314 : if(.not. use_psolver) then
315 : ! T.Rangel:
316 : ! Use this approach for PAW and sometimes for NC since
317 : ! in psolver() the core density is not added.
318 : !
319 : ! PAW case:
320 : ! It is important to call H_potential before XC_potential:
321 : ! In XC_potential, if test_nhat, we do:
322 : ! 1) rhor=rhor-rhohat,
323 : ! 2) makepositive(rhor,tol20)
324 : ! 3) after Psolver, we do rhor=rhor+rhohat,
325 : ! I found that rhor at input and output are slightly different,
326 : ! These differences lead to a difference of ~0.01 hartree in V_hartree.
327 : ! If PAW, substract compensation density from effective density:
328 : ! - if GGA, because nhat gradients are computed separately
329 : ! - if nhat does not have to be included in XC
330 :
331 : ! save rhor in rhonow to avoid modifying it.
332 : ABI_MALLOC(rhonow,(nfft,nspden))
333 : ! copy rhor into rhonow:
334 : ! ABINIT convention is followed: (ispin=1: for spin up + spin down)
335 : rhonow(1:nfft,1:nspden)=abs(rhor(1:nfft,1:nspden))+1.0d-20
336 :
337 : if(usewvl==1) then
338 : call H_potential(datacode,&
339 : & kernel,rhonow,pot_ion,enhartr,&
340 : & zero,sumpion)
341 : else
342 : ! Vxc is passed as a dummy argument
343 : call H_potential(datacode,&
344 : & kernel,rhonow,vxc,enhartr,&
345 : & zero,sumpion)
346 : end if
347 : !
348 : vhartr(1:nfft)=rhonow(1:nfft,1)
349 : ! write(*,*)'erase me psolver_rhohxc l350, set vhartr=0'
350 : ! vhartr=zero ; enhartr=zero
351 : !
352 : ! Since rhonow was modified inside H_potential:
353 : ! copy rhor again into rhonow following the BigDFT convention:
354 : call wvl_rhov_abi2big(1,rhor,rhonow)
355 :
356 : ! Add n_c here:
357 : if(n3xccc>0 .and. add_n_c_here) then
358 : do ispin=1,nspden
359 : rhonow(:,ispin)=rhonow(:,ispin)+xccc3d(:)
360 : end do
361 : end if
362 : ! Remove \hat{n} here:
363 : if(test_nhat .and. rest_hat_n_here) then
364 : rhonow(1:nfft,1:nspden)=rhonow(1:nfft,1:nspden)-nhat(1:nfft,1:nspden)
365 : end if
366 :
367 : ! Make the density positive everywhere (but do not care about gradients)
368 : iwarn=0 ; opt_mkdenpos=0
369 : call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhonow, xc_denpos)
370 : ! do ispin=1,nspden
371 : ! do ifft=1,nfft
372 : ! rhonow(ifft,ispin)=abs(rhonow(ifft,ispin))+1.0d-20
373 : ! end do
374 : ! end do
375 :
376 : ! If PAW, substract compensation density from effective density:
377 : ! - if GGA, because nhat gradients are computed separately
378 : ! - if nhat does not have to be included in XC
379 : if (test_nhat .and. .not. rest_hat_n_here) then
380 :
381 : call XC_potential(bndcode,datacode,me,nproc,comm,&
382 : & n1i,n2i,n3i,&
383 : & wvl_den%denspot%xc,hgrid(1),hgrid(2),hgrid(3),&
384 : & rhonow,enxc,envxc,nspden,rhocore,&
385 : & vxc,xcstr,rhohat=rhohat)
386 :
387 : else
388 :
389 : call XC_potential(bndcode,datacode,me,nproc,comm,&
390 : & n1i,n2i,n3i,&
391 : & wvl_den%denspot%xc,hgrid(1),hgrid(2),hgrid(3),&
392 : & rhonow,enxc,envxc,nspden,rhocore,&
393 : & vxc,xcstr)
394 :
395 : end if
396 :
397 : ! write(*,*)'psolver_rhohxc: erase me, set vxc=0'
398 : ! vxc=zero
399 : ! enxc=zero
400 : ! envxc=zero
401 :
402 : ! deallocate temporary array
403 : ABI_FREE(rhonow)
404 :
405 : else
406 : ! NC case: here we optimize memory, and we reuse vhartree to store rhor:
407 :
408 : ! We save total rhor in vhartr
409 : vhartr(1:nfft) = rhor(1:nfft, 1)
410 :
411 : ! In non-wavelet case, we change the rhor values.
412 : if (nspden == 2) then
413 : do ifft = 1, nfft
414 : ! We change rhor for psolver call.
415 : tmpDown = rhor(ifft, 1) - rhonow(ifft, 2)
416 : tmpUp = rhor(ifft, 2)
417 : rhor(ifft, 1) = tmpUp
418 : rhor(ifft, 2) = tmpDown
419 : end do
420 : end if
421 : ! Make the density positive everywhere (but do not care about gradients)
422 : iwarn=0 ; opt_mkdenpos=0
423 : call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhor, xc_denpos)
424 : ! do ispin=1,nspden
425 : ! do ifft=1,nfft
426 : ! rhor(ifft,ispin)=abs(rhor(ifft,ispin))+1.0d-20
427 : ! end do
428 : ! end do
429 :
430 : ! Call Poisson solver, here rhor(:,1) will contain Vhartree at output
431 : ! This does not compile, check mklocl_realspace where it do work.
432 : ! call psolver(bndcode, datacode, me, nproc, n1i, &
433 : !& n2i,n3i, ixc, hgrid(1), hgrid(2), hgrid(3), &
434 : !& rhor, kernel, vxc, enhartr, enxc, envxc, 0.d0, .false., nspden)
435 :
436 : ! PSolver work in place, we set back the rhor values.
437 : do ifft = 1, nfft, 1
438 : tmpPot = rhor(ifft, 1)
439 : ! Rhor total was saved in vhartr and current rhor(:,2) is down spin
440 : rhor(ifft, 1) = vhartr(ifft)
441 : if (nspden == 2) rhor(ifft, 2) = rhor(ifft, 1) - rhor(ifft, 2)
442 : vhartr(ifft) = tmpPot
443 : end do
444 : end if
445 :
446 : !Pass vhartr and vxc to BigDFT objects (useless?)
447 : !if(usewvl==1) then
448 : ! write(message, '(a,a,a,a)' ) ch10, ' rhotoxc_wvlpaw : but why are you copying me :..o('
449 : ! call wvl_vhartr_abi2big(1,vhartr,wvl_den)
450 : ! (this can be commented out, since we do not use denspot%v_xc
451 : ! call wvl_vxc_abi2big(1,vxc,wvl_den)
452 : !end if
453 :
454 : !Compute vxcavg
455 : call mean_fftr(vxc, vxcmean, nfft, nfftot, nspden,mpi_comm_sphgrid=comm)
456 : vxcavg = vxcmean(1)
457 :
458 : !Pass energies to wvl object:
459 : if(usewvl==1) then
460 : wvl_e%energs%eh = enhartr
461 : wvl_e%energs%exc = enxc
462 : wvl_e%energs%evxc= envxc
463 : end if
464 :
465 : !Nullify pointers and deallocate arrays
466 : if(test_nhat .and. .not. rest_hat_n_here) then
467 : ! if(nspden==2) ABI_FREE(rhohat)
468 : ABI_FREE(rhohat)
469 : if(associated(rhohat)) nullify(rhohat)
470 : end if
471 : if( n3xccc>0 .and. .not. add_n_c_here) then
472 : if(usepaw==1) then
473 : ABI_FREE(rhocore)
474 : end if
475 : end if
476 : if(associated(rhocore)) nullify(rhocore)
477 : if(associated(pot_ion)) nullify(pot_ion)
478 :
479 : #else
480 0 : BIGDFT_NOTENABLED_ERROR()
481 : if (.false.) write(std_out,*) nhatdim,nspden,n3xccc,nfft,icoulomb,ixc,nscforder,usewvl,&
482 : & usexcnhat,usepaw,xclevel,rprimd(1,1),xc_denpos,enxc,envxc,enhartr,vxcavg,mpi_enreg%nproc,&
483 : & wvl%h(1),wvl_den%symObj,wvl_e%energs,ngfft(1),xccc3d(1),nhat(1,1),rhor(1,1),vhartr(1),vxc(1,1)
484 : #endif
485 :
486 : DBG_EXIT("COLL")
487 :
488 0 : end subroutine psolver_rhohxc
489 : !!***
490 :
491 : !!****f* ABINIT/Psolver_hartree
492 : !! NAME
493 : !! Psolver_hartree
494 : !!
495 : !! FUNCTION
496 : !! Given rho(r), compute Hartree potential considering the system as
497 : !! an isolated one. This potential is obtained from the convolution
498 : !! of 1/r and rho(r), treated in Fourier space. This method is a wrapper around
499 : !! Psolver() developped for BigDFT.
500 : !! It does not compute the xc energy nor potential. See psolver_rhohxc() to do it.
501 : !! WARNING : the XC energy and potential computation capability has been
502 : !! for spin-polarized case, as everything is done as if nspden=1
503 : !!
504 : !! INPUTS
505 : !! dtset <type(dataset_type)>=all input variables in this dataset
506 : !! mpi_enreg=MPI-parallelisation information.
507 : !! rhor(nfft,nspden)=electron density in real space in electrons/bohr**3
508 : !!
509 : !! OUTPUT
510 : !! enhartr=returned Hartree energy (hartree).
511 : !! vhartr(nfft)=Hartree potential.
512 : !!
513 : !! NOTE
514 : !! In PSolver, with nspden == 2, rhor(:,1) = density up and
515 : !! rhor(:,2) = density down.
516 : !! But in ABINIT (dtset%usewvl != 1) rhor(:,1) = total density and
517 : !! rhor(:,2) = density up .
518 : !! In ABINIT (dtset%usewvl != 1), the same convention is used as in PSolver.
519 : !!
520 : !! SOURCE
521 :
522 0 : subroutine psolver_hartree(enhartr, hgrid, icoulomb, me, mpi_comm, nfft, ngfft, nproc, &
523 0 : & nscforder, nspden, rhor, vhartr, usewvl)
524 :
525 : #if defined HAVE_BIGDFT
526 : use BigDFT_API, only : coulomb_operator
527 : use poisson_solver, only : H_potential
528 : #endif
529 :
530 : !Arguments ------------------------------------
531 : !scalars
532 : integer, intent(in) :: nfft, nspden, icoulomb, usewvl, mpi_comm, me, nproc, nscforder
533 : real(dp), intent(out) :: enhartr
534 : !arrays
535 : integer, intent(in) :: ngfft(3)
536 : real(dp),intent(in) :: hgrid(3)
537 : real(dp),intent(in) :: rhor(nfft,nspden)
538 : real(dp),intent(out) :: vhartr(nfft)
539 :
540 : !Local variables-------------------------------
541 : #if defined HAVE_BIGDFT
542 : !scalars
543 : character(len=500) :: message
544 : character(len = 1) :: datacode, bndcode
545 : !arrays
546 : real(dp), dimension(1) :: pot_ion_dummy
547 : type(coulomb_operator):: kernel
548 : #endif
549 :
550 : ! *********************************************************************
551 :
552 : #if defined HAVE_BIGDFT
553 :
554 : if (icoulomb == 0) then
555 : ! The kernel is built with 'P'eriodic boundary counditions.
556 : bndcode = 'P'
557 : else if (icoulomb == 1) then
558 : ! The kernel is built with 'F'ree boundary counditions.
559 : bndcode = 'F'
560 : else if (icoulomb == 2) then
561 : ! The kernel is built with 'S'urface boundary counditions.
562 : bndcode = 'S'
563 : end if
564 :
565 : if(nspden > 2 .and. usewvl/=0 )then
566 : write(message, '(a,a,a,i0)' )&
567 : & 'Only non-spin-polarised or collinear spin is allowed for wavelets,',ch10,&
568 : & 'while the argument nspden = ', nspden
569 : ABI_BUG(message)
570 : end if
571 :
572 : !We do the computation.
573 : write(message, "(A,A,A,3I6)") "Psolver_hartree(): compute potential (Vhartree)...", ch10, &
574 : & " | dimension:", ngfft(1:3)
575 : call wrtout(std_out, message,'COLL')
576 :
577 : if (usewvl == 0) then
578 : vhartr(:) = rhor(:, 1)
579 :
580 : datacode = 'G'
581 : ! This may not work with MPI in the planewave code...
582 : else
583 : if(nspden==1)vhartr(:) = rhor(:, 1)
584 : if(nspden==2)vhartr(:) = rhor(:, 1) + rhor(:, 2)
585 : ! The data are 'D'istributed in the wavelet case or 'G'lobal otherwise.
586 : if (nproc > 1) then
587 : datacode = 'D'
588 : else
589 : datacode = 'G'
590 : end if
591 : end if
592 :
593 : !We get the kernel.
594 : call psolver_kernel( hgrid, 2, icoulomb, me, kernel, mpi_comm, ngfft, nproc, nscforder)
595 :
596 :
597 : !We attack PSolver with the total density contained in vhartr.
598 : !This is also valid for spin-polarized (collinear and non-collinear)
599 : !systems. Thus we enter nspden (last arg of PSolver) as being 1.
600 : !Warning : enxc and evxc are meaningless.
601 : ! call psolver(bndcode, datacode, me, nproc, ngfft(1), ngfft(2), ngfft(3),&
602 : !& 0, hgrid(1), hgrid(2), hgrid(3), vhartr, kernel%co%kernel, pot_ion_dummy, &
603 : !& enhartr, enxc, evxc, 0.d0, .false., 1)
604 :
605 : call H_potential(datacode,kernel,vhartr,pot_ion_dummy,&
606 : & enhartr,zero,.false.)
607 :
608 :
609 : #else
610 0 : BIGDFT_NOTENABLED_ERROR()
611 : if (.false.) write(std_out,*) nfft,nspden,icoulomb,usewvl,mpi_comm,me,nproc,nscforder,enhartr,&
612 : & ngfft(1),hgrid(1),rhor(1,1),vhartr(1)
613 : #endif
614 :
615 0 : end subroutine psolver_hartree
616 : !!***
617 :
618 : !!****f* ABINIT/psolver_kernel
619 : !! NAME
620 : !! psolver_kernel
621 : !!
622 : !! FUNCTION
623 : !! Build, get or free the kernel matrix used by the Poisson solver to compute the
624 : !! the convolution between 1/r and rho. The kernel is a saved variable. If
625 : !! this routine is called for building while a kernel already exists, it is not
626 : !! recomputed if all parameters (grid step and data size) are unchanged. Otherwise
627 : !! the kernel is freed and recompute again. The build action has a returned variable
628 : !! which is a pointer on the kernel. The get action also returns the kernel, or
629 : !! NULL if none has been associated.
630 : !!
631 : !! INPUTS
632 : !! iaction=0 to free all kernel allocated array,
633 : !! 1 to compute the kernel (parallel case),
634 : !! 2 to get it (parallel case),
635 : !! 3 to compute the kernel (sequential),
636 : !! 4 to get the sequential kernel.
637 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
638 : !!
639 : !! OUTPUT
640 : !! kernel= associated kernel on build (iaction = 1) and get action (iaction = 2).
641 : !!
642 : !! SOURCE
643 :
644 0 : subroutine psolver_kernel(hgrid, iaction, icoulomb, &
645 : & iproc, kernel, mpi_comm, ngfft, nproc, nscforder)
646 :
647 : #if defined HAVE_BIGDFT
648 : use BigDFT_API, only : coulomb_operator,nullify_coulomb_operator, &
649 : & deallocate_coulomb_operator,mpi_environment
650 : use poisson_solver, only : pkernel_init,pkernel_set
651 : #else
652 : use defs_wvltypes, only : coulomb_operator
653 : #endif
654 :
655 : !Arguments ------------------------------------
656 : !scalars
657 : integer,intent(in) :: iaction, icoulomb, mpi_comm, nscforder, iproc, nproc
658 : !arrays
659 : integer, intent(in) :: ngfft(3)
660 : type(coulomb_operator),intent(inout)::kernel
661 : real(dp),intent(in) :: hgrid(3)
662 :
663 : !Local variables-------------------------
664 : #if defined HAVE_BIGDFT
665 : !scalars
666 : integer,parameter :: igpu=0 !no GPUs
667 : !arrays
668 : integer, save :: kernel_scfOrder
669 : integer, save :: kernel_icoulomb
670 : integer, save :: data_size(3) = (/ -2, -2, -2 /)
671 : real(dp), save :: kernel_hgrid(3) ! Grid step used when creating the kernel.
672 : character(len = 1) :: geocode
673 : character(len=500) :: message
674 : integer :: current_size(3)
675 : type(coulomb_operator),save :: pkernel, kernelseq
676 : type(mpi_environment) :: mpi_env
677 : #endif
678 :
679 : ! *************************************************************************
680 :
681 : #if defined HAVE_BIGDFT
682 :
683 : if (icoulomb == 0) then
684 : ! The kernel is built with 'P'eriodic boundary counditions.
685 : geocode = 'P'
686 : else if (icoulomb == 1) then
687 : ! The kernel is built with 'F'ree boundary counditions.
688 : geocode = 'F'
689 : else if (icoulomb == 2) then
690 : ! The kernel is built with 'S'urface boundary counditions.
691 : geocode = 'S'
692 : end if
693 : current_size(:) = ngfft(1:3)
694 :
695 : !Initialise kernel_array pointer.
696 : if (maxval(data_size) == -2) then
697 : call nullify_coulomb_operator(pkernel)
698 : call nullify_coulomb_operator(kernelseq)
699 : end if
700 :
701 : !If iaction == 0, we free the kernel.
702 : if (iaction == 0) then
703 : if (associated(pkernel%kernel)) then
704 : write(message, "(A)") "Psolver_kernel() : deallocating pkernel..."
705 : call wrtout(std_out, message,'COLL')
706 :
707 : call deallocate_coulomb_operator(pkernel)
708 : end if
709 : if (associated(kernelseq%kernel)) then
710 : write(message, "(A)") "Psolver_kernel() : deallocating kernelseq..."
711 : call wrtout(std_out, message,'COLL')
712 :
713 : call deallocate_coulomb_operator(kernelseq)
714 : end if
715 : data_size = (/ -1, -1, -1 /)
716 : return
717 : end if
718 :
719 :
720 : !Action is build or get. We check the sizes before doing anything else.
721 :
722 : !!$!Get the size depending on wavelets calculations or not
723 : !!$ if (dtset%usewvl == 0) then
724 : !!$ hgrid(1) = rprimd(1, 1) / ngfft(1)
725 : !!$ hgrid(2) = rprimd(2, 2) / ngfft(2)
726 : !!$ hgrid(3) = rprimd(3, 3) / ngfft(3)
727 : !!$
728 : !!$ else
729 : !!$ hgrid(:) = 0.5d0 * wvl%h(:)
730 : !!$ current_size(1:3) = (/ wvl%Glr%d%n1i, wvl%Glr%d%n2i, wvl%Glr%d%n3i /)
731 : !!$ end if
732 :
733 : !Compute a new kernel if grid size has changed or if the kernel
734 : !has never been computed.
735 : if ((iaction == 1 .and. .not. associated(pkernel%kernel)) .or. &
736 : & (iaction == 3 .and. .not. associated(kernelseq%kernel)) .or. &
737 : & kernel_icoulomb /= icoulomb .or. &
738 : & data_size(1) /= current_size(1) .or. &
739 : & data_size(2) /= current_size(2) .or. &
740 : & data_size(3) /= current_size(3) .or. &
741 : & kernel_hgrid(1) /= hgrid(1) .or. &
742 : & kernel_hgrid(2) /= hgrid(2) .or. &
743 : & kernel_hgrid(3) /= hgrid(3) .or. &
744 : & kernel_scfOrder /= nscforder) then
745 : write(message, "(A,A,A,3I6)") "Psolver_kernel() : building kernel...", ch10, &
746 : & " | data dimensions:", current_size
747 : call wrtout(std_out, message, 'COLL')
748 :
749 : if (iaction == 1 .or. iaction == 2) then
750 : if (associated(pkernel%kernel)) then
751 : call deallocate_coulomb_operator(pkernel)
752 : end if
753 : mpi_env%mpi_comm = mpi_comm
754 : mpi_env%iproc = iproc
755 : mpi_env%nproc = nproc
756 : mpi_env%igroup = 0 ! no task groups
757 : mpi_env%ngroup = 1 ! no task groups
758 : pkernel= pkernel_init(.True.,iproc,nproc,igpu,geocode,&
759 : & current_size,hgrid,nscforder, mpi_env = mpi_env)
760 : call pkernel_set(pkernel,.True.)
761 : end if
762 :
763 : if (iaction == 3 .or. iaction == 4) then
764 : if (associated(kernelseq%kernel)) then
765 : call deallocate_coulomb_operator(kernelseq)
766 : end if
767 : mpi_env%mpi_comm = mpi_comm
768 : mpi_env%iproc = 0
769 : mpi_env%nproc = 1
770 : mpi_env%igroup = 0 ! no task groups
771 : mpi_env%ngroup = 1 ! no task groups
772 : kernelseq= pkernel_init(.True.,iproc,nproc,igpu,geocode,&
773 : & current_size,hgrid,nscforder, mpi_env = mpi_env)
774 : call pkernel_set(kernelseq,.True.)
775 : end if
776 :
777 : ! Storing variables which were used to make the kernel
778 : kernel_icoulomb = icoulomb
779 : data_size(:) = current_size(:)
780 : kernel_hgrid(:) = hgrid(:)
781 : kernel_scfOrder = nscforder
782 : end if
783 :
784 : ! Shallow copy if kernel has been associated.
785 : if (iaction == 1 .or. iaction == 2) then
786 : kernel = pkernel
787 : end if
788 : if (iaction == 3 .or. iaction == 4) then
789 : kernel = kernelseq
790 : end if
791 :
792 : #else
793 0 : BIGDFT_NOTENABLED_ERROR()
794 : if (.false.) write(std_out,*) iaction,icoulomb,mpi_comm,nscforder,iproc,nproc,ngfft(1),kernel%co,hgrid(1)
795 : #endif
796 :
797 0 : end subroutine psolver_kernel
798 : !!***
799 :
800 : end module m_psolver
801 : !!***
|
