TABLE OF CONTENTS


ABINIT/m_entropyDMFT [ Modules ]

[ Top ] [ Modules ]

NAME

  m_entropyDMFT

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

NOTES

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 module m_entropyDMFT
25 
26   use defs_basis
27   use m_errors
28   use m_abicore
29   use m_xmpi
30   use m_dtset
31 
32   use m_energies, only : energies_type, energies_eval_eint
33   use m_splines, only : spline_integrate, spline, splint
34   use m_pawang, only : pawang_type
35   use m_pawrad, only : pawrad_type, simp_gen, poisson
36   use m_pawtab, only : pawtab_type
37   use m_paw_correlations,only : pawpuxinit
38   use m_io_tools, only : get_unit
39   use m_data4entropyDMFT
40 
41   implicit none
42 
43   private
44 
45   public :: entropyDMFT_init
46   public :: entropyDMFT_destroy
47   public :: entropyDMFT_nextLambda
48   public :: entropyDMFT_addIntegrand
49   public :: entropyDMFT_computeEntropy
50 
51   integer, parameter :: E_U0       = 1
52   integer, parameter :: E_UU       = 2
53   integer, parameter :: E_DIRECT   = 1
54   integer, parameter :: E_DC       = 2
55   integer, parameter :: AC_NOTHING = 0
56   integer, parameter :: AC_ETOT    = 1
57   character(len=21), parameter :: HDR_NAME = "DATA FOR ETOT DMFT v="

ABINIT/m_entropyDMFT/entropyDMFT_addIntegrand [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_addIntegrand

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

723   subroutine entropyDMFT_addIntegrand(e_t,dt,energies,data4etot)
724 
725 !Arguments ------------------------------------
726     type(entropyDMFT_t)   , intent(inout) :: e_t
727     type(dataset_type) , intent(in   ) :: dt
728     type(energies_type), intent(in   ) :: energies
729     type(data4entropyDMFT_t)  , intent(in   ) :: data4etot
730 !Local variables ------------------------------
731     integer :: optdc
732     integer :: ictypat
733     integer :: iatom, icatom
734     integer :: iflavor1, iflavor2, ndim
735     integer :: icouple
736 
737     if ( e_t%action == AC_NOTHING ) return
738 
739     ! Write lambda for restart
740     if ( e_t%rank == 0 ) then
741       open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
742       write(e_t%ofile) e_t%lambda(e_t%mylambda)
743       close(e_t%ofile)
744     end if
745 
746     if ( e_t%mylambda == 1 ) then
747       ! Save entropy and internal energy for U=0
748       e_t%entropy0 = energies%entropy
749       ! 1 is for usepaw that is 1 in DMFT, optdc is to know if the DC scheme is
750       ! calculated.
751       call energies_eval_eint(energies,dt,1,optdc,e_t%energies(E_DIRECT,E_U0),e_t%energies(E_DC,E_U0))
752       if ( e_t%rank == 0 ) then
753         open(unit=e_t%ofile,file=e_t%filename,position="append")
754         write(e_t%ofile,'(a)') "# Temperature [Ha]:"
755         write(e_t%ofile,'(es22.14)') e_t%temp
756         write(e_t%ofile,'(a)') "# Entropy for lambda=0 [kb]:"
757         write(e_t%ofile,'(es22.14)') e_t%entropy0
758         write(e_t%ofile,'(a)') "# Internal energy for lambda=0 [Ha]:"
759         write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_U0)
760         close(e_t%ofile)
761         open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
762         write(e_t%ofile) e_t%temp
763         write(e_t%ofile) e_t%entropy0
764         write(e_t%ofile) e_t%energies(E_DC,E_U0)
765         close(e_t%ofile)
766       end if
767     else if ( e_t%mylambda == e_t%nlambda ) then
768       ! Save internal energy for U=Umax
769       call energies_eval_eint(energies,dt,1,optdc,e_t%energies(E_DIRECT,E_UU),e_t%energies(E_DC,E_UU))
770       if ( e_t%rank == 0 ) then
771         open(unit=e_t%ofile,file=e_t%filename,position="append")
772         write(e_t%ofile,'(a)') "# Internal energy for lambda=1 [Ha]:"
773         write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_UU)
774       end if
775       ! Compute U_max,  and uij
776       do ictypat = 1, e_t%nctypat
777         ndim = 2*(2*e_t%lpawu(ictypat)+1)
778         icouple = 0
779         !write(*,*) matU
780         if ( e_t%rank == 0 ) &
781            write(e_t%ofile,'(a,f7.5,1x,a,i4)') "# Interaction Matrix normalized by U=",e_t%U_input(ictypat) , &
782            "[Ha] for atom type", e_t%index_typat(ictypat)
783         do iflavor1 = 1, ndim
784           if ( e_t%rank == 0 ) &
785             write(e_t%ofile,'(14(a21,2x))',advance="no") (/ ( "...", iflavor2=1,iflavor1 ) /)
786           do iflavor2 = iflavor1+1, ndim
787             icouple = icouple + 1
788             e_t%uij(icouple,ictypat) = data4etot%hu_dens(iflavor1,iflavor2,e_t%index_typat(ictypat)) &
789                                       /e_t%U_input(ictypat)   ! to modify, this is just the idea
790             if ( e_t%rank == 0 ) &
791               write(e_t%ofile,'(14(es21.14,2x))',advance="no") e_t%uij(icouple,ictypat)
792           end do
793           if ( e_t%rank == 0 ) write(e_t%ofile,*)
794         end do
795       end do
796 
797       if ( e_t%rank == 0 ) then
798         close(e_t%ofile)
799         open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
800         write(e_t%ofile) e_t%energies(E_DC,E_UU)
801         do ictypat = 1, e_t%nctypat
802           ndim = 2*(2*e_t%lpawu(ictypat)+1)
803           icouple = 0
804           do iflavor1 = 1, ndim
805             do iflavor2 = iflavor1+1, ndim
806               icouple = icouple + 1
807               write(e_t%ofile) e_t%uij(icouple,ictypat)
808             end do
809           end do
810         end do
811         close(e_t%ofile)
812       end if
813     endif
814 
815     ! For all lambda
816     ! Save Docc, Nup and Ndwn
817     if ( e_t%rank == 0 ) then
818       open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
819     end if
820     do icatom = 1, e_t%ncatom
821       iatom = e_t%index_atom(icatom)
822       ndim = 2*(2*e_t%lpawu(e_t%typat(icatom))+1)
823       icouple = 0
824       e_t%e_dc(icatom,e_t%mylambda) = data4etot%e_dc(iatom)
825       if ( e_t%rank == 0 ) then
826         write(e_t%ofile) e_t%e_dc(icatom,e_t%mylambda)
827       end if
828       do iflavor1 = 1, ndim
829         do iflavor2 = iflavor1+1, ndim
830           icouple = icouple + 1
831           e_t%docc(icouple,icatom,e_t%mylambda) = data4etot%Docc(iflavor1,iflavor2,iatom)
832           if ( e_t%rank == 0 ) then
833             write(e_t%ofile) e_t%docc(icouple,icatom,e_t%mylambda)
834           end if
835         end do
836       end do
837     end do
838 
839     if ( e_t%rank == 0 ) then
840       close(e_t%ofile)
841     end if
842 
843     !call entropyDMFT_computeIntegrand(e_t)
844 
845   end subroutine entropyDMFT_addIntegrand

