Branch data 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 : 13301 : subroutine prcref(atindx,dielar,dielinv,&
192 : 13301 : & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,&
193 : 13301 : & istep,kg_diel,kxc,&
194 : : & mgfft,moved_atm_inside,mpi_enreg,my_natom,&
195 : 13301 : & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
196 : 13301 : & optreal,optres,pawrhoij,pawtab,ph1d,precon,psps,rhog,rhoijrespc,rhor,rprimd,&
197 : 13301 : & 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 : 13301 : integer,allocatable :: indpw_prc(:)
244 : : real(dp) :: dummy6(6),dummy7(6),gprimd(3,3),qphon(3),rmet(3,3)
245 : 26602 : real(dp) :: vmean(dtset%nspden),vprtrb(2)
246 : 13301 : real(dp),allocatable :: dummy(:),dummy1(:),dummy2(:),dummy3(:),dummy4(:),dummy5(:),dummy8(:),dummy9(:)
247 : 13301 : real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:)
248 : 13301 : real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:)
249 : 13301 : real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:),rhor_new(:,:)
250 : 13301 : real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:)
251 : 13301 : real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:)
252 : 13301 : real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:)
253 : 13301 : logical,allocatable :: mask(:)
254 : : ! *************************************************************************
255 : :
256 : : !Compute different geometric tensor, as well as ucvol, from rprimd
257 : 13301 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
258 : :
259 : : !1) Eventually take care of the forces
260 : :
261 [ - + ]: 13301 : 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 [ - + ]: 13301 : 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 [ + + - + : 13301 : if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then
- - ]
336 : 13283 : cplex=optreal
337 : 13283 : qphon(:)=zero
338 : : ! Simple scalar multiplication, or model dielectric function
339 : 13283 : 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 [ + - + + : 13301 : if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then
+ - ]
481 : : mixfac=one;mixfacmag=one
482 : : else
483 : 13301 : mixfac=dielar(4);mixfacmag=abs(dielar(7))
484 : : end if
485 [ + + + + ]: 13301 : if (psps%usepaw==1.and.my_natom>0) then
486 [ + + ]: 11640 : if (pawrhoij(1)%cplex_rhoij==1) then
487 : : index=0
488 [ + + ]: 37830 : do iatom=1,my_natom
489 [ + + ]: 65863 : do iq=1,pawrhoij(iatom)%qphase
490 [ - + ]: 28033 : iq0=merge(0,pawrhoij(iatom)%lmn2_size,iq==1)
491 [ + + ]: 89316 : do ispden=1,pawrhoij(iatom)%nspden
492 [ + + ]: 33250 : mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
493 [ + + ]: 2170027 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
494 : 2108744 : index=index+1;klmn=iq0+pawrhoij(iatom)%kpawmix(kmix)
495 : 2141994 : 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 [ + + ]: 4446 : do iatom=1,my_natom
503 [ + + ]: 7049 : do iq=1,pawrhoij(iatom)%qphase
504 [ - + ]: 2603 : iq0=merge(0,2*pawrhoij(iatom)%lmn2_size,iq==1)
505 [ + + ]: 14289 : do ispden=1,pawrhoij(iatom)%nspden
506 [ + + ]: 9083 : mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
507 [ + + ]: 925950 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
508 : 914264 : index=index+2;klmn=iq0+2*pawrhoij(iatom)%kpawmix(kmix)-1
509 [ + + ]: 2751875 : 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 [ + + + + ]: 13301 : 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 [ + + ]: 13301 : 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 [ + - - - ]: 13301 : 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 [ - + ]: 13301 : 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 : 28371 : subroutine prcref_PMA(atindx,dielar,dielinv,&
851 : 28371 : & dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,&
852 : 28371 : & istep,kg_diel,kxc,&
853 : : & mgfft,moved_atm_inside,mpi_enreg,my_natom,&
854 : 28371 : & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
855 : 28371 : & optreal,optres,pawrhoij,ph1d,precon,psps,rhog, rhoijrespc,rhor,rprimd,&
856 : 28371 : & susmat,vhartr,vpsp,vresid,vrespc,vxc,xred,&
857 : 28371 : & 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 : 28371 : integer,allocatable :: indpw_prc(:)
904 : : real(dp) :: dummy6(6),gprimd(3,3),qphon(3),rmet(3,3)
905 : 56742 : real(dp) :: vmean(dtset%nspden),vprtrb(2)
906 : 28371 : 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 : 28371 : real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:)
909 : 28371 : real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:)
910 : 28371 : real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:)
911 : 28371 : real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:)
912 : 28371 : real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:)
913 : 28371 : real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:)
914 : 28371 : logical,allocatable :: mask(:)
915 : : ! *************************************************************************
916 : :
917 [ - + ]: 28371 : 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 : 28371 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
923 : :
924 : : !1) Eventually take care of the forces
925 : :
926 [ + + ]: 28371 : 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 [ + + ]: 28371 : 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 [ + + + + : 28360 : if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then
+ + ]
1001 : 28014 : cplex=optreal
1002 : 28014 : qphon(:)=zero
1003 : : ! Simple scalar multiplication, or model dielectric function
1004 : 28014 : 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 [ + + + - ]: 28371 : 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 [ + + + + ]: 28371 : 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 [ - + ]: 28371 : 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 : 54878 : 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 : 54878 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1395 : 54878 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1396 : : real(dp) :: gmet(3,3),gprimd(3,3),potg0(4),rmet(3,3)
1397 : 54878 : real(dp),allocatable :: gq(:,:),work1(:,:),work2(:)
1398 : : ! *************************************************************************
1399 : :
1400 : : !Check that cplex has an allowed value
1401 : 54878 : 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 [ + + - + ]: 54878 : 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 : 54878 : qeq0=0
1414 [ + + ]: 54878 : 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 [ + + - + ]: 54878 : 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 : 54878 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1423 : 54878 : me_fft=ngfft(11)
1424 : 54878 : nproc_fft=ngfft(10)
1425 : :
1426 : : !Get the distrib associated with this fft_grid
1427 : 54878 : 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 : 54878 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1431 : :
1432 : 54878 : dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ; diemixmag=dielar(7)
1433 : :
1434 : 54878 : magn_precon=(diemixmag>=zero) ! Set to true if magnetization has to be preconditionned
1435 : 54878 : diemixmag=abs(diemixmag)
1436 : :
1437 : : !write(std_out,*)' moddiel : diemac, diemix, diemixmag =',diemac,diemix,diemixmag
1438 : :
1439 [ + + ]: 54878 : 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 [ + + + + : 1753168 : 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 : 48060 : length2=(two_pi*dielng)**2
1453 : 48060 : diemac_inv=1.0_dp/diemac
1454 [ - + + - : 144180 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
1455 [ + + ]: 48060 : if (optreal==1) then
1456 [ - + + - : 115314 : ABI_MALLOC(work2,(cplex*nfft))
+ - + - -
+ - + ]
1457 : : end if
1458 : :
1459 : : ! In order to speed the routine, precompute the components of g
1460 [ + + ]: 913140 : mg=maxval(ngfft)
1461 [ - + + - : 144180 : ABI_MALLOC(gq,(3,mg))
+ - + - -
+ + - ]
1462 [ + + ]: 192240 : do ii=1,3
1463 : 144180 : id(ii)=ngfft(ii)/2+2
1464 [ + + ]: 3057395 : do ing=1,ngfft(ii)
1465 : 2865155 : ig=ing-(ing/id(ii))*ngfft(ii)-1
1466 : 3009335 : 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 [ + + + + ]: 48060 : nspden_eff=nspden;if (optres==1.and.(.not.magn_precon)) nspden_eff=1
1473 [ + + ]: 107804 : do ispden=1,nspden_eff
1474 : :
1475 [ + + ]: 59744 : diemix_eff=diemix;if (ispden>1) diemix_eff=diemixmag
1476 : :
1477 : : ! Do fft from real space (work2) to G space (work1)
1478 [ + + ]: 59744 : if (optreal==1) then
1479 [ + + ]: 617997607 : work2(:)=vresid(:,ispden)
1480 : 40415 : 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 [ + + ]: 240385048 : do ifft=1,nfft
1485 : 240365719 : work1(1,ifft)=vresid(2*ifft-1,ispden)
1486 : 240385048 : work1(2,ifft)=vresid(2*ifft ,ispden)
1487 : : end do
1488 : : end if
1489 : :
1490 : : ! Store G=0 value
1491 : 59744 : potg0(ispden)=work1(re,1)
1492 : :
1493 : : ! Triple loop, for the three dimensions
1494 [ + + ]: 1300700 : do i3=1,n3
1495 : : ! Precompute some products that do not depend on i2 and i1
1496 : 1240956 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
1497 : 1240956 : gqgm23=gq(3,i3)*gmet(2,3)*2
1498 : 1240956 : gqgm13=gq(3,i3)*gmet(1,3)*2
1499 [ + + ]: 30174557 : do i2=1,n2
1500 [ + + ]: 30114813 : if (fftn2_distrib(i2)==me_fft) then
1501 : 27282521 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
1502 : 27282521 : gqgm12=gq(2,i2)*gmet(1,2)*2
1503 : 27282521 : gqg2p3=gqgm13+gqgm12
1504 : 27282521 : 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 : 27282521 : ii1=1
1509 [ + + + + ]: 27282521 : 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 : 27282521 : ii1=2
1513 : : end if
1514 : :
1515 : : ! Here, unlike in hartre.f, the G=0 term is not eliminated, albeit
1516 : : ! not changed.
1517 [ + + ]: 856183872 : do i1=ii1,n1
1518 : :
1519 : : ! One obtains the square of the norm of q+G (defined by i1,i2,i3)
1520 : 828901351 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
1521 : 828901351 : ifft=i1+i23
1522 : :
1523 : 828901351 : l2g2=length2*gs
1524 : : ! The model dielectric function is now computed
1525 : 828901351 : factor = (l2g2+diemac_inv)/(l2g2+1.0_dp) * diemix_eff
1526 : 828901351 : work1(re,ifft)=work1(re,ifft)*factor
1527 : 856183872 : 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 [ + + ]: 107804 : if (optreal==1) then
1542 : 40415 : call fourdp(cplex,work1,work2,1,mpi_enreg,nfft,1,ngfft,0)
1543 [ + + ]: 617997607 : 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 [ + + ]: 240385048 : do ifft=1,nfft
1548 : 240365719 : vrespc(2*ifft-1,ispden)=work1(1,ifft)
1549 : 240385048 : 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 [ - + ]: 48060 : ABI_FREE(gq)
1557 [ - + ]: 48060 : ABI_FREE(work1)
1558 [ + + ]: 48060 : if (optreal==1) then
1559 [ - + ]: 33627 : ABI_FREE(work2)
1560 : : end if
1561 : :
1562 : : ! End condition diemac/=1.0
1563 : : end if
1564 : :
1565 [ - + - + ]: 109756 : 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 : : !!***
|
