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