Branch data 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 : : !!***
|