ABINIT/m_entropyDMFT/entropyDMFT_allocateAll [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_allocateAll

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

317   subroutine entropyDMFT_allocateAll(e_t)
318 
319 !Arguments ------------------------------------
320     type(entropyDMFT_t), intent(inout) :: e_t
321 
322     ABI_MALLOC(e_t%index_atom, (1:e_t%natom))
323     e_t%index_atom = 0
324     ABI_MALLOC(e_t%index_typat,(1:e_t%nctypat))
325     e_t%index_atom = 0
326     ABI_MALLOC(e_t%typat,      (1:e_t%ncatom))
327     e_t%typat      = 0
328     ABI_MALLOC(e_t%lpawu,      (1:e_t%nctypat))
329     e_t%lpawu      = 0
330     ABI_MALLOC(e_t%U_input,    (1:e_t%nctypat))
331     e_t%U_input    = zero
332     ABI_MALLOC(e_t%J_input,    (1:e_t%nctypat))
333     e_t%J_input    = zero
334     ABI_MALLOC(e_t%lambda,     (1:e_t%nlambda))
335     e_t%lambda     = 0
336     ABI_MALLOC(e_t%docc,       (1:(14*13)/2,1:e_t%ncatom,1:e_t%nlambda))
337     e_t%docc       = zero
338     ABI_MALLOC(e_t%e_dc,       (1:e_t%ncatom,1:e_t%nlambda))
339     e_t%e_dc       = zero
340     ABI_MALLOC(e_t%uij,        (1:(14*13)/2,1:e_t%nctypat))
341     e_t%uij        = zero
342   end subroutine entropyDMFT_allocateAll

ABINIT/m_entropyDMFT/entropyDMFT_computeEntropy [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_computeEntropy

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

873   subroutine entropyDMFT_computeEntropy(e_t,entropy)
874 
875 !Arguments ------------------------------------
876     type(entropyDMFT_t), intent(inout) :: e_t
877     real(dp)        , intent(inout) :: entropy !inout in case we do nothing, avoid to change entropy in NaN
878 !Local variables ------------------------------
879     integer :: ilambda
880     integer :: nlambda
881     integer :: icatom
882     integer :: ncatom
883     integer :: icouple
884     integer :: ndim
885     real(dp), allocatable :: integrand(:,:)
886     real(dp) :: docc
887     real(dp) :: integral
888     real(dp) :: entropyDirect
889     character(len=500) :: string
890     character(len=50) :: i2str
891 
892     if ( e_t%action == AC_NOTHING ) return
893 
894     nlambda = e_t%nlambda
895     ncatom  = e_t%ncatom
896 
897     ABI_MALLOC(integrand,(1:nlambda,1:ncatom))
898     integrand(1:nlambda,1:ncatom) = zero
899 
900     if ( e_t%rank == 0 ) then
901       open(unit=e_t%ofile,file=e_t%filename,position="append")
902       write(e_t%ofile,'(2a)') ch10,"#Decomposition for each lambda per atom"
903     end if
904     do ilambda = 1, nlambda
905       if ( e_t%rank == 0 ) then
906         write(e_t%ofile,'(a,i2,a)') "# ------| lambda ",ilambda," |------ #"
907         write(e_t%ofile,'(a5,a23,a23)') "#Atom","DC","<uij*n_i*n_j>"
908       end if
909       do icatom = 1, ncatom
910         integrand(ilambda,icatom) = -e_t%e_dc(icatom,ilambda)
911         ndim = 2*(2*e_t%lpawu(e_t%typat(e_t%index_atom(icatom)))+1) ! nflavors
912         ndim = ndim*(ndim-1)/2 ! number of couples
913         docc = zero
914         do icouple = 1, ndim
915           docc = docc + e_t%uij(icouple,e_t%typat(icatom)) * e_t%docc(icouple,icatom,ilambda)
916         end do
917         integrand(ilambda,icatom) = integrand(ilambda,icatom) + docc
918         if ( e_t%rank == 0 ) &
919           write(e_t%ofile,'(i5,2es23.14)') e_t%index_atom(icatom),e_t%e_dc(icatom,ilambda),docc
920       end do
921     end do
922 
923     ! print for test purpose
924     if ( e_t%rank == 0 ) then
925       write(i2str,'(i6)') e_t%ncatom
926       string='(a,a7,'//TRIM(ADJUSTL(i2str))//'(12x,"Atom",i6))'
927       !open(unit=e_t%ofile,file=e_t%filename,position="append")
928       write(e_t%ofile,string) ch10,"#lambda", (/ (e_t%index_atom(icatom),icatom=1,ncatom) /)
929       string='(1x,f6.4,'//TRIM(ADJUSTL(i2str))//'(1x,e21.14))'
930       do ilambda = 1, nlambda
931         write(e_t%ofile,string) e_t%lambda(ilambda), (/ (integrand(ilambda,icatom),icatom=1,ncatom) /)
932       end do
933       close(e_t%ofile)
934     end if
935 
936     ! Integrate the integrand for all correlated atoms
937     call entropyDMFT_integrate(e_t,integrand,integral)
938     ABI_FREE(integrand)
939 
940     write(string,'(a,1x,78a)') ch10,"+",(/ ("-",ilambda=1,76) /), "+"
941     call wrtout(std_out,string,"COLL")
942     call wrtout(ab_out,string,"COLL")
943     write(string,'(1x,a)') "|             Calculation of entropy within  the DMFT Framework              |"
944     call wrtout(std_out,string,"COLL")
945     call wrtout(ab_out,string,"COLL")
946     write(string,'(1x,40a)') "+",(/ ("- ",ilambda=1,38) /), "+"
947     call wrtout(std_out,string,"COLL")
948 
949     entropyDirect = e_t%entropy0 + ( e_t%energies(E_DIRECT,E_UU) - e_t%energies(E_DIRECT,E_U0) - integral ) /e_t%temp
950     entropy = e_t%entropy0 + ( e_t%energies(E_DC,E_UU) - e_t%energies(E_DC,E_U0) - integral ) /e_t%temp
951 
952     write(i2str,'(a)') '(1x,a,19x,a11,es21.14,1x,a4,20x,a)'
953     write(string,i2str) "|","Integral = ", integral, "[Ha]", "|"
954     call wrtout(std_out,string,"COLL")
955     write(string,i2str) "|","E(0)     = ", e_t%energies(E_DC,E_U0), "[Ha]", "|"
956     call wrtout(std_out,string,"COLL")
957     write(string,i2str) "|","E(U)     = ", e_t%energies(E_DC,E_UU), "[Ha]", "|"
958     call wrtout(std_out,string,"COLL")
959     write(string,i2str) "|","S(0)     = ", e_t%entropy0, "[kb]", "|"
960     call wrtout(std_out,string,"COLL")
961     write(string,i2str) "|","S(U)     = ", entropy, "[kb]", "|"
962     call wrtout(std_out,string,"COLL")
963     write(string,'(1x,40a)') "+",(/ ("- ",ilambda=1,38) /), "+"
964     call wrtout(std_out,string,"COLL")
965 
966     write(string,'(1x,a,16x,a22,es21.14,1x,a4,12x,a)') "|","-kT*entropy is set to ", -e_t%temp*entropy, "[Ha]", "|"
967     call wrtout(std_out,string,"COLL")
968     call wrtout(ab_out,string,"COLL")
969 
970     write(string,'(1x,78a)') "+",(/ ("-",ilambda=1,76) /), "+"
971     call wrtout(std_out,string,"COLL")
972     call wrtout(ab_out,string,"COLL")
973 
974     if ( entropy < zero ) then
975       write(string,'(3a)') "Entropy is negative !!!!",ch10,&
976       "It does not make any sense"
977       ABI_WARNING(string)
978     end if
979 
980     if ( abs(entropy-entropyDirect) >= tol3 ) then
981       write(string,'(1x,a,1x,f8.3,1x,a,1x,f8.3,2a)') "Difference between Direct and DC entropies is", abs(entropy-entropyDirect), &
982         "which is greater than", tol3,ch10,"Action : converge better the DMFT and/or DFT loops"
983       ABI_WARNING(string)
984     end if
985 
986 
987 
988   end subroutine entropyDMFT_computeEntropy

ABINIT/m_entropyDMFT/entropyDMFT_destroy [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_destroy

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

1107   subroutine entropyDMFT_destroy(e_t)
1108 
1109 !Arguments ------------------------------------
1110     type(entropyDMFT_t), intent(inout) :: e_t
1111 
1112     if ( allocated(e_t%index_atom) ) then
1113       ABI_FREE(e_t%index_atom)
1114     endif
1115     if ( allocated(e_t%index_typat)) then
1116       ABI_FREE(e_t%index_typat)
1117     endif
1118     if ( allocated(e_t%typat     ) ) then
1119       ABI_FREE(e_t%typat     )
1120     endif
1121     if ( allocated(e_t%lpawu     ) ) then
1122       ABI_FREE(e_t%lpawu)
1123     endif
1124     if ( allocated(e_t%U_input   ) ) then
1125       ABI_FREE(e_t%U_input)
1126     endif
1127     if ( allocated(e_t%J_input   ) ) then
1128       ABI_FREE(e_t%J_input)
1129     endif
1130     if ( allocated(e_t%lambda    ) ) then
1131       ABI_FREE(e_t%lambda)
1132     endif
1133     if ( allocated(e_t%docc      ) ) then
1134       ABI_FREE(e_t%docc)
1135     endif
1136     if ( allocated(e_t%e_dc      ) ) then
1137       ABI_FREE(e_t%e_dc)
1138     endif
1139     if ( allocated(e_t%uij       ) ) then
1140       ABI_FREE(e_t%uij)
1141     endif
1142     e_t%isset = .FALSE.
1143   end subroutine entropyDMFT_destroy

ABINIT/m_entropyDMFT/entropyDMFT_dump [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_dump

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

512   subroutine entropyDMFT_dump(e_t)
513 
514 !Arguments ------------------------------------
515     type(entropyDMFT_t), intent(inout) :: e_t
516 !Local variables-------------------------------
517     integer            :: ilambda
518     integer            :: ictypat
519     integer            :: ndim
520     integer            :: icouple
521     integer            :: iflavor1
522     integer            :: iflavor2
523     integer            :: icatom
524     integer            :: iatom
525 
526     if ( e_t%rank /= 0 ) return
527 
528     ! Dump _data4EtotDMFT
529     open(unit=e_t%ofile,file=e_t%filedata,action="write",form="unformatted")
530     write(e_t%ofile) HDR_NAME, 1
531 
532     do ilambda = 1, e_t%mylambda
533       write(e_t%ofile) e_t%lambda(ilambda)
534 
535       if ( ilambda == 1 ) then
536         write(e_t%ofile) e_t%temp
537         write(e_t%ofile) e_t%entropy0
538         write(e_t%ofile) e_t%energies(E_DC,E_U0)
539       else if ( ilambda == e_t%nlambda ) then ! should never happend ?!
540         write(e_t%ofile) e_t%energies(E_DC,E_UU)
541         do ictypat = 1, e_t%nctypat
542           ndim = 2*(2*e_t%lpawu(ictypat)+1)
543           icouple = 0
544           do iflavor1 = 1, ndim
545             do iflavor2 = iflavor1+1, ndim
546               icouple = icouple + 1
547               write(e_t%ofile) e_t%uij(icouple,ictypat)
548             end do
549           end do
550         end do
551       end if
552 
553       do icatom = 1, e_t%ncatom
554         iatom = e_t%index_atom(icatom)
555         ndim = 2*(2*e_t%lpawu(e_t%typat(icatom))+1)
556         icouple = 0
557         write(e_t%ofile) e_t%e_dc(icatom,ilambda)
558         do iflavor1 = 1, ndim
559           do iflavor2 = iflavor1+1, ndim
560             icouple = icouple + 1
561             write(e_t%ofile) e_t%docc(icouple,icatom,ilambda)
562           end do
563         end do
564       end do
565     end do
566     close(e_t%ofile)
567 
568     ! Dump _EtotDMFT
569     open(unit=e_t%ofile,file=e_t%filename)
570     write(e_t%ofile,'(2a)') "# Data for entropy calculation in DMFT",ch10
571     do ilambda = 1, e_t%mylambda
572       if ( ilambda == 1 ) then
573         write(e_t%ofile,'(a)') "# Temperature [Ha]:"
574         write(e_t%ofile,'(es22.14)') e_t%temp
575         write(e_t%ofile,'(a)') "# Entropy for lambda=0 [kb]:"
576         write(e_t%ofile,'(es22.14)') e_t%entropy0
577         write(e_t%ofile,'(a)') "# Internal energy for lambda=0 [Ha]:"
578         write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_U0)
579       else if ( e_t%mylambda == e_t%nlambda ) then
580         write(e_t%ofile,'(a)') "# Internal energy for lambda=1 [Ha]:"
581         write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_UU)
582         do ictypat = 1, e_t%nctypat
583           ndim = 2*(2*e_t%lpawu(ictypat)+1)
584           icouple = 0
585           write(e_t%ofile,'(a,f7.5,1x,a,i4)') "# Interaction Matrix normalized by U=",e_t%U_input(ictypat) , &
586           "[Ha] for atom type", e_t%index_typat(ictypat)
587           do iflavor1 = 1, ndim
588             write(e_t%ofile,'(14(a21,2x))',advance="no") (/ ( "...", iflavor2=1,iflavor1 ) /)
589             do iflavor2 = iflavor1+1, ndim
590               icouple = icouple + 1
591               write(e_t%ofile,'(14(es21.14,2x))',advance="no") e_t%uij(icouple,ictypat)
592             end do
593             write(e_t%ofile,*)
594           end do
595         end do
596       end if
597     end do
598     close(e_t%ofile)
599   end subroutine entropyDMFT_dump

ABINIT/m_entropyDMFT/entropyDMFT_init [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_init

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

138 subroutine entropyDMFT_init(e_t,dt,pawtab,spacecomm,ifilename,ofilename)
139 
140 !Arguments ------------------------------------
141     type(entropyDMFT_t) , intent(inout) :: e_t
142     type(dataset_type)  , intent(in   ) :: dt
143     type(pawtab_type)   , intent(in   ) :: pawtab(:)
144     integer             , intent(in   ) :: spacecomm
145     character(len=fnlen), intent(in   ) :: ifilename
146     character(len=fnlen), intent(in   ) :: ofilename
147 !Local variables ------------------------------
148     logical :: doRestart
149     integer :: natom, ncatom
150     integer :: iatom, icatom
151     integer :: nctypat,itypat, ictypat
152     integer :: ilambda
153     integer, allocatable :: maptypat(:)
154     character(len=500) :: message
155 
156     e_t%action = dt%dmft_entropy
157     e_t%mylambda = e_t%action-1   ! should be usefull to start form a given value of lambda
158                                   ! -1 added to start from 1 at the first call
159                                   ! of nextLambda
160 
161     if ( e_t%action == AC_NOTHING ) then
162       e_t%isset = .TRUE.
163       return
164     endif
165 
166     e_t%action = AC_ETOT
167 
168     natom = dt%natom
169     ncatom = dt%natpawu
170 
171     icatom = 0
172     do iatom=1,natom
173       if ( pawtab(dt%typat(iatom))%lpawu /= -1 ) icatom = icatom + 1
174     end do
175 
176     if ( icatom /= ncatom ) ABI_ERROR("Inconsistent number of correlated atoms")
177 
178     nctypat = 0
179     do itypat=1,dt%ntypat
180       if ( pawtab(itypat)%lpawu /= -1 ) nctypat = nctypat + 1
181     end do
182 
183     e_t%natom  = natom
184     e_t%ncatom = ncatom
185     e_t%ntypat = dt%ntypat
186     e_t%nctypat = nctypat
187 
188     if ( dt%dmft_nlambda < 3 ) then
189       write(message,'(2a,i4,2a)') "DMFT must have dmft_nlamda >= 3 to compute entropy", &
190         "whereas its value is dmft_nlambda = ",dt%dmft_nlambda,ch10,&
191         "Action : check you input variable dmft_nlambda"
192       ABI_ERROR(message)
193     end if
194 
195     e_t%nlambda = dt%dmft_nlambda
196 
197     doRestart = .false.
198     if ( e_t%nlambda < (e_t%mylambda+1) ) then
199       ABI_ERROR("Restart calculation of DMFT entropy with a value of dmft_entropy greater than dmft_nlambda")
200     else if ( e_t%mylambda > 0 ) then
201       doRestart = .true.
202     end if
203 
204     call entropyDMFT_allocateAll(e_t)
205 
206     e_t%entropy0      = zero
207     e_t%energies(:,:) = zero
208 
209     e_t%lambda(:)     = (/ (DBLE(ilambda-1)/DBLE(e_t%nlambda-1), ilambda=1,e_t%nlambda) /)
210 
211     e_t%temp          = dt%tsmear
212 
213     ! Save each value of U and J for each correlated type
214     ictypat = 1
215     ABI_MALLOC(maptypat,(1:dt%ntypat))
216     do itypat=1,dt%ntypat
217       if ( pawtab(itypat)%lpawu /= -1 ) then
218         e_t%lpawu(ictypat) = pawtab(itypat)%lpawu
219         if(dt%dmft_t2g==1.and.e_t%lpawu(ictypat)==2) then
220           e_t%lpawu(ictypat)=1
221         end if
222         e_t%U_input(ictypat) = pawtab(itypat)%upawu
223         e_t%J_input(ictypat) = pawtab(itypat)%jpawu
224         e_t%index_typat(ictypat) = itypat
225         maptypat(itypat) = ictypat
226         ictypat = ictypat + 1
227       end if
228     end do
229 
230     ! Save type local and global of each correlated atom
231     ! Get the correct value of lpawu for each type
232     icatom = 1
233     do iatom=1,e_t%ncatom
234       if ( pawtab(dt%typat(iatom))%lpawu /= -1 ) then
235         e_t%typat(icatom) = maptypat(dt%typat(iatom))
236         e_t%index_atom(icatom) = iatom
237         icatom = icatom + 1
238       end if
239     end do
240     ABI_FREE(maptypat)
241 
242     write(message,'(a,1x,78a)') ch10,"+",(/ ("-",ilambda=1,76) /), "+"
243     call wrtout(std_out,message,"COLL")
244     call wrtout(ab_out,message,"COLL")
245     write(message,'(1x,a)') "|             Calculation of entropy within  the DMFT Framework              |"
246     call wrtout(std_out,message,"COLL")
247     call wrtout(ab_out,message,"COLL")
248     write(message,'(1x,40a)') "+",(/ ("- ",ilambda=1,38) /), "+"
249     call wrtout(std_out,message,"COLL")
250     do ilambda = 1, e_t%nlambda
251       write(message,'(1x,a,i4,a11,f6.4,55x,a)') &
252         "|", ilambda, ") lambda = ", e_t%lambda(ilambda), "|"
253       call wrtout(std_out,message,"COLL")
254       do ictypat=1,e_t%nctypat
255         write(message,'(1x,a,6x,a12,i4,4x,a3,5x,2(3x,a4,f6.4,1x,a2),10x,a)') "|", &
256           "- Atom type ", e_t%index_typat(ictypat), "->", &
257           "U = ",e_t%U_input(ictypat)*e_t%lambda(ilambda), "Ha", &
258           "J = ",e_t%J_input(ictypat)*e_t%lambda(ilambda), "Ha","|"
259         call wrtout(std_out,message,"COLL")
260       end do
261     end do
262     write(message,'(1x,78a)') "+",(/ ("-",ilambda=1,76) /), "+"
263     call wrtout(std_out,message,"COLL")
264     call wrtout(ab_out,message,"COLL")
265 
266     ! Set up MPI
267     e_t%spacecomm = spacecomm
268     e_t%rank      = xmpi_comm_rank(spacecomm)
269     e_t%comm_size = xmpi_comm_size(spacecomm)
270 
271     e_t%ofile = get_unit()
272     e_t%ofilename = ofilename
273     e_t%ifilename = ifilename
274     e_t%filename = TRIM(e_t%ofilename)//"_EntropyDMFT"
275     e_t%filedata = TRIM(e_t%ofilename)//"_data4EntropyDMFT"
276 
277     if ( doRestart .eqv. .true. ) then
278       ! If restart fails, then nothing changes and the full calculation is
279       ! perform. Otherwise, we complete as much a possible the structure.
280       call entropyDMFT_restart(e_t)
281     end if
282 
283     ! Rewrite the files with the previous data
284     call entropyDMFT_dump(e_t)
285 
286     e_t%isset = .TRUE.
287 
288 
289   end subroutine entropyDMFT_init

ABINIT/m_entropyDMFT/entropyDMFT_integrate [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_integrate

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

1016   subroutine entropyDMFT_integrate(e_t,integrand,integral)
1017 
1018 !Arguments ------------------------------------
1019     type(entropyDMFT_t), intent(inout) :: e_t
1020     real(dp)        , intent(in   ) :: integrand(:,:)
1021     real(dp)        , intent(  out) :: integral
1022 !Local variables ------------------------------
1023     real(dp), allocatable :: ypp(:)
1024     real(dp), allocatable :: fitx(:)
1025     real(dp), allocatable :: fity(:)
1026     real(dp) :: integral1
1027     !real(dp) :: integral2
1028     !real(dp) :: dx
1029     integer :: unit
1030     integer :: Nfit
1031     integer :: icatom
1032     integer :: i
1033 
1034     Nfit = 1000
1035 
1036     ABI_MALLOC(ypp,(1:e_t%nlambda))
1037     ABI_MALLOC(fitx,(1:Nfit))
1038     ABI_MALLOC(fity,(1:Nfit))
1039 
1040     integral = zero
1041 
1042     unit = get_unit()
1043     if ( e_t%rank == 0 ) open(unit=unit,file="fit.log")
1044     do icatom = 1, e_t%ncatom
1045       integral1 = zero
1046       !integral2 = zero
1047       !dx = e_t%U_input(icatom)/dble(e_t%nlambda-1)
1048 
1049       call spline(e_t%lambda,integrand(:,icatom),e_t%nlambda,&
1050         (integrand(2,icatom)-integrand(1,icatom))/e_t%lambda(2),&   ! first derivative left
1051         (integrand(e_t%nlambda,icatom)-integrand(e_t%nlambda-1,icatom))/e_t%lambda(2),& ! first derivative right
1052         ypp)
1053       fitx(1:Nfit)=(/ (dble(i-1)/dble(Nfit-1),i=1,Nfit) /)
1054       call splint(e_t%nlambda,e_t%lambda,integrand(:,icatom),ypp,Nfit,fitx,fity)
1055       integral1 = (sum(fity(1:Nfit))-(fity(1)+fity(Nfit))*half)*e_t%U_input(e_t%typat(icatom))/dble(Nfit-1)
1056       if ( e_t%rank == 0 ) then
1057         do i=1,Nfit
1058           write(unit,'(2ES21.14)') fitx(i),fity(i)
1059         end do
1060       end if
1061       !Unfortunately I don't trust this function and I don't have time to understand it.
1062       !call spline_integrate(integral2,Nfit,dx,integrand(:,icatom))
1063       !write(*,*) integral1, integral2
1064 
1065       !if ( abs(integral2-integral1) >= tol6 ) then
1066       !  write(msg,'(1x,a,1x,f8.6,1x,a,1x,f8.6,1x,a,i4)') "Difference between two different ways of integration is", abs(integral2-integral1), &
1067       !    "which is greater than", tol6, "for correlated atom",e_t%index_atom(icatom)
1068       !  ABI_WARNING(msg)
1069       !end if
1070 
1071       integral = integral+integral1
1072     end do
1073     if ( e_t%rank == 0 ) close(unit)
1074 
1075     ABI_FREE(ypp)
1076     ABI_FREE(fitx)
1077     ABI_FREE(fity)
1078 
1079   end subroutine entropyDMFT_integrate

ABINIT/m_entropyDMFT/entropyDMFT_nextLambda [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_nextLambda

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

627   function entropyDMFT_nextLambda(e_t,dt,pawtab,pawang,pawrad) result(nextstep)
628 
629 !Arguments ------------------------------------
630     type(entropyDMFT_t) , intent(inout) :: e_t
631     type(dataset_type)  , intent(in) :: dt
632     type(pawtab_type)   , intent(inout) :: pawtab(:)
633     type(pawang_type)   , intent(in   ) :: pawang
634     type(pawrad_type)   , intent(inout) :: pawrad(:)
635 !Local variables ------------------------------
636     logical :: nextstep
637     integer :: itypat
638     integer :: mylambda
639     logical :: is_dfpt=.false.
640     real(dp),allocatable :: upawu(:),jpawu(:)
641     character(len=100) :: message
642 
643     if ( e_t%isset .eqv. .FALSE. ) &
644       ABI_ERROR("entropyDMFT is not initialized")
645 
646     ! go to next lambda
647     mylambda = e_t%mylambda + 1
648     !if ( present(ilambda) ) mylambda = ilambda
649     e_t%mylambda = mylambda
650 
651     if ( e_t%action == AC_NOTHING .and. mylambda >= 1 ) then
652       nextstep = .FALSE.
653       ! We do nothing and return false, one scfvc has already performed
654     else if ( e_t%action == AC_NOTHING .and. mylambda < 1 ) then
655       nextstep = .TRUE.
656       ! We do nothing but return true to perform scfcv a least once
657     else if ( e_t%action == AC_ETOT .and. mylambda <= e_t%nlambda ) then ! we iterate over lambda
658       nextstep = .TRUE.
659       ABI_MALLOC(upawu,(dt%ntypat))
660       ABI_MALLOC(jpawu,(dt%ntypat))
661       do itypat = 1, e_t%nctypat
662         upawu(e_t%index_typat(itypat)) = e_t%lambda(mylambda) * e_t%U_input(itypat)
663         jpawu(e_t%index_typat(itypat)) = e_t%lambda(mylambda) * e_t%J_input(itypat)
664       end do
665     else ! we did all lambda values
666       nextstep = .FALSE.
667     endif
668 
669     if ( e_t%action == AC_ETOT .and. nextstep .eqv. .true. ) then
670       write(message,'(a,1x,78a)') ch10,"+",(/ ("-",itypat=1,76) /), "+"
671       call wrtout(std_out,message,"COLL")
672       call wrtout(ab_out,message,"COLL")
673       write(message,'(1x,a,i4,a11,f6.4,55x,a)') &
674         "|", mylambda, ") lambda = ", e_t%lambda(mylambda), "|"
675       call wrtout(std_out,message,"COLL")
676       call wrtout(ab_out,message,"COLL")
677       do itypat=1,e_t%nctypat
678         write(message,'(1x,a,6x,a12,i4,4x,a3,5x,2(3x,a4,f6.4,1x,a2),10x,a)') "|",&
679           "- Atom type ", e_t%index_typat(itypat), "->", &
680           "U = ",e_t%U_input(itypat)*e_t%lambda(mylambda), "Ha", &
681           "J = ",e_t%J_input(itypat)*e_t%lambda(mylambda), "Ha","|"
682           call wrtout(std_out,message,"COLL")
683           call wrtout(ab_out,message,"COLL")
684       end do
685       write(message,'(1x,78a)') "+",(/ ("-",itypat=1,76) /), "+"
686       call wrtout(std_out,message,"COLL")
687       call wrtout(ab_out,message,"COLL")
688       call pawpuxinit(dt%dmatpuopt,dt%exchmix,dt%f4of2_sla,dt%f6of2_sla,&
689 &        is_dfpt,jpawu,dt%lexexch,dt%lpawu,dt%nspinor,dt%ntypat,dt%optdcmagpawu,pawang,dt%pawprtvol,&
690 &        pawrad,pawtab,upawu,dt%usedmft,dt%useexexch,dt%usepawu)
691       ABI_FREE(upawu)
692       ABI_FREE(jpawu)
693     end if
694 
695   end function entropyDMFT_nextLambda

ABINIT/m_entropyDMFT/entropyDMFT_restart [ Functions ]

[ Top ] [ Functions ]

NAME

  entropyDMFT_restart

FUNCTION

  FIXME: add description.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  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

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

370   subroutine entropyDMFT_restart(e_t)
371 
372 !Arguments ------------------------------------
373     type(entropyDMFT_t), intent(inout) :: e_t
374 !Local variables-------------------------------
375     logical            :: doBcast
376     character(len=200) :: msg
377     character(len=21 ) :: hdr
378     integer            :: iostate
379     integer            :: ilambda
380     integer            :: tlambda
381     integer            :: ictypat
382     integer            :: ndim
383     integer            :: icouple
384     integer            :: iflavor1
385     integer            :: iflavor2
386     integer            :: icatom
387     integer            :: iatom
388     real(dp)           :: lambda
389 
390     write(msg,'(a,1x,2a,1x,a,f6.4,a)') ch10,"EtotDMFT will try to restart calculation from a previous run", ch10, &
391     "and this calculation should start at lambda = ",e_t%lambda(e_t%mylambda+1), ch10
392     call wrtout(std_out,msg,"COLL")
393     call wrtout(ab_out,msg,"COLL")
394 
395     doBcast = .true.
396 
397     if ( e_t%rank == 0 ) then
398       tlambda = 1
399       inquire(file=e_t%filedata,iostat=iostate)
400 
401       if ( iostate /= 0 ) then
402         write(msg,'(5a)') "File ", trim(e_t%filedata), " does not exist or is not accessible", ch10, &
403         "-> No restart performed but full calculation."
404         ABI_WARNING(msg)
405         e_t%mylambda = 0
406       end if
407 
408       open(unit=e_t%ofile,file=e_t%filedata,action="read",form="unformatted")
409       read(e_t%ofile,end=42) hdr, iostate
410 
411       if ( hdr /= HDR_NAME .or. iostate /= 1 ) then
412         write(msg,'(5a)') "File ", trim(e_t%filedata), " does not contain the proper header", ch10, &
413         "-> No restart performed but full calculation."
414         ABI_WARNING(msg)
415         e_t%mylambda = 0
416       end if
417 
418       do ilambda = 1, e_t%mylambda
419         tlambda = ilambda
420         read(e_t%ofile,end=42) lambda
421         if ( ABS(lambda - e_t%lambda(ilambda)) >= tol9 ) then
422           write(msg,'(5a,f6.4,a,f6.4)') "File ", trim(e_t%filedata), " is wrong:", ch10, &
423           "Lambda values are differente: in file ", lambda, " instead of ", e_t%lambda(ilambda)
424           ABI_WARNING(msg)
425           goto 42
426         end if
427 
428         if ( ilambda == 1 ) then
429           read(e_t%ofile,end=42) lambda
430           if ( lambda /= e_t%temp ) then
431             write(msg,'(7a,f6.4)') "File ", trim(e_t%filedata), " is wrong:", ch10, &
432             "Temperature is different than the value of tsmear", ch10, &
433             "-> No restart performed but full calculation."
434             ABI_WARNING(msg)
435             goto 42
436           end if
437           read(e_t%ofile,end=42) e_t%entropy0
438           read(e_t%ofile,end=42) e_t%energies(E_DC,E_U0)
439         else if ( ilambda == e_t%nlambda ) then ! should never happend ?!
440           read(e_t%ofile,end=42) e_t%energies(E_DC,E_UU)
441           do ictypat = 1, e_t%nctypat
442             ndim = 2*(2*e_t%lpawu(ictypat)+1)
443             icouple = 0
444             do iflavor1 = 1, ndim
445               do iflavor2 = iflavor1+1, ndim
446                 icouple = icouple + 1
447                 read(e_t%ofile,end=42) e_t%uij(icouple,ictypat)
448               end do
449             end do
450           end do
451         end if
452 
453         do icatom = 1, e_t%ncatom
454           iatom = e_t%index_atom(icatom)
455           ndim = 2*(2*e_t%lpawu(e_t%typat(icatom))+1)
456           icouple = 0
457           read(e_t%ofile,end=42) e_t%e_dc(icatom,ilambda)
458           do iflavor1 = 1, ndim
459             do iflavor2 = iflavor1+1, ndim
460               icouple = icouple + 1
461               read(e_t%ofile,end=42) e_t%docc(icouple,icatom,ilambda)
462             end do
463           end do
464         end do
465       end do
466       close(e_t%ofile)
467       goto 43
468 42    write(msg,'(5a,f6.4)') "File ", trim(e_t%filedata), " is wrong or incomplete", ch10, &
469           "-> Restart calculation will restart at lambda = ",e_t%lambda(tlambda)
470       ABI_WARNING(msg)
471       close(e_t%ofile)
472       e_t%mylambda = tlambda-1 ! -1 to go to previous lambda
473     end if
474     ! MPI BDCAST
475 43  call xmpi_bcast(e_t%mylambda,0, e_t%spacecomm, ictypat)
476     call xmpi_bcast(e_t%entropy0,0, e_t%spacecomm, ictypat)
477     call xmpi_bcast(e_t%energies,0, e_t%spacecomm, ictypat)
478     !call xmpi_bcast(e_t%uij,0, e_t%spacecomm, ictypat) ! No need since it
479     ! restart always perform the last lambda and this is calculation at the end
480     ! of last lambda
481     call xmpi_bcast(e_t%e_dc,0, e_t%spacecomm, ictypat)
482     call xmpi_bcast(e_t%docc,0, e_t%spacecomm, ictypat)
483 
484   end subroutine entropyDMFT_restart

m_entropyDMFT/entropyDMFT [ Types ]

[ Top ] [ m_entropyDMFT ] [ Types ]

NAME

  entropyDMFT

FUNCTION

  This structured datatype contains the necessary data

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (J. Bieder)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

 76   type, public :: entropyDMFT_t
 77     logical               :: isset=.FALSE.  ! flag to be sure we are initialize
 78     integer               :: spacecomm      ! MPI comm
 79     integer               :: rank           ! rank in the comm
 80     integer               :: comm_size      ! Number of cpus in the comm
 81     integer               :: action         ! what to do in gstate
 82     integer               :: mylambda       ! Current lambda
 83     integer               :: natom          ! number of atoms
 84     integer               :: ncatom         ! number of correlated atoms
 85     integer               :: ntypat         ! number of type of atoms
 86     integer               :: nctypat        ! number of type of correlated atoms
 87     integer               :: nlambda        ! number of integration points
 88     integer               :: ofile          ! unit for file output data
 89     character(len=fnlen)  :: filename       ! name fo the file user readable
 90     character(len=fnlen)  :: filedata       ! name fo the file for restart purposes
 91     character(len=fnlen)  :: ofilename      ! ofstream prefix
 92     character(len=fnlen)  :: ifilename      ! ifstream prefix
 93     real(dp)              :: temp           ! temperature
 94     real(dp)              :: entropy0       ! entropy for lambda=0
 95     real(dp)              :: energies(2,2)  ! internal energy for lambda=0 and 1
 96     integer , allocatable :: index_atom(:)  ! index for correlated atoms
 97     integer , allocatable :: index_typat(:) ! index for correlated types
 98     integer , allocatable :: lpawu(:)       ! orbital moment to treat (ntypat,2)
 99     integer , allocatable :: typat(:)       ! type of each correlated atom
100     real(dp), allocatable :: U_input(:)     ! U from input file (ntypat)
101     real(dp), allocatable :: J_input(:)     ! J from input file (ntypat)
102     real(dp), allocatable :: lambda(:)      ! ilamda
103     real(dp), allocatable :: docc(:,:,:)    ! n_In_j, natom, ilamda
104     real(dp), allocatable :: e_dc(:,:)      ! natom, ilamda
105     real(dp), allocatable :: uij(:,:)       ! uij to compute <uij n_i n_j>, ntypat
106   end type entropyDMFT_t