Line data Source code
1 : !!****m* ABINIT/m_calc_ucrpa
2 : !! NAME
3 : !! m_calc_ucrpa
4 : !!
5 : !! FUNCTION
6 : !! Calculate the effective interaction in the correlated orbitals
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon,ROuterovitch)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! INPUTS
15 : !!
16 : !! OUTPUT
17 : !!
18 : !! SOURCE
19 :
20 : #if defined HAVE_CONFIG_H
21 : #include "config.h"
22 : #endif
23 :
24 :
25 : #include "abi_common.h"
26 :
27 : MODULE m_calc_ucrpa
28 :
29 : #ifndef HAVE_CRPA_OPTIM
30 : #ifdef FC_INTEL
31 : #warning "optimization of m_calc_ucrpa is deactivated on intel fortran"
32 : !DEC$ NOOPTIMIZE
33 : #endif
34 : #endif
35 :
36 : use defs_basis
37 : implicit none
38 :
39 : private
40 :
41 : public :: calc_ucrpa
42 : !!***
43 :
44 : contains
45 : !!****f* ABINIT/calc_ucrpa
46 : !! NAME
47 : !! calc_ucrpa
48 : !!
49 : !! FUNCTION
50 : !! Calculate the effective interaction in the correlated orbitals
51 : !!
52 : !! COPYRIGHT
53 : !! Copyright (C) 1999-2026 ABINIT group (TApplencourt,BA)
54 : !! This file is distributed under the terms of the
55 : !! GNU General Public License, see ~abinit/COPYING
56 : !! or http://www.gnu.org/copyleft/gpl.txt .
57 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
58 : !!
59 : !! INPUTS
60 : !! npwe : number of plane wave for the dielectric constant
61 : !! npw : number of plane wave
62 : !! nomega : number of frequencis
63 : !! bandinf,bandsup : kohn sham band
64 : !! optimisation : string for the optimisation
65 : !! Wfd:: MPI communicator
66 : !! mesh <kmesh_t>
67 : !! %nbz=Number of points in the BZ
68 : !! %nibz=Number of points in IBZ
69 : !! %kibz(,nibz)=k-point coordinates, irreducible Brillouin zone
70 : !! %kbz(3,nbz)=k-point coordinates, full Brillouin zone
71 : !! %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
72 : !! %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
73 : !! %ktabp(nbz)= phase factor associated to tnons
74 : !! M1_q_m(bandinf:bandsup,bandinf:bandsup,npw,Qmesh%nibz): Oscillator strengh in Wannier basis
75 : !! rhot1_q_m(bandinf:bandsup,bandinf:bandsup,npw,Qmesh%nibz): Oscillator strengh
76 : !! multiplied by coulomb potential in Wannier basis
77 : !! OUTPUT
78 : !!
79 : !! NOTES
80 : !!
81 : !! SOURCE
82 :
83 0 : subroutine calc_ucrpa(itypatcor,cryst,Kmesh,lpawu,M1_q_m,Qmesh,npwe,&
84 0 : & npw,nsym,nomega,omegamin,omegamax,bandinf,bandsup,optimisation,ucvol,Wfd,fname,plowan_compute,rhot1,wanbz)
85 :
86 : use defs_basis
87 : use m_abicore
88 : use m_xmpi
89 : use m_errors
90 :
91 : use m_io_tools, only : open_file
92 : use m_wfd, only : wfd_t
93 : use m_io_screening, only : read_screening, em1_ncname
94 : use m_bz_mesh, only : kmesh_t
95 : use m_crystal, only : crystal_t
96 : use m_plowannier, only : operwan_realspace_type,plowannier_type
97 : ! _____ _
98 : ! |_ _| | |
99 : ! | | _ __ _ __ _ _| |_
100 : ! | | | '_ \| '_ \| | | | __|
101 : ! _| |_| | | | |_) | |_| | |_
102 : ! |_____|_| |_| .__/ \__,_|\__|
103 : ! | |
104 : ! |_|
105 :
106 : !Arguments ------------------------------------
107 : integer, intent(in) :: itypatcor,lpawu,npw,npwe,nsym
108 : integer, intent(in) :: nomega
109 : integer, intent(in) :: bandinf
110 : integer, intent(in) :: bandsup
111 : integer, intent(in) :: plowan_compute
112 : character(len=fnlen), intent(in) :: fname
113 : character(len=*), intent(in) :: optimisation
114 : real(dp), intent(in) :: ucvol,omegamin,omegamax
115 :
116 : class(wfd_t),intent(inout) :: Wfd
117 : type(kmesh_t),intent(in) :: Kmesh,Qmesh
118 : type(crystal_t),intent(in) :: Cryst
119 : type(operwan_realspace_type),intent(in) :: rhot1(npw,Qmesh%nibz)
120 : type(plowannier_type),intent(in) :: wanbz
121 : complex(dp), intent(in) :: M1_q_m(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,2*lpawu+1,2*lpawu+1,npw,Qmesh%nibz)
122 :
123 : !Local variables ------------------------------
124 : !scalars
125 : real(dp) :: x
126 : real(dp) :: t1,t2
127 : real(dp):: tol
128 : complex(dp) :: uu,jj
129 :
130 : complex :: nC,ualter
131 :
132 : integer :: iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,ispin,one_orbital
133 : integer :: im_paral,iqalloc,ib1,ib2,m1,m2,m3,m4,iqibz,mbband1,mbband2,mbband3,mbband4,spin1,spin2
134 : integer :: ierr,ik_bz,ik_ibz,iq_ibz,i,iG1,iG2,iG,iiG,iomega,iomega1,ispinor1,ispinor2,ispinor3,ispinor4
135 : integer :: lpawu_read,nkibz,nbband,nkbz,nprocs,nqalloc,nqibz,ms1,ms2,ms3,ms4,mbband,nspinor
136 : integer :: isym_kgw,iik,unt, cp_paral
137 : complex(dp) ::ph_mkt,cplx1,cplx2
138 :
139 : logical :: wannier=.TRUE.
140 : logical :: verbose=.FALSE.
141 : logical :: bug=.FALSE.
142 : logical :: lscr_one
143 :
144 : character(len=500) :: message
145 :
146 : !arrays
147 0 : complex(dp), allocatable :: V_m(:,:,:,:)
148 0 : complex(dp), allocatable :: U_m(:,:,:,:)
149 0 : complex(dp),allocatable :: uspin(:,:),jspin(:,:)
150 : ! complex(dp), allocatable :: coeffW_BZ(:,:,:),coeffW_IBZ(:,:,:)
151 0 : complex(dp), allocatable :: rhot_q_m1m3(:,:,:,:,:,:),rhot_q_m2m4(:,:,:,:,:,:)
152 0 : complex(dp), allocatable :: rhot_q_m1m3_npwe(:,:,:,:,:,:),rhot_q_m2m4_npwe(:,:,:,:,:,:)
153 0 : complex(dp),allocatable :: trrho(:,:),sumrhorhoeps(:)
154 0 : complex(gwp), allocatable :: scr(:,:,:,:)
155 :
156 0 : real(dp),allocatable :: k_coord(:,:)!,k_coordIBZ(:,:)
157 0 : real(dp),allocatable :: q_coord(:,:)
158 0 : real(dp),allocatable:: normG(:)
159 0 : complex(dp),allocatable:: uomega(:),jomega(:)
160 0 : real(dp),allocatable:: omega(:)
161 0 : complex(dp),allocatable:: eiqr(:)
162 :
163 0 : integer,allocatable:: ikmq_bz_t(:,:)
164 :
165 0 : logical,allocatable :: bijection(:)
166 : !************************************************************************
167 :
168 0 : write(message,*) ch10, '==== Calculation of the screened interaction ===='
169 0 : call wrtout(std_out,message,'COLL')
170 0 : call wrtout(ab_out,message,'COLL')
171 0 : write(message,*) ""
172 0 : call wrtout(std_out,message,'COLL')
173 0 : call wrtout(ab_out,message,'COLL')
174 0 : nkbz = Kmesh%nbz
175 0 : nqibz= Qmesh%nibz
176 0 : nspinor=Wfd%nspinor
177 :
178 0 : nbband=1+bandsup-bandinf
179 : ! _ __ ____
180 : ! | |/ / ___ / __ \
181 : ! | ' / ( _ ) | | | |
182 : ! | < / _ \/\ | | | |
183 : ! | . \ | (_> < | |__| |
184 : ! |_|\_\ \___/\/ \___\_\
185 :
186 0 : write(message,*) "Read K and Q mesh"
187 0 : call wrtout(std_out,message,'COLL')
188 0 : call wrtout(ab_out,message,'COLL')
189 :
190 0 : ABI_MALLOC(k_coord,(nkbz,3))
191 0 : ABI_MALLOC(q_coord,(nqibz,4))
192 0 : ABI_MALLOC(eiqr,(nqibz))
193 0 : eiqr=czero
194 : !==Read k and q==!
195 : !open(unit=2012,file='ikbz_COORD',form='formatted',status='unknown')
196 : !read(2012,*) (ik_bz,k_coord(ik_bz,:),i=1,nkbz)
197 : !close(2012)
198 :
199 0 : do ik_bz=1,nkbz
200 0 : call kmesh%get_BZ_item(ik_bz,k_coord(ik_bz,:),ik_ibz,isym_kgw,iik,ph_mkt)
201 : end do
202 :
203 : ! open(unit=2012,file='iqbz_COORD',form='formatted',status='unknown')
204 : ! read(2012,*)
205 0 : do i=1,nqibz
206 : ! read(2012,*) iq_ibz,q_coord(iq_ibz,:)
207 0 : q_coord(i,1)=Qmesh%ibz(1,i)
208 0 : q_coord(i,2)=Qmesh%ibz(2,i)
209 0 : q_coord(i,3)=Qmesh%ibz(3,i)
210 0 : q_coord(i,4)=Qmesh%wt(i)
211 : ! if (iq_ibz > nqibz) then
212 : ! write(message,*) iq_ibz,nqibz," Error on line",i,"Are you in iBZ ?"
213 : ! call wrtout(std_out,message,'COLL')
214 : ! end if
215 : end do
216 : ! close(2012)
217 :
218 : !==Bijection and array for k-q==!
219 0 : ABI_MALLOC(bijection,(nkbz))
220 0 : ABI_MALLOC(ikmq_bz_t,(nkbz,nqibz))
221 0 : bijection(:)=.FALSE.
222 0 : if (nsym==1) then
223 0 : do ik_bz=1,nkbz
224 0 : do iq_ibz=1,nqibz
225 0 : ikmq_bz_t(ik_bz,iq_ibz)=findkmq(ik_bz,k_coord,q_coord(iq_ibz,:),nkbz)
226 0 : if (ikmq_bz_t(ik_bz,iq_ibz)>nkbz.and.nsym==1) then
227 0 : BUG=.TRUE.
228 0 : write(message,*) "No K-Q for K/Q =",ik_bz,iq_ibz
229 0 : ABI_ERROR(message)
230 : end if
231 0 : bijection(ikmq_bz_t(ik_bz,iq_ibz))=.TRUE.
232 : end do
233 :
234 0 : if (count(bijection).NE.nqibz.and.nsym==1) then
235 0 : BUG=.TRUE.
236 0 : write(message,*) 'No bijection ',ik_bz
237 0 : ABI_ERROR(message)
238 : end if
239 :
240 0 : bijection(:)=.FALSE.
241 : end do
242 :
243 0 : if (.NOT.BUG.and.nsym==1) then
244 0 : write(message,*) "Bijection Ok."
245 0 : call wrtout(std_out,message,'COLL')
246 : end if
247 : endif
248 : ! _____
249 : ! / ____|
250 : ! _ __ ___ _ __ _ __ ___ ___ | | __
251 : ! | '_ \ / _ \| '__| '_ ` _ \ / _ \ | | |_ |
252 : ! | | | | (_) | | | | | | | | __/ | |__| |
253 : ! |_| |_|\___/|_| |_| |_| |_|\___| \_____|
254 :
255 0 : ABI_MALLOC(normG,(npw))
256 0 : if (verbose) then
257 0 : write(message,*) 'Read the potential and G norm'
258 0 : call wrtout(std_out,message,'COLL')
259 0 : if (open_file('normeG',message,newunit=unt,form='formatted',status='unknown') /= 0) then
260 0 : ABI_ERROR(message)
261 : end if
262 0 : read(unt,*) (iiG,x,normG(iiG),iG=1,npw)
263 0 : close(unt)
264 : !!False norme for G=0 idd G is the inverse of the potential inverse du potentiel (q=0)
265 0 : normG(1)=0
266 : end if
267 :
268 : !========================================================================
269 : !------------------------------------------------------------------------
270 :
271 : ! FIRST PART OF THE ROUTINE: USE M_G^(nn')(q,k) to do formal checks.
272 : ! USE rhot_q_n to do compute bare interaction
273 :
274 : !------------------------------------------------------------------------
275 : !========================================================================
276 :
277 : ! _
278 : ! ( )
279 : ! _ __ ___ _ __ _ __ |/
280 : ! | '_ ` _ \ | '_ \| '_ \
281 : ! | | | | | | | | | | | | |
282 : ! |_| |_| |_| |_| |_|_| |_|
283 :
284 : !==========================================================
285 : !==========================================================
286 : ! tol=1E-1
287 : ! tolerance for the normalization of wfc: should be around 0.01.
288 0 : tol = 1 ! very large for test.
289 0 : write(message,*) 'Check the norm of M'
290 0 : call wrtout(std_out,message,'COLL')
291 0 : write(message,*) 'Tolerance :',tol
292 0 : call wrtout(std_out,message,'COLL')
293 : ! __ __
294 : ! \ \ / /
295 : ! \ \ / / _ __
296 : ! \ \/ / | '_ \
297 : ! \ / | | | |
298 : ! \/ |_| |_|
299 : !==========================================================================
300 : !==Compute V_{n,n'}: bare interaction in the KS basis <- rhot_q_n -> V_n
301 : !==========================================================================
302 : !==========================================================
303 : !==Compute V_{n,n'}
304 : !==========================================================
305 0 : if(verbose) then
306 0 : write(message,*) ""
307 0 : call wrtout(std_out,message,'COLL')
308 0 : call wrtout(ab_out,message,'COLL')
309 0 : write(message,*) "==Calcul of the bare kohn-sham interaction V n=="
310 0 : call wrtout(std_out,message,'COLL')
311 : endif
312 0 : tol=1E+1
313 :
314 0 : if (.NOT.wannier) RETURN
315 :
316 : !========================================================================
317 : !------------------------------------------------------------------------
318 :
319 : ! SECOND PART OF THE ROUTINE: Read Wannier projections
320 :
321 : !------------------------------------------------------------------------
322 : !========================================================================
323 : !
324 : ! \ \ / / (_)
325 : ! \ \ /\ / /_ _ _ __ _ __ _ ___ _ __
326 : ! \ \/ \/ / _` | '_ \| '_ \| |/ _ \ '__|
327 : ! \ /\ / (_| | | | | | | | | __/ |
328 : ! \/ \/ \__,_|_| |_|_| |_|_|\___|_|
329 :
330 : !==========================================================
331 : !==========================================================
332 : !== Read Wannier coefficient in forlb.ovlp
333 : !==========================================================
334 0 : nkibz=Kmesh%nibz
335 :
336 : !Read "l"
337 0 : if (plowan_compute<10) then
338 0 : write(message,*) ""
339 0 : call wrtout(std_out,message,'COLL')
340 0 : call wrtout(ab_out,message,'COLL')
341 0 : write(message,*) "Read wannier in iBZ"
342 0 : call wrtout(ab_out,message,'COLL')
343 0 : call wrtout(std_out,message,'COLL')
344 0 : if (open_file('forlb.ovlp',message,newunit=unt,form='formatted',status='unknown') /= 0) then
345 0 : ABI_ERROR(message)
346 : end if
347 0 : rewind(unt)
348 0 : read(unt,*) message
349 0 : read(unt,*) message, lpawu_read
350 0 : read(unt,*) message, ib1, ib2
351 0 : close(unt)
352 0 : mbband=2*lpawu_read+1
353 : else
354 0 : ib1=wanbz%bandi_wan
355 0 : ib2=wanbz%bandf_wan
356 0 : mbband=2*wanbz%latom_wan(1)%lcalc(1)+1
357 0 : write(message,*)"Read l and bands from wanbz",ib1,ib2,mbband
358 0 : call wrtout(ab_out,message,'COLL')
359 0 : call wrtout(std_out,message,'COLL')
360 : endif
361 0 : if(ib1/=bandinf.and.ib2/=bandsup) then
362 0 : write(message,*) "Error with bands",ib1,bandinf,ib2,bandsup
363 0 : ABI_ERROR(message)
364 : endif
365 : !!Read the bandinf, bandinf redondance information
366 :
367 :
368 : !*******************************************************
369 : !if (3==4) then
370 : !! USELESS START
371 : !!*******************************************************
372 : ! ABI_MALLOC(coeffW_IBZ,(bandinf:bandsup,nkibz,mbband))
373 : ! coeffW_IBZ=czero
374 : !
375 : !
376 : ! if(Wfd%my_rank==0) then
377 : ! do ik_ibz=1,nkibz
378 : ! !read k
379 : ! read(2012,*)
380 : ! do iband=bandinf,bandsup
381 : ! !read band
382 : ! read(2012,*)
383 : ! !read projection
384 : ! do m1=1,mbband
385 : ! read(2012,*) binR,binR,binR,x,y
386 : ! !write(message,*) binR,binR,binR,x,y
387 : ! coeffW_IBZ(iband,ik_ibz,m1)=cmplx(x,y)
388 : ! end do
389 : ! end do
390 : ! end do
391 : ! endif
392 : ! call xmpi_barrier(Wfd%comm)
393 : ! call xcast_mpi(coeffW_IBZ,0,Wfd%comm,ierr)
394 : ! call xmpi_barrier(Wfd%comm)
395 : ! close(2012)
396 : !
397 : ! ABI_MALLOC(coeffW_BZ,(bandinf:bandsup,nkbz,mbband))
398 : !
399 : ! if (nkbz==nkibz) then
400 : ! coeffW_BZ=coeffW_IBZ
401 : ! else
402 : ! write(message,*) "Reconstruct in full BZ"
403 : ! call wrtout(std_out,message,'COLL')
404 : ! call wrtout(ab_out,message,'COLL')
405 : ! ABI_MALLOC(k_coordIBZ,(nkibz,3))
406 : !
407 : !! k_coordIBZ(:,:)=q_coord(:,1:3)
408 : !
409 : ! bijection(:)=.FALSE.
410 : ! write(message,*) "Indice in iBZ | Indice in BZ | Inverse in BZ"
411 : ! call wrtout(std_out,message,'COLL')
412 : ! do ik_bz=1,Kmesh%nbz
413 : ! write(6,*) "ik",ik_bz,Kmesh%tab(ik_bz),Kmesh%tabi(ik_bz),Kmesh%tabo(ik_bz)
414 : ! if(Kmesh%tabi(ik_bz)==1) then
415 : ! coeffW_BZ(:,ik_bz,:)=coeffW_IBZ(:,Kmesh%tab(ik_bz),:)
416 : ! bijection(ik_bz)=.TRUE.
417 : ! else if(Kmesh%tabi(ik_bz)==-1) then
418 : ! coeffW_BZ(:,ik_bz,:)=conjg(coeffW_IBZ(:,Kmesh%tab(ik_bz),:))
419 : ! write(message,*) Kmesh%tab(ik_bz),ik_bz
420 : ! inverse_ik_bz=Kmesh%tab(ik_bz)
421 : ! bijection(ik_bz)=.TRUE.
422 : !! bijection(inverse_ik_bz)=.TRUE.
423 : ! endif
424 : ! enddo
425 : !
426 : ! if (count(bijection).NE.nkbz) then
427 : ! BUG=.TRUE.
428 : ! write(message,*) 'Miss somme K point for the Wannier',count(bijection),"/",nkbz
429 : ! ABI_ERROR(message)
430 : ! end if
431 : !
432 : ! if (.NOT.BUG) then
433 : ! write(message,*) "Reconstruction Success"
434 : ! call wrtout(std_out,message,'COLL')
435 : ! call wrtout(ab_out,message,'COLL')
436 : ! end if
437 : ! ABI_FREE(k_coordIBZ)
438 : ! end if
439 : !
440 : ! ABI_FREE(coeffW_IBZ)
441 : !
442 : ! wk=1.0/nkbz
443 : !
444 : ! write(message,*) 'Orthogonality check'
445 : ! call wrtout(std_out,message,'COLL')
446 : ! write(message,*) 'Sum on all the k point ,on all the Kohn-Sham band of C_(m1)*C_m(2)'
447 : ! call wrtout(std_out,message,'COLL')
448 : !
449 : !! tolerance for the sum over k-points of the Wannier functions (orthogonality)
450 : ! tol=1E-5
451 : !
452 : !! Sum for one k-point (should be around 0.1).
453 : ! tol2=10E0
454 : !
455 : ! write(message,*) 'Tolerance : k',tol2,'m',tol
456 : ! call wrtout(std_out,message,'COLL')
457 : !
458 : !
459 : ! ! Here checks on Wannier coeff are done (ortho, norm..)
460 : ! nC=cmplx(0,0)
461 : ! BUG=.FALSE.
462 : !
463 : ! do m1=1,mbband
464 : ! do m2=1,mbband
465 : ! do ik_bz=1,nkbz
466 : ! nCt=sum(conjg(coeffW_BZ(:,ik_bz,m1))*coeffW_BZ(:,ik_bz,m2))
467 : ! if ( ((m1==m2).and.(abs(abs(ncT)-1)>tol2)).OR.&
468 : ! ((m1.NE.m2).and.(abs(ncT)>tol2)) ) then
469 : ! BUG=.TRUE.
470 : ! write(message,*) "No orthogonality for m1,m2",m1,m2,"kpt",ik_bz,abs(nCt)
471 : ! ABI_ERROR(message)
472 : ! end if
473 : ! nC=nC+wk*nCt
474 : ! end do
475 : !
476 : ! if ( ((m1==m2).and.(abs(abs(nC)-1)>tol)).OR.&
477 : ! ((m1.NE.m2).and.(abs(nC)>tol)) ) then
478 : ! bug=.TRUE.
479 : ! write(message,*) "No orthogonality for",m1,m2,abs(nC)
480 : ! ABI_ERROR(message)
481 : ! end if
482 : ! write(message,*) m1,m2,abs(nC)
483 : ! nC=cmplx(0,0)
484 : ! end do
485 : ! end do
486 : ! if (.NOT.bug) then
487 : ! write(message,*) "Orthogonality check"
488 : ! call wrtout(std_out,message,'COLL')
489 : ! end if
490 : !
491 : !!*******************************************************
492 : !end if
493 : ! USELESS END
494 : !*******************************************************
495 :
496 : ! do iq_ibz=1,nqibz
497 : ! do m1=1,mbband
498 : ! do m2=1,mbband
499 : ! write(6,*) "M1_q_m",M1_q_m(m1,m2,1,iq_ibz)
500 : ! end do
501 : ! end do
502 : ! end do
503 :
504 : if(real(M1_q_m(1,1,1,1,1,1,1))>0) then
505 : endif
506 :
507 : ! _
508 : ! ( )
509 : ! _ __ ___ _ __ ___ _ __ ___ |/
510 : ! | '_ ` _ \ | '_ ` _ \| '_ ` _ \
511 : ! | | | | | | | | | | | | | | | | |
512 : ! |_| |_| |_| |_| |_| |_|_| |_| |_|
513 :
514 :
515 : !!!==Calculation of M_G^(mm')(q)==!
516 :
517 0 : write(message,*) 'Calculation of M m'
518 0 : call wrtout(std_out,message,'COLL')
519 :
520 : !=================================================================!
521 : !==Compute V_{m,m'}(q,z): Oscillator strengh in the Wannier basis
522 : !=================================================================!
523 :
524 : ! Sum over k-points for the oscillator strengh in the Wannier basis
525 :
526 : if (verbose) then
527 : ! call Sauvegarde_M_q_m(M1_q_m,normG,nqibz,npw,mbband)
528 : end if
529 :
530 : ! write(message,*) "M1_q_m,m for iG=53 and iq=8"
531 : ! call wrtout(std_out,message,'COLL')
532 :
533 : ! _ _
534 : ! | | ( )
535 : ! _ __| |__ ___ _ __ ___ _ __ ___ |/
536 : ! | '__| '_ \ / _ \ | '_ ` _ \| '_ ` _ \
537 : ! | | | | | | (_) | | | | | | | | | | | |
538 : ! |_| |_| |_|\___/ |_| |_| |_|_| |_| |_|
539 :
540 0 : write(message,*) ""
541 0 : call wrtout(std_out,message,'COLL')
542 0 : call wrtout(ab_out,message,'COLL')
543 : !==Calcul de Ro_G^(mm')(q)==!
544 0 : write(message,*) 'Calculation of rhotwilde m'
545 0 : call wrtout(std_out,message,'COLL')
546 :
547 : !===============================================================!
548 : !===============================================================!
549 : !==Compute V_{m,m'}(q,z): bare interaction in the Wannier basis
550 : !===============================================================!
551 : !===============================================================!
552 0 : ABI_FREE(ikmq_bz_t)
553 :
554 : !ABI_MALLOC(rhot_q_m,(nspinor,nspinor,mbband,mbband,npw,nqibz))
555 :
556 : ! if (plowan_compute>=10)then
557 : ! write(message,*)" cRPA calculation using plowan module"
558 : ! call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
559 : ! do iG=1,npw
560 : ! do iqibz=1,nqibz
561 : ! do ispinor1=1,nspinor
562 : ! do ispinor2=1,nspinor
563 : ! do m1=1,2*wanbz%latom_wan(1)%lcalc(1)+1
564 : ! do m2=1,2*wanbz%latom_wan(1)%lcalc(1)+1
565 : ! rhot_q_m(ispinor1,ispinor2,m1,m2,iG,iqibz)=&
566 : ! &rhot1(iG,iqibz)%atom_index(1,1)%position(1,1)%atom(1,1)%matl(m1,m2,1,ispinor1,ispinor2)
567 : ! !write(6,*) rhot_q_m(ispinor1,ispinor2,m1,m2,iG,iqbz)
568 : ! enddo
569 : ! enddo
570 : ! enddo
571 : ! enddo
572 : ! enddo
573 : ! enddo
574 : ! endif
575 0 : do iatom1=1,wanbz%natom_wan
576 0 : do iatom2=1,wanbz%natom_wan
577 0 : do iatom3=1,wanbz%natom_wan
578 0 : do iatom4=1,wanbz%natom_wan
579 0 : if (iatom1/=iatom3 .or. iatom2/=iatom4)cycle
580 0 : do pos1=1,size(wanbz%nposition(iatom1)%pos,1)
581 0 : do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
582 0 : do pos3=1,size(wanbz%nposition(iatom3)%pos,1)
583 0 : do pos4=1,size(wanbz%nposition(iatom4)%pos,1)
584 0 : if (pos1/=pos3 .or. pos2/=pos4)cycle
585 0 : do il1=1,wanbz%nbl_atom_wan(iatom1)
586 0 : do il2=1,wanbz%nbl_atom_wan(iatom2)
587 0 : do il3=1,wanbz%nbl_atom_wan(iatom3)
588 0 : do il4=1,wanbz%nbl_atom_wan(iatom4)
589 0 : if(il1/=il3 .or. il2/=il4)cycle
590 0 : if (wanbz%nsppol/=1)then
591 0 : ABI_MALLOC(uspin,(4,nomega))
592 0 : ABI_MALLOC(jspin,(4,nomega))
593 : endif
594 0 : if (iatom1==iatom2 .and. pos1==pos2 .and. il1 == il2)then
595 0 : one_orbital=1
596 : else
597 0 : one_orbital=0
598 : endif
599 0 : ABI_MALLOC(omega,(nomega))
600 0 : do spin1=1,wanbz%nsppol
601 0 : do spin2=1,wanbz%nsppol
602 0 : cp_paral=0
603 0 : mbband1=2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
604 0 : mbband2=2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
605 0 : mbband3=2*wanbz%latom_wan(iatom3)%lcalc(il3)+1
606 0 : mbband4=2*wanbz%latom_wan(iatom4)%lcalc(il4)+1
607 0 : ABI_MALLOC(rhot_q_m1m3,(npw,nqibz,mbband1,mbband3,nspinor,nspinor))
608 0 : ABI_MALLOC(rhot_q_m2m4,(npw,nqibz,mbband2,mbband4,nspinor,nspinor))
609 0 : rhot_q_m1m3=czero
610 0 : rhot_q_m2m4=czero
611 0 : do ispinor1=1,wanbz%nspinor
612 0 : do ispinor2=1,wanbz%nspinor
613 0 : do ispinor3=1,wanbz%nspinor
614 0 : do ispinor4=1,wanbz%nspinor
615 0 : do m1=1,mbband1
616 0 : do m2=1,mbband2
617 0 : do m3=1,mbband3
618 0 : do m4=1,mbband4
619 0 : do iqibz=1,nqibz
620 : !Loig Vaugier PhD eq. 5.15
621 : eiqr(iqibz)=exp(cmplx(0.0, 1.0 ) * two_pi * ( &
622 : q_coord(iqibz,1)* ( cryst%xred(1,wanbz%iatom_wan(iatom3)) - cryst%xred(1,wanbz%iatom_wan(iatom2)) )+ &
623 : q_coord(iqibz,2)* ( cryst%xred(2,wanbz%iatom_wan(iatom3)) - cryst%xred(2,wanbz%iatom_wan(iatom2)) )+ &
624 0 : q_coord(iqibz,3)* ( cryst%xred(3,wanbz%iatom_wan(iatom3)) - cryst%xred(3,wanbz%iatom_wan(iatom2)) )))
625 0 : do iG=1,npw
626 : ! cp_paral=cp_paral+1
627 : ! if(mod(cp_paral-1,nprocs)==Wfd%my_rank) then
628 0 : cplx1=rhot1(iG,iqibz)%atom_index(iatom1,iatom3)%position(pos1,pos3)%atom(il1,il3)%matl(m1,m3,spin1,ispinor1,ispinor3)
629 0 : rhot_q_m1m3(iG,iqibz,m1,m3,ispinor1,ispinor3)=cplx1
630 :
631 0 : cplx2=rhot1(iG,iqibz)%atom_index(iatom2,iatom4)%position(pos2,pos4)%atom(il2,il4)%matl(m2,m4,spin2,ispinor2,ispinor4)
632 0 : rhot_q_m2m4(iG,iqibz,m2,m4,ispinor2,ispinor4)=cplx2
633 :
634 : !endif
635 : enddo!iG
636 : enddo!iqibz
637 : enddo!m4
638 : enddo!m3
639 : enddo!m2
640 : enddo!m1
641 : enddo!ispinor4
642 : enddo!ispinor3
643 : enddo!ispinor2
644 : enddo!ispinor1
645 : ! call xmpi_barrier(Wfd%comm) ! First make sure that all processors are here
646 : ! call xmpi_sum(rhot_q_m1m3,Wfd%comm,ierr)
647 : ! call xmpi_sum(rhot_q_m2m4,Wfd%comm,ierr)
648 : ! call xmpi_barrier(Wfd%comm) ! First make sure that all processors are here
649 :
650 : ! __ __
651 : ! \ \ / /
652 : ! \ \ / / _ __ ___
653 : ! \ \/ / | '_ ` _ \
654 : ! \ / | | | | | |
655 : ! \/ |_| |_| |_|
656 :
657 :
658 0 : ABI_MALLOC(V_m,(mbband1*nspinor,mbband2*nspinor,mbband3*nspinor,mbband4*nspinor))
659 0 : write(message,*) "";call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
660 0 : write(message,*) "==Calculation of the bare interaction V m=="
661 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
662 :
663 0 : call print_orbitals(spin1,spin2,iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,wanbz,1)
664 0 : V_m=czero
665 0 : call cpu_time ( t1 )
666 0 : im_paral=0
667 0 : nprocs = xmpi_comm_size(Wfd%comm)
668 :
669 :
670 :
671 :
672 :
673 :
674 0 : do ispinor1=1,nspinor
675 0 : do ispinor2=1,nspinor
676 0 : do ispinor3=1,nspinor
677 0 : do ispinor4=1,nspinor
678 0 : do m1=1,mbband1
679 0 : do m2=1,mbband2
680 0 : do m3=1,mbband3
681 0 : do m4=1,mbband4
682 0 : ms1=m1+(ispinor1-1)*mbband1
683 0 : ms2=m2+(ispinor2-1)*mbband2
684 0 : ms3=m3+(ispinor3-1)*mbband3
685 0 : ms4=m4+(ispinor4-1)*mbband4
686 0 : im_paral=im_paral+1
687 0 : if(mod(im_paral-1,nprocs)==Wfd%my_rank) then
688 : !!somme interne sur iG, puis somme externe sur iq_ibz
689 : !! Sum_(iq_ibz) wi(iq_ibz)*Sum_ig Rho(m3,m1,iG,iq)cong*Rho(m2,m4,ig,iq)
690 : V_m(ms1,ms2,ms3,ms4)=sum(eiqr(:)*q_coord(:,4)*sum(conjg(rhot_q_m1m3(:,:,m3,m1,ispinor3,ispinor1))* &
691 0 : &rhot_q_m2m4(:,:,m2,m4,ispinor2,ispinor4),dim = 1))*Ha_eV/(ucvol)
692 : !end do
693 : endif!paral
694 : end do!m4
695 : end do!m3
696 : end do!m2
697 : end do!m1
698 : end do!ispinor4
699 : end do!ispinor3
700 : end do!ispinor2
701 : end do!ispinor1
702 0 : call xmpi_barrier(Wfd%comm) ! First make sure that all processors are here
703 0 : call xmpi_sum(V_m,Wfd%comm,ierr)
704 0 : call xmpi_barrier(Wfd%comm) ! First make sure that all processors are here
705 0 : call cpu_time ( t2 )
706 0 : write(message,*) "in ",t2-t1,"sec"
707 0 : call wrtout(std_out,message,'COLL')
708 : ! !==Check if calculation is correct
709 0 : tol=1E-2
710 :
711 0 : write(message,*) "BARE INTERACTION"
712 0 : call wrtout(std_out,message,'COLL')
713 0 : call checkk(V_m,1,mbband*nspinor,1,0,uu,jj,"bare interaction",mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital)
714 : !call print_U(mbband1,mbband2,mbband3,mbband4,nspinor,V_m)
715 : ! ========================================================================
716 : ! ------------------------------------------------------------------------
717 :
718 : ! THIRD PART OF THE ROUTINE: Read dielectric function
719 :
720 : ! ------------------------------------------------------------------------
721 : ! ========================================================================
722 :
723 : ! _____ _
724 : ! / ____| (_)
725 : ! | (___ ___ _ __ ___ ___ _ __ _ _ __ __ _
726 : ! \___ \ / __| '__/ _ \/ _ \ '_ \| | '_ \ / _` |
727 : ! ____) | (__| | | __/ __/ | | | | | | | (_| |
728 : ! |_____/ \___|_| \___|\___|_| |_|_|_| |_|\__, |
729 : ! __/ |
730 : ! |___/
731 : ! ==========================================================
732 : ! == Read Dielectric Matrix for _SCR file
733 : ! ==========================================================
734 0 : write(message,*) "";call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
735 0 : call wrtout(std_out,message,'COLL')
736 0 : call wrtout(ab_out,message,'COLL')
737 : ! write(message,*) "==Read the dielectric matrix=="
738 : ! call wrtout(ab_out,message,'COLL'); call wrtout(std_out,message,'COLL')
739 :
740 0 : lscr_one=.true.
741 :
742 : ! if symetry is not enabled then a large number of q-point are used
743 : ! and they are read with direct access to avoid having too large
744 : ! memory.
745 : !-------------------------------------------------------------------
746 : if(nsym>1.and..not.lscr_one) nqalloc=nqibz
747 : if(nsym==1.or.lscr_one) nqalloc=1
748 :
749 0 : ABI_MALLOC(scr,(npwe,npwe,nomega,nqalloc))
750 0 : scr=czero
751 : if(nsym>1.and..not.lscr_one) then
752 : write(message,*) "==Read the dielectric matrix=="
753 : call wrtout(ab_out,message,'COLL'); call wrtout(std_out,message,'COLL')
754 : call read_screening(em1_ncname,fname,npwe,nqibz,nomega,scr,IO_MODE_MPI,Wfd%comm)
755 : endif
756 :
757 : ! if (verbose) then
758 : ! open(unit=2211,file='Screening',form='formatted',status='unknown')
759 : ! do iG1=1,npwe
760 : ! do iG2=1,npwe
761 : ! write(2211,*) iG1,iG2,normG(iG1),normG(iG2),abs(scr(iG1,iG2,1,1)),abs(scr(iG1,iG2,2,1))
762 : ! end do
763 : ! end do
764 : ! close(2211)
765 : ! end if
766 : !
767 : if(nsym>1.and..not.lscr_one) then
768 : write(message,*) "Check the hermiticity"
769 : call wrtout(ab_out,message,'COLL')
770 : call wrtout(std_out,message,'COLL')
771 : tol = 1E-2
772 : do iG1=1,npwe
773 : if (modulo(iG1,100).EQ.1) then
774 : write(message,*) iG1,"/",npw
775 : call wrtout(std_out,message,'COLL')
776 : end if
777 : do iG2=iG1,npwe
778 : if (ANY(abs(scr(iG1,iG2,:,:)-scr(iG2,iG1,:,:))>tol)) then
779 : write(message,*) iG1,iG2,"False"
780 : ABI_ERROR(message)
781 : do iomega1=1,nomega
782 : if(abs(scr(iG1,iG2,iomega1,1)-scr(iG2,iG1,iomega1,1))>tol) then
783 : write(message,*) iG1,iG2,"False",scr(iG1,iG2,iomega1,1),scr(iG2,iG1,iomega1,1)
784 : call wrtout(std_out,message,'COLL')
785 : endif
786 : enddo
787 : end if
788 : if(iG1==iG2) then
789 : scr(iG1,iG1,:,:)=scr(iG1,iG1,:,:)-one ! unscreened part of
790 : endif
791 : end do
792 : end do
793 : write(message,*) "Done: Hermiticity of dielectric matrix checked"
794 : call wrtout(std_out,message,'COLL')
795 : call wrtout(ab_out,message,'COLL')
796 : endif
797 :
798 : ! ========================================================================
799 : ! ------------------------------------------------------------------------
800 :
801 : ! FOURTH PART OF THE ROUTINE: Use dielectric function and oscillator
802 : ! strengh in Wannier basis to compute Screened cRPA interactions.
803 :
804 : ! ------------------------------------------------------------------------
805 : ! ========================================================================
806 : ! _ _
807 : ! | | | |
808 : ! | | | | _ __ ___
809 : ! | | | | | '_ ` _ \
810 : ! | |__| | | | | | | |
811 : ! \____/ |_| |_| |_|
812 : !
813 0 : write(message,*) "";call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
814 0 : write(message,*) "== Calculation of the screened interaction on the correlated orbital U m =="
815 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
816 0 : write(message,*)ch10, " = Start loop over frequency "
817 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
818 :
819 0 : ABI_MALLOC(U_m,(mbband1*nspinor,mbband2*nspinor,mbband3*nspinor,mbband4*nspinor))
820 0 : ABI_MALLOC(rhot_q_m1m3_npwe,(npwe,nqibz,mbband1,mbband3,nspinor,nspinor))
821 0 : ABI_MALLOC(rhot_q_m2m4_npwe,(npwe,nqibz,mbband2,mbband4,nspinor,nspinor))
822 0 : ABI_MALLOC(trrho,(npwe,nqibz))
823 0 : ABI_MALLOC(sumrhorhoeps,(nqibz))
824 :
825 : ! Following lines are only for debug
826 : !--------------------------------------------
827 0 : trrho(:,:)=czero
828 0 : do iG1=1,npwe
829 0 : rhot_q_m1m3_npwe(iG1,:,:,:,:,:)=rhot_q_m1m3(iG1,:,:,:,:,:)
830 0 : rhot_q_m2m4_npwe(iG1,:,:,:,:,:)=rhot_q_m2m4(iG1,:,:,:,:,:)
831 : enddo
832 : ! do ispinor1=1,nspinor
833 : ! do m1=1,mbband
834 : ! trrho(iG1,:)= trrho(iG1,:)+rhot_q_m(ispinor1,ispinor1,m1,m1,iG1,:)
835 : ! enddo
836 : ! enddo
837 : ! enddo
838 : ! write(6,*) "trrho"
839 : ! do iG1=1,npwe
840 : ! do iq_ibz=1,nqibz
841 : ! write(6,*) iG1,iq_ibz,trrho(iG1,iq_ibz)
842 : ! enddo
843 : ! enddo
844 : ! write(6,*)
845 :
846 : ! Loop over frequencies to compute cRPA U(w)
847 : !--------------------------------------------
848 0 : ABI_MALLOC(uomega,(nomega))
849 0 : ABI_MALLOC(jomega,(nomega))
850 :
851 0 : iomega=1
852 0 : do iomega=1,nomega
853 0 : write(message,'(2a,i4,2a)')ch10, " --- For frequency w =",iomega, " -------------",ch10
854 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
855 0 : sumrhorhoeps(:)=czero
856 0 : ualter=czero
857 :
858 0 : write(std_out,*) Optimisation
859 0 : U_m=cmplx(0,0)
860 0 : SELECT CASE(trim(Optimisation))
861 : CASE("naif")
862 : CASE("onlyG")
863 : CASE("Gsum2")
864 : CASE("Gsum")
865 : ! write(message,*) "Optimisation on G and sum"
866 : ! call wrtout(std_out,message,'COLL')
867 0 : U_m=cmplx(0,0)
868 0 : nc=cmplx(0,0)
869 0 : im_paral=0
870 0 : call cpu_time ( t1 )
871 0 : nprocs = xmpi_comm_size(Wfd%comm)
872 0 : do iq_ibz=1,nqibz
873 : ! eiqr=exp(cmplx(0.0, 1.0 ) * two_pi * ( &
874 : ! q_coord(iq_ibz,1)* ( cryst%xred(1,wanbz%iatom_wan(iatom3)) - cryst%xred(1,wanbz%iatom_wan(iatom2)) )+ &
875 : ! q_coord(iq_ibz,2)* ( cryst%xred(2,wanbz%iatom_wan(iatom3)) - cryst%xred(2,wanbz%iatom_wan(iatom2)) )+ &
876 : ! q_coord(iq_ibz,3)* ( cryst%xred(3,wanbz%iatom_wan(iatom3)) - cryst%xred(3,wanbz%iatom_wan(iatom2)) )))
877 : if(nsym>1.and..not.lscr_one) iqalloc=iq_ibz
878 : if(nsym==1.or.lscr_one) iqalloc=1
879 : if(nsym==1.or.lscr_one) then
880 : #ifdef HAVE_MPI_IO
881 0 : call read_screening(em1_ncname,fname,npwe,1,nomega,scr,IO_MODE_MPI,Wfd%comm,iqiA=iq_ibz)
882 : #else
883 : call read_screening(em1_ncname,fname,npwe,1,nomega,scr,IO_MODE_FORTRAN,Wfd%comm,iqiA=iq_ibz)
884 : #endif
885 0 : write(message,*) "Check the hermiticity"
886 0 : call wrtout(std_out,message,'COLL')
887 0 : tol = 0.01_dp
888 0 : do iG1=1,npwe
889 0 : if (modulo(iG1,100).EQ.1) then
890 0 : write(message,*) iG1,"/",npw
891 0 : call wrtout(std_out,message,'COLL')
892 : end if
893 0 : do iG2=iG1,npwe
894 0 : if (ANY(abs(scr(iG1,iG2,:,:)-scr(iG2,iG1,:,:))>tol)) then
895 0 : do iomega1=1,nomega
896 0 : if(abs(scr(iG1,iG2,iomega1,1)-scr(iG2,iG1,iomega1,1))>tol) then
897 0 : write(message,*) iG1,iG2,"False",scr(iG1,iG2,iomega1,1),scr(iG2,iG1,iomega1,1)
898 0 : call wrtout(std_out,message,'COLL')
899 : endif
900 : enddo
901 0 : write(message,*) "CHECK THE HERMITICITY"
902 0 : ABI_WARNING(message)
903 : end if
904 0 : if(iG1==iG2) then
905 0 : scr(iG1,iG1,:,:)=scr(iG1,iG1,:,:)-one ! unscreened part of
906 : ! interaction is computed before
907 : ! scr(iG1,iG1,:,:)=one
908 : ! else
909 : ! scr(iG1,iG2,:,:)=zero
910 : endif
911 : end do
912 : end do
913 : endif
914 0 : write(message,*) "Done: Hermiticity of dielectric matrix checked"
915 :
916 : ! do iG1=1,npwe
917 : ! do iG2=1,npwe
918 : ! sumrhorhoeps(iq_ibz)=sumrhorhoeps(iq_ibz)+conjg(trrho(iG1,iq_ibz))*trrho(iG2,iq_ibz)*scr(iG2,iG1,iomega,1)
919 : ! sumrhorhoeps(iq_ibz)=sumrhorhoeps(iq_ibz) + scr(iG2,iG1,iomega,iqalloc)
920 : ! enddo
921 : ! enddo
922 : ! ualter=ualter+sumrhorhoeps(iq_ibz)*q_coord(iq_ibz,4)*Ha_eV/ucvol/(mbband)**2
923 : ! write(6,*) "sumrhorhoeps",iq_ibz,sumrhorhoeps(iq_ibz)
924 : ! write(6,*) "ualter",iq_ibz,ualter
925 0 : do ispinor1=1,nspinor
926 0 : do ispinor2=1,nspinor
927 0 : do ispinor3=1,nspinor
928 0 : do ispinor4=1,nspinor
929 0 : do m1=1,mbband1
930 0 : do m2=1,mbband2
931 0 : do m3=1,mbband3
932 0 : do m4=1,mbband4
933 0 : ms1=m1+(ispinor1-1)*mbband1
934 0 : ms2=m2+(ispinor2-1)*mbband2
935 0 : ms3=m3+(ispinor3-1)*mbband3
936 0 : ms4=m4+(ispinor4-1)*mbband4
937 0 : im_paral=im_paral+1
938 0 : if(mod(im_paral-1,nprocs)==Wfd%my_rank) then
939 : ! sum q sum G1 sum G2 f(G1,q)f(G2,q)G(G1,G2,q)
940 : ! sum q sum G1 f(g1,q) sum G2 f(G2,q)G(G1,G2,q)
941 0 : do iG1=1,npwe
942 : nc=nc+conjg(rhot_q_m1m3_npwe(iG1,iq_ibz,m1,m3,ispinor1,ispinor3))*&
943 0 : &sum(rhot_q_m2m4_npwe(:,iq_ibz,m2,m4,ispinor2,ispinor4)*scr(:,iG1,iomega,iqalloc))
944 : end do
945 :
946 0 : U_m(ms1,ms2,ms3,ms4)=U_m(ms1,ms2,ms3,ms4)+nc*q_coord(iq_ibz,4)*eiqr(iq_ibz)
947 : ! if(m1==1.and.m2==1.and.m3==1.and.m4==1) then
948 : ! write(6,*) "TEST2"
949 : ! write(6,*) iq_ibz,nc
950 : ! write(6,*) q_coord(iq_ibz,4)
951 : ! endif
952 0 : nc=cmplx(0,0)
953 : endif!paral
954 : end do!m4
955 : end do!m3
956 : end do!m2
957 : end do!m1
958 : end do!ispinor4
959 : end do!ispinor3
960 : end do!ispinor2
961 : end do!ispinor1
962 : end do!iq_ibz
963 0 : U_m(:,:,:,:)=U_m(:,:,:,:)*Ha_eV/(ucvol)
964 0 : call xmpi_barrier(Wfd%comm) ! First make sure that all processors are here
965 0 : call xmpi_sum(U_m,Wfd%comm,ierr)
966 0 : call xmpi_barrier(Wfd%comm) ! First make sure that all processors are here
967 0 : call cpu_time ( t2 )
968 0 : write(message,*) "in ",t2-t1, "sec"
969 0 : call wrtout(std_out,message,'COLL')
970 : END SELECT
971 : ! tolerance of the symetry of screened U.
972 0 : tol=1E-2
973 0 : call checkk(U_m,1,mbband*nspinor,0,iomega,uu,jj,"UminusVbare",mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital)
974 0 : U_m=V_m+U_m
975 0 : write(message,*) "UCRPA interaction"
976 0 : call wrtout(std_out,message,'COLL')
977 : call checkk(U_m,1,mbband*nspinor,1,iomega,uomega(iomega),jomega(iomega),&
978 0 : &"cRPA interaction",mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital)
979 0 : if (spin1==1 .and. spin2==1) then
980 : ispin=1
981 0 : else if(spin1==1 .and. spin2==2) then
982 : ispin=2
983 0 : else if (spin1==2 .and. spin2==1)then
984 : ispin=3
985 : else
986 0 : ispin=4
987 : endif
988 0 : if (wanbz%nsppol/=1)then
989 0 : uomega(iomega)=4*uomega(iomega)
990 0 : jomega(iomega)=4*jomega(iomega)
991 0 : uspin(ispin,iomega)=uomega(iomega)
992 0 : jspin(ispin,iomega)=jomega(iomega)
993 : endif
994 : !write(message,*)ch10,"SCREENED INTERACTION"
995 : !call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
996 : !call print_U(mbband1,mbband2,mbband3,mbband4,nspinor,U_m)
997 : enddo
998 :
999 0 : call print_orbitals(spin1,spin2,iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,wanbz,0)
1000 0 : write(message,*)ch10," -------------------------------------------------------------"
1001 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1002 0 : write(message,*)" Average U and J as a function of frequency "
1003 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1004 0 : write(message,*)" -------------------------------------------------------------"
1005 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1006 0 : write(message,*)" omega U(omega) J(omega)"
1007 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1008 0 : do iomega=1,nomega
1009 0 : if(nomega==1) then
1010 0 : omega(iomega)=omegamin
1011 : else
1012 0 : omega(iomega)=(omegamax-omegamin)/(nomega-1)*(iomega-1)+omegamin
1013 : endif
1014 0 : if(nomega==1) omega(iomega)=omegamin
1015 0 : write(message,'(2x,f11.3,2x,2f10.4,2x,2f10.4)') omega(iomega)*Ha_eV, uomega(iomega),jomega(iomega)
1016 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1017 : enddo
1018 0 : write(message,*)" -------------------------------------------------------------"
1019 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1020 :
1021 : ! END OF LOOP ON ATOMS !!
1022 :
1023 0 : ABI_FREE(uomega)
1024 0 : ABI_FREE(jomega)
1025 0 : ABI_FREE(V_m)
1026 0 : ABI_FREE(rhot_q_m1m3_npwe)
1027 0 : ABI_FREE(rhot_q_m2m4_npwe)
1028 0 : ABI_FREE(trrho)
1029 0 : ABI_FREE(sumrhorhoeps)
1030 0 : ABI_FREE(scr)
1031 0 : ABI_FREE(U_m)
1032 0 : ABI_FREE(rhot_q_m1m3)
1033 0 : ABI_FREE(rhot_q_m2m4)
1034 : ! Print dielectric matrix
1035 : ! do iq_ibz=1,nqibz
1036 : ! call read_screening(em1_ncname,fname,npwe,1,nomega,scr,IO_MODE_FORTRAN,Wfd%comm,iqiA=iq_ibz)
1037 : ! enddo
1038 : enddo!spin2
1039 : enddo!spin1
1040 0 : if (wanbz%nsppol/=1)then
1041 0 : do iomega=1,nomega
1042 0 : if(nomega==1) then
1043 0 : omega(iomega)=omegamin
1044 : else
1045 0 : omega(iomega)=(omegamax-omegamin)/(nomega-1)*(iomega-1)+omegamin
1046 : endif
1047 0 : if(nomega==1) omega(iomega)=omegamin
1048 : enddo
1049 : !call print_orbitals(1,1,iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,wanbz,2)
1050 0 : call print_uj_spin(nomega,uspin,jspin,omega,one_orbital)
1051 : endif
1052 0 : if (wanbz%nsppol/=1)then
1053 0 : ABI_FREE(uspin)
1054 0 : ABI_FREE(jspin)
1055 : endif
1056 0 : ABI_FREE(omega)
1057 : enddo!il4
1058 : enddo!il3
1059 : enddo!il2
1060 : enddo!il1
1061 : enddo!pos4
1062 : enddo!pos3
1063 : enddo!pos2
1064 : enddo!pos1
1065 : enddo!iatom4
1066 : enddo!iatom3
1067 : enddo!iatom2
1068 : enddo!iatom1
1069 0 : ABI_FREE(k_coord)
1070 0 : ABI_FREE(q_coord)
1071 0 : ABI_FREE(eiqr)
1072 0 : ABI_FREE(bijection)
1073 0 : ABI_FREE(normG)
1074 : ! ABI_FREE(coeffW_BZ)
1075 : !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1076 : !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1077 :
1078 : CONTAINS
1079 :
1080 : !! __ _ _
1081 : !! / _| | | (_)
1082 : !! | |_ ___ _ __ ___| |_ _ ___ _ __
1083 : !! | _/ _ \| '_ \ / __| __| |/ _ \| '_ \
1084 : !! | || (_) | | | | (__| |_| | (_) | | | |
1085 : !! |_| \___/|_| |_|\___|\__|_|\___/|_| |_|
1086 : !
1087 0 : integer FUNCTION fi(nkbz,k_coord,kprime_coord)
1088 :
1089 : integer,intent(in) :: nkbz
1090 : real(dp),dimension(nkbz,3),intent(in) ::k_coord
1091 : real(dp),dimension(3),intent(in) :: kprime_coord(3)
1092 :
1093 0 : do fi=1,nkbz
1094 0 : if (ALL(abs(kprime_coord(:)-k_coord(fi,:))<0.001)) then
1095 : exit
1096 : end if
1097 : end do
1098 0 : END FUNCTION fi
1099 :
1100 0 : integer FUNCTION findkmq(ik_bz,k_coord,q_coord,nkbz)
1101 :
1102 : integer,intent(in) :: ik_bz,nkbz
1103 : real(dp),dimension(nkbz,3),intent(in) ::k_coord
1104 : real(dp),dimension(4),intent(in) ::q_coord
1105 : real(dp),dimension(3) :: kprime_coord
1106 : integer :: i,j,k
1107 :
1108 :
1109 0 : kprime_coord(:)=k_coord(ik_bz,:)-q_coord(1:3)
1110 :
1111 0 : where (kprime_coord > 0.5)
1112 : kprime_coord(:)= kprime_coord(:)-1
1113 : elsewhere(kprime_coord < -0.5)
1114 : kprime_coord(:)= kprime_coord(:)+1
1115 : end where
1116 :
1117 : !! indice of k -q
1118 0 : findkmq=fi(nkbz,k_coord,kprime_coord)
1119 :
1120 : !!Test if k-q exists
1121 0 : if (findkmq.EQ.(nkbz+1)) then
1122 : !!The prb comes from PBC included born_inf,born_sup]
1123 : !! One test all combination over boundaries
1124 0 : do i=1,2
1125 0 : if (abs(abs(kprime_coord(1))-0.5)<0.001) kprime_coord(1)=(-1)**i*0.5
1126 0 : do j=1,2
1127 0 : if (abs(abs(kprime_coord(2))-0.5)<0.001) kprime_coord(2)=(-1)**j*0.5
1128 0 : do k=1,2
1129 0 : if (abs(abs(kprime_coord(3))-0.5)<0.01) kprime_coord(3)=(-1)**k*0.5
1130 0 : findkmq=fi(nkbz,k_coord,kprime_coord)
1131 : !!Quand on a trouver la bonne valeur on part
1132 0 : if (findkmq.NE.(nkbz+1)) return
1133 : end do
1134 : end do
1135 : end do
1136 : end if
1137 : END FUNCTION findkmq
1138 :
1139 0 : SUBROUTINE checkk(Interaction,m_inf,m_sup,prtopt,ifreq,uu,jj,utype,mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital)
1140 :
1141 : integer, intent(in) :: m_inf,m_sup,ifreq,mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital
1142 : complex(dp), intent(in) :: Interaction(mbband1*nspinor,mbband2*nspinor,mbband3*nspinor,mbband4*nspinor)
1143 : complex(dp), intent(out) :: uu,jj
1144 : character(len=*), intent(in) :: utype
1145 : integer :: prtopt
1146 :
1147 :
1148 : !==Check correctness
1149 : ! write(message,*) "== Check == "
1150 : ! call wrtout(std_out,message,'COLL')
1151 : ! write(message,*) 'Tolerance :',tol
1152 : ! call wrtout(std_out,message,'COLL')
1153 : ! do i=m_inf,m_sup
1154 : ! do j=i+1,m_sup
1155 : !if (abs(abs(Interaction(i,i,i,i)-Interaction(j,j,j,j)))>tol) then
1156 : ! BUG=.TRUE.
1157 : ! write(message,*) "Problem in the interband calculation"&
1158 : !& ,i,j,abs(Interaction(i,i,i,i)),abs(Interaction(j,j,j,j))
1159 : ! call wrtout(std_out,message,'COLL')
1160 : !end if
1161 :
1162 : ! if (abs(Interaction(i,j,i,j)-Interaction(j,i,j,i))>tol) then
1163 : ! BUG=.TRUE.
1164 : ! write(message,*) "Warning in the symetry of U'",i,j,abs(Interaction(i,j,i,j)),&
1165 : ! & abs(Interaction(j,i,j,i)),abs(Interaction(i,j,i,j)-Interaction(j,i,j,i))
1166 : ! call wrtout(std_out,message,'COLL')
1167 : ! end if
1168 :
1169 : ! if (abs(Interaction(i,i,j,j)-Interaction(j,j,i,i))>tol) then
1170 : ! BUG=.TRUE.
1171 : ! write(message,*) "Warning in the symetry of J'",i,j,abs(Interaction(i,i,j,j)),&
1172 : ! & abs(Interaction(j,j,i,i)),abs(Interaction(j,j,i,i)-Interaction(i,i,j,j))
1173 : ! call wrtout(std_out,message,'COLL')
1174 : ! end if
1175 : ! end do
1176 : ! end do
1177 :
1178 :
1179 : ! do i=m_inf,m_sup
1180 : ! do j=m_inf,m_sup
1181 : ! if (i.EQ.j) cycle
1182 : ! if (abs(Interaction(i,j,j,j))>tol) then
1183 : ! BUG=.TRUE.
1184 : ! write(message,*) "Warning in the symetry U(i,j,j,j) should vanish (in the Ylm basis) for",&
1185 : !& i,j,abs(Interaction(i,j,j,j))
1186 : ! call wrtout(std_out,message,'COLL')
1187 : ! end if
1188 : ! end do
1189 : ! end do
1190 :
1191 : ! if (.not.BUG) then
1192 : ! call wrtout(std_out,'Calcul is possibly correct','COLL')
1193 : ! call Affichage(Interaction,m_inf,m_sup,2)
1194 : ! else
1195 : ! call wrtout(std_out,'Maybe somme error','COLL')
1196 : ! call Affichage(Interaction,m_inf,m_sup,1)
1197 : ! end if
1198 :
1199 0 : if(prtopt>0) call Affichage(Interaction,m_inf,m_sup,1,ifreq,uu,jj,utype,mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital)
1200 :
1201 : END SUBROUTINE checkk
1202 :
1203 0 : SUBROUTINE Affichage(Interaction,m_inf,m_sup,option,ifreq,uu,jj,utype,mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital)
1204 :
1205 : integer, intent(in) :: m_inf,m_sup,option,ifreq,mbband1,mbband2,mbband3,mbband4,nspinor,one_orbital
1206 : complex(dp), intent(in) :: Interaction(mbband1*nspinor,mbband2*nspinor,mbband3*nspinor,mbband4*nspinor)
1207 : complex(dp),intent(out) :: UU,JJ
1208 : character(len=*), intent(in) :: utype
1209 : complex(dp) :: UU1,UUmJJ,JJ1,JJ2
1210 : integer :: m1,m2
1211 : logical :: lprint
1212 : character(len=500) :: message
1213 :
1214 :
1215 0 : if(utype=="cRPA interaction".or.utype=="bare interaction") then
1216 : lprint=.true.
1217 : else
1218 0 : lprint=.false.
1219 : endif
1220 0 : write(message,*) "";call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1221 :
1222 0 : if (one_orbital==1) then
1223 0 : if(lprint) then
1224 0 : write(message,*)" Diagonal ",utype
1225 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1226 : endif
1227 :
1228 0 : do m1=m_inf,m_sup
1229 0 : if (option.EQ.1) then
1230 0 : write(message,'(a,i3,14f7.3)') " ",m1,real(Interaction(m1,m1,m1,m1))
1231 0 : call wrtout(std_out,message,'COLL')
1232 0 : call wrtout(ab_out,message,'COLL')
1233 : end if
1234 : end do
1235 : endif
1236 :
1237 0 : if(lprint) then
1238 0 : write(message,*) "";call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1239 0 : write(message,*)" U'=U(m1,m2,m1,m2) for the ",utype
1240 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1241 :
1242 0 : write(message,'(a,14i7)') " -",(m2,m2=1,mbband2)
1243 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1244 0 : do m1=m_inf,m_sup
1245 0 : if (option.EQ.1) then
1246 0 : write(message,'(a,i3,14f7.3)') " ",m1,(real(Interaction(m1,m2,m1,m2)),m2=1,mbband2)
1247 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1248 : end if
1249 : end do
1250 : endif
1251 :
1252 0 : write(message,*) ""
1253 0 : call wrtout(std_out,message,'COLL')
1254 0 : call wrtout(ab_out,message,'COLL')
1255 :
1256 0 : UU=czero
1257 0 : do m1=1,mbband1
1258 0 : do m2=1,mbband2
1259 0 : UU=UU+Interaction(m1,m2,m1,m2)
1260 : enddo
1261 : enddo
1262 0 : UU=UU/float(mbband1*mbband2)
1263 0 : if(ifreq/=0) write(message,'(3a,i3,a,2f10.4,a)')' Hubbard ',utype,' for w =',ifreq,', U=1/(2l+1)**2 \sum U(m1,m2,m1,m2)=',UU,ch10
1264 0 : if(ifreq==0) write(message,'(3a,2f10.4,a)') ' Hubbard ',utype, ' U=1/(2l+1)**2 \sum U(m1,m2,m1,m2)=',UU,ch10
1265 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1266 :
1267 :
1268 0 : if (one_orbital==1)then
1269 0 : UU1=czero
1270 0 : do m1=1,mbband1
1271 0 : UU1=UU1+Interaction(m1,m1,m1,m1)
1272 : enddo
1273 0 : UU1=UU1/((mbband1))
1274 0 : if(ifreq/=0) write(message,'(3a,i4,a,2f10.4,2a)')&
1275 0 : & ' (Hubbard ',utype,' for w =',ifreq,', U=1/(2l+1) \sum U(m1,m1,m1,m1)=',UU1,')',ch10
1276 0 : if(ifreq==0) write(message,'(3a,2f10.4,2a)')' (Hubbard ',utype,' U=1/(2l+1) \sum U(m1,m1,m1,m1)=',UU1,')',ch10
1277 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1278 : endif
1279 :
1280 0 : if(lprint .and. one_orbital==1) then
1281 0 : write(message,*)' Hund coupling J=U(m1,m1,m2,m2) for the ', utype
1282 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1283 :
1284 0 : write(message,'(a,14i7)') " -",(m2,m2=1,mbband1)
1285 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1286 0 : do m1=1,mbband1
1287 0 : if (option.EQ.1) then
1288 0 : write(message,'(a,i3,14f7.3)') " ",m1,(real(Interaction(m1,m1,m2,m2)),m2=1,mbband1)
1289 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1290 : end if
1291 : end do
1292 : endif
1293 0 : if (one_orbital==1)then
1294 0 : UUmJJ=czero
1295 0 : do m1=1,mbband1
1296 0 : do m2=1,mbband1
1297 0 : UUmJJ=UUmJJ+Interaction(m1,m2,m1,m2)-Interaction(m1,m2,m2,m1)
1298 : enddo
1299 : enddo
1300 0 : if (mbband1/=1) then
1301 0 : UUmJJ=UUmJJ/float((mbband1)*(mbband1-1))
1302 0 : JJ1=UU-UUmJJ
1303 : endif
1304 :
1305 :
1306 0 : JJ=czero
1307 0 : do m1=1,mbband1
1308 0 : do m2=1,mbband1
1309 0 : if(m1/=m2) JJ=JJ+Interaction(m1,m2,m2,m1)
1310 : enddo
1311 : enddo
1312 0 : if (mbband1/=1) then
1313 0 : JJ=JJ/float((mbband1)*(mbband1-1))
1314 : endif
1315 :
1316 0 : write(message,'(a,3x,2a,2f10.4,a)')ch10,utype,&
1317 0 : & ' value of J=U-1/((2l+1)(2l)) \sum_{m1,m2} (U(m1,m2,m1,m2)-U(m1,m2,m2,m1))=',JJ1,ch10
1318 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1319 :
1320 0 : if (mbband1/=1) then
1321 0 : UUmJJ=czero
1322 0 : do m1=1,mbband1
1323 0 : do m2=1,mbband1
1324 0 : UUmJJ=UUmJJ+Interaction(m1,m2,m1,m2)-Interaction(m1,m1,m2,m2)
1325 : enddo
1326 : enddo
1327 0 : UUmJJ=UUmJJ/float((mbband1)*(mbband1-1))
1328 0 : JJ2=UU-UUmJJ
1329 0 : if(abs(JJ1-JJ2)<0.0001) then
1330 0 : JJ=JJ1
1331 : else
1332 :
1333 0 : write(message,'(a,3x,2a,2f10.4,a)')ch10,utype,&
1334 0 : & ' value of J=U-1/((2l+1)(2l)) \sum_{m1,m2} (U(m1,m2,m1,m2)-U(m1,m1,m2,m2))=',JJ2,ch10
1335 0 : call wrtout(std_out,message,'COLL')
1336 : ! stop
1337 : endif
1338 : endif
1339 : endif
1340 :
1341 0 : if(lprint .and. one_orbital==1) then
1342 0 : write(message,*)ch10,' Hund coupling J2=U(m1,m2,m2,m1) for the ', utype
1343 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1344 :
1345 0 : write(message,'(a,14i7)') " -",(m2,m2=1,mbband1)
1346 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1347 0 : do m1=1,mbband1
1348 0 : if (option.EQ.1) then
1349 0 : write(message,'(a,i3,14f7.3)') " ",m1,(real(Interaction(m1,m2,m2,m1)),m2=1,mbband1)
1350 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1351 : end if
1352 : end do
1353 : ! write(message,*) "";call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1354 : ! write(message,*) "U'=U-2J for the t2g should be checked"
1355 : ! call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1356 : endif
1357 :
1358 0 : END SUBROUTINE Affichage
1359 :
1360 : SUBROUTINE Sauvegarde_M_q_m(M_q_m,normG,nqibz,npw,mbband)
1361 :
1362 : integer, intent(in) :: nqibz,npw,mbband
1363 : complex(dp), intent(in) :: M_q_m(mbband,mbband,npw,nqibz)
1364 : real(dp), intent(in) :: normG(npw)
1365 : integer :: i,j,iq_ibz,iG,unt
1366 : character(len=500) :: msg
1367 :
1368 : !==Ecriture de M_(G=0)^(mm')(q) ==!
1369 : if (open_file('M_mimj(n=1_2_3)(q,G=0)',msg,newunit=unt,form='formatted',status='unknown') /=0) then
1370 : ABI_ERROR(msg)
1371 : end if
1372 : do iq_ibz=1,nqibz
1373 : write(unt,*) iq_ibz,((abs(M_q_m(i,j,1,iq_ibz)),i=1,mbband),j=1,mbband)
1374 : end do
1375 : close(unt)
1376 :
1377 : !==Ecriture de M_G^(mm')(q=0) ==!
1378 : if (open_file('M_mm(m=1..mbband)(q=0)',msg,newunit=unt,form='formatted',status='unknown') /= 0) then
1379 : ABI_ERROR(msg)
1380 : end if
1381 : do iG=1,npw
1382 : write(unt,*) normG(iG),(abs(M_q_m(i,i,iG,1)),i=1,mbband)
1383 : end do
1384 : close(unt)
1385 : END SUBROUTINE Sauvegarde_M_q_m
1386 :
1387 0 : subroutine print_orbitals(spin1,spin2,iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,wanbz,opt)
1388 : integer, intent(in) :: spin1,spin2,iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,opt
1389 : type(plowannier_type),intent(in) :: wanbz
1390 : character(len=5000) ::message
1391 : character(len=10):: print_spin
1392 :
1393 0 : if (spin1==spin2 .and. spin1==1)then
1394 0 : print_spin=" Up-Up"
1395 0 : else if(spin1==1 .and. spin2==2)then
1396 0 : print_spin=" Up-Down"
1397 0 : else if(spin1==2 .and. spin2==1) then
1398 0 : print_spin=" Down-Up"
1399 : else
1400 0 : print_spin=" Down-Down"
1401 : endif
1402 0 : if (opt==1)then
1403 0 : write(message,*)ch10,"==Definition of the orbitals=="
1404 0 : else if (opt==0) then
1405 0 : write(message,*)ch10,"==Reminder of the orbitals=="
1406 0 : else if (opt==2) then
1407 0 : write(message,*)ch10,"==Reminder of the orbitals=="
1408 0 : write(print_spin,*)" summary"
1409 0 : else if (opt==3) then
1410 0 : write(message,*)ch10,"==Reminder of the orbitals=="
1411 0 : write(print_spin,*)" average"
1412 : endif
1413 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1414 0 : if (iatom1==iatom2.and. iatom3==iatom4 .and. iatom3==iatom2)then
1415 0 : if(pos1==pos2 .and. pos3==pos4 .and. pos1==pos3)then
1416 0 : if (il1==il2 .and. il3==il4 .and. il3==il1)then
1417 0 : write(message,*)" Only one orbital"
1418 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1419 0 : write(message,*)" Orbital with l=",wanbz%latom_wan(iatom1)%lcalc(il1),&
1420 0 : &"on atom",wanbz%iatom_wan(iatom1),"with spin's orientations",print_spin
1421 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1422 : else
1423 0 : write(message,*)"Different orbitals on atom",wanbz%iatom_wan(iatom1),"with spin's orientations",print_spin
1424 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1425 0 : write(message,*)" orbital 1 with l=",wanbz%latom_wan(iatom1)%lcalc(il1)
1426 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1427 0 : write(message,*)" orbital 2 with l=",wanbz%latom_wan(iatom1)%lcalc(il2)
1428 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1429 0 : write(message,*)" orbital 3 with l=",wanbz%latom_wan(iatom1)%lcalc(il3)
1430 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1431 0 : write(message,*)" orbital 4 with l=",wanbz%latom_wan(iatom1)%lcalc(il4)
1432 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1433 : endif
1434 : else
1435 0 : if(il1==il2 .and. il3==il4 .and. il1==il3)then
1436 0 : write(message,*)" Different position of the same orbital"
1437 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1438 0 : write(message,*)" Orbitals with l=",wanbz%latom_wan(iatom1)%lcalc(il1),&
1439 0 : &"on atom",wanbz%iatom_wan(iatom1),"with spin's orientations",print_spin
1440 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1441 0 : write(message,*)" orbital 1 at postion ",wanbz%nposition(iatom1)%pos(pos1,:)
1442 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1443 0 : write(message,*)" orbital 2 at postion ",wanbz%nposition(iatom1)%pos(pos2,:)
1444 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1445 0 : write(message,*)" orbital 3 at postion ",wanbz%nposition(iatom1)%pos(pos3,:)
1446 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1447 0 : write(message,*)" orbital 4 at postion ",wanbz%nposition(iatom1)%pos(pos4,:)
1448 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1449 : else
1450 0 : write(message,*)" Different orbitals on the same atom, with different positions"
1451 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1452 0 : write(message,*)" Orbitals on atom",wanbz%iatom_wan(iatom1),"with spin's orientations",print_spin
1453 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1454 0 : write(message,*)" orbital 1 with l=",wanbz%latom_wan(iatom1)%lcalc(il1),"at position ",wanbz%nposition(iatom1)%pos(pos1,:)
1455 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1456 0 : write(message,*)" orbital 2 with l=",wanbz%latom_wan(iatom1)%lcalc(il2),"at position ",wanbz%nposition(iatom1)%pos(pos2,:)
1457 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1458 0 : write(message,*)" orbital 3 with l=",wanbz%latom_wan(iatom1)%lcalc(il3),"at position ",wanbz%nposition(iatom1)%pos(pos3,:)
1459 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1460 0 : write(message,*)" orbital 4 with l=",wanbz%latom_wan(iatom1)%lcalc(il4),"at position ",wanbz%nposition(iatom1)%pos(pos4,:)
1461 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1462 : endif
1463 : endif
1464 : else
1465 0 : if(pos1==pos2 .and. pos3==pos4 .and. pos1==pos3)then
1466 0 : write(message,*)"Different atoms, with spin orientation",print_spin
1467 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1468 0 : write(message,*)" orbital 1 with l=",wanbz%latom_wan(iatom1)%lcalc(il1),"on atom",wanbz%iatom_wan(iatom1)
1469 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1470 0 : write(message,*)" orbital 2 with l=",wanbz%latom_wan(iatom2)%lcalc(il2),"on atom",wanbz%iatom_wan(iatom2)
1471 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1472 0 : write(message,*)" orbital 3 with l=",wanbz%latom_wan(iatom3)%lcalc(il3),"on atom",wanbz%iatom_wan(iatom3)
1473 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1474 0 : write(message,*)" orbital 4 with l=",wanbz%latom_wan(iatom4)%lcalc(il4),"on atom",wanbz%iatom_wan(iatom4)
1475 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1476 : else
1477 0 : write(message,*)"Different atoms, in different postions with spin's orientations",print_spin
1478 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1479 0 : write(message,*)" orbital 1 with l=",wanbz%latom_wan(iatom1)%lcalc(il1),&
1480 0 : &"on atom",wanbz%iatom_wan(iatom1),"at position ",wanbz%nposition(iatom1)%pos(pos1,:)
1481 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1482 0 : write(message,*)" orbital 2 with l=",wanbz%latom_wan(iatom2)%lcalc(il2),&
1483 0 : &"on atom",wanbz%iatom_wan(iatom2),"at position ",wanbz%nposition(iatom2)%pos(pos2,:)
1484 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1485 0 : write(message,*)" orbital 3 with l=",wanbz%latom_wan(iatom3)%lcalc(il3),&
1486 0 : &"on atom",wanbz%iatom_wan(iatom3),"at position ",wanbz%nposition(iatom3)%pos(pos3,:)
1487 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1488 0 : write(message,*)" orbital 4 with l=",wanbz%latom_wan(iatom4)%lcalc(il4),&
1489 0 : &"on atom",wanbz%iatom_wan(iatom4),"at position ",wanbz%nposition(iatom4)%pos(pos4,:)
1490 0 : call wrtout(std_out,message,'COLL');call wrtout(ab_out,message,'COLL')
1491 : endif
1492 : endif
1493 0 : end subroutine print_orbitals
1494 :
1495 :
1496 0 : subroutine print_uj_spin(nomega,uspin,jspin,omega,one_orbital)
1497 : integer,intent(in) :: nomega,one_orbital
1498 : complex(dp),intent(in) :: uspin(4,nomega)
1499 : complex(dp),intent(in) :: jspin(4,nomega)
1500 : real(dp),intent(in) :: omega(nomega)
1501 : integer :: iomega,ispin
1502 0 : complex(dp) :: uomega(nomega),jomega(nomega)
1503 : character(len=500)::message
1504 :
1505 0 : call print_orbitals(spin1,spin2,iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,wanbz,3)
1506 0 : write(message,*)ch10," --------------------------------------------------------------------"
1507 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1508 0 : write(message,*)" Sum up of the calcul for different spin polarization and frequencies"
1509 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1510 0 : write(message,*)" --------------------------------------------------------------------"
1511 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1512 :
1513 :
1514 0 : write(message,*)ch10,"Sum up of U",ch10
1515 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1516 :
1517 0 : write(message,'(6a)')" -omega (eV)- "," Up-Up "," Up-Down "," Down-Up "," Down-Down "," Average "
1518 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1519 0 : do iomega=1,nomega
1520 0 : write(message,'(a,f7.2,5f9.3)')" ",omega(iomega)*Ha_eV,(real(uspin(ispin,iomega)),ispin=1,4),sum(real(uspin(:,iomega)))/4
1521 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1522 : end do
1523 :
1524 0 : if (one_orbital==1)then
1525 0 : write(message,*)ch10,"Sum up of J",ch10
1526 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1527 :
1528 0 : write(message,'(6a)')" -omega (eV)- "," Up-Up "," Up-Down "," Down-Up "," Down-Down "," Average "
1529 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1530 0 : do iomega=1,nomega
1531 0 : write(message,'(a,f7.2,5f9.3)')" ",omega(iomega)*Ha_eV,(real(jspin(ispin,iomega)),ispin=1,4),sum(real(jspin(:,iomega)))/4
1532 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1533 : end do
1534 : endif
1535 0 : uomega(:)=sum(uspin,dim=1)/4
1536 0 : jomega(:)=sum(jspin,dim=1)/4
1537 : !call print_orbitals(spin1,spin2,iatom1,iatom2,iatom3,iatom4,pos1,pos2,pos3,pos4,il1,il2,il3,il4,wanbz,3)
1538 0 : write(message,*)ch10," -------------------------------------------------------------"
1539 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1540 0 : write(message,*)" Average U and J as a function of frequency "
1541 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1542 0 : write(message,*)" -------------------------------------------------------------"
1543 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1544 0 : write(message,*)" omega U(omega) J(omega)"
1545 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1546 0 : do iomega=1,nomega
1547 0 : write(message,'(2x,f11.3,2x,2f10.4,2x,2f10.4)') omega(iomega)*Ha_eV, uomega(iomega),jomega(iomega)
1548 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1549 : enddo
1550 0 : write(message,*)" -------------------------------------------------------------"
1551 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1552 0 : end subroutine print_uj_spin
1553 :
1554 : end subroutine calc_ucrpa
1555 : !!***
1556 :
1557 : END MODULE m_calc_ucrpa
1558 : !!***
|
