Line data Source code
1 : !!****m* ABINIT/m_paw_uj
2 : !! NAME
3 : !! m_paw_uj
4 : !!
5 : !! FUNCTION
6 : !! This module contains several routines relevant only for automatic determination of U
7 : !! in PAW+U context (linear response method according to Phys. Rev. B 71, 035105)
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2018-2026 ABINIT group (DJA)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_paw_uj
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use m_linalg_interfaces
29 : use m_xmpi
30 : use m_dtset
31 : use netcdf
32 : use m_nctk
33 :
34 : use m_fstrings, only : strcat
35 : use m_pptools, only : prmat
36 : use m_special_funcs, only : iradfnh
37 : use m_geometry, only : shellstruct,ioniondist
38 : use m_parser, only : prttagm
39 : use m_supercell, only : mksupercell
40 : use m_pawrad, only : pawrad_type
41 : use m_pawtab, only : pawtab_type
42 : use m_paw_ij, only : paw_ij_type
43 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
44 : use m_dtfil, only : datafiles_type
45 : use m_crystal, only : crystal_t
46 :
47 : implicit none
48 :
49 : private
50 : !!***
51 :
52 : !----------------------------------------------------------------------
53 :
54 : !!****t* m_paw_uj/macro_uj_type
55 : !! NAME
56 : !! dtmacro_uj
57 : !!
58 : !! FUNCTION
59 : !! This data type contains the potential shifts and the occupations
60 : !! for the determination of U and J for the DFT+U calculations.
61 : !! iuj=1,3: non-selfconsistent calculations. iuj=2,4 selfconsistent calculations.
62 : !!
63 : !! SOURCE
64 :
65 : type, public :: macro_uj_type
66 :
67 : ! Integer
68 : integer :: iuj ! dataset treated
69 : integer :: nat ! number of atoms U (J) is determined on
70 : integer :: ndtset ! total number of datasets
71 : integer :: nspden ! number of densities treated
72 : integer :: macro_uj ! which mode the determination runs in
73 : integer :: pawujat ! which atom U (J) is determined on
74 : integer :: pawprtvol ! controlling amount of output
75 : integer :: option ! controls the determination of U (1 with compensating charge bath, 2 without)
76 : integer :: dmatpuopt ! controls the renormalisation of the PAW projectors
77 :
78 : ! Real
79 : real(dp) :: diemix ! mixing parameter
80 : real(dp) :: diemixmag ! magnetic mixing parameter
81 : real(dp) :: mdist ! maximal distance of ions
82 : real(dp) :: pawujga ! gamma for inversion of singular matrices
83 : real(dp) :: ph0phiint ! integral of phi(:,1)*phi(:,1)
84 : real(dp) :: pawujrad ! radius to which U should be extrapolated.
85 : real(dp) :: pawrad ! radius of the paw atomic data
86 :
87 : ! Integer arrays
88 : integer , allocatable :: scdim(:)
89 : ! size of supercell
90 :
91 : ! Real arrays
92 : real(dp) , allocatable :: occ(:,:)
93 : ! occupancies after a potential shift: occ(ispden,nat)
94 :
95 : real(dp) , allocatable :: rprimd(:,:)
96 : ! unit cell for symmetrization
97 :
98 : real(dp) , allocatable :: vsh(:,:)
99 : ! potential shifts on atoms, dimensions: nspden,nat
100 :
101 : real(dp) , allocatable :: xred(:,:)
102 : ! atomic position for symmetrization
103 :
104 : real(dp) , allocatable :: wfchr(:)
105 : ! wfchr(1:3): zion, n and l of atom on which projection is done
106 : ! wfchr(4:6): coefficients ai of a0+a1*r+a2*r^2, fit to the wfc for r< r_paw
107 :
108 : real(dp), allocatable :: zioneff(:)
109 : ! zioneff(ij_proj), "Effective charge"*n "seen" at r_paw, deduced from Phi at r_paw, n:
110 : ! pricipal quantum number; good approximation to model wave function outside PAW-sphere through
111 :
112 : end type macro_uj_type
113 : !!***
114 :
115 : !public procedures.
116 : public :: pawuj_ini ! Initialize dtpawuj datastructure
117 : public :: pawuj_free ! Deallocate dtpawuj datastructure
118 : public :: pawuj_det ! Determine U (or J) parameter
119 : public :: pawuj_red ! Store atomic occupancies, potential shift, positions
120 :
121 : !private procedures.
122 : ! chiscwrt: distributes values of chi_org on chi_sc according to ion-ion distances
123 : ! linvmat: inverts real matrix inmat
124 : ! lprtmat: prints out the real matrix mmat
125 : ! lcalcu: prints out real the real matrice mmat
126 : ! blow_pawuj: reads a real nxn matrice and appends lines n+1 and clumn n+1
127 :
128 : CONTAINS !========================================================================================
129 : !!***
130 :
131 : !----------------------------------------------------------------------
132 :
133 : !!****f* m_paw_uj/pawuj_ini
134 : !! NAME
135 : !! pawuj_ini
136 : !!
137 : !! FUNCTION
138 : !! Initialize dtpawuj datastructure for one SCF run
139 : !! Relevant only for automatic determination of U in PAW+U context
140 : !!
141 : !! INPUTS
142 : !!
143 : !! OUTPUT
144 : !!
145 : !! SIDE EFFECTS
146 : !! dtpawuj(0:ndtpawuj) (initialization of fields vsh, occ, iuj,nnat)
147 : !!
148 : !! SOURCE
149 :
150 10 : subroutine pawuj_ini(dtpawuj,ndtset)
151 :
152 : !Arguments ------------------------------------
153 : !scalars
154 : integer :: ndtset
155 : type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtset)
156 :
157 : !Local variables -------------------------
158 : !Variables for partial dos calculation
159 : !scalars
160 : integer, parameter :: nwfchr=6
161 : integer :: iuj,im1
162 : ! *********************************************************************
163 :
164 : DBG_ENTER("COLL")
165 60 : do iuj=0,ndtset
166 : !write(std_out,*)'pawuj_ini iuj ',iuj
167 50 : dtpawuj(iuj)%diemix=0.45_dp
168 50 : dtpawuj(iuj)%diemixmag=0.45_dp
169 50 : dtpawuj(iuj)%iuj=0
170 50 : dtpawuj(iuj)%nat=0
171 50 : dtpawuj(iuj)%ndtset=1
172 50 : dtpawuj(iuj)%nspden=1
173 50 : dtpawuj(iuj)%macro_uj=0
174 50 : dtpawuj(iuj)%option=1
175 50 : dtpawuj(iuj)%pawujat=1
176 50 : dtpawuj(iuj)%pawujga=one
177 50 : dtpawuj(iuj)%pawprtvol=1
178 50 : dtpawuj(iuj)%ph0phiint=one
179 50 : dtpawuj(iuj)%dmatpuopt=2
180 50 : dtpawuj(iuj)%pawujrad=3.0_dp
181 50 : dtpawuj(iuj)%pawrad=20.0_dp
182 : !Allocate arrays
183 : !write(std_out,*)'pawuj_ini before arrays'
184 50 : ABI_MALLOC(dtpawuj(iuj)%rprimd,(3,3))
185 50 : ABI_MALLOC(dtpawuj(iuj)%scdim,(3))
186 50 : ABI_MALLOC(dtpawuj(iuj)%wfchr,(nwfchr))
187 700 : dtpawuj(iuj)%rprimd=reshape((/ 1,0,0,0,1,0,0,0,1/),(/ 3,3 /))
188 250 : dtpawuj(iuj)%scdim=reshape((/ 250,0,0 /),(/3 /))
189 400 : dtpawuj(iuj)%wfchr=(/ (0,im1=1,nwfchr) /)
190 60 : if (iuj>0) then
191 40 : dtpawuj(iuj)%iuj=-1
192 : end if
193 : end do
194 :
195 : DBG_EXIT("COLL")
196 :
197 10 : end subroutine pawuj_ini
198 : !!***
199 :
200 : !----------------------------------------------------------------------
201 :
202 : !!****f* m_paw_uj/pawuj_free
203 : !! NAME
204 : !! pawuj_free
205 : !!
206 : !! FUNCTION
207 : !! deallocate pawuj stuff
208 : !!
209 : !! INPUTS
210 : !!
211 : !! OUTPUT
212 : !!
213 : !! SOURCE
214 :
215 50 : subroutine pawuj_free(dtpawuj)
216 :
217 : !Arguments -------------------------------
218 : type(macro_uj_type),intent(inout) :: dtpawuj
219 :
220 : ! *********************************************************************
221 :
222 50 : if (allocated(dtpawuj%scdim)) then
223 50 : ABI_FREE(dtpawuj%scdim)
224 : end if
225 50 : if (allocated(dtpawuj%occ)) then
226 50 : ABI_FREE(dtpawuj%occ)
227 : end if
228 50 : if (allocated(dtpawuj%rprimd)) then
229 50 : ABI_FREE(dtpawuj%rprimd)
230 : end if
231 50 : if (allocated(dtpawuj%vsh)) then
232 50 : ABI_FREE(dtpawuj%vsh)
233 : end if
234 50 : if (allocated(dtpawuj%xred)) then
235 50 : ABI_FREE(dtpawuj%xred)
236 : end if
237 50 : if (allocated(dtpawuj%wfchr)) then
238 50 : ABI_FREE(dtpawuj%wfchr)
239 : end if
240 50 : if (allocated(dtpawuj%zioneff)) then
241 0 : ABI_FREE(dtpawuj%zioneff)
242 : end if
243 :
244 50 : end subroutine pawuj_free
245 : !!***
246 :
247 : !----------------------------------------------------------------------
248 : !!****f* m_paw_uj/pawuj_det
249 : !! NAME
250 : !! pawuj_det
251 : !!
252 : !! FUNCTION
253 : !! From the complete dtpawuj-dataset determines U (or J) parameter for
254 : !! PAW+U calculations
255 : !! Relevant only for automatic determination of Hubbard Parameters in
256 : !! PAW+U context.
257 : !! Hubbard U = diemix/chi0-1/chi
258 : !! Hund's J = 1/chi-diemixmag/chi0
259 : !!
260 : !! INPUTS
261 : !! dtpawuj=potential shifts (vsh) and atomic occupations (occ)
262 : !!
263 : !! OUTPUT
264 : !! only printing
265 : !! (among other things a section in the ab.out that can be used for input in ujdet)
266 : !!
267 : !! SOURCE
268 :
269 10 : subroutine pawuj_det(dtpawuj,ndtpawuj,dtset,dtfil,ures,comm)
270 :
271 : !Arguments ------------------------------------
272 : !scalars
273 : !arrays
274 : integer :: ndtpawuj, comm
275 : type(macro_uj_type),intent(in) :: dtpawuj(0:ndtpawuj)
276 : type(dataset_type),intent(in) :: dtset
277 : type(datafiles_type),intent(in) :: dtfil
278 : real(dp),intent(out) :: ures
279 :
280 : !Local variables-------------------------------
281 : !scalars
282 : integer,parameter :: natmax=2,nwfchr=6
283 : integer :: ii,jj,nat_org,natom,jdtset,nspden,macro_uj,marr
284 : integer :: im1,ndtuj,nsh_org, nsh_sc,nat_sc,maxnat
285 : integer :: pawujat,pawprtvol,pawujoption
286 : integer :: dmatpuopt,invopt,ipert
287 : integer :: my_rank, ncid, ncerr
288 : real(dp) :: pawujga,ph0phiint,intg,fcorr,eyp
289 : real(dp) :: diem,signum,scalarHP !LMac quantities
290 :
291 : character(len=500) :: message
292 : !character(len=2) :: hstr
293 : character(len=5) :: pertname
294 : character(len=1) :: parname
295 : character(len=14) :: occmag
296 520 : type(crystal_t) :: cryst
297 : !arrays
298 : integer :: ext(3)
299 20 : real(dp) :: rprimd_sc(3,3),vsh(ndtpawuj),a(5),b(5)
300 10 : integer,allocatable :: jdtset_(:),smult_org(:),smult_sc(:),intmagv_org(:),intmagv_sc(:)
301 10 : real(dp),allocatable :: luocc(:,:),xred_org(:,:)
302 10 : real(dp),allocatable :: magv_org(:),magv_sc(:),chi(:),chi0(:),chi0_sc(:), chi_sc(:), xred_sc(:,:)
303 10 : real(dp),allocatable :: sdistv_org(:),sdistv_sc(:),distv_org(:),distv_sc(:)
304 : integer,parameter :: ncid0 = 0, master = 0
305 : integer :: units(2)
306 : ! *********************************************************************
307 :
308 : DBG_ENTER("COLL")
309 :
310 10 : my_rank = xmpi_comm_rank(comm)
311 30 : units = [std_out, ab_out]
312 :
313 : !write(std_out,*) 'pawuj 01'
314 : !###########################################################
315 : !### 01. Allocations
316 :
317 : !Initializations
318 60 : ndtuj=count(dtpawuj(:)%iuj/=-1)-1 ! number of datasets initialized by pawuj_red
319 30 : ABI_MALLOC(jdtset_,(0:ndtuj))
320 110 : jdtset_(0:ndtuj)=pack(dtpawuj(:)%iuj,dtpawuj(:)%iuj/=-1)
321 60 : jdtset=maxval(dtpawuj(:)%iuj)
322 :
323 : !DEBUG
324 10 : write(message,'(10(a,i3))')'pawuj_det jdtset ',jdtset,&
325 20 : & ' ndtuj ', ndtuj,' ndtpawuj ',ndtpawuj
326 10 : call wrtout(std_out,message,'COLL')
327 : !call flush_unit(6)
328 : !END DEBUG
329 :
330 10 : nspden=dtpawuj(jdtset)%nspden
331 10 : nat_org=dtpawuj(jdtset)%nat
332 10 : natom=dtset%natom
333 10 : write(message,'(a,i3)')'natom,nat_org',nat_org
334 10 : call wrtout(std_out,message,'COLL')
335 : ! ABI_CHECK(nat_org == 1, "MG: I'm not sure we access data with the right index if nat_org > 1")
336 :
337 10 : macro_uj=dtpawuj(jdtset)%macro_uj
338 10 : pawujat=dtpawuj(jdtset)%pawujat
339 10 : pawprtvol=dtpawuj(jdtset)%pawprtvol
340 10 : pawujga=dtpawuj(jdtset)%pawujga
341 10 : pawujoption=dtpawuj(jdtset)%option
342 10 : ph0phiint=dtpawuj(jdtset)%ph0phiint
343 10 : dmatpuopt=dtpawuj(jdtset)%dmatpuopt
344 10 : marr=maxval((/ 9, nspden*nat_org ,nat_org*3 /))
345 10 : eyp=2.5_dp ! for dmatpuopt==2 and 3
346 10 : if (dmatpuopt==1) eyp=eyp+3.0_dp
347 10 : if (dmatpuopt>=3) eyp=(eyp+3.0_dp-dmatpuopt)
348 :
349 : !DEBUG
350 10 : write(message,'(a,i3)')'pawujdet dmatpuopt',dmatpuopt
351 10 : call wrtout(std_out,message,'COLL')
352 : !END DEBUG
353 :
354 40 : ABI_MALLOC(luocc,(ndtpawuj,nat_org))
355 30 : ABI_MALLOC(magv_org,(nat_org))
356 30 : ABI_MALLOC(intmagv_org,(nat_org))
357 30 : ABI_MALLOC(xred_org,(3,nat_org))
358 20 : ABI_MALLOC(chi0,(nat_org))
359 20 : ABI_MALLOC(chi,(nat_org))
360 20 : ABI_MALLOC(distv_org,(nat_org))
361 : !DEBUG
362 : !write(message,fmt='((a,i3,a))')'pawuj_det init sg'
363 : !call wrtout(std_out,message,'COLL')
364 : !END DEBUG
365 85 : luocc=zero
366 :
367 : !###########################################################
368 : !### 02. Testing consistency of parameters and outputting info
369 :
370 : !LMac
371 : ! if (ndtuj/=4) return
372 :
373 10 : write(message, '(3a)' ) ch10,' ---------- calculate U, (J) start ---------- ',ch10
374 10 : call wrtout(ab_out,message,'COLL')
375 :
376 : !Tests if perturbed atom is consistent with ujdet atom
377 50 : if (all(dtpawuj(1:ndtpawuj)%pawujat==pawujat)) then
378 10 : write (message,fmt='(a,i3)') ' All pawujat ok and equal to ',pawujat
379 10 : call wrtout(ab_out,message,'COLL')
380 : else
381 0 : write (message,fmt='(a,4i3,2a)') ' Differing values of pawujat were found: ',dtpawuj(1:ndtuj)%pawujat,ch10,&
382 0 : & 'No determination of U.'
383 0 : call wrtout(ab_out,message,'COLL')
384 0 : return
385 : end if
386 :
387 : !Tests consistency of macro_uj, then writes message about macro_uj procedure selected. LMac
388 10 : if (nspden==1) then
389 0 : write(message,fmt='(2a)') ' nspden==1, determination',&
390 0 : & ' of U-parameter for unpolarized structure (non standard)'
391 10 : else if (macro_uj==1.and.nspden==2) then
392 6 : write(message,fmt='(2a)') ' macro_uj=1 and nspden=2:',&
393 12 : & ' standard determination of Hubbard U-parameter'
394 4 : else if (macro_uj==2.and.nspden==2) then
395 0 : write(message,fmt='(2a)') ' macro_uj=2 and nspden=2:',&
396 0 : & ' determination of parameter on single spin channel (experimental)'
397 4 : else if (macro_uj==3.and.nspden==2) then
398 0 : write(message,fmt='(2a)') ' macro_uj=3 and nspden=2,',&
399 0 : & ' determination of (not Hunds) J-parameter on single spin channel (experimental)'
400 4 : else if (macro_uj==4.and.nspden==2) then
401 4 : write(message,fmt='(2a)') ' macro_uj=4 and nspden=2,',&
402 8 : & ' Hunds J determination -- L. MacEnulty August 2021'
403 : end if
404 10 : call wrtout(ab_out,message,'COLL')
405 :
406 : !Tests compatibility of nspden and macro_uj
407 10 : if (macro_uj>1.and.nspden==1) then
408 0 : write (message,'(4a,2a)') ' U on a single spin channel (or J) can only be determined for nspden=2 ,',ch10,&
409 0 : & 'No determination of U.'
410 0 : call wrtout(ab_out,message,'COLL')
411 0 : return
412 : end if
413 :
414 : !Calculation of response matrix
415 10 : diem=dtpawuj(1)%diemix !Unscreened response in Hubbard U impacted by diemix
416 10 : pertname='alpha' !Hubbard U perturbation; applied equally to spins up and down
417 10 : parname='U'
418 10 : occmag=' Occupations'
419 10 : signum=1.0d0 !Hubbard U is signum*(1/chi0-1/chi)
420 50 : do jdtset=1,4
421 40 : if (nspden==1) then
422 0 : luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)
423 : !Hubbard U: uses sum of spin up and spin down occupancies
424 40 : else if (macro_uj==1.and.nspden==2) then
425 52 : luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)+dtpawuj(jdtset)%occ(2,:) !Total occupation
426 16 : else if (macro_uj==2.and.nspden==2) then
427 0 : luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)
428 16 : else if (macro_uj==3.and.nspden==2) then
429 0 : luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(2,:)
430 0 : parname='J'
431 : !Hund's J: uses difference of spin up and spin down occupancies
432 16 : else if (macro_uj==4.and.nspden==2) then
433 48 : luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)-dtpawuj(jdtset)%occ(2,:) !Magnetization
434 16 : diem=dtpawuj(1)%diemixmag !Unscreened response in Hund's J impacted by diemixmag
435 16 : pertname='beta ' !Hund's J perturbation: +beta to spin up, -beta to down
436 16 : parname='J'
437 16 : occmag='Magnetizations'
438 16 : signum=-1.0d0 !Hund's J is -1*(1/chi0-1/chi)
439 : end if
440 : ! MG: Cannot use pawujat index to extract vsh as we have used:
441 : !
442 : ! dtpawuj(icyc)%vsh=reshape(pack(atvshift,atvshmusk),(/ nspden,nnat /))
443 :
444 : ! thus the pawujat atom in the atvshift array becomes the first one in %vsh
445 : !
446 : !vsh(jdtset)=dtpawuj(jdtset)%vsh(1,pawujat)
447 40 : vsh(jdtset)=dtpawuj(jdtset)%vsh(1,1)
448 50 : if (pawprtvol==-3) then
449 0 : write(message,fmt='(2a,i3,a,f15.12)') ch10,' Potential shift vsh(',jdtset,') =',vsh(jdtset)
450 0 : call wrtout(std_out,message,'COLL')
451 0 : write(message,fmt='( a,i3,a,120f15.9)') ' Occupations occ(',jdtset,') ',( luocc(jdtset,ii),ii=1,nat_org )
452 0 : call wrtout(std_out,message,'COLL')
453 : end if
454 : end do
455 :
456 10 : write(message,fmt='( a)') 'Occupations assigned.'
457 10 : call wrtout(std_out,message,'COLL')
458 :
459 : !Two-point linear regression of response matrices.
460 35 : chi0=(luocc(1,1:nat_org)-luocc(3,1:nat_org))/(vsh(1)-vsh(3))/diem
461 35 : chi=(luocc(2,1:nat_org)-luocc(4,1:nat_org))/(vsh(2)-vsh(4))
462 :
463 : ! MG: pawujat replaced by 1 because arrays are dimensioned with nat_org (usually 1) and not natom!
464 10 : if ((abs(chi0(1))<0.0000001).or.(abs(chi(1))<0.0000001)) then
465 0 : write(message, '(2a,2f12.5,a)' ) ch10,'Chi0 or Chi is too small for inversion.',&
466 0 : &chi0(1),chi(1),ch10
467 0 : call wrtout(ab_out,message,'COLL')
468 0 : return
469 : end if
470 :
471 : !LMac: Scalar Hubbard Parameter
472 10 : scalarHP=signum*(1.0d0/chi0(1)-1.0d0/chi(1))*Ha_eV
473 :
474 10 : write(message,fmt='(a)')': '
475 10 : if (nspden==2) then
476 35 : magv_org=dtpawuj(1)%occ(1,:)-dtpawuj(1)%occ(2,:)
477 10 : if (all(abs(magv_org)<0.001)) then
478 0 : magv_org=(/(1,im1=1,nat_org)/)
479 : else
480 25 : magv_org=abs(magv_org)/magv_org
481 10 : if (magv_org(1).lt.0) magv_org=magv_org*(-1_dp)
482 10 : if (all(magv_org(2:nat_org).lt.0)) then
483 10 : magv_org=abs(magv_org)
484 5 : write(message,'(a)')', (reset to fm): '
485 : end if
486 : end if
487 : else
488 0 : magv_org=(/(1,im1=1,nat_org)/)
489 : end if
490 65 : intmagv_org=(/(int(magv_org(im1)),im1=1,nat_org)/)
491 :
492 : ! if (pawprtvol==-3) then
493 10 : write(message,fmt='(3a, 150f4.1)') ch10,' Magnetization',trim(message),magv_org
494 10 : call wrtout(std_out,message,'COLL')
495 : ! end if
496 :
497 : !Case of extrapolation to larger r_paw: calculate intg
498 70 : if (all(dtpawuj(1)%wfchr(:)/=0).and.ph0phiint/=1) then
499 10 : if (dtpawuj(1)%pawujrad<20.and.dtpawuj(1)%pawujrad>dtpawuj(1)%pawrad) then
500 : fcorr=(1-ph0phiint)/(IRadFnH(dtpawuj(1)%pawrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),&
501 10 : & nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1)))
502 : intg=ph0phiint/(1-fcorr*IRadFnH(dtpawuj(1)%pawujrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),&
503 10 : & nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1)))
504 10 : write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met2 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg
505 10 : call wrtout(std_out,message,'COLL')
506 0 : else if (dtpawuj(1)%pawujrad<dtpawuj(1)%pawrad) then
507 0 : a=0 ; a(1:3)=dtpawuj(1)%wfchr(4:6)
508 0 : a=(/a(1)**2,a(1)*a(2),a(2)**2/3.0_dp+2.0_dp/3.0_dp*a(1)*a(3),a(2)*a(3)/2.0_dp,a(3)**2/5.0_dp/)
509 0 : b=(/(dtpawuj(1)%pawujrad**(im1)-dtpawuj(1)%pawrad**(im1),im1=1,5)/)
510 0 : intg=dot_product(a,b)
511 0 : intg=ph0phiint/(ph0phiint+intg)
512 0 : write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met1 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg
513 0 : call wrtout(std_out,message,'COLL')
514 0 : else if (dtpawuj(1)%pawujrad==dtpawuj(1)%pawrad) then
515 0 : intg=one
516 0 : write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (pawujrad=pawrad)'
517 0 : call wrtout(std_out,message,'COLL')
518 : else
519 0 : intg=ph0phiint
520 0 : write(message, fmt='(a,3f12.5)') ' pawuj_det: extrapolation to r->\inf using factor ',intg
521 0 : call wrtout(std_out,message,'COLL')
522 : end if
523 : else
524 0 : write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (ph0phiint=1 or wfchr=0)'
525 0 : call wrtout(std_out,message,'COLL')
526 0 : intg=one
527 : end if
528 :
529 : !Determine U in primitive cell
530 10 : write(message,fmt='(a)')' pawuj_det: determine U in primitive cell'
531 10 : call wrtout(std_out,message,'COLL')
532 :
533 10 : call lcalcu(intmagv_org,nat_org,dtpawuj(1)%rprimd,dtpawuj(1)%xred,chi,chi0,pawujat,ures,pawprtvol,pawujga,pawujoption)
534 : !Begin calculate U in supercell
535 :
536 : !Analyze shell structure of primitive cell
537 : !and atomic distances in primitive cell
538 20 : ABI_MALLOC(smult_org,(nat_org))
539 20 : ABI_MALLOC(sdistv_org,(nat_org))
540 : call shellstruct(dtpawuj(1)%xred,dtpawuj(1)%rprimd,nat_org,&
541 10 : & intmagv_org,distv_org,smult_org,sdistv_org,nsh_org,pawujat,pawprtvol)
542 :
543 :
544 : !LMac: Printing relevant information about the Hubbard parameter just calculated.
545 10 : write(message,'(3a)') ch10,ch10,'*********************************************************************'
546 10 : call wrtout(units,message)
547 10 : write(message,'(4a)') '************************ Linear Response ',parname,' ************************',ch10
548 10 : call wrtout(units,message)
549 10 : write(message, '(a,i4,a)' ) ' Info printed for perturbed atom: ',pawujat,ch10
550 10 : call wrtout(units,message)
551 10 : write(message, fmt='(10a)')' Perturbations ',occmag,ch10,&
552 10 : ' --------------- -----------------------------',ch10,&
553 10 : ' ',pertname,' [eV] Unscreened Screened',ch10,&
554 20 : ' --------------- -----------------------------'
555 10 : call wrtout(units,message)
556 : ! MG: pawujat --> 1.
557 30 : do ipert=1,2
558 20 : write(message, fmt='(3f15.10)') vsh(ipert*2-1)*Ha_eV,luocc(ipert*2-1,1),luocc(ipert*2,1)
559 30 : call wrtout(units,message)
560 : end do
561 10 : write(message,'(2a)') ch10,' Scalar response functions:'
562 10 : call wrtout(units,message)
563 10 : write(message,fmt='(a,f12.5)') ' Chi0 [eV^-1]: ',chi0(1)/Ha_eV
564 10 : call wrtout(units,message)
565 10 : write(message,fmt='(a,f12.5)') ' Chi [eV^-1]: ',chi(1)/Ha_eV
566 10 : call wrtout(units,message)
567 10 : write(message,'(4a,f9.5,a)') ch10,' The scalar ',parname,' from the two-point regression scheme is ',scalarHP,' eV.'
568 10 : call wrtout(units,message)
569 10 : write(message,'(3a)') '*********************************************************************',ch10,&
570 20 : '*********************************************************************'
571 10 : call wrtout(units,message)
572 10 : write(message,'(7a)') 'Note: For more reliable linear regressions of the response',ch10,&
573 10 : 'matrices, it is advised that you have more than two points.',ch10,&
574 20 : 'See the LRUJ protocol for more information.',ch10,ch10
575 10 : call wrtout(std_out,message,'COLL')
576 10 : call wrtout(ab_out,message,'COLL')
577 :
578 : !###########################################################
579 : !### Create the file LRUJ.nc for the LRUJ ujdet utility
580 :
581 10 : if (my_rank == master) then
582 : ! First call:
583 : ! - Create NC file, define dimensions, scalars and arrays.
584 : ! - Add crystal structure and metadata required to post-process the data.
585 10 : NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), "_LRUJ.nc"), xmpi_comm_self))
586 :
587 10 : cryst = dtset%get_crystal(1)
588 10 : NCF_CHECK(cryst%ncwrite(ncid))
589 10 : call cryst%free()
590 :
591 : !Define dimensions.
592 : ncerr = nctk_def_dims(ncid, [ &
593 : nctkdim_t("natom", natom), &
594 : nctkdim_t("nnat", nat_org), &
595 : nctkdim_t("ndtpawuj", ndtpawuj), &
596 : nctkdim_t("nspden", dtset%nspden), &
597 : nctkdim_t("nsppol", dtset%nsppol) ], &
598 60 : defmode=.True.)
599 10 : NCF_CHECK(ncerr)
600 :
601 : ! Define integer scalars
602 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
603 : "usepaw", "macro_uj", "pawujat", "nspden", "dmatpuopt" &
604 60 : ])
605 :
606 : ! Define double precision scalars
607 : ! @lmacenul: Here I write pawujv so that one can order the results by alpha in the post-processing tool.
608 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
609 : "diemix", "diemixmag", "ph0phiint", "uj_pert" &
610 50 : ])
611 10 : NCF_CHECK(ncerr)
612 :
613 : ! Define arrays with results.
614 : ! TODO: Here I need an extra dimension to store results with different iuj and/or different names.
615 : ncerr = nctk_def_arrays(ncid, [ &
616 : nctkarr_t("luocc", "dp", "ndtpawuj, nnat") &
617 : !nctkarr_t("vsh", "dp", "nspden, nnat") &
618 20 : ])
619 10 : NCF_CHECK(ncerr)
620 :
621 : ! ===========================================
622 : ! Write metadata that does not depend on icyc
623 : ! ===========================================
624 10 : NCF_CHECK(nctk_set_datamode(ncid))
625 :
626 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
627 : "usepaw", "macro_uj", "pawujat", "nspden", "dmatpuopt"], &
628 110 : [dtset%usepaw, macro_uj, pawujat, nspden, dmatpuopt])
629 10 : NCF_CHECK(ncerr)
630 :
631 : associate (dt => dtpawuj(1))
632 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
633 : "diemix", "diemixmag", "ph0phiint", "uj_pert" ], &
634 90 : [dt%diemix, dt%diemixmag, ph0phiint, vsh(3)])
635 20 : NCF_CHECK(ncerr)
636 : end associate
637 :
638 : ! Write arrays
639 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "luocc"), luocc))
640 :
641 10 : NCF_CHECK(nf90_close(ncid))
642 : end if
643 :
644 10 : ii=1
645 10 : write(message, fmt='(8a)') ' URES ',' ii',' nat',' r_max',' U(J)[eV]',' U_ASA[eV]',' U_inf[eV]',ch10
646 : write(message, fmt='(a,2i7,4f12.5)') &
647 35 : trim(message)//' URES ',ii,nat_org,maxval(abs(distv_org)),signum*ures,signum*ures*exp(log(intg)*eyp),&
648 20 : signum*ures*exp(log(ph0phiint)*eyp)
649 10 : call wrtout(units,message)
650 :
651 10 : if (pawprtvol>1) then
652 0 : write(message,fmt='(a, 150f10.5)')' pawuj_det: ionic distances in original cell (distv_org) ', distv_org
653 0 : call wrtout(std_out,message,'COLL')
654 : end if
655 :
656 : !Construct supercell, calculate limit dimensions of supercell
657 : ii=1
658 40 : maxnat=product(dtpawuj(1)%scdim(:))*nat_org
659 10 : if (maxnat==0) then
660 10 : maxnat=dtpawuj(1)%scdim(1)
661 40 : ext=(/ii, ii, ii/)
662 : else
663 : jj=1
664 0 : do while (jj<=3)
665 0 : ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) )
666 0 : jj=jj+1
667 : end do
668 : end if
669 40 : ext=ext+(/ 1,1,1 /)
670 10 : ii=ii+1
671 :
672 :
673 40 : nat_sc=product(ext)*nat_org
674 :
675 : !DEBUG
676 : !write(message,fmt='(a,3i4)')'pawuj_det: ext ',ext
677 : !call wrtout(std_out,message,'COLL')
678 : !END DEBUG
679 :
680 55 : do while (nat_sc<=maxnat)
681 135 : ABI_MALLOC(chi0_sc,(nat_sc))
682 90 : ABI_MALLOC(chi_sc,(nat_sc))
683 90 : ABI_MALLOC(distv_sc,(nat_sc))
684 90 : ABI_MALLOC(magv_sc,(nat_sc))
685 135 : ABI_MALLOC(intmagv_sc,(nat_sc))
686 90 : ABI_MALLOC(sdistv_sc,(nat_sc))
687 90 : ABI_MALLOC(smult_sc,(nat_sc))
688 135 : ABI_MALLOC(xred_sc,(3,nat_sc))
689 :
690 : ! Determine positions=xred_sc and supercell dimensions=rpimd_sc
691 : call mksupercell(dtpawuj(1)%xred,intmagv_org,dtpawuj(1)%rprimd,nat_org,&
692 45 : & nat_sc,xred_sc,magv_sc,rprimd_sc,ext,pawprtvol)
693 :
694 : ! Determine shell structure of supercell: multiplicities (smult_sc), radii of shells (sdistv_sc)
695 : ! number of shells (nsh_sc) and atomic distances in supercell (distv_sc)
696 :
697 45 : write(message,fmt='(a,3i3,a)')' pawuj_det: determine shell structure of ',ext(1:3),' supercell'
698 45 : call wrtout(std_out,message,'COLL')
699 :
700 13410 : intmagv_sc=(/(int(magv_sc(im1)),im1=1,nat_sc)/)
701 45 : call shellstruct(xred_sc,rprimd_sc,nat_sc,intmagv_sc,distv_sc,smult_sc,sdistv_sc,nsh_sc,pawujat,pawprtvol)
702 :
703 45 : if (pawprtvol>=2) then
704 0 : write(message,fmt='(a)')' pawuj_det: ionic distances in supercell (distv_sc) '
705 0 : call wrtout(std_out,message,'COLL')
706 0 : call prmat(distv_sc,1,nat_sc,1,std_out)
707 : end if
708 :
709 : ! Determine chi and chi0 in supercell (chi0_sc, chi_sc)
710 : ! DEBUG
711 : ! write(message,fmt='(a)')'pawuj_det: chi and chi0 in supercell'
712 : ! call wrtout(std_out,message,'COLL')
713 : ! END DEBUG
714 :
715 45 : if (pawujoption>2) then
716 0 : invopt=2
717 : else
718 45 : invopt=1
719 : end if
720 :
721 : call chiscwrt(chi0,distv_org,nat_org,sdistv_org,smult_org,nsh_org,&
722 45 : & chi0_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol)
723 : call chiscwrt(chi,distv_org,nat_org,sdistv_org,smult_org,nsh_org,&
724 45 : & chi_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol)
725 :
726 : ! Calculate U in supercell
727 : ! DEBUG
728 : ! write(message,fmt='(a)')'pawuj_det: U in supercell'
729 : ! call wrtout(std_out,message,'COLL')
730 : ! END DEBUG
731 45 : call lcalcu(intmagv_sc,nat_sc,rprimd_sc,xred_sc,chi_sc,chi0_sc,pawujat,ures,pawprtvol,pawujga,pawujoption)
732 4530 : write(message, fmt='(a,2i7,4f12.5)') ' URES ',ii,nat_sc,maxval(abs(distv_sc)),signum*ures,signum*ures*exp(log(intg)*eyp),&
733 90 : & signum*ures*exp(log(ph0phiint)*eyp)
734 45 : call wrtout(units,message)
735 :
736 45 : ABI_FREE(chi0_sc)
737 45 : ABI_FREE(chi_sc)
738 45 : ABI_FREE(distv_sc)
739 45 : ABI_FREE(magv_sc)
740 45 : ABI_FREE(intmagv_sc)
741 45 : ABI_FREE(sdistv_sc)
742 45 : ABI_FREE(smult_sc)
743 45 : ABI_FREE(xred_sc)
744 :
745 180 : if (product(dtpawuj(1)%scdim(:))==0) then
746 180 : ext=(/ii, ii, ii/)
747 : else
748 : jj=1
749 0 : do while (jj<=3)
750 0 : ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) )
751 0 : jj=jj+1
752 : end do
753 : end if
754 180 : ext=ext+(/ 1,1,1 /)
755 45 : ii=ii+1
756 :
757 280 : nat_sc=product(ext)*nat_org
758 :
759 : end do
760 :
761 10 : write(message, '(3a)' )ch10,' ---------- calculate U, (J) end -------------- ',ch10
762 10 : call wrtout(ab_out,message,'COLL')
763 :
764 10 : ABI_FREE(jdtset_)
765 10 : ABI_FREE(chi)
766 10 : ABI_FREE(chi0)
767 10 : ABI_FREE(smult_org)
768 10 : ABI_FREE(sdistv_org)
769 10 : ABI_FREE(luocc)
770 10 : ABI_FREE(magv_org)
771 10 : ABI_FREE(intmagv_org)
772 10 : ABI_FREE(xred_org)
773 10 : ABI_FREE(distv_org)
774 :
775 : DBG_EXIT("COLL")
776 :
777 30 : end subroutine pawuj_det
778 : !!***
779 :
780 : !----------------------------------------------------------------------
781 : !!****f* m_paw_uj/pawuj_red
782 : !! NAME
783 : !! pawuj_red
784 : !!
785 : !! FUNCTION
786 : !! Store atomic occupancies, potential shift, positions in dtpawuj datastructure.
787 : !!
788 : !! INPUTS
789 : !! istep: SCF iteration step
790 : !! pert_state: 0 if the routine is called with unperturbed occupancies
791 : !! 1 if the routine is called after having applied the perturbation.
792 : !! fatvshift=factor that multiplies atvshift
793 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
794 : !! comm_atom=--optional-- MPI communicator over atoms
795 : !! my_natom=number of atoms treated by current processor
796 : !! natom=number of atoms in cell
797 : !! ntypat = number of atom types
798 : !! paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
799 : !! pawprtvol= printing volume
800 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
801 : !! comm=MPI commuicator
802 : !!
803 : !! OUTPUT
804 : !! dtpawuj(0:ndtpawuj) (initialization of fields vsh, occ, occ0, iuj,nnat)
805 : !!
806 : !! SOURCE
807 :
808 146 : subroutine pawuj_red(istep, pert_state, dtfil, &
809 292 : dtset,dtpawuj,fatvshift,my_natom,natom,ntypat,paw_ij,pawrad,pawtab,ndtpawuj, comm, &
810 146 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
811 :
812 : !Arguments ------------------------------------
813 : !scalars
814 : integer,intent(in) :: istep, pert_state, my_natom,natom,ntypat,ndtpawuj, comm
815 : integer,optional,intent(in) :: comm_atom
816 : real(dp),intent(in) :: fatvshift
817 : !arrays
818 : integer,optional,target,intent(in) :: mpi_atmtab(:)
819 : type(paw_ij_type),intent(in) :: paw_ij(my_natom)
820 : type(pawtab_type),intent(in) :: pawtab(ntypat)
821 : type(pawrad_type),intent(in) :: pawrad(ntypat)
822 : type(dataset_type),intent(in) :: dtset
823 : type(datafiles_type),intent(in) :: dtfil
824 : type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtpawuj)
825 :
826 : !Local variables-------------------------------
827 : !scalars
828 : integer,parameter :: natmax=2,ncoeff=3, master = 0
829 : integer :: iatom,iatom_tot,ierr,im1,im2,ispden,itypat,ll,nspden,nsppol,iuj,ncyc,icyc
830 : integer :: my_comm_atom,nnat,natpawu,natvshift,pawujat,ndtset,typawujat
831 : !integer :: my_rank, ncid, ncerr
832 : logical :: usepawu !antiferro,
833 : logical :: my_atmtab_allocated,paral_atom
834 : character(len=1000) :: message,hstr
835 : character(len=500) :: messg
836 : !arrays
837 292 : logical :: musk(3,natom)
838 146 : integer,pointer :: my_atmtab(:)
839 : real(dp) :: rrtab(ncoeff),wftab(ncoeff),a(ncoeff,ncoeff),b(ncoeff,ncoeff)! ,a(ncoeff,ncoeff)
840 : real(dp),allocatable :: nnocctot(:,:) !,magv(:)
841 : real(dp),allocatable :: atvshift(:,:,:) ! atvshift(natvshift,2,natom)
842 : logical,allocatable :: dmusk(:,:),atvshmusk(:,:,:) !atvshmusk(natvshift,2,natom)
843 :
844 : ! *********************************************************************
845 :
846 : ABI_UNUSED((/dtfil%ireadwf, istep, pert_state, comm/))
847 :
848 : !Initializations
849 146 : nspden=1;nsppol=1
850 146 : if (my_natom>0) then
851 146 : nspden=paw_ij(1)%nspden ; nsppol=paw_ij(1)%nsppol
852 : end if
853 :
854 : !my_rank = xmpi_comm_rank(comm)
855 :
856 146 : natvshift=dtset%natvshift
857 146 : pawujat=dtset%pawujat
858 146 : natpawu=dtset%natpawu ; ndtset=dtset%ndtset
859 730 : ABI_MALLOC(atvshift,(natvshift,nspden,natom))
860 730 : ABI_MALLOC(atvshmusk,(natvshift,nspden,natom))
861 584 : ABI_MALLOC(dmusk,(nspden,natom))
862 584 : ABI_MALLOC(nnocctot,(nspden,natom))
863 3366 : musk=.false.; dmusk=.false.
864 6272 : atvshift=fatvshift*dtset%atvshift
865 6126 : atvshmusk=.false.
866 1526 : nnocctot=zero
867 146 : typawujat=dtset%typat(pawujat)
868 460 : usepawu=(count(pawtab(:)%usepawu/=0)>0)
869 :
870 : !Set up parallelism over atoms
871 146 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
872 146 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
873 146 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
874 146 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
875 :
876 146 : nnat=0
877 146 : if (usepawu) then
878 146 : write(message,'(3a)') ch10, '---------- pawuj_red ------ ',ch10
879 146 : call wrtout(std_out, message,'COLL');
880 606 : do iatom=1,my_natom
881 460 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
882 460 : itypat=paw_ij(iatom)%itypat;ll=pawtab(itypat)%lpawu
883 606 : if ((ll>=0).and.(pawtab(itypat)%usepawu/=0).and.itypat==typawujat) then
884 832 : musk(:,iatom_tot)=(/.true., .true., .true. /)
885 5200 : atvshmusk(:,:,iatom_tot)=reshape((/ (( (im1==1), im1=1,natvshift) ,im2=1,nspden ) /),(/natvshift,nspden/))
886 624 : do ispden=1,nspden
887 416 : nnocctot(ispden,iatom_tot)=paw_ij(iatom)%nocctot(ispden)
888 624 : dmusk(ispden,iatom_tot)=.true.
889 : end do
890 208 : nnat=nnat+1
891 : end if
892 : end do
893 :
894 : ! Reduction in case of parallelism
895 146 : if (paral_atom) then
896 0 : call xmpi_sum(nnocctot ,my_comm_atom,ierr)
897 0 : call xmpi_lor(atvshmusk,my_comm_atom) ! dim=natpawu ???
898 0 : call xmpi_lor(dmusk ,my_comm_atom)
899 0 : call xmpi_lor(musk ,my_comm_atom)
900 0 : call xmpi_sum(nnat ,my_comm_atom,ierr)
901 : end if
902 :
903 876 : iuj=maxval(dtpawuj(:)%iuj)
904 146 : write(std_out,*)' pawuj_red: iuj',iuj
905 : ! write(std_out,*)' pawuj_red: dtpawuj(:)%iuj ',(dtpawuj(ii)%iuj,ii=1,ndtpawuj)
906 :
907 : !If this is the unperturbed state, then unscreened and screened occupancies
908 : !are the same. Also set perturbation vsh to zero.
909 : !If this is the unperturbed case, then do everything twice; once for iuj=1
910 : !and the same for iuj=2. If this is the perturbed case, do everything only
911 : !once. LMac
912 146 : if (iuj==1) then
913 40 : ABI_MALLOC(dtpawuj(0)%vsh,(nspden,nnat))
914 30 : ABI_MALLOC(dtpawuj(0)%occ,(nspden,nnat))
915 30 : ABI_MALLOC(dtpawuj(0)%xred,(3,nnat))
916 55 : dtpawuj(0)%vsh=0
917 55 : dtpawuj(0)%occ=0
918 70 : dtpawuj(0)%xred=0
919 : ncyc=2
920 : else
921 : ncyc=iuj
922 : end if
923 :
924 302 : do icyc=iuj,ncyc
925 :
926 156 : if (icyc==1) then ! 1 and 3: non-scf steps
927 10 : dtpawuj(icyc+1)%iuj=icyc+1
928 10 : dtpawuj(icyc+2)%iuj=icyc+2
929 146 : else if (icyc==3) then
930 10 : dtpawuj(icyc+1)%iuj=icyc+1
931 : end if
932 :
933 : ! TODO: check that this is correct: this point is passed several times
934 : ! for a given value of iuj - should the stuff be accumulated instead of replaced?
935 156 : if(allocated(dtpawuj(icyc)%vsh)) then
936 116 : ABI_FREE(dtpawuj(icyc)%vsh)
937 : end if
938 624 : ABI_MALLOC(dtpawuj(icyc)%vsh,(nspden,nnat))
939 : !Allocate screened occupancy array: LMac
940 156 : if(allocated(dtpawuj(icyc)%occ)) then
941 116 : ABI_FREE(dtpawuj(icyc)%occ)
942 : end if
943 468 : ABI_MALLOC(dtpawuj(icyc)%occ,(nspden,nnat))
944 156 : if(allocated(dtpawuj(icyc)%xred)) then
945 116 : ABI_FREE(dtpawuj(icyc)%xred)
946 : end if
947 468 : ABI_MALLOC(dtpawuj(icyc)%xred,(3,nnat))
948 :
949 :
950 624 : rrtab=(/0.75_dp,0.815_dp,1.0_dp/)*pawtab(typawujat)%rpaw
951 624 : wftab=pawtab(typawujat)%phi(pawtab(typawujat)%mesh_size,pawtab(typawujat)%lnproju(1))
952 :
953 624 : do im1=1,ncoeff
954 468 : if (pawrad(typawujat)%mesh_type==1) then
955 0 : im2=nint(rrtab(im1)/pawrad(typawujat)%rstep+1)
956 468 : else if (pawrad(typawujat)%mesh_type==2) then
957 468 : im2=nint(log(rrtab(im1)/pawrad(typawujat)%rstep+1)/pawrad(typawujat)%lstep+1)
958 0 : else if (pawrad(typawujat)%mesh_type==3) then
959 0 : im2=nint(log(rrtab(im1)/pawrad(typawujat)%rstep)/pawrad(typawujat)%lstep+1)
960 0 : else if (pawrad(typawujat)%mesh_type==4) then
961 0 : im2=nint(pawtab(typawujat)%mesh_size*(1-exp((-one)*rrtab(im1)/pawrad(typawujat)%rstep))+1)
962 : end if
963 :
964 468 : rrtab(im1)=pawrad(typawujat)%rad(im2)
965 624 : wftab(im1)=pawtab(typawujat)%phi(im2,pawtab(typawujat)%lnproju(1))
966 : end do
967 156 : write(message,fmt='(a,i3,a,10f10.5)')' pawuj_red: mesh_type',pawrad(typawujat)%mesh_type,' rrtab:', rrtab
968 156 : call wrtout(std_out,message,'COLL')
969 156 : write(message,fmt='(a,10f10.5)')' pawuj_red: wftab', wftab
970 156 : call wrtout(std_out,message,'COLL')
971 3432 : a=reshape((/ (( rrtab(im2)**(im1-1), im1=1,3) ,im2=1,3 )/),(/ncoeff,ncoeff/))
972 156 : write(messg,fmt='(a)')'A'
973 156 : call linvmat(a,b,ncoeff,messg,2,0.0_dp,3) ! linvmat(inmat,oumat,nat,nam,option,gam,prtvol)
974 156 : write(std_out,*) 'pawuj_red: a,b ', a,b
975 2028 : wftab=matmul(wftab,b)
976 156 : write(std_out,*) 'pawuj_red: wftab ', wftab
977 624 : dtpawuj(icyc)%wfchr(4:6)=wftab
978 :
979 156 : dtpawuj(icyc)%nat=nnat
980 156 : write(std_out,*) 'pawuj_red: m1'
981 1293 : dtpawuj(icyc)%vsh=reshape(pack(atvshift,atvshmusk),(/ nspden,nnat /))
982 : !factor in next line to compensate nocctot contains just occ of 1 spin channel for nspden=1
983 156 : write(std_out,*) 'pawuj_red: m2'
984 1962 : dtpawuj(icyc)%occ=reshape(pack(nnocctot,dmusk),(/nspden,nnat/))*(3-nspden)
985 156 : write(std_out,*) 'pawuj_red: m3'
986 :
987 825 : write(std_out,*) 'pawuj_red: occ ', dtpawuj(icyc)%occ
988 :
989 1516 : dtpawuj(icyc)%xred=reshape(pack(dtset%xred_orig(:,:,1),musk),(/3,nnat/))
990 156 : dtpawuj(icyc)%ph0phiint=pawtab(typawujat)%ph0phiint(1)
991 : dtpawuj(icyc)%wfchr(1:3)=(/ pawtab(typawujat)%zioneff(1)*(dtset%lpawu(typawujat)+2),&
992 624 : & one*(dtset%lpawu(typawujat)+1),one*(dtset%lpawu(typawujat))/)
993 156 : dtpawuj(icyc)%pawrad=pawtab(typawujat)%rpaw
994 :
995 1092 : write(std_out,*) 'pawuj_red: wfchr ',dtpawuj(icyc)%wfchr
996 :
997 246 : if (icyc.le.2) dtpawuj(icyc)%vsh=0.0d0
998 :
999 156 : write (hstr,'(I0)') icyc
1000 156 : write(message,'(a,a,I3,a)') ch10, '---------- MARK ------ ',icyc,ch10
1001 156 : call wrtout(std_out,message,'COLL')
1002 156 : write(message,fmt='(a)') 'vsh'//trim(hstr)
1003 156 : call wrtout(std_out,message,'COLL')
1004 156 : call prmat(dtpawuj(icyc)%vsh(:,:),1,nnat*nspden,1)
1005 156 : write(message,fmt='(a)') 'occ'//trim(hstr)
1006 156 : call wrtout(std_out,message,'COLL')
1007 156 : call prmat(dtpawuj(icyc)%occ(:,:),1,nnat*nspden,1)
1008 156 : write(message, '(3a)' )'---------- MARK ---------- ',ch10
1009 302 : call wrtout(std_out,message,'COLL')
1010 : end do
1011 : end if !usepawu
1012 :
1013 146 : ABI_FREE(nnocctot)
1014 146 : ABI_FREE(dmusk)
1015 146 : ABI_FREE(atvshift)
1016 146 : ABI_FREE(atvshmusk)
1017 :
1018 : !Destroy atom table used for parallelism
1019 146 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1020 :
1021 146 : end subroutine pawuj_red
1022 : !!***
1023 :
1024 : !----------------------------------------------------------------------
1025 :
1026 : !!****f* m_paw_uj/chiscwrt
1027 : !! NAME
1028 : !! chiscwrt
1029 : !!
1030 : !! FUNCTION
1031 : !! PRIVATE ROUTINE
1032 : !! Distributes values of chi_org on chi_sc according to ion-ion distances and multiplicities in shells
1033 : !!
1034 : !! INPUTS
1035 : !! chi_org = response matrix in original unit cell
1036 : !! disv_org = distances (multiplied by magntization) in original unit cell
1037 : !! nat_org = number of atoms in original unit cell
1038 : !! sdisv_org = radii of shells in original unit cell
1039 : !! smult_org = multiplicities of shells in original unit cell
1040 : !! nsh_org = number of shells in original unit cell
1041 : !! disv_sc = distances (multiplied by magntization) in super-cell
1042 : !! nat_sc = number of atoms in super-cell
1043 : !! sdisv_sc = radii of shells in super-cell (was unused, so suppressed)
1044 : !! smult_sc = multiplicities of shells in super-cell
1045 : !! nsh_sc = number of shells in super-cell
1046 : !! opt =
1047 : !!
1048 : !! OUTPUT
1049 : !! chi_sc = first line of reponse matrix in super-cell
1050 : !!
1051 : !! SIDE EFFECTS
1052 : !!
1053 : !! SOURCE
1054 :
1055 90 : subroutine chiscwrt(chi_org,disv_org,nat_org,sdisv_org,smult_org,nsh_org,chi_sc,&
1056 90 : & disv_sc,nat_sc,smult_sc,nsh_sc,opt,prtvol)
1057 :
1058 : !Arguments ------------------------------------
1059 : !scalars
1060 : integer,intent(in) :: nat_org,nsh_org,nsh_sc
1061 : integer,intent(in) :: nat_sc
1062 : integer,intent(in),optional :: prtvol
1063 : integer,intent(in),optional :: opt
1064 : !arrays
1065 : real(dp),intent(in) :: chi_org(nat_org),disv_org(nat_org),sdisv_org(nsh_org)
1066 : integer,intent(in) :: smult_org(nsh_org), smult_sc(nsh_sc)
1067 : real(dp),intent(in) :: disv_sc(nat_sc)
1068 : real(dp),intent(out) :: chi_sc(nat_sc)
1069 :
1070 : !Local variables-------------------------------
1071 : !scalars
1072 : integer :: iatom,jatom,jsh,optt
1073 : character(len=500) :: message
1074 : !arrays
1075 90 : real(dp) :: chi_orgl(nat_org)
1076 :
1077 : ! *************************************************************************
1078 :
1079 : if (present(opt)) then
1080 : optt=opt
1081 : else
1082 : optt=1
1083 : end if
1084 :
1085 :
1086 220 : do iatom=1,nat_org
1087 260 : do jsh=1,nsh_org
1088 170 : if (disv_org(iatom)==sdisv_org(jsh)) then
1089 130 : if (opt==2) then
1090 0 : chi_orgl(iatom)=chi_org(iatom)
1091 130 : else if (opt==1) then
1092 130 : chi_orgl(iatom)=chi_org(iatom)*smult_org(jsh)/smult_sc(jsh)
1093 : end if
1094 : exit
1095 : end if
1096 : end do !iatom
1097 : end do !jsh
1098 :
1099 90 : if (prtvol>1) then
1100 0 : write(message,fmt='(a,150f10.5)')' chiscwrt: chi at input ',chi_org
1101 0 : call wrtout(std_out,message,'COLL')
1102 0 : write(message,fmt='(a,150f10.5)')' chiscwrt: chi after division ',chi_orgl
1103 0 : call wrtout(std_out,message,'COLL')
1104 : end if
1105 :
1106 8970 : do iatom=1,nat_sc
1107 21880 : do jatom=1,nat_org
1108 21790 : if (disv_org(jatom)==disv_sc(iatom)) then
1109 410 : chi_sc(iatom)=chi_orgl(jatom)
1110 410 : exit
1111 12910 : else if (jatom==nat_org) then
1112 8470 : chi_sc(iatom)=0_dp
1113 : end if
1114 : end do
1115 : end do
1116 :
1117 90 : if (prtvol>1) then
1118 0 : write(message,'(a)')' chiscwrt, chi_sc '
1119 0 : call wrtout(std_out,message,'COLL')
1120 0 : call prmat(chi_sc,1,nat_sc,1,std_out)
1121 : end if
1122 :
1123 90 : end subroutine chiscwrt
1124 : !!***
1125 :
1126 : !----------------------------------------------------------------------
1127 :
1128 : !!****f* m_paw_uj/linvmat
1129 : !! NAME
1130 : !! linvmat
1131 : !!
1132 : !! FUNCTION
1133 : !! inverts real matrix inmat
1134 : !!
1135 : !! INPUTS
1136 : !! inmat(1:nat,1:nat)=matrix to be inverted
1137 : !! nat=dimension of inmat
1138 : !! nam=comment specifiying the input matrix (to be printed in output)
1139 : !! option=how to invert inmat
1140 : !! option=1 or 3 add charge bath to matrix and add gam for inversion
1141 : !! option=2 simply invert matrix
1142 : !! gam=gamma added to inmat before inversion in case charge bath is used (allows inversion of otherwise
1143 : !! singular matrix)
1144 : !! prtvol=controls output to files (see subroutine lprtmat)
1145 : !!
1146 : !! OUTPUT
1147 : !! oumat(nnat,nnat)=inverse of inmat, nnat=nat+1 for option=1 or option=3; nnat=nat for option=2
1148 : !!
1149 : !! SOURCE
1150 :
1151 266 : subroutine linvmat(inmat,oumat,nat,nam,option,gam,prtvol)
1152 :
1153 : !Arguments -------------------------------
1154 :
1155 : integer,intent(in) :: nat
1156 : real(dp),intent(in) :: gam,inmat(nat,nat)
1157 : real(dp),intent(inout) :: oumat(:,:) ! nat+1,nat+1 for option=1 or 2
1158 : ! nat,nat for option=2
1159 : character(len=500),intent(in) :: nam
1160 : integer,intent(in),optional :: prtvol,option
1161 :
1162 : !Local variables -------------------------
1163 : character(len=500) :: message
1164 : character(len=500) :: bastrin,gastrin
1165 : integer :: info,nnat,optionn !,ii,jj
1166 266 : integer,allocatable :: ipvt(:)
1167 266 : real(dp),allocatable :: hma(:,:),work(:)
1168 :
1169 : ! *********************************************************************
1170 :
1171 : if (present(option)) then
1172 : optionn=option
1173 : else
1174 : optionn=1
1175 : end if
1176 :
1177 266 : if (option==1.or.option==3) then
1178 110 : write(bastrin,'(a)')'+ charge bath '
1179 110 : write(gastrin,'(a,d10.2,a)')'+ gamma (=',gam,') '
1180 : else
1181 156 : write(bastrin,'(a)')''
1182 156 : write(gastrin,'(a)')''
1183 : end if
1184 :
1185 266 : write(message,fmt='(a)')' matrix '//trim(nam)
1186 266 : call lprtmat(message,1,prtvol,inmat,nat)
1187 :
1188 266 : if (option==1.or.option==3) then
1189 110 : call blow_pawuj(inmat,nat,oumat)
1190 110 : write(message,fmt='(a,a)')' ',trim(nam)//trim(bastrin)
1191 110 : call lprtmat(message,1,prtvol,oumat,nat+1)
1192 1519370 : oumat=oumat+gam
1193 110 : nnat=nat+1
1194 : else
1195 156 : nnat=nat
1196 2028 : oumat=inmat
1197 156 : oumat(1,1)=inmat(1,1)
1198 : end if
1199 :
1200 :
1201 1064 : ABI_MALLOC(hma,(nnat,nnat))
1202 798 : ABI_MALLOC(work,(nnat))
1203 798 : ABI_MALLOC(ipvt,(nnat))
1204 9754 : work=0_dp
1205 1521398 : hma(:,:)=oumat
1206 :
1207 266 : call dgetrf(nnat,nnat,hma,nnat,ipvt,info)
1208 266 : if (.not.info==0) then
1209 0 : write(message, '(3a)' ) 'Matrix '//trim(nam)//' is singular',ch10,'Probably too many symmetries kept'
1210 0 : call wrtout(ab_out,message,'COLL')
1211 0 : return
1212 : end if
1213 :
1214 266 : call dgetri(nnat,hma,nnat,ipvt,work,nnat,info)
1215 1521398 : oumat=hma(:,:)
1216 :
1217 266 : write(message,fmt='(2a,a)')' ('//trim(nam)//trim(bastrin)//trim(gastrin)//')^(-1)'
1218 266 : call lprtmat(message,1,prtvol,oumat,nnat)
1219 :
1220 266 : ABI_FREE(hma)
1221 266 : ABI_FREE(work)
1222 266 : ABI_FREE(ipvt)
1223 :
1224 0 : end subroutine linvmat
1225 : !!***
1226 :
1227 : !----------------------------------------------------------------------
1228 :
1229 : !!****f* m_paw_uj/lprtmat
1230 : !! NAME
1231 : !! lprtmat
1232 : !!
1233 : !! FUNCTION
1234 : !! prints out the real matrix mmat
1235 : !!
1236 : !! INPUTS
1237 : !! mmat(nat,nat)=matrix to be printed
1238 : !! nat=dimension of mmat
1239 : !! prtvol specifies the volume of printing
1240 : !! 3: print the whole matrix
1241 : !! 2: print the first line
1242 : !! 1: do not print anything
1243 : !! chan specifies the output files
1244 : !! 1: output only to std_out
1245 : !! 2: output also to ab_out
1246 : !! commnt=comment specifying matirix
1247 : !!
1248 : !! OUTPUT
1249 : !! oumat(nat+1,nat+1)=inverse of inmat
1250 : !!
1251 : !! SOURCE
1252 :
1253 697 : subroutine lprtmat(commnt,chan,prtvol,mmat,nat)
1254 :
1255 : !Arguments -------------------------------
1256 : integer,intent(in) :: nat,chan,prtvol
1257 : real(dp),intent(in) :: mmat(nat,nat)
1258 : character(len=500),intent(in) :: commnt
1259 :
1260 : !Local variables -------------------------
1261 : character(len=500) :: message
1262 : ! *********************************************************************
1263 :
1264 697 : if (prtvol==3) then
1265 312 : write(message,fmt='(a)') trim(commnt)
1266 312 : call wrtout(std_out,message,'COLL')
1267 312 : call prmat(mmat,nat,nat,nat,std_out)
1268 312 : if (chan==2) then
1269 0 : call wrtout(ab_out,message,'COLL')
1270 0 : call prmat(mmat,nat,nat,nat,ab_out)
1271 : end if
1272 312 : write(message,*)ch10
1273 312 : call wrtout(std_out,message,'COLL')
1274 312 : if (chan==2) then
1275 0 : call wrtout(ab_out,message,'COLL')
1276 : end if
1277 : end if
1278 :
1279 697 : if (prtvol==2) then
1280 0 : write(message,fmt='(a)') trim(commnt)
1281 0 : call wrtout(std_out,message,'COLL')
1282 0 : call prmat(mmat,1,nat,nat,std_out)
1283 0 : if (chan==2) then
1284 0 : call wrtout(ab_out,message,'COLL')
1285 0 : call prmat(mmat,1,nat,nat,ab_out)
1286 : end if
1287 0 : write(message,*)ch10
1288 0 : call wrtout(std_out,message,'COLL')
1289 0 : if (chan==2) then
1290 0 : call wrtout(ab_out,message,'COLL')
1291 : end if
1292 : end if
1293 :
1294 697 : end subroutine lprtmat
1295 : !!***
1296 :
1297 : !----------------------------------------------------------------------
1298 :
1299 : !!****f* m_paw_uj/lcalcu
1300 : !! NAME
1301 : !! lcalcu
1302 : !!
1303 : !! FUNCTION
1304 : !! prints out real the real matrice mmat
1305 : !!
1306 : !! INPUTS
1307 : !! magv=magnetic ordering of the ions (-1 of down/1 for up)
1308 : !! natom=number of atoms
1309 : !! rprimd(3,3)=lattic vectors of unit cell
1310 : !! xred(3,natom)=positions of atoms
1311 : !! chi(natom)=full response of atoms due to shift on atom pawujat
1312 : !! chi0(natom)= response of atoms due to shift on atom pawujat
1313 : !! pawujat=specifies on which atom the potential shift was done
1314 : !! prtvol=controls output to files (see subroutine lprtmat)
1315 : !! gam=gamma to be used for inversion of matrices (see subroutine livmat)
1316 : !! opt=wether to use charge bath (1 or 3) or not (else)
1317 : !!
1318 : !! OUTPUT
1319 : !! ures=resulting U (in eV) on atom pawujat
1320 : !!
1321 : !! SOURCE
1322 :
1323 55 : subroutine lcalcu(magv,natom,rprimd,xred,chi,chi0,pawujat,ures,prtvol,gam,opt)
1324 :
1325 : !Arguments ------------------------------------
1326 : !scalars
1327 : integer,intent(in) :: natom
1328 : integer,intent(in),optional :: opt,pawujat,prtvol
1329 : real(dp),intent(in),optional:: gam
1330 : real(dp),intent(out) :: ures
1331 : !arrays
1332 : integer,intent(in) :: magv(natom)
1333 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom),chi(natom),chi0(natom)
1334 :
1335 : !Local variables-------------------------------
1336 : !scalars
1337 : character(len=500) :: message
1338 : integer :: optt,prtvoll,nnatom
1339 : real(dp) :: gamm
1340 : !arrays
1341 55 : real(dp),allocatable :: chi0matrix1(:,:),chimatrix2(:,:),auginvchi0matrix3(:,:)
1342 55 : real(dp),allocatable :: auginvchimatrix4(:,:),uresmatrix(:,:)
1343 :
1344 : ! *********************************************************************
1345 :
1346 55 : if (present(opt)) then
1347 55 : optt=opt
1348 : else
1349 0 : optt=1
1350 : end if
1351 :
1352 55 : if (present(prtvol)) then
1353 55 : prtvoll=prtvol
1354 : else
1355 0 : prtvoll=1
1356 : end if
1357 :
1358 55 : if (present(gam)) then
1359 55 : gamm=gam
1360 : else
1361 0 : gamm=1_dp
1362 : end if
1363 :
1364 55 : if (optt==1.or.optt==3) then
1365 55 : nnatom=natom+1
1366 : else
1367 0 : nnatom=natom
1368 : end if
1369 :
1370 220 : ABI_MALLOC(chi0matrix1,(natom,natom))
1371 165 : ABI_MALLOC(chimatrix2,(natom,natom))
1372 220 : ABI_MALLOC(auginvchi0matrix3,(nnatom,nnatom))
1373 165 : ABI_MALLOC(auginvchimatrix4,(nnatom,nnatom))
1374 165 : ABI_MALLOC(uresmatrix,(nnatom,nnatom))
1375 :
1376 55 : call ioniondist(natom,rprimd,xred,chi0matrix1,3,chi0,magv,pawujat,prtvoll)
1377 55 : call ioniondist(natom,rprimd,xred,chimatrix2,3,chi,magv,pawujat,prtvoll)
1378 :
1379 55 : write(message,fmt='(a)')'response chi_0'
1380 55 : call linvmat(chi0matrix1,auginvchi0matrix3,natom,message,optt,gamm,prtvoll)
1381 55 : call wrtout(std_out,message,'COLL')
1382 :
1383 55 : write(message,fmt='(a)')'response chi'
1384 55 : call linvmat(chimatrix2,auginvchimatrix4,natom,message,optt,gamm,prtvoll)
1385 55 : call wrtout(std_out,message,'COLL')
1386 :
1387 759740 : uresmatrix=(auginvchi0matrix3-auginvchimatrix4)*Ha_eV
1388 :
1389 55 : write(message,fmt='(a,i3,a)')' (chi_0)^(-1)-(chi)^(-1) (eV)'
1390 55 : call lprtmat(message,2,prtvoll,uresmatrix,nnatom)
1391 55 : call wrtout(std_out,message,'COLL')
1392 :
1393 55 : ures=uresmatrix(1,pawujat)
1394 :
1395 55 : ABI_FREE(chi0matrix1)
1396 55 : ABI_FREE(chimatrix2)
1397 55 : ABI_FREE(auginvchi0matrix3)
1398 55 : ABI_FREE(auginvchimatrix4)
1399 55 : ABI_FREE(uresmatrix)
1400 :
1401 :
1402 55 : end subroutine lcalcu
1403 : !!***
1404 :
1405 : !----------------------------------------------------------------------
1406 :
1407 : !!****f* m_paw_uj/blow_pawuj
1408 : !!
1409 : !! NAME
1410 : !! blow_pawuj
1411 : !!
1412 : !! FUNCTION
1413 : !! This subroutine reads a real nxn matrice and appends lines n+1 and clumn n+1 containing
1414 : !! the sum of the lines
1415 : !!
1416 : !! INPUTS
1417 : !! mat(nj,nj) matrix to be completed
1418 : !!
1419 : !! OUTPUT
1420 : !! matt(nj+1,nj+1) completed matrix
1421 : !!
1422 : !! SOURCE
1423 :
1424 110 : subroutine blow_pawuj(mat,nj,matt)
1425 :
1426 : !Arguments ------------------------------------
1427 : !scalars
1428 : integer,intent(in) :: nj
1429 : !arrays
1430 : real(dp),intent(in) :: mat(nj,nj)
1431 : real(dp),intent(out) :: matt(nj+1,nj+1)
1432 :
1433 : !Local variables-------------------------------
1434 : !scalars
1435 : integer :: ii
1436 : !arrays
1437 :
1438 : ! *************************************************************************
1439 :
1440 1501330 : matt(1:nj,1:nj)=mat
1441 9020 : do ii = 1,nj
1442 1501330 : matt(ii,nj+1)=-sum(matt(ii,1:nj))
1443 : end do
1444 :
1445 9130 : do ii = 1,nj+1
1446 1510350 : matt(nj+1,ii)=-sum(matt(1:nj,ii))
1447 : end do
1448 :
1449 110 : end subroutine blow_pawuj
1450 : !!***
1451 :
1452 : !----------------------------------------------------------------------
1453 :
1454 0 : END MODULE m_paw_uj
1455 : !!***
|
