TABLE OF CONTENTS
ABINIT/qmc_prep_ctqmc [ Functions ]
NAME
qmc_prep_ctqmc
FUNCTION
Prepare and call the qmc subroutines
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (BAmadon,VPlanes) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cryst_struc <type(crystal_t)>=crystal structure data hu <type(hu_type)>= U interaction paw_dmft <type(paw_dmft_type)>= DMFT data structure pawang <type(pawang)>=paw angular mesh and related data pawprtvol = drive the amount of writed data. weiss <type(green_type)>= weiss function
OUTPUT
green <type(green_type)>= green function
NOTES
PARENTS
impurity_solve
CHILDREN
add_matlu,checkreal_matlu,compute_levels,copy_green,copy_matlu ctqmc_triqs_run,ctqmcinterface_finalize,ctqmcinterface_init ctqmcinterface_run,ctqmcinterface_setopts,data4entropydmft_setdocc destroy_green,destroy_matlu,destroy_oper,diag_matlu,diff_matlu fac_matlu,flush_unit,fourier_green,hybridization_asymptotic_coefficient identity_matlu,init_green,init_matlu,init_oper,int_fct,inverse_oper jbessel,occup_green_tau,print_green,print_matlu,printocc_green printplot_matlu,prod_matlu,rotate_matlu,rotatevee_hu,sbf8,shift_matlu slm2ylm_matlu,sym_matlu,testcode_ctqmc,vee_ndim2tndim_hu_r,wrtout,xginv xmpi_barrier,xmpi_bcast
SOURCE
46 #if defined HAVE_CONFIG_H 47 #include "config.h" 48 #endif 49 50 #include "abi_common.h" 51 52 subroutine qmc_prep_ctqmc(cryst_struc,green,self,hu,paw_dmft,pawang,pawprtvol,weiss) 53 54 use defs_basis 55 use defs_datatypes 56 use defs_abitypes 57 use m_errors 58 use m_xmpi 59 60 use m_special_funcs, only : sbf8 61 use m_pawang, only : pawang_type 62 use m_crystal, only : crystal_t 63 use m_green, only : green_type,occup_green_tau,print_green,printocc_green,spline_fct,copy_green,init_green,destroy_green,& 64 & int_fct,greenldacompute_green,fourier_green 65 use m_paw_dmft, only : paw_dmft_type 66 use m_abilasi, only : xginv 67 use m_oper, only : oper_type,destroy_oper,init_oper,inverse_oper 68 use m_self, only : self_type 69 use m_matlu, only : matlu_type,sym_matlu, print_matlu, & 70 & diag_matlu,init_matlu,destroy_matlu,rotate_matlu,checkdiag_matlu,checkreal_matlu, & 71 & copy_matlu, diff_matlu, slm2ylm_matlu, shift_matlu, prod_matlu,fac_matlu,& 72 & add_matlu,printplot_matlu,identity_matlu 73 use m_hu, only : hu_type,rotatevee_hu,vee_ndim2tndim_hu_r 74 use m_Ctqmc 75 use m_CtqmcInterface 76 use m_GreenHyb 77 use m_data4entropyDMFT 78 !use m_self, only : self_type,initialize_self,destroy_self,print_self,rw_self 79 use m_io_tools, only : flush_unit, open_file 80 use m_paw_numeric, only : jbessel=>paw_jbessel 81 82 #if defined HAVE_TRIQS 83 use TRIQS_CTQMC !Triqs module 84 #endif 85 use ISO_C_BINDING 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'qmc_prep_ctqmc' 91 use interfaces_14_hidewrite 92 use interfaces_68_dmft, except_this_one => qmc_prep_ctqmc 93 !End of the abilint section 94 95 implicit none 96 97 !Arguments ------------------------------------ 98 !scalars 99 ! type(pawang_type), intent(in) :: pawang 100 type(crystal_t),intent(in) :: cryst_struc 101 type(green_type), intent(inout) :: green ! MGNAG: This fix the problem with v7[27:29] on nag@petrus 102 type(hu_type), intent(in) :: hu(cryst_struc%ntypat) 103 type(paw_dmft_type), intent(inout) :: paw_dmft 104 type(pawang_type), intent(in) :: pawang 105 integer, intent(in) :: pawprtvol 106 type(green_type), intent(inout) :: weiss 107 type(self_type), intent(in) :: self 108 109 !Local variables ------------------------------ 110 character(len=500) :: message 111 character(len=2) :: gtau_iter,iatomnb 112 integer :: iatom,ierr,if1,if2,iflavor,iflavor1,iflavor2,iflavor3,ifreq,im,im1,ispinor,ispinor1,isppol,itau,itypat,im2,ispinor2 113 integer :: lpawu,master,mbandc,natom,nflavor,nkpt,nspinor,nsppol,nsppol_imp,tndim,ispa,ispb,ima,imb 114 integer :: nproc,opt_diag,opt_nondiag,testcode,testrot,dmft_nwlo,opt_fk,useylm,nomega,opt_rot 115 integer :: ier,rot_type_vee 116 complex(dpc) :: omega_current,integral(2,2) 117 real(dp) :: omega 118 real(dp) :: facd,facnd 119 logical :: nondiaglevels 120 ! arrays 121 real(dp), allocatable :: docc(:,:) 122 real(dp), allocatable, target :: gtmp(:,:), levels_ctqmc(:) !modif 123 complex(dpc), allocatable :: levels_ctqmc_nd(:,:) 124 complex(dpc), allocatable :: hybri_limit(:,:) 125 real(dp), allocatable, target :: gtmp_nd(:,:,:) 126 real(dp) :: umod(2,2) 127 character(len=30) :: tmpfil 128 complex(dpc), allocatable :: fw1(:,:),gw_tmp(:,:) 129 complex(dpc), allocatable, target :: gw_tmp_nd(:,:,:) !modif 130 complex(dpc), allocatable, target :: fw1_nd(:,:,:) !modif 131 complex(dpc), allocatable :: gw1_nd(:,:,:) 132 complex(dpc), allocatable :: shift(:) 133 integer,parameter :: optdb=0 134 type(coeff2_type), allocatable :: udens_atoms(:) 135 ! Type ----------------------------------------- 136 type(coeff2c_type), allocatable :: eigvectmatlu(:,:) 137 type(green_type) :: weiss_for_rot 138 type(matlu_type), allocatable :: dmat_diag(:) 139 type(matlu_type), allocatable :: matlu1(:) 140 type(matlu_type), allocatable :: matlu2(:) 141 type(matlu_type), allocatable :: matlu3(:) 142 type(matlu_type), allocatable :: matlu4(:) 143 type(matlu_type), allocatable :: identity(:) 144 type(matlu_type), allocatable :: level_diag(:) 145 type(oper_type) :: energy_level 146 !type(self_type) :: self 147 ! type(green_type) :: gw_loc 148 type(CtqmcInterface) :: hybrid !!! WARNING THIS IS A BACKUP PLAN 149 type(green_type) :: greenlda 150 type(matlu_type), allocatable :: hybri_coeff(:) 151 ! Var added to the code for TRIQS_CTQMC test and default value ----------------------------------------------------------- 152 logical(kind=1) :: rot_inv = .false. 153 logical(kind=1) :: leg_measure = .true. 154 #if defined HAVE_TRIQS 155 logical(kind=1) :: hist = .false. 156 logical(kind=1) :: wrt_files = .true. 157 logical(kind=1) :: tot_not = .true. 158 #endif 159 160 integer :: nfreq,unt,unt2 161 integer :: ntau ! >= 2*nfreq + 1 162 integer :: nleg 163 integer :: ileg 164 integer :: verbosity_solver ! min 0 -> max 3 165 166 real(dp) :: beta,besp,bespp,xx 167 complex(dpc) :: u_nl 168 169 complex(dpc), allocatable, target ::fw1_nd_tmp(:,:,:) 170 complex(dpc), allocatable, target :: g_iw(:,:,:) 171 real(dp), allocatable, target :: u_mat_ij(:,:) 172 real(dp), allocatable, target :: u_mat_ijkl(:,:,:,:) 173 real(dp), allocatable, target :: u_mat_ijkl_tmp(:,:,:,:) 174 ! real(dp), allocatable, target :: gtau(:,:,:) 175 real(dp), allocatable, target :: gl_nd(:,:,:) 176 real(dp), allocatable :: jbes(:) 177 178 type(c_ptr) :: levels_ptr, fw1_nd_ptr, u_mat_ij_ptr, u_mat_ijkl_ptr, g_iw_ptr, gtau_ptr, gl_ptr 179 180 ! ************************************************************************ 181 mbandc=paw_dmft%mbandc 182 nkpt=paw_dmft%nkpt 183 nsppol=paw_dmft%nsppol 184 natom=paw_dmft%natom 185 nspinor=paw_dmft%nspinor 186 greenlda%whichgreen="LDA" 187 188 call init_green(weiss_for_rot,paw_dmft,opt_oper_ksloc=2) 189 ! weiss_for_rot=>weiss 190 ! call init_green(gw_loc,paw_dmft) 191 call copy_green(weiss,weiss_for_rot,opt_tw=2) 192 !======================================================================= 193 !== Use one QMC solver =============================================== 194 !======================================================================= 195 write(message,'(2a)') ch10,' === CT-QMC solver === ' 196 call wrtout(std_out,message,'COLL') 197 198 ! Initialise for compiler 199 omega_current=czero 200 201 ! Initialise nproc 202 nproc=paw_dmft%nproc 203 204 ! ====================================== 205 ! Allocations: diagonalization and eigenvectors 206 ! ====================================== 207 ABI_DATATYPE_ALLOCATE(udens_atoms,(natom)) 208 ABI_DATATYPE_ALLOCATE(eigvectmatlu,(natom,nsppol)) 209 ABI_DATATYPE_ALLOCATE(dmat_diag,(natom)) 210 ABI_DATATYPE_ALLOCATE(identity,(natom)) 211 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,dmat_diag) 212 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,identity) 213 call identity_matlu(identity,natom) 214 do iatom=1,cryst_struc%natom 215 lpawu=paw_dmft%lpawu(iatom) 216 if(lpawu/=-1) then 217 tndim=nspinor*(2*lpawu+1) 218 do isppol=1,nsppol 219 ABI_ALLOCATE(eigvectmatlu(iatom,isppol)%value,(tndim,tndim)) 220 end do 221 ABI_ALLOCATE(udens_atoms(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1))) 222 dmat_diag(iatom)%mat=czero 223 end if 224 end do 225 226 ! ___________________________________________________________________________________ 227 ! 228 ! FIRST PART: DIAGONALISATION AND ROTATIONS. 229 ! ___________________________________________________________________________________ 230 ! 231 232 ! ================================================================= 233 ! Choose to diagonalize and how to do it 234 ! ================================================================= 235 236 ! ================================================================= 237 ! Impose diago of density matrix 238 ! ================================================================= 239 240 ! ================================================================= 241 ! Impose diago of levels and Ylm basis if opt_nondiag=1 242 ! ================================================================= 243 ! opt_diag=1 ! 1: diago the levels (The best choice). 244 ! opt_diag=2 ! 2: diago density matrix (can be used for historical reasons) 245 246 ! Need in the general case of two input variable for opt_diag and 247 ! opt_nondiag! 248 ! The default value of opt_diag should be 2 for historical reasons (or 249 ! we decide to change the automatic tests) 250 ! opt_nondiag should be 0 by default 251 opt_diag = 1 252 if(paw_dmft%dmft_solv>=6) then 253 opt_nondiag = 1 ! Use cthyb in triqs 254 else 255 opt_nondiag = 0 ! use fast ctqmc in ABINIT without non diagonal terms. 256 end if 257 258 useylm=0 259 if(nspinor==2) then 260 useylm=1 ! to avoid complex G(tau) 261 end if 262 263 !write(6,*) "nspinor,useylm",nspinor,useylm 264 if(useylm==0) then 265 write(std_out,*) " Slm basis is used (before rotation)" 266 rot_type_vee=1 ! for rotatevee_hu 267 else if(useylm==1) then 268 write(std_out,*) " Ylm basis is used (before rotation)" 269 rot_type_vee=4 ! for rotatevee_hu 270 end if 271 272 273 ! if(useylm==1.and.opt_diag/=1) MSG_ERROR("useylm==1 and opt_diag/=0 is not possible") 274 if(hu(1)%jpawu_zero.and.nsppol==2) nsppol_imp=2 ! J=0 and nsppol=2 275 if(.not.hu(1)%jpawu_zero.or.nsppol/=2) nsppol_imp=1 ! J/=0 ou nsppol=1 276 ! ================================================================= 277 ! Compute LDA Green's function to compare to weiss_for_rot (check) 278 ! ================================================================= 279 ! call init_green(greenlda,paw_dmft,opt_oper_ksloc=3) 280 ! call greenldacompute_green(cryst_struc,greenlda,pawang,paw_dmft) 281 !! call copy_green(greenlda,weiss_for_rot,2) 282 283 ! ================================================================= 284 ! Compute atomic levels 285 ! ================================================================= 286 call init_oper(paw_dmft,energy_level,opt_ksloc=3) 287 288 ! Compute atomic levels in Slm basis 289 ! ---------------------------------- 290 call compute_levels(cryst_struc,energy_level,self%hdc,pawang,paw_dmft,nondiag=nondiaglevels) 291 292 ! If levels are not diagonal, then diagonalize it (according to 293 ! dmftctqmc_basis) 294 ! ------------------------------------------------ 295 if(paw_dmft%dmftctqmc_basis==1) then 296 if(nondiaglevels.or.useylm==1) then 297 opt_diag=1 298 write(message,'(3a)') ch10, " == Hamiltonian in local basis is non diagonal: diagonalise it",ch10 299 else 300 opt_diag=0 301 write(message,'(5a)') ch10, " == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 & 302 & ," CTQMC will use this basis",ch10 303 end if 304 else if (paw_dmft%dmftctqmc_basis==2) then 305 if(nondiaglevels.or.useylm==1) then 306 write(message,'(7a)') ch10, " == Hamiltonian in local basis is non diagonal",ch10, & 307 & " == According to dmftctqmc_basis: diagonalise density matrix",ch10, & 308 & " == Warning : Check that the Hamiltonian is diagonal !",ch10 309 opt_diag=2 310 else 311 write(message,'(5a)') ch10, " == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 & 312 & ," CTQMC will use this basis",ch10 313 opt_diag=0 314 end if 315 else if (paw_dmft%dmftctqmc_basis==0) then 316 if(nondiaglevels) then 317 write(message,'(4a)') ch10, " == Hamiltonian in local basis is non diagonal",ch10, & 318 & " == According to dmftctqmc_basis: keep this non diagonal basis for the calculation" 319 else 320 write(message,'(5a)') ch10, " == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 & 321 & ," CTQMC will use this basis",ch10 322 end if 323 opt_diag=0 324 end if 325 call wrtout(std_out,message,'COLL') 326 327 if(opt_diag==1) then 328 write(std_out,*) " == The atomic levels are diagonalized" 329 else if(opt_diag==2) then 330 write(std_out,*) " == The correlated occupation matrix is diagonalized" 331 end if 332 333 ! ================================================================= 334 ! Now, check if diagonalisation is necessary 335 ! ================================================================= 336 337 338 ! ================================================================= 339 ! First rotate to Ylm basis the atomic levels 340 ! ================================================================= 341 342 if(useylm==1) then 343 344 ! Rotate from Slm to Ylm the atomic levels 345 ! ---------------------------------------- 346 call slm2ylm_matlu(energy_level%matlu,natom,1,pawprtvol) 347 348 ! Print atomic energy levels in Ylm basis 349 ! -------------------------------- 350 if(pawprtvol>=3) then 351 write(message,'(a,a)') ch10, " == Print Energy levels in Ylm basis" 352 call wrtout(std_out,message,'COLL') 353 call print_matlu(energy_level%matlu,natom,1) 354 end if 355 356 end if ! useylm 357 358 ! =========================================================================================== 359 ! Start for diagonalization of levels/density matrix according to opt_diag 360 ! =========================================================================================== 361 !opt_rot=2 ! do it one time before CTQMC 362 opt_rot=1 ! do all the rotations successively on all different quantities. 363 if(opt_diag==1.or.opt_diag==0) then 364 365 366 if(opt_diag==1) then 367 ! ================================================================= 368 ! Diagonalize atomic levels 369 ! ================================================================= 370 ABI_DATATYPE_ALLOCATE(level_diag,(natom)) 371 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,level_diag) 372 373 ! Diagonalise atomic levels (opt_real is necessary, because 374 ! rotation must be real in order for the occupations and Green's 375 ! function to be real) 376 ! --------------------------------------------------------------- 377 call diag_matlu(energy_level%matlu,level_diag,natom,& 378 & prtopt=pawprtvol,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,optreal=1) 379 380 ! call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1) 381 ! write(message,'(a,2x,a,f13.5)') ch10,& 382 !& " == Print first Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie 383 ! call wrtout(std_out,message,'COLL') 384 ! call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1) 385 386 if(opt_rot==1) call copy_matlu(level_diag,energy_level%matlu,natom) 387 388 389 call destroy_matlu(level_diag,natom) 390 ABI_DATATYPE_DEALLOCATE(level_diag) 391 392 ! Print diagonalized levels 393 ! -------------------------- 394 if(pawprtvol>=3) then 395 write(message,'(a,2x,a,f13.5)') ch10,& 396 & " == Print Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie 397 call wrtout(std_out,message,'COLL') 398 call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1) 399 else 400 write(message,'(a,2x,a,f13.5)') ch10,& 401 & " == Energy levels Diagonalized for Fermi Level=",paw_dmft%fermie 402 call wrtout(std_out,message,'COLL') 403 end if 404 call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee) 405 else if (opt_diag==0) then 406 do iatom=1,cryst_struc%natom 407 lpawu=paw_dmft%lpawu(iatom) 408 itypat=cryst_struc%typat(iatom) 409 if(lpawu/=-1) then 410 ! write(6,*) size(udens_atoms(iatom)%value) 411 ! write(6,*) size(hu(itypat)%udens) 412 ! write(6,*) udens_atoms(iatom)%value 413 ! write(6,*) hu(itypat)%udens 414 udens_atoms(iatom)%value=hu(itypat)%udens 415 end if 416 end do 417 end if 418 ! call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms) 419 420 else if(opt_diag==2) then 421 ! ================================================================= 422 ! Diagonalizes density matrix and keep eigenvectors in eigvectmatlu 423 ! ================================================================= 424 425 ! Print density matrix before diagonalization 426 ! ------------------------------------------- 427 if(pawprtvol>=3) then 428 write(message,'(a,2x,a)') ch10, " == Density Matrix before diagonalisation =" 429 call wrtout(std_out,message,'COLL') 430 !MGNAG: This call is wrong if green has intent(out), now we use intent(inout) 431 call print_matlu(green%occup%matlu,natom,1) 432 end if 433 434 !! checkstop: we can have two different diagonalisation basis for the up and dn 435 !! but one use the same basis, unless the error is really to large(>0.1) 436 437 ! Diagonalize density matrix 438 ! --------------------------- 439 call diag_matlu(green%occup%matlu,dmat_diag,natom,& 440 & prtopt=4,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,checkstop=.false.) 441 442 ! Print diagonalized density matrix 443 ! ---------------------------------- 444 if(pawprtvol>=3) then 445 write(message,'(a,2x,a)') ch10,& 446 & " == Diagonalized Density Matrix in the basis used for QMC =" 447 call wrtout(std_out,message,'COLL') 448 call print_matlu(dmat_diag,natom,1) 449 450 !write(message,'(2a,i3,13x,a)') ch10,' == Rotation of interaction matrix ==' 451 !call wrtout(std_out,message,'COLL') 452 end if 453 454 !if (.not.hu(1)%jpawu_zero) & 455 !MSG_WARNING("In qmc_prep_ctqmc J/=0 and rotation matrix not rotated") 456 ! Rotate interaction. 457 ! call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms) 458 call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee) 459 460 end if 461 ! =========================================================================================== 462 ! END Of diagonalization 463 ! =========================================================================================== 464 465 call flush_unit(std_out) 466 467 ! =========================================================================================== 468 ! Broadcast matrix of rotation from processor 0 to the other 469 ! In case of degenerate levels, severals rotations are possible. Here we 470 ! choose the rotation of proc 0. It is arbitrary. 471 ! =========================================================================================== 472 do iatom=1,cryst_struc%natom 473 lpawu=paw_dmft%lpawu(iatom) 474 if(lpawu/=-1) then 475 tndim=nspinor*(2*lpawu+1) 476 do isppol=1,nsppol 477 call xmpi_bcast(eigvectmatlu(iatom,isppol)%value,0,paw_dmft%spacecomm,ier) 478 end do 479 end if 480 end do 481 482 483 !unitnb=300000+paw_dmft%myproc 484 !call int2char4(paw_dmft%myproc,tag_proc) 485 !tmpfil = 'eigvectmatluaftermpi'//tag_proc 486 !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted') 487 !do iflavor1=1,14 488 ! do iflavor2=1,14 489 ! write(unitnb,*) iflavor1,iflavor2,eigvectmatlu(1,1)%value(iflavor1,iflavor2) 490 ! enddo 491 !enddo 492 493 ! =========================================================================================== 494 ! Now rotate various quantities in the new basis 495 ! =========================================================================================== 496 497 !======================================================= 498 ! Allocate, Compute, and Rotate atomic levels for CTQMC 499 !======================================================= 500 501 ! If levels not rotated, rotate them 502 ! ----------------------------------- 503 if(opt_diag==2.and.opt_rot==1) call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1) 504 505 ! Print atomic levels 506 ! ------------------- 507 if(pawprtvol>=3) then 508 write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels after rotation" 509 call wrtout(std_out,message,'COLL') 510 call print_matlu(energy_level%matlu,natom,1) 511 else 512 write(message,'(a,2x,a,f13.5)') ch10," == CT-QMC Energy levels rotated" 513 call wrtout(std_out,message,'COLL') 514 end if 515 516 !==================================================================== 517 ! If levels were diagonalized before, then rotate density matrix for 518 ! information. 519 !==================================================================== 520 if(opt_diag==1) then 521 ABI_DATATYPE_ALLOCATE(matlu1,(natom)) 522 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 523 call copy_matlu(green%occup%matlu,matlu1,natom) 524 if(pawprtvol>=3) then 525 write(message,'(a,2x,a)') ch10,& 526 & " == Occupations before rotations" 527 call wrtout(std_out,message,'COLL') 528 call print_matlu(green%occup%matlu,natom,1) 529 end if 530 531 ! 1) rotate density matrix to Ylm basis 532 ! -------------------------------------- 533 if(useylm==1) then 534 call slm2ylm_matlu(matlu1,natom,1,pawprtvol) 535 if(pawprtvol>=3) then 536 write(message,'(a,a)') ch10, " == Print occupations in Ylm basis" 537 call wrtout(std_out,message,'COLL') 538 call print_matlu(matlu1,natom,1) 539 end if 540 end if 541 542 ! 2) rotate density matrix to rotated basis 543 ! ------------------------------------------- 544 if(opt_rot==1.or.opt_rot==2) call rotate_matlu(matlu1,eigvectmatlu,natom,3,1) 545 write(message,'(a,2x,a,f13.5)') ch10," == Rotated occupations (for information)" 546 call wrtout(std_out,message,'COLL') 547 call print_matlu(matlu1,natom,1,compl=1) 548 call checkreal_matlu(matlu1,natom,tol10) 549 call destroy_matlu(matlu1,natom) 550 ABI_DATATYPE_DEALLOCATE(matlu1) 551 552 end if 553 554 call flush_unit(std_out) 555 556 557 ! ================================================================= 558 ! Rotate weiss function according to eigenvectors. 559 ! ================================================================= 560 !!!stop 561 ! Rotate Weiss function first in Ylm basis 562 ! ------------------------------------------------------------------- 563 if(useylm==1) then 564 write(message,'(a,2x,a)') ch10, " == Rotation of weiss and greenlda in the Ylm Basis=" 565 call wrtout(std_out,message,'COLL') 566 do ifreq=1,paw_dmft%dmft_nwlo 567 call slm2ylm_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,1,0) 568 call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,1,0) 569 ! call slm2ylm_matlu(greenlda%oper(ifreq)%matlu,natom,1,0) 570 end do 571 end if 572 573 if(pawprtvol>=3) then 574 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 575 ! " == Print weiss for small freq 1 before rot" ! debug 576 ! call wrtout(std_out,message,'COLL') ! debug 577 ! call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) ! debug 578 579 ! Print Weiss function 580 ! -------------------- 581 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 582 & " == Print weiss for 1st freq before rot" ! debug 583 call wrtout(std_out,message,'COLL') ! debug 584 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) ! debug 585 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 586 & " == Print weiss for last freq before rot" ! debug 587 call wrtout(std_out,message,'COLL') ! debug 588 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 589 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 590 !& " == Print LDA G for 1st freq before rot" ! debug 591 ! call wrtout(std_out,message,'COLL') ! debug 592 ! call print_matlu(greenlda%oper(1)%matlu,natom,1,compl=1,opt_exp=2) ! debug 593 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 594 !& " == Print LDA G for last freq before rot" ! debug 595 ! call wrtout(std_out,message,'COLL') ! debug 596 ! call print_matlu(greenlda%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) ! debug 597 end if 598 599 if(opt_diag/=0) then 600 ! Rotate Weiss function first in Ylm basis then in the rotated basis. 601 ! ------------------------------------------------------------------- 602 write(message,'(a,2x,a)') ch10, " == Rotation of weiss =" 603 call wrtout(std_out,message,'COLL') 604 do ifreq=1,paw_dmft%dmft_nwlo 605 if(opt_rot==1) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 606 if(opt_rot==1) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 607 ! call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6) 608 end do 609 610 if(paw_dmft%myproc .eq. mod(nproc+1,nproc)) then 611 if (open_file(trim(paw_dmft%filapp)//"_atom__G0w_.dat", message, newunit=unt) /= 0) then 612 MSG_ERROR(message) 613 end if 614 do ifreq=1,paw_dmft%dmft_nwlo 615 write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),& 616 & (((weiss_for_rot%oper(ifreq)%matlu(1)%mat(im1,im1,isppol,ispinor,ispinor),& 617 & im1=1,3),ispinor=1,nspinor),isppol=1,nsppol) 618 end do 619 close(unt) 620 end if 621 622 call flush_unit(std_out) 623 if(pawprtvol>=3) then 624 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 625 & " == Print weiss for small freq 1 after rot" ! debug 626 call wrtout(std_out,message,'COLL') ! debug 627 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) ! debug 628 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 629 & " == Print weiss for last freq after rot" ! debug 630 call wrtout(std_out,message,'COLL') ! debug 631 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 632 end if 633 634 ! ! Rotate LDA Green's function first in Ylm basis then in the rotated basis and compare to weiss_for_rot 635 ! ! ----------------------------------------------------------------------------------------------------- 636 ! write(message,'(a,2x,a)') ch10, " == Rotation of greenlda =" 637 ! call wrtout(std_out,message,'COLL') 638 ! do ifreq=1,paw_dmft%dmft_nwlo 639 ! if(opt_rot==1) call rotate_matlu(greenlda%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 640 ! call diff_matlu("Weiss_for_rot","greenlda",weiss_for_rot%oper(ifreq)%matlu,greenlda%oper(ifreq)%matlu,natom,1,tol14) 641 !! call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6) 642 ! end do 643 ! if(pawprtvol>=3) then 644 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 645 !& " == Print greenlda for small freq 1 after rot" ! debug 646 ! call wrtout(std_out,message,'COLL') ! debug 647 ! call print_matlu(greenlda%oper(1)%matlu,natom,1,compl=1,opt_exp=2) ! debug 648 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 649 !& " == Print greenlda for last freq after rot" ! debug 650 ! call wrtout(std_out,message,'COLL') ! debug 651 ! call print_matlu(greenlda%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) ! debug 652 ! end if 653 ! call flush_unit(std_out) 654 end if 655 656 ! ================================================================= 657 ! Compute analytic limit of hybridization and rotate it 658 ! ================================================================= 659 ABI_DATATYPE_ALLOCATE(hybri_coeff,(paw_dmft%natom)) 660 call init_matlu(paw_dmft%natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,hybri_coeff) 661 !write(6,*)"hybri1",hybri_coeff(1)%mat(1,1,1,1,1),paw_dmft%natom,cryst_struc%natom 662 663 ! Compute analytical C_ij such that F_ij -> C_ij/iw_n 664 ! --------------------------------------- 665 call hybridization_asymptotic_coefficient(cryst_struc,paw_dmft,pawang,hybri_coeff) 666 write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n for large frequency" 667 call wrtout(std_out,message,'COLL') 668 669 ! Print analytical C_ij (not rotated) 670 ! --------------------------------------- 671 call print_matlu(hybri_coeff,natom,1) 672 673 ! Rotate analytical C_ij in Ylm basis 674 ! --------------------------------------- 675 if(useylm==1) call slm2ylm_matlu(hybri_coeff,natom,1,pawprtvol) 676 if(opt_diag/=0) then 677 678 ! Rotate analytical C_ij in rotated basis 679 ! --------------------------------------- 680 if(opt_rot==1.or.opt_rot==2) call rotate_matlu(hybri_coeff,eigvectmatlu,natom,3,1) 681 682 ! Print analytical C_ij (rotated) 683 ! --------------------------------------- 684 write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n after rotation" 685 call wrtout(std_out,message,'COLL') 686 call print_matlu(hybri_coeff,natom,1,compl=1,opt_exp=1) 687 end if 688 689 ! ================================================================= 690 ! Check if rotation is properly done. 691 ! ================================================================= 692 if(3==4) then 693 write(message,'(a,2x,a)') ch10,& 694 & " == Print dmat before rot" 695 call wrtout(std_out,message,'COLL') 696 call print_matlu(green%occup%matlu,natom,1) 697 if(useylm==1) call slm2ylm_matlu(green%occup%matlu,natom,1,pawprtvol) 698 if(opt_rot==1) call rotate_matlu(green%occup%matlu,eigvectmatlu,natom,3,1) 699 write(message,'(a,2x,a)') ch10,& 700 & " == Print dmat after rot" 701 call wrtout(std_out,message,'COLL') 702 call print_matlu(green%occup%matlu,natom,1) 703 704 write(message,'(2a)') ch10,' QMC STOP: DEBUG' 705 call wrtout(std_out,message,'COLL') 706 MSG_ERROR(message) 707 end if 708 ! ================================================================= 709 ! Check 710 ! ================================================================= 711 712 ! write(message,'(a,2x,a,f13.5)') ch10,& 713 !& " == Print weiss for small tau" 714 ! call wrtout(std_out,message,'COLL') 715 ! call print_matlu(weiss%oper(1)%matlu,natom,1) 716 ! write(message,'(a,2x,a,f13.5)') ch10,& 717 !& " == Print weiss for large tau" 718 ! call wrtout(std_out,message,'COLL') 719 ! call print_matlu(weiss%oper(paw_dmft%dmft_nwlo)%matlu,natom,1) 720 ! call flush_unit(std_out) 721 ! write(message,'(2a)') ch10,' Check weiss_for_rot(last freq)' 722 ! call wrtout(std_out,message,'COLL') 723 ! call checkdiag_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,tol6,opt=nspinor) 724 ! call flush_unit(std_out) 725 ! write(message,'(2a)') ch10,' Check weiss_for_rot(ifreq=1)' 726 ! call wrtout(std_out,message,'COLL') 727 ! call checkdiag_matlu(weiss_for_rot%oper(1)%matlu,natom,tol6,opt=nspinor) 728 ! call flush_unit(std_out) 729 730 master=0 731 732 ! ================================================================= 733 ! Print out 734 ! ================================================================= 735 736 ! Print Weiss 737 ! ------------- 738 if(paw_dmft%dmft_prgn==1) then 739 call print_green('Weiss_diag',weiss_for_rot,1,paw_dmft,pawprtvol=1,opt_wt=1,opt_decim=1) 740 end if 741 742 write(message,'(a,2x,a,f13.5)') ch10,& 743 & " == Preparing data for CTQMC" 744 call wrtout(std_out,message,'COLL') 745 746 ! Print Rotate Weiss for 1st and last frequencies 747 ! ------------------------------------------------ 748 if (pawprtvol>=3) then 749 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 750 & " == Print rotated weiss function for small freq in the rotated basis" ! debug 751 call wrtout(std_out,message,'COLL') ! debug 752 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) ! debug 753 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 754 & " == Print rotated weiss function for largest freq in the rotated basis" ! debug 755 call wrtout(std_out,message,'COLL') ! debug 756 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 757 end if 758 759 ! ================================================================= 760 ! VARIABLES FOR CTQMC TESTS 761 testcode = 0 762 testrot = 0 763 opt_fk=0 ! for developpers to check Fourier transform and computes G0(tau) 764 opt_fk=1 ! usual case: for real calculations 765 ! ================================================================= 766 767 ! ___________________________________________________________________________________ 768 ! 769 ! SECOND PART : BUILT HYBRIDIZATION FROM G0 770 ! ___________________________________________________________________________________ 771 ! 772 ! =========================================================================================== 773 ! Compute inverse of weiss and compute hybridization 774 ! =========================================================================================== 775 776 ! Compute inverse of weiss for each Frequency 777 ! ---------------------------------------------- 778 do ifreq=1,paw_dmft%dmft_nwlo 779 ABI_DATATYPE_ALLOCATE(matlu1,(natom)) 780 ABI_DATATYPE_ALLOCATE(matlu2,(natom)) 781 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 782 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2) 783 784 call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,natom) 785 786 ! Print G_0(iw_n) 787 ! ---------------- 788 if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"go",60000,imre=1) 789 790 ! Compute G_0^-1 791 ! ------------------------------------------- 792 ! if opt_fk=1 or testcode/=0 Do the inversion 793 ! if opt_fk=0 Do not inverse. 794 ! If testcode=2 and opt_fk=0 Do the inversion 795 ! If testcode=1 and opt_fk=0 Do the inversion but no effect, because it will nevertheless be erased 796 ! If opt_fk=1 Do the inversion 797 ! ------------------------------------------- 798 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"weiss",12000,imre=1) 799 if(opt_fk==1.or.testcode/=0) call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1) 800 801 ! Print G_0^-1(iw_n) 802 ! ---------------- 803 if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"goinv",70000,imre=1) 804 805 if(pawprtvol>=4.or.ifreq==paw_dmft%dmft_nwlo) then 806 if(opt_fk==1.or.testcode/=0) then 807 ! Check inversion : do the product 808 ! ---------------------------------------------- 809 call prod_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,matlu2,natom) 810 write(message,'(a,2x,a,i7)') ch10,& ! debug 811 & " == Print product of weiss times invers for freq",ifreq 812 call wrtout(std_out,message,'COLL') ! debug 813 call print_matlu(matlu2,natom,1) ! debug 814 end if 815 end if 816 817 call destroy_matlu(matlu1,natom) 818 call destroy_matlu(matlu2,natom) 819 ABI_DATATYPE_DEALLOCATE(matlu1) 820 ABI_DATATYPE_DEALLOCATE(matlu2) 821 end do 822 823 ! Copy weiss_for_rot into weiss 824 ! ------------------------------- 825 !call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,weiss%oper(ifreq)%matlu,natom) 826 827 828 ! Print G_0^-1 for 1st and last frequencies. 829 ! ----------------------------------------- 830 if(pawprtvol>=3) then 831 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 832 & " == Print G_0^-1 for small freq in the rotated basis" ! debug 833 call wrtout(std_out,message,'COLL') ! debug 834 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) ! debug 835 write(message,'(a,2x,a,e18.10,a)') ch10,& ! debug 836 & " == Print G_0^-1 for last freq in the rotated basis (last freq=", paw_dmft%omega_lo(paw_dmft%dmft_nwlo),")" ! debug 837 call wrtout(std_out,message,'COLL') ! debug 838 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 839 end if 840 841 ! Substract frequency from diagonal part 842 ! ====================================== 843 844 ABI_ALLOCATE(shift,(natom)) 845 do ifreq=1,paw_dmft%dmft_nwlo 846 shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 847 848 ! write(5555,'(400e17.4)') paw_dmft%omega_lo(ifreq),((((((weiss_for_rot%oper(ifreq)%matlu(1)%mat& 849 ! & (im,im1,isppol,ispinor,ispinor1)-cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)),im=1,2*3+1),& 850 !& im1=1,2*3+1),isppol=1,nsppol),ispinor=1,nspinor),ispinor1=1,nspinor) 851 852 ! Compute G_0^-1-iw_n 853 ! -------------------- 854 if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift) 855 856 ! Compute -G_0^-1+iw_n 857 ! -------------------- 858 if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone) 859 860 ! Print -G_0^-1+iw_n 861 ! -------------------- 862 if(optdb==1) then 863 call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"G0inv_minus_omega",20000,imre=1) 864 end if 865 end do 866 867 ! Print -G_0^+1-iw_n=(F-levels) for last freq in the rotated basis" 868 ! ------------------------------------------------------------------ 869 ABI_DEALLOCATE(shift) 870 if(pawprtvol>=3) then 871 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 872 & " == Print G_0^-1-iw_n=-(F-levels) for last freq in the rotated basis" ! debug 873 call wrtout(std_out,message,'COLL') ! debug 874 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 875 end if 876 877 ! Check numerical limit of F(i_wn)*iw_n (can be used also to compute F ) 878 ! ====================================== 879 880 if(opt_nondiag==1) then 881 ABI_DATATYPE_ALLOCATE(matlu1,(natom)) 882 ABI_DATATYPE_ALLOCATE(matlu2,(natom)) 883 ABI_DATATYPE_ALLOCATE(matlu3,(natom)) 884 ABI_DATATYPE_ALLOCATE(matlu4,(natom)) 885 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 886 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2) 887 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu3) 888 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu4) 889 890 write(message,'(a,2x,a,f13.5)') ch10, " == energy_levels" 891 call wrtout(std_out,message,'COLL') 892 call print_matlu(energy_level%matlu,natom,1,opt_exp=2,compl=1) 893 894 do ifreq=paw_dmft%dmft_nwlo,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency. 895 !do ifreq=paw_dmft%dmftqmc_l,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency. 896 ! Compute F (substract levels) for max frequency 897 ! ----------------------------------------------- 898 call add_matlu(weiss_for_rot%oper(ifreq)%matlu,energy_level%matlu,matlu1,natom,-1) 899 900 ! Print F(iw_n)=-(G_0^-1-iw_n+levels) for last frequency. 901 ! -------------------------------------------------------- 902 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 903 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 904 & " == Print F(iw_n)=-(G_0^-1-iw_n+levels) for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 905 call wrtout(std_out,message,'COLL') 906 call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 907 end if 908 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"Hybridization",10000,imre=1) 909 910 ! Put F in weiss_for_rot -> CTQMC 911 ! ------------------------------- 912 if(opt_rot==2) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 913 ! The following line will produce directly the weiss function for the CTQMC code 914 if(opt_fk==1) call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom) 915 916 917 ! Multiply F by frequency 918 ! ------------------------ 919 call copy_matlu(matlu1,matlu2,natom) 920 call fac_matlu(matlu1,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)) 921 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 922 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 923 & " == Print numerical C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 924 call wrtout(std_out,message,'COLL') 925 call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 926 end if 927 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij",72800,imre=1) 928 !call rotate_matlu(matlu1,eigvectmatlu,natom,3,1) 929 930 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 931 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 932 & " == Print numerical after back rotation C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 933 call wrtout(std_out,message,'COLL') 934 call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 935 end if 936 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_rotated",72900,imre=1) 937 938 ! Built C_ij/iw_n 939 ! ------------------------ 940 call copy_matlu(hybri_coeff,matlu1,natom) 941 call fac_matlu(matlu1,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)) 942 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_over_omega",72000) 943 ! if(ifreq==paw_dmft%dmft_nwlo) then 944 ! write(message,'(a,2x,a,f13.5)') ch10, " == Print numerical C_ij/iw_n for frequency",paw_dmft%omega_lo(ifreq) 945 ! call wrtout(std_out,message,'COLL') 946 ! call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 947 ! endif 948 949 ! For test: put C_ij/i_wn into weiss_for_rot 950 ! -------------------------------------------- 951 !call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1) 952 953 ! Compute Hybri - C_ij/iw_n 954 ! ------------------------ 955 call add_matlu(matlu2,matlu1,matlu3,natom,-1) 956 957 ! Print Hybri - C_ij/iw_n 958 ! ------------------------ 959 if(optdb==1) call printplot_matlu(matlu3,natom,paw_dmft%omega_lo(ifreq),"hybri_minus_asymp",74000,imre=1) 960 961 ! Multiply (F-C_ij/i_wn) by (iw_n)**2 to find D_ij such that (F-C_ij/i_wn) -> D_ij/(iw_n)^2 only for last frequency. 962 ! ------------------------------------------------------------------------------------------------------------------ 963 call copy_matlu(matlu3,matlu2,natom) 964 call fac_matlu(matlu2,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2) 965 if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"fminuscijtimesw2",75000,imre=1) 966 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 967 call copy_matlu(matlu2,matlu4,natom) 968 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 969 & " == Print numerical (F(iw_n)-C_ij/iw_n)%iw_n^2 for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 970 call wrtout(std_out,message,'COLL') 971 call print_matlu(matlu4,natom,1) 972 end if 973 974 ! Built C_ij/iw_n+D_ij/(iw_n)^2 975 ! ------------------------ 976 call copy_matlu(matlu4,matlu3,natom,opt_re=1) 977 call fac_matlu(matlu3,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2) 978 call add_matlu(matlu1,matlu3,matlu2,natom,1) 979 if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"cij_w_plus_dij_w2",72700,imre=1) 980 ! For test: put C_ij/i_wn +D_ij/(iw_n)^2 into weiss_for_rot 981 ! -------------------------------------------- 982 !call copy_matlu(matlu2,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1) 983 984 985 end do 986 987 call destroy_matlu(matlu1,natom) 988 call destroy_matlu(matlu2,natom) 989 call destroy_matlu(matlu3,natom) 990 call destroy_matlu(matlu4,natom) 991 ABI_DATATYPE_DEALLOCATE(matlu1) 992 ABI_DATATYPE_DEALLOCATE(matlu2) 993 ABI_DATATYPE_DEALLOCATE(matlu3) 994 ABI_DATATYPE_DEALLOCATE(matlu4) 995 end if 996 997 ! ========================================================================================= 998 ! Start big loop to compute hybridization 999 ! ========================================================================================= 1000 do iatom=1,cryst_struc%natom 1001 green%ecorr_qmc(iatom)=zero 1002 itypat=cryst_struc%typat(iatom) 1003 lpawu=paw_dmft%lpawu(iatom) 1004 tndim=2*lpawu+1 1005 if(lpawu/=-1) then 1006 1007 nflavor=2*(tndim) 1008 if(testcode>=1) then 1009 nflavor=2 1010 if(testcode==2) then 1011 ispa=1 1012 ispb=2 1013 if(nspinor==1) ispb=1 1014 ima=1 1015 imb=1 1016 if(tndim>4) then 1017 ima=5 ! row 1018 imb=4 ! column 1019 end if 1020 end if 1021 end if 1022 1023 !allocate(correl_loc(nflavor,nflavor)) 1024 !ABI_ALLOCATE(f_with_k,(MIN(paw_dmft%dmft_nwli,paw_dmft%dmftqmc_l),nflavor)) 1025 ABI_ALLOCATE(fw1,(paw_dmft%dmft_nwlo,nflavor)) 1026 ABI_ALLOCATE(fw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor)) 1027 ABI_ALLOCATE(levels_ctqmc,(nflavor)) 1028 ABI_ALLOCATE(levels_ctqmc_nd,(nflavor,nflavor)) 1029 levels_ctqmc_nd=czero 1030 ABI_ALLOCATE(hybri_limit,(nflavor,nflavor)) 1031 hybri_limit=czero 1032 fw1_nd=czero 1033 fw1=czero 1034 !allocate(fw2(paw_dmft%dmft_nwli)) 1035 ! ================================================================= 1036 ! Compute Hybridization 1037 ! ================================================================= 1038 1039 if (testcode==0) then 1040 iflavor1=0 1041 iflavor2=0 1042 do isppol=1,nsppol 1043 do ispinor1=1,nspinor 1044 do ispinor2=1,nspinor 1045 do im1=1,tndim 1046 do im2=1,tndim 1047 ! first diagonal terms whatever opt_nondiag 1048 iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1) 1049 iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1) 1050 !write(6,*) isppol,ispinor1,ispinor2,im1,im2 1051 !write(6,*) iflavor1,iflavor2 1052 1053 if ( iflavor1==iflavor2 ) then 1054 1055 ! Do spline of weiss function for im and isppol 1056 ! Construction of fw1 1057 ! ================================================================= 1058 1059 do ifreq=1,paw_dmft%dmft_nwlo 1060 if(opt_fk==1) then 1061 fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1) 1062 else if (opt_fk==0) then 1063 fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1) 1064 end if 1065 end do 1066 fw1_nd(:,iflavor1,iflavor1)=fw1(:,iflavor1) 1067 1068 levels_ctqmc(iflavor1)=real(energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1),kind=dp) 1069 hybri_limit(iflavor1,iflavor1)=hybri_coeff(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1) 1070 1071 1072 if(nsppol==1.and.nspinor==1) then 1073 !f_with_k(:,iflavor+tndim)=f_with_k(:,iflavor) 1074 fw1(:,iflavor1+tndim)=fw1(:,iflavor1) 1075 fw1_nd(:,iflavor1+tndim,iflavor1+tndim)=fw1(:,iflavor1) 1076 levels_ctqmc(iflavor1+tndim)=levels_ctqmc(iflavor1) 1077 hybri_limit(iflavor1+tndim,iflavor1+tndim)=hybri_limit(iflavor1,iflavor1) 1078 end if 1079 1080 else 1081 1082 do ifreq=1,paw_dmft%dmft_nwlo 1083 if(opt_fk==1) then 1084 fw1_nd(ifreq,iflavor1,iflavor2)= & 1085 & weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 1086 else if (opt_fk==0) then 1087 fw1_nd(ifreq,iflavor1,iflavor2)= & 1088 & weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 1089 end if 1090 1091 ! omega=pi*paw_dmft%temp*(two*float(ifreq)-1) 1092 !write(3333,*) omega,weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 1093 !write(4444,*) omega,fw1_nd(ifreq,iflavor1,iflavor2),"#",iflavor1,iflavor2 1094 ! if(iflavor1/=iflavor2)write(5555,*) omega,imag(fw1_nd(ifreq,iflavor1,iflavor2)),"#",iflavor1,iflavor2 1095 1096 end do 1097 hybri_limit(iflavor1,iflavor2)=hybri_coeff(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 1098 1099 ! write(3333,*) 1100 ! write(4444,*) 1101 ! write(5555,*) 1102 1103 if(nsppol==1.and.nspinor==1) then 1104 fw1_nd(:,iflavor1+tndim,iflavor2+tndim) = fw1_nd(:,iflavor1,iflavor2) 1105 hybri_limit(iflavor1+tndim,iflavor2+tndim)=hybri_limit(iflavor1,iflavor2) 1106 end if 1107 1108 end if 1109 1110 ! < / HACK > 1111 end do !im2 1112 end do !im1 1113 end do !ispinor2 1114 end do !ispinor1 1115 end do !isppol 1116 ! < HACK > 1117 ! JB. On 1000 cpus this can not work since all CPU try to open/write the files 1118 ! Action : Don't print it or check only one cpu does it. 1119 if(pawprtvol>=10000000) then 1120 write(message,'(a,2x,a)') ch10, " == Hybri for all flavors for CTQMC " 1121 call wrtout(std_out,message,'COLL') 1122 do iflavor1=1,nflavor 1123 write(message,'(4x,14(2e14.5,2x))') (hybri_limit(iflavor1,iflavor2),iflavor2=1,nflavor) 1124 call wrtout(std_out,message,'COLL') 1125 end do 1126 1127 if (open_file('Hybri_cijoveromega',message, newunit=unt, status='unknown', form='formatted') /= 0) then 1128 MSG_ERROR(message) 1129 end if 1130 if (open_file('Hybri',message,newunit=unt2,status='unknown',form='formatted') /= 0) then 1131 MSG_ERROR(message) 1132 end if 1133 do ifreq=1,paw_dmft%dmft_nwlo 1134 ! weiss_for_rot is G_0^-1-iw_n=-(F-levels) 1135 if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"weissbefore112",30000) 1136 end do 1137 do iflavor1=1,nflavor 1138 do iflavor2=1,nflavor 1139 do ifreq=1,paw_dmft%dmft_nwlo 1140 omega=pi*paw_dmft%temp*(two*float(ifreq)-1) 1141 ! fw1_nd is -G_0^+1-iw_n=(F-levels) 1142 write(unt,'(300e16.5)') paw_dmft%omega_lo(ifreq)& 1143 & ,fw1_nd(ifreq,iflavor1,iflavor2)-hybri_limit(iflavor1,iflavor2)/cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp) 1144 write(unt2,'(300e16.5)') paw_dmft%omega_lo(ifreq),fw1_nd(ifreq,iflavor1,iflavor2) 1145 !write(1111,*) omega,real(fw1_nd(ifreq,iflavor1,iflavor2)) 1146 !write(1112,*) omega,imag(fw1_nd(ifreq,iflavor1,iflavor2)) 1147 end do 1148 write(unt,*) 1149 write(unt2,*) 1150 ! write(1111,*) 1151 ! write(1112,*) 1152 end do 1153 end do 1154 close(unt) 1155 close(unt2) 1156 end if 1157 end if ! testcode 1158 ! </ HACK > 1159 1160 ! ==================================================================================== 1161 ! TEST 1162 ! For testing purpose, built ultra simple hybridization (constant in 1163 ! imaginary time or very simple) or extract some part of the calculated hybridization 1164 ! ==================================================================================== 1165 1166 if(testcode>=1) then 1167 dmft_nwlo=paw_dmft%dmft_nwlo 1168 paw_dmft%dmft_nwlo=paw_dmft%dmftqmc_l 1169 ABI_ALLOCATE(gw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor)) 1170 gw1_nd=czero 1171 1172 ! Call testcode_ctqmc: built simple hybridization 1173 !-------------------------------------------------- 1174 if (testcode==1) then 1175 call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,& 1176 & levels_ctqmc,hybri_limit,nflavor,1,paw_dmft%temp,testrot,testcode,umod) 1177 ! Select 2x2 hybridization matrix from the current larger matrix 1178 ! ima and imb are defined above. 1179 !---------------------------------------------------------------- 1180 else if (testcode==2) then 1181 facnd=0.8d0 1182 facd=1.0d0 1183 !write(6,*) "fac",facnd,facd 1184 levels_ctqmc_nd(2,2) = energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb) 1185 levels_ctqmc_nd(1,1) = energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa) 1186 levels_ctqmc(2) = real(energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb),kind=dp) 1187 levels_ctqmc(1) = real(energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa),kind=dp) 1188 if(opt_diag/=1) then 1189 levels_ctqmc_nd(1,2) = energy_level%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1190 levels_ctqmc_nd(2,1) = energy_level%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1191 end if 1192 hybri_limit(1,1) = facd*hybri_coeff(iatom)%mat(ima,ima,1,ispa,ispa) 1193 hybri_limit(2,2) = facd*hybri_coeff(iatom)%mat(imb,imb,1,ispb,ispb) 1194 hybri_limit(1,2) = facnd*hybri_coeff(iatom)%mat(ima,imb,1,ispa,ispb) 1195 hybri_limit(2,1) = facnd*hybri_coeff(iatom)%mat(imb,ima,1,ispb,ispa) 1196 !write(6,*) "hybri_limit",hybri_limit 1197 !write(6,*) "levels_ctqmc",levels_ctqmc 1198 umod=zero 1199 1200 tmpfil = 'fw1_nd_re' 1201 !if (open_file(newunit=unt,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then 1202 ! MSG_ERROR(message) 1203 !end if 1204 tmpfil = 'fw1_nd_im' 1205 !if (open_file(newunit=unt2,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then 1206 ! MSG_ERROR(message) 1207 !end if 1208 write(std_out,*) "testcode==2",ispa,ispb,ima,imb 1209 write(std_out,*) "opt_fk==",opt_fk 1210 do ifreq=1,paw_dmft%dmftqmc_l 1211 if (opt_fk==1) then 1212 fw1_nd(ifreq,1,1) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa) 1213 fw1_nd(ifreq,2,2) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb) 1214 !fw1_nd(ifreq,1,2) = weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1215 !fw1_nd(ifreq,2,1) = weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1216 fw1_nd(ifreq,1,2) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1217 fw1_nd(ifreq,2,1) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1218 omega=pi*paw_dmft%temp*(two*float(ifreq)-1) 1219 else if (opt_fk==0) then 1220 fw1_nd(ifreq,1,1) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa) 1221 fw1_nd(ifreq,2,2) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb) 1222 fw1_nd(ifreq,1,2) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1223 fw1_nd(ifreq,2,1) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1224 call xginv(fw1_nd(ifreq,:,:),2) 1225 end if 1226 end do 1227 !close(unt) 1228 !close(unt2) 1229 end if 1230 1231 ! Calculation of Inverse Green's function from hybridization 1232 !------------------------------------------------------------- 1233 do if1=1,2 1234 do if2=1,2 1235 do ifreq=1,paw_dmft%dmftqmc_l 1236 omega=pi*paw_dmft%temp*(two*float(ifreq)-1) 1237 if(if1==if2) then 1238 gw1_nd(ifreq,if1,if2) = (cmplx(0.d0,omega,kind=dp)-fw1_nd(ifreq,if1,if2)) 1239 else 1240 gw1_nd(ifreq,if1,if2) = (-fw1_nd(ifreq,if1,if2)) 1241 end if 1242 end do 1243 end do 1244 end do 1245 ! Calculation of Green's function (call inverse) 1246 !------------------------------------------------------------- 1247 do ifreq=1,paw_dmft%dmftqmc_l 1248 call xginv(gw1_nd(ifreq,:,:),2) 1249 end do 1250 write(std_out,*) " testctqmc high frequency limit of hybridization",fw1_nd(paw_dmft%dmftqmc_l,:,:) 1251 1252 ! Integrate Green's function 1253 !------------------------------------------------------------- 1254 do if1=1,2 1255 do if2=1,2 1256 do ifreq=1,paw_dmft%dmftqmc_l 1257 omega=pi*paw_dmft%temp*(two*float(ifreq)-1) 1258 !write(999,*) omega,gw1_nd(ifreq,if1,if2) 1259 end do 1260 !write(999,*) 1261 call int_fct(gw1_nd(:,if1,if2),(if1==if2),2,paw_dmft,integral(if1,if2)) ! test_1 1262 ! write(std_out,*) "testctqmc occupations of input Green's function",(integral(if1,if2)) 1263 end do 1264 end do 1265 ! Write Occupations 1266 !------------------------------------------------------------- 1267 do if1=1,2 1268 do if2=1,2 1269 ! write(std_out,*) "testctqmc occupations of input Green's function",(integral(if1,if2)+conjg(integral(if2,if1)))/two 1270 end do 1271 end do 1272 write(std_out,*) "Occupation of model in matrix form" 1273 do if1=1,2 1274 write(std_out,'(2(2f13.5,3x))') ((integral(if1,if2)+conjg(integral(if2,if1)))/two,if2=1,2) 1275 end do 1276 write(std_out,*) "Limit of hybridization " 1277 do if1=1,2 1278 write(std_out,'(2(2f13.5,3x))') (hybri_limit(if1,if2),if2=1,2) 1279 end do 1280 1281 ! If opt_fk=0, give Green's function to CTQMC code instead of 1282 ! hybridization 1283 !------------------------------------------------------------- 1284 if(opt_fk==0) then 1285 fw1_nd=gw1_nd 1286 end if 1287 1288 ABI_DEALLOCATE(gw1_nd) 1289 paw_dmft%dmft_nwlo=dmft_nwlo 1290 end if 1291 1292 call flush_unit(std_out) 1293 ! ================================================================= 1294 1295 ! ___________________________________________________________________________________ 1296 ! 1297 ! THIRD PART : CALL CTQMC 1298 ! ___________________________________________________________________________________ 1299 1300 ! ================================================================= 1301 ! Main calls to CTQMC code in ABINIT 1302 ! ================================================================= 1303 if(paw_dmft%dmft_solv==5) then 1304 1305 write(message,'(a,2x,a)') ch10,& 1306 & " == Initializing CTQMC" 1307 call wrtout(std_out,message,'COLL') 1308 1309 ! Initialisation 1310 ! ================================================================= 1311 nomega=paw_dmft%dmftqmc_l 1312 call CtqmcInterface_init(hybrid,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, & 1313 & paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,& 1314 & std_out,paw_dmft%spacecomm) 1315 1316 ! for non diagonal code 1317 ! call CtqmcInterface_init(hybrid,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, & 1318 !& paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,& 1319 !& paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,& 1320 !& std_out,paw_dmft%spacecomm,opt_nondiag) 1321 1322 ! options 1323 ! ================================================================= 1324 call CtqmcInterface_setOpts(hybrid,& 1325 opt_Fk =opt_fk,& 1326 & opt_order =paw_dmft%dmftctqmc_order ,& 1327 & opt_movie =paw_dmft%dmftctqmc_mov ,& 1328 & opt_analysis=paw_dmft%dmftctqmc_correl,& 1329 & opt_check =paw_dmft%dmftctqmc_check ,& 1330 & opt_noise =paw_dmft%dmftctqmc_grnns ,& 1331 & opt_spectra =paw_dmft%dmftctqmc_mrka ,& 1332 & opt_gmove =paw_dmft%dmftctqmc_gmove ) 1333 write(message,'(a,2x,2a)') ch10,& 1334 & " == Initialization CTQMC done", ch10 1335 call wrtout(std_out,message,'COLL') 1336 end if 1337 1338 if(paw_dmft%dmft_solv>=6) then 1339 ABI_ALLOCATE(gw_tmp_nd,(paw_dmft%dmft_nwli,nflavor,nflavor)) !because size allocation problem with TRIQS paw_dmft%dmft_nwlo must be >= paw_dmft%dmft_nwli 1340 open(unit=505,file=trim(paw_dmft%filapp)//"_Legendre_coefficients.dat", status='unknown',form='formatted') 1341 else 1342 ABI_ALLOCATE(gw_tmp,(paw_dmft%dmft_nwlo,nflavor+1)) 1343 ABI_ALLOCATE(gw_tmp_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor+1)) 1344 ! use gw_tmp to put freq 1345 do ifreq=1,paw_dmft%dmft_nwlo 1346 gw_tmp(ifreq,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 1347 gw_tmp_nd(ifreq,nflavor,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 1348 end do 1349 end if 1350 1351 ABI_ALLOCATE(gtmp,(paw_dmft%dmftqmc_l,nflavor)) 1352 ! THIS IS A BACKUP PLAN. USING paw_dmft%hybrid makes a segfault on TIKAL 1353 ! PSC with MPI only (and max2_open64). paw_dmf%hybrid is corrupted 1354 ! somewhere but I could not find the place in all DMFT routines 1355 ABI_ALLOCATE(gtmp_nd,(paw_dmft%dmftqmc_l,nflavor,nflavor)) 1356 call flush_unit(std_out) 1357 1358 if(testcode==0) then 1359 !unitnb=100000+paw_dmft%myproc 1360 !call int2char4(paw_dmft%myproc,tag_proc) 1361 !tmpfil = 'hybrilimit'//tag_proc 1362 !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted') 1363 !do iflavor1=1,nflavor 1364 ! do iflavor2=1,nflavor 1365 ! write(unitnb,*) iflavor1,iflavor2,hybri_limit(iflavor1,iflavor2) 1366 ! enddo 1367 !enddo 1368 1369 ! ================================================================= 1370 ! CTQMC run Abinit 1371 ! ================================================================= 1372 if(paw_dmft%dmft_solv==5) then 1373 1374 ABI_ALLOCATE(docc,(1:nflavor,1:nflavor)) 1375 docc(:,:) = zero 1376 call CtqmcInterface_run(hybrid,fw1(1:paw_dmft%dmftqmc_l,:),Gtau=gtmp,& 1377 & Gw=gw_tmp,D=docc(:,:),E=green%ecorr_qmc(iatom),& 1378 !& matU=hu(itypat)%udens,opt_levels=levels_ctqmc) 1379 & matU=udens_atoms(iatom)%value,opt_levels=levels_ctqmc) 1380 call data4entropyDMFT_setDocc(paw_dmft%forentropyDMFT,iatom,docc) 1381 ABI_DEALLOCATE(docc) 1382 1383 ! ================================================================= 1384 ! CTQMC run TRIQS 1385 ! ================================================================= 1386 else if (paw_dmft%dmft_solv>=6) then 1387 1388 ! fw1_nd: Hybridation 1389 ! levels_ctqmc: niveaux 1390 ! hu(itypat)%udens(:,:) : U_ij 1391 ! hu(itypat)%u(:,:,:,:) : uijkl 1392 ! temperature : paw_dmft%temp 1393 ! paw_dmft%dmftqmc_l: nombre de points en temps -1 1394 ! paw_dmft%dmftqmc_n: nombre de cycles 1395 ! ?? Quelles sorties: Les fonctions de Green 1396 ! frequence/temps/Legendre. 1397 ! Double occupations ?? <n_i n_j> 1398 ! test n_tau > 2*nfreq => ntau = 2*nfreq + 1 1399 ! for non diagonal code: 1400 ! call CtqmcInterface_run(hybrid,fw1_nd(1:paw_dmft%dmftqmc_l,:,:),Gtau=gtmp_nd,& 1401 !& Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),& 1402 !& Noise=Noise,matU=hu(itypat)%udens,opt_levels=levels_ctqmc,hybri_limit=hybri_limit) 1403 !Check choice of user to fix model bool var for the solver 1404 if (paw_dmft%dmft_solv==6) then 1405 rot_inv = .false. 1406 else !obviously paw_dmft%dmft_solv==7 with rot invariant terms 1407 rot_inv = .true. 1408 end if 1409 1410 nfreq = paw_dmft%dmft_nwli 1411 !paw_dmft%dmft_nwlo = paw_dmft%dmft_nwli !transparent for user 1412 ntau = paw_dmft%dmftqmc_l !(2*paw_dmft%dmftqmc_l)+1 !nfreq=paw_dmft%dmft_nwli 1413 nleg = paw_dmft%dmftctqmc_triqs_nleg 1414 1415 if ( ntau >= (2*nfreq)+1 ) then 1416 1417 verbosity_solver = paw_dmft%prtvol 1418 beta = 1.0/(paw_dmft%temp*Ha_eV) 1419 1420 !Allocation in/output array phase: 1421 ABI_ALLOCATE(fw1_nd_tmp,(1:nflavor,1:nflavor,1:nfreq)) !column major 1422 ABI_ALLOCATE(g_iw,(1:nflavor,1:nflavor,1:nfreq)) !column major 1423 ! ABI_ALLOCATE(gtau,(1:nflavor,1:nflavor,1:ntau)) !column major 1424 ABI_ALLOCATE(u_mat_ij,(1:nflavor,1:nflavor)) !column major 1425 ABI_ALLOCATE(u_mat_ijkl,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major 1426 ABI_ALLOCATE(u_mat_ijkl_tmp,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major 1427 1428 if ( leg_measure ) then !only if functionality is enabled 1429 ABI_ALLOCATE(gl_nd,(1:nleg,1:nflavor,1:nflavor)) !column major !nl = 30 by default 1430 end if 1431 1432 !Conversion datas Ha -> eV (some duplications for test...) 1433 !fw1_nd_tmp = fw1_nd(1:paw_dmft%dmftqmc_l,:,:) * Ha_eV !fw1_nd = fw1_nd * Ha_eV !Ok? 1434 1435 do iflavor=1,nflavor 1436 do iflavor1=1,nflavor 1437 do ifreq=1,nfreq 1438 fw1_nd_tmp(iflavor,iflavor1,ifreq) = fw1_nd(ifreq,iflavor,iflavor1) * Ha_eV 1439 ! WRITE(500,*) "[IN Fortran] F[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",fw1_nd(ifreq,iflavor,iflavor1) 1440 end do 1441 end do 1442 end do 1443 1444 !Report test 1445 ! WRITE(502,*) hu(itypat)%udens 1446 ! do ifreq=1,paw_dmft%dmftqmc_l 1447 ! write(501,*) ((fw1_nd(ifreq,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor) 1448 ! enddo 1449 !write(866,*)paw_dmft%dmft_nwlo,paw_dmft%dmftqmc_l 1450 !write(866,*) u_mat_ij 1451 ! do iflavor=1,nflavor+1 1452 ! do iflavor1=1,nflavor+1 1453 ! WRITE(502,*) "[OUT Fortran] U(i,j)[ l= ",iflavor," l_= ",iflavor1,"] = ",hu(itypat)%udens(iflavor,iflavor1) 1454 ! enddo 1455 ! enddo 1456 1457 ! if(paw_dmft%myproc==0) then 1458 ! do iflavor=1,nflavor 1459 ! do iflavor1=1,nflavor 1460 ! do iflavor2=1,nflavor 1461 ! do iflavor3=1,nflavor 1462 ! write(490,*), hu(itypat)%vee(iflavor,iflavor1,iflavor2,iflavor3) 1463 ! enddo 1464 ! enddo 1465 ! enddo 1466 ! enddo 1467 ! endif 1468 1469 ! if(paw_dmft%myproc==0) then 1470 ! do iflavor=1,nflavor 1471 ! do iflavor1=1,nflavor 1472 ! write(491,*), hu(itypat)%udens(iflavor,iflavor1) !(1,1,1,1) 1473 ! enddo 1474 ! enddo 1475 ! endif 1476 1477 ! do iflavor=1,nflavor 1478 ! do iflavor1=1,nflavor 1479 ! do iflavor2=1,nflavor 1480 ! do iflavor3=1,nflavor 1481 ! WRITE(552,*), hu(itypat)%vee!(iflavor,iflavor1,iflavor2,iflavor3) 1482 ! enddo 1483 ! enddo 1484 ! enddo 1485 ! enddo 1486 1487 call vee_ndim2tndim_hu_r(lpawu,hu(itypat)%vee,u_mat_ijkl_tmp,1) 1488 do iflavor=1,nflavor 1489 do iflavor1=1,nflavor 1490 do iflavor2=1,nflavor 1491 do iflavor3=1,nflavor 1492 u_mat_ijkl(iflavor,iflavor1,iflavor2,iflavor3) = Ha_eV * u_mat_ijkl_tmp(iflavor,iflavor1,iflavor2,iflavor3) 1493 end do 1494 end do 1495 end do 1496 end do 1497 1498 !u_mat_ijkl = Ha_eV * reshape( u_mat_ijkl , [nflavor,nflavor,nflavor,nflavor] ) !column -> row major + conversion 1499 u_mat_ij = transpose( hu(itypat)%udens ) * Ha_eV !column -> row major + conversion 1500 levels_ctqmc = levels_ctqmc * Ha_eV 1501 1502 !Location array in memory for C++ pointer args to pass 1503 g_iw_ptr = C_LOC( gw_tmp_nd ) !C_LOC( g_iw ) 1504 gtau_ptr = C_LOC( gtmp_nd ) !C_LOC( gtau ) 1505 gl_ptr = C_LOC( gl_nd ) 1506 fw1_nd_ptr = C_LOC( fw1_nd_tmp ) 1507 u_mat_ij_ptr = C_LOC( u_mat_ij ) 1508 u_mat_ijkl_ptr = C_LOC( u_mat_ijkl ) 1509 levels_ptr = C_LOC( levels_ctqmc ) 1510 1511 !Calling interfaced TRIQS solver subroutine from src/01_triqs_ext package 1512 #if defined HAVE_TRIQS 1513 call Ctqmc_triqs_run ( rot_inv, leg_measure, hist, wrt_files, tot_not, & 1514 & nflavor, nfreq, ntau , nleg, int(paw_dmft%dmftqmc_n/paw_dmft%nproc), & 1515 & paw_dmft%dmftctqmc_meas*2*2*nflavor, paw_dmft%dmftqmc_therm, & 1516 & verbosity_solver, paw_dmft%dmftqmc_seed,beta, & 1517 & levels_ptr, u_mat_ij_ptr, u_mat_ijkl_ptr, fw1_nd_ptr, & 1518 & g_iw_ptr, gtau_ptr, gl_ptr, paw_dmft%spacecomm ) 1519 #endif 1520 1521 !WRITE(*,*) "Hello Debug" 1522 !call xmpi_barrier(paw_dmft%spacecomm) !Resynch all processus after calling Impurity solver from TRIQS 1523 1524 !Report output datas from TRIQS to Abinit 1525 !Interacting G(iw) 1526 do ifreq=1,nfreq 1527 do iflavor1=1,nflavor 1528 do iflavor=1,nflavor 1529 ! gw_tmp_nd(ifreq,iflavor,iflavor1) = g_iw(iflavor,iflavor1,ifreq) !* Ha_eV !because 1/ G0(eV) 1530 ! WRITE(503,*) "[OUT Fortran] G(iw)[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",gw_tmp_nd(ifreq,iflavor,iflavor1)!g_iw(iflavor,iflavor1,ifreq) 1531 end do 1532 end do 1533 end do 1534 1535 ! Convert in Ha 1536 gw_tmp_nd = gw_tmp_nd*Ha_eV 1537 1538 ! do iflavor1=1,nflavor 1539 ! do iflavor=1,nflavor 1540 ! 1541 ! WRITE(510,*) "[OUT Fortran] U[ l= ",iflavor," l_= ",iflavor1,"] = ",u_mat_ij(iflavor,iflavor1) 1542 ! enddo 1543 ! enddo 1544 1545 ! if(paw_dmft%myproc==0) write(6,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1) 1546 ! if(paw_dmft%myproc==1) write(6,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1) 1547 ! if(paw_dmft%myproc==0) write(621,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1) 1548 ! if(paw_dmft%myproc==1) write(622,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1) 1549 ! call flush_unit(621) 1550 ! call flush_unit(622) 1551 ! write(message,*) ch10, "essai",paw_dmft%myproc, paw_dmft%myproc,paw_dmft%dmftqmc_seed!gw_tmp_nd(2,1,1) 1552 ! call wrtout(555,message,'PERS',.true.) 1553 ! if(paw_dmft%myproc==0) write(499,*) "essai",paw_dmft%myproc, paw_dmft%dmftqmc_seed 1554 ! if(paw_dmft%myproc==1) write(498,*) "essai",paw_dmft%myproc,paw_dmft%dmftqmc_seed 1555 1556 !Its associated G(tau): Problem of compatibility => paw_dmft%dmftqmc_l < (2*paw_dmft%dmftqmc_l)+1 => We report only paw_dmft%dmftqmc_l = first values of G(tau)... 1557 ! do iflavor=1,nflavor 1558 ! do iflavor1=1,nflavor 1559 ! do itau=1,ntau 1560 ! if ( modulo(itau,2) == 1 ) then !Problem of binding: paw_dmft%dmftqmc_l =! ntau => We take one value by 2 and Write in file all the G(tau) out function from TRIQS 1561 !gtmp_nd(itau,iflavor,iflavor1) = gtau(iflavor,iflavor1,itau) 1562 ! endif 1563 ! if(paw_dmft%myproc==0) then 1564 ! WRITE(504,*) "[OUT Fortran] G[ tau= ",itau," l= ",iflavor," l_= ",iflavor1,"] = ",gtmp_nd(itau,iflavor,iflavor1) !gtmp_nd(itau,iflavor,iflavor1) !passage ok avec ntau/iflavor1/iflavor (iflavor,iflavor1,ntau) 1565 ! endif 1566 ! enddo 1567 ! enddo 1568 ! enddo 1569 1570 !Legendre Polynoms G(L) for extrapolation of Interacting G(iw) by FT, only if leg_measure == TRUE 1571 if (leg_measure) then 1572 do ileg=1,nleg 1573 WRITE(505,*) ileg,((gl_nd(ileg,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor) 1574 end do 1575 close(505) 1576 end if 1577 ! if(paw_dmft%myproc==0) then 1578 ! do itau=1,paw_dmft%dmftqmc_l 1579 ! write(490,*) ((gtmp_nd(itau,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor) 1580 ! enddo 1581 ! endif 1582 ABI_DEALLOCATE( fw1_nd_tmp ) 1583 ABI_DEALLOCATE( g_iw ) 1584 ABI_DEALLOCATE( u_mat_ijkl ) 1585 ABI_DEALLOCATE( u_mat_ijkl_tmp ) 1586 ABI_DEALLOCATE( u_mat_ij ) 1587 ! ========================================================================= 1588 ! Compute Green's function in imaginary freq using Legendre coefficients 1589 ! ========================================================================= 1590 if (leg_measure) then 1591 call xmpi_barrier(paw_dmft%spacecomm) 1592 call flush_unit(std_out) 1593 write(message,'(2a)') ch10," == Compute G(iw_n) from Legendre coefficients" 1594 call wrtout(std_out,message,'COLL') 1595 ABI_ALLOCATE( jbes, (nleg)) 1596 gw_tmp_nd=czero 1597 1598 ! write(77,*) " TEST OF BESSEL S ROUTINES 0 0" 1599 1600 ! xx=0_dp 1601 ! ileg=0 1602 ! call sbf8(ileg+1,xx,jbes) 1603 ! write(77,*) "T0 A",jbes(ileg+1) 1604 ! call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx) 1605 ! write(77,*) "T0 B",jbes(ileg+1) 1606 ! write(77,*) "T0 C",bessel_jn(ileg,xx) 1607 1608 ! write(77,*) " TEST OF BESSEL S ROUTINES 1.5 0" 1609 1610 ! xx=1.5_dp 1611 ! ileg=0 1612 ! call sbf8(ileg+1,xx,jbes) 1613 ! write(77,*) "T1 A",jbes(ileg+1) 1614 ! call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx) 1615 ! write(77,*) "T1 B",jbes(ileg+1) 1616 ! write(77,*) "T1 C",bessel_jn(ileg,xx) 1617 1618 ! write(77,*) " TEST OF BESSEL S ROUTINES 1.5 1" 1619 1620 ! xx=1.5_dp 1621 ! ileg=1 1622 ! call sbf8(ileg+1,xx,jbes) 1623 ! write(77,*) "T2 A",jbes(ileg+1) 1624 ! call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx) 1625 ! write(77,*) "T2 B",jbes(ileg+1) 1626 ! write(77,*) "T2 C",bessel_jn(ileg,xx) 1627 1628 1629 do ifreq=1,paw_dmft%dmft_nwli 1630 xx=real(2*ifreq-1,kind=dp)*pi/two 1631 if(xx<=100_dp) call sbf8(nleg,xx,jbes) 1632 do ileg=1,nleg 1633 ! write(77,*) "A",ifreq,jbes(ileg),xx 1634 1635 if(xx>=99) call jbessel(jbes(ileg),besp,bespp,ileg-1,1,xx) 1636 ! write(77,*) "B",ifreq,jbes(ileg),xx 1637 1638 !write(77,*) "C",ifreq,jbes(ileg),xx 1639 1640 u_nl=sqrt(float(2*ileg-1))*(-1)**(ifreq-1)*cmplx(0_dp,one)**(ileg)*jbes(ileg) 1641 ! write(77,*) "----------",ileg,jbes(ileg), u_nl,gl_nd(ileg,1,1) 1642 1643 do iflavor=1,nflavor 1644 do iflavor1=1,nflavor 1645 gw_tmp_nd(ifreq,iflavor,iflavor1)= gw_tmp_nd(ifreq,iflavor,iflavor1) + & 1646 & u_nl*gl_nd(ileg,iflavor,iflavor1) 1647 end do 1648 end do 1649 1650 ! write(77,*) "------------------", gw_tmp_nd(ifreq,1,1) 1651 1652 end do 1653 ! write(77,*) "------------------ sum ", gw_tmp_nd(ifreq,1,1) 1654 end do 1655 ABI_DEALLOCATE( jbes ) 1656 call xmpi_barrier(paw_dmft%spacecomm) 1657 call flush_unit(std_out) 1658 end if 1659 gw_tmp_nd = gw_tmp_nd*Ha_eV 1660 1661 1662 1663 ! ========================================================================= 1664 1665 if ( leg_measure ) then !only if functionality is enabled 1666 ABI_DEALLOCATE(gl_nd) 1667 end if 1668 1669 else 1670 write(message,'(2a)') ch10," Can't launch TRIQS CTHYB solver because dmftqmc_l must be >= 2*dmft_nwli + 1" 1671 MSG_ERROR(message) 1672 end if 1673 1674 end if 1675 1676 ! ================================================================= 1677 ! END CTQMC run TRIQS 1678 ! ================================================================= 1679 1680 else if (testcode>=1) then 1681 1682 1683 ! ================================================================= 1684 ! CTQMC run for tests 1685 ! ================================================================= 1686 write(std_out,*) "nomega,dmftqmc_l",nomega,paw_dmft%dmftqmc_l 1687 call CtqmcInterface_run(hybrid,fw1(1:nomega,:),Gtau=gtmp,& 1688 & Gw=gw_tmp,E=green%ecorr_qmc(iatom),& 1689 & matU=umod,opt_levels=levels_ctqmc) 1690 1691 ! for non diagonal code 1692 ! call CtqmcInterface_run(hybrid,fw1_nd(1:nomega,:,:),Gtau=gtmp_nd,& 1693 !& Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),& 1694 !& Noise=Noise,matU=umod,opt_levels=levels_ctqmc,hybri_limit=hybri_limit) 1695 1696 end if 1697 1698 ! ================================================================= 1699 ! TEST 1700 ! If test of the code is activated, and testrot =1 rotate back green's function 1701 ! and stop the code. 1702 ! ================================================================= 1703 if(testcode==1) then 1704 1705 call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,& 1706 & levels_ctqmc,hybri_limit,nflavor,2,paw_dmft%temp,testrot,testcode,umod) 1707 1708 write(message,'(2a)') ch10,' testcode end of test calculation' 1709 MSG_ERROR(message) 1710 end if 1711 if(testcode==2) then 1712 write(message,'(2a)') ch10,' testcode 2 end of test calculation' 1713 MSG_ERROR(message) 1714 end if 1715 1716 ! ================================================================= 1717 ABI_DEALLOCATE(fw1) 1718 ABI_DEALLOCATE(fw1_nd) 1719 1720 !---------------------------------------- 1721 ! <DEBUG> 1722 !---------------------------------------- 1723 ! Construct UNIT 1724 if(paw_dmft%idmftloop < 10) then 1725 write(gtau_iter,'("0",i1)') paw_dmft%idmftloop 1726 elseif(paw_dmft%idmftloop >= 10 .and. paw_dmft%idmftloop < 100) then 1727 write(gtau_iter,'(i2)') paw_dmft%idmftloop 1728 else 1729 gtau_iter="xx" 1730 end if 1731 if(iatom < 10) then 1732 write(iatomnb,'("0",i1)') iatom 1733 elseif(iatom >= 10 .and. iatom < 100) then 1734 write(iatomnb,'(i2)') iatom 1735 else 1736 iatomnb='xx' 1737 end if 1738 1739 if(paw_dmft%myproc .eq. mod(nproc+1,nproc)) then 1740 ! < HACK > 1741 if(paw_dmft%dmft_solv>=6) then 1742 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /=0) then 1743 MSG_ERROR(message) 1744 end if 1745 do ifreq=1,paw_dmft%dmft_nwli 1746 write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),((gw_tmp_nd(ifreq,iflavor,iflavor)), iflavor=1, nflavor) 1747 end do 1748 close(unt) 1749 else 1750 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_"//gtau_iter//".dat", message, newunit=unt) /= 0) then 1751 MSG_ERROR(message) 1752 end if 1753 !call Ctqmc_printGreen(paw_dmft%hybrid(iatom)%hybrid,unt) !Problem here 1754 do itau=1,paw_dmft%dmftqmc_l 1755 write(unt,'(29f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,& 1756 (gtmp(itau,iflavor), iflavor=1, nflavor) 1757 end do 1758 write(unt,'(29f21.14)') 1/paw_dmft%temp, (-1_dp-gtmp(1,iflavor), iflavor=1, nflavor) 1759 close(unt) 1760 !open(unit=4243, file=trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_F_"//gtau_iter//".dat") 1761 !call BathOperator_printF(paw_dmft%hybrid(iatom)%hybrid%bath,4243) !Already comment here 1762 !close(4243) 1763 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /= 0) then 1764 MSG_ERROR(message) 1765 end if 1766 do ifreq=1,paw_dmft%dmft_nwlo 1767 write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq), & 1768 & (gw_tmp(ifreq,iflavor), iflavor=1, nflavor) 1769 end do 1770 close(unt) 1771 end if 1772 ! </ HACK > 1773 end if 1774 1775 !---------------------------------------- 1776 ! </DEBUG> 1777 !---------------------------------------- 1778 if(paw_dmft%dmft_solv>=6) then 1779 !Nothing just hybrid var problem 1780 else 1781 write(message,'(a,2x,a)') ch10,& 1782 & " == Destroy CTQMC" 1783 call wrtout(std_out,message,'COLL') 1784 call CtqmcInterface_finalize(hybrid) 1785 write(message,'(a,2x,a)') ch10,& 1786 & " == Destroy CTQMC done" 1787 call wrtout(std_out,message,'COLL') 1788 end if 1789 !ABI_DEALLOCATE(f_with_k) 1790 ABI_DEALLOCATE(hybri_limit) 1791 ABI_DEALLOCATE(levels_ctqmc_nd) 1792 ABI_DEALLOCATE(levels_ctqmc) 1793 1794 ! ___________________________________________________________________________________ 1795 ! 1796 ! FOURTH PART : USE OUTPUT OF CTQMC AND DO BACK ROTATION 1797 ! ___________________________________________________________________________________ 1798 ! 1799 do itau=1,paw_dmft%dmftqmc_l 1800 green%oper_tau(itau)%matlu(iatom)%mat(:,:,:,:,:)=czero 1801 end do 1802 green%occup_tau%matlu(iatom)%mat(nflavor:,:,:,:,:)=czero 1803 1804 do ifreq=1,paw_dmft%dmft_nwlo 1805 green%oper(ifreq)%matlu(iatom)%mat(:,:,:,:,:)=czero 1806 end do 1807 green%occup%matlu(iatom)%mat(:,:,:,:,:)=czero 1808 1809 ! built time and frequency green's function from output of CTQMC 1810 ! ================================================================= 1811 if(opt_nondiag==1) then 1812 do isppol=1,nsppol 1813 do ispinor1=1,nspinor 1814 do im1=1,tndim 1815 iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1) 1816 do ispinor2=1,nspinor 1817 do im2=1,tndim 1818 iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1) 1819 ! NNtodo here: nsppol=1 above, only, with symetrisation 1820 ! automatic 1821 ! iflavor1=im+(isppol-1)*tndim 1822 do itau=1,paw_dmft%dmftqmc_l 1823 green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 1824 & gtmp_nd(itau,iflavor1,iflavor2) 1825 if(nsppol==1.and.nspinor==1) then 1826 green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 1827 & (gtmp_nd(itau,iflavor1,iflavor2)+gtmp_nd(itau,iflavor1+tndim,iflavor2+tndim))/two 1828 end if 1829 ! NNtodo here: isppol above should be one and symetrized 1830 ! gtmp 1831 end do !itau 1832 ! ifreq2=0 1833 if(paw_dmft%dmft_solv<6.or.leg_measure) then 1834 do ifreq=1,paw_dmft%dmft_nwlo 1835 ! if(paw_dmft%select_log(ifreq)==1) then 1836 ! ifreq2=ifreq2+1 1837 green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 1838 & gw_tmp_nd(ifreq,iflavor1,iflavor2) 1839 if(nsppol==1.and.nspinor==1) then 1840 green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 1841 & (gw_tmp_nd(ifreq,iflavor1,iflavor2)+& 1842 & gw_tmp_nd(ifreq,iflavor1+tndim,iflavor2+tndim))/two 1843 end if 1844 ! NNtodo here: isppol above should be one and symetrized 1845 ! gw_tmp 1846 ! endif 1847 end do ! ifreq 1848 end if 1849 end do ! im2 1850 end do ! ispinor2 1851 end do ! im1 1852 end do ! ispinor 1853 end do ! isppol 1854 if(paw_dmft%dmft_solv>=6.and..not.leg_measure) then 1855 write(message,'(2a,i3,13x,a)') ch10,' === Direct Fourier Transform t->w of Weiss Field' 1856 call wrtout(std_out,message,'COLL') 1857 call fourier_green(cryst_struc,green,paw_dmft,& 1858 & pawang,opt_ksloc=2,opt_tw=1) 1859 end if 1860 else 1861 iflavor=0 1862 do isppol=1,nsppol 1863 do ispinor=1,nspinor 1864 do im=1,tndim 1865 ! NNtodo here: nsppol=1 above, only, with symetrisation 1866 ! automatic 1867 ! iflavor=im+(isppol-1)*tndim 1868 iflavor=iflavor+1 1869 do itau=1,paw_dmft%dmftqmc_l 1870 green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gtmp(itau,iflavor) 1871 if(nsppol==1.and.nspinor==1) then 1872 green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=& 1873 & (gtmp(itau,iflavor)+gtmp(itau,iflavor+tndim))/two 1874 end if 1875 ! NNtodo here: isppol above should be one and symetrized 1876 ! gtmp 1877 end do 1878 ! ifreq2=0 1879 do ifreq=1,paw_dmft%dmft_nwlo 1880 ! if(paw_dmft%select_log(ifreq)==1) then 1881 ! ifreq2=ifreq2+1 1882 green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gw_tmp(ifreq,iflavor) 1883 if(nsppol==1.and.nspinor==1) then 1884 green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=& 1885 & (gw_tmp(ifreq,iflavor)+gw_tmp(ifreq,iflavor+tndim))/two 1886 end if 1887 ! if opt_nondiag=0, then: 1888 ! As Green's function is diagonal, one suppress off diag terms in Weiss, if any. 1889 ! (If off diag are non zero in the density matrix and thus in the Green's function, 1890 ! there is a warning in checkreal_matlu above). 1891 do im1=1,tndim 1892 do ispinor1=1,nspinor 1893 if(im/=im1.or.ispinor/=ispinor1) then 1894 weiss%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=czero 1895 end if 1896 end do 1897 end do 1898 ! NNtodo here: isppol above should be one and symetrized 1899 ! gw_tmp 1900 ! endif 1901 end do 1902 end do 1903 end do 1904 end do 1905 end if 1906 if(paw_dmft%dmft_solv<6) ABI_DEALLOCATE(gw_tmp) 1907 ABI_DEALLOCATE(gw_tmp_nd) 1908 ABI_DEALLOCATE(gtmp) 1909 ABI_DEALLOCATE(gtmp_nd) 1910 if(nsppol==1.and.nspinor==1) then 1911 write(message,'(a,2x,a,f13.5)') ch10,& 1912 & " == nsppol==1 and nspden==1: Green functions from CTQMC have been symetrized over spin" 1913 call wrtout(std_out,message,'COLL') 1914 end if 1915 !write(message,'(i3,4x,2e21.14)') 5,weiss_for_rot%oper(1)%matlu(1)%mat(1,1,1,1,1) 1916 !call wrtout(std_out,message,'COLL') ! debug 1917 ! do im=1,tndim 1918 ! do itau=1,paw_dmft%dmftqmc_l 1919 ! gtt=(green%oper_tau(itau)%matlu(iatom)%mat(im,im,1,1,1)+& 1920 !& green%oper_tau(itau)%matlu(iatom)%mat(im,im,2,1,1))/two 1921 ! green%oper_tau(itau)%matlu(iatom)%mat(im,im,1,1,1)=gtt 1922 ! green%oper_tau(itau)%matlu(iatom)%mat(im,im,2,1,1)=gtt 1923 ! enddo 1924 ! enddo 1925 ! write(6,*)" SYMETRISATION OVER ISPPOL" 1926 !deallocate(correl_loc,f_with_k,gtmp,fw1,fw2) 1927 end if 1928 end do ! iatom 1929 1930 1931 if(paw_dmft%dmft_prgn==1) then 1932 call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=2) 1933 call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 1934 end if 1935 !write(message,'(i3,4x,2e21.14)') 6,weiss_for_rot%oper(1)%matlu(1)%mat(1,1,1,1,1) 1936 !call wrtout(std_out,message,'COLL') ! debug 1937 ! ================================================================= 1938 ! Inverse Weiss, then 1939 ! Copy Weiss_for_rot into weiss and rotate back weiss to the original basis 1940 ! ================================================================= 1941 1942 ! ABI_ALLOCATE(shift,(natom)) 1943 ! do ifreq=1,paw_dmft%dmft_nwlo 1944 ! ! First weiss_for_rot contains -G_0^-1+iw_n 1945 ! ! ------------------------------------------- 1946 ! ! Compute G_0^-1-iw_n 1947 ! ! -------------------- 1948 ! write(6,*) "1" 1949 ! if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone) 1950 ! 1951 ! 1952 ! write(6,*) "2" 1953 ! ! Compute G_0^-1 1954 ! ! -------------------- 1955 ! shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 1956 ! if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift,signe=1) 1957 ! 1958 ! write(6,*) "3" 1959 ! ! Compute G_0 1960 ! ! -------------------- 1961 ! call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1) 1962 ! ! No need to copy if weiss_for_rot is a pointer to weiss ... 1963 !! if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0) 1964 !! if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 1965 ! 1966 ! ! Compute G_0 in the original basis 1967 ! ! -------------------- 1968 ! call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 1969 ! end do 1970 ! ABI_DEALLOCATE(shift) 1971 1972 ! ================================================================= 1973 ! Here compute Self energy from Dyson and print it 1974 ! Warning : Weiss_for_rot is inversed inside dyson 1975 ! ================================================================= 1976 ! call initialize_self(self,paw_dmft) 1977 ! call dyson(green,paw_dmft,self,weiss_for_rot,opt_weissself=2) 1978 ! call rw_self(self,mpi_enreg,paw_dmft,prtopt=2,opt_rw=2,opt_char="diag") 1979 ! call destroy_self(self) 1980 !write(message,'(i3,4x,2e21.14)') 7,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1) 1981 !call wrtout(std_out,message,'COLL') ! debug 1982 1983 ! ================================================================= 1984 ! Rotate back green function to original basis (non-diagonal) 1985 ! (and Weiss for further use: might be useful if an back Fourier 1986 ! transformation is done). 1987 ! ================================================================= 1988 if(pawprtvol>=3) then 1989 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1990 & " == Print green function for small tau after CTQMC" ! debug 1991 call wrtout(std_out,message,'COLL') ! debug 1992 call print_matlu(green%oper_tau(1)%matlu,natom,1) ! debug 1993 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1994 & " == Print green function for small freq after CTQMC" ! debug 1995 call wrtout(std_out,message,'COLL') ! debug 1996 call print_matlu(green%oper(1)%matlu,natom,1) ! debug 1997 end if 1998 1999 2000 ! === Compute rotated Occupations in green%occup_tau 2001 call occup_green_tau(green) 2002 2003 if(pawprtvol>=3) then 2004 ! === Compute non rotated Occupations in green%occup_tau 2005 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the ctqmc basis" 2006 call wrtout(std_out,message,'COLL') 2007 call print_matlu(green%occup_tau%matlu,natom,1) 2008 end if 2009 2010 write(message,'(a,2x,a,f13.5)') ch10,& 2011 & " == Rotate Green function to original basis " 2012 call wrtout(std_out,message,'COLL') 2013 !write(message,'(i3,4x,2e21.14)') 8,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1) 2014 !call wrtout(std_out,message,'COLL') ! debug 2015 2016 ! Rotate oper_tau into Ylm basis and then Slm basis 2017 !------------------------------------------------------------- 2018 do itau=1,paw_dmft%dmftqmc_l 2019 if(opt_diag/=0) call rotate_matlu(green%oper_tau(itau)%matlu,eigvectmatlu,natom,3,0) 2020 if(useylm==1) call slm2ylm_matlu(green%oper_tau(itau)%matlu,natom,2,0) 2021 end do 2022 2023 ! Rotate occup_tau into Ylm basis and then Slm basis 2024 2025 ! Rotate occup_tau into Ylm basis and then Slm basis 2026 !------------------------------------------------------------- 2027 if(opt_diag/=0) call rotate_matlu(green%occup_tau%matlu,eigvectmatlu,natom,3,0) 2028 if(useylm==1) then 2029 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Ylm basis" 2030 call wrtout(std_out,message,'COLL') 2031 call print_matlu(green%occup_tau%matlu,natom,1) 2032 end if 2033 2034 if(useylm==1) call slm2ylm_matlu(green%occup_tau%matlu,natom,2,0) 2035 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Slm basis" 2036 call wrtout(std_out,message,'COLL') 2037 call print_matlu(green%occup_tau%matlu,natom,1) 2038 2039 ! Rotate Green's and Weiss functions into Ylm basis and then Slm basis 2040 !------------------------------------------------------------- 2041 do ifreq=1,paw_dmft%dmft_nwlo 2042 if(opt_diag/=0) call rotate_matlu(green%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 2043 if(useylm==1) call slm2ylm_matlu(green%oper(ifreq)%matlu,natom,2,0) 2044 if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 2045 if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0) 2046 end do 2047 !write(message,'(i3,4x,2e21.14)') 10,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1) 2048 !call wrtout(std_out,message,'COLL') ! debug 2049 2050 if(pawprtvol>=3) then 2051 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 2052 & " == Print green function for small time after rotation (in the original basis)" ! debug 2053 call wrtout(std_out,message,'COLL') ! debug 2054 call print_matlu(green%oper_tau(1)%matlu,natom,1) ! debug 2055 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 2056 & " == Print green function for small freq after rotation (in the original basis)" ! debug 2057 call wrtout(std_out,message,'COLL') ! debug 2058 call print_matlu(green%oper(1)%matlu,natom,1) ! debug 2059 !< HACK > 2060 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 2061 & " == Print diagonalized weiss_for_rot function after rotation for small freq in the ctqmc basis" ! debug 2062 call wrtout(std_out,message,'COLL') ! debug 2063 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) ! debug 2064 !</ HACK > 2065 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 2066 & " == Print weiss function for small freq in the original basis" ! debug 2067 call wrtout(std_out,message,'COLL') ! debug 2068 call print_matlu(weiss%oper(1)%matlu,natom,1) ! debug 2069 2070 do ifreq=1,paw_dmft%dmft_nwlo 2071 call sym_matlu(cryst_struc,weiss%oper(ifreq)%matlu,pawang) 2072 end do 2073 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 2074 & " == Print symetrized weiss function for small freq in the original basis" ! debug 2075 call wrtout(std_out,message,'COLL') ! debug 2076 call print_matlu(weiss%oper(1)%matlu,natom,1) ! debug 2077 end if 2078 2079 2080 ABI_DATATYPE_ALLOCATE(matlu1,(natom)) 2081 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 2082 call copy_matlu(green%occup_tau%matlu,matlu1,natom) 2083 call sym_matlu(cryst_struc,matlu1,pawang) 2084 2085 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the original basis" 2086 call wrtout(std_out,message,'COLL') 2087 call print_matlu(green%occup_tau%matlu,natom,1) 2088 2089 write(message,'(a,2x,a,f13.5)') ch10," == Symetrized occupations" 2090 call wrtout(std_out,message,'COLL') 2091 call print_matlu(matlu1,natom,1) 2092 2093 call diff_matlu("CTQMC Occup","CTQMC Occup symetrized",green%occup_tau%matlu,matlu1,natom,0,tol4,ierr) 2094 call destroy_matlu(matlu1,natom) 2095 ABI_DATATYPE_DEALLOCATE(matlu1) 2096 2097 ! ================================================================= 2098 ! Symetrise green function G(tau) and G(ifreq) to recover symetry 2099 ! artificially broken by QMC 2100 ! ================================================================= 2101 write(message,'(a,2x,a,f13.5)') ch10,& 2102 & " == Symetrise green function after QMC " 2103 call wrtout(std_out,message,'COLL') 2104 do itau=1,paw_dmft%dmftqmc_l 2105 call sym_matlu(cryst_struc,green%oper_tau(itau)%matlu,pawang) 2106 end do 2107 do ifreq=1,paw_dmft%dmft_nwlo 2108 call sym_matlu(cryst_struc,green%oper(ifreq)%matlu,pawang) 2109 end do 2110 if(pawprtvol>=3) then 2111 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 2112 & " == Print green function for small time after symetrisation" ! debug 2113 call wrtout(std_out,message,'COLL') ! debug 2114 call print_matlu(green%oper_tau(1)%matlu,natom,1) ! debug 2115 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 2116 & " == Print green function for small freq after symetrisation" ! debug 2117 call wrtout(std_out,message,'COLL') ! debug 2118 call print_matlu(green%oper(1)%matlu,natom,1) ! debug 2119 end if 2120 if(paw_dmft%dmft_prgn==1) then 2121 call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=2) 2122 call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 2123 end if 2124 2125 ! === Compute Occupations (Symetrized from oper_tau) 2126 call occup_green_tau(green) 2127 2128 2129 ! === Print occupations 2130 call printocc_green(green,6,paw_dmft,3) 2131 2132 call destroy_oper(energy_level) 2133 call destroy_matlu(dmat_diag,natom) 2134 ABI_DATATYPE_DEALLOCATE(dmat_diag) 2135 call destroy_matlu(identity,natom) 2136 ABI_DATATYPE_DEALLOCATE(identity) 2137 do iatom=1,cryst_struc%natom 2138 lpawu=paw_dmft%lpawu(iatom) 2139 if(lpawu/=-1) then 2140 do isppol=1,nsppol 2141 ABI_DEALLOCATE(eigvectmatlu(iatom,isppol)%value) 2142 !ABI_DEALLOCATE(udens_atoms(iatom)) 2143 end do 2144 ABI_DEALLOCATE(udens_atoms(iatom)%value) 2145 end if 2146 end do 2147 ABI_DATATYPE_DEALLOCATE(udens_atoms) 2148 ABI_DATATYPE_DEALLOCATE(eigvectmatlu) 2149 call destroy_green(weiss_for_rot) 2150 ! call destroy_green(gw_loc) 2151 ! call destroy_green(greenlda) 2152 2153 ! destroy limit of hybridization 2154 call destroy_matlu(hybri_coeff,paw_dmft%natom) 2155 ABI_DATATYPE_DEALLOCATE(hybri_coeff) 2156 2157 end subroutine qmc_prep_ctqmc