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