TABLE OF CONTENTS
- ABINIT/m_forctqmc
- m_forctqmc/ctqmc_calltriqs
- m_forctqmc/ctqmcoutput_printgreen
- m_forctqmc/ctqmcoutput_to_green
- m_forctqmc/qmc_prep_ctqmc
- m_forctqmc/testcode_ctqmc
- m_forctqmc/testcode_ctqmc_b
ABINIT/m_forctqmc [ Modules ]
NAME
m_forctqmc
FUNCTION
Prepare CTQMC and call CTQMC
COPYRIGHT
Copyright (C) 2006-2022 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 .
INPUTS
OUTPUT
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 MODULE m_forctqmc 27 28 use defs_basis 29 use m_errors 30 use m_xmpi 31 use m_abicore 32 use m_Ctqmc 33 use m_CtqmcInterface 34 use m_Ctqmcoffdiag 35 use m_CtqmcoffdiagInterface 36 use m_GreenHyb 37 use m_data4entropyDMFT 38 39 use m_pawang, only : pawang_type 40 use m_crystal, only : crystal_t 41 use m_green, only : green_type,occup_green_tau,print_green,printocc_green,spline_fct,copy_green,init_green,destroy_green,& 42 & int_fct,greendftcompute_green,fourier_green 43 use m_paw_dmft, only : paw_dmft_type 44 use m_hide_lapack, only : xginv 45 use m_oper, only : oper_type,destroy_oper,init_oper,inverse_oper 46 use m_self, only : self_type 47 use m_matlu, only : matlu_type,sym_matlu, print_matlu, & 48 & diag_matlu,init_matlu,destroy_matlu,rotate_matlu,checkdiag_matlu,checkreal_matlu, & 49 & copy_matlu, diff_matlu, slm2ylm_matlu, shift_matlu, prod_matlu,fac_matlu,& 50 & add_matlu,printplot_matlu,identity_matlu,zero_matlu 51 use m_hu, only : hu_type,rotatevee_hu,vee_ndim2tndim_hu_r,copy_hu,destroy_hu 52 use m_io_tools, only : flush_unit, open_file 53 use m_datafordmft, only : hybridization_asymptotic_coefficient,compute_levels 54 use m_special_funcs, only : sbf8 55 use m_paw_numeric, only : jbessel=>paw_jbessel 56 use m_paw_correlations, only : calc_vee 57 #ifdef HAVE_NETCDF 58 use netcdf !If calling TRIQS via python invocation, write a .nc file 59 #endif 60 61 implicit none 62 63 private 64 65 public :: qmc_prep_ctqmc 66 public :: testcode_ctqmc 67 public :: testcode_ctqmc_b 68 public :: ctqmcoutput_to_green 69 public :: ctqmcoutput_printgreen 70 public :: ctqmc_calltriqs
m_forctqmc/ctqmc_calltriqs [ Functions ]
[ Top ] [ m_forctqmc ] [ Functions ]
NAME
ctqmc_calltriqs
FUNCTION
Call TRIQS solver and perform calculation of Green's function using Legendre coefficients.
INPUTS
paw_dmft <type(paw_dmft_type)>= DMFT data structure cryst_struc <type(crystal_t)>=crystal structure data hu <type(hu_type)>= U interaction levels_ctqmc(nflavor) = atomic levels gw_tmp_nd(nb_of_frequency,nflavor,nflavor) = Green's fct in imag freq (with off diag terms) gtmp_nd(dmftqmc_l,nflavor,nflavor) = Green's fct in imag time (with off diag terms) fw1_nd(dmft_nwlo,nflavor,nflavor) = Hybridization fct in imag time (with off diag terms) leg_measure = logical, true is legendre measurement is activated iatom= index of atom
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
2432 subroutine ctqmc_calltriqs(paw_dmft,cryst_struc,hu,levels_ctqmc,gtmp_nd,gw_tmp_nd,fw1_nd,leg_measure,iatom) 2433 2434 #if defined HAVE_TRIQS_v2_0 || defined HAVE_TRIQS_v1_4 2435 use TRIQS_CTQMC !Triqs module 2436 #endif 2437 #if defined HAVE_PYTHON_INVOCATION 2438 use INVOKE_PYTHON 2439 #endif 2440 use, intrinsic :: iso_c_binding 2441 2442 !Arguments ------------------------------------ 2443 !scalars 2444 type(paw_dmft_type), intent(in) :: paw_dmft 2445 type(crystal_t),intent(in) :: cryst_struc 2446 type(hu_type), intent(in) :: hu(cryst_struc%ntypat) 2447 real(dp), allocatable, target, intent(inout) :: gtmp_nd(:,:,:) 2448 complex(dpc), allocatable, target, intent(inout) :: gw_tmp_nd(:,:,:) 2449 complex(dpc), allocatable, target, intent(in) :: fw1_nd(:,:,:) 2450 real(dp), allocatable, target, intent(inout) :: levels_ctqmc(:) 2451 logical(kind=1), intent(in) :: leg_measure 2452 integer, intent(in) :: iatom 2453 2454 !Local variables ------------------------------ 2455 complex(dpc), allocatable, target ::fw1_nd_tmp(:,:,:) 2456 complex(dpc), allocatable, target :: g_iw(:,:,:) 2457 real(dp), allocatable, target :: u_mat_ij(:,:) 2458 real(dp), allocatable, target :: u_mat_ijkl(:,:,:,:) 2459 real(dp), allocatable, target :: u_mat_ijkl_tmp(:,:,:,:) 2460 real(dp), allocatable, target :: gl_nd(:,:,:) 2461 type(c_ptr) :: levels_ptr, fw1_nd_ptr, u_mat_ij_ptr, u_mat_ijkl_ptr, g_iw_ptr, gtau_ptr, gl_ptr 2462 real(dp), allocatable :: jbes(:) 2463 character(len=500) :: message 2464 integer :: itau, ifreq, iflavor1 2465 integer :: iflavor2,iflavor,nflavor,iflavor3,itypat 2466 integer :: nfreq,ntau,nleg,ileg 2467 integer :: verbosity_solver ! min 0 -> max 3 2468 logical(kind=1) :: rot_inv = .false. 2469 logical(kind=1) :: hist = .false. 2470 logical(kind=1) :: wrt_files = .true. 2471 logical(kind=1) :: tot_not = .true. 2472 real(dp) :: beta,besp,bespp,xx 2473 complex(dpc) :: u_nl 2474 2475 !---------- 2476 !Variables for writing out the NETCDF file 2477 !---------- 2478 integer(kind=4) :: ncid 2479 integer(kind=4) :: dim_one_id, dim_nflavor_id, dim_nwlo_id, dim_nwli_id 2480 integer(kind=4) :: dim_qmc_l_id, dim_nleg_id 2481 integer(kind=4), dimension(2) :: dim_u_mat_ij_id 2482 integer(kind=4), dimension(3) :: dim_fw1_id, dim_g_iw_id, dim_gl_id, dim_gtau_id 2483 integer(kind=4), dimension(4) :: dim_u_mat_ijkl_id 2484 integer(kind=4) :: var_rot_inv_id, var_leg_measure_id, var_hist_id, var_wrt_files_id 2485 integer(kind=4) :: var_tot_not_id, var_n_orbitals_id, var_n_freq_id, var_n_tau_id, var_n_l_id, var_n_cycles_id 2486 integer(kind=4) :: var_cycle_length_id, var_ntherm_id, var_verbo_id, var_seed_id, var_beta_id 2487 integer(kind=4) :: var_levels_id, var_u_mat_ij_id, var_u_mat_ijkl_id, var_real_fw1_nd_id, var_imag_fw1_nd_id 2488 integer(kind=4) :: var_real_g_iw_id, var_imag_g_iw_id, var_gtau_id, var_gl_id, var_spacecomm_id 2489 2490 integer(kind=4) :: varid 2491 logical :: file_exists 2492 complex :: i 2493 character(len=100) :: filename 2494 2495 real(dp), allocatable, target :: new_re_g_iw(:,:,:), new_im_g_iw(:,:,:) 2496 real(dp), allocatable, target :: new_g_tau(:,:,:), new_gl(:,:,:) 2497 !---------- 2498 ! ************************************************************************ 2499 2500 ! fw1_nd: Hybridation 2501 ! levels_ctqmc: niveaux 2502 ! hu(itypat)%udens(:,:) : U_ij 2503 ! hu(itypat)%u(:,:,:,:) : uijkl 2504 ! temperature : paw_dmft%temp 2505 ! paw_dmft%dmftqmc_l: nombre de points en temps -1 2506 ! paw_dmft%dmftqmc_n: nombre de cycles 2507 ! ?? Quelles sorties: Les fonctions de Green 2508 ! frequence/temps/Legendre. 2509 ! Double occupations ?? <n_i n_j> 2510 ! test n_tau > 2*nfreq => ntau = 2*nfreq + 1 2511 ! for non diagonal code: 2512 ! call CtqmcInterface_run(hybrid,fw1_nd(1:paw_dmft%dmftqmc_l,:,:),Gtau=gtmp_nd,& 2513 !& Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),& 2514 !& Noise=Noise,matU=hu(itypat)%udens,opt_levels=levels_ctqmc,hybri_limit=hybri_limit) 2515 !Check choice of user to fix model bool var for the solver 2516 if (paw_dmft%dmft_solv==6) then 2517 rot_inv = .false. 2518 else !obviously paw_dmft%dmft_solv==7 with rot invariant terms 2519 rot_inv = .true. 2520 end if 2521 2522 nfreq = paw_dmft%dmft_nwli 2523 !paw_dmft%dmft_nwlo = paw_dmft%dmft_nwli !transparent for user 2524 ntau = paw_dmft%dmftqmc_l !(2*paw_dmft%dmftqmc_l)+1 !nfreq=paw_dmft%dmft_nwli 2525 nleg = paw_dmft%dmftctqmc_triqs_nleg 2526 nflavor=2*(2*paw_dmft%lpawu(iatom)+1) 2527 itypat=cryst_struc%typat(iatom) 2528 2529 2530 verbosity_solver = paw_dmft%prtvol 2531 beta = 1.0/(paw_dmft%temp*Ha_eV) 2532 2533 !Allocation in/output array phase: 2534 ABI_MALLOC(fw1_nd_tmp,(1:nflavor,1:nflavor,1:nfreq)) !column major 2535 ABI_MALLOC(g_iw,(1:nflavor,1:nflavor,1:nfreq)) !column major 2536 ABI_MALLOC(u_mat_ij,(1:nflavor,1:nflavor)) !column major 2537 ABI_MALLOC(u_mat_ijkl,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major 2538 ABI_MALLOC(u_mat_ijkl_tmp,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major 2539 2540 if ( leg_measure ) then !only if functionality is enabled 2541 ABI_MALLOC(gl_nd,(1:nleg,1:nflavor,1:nflavor)) !column major !nl = 30 by default 2542 end if 2543 2544 !Conversion datas Ha -> eV (some duplications for test...) 2545 !fw1_nd_tmp = fw1_nd(1:paw_dmft%dmftqmc_l,:,:) * Ha_eV !fw1_nd = fw1_nd * Ha_eV !Ok? 2546 2547 do iflavor=1,nflavor 2548 do iflavor1=1,nflavor 2549 do ifreq=1,nfreq 2550 fw1_nd_tmp(iflavor,iflavor1,ifreq) = fw1_nd(ifreq,iflavor,iflavor1) * Ha_eV 2551 ! WRITE(500,*) "[IN Fortran] F[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",fw1_nd(ifreq,iflavor,iflavor1) 2552 end do 2553 end do 2554 end do 2555 2556 !Report test 2557 ! WRITE(502,*) hu(itypat)%udens 2558 ! do ifreq=1,paw_dmft%dmftqmc_l 2559 ! write(501,*) ((fw1_nd(ifreq,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor) 2560 ! enddo 2561 !write(866,*)paw_dmft%dmft_nwlo,paw_dmft%dmftqmc_l 2562 !write(866,*) u_mat_ij 2563 ! do iflavor=1,nflavor+1 2564 ! do iflavor1=1,nflavor+1 2565 ! WRITE(502,*) "[OUT Fortran] U(i,j)[ l= ",iflavor," l_= ",iflavor1,"] = ",hu(itypat)%udens(iflavor,iflavor1) 2566 ! enddo 2567 ! enddo 2568 2569 ! if(paw_dmft%myproc==0) then 2570 ! do iflavor=1,nflavor 2571 ! do iflavor1=1,nflavor 2572 ! do iflavor2=1,nflavor 2573 ! do iflavor3=1,nflavor 2574 ! write(490,*), hu(itypat)%vee(iflavor,iflavor1,iflavor2,iflavor3) 2575 ! enddo 2576 ! enddo 2577 ! enddo 2578 ! enddo 2579 ! endif 2580 2581 ! if(paw_dmft%myproc==0) then 2582 ! do iflavor=1,nflavor 2583 ! do iflavor1=1,nflavor 2584 ! write(491,*), hu(itypat)%udens(iflavor,iflavor1) !(1,1,1,1) 2585 ! enddo 2586 ! enddo 2587 ! endif 2588 2589 ! do iflavor=1,nflavor 2590 ! do iflavor1=1,nflavor 2591 ! do iflavor2=1,nflavor 2592 ! do iflavor3=1,nflavor 2593 ! WRITE(552,*), hu(itypat)%vee!(iflavor,iflavor1,iflavor2,iflavor3) 2594 ! enddo 2595 ! enddo 2596 ! enddo 2597 ! enddo 2598 2599 call vee_ndim2tndim_hu_r(paw_dmft%lpawu(iatom),hu(itypat)%vee,u_mat_ijkl_tmp,1) 2600 do iflavor=1,nflavor 2601 do iflavor1=1,nflavor 2602 do iflavor2=1,nflavor 2603 do iflavor3=1,nflavor 2604 u_mat_ijkl(iflavor,iflavor1,iflavor2,iflavor3) = Ha_eV * u_mat_ijkl_tmp(iflavor,iflavor1,iflavor2,iflavor3) 2605 end do 2606 end do 2607 end do 2608 end do 2609 2610 !u_mat_ijkl = Ha_eV * reshape( u_mat_ijkl , [nflavor,nflavor,nflavor,nflavor] ) !column -> row major + conversion 2611 u_mat_ij = transpose( hu(itypat)%udens ) * Ha_eV !column -> row major + conversion 2612 levels_ctqmc = levels_ctqmc * Ha_eV 2613 2614 !Location array in memory for C++ pointer args to pass 2615 !---------------------------------------------------- 2616 g_iw_ptr = C_LOC( gw_tmp_nd ) !C_LOC( g_iw ) 2617 gtau_ptr = C_LOC( gtmp_nd ) !C_LOC( gtau ) 2618 gl_ptr = C_LOC( gl_nd ) 2619 fw1_nd_ptr = C_LOC( fw1_nd_tmp ) 2620 u_mat_ij_ptr = C_LOC( u_mat_ij ) 2621 u_mat_ijkl_ptr = C_LOC( u_mat_ijkl ) 2622 levels_ptr = C_LOC( levels_ctqmc ) 2623 2624 !Calling interfaced TRIQS solver subroutine from src/67_triqs_ext package 2625 if (paw_dmft%dmft_solv==9) then 2626 #ifndef HAVE_NETCDF 2627 write(message,'(2a)') ch10,' NETCDF requiered! ABINIT communicates with the python script through netcdf.' 2628 call wrtout(std_out,message,'COLL') 2629 ABI_ERROR(message) 2630 #else 2631 #ifndef HAVE_PYTHON_INVOCATION 2632 write(message,'(23a)') ch10,' Python invocation flag requiered! You need to install ABINIT with ',& 2633 'enable_python_invocation = yes" in your "configure.ac" file.' 2634 call wrtout(std_out,message,'COLL') 2635 ABI_ERROR(message) 2636 #else 2637 ! Creating the NETCDF file 2638 ! write(std_out, "(a)") trim(paw_dmft%filapp) 2639 write(filename, '(a, a)') trim(paw_dmft%filnamei), "_abinit_output_for_py.nc" 2640 write(std_out, '(3a)') ch10, " Creating NETCDF file: ", trim(filename) 2641 NCF_CHECK(nf90_create(filename, NF90_CLOBBER, ncid)) 2642 2643 ! Defining the dimensions of the variables to write in the NETCDF file 2644 NCF_CHECK(nf90_def_dim(ncid, "one", 1, dim_one_id)) 2645 NCF_CHECK(nf90_def_dim(ncid, "nflavor", nflavor, dim_nflavor_id)) 2646 NCF_CHECK(nf90_def_dim(ncid, "nwlo", paw_dmft%dmft_nwlo, dim_nwlo_id)) 2647 NCF_CHECK(nf90_def_dim(ncid, "nwli", paw_dmft%dmft_nwli, dim_nwli_id)) 2648 NCF_CHECK(nf90_def_dim(ncid, "qmc_l", paw_dmft%dmftqmc_l, dim_qmc_l_id)) 2649 NCF_CHECK(nf90_def_dim(ncid, "nleg", nleg, dim_nleg_id)) 2650 2651 dim_u_mat_ij_id = (/ dim_nflavor_id, dim_nflavor_id /) 2652 dim_u_mat_ijkl_id = (/ dim_nflavor_id, dim_nflavor_id, dim_nflavor_id, dim_nflavor_id /) 2653 dim_fw1_id = (/ dim_nflavor_id, dim_nflavor_id, dim_nwli_id /) 2654 dim_g_iw_id = (/ dim_nwli_id, dim_nflavor_id, dim_nflavor_id /) 2655 dim_gtau_id = (/ dim_qmc_l_id, dim_nflavor_id, dim_nflavor_id /) 2656 dim_gl_id = (/ dim_nleg_id, dim_nflavor_id, dim_nflavor_id /) 2657 2658 ! Defining the variables 2659 NCF_CHECK(nf90_def_var(ncid, "rot_inv", NF90_INT, dim_one_id, var_rot_inv_id)) 2660 NCF_CHECK(nf90_def_var(ncid, "leg_measure", NF90_INT, dim_one_id, var_leg_measure_id)) 2661 NCF_CHECK(nf90_def_var(ncid, "hist", NF90_INT, dim_one_id, var_hist_id)) 2662 NCF_CHECK(nf90_def_var(ncid, "wrt_files", NF90_INT, dim_one_id, var_wrt_files_id)) 2663 NCF_CHECK(nf90_def_var(ncid, "tot_not", NF90_INT, dim_one_id, var_tot_not_id)) 2664 NCF_CHECK(nf90_def_var(ncid, "n_orbitals", NF90_INT, dim_one_id, var_n_orbitals_id)) 2665 NCF_CHECK(nf90_def_var(ncid, "n_freq", NF90_INT, dim_one_id, var_n_freq_id)) 2666 NCF_CHECK(nf90_def_var(ncid, "n_tau", NF90_INT, dim_one_id, var_n_tau_id)) 2667 NCF_CHECK(nf90_def_var(ncid, "n_l", NF90_INT, dim_one_id, var_n_l_id)) 2668 NCF_CHECK(nf90_def_var(ncid, "n_cycles", NF90_INT, dim_one_id, var_n_cycles_id)) 2669 NCF_CHECK(nf90_def_var(ncid, "cycle_length", NF90_INT, dim_one_id, var_cycle_length_id)) 2670 NCF_CHECK(nf90_def_var(ncid, "ntherm", NF90_INT, dim_one_id, var_ntherm_id)) 2671 NCF_CHECK(nf90_def_var(ncid, "verbo", NF90_INT, dim_one_id, var_verbo_id)) 2672 NCF_CHECK(nf90_def_var(ncid, "seed", NF90_INT, dim_one_id, var_seed_id)) 2673 NCF_CHECK(nf90_def_var(ncid, "beta", NF90_FLOAT, dim_one_id, var_beta_id)) 2674 NCF_CHECK(nf90_def_var(ncid, "levels", NF90_DOUBLE, dim_nflavor_id, var_levels_id)) 2675 NCF_CHECK(nf90_def_var(ncid, "u_mat_ij", NF90_DOUBLE, dim_u_mat_ij_id, var_u_mat_ij_id)) 2676 NCF_CHECK(nf90_def_var(ncid, "u_mat_ijkl", NF90_DOUBLE, dim_u_mat_ijkl_id, var_u_mat_ijkl_id)) 2677 NCF_CHECK(nf90_def_var(ncid, "real_fw1_nd", NF90_DOUBLE, dim_fw1_id, var_real_fw1_nd_id)) 2678 NCF_CHECK(nf90_def_var(ncid, "imag_fw1_nd", NF90_DOUBLE, dim_fw1_id, var_imag_fw1_nd_id)) 2679 NCF_CHECK(nf90_def_var(ncid, "real_g_iw", NF90_DOUBLE, dim_g_iw_id, var_real_g_iw_id)) 2680 NCF_CHECK(nf90_def_var(ncid, "imag_g_iw", NF90_DOUBLE, dim_g_iw_id, var_imag_g_iw_id)) 2681 NCF_CHECK(nf90_def_var(ncid, "gtau", NF90_DOUBLE, dim_gtau_id, var_gtau_id)) 2682 NCF_CHECK(nf90_def_var(ncid, "gl", NF90_DOUBLE, dim_gl_id, var_gl_id)) 2683 NCF_CHECK(nf90_def_var(ncid, "spacecomm", NF90_INT, dim_one_id, var_spacecomm_id)) 2684 NCF_CHECK(nf90_enddef(ncid)) 2685 2686 ! Filling the variables with actual data 2687 if (rot_inv) then 2688 NCF_CHECK(nf90_put_var(ncid, var_rot_inv_id, 1)) 2689 else 2690 NCF_CHECK(nf90_put_var(ncid, var_rot_inv_id, 0)) 2691 end if 2692 if (leg_measure) then 2693 NCF_CHECK(nf90_put_var(ncid, var_leg_measure_id, 1)) 2694 else 2695 NCF_CHECK(nf90_put_var(ncid, var_leg_measure_id, 0)) 2696 end if 2697 if (hist) then 2698 NCF_CHECK(nf90_put_var(ncid, var_hist_id, 1)) 2699 else 2700 NCF_CHECK(nf90_put_var(ncid, var_hist_id, 0)) 2701 end if 2702 if (wrt_files) then 2703 NCF_CHECK(nf90_put_var(ncid, var_wrt_files_id, 1)) 2704 else 2705 NCF_CHECK(nf90_put_var(ncid, var_wrt_files_id, 0)) 2706 end if 2707 if (tot_not) then 2708 NCF_CHECK(nf90_put_var(ncid, var_tot_not_id, 1)) 2709 else 2710 NCF_CHECK(nf90_put_var(ncid, var_tot_not_id, 0)) 2711 end if 2712 NCF_CHECK(nf90_put_var(ncid, var_n_orbitals_id, nflavor)) 2713 NCF_CHECK(nf90_put_var(ncid, var_n_freq_id, nfreq)) 2714 NCF_CHECK(nf90_put_var(ncid, var_n_tau_id, ntau)) 2715 NCF_CHECK(nf90_put_var(ncid, var_n_l_id, nleg)) 2716 NCF_CHECK(nf90_put_var(ncid, var_n_cycles_id, int(paw_dmft%dmftqmc_n/paw_dmft%nproc))) 2717 NCF_CHECK(nf90_put_var(ncid, var_cycle_length_id, paw_dmft%dmftctqmc_meas*2*2*nflavor)) 2718 NCF_CHECK(nf90_put_var(ncid, var_ntherm_id, paw_dmft%dmftqmc_therm)) 2719 NCF_CHECK(nf90_put_var(ncid, var_verbo_id, verbosity_solver)) 2720 NCF_CHECK(nf90_put_var(ncid, var_seed_id, paw_dmft%dmftqmc_seed)) 2721 NCF_CHECK(nf90_put_var(ncid, var_beta_id, beta)) 2722 NCF_CHECK(nf90_put_var(ncid, var_levels_id, levels_ctqmc)) 2723 NCF_CHECK(nf90_put_var(ncid, var_u_mat_ij_id, u_mat_ij)) 2724 NCF_CHECK(nf90_put_var(ncid, var_u_mat_ijkl_id, u_mat_ijkl)) 2725 NCF_CHECK(nf90_put_var(ncid, var_real_fw1_nd_id, real(fw1_nd_tmp))) 2726 NCF_CHECK(nf90_put_var(ncid, var_imag_fw1_nd_id, aimag(fw1_nd_tmp))) 2727 NCF_CHECK(nf90_put_var(ncid, var_real_g_iw_id, real(gw_tmp_nd))) 2728 NCF_CHECK(nf90_put_var(ncid, var_imag_g_iw_id, aimag(gw_tmp_nd))) 2729 NCF_CHECK(nf90_put_var(ncid, var_gtau_id, gtmp_nd)) 2730 NCF_CHECK(nf90_put_var(ncid, var_gl_id, gl_nd)) 2731 NCF_CHECK(nf90_put_var(ncid, var_spacecomm_id, paw_dmft%spacecomm)) 2732 NCF_CHECK(nf90_close(ncid)) 2733 2734 write(std_out, '(4a)') ch10, " NETCDF file ", trim(filename), " written; Launching python invocation" 2735 2736 ! Invoking python to execute the script 2737 ! call Invoke_python_triqs (paw_dmft%myproc, trim(paw_dmft%filnamei)//c_null_char, paw_dmft%spacecomm) 2738 call Invoke_python_triqs (paw_dmft%myproc, trim(paw_dmft%filnamei)//c_null_char) 2739 call xmpi_barrier(paw_dmft%spacecomm) 2740 call flush_unit(std_out) 2741 2742 ! Allocating the fortran variables for the results 2743 ABI_MALLOC(new_re_g_iw,(nflavor,nflavor, paw_dmft%dmft_nwli)) 2744 ABI_MALLOC(new_im_g_iw,(nflavor,nflavor, paw_dmft%dmft_nwli)) 2745 ABI_MALLOC(new_g_tau,(nflavor,nflavor, paw_dmft%dmftqmc_l)) 2746 ABI_MALLOC(new_gl,(nflavor,nflavor, nleg)) 2747 i = (0, 1) 2748 2749 ! Check if file exists 2750 write(filename, '(a, a)') trim(paw_dmft%filnamei), "_py_output_for_abinit.nc" 2751 2752 INQUIRE(FILE=filename, EXIST=file_exists) 2753 if(.not. file_exists) then 2754 write(message,'(4a)') ch10,' Cannot find file ', trim(filename), '! Make sure the python script writes it with the right name and at the right place!' 2755 call wrtout(std_out,message,'COLL') 2756 ABI_ERROR(message) 2757 endif 2758 2759 write(std_out, '(3a)') ch10, " Reading NETCDF file ", trim(filename) 2760 2761 ! Opening the NETCDF file 2762 NCF_CHECK(nf90_open(filename, nf90_nowrite, ncid)) 2763 2764 ! Read from the file 2765 ! Re{G_iw} 2766 write(std_out, '(2a)') ch10, " -- Re[G(iw_n)]" 2767 NCF_CHECK(nf90_inq_varid(ncid, "re_g_iw", varid)) 2768 NCF_CHECK(nf90_get_var(ncid, varid, new_re_g_iw)) 2769 ! Im{G_iw} 2770 write(std_out, '(2a)') ch10, " -- Im[G(iw_n)]" 2771 NCF_CHECK(nf90_inq_varid(ncid, "im_g_iw", varid)) 2772 NCF_CHECK(nf90_get_var(ncid, varid, new_im_g_iw)) 2773 ! G_tau 2774 write(std_out, '(2a)') ch10, " -- G(tau)" 2775 NCF_CHECK(nf90_inq_varid(ncid, "g_tau", varid)) 2776 NCF_CHECK(nf90_get_var(ncid, varid, new_g_tau)) 2777 ! G_l 2778 write(std_out, '(2a)') ch10, " -- G_l" 2779 NCF_CHECK(nf90_inq_varid(ncid, "gl", varid)) 2780 NCF_CHECK(nf90_get_var(ncid, varid, new_gl)) 2781 2782 ! Assigning data 2783 do iflavor1=1, nflavor 2784 do iflavor2=1, nflavor 2785 do ifreq=1, paw_dmft%dmft_nwli 2786 gw_tmp_nd(ifreq, iflavor1, iflavor2) = new_re_g_iw(iflavor1, iflavor2, ifreq) & 2787 & + i*new_im_g_iw(iflavor1, iflavor2, ifreq) 2788 end do 2789 do itau=1, paw_dmft%dmftqmc_l 2790 gtmp_nd(itau, iflavor1, iflavor2) = new_g_tau(iflavor1, iflavor2, itau) 2791 end do 2792 do ileg=1, nleg 2793 gl_nd(ileg, iflavor1, iflavor2) = new_gl(iflavor1, iflavor2, ileg) 2794 end do 2795 end do 2796 end do 2797 2798 ! Deallocating 2799 ABI_FREE(new_re_g_iw) 2800 ABI_FREE(new_im_g_iw) 2801 ABI_FREE(new_g_tau) 2802 ABI_FREE(new_gl) 2803 #endif 2804 #endif 2805 elseif(paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then 2806 !Calling interfaced TRIQS solver subroutine from src/01_triqs_ext package 2807 !---------------------------------------------------- 2808 #if defined HAVE_TRIQS_v2_0 || defined HAVE_TRIQS_v1_4 2809 call Ctqmc_triqs_run ( rot_inv, leg_measure, hist, wrt_files, tot_not, & 2810 & nflavor, nfreq, ntau , nleg, int(paw_dmft%dmftqmc_n/paw_dmft%nproc), & 2811 & paw_dmft%dmftctqmc_meas*2*2*nflavor, paw_dmft%dmftqmc_therm, & 2812 & verbosity_solver, paw_dmft%dmftqmc_seed,beta, & 2813 & levels_ptr, u_mat_ij_ptr, u_mat_ijkl_ptr, fw1_nd_ptr, & 2814 !& g_iw_ptr, gtau_ptr, gl_ptr, paw_dmft%spacecomm ) 2815 & g_iw_ptr, gtau_ptr, gl_ptr, paw_dmft%myproc ) 2816 #endif 2817 endif 2818 2819 !WRITE(*,*) "Hello Debug" 2820 !call xmpi_barrier(paw_dmft%spacecomm) !Resynch all processus after calling Impurity solver from TRIQS 2821 2822 !Report output datas from TRIQS to Abinit 2823 !Interacting G(iw) 2824 do ifreq=1,nfreq 2825 do iflavor1=1,nflavor 2826 do iflavor=1,nflavor 2827 ! gw_tmp_nd(ifreq,iflavor,iflavor1) = g_iw(iflavor,iflavor1,ifreq) !* Ha_eV !because 1/ G0(eV) 2828 ! WRITE(503,*) "[OUT Fortran] G(iw)[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",gw_tmp_nd(ifreq,iflavor,iflavor1)!g_iw(iflavor,iflavor1,ifreq) 2829 end do 2830 end do 2831 end do 2832 2833 ! Convert in Ha 2834 gw_tmp_nd = gw_tmp_nd*Ha_eV 2835 2836 ! do iflavor1=1,nflavor 2837 ! do iflavor=1,nflavor 2838 ! 2839 ! WRITE(510,*) "[OUT Fortran] U[ l= ",iflavor," l_= ",iflavor1,"] = ",u_mat_ij(iflavor,iflavor1) 2840 ! enddo 2841 ! enddo 2842 2843 ! if(paw_dmft%myproc==0) write(6,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1) 2844 ! if(paw_dmft%myproc==1) write(6,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1) 2845 ! if(paw_dmft%myproc==0) write(621,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1) 2846 ! if(paw_dmft%myproc==1) write(622,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1) 2847 ! call flush_unit(621) 2848 ! call flush_unit(622) 2849 ! write(message,*) ch10, "essai",paw_dmft%myproc, paw_dmft%myproc,paw_dmft%dmftqmc_seed!gw_tmp_nd(2,1,1) 2850 ! call wrtout(555,message,'PERS',.true.) 2851 ! if(paw_dmft%myproc==0) write(499,*) "essai",paw_dmft%myproc, paw_dmft%dmftqmc_seed 2852 ! if(paw_dmft%myproc==1) write(498,*) "essai",paw_dmft%myproc,paw_dmft%dmftqmc_seed 2853 2854 !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)... 2855 ! do iflavor=1,nflavor 2856 ! do iflavor1=1,nflavor 2857 ! do itau=1,ntau 2858 ! 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 2859 !gtmp_nd(itau,iflavor,iflavor1) = gtau(iflavor,iflavor1,itau) 2860 ! endif 2861 ! if(paw_dmft%myproc==0) then 2862 ! 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) 2863 ! endif 2864 ! enddo 2865 ! enddo 2866 ! enddo 2867 2868 ! Write Legendre Polynoms G(L) for extrapolation of Interacting G(iw) by FT, only if leg_measure == TRUE 2869 ! ------------------------------------------------------------------------------------------- 2870 if (leg_measure) then 2871 do ileg=1,nleg 2872 WRITE(505,*) ileg,((gl_nd(ileg,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor) 2873 end do 2874 close(505) 2875 end if 2876 ! f(paw_dmft%myproc==0) then 2877 ! do itau=1,paw_dmft%dmftqmc_l 2878 ! write(490,*) ((gtmp_nd(itau,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor) 2879 ! enddo 2880 ! ndif 2881 ABI_FREE( fw1_nd_tmp ) 2882 ABI_FREE( g_iw ) 2883 ABI_FREE( u_mat_ijkl ) 2884 ABI_FREE( u_mat_ijkl_tmp ) 2885 ABI_FREE( u_mat_ij ) 2886 2887 2888 ! Compute Green's function in imaginary freq using Legendre coefficients 2889 ! ----------------------------------------------------------------------- 2890 if (leg_measure) then 2891 call xmpi_barrier(paw_dmft%spacecomm) 2892 call flush_unit(std_out) 2893 write(message,'(2a)') ch10," == Compute G(iw_n) from Legendre coefficients" 2894 call wrtout(std_out,message,'COLL') 2895 ABI_MALLOC( jbes, (nleg)) 2896 gw_tmp_nd=czero 2897 2898 ! write(77,*) " TEST OF BESSEL S ROUTINES 0 0" 2899 2900 ! xx=0_dp 2901 ! ileg=0 2902 ! call sbf8(ileg+1,xx,jbes) 2903 ! write(77,*) "T0 A",jbes(ileg+1) 2904 ! call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx) 2905 ! write(77,*) "T0 B",jbes(ileg+1) 2906 ! write(77,*) "T0 C",bessel_jn(ileg,xx) 2907 2908 ! write(77,*) " TEST OF BESSEL S ROUTINES 1.5 0" 2909 2910 ! xx=1.5_dp 2911 ! ileg=0 2912 ! call sbf8(ileg+1,xx,jbes) 2913 ! write(77,*) "T1 A",jbes(ileg+1) 2914 ! call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx) 2915 ! write(77,*) "T1 B",jbes(ileg+1) 2916 ! write(77,*) "T1 C",bessel_jn(ileg,xx) 2917 2918 ! write(77,*) " TEST OF BESSEL S ROUTINES 1.5 1" 2919 2920 ! xx=1.5_dp 2921 ! ileg=1 2922 ! call sbf8(ileg+1,xx,jbes) 2923 ! write(77,*) "T2 A",jbes(ileg+1) 2924 ! call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx) 2925 ! write(77,*) "T2 B",jbes(ileg+1) 2926 ! write(77,*) "T2 C",bessel_jn(ileg,xx) 2927 2928 2929 do ifreq=1,paw_dmft%dmft_nwli 2930 xx=real(2*ifreq-1,kind=dp)*pi/two 2931 if(xx<=100_dp) call sbf8(nleg,xx,jbes) 2932 do ileg=1,nleg 2933 ! write(77,*) "A",ifreq,jbes(ileg),xx 2934 2935 if(xx>=99) call jbessel(jbes(ileg),besp,bespp,ileg-1,1,xx) 2936 ! write(77,*) "B",ifreq,jbes(ileg),xx 2937 2938 !write(77,*) "C",ifreq,jbes(ileg),xx 2939 2940 u_nl=sqrt(float(2*ileg-1))*(-1)**(ifreq-1)*cmplx(0_dp,one)**(ileg)*jbes(ileg) 2941 write(77,*) "----------",ileg,jbes(ileg), u_nl,gl_nd(ileg,1,1) 2942 2943 do iflavor=1,nflavor 2944 do iflavor1=1,nflavor 2945 gw_tmp_nd(ifreq,iflavor,iflavor1)= gw_tmp_nd(ifreq,iflavor,iflavor1) + & 2946 & u_nl*gl_nd(ileg,iflavor,iflavor1) 2947 end do 2948 end do 2949 2950 ! write(77,*) "------------------", gw_tmp_nd(ifreq,1,1) 2951 2952 end do 2953 ! write(77,*) "------------------ sum ", gw_tmp_nd(ifreq,1,1) 2954 end do 2955 ABI_FREE( jbes ) 2956 call xmpi_barrier(paw_dmft%spacecomm) 2957 call flush_unit(std_out) 2958 end if 2959 gw_tmp_nd = gw_tmp_nd*Ha_eV 2960 2961 2962 if ( leg_measure ) then !only if functionality is enabled 2963 ABI_FREE(gl_nd) 2964 end if 2965 2966 2967 end subroutine ctqmc_calltriqs
m_forctqmc/ctqmcoutput_printgreen [ Functions ]
[ Top ] [ m_forctqmc ] [ Functions ]
NAME
ctqmcoutput_printgreen
FUNCTION
Print values of green function in files. Symetrize imaginary time Green's function in a peculiar case (dmft_solv=8 and natom=1). Should be moved later.
INPUTS
paw_dmft <type(paw_dmft_type)>= DMFT data structure gtmp_nd(dmftqmc_l,nflavor,nflavor) = Green's fct in imag time (with off diag terms) gw_tmp_nd(nb_of_frequency,nflavor,nflavor) = Green's fct in imag freq (with off diag terms) gtmp(dmftqmc_l,nflavor) = Green's fct in imag time (diag) gw_tmp(nb_of_frequency,nflavor+1) =Green's fct in imag freq (diag) iatom = atoms on which the calculation has been done
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
2261 subroutine ctqmcoutput_printgreen(paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom) 2262 2263 !Arguments ------------------------------------ 2264 !scalars 2265 type(paw_dmft_type), intent(in) :: paw_dmft 2266 real(dp), allocatable, intent(inout) :: gtmp_nd(:,:,:) 2267 complex(dpc), allocatable, intent(in) :: gw_tmp(:,:) 2268 complex(dpc), allocatable, intent(in) :: gw_tmp_nd(:,:,:) 2269 real(dp), allocatable, intent(in) :: gtmp(:,:) 2270 integer, intent(in) :: iatom 2271 2272 !Local variables ------------------------------ 2273 character(len=500) :: message 2274 integer :: ifreq, itau,iflavor1 2275 integer :: tndim,iflavor,nflavor 2276 character(len=2) :: gtau_iter,iatomnb 2277 integer :: unt 2278 ! ************************************************************************ 2279 tndim=2*paw_dmft%lpawu(iatom)+1 2280 nflavor=2*(tndim) 2281 !---------------------------------------- 2282 ! <DEBUG> 2283 !---------------------------------------- 2284 ! Construct UNIT 2285 if(paw_dmft%idmftloop < 10) then 2286 write(gtau_iter,'("0",i1)') paw_dmft%idmftloop 2287 elseif(paw_dmft%idmftloop >= 10 .and. paw_dmft%idmftloop < 100) then 2288 write(gtau_iter,'(i2)') paw_dmft%idmftloop 2289 else 2290 gtau_iter="xx" 2291 end if 2292 if(iatom < 10) then 2293 write(iatomnb,'("0",i1)') iatom 2294 elseif(iatom >= 10 .and. iatom < 100) then 2295 write(iatomnb,'(i2)') iatom 2296 else 2297 iatomnb='xx' 2298 end if 2299 2300 if(paw_dmft%myproc .eq. mod(paw_dmft%nproc+1,paw_dmft%nproc)) then 2301 ! < HACK > 2302 if(paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7) then 2303 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /=0) then 2304 ABI_ERROR(message) 2305 end if 2306 do ifreq=1,paw_dmft%dmft_nwli 2307 write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),((gw_tmp_nd(ifreq,iflavor,iflavor)), iflavor=1, nflavor) 2308 end do 2309 close(unt) 2310 else 2311 if(paw_dmft%dmft_solv==5) then 2312 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_"//gtau_iter//".dat", message, newunit=unt) /= 0) then 2313 ABI_ERROR(message) 2314 end if 2315 do itau=1,paw_dmft%dmftqmc_l 2316 write(unt,'(29f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,& 2317 (gtmp(itau,iflavor), iflavor=1, nflavor) 2318 end do 2319 write(unt,'(29f21.14)') 1/paw_dmft%temp, (-1_dp-gtmp(1,iflavor), iflavor=1, nflavor) 2320 close(unt) 2321 endif 2322 if(paw_dmft%dmft_solv==8) then 2323 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_offdiag_unsym_"//gtau_iter//".dat",& 2324 & message, newunit=unt) /= 0) then 2325 ABI_ERROR(message) 2326 end if 2327 do itau=1,paw_dmft%dmftqmc_l 2328 write(unt,'(196f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,& 2329 ((gtmp_nd(itau,iflavor,iflavor1), iflavor=1, nflavor),iflavor1=1, nflavor) 2330 end do 2331 close(unt) 2332 ! if(paw_dmft%natom==1) then ! If natom>1, it should be moved outside the loop over atoms 2333 ! ABI_MALLOC(matlu1,(paw_dmft%natom)) 2334 ! call init_matlu(paw_dmft%natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,matlu1) 2335 ! do itau=1,paw_dmft%dmftqmc_l 2336 ! do isppol=1,paw_dmft%nsppol 2337 ! do ispinor1=1,paw_dmft%nspinor 2338 ! do im1=1,tndim 2339 ! iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1) 2340 ! do ispinor2=1,paw_dmft%nspinor 2341 ! do im2=1,tndim 2342 ! iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1) 2343 ! matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 2344 !& gtmp_nd(itau,iflavor1,iflavor2) 2345 ! end do ! im2 2346 ! end do ! ispinor2 2347 ! end do ! im1 2348 ! end do ! ispinor 2349 ! end do ! isppol 2350 ! call rotate_matlu(matlu1,eigvectmatlu,paw_dmft%natom,3,0) 2351 ! call slm2ylm_matlu(matlu1,paw_dmft%natom,2,0) 2352 ! call sym_matlu(cryst_struc,matlu1,pawang,paw_dmft) 2353 ! call slm2ylm_matlu(matlu1,paw_dmft%natom,1,0) 2354 ! call rotate_matlu(matlu1,eigvectmatlu,paw_dmft%natom,3,1) 2355 ! do isppol=1,paw_dmft%nsppol 2356 ! do ispinor1=1,paw_dmft%nspinor 2357 ! do im1=1,tndim 2358 ! iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1) 2359 ! do ispinor2=1,paw_dmft%nspinor 2360 ! do im2=1,tndim 2361 ! iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1) 2362 ! gtmp_nd(itau,iflavor1,iflavor2)=& 2363 ! matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 2364 ! end do ! im2 2365 ! end do ! ispinor2 2366 ! end do ! im1 2367 ! end do ! ispinor 2368 ! end do ! isppol 2369 ! end do !itau 2370 ! call destroy_matlu(matlu1,paw_dmft%natom) 2371 ! ABI_FREE(matlu1) 2372 ! endif ! if natom=1 2373 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_offdiag_"//gtau_iter//".dat",& 2374 & message, newunit=unt) /= 0) then 2375 ABI_ERROR(message) 2376 end if 2377 do itau=1,paw_dmft%dmftqmc_l 2378 write(unt,'(196f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,& 2379 ((gtmp_nd(itau,iflavor,iflavor1), iflavor=1, nflavor),iflavor1=1, nflavor) 2380 end do 2381 close(unt) 2382 endif 2383 !open(unit=4243, file=trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_F_"//gtau_iter//".dat") 2384 !call BathOperator_printF(paw_dmft%hybrid(iatom)%hybrid%bath,4243) !Already comment here 2385 !close(4243) 2386 if(paw_dmft%dmft_solv==5) then 2387 if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /= 0) then 2388 ABI_ERROR(message) 2389 end if 2390 do ifreq=1,paw_dmft%dmft_nwlo 2391 write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq), & 2392 & (gw_tmp(ifreq,iflavor), iflavor=1, nflavor) 2393 end do 2394 endif 2395 close(unt) 2396 end if 2397 ! </ HACK > 2398 end if 2399 2400 2401 end subroutine ctqmcoutput_printgreen
m_forctqmc/ctqmcoutput_to_green [ Functions ]
[ Top ] [ m_forctqmc ] [ Functions ]
NAME
ctqmcoutput_to_green
FUNCTION
Put values of green function from ctqmc into green datatype Symetrize over spin if calculation is non magnetic
INPUTS
paw_dmft <type(paw_dmft_type)>= DMFT data structure gtmp_nd(dmftqmc_l,nflavor,nflavor) = Green's fct in imag time (with off diag terms) gw_tmp_nd(nb_of_frequency,nflavor,nflavor) = Green's fct in imag freq (with off diag terms) gtmp(dmftqmc_l,nflavor) = Green's fct in imag time (diag) gw_tmp(nb_of_frequency,nflavor+1) =Green's fct in imag freq (diag) iatom = atoms on which the calculation has been done leg_measure = logical, to Legendre Measurement or not (if done Green function is frequency is computed) opt_nondiag = integer, it activated, then
OUTPUT
green <type(green_type)>= green's function
SIDE EFFECTS
NOTES
SOURCE
2130 subroutine ctqmcoutput_to_green(green,paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom,leg_measure,opt_nondiag) 2131 2132 !Arguments ------------------------------------ 2133 !scalars 2134 type(paw_dmft_type), intent(in) :: paw_dmft 2135 type(green_type), intent(inout) :: green 2136 real(dp), allocatable, intent(in) :: gtmp_nd(:,:,:) 2137 complex(dpc), allocatable, intent(in) :: gw_tmp(:,:) 2138 complex(dpc), allocatable, intent(in) :: gw_tmp_nd(:,:,:) 2139 real(dp), allocatable, intent(in) :: gtmp(:,:) 2140 integer, intent(in) :: iatom,opt_nondiag 2141 logical(kind=1), intent(in) :: leg_measure 2142 character(len=500) :: message 2143 2144 !Local variables ------------------------------ 2145 integer :: ifreq, itau,im1,im2,isppol,ispinor1,ispinor2,iflavor1 2146 integer :: iflavor2,tndim,ispinor,iflavor,im,nflavor 2147 ! ************************************************************************ 2148 tndim=2*paw_dmft%lpawu(iatom)+1 2149 nflavor=2*(tndim) 2150 2151 do itau=1,paw_dmft%dmftqmc_l 2152 green%oper_tau(itau)%matlu(iatom)%mat(:,:,:,:,:)=czero 2153 end do 2154 green%occup_tau%matlu(iatom)%mat(nflavor:,:,:,:,:)=czero 2155 2156 do ifreq=1,paw_dmft%dmft_nwlo 2157 green%oper(ifreq)%matlu(iatom)%mat(:,:,:,:,:)=czero 2158 end do 2159 green%occup%matlu(iatom)%mat(:,:,:,:,:)=czero 2160 2161 ! built time and frequency green's function from output of CTQMC 2162 ! ================================================================= 2163 if(opt_nondiag==1) then 2164 do isppol=1,paw_dmft%nsppol 2165 do ispinor1=1,paw_dmft%nspinor 2166 do im1=1,tndim 2167 iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1) 2168 do ispinor2=1,paw_dmft%nspinor 2169 do im2=1,tndim 2170 iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1) 2171 do itau=1,paw_dmft%dmftqmc_l 2172 green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 2173 & gtmp_nd(itau,iflavor1,iflavor2) 2174 ! symetrize over spin if nsppol=nspinor=1 2175 if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then 2176 green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 2177 & (gtmp_nd(itau,iflavor1,iflavor2)+gtmp_nd(itau,iflavor1+tndim,iflavor2+tndim))/two 2178 end if 2179 end do !itau 2180 if(paw_dmft%dmft_solv<6.or.leg_measure) then 2181 do ifreq=1,paw_dmft%dmft_nwlo 2182 green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 2183 & gw_tmp_nd(ifreq,iflavor1,iflavor2) 2184 ! symetrize over spin if nsppol=nspinor=1 2185 if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then 2186 green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=& 2187 & (gw_tmp_nd(ifreq,iflavor1,iflavor2)+& 2188 & gw_tmp_nd(ifreq,iflavor1+tndim,iflavor2+tndim))/two 2189 end if 2190 end do ! ifreq 2191 end if 2192 end do ! im2 2193 end do ! ispinor2 2194 end do ! im1 2195 end do ! ispinor 2196 end do ! isppol 2197 else 2198 iflavor=0 2199 do isppol=1,paw_dmft%nsppol 2200 do ispinor=1,paw_dmft%nspinor 2201 do im=1,tndim 2202 iflavor=iflavor+1 2203 do itau=1,paw_dmft%dmftqmc_l 2204 green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gtmp(itau,iflavor) 2205 ! symetrize over spin if nsppol=paw_dmft%nspinor=1 2206 if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then 2207 green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=& 2208 & (gtmp(itau,iflavor)+gtmp(itau,iflavor+tndim))/two 2209 end if 2210 end do 2211 ! ifreq2=0 2212 do ifreq=1,paw_dmft%dmft_nwlo 2213 ! if(paw_dmft%select_log(ifreq)==1) then 2214 ! ifreq2=ifreq2+1 2215 green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gw_tmp(ifreq,iflavor) 2216 ! symetrize over spin if nsppol=paw_dmft%nspinor=1 2217 if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then 2218 green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=& 2219 & (gw_tmp(ifreq,iflavor)+gw_tmp(ifreq,iflavor+tndim))/two 2220 end if 2221 end do 2222 end do 2223 end do 2224 end do 2225 end if 2226 if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then 2227 write(message,'(a,2x,a,f13.5)') ch10,& 2228 & " == nsppol==1 and nspden==1: Green functions from CTQMC have been symetrized over spin" 2229 call wrtout(std_out,message,'COLL') 2230 end if 2231 2232 end subroutine ctqmcoutput_to_green
m_forctqmc/qmc_prep_ctqmc [ Functions ]
[ Top ] [ m_forctqmc ] [ Functions ]
NAME
qmc_prep_ctqmc
FUNCTION
Prepare and call the qmc subroutines
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
SOURCE
97 subroutine qmc_prep_ctqmc(cryst_struc,green,self,hu,paw_dmft,pawang,pawprtvol,weiss) 98 99 !Arguments ------------------------------------ 100 !scalars 101 ! type(pawang_type), intent(in) :: pawang 102 type(crystal_t),intent(in) :: cryst_struc 103 type(green_type), intent(inout) :: green ! MGNAG: This fix the problem with v7[27:29] on nag@petrus 104 type(hu_type), intent(in) :: hu(cryst_struc%ntypat) 105 type(paw_dmft_type), intent(inout) :: paw_dmft 106 type(pawang_type), intent(in) :: pawang 107 integer, intent(in) :: pawprtvol 108 type(green_type), intent(inout) :: weiss 109 type(self_type), intent(in) :: self 110 111 !Local variables ------------------------------ 112 character(len=500) :: message 113 integer :: iatom,ierr,if1,if2,iflavor1,iflavor2,ifreq,im1,ispinor,ispinor1,isppol,itau,itypat,im2,ispinor2 114 integer :: lpawu,master,mbandc,natom,nflavor,nkpt,nspinor,nsppol,nsppol_imp,tndim,ispa,ispb,ima,imb 115 integer :: nproc,opt_diag,opt_nondiag,testcode,testrot,dmft_nwlo,opt_fk,useylm,nomega,opt_rot 116 integer :: ier,rot_type_vee 117 real(dp) :: f4of2_sla,f6of2_sla 118 complex(dpc) :: omega_current,integral(2,2) 119 real(dp) :: doccsum,noise,omega,EE 120 logical :: nondiaglevels 121 ! arrays 122 real(dp), allocatable :: docc(:,:) 123 real(dp), allocatable, target :: gtmp(:,:), levels_ctqmc(:) !modif 124 complex(dpc), allocatable :: levels_ctqmc_nd(:,:) 125 complex(dpc), allocatable :: hybri_limit(:,:) 126 real(dp), allocatable, target :: gtmp_nd(:,:,:) 127 real(dp) :: umod(2,2) 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(coeff2_type), allocatable :: udens_atoms_for_s(:) 136 ! Type ----------------------------------------- 137 type(coeff2c_type), allocatable :: eigvectmatlu(:,:) 138 type(green_type) :: weiss_for_rot 139 type(matlu_type), allocatable :: dmat_diag(:) 140 type(matlu_type), allocatable :: matlu1(:) 141 type(matlu_type), allocatable :: matlu2(:) 142 type(matlu_type), allocatable :: matlu3(:) 143 type(matlu_type), allocatable :: matlu4(:) 144 type(matlu_type), allocatable :: identity(:) 145 type(matlu_type), allocatable :: level_diag(:) 146 type(oper_type) :: energy_level 147 type(hu_type), allocatable :: hu_for_s(:) 148 !type(self_type) :: self 149 ! type(green_type) :: gw_loc 150 type(CtqmcInterface) :: hybrid !!! WARNING THIS IS A BACKUP PLAN 151 type(CtqmcoffdiagInterface) :: hybridoffdiag !!! WARNING THIS IS A BACKUP PLAN 152 type(green_type) :: greendft 153 type(matlu_type), allocatable :: hybri_coeff(:) 154 integer :: unt,unt2 155 ! Var added to the code for TRIQS_CTQMC test and default value ----------------------------------------------------------- 156 logical(kind=1) :: leg_measure = .true. 157 ! ************************************************************************ 158 mbandc=paw_dmft%mbandc 159 nkpt=paw_dmft%nkpt 160 nsppol=paw_dmft%nsppol 161 natom=paw_dmft%natom 162 nspinor=paw_dmft%nspinor 163 greendft%whichgreen="DFT" 164 165 call init_green(weiss_for_rot,paw_dmft,opt_oper_ksloc=2) 166 ! weiss_for_rot=>weiss 167 ! call init_green(gw_loc,paw_dmft) 168 call copy_green(weiss,weiss_for_rot,opt_tw=2) 169 !======================================================================= 170 !== Use one QMC solver =============================================== 171 !======================================================================= 172 write(message,'(2a)') ch10,' === CT-QMC solver === ' 173 call wrtout(std_out,message,'COLL') 174 175 ! Initialise for compiler 176 omega_current=czero 177 178 ! Initialise nproc 179 nproc=paw_dmft%nproc 180 181 ! ====================================== 182 ! Allocations: diagonalization and eigenvectors 183 ! ====================================== 184 ABI_MALLOC(udens_atoms,(natom)) 185 ABI_MALLOC(eigvectmatlu,(natom,nsppol)) 186 if(paw_dmft%ientropy==1) then 187 ABI_MALLOC(udens_atoms_for_s,(natom)) 188 endif 189 ABI_MALLOC(dmat_diag,(natom)) 190 ABI_MALLOC(identity,(natom)) 191 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,dmat_diag) 192 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,identity) 193 call identity_matlu(identity,natom) 194 do iatom=1,cryst_struc%natom 195 lpawu=paw_dmft%lpawu(iatom) 196 if(lpawu/=-1) then 197 tndim=nspinor*(2*lpawu+1) 198 do isppol=1,nsppol 199 ABI_MALLOC(eigvectmatlu(iatom,isppol)%value,(tndim,tndim)) 200 end do 201 ABI_MALLOC(udens_atoms(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1))) 202 if(paw_dmft%ientropy==1) then 203 ABI_MALLOC(udens_atoms_for_s(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1))) 204 endif 205 dmat_diag(iatom)%mat=czero 206 end if 207 end do 208 209 ! ___________________________________________________________________________________ 210 ! 211 ! FIRST PART: DIAGONALISATION AND ROTATIONS. 212 ! ___________________________________________________________________________________ 213 ! 214 215 ! ================================================================= 216 ! Choose to diagonalize and how to do it 217 ! ================================================================= 218 219 ! ================================================================= 220 ! Impose diago of density matrix 221 ! ================================================================= 222 223 ! ================================================================= 224 ! Impose diago of levels and Ylm basis if opt_nondiag=1 225 ! ================================================================= 226 ! opt_diag=1 ! 1: diago the levels (The best choice). 227 ! opt_diag=2 ! 2: diago density matrix (can be used for historical reasons) 228 229 ! Need in the general case of two input variable for opt_diag and 230 ! opt_nondiag! 231 ! The default value of opt_diag should be 2 for historical reasons (or 232 ! we decide to change the automatic tests) 233 ! opt_nondiag should be 0 by default 234 opt_diag = 1 235 if(paw_dmft%dmft_solv>=6) then 236 opt_nondiag = 1 ! Use cthyb in triqs or ctqmc in abinit with offdiag terms in F 237 else 238 opt_nondiag = 0 ! use fast ctqmc in ABINIT without off diagonal terms in F 239 end if 240 241 useylm=0 242 if(nspinor==2) then 243 useylm=1 ! to avoid complex G(tau) 244 end if 245 246 !write(6,*) "nspinor,useylm",nspinor,useylm 247 if(useylm==0) then 248 write(std_out,*) " Slm basis is used (before rotation)" 249 rot_type_vee=1 ! for rotatevee_hu 250 else if(useylm==1) then 251 write(std_out,*) " Ylm basis is used (before rotation)" 252 rot_type_vee=4 ! for rotatevee_hu 253 end if 254 255 256 ! if(useylm==1.and.opt_diag/=1) ABI_ERROR("useylm==1 and opt_diag/=0 is not possible") 257 if(hu(1)%jpawu_zero.and.nsppol==2) nsppol_imp=2 ! J=0 and nsppol=2 258 if(.not.hu(1)%jpawu_zero.or.nsppol/=2) nsppol_imp=1 ! J/=0 ou nsppol=1 259 ! ================================================================= 260 ! Compute DFT Green's function to compare to weiss_for_rot (check) 261 ! ================================================================= 262 ! call init_green(greendft,paw_dmft,opt_oper_ksloc=3) 263 ! call greendftcompute_green(cryst_struc,greendft,pawang,paw_dmft) 264 !! call copy_green(greendft,weiss_for_rot,2) 265 266 ! ================================================================= 267 ! Compute atomic levels 268 ! ================================================================= 269 call init_oper(paw_dmft,energy_level,opt_ksloc=3) 270 271 ! Compute atomic levels in Slm basis 272 ! ---------------------------------- 273 call compute_levels(cryst_struc,energy_level,self%hdc,pawang,paw_dmft,nondiag=nondiaglevels) 274 275 ! If levels are not diagonal, then diagonalize it (according to 276 ! dmftctqmc_basis) 277 ! ------------------------------------------------ 278 if(paw_dmft%dmftctqmc_basis==1) then 279 if(nondiaglevels.or.useylm==1) then 280 opt_diag=1 281 write(message,'(3a)') ch10, " == Hamiltonian in local basis is non diagonal: diagonalise it",ch10 282 else 283 opt_diag=0 284 write(message,'(5a)') ch10, " == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 & 285 & ," CTQMC will use this basis",ch10 286 end if 287 else if (paw_dmft%dmftctqmc_basis==2) then 288 if(nondiaglevels.or.useylm==1) then 289 write(message,'(7a)') ch10, " == Hamiltonian in local basis is non diagonal",ch10, & 290 & " == According to dmftctqmc_basis: diagonalise density matrix",ch10, & 291 & " == Warning : Check that the Hamiltonian is diagonal !",ch10 292 opt_diag=2 293 else 294 write(message,'(5a)') ch10, " == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 & 295 & ," CTQMC will use this basis",ch10 296 opt_diag=0 297 end if 298 else if (paw_dmft%dmftctqmc_basis==0) then 299 if(nondiaglevels) then 300 write(message,'(4a)') ch10, " == Hamiltonian in local basis is non diagonal",ch10, & 301 & " == According to dmftctqmc_basis: keep this non diagonal basis for the calculation" 302 else 303 write(message,'(5a)') ch10, " == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 & 304 & ," CTQMC will use this basis",ch10 305 end if 306 opt_diag=0 307 end if 308 call wrtout(std_out,message,'COLL') 309 if(opt_diag==1) then 310 write(std_out,*) " == The atomic levels are diagonalized" 311 else if(opt_diag==2) then 312 write(std_out,*) " == The correlated occupation matrix is diagonalized" 313 end if 314 315 ! ================================================================= 316 ! Now, check if diagonalisation is necessary 317 ! ================================================================= 318 319 320 ! ================================================================= 321 ! First rotate to Ylm basis the atomic levels 322 ! ================================================================= 323 324 if(useylm==1) then 325 326 ! Rotate from Slm to Ylm the atomic levels 327 ! ---------------------------------------- 328 call slm2ylm_matlu(energy_level%matlu,natom,1,pawprtvol) 329 330 ! Print atomic energy levels in Ylm basis 331 ! -------------------------------- 332 if(pawprtvol>=3) then 333 write(message,'(a,a)') ch10, " == Print Energy levels in Ylm basis" 334 call wrtout(std_out,message,'COLL') 335 call print_matlu(energy_level%matlu,natom,1) 336 end if 337 338 end if ! useylm 339 340 ! =========================================================================================== 341 ! Start for diagonalization of levels/density matrix according to opt_diag 342 ! =========================================================================================== 343 !opt_rot=2 ! do it one time before CTQMC 344 opt_rot=1 ! do all the rotations successively on all different quantities. 345 if(opt_diag==1.or.opt_diag==0) then 346 347 348 if(opt_diag==1) then 349 ! ================================================================= 350 ! Diagonalize atomic levels 351 ! ================================================================= 352 ABI_MALLOC(level_diag,(natom)) 353 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,level_diag) 354 355 ! Diagonalise atomic levels (opt_real is necessary, because 356 ! rotation must be real in order for the occupations and Green's 357 ! function to be real) 358 ! --------------------------------------------------------------- 359 call diag_matlu(energy_level%matlu,level_diag,natom,& 360 & prtopt=pawprtvol,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,optreal=1,& 361 & test=paw_dmft%dmft_solv) ! temporary: test should be extended to all cases. 362 363 ! call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1) 364 ! write(message,'(a,2x,a,f13.5)') ch10,& 365 !& " == Print first Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie 366 ! call wrtout(std_out,message,'COLL') 367 ! call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1) 368 369 if(opt_rot==1) call copy_matlu(level_diag,energy_level%matlu,natom) 370 371 372 call destroy_matlu(level_diag,natom) 373 ABI_FREE(level_diag) 374 375 ! Print diagonalized levels 376 ! -------------------------- 377 if(pawprtvol>=3) then 378 write(message,'(a,2x,a,f13.5)') ch10,& 379 & " == Print Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie 380 call wrtout(std_out,message,'COLL') 381 call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1) 382 else 383 write(message,'(a,2x,a,f13.5)') ch10,& 384 & " == Energy levels Diagonalized for Fermi Level=",paw_dmft%fermie 385 call wrtout(std_out,message,'COLL') 386 end if 387 388 call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee) 389 390 391 else if (opt_diag==0) then 392 do iatom=1,cryst_struc%natom 393 lpawu=paw_dmft%lpawu(iatom) 394 itypat=cryst_struc%typat(iatom) 395 if(lpawu/=-1) then 396 ! write(6,*) size(udens_atoms(iatom)%value) 397 ! write(6,*) size(hu(itypat)%udens) 398 ! write(6,*) udens_atoms(iatom)%value 399 ! write(6,*) hu(itypat)%udens 400 udens_atoms(iatom)%value=hu(itypat)%udens 401 end if 402 end do 403 end if 404 ! call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms) 405 406 else if(opt_diag==2) then 407 ! ================================================================= 408 ! Diagonalizes density matrix and keep eigenvectors in eigvectmatlu 409 ! ================================================================= 410 411 ! Print density matrix before diagonalization 412 ! ------------------------------------------- 413 if(pawprtvol>=3) then 414 write(message,'(a,2x,a)') ch10, " == Density Matrix before diagonalisation =" 415 call wrtout(std_out,message,'COLL') 416 !MGNAG: This call is wrong if green has intent(out), now we use intent(inout) 417 call print_matlu(green%occup%matlu,natom,1) 418 end if 419 420 !! checkstop: we can have two different diagonalisation basis for the up and dn 421 !! but one use the same basis, unless the error is really to large(>0.1) 422 423 ! Diagonalize density matrix 424 ! --------------------------- 425 call diag_matlu(green%occup%matlu,dmat_diag,natom,& 426 & prtopt=4,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,checkstop=.false.) 427 428 ! Print diagonalized density matrix 429 ! ---------------------------------- 430 if(pawprtvol>=3) then 431 write(message,'(a,2x,a)') ch10,& 432 & " == Diagonalized Density Matrix in the basis used for QMC =" 433 call wrtout(std_out,message,'COLL') 434 call print_matlu(dmat_diag,natom,1) 435 436 !write(message,'(2a,i3,13x,a)') ch10,' == Rotation of interaction matrix ==' 437 !call wrtout(std_out,message,'COLL') 438 end if 439 440 !if (.not.hu(1)%jpawu_zero) & 441 !ABI_WARNING("In qmc_prep_ctqmc J/=0 and rotation matrix not rotated") 442 ! Rotate interaction. 443 ! call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms) 444 call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee) 445 446 447 end if 448 ! =========================================================================================== 449 ! END Of diagonalization 450 ! =========================================================================================== 451 if(paw_dmft%ientropy==1) then 452 ABI_MALLOC(hu_for_s,(cryst_struc%ntypat)) 453 ! Usefull to compute interaction energy for U=1 J=J/U when U=0. 454 call copy_hu(cryst_struc%ntypat,hu,hu_for_s) 455 f4of2_sla=-1_dp 456 f6of2_sla=-1_dp 457 do itypat=1,cryst_struc%ntypat 458 call calc_vee(f4of2_sla,f6of2_sla,paw_dmft%j_for_s/paw_dmft%u_for_s,hu_for_s(itypat)%lpawu,pawang,one,hu_for_s(itypat)%vee) 459 enddo 460 call rotatevee_hu(cryst_struc,hu_for_s,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms_for_s,rot_type_vee) 461 call destroy_hu(hu_for_s,cryst_struc%ntypat,paw_dmft%dmftqmc_t2g,paw_dmft%dmftqmc_x2my2d) 462 ! udens_atoms_for_s will be used later. 463 ABI_FREE(hu_for_s) 464 endif 465 466 467 call flush_unit(std_out) 468 469 ! =========================================================================================== 470 ! Broadcast matrix of rotation from processor 0 to the other 471 ! In case of degenerate levels, severals rotations are possible. Here we 472 ! choose the rotation of proc 0. It is arbitrary. 473 ! =========================================================================================== 474 do iatom=1,cryst_struc%natom 475 lpawu=paw_dmft%lpawu(iatom) 476 if(lpawu/=-1) then 477 tndim=nspinor*(2*lpawu+1) 478 do isppol=1,nsppol 479 call xmpi_bcast(eigvectmatlu(iatom,isppol)%value,0,paw_dmft%spacecomm,ier) 480 end do 481 end if 482 end do 483 484 485 !unitnb=300000+paw_dmft%myproc 486 !call int2char4(paw_dmft%myproc,tag_proc) 487 !tmpfil = 'eigvectmatluaftermpi'//tag_proc 488 !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted') 489 !do iflavor1=1,14 490 ! do iflavor2=1,14 491 ! write(unitnb,*) iflavor1,iflavor2,eigvectmatlu(1,1)%value(iflavor1,iflavor2) 492 ! enddo 493 !enddo 494 495 ! =========================================================================================== 496 ! Now rotate various quantities in the new basis 497 ! =========================================================================================== 498 499 !======================================================= 500 ! Allocate, Compute, and Rotate atomic levels for CTQMC 501 !======================================================= 502 503 ! If levels not rotated, rotate them 504 ! ----------------------------------- 505 if(opt_diag==2.and.opt_rot==1) call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1) 506 507 ! Print atomic levels 508 ! ------------------- 509 if(pawprtvol>=3) then 510 write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels after rotation" 511 call wrtout(std_out,message,'COLL') 512 call print_matlu(energy_level%matlu,natom,1) 513 else 514 write(message,'(a,2x,a,f13.5)') ch10," == CT-QMC Energy levels rotated" 515 call wrtout(std_out,message,'COLL') 516 end if 517 518 !==================================================================== 519 ! If levels were diagonalized before, then rotate density matrix for 520 ! information. 521 !==================================================================== 522 if(opt_diag==1) then 523 ABI_MALLOC(matlu1,(natom)) 524 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 525 call copy_matlu(green%occup%matlu,matlu1,natom) 526 if(pawprtvol>=3) then 527 write(message,'(a,2x,a)') ch10,& 528 & " == Occupations before rotations" 529 call wrtout(std_out,message,'COLL') 530 call print_matlu(green%occup%matlu,natom,1) 531 end if 532 533 ! 1) rotate density matrix to Ylm basis 534 ! -------------------------------------- 535 if(useylm==1) then 536 call slm2ylm_matlu(matlu1,natom,1,pawprtvol) 537 if(pawprtvol>=3) then 538 write(message,'(a,a)') ch10, " == Print occupations in Ylm basis" 539 call wrtout(std_out,message,'COLL') 540 call print_matlu(matlu1,natom,1) 541 end if 542 end if 543 544 ! 2) rotate density matrix to rotated basis 545 ! ------------------------------------------- 546 if(opt_rot==1.or.opt_rot==2) call rotate_matlu(matlu1,eigvectmatlu,natom,3,1) 547 write(message,'(a,2x,a,f13.5)') ch10," == Rotated occupations (for information)" 548 call wrtout(std_out,message,'COLL') 549 call print_matlu(matlu1,natom,1,compl=1) 550 call checkreal_matlu(matlu1,natom,tol10) 551 call destroy_matlu(matlu1,natom) 552 ABI_FREE(matlu1) 553 554 end if 555 556 call flush_unit(std_out) 557 558 559 ! ================================================================= 560 ! Rotate weiss function according to eigenvectors. 561 ! ================================================================= 562 !!!stop 563 ! Rotate Weiss function first in Ylm basis 564 ! ------------------------------------------------------------------- 565 if(useylm==1) then 566 write(message,'(a,2x,a)') ch10, " == Rotation of weiss and greendft in the Ylm Basis=" 567 call wrtout(std_out,message,'COLL') 568 do ifreq=1,paw_dmft%dmft_nwlo 569 call slm2ylm_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,1,0) 570 call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,1,0) 571 ! call slm2ylm_matlu(greendft%oper(ifreq)%matlu,natom,1,0) 572 end do 573 end if 574 575 if(pawprtvol>=3) then 576 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 577 ! " == Print weiss for small freq 1 before rot" ! debug 578 ! call wrtout(std_out,message,'COLL') ! debug 579 ! call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) ! debug 580 581 ! Print Weiss function 582 ! -------------------- 583 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 584 & " == Print weiss for 1st freq before rot" ! debug 585 call wrtout(std_out,message,'COLL') ! debug 586 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) ! debug 587 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 588 & " == Print weiss for last freq before rot" ! debug 589 call wrtout(std_out,message,'COLL') ! debug 590 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 591 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 592 !& " == Print DFT G for 1st freq before rot" ! debug 593 ! call wrtout(std_out,message,'COLL') ! debug 594 ! call print_matlu(greendft%oper(1)%matlu,natom,1,compl=1,opt_exp=2) ! debug 595 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 596 !& " == Print DFT G for last freq before rot" ! debug 597 ! call wrtout(std_out,message,'COLL') ! debug 598 ! call print_matlu(greendft%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) ! debug 599 end if 600 601 if(opt_diag/=0) then 602 ! Rotate Weiss function from the Slm (or Ylm) to the basis of diagonalisation 603 ! ------------------------------------------------------------------- 604 write(message,'(a,2x,a)') ch10, " == Rotation of weiss =" 605 call wrtout(std_out,message,'COLL') 606 do ifreq=1,paw_dmft%dmft_nwlo 607 if(opt_rot==1) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 608 if(opt_rot==1) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 609 ! call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6) 610 end do 611 612 if(paw_dmft%myproc .eq. mod(nproc+1,nproc)) then 613 if (open_file(trim(paw_dmft%filapp)//"_atom__G0w_.dat", message, newunit=unt) /= 0) then 614 ABI_ERROR(message) 615 end if 616 do ifreq=1,paw_dmft%dmft_nwlo 617 write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),& 618 & (((weiss_for_rot%oper(ifreq)%matlu(1)%mat(im1,im1,isppol,ispinor,ispinor),& 619 & im1=1,3),ispinor=1,nspinor),isppol=1,nsppol) 620 end do 621 close(unt) 622 end if 623 624 call flush_unit(std_out) 625 if(pawprtvol>=3) then 626 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 627 & " == Print weiss for small freq 1 after rot" ! debug 628 call wrtout(std_out,message,'COLL') ! debug 629 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) ! debug 630 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 631 & " == Print weiss for last freq after rot" ! debug 632 call wrtout(std_out,message,'COLL') ! debug 633 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 634 end if 635 636 ! ! Rotate DFT Green's function first in Ylm basis then in the rotated basis and compare to weiss_for_rot 637 ! ! ----------------------------------------------------------------------------------------------------- 638 ! write(message,'(a,2x,a)') ch10, " == Rotation of greendft =" 639 ! call wrtout(std_out,message,'COLL') 640 ! do ifreq=1,paw_dmft%dmft_nwlo 641 ! if(opt_rot==1) call rotate_matlu(greendft%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 642 ! call diff_matlu("Weiss_for_rot","greendft",weiss_for_rot%oper(ifreq)%matlu,greendft%oper(ifreq)%matlu,natom,1,tol14) 643 !! call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6) 644 ! end do 645 ! if(pawprtvol>=3) then 646 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 647 !& " == Print greendft for small freq 1 after rot" ! debug 648 ! call wrtout(std_out,message,'COLL') ! debug 649 ! call print_matlu(greendft%oper(1)%matlu,natom,1,compl=1,opt_exp=2) ! debug 650 ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug 651 !& " == Print greendft for last freq after rot" ! debug 652 ! call wrtout(std_out,message,'COLL') ! debug 653 ! call print_matlu(greendft%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) ! debug 654 ! end if 655 ! call flush_unit(std_out) 656 end if 657 658 ! ================================================================= 659 ! Compute analytic limit of hybridization and rotate it 660 ! ================================================================= 661 ABI_MALLOC(hybri_coeff,(paw_dmft%natom)) 662 call init_matlu(paw_dmft%natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,hybri_coeff) 663 !write(6,*)"hybri1",hybri_coeff(1)%mat(1,1,1,1,1),paw_dmft%natom,cryst_struc%natom 664 665 ! Compute analytical C_ij such that F_ij -> C_ij/iw_n 666 ! --------------------------------------- 667 call hybridization_asymptotic_coefficient(cryst_struc,paw_dmft,pawang,hybri_coeff) 668 write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n for large frequency" 669 call wrtout(std_out,message,'COLL') 670 671 ! Print analytical C_ij (not rotated) 672 ! --------------------------------------- 673 call print_matlu(hybri_coeff,natom,1) 674 675 ! Rotate analytical C_ij in Ylm basis 676 ! --------------------------------------- 677 if(useylm==1) call slm2ylm_matlu(hybri_coeff,natom,1,pawprtvol) 678 if(opt_diag/=0) then 679 680 ! Rotate analytical C_ij in rotated basis 681 ! --------------------------------------- 682 if(opt_rot==1.or.opt_rot==2) call rotate_matlu(hybri_coeff,eigvectmatlu,natom,3,1) 683 684 ! Print analytical C_ij (rotated) 685 ! --------------------------------------- 686 write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n after rotation" 687 call wrtout(std_out,message,'COLL') 688 call print_matlu(hybri_coeff,natom,1,compl=1,opt_exp=1) 689 end if 690 691 ! ================================================================= 692 ! Check if rotation is properly done. 693 ! ================================================================= 694 if(3==4) then 695 write(message,'(a,2x,a)') ch10,& 696 & " == Print dmat before rot" 697 call wrtout(std_out,message,'COLL') 698 call print_matlu(green%occup%matlu,natom,1) 699 if(useylm==1) call slm2ylm_matlu(green%occup%matlu,natom,1,pawprtvol) 700 if(opt_rot==1) call rotate_matlu(green%occup%matlu,eigvectmatlu,natom,3,1) 701 write(message,'(a,2x,a)') ch10,& 702 & " == Print dmat after rot" 703 call wrtout(std_out,message,'COLL') 704 call print_matlu(green%occup%matlu,natom,1) 705 706 write(message,'(2a)') ch10,' QMC STOP: DEBUG' 707 call wrtout(std_out,message,'COLL') 708 ABI_ERROR(message) 709 end if 710 ! ================================================================= 711 ! Check 712 ! ================================================================= 713 714 ! write(message,'(a,2x,a,f13.5)') ch10,& 715 !& " == Print weiss for small tau" 716 ! call wrtout(std_out,message,'COLL') 717 ! call print_matlu(weiss%oper(1)%matlu,natom,1) 718 ! write(message,'(a,2x,a,f13.5)') ch10,& 719 !& " == Print weiss for large tau" 720 ! call wrtout(std_out,message,'COLL') 721 ! call print_matlu(weiss%oper(paw_dmft%dmft_nwlo)%matlu,natom,1) 722 ! call flush_unit(std_out) 723 ! write(message,'(2a)') ch10,' Check weiss_for_rot(last freq)' 724 ! call wrtout(std_out,message,'COLL') 725 ! call checkdiag_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,tol6,opt=nspinor) 726 ! call flush_unit(std_out) 727 ! write(message,'(2a)') ch10,' Check weiss_for_rot(ifreq=1)' 728 ! call wrtout(std_out,message,'COLL') 729 ! call checkdiag_matlu(weiss_for_rot%oper(1)%matlu,natom,tol6,opt=nspinor) 730 ! call flush_unit(std_out) 731 732 master=0 733 734 ! ================================================================= 735 ! Print out 736 ! ================================================================= 737 738 ! Print Weiss 739 ! ------------- 740 if(paw_dmft%dmft_prgn==1) then 741 call print_green('Weiss_diag',weiss_for_rot,1,paw_dmft,pawprtvol=1,opt_wt=1,opt_decim=1) 742 end if 743 744 write(message,'(a,2x,a,f13.5)') ch10,& 745 & " == Preparing data for CTQMC" 746 call wrtout(std_out,message,'COLL') 747 748 ! Print Rotate Weiss for 1st and last frequencies 749 ! ------------------------------------------------ 750 if (pawprtvol>=3) then 751 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 752 & " == Print rotated weiss function for small freq in the rotated basis" ! debug 753 call wrtout(std_out,message,'COLL') ! debug 754 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) ! debug 755 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 756 & " == Print rotated weiss function for largest freq in the rotated basis" ! debug 757 call wrtout(std_out,message,'COLL') ! debug 758 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 759 end if 760 761 ! ================================================================= 762 ! VARIABLES FOR CTQMC TESTS 763 testcode = 0 764 testrot = 0 765 opt_fk=0 ! for developpers to check Fourier transform and computes G0(tau) 766 opt_fk=1 ! usual case: for real calculations 767 ! ================================================================= 768 769 ! ___________________________________________________________________________________ 770 ! 771 ! SECOND PART : BUILT HYBRIDIZATION FROM G0 772 ! ___________________________________________________________________________________ 773 ! 774 ! =========================================================================================== 775 ! Compute inverse of weiss and compute hybridization 776 ! =========================================================================================== 777 778 ! Compute inverse of weiss for each Frequency 779 ! ---------------------------------------------- 780 do ifreq=1,paw_dmft%dmft_nwlo 781 ABI_MALLOC(matlu1,(natom)) 782 ABI_MALLOC(matlu2,(natom)) 783 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 784 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2) 785 786 call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,natom) 787 788 ! Print G_0(iw_n) 789 ! ---------------- 790 if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"go",60000,imre=1) 791 792 ! Compute G_0^-1 793 ! ------------------------------------------- 794 ! if opt_fk=1 or testcode/=0 Do the inversion 795 ! if opt_fk=0 Do not inverse. 796 ! If testcode=2 and opt_fk=0 Do the inversion 797 ! If testcode=1 and opt_fk=0 Do the inversion but no effect, because it will nevertheless be erased 798 ! If opt_fk=1 Do the inversion 799 ! ------------------------------------------- 800 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"weiss",12000,imre=1) 801 if(opt_fk==1.or.testcode/=0) call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1) 802 803 ! Print G_0^-1(iw_n) 804 ! ---------------- 805 if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"goinv",70000,imre=1) 806 807 if(pawprtvol>=4.or.ifreq==paw_dmft%dmft_nwlo) then 808 if(opt_fk==1.or.testcode/=0) then 809 ! Check inversion : do the product 810 ! ---------------------------------------------- 811 call prod_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,matlu2,natom) 812 write(message,'(a,2x,a,i7)') ch10,& ! debug 813 & " == Print product of weiss times invers for freq",ifreq 814 call wrtout(std_out,message,'COLL') ! debug 815 call print_matlu(matlu2,natom,1) ! debug 816 end if 817 end if 818 819 call destroy_matlu(matlu1,natom) 820 call destroy_matlu(matlu2,natom) 821 ABI_FREE(matlu1) 822 ABI_FREE(matlu2) 823 end do 824 825 ! Copy weiss_for_rot into weiss 826 ! ------------------------------- 827 !call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,weiss%oper(ifreq)%matlu,natom) 828 829 830 ! Print G_0^-1 for 1st and last frequencies. 831 ! ----------------------------------------- 832 if(pawprtvol>=3) then 833 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 834 & " == Print G_0^-1 for small freq in the rotated basis" ! debug 835 call wrtout(std_out,message,'COLL') ! debug 836 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) ! debug 837 write(message,'(a,2x,a,e18.10,a)') ch10,& ! debug 838 & " == Print G_0^-1 for last freq in the rotated basis (last freq=", paw_dmft%omega_lo(paw_dmft%dmft_nwlo),")" ! debug 839 call wrtout(std_out,message,'COLL') ! debug 840 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 841 end if 842 843 ! Substract frequency from diagonal part 844 ! ====================================== 845 846 ABI_MALLOC(shift,(natom)) 847 do ifreq=1,paw_dmft%dmft_nwlo 848 shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 849 850 ! write(5555,'(400e17.4)') paw_dmft%omega_lo(ifreq),((((((weiss_for_rot%oper(ifreq)%matlu(1)%mat& 851 ! & (im,im1,isppol,ispinor,ispinor1)-cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)),im=1,2*3+1),& 852 !& im1=1,2*3+1),isppol=1,nsppol),ispinor=1,nspinor),ispinor1=1,nspinor) 853 854 ! Compute G_0^-1-iw_n 855 ! -------------------- 856 if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift) 857 858 ! Compute -G_0^-1+iw_n 859 ! -------------------- 860 if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone) 861 862 ! Print -G_0^-1+iw_n 863 ! -------------------- 864 if(optdb==1) then 865 call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"G0inv_minus_omega",20000,imre=1) 866 end if 867 end do 868 869 ! Print -G_0^+1-iw_n=(F-levels) for last freq in the rotated basis" 870 ! ------------------------------------------------------------------ 871 ABI_FREE(shift) 872 if(pawprtvol>=3) then 873 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 874 & " == Print G_0^-1-iw_n=-(F-levels) for last freq in the rotated basis" ! debug 875 call wrtout(std_out,message,'COLL') ! debug 876 call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug 877 end if 878 879 ! Check numerical limit of F(i_wn)*iw_n (can be used also to compute F ) 880 ! ====================================== 881 882 if(opt_nondiag==1) then 883 ABI_MALLOC(matlu1,(natom)) 884 ABI_MALLOC(matlu2,(natom)) 885 ABI_MALLOC(matlu3,(natom)) 886 ABI_MALLOC(matlu4,(natom)) 887 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 888 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2) 889 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu3) 890 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu4) 891 892 write(message,'(a,2x,a,f13.5)') ch10, " == energy_levels" 893 call wrtout(std_out,message,'COLL') 894 call print_matlu(energy_level%matlu,natom,1,opt_exp=2,compl=1) 895 896 do ifreq=paw_dmft%dmft_nwlo,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency. 897 !do ifreq=paw_dmft%dmftqmc_l,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency. 898 ! Compute F (substract levels) for max frequency 899 ! ----------------------------------------------- 900 call add_matlu(weiss_for_rot%oper(ifreq)%matlu,energy_level%matlu,matlu1,natom,-1) 901 902 ! Print F(iw_n)=-(G_0^-1-iw_n+levels) for last frequency. 903 ! -------------------------------------------------------- 904 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 905 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 906 & " == Print F(iw_n)=-(G_0^-1-iw_n+levels) for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 907 call wrtout(std_out,message,'COLL') 908 call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 909 end if 910 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"Hybridization",10000,imre=1) 911 912 ! Put F in weiss_for_rot -> CTQMC 913 ! ------------------------------- 914 if(opt_rot==2) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1) 915 ! The following line will produce directly the weiss function for the CTQMC code 916 if(opt_fk==1) call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom) 917 918 919 ! Multiply F by frequency 920 ! ------------------------ 921 call copy_matlu(matlu1,matlu2,natom) 922 call fac_matlu(matlu1,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)) 923 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 924 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 925 & " == Print numerical C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 926 call wrtout(std_out,message,'COLL') 927 call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 928 end if 929 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij",72800,imre=1) 930 !call rotate_matlu(matlu1,eigvectmatlu,natom,3,1) 931 932 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 933 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 934 & " == Print numerical after back rotation C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 935 call wrtout(std_out,message,'COLL') 936 call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 937 end if 938 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_rotated",72900,imre=1) 939 940 ! Built C_ij/iw_n 941 ! ------------------------ 942 call copy_matlu(hybri_coeff,matlu1,natom) 943 call fac_matlu(matlu1,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)) 944 if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_over_omega",72000) 945 ! if(ifreq==paw_dmft%dmft_nwlo) then 946 ! write(message,'(a,2x,a,f13.5)') ch10, " == Print numerical C_ij/iw_n for frequency",paw_dmft%omega_lo(ifreq) 947 ! call wrtout(std_out,message,'COLL') 948 ! call print_matlu(matlu1,natom,1,opt_exp=1,compl=1) 949 ! endif 950 951 ! For test: put C_ij/i_wn into weiss_for_rot 952 ! -------------------------------------------- 953 !call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1) 954 955 ! Compute Hybri - C_ij/iw_n 956 ! ------------------------ 957 call add_matlu(matlu2,matlu1,matlu3,natom,-1) 958 959 ! Print Hybri - C_ij/iw_n 960 ! ------------------------ 961 if(optdb==1) call printplot_matlu(matlu3,natom,paw_dmft%omega_lo(ifreq),"hybri_minus_asymp",74000,imre=1) 962 963 ! 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. 964 ! ------------------------------------------------------------------------------------------------------------------ 965 call copy_matlu(matlu3,matlu2,natom) 966 call fac_matlu(matlu2,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2) 967 if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"fminuscijtimesw2",75000,imre=1) 968 if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then 969 call copy_matlu(matlu2,matlu4,natom) 970 write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, & 971 & " == Print numerical (F(iw_n)-C_ij/iw_n)%iw_n^2 for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")" 972 call wrtout(std_out,message,'COLL') 973 call print_matlu(matlu4,natom,1) 974 end if 975 976 ! Built C_ij/iw_n+D_ij/(iw_n)^2 977 ! ------------------------ 978 call copy_matlu(matlu4,matlu3,natom,opt_re=1) 979 call fac_matlu(matlu3,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2) 980 call add_matlu(matlu1,matlu3,matlu2,natom,1) 981 if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"cij_w_plus_dij_w2",72700,imre=1) 982 ! For test: put C_ij/i_wn +D_ij/(iw_n)^2 into weiss_for_rot 983 ! -------------------------------------------- 984 !call copy_matlu(matlu2,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1) 985 986 987 end do 988 989 call destroy_matlu(matlu1,natom) 990 call destroy_matlu(matlu2,natom) 991 call destroy_matlu(matlu3,natom) 992 call destroy_matlu(matlu4,natom) 993 ABI_FREE(matlu1) 994 ABI_FREE(matlu2) 995 ABI_FREE(matlu3) 996 ABI_FREE(matlu4) 997 end if ! if opt_nondiag=1 998 999 ! ========================================================================================= 1000 ! Start big loop over atoms to compute hybridization and do the CTQMC 1001 ! ========================================================================================= 1002 do iatom=1,cryst_struc%natom 1003 green%ecorr_qmc(iatom)=zero 1004 itypat=cryst_struc%typat(iatom) 1005 lpawu=paw_dmft%lpawu(iatom) 1006 tndim=2*lpawu+1 1007 if(lpawu/=-1) then 1008 1009 nflavor=2*(tndim) 1010 if(testcode>=1) then 1011 nflavor=2 1012 if(testcode==2) then 1013 ispa=1 1014 ispb=2 1015 if(nspinor==1) ispb=1 1016 ima=1 1017 imb=1 1018 if(tndim>4) then 1019 ima=5 ! row 1020 imb=4 ! column 1021 end if 1022 end if 1023 end if 1024 1025 ABI_MALLOC(fw1,(paw_dmft%dmft_nwlo,nflavor)) 1026 ABI_MALLOC(fw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor)) 1027 ABI_MALLOC(levels_ctqmc,(nflavor)) 1028 ABI_MALLOC(levels_ctqmc_nd,(nflavor,nflavor)) 1029 levels_ctqmc_nd=czero 1030 ABI_MALLOC(hybri_limit,(nflavor,nflavor)) 1031 hybri_limit=czero 1032 fw1_nd=czero 1033 fw1=czero 1034 1035 ! ================================================================= 1036 ! Put hybridization in new arrays for CTQMC 1037 ! ================================================================= 1038 if (testcode==0) then 1039 iflavor1=0 1040 iflavor2=0 1041 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 1048 ! first diagonal terms whatever opt_nondiag 1049 iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1) 1050 iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1) 1051 1052 if ( iflavor1==iflavor2 ) then 1053 1054 ! Put weiss_for_rot in fw1 1055 do ifreq=1,paw_dmft%dmft_nwlo 1056 if(opt_fk==1) then 1057 fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1) 1058 else if (opt_fk==0) then 1059 fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1) 1060 end if 1061 end do 1062 fw1_nd(:,iflavor1,iflavor1)=fw1(:,iflavor1) 1063 1064 levels_ctqmc(iflavor1)=real(energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1),kind=dp) 1065 hybri_limit(iflavor1,iflavor1)=hybri_coeff(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1) 1066 1067 1068 ! case nsppol=nspinor=1 1069 if(nsppol==1.and.nspinor==1) then 1070 fw1(:,iflavor1+tndim)=fw1(:,iflavor1) 1071 fw1_nd(:,iflavor1+tndim,iflavor1+tndim)=fw1(:,iflavor1) 1072 levels_ctqmc(iflavor1+tndim)=levels_ctqmc(iflavor1) 1073 hybri_limit(iflavor1+tndim,iflavor1+tndim)=hybri_limit(iflavor1,iflavor1) 1074 end if 1075 1076 ! off diagonal terms 1077 else 1078 1079 ! Put weiss_for_rot in fw1_nd 1080 do ifreq=1,paw_dmft%dmft_nwlo 1081 if(opt_fk==1) then 1082 fw1_nd(ifreq,iflavor1,iflavor2)= & 1083 & weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 1084 else if (opt_fk==0) then 1085 fw1_nd(ifreq,iflavor1,iflavor2)= & 1086 & weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 1087 end if 1088 end do 1089 hybri_limit(iflavor1,iflavor2)=hybri_coeff(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2) 1090 1091 ! case nsppol=nspinor=1 1092 if(nsppol==1.and.nspinor==1) then 1093 fw1_nd(:,iflavor1+tndim,iflavor2+tndim) = fw1_nd(:,iflavor1,iflavor2) 1094 hybri_limit(iflavor1+tndim,iflavor2+tndim)=hybri_limit(iflavor1,iflavor2) 1095 end if 1096 1097 end if 1098 1099 ! < / HACK > 1100 end do !im2 1101 end do !im1 1102 end do !ispinor2 1103 end do !ispinor1 1104 end do !isppol 1105 ! < HACK > 1106 ! JB. On 1000 cpus this can not work since all CPU try to open/write the files 1107 ! Action : Don't print it or check only one cpu does it. 1108 1109 if(pawprtvol>=10000000) then 1110 write(message,'(a,2x,a)') ch10, " == Hybri for all flavors for CTQMC " 1111 call wrtout(std_out,message,'COLL') 1112 do iflavor1=1,nflavor 1113 write(message,'(4x,14(2e14.5,2x))') (hybri_limit(iflavor1,iflavor2),iflavor2=1,nflavor) 1114 call wrtout(std_out,message,'COLL') 1115 end do 1116 1117 if (open_file('Hybri_cijoveromega',message, newunit=unt, status='unknown', form='formatted') /= 0) then 1118 ABI_ERROR(message) 1119 end if 1120 if (open_file('Hybri',message,newunit=unt2,status='unknown',form='formatted') /= 0) then 1121 ABI_ERROR(message) 1122 end if 1123 do ifreq=1,paw_dmft%dmft_nwlo 1124 ! weiss_for_rot is G_0^-1-iw_n=-(F-levels) 1125 if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"weissbefore112",30000) 1126 end do 1127 do iflavor1=1,nflavor 1128 do iflavor2=1,nflavor 1129 do ifreq=1,paw_dmft%dmft_nwlo 1130 omega=pi*paw_dmft%temp*(two*float(ifreq)-1) 1131 ! fw1_nd is -G_0^+1-iw_n=(F-levels) 1132 write(unt,'(300e16.5)') paw_dmft%omega_lo(ifreq)& 1133 & ,fw1_nd(ifreq,iflavor1,iflavor2)-hybri_limit(iflavor1,iflavor2)/cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp) 1134 write(unt2,'(300e16.5)') paw_dmft%omega_lo(ifreq),fw1_nd(ifreq,iflavor1,iflavor2) 1135 end do 1136 write(unt,*) 1137 write(unt2,*) 1138 end do 1139 end do 1140 close(unt) 1141 close(unt2) 1142 end if 1143 end if ! testcode 1144 ! </ HACK > 1145 1146 ! ==================================================================================== 1147 ! TEST 1148 ! For testing purpose, built ultra simple hybridization (constant in 1149 ! imaginary time or very simple) or extract some part of the calculated hybridization 1150 ! ==================================================================================== 1151 if(testcode>=1) then 1152 dmft_nwlo=paw_dmft%dmft_nwlo 1153 paw_dmft%dmft_nwlo=paw_dmft%dmftqmc_l 1154 ABI_MALLOC(gw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor)) 1155 gw1_nd=czero 1156 1157 ! Call testcode_ctqmc: built simple hybridization 1158 !-------------------------------------------------- 1159 if (testcode==1) then 1160 call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,& 1161 & levels_ctqmc,hybri_limit,nflavor,1,paw_dmft%temp,testrot,testcode,umod) 1162 ! Select 2x2 hybridization matrix from the current larger matrix 1163 ! ima and imb are defined above. 1164 !---------------------------------------------------------------- 1165 else if (testcode==2) then 1166 !close(unt) 1167 !close(unt2) 1168 call testcode_ctqmc_b(energy_level,hybri_coeff,weiss_for_rot,paw_dmft%dmftqmc_l,fw1_nd,& 1169 & levels_ctqmc,levels_ctqmc_nd,hybri_limit,paw_dmft%temp,umod,opt_diag,opt_fk) 1170 end if 1171 1172 ! Calculation of Inverse Green's function from hybridization 1173 !------------------------------------------------------------- 1174 do if1=1,2 1175 do if2=1,2 1176 do ifreq=1,paw_dmft%dmftqmc_l 1177 omega=pi*paw_dmft%temp*(two*float(ifreq)-1) 1178 if(if1==if2) then 1179 gw1_nd(ifreq,if1,if2) = (cmplx(0.d0,omega,kind=dp)-fw1_nd(ifreq,if1,if2)) 1180 else 1181 gw1_nd(ifreq,if1,if2) = (-fw1_nd(ifreq,if1,if2)) 1182 end if 1183 end do 1184 end do 1185 end do 1186 ! Calculation of Green's function (call inverse) 1187 !------------------------------------------------------------- 1188 do ifreq=1,paw_dmft%dmftqmc_l 1189 call xginv(gw1_nd(ifreq,:,:),2) 1190 end do 1191 write(std_out,*) " testctqmc high frequency limit of hybridization",fw1_nd(paw_dmft%dmftqmc_l,:,:) 1192 1193 ! Integrate Green's function 1194 !------------------------------------------------------------- 1195 do if1=1,2 1196 do if2=1,2 1197 call int_fct(gw1_nd(:,if1,if2),(if1==if2),2,paw_dmft,integral(if1,if2)) ! test_1 1198 end do 1199 end do 1200 ! Write Occupations 1201 write(std_out,*) "Occupation of model in matrix form" 1202 do if1=1,2 1203 write(std_out,'(2(2f13.5,3x))') ((integral(if1,if2)+conjg(integral(if2,if1)))/two,if2=1,2) 1204 end do 1205 write(std_out,*) "Limit of hybridization " 1206 do if1=1,2 1207 write(std_out,'(2(2f13.5,3x))') (hybri_limit(if1,if2),if2=1,2) 1208 end do 1209 1210 ! If opt_fk=0, give Green's function to CTQMC code instead of 1211 ! hybridization 1212 !------------------------------------------------------------- 1213 if(opt_fk==0) then 1214 fw1_nd=gw1_nd 1215 end if 1216 1217 ABI_FREE(gw1_nd) 1218 paw_dmft%dmft_nwlo=dmft_nwlo 1219 1220 ! and testcode>1 1221 end if 1222 1223 1224 call flush_unit(std_out) 1225 ! ================================================================= 1226 1227 ! ___________________________________________________________________________________ 1228 ! 1229 ! THIRD PART : CALL CTQMC 1230 ! ___________________________________________________________________________________ 1231 1232 ! ================================================================= 1233 ! Main calls to CTQMC code in ABINIT (INITIALIZATION and OPTIONS) 1234 ! ================================================================= 1235 if(paw_dmft%dmft_solv==5.or.paw_dmft%dmft_solv==8) then 1236 1237 write(message,'(a,2x,a)') ch10,& 1238 & " == Initializing CTQMC" 1239 call wrtout(std_out,message,'COLL') 1240 1241 ! Initialisation 1242 ! ================================================================= 1243 if(paw_dmft%dmft_solv==5) then 1244 nomega=paw_dmft%dmftqmc_l 1245 call CtqmcInterface_init(hybrid,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, & 1246 & paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,& 1247 & std_out,paw_dmft%spacecomm) 1248 ! options 1249 ! ================================================================= 1250 call CtqmcInterface_setOpts(hybrid,& 1251 opt_Fk =opt_fk,& 1252 & opt_order =paw_dmft%dmftctqmc_order ,& 1253 & opt_histo =paw_dmft%dmftctqmc_config ,& 1254 & opt_movie =paw_dmft%dmftctqmc_mov ,& 1255 & opt_analysis=paw_dmft%dmftctqmc_correl,& 1256 & opt_check =paw_dmft%dmftctqmc_check ,& 1257 & opt_noise =paw_dmft%dmftctqmc_grnns ,& 1258 & opt_spectra =paw_dmft%dmftctqmc_mrka ,& 1259 & opt_gmove =paw_dmft%dmftctqmc_gmove ) 1260 endif 1261 1262 if(paw_dmft%dmft_solv==8) then 1263 nomega=paw_dmft%dmftqmc_l 1264 call CtqmcoffdiagInterface_init(hybridoffdiag,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, & 1265 & paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,& 1266 & paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,& 1267 & std_out,paw_dmft%spacecomm,opt_nondiag) 1268 ! options 1269 ! ================================================================= 1270 call CtqmcoffdiagInterface_setOpts(hybridoffdiag,& 1271 opt_Fk =opt_fk,& 1272 & opt_order =paw_dmft%dmftctqmc_order ,& 1273 & opt_histo =paw_dmft%dmftctqmc_config ,& 1274 & opt_movie =paw_dmft%dmftctqmc_mov ,& 1275 & opt_analysis=paw_dmft%dmftctqmc_correl,& 1276 & opt_check =paw_dmft%dmftctqmc_check ,& 1277 & opt_noise =paw_dmft%dmftctqmc_grnns ,& 1278 & opt_spectra =paw_dmft%dmftctqmc_mrka ,& 1279 & opt_gmove =paw_dmft%dmftctqmc_gmove ) 1280 endif 1281 1282 write(message,'(a,2x,2a)') ch10, " == Initialization CTQMC done", ch10 1283 call wrtout(std_out,message,'COLL') 1284 end if 1285 1286 if(paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7.or.paw_dmft%dmft_solv==9) then 1287 ABI_MALLOC(gw_tmp_nd,(paw_dmft%dmft_nwli,nflavor,nflavor)) 1288 !because size allocation problem with TRIQS paw_dmft%dmft_nwlo must be >= paw_dmft%dmft_nwli 1289 open(unit=505,file=trim(paw_dmft%filapp)//"_Legendre_coefficients.dat", status='unknown',form='formatted') 1290 else 1291 if(paw_dmft%dmft_solv==5) then 1292 ABI_MALLOC(gw_tmp,(paw_dmft%dmft_nwlo,nflavor+1)) 1293 end if 1294 ABI_MALLOC(gw_tmp_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor+1)) 1295 !use gw_tmp to put freq 1296 do ifreq=1,paw_dmft%dmft_nwlo 1297 if(paw_dmft%dmft_solv==5) gw_tmp(ifreq,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 1298 gw_tmp_nd(ifreq,nflavor,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 1299 end do 1300 end if 1301 1302 ABI_MALLOC(gtmp,(paw_dmft%dmftqmc_l,nflavor)) 1303 ! THIS IS A BACKUP PLAN. USING paw_dmft%hybrid makes a segfault on TIKAL 1304 ! PSC with MPI only (and max2_open64). paw_dmf%hybrid is corrupted 1305 ! somewhere but I could not find the place in all DMFT routines 1306 ABI_MALLOC(gtmp_nd,(paw_dmft%dmftqmc_l,nflavor,nflavor)) 1307 call flush_unit(std_out) 1308 1309 ! ================================================================= 1310 ! BEGIN CALL TO CTQMC SOLVERS 1311 ! ================================================================= 1312 1313 if(testcode==0) then 1314 1315 ! ================================================================= 1316 ! CTQMC run Abinit 1317 ! ================================================================= 1318 if(paw_dmft%dmft_solv==5) then 1319 1320 ABI_MALLOC(docc,(1:nflavor,1:nflavor)) 1321 docc(:,:) = zero 1322 call CtqmcInterface_run(hybrid,fw1(1:paw_dmft%dmftqmc_l,:),Gtau=gtmp,& 1323 & Gw=gw_tmp,D=docc(:,:),E=green%ecorr_qmc(iatom),& 1324 !& matU=hu(itypat)%udens,opt_levels=levels_ctqmc) 1325 & matU=udens_atoms(iatom)%value,opt_levels=levels_ctqmc) 1326 call data4entropyDMFT_setDocc(paw_dmft%forentropyDMFT,iatom,docc) 1327 ABI_FREE(docc) 1328 !DO iflavor = 1, nflavor 1329 ! hybrid%Hybrid%Greens(iflavor)%oper(1:this%samples) = gtmp(1:this%samples,iflavor) 1330 ! CALL GreenHyb_forFourier(this%Greens(iflavor), Gomega=Gw(:,iflavor), omega=Gw(:,this%flavors+1)) 1331 !END DO 1332 1333 ! ================================================================= 1334 ! CTQMC run Abinit off diagonal terms in hybridization 1335 ! ================================================================= 1336 else if (paw_dmft%dmft_solv==8) then 1337 ! ================================================================= 1338 1339 ABI_MALLOC(docc,(1:nflavor,1:nflavor)) 1340 docc(:,:) = zero 1341 call CtqmcoffdiagInterface_run(hybridoffdiag,fw1_nd(1:paw_dmft%dmftqmc_l,:,:),Gtau=gtmp_nd,& 1342 & Gw=gw_tmp_nd,D=doccsum,E=green%ecorr_qmc(iatom),& 1343 & Noise=noise,matU=udens_atoms(iatom)%value,Docc=docc,opt_levels=levels_ctqmc,hybri_limit=hybri_limit) 1344 ! For entropy (alternative formulation) 1345 if(paw_dmft%ientropy==1) then 1346 EE=zero 1347 do if1=1,nflavor 1348 do if2=if1+1,nflavor 1349 EE=EE+docc(if1,if2)*udens_atoms_for_s(iatom)%value(if1,if2) 1350 ! write(std_out,*) udens_atoms_for_s(iatom)%value(if1,if2),docc(if1,if2) 1351 enddo 1352 enddo 1353 ! Here in udens U=1, J=J/U, so we need to multiply bu U/Ha_eV 1354 write(message,'(a,3(f14.10,3x))') "For entropy calculation E_corr_qmc, u_for_s, j_for,s", & 1355 & paw_dmft%u_for_s*EE/Ha_eV,paw_dmft%u_for_s,paw_dmft%j_for_s 1356 call wrtout(std_out,message,'COLL') 1357 EE=zero 1358 do if1=1,nflavor 1359 do if2=if1+1,nflavor 1360 EE=EE+docc(if1,if2)*udens_atoms(iatom)%value(if1,if2) 1361 ! write(std_out,*) udens_atoms(iatom)%value(if1,if2),docc(if1,if2) 1362 enddo 1363 enddo 1364 ! Here in udens U=U, J=J, so we obtain directly the results 1365 write(message,'(a,3(f14.10,3x))') "Reference calculation E_corr_qmc, upawu , jpawu ", & 1366 & EE,hu(itypat)%upawu*Ha_eV,hu(itypat)%jpawu*Ha_eV 1367 call wrtout(std_out,message,'COLL') 1368 endif 1369 ABI_FREE(docc) 1370 ! TODO: Handle de luj0 case for entropy 1371 1372 ! ================================================================= 1373 ! CTQMC run TRIQS 1374 ! ================================================================= 1375 else if ((paw_dmft%dmft_solv>=6.and.paw_dmft%dmft_solv<=7).or.paw_dmft%dmft_solv==9) then 1376 ! ================================================================= 1377 1378 if ( paw_dmft%dmft_solv==9.or.(paw_dmft%dmftqmc_l >= (2*paw_dmft%dmft_nwli)+1) ) then 1379 1380 call ctqmc_calltriqs(paw_dmft,cryst_struc,hu,levels_ctqmc,gtmp_nd,gw_tmp_nd,fw1_nd,leg_measure,iatom) 1381 1382 else 1383 write(message,'(2a)') ch10," Can't launch TRIQS CTHYB solver because dmftqmc_l must be >= 2*dmft_nwli + 1" 1384 ABI_ERROR(message) 1385 end if 1386 1387 end if 1388 1389 ! ================================================================= 1390 ! CTQMC run for tests 1391 ! ================================================================= 1392 else if (testcode>=1) then 1393 call CtqmcInterface_run(hybrid,fw1(1:nomega,:),Gtau=gtmp,& 1394 & Gw=gw_tmp,E=green%ecorr_qmc(iatom),& 1395 & matU=umod,opt_levels=levels_ctqmc) 1396 1397 ! for non diagonal code 1398 ! call CtqmcInterface_run(hybrid,fw1_nd(1:nomega,:,:),Gtau=gtmp_nd,& 1399 !& Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),& 1400 !& Noise=Noise,matU=umod,opt_levels=levels_ctqmc,hybri_limit=hybri_limit) 1401 1402 ! If test of the code is activated, and testrot =1 rotate back green's function and stop the code. 1403 ! -------------------------------------------------------------------------------------------------- 1404 if(testcode==1) then 1405 1406 call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,& 1407 & levels_ctqmc,hybri_limit,nflavor,2,paw_dmft%temp,testrot,testcode,umod) 1408 1409 write(message,'(2a)') ch10,' testcode end of test calculation' 1410 ABI_ERROR(message) 1411 end if 1412 if(testcode==2) then 1413 write(message,'(2a)') ch10,' testcode 2 end of test calculation' 1414 ABI_ERROR(message) 1415 end if 1416 1417 end if 1418 ! ================================================================= 1419 ! END CALL TO CTQMC SOLVERS 1420 ! ================================================================= 1421 1422 1423 ! Print green function is files directly from CTQMC 1424 ! -------------------------------------------------- 1425 call ctqmcoutput_printgreen(paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom) 1426 1427 1428 ! If the CTQMC code in ABINIT was used, then destroy it and deallocate arrays 1429 ! ---------------------------------------------------------------------------- 1430 if(paw_dmft%dmft_solv<6.and.paw_dmft%dmft_solv>7) then 1431 !Nothing just hybrid var problem 1432 else 1433 write(message,'(a,2x,a)') ch10,& 1434 & " == Destroy CTQMC" 1435 call wrtout(std_out,message,'COLL') 1436 if(paw_dmft%dmft_solv==5) call CtqmcInterface_finalize(hybrid) 1437 if(paw_dmft%dmft_solv==8) call CtqmcoffdiagInterface_finalize(hybridoffdiag) 1438 write(message,'(a,2x,a)') ch10,& 1439 & " == Destroy CTQMC done" 1440 call wrtout(std_out,message,'COLL') 1441 end if 1442 ABI_FREE(hybri_limit) 1443 ABI_FREE(levels_ctqmc_nd) 1444 ABI_FREE(levels_ctqmc) 1445 ABI_FREE(fw1) 1446 ABI_FREE(fw1_nd) 1447 1448 ! ___________________________________________________________________________________ 1449 ! 1450 ! FOURTH PART : USE OUTPUT OF CTQMC AND THEN DO BACK ROTATION 1451 ! ___________________________________________________________________________________ 1452 ! 1453 1454 ! Put green's function values from CTQMC into green structure 1455 !------------------------------------------------------------- 1456 call ctqmcoutput_to_green(green,paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom,leg_measure,opt_nondiag) 1457 1458 ! Deallocate arrays for CTQMC 1459 !----------------------------- 1460 if(paw_dmft%dmft_solv<6) then 1461 ABI_FREE(gw_tmp) 1462 endif 1463 ABI_FREE(gw_tmp_nd) 1464 ABI_FREE(gtmp) 1465 ABI_FREE(gtmp_nd) 1466 1467 1468 1469 ! Do Fourier transform if it was not done (ie if TRIQS is used without legendre measurement) 1470 !---------------------------------------------------------------------------------------------- 1471 if(opt_nondiag==1) then ! (As leg_measure is activated by defautl, this fourier is never done). 1472 if(paw_dmft%dmft_solv>=6.and..not.leg_measure.and.paw_dmft%dmft_solv<=7) then 1473 write(message,'(2a,i3,13x,a)') ch10,' === Direct Fourier Transform t->w of Weiss Field' 1474 call wrtout(std_out,message,'COLL') 1475 call fourier_green(cryst_struc,green,paw_dmft,& 1476 & pawang,opt_ksloc=2,opt_tw=1) 1477 end if 1478 endif 1479 1480 1481 end if 1482 1483 end do ! iatom 1484 ! ========================================================================================= 1485 ! End big loop over atoms to compute hybridization and do the CTQMC 1486 ! ========================================================================================= 1487 1488 if(paw_dmft%dmft_prgn==1) then 1489 call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=2) 1490 call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 1491 end if 1492 !write(message,'(i3,4x,2e21.14)') 6,weiss_for_rot%oper(1)%matlu(1)%mat(1,1,1,1,1) 1493 !call wrtout(std_out,message,'COLL') ! debug 1494 ! ================================================================= 1495 ! Inverse Weiss, then 1496 ! Copy Weiss_for_rot into weiss and rotate back weiss to the original basis 1497 ! ================================================================= 1498 1499 ! ABI_MALLOC(shift,(natom)) 1500 ! do ifreq=1,paw_dmft%dmft_nwlo 1501 ! ! First weiss_for_rot contains -G_0^-1+iw_n 1502 ! ! ------------------------------------------- 1503 ! ! Compute G_0^-1-iw_n 1504 ! ! -------------------- 1505 ! write(6,*) "1" 1506 ! if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone) 1507 ! 1508 ! 1509 ! write(6,*) "2" 1510 ! ! Compute G_0^-1 1511 ! ! -------------------- 1512 ! shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp) 1513 ! if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift,signe=1) 1514 ! 1515 ! write(6,*) "3" 1516 ! ! Compute G_0 1517 ! ! -------------------- 1518 ! call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1) 1519 ! ! No need to copy if weiss_for_rot is a pointer to weiss ... 1520 !! if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0) 1521 !! if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 1522 ! 1523 ! ! Compute G_0 in the original basis 1524 ! ! -------------------- 1525 ! call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 1526 ! end do 1527 ! ABI_FREE(shift) 1528 1529 ! ================================================================= 1530 ! Here compute Self energy from Dyson and print it 1531 ! Warning : Weiss_for_rot is inversed inside dyson 1532 ! ================================================================= 1533 ! call initialize_self(self,paw_dmft) 1534 ! call dyson(green,paw_dmft,self,weiss_for_rot,opt_weissself=2) 1535 ! call rw_self(self,mpi_enreg,paw_dmft,prtopt=2,opt_rw=2,opt_char="diag") 1536 ! call destroy_self(self) 1537 !write(message,'(i3,4x,2e21.14)') 7,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1) 1538 !call wrtout(std_out,message,'COLL') ! debug 1539 1540 ! ================================================================= 1541 ! Rotate back green function to original basis (non-diagonal) 1542 ! (and Weiss for further use: might be useful if an back Fourier 1543 ! transformation is done). 1544 ! ================================================================= 1545 if(pawprtvol>=3) then 1546 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1547 & " == Print green function for small tau after CTQMC" ! debug 1548 call wrtout(std_out,message,'COLL') ! debug 1549 call print_matlu(green%oper_tau(1)%matlu,natom,1) ! debug 1550 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1551 & " == Print green function for small freq after CTQMC" ! debug 1552 call wrtout(std_out,message,'COLL') ! debug 1553 call print_matlu(green%oper(1)%matlu,natom,1) ! debug 1554 end if 1555 1556 1557 ! === Compute rotated Occupations in green%occup_tau 1558 call occup_green_tau(green) 1559 1560 if(pawprtvol>=3) then 1561 ! === Compute non rotated Occupations in green%occup_tau 1562 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the ctqmc basis" 1563 call wrtout(std_out,message,'COLL') 1564 call print_matlu(green%occup_tau%matlu,natom,1) 1565 end if 1566 1567 write(message,'(a,2x,a,f13.5)') ch10,& 1568 & " == Rotate Green function to original basis " 1569 call wrtout(std_out,message,'COLL') 1570 !write(message,'(i3,4x,2e21.14)') 8,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1) 1571 !call wrtout(std_out,message,'COLL') ! debug 1572 1573 ! Rotate oper_tau into Ylm basis and then Slm basis 1574 !------------------------------------------------------------- 1575 do itau=1,paw_dmft%dmftqmc_l 1576 if(opt_diag/=0) call rotate_matlu(green%oper_tau(itau)%matlu,eigvectmatlu,natom,3,0) 1577 if(useylm==1) call slm2ylm_matlu(green%oper_tau(itau)%matlu,natom,2,0) 1578 end do 1579 1580 ! Rotate occup_tau into Ylm basis and then Slm basis 1581 1582 ! Rotate occup_tau into Ylm basis and then Slm basis 1583 !------------------------------------------------------------- 1584 if(opt_diag/=0) call rotate_matlu(green%occup_tau%matlu,eigvectmatlu,natom,3,0) 1585 if(useylm==1) then 1586 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Ylm basis" 1587 call wrtout(std_out,message,'COLL') 1588 call print_matlu(green%occup_tau%matlu,natom,1) 1589 end if 1590 1591 if(useylm==1) call slm2ylm_matlu(green%occup_tau%matlu,natom,2,0) 1592 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Slm basis" 1593 call wrtout(std_out,message,'COLL') 1594 call print_matlu(green%occup_tau%matlu,natom,1) 1595 1596 ! Put Weiss off diagonal terms to zero because Green function will not have any offdiag terms 1597 !------------------------------------------------------------------------------ 1598 ! (if opt_nondiag=0 ie dmft_solv=5) 1599 do ifreq=1,paw_dmft%dmft_nwlo 1600 if(opt_nondiag==0) call zero_matlu(weiss%oper(ifreq)%matlu,natom,onlynondiag=1) 1601 end do 1602 ! ( if opt_nondiag=0, then: 1603 ! As Green's function is diagonal, one suppress off diag terms in Weiss, if any. 1604 ! (If off diag are non zero in the density matrix and thus in the Green's function, 1605 ! there is a warning in checkreal_matlu above).) 1606 1607 ! Rotate Green's and Weiss functions into Ylm basis and then Slm basis 1608 !------------------------------------------------------------- 1609 do ifreq=1,paw_dmft%dmft_nwlo 1610 if(opt_diag/=0) call rotate_matlu(green%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 1611 if(useylm==1) call slm2ylm_matlu(green%oper(ifreq)%matlu,natom,2,0) 1612 if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 1613 if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0) 1614 end do 1615 !write(message,'(i3,4x,2e21.14)') 10,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1) 1616 !call wrtout(std_out,message,'COLL') ! debug 1617 1618 if(pawprtvol>=3) then 1619 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1620 & " == Print green function for small time after rotation (in the original basis)" ! debug 1621 call wrtout(std_out,message,'COLL') ! debug 1622 call print_matlu(green%oper_tau(1)%matlu,natom,1) ! debug 1623 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1624 & " == Print green function for small freq after rotation (in the original basis)" ! debug 1625 call wrtout(std_out,message,'COLL') ! debug 1626 call print_matlu(green%oper(1)%matlu,natom,1) ! debug 1627 !< HACK > 1628 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1629 & " == Print diagonalized weiss_for_rot function after rotation for small freq in the ctqmc basis" ! debug 1630 call wrtout(std_out,message,'COLL') ! debug 1631 call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) ! debug 1632 !</ HACK > 1633 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1634 & " == Print weiss function for small freq in the original basis" ! debug 1635 call wrtout(std_out,message,'COLL') ! debug 1636 call print_matlu(weiss%oper(1)%matlu,natom,1) ! debug 1637 1638 do ifreq=1,paw_dmft%dmft_nwlo 1639 call sym_matlu(cryst_struc,weiss%oper(ifreq)%matlu,pawang,paw_dmft) 1640 end do 1641 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1642 & " == Print symetrized weiss function for small freq in the original basis" ! debug 1643 call wrtout(std_out,message,'COLL') ! debug 1644 call print_matlu(weiss%oper(1)%matlu,natom,1) ! debug 1645 end if 1646 1647 1648 ABI_MALLOC(matlu1,(natom)) 1649 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1) 1650 call copy_matlu(green%occup_tau%matlu,matlu1,natom) 1651 call sym_matlu(cryst_struc,matlu1,pawang,paw_dmft) 1652 1653 write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the original basis" 1654 call wrtout(std_out,message,'COLL') 1655 call print_matlu(green%occup_tau%matlu,natom,1) 1656 1657 write(message,'(a,2x,a,f13.5)') ch10," == Symetrized occupations" 1658 call wrtout(std_out,message,'COLL') 1659 call print_matlu(matlu1,natom,1) 1660 1661 call diff_matlu("CTQMC Occup","CTQMC Occup symetrized",green%occup_tau%matlu,matlu1,natom,0,tol4,ierr) 1662 call destroy_matlu(matlu1,natom) 1663 ABI_FREE(matlu1) 1664 1665 ! ================================================================= 1666 ! Symetrise green function G(tau) and G(ifreq) to recover symetry 1667 ! artificially broken by QMC 1668 ! ================================================================= 1669 write(message,'(a,2x,a,f13.5)') ch10,& 1670 & " == Symetrise green function after QMC " 1671 call wrtout(std_out,message,'COLL') 1672 do itau=1,paw_dmft%dmftqmc_l 1673 call sym_matlu(cryst_struc,green%oper_tau(itau)%matlu,pawang,paw_dmft) 1674 end do 1675 do ifreq=1,paw_dmft%dmft_nwlo 1676 call sym_matlu(cryst_struc,green%oper(ifreq)%matlu,pawang,paw_dmft) 1677 end do 1678 if(pawprtvol>=3) then 1679 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1680 & " == Print green function for small time after symetrisation" ! debug 1681 call wrtout(std_out,message,'COLL') ! debug 1682 call print_matlu(green%oper_tau(1)%matlu,natom,1) ! debug 1683 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 1684 & " == Print green function for small freq after symetrisation" ! debug 1685 call wrtout(std_out,message,'COLL') ! debug 1686 call print_matlu(green%oper(1)%matlu,natom,1) ! debug 1687 end if 1688 if(paw_dmft%dmft_prgn==1) then 1689 call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=2) 1690 call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 1691 end if 1692 1693 ! === Compute Occupations (Symetrized from oper_tau) 1694 call occup_green_tau(green) 1695 1696 1697 ! === Print occupations 1698 call printocc_green(green,6,paw_dmft,3) 1699 1700 call destroy_oper(energy_level) 1701 call destroy_matlu(dmat_diag,natom) 1702 ABI_FREE(dmat_diag) 1703 call destroy_matlu(identity,natom) 1704 ABI_FREE(identity) 1705 do iatom=1,cryst_struc%natom 1706 lpawu=paw_dmft%lpawu(iatom) 1707 if(lpawu/=-1) then 1708 do isppol=1,nsppol 1709 ABI_FREE(eigvectmatlu(iatom,isppol)%value) 1710 !ABI_FREE(udens_atoms(iatom)) 1711 end do 1712 ABI_FREE(udens_atoms(iatom)%value) 1713 end if 1714 end do 1715 ABI_FREE(udens_atoms) 1716 ABI_FREE(eigvectmatlu) 1717 call destroy_green(weiss_for_rot) 1718 ! call destroy_green(gw_loc) 1719 ! call destroy_green(greendft) 1720 1721 ! destroy limit of hybridization 1722 call destroy_matlu(hybri_coeff,paw_dmft%natom) 1723 ABI_FREE(hybri_coeff) 1724 1725 end subroutine qmc_prep_ctqmc
m_forctqmc/testcode_ctqmc [ Functions ]
[ Top ] [ m_forctqmc ] [ Functions ]
NAME
testcode_ctqmc
FUNCTION
Setup ultra simple hybridization to test CTQMC in simple situations.
INPUTS
gtmp_nd gw_tmp_nd temp = temperature dmftqmc_l = number of times slices nflavor = number of flavor testrot = 0/1 if rotation of hybridization is tested or not testcode = 1 if tests are activated. opt = 1/2 if pre or postprocessing of CTQMC data.
OUTPUT
fw1_nd = non diagonal hybridization fw1 = hybridization umod = value of U
SIDE EFFECTS
gtmp_nd gw_tmp_nd
NOTES
SOURCE
1856 subroutine testcode_ctqmc(dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,levels_ctqmc,hybri_limit,& 1857 & nflavor,opt,temp,testrot,testcode,umod) 1858 1859 1860 !Arguments ------------------------------------ 1861 !scalars 1862 integer, intent(in) :: dmftqmc_l,nflavor,testrot,testcode,opt 1863 real(dp), intent(in) :: temp 1864 real(dp), intent(out) :: umod(2,2) 1865 complex(dpc), intent(inout) :: gw_tmp_nd(:,:,:) 1866 real(dp), intent(inout) :: gtmp_nd(:,:,:) 1867 complex(dpc), intent(out) :: fw1(:,:) 1868 complex(dpc), intent(out) :: fw1_nd(:,:,:) 1869 real(dp), intent(inout) :: levels_ctqmc(:) 1870 complex(dpc), intent(inout) :: hybri_limit(:,:) 1871 1872 !Local variables ------------------------------ 1873 character(len=500) :: message 1874 integer :: ifreq, itau,realrot,simplehyb 1875 real(dp) :: omega 1876 real(dp) :: tbi1,tbi2,e2,tbi3,tbi4,e3,e4,tbi21,tbi12,e3b,e4b,tbi21b,tbi12b 1877 complex(dpc) :: e1 1878 ! arrays 1879 complex(dpc) :: RR(2,2) 1880 complex(dpc) :: RR1(2,2) 1881 complex(dpc) :: RRi(2,2) 1882 complex(dpc) :: RRt(2,2) 1883 ! ************************************************************************ 1884 if (testcode==0) return 1885 if (nflavor/=2) then 1886 write(message,'(2a)') ch10,' testcode nflavor.ne.2' 1887 ABI_ERROR(message) 1888 end if 1889 1890 simplehyb=2 1891 simplehyb=1 1892 simplehyb=3 1893 !========================= 1894 ! Built rotation matrix 1895 !========================= 1896 realrot=0 1897 realrot=2 1898 if (realrot==1) then 1899 ! Real rotation 1900 !========================= 1901 RR(1,1) = SQRT(3.d0)/2.d0 1902 RR(1,2) = -1.d0/2.d0 1903 RR(2,1) = 1.d0/2.d0 1904 RR(2,2) = SQRT(3.d0)/2.d0 1905 else if (realrot==2) then 1906 ! Real rotation 1907 !========================= 1908 RR(1,1) = SQRT(1.d0/2.d0) 1909 RR(1,2) = -SQRT(1.d0/2.d0) 1910 RR(2,1) = SQRT(1.d0/2.d0) 1911 RR(2,2) = SQRT(1.d0/2.d0) 1912 else 1913 ! Complex rotation 1914 !========================= 1915 RR(1,1) = CMPLX(one,two) 1916 RR(1,2) = CMPLX(one,one) 1917 RR(2,1) = CMPLX(one,-one) 1918 RR(2,2) = CMPLX(-one,two) 1919 RR=RR/sqrt(seven) 1920 end if 1921 ! Check rotation is unitary 1922 !========================== 1923 RRi(1,1) = conjg(RR(1,1)) 1924 RRi(1,2) = conjg(RR(2,1)) 1925 RRi(2,1) = conjg(RR(1,2)) 1926 RRi(2,2) = conjg(RR(2,2)) 1927 RR1(:,:) = MATMUL ( RR(:,:) , RRi(:,:) ) 1928 !write(6,*) "RR1",RR1 1929 if(abs(RR1(1,1)-one).gt.tol7.or.abs(RR1(1,2)).gt.tol7.or.abs(RR1(2,2)-one).gt.tol7.or.abs(RR1(2,1)).gt.tol7) then 1930 write(message,'(2a)') ch10,' testcode error in rotation matrix' 1931 ABI_ERROR(message) 1932 end if 1933 1934 1935 !================================= 1936 ! Built hybridization for CTQMC 1937 !================================= 1938 if (opt==1) then 1939 1940 ! Parameters: tight-binding + U 1941 ! firt test of the code try umod=0, and (tbi1,tbi2,e1,e2)=(2,1,0.5,0.0) testrot=1 1942 ! second test of the code try umod=four, and (tbi1,tbi2,e1,e2)=(2,1,0.0,0.0) testrot=1 1943 !======================================================================================= 1944 fw1_nd(:,:,:)= czero 1945 tbi1=2.0_dp 1946 tbi2=1.0_dp 1947 tbi3=1.0_dp 1948 tbi4=1.0_dp 1949 tbi12=2.5_dp 1950 tbi12b=2.5_dp 1951 tbi21=2.5_dp 1952 tbi21b=2.5_dp 1953 e1=cmplx(0.0,0.0,8) 1954 e2=zero 1955 e3=0.2 1956 e4=0.3 1957 e3b=0.3 1958 e4b=-0.2 1959 umod(:,:)=0.d0 1960 1961 if(testrot==1.and.(abs(tbi1-tbi2)<tol6)) then 1962 write(message,'(3a)') ch10,' testrot=1 with tbi1=tbi2 is equivalent' & 1963 ,'to testrot=0: change testrot' 1964 ABI_WARNING(message) 1965 end if 1966 ! Built fw1_nd 1967 !============== 1968 do ifreq=1,dmftqmc_l 1969 1970 omega=pi*temp*(two*float(ifreq)-1) 1971 1972 if(simplehyb==1) then 1973 fw1_nd(ifreq,1,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 1974 fw1_nd(ifreq,2,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 1975 fw1(ifreq,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 1976 fw1(ifreq,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 1977 hybri_limit(1,1)=tbi1**2 1978 hybri_limit(2,2)=tbi2**2 1979 hybri_limit(1,2)=0.d0 1980 hybri_limit(2,1)=0.d0 1981 else if(simplehyb==2) then 1982 fw1_nd(ifreq,1,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)+tbi3**2/(dcmplx(0.d0,omega)-e3) 1983 fw1_nd(ifreq,2,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)+tbi4**2/(dcmplx(0.d0,omega)-e4) 1984 fw1(ifreq,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 1985 fw1(ifreq,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 1986 else if(simplehyb==3) then 1987 fw1_nd(ifreq,1,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 1988 fw1_nd(ifreq,2,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 1989 fw1_nd(ifreq,1,2) = tbi12**2/(dcmplx(0.d0,omega)-e3)+tbi12b**2/(dcmplx(0.d0,omega)-e3b) 1990 fw1_nd(ifreq,2,1) = tbi21**2/(dcmplx(0.d0,omega)-e4)+tbi21b**2/(dcmplx(0.d0,omega)-e4b) 1991 fw1(ifreq,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 1992 fw1(ifreq,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 1993 hybri_limit(1,1)=tbi1**2 1994 hybri_limit(2,2)=tbi2**2 1995 hybri_limit(1,2)=tbi12**2+tbi12b**2 1996 hybri_limit(2,1)=tbi21**2+tbi21b**2 1997 end if 1998 write(132,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1)) 1999 write(133,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2)) 2000 write(134,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1)) 2001 write(135,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2)) 2002 write(1234,*) omega, real(fw1(ifreq,1)),aimag(fw1(ifreq,1)) 2003 end do 2004 ! Built level and limit of hybridization 2005 !======================================= 2006 levels_ctqmc(1:nflavor)=-umod(1,1)/two 2007 2008 write(std_out,*) "fw1_nd" 2009 write(std_out,*) fw1_nd(1,1,1), fw1_nd(1,1,2) 2010 write(std_out,*) fw1_nd(1,2,1), fw1_nd(1,2,2) 2011 write(std_out,*) "fw1" 2012 write(std_out,*) fw1(1,1), fw1(1,2) 2013 write(std_out,*) fw1(2,1), fw1(2,2) 2014 2015 ! Rotate hybridization if testrot=1 2016 !================================== 2017 if(testrot==1) then 2018 2019 do ifreq=1,dmftqmc_l 2020 RRt(:,:) = MATMUL ( RR(:,:) , fw1_nd(ifreq,:,:) ) 2021 !write(6,*) "RRt" 2022 !write(6,*) RRt(1,1), RRt(1,2) 2023 !write(6,*) RRt(2,1), RRt(2,2) 2024 RR1(:,:) = MATMUL ( RRt(:,:) , RRi(:,:) ) 2025 !write(6,*) "RR1" 2026 !write(6,*) RR1(1,1), RR1(1,2) 2027 !write(6,*) RR1(2,1), RR1(2,2) 2028 fw1_nd(ifreq,:,:)=RR1(:,:) 2029 omega=pi*temp*(two*float(ifreq)+1) 2030 write(3322,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1)) 2031 write(232,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1)) 2032 write(233,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2)) 2033 write(234,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1)) 2034 write(235,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2)) 2035 end do 2036 2037 ! Rotate limit of hybridization 2038 !======================================= 2039 RRt(:,:) = MATMUL ( RR(:,:) , hybri_limit(:,:) ) 2040 RR1(:,:) = MATMUL ( RRt(:,:) , RRi(:,:) ) 2041 hybri_limit(:,:)=RR1(:,:) 2042 2043 end if 2044 ! rajouter test real(fw1_nd(1,:,:)) doit etre diagonale 2045 2046 !====================================== 2047 ! Rotate Green's function from CTQMC 2048 !====================================== 2049 else if(opt==2) then 2050 2051 write(std_out,*) "gw_tmp_nd" 2052 write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2) 2053 write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2) 2054 ! Rotate Green's function back 2055 !============================== 2056 if(testrot==1) then 2057 do ifreq=1,dmftqmc_l 2058 RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gw_tmp_nd(ifreq,1:nflavor,1:nflavor) ) 2059 RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) ) 2060 gw_tmp_nd(ifreq,1:nflavor,1:nflavor)=RR1(1:nflavor,1:nflavor) 2061 end do 2062 2063 write(std_out,*) "gw_tmp_nd after rotation" 2064 write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2) 2065 write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2) 2066 2067 do itau=1,dmftqmc_l 2068 RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) ) 2069 RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) ) 2070 gtmp_nd(itau,1:nflavor,1:nflavor)=real(RR1(1:nflavor,1:nflavor)) 2071 end do 2072 2073 ! Rotate Green's function for comparison with testrot=1 2074 !====================================================== 2075 else if (testrot==0) then ! produce rotated green's function to compare to testrot=1 case 2076 2077 do itau=1,dmftqmc_l 2078 RRt(1:nflavor,1:nflavor) = MATMUL ( RR(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) ) 2079 RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RRi(1:nflavor,1:nflavor) ) 2080 write(444,*) real(itau-1)/(temp*real(dmftqmc_l)),real(RR1(1,1)),real(RR1(2,2)),real(RR1(1,2)),real(RR1(2,1)) 2081 end do 2082 2083 end if 2084 2085 ! Print out rotated Green's function 2086 !===================================== 2087 do itau=1,dmftqmc_l 2088 write(555,'(e14.5,4(2e14.5,3x))') real(itau-1)/(temp*real(dmftqmc_l)),gtmp_nd(itau,1,1),& 2089 & gtmp_nd(itau,2,2),gtmp_nd(itau,1,2),gtmp_nd(itau,2,1) 2090 end do 2091 2092 write(message,'(2a)') ch10,' testcode end of test calculation' 2093 ABI_ERROR(message) 2094 2095 end if 2096 close(444) 2097 close(555) 2098 2099 end subroutine testcode_ctqmc
m_forctqmc/testcode_ctqmc_b [ Functions ]
[ Top ] [ m_forctqmc ] [ Functions ]
NAME
testcode_ctqmc_b
FUNCTION
Setup ultra simple hybridization to test CTQMC in simple situations.
INPUTS
temp = temperature dmftqmc_l = number of times slices levels_ctqmc_nd=level matrix
OUTPUT
fw1_nd=hybridization matrix umod = value of U hybri_limit= limit of F weiss_for_rot= weiss function hybri_coeff
SIDE EFFECTS
NOTES
SOURCE
1753 subroutine testcode_ctqmc_b(energy_level,hybri_coeff,weiss_for_rot,dmftqmc_l,fw1_nd,levels_ctqmc,& 1754 & levels_ctqmc_nd,hybri_limit,temp,umod,opt_diag,opt_fk) 1755 1756 !Arguments ------------------------------------ 1757 !scalars 1758 integer, intent(in) :: dmftqmc_l,opt_diag,opt_fk 1759 real(dp), intent(in) :: temp 1760 real(dp), intent(out) :: umod(2,2) 1761 real(dp), intent(inout) :: levels_ctqmc(:) 1762 complex(dpc), intent(out) :: fw1_nd(:,:,:) 1763 complex(dpc), intent(inout) :: levels_ctqmc_nd(:,:) 1764 complex(dpc), intent(inout) :: hybri_limit(:,:) 1765 type(oper_type) :: energy_level 1766 type(matlu_type), allocatable :: hybri_coeff(:) 1767 type(green_type) :: weiss_for_rot 1768 1769 !Local variables ------------------------------ 1770 integer :: ifreq,iatom,ima,imb,ispa,ispb 1771 real(dp) :: omega 1772 real(dp) :: facnd, facd 1773 character(len=30) :: tmpfil 1774 ! ************************************************************************ 1775 facnd=0.8d0 1776 facd=1.0d0 1777 !write(6,*) "fac",facnd,facd 1778 levels_ctqmc_nd(2,2) = energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb) 1779 levels_ctqmc_nd(1,1) = energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa) 1780 levels_ctqmc(2) = real(energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb),kind=dp) 1781 levels_ctqmc(1) = real(energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa),kind=dp) 1782 if(opt_diag/=1) then 1783 levels_ctqmc_nd(1,2) = energy_level%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1784 levels_ctqmc_nd(2,1) = energy_level%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1785 end if 1786 hybri_limit(1,1) = facd*hybri_coeff(iatom)%mat(ima,ima,1,ispa,ispa) 1787 hybri_limit(2,2) = facd*hybri_coeff(iatom)%mat(imb,imb,1,ispb,ispb) 1788 hybri_limit(1,2) = facnd*hybri_coeff(iatom)%mat(ima,imb,1,ispa,ispb) 1789 hybri_limit(2,1) = facnd*hybri_coeff(iatom)%mat(imb,ima,1,ispb,ispa) 1790 !write(6,*) "hybri_limit",hybri_limit 1791 !write(6,*) "levels_ctqmc",levels_ctqmc 1792 umod=zero 1793 1794 tmpfil = 'fw1_nd_re' 1795 !if (open_file(newunit=unt,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then 1796 ! ABI_ERROR(message) 1797 !end if 1798 tmpfil = 'fw1_nd_im' 1799 !if (open_file(newunit=unt2,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then 1800 ! ABI_ERROR(message) 1801 !end if 1802 write(std_out,*) "testcode==2",ispa,ispb,ima,imb 1803 write(std_out,*) "opt_fk==",opt_fk 1804 do ifreq=1,dmftqmc_l 1805 if (opt_fk==1) then 1806 fw1_nd(ifreq,1,1) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa) 1807 fw1_nd(ifreq,2,2) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb) 1808 !fw1_nd(ifreq,1,2) = weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1809 !fw1_nd(ifreq,2,1) = weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1810 fw1_nd(ifreq,1,2) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1811 fw1_nd(ifreq,2,1) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1812 omega=pi*temp*(two*float(ifreq)-1) 1813 else if (opt_fk==0) then 1814 fw1_nd(ifreq,1,1) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa) 1815 fw1_nd(ifreq,2,2) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb) 1816 fw1_nd(ifreq,1,2) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb) 1817 fw1_nd(ifreq,2,1) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa) 1818 call xginv(fw1_nd(ifreq,:,:),2) 1819 end if 1820 end do 1821 end subroutine testcode_ctqmc_b