Line data Source code
1 : !!****m* ABINIT/m_prcref
2 : !! NAME
3 : !! m_prcref
4 : !!
5 : !! FUNCTION
6 : !! Routines to precondition residual potential (or density) and forces.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, MT, PMA)
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 : !Include and generate MKL_RCI module
23 : !#if defined HAVE_LINALG_MKL_OMATCOPY
24 : !#include "mkl_rci.f90"
25 : !#endif
26 :
27 : module m_prcref
28 :
29 : use iso_c_binding
30 : use defs_basis
31 : use defs_wvltypes
32 : use m_errors
33 : use m_abicore
34 : use m_xmpi
35 : use m_xcdata
36 : use m_frskerker1
37 : use m_frskerker2
38 : use mod_prc_memory
39 : use m_dtset
40 : use m_precon
41 :
42 : use defs_datatypes, only : pseudopotential_type
43 : use defs_abitypes, only : MPI_type
44 : use m_time, only : timab
45 : use m_numeric_tools, only : dotproduct
46 : use m_geometry, only : xcart2xred, metric
47 : use m_cgtools, only : dotprod_vn, mean_fftr
48 : use m_mpinfo, only : ptabs_fourdp, destroy_mpi_enreg, initmpi_seq
49 : use m_pawtab, only : pawtab_type
50 : use m_pawrhoij, only : pawrhoij_type
51 : use m_rcpaw, only : rcpaw_type
52 : use m_extfpmd, only : extfpmd_type
53 : use m_fftcore, only : kgindex
54 : use m_fft, only : zerosym, indirect_parallel_fourier, fourdp
55 : use m_kg, only : getph
56 : use m_spacepar, only : hartre, laplacian
57 : use m_forces, only : fresid
58 : use m_atm2fft, only : atm2fft
59 : use m_rhotoxc, only : rhotoxc
60 : use m_mklocl, only : mklocl
61 : use m_mkcore, only : mkcore
62 : !use m_iterative_solvers, only : linsolve
63 :
64 : implicit none
65 :
66 : private
67 : !!***
68 :
69 : public :: prcref
70 : public :: prcref_PMA
71 : public :: moddiel ! Precondition the residual, using a model dielectric function.
72 : private :: chi0diel ! Precondition the residual, using a model chi0 operator defined in the object precon.
73 :
74 : !!***
75 :
76 : contains
77 : !!***
78 :
79 : !!****f* ABINIT/prcref
80 : !!
81 : !! NAME
82 : !! prcref
83 : !!
84 : !! FUNCTION
85 : !! Compute preconditioned residual potential (or density) and forces.
86 : !! iprcel, densfor_pred and iprcfc govern the choice of the preconditioner.
87 : !! Three tasks are done:
88 : !! 1) Preconditioning of the forces (residual has already been included)
89 : !! using the approximate force constant matrix. Get proposed
90 : !! change of atomic positions.
91 : !! 2) Precondition the residual, get first part of proposed trial
92 : !! potential change.
93 : !! 3) PAW only: precondition the rhoij residuals (simple preconditionning)
94 : !! 4) Take into account the proposed change of atomic positions to
95 : !! modify the proposed trial potential change.
96 : !!
97 : !! NOTE
98 : !! This routine is almost similar to prcref_PMA.F90 which is employed in
99 : !! case of potential mixing. Yet it has undergone strong changes simultaneously
100 : !! from two different sources at the same time which resulted in a splitting.
101 : !!
102 : !! INPUTS
103 : !! atindx(natom)=index table for atoms (see gstate.f)
104 : !! dielar(7)=input parameters for dielectric matrix:
105 : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
106 : !! dielstrt=number of the step at which the dielectric preconditioning begins.
107 : !! dtset <type(dataset_type)>=all input variables in this dataset
108 : !! | intxc=control xc quadrature
109 : !! | densfor_pred= not yet used here
110 : !! | iprcel= governs the preconditioning of the potential residual
111 : !! | 0 => simple model dielectric matrix, described by the
112 : !! | parameters dielng, diemac, diemix and diemixmag contained in dielar.
113 : !! | between 21 and 39 => until istep=dielstart, same as iprcel=0, then uses
114 : !! | the RPA dielectric matrix (routine dielmt)
115 : !! | between 41 and 49 => uses the RPA dielectric matrix (routine dielmt).
116 : !! | between 51 and 59 => uses the RPA dielectric matrix (routine dieltcel).
117 : !! | between 61 and 69 => uses the electronic dielectric matr (routine dieltcel).
118 : !! | between 71 and 79 => uses the real-space preconditioner based on Kerker prc (prcrskerkerN)
119 : !! | between 81 and 99 => reserved for futur version of the real-space preconditioner
120 : !! | between 141 and 169 -> same as between 41 and 69 but with a different periodicity: modulo(iprcel modulo (10))
121 : !! | iprcfc= governs the preconditioning of the forces
122 : !! | 0 => hessian is the identity matrix
123 : !! | 1 => hessian is 0.5 times the identity matrix
124 : !! | 2 => hessian is 0.25 times the identity matrix
125 : !! | ixc=exchange-correlation choice parameter.
126 : !! | natom=number of atoms
127 : !! | nspden=number of spin-density components
128 : !! | occopt=option for occupancies
129 : !! | prtvol=control print volume and debugging
130 : !! | typat(natom)=integer type for each atom in cell
131 : !! etotal=total ennergy
132 : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
133 : !! fcart(3,natom)=cartesian forces (hartree/bohr)
134 : !! ffttomix(nfft*(1-nfftprc/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
135 : !! gmet(3,3)=metric tensor in G space in Bohr**-2.
136 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
137 : !! istep= number of the step in the SCF cycle
138 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
139 : !! mgfft=maximum size of 1D FFTs
140 : !! moved_atm_inside= if 1, then the preconditioned forces
141 : !! as well as the preconditioned potential residual must be computed;
142 : !! otherwise, compute only the preconditioned potential residual.
143 : !! mpi_enreg=information about MPI parallelization
144 : !! my_natom=number of atoms treated by current processor
145 : !! nattyp(ntypat)=number of atoms of each type in cell.
146 : !! nfft=number of fft grid points
147 : !! nfftprc=size of FFT grid on which the potential residual will be preconditionned
148 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
149 : !! ngfftprc(18)=contain all needed information about 3D FFT for the grid corresponding to nfftprc
150 : !! nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description
151 : !! npawmix=-PAW only- number of spherical part elements to be mixed
152 : !! npwdiel=number of planewaves for dielectric matrix
153 : !! ntypat=number of types of atoms in cell.
154 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
155 : !! optreal=1 if residual potential is is REAL space, 2 if it is in RECIPROCAL SPACE
156 : !! optres=0: the array vresid contains a potential residual
157 : !! 1: the array vresid contains a density residual
158 : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
159 : !! Use here rhoij residuals (and gradients)
160 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
161 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
162 : !! rhog(2,nfft)=array for electron density in reciprocal space
163 : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
164 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
165 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
166 : !! the susceptibility (or density-density response) matrix in reciprocal space
167 : !! vresid(optreal*nfftprc,nspden)=residual potential
168 : !! vxc(nfft,nspden)=exchange-correlation potential (hartree)
169 : !! vhartr(nfft)=array for holding Hartree potential
170 : !! vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
171 : !! vpsp(nfft)=array for holding local psp
172 : !! xred(3,natom)=reduced dimensionless atomic coordinates
173 : !!
174 : !! OUTPUT
175 : !! dtn_pc(3,natom)=preconditioned change of atomic position,
176 : !! in reduced coordinates
177 : !! vrespc(optreal*nfftprc,nspden)=preconditioned residual of the potential
178 : !! ==== if psps%usepaw==1
179 : !! rhoijrespc(npawmix)= preconditionned rhoij residuals at output
180 : !!
181 : !! SIDE EFFECT
182 : !! dielinv(2,npwdiel,nspden,npwdiel,nspden)=
183 : !! inverse of the dielectric matrix in rec. space
184 : !! kxc(nfft,nkxc)=exchange-correlation kernel,
185 : !! needed if the electronic dielectric matrix is computed
186 : !! ===== if densfor_pred==3 .and. moved_atm_inside==1 =====
187 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
188 : !!
189 : !! SOURCE
190 :
191 12922 : subroutine prcref(atindx,dielar,dielinv,&
192 12922 : & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,&
193 12922 : & istep,kg_diel,kxc,&
194 : & mgfft,moved_atm_inside,mpi_enreg,my_natom,&
195 12922 : & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
196 12922 : & optreal,optres,pawrhoij,pawtab,ph1d,precon,psps,rhog,rhoijrespc,rhor,rprimd,&
197 12922 : & susmat,vhartr,vpsp,vresid,vrespc,vxc,wvl,wvl_den,xred,rcpaw,extfpmd)
198 :
199 : !Arguments-------------------------------
200 : !scalars
201 : integer,intent(in) :: dielstrt,istep,my_natom,mgfft,moved_atm_inside,n1xccc
202 : integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat,optreal,optres
203 : real(dp),intent(in) :: etotal,gsqcut
204 : type(MPI_type),intent(in) :: mpi_enreg
205 : type(dataset_type),intent(in) :: dtset
206 : type(pseudopotential_type),intent(in) :: psps
207 : type(wvl_internal_type), intent(in) :: wvl
208 : type(wvl_denspot_type), intent(inout) :: wvl_den
209 : type(rcpaw_type),intent(inout),pointer :: rcpaw
210 : type(extfpmd_type),intent(inout),pointer :: extfpmd
211 : type(precon_object), intent(inout) :: precon
212 :
213 : !arrays
214 : integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft))
215 : integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18)
216 : real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),rhog(2,nfft)
217 : real(dp),intent(in) :: rhor(nfft,dtset%nspden)
218 : real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
219 : real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden)
220 : real(dp),intent(in) :: vxc(nfft,dtset%nspden)
221 : real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
222 : real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc)
223 : real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft)
224 : real(dp),intent(inout) :: xred(3,dtset%natom)
225 : real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix),rprimd(3,3)
226 : real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden)
227 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
228 : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
229 :
230 : !Local variables-------------------------------
231 : !scalars
232 : integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1
233 : integer :: ipw2,iq,iq0,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm
234 : integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method
235 : real(dp) :: ai,ar,diemix,diemixmag,eei,bigexc,bigsxc
236 : real(dp) :: mixfac
237 : real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg
238 : logical :: computediel,non_magnetic_xc
239 : character(len=500) :: message
240 : type(xcdata_type) :: xcdata
241 : !arrays
242 : integer :: qprtrb(3)
243 12922 : integer,allocatable :: indpw_prc(:)
244 : real(dp) :: dummy6(6),dummy7(6),gprimd(3,3),qphon(3),rmet(3,3)
245 25844 : real(dp) :: vmean(dtset%nspden),vprtrb(2)
246 12922 : real(dp),allocatable :: dummy(:),dummy1(:),dummy2(:),dummy3(:),dummy4(:),dummy5(:),dummy8(:),dummy9(:)
247 12922 : real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:)
248 12922 : real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:)
249 12922 : real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:),rhor_new(:,:)
250 12922 : real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:)
251 12922 : real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:)
252 12922 : real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:)
253 12922 : logical,allocatable :: mask(:)
254 : ! *************************************************************************
255 :
256 : !Compute different geometric tensor, as well as ucvol, from rprimd
257 12922 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
258 :
259 : !1) Eventually take care of the forces
260 :
261 12922 : if(moved_atm_inside==1)then
262 0 : ABI_MALLOC(fcart_pc,(3,dtset%natom))
263 :
264 0 : if(dtset%iprcfc==0)then
265 0 : fcart_pc(:,:)=fcart(:,:)
266 : else
267 0 : fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:)
268 : end if
269 :
270 : ! Compute preconditioned delta xred from preconditioned fcart and rprimd
271 0 : call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc)
272 :
273 0 : ABI_FREE(fcart_pc)
274 : end if
275 :
276 : !#######################################################################
277 :
278 : !2) Take care of the potential residual
279 :
280 : !Compute the residuals corresponding to the solution
281 : !of an approximate realspace dielectric function according
282 : !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]]
283 12922 : if(dtset%iprcel>=71.and.dtset%iprcel<=79) then
284 0 : if (nfft==nfftprc) then
285 0 : if (dtset%iprcel<=78) then
286 : call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, &
287 0 : & gprimd,vresid,vrespc,rhor(:,1))
288 : else
289 : call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, &
290 0 : & vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol)
291 : end if
292 : else
293 : ! If preconditionning has to be done on a coarse grid,
294 : ! has to transfer several arrays
295 0 : ABI_MALLOC(work1,(nfftprc,dtset%nspden))
296 0 : ABI_MALLOC(work3,(nfftprc,dtset%nspden))
297 0 : ABI_MALLOC(work,(2*nfftprc))
298 0 : do ispden=1,dtset%nspden
299 0 : work(:)=vresid(:,ispden)
300 0 : call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,1,ngfftprc,0)
301 : end do
302 0 : ABI_FREE(work)
303 0 : if (dtset%iprcel<=78) then
304 0 : ABI_MALLOC(rhog_wk,(2,nfftprc))
305 0 : rhog_wk(:,:)=zero
306 0 : if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then
307 0 : nfftot=PRODUCT(ngfft(1:3))
308 : call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,&
309 0 : & ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot)
310 : else
311 0 : do ii=1,nfft
312 0 : if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii)
313 : end do
314 : end if
315 : call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),&
316 0 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
317 0 : ABI_MALLOC(work,(nfftprc))
318 0 : call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,1,ngfftprc,0)
319 : call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, &
320 0 : & gprimd,work1,work3,work)
321 0 : ABI_FREE(work)
322 : else
323 : call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, &
324 0 : & work1,work3,dtset%natom,xred,mpi_enreg,ucvol)
325 : end if
326 0 : do ispden=1,dtset%nspden
327 0 : call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,1,ngfftprc,0)
328 : end do
329 0 : ABI_FREE(work1)
330 0 : ABI_FREE(work3)
331 : end if
332 :
333 : else
334 :
335 12922 : if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then
336 12904 : cplex=optreal
337 12904 : qphon(:)=zero
338 : ! Simple scalar multiplication, or model dielectric function
339 12904 : call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,qphon,rprimd,vresid,vrespc)
340 :
341 : ! Use the inverse dielectric matrix in a small G sphere
342 18 : else if( (istep>=dielstrt .and. dtset%iprcel>=21 .and. dtset%iprcel<200) .or. modulo(dtset%iprcel,100)>=41 )then
343 :
344 : ! With dielop=1, the matrices will be computed when istep=dielstrt
345 : ! With dielop=2, the matrices will be computed when istep=dielstrt and 1
346 0 : dielop=1
347 0 : if(modulo(dtset%iprcel,100)>=41)dielop=2
348 0 : computediel = dtset%testsusmat(dielop, dielstrt, istep) !test if the matrix is to be computed
349 0 : if(computediel) then
350 : ! Compute the inverse dielectric matrix from the susceptibility matrix
351 : ! There are two routines for the RPA matrix, while for the electronic
352 : ! dielectric matrix, only dieltcel will do the work
353 0 : if(modulo(dtset%iprcel,100)<=49)then
354 : call dielmt(dielinv,gmet,kg_diel,&
355 0 : & npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat)
356 : else
357 0 : option=1
358 0 : if(modulo(dtset%iprcel,100)>=61)option=2
359 : call dieltcel(dielinv,gmet,kg_diel,kxc,&
360 0 : & nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%prtvol,susmat)
361 : end if
362 : end if
363 :
364 0 : ABI_MALLOC(work1,(2,nfftprc))
365 0 : ABI_MALLOC(work2,(optreal*nfftprc))
366 :
367 : ! Presently, one uses the inverse of the RPA dielectric matrix,
368 : ! for which spin must be averaged.
369 :
370 : ! Do fft from real space (work2) to G space (work1)
371 0 : if (optreal==1) then
372 0 : work2(:)=vresid(:,1)
373 : ! Must average over spins in the case of a potential residual
374 0 : if(dtset%nspden/=1.and.optres==0)work2(:)=(work2(:)+vresid(:,2))*half
375 0 : call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,1,ngfftprc,0)
376 : else
377 0 : work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/))
378 0 : if(dtset%nspden/=1.and.optres==0)work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half
379 : end if
380 :
381 : ! Multiply by restricted inverse of dielectric matrix.
382 : ! Must first copy relevant elements of work1 to a npwdiel-dimensioned array,
383 : ! then zero work1, operate with the dielinv matrix, and store in work1.
384 :
385 0 : ABI_MALLOC(vres_diel,(2,npwdiel))
386 0 : ABI_MALLOC(indpw_prc,(npwdiel))
387 0 : ABI_MALLOC(mask,(npwdiel))
388 0 : mask(:)=.true.
389 0 : call kgindex(indpw_prc,kg_diel,mask,mpi_enreg,ngfftprc,npwdiel)
390 :
391 0 : do ipw1=1,npwdiel
392 0 : if(mask(ipw1)) then
393 0 : vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1))
394 0 : vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1))
395 : end if
396 : end do
397 0 : work1(:,:)=zero
398 0 : do ipw1=1,npwdiel
399 0 : ar=zero ; ai=zero
400 :
401 : ! Use inverse of dielectric matrix (potential mixing)
402 0 : if (optres==0) then
403 0 : do ipw2=1,npwdiel
404 0 : if(mask(ipw2))then
405 : ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
406 0 : & -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
407 : ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
408 0 : & +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
409 : end if
410 : end do
411 : else
412 : ! Use symetric of inverse of dielectric matrix (density mixing)
413 0 : do ipw2=1,npwdiel
414 0 : if(mask(ipw2))then
415 : ar=ar+dielinv(1,ipw2,1,ipw1,1)*vres_diel(1,ipw2) &
416 0 : & +dielinv(2,ipw2,1,ipw1,1)*vres_diel(2,ipw2)
417 : ai=ai-dielinv(2,ipw2,1,ipw1,1)*vres_diel(1,ipw2) &
418 0 : & +dielinv(1,ipw2,1,ipw1,1)*vres_diel(2,ipw2)
419 : end if
420 : end do
421 : end if
422 : ! Must be careful not to count the diagonal 1 twice : it is added later,
423 : ! so must be subtracted now.
424 0 : call xmpi_sum(ar,mpi_enreg%comm_fft,ier)
425 0 : call xmpi_sum(ai,mpi_enreg%comm_fft,ier)
426 0 : if(mask(ipw1)) then
427 0 : work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1)
428 0 : work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1)
429 : end if !mask(ipw1)
430 : end do ! ipw1
431 0 : ABI_FREE(vres_diel)
432 0 : ABI_FREE(indpw_prc)
433 0 : ABI_FREE(mask)
434 : ! Fourier transform
435 0 : if (optreal==1) then
436 0 : call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,1,ngfftprc,0)
437 : else
438 0 : work2(:)=reshape(work1(:,:),(/nfftprc*2/))
439 : end if
440 :
441 : ! Add to get the preconditioned vresid, must be careful about spins.
442 0 : if(dtset%iprcel>=30)then
443 0 : diemix=dielar(4);diemixmag=abs(dielar(7))
444 0 : vrespc(:,1)=diemix*(vresid(:,1)+work2(:))
445 0 : if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:))
446 0 : if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=diemixmag*vresid(:,3:4)
447 0 : if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=diemixmag*vresid(:,2:dtset%nspden)
448 : else
449 0 : vrespc(:,1)=vresid(:,1)+work2(:)
450 0 : if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=vresid(:,2)+work2(:)
451 0 : if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=vresid(:,3:4)
452 0 : if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=vresid(:,2:dtset%nspden)
453 : end if
454 :
455 0 : ABI_FREE(work1)
456 0 : ABI_FREE(work2)
457 :
458 18 : else if (dtset%iprcel>=200 .and. dtset%iprcel<300) then
459 18 : call chi0diel(precon, dtset, mpi_enreg, optreal, optres, vresid, vrespc)
460 : ! Other choice ?
461 :
462 : else
463 : write(message, '(a,i3,a,a,a,a)' )&
464 0 : & 'From the calling routine, iprcel=',dtset%iprcel,ch10,&
465 0 : & 'The only allowed values are 0 or larger than 20.',ch10,&
466 0 : & 'Action: correct your input file.'
467 0 : ABI_ERROR(message)
468 : end if
469 : end if
470 : !#######################################################################
471 :
472 : !3) PAW : precondition the rhoij quantities (augmentation
473 : !occupancies) residuals. Use a simple preconditionning
474 : !with the same mixing factor as the model dielectric function.
475 : ! RCPAW : precondition the core occupations residuals
476 : !with the same mixing factor as the model dielectric function
477 : ! Extfpmd : precondition the extfpmd number of electrons
478 : !with the same mixing factor as the model dielectric function
479 :
480 12922 : if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then
481 : mixfac=one;mixfacmag=one
482 : else
483 12922 : mixfac=dielar(4);mixfacmag=abs(dielar(7))
484 : end if
485 12922 : if (psps%usepaw==1.and.my_natom>0) then
486 11261 : if (pawrhoij(1)%cplex_rhoij==1) then
487 : index=0
488 36945 : do iatom=1,my_natom
489 64362 : do iq=1,pawrhoij(iatom)%qphase
490 27417 : iq0=merge(0,pawrhoij(iatom)%lmn2_size,iq==1)
491 87312 : do ispden=1,pawrhoij(iatom)%nspden
492 32478 : mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
493 2130107 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
494 2070212 : index=index+1;klmn=iq0+pawrhoij(iatom)%kpawmix(kmix)
495 2102690 : rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden)
496 : end do
497 : end do
498 : end do
499 : end do
500 : else
501 : index=-1
502 4212 : do iatom=1,my_natom
503 6691 : do iq=1,pawrhoij(iatom)%qphase
504 2479 : iq0=merge(0,2*pawrhoij(iatom)%lmn2_size,iq==1)
505 13833 : do ispden=1,pawrhoij(iatom)%nspden
506 8875 : mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
507 918130 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
508 906776 : index=index+2;klmn=iq0+2*pawrhoij(iatom)%kpawmix(kmix)-1
509 2729203 : rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
510 : end do
511 : end do
512 : end do
513 : end do
514 : end if
515 : end if
516 :
517 12922 : if (psps%usepaw==1.and.associated(rcpaw)) then
518 18 : mixfac_eff=mixfac!;if (ispden>1) mixfac_eff=mixfacmag
519 36 : do iatom=1,rcpaw%ntypat
520 36 : if(rcpaw%atm(iatom)%zcore_orig>zero) then
521 108 : rcpaw%atm(iatom)%occ_respc=mixfac_eff*rcpaw%atm(iatom)%occ_res
522 : endif
523 : enddo
524 : endif
525 :
526 12922 : if(associated(extfpmd)) then
527 83 : mixfac_eff=mixfac!;if (ispden>1) mixfac_eff=mixfacmag
528 83 : extfpmd%nelect_respc=mixfac_eff*extfpmd%nelect_res
529 : endif
530 :
531 :
532 : !#######################################################################
533 :
534 : !4) Take care of the change of atomic positions
535 : !Note : this part is very demanding on memory...
536 : !however, since this algorithm is still in development,
537 : !it was NOT included in the estimation provided by memory.f
538 12922 : if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then
539 :
540 : ! Not yet compatible with resid given in reciprocal space
541 0 : if (optreal/=1) then
542 : write(message, '(5a)' )&
543 0 : & 'From the calling routine, densfor_pred=3',ch10,&
544 0 : & 'You cannot use residuals in reciprocal space.',ch10,&
545 0 : & 'Action: correct your input file.'
546 0 : ABI_ERROR(message)
547 : end if
548 : ! Not compatible with non-collinear magnetism
549 0 : if(dtset%nspden==4)then
550 0 : ABI_ERROR('densfor_pred=3 does not work for nspden=4 !')
551 : end if
552 :
553 0 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
554 0 : nfftot=PRODUCT(ngfft(1:3))
555 :
556 0 : if (optres==0) then ! Array vresid contains a potential residual
557 : ! -----------------------------------------------------------------
558 :
559 : ! First subtract the current local, hartree and exchange correlation potentials
560 0 : do ispden=1,min(dtset%nspden,2)
561 0 : vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden)
562 : end do
563 0 : if (dtset%nspden==4) then
564 0 : do ispden=3,4
565 0 : vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden)
566 : end do
567 : end if
568 :
569 : ! Compute the modified density, in rhor_wk
570 0 : option=2
571 0 : ABI_MALLOC(gresid,(3,dtset%natom))
572 0 : ABI_MALLOC(grxc,(3,dtset%natom))
573 0 : ABI_MALLOC(rhor_wk,(nfft,dtset%nspden))
574 0 : ABI_MALLOC(rhor_wk0,(nfft,dtset%nspden))
575 0 : ABI_MALLOC(xred_wk,(3,dtset%natom))
576 0 : xred_wk(:,:)=xred(:,:)+dtn_pc(:,:)
577 : call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
578 : & ntypat,option,pawtab,rhor,rprimd,&
579 0 : & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp)
580 :
581 : ! Compute up+down rhog_wk(G) by fft
582 0 : ABI_MALLOC(work,(nfft))
583 0 : ABI_MALLOC(rhog_wk,(2,nfft))
584 0 : work(:)=rhor_wk(:,1)
585 0 : call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,1,ngfft,0)
586 0 : ABI_FREE(work)
587 :
588 : ! Compute structure factor phases for new atomic pos:
589 0 : call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk)
590 :
591 : ! Compute local ionic pseudopotential vpsp:
592 : ! and core electron density xccc3d, if needed.
593 0 : n3xccc=0;if (n1xccc/=0) n3xccc=nfft
594 0 : ABI_MALLOC(xccc3d,(n3xccc))
595 0 : ABI_MALLOC(vpsp_wk,(nfft))
596 0 : vprtrb(1:2)=zero
597 :
598 : ! Determine by which method the local ionic potential and/or
599 : ! the pseudo core charge density contributions have to be computed
600 : ! Local ionic potential:
601 : ! Method 1: PAW
602 : ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
603 0 : vloc_method=1;if (psps%usepaw==0) vloc_method=2
604 0 : if (dtset%icoulomb>0) vloc_method=2
605 0 : if (psps%usewvl==1) vloc_method=2
606 : ! Pseudo core charge density:
607 : ! Method 1: PAW, nc_xccc_gspace
608 : ! Method 2: Norm-conserving PP, wavelets
609 : coredens_method=1;if (psps%usepaw==0) coredens_method=2
610 0 : if (psps%nc_xccc_gspace==1) coredens_method=1
611 0 : if (psps%nc_xccc_gspace==0) coredens_method=2
612 0 : if (psps%usewvl==1) coredens_method=2
613 :
614 : ! Local ionic potential and/or pseudo core charge by method 1
615 0 : if (vloc_method==1.or.coredens_method==1) then
616 0 : optv=0;if (vloc_method==1) optv=1
617 0 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
618 0 : optatm=1;optdyfr=0;optgr=0;optstr=0;optn2=1;opteltfr=0
619 : ! Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused...
620 : call atm2fft(atindx,xccc3d,vpsp,dummy,dummy2,dummy9,dummy1,gmet,gprimd,dummy3,dummy4,gsqcut,&
621 : & mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,&
622 : & optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dtset%rcut,dummy5,rprimd,dummy6,dummy7,&
623 : & ucvol,psps%usepaw,dummy8,dummy8,dummy8,vprtrb,psps%vlspl,&
624 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
625 0 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
626 : end if
627 :
628 : ! Local ionic potential by method 2
629 0 : if (vloc_method==2) then
630 0 : option=1
631 0 : ABI_MALLOC(dyfrlo_indx,(3,3,dtset%natom))
632 0 : ABI_MALLOC(grtn_indx,(3,dtset%natom))
633 : call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,&
634 : & mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,&
635 : & ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,&
636 0 : & ucvol,vprtrb,vpsp_wk,wvl,wvl_den,xred)
637 0 : ABI_FREE(dyfrlo_indx)
638 0 : ABI_FREE(grtn_indx)
639 : end if
640 :
641 : ! Pseudo core electron density by method 2
642 0 : if (coredens_method==2.and.n1xccc/=0) then
643 0 : option=1
644 0 : ABI_MALLOC(dyfrx2,(3,3,dtset%natom))
645 0 : ABI_MALLOC(grxc_indx,(3,dtset%natom))
646 : call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
647 : & n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,&
648 0 : & psps%xccc1d,xccc3d,xred_wk)
649 0 : ABI_FREE(dyfrx2)
650 0 : ABI_FREE(grxc_indx)
651 : end if
652 :
653 : ! Compute Hartree+xc potentials
654 0 : ABI_MALLOC(vxc_wk,(nfft,dtset%nspden))
655 0 : ABI_MALLOC(vhartr_wk,(nfft))
656 0 : option=1
657 :
658 : call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfft,ngfft,&
659 0 : &dtset%nkpt,dtset%rcut,rhog_wk,rprimd,dtset%vcutgeo,vhartr_wk)
660 :
661 : ! Prepare the call to rhotoxc
662 0 : call xcdata_init(xcdata,dtset=dtset)
663 0 : nk3xc=1 ; non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
664 : call rhotoxc(bigexc,bigsxc,kxc,mpi_enreg,nfft,ngfft,&
665 : & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor_wk,rprimd,1,&
666 0 : & vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk)
667 0 : ABI_FREE(xccc3d)
668 :
669 : ! Sum all contributions
670 0 : do ispden=1,min(dtset%nspden,2)
671 0 : do ifft=1,nfft
672 0 : vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden)
673 : end do
674 : end do
675 0 : if (dtset%nspden==4) then
676 0 : do ispden=3,4
677 0 : do ifft=1,nfft
678 0 : vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden)
679 : end do
680 : end do
681 : end if
682 0 : call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
683 0 : if(dtset%nspden==2) then
684 0 : vmean(1)=half*(vmean(1)+vmean(2))
685 0 : vmean(2)=vmean(1)
686 : end if
687 0 : do ispden=1,dtset%nspden
688 0 : vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden)
689 : end do
690 0 : ABI_FREE(gresid)
691 0 : ABI_FREE(grxc)
692 0 : ABI_FREE(rhog_wk)
693 0 : ABI_FREE(rhor_wk)
694 0 : ABI_FREE(rhor_wk0)
695 0 : ABI_FREE(xred_wk)
696 0 : ABI_FREE(vhartr_wk)
697 0 : ABI_FREE(vpsp_wk)
698 0 : ABI_FREE(vxc_wk)
699 :
700 : else ! Array vresid contains a density residual
701 : ! -----------------------------------------------------------------
702 :
703 : ! Only have to compute the modified preconditionned density residual
704 0 : option=2
705 0 : ABI_MALLOC(gresid,(3,dtset%natom))
706 0 : ABI_MALLOC(grxc,(3,dtset%natom))
707 0 : ABI_MALLOC(rhor_new,(nfft,dtset%nspden))
708 0 : ABI_MALLOC(rhor_wk,(nfft,dtset%nspden))
709 0 : ABI_MALLOC(rhor_wk0,(nfft,dtset%nspden))
710 0 : ABI_MALLOC(xred_wk,(3,dtset%natom))
711 0 : xred_wk(:,:)=xred(:,:)+dtn_pc(:,:)
712 0 : rhor_new(:,1)=rhor(:,1)+vrespc(:,1)
713 0 : if (dtset%nspden==2) then
714 0 : rhor_new(:,1)=rhor_new(:,1)+vrespc(:,2)
715 0 : rhor_new(:,2)=rhor(:,2)+vrespc(:,1)
716 : end if
717 : call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
718 : & ntypat,option,pawtab,rhor,rprimd,&
719 0 : & ucvol,rhor_wk0,xred_wk,xred,psps%znuclpsp)
720 : call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
721 : & ntypat,option,pawtab,rhor_new,rprimd,&
722 0 : & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp)
723 0 : vrespc(:,1)=rhor_wk(:,dtset%nspden)-rhor_wk0(:,dtset%nspden)
724 0 : if (dtset%nspden==2) vrespc(:,2)=rhor_wk(:,1)-rhor_wk0(:,1)-vrespc(:,1)
725 0 : ABI_FREE(gresid)
726 0 : ABI_FREE(grxc)
727 0 : ABI_FREE(rhor_new)
728 0 : ABI_FREE(rhor_wk)
729 0 : ABI_FREE(rhor_wk0)
730 0 : ABI_FREE(xred_wk)
731 : end if
732 :
733 : end if
734 :
735 12922 : end subroutine prcref
736 : !!***
737 :
738 : !!****f* ABINIT/prcref_PMA
739 : !!
740 : !! NAME
741 : !! prcref_PMA
742 : !!
743 : !! FUNCTION
744 : !! Compute preconditioned residual potential (or density) and forces.
745 : !! iprcel, densfor_pred and iprcfc govern the choice of the preconditioner.
746 : !! Three tasks are done:
747 : !! 1) Preconditioning of the forces (residual has already been included)
748 : !! using the approximate force constant matrix. Get proposed
749 : !! change of atomic positions.
750 : !! 2) Precondition the residual, get first part of proposed trial
751 : !! potential change.
752 : !! 3) PAW only: precondition the rhoij residuals (simple preconditionning)
753 : !! 4) Take into account the proposed change of atomic positions to
754 : !! modify the proposed trial potential change.
755 : !!
756 : !! NOTE
757 : !! This routine is almost similar to prcref.F90 which is employed in
758 : !! case of density mixing. Yet it has undergone strong changes simultaneously
759 : !! from two different sources at the same time which resulted in a splitting.
760 : !!
761 : !! INPUTS
762 : !! atindx(natom)=index table for atoms (see gstate.f)
763 : !! dielar(7)=input parameters for dielectric matrix:
764 : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
765 : !! dielstrt=number of the step at which the dielectric preconditioning begins.
766 : !! dtset <type(dataset_type)>=all input variables in this dataset
767 : !! | intxc=control xc quadrature
768 : !! | densfor_pred= not yet used here
769 : !! | iprcel= governs the preconditioning of the potential residual
770 : !! | 0 => simple model dielectric matrix, described by the
771 : !! | parameters dielng, diemac, diemix and diemixmag contained in dielar.
772 : !! | between 21 and 39 => until istep=dielstart, same as iprcel=0, then uses
773 : !! | the RPA dielectric matrix (routine dielmt)
774 : !! | between 41 and 49 => uses the RPA dielectric matrix (routine dielmt).
775 : !! | between 51 and 59 => uses the RPA dielectric matrix (routine dieltcel).
776 : !! | between 61 and 69 => uses the electronic dielectric matr (routine dieltcel).
777 : !! | between 71 and 79 => uses the real-space preconditioner based on Kerker prc (prcrskerkerN)
778 : !! | between 81 and 99 => reserved for futur version of the real-space preconditioner
779 : !! | between 141 and 169 -> same as between 41 and 69 but with a different periodicity: modulo(iprcel modulo (10))
780 : !! | iprcfc= governs the preconditioning of the forces
781 : !! | 0 => hessian is the identity matrix
782 : !! | 1 => hessian is 0.5 times the identity matrix
783 : !! | 2 => hessian is 0.25 times the identity matrix
784 : !! | ixc=exchange-correlation choice parameter.
785 : !! | natom=number of atoms
786 : !! | nspden=number of spin-density components
787 : !! | occopt=option for occupancies
788 : !! | prtvol=control print volume and debugging
789 : !! | typat(natom)=integer type for each atom in cell
790 : !! fcart(3,natom)=cartesian forces (hartree/bohr)
791 : !! ffttomix(nfft*(1-nfftprc/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
792 : !! gmet(3,3)=metric tensor in G space in Bohr**-2.
793 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
794 : !! istep= number of the step in the SCF cycle
795 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
796 : !! mgfft=maximum size of 1D FFTs
797 : !! moved_atm_inside= if 1, then the preconditioned forces
798 : !! as well as the preconditioned potential residual must be computed;
799 : !! otherwise, compute only the preconditioned potential residual.
800 : !! mpi_enreg=information about MPI parallelization
801 : !! my_natom=number of atoms treated by current processor
802 : !! nattyp(ntypat)=number of atoms of each type in cell.
803 : !! nfft=number of fft grid points
804 : !! nfftprc=size of FFT grid on which the potential residual will be preconditionned
805 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
806 : !! ngfftprc(18)=contain all needed information about 3D FFT for the grid corresponding to nfftprc
807 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
808 : !! npawmix=-PAW only- number of spherical part elements to be mixed
809 : !! npwdiel=number of planewaves for dielectric matrix
810 : !! ntypat=number of types of atoms in cell.
811 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
812 : !! optreal=1 if residual potential is is REAL space, 2 if it is in RECIPROCAL SPACE
813 : !! optres=0: the array vresid contains a potential residual
814 : !! 1: the array vresid contains a density residual
815 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
816 : !! Use here rhoij residuals (and gradients)
817 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
818 : !! rhog(2,nfft)=array for electron density in reciprocal space
819 : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
820 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
821 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
822 : !! the susceptibility (or density-density response) matrix in reciprocal space
823 : !! vresid(optreal*nfftprc,nspden)=residual potential
824 : !! vxc(nfft,nspden)=exchange-correlation potential (hartree)
825 : !! vhartr(nfft)=array for holding Hartree potential
826 : !! vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
827 : !! vpsp(nfft)=array for holding local psp
828 : !! xred(3,natom)=reduced dimensionless atomic coordinates
829 : !!
830 : !! etotal
831 : !! pawtab
832 : !!
833 : !! OUTPUT
834 : !! dtn_pc(3,natom)=preconditioned change of atomic position,
835 : !! in reduced coordinates
836 : !! vrespc(optreal*nfftprc,nspden)=preconditioned residual of the potential
837 : !! ==== if psps%usepaw==1
838 : !! rhoijrespc(npawmix)= preconditionned rhoij residuals at output
839 : !!
840 : !! SIDE EFFECT
841 : !! dielinv(2,npwdiel,nspden,npwdiel,nspden)=
842 : !! inverse of the dielectric matrix in rec. space
843 : !! kxc(nfft,nkxc)=exchange-correlation kernel,
844 : !! needed if the electronic dielectric matrix is computed
845 : !! ===== if densfor_pred==3 .and. moved_atm_inside==1 =====
846 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
847 : !!
848 : !! SOURCE
849 :
850 28280 : subroutine prcref_PMA(atindx,dielar,dielinv,&
851 28280 : & dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,&
852 28280 : & istep,kg_diel,kxc,&
853 : & mgfft,moved_atm_inside,mpi_enreg,my_natom,&
854 28280 : & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
855 28280 : & optreal,optres,pawrhoij,ph1d,precon,psps,rhog, rhoijrespc,rhor,rprimd,&
856 28280 : & susmat,vhartr,vpsp,vresid,vrespc,vxc,xred,&
857 28280 : & etotal,pawtab,wvl)
858 :
859 : !Arguments-------------------------------
860 : !variables used for tfvw
861 : !scalars
862 : integer,intent(in) :: dielstrt,istep,mgfft,moved_atm_inside,my_natom,n1xccc
863 : integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat
864 : integer,intent(in) :: optreal,optres
865 : real(dp),intent(in) :: etotal,gsqcut
866 : type(MPI_type),intent(in) :: mpi_enreg
867 : type(dataset_type),intent(in) :: dtset
868 : type(pseudopotential_type),intent(in) :: psps
869 : type(wvl_data), intent(inout) :: wvl
870 : type(precon_object), intent(inout) :: precon
871 :
872 : !arrays
873 : integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft))
874 : integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18)
875 : real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom)
876 : real(dp),intent(in) :: rhog(2,nfft)
877 : real(dp),intent(in) :: rhor(nfft,dtset%nspden)
878 : real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
879 : real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden)
880 : real(dp),intent(in) :: vxc(nfft,dtset%nspden)
881 : real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
882 : real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc)
883 : real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft)
884 : real(dp),intent(inout) :: xred(3,dtset%natom),rprimd(3,3)
885 : real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix)
886 : real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden)
887 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
888 : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
889 :
890 : !Local variables-------------------------------
891 : !scalars
892 : integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1
893 : integer :: ipw2,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm
894 : integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method
895 : real(dp) :: ai,ar,diemix,diemixmag,eei,bigexc,bigsxc
896 : real(dp) :: mixfac
897 : real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg
898 : logical :: computediel,non_magnetic_xc
899 : character(len=500) :: message
900 : type(xcdata_type) :: xcdata
901 : !arrays
902 : integer :: qprtrb(3)
903 28280 : integer,allocatable :: indpw_prc(:)
904 : real(dp) :: dummy6(6),gprimd(3,3),qphon(3),rmet(3,3)
905 56560 : real(dp) :: vmean(dtset%nspden),vprtrb(2)
906 28280 : real(dp),allocatable :: dummy_in(:)
907 : real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0),dummy_out7(0)
908 28280 : real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:)
909 28280 : real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:)
910 28280 : real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:)
911 28280 : real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:)
912 28280 : real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:)
913 28280 : real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:)
914 28280 : logical,allocatable :: mask(:)
915 : ! *************************************************************************
916 :
917 28280 : if(optres==1)then
918 0 : ABI_ERROR('density mixing (optres=1) not admitted!')
919 : end if
920 :
921 : !Compute different geometric tensor, as well as ucvol, from rprimd
922 28280 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
923 :
924 : !1) Eventually take care of the forces
925 :
926 28280 : if(moved_atm_inside==1)then
927 330 : ABI_MALLOC(fcart_pc,(3,dtset%natom))
928 :
929 110 : if(dtset%iprcfc==0)then
930 288 : fcart_pc(:,:)=fcart(:,:)
931 : else
932 1158 : fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:)
933 : end if
934 :
935 : ! Compute preconditioned delta xred from preconditioned fcart and rprimd
936 110 : call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc)
937 :
938 110 : ABI_FREE(fcart_pc)
939 : end if
940 :
941 : !#######################################################################
942 :
943 : !2) Take care of the potential residual
944 :
945 : !Compute the residuals corresponding to the solution
946 : !of an approximate realspace dielectric function according
947 : !to X. Gonze PRB vol54 nb7 p4383 (1996) [[cite:Gonze1996]]
948 28280 : if(dtset%iprcel>=71.and.dtset%iprcel<=79) then
949 11 : if (nfft==nfftprc) then
950 11 : if (dtset%iprcel<=78) then
951 : call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, &
952 7 : & gprimd,vresid,vrespc,rhor(:,1))
953 : else
954 : call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, &
955 4 : & vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol)
956 : end if
957 : else
958 : ! If preconditionning has to be done on a coarse grid,
959 : ! has to transfer several arrays
960 0 : ABI_MALLOC(work1,(nfftprc,dtset%nspden))
961 0 : ABI_MALLOC(work3,(nfftprc,dtset%nspden))
962 0 : ABI_MALLOC(work,(2*nfftprc))
963 0 : do ispden=1,dtset%nspden
964 0 : work(:)=vresid(:,ispden)
965 0 : call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,1,ngfftprc,0)
966 : end do
967 0 : ABI_FREE(work)
968 0 : if (dtset%iprcel<=78) then
969 0 : ABI_MALLOC(rhog_wk,(2,nfftprc))
970 0 : rhog_wk(:,:)=zero
971 0 : if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then
972 0 : nfftot=PRODUCT(ngfft(1:3))
973 : call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,&
974 0 : & ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot)
975 : else
976 0 : do ii=1,nfft
977 0 : if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii)
978 : end do
979 : end if
980 : call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),&
981 0 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
982 0 : ABI_MALLOC(work,(nfftprc))
983 0 : call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,1,ngfftprc,0)
984 : call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, &
985 0 : & gprimd,work1,work3,work)
986 0 : ABI_FREE(work)
987 : else
988 : call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, &
989 0 : & work1,work3,dtset%natom,xred,mpi_enreg,ucvol)
990 : end if
991 0 : do ispden=1,dtset%nspden
992 0 : call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,1,ngfftprc,0)
993 : end do
994 0 : ABI_FREE(work1)
995 0 : ABI_FREE(work3)
996 : end if
997 :
998 : else
999 :
1000 28269 : if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then
1001 27923 : cplex=optreal
1002 27923 : qphon(:)=zero
1003 : ! Simple scalar multiplication, or model dielectric function
1004 27923 : call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,qphon,rprimd,vresid,vrespc)
1005 :
1006 : ! Use the inverse dielectric matrix in a small G sphere
1007 346 : else if( (istep>=dielstrt .and. dtset%iprcel>=21 .and. dtset%iprcel<200) .or. modulo(dtset%iprcel,100)>=41 )then
1008 :
1009 : ! Wnith dielop=1, the matrices will be computed when istep=dielstrt
1010 : ! With dielop=2, the matrices will be computed when istep=dielstrt and 1
1011 346 : dielop=1
1012 346 : if(modulo(dtset%iprcel,100)>=41)dielop=2
1013 346 : computediel = dtset%testsusmat(dielop, dielstrt, istep) !test if the matrix is to be computed
1014 346 : if(computediel) then
1015 : ! Compute the inverse dielectric matrix from the susceptibility matrix
1016 : ! There are two routines for the RPA matrix, while for the electronic
1017 : ! dielectric matrix, only dieltcel will do the work
1018 74 : if(modulo(dtset%iprcel,100)<=49)then
1019 : call dielmt(dielinv,gmet,kg_diel,&
1020 68 : & npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat)
1021 : else
1022 6 : option=1
1023 6 : if(modulo(dtset%iprcel,100)>=61)option=2
1024 : call dieltcel(dielinv,gmet,kg_diel,kxc,&
1025 6 : & nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%prtvol,susmat)
1026 : end if
1027 : end if
1028 :
1029 1038 : ABI_MALLOC(work1,(2,nfftprc))
1030 1038 : ABI_MALLOC(work2,(optreal*nfftprc))
1031 :
1032 : ! Presently, one uses the inverse of the RPA dielectric matrix,
1033 : ! for which spin must be averaged.
1034 :
1035 : ! Do fft from real space (work2) to G space (work1)
1036 346 : if (optreal==1) then
1037 3628122 : work2(:)=vresid(:,1)
1038 : ! Must average over spins if needed.
1039 228346 : if(dtset%nspden/=1)work2(:)=(work2(:)+vresid(:,2))*half
1040 346 : call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,1,ngfftprc,0)
1041 : else
1042 0 : work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/))
1043 0 : if (dtset%nspden/=1) work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half
1044 : end if
1045 :
1046 : ! Multiply by restricted inverse of dielectric matrix.
1047 : ! Must first copy relevant elements of work1 to a npwdiel-dimensioned array,
1048 : ! then zero work1, operate with the dielinv matrix, and store in work1.
1049 :
1050 1038 : ABI_MALLOC(vres_diel,(2,npwdiel))
1051 1038 : ABI_MALLOC(indpw_prc,(npwdiel))
1052 692 : ABI_MALLOC(mask,(npwdiel))
1053 64568 : mask(:)=.true.
1054 346 : call kgindex(indpw_prc,kg_diel,mask,mpi_enreg,ngfftprc,npwdiel)
1055 64568 : do ipw1=1,npwdiel
1056 64568 : if(mask(ipw1)) then
1057 64222 : vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1))
1058 64222 : vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1))
1059 : end if
1060 : end do
1061 :
1062 10883674 : work1(:,:)=zero
1063 64568 : do ipw1=1,npwdiel
1064 64222 : ar=zero ; ai=zero
1065 :
1066 : ! Use inverse of dielectric matrix (potential mixing)
1067 17247856 : do ipw2=1,npwdiel
1068 17247856 : if(mask(ipw2))then
1069 : ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
1070 17183634 : & -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
1071 : ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) &
1072 17183634 : & +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2)
1073 : end if
1074 : end do
1075 : ! Must be careful not to count the diagonal 1 twice : it is added later,
1076 : ! so must be subtracted now.
1077 64222 : call xmpi_sum(ar,mpi_enreg%comm_fft,ier)
1078 64222 : call xmpi_sum(ai,mpi_enreg%comm_fft,ier)
1079 64568 : if(mask(ipw1)) then
1080 64222 : work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1)
1081 64222 : work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1)
1082 : end if !mask(ipw1)
1083 : end do ! ipw1
1084 346 : ABI_FREE(vres_diel)
1085 346 : ABI_FREE(indpw_prc)
1086 346 : ABI_FREE(mask)
1087 :
1088 : ! Fourier transform
1089 346 : if (optreal==1) then
1090 346 : call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,1,ngfftprc,0)
1091 : else
1092 0 : work2(:)=reshape(work1(:,:),(/nfftprc*2/))
1093 : end if
1094 :
1095 : ! Add to get the preconditioned vresid, must be careful about spins.
1096 346 : if(dtset%iprcel>=30)then
1097 309 : diemix=dielar(4);diemixmag=abs(dielar(7))
1098 3286133 : vrespc(:,1)=diemix*(vresid(:,1)+work2(:))
1099 309 : if(dtset%nspden/=1)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:))
1100 309 : if(dtset%nspden==4)vrespc(:,3:4)=diemixmag*vresid(:,3:4)
1101 : else
1102 341989 : vrespc(:,1)=vresid(:,1)+work2(:)
1103 228037 : if(dtset%nspden/=1)vrespc(:,2)=vresid(:,2)+work2(:)
1104 37 : if(dtset%nspden==4)vrespc(:,3:4)=vresid(:,3:4)
1105 : end if
1106 :
1107 346 : ABI_FREE(work1)
1108 346 : ABI_FREE(work2)
1109 :
1110 0 : else if (dtset%iprcel>=200 .and. dtset%iprcel<300) then
1111 0 : call chi0diel(precon, dtset, mpi_enreg, optreal, optres, vresid, vrespc)
1112 :
1113 : ! Other choice ?
1114 : else
1115 : write(message, '(a,i0,a,a,a,a)' )&
1116 0 : & 'From the calling routine, iprcel= ',dtset%iprcel,ch10,&
1117 0 : & 'The only allowed values are 0 or larger than 20.',ch10,&
1118 0 : & 'Action: correct your input file.'
1119 0 : ABI_ERROR(message)
1120 : end if
1121 : end if
1122 : !#######################################################################
1123 :
1124 : !3) PAW only : precondition the rhoij quantities (augmentation
1125 : !occupancies) residuals. Use a simple preconditionning
1126 : !with the same mixing factor as the model dielectric function.
1127 :
1128 28280 : if (psps%usepaw==1.and.my_natom>0) then
1129 782 : if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then
1130 : mixfac=one;mixfacmag=one
1131 : else
1132 782 : mixfac=dielar(4);mixfacmag=abs(dielar(7))
1133 : end if
1134 782 : if (pawrhoij(1)%cplex_rhoij==1) then
1135 : index=0
1136 1921 : do iatom=1,my_natom
1137 3675 : do ispden=1,pawrhoij(iatom)%nspden
1138 1754 : mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
1139 243057 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
1140 240164 : index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix)
1141 241918 : rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden)
1142 : end do
1143 : end do
1144 : end do
1145 : else
1146 : index=-1
1147 0 : do iatom=1,my_natom
1148 0 : do ispden=1,pawrhoij(iatom)%nspden
1149 0 : mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
1150 0 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
1151 0 : index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
1152 0 : rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
1153 : end do
1154 : end do
1155 : end do
1156 : end if
1157 : end if
1158 : !#######################################################################
1159 :
1160 : !4) Take care of the change of atomic positions
1161 : !Note : this part is very demanding on memory...
1162 : !however, since this algorithm is still in development,
1163 : !it was NOT included in the estimation provided by memory.f
1164 28280 : if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then
1165 :
1166 : ! Not yet compatible with resid given in reciprocal space
1167 41 : if (optreal/=1) then
1168 : write(message, '(5a)' )&
1169 0 : & 'From the calling routine, densfor_pred=3',ch10,&
1170 0 : & 'You cannot use residuals in reciprocal space.',ch10,&
1171 0 : & 'Action: correct your input file.'
1172 0 : ABI_ERROR(message)
1173 : end if
1174 :
1175 : ! Not compatible with non-collinear magnetism
1176 41 : if(dtset%nspden==4)then
1177 0 : ABI_ERROR('densfor_pred=3 does not work for nspden=4!')
1178 : end if
1179 :
1180 41 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1181 164 : nfftot=PRODUCT(ngfft(1:3))
1182 :
1183 : ! First subtract the current local, hartree and exchange correlation potentials
1184 123 : do ispden=1,min(dtset%nspden,2)
1185 984123 : vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden)
1186 : end do
1187 41 : if (dtset%nspden==4) then
1188 0 : do ispden=3,4
1189 0 : vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden)
1190 : end do
1191 : end if
1192 :
1193 : ! Compute the modified density, in rhor_wk
1194 41 : option=2
1195 123 : ABI_MALLOC(gresid,(3,dtset%natom))
1196 82 : ABI_MALLOC(grxc,(3,dtset%natom))
1197 164 : ABI_MALLOC(rhor_wk,(nfft,dtset%nspden))
1198 123 : ABI_MALLOC(rhor_wk0,(nfft,dtset%nspden))
1199 82 : ABI_MALLOC(xred_wk,(3,dtset%natom))
1200 369 : xred_wk(:,:)=xred(:,:)+dtn_pc(:,:)
1201 :
1202 : call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,&
1203 : & ntypat,option,pawtab,rhor,rprimd,&
1204 41 : & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp)
1205 :
1206 : ! Compute up+down rhog_wk(G) by fft
1207 123 : ABI_MALLOC(work,(nfft))
1208 123 : ABI_MALLOC(rhog_wk,(2,nfft))
1209 492041 : work(:)=rhor_wk(:,1)
1210 41 : call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,1,ngfft,0)
1211 41 : ABI_FREE(work)
1212 :
1213 : ! Compute structure factor phases for new atomic pos:
1214 41 : call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk)
1215 :
1216 : ! Compute local ionic pseudopotential vpsp:
1217 : ! and core electron density xccc3d, if needed.
1218 41 : n3xccc=0;if (n1xccc/=0) n3xccc=nfft
1219 123 : ABI_MALLOC(xccc3d,(n3xccc))
1220 82 : ABI_MALLOC(vpsp_wk,(nfft))
1221 41 : vprtrb(1:2)=zero
1222 :
1223 : ! Determine by which method the local ionic potential and/or
1224 : ! the pseudo core charge density have to be computed
1225 : ! Local ionic potential:
1226 : ! Method 1: PAW
1227 : ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
1228 41 : vloc_method=1;if (psps%usepaw==0) vloc_method=2
1229 41 : if (dtset%icoulomb>0) vloc_method=2
1230 41 : if (psps%usewvl==1) vloc_method=2
1231 : ! Pseudo core charge density:
1232 : ! Method 1: PAW, nc_xccc_gspace
1233 : ! Method 2: Norm-conserving PP, wavelets
1234 : coredens_method=1;if (psps%usepaw==0) coredens_method=2
1235 41 : if (psps%nc_xccc_gspace==1) coredens_method=1
1236 41 : if (psps%nc_xccc_gspace==0) coredens_method=2
1237 41 : if (psps%usewvl==1) coredens_method=2
1238 :
1239 : ! Local ionic potential and/or pseudo core charge by method 1
1240 41 : if (vloc_method==1.or.coredens_method==1) then
1241 0 : optv=0;if (vloc_method==1) optv=1
1242 0 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
1243 0 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1
1244 : ! Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused...
1245 : call atm2fft(atindx,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,gmet,&
1246 : & gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,&
1247 : & nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
1248 : & psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dtset%rcut,dummy_in,rprimd,dummy_out6,dummy_out7,ucvol,&
1249 : & psps%usepaw,dummy_in,dummy_in,dummy_in,vprtrb,psps%vlspl,&
1250 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
1251 0 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
1252 : end if
1253 :
1254 : ! Local ionic potential by method 2
1255 41 : if (vloc_method==2) then
1256 41 : option=1
1257 123 : ABI_MALLOC(dyfrlo_indx,(3,3,dtset%natom))
1258 82 : ABI_MALLOC(grtn_indx,(3,dtset%natom))
1259 : call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,&
1260 : & mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,&
1261 : & ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,&
1262 41 : & ucvol,vprtrb,vpsp_wk,wvl%descr,wvl%den,xred)
1263 41 : ABI_FREE(dyfrlo_indx)
1264 41 : ABI_FREE(grtn_indx)
1265 : end if
1266 :
1267 : ! Pseudo core electron density by method 2
1268 41 : if (coredens_method==2.and.n1xccc/=0) then
1269 41 : option=1
1270 123 : ABI_MALLOC(dyfrx2,(3,3,dtset%natom))
1271 82 : ABI_MALLOC(grxc_indx,(3,dtset%natom))
1272 : call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
1273 : & n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,&
1274 41 : & psps%xccc1d,xccc3d,xred_wk)
1275 41 : ABI_FREE(dyfrx2)
1276 41 : ABI_FREE(grxc_indx)
1277 : end if
1278 :
1279 : ! Compute Hartree+xc potentials
1280 123 : ABI_MALLOC(vxc_wk,(nfft,dtset%nspden))
1281 82 : ABI_MALLOC(vhartr_wk,(nfft))
1282 41 : option=1
1283 :
1284 : call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfft,ngfft,&
1285 41 : &dtset%nkpt,dtset%rcut,rhog_wk,rprimd,dtset%vcutgeo,vhartr_wk)
1286 :
1287 : ! Prepare the call to rhotoxc
1288 41 : call xcdata_init(xcdata,dtset=dtset)
1289 41 : nk3xc=1 ; non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
1290 41 : ABI_MALLOC(work,(0))
1291 : call rhotoxc(bigexc,bigsxc,kxc,mpi_enreg,nfft,ngfft,&
1292 : & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor_wk,rprimd,1,&
1293 41 : & vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk)
1294 41 : ABI_FREE(work)
1295 41 : ABI_FREE(xccc3d)
1296 :
1297 : ! Sum all contributions
1298 123 : do ispden=1,min(dtset%nspden,2)
1299 984123 : do ifft=1,nfft
1300 984082 : vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden)
1301 : end do
1302 : end do
1303 41 : if (dtset%nspden==4) then
1304 0 : do ispden=3,4
1305 0 : do ifft=1,nfft
1306 0 : vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden)
1307 : end do
1308 : end do
1309 : end if
1310 41 : call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1311 41 : if(dtset%nspden==2) then
1312 41 : vmean(1)=half*(vmean(1)+vmean(2))
1313 41 : vmean(2)=vmean(1)
1314 : end if
1315 123 : do ispden=1,dtset%nspden
1316 984123 : vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden)
1317 : end do
1318 41 : ABI_FREE(gresid)
1319 41 : ABI_FREE(grxc)
1320 41 : ABI_FREE(rhog_wk)
1321 41 : ABI_FREE(rhor_wk)
1322 41 : ABI_FREE(rhor_wk0)
1323 41 : ABI_FREE(xred_wk)
1324 41 : ABI_FREE(vhartr_wk)
1325 41 : ABI_FREE(vpsp_wk)
1326 82 : ABI_FREE(vxc_wk)
1327 :
1328 : end if
1329 :
1330 28280 : end subroutine prcref_PMA
1331 : !!***
1332 :
1333 :
1334 : !!****f* ABINIT/moddiel
1335 : !! NAME
1336 : !! moddiel
1337 : !!
1338 : !! FUNCTION
1339 : !! Precondition the residual, using a model dielectric function.
1340 : !! When cplex=1, assume q=(0 0 0), and vresid and vrespc will be REAL
1341 : !! When cplex=2, q must be taken into account, and vresid and vrespc will be COMPLEX
1342 : !!
1343 : !! INPUTS
1344 : !! cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX
1345 : !! dielar(7)=input parameters for dielectric matrix:
1346 : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
1347 : !! mpi_enreg=information about MPI parallelization
1348 : !! nfft=(effective) number of FFT grid points (for this processor)
1349 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1350 : !! nspden=number of spin-density components
1351 : !! optreal=1 if residual potential is in REAL space, 2 if it is in RECIPROCAL SPACE
1352 : !! optres=0: the array vresid contains a potential residual
1353 : !! 1: the array vresid contains a density residual
1354 : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
1355 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1356 : !! vresid(cplex*nfft,nspden)=residual density/potential in REAL space (if optreal==1)
1357 : !! residual density/potential in RECIPROCAL space (if optreal==2)
1358 : !!
1359 : !! OUTPUT
1360 : !! vrespc(cplex*nfft,nspden)=preconditioned residual of the density/potential
1361 : !! in REAL space if optreal==1
1362 : !! in RECIPROCAL space if optreal==2
1363 : !!
1364 : !! SIDE EFFECTS
1365 : !!
1366 : !! NOTES
1367 : !! optreal==2 is not compatible with cplex==1
1368 : !!
1369 : !! SOURCE
1370 :
1371 54368 : subroutine moddiel(cplex,dielar,mpi_enreg,nfft,ngfft,nspden,optreal,optres,qphon,rprimd,vresid,vrespc)
1372 :
1373 : !Arguments-------------------------------
1374 : !scalars
1375 : integer,intent(in) :: cplex,nfft,nspden,optreal,optres
1376 : type(MPI_type),intent(in) :: mpi_enreg
1377 : !arrays
1378 : integer,intent(in) :: ngfft(18)
1379 : real(dp),intent(in) :: dielar(7),qphon(3),rprimd(3,3)
1380 : real(dp),intent(in) :: vresid(cplex*nfft,nspden)
1381 : real(dp),intent(out) :: vrespc(cplex*nfft,nspden)
1382 :
1383 : !Local variables-------------------------------
1384 : !scalars
1385 : integer,parameter :: im=2,re=1
1386 : integer :: i1,i2,i23,i3,ifft,ig,ii,ii1,ing,ispden,me_fft,mg,n1,n2,n3,nproc_fft
1387 : integer :: nspden_eff,qeq0
1388 : logical :: magn_precon
1389 : real(dp) :: dielng,diemac,diemac_inv,diemix,diemixmag,diemix_eff,factor,gqg2p3,gqgm12,gqgm13
1390 : real(dp) :: gqgm23,gs,gs2,gs3,l2g2,length2,ucvol
1391 : character(len=500) :: message
1392 : !arrays
1393 : integer :: id(3)
1394 54368 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1395 54368 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1396 : real(dp) :: gmet(3,3),gprimd(3,3),potg0(4),rmet(3,3)
1397 54368 : real(dp),allocatable :: gq(:,:),work1(:,:),work2(:)
1398 : ! *************************************************************************
1399 :
1400 : !Check that cplex has an allowed value
1401 54368 : if(cplex/=1 .and. cplex/=2)then
1402 : write(message,'(a,i0,a,a)')&
1403 0 : & ' From the calling routine, cplex=',cplex,ch10,&
1404 0 : & ' but the only value allowed are 1 and 2.'
1405 0 : ABI_BUG(message)
1406 : end if
1407 :
1408 54368 : if(cplex==1.and.optreal==2)then
1409 0 : ABI_BUG('When optreal=2, cplex must be 2.')
1410 : end if
1411 :
1412 : !This is to allow q=0
1413 54368 : qeq0=0
1414 54368 : if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
1415 :
1416 : !If cplex=1 then qphon should be 0 0 0
1417 54368 : if (cplex==1.and. qeq0/=1) then
1418 0 : write(message,'(a,3e12.4,a)' )' cplex=1 but qphon=',qphon,' qphon should be 0 0 0.'
1419 0 : ABI_BUG(message)
1420 : end if
1421 :
1422 54368 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1423 54368 : me_fft=ngfft(11)
1424 54368 : nproc_fft=ngfft(10)
1425 :
1426 : !Get the distrib associated with this fft_grid
1427 54368 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1428 :
1429 : !Compute different geometric tensor, as well as ucvol, from rprimd
1430 54368 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1431 :
1432 54368 : dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ; diemixmag=dielar(7)
1433 :
1434 54368 : magn_precon=(diemixmag>=zero) ! Set to true if magnetization has to be preconditionned
1435 54368 : diemixmag=abs(diemixmag)
1436 :
1437 : !write(std_out,*)' moddiel : diemac, diemix, diemixmag =',diemac,diemix,diemixmag
1438 :
1439 54368 : if(abs(diemac-1.0_dp)<1.0d-6)then
1440 :
1441 : ! Here, simple mixing is required, through macroscopic
1442 : ! dielectric constant set to 1.0_dp .
1443 149078034 : vrespc(:,1)=diemix*vresid(:,1)
1444 52532569 : if (nspden/=1) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden)
1445 : else
1446 :
1447 : ! Magnetization is not preconditionned
1448 1752658 : if (optres==1.and.nspden>1.and.(.not.magn_precon)) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden)
1449 :
1450 : ! Here, model dielectric function (G-diagonal operator)
1451 :
1452 47550 : length2=(two_pi*dielng)**2
1453 47550 : diemac_inv=1.0_dp/diemac
1454 142650 : ABI_MALLOC(work1,(2,nfft))
1455 47550 : if (optreal==1) then
1456 114430 : ABI_MALLOC(work2,(cplex*nfft))
1457 : end if
1458 :
1459 : ! In order to speed the routine, precompute the components of g
1460 903450 : mg=maxval(ngfft)
1461 142650 : ABI_MALLOC(gq,(3,mg))
1462 190200 : do ii=1,3
1463 142650 : id(ii)=ngfft(ii)/2+2
1464 3000467 : do ing=1,ngfft(ii)
1465 2810267 : ig=ing-(ing/id(ii))*ngfft(ii)-1
1466 2952917 : gq(ii,ing)=ig+qphon(ii)
1467 : end do
1468 : end do
1469 :
1470 : ! Do-loop on spins
1471 : ! Note XG 010922 : I doubt the preconditioner is OK for the magnetization
1472 47550 : nspden_eff=nspden;if (optres==1.and.(.not.magn_precon)) nspden_eff=1
1473 106525 : do ispden=1,nspden_eff
1474 :
1475 58975 : diemix_eff=diemix;if (ispden>1) diemix_eff=diemixmag
1476 :
1477 : ! Do fft from real space (work2) to G space (work1)
1478 58975 : if (optreal==1) then
1479 565062320 : work2(:)=vresid(:,ispden)
1480 40048 : call fourdp(cplex,work1,work2,-1,mpi_enreg,nfft,1,ngfft,0)
1481 : else
1482 : ! work1(:,:)=reshape(vresid(:,ispden),(/2,nfft/))
1483 : ! Reshape function does not work with big arrays for some compilers
1484 234145414 : do ifft=1,nfft
1485 234126487 : work1(1,ifft)=vresid(2*ifft-1,ispden)
1486 234145414 : work1(2,ifft)=vresid(2*ifft ,ispden)
1487 : end do
1488 : end if
1489 :
1490 : ! Store G=0 value
1491 58975 : potg0(ispden)=work1(re,1)
1492 :
1493 : ! Triple loop, for the three dimensions
1494 1274247 : do i3=1,n3
1495 : ! Precompute some products that do not depend on i2 and i1
1496 1215272 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
1497 1215272 : gqgm23=gq(3,i3)*gmet(2,3)*2
1498 1215272 : gqgm13=gq(3,i3)*gmet(1,3)*2
1499 28775584 : do i2=1,n2
1500 28716609 : if (fftn2_distrib(i2)==me_fft) then
1501 26390001 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
1502 26390001 : gqgm12=gq(2,i2)*gmet(1,2)*2
1503 26390001 : gqg2p3=gqgm13+gqgm12
1504 26390001 : i23=n1*( ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
1505 :
1506 : ! Do the test that eliminates the Gamma point outside
1507 : ! of the inner loop
1508 26390001 : ii1=1
1509 26390001 : if(i2 == 1 .and. i3 == 1 .and. qeq0==1)then
1510 : ! if(i23==0 .and. qeq0==1)then: this changes with the number of fft procs...
1511 : ! and seems to be wrong.Pls check
1512 26390001 : ii1=2
1513 : end if
1514 :
1515 : ! Here, unlike in hartre.f, the G=0 term is not eliminated, albeit
1516 : ! not changed.
1517 796117921 : do i1=ii1,n1
1518 :
1519 : ! One obtains the square of the norm of q+G (defined by i1,i2,i3)
1520 769727920 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
1521 769727920 : ifft=i1+i23
1522 :
1523 769727920 : l2g2=length2*gs
1524 : ! The model dielectric function is now computed
1525 769727920 : factor = (l2g2+diemac_inv)/(l2g2+1.0_dp) * diemix_eff
1526 769727920 : work1(re,ifft)=work1(re,ifft)*factor
1527 796117921 : work1(im,ifft)=work1(im,ifft)*factor
1528 :
1529 : end do
1530 : end if
1531 : end do
1532 : end do
1533 :
1534 : ! Might get rid of the G=0 term
1535 : ! if(qeq0==1)then
1536 : ! work1(re,1)=0.0_dp
1537 : ! work1(im,1)=0.0_dp
1538 : ! end if
1539 :
1540 : ! Fourier transform
1541 106525 : if (optreal==1) then
1542 40048 : call fourdp(cplex,work1,work2,1,mpi_enreg,nfft,1,ngfft,0)
1543 565062320 : vrespc(:,ispden)=work2(:)
1544 : else
1545 : ! vrespc(:,ispden)=reshape(work1(:,:),(/nfft*2/))
1546 : ! Reshape function does not work with big arrays for some compilers
1547 234145414 : do ifft=1,nfft
1548 234126487 : vrespc(2*ifft-1,ispden)=work1(1,ifft)
1549 234145414 : vrespc(2*ifft ,ispden)=work1(2,ifft)
1550 : end do
1551 : end if
1552 :
1553 : ! End of loop on spin polarizations
1554 : end do
1555 :
1556 47550 : ABI_FREE(gq)
1557 47550 : ABI_FREE(work1)
1558 47550 : if (optreal==1) then
1559 33440 : ABI_FREE(work2)
1560 : end if
1561 :
1562 : ! End condition diemac/=1.0
1563 : end if
1564 :
1565 108736 : end subroutine moddiel
1566 : !!***
1567 :
1568 : !!****f* ABINIT/dielmt
1569 : !! NAME
1570 : !! dielmt
1571 : !!
1572 : !! FUNCTION
1573 : !! Compute dielectric matrix from susceptibility matrix
1574 : !! Diagonalize it, then invert it.
1575 : !!
1576 : !! INPUTS
1577 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
1578 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
1579 : !! npwdiel=size of the dielinv and susmat arrays.
1580 : !! nspden=number of spin-density components
1581 : !! occopt=option for occupancies
1582 : !! prtvol=control print volume and debugging output
1583 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
1584 : !! the susceptibility (or density-density response) matrix in reciprocal space
1585 : !!
1586 : !! OUTPUT
1587 : !! dielinv(2,npwdiel,(nspden+4)/3,npwdiel,(nspden+4)/3)=inverse of the (non-hermitian)
1588 : !! TC dielectric matrix in reciprocal space.
1589 : !!
1590 : !! NOTES
1591 : !! Warning : will not work in the spin-polarized, metallic case.
1592 : !! Output (not cleaned)
1593 : !! !!! Spin behaviour is not obvious !!!
1594 : !!
1595 : !! TODO
1596 : !! Write equation below (hermitian matrix)
1597 : !!
1598 : !! SOURCE
1599 :
1600 68 : subroutine dielmt(dielinv,gmet,kg_diel,npwdiel,nspden,occopt,prtvol,susmat)
1601 :
1602 : !Arguments ------------------------------------
1603 : !scalars
1604 : integer,intent(in) :: npwdiel,nspden,occopt,prtvol
1605 : !arrays
1606 : integer,intent(in) :: kg_diel(3,npwdiel)
1607 : real(dp),intent(in) :: gmet(3,3)
1608 : real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
1609 : real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden)
1610 :
1611 : !Local variables-------------------------------
1612 : !scalars
1613 : integer :: ieig,ier,ii,index,ipw,ipw1,ipw2,isp,jj,npwsp
1614 : real(dp) :: ai1,ai2,ar1,ar2,eiginv,gfact,gfactinv,kg_red1,kg_red2,kg_red3,gsquar
1615 : real(dp) :: tpisq
1616 : character(len=500) :: message
1617 : !arrays
1618 : real(dp) :: tsec(2)
1619 68 : real(dp),allocatable :: dielh(:),dielmat(:,:,:,:,:),dielvec(:,:,:)
1620 68 : real(dp),allocatable :: eig_diel(:),zhpev1(:,:),zhpev2(:)
1621 : !integer :: ipw3
1622 : !real(dp) :: elementi,elementr
1623 : ! *************************************************************************
1624 :
1625 : !DEBUG
1626 : !write(std_out,*)' dielmt : enter '
1627 : !ENDDEBUG
1628 :
1629 : !tpisq is (2 Pi) **2:
1630 68 : tpisq=(two_pi)**2
1631 :
1632 68 : call timab(90,1,tsec)
1633 :
1634 : !-Compute now the hermitian dielectric matrix------------------------------
1635 : !Following remarks are only valid within RPA approximation (Kxc=0):
1636 :
1637 : !for the spin-unpolarized case, 1 - 4pi (1/G) chi0(G,Gp) (1/Gp)
1638 :
1639 : !for the spin-polarized case,
1640 : !( 1 0 ) - 4pi ( 1/G 1/G ) ( chi0 upup chi0 updn ) ( 1/Gp 1/Gp )
1641 : !( 0 1 ) ( 1/G 1/G ) ( chi0 dnup chi0 dndn ) ( 1/Gp 1/Gp )
1642 : !which is equal to
1643 : !( 1 0 ) - 4pi (1/G 0 ) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) (1/Gp 0 )
1644 : !( 0 1 ) ( 0 1/G) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) ( 0 1/Gp)
1645 : !So, if spin-polarized, sum all spin contributions
1646 : !Note: chi0 updn = chi0 dnup = zero for non-metallic systems
1647 :
1648 : !In the case of non-collinear magnetism, within RPA, this is the same because:
1649 : !chi0_(s1,s2),(s3,s4) = delta_s1,s2 * delta_s3,s4 * chi0_(s1,s1),(s3,s3)
1650 : !Only chi_upup,upup, chi_dndn,dndn, chi_upup,dndn and chi_dndn,upup
1651 : !have to be taken into account (stored, susmat(:,ipw1,1:2,ipw2,1:2)
1652 :
1653 408 : ABI_MALLOC(dielmat,(2,npwdiel,min(nspden,2),npwdiel,min(nspden,2)))
1654 :
1655 68 : if(nspden/=1)then
1656 2 : if (occopt<3) then
1657 120 : do ipw2=1,npwdiel
1658 7082 : do ipw1=1,npwdiel
1659 6962 : dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2)
1660 7080 : dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2)
1661 : end do
1662 : end do
1663 : else
1664 0 : do ipw2=1,npwdiel
1665 0 : do ipw1=1,npwdiel
1666 0 : dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2)+susmat(1,ipw1,1,ipw2,2)+susmat(1,ipw1,2,ipw2,1)
1667 0 : dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2)+susmat(2,ipw1,1,ipw2,2)+susmat(2,ipw1,2,ipw2,1)
1668 : end do
1669 : end do
1670 : end if
1671 : else
1672 11080 : do ipw2=1,npwdiel
1673 2602266 : do ipw1=1,npwdiel
1674 2591186 : dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)
1675 2602200 : dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)
1676 : end do
1677 : end do
1678 : end if
1679 : !Compute 1/G factors and include them in the dielectric matrix
1680 11200 : do ipw1=1,npwdiel
1681 11132 : kg_red1=dble(kg_diel(1,ipw1))
1682 11132 : kg_red2=dble(kg_diel(2,ipw1))
1683 11132 : kg_red3=dble(kg_diel(3,ipw1))
1684 : gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 &
1685 : & +two*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + &
1686 11132 : & gmet(2,3)*kg_red2*kg_red3) )
1687 : ! Distinguish G=0 from other elements
1688 11200 : if(gsquar>tol12)then
1689 : ! !$ gfact=\sqrt (4.0_dp \pi/gsquar/dble(nspden))$
1690 11064 : gfact=sqrt(four_pi/gsquar)
1691 2598080 : do ipw2=1,npwdiel
1692 : ! Must multiply both rows and columns, and also changes the sign
1693 2587016 : dielmat(1,ipw2,1,ipw1,1)=-dielmat(1,ipw2,1,ipw1,1)*gfact
1694 2587016 : dielmat(2,ipw2,1,ipw1,1)=-dielmat(2,ipw2,1,ipw1,1)*gfact
1695 2587016 : dielmat(1,ipw1,1,ipw2,1)= dielmat(1,ipw1,1,ipw2,1)*gfact
1696 2598080 : dielmat(2,ipw1,1,ipw2,1)= dielmat(2,ipw1,1,ipw2,1)*gfact
1697 : end do
1698 : else
1699 : ! Zero the G=0 elements, head and wings
1700 11200 : do ipw2=1,npwdiel
1701 11132 : dielmat(1,ipw2,1,ipw1,1)=zero
1702 11132 : dielmat(2,ipw2,1,ipw1,1)=zero
1703 11132 : dielmat(1,ipw1,1,ipw2,1)=zero
1704 11200 : dielmat(2,ipw1,1,ipw2,1)=zero
1705 : end do
1706 : end if
1707 : end do
1708 :
1709 : !Complete the matrix in the spin-polarized case
1710 : !should this be nspden==2??
1711 68 : if(nspden/=1)then
1712 120 : do ipw1=1,npwdiel
1713 7082 : do ipw2=1,npwdiel
1714 6962 : dielmat(1,ipw1,1,ipw2,2)=dielmat(1,ipw1,1,ipw2,1)
1715 6962 : dielmat(2,ipw1,1,ipw2,2)=dielmat(2,ipw1,1,ipw2,1)
1716 6962 : dielmat(1,ipw1,2,ipw2,1)=dielmat(1,ipw1,1,ipw2,1)
1717 6962 : dielmat(2,ipw1,2,ipw2,1)=dielmat(2,ipw1,1,ipw2,1)
1718 6962 : dielmat(1,ipw1,2,ipw2,2)=dielmat(1,ipw1,1,ipw2,1)
1719 7080 : dielmat(2,ipw1,2,ipw2,2)=dielmat(2,ipw1,1,ipw2,1)
1720 : end do
1721 : end do
1722 : end if
1723 :
1724 : !DEBUG
1725 : !write(std_out,*)' dielmt : make dielmat equal to identity matrix '
1726 : !do ipw1=1,npwdiel
1727 : !do ipw2=1,npwdiel
1728 : !dielmat(1,ipw1,1,ipw2,1)=0.0_dp
1729 : !dielmat(2,ipw1,1,ipw2,1)=0.0_dp
1730 : !end do
1731 : !end do
1732 : !ENDDEBUG
1733 :
1734 : !Add the diagonal part
1735 138 : do isp=1,min(nspden,2)
1736 11388 : do ipw=1,npwdiel
1737 11320 : dielmat(1,ipw,isp,ipw,isp)=one+dielmat(1,ipw,isp,ipw,isp)
1738 : end do
1739 : end do
1740 :
1741 : !-The hermitian dielectric matrix is computed ------------------------------
1742 : !-Now, diagonalize it ------------------------------------------------------
1743 :
1744 : !In RPA, everything is projected on the spin-symmetrized
1745 : !space. This was coded here (for the time being).
1746 :
1747 : !Diagonalize the hermitian dielectric matrix
1748 :
1749 : !npwsp=npwdiel*nspden
1750 68 : npwsp=npwdiel
1751 :
1752 204 : ABI_MALLOC(dielh,(npwsp*(npwsp+1)))
1753 272 : ABI_MALLOC(dielvec,(2,npwsp,npwsp))
1754 204 : ABI_MALLOC(eig_diel,(npwsp))
1755 204 : ABI_MALLOC(zhpev1,(2,2*npwsp-1))
1756 204 : ABI_MALLOC(zhpev2,(3*npwsp-2))
1757 68 : ier=0
1758 : !Store the dielectric matrix in proper mode before calling zhpev
1759 68 : index=1
1760 11200 : do ii=1,npwdiel
1761 1315840 : do jj=1,ii
1762 1304640 : dielh(index )=dielmat(1,jj,1,ii,1)
1763 1304640 : dielh(index+1)=dielmat(2,jj,1,ii,1)
1764 1315772 : index=index+2
1765 : end do
1766 : end do
1767 : !If spin-polarized and non RPA, need to store other parts of the matrix
1768 : !if(nspden/=1)then
1769 : !do ii=1,npwdiel
1770 : !Here, spin-flip contribution
1771 : !do jj=1,npwdiel
1772 : !dielh(index )=dielmat(1,jj,1,ii,2)
1773 : !dielh(index+1)=dielmat(2,jj,1,ii,2)
1774 : !index=index+2
1775 : !end do
1776 : !Here spin down-spin down upper matrix
1777 : !do jj=1,ii
1778 : !dielh(index )=dielmat(1,jj,2,ii,2)
1779 : !dielh(index+1)=dielmat(2,jj,2,ii,2)
1780 : !index=index+2
1781 : !end do
1782 : !end do
1783 : !end if
1784 :
1785 : call ZHPEV ('V','U',npwsp,dielh,eig_diel,dielvec,npwdiel,zhpev1,&
1786 68 : & zhpev2,ier)
1787 68 : ABI_FREE(zhpev1)
1788 68 : ABI_FREE(zhpev2)
1789 :
1790 68 : if(prtvol>=10)then
1791 15 : write(message, '(a,a,a,5es12.4)' )ch10,&
1792 15 : & ' Five largest eigenvalues of the hermitian RPA dielectric matrix:',&
1793 30 : & ch10,eig_diel(npwdiel:npwdiel-4:-1)
1794 15 : call wrtout(ab_out,message,'COLL')
1795 : end if
1796 :
1797 68 : write(message, '(a,a)' )ch10,&
1798 136 : & ' dielmt : 15 largest eigenvalues of the hermitian RPA dielectric matrix'
1799 68 : call wrtout(std_out,message,'COLL')
1800 68 : write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(npwdiel:npwdiel-4:-1)
1801 68 : call wrtout(std_out,message,'COLL')
1802 68 : write(message, '(a,5es12.5)' )' 6-10 :',eig_diel(npwdiel-5:npwdiel-9:-1)
1803 68 : call wrtout(std_out,message,'COLL')
1804 68 : write(message, '(a,5es12.5)' )' 11-15:',eig_diel(npwdiel-10:npwdiel-14:-1)
1805 68 : call wrtout(std_out,message,'COLL')
1806 68 : write(message, '(a,a)' )ch10,&
1807 136 : & ' dielmt : 5 smallest eigenvalues of the hermitian RPA dielectric matrix'
1808 68 : call wrtout(std_out,message,'COLL')
1809 68 : write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(1:5)
1810 68 : call wrtout(std_out,message,'COLL')
1811 :
1812 : !Invert the hermitian dielectric matrix,
1813 : !Should use a BLAS call !
1814 11200 : do ipw2=1,npwdiel
1815 1315840 : do ipw1=ipw2,npwdiel
1816 1304640 : dielinv(1,ipw1,1,ipw2,1)=zero
1817 1315772 : dielinv(2,ipw1,1,ipw2,1)=zero
1818 : end do
1819 : end do
1820 11200 : do ieig=1,npwdiel
1821 11132 : eiginv=one/eig_diel(ieig)
1822 2609348 : do ipw2=1,npwdiel
1823 384025872 : do ipw1=ipw2,npwdiel
1824 381416592 : ar1=dielvec(1,ipw1,ieig)
1825 381416592 : ai1=dielvec(2,ipw1,ieig)
1826 381416592 : ar2=dielvec(1,ipw2,ieig)
1827 381416592 : ai2=dielvec(2,ipw2,ieig)
1828 : dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)+&
1829 381416592 : & (ar1*ar2+ai1*ai2)*eiginv
1830 : dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)+&
1831 384014740 : & (ai1*ar2-ar1*ai2)*eiginv
1832 : end do
1833 : end do
1834 : end do
1835 11132 : do ipw2=1,npwdiel-1
1836 1304640 : do ipw1=ipw2+1,npwdiel
1837 1293508 : dielinv(1,ipw2,1,ipw1,1)= dielinv(1,ipw1,1,ipw2,1)
1838 1304572 : dielinv(2,ipw2,1,ipw1,1)=-dielinv(2,ipw1,1,ipw2,1)
1839 : end do
1840 : end do
1841 :
1842 68 : ABI_FREE(dielh)
1843 68 : ABI_FREE(dielvec)
1844 68 : ABI_FREE(eig_diel)
1845 :
1846 : !DEBUG
1847 : !Checks whether the inverse of the hermitian dielectric matrix
1848 : !has been correctly generated
1849 : !do ipw1=1,npwdiel
1850 : !do ipw2=1,npwdiel
1851 : !elementr=0.0_dp
1852 : !elementi=0.0_dp
1853 : !do ipw3=1,npwdiel
1854 : !elementr=elementr+dielinv(1,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1)&
1855 : !& -dielinv(2,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1)
1856 : !elementi=elementi+dielinv(1,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1)&
1857 : !& +dielinv(2,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1)
1858 : !end do
1859 : !if(elementr**2+elementi**2 > 1.0d-12)then
1860 : !if( ipw1 /= ipw2 .or. &
1861 : !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then
1862 : !write(std_out,*)' dielmt : the inversion procedure is not correct '
1863 : !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2
1864 : !write(std_out,*)' elementr,elementi=',elementr,elementi
1865 : !stop
1866 : !end if
1867 : !end if
1868 : !end do
1869 : !end do
1870 : !write(std_out,*)'dielmt : matrix has been inverted successfully '
1871 : !stop
1872 : !ENDDEBUG
1873 :
1874 : !Then get the inverse of the asymmetric
1875 : !dielectric matrix, as required for the preconditioning.
1876 :
1877 : !Inverse of the dielectric matrix : ( 1 - 4pi (1/G^2) chi0(G,Gp) )^(-1)
1878 : !In dielinv there is now (1 - 4pi (1/G) chi0(G,Gp) (1/Gp) )^(-1)
1879 : !So, evaluate dielinv_after(G,Gp) =
1880 : !(4pi/G^2)^(1/2) dielinv_before(G,Gp) (4pi/Gp^2)^(-1/2)
1881 : !In RPA, can focus on the spin-averaged quantities
1882 11200 : do ipw1=1,npwdiel
1883 11132 : kg_red1=dble(kg_diel(1,ipw1))
1884 11132 : kg_red2=dble(kg_diel(2,ipw1))
1885 11132 : kg_red3=dble(kg_diel(3,ipw1))
1886 : gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 &
1887 : & +two*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + &
1888 11132 : & gmet(2,3)*kg_red2*kg_red3) )
1889 : ! Distinguish G=0 from other elements
1890 11200 : if(gsquar>tol12)then
1891 11064 : gfact=sqrt(four_pi/gsquar)
1892 11064 : gfactinv=one/gfact
1893 2598080 : do ipw2=1,npwdiel
1894 : ! Must multiply both rows and columns
1895 2587016 : dielinv(1,ipw2,1,ipw1,1)=dielinv(1,ipw2,1,ipw1,1)*gfactinv
1896 2587016 : dielinv(2,ipw2,1,ipw1,1)=dielinv(2,ipw2,1,ipw1,1)*gfactinv
1897 2587016 : dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)*gfact
1898 2598080 : dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)*gfact
1899 : end do
1900 : else
1901 : ! Zero the G=0 elements, head
1902 11200 : do ipw2=1,npwdiel
1903 33328 : if (ipw2/=ipw1) dielinv(1:2,ipw1,1,ipw2,1)=zero
1904 : end do
1905 : end if
1906 : end do
1907 :
1908 68 : ABI_FREE(dielmat)
1909 :
1910 68 : call timab(90,2,tsec)
1911 :
1912 68 : end subroutine dielmt
1913 : !!***
1914 :
1915 :
1916 : !!****f* ABINIT/dieltcel
1917 : !! NAME
1918 : !! dieltcel
1919 : !!
1920 : !! FUNCTION
1921 : !! Compute either test charge or electronic dielectric matrices
1922 : !! from susceptibility matrix
1923 : !! Diagonalize it, then invert it.
1924 : !!
1925 : !! INPUTS
1926 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
1927 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
1928 : !! kxc(nfft,nkxc)=exchange-correlation kernel,
1929 : !! needed if the electronic dielectric matrix is computed
1930 : !! nfft=(effective) number of FFT grid points (for this processor)
1931 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1932 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
1933 : !! npwdiel=size of the dielinv and susmat arrays.
1934 : !! nspden=number of spin-density components
1935 : !! occopt=option for occupancies
1936 : !! option=1 for Test Charge dielectric matrix, 2 for electronic dielectric matrix
1937 : !! prtvol=control print volume and debugging output
1938 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
1939 : !! the susceptibility (or density-density response) matrix in reciprocal space
1940 : !!
1941 : !! OUTPUT
1942 : !! dielinv(2,npwdiel,nspden,npwdiel,nspden)=inverse of the (non-hermitian)
1943 : !! TC dielectric matrix in reciprocal space.
1944 : !!
1945 : !! NOTES
1946 : !! Output (not cleaned)
1947 : !! !!! Spin behaviour is not obvious !!!
1948 : !! Will not work in the spin-polarized, metallic case.
1949 : !!
1950 : !! SOURCE
1951 :
1952 6 : subroutine dieltcel(dielinv,gmet,kg_diel,kxc,nfft,ngfft,nkxc,npwdiel,nspden,occopt,option,prtvol,susmat)
1953 :
1954 : !Arguments ------------------------------------
1955 : !scalars
1956 : integer,intent(in) :: nfft,nkxc,npwdiel,nspden,occopt,option
1957 : integer,intent(in) :: prtvol
1958 : !arrays
1959 : integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18)
1960 : real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc)
1961 : real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
1962 : real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden)
1963 :
1964 : !Local variables-------------------------------
1965 : !scalars
1966 : integer :: i1,i2,i3,ieig,ier,ifft,ii,index,ipw0,ipw1,ipw2,ispden,j1
1967 : integer :: j2,j3,jj,k1,k2,k3,n1,n2,n3
1968 : real(dp) :: ai,ai2,ar,ar2,eiginv,kg_red1,kg_red2,kg_red3,gsquar,si
1969 : real(dp) :: sr,tpisq
1970 : character(len=500) :: message
1971 6 : type(MPI_type) :: mpi_enreg_seq
1972 : !arrays
1973 : real(dp) :: tsec(2)
1974 6 : real(dp),allocatable :: eig_msusinvsqr(:),eig_msussqr(:)
1975 6 : real(dp),allocatable :: eig_sus(:),eig_sym(:),invsqrsus(:,:,:,:,:)
1976 6 : real(dp),allocatable :: khxc(:,:,:,:,:),kxcg(:,:),sqrsus(:,:,:,:,:),sush(:)
1977 6 : real(dp),allocatable :: susvec(:,:,:),symdielmat(:,:,:,:,:),symh(:)
1978 6 : real(dp),allocatable :: symvec(:,:,:,:,:),wkxc(:),work(:,:,:,:,:)
1979 6 : real(dp),allocatable :: work2(:,:,:,:,:),zhpev1(:,:),zhpev2(:)
1980 : !no_abirules
1981 : !integer :: ipw3
1982 : !real(dp) :: elementi,elementr
1983 : !DEBUG
1984 : !Used to moderate divergence effect near rho=0 of the Kxc
1985 : !this limit value is truly empirical (exprmt on small Sr cell).
1986 : !real(dp) :: kxc_min=-200.0
1987 : !ENDDEBUG
1988 : ! *************************************************************************
1989 :
1990 6 : call timab(96,1,tsec)
1991 :
1992 : !tpisq is (2 Pi) **2:
1993 6 : tpisq=(two_pi)**2
1994 :
1995 6 : if(nspden/=1 .and. (occopt>=3 .and. occopt<=8) )then
1996 : write(message, '(a,a,a)' )&
1997 0 : & 'In the present version of the code, one cannot produce',ch10,&
1998 0 : & 'the dielectric matrix in the metallic, spin-polarized case.'
1999 0 : ABI_BUG(message)
2000 : end if
2001 :
2002 6 : if(nspden==4)then
2003 : write(message,'(a,a,a)')&
2004 0 : & 'In the present version of the code, one cannot produce',ch10,&
2005 0 : & 'the dielectric matrix in the non-collinear spin-polarized case.'
2006 0 : ABI_ERROR(message)
2007 : end if
2008 :
2009 :
2010 : !-Diagonalize the susceptibility matrix
2011 :
2012 18 : ABI_MALLOC(sush,(npwdiel*(npwdiel+1)))
2013 24 : ABI_MALLOC(susvec,(2,npwdiel,npwdiel))
2014 18 : ABI_MALLOC(eig_msusinvsqr,(npwdiel))
2015 12 : ABI_MALLOC(eig_msussqr,(npwdiel))
2016 12 : ABI_MALLOC(eig_sus,(npwdiel))
2017 18 : ABI_MALLOC(zhpev1,(2,2*npwdiel-1))
2018 18 : ABI_MALLOC(zhpev2,(3*npwdiel-2))
2019 36 : ABI_MALLOC(work,(2,npwdiel,nspden,npwdiel,nspden))
2020 30 : ABI_MALLOC(work2,(2,npwdiel,nspden,npwdiel,nspden))
2021 30 : ABI_MALLOC(sqrsus,(2,npwdiel,nspden,npwdiel,nspden))
2022 30 : ABI_MALLOC(invsqrsus,(2,npwdiel,nspden,npwdiel,nspden))
2023 :
2024 : !At some time, should take care of different spin channels
2025 12 : do ispden=1,nspden
2026 :
2027 6 : if(nspden/=1)then
2028 0 : ABI_ERROR('dieltcel : stop, nspden/=1')
2029 : end if
2030 :
2031 : ! Store the susceptibility matrix in proper mode before calling zhpev
2032 : index=1
2033 240 : do ii=1,npwdiel
2034 4920 : do jj=1,ii
2035 4680 : sush(index )=susmat(1,jj,1,ii,1)
2036 4680 : sush(index+1)=susmat(2,jj,1,ii,1)
2037 4914 : index=index+2
2038 : end do
2039 : end do
2040 :
2041 6 : ier=0
2042 6 : call ZHPEV ('V','U',npwdiel,sush,eig_sus,susvec,npwdiel,zhpev1,zhpev2,ier)
2043 :
2044 : ! DEBUG
2045 : ! write(std_out,*)' dieltcel : print eigenvalues of the susceptibility matrix'
2046 : ! do ii=1,npwdiel
2047 : ! write(std_out,'(i5,es16.6)' )ii,eig_sus(ii)
2048 : ! end do
2049 : ! ENDDEBUG
2050 :
2051 240 : do ii=1,npwdiel
2052 240 : if(-eig_sus(ii)>1.d-12)then
2053 228 : eig_msussqr(ii)=sqrt(-eig_sus(ii))
2054 228 : eig_msusinvsqr(ii)=1._dp/eig_msussqr(ii)
2055 6 : else if(-eig_sus(ii)< -1.d-12)then
2056 0 : message = "Found positive eigenvalue of susceptibility matrix."
2057 0 : ABI_BUG(message)
2058 : else
2059 : ! Set the eigenvalue corresponding to a constant potential change to 1,
2060 : ! while it will be set to zero in Khx.
2061 6 : eig_msussqr(ii)=1._dp
2062 6 : eig_msusinvsqr(ii)=1._dp
2063 : end if
2064 : end do
2065 :
2066 : ! Compute square root of minus susceptibility matrix
2067 : ! and inverse square root of minus susceptibility matrix
2068 240 : do ii=1,npwdiel
2069 27612 : work(:,:,1,ii,1)=susvec(:,:,ii)*eig_msussqr(ii)
2070 27618 : work2(:,:,1,ii,1)=susvec(:,:,ii)*eig_msusinvsqr(ii)
2071 : end do
2072 246 : do ipw2=1,npwdiel
2073 4920 : do ipw1=ipw2,npwdiel
2074 : ar=0._dp ; ai=0._dp ; ar2=0._dp ; ai2=0._dp
2075 187200 : do ii=1,npwdiel
2076 182520 : sr=susvec(1,ipw2,ii) ; si=susvec(2,ipw2,ii)
2077 182520 : ar =ar +work(1,ipw1,1,ii,1)*sr +work(2,ipw1,1,ii,1)*si
2078 182520 : ai =ai +work(2,ipw1,1,ii,1)*sr -work(1,ipw1,1,ii,1)*si
2079 182520 : ar2=ar2 +work2(1,ipw1,1,ii,1)*sr +work2(2,ipw1,1,ii,1)*si
2080 187200 : ai2=ai2 +work2(2,ipw1,1,ii,1)*sr -work2(1,ipw1,1,ii,1)*si
2081 : end do
2082 4680 : sqrsus(1,ipw1,1,ipw2,1)=ar
2083 4680 : sqrsus(2,ipw1,1,ipw2,1)=ai
2084 4680 : invsqrsus(1,ipw1,1,ipw2,1)=ar2
2085 4680 : invsqrsus(2,ipw1,1,ipw2,1)=ai2
2086 4914 : if(ipw1/=ipw2)then
2087 4446 : sqrsus(1,ipw2,1,ipw1,1)=ar
2088 4446 : sqrsus(2,ipw2,1,ipw1,1)=-ai
2089 4446 : invsqrsus(1,ipw2,1,ipw1,1)=ar2
2090 4446 : invsqrsus(2,ipw2,1,ipw1,1)=-ai2
2091 : end if
2092 : end do
2093 : end do
2094 :
2095 : ! DEBUG
2096 : ! Checks whether sqrsus and invsqrsus are inverse of each other.
2097 : ! do ipw1=1,npwdiel
2098 : ! do ipw2=1,npwdiel
2099 : ! elementr=0.0_dp
2100 : ! elementi=0.0_dp
2101 : ! do ipw3=1,npwdiel
2102 : ! elementr=elementr+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1)&
2103 : ! & -sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1)
2104 : ! elementi=elementi+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1)&
2105 : ! & +sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1)
2106 : ! end do
2107 : ! if(elementr**2+elementi**2 > 1.0d-12)then
2108 : ! if( ipw1 /= ipw2 .or. &
2109 : ! & ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then
2110 : ! write(std_out,*)' dieltcel : sqrsus and invsqrsus are not (pseudo)',&
2111 : ! & 'inverse of each other'
2112 : ! write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2
2113 : ! write(std_out,*)' elementr,elementi=',elementr,elementi
2114 : ! stop
2115 : ! end if
2116 : ! end if
2117 : ! end do
2118 : ! end do
2119 : ! ENDDEBUG
2120 :
2121 : ! End loop over spins
2122 : end do
2123 :
2124 6 : ABI_FREE(eig_msusinvsqr)
2125 6 : ABI_FREE(eig_msussqr)
2126 6 : ABI_FREE(eig_sus)
2127 6 : ABI_FREE(sush)
2128 6 : ABI_FREE(susvec)
2129 :
2130 : !-Compute the Hxc kernel
2131 :
2132 30 : ABI_MALLOC(khxc,(2,npwdiel,nspden,npwdiel,nspden))
2133 30 : ABI_MALLOC(symdielmat,(2,npwdiel,nspden,npwdiel,nspden))
2134 :
2135 27858 : khxc(:,:,:,:,:)=0.0_dp
2136 :
2137 : !Compute Hartree kernel
2138 240 : do ipw1=1,npwdiel
2139 234 : kg_red1=dble(kg_diel(1,ipw1))
2140 234 : kg_red2=dble(kg_diel(2,ipw1))
2141 234 : kg_red3=dble(kg_diel(3,ipw1))
2142 : gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 &
2143 : & +2.0_dp*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + &
2144 234 : & gmet(2,3)*kg_red2*kg_red3) )
2145 : ! Distinguish G=0 from other elements
2146 240 : if(gsquar>1.0d-12)then
2147 228 : khxc(1,ipw1,1,ipw1,1)= 4.0_dp*pi/gsquar
2148 : else
2149 : ! G=0
2150 : ipw0=ipw1
2151 : end if
2152 : end do
2153 :
2154 : !Eventually add the xc part
2155 6 : if(option>=2)then
2156 :
2157 12 : ABI_MALLOC(wkxc,(nfft))
2158 12 : ABI_MALLOC(kxcg,(2,nfft))
2159 25924 : wkxc(:)=kxc(:,1)
2160 : ! DEBUG
2161 : ! Used to moderate divergenc effect near rho=0 of the Kxc (see above).
2162 : ! wkxc(:)=merge(kxc(:,1), kxc_min, kxc(:,1) > kxc_min)
2163 : ! ENDDEBUG
2164 4 : call initmpi_seq(mpi_enreg_seq)
2165 4 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft(2),ngfft(3),'all')
2166 4 : call fourdp(1,kxcg,wkxc,-1,mpi_enreg_seq,nfft,1,ngfft,0) ! trsfrm R to G
2167 4 : call destroy_mpi_enreg(mpi_enreg_seq)
2168 :
2169 : ! Compute difference in G vectors
2170 4 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
2171 160 : do ipw2=1,npwdiel
2172 160 : if(ipw2/=ipw0)then
2173 :
2174 152 : j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2)
2175 : ! Fills diagonal
2176 152 : khxc(1,ipw2,1,ipw2,1)=khxc(1,ipw2,1,ipw2,1)+kxcg(1,1)
2177 152 : khxc(2,ipw2,1,ipw2,1)=khxc(2,ipw2,1,ipw2,1)+kxcg(2,1)
2178 :
2179 152 : if(ipw2/=npwdiel)then
2180 : ! Fills off-diagonal part of the matrix, except G=0
2181 2960 : do ipw1=ipw2+1,npwdiel
2182 2960 : if(ipw1/=ipw0)then
2183 2812 : i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1)
2184 : ! Use of two mod calls handles both i1-j1>=ndiel1 AND i1-j1<0
2185 2812 : k1=mod(n1+mod(i1-j1,n1),n1)
2186 2812 : k2=mod(n2+mod(i2-j2,n2),n2)
2187 2812 : k3=mod(n3+mod(i3-j3,n3),n3)
2188 2812 : ifft=k1+1+n1*(k2+n2*k3)
2189 : ! The signs of imaginary contributions have been checked
2190 2812 : khxc(1,ipw1,1,ipw2,1)=kxcg(1,ifft)
2191 2812 : khxc(2,ipw1,1,ipw2,1)=kxcg(2,ifft)
2192 2812 : khxc(1,ipw2,1,ipw1,1)=kxcg(1,ifft)
2193 2812 : khxc(2,ipw2,1,ipw1,1)=-kxcg(2,ifft)
2194 : end if
2195 : end do
2196 : end if
2197 :
2198 : end if
2199 : end do
2200 :
2201 4 : ABI_FREE(wkxc)
2202 4 : ABI_FREE(kxcg)
2203 :
2204 : ! Endif option 2
2205 : end if
2206 :
2207 : !Now, get the symmetric dielectric matrix
2208 : !Premultiplication by square root of minus susceptibility matrix
2209 240 : do ipw2=1,npwdiel
2210 9366 : do ipw1=1,npwdiel
2211 : ar=0._dp ; ai=0._dp
2212 365040 : do ii=1,npwdiel
2213 : ar=ar+sqrsus(1,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) &
2214 355914 : & -sqrsus(2,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1)
2215 : ai=ai+sqrsus(2,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) &
2216 365040 : & +sqrsus(1,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1)
2217 : end do
2218 9126 : work(1,ipw1,1,ipw2,1)=ar
2219 9360 : work(2,ipw1,1,ipw2,1)=ai
2220 : end do
2221 : end do
2222 : !Postmultiplication by square root of minus susceptibility matrix
2223 240 : do ipw2=1,npwdiel
2224 : ! do ipw1=ipw2,npwdiel
2225 9360 : do ipw1=1,npwdiel
2226 : ar=0._dp ; ai=0._dp
2227 365040 : do ii=1,npwdiel
2228 : ar=ar+work(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2229 355914 : & -work(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2230 : ai=ai+work(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2231 365040 : & +work(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2232 : end do
2233 9126 : symdielmat(1,ipw1,1,ipw2,1)=ar
2234 9360 : symdielmat(2,ipw1,1,ipw2,1)=ai
2235 : ! if(ipw1/=ipw2)then
2236 : ! symdielmat(1,ipw2,1,ipw1,1)=ar
2237 : ! symdielmat(2,ipw2,1,ipw1,1)=-ai
2238 : ! end if
2239 : end do
2240 : ! Add the unity matrix
2241 240 : symdielmat(1,ipw2,1,ipw2,1)=1._dp+symdielmat(1,ipw2,1,ipw2,1)
2242 : end do
2243 :
2244 6 : ABI_FREE(khxc)
2245 :
2246 12 : ABI_MALLOC(symh,(npwdiel*(npwdiel+1)))
2247 30 : ABI_MALLOC(symvec,(2,npwdiel,nspden,npwdiel,nspden))
2248 12 : ABI_MALLOC(eig_sym,(npwdiel))
2249 :
2250 : !Store the symmetrized dielectric matrix in proper mode before calling zhpev
2251 240 : index=1
2252 240 : do ii=1,npwdiel
2253 4920 : do jj=1,ii
2254 4680 : symh(index )=symdielmat(1,jj,1,ii,1)
2255 4680 : symh(index+1)=symdielmat(2,jj,1,ii,1)
2256 4914 : index=index+2
2257 : end do
2258 : end do
2259 :
2260 6 : ier=0
2261 : call ZHPEV ('V','U',npwdiel,symh,eig_sym,symvec,npwdiel,zhpev1,&
2262 6 : & zhpev2,ier)
2263 :
2264 6 : if(prtvol>=10)then
2265 6 : write(message, '(a,a,a,5es12.4)' )ch10,&
2266 6 : & ' Five largest eigenvalues of the symmetrized dielectric matrix:',&
2267 12 : & ch10,eig_sym(npwdiel:npwdiel-4:-1)
2268 6 : call wrtout(ab_out,message,'COLL')
2269 : end if
2270 :
2271 6 : write(message,'(2a)')ch10,' dieltcel : 15 largest eigenvalues of the symmetrized dielectric matrix'
2272 6 : call wrtout(std_out,message,'COLL')
2273 6 : write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(npwdiel:npwdiel-4:-1)
2274 6 : call wrtout(std_out,message,'COLL')
2275 6 : write(message, '(a,5es12.5)' )' 6-10 :',eig_sym(npwdiel-5:npwdiel-9:-1)
2276 6 : call wrtout(std_out,message,'COLL')
2277 6 : write(message, '(a,5es12.5)' )' 11-15:',eig_sym(npwdiel-10:npwdiel-14:-1)
2278 6 : call wrtout(std_out,message,'COLL')
2279 6 : write(message, '(2a)' )ch10,' dieltcel : 5 smallest eigenvalues of the symmetrized dielectric matrix'
2280 6 : call wrtout(std_out,message,'COLL')
2281 6 : write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(1:5)
2282 6 : call wrtout(std_out,message,'COLL')
2283 :
2284 : !Invert the hermitian dielectric matrix,
2285 27858 : work(:,:,:,:,:)=0.0_dp
2286 240 : do ieig=1,npwdiel
2287 234 : eiginv=1.0_dp/eig_sym(ieig)
2288 9366 : do ipw2=1,npwdiel
2289 : ! do ipw1=ipw2,npwdiel
2290 365274 : do ipw1=1,npwdiel
2291 : work(1,ipw1,1,ipw2,1)=work(1,ipw1,1,ipw2,1)+&
2292 : & (symvec(1,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)+ &
2293 355914 : & symvec(2,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv
2294 : work(2,ipw1,1,ipw2,1)=work(2,ipw1,1,ipw2,1)+&
2295 : & (symvec(2,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)- &
2296 365040 : & symvec(1,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv
2297 : end do
2298 : end do
2299 : end do
2300 : !if(npwdiel>1)then
2301 : !do ipw2=2,npwdiel
2302 : !do ipw1=1,ipw2-1
2303 : !work(1,ipw1,1,ipw2,1)= work(1,ipw2,1,ipw1,1)
2304 : !work(2,ipw1,1,ipw2,1)=-work(2,ipw2,1,ipw1,1)
2305 : !end do
2306 : !end do
2307 : !end if
2308 :
2309 6 : ABI_FREE(eig_sym)
2310 6 : ABI_FREE(symh)
2311 6 : ABI_FREE(symvec)
2312 :
2313 : !DEBUG
2314 : !Checks whether the inverse of the symmetric dielectric matrix
2315 : !has been correctly generated
2316 : !do ipw1=1,npwdiel
2317 : !do ipw2=1,npwdiel
2318 : !elementr=0.0_dp
2319 : !elementi=0.0_dp
2320 : !do ipw3=1,npwdiel
2321 : !elementr=elementr+work(1,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1)&
2322 : !& -work(2,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1)
2323 : !elementi=elementi+work(1,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1)&
2324 : !& +work(2,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1)
2325 : !end do
2326 : !if(elementr**2+elementi**2 > 1.0d-12)then
2327 : !if( ipw1 /= ipw2 .or. &
2328 : !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then
2329 : !write(std_out,*)' dieltcel : the inversion procedure is not correct '
2330 : !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2
2331 : !write(std_out,*)' elementr,elementi=',elementr,elementi
2332 : !stop
2333 : !end if
2334 : !end if
2335 : !end do
2336 : !end do
2337 : !write(std_out,*)'dieltcel : matrix has been inverted successfully '
2338 : !ENDDEBUG
2339 :
2340 6 : ABI_FREE(symdielmat)
2341 :
2342 : !Then get the inverse of the asymmetric
2343 : !dielectric matrix, as required for the preconditioning.
2344 : !Premultiplication by square root of minus susceptibility matrix
2345 240 : do ipw2=1,npwdiel
2346 9366 : do ipw1=1,npwdiel
2347 : ar=0._dp ; ai=0._dp
2348 365040 : do ii=1,npwdiel
2349 : ar=ar+invsqrsus(1,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) &
2350 355914 : & -invsqrsus(2,ipw1,1,ii,1)*work(2,ii,1,ipw2,1)
2351 : ai=ai+invsqrsus(2,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) &
2352 365040 : & +invsqrsus(1,ipw1,1,ii,1)*work(2,ii,1,ipw2,1)
2353 : end do
2354 9126 : work2(1,ipw1,1,ipw2,1)=ar
2355 9360 : work2(2,ipw1,1,ipw2,1)=ai
2356 : end do
2357 : end do
2358 : !Postmultiplication by square root of minus susceptibility matrix
2359 240 : do ipw2=1,npwdiel
2360 9366 : do ipw1=1,npwdiel
2361 : ar=0._dp ; ai=0._dp
2362 365040 : do ii=1,npwdiel
2363 : ar=ar+work2(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2364 355914 : & -work2(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2365 : ai=ai+work2(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) &
2366 365040 : & +work2(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1)
2367 : end do
2368 9126 : dielinv(1,ipw1,1,ipw2,1)=ar
2369 9360 : dielinv(2,ipw1,1,ipw2,1)=ai
2370 : end do
2371 : end do
2372 :
2373 6 : ABI_FREE(invsqrsus)
2374 6 : ABI_FREE(sqrsus)
2375 6 : ABI_FREE(work)
2376 6 : ABI_FREE(work2)
2377 6 : ABI_FREE(zhpev1)
2378 6 : ABI_FREE(zhpev2)
2379 :
2380 6 : call timab(96,2,tsec)
2381 :
2382 6 : end subroutine dieltcel
2383 : !!***
2384 :
2385 : !!****f* ABINIT/chi0diel
2386 : !! NAME
2387 : !! chi0diel
2388 : !!
2389 : !! FUNCTION
2390 : !! Computes the preconditioned residual vrespc = P^-1 vresid where P is an approximation of the
2391 : !! dielectric matrix (if 'optres'=0) or its adjoint (if 'optres'=1) based of a model of the
2392 : !! non-interacting susceptibility chi0. The approximation is described by the object 'precon'
2393 : !! (see the abinit documentation of the input variable 'iprcel').
2394 : !!
2395 : !! INPUTS
2396 : !! precon = precon_object that contain the model chi0 operator.
2397 : !! dtset = All input variables for this dataset.
2398 : !! mpi_enreg = Information about MPI parallelization.
2399 : !! optreal = 1: vresid is given in the REAL space.
2400 : !! 2: vresid is given in the RECIPROCAL space.
2401 : !! optres = 0: the array vresid contains a potential residual.
2402 : !! 1: the array vresid contains a density residual.
2403 : !! vresid (:, nspden) = residual density/potential in REAL space (if optreal==1)
2404 : !! residual density/potential in RECIPROCAL space (if optreal==2)
2405 : !!
2406 : !! OUTPUT
2407 : !! vrespc (:, nspden) = preconditioned residual of the density/potential in REAL space
2408 : !! in REAL space (optreal==1) or RECIPROCAL space (optreal==2)
2409 : !!
2410 : !! SOURCE
2411 :
2412 18 : subroutine chi0diel(precon, dtset, mpi_enreg, optreal, optres, vresid, vrespc)
2413 :
2414 : !Arguments ------------------------------------
2415 : type(precon_object) :: precon
2416 : !scalars
2417 : integer,intent(in) :: optreal, optres
2418 : type(MPI_type),intent(in) :: mpi_enreg
2419 : type(dataset_type),intent(in) :: dtset
2420 : !arrays
2421 : real(dp),intent(in) :: vresid(optreal*precon%nfftprc, dtset%nspden)
2422 : real(dp),intent(out) :: vrespc(optreal*precon%nfftprc, dtset%nspden)
2423 :
2424 : ! *************************************************************************
2425 :
2426 18 : call precon%apply_precon(dtset, mpi_enreg, optreal, optres, vresid, vrespc)
2427 :
2428 : !Simple mixing
2429 121554 : vrespc = precon%diemix * vrespc
2430 :
2431 18 : end subroutine chi0diel
2432 : !!***
2433 :
2434 : !!****f* ABINIT/prcrskerker1
2435 : !! NAME
2436 : !! prcrskerker1
2437 : !!
2438 : !! FUNCTION
2439 : !! preconditionning by a real-space conjugate gradient on residual
2440 : !! using a model dielectric function in real space
2441 : !!
2442 : !! INPUTS
2443 : !! nfft=number of fft grid points
2444 : !! nspden=number of spin-density components
2445 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2446 : !! dielar(7)=input parameters for dielectric matrix:
2447 : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
2448 : !! gprimd(3,3)=dimensional primitive translations in fourier space (bohr**-1)
2449 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2450 : !! vresid(nfft,nspden)=residual potential
2451 : !! base(nfft) = real space function used as a basis to guess a fine dielectric function
2452 : !! see the calling routine to know the content
2453 : !!
2454 : !! OUTPUT
2455 : !! vrespc(nfft,nspden)=preconditioned residual of the potential
2456 : !!
2457 : !! WARNINGS
2458 : !! This is experimental code : input, ouptput, results and any other feature may vary greatly.
2459 : !!
2460 : !! NOTES
2461 : !! needs severe cleaning and this is abuse of modules as common blocks...
2462 : !!
2463 : !! SOURCE
2464 :
2465 7 : subroutine prcrskerker1(dtset,mpi_enreg,nfft,nspden,ngfft,dielar,etotal,gprimd,vresid,vrespc,base)
2466 :
2467 : !Arguments ------------------------------------
2468 : !scalars
2469 : integer,intent(in) :: nfft,nspden
2470 : real(dp) :: etotal
2471 : type(MPI_type),intent(in) :: mpi_enreg
2472 : type(dataset_type),intent(in) :: dtset
2473 : !arrays
2474 : integer,intent(in) :: ngfft(18)
2475 : real(dp),intent(in) :: base(nfft),dielar(7),gprimd(3,3)
2476 : real(dp),intent(in) :: vresid(nfft,nspden)
2477 : real(dp),intent(out) :: vrespc(nfft,nspden)
2478 :
2479 : !Local variables-------------------------------
2480 : !scalars
2481 : integer :: ifft,ispden,n1,n2,n3
2482 : real(dp) :: base_delta,base_max,base_min,dielng,diemac,diemix
2483 : real(dp) :: diemixmag
2484 : real(dp) :: rdummy1,rdummy2
2485 : logical :: new_prc_func
2486 : !arrays
2487 14 : real(dp) :: deltaW(nfft,nspden)
2488 14 : real(dp) :: g2cart(nfft)
2489 14 : real(dp) :: mat(nfft,nspden)
2490 : ! *************************************************************************
2491 :
2492 : !DEBUG
2493 : !write(std_out,*)' prckerker1 : enter '
2494 : !ENDDEBUG
2495 : !if(cycle==0) then
2496 7 : call prc_mem_init(nfft)
2497 :
2498 7 : if(cycle==0) then
2499 1 : new_prc_func=.TRUE.
2500 1 : energy_min=etotal
2501 6 : else if(etotal < energy_min) then
2502 5 : new_prc_func=.TRUE.
2503 5 : energy_min=etotal
2504 : else
2505 : new_prc_func=.FALSE.
2506 : end if
2507 :
2508 :
2509 7 : dielng=dielar(2)
2510 7 : diemac=dielar(3)
2511 7 : diemix=dielar(4)
2512 7 : diemixmag=dielar(7)
2513 : !******************************************************************
2514 : !compute the diemac(r) **
2515 : !******************************************************************
2516 : !this task will be devoted to a general function later
2517 7 : n1=ngfft(1)
2518 7 : n2=ngfft(2)
2519 7 : n3=ngfft(3)
2520 : !base_cp=base
2521 7 : if(new_prc_func) then
2522 6 : base_min=base(1)
2523 6 : base_max=base(1)
2524 60006 : do ifft=1,nfft
2525 60000 : base_min = min(base_min,base(ifft))
2526 60006 : base_max = max(base_max,base(ifft))
2527 : end do
2528 6 : base_delta = base_max - base_min
2529 : ! if(cycle.lt.2) then
2530 60006 : rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one)
2531 : ! else
2532 : ! rdiemac(:) = rdiemac(:)*0.5_dp+0.5_dp*(((base(:)-base_min) / (base_delta) ) *(diemac-one) + one)
2533 : ! end if
2534 : ! if(cycle==0) rdiemac(:) = (((base(:)-base_min) / (base_delta) ) *(diemac-one) + one)
2535 : ! rdiemac(:) = exp(((base(:)-base_min) / (base_delta) *log(diemac)))
2536 : end if
2537 7 : cycle=cycle+1
2538 : !if(cycle==5) cycle=0
2539 : !end if
2540 : !******************************************************************
2541 : !compute deltaW **
2542 : !******************************************************************
2543 70014 : vrespc=vresid !starting point
2544 : ! put the laplacian of the residuals into deltaW
2545 7 : call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,rdfuncr=vrespc,laplacerdfuncr=deltaW,g2cart_out=g2cart)
2546 :
2547 : !call laplacian(vrespc,buffer,ngfft,gprimd) ! put the laplacian of the residuals into deltaW
2548 : !do ifft=1,nfft
2549 : !if (buffer(ifft,1)/=deltaW(ifft,1)) then
2550 : !stop
2551 : !end if
2552 : !end do
2553 70007 : deltaW(:,1)= diemix*(((one/rdiemac(:))*vresid(:,1))-(((dielng)**2)*deltaW(:,1)))
2554 7 : if (nspden>1.and.(diemixmag>=zero)) then
2555 0 : do ispden=2,nspden
2556 0 : deltaW(:,ispden)= abs(diemixmag)*(((one/rdiemac(:))*vresid(:,ispden))-(((dielng)**2)*deltaW(:,ispden)))
2557 : end do
2558 : end if
2559 : !call random_number(deltaW)
2560 : !call random_number(vrespc)
2561 : !******************************************************************
2562 : !Finding the preconditionned residuals which minimizes **
2563 : !half*(vrespc*(1-dielng2/4pi2 nabla2) vrespc) - vrespc * deltaW **
2564 : !***********************************************************************
2565 70007 : vrespc(:,1)=diemix*vrespc(:,1)
2566 7 : if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*vrespc(:,2:nspden)
2567 : !buffer=vrespc
2568 :
2569 :
2570 : !==============================================================================
2571 : !==============================================================================
2572 : !! Original loop
2573 : !==============================================================================
2574 : !==============================================================================
2575 :
2576 7 : call frskerker1__init(dtset,mpi_enreg,nfft,ngfft,nspden,dielng,deltaW,gprimd,mat,g2cart)
2577 :
2578 : !call cgpr(pf_rscgres,dpf_rscgres,newvres,real(1e-40,dp),700,vrespc,rdummy1,rdummy2)
2579 : !rdummy1 = pf_rscgres(nfft,nspden,vrespc)
2580 : call cgpr(nfft,nspden,frskerker1__pf,frskerker1__dpf,frskerker1__newvres,&
2581 7 : & real(1e-10,dp),700,vrespc,rdummy1,rdummy2)
2582 7 : call frskerker1__end()
2583 :
2584 : !==============================================================================
2585 : !==============================================================================
2586 : !! Original loop end
2587 : !==============================================================================
2588 : !==============================================================================
2589 :
2590 :
2591 : !cplex=1
2592 : !qphon(:)=zero
2593 : !call moddiel(cplex,dielar,nfft,ngfft,nspden,1,0,qphon,rprimd,vresid,buffer)
2594 : !c1=0
2595 : !do ifft=1,nfft,1
2596 : !if((abs(buffer(ifft,1)-vrespc(ifft,1))/(abs(buffer(ifft,1)+vrespc(ifft,1))*half)) > 5e-3) then
2597 : !c1=c1+1
2598 : !end if
2599 : !end do
2600 : !call laplacian(vrespc,buffer,ngfft,gprimd)
2601 : !buffer=vrespc(:,:)-buffer(:,:)*dielng**2
2602 : !c2=0
2603 : !do ifft=1,nfft,1
2604 : !if((abs(buffer(ifft,1)-deltaW(ifft,1))/(abs(buffer(ifft,1)+deltaW(ifft,1))*half)) > 5e-3) then
2605 : !c2=c2+1
2606 : !end if
2607 : !end do
2608 : !!! !stop
2609 : !call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,&
2610 : !& g2cart_out=g2cart)
2611 :
2612 : !vrespc=vresid
2613 : !do ispden=1,nspden
2614 : !call fourdp(1, gvrespc(:,:,ispden), vrespc(:,ispden),-1,mpi_enreg,nfft,ngfft,0)
2615 : !end do
2616 : !filtering
2617 : !do ispden=1,nspden
2618 : !do ifft=1,nfft
2619 : !! gvrespc(:,ifft,ispden)=(one-exp(-g2cart(ifft)*15.0_dp))*gvrespc(:,ifft,ispden)
2620 : !! gvrespc(:,ifft,ispden)=(exp(-g2cart(ifft)*10.0_dp))*gvrespc(:,ifft,ispden)
2621 : !! gvrespc(:,ifft,ispden)=(one-one/(exp(-0.002_dp/g2cart(ifft)**2)+one))*gvrespc(:,ifft,ispden)
2622 : !gvrespc(:,ifft,ispden)=(two-2_dp/(exp(-0.008_dp/(g2cart(ifft)+0.0012_dp))+one))*gvrespc(:,ifft,ispden)
2623 : !gvrespc(:,ifft,ispden)=min(one,(sqrt(g2cart(ifft)/0.006_dp))**(one))*gvrespc(:,ifft,ispden)
2624 : !end do
2625 : !end do
2626 : !change resulting potential to real space
2627 : !do ispden=1,nspden
2628 : !call fourdp(1,gvrespc(:,:,ispden),vrespc(:,ispden),1,mpi_enreg,nfft,ngfft,0)
2629 : !end do
2630 : !vrespc=vrespc*diemix
2631 : !maxg2=g2cart(1)
2632 : !ming2=g2cart(5)
2633 : !do ifft=1,nfft
2634 : !maxg2=max(g2cart(ifft),maxg2)
2635 : !if(g2cart(ifft) .gt. zero) ming2=min(g2cart(ifft),ming2)
2636 : !end do
2637 : !stop
2638 :
2639 : !DEBUG
2640 : !write(std_out,*)' prckerker1 : exit '
2641 : !ENDDEBUG
2642 :
2643 7 : end subroutine prcrskerker1
2644 : !!***
2645 :
2646 : !!****f* ABINIT/prcrskerker2
2647 : !! NAME
2648 : !! prcrskerker2
2649 : !!
2650 : !! FUNCTION
2651 : !! preconditionning by a real-space conjugate gradient on residual
2652 : !! using a model dielectric function in real space
2653 : !! differing from prcrskerker1 by the
2654 : !! use of a linear response approach
2655 : !!
2656 : !! INPUTS
2657 : !! nfft=number of fft grid points
2658 : !! nspden=number of spin-density components
2659 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2660 : !! dielar(7)=input parameters for dielectric matrix:
2661 : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
2662 : !! gprimd(3,3)=dimensional primitive translations in fourier space (bohr**-1)
2663 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2664 : !! vresid(nfft,nspden)=residual potential
2665 : !!
2666 : !! OUTPUT
2667 : !! vrespc(nfft,nspden)=preconditioned residual of the potential
2668 : !!
2669 : !! WARNINGS
2670 : !! This is experimental code : input, ouptput, results and any other feature may vary greatly.
2671 : !!
2672 : !! NOTES
2673 : !!
2674 : !! SOURCE
2675 :
2676 4 : subroutine prcrskerker2(dtset,nfft,nspden,ngfft,dielar,gprimd,rprimd,vresid,vrespc,natom,xred,mpi_enreg,ucvol)
2677 :
2678 : !Arguments ------------------------------------
2679 : !scalars
2680 : integer,intent(in) :: natom,nfft,nspden
2681 : real(dp),intent(in) :: ucvol
2682 : type(MPI_type),intent(in) :: mpi_enreg
2683 : type(dataset_type),intent(in) :: dtset
2684 : !arrays
2685 : integer,intent(in) :: ngfft(18)
2686 : real(dp),intent(in) :: dielar(7),gprimd(3,3),rprimd(3,3),vresid(nfft,nspden)
2687 : real(dp),intent(in) :: xred(3,natom)
2688 : real(dp),intent(out) :: vrespc(nfft,nspden)
2689 :
2690 : !Local variables-------------------------------
2691 : !logical,save ::ok=.FALSE.
2692 : !scalars
2693 : integer :: cplex,i1,i2,i3,iatom,iatom27,ifft,ispden,n1,n2,n3,natom27,nfftotf
2694 : integer :: option
2695 : real(dp),save :: lastp1=one,lastp2=one
2696 : real(dp) :: C1,C2,DE,core,dielng,diemac,diemix,diemixmag,doti,dr,l1,l2,l3,l4,r
2697 : real(dp) :: rdummy1,rdummy2,rmin,xr,y,yr,zr
2698 : !arrays
2699 4 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
2700 4 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
2701 8 : real(dp) :: V1(nfft,nspden),V2(nfft,nspden),buffer(nfft,nspden)
2702 8 : real(dp) :: deltaW(nfft,nspden)
2703 8 : real(dp) :: mat(nfft,nspden)
2704 8 : real(dp) :: rdielng(nfft),rdiemac(nfft),xcart(3,natom)
2705 8 : real(dp) :: xcart27(3,natom*27)
2706 : ! *************************************************************************
2707 :
2708 4 : dielng=dielar(2)
2709 4 : diemac=dielar(3)
2710 4 : diemix=dielar(4)
2711 4 : diemixmag=dielar(7)
2712 : !******************************************************************
2713 : !compute the diemac(r) **
2714 : !******************************************************************
2715 : !this task will be devoted to a general function later
2716 4 : n1=ngfft(1)
2717 4 : n2=ngfft(2)
2718 4 : n3=ngfft(3)
2719 4 : nfftotf=n1*n2*n3
2720 : !if(.not.ok) then
2721 36 : xcart(1,:)=xred(1,:)*rprimd(1,1)+xred(2,:)*rprimd(1,2)+xred(3,:)*rprimd(1,3)
2722 36 : xcart(2,:)=xred(1,:)*rprimd(2,1)+xred(2,:)*rprimd(2,2)+xred(3,:)*rprimd(2,3)
2723 36 : xcart(3,:)=xred(1,:)*rprimd(3,1)+xred(2,:)*rprimd(3,2)+xred(3,:)*rprimd(3,3)
2724 :
2725 : iatom27=0
2726 16 : do i1=-1,1
2727 52 : do i2=-1,1
2728 156 : do i3=-1,1
2729 1008 : do iatom=1,natom
2730 864 : iatom27=iatom27+1
2731 3564 : xcart27(:,iatom27)=xcart(:,iatom)+rprimd(:,1)*i1+rprimd(:,2)*i2+rprimd(:,3)*i3
2732 : end do
2733 : end do
2734 : end do
2735 : end do
2736 :
2737 : !stop
2738 : natom27=27*natom
2739 :
2740 : l1=0.34580850339844665
2741 : !l2=0.5123510203906797 !0.41242551019533985
2742 : !l3=0.8001489796093203 !0.90007448980466009
2743 :
2744 : l2=0.41242551019533985
2745 : l3=0.90007448980466009
2746 : l4=0.9666914966015534
2747 :
2748 :
2749 : l1=0.31387233559896449
2750 : l2=0.35828367346355994
2751 : l3=0.9333829932031068
2752 : l4=0.9777943310677023
2753 :
2754 : l1=3.5
2755 : l2=11.5
2756 : l3=2.5
2757 : l4=6.5
2758 : !l1=30. !cellules pleines
2759 :
2760 40004 : rdielng=zero
2761 4 : core=1. !(value of Er at the core of atoms)
2762 4 : dr=2.65 ! radius of atoms=2.65165
2763 4 : y=1. ! value of Er in the empty region
2764 :
2765 : !Get the distrib associated with this fft_grid
2766 4 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2767 :
2768 404 : do i3=1,n3
2769 400 : ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft)
2770 4404 : do i2=1,n2
2771 4400 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
2772 44000 : do i1=1,n1
2773 40000 : ifft=ifft+1
2774 : ! !!!!!!!!!!!!!!!!!!!!!!!!!
2775 : ! ! calculation of the simplest part void/metal
2776 : ! !! x=real(real(i3,dp)/real(n3,dp),dp)
2777 : ! !! !x=i3/n3
2778 : ! !! if(x < l1) then
2779 : ! !! rdiemac(ifft)=diemac
2780 : ! !! rdielng(ifft)=dielng
2781 : ! !! else if(x < l2) then
2782 : ! !! xp=(l2-x)/(l2-l1)
2783 : ! !! rdiemac(ifft)=y+(diemac-y)&
2784 : ! !! & * (1.-(1.-xp)**4)**4
2785 : ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4
2786 : ! !! else if(x < l3) then
2787 : ! !! rdiemac(ifft)=y
2788 : ! !! rdielng(ifft)=zero
2789 : ! !! else if(x < l4) then
2790 : ! !! xp=(l3-x)/(l3-l4)
2791 : ! !! rdiemac(ifft)=y+(diemac-y)&
2792 : ! !! & * (1.-(1.-xp)**4)**4
2793 : ! !! rdielng(ifft)=dielng*(1.-(1.-xp)**4)**4
2794 : ! !! else
2795 : ! !! rdiemac(ifft)=diemac
2796 : ! !! rdielng(ifft)=dielng
2797 : ! !! end if
2798 : ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2799 : ! !!!! calculation of atomic core dielectric
2800 : ! !! rmin=1e16
2801 : ! !! xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)&
2802 : ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(1,3)
2803 : ! !! yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)&
2804 : ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(2,3)
2805 : ! !! zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)&
2806 : ! !! &+real((i3-1),dp)/real(n3,dp)*rprimd(3,3)
2807 : ! !! do iatom=1,natom27
2808 : ! !! r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2
2809 : ! !! if (r<rmin) then
2810 : ! !! rmin=r
2811 : ! !! end if
2812 : ! !! end do
2813 : ! !! if(rmin < dr**2) then
2814 : ! !! rdiemac(ifft)=min(rdiemac(ifft),core+(diemac-core)*(1.-(1.-sqrt(rmin)/dr)**2)**2)
2815 : ! !! rdielng(ifft)=dielng-dielng*(1.-(1.-sqrt(rmin)/dr)**4)**4
2816 : ! !! else
2817 : ! !! rdiemac(ifft)=min(rdiemac(ifft),diemac)
2818 : ! !! end if
2819 40000 : rmin=1e16
2820 : xr=real(real((i1-1),dp)/n1,dp)*rprimd(1,1)+real(real((i2-1),dp)/n2,dp)*rprimd(1,2)&
2821 40000 : & +real((i3-1),dp)/real(n3,dp)*rprimd(1,3)
2822 : yr=real(real((i1-1),dp)/n1,dp)*rprimd(2,1)+real(real((i2-1),dp)/n2,dp)*rprimd(2,2)&
2823 40000 : & +real((i3-1),dp)/real(n3,dp)*rprimd(2,3)
2824 : zr=real(real((i1-1),dp)/n1,dp)*rprimd(3,1)+real(real((i2-1),dp)/n2,dp)*rprimd(3,2)&
2825 40000 : & +real((i3-1),dp)/real(n3,dp)*rprimd(3,3)
2826 :
2827 40000 : rdiemac(ifft)=y
2828 40000 : rdielng(ifft)=zero
2829 8680000 : do iatom=1,natom27
2830 8640000 : r=(xr-xcart27(1,iatom))**2+(yr-xcart27(2,iatom))**2+(zr-xcart27(3,iatom))**2
2831 : if (r<rmin) then
2832 : rmin=r
2833 : end if
2834 8640000 : if(r < l1) then
2835 4840 : rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y)
2836 8635160 : else if(r < l2) then
2837 23680 : rdiemac(ifft)= rdiemac(ifft) + 0.7_dp * (diemac-y)*(one-((sqrt(r)-l1)/(l2-l1))**2)**2
2838 : else
2839 : rdiemac(ifft)=rdiemac(ifft)
2840 : end if
2841 8680000 : if(r < l3) then
2842 2792 : rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng)
2843 8637208 : else if(r < l4) then
2844 9408 : rdielng(ifft)= rdielng(ifft) + 0.5_dp * (dielng) *(one-((sqrt(r)-l3)/(l4-l3))**2)**2
2845 : end if
2846 : end do
2847 :
2848 40000 : rdielng(ifft)=min(rdielng(ifft),dielng)
2849 : ! rdielng(ifft)=dielng
2850 : rdiemac(ifft)=min(rdiemac(ifft),diemac)
2851 44000 : rdiemac(ifft)=diemac
2852 : end do
2853 : end if
2854 : end do
2855 : end do
2856 : !rdielng(:)=dielng
2857 :
2858 : !****************************************************************************************
2859 : !****************************************************************************************
2860 : !****************************************************************************************
2861 : !****************************************************************************************
2862 : !******************************************************************
2863 : !compute V1
2864 : !******************************************************************
2865 40008 : V1=vresid
2866 4 : call laplacian(gprimd,mpi_enreg,nfft,nspden,ngfft,rdfuncr=V1,laplacerdfuncr=deltaW)
2867 40004 : deltaW(:,1)= (((one/rdiemac(:))*V1(:,1))-(((rdielng(:))**2)*deltaW(:,1)))
2868 : !deltaW(:,1)= -diemix*(((rdielng(:))**2)*deltaW(:,ispden))
2869 4 : if (nspden>1) then
2870 0 : do ispden=2,nspden
2871 0 : deltaW(:,ispden)= (((one/rdiemac(:))*V1(:,ispden))-(((rdielng(:))**2)*deltaW(:,ispden)))
2872 : ! deltaW(:,ispden)= -abs(diemixmag)*(((rdielng(:))**2)*deltaW(:,ispden))
2873 : end do
2874 : end if
2875 4 : call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat)
2876 : call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,&
2877 4 : & frskerker2__newvres2,lastp1*real(1e-6 ,dp),700,V1,rdummy1,rdummy2)
2878 4 : lastp1=min(abs(rdummy1),1e-6_dp)
2879 4 : call frskerker2__end()
2880 :
2881 : !******************************************************************
2882 : !compute V2
2883 : !******************************************************************
2884 40008 : V2=vresid
2885 8 : do ispden=1,nspden
2886 40008 : deltaW(:,ispden)= (rdielng(:)**2)
2887 : end do
2888 4 : call frskerker2__init(dtset,mpi_enreg,nfft,ngfft,nspden,rdielng,deltaW,gprimd,mat)
2889 : call cgpr(nfft,nspden,frskerker2__pf,frskerker2__dpf,&
2890 4 : & frskerker2__newvres2,lastp2*real(1e-6,dp),700,V2,rdummy1,rdummy2)
2891 4 : lastp2=min(abs(rdummy1),1e-6_dp)
2892 4 : call frskerker2__end()
2893 :
2894 :
2895 : !******************************************************************
2896 : !compute C1, C2 & DE
2897 : !******************************************************************
2898 4 : cplex=1;
2899 4 : option=1;
2900 : call dotprod_vn(cplex,& !complex density/pot
2901 : &rdielng,& !the density
2902 : &DE,& !resulting dorproduct integrated over r ! here DE is used has a buffer
2903 : &doti,& !imaginary part of the integral
2904 : &size(rdielng,1),& !number of locally(cpu) attributed grid point
2905 : &nfftotf,& !real total number of grid point
2906 : &nspden,& !nspden
2907 : &option,& !1=compute only the real part 2=compute also the imaginary part
2908 : &rdielng,& !the potential
2909 : &ucvol,& !cell volume
2910 4 : &mpi_comm_sphgrid=mpi_enreg%comm_fft)
2911 8 : do ispden=1,nspden
2912 40008 : buffer(:,ispden)=rdielng(:)*V1(:,ispden)
2913 : end do
2914 : call dotprod_vn(cplex,& !complex density/pot
2915 : &rdielng,& !the density
2916 : &C1,& !resulting dorproduct integrated over r ! here DE is used has a buffer
2917 : &doti,& !imaginary part of the integral
2918 : &size(rdielng,1),& !number of locally(cpu) attributed grid point
2919 : &nfftotf,& !real total number of grid point
2920 : &nspden,& !nspden
2921 : &option,& !1=compute only the real part 2=compute also the imaginary part
2922 : &buffer,& !the potential
2923 : &ucvol,& !cell volume
2924 4 : &mpi_comm_sphgrid=mpi_enreg%comm_fft)
2925 8 : do ispden=1,nspden
2926 40008 : buffer(:,ispden)=rdielng(:)*V2(:,ispden)
2927 : end do
2928 : call dotprod_vn(cplex,& !complex density/pot
2929 : &rdielng,& !the density
2930 : &C2,& !resulting dorproduct integrated over r ! here DE is used has a buffer
2931 : &doti,& !imaginary part of the integral
2932 : &size(rdielng,1),& !number of locally(cpu) attributed grid point
2933 : &nfftotf,& !real total number of grid point
2934 : &nspden,& !nspden
2935 : &option,& !1=compute only the real part 2=compute also the imaginary part
2936 : &buffer,& !the potential
2937 : &ucvol,& !cell volume
2938 4 : &mpi_comm_sphgrid=mpi_enreg%comm_fft)
2939 4 : C1=C1/DE
2940 4 : C2=C2/DE
2941 4 : DE=C1/(one-C2)
2942 :
2943 : !******************************************************************
2944 : !compute the new preconditionned residuals
2945 : !******************************************************************
2946 40004 : vrespc(:,1)=diemix*(V1(:,1)+DE*V2(:,1))
2947 4 : if (nspden>1) vrespc(:,2:nspden)=abs(diemixmag)*(V1(:,2:nspden)+DE*V2(:,2:nspden))
2948 :
2949 4 : end subroutine prcrskerker2
2950 : !!***
2951 :
2952 : !!****f* ABINIT/cgpr
2953 : !! NAME
2954 : !! cgpr
2955 : !!
2956 : !! FUNCTION
2957 : !! perform Polak-Ribiere conjugate gradient on a function f
2958 : !! implementation based on the cg recipe of "numerical recipe"
2959 : !!
2960 : !! INPUTS
2961 : !! dp_dum_vdp: function to be minimized (return a dp from a vector of dp)
2962 : !! vdp_dum_vdp: derivative of f
2963 : !! dtol: precision precision required for the minimization
2964 : !! itmax: number of iterations allowed (each linmin will be done with at max 10 times
2965 : !! this number
2966 : !!
2967 : !! OUTPUT
2968 : !! fmin: value of f at the minimum
2969 : !! lastdelta: absolute value of the last delta between steps
2970 : !! SIDE EFFECTS
2971 : !! v: vector on which minimization is to be performed, starting point
2972 : !! and resulting min
2973 : !!
2974 : !! SOURCE
2975 :
2976 15 : subroutine cgpr(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,dtol,itmax,v,fmin,delta)
2977 :
2978 : !Arguments ------------------------------------
2979 : include "dummy_functions.inc"
2980 : !scalars
2981 : integer,intent(in) :: itmax,nv1,nv2
2982 : real(dp),intent(in) :: dtol
2983 : real(dp),intent(out) :: delta,fmin
2984 : !arrays
2985 : real(dp),intent(inout) :: v(nv1,nv2)
2986 :
2987 : !Local variables-------------------------------
2988 : !scalars
2989 : integer :: iiter
2990 : real(dp) :: fv,gam,gscal,gscal2,sto
2991 : !arrays
2992 15 : real(dp) :: grad0(nv1,nv2),grad1(nv1,nv2),grad2(nv1,nv2),grad3(nv1,nv2)
2993 : !no_abirules
2994 :
2995 : !************************************************************************
2996 15 : fv = dp_dum_v2dp(nv1,nv2,v(:,:))
2997 150030 : grad0(:,:) = -v2dp_dum_v2dp(nv1,nv2,v(:,:))
2998 150030 : grad1(:,:) = grad0(:,:)
2999 150030 : grad2(:,:) = grad0(:,:)
3000 448 : do iiter=1,itmax
3001 448 : call linmin(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,v,grad0,fmin)
3002 : ! return if the min is reached
3003 448 : sto=dtol*(abs(fmin)+abs(fv)+tol14)
3004 448 : delta=abs(fv-fmin)
3005 : delta=abs(delta)
3006 448 : if((delta.lt.sto).or.(iiter==itmax)) then
3007 : ! DEBUG
3008 : ! write(std_out,*) 'cgpr (01cg) : stop cond for cgpr:',sto,'delta:',delta,'fv:',fv,'fmin:',fmin
3009 : ! ENDDEBUG
3010 : return
3011 : end if
3012 : ! a new step
3013 433 : fv=fmin
3014 433 : grad0(:,:)=v2dp_dum_v2dp(nv1,nv2,v(:,:))
3015 433 : gscal=dotproduct(nv1,nv2,grad1(:,:),grad1(:,:))
3016 4330866 : grad3(:,:)=grad0(:,:)+grad1(:,:)
3017 433 : gscal2=dotproduct(nv1,nv2,grad3(:,:),grad0(:,:))
3018 433 : gam=gscal2/gscal
3019 4330866 : grad1(:,:)=-grad0(:,:)
3020 4330866 : grad2(:,:)=grad1(:,:)+gam*grad2(:,:)
3021 4331314 : grad0(:,:)=grad2(:,:)
3022 : ! DEBUG
3023 : ! write(std_out,*) 'cgpr (01cg) :================================================================================='
3024 : ! write(std_out,*) 'cgpr (01cg) : step',iiter,'delta:',delta ,'fv',fv,'fmin',fmin
3025 : ! write(std_out,*) 'cgpr (01cg) :================================================================================='
3026 : ! ENDDEBUG
3027 : end do
3028 :
3029 0 : end subroutine cgpr
3030 : !!***
3031 :
3032 : !!****f* ABINIT/linmin
3033 : !! NAME
3034 : !! linmin
3035 : !!
3036 : !! FUNCTION
3037 : !! minimizes a function along a gradient line:
3038 : !! first bracket the minimum then perform the minimization
3039 : !!
3040 : !! COPYRIGHT
3041 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, MT)
3042 : !! This file is distributed under the terms of the
3043 : !! GNU General Public License, see ~ABINIT/COPYING
3044 : !! or http://www.gnu.org/copyleft/gpl.txt .
3045 : !! For the initials of contributors, see ~ABINIT/Infos/contributors .
3046 : !!
3047 : !! INPUTS
3048 : !! dp_dum_vdp: function to be minimized (return a dp from a vector of dp)
3049 : !! vdp_dum_vdp: derivative of f
3050 : !!
3051 : !! OUTPUT
3052 : !! fmin: minimum value reached for dp_dum_vdp
3053 : !!
3054 : !! SIDE EFFECTS
3055 : !! grad: the gradient line along which the minimization is performed (not changed)
3056 : !! v: the starting and then ending point of the minimization
3057 : !!
3058 : !! SOURCE
3059 :
3060 448 : subroutine linmin(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,v,grad,fmin)
3061 :
3062 : !Arguments ------------------------------------
3063 : include "dummy_functions.inc"
3064 : !scalars
3065 : integer,intent(in) :: nv1,nv2
3066 : real(dp),intent(out) :: fmin
3067 : !arrays
3068 : real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2)
3069 :
3070 : !Local variables-------------------------------
3071 : !scalars
3072 : real(dp),parameter :: maglimit=10000.0_dp,tol=tol8*tol8*tol8
3073 : real(dp) :: a,b,fa,fb,fx,x,xmin
3074 : !no_abirules
3075 :
3076 : !************************************************************************
3077 448 : a=zero
3078 448 : x=ninth*real(1e-4,dp)
3079 448 : call bracketing (nv1,nv2,dp_dum_v2dp,v,grad,a,x,b,fa,fx,fb)
3080 : !DEBUG
3081 : !write(std_out,*) 'linmin (01cg) : linmin after bracketing'
3082 : !write(std_out,*) 'linmin (01cg) : point',a,x,b,'value',fa,fx,fb
3083 : !ENDDEBUG
3084 448 : fmin =brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,6,v,grad,a,x,b,tol,xmin)
3085 :
3086 448 : end subroutine linmin
3087 : !!***
3088 :
3089 : !!****f* ABINIT/bracketing
3090 : !! NAME
3091 : !! bracketing
3092 : !!
3093 : !! FUNCTION
3094 : !! bracket a minimum of a function f
3095 : !!
3096 : !! INPUTS
3097 : !! dp_dum_vdp: the function of which the mimimum should be bracketted
3098 : !!
3099 : !! OUTPUT
3100 : !! b= last member of the bracketing triplet a < x < b
3101 : !! fa,fx,fb= value of the function at dp_dum_vdp(v(:)+y*grad(:))
3102 : !!
3103 : !! SIDE EFFECTS
3104 : !! v: the initial vector for the function (return unchanged)
3105 : !! grad: the direction on which the bracketting is to be performed (return unchanged)
3106 : !! a,x: two members of the bracketing triplet (see b)
3107 : !!
3108 : !! SOURCE
3109 :
3110 448 : subroutine bracketing (nv1,nv2,dp_dum_v2dp,v,grad,a,x,b,fa,fx,fb)
3111 :
3112 : !Arguments ------------------------------------
3113 : include "dummy_functions.inc"
3114 : !scalars
3115 : integer,intent(in) :: nv1,nv2
3116 : real(dp),intent(inout) :: a,x
3117 : real(dp),intent(out) :: b,fa,fb,fx
3118 : !arrays
3119 : real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2)
3120 :
3121 : !Local variables-------------------------------
3122 : !scalars
3123 : real(dp),parameter :: maglimit=10000.0_dp
3124 : real(dp) :: c,fu,q,r,u,ulim
3125 : ! *************************************************************************
3126 :
3127 4480896 : fa=dp_dum_v2dp(nv1,nv2,v(:,:)+(a*grad(:,:)))
3128 4480896 : fx=dp_dum_v2dp(nv1,nv2,(x*grad(:,:))+v(:,:))
3129 448 : if(fx > fa) then
3130 26 : c=a
3131 26 : a=x
3132 26 : x=c
3133 26 : c=fa
3134 26 : fa=fx
3135 26 : fx=c
3136 : end if
3137 448 : b=x+gold*(x-a)
3138 4480896 : fb=dp_dum_v2dp(nv1,nv2,(b*grad(:,:))+v(:,:))
3139 574 : do
3140 1022 : if (fx <= fb) return
3141 589 : r=(x-a)*(fx-fb)
3142 589 : q=(x-b)*(fx-fa)
3143 589 : u=x-((x-b)*q-(x-a)*r)/(two*sign(max(abs(q-r),smallest_real),q-r))
3144 589 : ulim=x+maglimit*(b-x)
3145 589 : if((x-u)*(u-b) > zero) then
3146 150030 : fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
3147 15 : if(fu < fb) then
3148 15 : a=x
3149 15 : fa=fx
3150 15 : x=u
3151 15 : fx=fu
3152 15 : return
3153 0 : else if (fx < fu) then
3154 0 : b=u
3155 0 : fb=fu
3156 0 : return
3157 : end if
3158 0 : u=b+gold*(b-x)
3159 0 : fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
3160 574 : else if((b-u)*(u-ulim) > zero) then
3161 4880976 : fu=dp_dum_v2dp(nv1,nv2,u*grad(:,:)+v(:,:))
3162 488 : if(fu<fb) then
3163 478 : x=b
3164 478 : b=u
3165 478 : u=b+gold*(b-x)
3166 478 : fx=fb
3167 478 : fb=fu
3168 4780956 : fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
3169 : end if
3170 86 : else if((u-ulim)*(ulim-b) >= zero) then
3171 27 : u=ulim
3172 270054 : fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
3173 : else
3174 59 : u=b+gold*(b-x)
3175 590118 : fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
3176 : end if
3177 574 : a=x
3178 574 : x=b
3179 574 : b=u
3180 574 : fa=fx
3181 574 : fx=fb
3182 574 : fb=fu
3183 : end do
3184 :
3185 : end subroutine bracketing
3186 : !!***
3187 :
3188 : !!****f* ABINIT/brent
3189 : !! NAME
3190 : !! brent
3191 : !!
3192 : !! FUNCTION
3193 : !! minimizes a function along a line
3194 : !!
3195 : !! INPUTS
3196 : !! dp_dum_vdp: function to be minimized (return a dp from a vector of dp)
3197 : !! vdp_dum_vdp: derivative of the function (return a vector of dp from a vector of dp)
3198 : !! itmax: number of iterations allowed
3199 : !! tol: tolerance on error. It depend on the precision of the numbers
3200 : !! (usually chosen as sqrt(max precision available with your floating point reresentation))
3201 : !! ax,xx,bx: a bracketing triplet around the minimum to be find
3202 : !! OUTPUT
3203 : !! xmin: value such that dp_dum_vdp(v(:)+xmin*grad(:)) is minimum
3204 : !! brent: dp_dum_vdp(v(:)+xmin*grad(:))
3205 : !!
3206 : !! SIDE EFFECTS
3207 : !! grad(:): direction along which the minimization is performed
3208 : !! v(:): starting and ending point of the minimization
3209 : !!
3210 : !! SOURCE
3211 :
3212 448 : function brent(nv1,nv2,dp_dum_v2dp,v2dp_dum_v2dp,sub_dum_dp_v2dp_v2dp,itmax,v,grad,ax,xx,bx,tol,xmin)
3213 :
3214 : !Arguments ------------------------------------
3215 : include "dummy_functions.inc"
3216 : !scalars
3217 : integer,intent(in) :: itmax,nv1,nv2
3218 : real(dp) :: brent
3219 : real(dp),intent(in) :: ax,bx,tol,xx
3220 : real(dp),intent(out) :: xmin
3221 : !arrays
3222 : real(dp),intent(inout) :: grad(nv1,nv2),v(nv1,nv2)
3223 :
3224 : !Local variables-------------------------------
3225 : !scalars
3226 : integer :: iter
3227 : real(dp) :: a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,vv,w
3228 : real(dp) :: x,xm,zeps
3229 : logical :: ok1,ok2,ok3,ok4
3230 :
3231 : !************************************************************************
3232 448 : zeps=epsilon(ax*real(1e-2,dp))
3233 448 : a=min(ax,bx)
3234 448 : b=max(ax,bx)
3235 448 : vv=xx
3236 448 : w=xx
3237 448 : x=xx
3238 448 : e=zero
3239 4480896 : fx=dp_dum_v2dp(nv1,nv2,x*grad(:,:)+v(:,:))
3240 448 : fv=fx
3241 448 : fw=fx
3242 : !the function sub_dum_dp_v2dp_v2dp must do the equivalent of
3243 : !v(:,:)=v(:,:)+(grad(:,:)*x)
3244 : !but for instance renormilizing the density if brent is used on a density...
3245 : !vp(:,:) = v(:,:)
3246 : !sub_dum_dp_v2dp_v2dp(x,grad(:,:),vp(:,:)
3247 : !dx=dotproduct(v2dp_dum_v2dp(vp(:,:)),grad(:,:))
3248 4480896 : dx=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,v(:,:)+x*grad(:,:)),grad(:,:))
3249 448 : dv=dx
3250 448 : dw=dx
3251 2562 : do iter=1,itmax
3252 2350 : xm=half*(a+b)
3253 2350 : tol1=tol*abs(x)+zeps
3254 2350 : tol2=two*tol1
3255 2350 : if(abs(x-xm) <= (tol2-half*(b-a))) then
3256 : exit
3257 : end if
3258 2253 : if(abs(e) > tol1) then
3259 1744 : d1=two*(b-a)
3260 1744 : d2=d1
3261 1744 : if(dw /= dx) d1=(w-x)*dx/(dx-dw)
3262 1744 : if(dv /= dx) d2=(vv-x)*dx/(dx-dv)
3263 1744 : u1=x+d1
3264 1744 : u2=x+d2
3265 1744 : ok1=((a-u1)*(u1-b)>zero).and.(dx*d1<=zero)
3266 1744 : ok2=((a-u2)*(u2-b)>zero).and.(dx*d2<=zero)
3267 2225 : olde=e
3268 2225 : e=d
3269 1170 : if(ok1.or.ok2) then
3270 1055 : if(ok1.and.ok2) then
3271 549 : d=merge(d1,d2,abs(d1)<abs(d2))
3272 : else
3273 506 : d=merge(d1,d2,ok1)
3274 : end if
3275 1055 : if(abs(d)<=abs(half*olde)) then
3276 1013 : u=x+d
3277 1013 : if(((u-a)<tol2).or.((b-u)<tol2)) d=sign(tol1,xm-x)
3278 : else
3279 42 : e=merge(a,b,dx>=zero)-x
3280 42 : d=half*e
3281 : end if
3282 : else
3283 689 : e=merge(a,b,dx>=zero)-x
3284 689 : d=half*e
3285 : end if
3286 : else
3287 509 : e=merge(a,b,dx>=zero)-x
3288 509 : d=half*e
3289 : end if
3290 :
3291 2253 : if(abs(d) >=tol1)then
3292 1950 : u=x+d
3293 19503900 : fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
3294 : else
3295 303 : u=x+sign(tol1,d)
3296 3030606 : fu=dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:))
3297 303 : if(fu>fx) then
3298 : exit
3299 : end if
3300 : end if
3301 21144228 : du=dotproduct(nv1,nv2,v2dp_dum_v2dp(nv1,nv2,(u*grad(:,:))+v(:,:)),grad(:,:))
3302 2562 : if(fu<=fx)then
3303 809 : if(u>=x)then
3304 : a=x
3305 : else
3306 412 : b=x
3307 : end if
3308 809 : vv=w
3309 809 : fv=fw
3310 809 : dv=dw
3311 809 : w=x
3312 809 : fw=fx
3313 809 : dw=dx
3314 809 : x=u
3315 809 : dx=du
3316 809 : fx=fu
3317 : else
3318 1305 : if(u<x) then
3319 : a=u
3320 : else
3321 733 : b=u
3322 : end if
3323 1305 : ok3=(w==x).or.(fu.le.fw)
3324 1305 : ok4=(vv==w).or.(vv==x).or.(fu.lt.fv)
3325 1305 : if(ok3) then
3326 : vv=w
3327 : fv=fw
3328 : dv=dw
3329 : w=u
3330 : fw=fu
3331 : dw=du
3332 160 : else if( ok4 ) then
3333 86 : vv=u
3334 86 : fv=fu
3335 86 : dv=du
3336 : end if
3337 : end if
3338 : end do
3339 448 : xmin=x
3340 : !the function sub_dum_dp_v2dp_v2dp must do the equivalent of
3341 : !v(:,:)=v(:,:)+(grad(:,:)*x)
3342 : !but for instance renormilizing the density if brent is used on a density...
3343 448 : call sub_dum_dp_v2dp_v2dp(nv1,nv2,x,grad(:,:),v(:,:))
3344 448 : brent=fx
3345 :
3346 448 : end function brent
3347 : !!***
3348 :
3349 : end module m_prcref
3350 : !!***
|
