TABLE OF CONTENTS


ABINIT/m_Ctqmcoffdiag [ Modules ]

[ Top ] [ Modules ]

NAME

  m_Ctqmcoffdiag

FUNCTION

  Manage and drive all the CTQMC
  Should not be used if you don't know what you do
  Please use CtqmcoffdiagInterface

COPYRIGHT

  Copyright (C) 2013-2024 ABINIT group (J. Bieder, B. Amadon, J. Denier)
  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

24 #include "defs.h"
25 #define CTQMC_SLICE1 100
26 ! Coupe Sweeps en 100
27 #define CTQMC_SLICE2 100
28 ! Coupe modNoise1 en 2000
29 #define CTQMC_SEGME  1
30 #define CTQMC_ANTIS -2
31 #define CTQMC_ADDED  3  
32 #define CTQMC_REMOV  4
33 #define CTQMC_DETSI  5
34 MODULE m_Ctqmcoffdiag
35 
36  USE m_Global
37  USE m_GreenHyboffdiag
38  USE m_BathOperatoroffdiag
39  USE m_ImpurityOperator
40  USE m_Stat
41  USE m_FFTHyb
42  USE m_OurRng
43  use defs_basis
44 #ifdef HAVE_MPI2
45  USE mpi
46 #endif
47  IMPLICIT NONE
48 
49  public :: Ctqmcoffdiag_init
50  public :: Ctqmcoffdiag_setParameters
51  public :: Ctqmcoffdiag_setSweeps
52  public :: Ctqmcoffdiag_setSeed
53  public :: Ctqmcoffdiag_allocateAll
54  public :: Ctqmcoffdiag_allocateOpt
55  public :: Ctqmcoffdiag_setG0wTab
56  public :: Ctqmcoffdiag_setU
57  public :: Ctqmcoffdiag_clear
58  public :: Ctqmcoffdiag_reset
59  public :: Ctqmcoffdiag_setMu
60  public :: Ctqmcoffdiag_computeF
61  public :: Ctqmcoffdiag_run
62  public :: Ctqmcoffdiag_tryAddRemove
63  public :: Ctqmcoffdiag_trySwap
64  public :: Ctqmcoffdiag_measN
65  public :: Ctqmcoffdiag_measCorrelation
66  public :: Ctqmcoffdiag_measPerturbation
67  public :: Ctqmcoffdiag_getResult
68  public :: Ctqmcoffdiag_symmetrizeGreen
69  public :: Ctqmcoffdiag_getGreen
70  public :: Ctqmcoffdiag_getD
71  public :: Ctqmcoffdiag_getE
72  public :: Ctqmcoffdiag_printAll
73  public :: Ctqmcoffdiag_printQMC
74  public :: Ctqmcoffdiag_printGreen
75  public :: Ctqmcoffdiag_printD
76  public :: Ctqmcoffdiag_printE
77  public :: Ctqmcoffdiag_printPerturbation
78  public :: Ctqmcoffdiag_printCorrelation
79  public :: Ctqmcoffdiag_printSpectra
80  public :: Ctqmcoffdiag_destroy
81  public :: Ctqmcoffdiag_setMagmom

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_allocateAll [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_allocateAll

FUNCTION

  Allocate all non option variables

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

666 SUBROUTINE Ctqmcoffdiag_allocateAll(op)
667 
668 !Arguments ------------------------------------
669   TYPE(Ctqmcoffdiag), INTENT(INOUT) :: op
670 !Local variables ------------------------------
671   INTEGER                  :: flavors
672 
673   IF ( .NOT. op%para ) &
674     CALL ERROR("Ctqmcoffdiag_allocateAll : Ctqmcoffdiag_setParameters never called  ")
675 
676   flavors = op%flavors
677 
678 
679 !  number of electrons
680   FREEIF(op%measN)
681   MALLOC(op%measN,(1:4,1:flavors))
682   op%measN = 0.d0
683 
684 !  double occupancies 
685   FREEIF(op%measDE)
686   MALLOC(op%measDE,(1:flavors,1:flavors) )
687   op%measDE = 0.d0
688 
689   FREEIF(op%mu)
690 #ifdef FC_LLVM
691   ! LLVM 16 doesn't recognize this macro here
692   MALLOC(op%mu, (1:flavors) )
693 #else
694   MALLOC(op%mu, (1:flavors))
695 #endif
696   op%mu = 0.d0
697   FREEIF(op%hybri_limit)
698 #ifdef FC_LLVM
699   ! LLVM 16 doesn't recognize this macro here
700   MALLOC(op%hybri_limit, (flavors,flavors) )
701 #else
702   MALLOC(op%hybri_limit, (flavors,flavors))
703 #endif
704   op%hybri_limit = czero
705 END SUBROUTINE Ctqmcoffdiag_allocateAll

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_allocateOpt [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_allocateOpt

FUNCTION

  allocate all option variables 

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

733 SUBROUTINE Ctqmcoffdiag_allocateOpt(op)
734 
735 !Arguments ------------------------------------
736   TYPE(Ctqmcoffdiag), INTENT(INOUT) :: op
737 !Local variables ------------------------------
738   INTEGER :: i
739   INTEGER :: j
740   INTEGER :: k
741 
742   IF ( .NOT. op%para ) &
743     CALL ERROR("Ctqmcoffdiag_allocateOpt : Ctqmcoffdiag_setParameters never called  ")
744 
745   IF ( op%opt_analysis .EQ. 1 ) THEN
746     FREEIF(op%measCorrelation)
747     MALLOC(op%measCorrelation,(1:op%samples+1,1:3,1:op%flavors))
748     op%measCorrelation = 0.d0
749   END IF
750 
751   IF ( op%opt_order .GT. 0 ) THEN
752     FREEIF(op%measPerturbation)
753     MALLOC(op%measPerturbation,(1:op%opt_order,1:op%flavors))
754     op%measPerturbation = 0.d0
755     FREEIF(op%meas_fullemptylines)
756     MALLOC(op%meas_fullemptylines,(2,1:op%flavors))
757     op%meas_fullemptylines = 0.d0
758   END IF
759 
760   IF ( op%opt_histo .GT. 0 ) THEN
761     FREEIF(op%occup_histo_time)
762     MALLOC(op%occup_histo_time,(1:op%flavors+1))
763     op%occup_histo_time= 0.d0
764     FREEIF(op%occupconfig)
765     MALLOC(op%occupconfig,(1:2**op%flavors))
766     op%occupconfig= 0.d0
767     FREEIF(op%suscep)
768     MALLOC(op%suscep,(1:3,1:op%samples))
769     op%suscep= 0.d0
770     FREEIF(op%chi)
771     MALLOC(op%chi,(1:3,1:op%samples))
772     op%chi= 0.d0
773     FREEIF(op%chicharge)
774     MALLOC(op%chicharge,(1:3,1:op%samples))
775     op%chicharge= 0.d0
776     FREEIF(op%ntot)
777     MALLOC(op%ntot,(1:3))
778     op%ntot= 0.d0
779   END IF
780 
781   IF ( op%opt_noise .EQ. 1 ) THEN
782     IF ( ALLOCATED(op%measNoiseG) ) THEN
783       DO i=1,2
784         DO j = 1, op%flavors
785           DO k= 1, op%samples+1
786             CALL Vector_destroy(op%measNoiseG(k,j,i))
787           END DO
788         END DO
789       END DO
790       DT_FREE(op%measNoiseG)
791     END IF
792     DT_MALLOC(op%measNoiseG,(1:op%samples+1,1:op%flavors,1:2))
793     !DO i=1,2
794       DO j = 1, op%flavors
795         DO k= 1, op%samples+1
796           CALL Vector_init(op%measNoiseG(k,j,1),CTQMC_SLICE1)
797         END DO
798       END DO
799       DO j = 1, op%flavors
800         DO k= 1, op%samples+1
801           CALL Vector_init(op%measNoiseG(k,j,2),CTQMC_SLICE1*CTQMC_SLICE2+1) ! +1 pour etre remplacer ceil
802         END DO
803       END DO
804     !END DO
805     FREEIF(op%abNoiseG)
806     MALLOC(op%aBNoiseG,(1:2,1:op%samples+1,op%flavors))
807     op%abNoiseG = 0.d0
808   END IF
809 
810   IF (op%opt_spectra .GE. 1 ) THEN
811     FREEIF(op%density)
812     !MALLOC(op%density,(1:op%thermalization,1:op%flavors))
813     i = CEILING(DBLE(op%thermalization+op%sweeps)/DBLE(op%measurements*op%opt_spectra))
814     MALLOC(op%density,(1:op%flavors+1,1:i))
815     op%density = 0.d0
816   END IF
817 !#endif
818 END SUBROUTINE Ctqmcoffdiag_allocateOpt

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_clear [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_clear

FUNCTION

  clear a ctqmc run

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

1029 SUBROUTINE Ctqmcoffdiag_clear(op)
1030 
1031 !Arguments ------------------------------------
1032   TYPE(Ctqmcoffdiag), INTENT(INOUT) :: op
1033 !Local variables ------------------------------
1034   INTEGER :: i
1035   INTEGER :: j
1036   INTEGER :: k
1037 
1038   op%measN(1,:) = 0.d0
1039   op%measN(2,:) = 0.d0
1040   !Do not set measN(3,:) to 0 to avoid erasing N between therm and ctqmc
1041   op%measN(4,:) = 0.d0
1042   op%measDE = 0.d0
1043 !  op%seg_added    = 0.d0
1044 !  op%anti_added   = 0.d0
1045 !  op%seg_removed  = 0.d0
1046 !  op%anti_removed = 0.d0
1047 !  op%seg_sign     = 0.d0
1048 !  op%anti_sign    = 0.d0
1049   op%stats(:)     = 0.d0
1050 !  op%signvaluecurrent    = 0.d0
1051 !  op%signvaluemeas    = 0.d0
1052   op%swap         = 0.d0
1053   op%runTime      = 0.d0
1054   op%modGlobalMove(2) = 0 
1055   CALL Vector_clear(op%measNoise(1))
1056   CALL Vector_clear(op%measNoise(2))
1057 !#ifdef CTCtqmcoffdiag_CHECK
1058   op%errorImpurity = 0.d0
1059   op%errorBath     = 0.d0
1060 !#endif
1061   CALL GreenHyboffdiag_clear(op%Greens)
1062 !#ifdef CTCtqmcoffdiag_ANALYSIS
1063   IF ( op%opt_analysis .EQ. 1 .AND. ALLOCATED(op%measCorrelation) ) &    
1064     op%measCorrelation = 0.d0 
1065   IF ( op%opt_order .GT. 0 .AND. ALLOCATED(op%measPerturbation) ) &
1066     op%measPerturbation = 0.d0
1067   IF ( op%opt_order .GT. 0 .AND. ALLOCATED(op%meas_fullemptylines) ) &
1068     op%meas_fullemptylines = 0.d0
1069   IF ( op%opt_noise .EQ. 1 .AND. ALLOCATED(op%measNoiseG) ) THEN
1070     DO i=1,2
1071       DO j = 1, op%flavors
1072         DO k= 1, op%samples+1
1073           CALL Vector_clear(op%measNoiseG(k,j,i))
1074         END DO
1075       END DO
1076     END DO
1077     !DO j = 1, op%flavors
1078     !  CALL GreenHyboffdiag_clear(op%Greens(j))
1079     !END DO
1080   END IF
1081 !#endif
1082 END SUBROUTINE Ctqmcoffdiag_clear

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_computeF [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_computeF

FUNCTION

  Compute the hybridization function

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  Gomega=G0 to compute F
  opt_fk=What is Gomega

OUTPUT

  F=hybridization function

SIDE EFFECTS

NOTES

SOURCE

1251 SUBROUTINE Ctqmcoffdiag_computeF(op, Gomega, F, opt_fk)
1252 
1253  use m_hide_lapack,  only : xginv
1254 !Arguments ------------------------------------
1255   TYPE(Ctqmcoffdiag)                       , INTENT(INOUT) :: op
1256   COMPLEX(KIND=8), DIMENSION(:,:,:), INTENT(IN   ) :: Gomega
1257   !INTEGER                         , INTENT(IN   ) :: Wmax
1258   DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(INOUT) :: F
1259   INTEGER                         , INTENT(IN   ) :: opt_fk
1260 !Local variables ------------------------------
1261   INTEGER                                         :: flavors
1262   INTEGER                                         :: samples
1263   INTEGER                                         :: iflavor,ifl
1264   INTEGER                                         :: iflavor2
1265   INTEGER                                         :: iomega
1266   INTEGER                                         :: itau
1267   DOUBLE PRECISION                                :: pi_invBeta
1268   DOUBLE PRECISION                                :: K
1269   !DOUBLE PRECISION                                :: re
1270   !DOUBLE PRECISION                                :: im
1271   !DOUBLE PRECISION                                :: det
1272   COMPLEX(KIND=8), DIMENSION(:,:,:), ALLOCATABLE   :: F_omega
1273   COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE   :: F_omega_inv
1274   COMPLEX(KIND=8), DIMENSION(:,:,:), ALLOCATABLE   :: Gomega_tmp
1275   TYPE(GreenHyboffdiag)                                     :: F_tmp
1276   !character(len=4) :: tag_proc
1277   !character(len=30) :: tmpfil
1278   !INTEGER :: unitnb
1279 
1280   ABI_UNUSED((/opt_fk/))
1281 
1282   flavors    = op%flavors
1283 
1284   samples    = op%samples
1285   pi_invBeta = ACOS(-1.d0) / op%beta
1286   op%Wmax=SIZE(Gomega,1)
1287 !sui!write(std_out,*) "op%Wmax",op%Wmax
1288   !=================================
1289   ! --- Initialize F_tmp 
1290   !=================================
1291   IF ( op%have_MPI .EQV. .TRUE. ) THEN
1292     CALL GreenHyboffdiag_init(F_tmp,samples,op%beta,flavors,MY_COMM=op%MY_COMM)
1293   ELSE
1294     CALL GreenHyboffdiag_init(F_tmp,samples,op%beta,flavors)
1295   END IF
1296 !  K = op%mu
1297 
1298   !=================================
1299   ! --- Allocate F_omega
1300   !=================================
1301   MALLOC(F_omega,(1:op%Wmax,1:flavors,1:flavors))
1302   MALLOC(F_omega_inv,(1:flavors,1:flavors))
1303   MALLOC(Gomega_tmp,(1:op%Wmax,1:flavors,1:flavors))
1304   !op%hybri_limit(2,2)=op%hybri_limit(1,1)
1305   !op%mu(1)=op%mu(1)/10
1306   !op%mu(2)=op%mu(1)
1307   DO iomega=1,op%Wmax
1308     do iflavor=1,flavors
1309       do iflavor2=1,flavors
1310        ! Gomega_tmp(iomega,iflavor,iflavor2)=op%hybri_limit(iflavor,iflavor2)/(cmplx(0.d0,(2.d0*DBLE(iomega)-1.d0) * pi_invBeta))/3.d0
1311       enddo
1312     enddo
1313   END DO
1314   Gomega_tmp=Gomega
1315 
1316   !IF ( op%rank .EQ. 0 ) &
1317     !OPEN(UNIT=9876,FILE="K.dat",POSITION="APPEND")
1318   
1319   !=============================================================================================
1320   ! --- Compute Bath Green's function from Hybridization function in imaginary time
1321   !=============================================================================================
1322   !IF ( opt_fk .EQ. 0 ) THEN
1323    IF ( op%rank .EQ. 0 ) THEN
1324    !  DO iflavor = 1, flavors
1325    !    DO iflavor2 = 1, flavors
1326    !        write(330,*) "#",iflavor,iflavor2
1327    !        write(331,*) "#",iflavor,iflavor2
1328    !      do  iomega=1,op%Wmax
1329    !        write(330,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(Gomega_tmp(iomega,iflavor,iflavor2))
1330    !        write(331,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(Gomega_tmp(iomega,iflavor,iflavor2))
1331    !      enddo
1332    !        write(330,*) 
1333    !        write(331,*) 
1334    !    END DO
1335    !  END DO
1336    ENDIF
1337      DO iomega=1,op%Wmax
1338      !  be careful...here 
1339      ! Gomega in input is Fomega and 
1340      ! F_omega is   Gomega.
1341      ! COMPUTE G0 FROM F
1342        do iflavor=1,flavors
1343          do iflavor2=1,flavors
1344            if (iflavor==iflavor2) then
1345              F_omega_inv(iflavor,iflavor2)= (cmplx(0.d0,(2.d0*DBLE(iomega)-1.d0) * pi_invBeta,kind=8) &
1346 &             + op%mu(iflavor)- Gomega_tmp(iomega,iflavor,iflavor2))
1347            else
1348              F_omega_inv(iflavor,iflavor2)= (- Gomega_tmp(iomega,iflavor,iflavor2))
1349            endif
1350          enddo
1351        enddo
1352   !   END DO
1353   ! IF ( op%rank .EQ. 0 ) THEN
1354   !   DO iflavor = 1, flavors
1355   !     DO iflavor2 = 1, flavors
1356   !         write(334,*) "#",iflavor,iflavor2
1357   !         write(335,*) "#",iflavor,iflavor2
1358   !       do  iomega=1,op%Wmax
1359   !         write(334,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_omega(iomega,iflavor,iflavor2))
1360   !         write(335,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(F_omega(iomega,iflavor,iflavor2))
1361   !       enddo
1362   !         write(334,*) 
1363   !         write(335,*) 
1364   !     END DO
1365   !   END DO
1366   ! ENDIF
1367 
1368   !   DO iomega=1,op%Wmax
1369        call xginv(F_omega_inv,flavors)
1370        do iflavor=1,flavors
1371          do iflavor2=1,flavors
1372            F_omega(iomega,iflavor,iflavor2) = F_omega_inv(iflavor,iflavor2)
1373          enddo
1374        enddo
1375      END DO
1376 
1377    !IF ( op%rank .EQ. 0 ) THEN
1378    !  DO iflavor = 1, flavors
1379    !    DO iflavor2 = 1, flavors
1380    !        write(332,*) "#",iflavor,iflavor2
1381    !        write(333,*) "#",iflavor,iflavor2
1382    !      do  iomega=1,op%Wmax
1383    !        write(332,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_omega(iomega,iflavor,iflavor2))
1384    !        write(333,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(F_omega(iomega,iflavor,iflavor2))
1385    !      enddo
1386    !        write(332,*) 
1387    !        write(333,*) 
1388    !    END DO
1389    !  END DO
1390    !ENDIF
1391      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392      !  for test: Fourier of G0(iwn)
1393      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1394    !sui!write(std_out,*) "opt_fk=0"
1395      CALL GreenHyboffdiag_setOperW(F_tmp,F_omega)
1396   ! IF ( op%rank .EQ. 0 ) THEN
1397   !    DO iflavor = 1, flavors
1398   !      DO iflavor2 = 1, flavors
1399   !          write(336,*) "#",iflavor,iflavor2
1400   !          write(337,*) "#",iflavor,iflavor2
1401   !        do  iomega=1,op%Wmax
1402   !          write(336,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_tmp%oper_w(iomega,iflavor,iflavor2))
1403   !          write(337,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(F_tmp%oper_w(iomega,iflavor,iflavor2))
1404   !        enddo
1405   !          write(336,*) 
1406   !          write(337,*) 
1407   !      END DO
1408   !    END DO
1409   !  ENDIF
1410      CALL GreenHyboffdiag_backFourier(F_tmp,func="green")
1411      ! --- Put the result in F
1412      DO iflavor = 1, flavors
1413        DO iflavor2 = 1, flavors
1414          DO itau=1,samples+1
1415          F(itau,iflavor,iflavor2) = F_tmp%oper(itau,iflavor,iflavor2)
1416          END DO
1417        END DO
1418      END DO
1419      !IF ( op%rank .EQ. 0 ) THEN
1420      !  DO iflavor = 1, flavors
1421      !    DO iflavor2 = 1, flavors
1422      !        write(346,*) "#",iflavor,iflavor2
1423      !      do  itau=1,op%samples+1
1424      !        write(346,*) (itau-1)*op%beta/(op%samples),real(F(itau,iflavor,iflavor2))
1425      !      enddo
1426      !        write(346,*) 
1427      !    END DO
1428      !  END DO
1429      !ENDIF
1430      DO iflavor = 1, flavors
1431        DO iflavor2 = 1, flavors
1432          DO itau=1,samples+1
1433 !         This symetrization is general and valid even with SOC
1434 !         Without SOC, it leads to zero.
1435          F(itau,iflavor,iflavor2) = (F_tmp%oper(itau,iflavor,iflavor2)+F_tmp%oper(itau,iflavor2,iflavor))/2.d0
1436          END DO
1437        END DO
1438      END DO
1439      open (unit=4367,file='G0tau_fromF',status='unknown',form='formatted')
1440      rewind(4367)
1441      IF ( op%rank .EQ. 0 ) THEN
1442        DO iflavor = 1, flavors
1443          DO iflavor2 = 1, flavors
1444              write(4367,*) "#",iflavor,iflavor2
1445            do  itau=1,op%samples+1
1446              write(4367,*) (itau-1)*op%beta/(op%samples),F(itau,iflavor,iflavor2)
1447            enddo
1448              write(4367,*) 
1449          END DO
1450        !sui!write(std_out,'(5x,14(2f9.5,2x))') (F(op%samples+1,iflavor,iflavor2),iflavor2=1,flavors)
1451        END DO
1452      ENDIF
1453      !call flush(437)
1454      close(4367)
1455      !call flush(6)
1456      
1457      call xmpi_barrier(op%MY_COMM)
1458      !CALL ERROR("END OF CALCULATION")
1459      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1460      !  END OF TEST
1461      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1462 
1463      !DO iomega=1,op%Wmax
1464      !  call xginv(F_omega(iomega,:,:),flavors)
1465      !END DO
1466     !F_omega = CMPLX(-1.d0,0,8)/Gomega_tmp
1467   !ELSE
1468   !=============================================================================================
1469   ! --- Restore Hybridization in F_omega
1470   !=============================================================================================
1471 
1472 !   Restore Hybridization in F_omega for the following operations
1473   F_omega = Gomega_tmp
1474   !END IF
1475 
1476   !==================================================================
1477   ! --- Full double loop on flavors to compute F (remove levels)
1478   !==================================================================
1479   DO iflavor = 1, flavors
1480     DO iflavor2 = 1, flavors
1481 
1482   ! --- Compute or use the levels for the diagonal hybridization (else K=0)
1483       IF(iflavor==iflavor2) THEN
1484         IF ( op%opt_levels .EQ. 1 ) THEN
1485           K = op%mu(iflavor)
1486         ELSE
1487           K = -REAL(F_omega(op%Wmax, iflavor,iflavor))
1488 !        op%mu = K
1489           op%mu(iflavor) = K 
1490         END IF
1491       ELSE
1492         K=0.d0
1493       ENDIF
1494       !IF ( op%rank .EQ. 0 ) &
1495       !WRITE(9876,'(I4,2E22.14)') iflavor, K, REAL(-F_omega(op%Wmax, iflavor))
1496      ! IF(op%rank .EQ.0) &
1497      ! WRITE(op%ostream,*) "CTQMC K, op%mu = ",K,op%mu(iflavor)
1498       !WRITE(op%ostream,*) "CTQMC beta     = ",op%beta
1499 
1500   ! --- Compute F (by removing the levels) if opt_fk==0
1501     !  IF ( opt_fk .EQ. 0 ) THEN
1502     !   ! DO iomega = 1, op%Wmax
1503     !   !   re = REAL(F_omega(iomega,iflavor,iflavor2))
1504     !   !   im = AIMAG(F_omega(iomega,iflavor,iflavor2))
1505     !   !   if (iflavor==iflavor2) then
1506     !   !     F_omega(iomega,iflavor,iflavor) = CMPLX(re + K, im + (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, 8)
1507     !   !   else
1508     !   !     F_omega(iomega,iflavor,iflavor2) = CMPLX(re , im  , 8)
1509     !   !   endif
1510     !   !   !if(iflavor==1.and.op%rank==0) then
1511     !   !     !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
1512     !   !     !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega_tmp(iomega, iflavor)),imag(Gomega_tmp(iomega, iflavor))
1513     !   !   !end if 
1514     !   ! END DO
1515     !  ELSE
1516     !    DO iomega = 1, op%Wmax
1517     !      !F_omega(iomega,iflavor,iflavor2) = F_omega(iomega,iflavor,iflavor2) + CMPLX(K, 0.d0, 8)
1518 
1519 
1520     !      !if(iflavor==1.and.op%rank==0) then
1521     !        !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
1522     !        !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega_tmp(iomega, iflavor)),imag(Gomega_tmp(iomega, iflavor))
1523     !      !end if 
1524     !    END DO
1525     !  END IF
1526   ! --- compute residual K (?)
1527       K = REAL(CMPLX(0,(2.d0*DBLE(op%Wmax)-1.d0)*pi_invBeta,8)*F_omega(op%Wmax,iflavor,iflavor2))
1528       CALL GreenHyboffdiag_setMuD1(op%Greens,iflavor,iflavor2,op%mu(iflavor),K)
1529     END DO
1530   END DO
1531 
1532   do  iomega=1,op%Wmax
1533    ! write(336,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_omega(iomega,1,1)),imag(F_omega(iomega,1,1))
1534   enddo
1535 
1536   ! --- Creates F_tmp%oper_w
1537   CALL GreenHyboffdiag_setOperW(F_tmp,F_omega)
1538  ! do  iflavor=1, flavors ; do  iflavor2=1, flavors ; write(337,*) "#",iflavor,iflavor2 ; do  iomega=1,op%Wmax
1539  !   write(337,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_tmp%oper_w(iomega,iflavor,iflavor2)),&
1540  !&   imag(F_tmp%oper_w(iomega,iflavor,iflavor2))
1541  ! enddo ; write(337,*) ; enddo ; enddo
1542 !  IF ( op%rank .EQ. 0 ) THEN
1543 !    DO iflavor = 1, flavors
1544 !      DO iflavor2 = 1, flavors
1545 !        write(336,*) "#",iflavor,iflavor2
1546 !        write(337,*) "#",iflavor,iflavor2
1547 !        do  iomega=1,op%Wmax
1548 !          write(336,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_tmp%oper_w(iomega,iflavor,iflavor2))
1549 !          write(337,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(F_tmp%oper_w(iomega,iflavor,iflavor2))
1550 !        enddo
1551 !        write(336,*) 
1552 !        write(337,*) 
1553 !        write(136,*) "#",iflavor,iflavor2
1554 !        write(137,*) "#",iflavor,iflavor2
1555 !        do  iomega=1,op%Wmax
1556 !        write(136,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_tmp%oper_w(iomega,iflavor,iflavor2)-op%hybri_limit(iflavor,iflavor2)/(cmplx(0.d0,(2.d0*DBLE(iomega)-1.d0) * pi_invBeta,kind=8)))
1557 !        write(137,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(F_tmp%oper_w(iomega,iflavor,iflavor2)-op%hybri_limit(iflavor,iflavor2)/(cmplx(0.d0,(2.d0*DBLE(iomega)-1.d0) * pi_invBeta,kind=8)))
1558 !        enddo
1559 !        write(136,*) 
1560 !        write(137,*) 
1561 !        write(836,*) "#",iflavor,iflavor2
1562 !        write(837,*) "#",iflavor,iflavor2
1563 !        do  iomega=1,op%Wmax
1564 !        write(836,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(op%hybri_limit(iflavor,iflavor2)/(cmplx(0.d0,(2.d0*DBLE(iomega)-1.d0) * pi_invBeta)))
1565 !        write(837,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,imag(op%hybri_limit(iflavor,iflavor2)/(cmplx(0.d0,(2.d0*DBLE(iomega)-1.d0) * pi_invBeta)))
1566 !        enddo
1567 !        write(836,*) 
1568 !        write(837,*) 
1569 !      END DO
1570 !    END DO
1571 !  ENDIF
1572   !CALL GreenHyboffdiag_backFourier(F_tmp,F_omega(:,iflavor))
1573    ! DO iflavor = 1, flavors
1574    !   DO iflavor2 = 1, flavors
1575    !   unitnb=80000+F_tmp%rank
1576    !   call int2char4(F_tmp%rank,tag_proc)
1577    !   tmpfil = 'oper_wavantFOURIER'//tag_proc
1578    !   open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1579    !   write(unitnb,*) "#",iflavor,iflavor2
1580    !   ! C_omega et oper_w differents Domega identique. Est ce du a des
1581    !   ! diago differentes   pour chaque procs dans qmc_prep_ctqmc
1582    !   do  iomega=1,F_tmp%Wmax
1583    !   write(unitnb,*)  (2.d0*DBLE(iomega)-1.d0) * pi_invBeta,real(F_tmp%oper_w(iomega,iflavor,iflavor2))
1584    !   enddo
1585    !   write(unitnb,*) 
1586    !   END DO
1587    ! END DO
1588 
1589   ! --- For all iflavor and iflavor2, do the Fourier transformation to
1590   ! --- have (F(\tau))
1591   !CALL GreenHyboffdiag_backFourier(F_tmp,hybri_limit=op%hybri_limit,opt_hybri_limit=op%opt_hybri_limit)
1592    write(std_out,*) "WARNING opt_hybri_limit==0"
1593   CALL GreenHyboffdiag_backFourier(F_tmp,hybri_limit=op%hybri_limit,opt_hybri_limit=0)
1594 !  CALL GreenHyboffdiag_backFourier(F_tmp,hybri_limit=op%hybri_limit,opt_hybri_limit=1)
1595  ! CALL GreenHyboffdiag_backFourier(F_tmp)
1596 
1597   ! --- Put the result in F
1598   DO iflavor = 1, flavors
1599     DO iflavor2 = 1, flavors
1600       DO itau=1,samples+1
1601       F(itau,iflavor,iflavor2) = -F_tmp%oper(samples+2-itau,iflavor,iflavor2)
1602       END DO
1603     END DO
1604   END DO
1605 ! IF ( op%rank .EQ. 0 ) THEN
1606 !   ifl=0
1607 !   DO iflavor = 1, flavors
1608 !     DO iflavor2 = 1, flavors
1609 !       ifl=ifl+1
1610 !       write(346,*) "#",iflavor,iflavor2,ifl
1611 !       do  itau=1,op%samples+1
1612 !         write(346,*) itau,real(F(itau,iflavor,iflavor2))
1613 !       enddo
1614 !       write(346,*) 
1615 !     END DO
1616 !   END DO
1617 ! ENDIF
1618 ! close(346)
1619   DO iflavor = 1, flavors
1620     DO iflavor2 = 1, flavors
1621       DO itau=1,samples+1
1622 !      This symetrization is general and valid even with SOC
1623 !      Without SOC, it leads to zero.
1624       F(itau,iflavor,iflavor2) = -(F_tmp%oper(samples+2-itau,iflavor,iflavor2)+F_tmp%oper(samples+2-itau,iflavor2,iflavor))/2.d0
1625       END DO
1626     END DO
1627   END DO
1628   !DO iflavor = 1, flavors
1629   !  DO iflavor2 = 1, flavors
1630   !    DO itau=1,samples+1
1631   !    F(itau,iflavor,iflavor2) = F(samples/2,iflavor,iflavor2)
1632   !    END DO
1633   !  END DO
1634   !END DO
1635 
1636  !  SOME TRY TO ADJUST F
1637   !DO iflavor = 1, flavors
1638   !  DO iflavor2 = 1, flavors
1639   !    do  itau=1,op%samples+1
1640   !    !if(iflavor/=iflavor2) F(itau,iflavor,iflavor2)=F((op%samples+1)/2,iflavor,iflavor2)
1641   !    !if(iflavor==iflavor2) F(itau,iflavor,iflavor2)=F((op%samples+1)/2,iflavor,iflavor2)
1642   !    enddo
1643   !  END DO
1644   !END DO
1645   !write(6,*) "QQQQ1",op%rank
1646 
1647   open (unit=436,file='Hybridization.dat',status='unknown',form='formatted')
1648   rewind(436)
1649   IF ( op%rank .EQ. 0 ) THEN
1650     ifl=0
1651     DO iflavor = 1, flavors
1652       DO iflavor2 = 1, flavors
1653         ifl=ifl+1
1654           write(436,*) "#",iflavor,iflavor2 !,ifl,op%hybri_limit(iflavor,iflavor2)
1655         do  itau=1,op%samples+1
1656           write(436,*) itau,F(itau,iflavor,iflavor2)
1657         enddo
1658           write(436,*) 
1659       END DO
1660     END DO
1661   ENDIF
1662   close(436)
1663   !   call xmpi_barrier(op%MY_COMM)
1664   !write(6,*) "QQQQ3"
1665   FREE(Gomega_tmp)
1666   FREE(F_omega)
1667   FREE(F_omega_inv)
1668   !write(6,*) "QQQQ4"
1669   CALL GreenHyboffdiag_destroy(F_tmp)
1670   !write(6,*) "QQQQ2"
1671 
1672 
1673 END SUBROUTINE Ctqmcoffdiag_computeF

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_destroy [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_destroy

FUNCTION

  destroy and deallocate all variables

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4776 SUBROUTINE Ctqmcoffdiag_destroy(op)
4777 
4778 !Arguments ------------------------------------
4779   TYPE(Ctqmcoffdiag), INTENT(INOUT) :: op
4780 !Local variables ------------------------------
4781   !INTEGER                  :: iflavor
4782   INTEGER                  :: flavors
4783   INTEGER                  :: i
4784   INTEGER                  :: j
4785   INTEGER                  :: k
4786 
4787   flavors = op%flavors
4788 
4789   CALL ImpurityOperator_destroy(op%Impurity)
4790   CALL BathOperatoroffdiag_destroy(op%Bath)
4791   CALL Vector_destroy(op%measNoise(1))
4792   CALL Vector_destroy(op%measNoise(2))
4793 
4794 !sui!write(6,*) "before greenhyb_destroy in ctmqc_destroy"
4795   CALL GreenHyboffdiag_destroy(op%Greens)
4796 !#ifdef CTCtqmcoffdiag_ANALYSIS
4797   FREEIF(op%measCorrelation)
4798   FREEIF(op%measPerturbation)
4799   FREEIF(op%meas_fullemptylines)
4800   FREEIF(op%measN)
4801   IF ( op%opt_histo .GT. 0 ) THEN
4802     FREEIF(op%occup_histo_time)
4803     FREEIF(op%suscep)
4804     FREEIF(op%occupconfig)
4805     FREEIF(op%chi)
4806     FREEIF(op%chicharge)
4807     FREEIF(op%ntot)
4808   END IF
4809   FREEIF(op%measDE)
4810   FREEIF(op%mu)
4811   FREEIF(op%hybri_limit)
4812   FREEIF(op%abNoiseG)
4813   IF ( ALLOCATED(op%measNoiseG) ) THEN
4814     DO i=1,2
4815       DO j = 1, op%flavors
4816         DO k= 1, op%samples+1
4817           CALL Vector_destroy(op%measNoiseG(k,j,i))
4818         END DO
4819       END DO
4820     END DO
4821     DT_FREE(op%measNoiseG)
4822   END IF
4823   FREEIF(op%density)
4824 !#endif
4825   op%ostream        = 0
4826   op%istream        = 0
4827  
4828   op%sweeps         = 0
4829   op%thermalization = 0
4830   op%flavors        = 0
4831   op%samples        = 0
4832   op%beta           = 0.d0
4833 !  op%seg_added      = 0.d0
4834 !  op%anti_added     = 0.d0
4835 !  op%seg_removed    = 0.d0
4836 !  op%anti_removed   = 0.d0
4837 !  op%seg_sign       = 0.d0
4838 !  op%anti_sign      = 0.d0
4839   op%stats          = 0.d0
4840   op%swap           = 0.d0
4841 
4842 
4843   op%set  = .FALSE.
4844   op%done = .FALSE.
4845   op%init = .FALSE.
4846 END SUBROUTINE Ctqmcoffdiag_destroy

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_getD [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_getD

FUNCTION

  get double occupation

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

  D=full double occupation

SIDE EFFECTS

NOTES

SOURCE

4040 SUBROUTINE Ctqmcoffdiag_getD(op, D)
4041 
4042 !Arguments ------------------------------------
4043   TYPE(Ctqmcoffdiag)       , INTENT(IN ) :: op
4044   DOUBLE PRECISION, INTENT(OUT) :: D
4045 !Local variables ------------------------------
4046   INTEGER                       :: iflavor1
4047   INTEGER                       :: iflavor2
4048 
4049   D = 0.d0
4050 
4051   DO iflavor1 = 1, op%flavors
4052     DO iflavor2 = iflavor1+1, op%flavors
4053       D = D + op%measDE(iflavor2,iflavor1)
4054     END DO
4055   END DO
4056   !IF ( op%rank .EQ. 0 ) THEN
4057   !  DO iflavor1 = 1, op%flavors
4058   !    DO iflavor2 = iflavor1+1, op%flavors
4059   !     write(4533,*) op%measDE(iflavor2,iflavor1)k
4060   !     write(4534,*) op%Impurity%mat_U(iflavor2,iflavor1)k
4061   !    END DO
4062   !  END DO
4063 
4064   !ENDIF
4065 
4066 END SUBROUTINE Ctqmcoffdiag_getD

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_getE [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_getE

FUNCTION

  get interaction energy and noise on it

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

  E=interaction energy
  noise=noise on this value

SIDE EFFECTS

NOTES

SOURCE

4095 SUBROUTINE Ctqmcoffdiag_getE(op,E,noise)
4096 
4097 !Arguments ------------------------------------
4098   TYPE(Ctqmcoffdiag)       , INTENT(IN ) :: op
4099   DOUBLE PRECISION, INTENT(OUT) :: E
4100   DOUBLE PRECISION, INTENT(OUT) :: Noise
4101 
4102   E = op%measDE(1,1)  
4103   Noise = op%a_Noise*(DBLE(op%sweeps)*DBLE(op%size))**op%b_Noise
4104 END SUBROUTINE Ctqmcoffdiag_getE

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_getGreen [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_getGreen

FUNCTION

  Get the full green functions in time and/or frequency

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

  Gtau=green function in time
  Gw=green function in frequency

SIDE EFFECTS

NOTES

SOURCE

3769 SUBROUTINE Ctqmcoffdiag_getGreen(op, Gtau, Gw)
3770 
3771 !Arguments ------------------------------------
3772  USE m_GreenHyboffdiag
3773   TYPE(Ctqmcoffdiag)          , INTENT(INOUT)    :: op
3774   DOUBLE PRECISION, DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: Gtau
3775   COMPLEX(KIND=8), DIMENSION(:,:,:), OPTIONAL, INTENT(INOUT) :: Gw
3776 !Local variables ------------------------------
3777   !INTEGER                            :: itime
3778   INTEGER                            :: iflavor1
3779   INTEGER                            :: iflavor1b !,iflavor,iflavorbis
3780   INTEGER                            :: iflavor2
3781   INTEGER                            :: iflavor3
3782   INTEGER                            :: flavors,tail
3783   INTEGER                            :: ifreq,itime
3784   DOUBLE PRECISION :: u1 
3785   DOUBLE PRECISION :: u2
3786   DOUBLE PRECISION :: u3
3787   DOUBLE PRECISION :: Un
3788   DOUBLE PRECISION :: UUnn,iw !omega,
3789   CHARACTER(LEN=4)                   :: cflavors
3790   CHARACTER(LEN=50)                  :: string
3791   TYPE(GreenHyboffdiag)                     :: F_tmp
3792 
3793   flavors = op%flavors
3794   DO iflavor1 = 1, flavors
3795     u1 = 0.d0
3796     u2 = 0.d0
3797     u3 = 0.d0
3798     DO iflavor2 = 1, flavors
3799       IF ( iflavor2 .EQ. iflavor1 ) CYCLE
3800       Un = op%Impurity%mat_U(iflavor2,iflavor1) * op%measN(1,iflavor2)
3801 !      Un = op%Impurity%mat_U(iflavor2,iflavor1) * (op%Greens%oper(1,iflavor2,iflavor2) + 1.d0)
3802       !write(6,*) "forsetmoments",iflavor1,iflavor2,(op%Greens%oper(1,iflavor2,iflavor2) + 1.d0), Un
3803       u1 = u1 + Un 
3804       u2 = u2 + Un*op%Impurity%mat_U(iflavor2,iflavor1) 
3805       DO iflavor3 = 1, flavors
3806         IF ( iflavor3 .EQ. iflavor2 .OR. iflavor3 .EQ. iflavor1 ) CYCLE
3807         UUnn = (op%Impurity%mat_U(iflavor2,iflavor1)*op%Impurity%mat_U(iflavor3,iflavor1)) * &
3808 &                                                    op%measDE(iflavor2,iflavor3) 
3809         u2 = u2 + UUnn 
3810       END DO
3811     END DO  
3812      ! write(6,*) "u1,u2",u1,u2
3813 
3814     DO iflavor1b = 1, flavors
3815       u3 =-(op%Impurity%mat_U(iflavor1,iflavor1b))*op%Greens%oper(1,iflavor1,iflavor1b)
3816       ! u3=U_{1,1b}*G_{1,1b}
3817       CALL GreenHyboffdiag_setMoments(op%Greens,iflavor1,iflavor1b,u1,u2,u3)
3818     END DO ! iflavor1b
3819 
3820   END DO ! iflavor1
3821 
3822   IF ( PRESENT( Gtau ) ) THEN
3823     DO iflavor1 = 1, flavors
3824       DO iflavor2 = 1, flavors
3825         Gtau(1:op%samples,iflavor1,iflavor2) = op%Greens%oper(1:op%samples,iflavor1,iflavor2)
3826       END DO  
3827     END DO ! iflavor1
3828   END IF
3829 ! !--------- Write Occupation matrix before Gtau
3830 !  write(ostream,'(17x,a)') "Occupation matrix"
3831 !  write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
3832 !  write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
3833 !  do iflavor=1, op%flavors
3834 !    write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(-op%Greens%oper(op%samples,iflavor,iflavorbis),iflavorbis=1,op%flavors)
3835 !  enddo
3836 !  write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
3837 ! !------------------------------------------------------------------------------------------
3838 ! !--------- Write Occupation matrix Gtau
3839 !  write(ostream,'(17x,a)') "Occupation matrix"
3840 !  write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
3841 !  write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
3842 !  do iflavor=1, op%flavors
3843 !    write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(Gtau(op%samples,iflavor,iflavorbis),iflavorbis=1,op%flavors)
3844 !  enddo
3845 !  write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
3846 ! !------------------------------------------------------------------------------------------
3847 
3848 !================================================
3849   if(3==4) then
3850 !================================================
3851     DO iflavor1 = 1, flavors
3852       DO iflavor1b = 1, flavors
3853        !call nfourier3(op%Greens%oper(1:op%samples,iflavor1,iflavor1b),Gw(1:op%samples,iflavor1,iflavor1b),iflavor1==iflavor1b,op%Greens%samples,op%Greens%samples-1,op%Greens%beta,1.d0,op%Greens%Mk(iflavor1,iflavor1b,1),op%Greens%Mk(iflavor1,iflavor1b,2),op%Greens%Mk(iflavor1,iflavor1b,3))
3854       END DO  
3855     END DO ! iflavor1
3856 !    ============== write Gomega_nd.dat
3857     if(op%rank==0) then
3858     OPEN(UNIT=44, FILE="Gomega_nd_nfourier2.dat")
3859     WRITE(cflavors,'(I4)') 2*(flavors*flavors+1)
3860     string = '(1x,'//TRIM(ADJUSTL(cflavors))//'E15.5)'
3861     !write(6,*) " op%Greens%Wmax", op%Greens%Wmax
3862     do  iflavor1=1, flavors 
3863       do  iflavor1b=1, flavors 
3864     write(44,*) "#op%Greens%Mk(iflavor1,iflavor2,1",op%Greens%Mk(iflavor1,iflavor1b,:)
3865         DO ifreq = 1, op%samples
3866 !      !write(6,string) (DBLE(ifreq)*2-1)*3.1415/op%Greens%beta, &
3867 !      (/ ((real(Gw(ifreq,iflavor1,iflavor1b)),imag(Gw(ifreq,iflavor1,iflavor1b)), iflavor1=1, flavors),iflavor1b=1,flavors) /)
3868 !      WRITE(44,string) (DBLE(ifreq)*2.d0-1.d0)*3.1415926/op%Greens%beta, &
3869         iw=aimag(Gw(ifreq,op%flavors,op%flavors+1))
3870          WRITE(44,string) aimag(Gw(ifreq,op%flavors,op%flavors+1)),&
3871          real(Gw(ifreq,iflavor1,iflavor1b)),aimag(Gw(ifreq,iflavor1,iflavor1b)),&
3872             ( -op%Greens%Mk(iflavor1,iflavor1b,2) )/(iw*iw) , (op%Greens%Mk(iflavor1,iflavor1b,1))/iw!-op%Greens%Mk(iflavor1,iflavor1b,3)/(iw*iw))/iw 
3873       !   WRITE(102,*) aimag(Gw(ifreq,op%flavors,op%flavors+1)), (op%Greens%Mk(iflavor1,iflavor1b,1))/iw,op%Greens%Mk(iflavor1,iflavor1b,1),iw
3874         END DO
3875          WRITE(44,*) 
3876       END DO
3877     END DO
3878     close(44)
3879     endif
3880 !================================================
3881   endif
3882 !================================================
3883        !!write(6,*) "present gw", present(gw)
3884   IF ( PRESENT( Gw ) ) THEN
3885      !!write(6,*) "size gw",SIZE(Gw,DIM=2) ,flavors+1 
3886     IF ( SIZE(Gw,DIM=3) .EQ. flavors+1 ) THEN
3887      ! CALL GreenHyboffdiag_forFourier(op%Greens, Gomega=Gw, omega=Gw(:,op%flavors,op%flavors+1))
3888       CALL GreenHyboffdiag_forFourier(op%Greens, Gomega=Gw, omega=Gw(:,op%flavors,op%flavors+1))
3889       !write(6,*) "1"
3890       !IF ( op%rank .EQ. 0 ) write(20,*) Gw(:,iflavor1)
3891     ELSE IF ( SIZE(Gw,DIM=3) .EQ. flavors ) THEN  
3892       CALL GreenHyboffdiag_forFourier(op%Greens,Gomega=Gw)
3893       !write(6,*) "2"
3894     ELSE
3895       CALL WARNALL("Ctqmcoffdiag_getGreen : Gw is not valid                    ")
3896       CALL GreenHyboffdiag_forFourier(op%Greens,Wmax=op%Wmax)
3897       !write(6,*) "3"
3898     END IF
3899   ELSE
3900     CALL GreenHyboffdiag_forFourier(op%Greens,Wmax=op%Wmax)
3901   END IF
3902 !  ============== write Gomega_nd.dat
3903 !================================================
3904 !  if(3==4) then
3905 !================================================
3906   if(op%rank==0.and.3==4) then
3907   OPEN(UNIT=44, FILE="Gomega_nd.dat")
3908   WRITE(cflavors,'(I4)') 2*(flavors*flavors+1)
3909   string = '(1x,'//TRIM(ADJUSTL(cflavors))//'E15.5)'
3910   !write(6,*) " op%Greens%Wmax", op%Greens%Wmax
3911   do  iflavor1=1, flavors 
3912     do  iflavor1b=1, flavors 
3913   write(44,*) "#op%Greens%Mk(iflavor1,iflavor2,1",op%Greens%Mk(iflavor1,iflavor1b,:)
3914       DO ifreq = 1, SIZE(Gw,1)    
3915 !    !write(6,string) (DBLE(ifreq)*2-1)*3.1415/op%Greens%beta, &
3916 !    (/ ((real(Gw(ifreq,iflavor1,iflavor1b)),imag(Gw(ifreq,iflavor1,iflavor1b)), iflavor1=1, flavors),iflavor1b=1,flavors) /)
3917 !    WRITE(44,string) (DBLE(ifreq)*2.d0-1.d0)*3.1415926/op%Greens%beta, &
3918       iw=aimag(Gw(ifreq,op%flavors,op%flavors+1))
3919        WRITE(44,string) aimag(Gw(ifreq,op%flavors,op%flavors+1)),&
3920        real(Gw(ifreq,iflavor1,iflavor1b)),aimag(Gw(ifreq,iflavor1,iflavor1b)),&
3921           ( -op%Greens%Mk(iflavor1,iflavor1b,2) )/(iw*iw) , &
3922 &          (op%Greens%Mk(iflavor1,iflavor1b,1)-op%Greens%Mk(iflavor1,iflavor1b,3)/(iw*iw))/iw 
3923       END DO
3924        WRITE(44,*) 
3925     END DO
3926   END DO
3927   endif
3928 !================================================
3929 !  endif
3930 !================================================
3931 
3932 
3933 !  ==============================
3934   ! --- Initialize F_tmp 
3935   !write(6,*) "10"
3936 
3937   IF ( op%have_MPI .EQV. .TRUE. ) THEN
3938     !CALL GreenHyboffdiag_init(F_tmp,op%samples,op%beta,op%flavors,MY_COMM=op%MY_COMM)
3939     CALL GreenHyboffdiag_init(F_tmp,op%samples,op%beta,flavors)
3940     !write(6,*) "10a"
3941   ELSE
3942     CALL GreenHyboffdiag_init(F_tmp,op%samples,op%beta,flavors)
3943     !write(6,*) "10b"
3944   END IF
3945 
3946   !write(6,*) "11"
3947 !  CALL GreenHyboffdiag_setOperW(F_tmp,Gw)
3948 
3949   tail = op%samples 
3950   F_tmp%Wmax=op%samples ! backFourier only works for linear freq: calculation of A and etc..
3951   MALLOC(F_tmp%oper_w,(1:tail,op%flavors,op%flavors))
3952   F_tmp%oper_w(1:tail,1:F_tmp%nflavors,1:F_tmp%nflavors) = Gw(1:tail,1:F_tmp%nflavors,1:F_tmp%nflavors)
3953   !write(6,*) "example",F_tmp%oper_w(1,1,1)
3954   !write(6,*) "example",Gw(1,1,1)
3955   F_tmp%setW = .TRUE.
3956   !write(6,*) size(F_tmp%oper_w,1)
3957   !write(6,*) size(F_tmp%oper_w,2)
3958   !write(6,*) size(F_tmp%oper_w,3)
3959   !write(6,*) size(Gw,1)
3960   !write(6,*) size(Gw,2)
3961   !write(6,*) size(Gw,3)
3962 
3963   !write(6,*) "eee", (2.d0*DBLE(ifreq)-1.d0) * 3.1415/op%beta,real(F_tmp%oper_w(1,1,1)),imag(F_tmp%oper_w(1,1,1))
3964 !================================================
3965   if(3==4) then
3966 !================================================
3967     OPEN(UNIT=3337, FILE="Gomega_nd2.dat")
3968     do  iflavor1=1, flavors 
3969       do  iflavor1b=1, flavors 
3970         do  ifreq=1, tail
3971 !         write(3337,*) (2.d0*DBLE(ifreq)-1.d0) * 3.1415/op%beta,real(F_tmp%oper_w(ifreq,iflavor1,iflavor1b)),&
3972 !   &     imag(F_tmp%oper_w(ifreq,iflavor1,iflavor1b))
3973           write(3337,*) aimag(Gw(ifreq,op%flavors,op%flavors+1)), real(F_tmp%oper_w(ifreq,iflavor1,iflavor1b)),&
3974  &        aimag(F_tmp%oper_w(ifreq,iflavor1,iflavor1b))
3975 
3976       !    omega=(2.d0*DBLE(ifreq)-1.d0) * 3.1415/op%beta
3977  !        F_tmp%oper_w(ifreq,iflavor1,iflavor1b)=0.1**2/Gw(ifreq,op%flavors,op%flavors+1)
3978         enddo 
3979         write(3337,*)
3980       enddo 
3981     enddo
3982     close(3337)
3983 !================================================
3984   endif
3985 !================================================
3986 
3987   !write(6,*) "12",F_tmp%Wmax
3988 
3989 !  CALL GreenHyboffdiag_backFourier(F_tmp,func="green")
3990 
3991   !write(6,*) "13"
3992 
3993 !================================================
3994   if(3==4) then
3995 !================================================
3996     OPEN(UNIT=48, FILE="Gtau_nd_2.dat")
3997 !    --- Print full non diagonal Gtau in Gtau_nd.dat
3998     WRITE(cflavors,'(I4)') flavors*flavors+1
3999     string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
4000     DO itime = 1, op%samples+1
4001       WRITE(48,string) DBLE(itime-1)*op%beta/DBLE(op%samples), &
4002  &    ((F_tmp%oper(itime,iflavor1,iflavor1b), iflavor1=1, flavors),iflavor1b=1,flavors)
4003     END DO
4004 !================================================
4005   endif
4006 !================================================
4007 
4008   CALL GreenHyboffdiag_destroy(F_tmp)
4009 
4010   !FREE(F_tmp%oper_w)
4011 !  ==============================
4012 END SUBROUTINE Ctqmcoffdiag_getGreen

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_getResult [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_getResult

FUNCTION

  reduce everything to get the result of the simulation

COPYRIGHT

  Copyright (C) 2013-2024 ABINIT group (J. Bieder,F. Gendron)
  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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

2992 SUBROUTINE Ctqmcoffdiag_getResult(op,Iatom,fname)
2993 
2994 
2995 #ifdef HAVE_MPI1
2996 include 'mpif.h'
2997 #endif
2998 !Arguments ------------------------------------
2999   TYPE(Ctqmcoffdiag)  , INTENT(INOUT)                    :: op
3000   INTEGER, INTENT(IN ) :: Iatom
3001   character(len=fnlen), INTENT(INOUT) :: fname
3002 !Local variables ------------------------------
3003   INTEGER                                       :: iflavor
3004 !  INTEGER                                       :: iflavor1
3005 !  INTEGER                                       :: iflavor2
3006   INTEGER                                       :: flavors
3007   INTEGER                                       :: itau
3008   INTEGER                                       :: endDensity
3009   character(len=2)                              :: atomnb
3010   DOUBLE PRECISION                              :: inv_flavors
3011   DOUBLE PRECISION                              :: a
3012   DOUBLE PRECISION                              :: b
3013   DOUBLE PRECISION                              :: r
3014   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: alpha
3015   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: beta
3016   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)   :: measN_1
3017   DOUBLE PRECISION,              DIMENSION(1:2) :: TabX
3018   DOUBLE PRECISION,              DIMENSION(1:2) :: TabY
3019   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)   :: freqs
3020   INTEGER, ALLOCATABLE, DIMENSION(:)   :: counts
3021   INTEGER, ALLOCATABLE, DIMENSION(:)   :: displs
3022   INTEGER, ALLOCATABLE, DIMENSION(:)   :: occtot
3023   INTEGER, ALLOCATABLE, DIMENSION(:)   :: spintot
3024   INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: occ
3025   INTEGER                                       :: sp1,spinmax,spinmin,dspin,nelec,spin
3026   INTEGER                                       :: spAll
3027   INTEGER                                       :: last
3028   INTEGER                                       :: n1
3029   INTEGER                                       :: n2,n3,quotient,remainder,signe
3030   INTEGER                                       :: debut
3031   DOUBLE PRECISION                                       :: signvaluemeassum
3032 !  INTEGER                                       :: fin
3033 #ifdef HAVE_MPI
3034   INTEGER                                       :: ierr
3035   DOUBLE PRECISION,              DIMENSION(1)   :: arr
3036 #endif
3037   INTEGER                                       :: sizeoper,nbprocs,myrank
3038   DOUBLE PRECISION                              :: inv_size,sumh,sumtot
3039   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: buffer 
3040   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: buffer2,buffer2s
3041   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: fullempty
3042   TYPE(FFTHyb) :: FFTmrka
3043 
3044 #if defined HAVE_MPI && !defined HAVE_MPI2_INPLACE
3045    DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:)   :: buffer1_out,freqs_buf
3046    DOUBLE PRECISION, ALLOCATABLE , DIMENSION(:,:) :: buffer2_out
3047 #endif
3048 
3049   IF ( .NOT. op%done ) &
3050     CALL ERROR("Ctqmcoffdiag_getResult : Simulation not run                ")
3051 
3052   flavors     =  op%flavors
3053   inv_flavors = 1.d0 / DBLE(flavors)
3054 
3055 
3056   inv_size = 1.d0 / DBLE(op%size)
3057   sp1 = 0
3058   spAll = 0
3059 
3060 !#ifdef CTCtqmcoffdiag_CHECK
3061   IF ( op%opt_check .GT. 0 ) THEN
3062     op%errorImpurity = ImpurityOperator_getError(op%Impurity) * inv_flavors 
3063     op%errorBath     = BathOperatoroffdiag_getError    (op%Bath    ) * inv_flavors 
3064   END IF
3065 !#endif
3066 
3067   MALLOC(alpha,(1,1))
3068   MALLOC(beta,(1,1))
3069   MALLOC(buffer,(1,1))
3070   IF ( op%opt_noise .EQ. 1) THEN
3071     FREEIF(alpha)
3072     MALLOC(alpha,(1:op%samples+1,1:flavors))
3073     FREEIF(beta)
3074     MALLOC(beta,(1:op%samples+1,1:flavors))
3075   END IF
3076 
3077   IF ( op%have_MPI .EQV. .TRUE.) THEN 
3078     sp1   = 0
3079     spAll = sp1 + flavors + 6 
3080 
3081 !#ifdef CTCtqmcoffdiag_ANALYSIS
3082     IF ( op%opt_analysis .EQ. 1 ) &
3083       spAll = spAll + 3*sp1 
3084     IF ( op%opt_order .GT. 0 ) &
3085       spAll = spAll + op%opt_order 
3086     IF ( op%opt_noise .EQ. 1 ) &
3087       spAll = spAll + 2*(op%samples + 1)
3088 !#endif
3089 
3090     FREEIF(buffer)
3091     MALLOC(buffer,(1:spAll,1:MAX(2,flavors)))
3092   END IF
3093 
3094 !  op%seg_added    = op%seg_added    * inv_flavors 
3095 !  op%seg_removed  = op%seg_removed  * inv_flavors
3096 !  op%seg_sign     = op%seg_sign     * inv_flavors
3097 !  op%anti_added   = op%anti_added   * inv_flavors
3098 !  op%anti_removed = op%anti_removed * inv_flavors
3099 !  op%anti_sign    = op%anti_sign    * inv_flavors
3100   op%stats(:) = op%stats(:) * inv_flavors
3101 
3102   DO iflavor = 1, flavors
3103     ! Accumulate last values of  N (see also ctqmc_measn)
3104     op%measN(1,iflavor) = op%measN(1,iflavor) + op%measN(3,iflavor)*op%measN(4,iflavor)
3105     op%measN(2,iflavor) = op%measN(2,iflavor) + op%measN(4,iflavor)
3106     ! Reduction
3107     op%measN(1,iflavor)  = op%measN(1,iflavor) / ( op%measN(2,iflavor) * op%beta )
3108     ! Correction
3109 !#ifdef CTCtqmcoffdiag_ANALYSIS
3110     IF ( op%opt_order .GT. 0 ) &
3111       op%measPerturbation(:   ,iflavor) = op%measPerturbation(:,iflavor) &
3112                                     / SUM(op%measPerturbation(:,iflavor))
3113     IF ( op%opt_order .GT. 0 ) &
3114       op%meas_fullemptylines(:   ,iflavor) = op%meas_fullemptylines(:,iflavor) &
3115                                     / SUM(op%meas_fullemptylines(:,iflavor))
3116     !write(6,*) "sum fullempty",iflavor,op%meas_fullemptylines(:,iflavor)
3117 
3118     IF ( op%opt_analysis .EQ. 1 ) THEN
3119       op%measCorrelation (:,1,iflavor) = op%measCorrelation  (:,1,iflavor) &
3120                                     / SUM(op%measCorrelation (:,1,iflavor)) &
3121                                     * op%inv_dt 
3122       op%measCorrelation (:,2,iflavor) = op%measCorrelation  (:,2,iflavor) &
3123                                     / SUM(op%measCorrelation (:,2,iflavor)) &
3124                                     * op%inv_dt 
3125       op%measCorrelation (:,3,iflavor) = op%measCorrelation  (:,3,iflavor) &
3126                                     / SUM(op%measCorrelation (:,3,iflavor)) &
3127                                     * op%inv_dt 
3128     END IF
3129 !#endif
3130     IF ( op%opt_noise .EQ. 1 ) THEN
3131       TabX(1) = DBLE(op%modNoise2)
3132       TabX(2) = DBLE(op%modNoise1)
3133       DO itau = 1, op%samples+1
3134         op%measNoiseG(itau,iflavor,2)%vec = -op%measNoiseG(itau,iflavor,2)%vec*op%inv_dt &  
3135                                            /(op%beta*DBLE(op%modNoise2))
3136         op%measNoiseG(itau,iflavor,1)%vec = -op%measNoiseG(itau,iflavor,1)%vec*op%inv_dt &  
3137                                            /(op%beta*DBLE(op%modNoise1))
3138         n2 = op%measNoiseG(itau,iflavor,2)%tail
3139         TabY(1) = Stat_deviation(op%measNoiseG(itau,iflavor,2)%vec(1:n2))!*SQRT(n2/(n2-1))
3140         n1 = op%measNoiseG(itau,iflavor,1)%tail
3141         TabY(2) = Stat_deviation(op%measNoiseG(itau,iflavor,1)%vec(1:n1))!*SQRT(n1/(n1-1))
3142         CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,alpha(itau,iflavor),beta(itau,iflavor),r)
3143         ! ecart type -> 60%
3144         ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma
3145       END DO
3146     END IF
3147 
3148   END DO
3149 !sui!write(6,*) "getresults"
3150   CALL GreenHyboffdiag_measHybrid(op%Greens, op%Bath%M, op%Impurity%Particles, .TRUE.,op%signvalue)
3151   CALL GreenHyboffdiag_getHybrid(op%Greens)
3152  ! write(6,*) "op%measN",op%measN(1,:)
3153   MALLOC(measN_1,(flavors))
3154   do iflavor=1,flavors
3155     measN_1(iflavor)=op%measN(1,iflavor)
3156   enddo
3157   CALL GreenHyboffdiag_setN(op%Greens, measN_1(:))
3158   FREE(measN_1)
3159 
3160 ! todoab case _nd and _d are not completely described.
3161   FREEIF(buffer2)
3162   FREEIF(buffer2s)
3163   sizeoper=size(op%Greens%oper,1)
3164   !write(6,*) "sss",size(op%Greens%oper,1),sizeoper
3165   !write(6,*) "sss",size(op%Greens%oper,2),flavors
3166   !write(6,*) "sss",size(op%Greens%oper,3),flavors
3167   MALLOC(buffer2,(1:sizeoper,flavors,flavors))
3168   MALLOC(buffer2s,(1:sizeoper,flavors,flavors))
3169   MALLOC(fullempty,(2,flavors))
3170       !sui!write(6,*) "greens1"
3171   IF ( op%have_MPI .EQV. .TRUE. ) THEN 
3172       !sui!write(6,*) "greens2"
3173     fullempty=0.d0
3174     buffer2 = op%Greens%oper
3175     !write(6,*) "buffer2",(op%Greens%oper(1,n1,n1),n1=1,flavors)
3176     buffer2s= 0.d0
3177     do iflavor=1,flavors
3178       do itau=1,sizeoper
3179     !sui!write(6,*) "greens",iflavor,itau,op%Greens%oper(itau,iflavor,iflavor)
3180       enddo
3181     enddo
3182    !write(6,*) "beforempi",op%Greens%oper(1,1,1) ,buffer2(1,1,1)
3183 #ifdef HAVE_MPI
3184    CALL MPI_COMM_SIZE(op%MY_COMM,nbprocs,ierr)
3185    CALL MPI_COMM_RANK(op%MY_COMM,myrank,ierr)
3186 #endif
3187   !write(6,*) "procs",nbprocs,myrank
3188   END IF
3189   last = sp1
3190 
3191   op%measDE(:,:) = op%measDE(:,:) * DBLE(op%measurements) /(DBLE(op%sweeps)*op%beta)
3192 
3193   IF ( op%opt_histo .GT. 0 ) THEN
3194     op%occup_histo_time(:) = op%occup_histo_time(:) / INT(op%sweeps/op%measurements)
3195     op%occupconfig(:) = op%occupconfig(:) / INT(op%sweeps/op%measurements)
3196     op%suscep(:,:) = op%suscep(:,:) / INT(op%sweeps/op%measurements)
3197     op%chi(:,:) = op%chi(:,:) / INT(op%sweeps/op%measurements)
3198     op%chicharge(:,:) = op%chicharge(:,:) / INT(op%sweeps/op%measurements)
3199     op%ntot(:) = op%ntot(:) / INT(op%sweeps/op%measurements)
3200   END IF
3201 ! HISTO before MPI_SUM
3202 !  write(6,*) "=== Histogram of occupations for complete simulation 3 ====",INT(op%sweeps/op%measurements)
3203 !  sumh=0
3204 !  do n1=1,op%flavors+1
3205 !     write(6,'(i4,f10.4)')  n1-1, op%occup_histo_time(n1)
3206 !     sumh=sumh+op%occup_histo_time(n1)
3207 !  enddo
3208 !  write(6,*) "=================================",sumh
3209 
3210   n1 = op%measNoise(1)%tail
3211   n2 = op%measNoise(2)%tail
3212 
3213   ! On utilise freqs comme tableau de regroupement
3214   ! Gather de Noise1
3215   IF ( op%have_MPI .EQV. .TRUE. ) THEN
3216     MALLOC(counts,(1:op%size))
3217     MALLOC(displs,(1:op%size))
3218     FREEIF(freqs)
3219     MALLOC(freqs,(1:op%size*n1))
3220     freqs = 0.d0
3221     counts(:) = n1
3222     displs(:) = (/ ( iflavor*n1, iflavor=0, op%size-1 ) /)
3223 #ifdef HAVE_MPI
3224 #if defined HAVE_MPI2_INPLACE
3225     freqs(n1*op%rank+1:n1*(op%rank+1)) = op%measNoise(1)%vec(1:n1) 
3226     CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
3227                         freqs, counts, displs, &
3228                         MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3229 #else
3230     MALLOC(freqs_buf,(n1))
3231     freqs_buf(1:n2)=op%measNoise(2)%vec(1:n1)
3232     CALL MPI_ALLGATHERV(freqs_buf, n1, MPI_DOUBLE_PRECISION, &
3233                         freqs, counts, displs, &
3234                         MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3235     FREE(freqs_buf)
3236 #endif
3237 #endif
3238     n1 = op%size*n1
3239     CALL Vector_setSize(op%measNoise(1),n1)
3240     op%measNoise(1)%vec(1:n1) = freqs(:)
3241     ! Gather de Noise2
3242     FREE(freqs)
3243     MALLOC(freqs,(1:op%size*n2))
3244     freqs = 0.d0
3245     counts(:) = n2
3246     displs(:) = (/ ( iflavor*n2, iflavor=0, op%size-1 ) /)
3247 #ifdef HAVE_MPI
3248 #if defined HAVE_MPI2_INPLACE
3249     freqs(n2*op%rank+1:n2*(op%rank+1)) = op%measNoise(2)%vec(1:n2) 
3250     CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
3251                         freqs, counts, displs, &
3252                         MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3253 #else
3254     MALLOC(freqs_buf,(n2))
3255     freqs_buf(1:n2)=op%measNoise(2)%vec(1:n2)
3256     CALL MPI_ALLGATHERV(freqs_buf, n2, MPI_DOUBLE_PRECISION, &
3257                         freqs, counts, displs, &
3258                         MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
3259     FREE(freqs_buf)
3260 #endif
3261 #endif
3262     n2 = op%size*n2
3263     CALL Vector_setSize(op%measNoise(2),n2)
3264     op%measNoise(2)%vec(1:n2) = freqs(:)
3265     FREE(counts)
3266     FREE(displs)
3267     FREE(freqs)
3268   END IF
3269   !n1 = op%measNoise(1)%tail
3270   !n2 = op%measNoise(2)%tail
3271 
3272   ! Transformation des paquets pour que ca fit a CTQMC_SLICE(1|2)
3273   IF ( n1 .GT. CTQMC_SLICE1 ) THEN
3274     itau = n1/CTQMC_SLICE1
3275     MALLOC(freqs,(1:n1/itau))
3276     DO debut=1, n1/itau
3277       freqs(debut)=SUM(op%measNoise(1)%vec((debut-1)*itau+1:itau*debut))
3278     END DO
3279     freqs(:) = freqs(:)/DBLE(itau)
3280     op%modNoise1 = op%modNoise1*itau
3281     n1 = n1/itau
3282     CALL Vector_setSize(op%measNoise(1),n1)
3283     op%measNoise(1)%vec(1:n1) = freqs(:)
3284     FREE(freqs)
3285   END IF
3286   IF ( n2 .GT. CTQMC_SLICE1*CTQMC_SLICE2 ) THEN
3287     itau = n2/(CTQMC_SLICE1*CTQMC_SLICE2)
3288     MALLOC(freqs,(1:n2/itau))
3289     DO debut=1, n2/itau
3290       freqs(debut)=SUM(op%measNoise(2)%vec((debut-1)*itau+1:itau*debut))
3291     END DO
3292     freqs(:) = freqs(:)/DBLE(itau)
3293     op%modNoise2 = op%modNoise2*itau
3294     n2 = n2/itau
3295     CALL Vector_setSize(op%measNoise(2),n2)
3296     op%measNoise(2)%vec(1:n2) = freqs(:)
3297     FREE(freqs)
3298   END IF
3299   ! On peut s'amuser avec nos valeur d'energies
3300   !MALLOC(TabX,(1:20))
3301   !MALLOC(TabY,(1:20))
3302 
3303   TabX(1) = DBLE(op%modNoise2)
3304   TabX(2) = DBLE(op%modNoise1)
3305 
3306   ! Il faut calculer pour chaque modulo 10 ecarts type sur les donnes acquises
3307   op%measNoise(1)%vec(1:n1) = op%measNoise(1)%vec(1:n1)/(op%beta*DBLE(op%modNoise1))*DBLE(op%measurements)
3308   op%measNoise(2)%vec(1:n2) = op%measNoise(2)%vec(1:n2)/(op%beta*DBLE(op%modNoise2))*DBLE(op%measurements)
3309 !  CALL Vector_print(op%measNoise(1),op%rank+70)
3310 !  CALL Vector_print(op%measNoise(2),op%rank+50)
3311 !  DO iflavor=1,10
3312 !    debut = (iflavor-1)*n2/10+1
3313 !    fin   = iflavor*n2/10
3314 !    TabY(iflavor) = Stat_deviation(op%measNoise(2)%vec(debut:fin))
3315 !    debut = (iflavor-1)*n1/10+1
3316 !    fin   = iflavor*n1/10
3317 !    TabY(10+iflavor) = Stat_deviation(op%measNoise(1)%vec(debut:fin))
3318 !  END DO
3319 !!  TabY(1:n) = (op%measNoise(2)%vec(1:n)   &
3320 !!              )
3321 !!             !/(op%beta*DBLE(op%modNoise2))*DBLE(op%measurements) &
3322 !!             !- op%measDE(1,1))
3323 !!  TabY(op%measNoise(2)%tail+1:n+op%measNoise(2)%tail) = (op%measNoise(1)%vec(1:n)   &
3324 !!               )
3325 !!             ! /(op%beta*DBLE(op%modNoise1))*DBLE(op%measurements) &
3326 !!             ! - op%measDE(1,1))
3327 !  IF ( op%rank .EQ. 0 ) THEN
3328 !    DO iflavor=1,20
3329 !      write(45,*) TabX(iflavor), TabY(iflavor)
3330 !    END DO
3331 !  END IF
3332 !
3333 
3334 
3335   TabY(1) = Stat_deviation(op%measNoise(2)%vec(1:n2))!*SQRT(n2/(n2-1))
3336 !!  write(op%rank+10,*) TabX(2)
3337 !!  write(op%rank+40,*) TabX(1)
3338 !!  CALL Vector_print(op%measNoise(1),op%rank+10)
3339 !!  CALL Vector_print(op%measNoise(2),op%rank+40)
3340 !!  CLOSE(op%rank+10)
3341 !!  CLOSE(op%rank+40)
3342   TabY(2) = Stat_deviation(op%measNoise(1)%vec(1:n1))!*SQRT(n1/(n1-1))
3343 !!  ! Ecart carre moyen ~ ecart type mais non biaise. Serait moins precis. Aucun
3344   ! impact sur la pente, juste sur l'ordonnee a l'origine.
3345 
3346   CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,a,b,r)
3347 !  FREE(TabX)
3348 !  FREE(TabY)
3349   ! ecart type -> 60%
3350   ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma
3351 
3352   !op%measDE(1,1) = SUM(op%measNoise(1)%vec(1:op%measNoise(1)%tail))/(DBLE(op%measNoise(1)%tail*op%modNoise1)*op%beta)
3353   !op%measDE(2:flavors,1:flavors) = op%measDE(2:flavors,1:flavors) /(DBLE(op%sweeps)*op%beta)
3354   CALL ImpurityOperator_getErrorOverlap(op%Impurity,op%measDE)
3355   ! Add the difference between true calculation and quick calculation of the
3356   ! last sweep overlap to measDE(2,2)
3357   !op%measDE = op%measDE * DBLE(op%measurements) 
3358   IF ( op%have_MPI .EQV. .TRUE. ) THEN 
3359     IF ( op%opt_analysis .EQ. 1 ) THEN
3360       buffer(last+1:last+sp1,:) = op%measCorrelation(:,1,:)
3361       last = last + sp1
3362       buffer(last+1:last+sp1,:) = op%measCorrelation(:,2,:)
3363       last = last + sp1
3364       buffer(last+1:last+sp1,:) = op%measCorrelation(:,3,:)
3365       last = last + sp1
3366     END IF
3367     IF ( op%opt_order .GT. 0 ) THEN
3368       buffer(last+1:last+op%opt_order, :) = op%measPerturbation(:,:)
3369       last = last + op%opt_order
3370     END IF
3371     IF ( op%opt_noise .EQ. 1 ) THEN
3372       buffer(last+1:last+op%samples+1,:) = alpha(:,:)
3373       last = last + op%samples + 1
3374       buffer(last+1:last+op%samples+1,:) = beta(:,:)
3375       last = last + op%samples + 1
3376     END IF
3377 !  op%measDE(2,2) = a*EXP(b*LOG(DBLE(op%sweeps*op%size)))
3378     buffer(spall-(flavors+5):spAll-6,:) = op%measDE(:,:)
3379 !    buffer(spAll  ,1) = op%seg_added   
3380 !    buffer(spAll-1,1) = op%seg_removed 
3381 !    buffer(spAll-2,1) = op%seg_sign    
3382 !    buffer(spAll  ,2) = op%anti_added  
3383 !    buffer(spAll-1,2) = op%anti_removed
3384 !    buffer(spAll-2,2) = op%anti_sign   
3385     buffer(spAll  ,1) = op%stats(1)
3386     buffer(spAll-1,1) = op%stats(2)
3387     buffer(spAll-2,1) = op%stats(3)
3388     buffer(spAll  ,2) = op%stats(4)
3389     buffer(spAll-1,2) = op%stats(5)
3390     buffer(spAll-2,2) = op%stats(6)
3391     buffer(spAll-3,1) = op%swap
3392     buffer(spAll-3,2) = DBLE(op%modGlobalMove(2))
3393     buffer(spAll-4,1) = a
3394     buffer(spAll-4,2) = b
3395 !#ifdef CTCtqmcoffdiag_CHECK
3396     buffer(spAll-5,1) = op%errorImpurity
3397     buffer(spAll-5,2) = op%errorBath 
3398     signvaluemeassum = 0
3399 !#endif
3400 
3401 #ifdef HAVE_MPI
3402 #if defined HAVE_MPI2_INPLACE
3403     CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, spAll*flavors, &
3404                      MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3405 #else
3406     MALLOC(buffer2_out,(spAll,flavors))
3407     CALL MPI_ALLREDUCE(buffer, buffer2_out, spAll*flavors, &
3408                      MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3409     buffer(1:spAll,1:flavors)=buffer2_out(1:spAll,1:flavors)
3410     FREE(buffer2_out)
3411 #endif
3412     CALL MPI_ALLREDUCE(buffer2, buffer2s, sizeoper*flavors*flavors, &
3413                      MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3414     CALL MPI_ALLREDUCE([op%runtime], arr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, op%MY_COMM, ierr)
3415     op%runtime=arr(1)
3416     CALL MPI_ALLREDUCE([op%Greens%signvaluemeas], arr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3417     signvaluemeassum=arr(1)
3418 #if defined HAVE_MPI2_INPLACE
3419     IF ( op%opt_histo .GT. 0 ) THEN
3420       CALL MPI_ALLREDUCE(MPI_IN_PLACE,op%occup_histo_time, flavors+1, MPI_DOUBLE_PRECISION, MPI_SUM, &
3421              op%MY_COMM, ierr)
3422       CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%occupconfig, 2**flavors, MPI_DOUBLE_PRECISION, MPI_SUM, &
3423                op%MY_COMM, ierr)
3424       CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%suscep, 3*op%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
3425                op%MY_COMM, ierr)
3426       CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%chi, 3*op%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
3427                op%MY_COMM, ierr)
3428       CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%chicharge, 3*op%samples, MPI_DOUBLE_PRECISION, MPI_SUM, &
3429                op%MY_COMM, ierr)
3430       CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%ntot, 3, MPI_DOUBLE_PRECISION, MPI_SUM, &
3431                op%MY_COMM, ierr)
3432     END IF
3433     CALL MPI_ALLREDUCE(MPI_IN_PLACE, sumh, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
3434               op%MY_COMM, ierr)     
3435 #else
3436       MALLOC(buffer1_out,(flavors+1))
3437       CALL MPI_ALLREDUCE(op%occup_histo_time, buffer1_out, flavors+1, MPI_DOUBLE_PRECISION, MPI_SUM, &
3438              op%MY_COMM, ierr)
3439       op%occup_histo_time(1:flavors+1) =buffer1_out(1:flavors+1)
3440       FREE(buffer1_out)
3441 #endif
3442     CALL MPI_ALLREDUCE([sumh], arr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
3443     sumh=arr(1)
3444     IF ( op%opt_order .GT. 0 ) THEN
3445       CALL MPI_ALLREDUCE(op%meas_fullemptylines, fullempty, 2*flavors, MPI_DOUBLE_PRECISION, MPI_SUM, &
3446                op%MY_COMM, ierr)
3447     ENDIF
3448 #endif
3449 
3450   
3451     buffer          = buffer * inv_size
3452     op%measDE(:,:)  = buffer(spall-(flavors+5):spAll-6,:)
3453 !    op%seg_added    = buffer(spAll  ,1)
3454 !    op%seg_removed  = buffer(spAll-1,1)
3455 !    op%seg_sign     = buffer(spAll-2,1)
3456 !    op%anti_added   = buffer(spAll  ,2)
3457 !    op%anti_removed = buffer(spAll-1,2)
3458 !    op%anti_sign    = buffer(spAll-2,2)
3459     op%stats(1)    = buffer(spAll  ,1)
3460     op%stats(2)    = buffer(spAll-1,1)
3461     op%stats(3)    = buffer(spAll-2,1)
3462     op%stats(4)    = buffer(spAll  ,2)
3463     op%stats(5)    = buffer(spAll-1,2)
3464     op%stats(6)    = buffer(spAll-2,2)
3465     op%swap         = buffer(spAll-3,1)
3466     op%modGlobalMove(2) = NINT(buffer(spAll-3,2))
3467     a               = buffer(spAll-4,1) 
3468     b               = buffer(spAll-4,2)
3469 !!#ifdef CTCtqmcoffdiag_CHECK
3470     op%errorImpurity= buffer(spAll-5,1) 
3471     op%errorBath    = buffer(spAll-5,2)   
3472 !#endif
3473 
3474    ! DO iflavor = 1, flavors
3475    !   op%Greens(iflavor)%oper          = buffer(1:sp1          , iflavor)
3476    ! END DO
3477     op%Greens%oper = buffer2s/float(nbprocs)
3478    ! write(6,*) "buffer2s",(op%Greens%oper(1,n1,n1),n1=1,flavors)
3479     op%Greens%signvaluemeas = signvaluemeassum/float(nbprocs)
3480     !sui!write(6,*) "nbprocs",nbprocs,op%Greens%signvaluemeas
3481     op%Greens%oper = op%Greens%oper / op%Greens%signvaluemeas
3482    ! write(6,*) "buffer3s",(op%Greens%oper(1,n1,n1),n1=1,flavors)
3483     IF ( op%opt_order .GT. 0 ) THEN
3484       op%meas_fullemptylines= fullempty/float(nbprocs)
3485     ENDIF
3486     do iflavor=1,flavors
3487       do itau=1,sizeoper
3488     !sui!write(6,*) "greens_av",iflavor,itau,op%Greens%oper(itau,iflavor,iflavor)
3489       enddo
3490     enddo
3491    !write(6,*) "aftermpi",op%Greens%oper(1,1,1) ,buffer2s(1,1,1)
3492     last = sp1
3493     IF ( op%opt_analysis .EQ. 1 ) THEN
3494       op%measCorrelation(:,1,:) = buffer(last+1:last+sp1,:) 
3495       last = last + sp1
3496       op%measCorrelation(:,2,:) = buffer(last+1:last+sp1,:) 
3497       last = last + sp1
3498       op%measCorrelation(:,3,:) = buffer(last+1:last+sp1,:) 
3499       last = last + sp1
3500     END IF
3501     IF ( op%opt_order .GT. 0 ) THEN
3502       op%measPerturbation(:,:) = buffer(last+1:last+op%opt_order, :)
3503       last = last + op%opt_order
3504     END IF
3505     IF ( op%opt_noise .EQ. 1 ) THEN
3506       alpha(:,:) = buffer(last+1:last+op%samples+1,:)
3507       last = last + op%samples + 1
3508       beta(:,:) = buffer(last+1:last+op%samples+1,:)
3509       last = last + op%samples + 1
3510     END IF
3511   END IF
3512   DO iflavor = 1, flavors
3513     ! complete DE matrix
3514     op%measDE(iflavor, iflavor+1:flavors) = op%measDE(iflavor+1:flavors,iflavor)
3515   END DO
3516   FREE(buffer)
3517   FREE(buffer2)
3518   FREE(buffer2s)
3519   FREE(fullempty)
3520 
3521   IF ( op%opt_spectra .GE. 1 ) THEN
3522     endDensity = SIZE(op%density,2)
3523     IF ( op%density(1,endDensity) .EQ. -1.d0 ) &
3524       endDensity = endDensity - 1
3525     CALL FFTHyb_init(FFTmrka,endDensity,DBLE(op%thermalization)/DBLE(op%measurements*op%opt_spectra))
3526     ! Not very Beauty 
3527     MALLOC(freqs,(1:FFTmrka%size/2))
3528     DO iflavor = 1, flavors
3529       ! mean value is removed to supress the continue composent 
3530       CALL FFTHyb_setData(FFTmrka,op%density(iflavor,1:endDensity)/op%beta+op%Greens%oper(op%samples+1,iflavor,iflavor))
3531       CALL FFTHyb_run(FFTmrka,1)
3532       CALL FFTHyb_getData(FFTmrka,endDensity,op%density(iflavor,:),freqs)
3533     END DO
3534     op%density(flavors+1,:) = -1.d0
3535     op%density(flavors+1,1:FFTmrka%size/2) = freqs
3536     CALL FFTHyb_destroy(FFTmrka)
3537     FREE(freqs)
3538   END IF
3539 
3540   op%a_Noise = a
3541   op%b_Noise = b
3542   IF ( op%opt_noise .EQ. 1 ) THEN
3543     op%abNoiseG(1,:,:) = alpha
3544     op%abNoiseG(2,:,:) = beta
3545   END IF
3546   FREE(alpha)
3547   FREE(beta)
3548   IF ( op%opt_histo .GT. 0 ) THEN
3549     write(op%ostream,*) "=== Histogram of occupations for complete simulation  ===="
3550  !   write(6,*) "sumh over procs", sumh
3551     sumh=0
3552     do n1=1,op%flavors+1
3553        write(op%ostream,'(i4,f10.4)')  n1-1, op%occup_histo_time(n1)/float(nbprocs)
3554        sumh=sumh+op%occup_histo_time(n1)/float(nbprocs)
3555     enddo
3556        write(op%ostream,'(a,f10.4)') " all" , sumh
3557     write(op%ostream,*) "================================="
3558 
3559     MALLOC(occ,(2**op%flavors,1:flavors))
3560     MALLOC(occtot,(2**op%flavors))
3561     MALLOC(spintot,(2**op%flavors))
3562     do n1=1,2**op%flavors
3563       ! Compute occupations of individual Orbitals
3564       n3=n1-1
3565       occtot(n1)=0
3566       spintot(n1)=0
3567       signe=1
3568       do n2=1,op%flavors
3569         remainder=modulo(n3,2)
3570         quotient=(n3-remainder)/2
3571         occ(n1,n2)=remainder
3572         n3=quotient
3573         occtot(n1)=occtot(n1)+occ(n1,n2)
3574         if(n2>=7) signe =-1
3575         !if(n2>=7) signe =0
3576         spintot(n1)=spintot(n1)+occ(n1,n2)*signe
3577       enddo
3578       op%occupconfig(n1)=op%occupconfig(n1)/float(nbprocs)
3579     enddo
3580 
3581     write(op%ostream,*) "=== Histogram of occupations of configurations for complete simulation  ===="
3582     sumh=0
3583     if(op%flavors==14) then
3584       do n1=1,2**op%flavors
3585          write(op%ostream,'(i4,14i2,f20.2)')  n1, (occ(n1,n2),n2=1,op%flavors),op%occupconfig(n1)
3586          sumh=sumh+op%occupconfig(n1)
3587       enddo
3588     else if(op%flavors==10) then
3589       do n1=1,2**op%flavors
3590          write(op%ostream,'(i4,10i2,f20.2)')  n1, (occ(n1,n2),n2=1,op%flavors),op%occupconfig(n1)
3591          sumh=sumh+op%occupconfig(n1)
3592       enddo
3593     end if
3594     write(op%ostream,'(a,f10.4)') " all" , sumh
3595     
3596     sumtot=0
3597     do nelec=0,10
3598       spinmin=modulo(nelec,2)
3599       if(nelec<=5) spinmax=nelec
3600       if(nelec>=6) spinmax=10-nelec
3601       spinmax=6
3602       dspin=1
3603       write(op%ostream,*) "=== Histogram of occupations of configurations for total number of electrons",nelec
3604       do spin=spinmin,spinmax,dspin
3605         sumh=0
3606         do n1=1,2**op%flavors
3607           if(occtot(n1)==nelec.and.abs(spintot(n1))==spin) then
3608             sumh=sumh+op%occupconfig(n1)
3609             if(op%flavors==14) then
3610               write(op%ostream,'(i8,14i2,a,i2,i3,f10.4)')  n1,(occ(n1,n2),n2=1,op%flavors),"  ",occtot(n1),spintot(n1),&
3611 & op%occupconfig(n1)
3612             else if(op%flavors==10) then
3613               write(op%ostream,'(i8,10i2,a,i2,i3,f10.4)')  n1,(occ(n1,n2),n2=1,op%flavors),"  ",occtot(n1),spintot(n1),&
3614  op%occupconfig(n1)
3615             end if
3616           endif
3617         enddo
3618         write(op%ostream,'(a,i4,a,i4,a,f10.4)') " === Sum of weights for",nelec," electrons and spin",spin," is ",sumh
3619         sumtot=sumtot+sumh
3620       enddo
3621     enddo
3622     write(op%ostream,'(a,f10.4)') "Full sum is",sumtot
3623 
3624    FREE(occ)
3625     FREE(occtot)
3626     FREE(spintot)
3627   !END IF
3628 
3629     !==============================
3630     ! Write Susceptibilities
3631     !==============================
3632     write(atomnb, '("0",i1)') Iatom
3633     ! Local Magnetic Susceptibility
3634     if(op%opt_histo .gt. 1) then
3635       ! Scalar
3636       if(op%nspinor .eq. 1) then
3637         open (unit=735,file=trim(fname)//'_LocalSpinSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
3638         write(735,*) '#Tau Total t2g eg'
3639         do n1=1,op%samples
3640           op%suscep(:,n1)=op%suscep(:,n1)/float(nbprocs)/float(op%samples)
3641           write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*op%beta/op%samples,(op%suscep(n2,n1),n2=1,3)
3642         enddo
3643 
3644       else
3645         ! SOC
3646         open (unit=735,file=trim(fname)//'_LocalMagnSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
3647         write(735,*) '#Tau Total Orbital Spin'
3648         do n1=1,op%samples
3649           op%chi(:,n1) = op%chi(:,n1)/float(nbprocs)/float(op%samples)
3650           write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*op%beta/op%samples,(op%chi(n2,n1),n2=1,3)
3651         end do
3652       endif
3653     close(unit=735)
3654     endif
3655 
3656     ! Local Charge Susceptiblity
3657     if(op%opt_histo .gt. 2) then 
3658       op%ntot(:) = op%ntot(:)/float(nbprocs)/float(op%samples)
3659       open (unit=735,file=trim(fname)//'_LocalChargeSuscept_atom_'//atomnb//'.dat',status='unknown',form='formatted')
3660       write(735,*) '#Tau Total <ntot>'
3661       do n1=1,op%samples
3662         op%chicharge(1,n1)=(op%chicharge(1,n1)/float(nbprocs)/float(op%samples))-(op%ntot(1)*op%ntot(1))
3663         !op%chicharge(2,n1)=(op%chicharge(2,n1)/float(nbprocs)/float(op%samples))-(op%ntot(2)*op%ntot(2))
3664         !op%chicharge(3,n1)=(op%chicharge(3,n1)/float(nbprocs)/float(op%samples))-(op%ntot(3)*op%ntot(3))
3665         write(735,'(1x,f14.8,2x,f12.8,2x,f12.8,2x,f12.8)') (n1-1)*op%beta/op%samples,(op%chicharge(1,n1)),op%ntot(1)
3666       enddo
3667       close(unit=735)
3668     endif
3669 
3670   END IF
3671 
3672 
3673 END SUBROUTINE Ctqmcoffdiag_getResult

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_init [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_init

FUNCTION

  Initialize the type Ctqmcoffdiag
  Allocate all the non optional variables

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  ostream=where to write
  istream=where to read the input parameters if so
  bFile=logical argument True if input is read from istream
  MY_COMM=mpi communicator for the CTQMC
  iBuffer=input parameters if bFile is false

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

344 SUBROUTINE Ctqmcoffdiag_init(op, ostream, istream, bFile, MY_COMM, iBuffer)
345 
346 
347 #ifdef HAVE_MPI1
348 include 'mpif.h'
349 #endif
350 !Arguments ------------------------------------
351   TYPE(Ctqmcoffdiag), INTENT(INOUT)                      :: op
352   INTEGER  , INTENT(IN   )                      :: ostream
353   INTEGER  , INTENT(IN   )                      :: istream
354   LOGICAL  , INTENT(IN   )                      :: bFile
355   DOUBLE PRECISION, DIMENSION(1:11), OPTIONAL, INTENT(IN) :: iBuffer
356   INTEGER  , OPTIONAL, INTENT(IN   )                      :: MY_COMM
357 !Local variables ------------------------------
358 #ifdef HAVE_MPI
359   INTEGER                                       :: ierr
360 #endif
361   !INTEGER                                       :: iflavor
362 #ifdef __GFORTRAN__
363 !  INTEGER                                       :: pid
364 !  CHARACTER(LEN=5)                              :: Cpid
365 !
366 #endif
367   DOUBLE PRECISION, DIMENSION(1:11)             :: buffer
368 
369   op%ostream = ostream
370   op%istream = istream
371   
372 ! --- RENICE ---
373 !#ifdef __GFORTRAN__
374 !  pid = GetPid()
375 !  WRITE(Cpid,'(I5)') pid
376 !  CALL SYSTEM('renice +19 '//TRIM(ADJUSTL(Cpid))//' > /dev/null')
377 !#endif
378 !! --- RENICE ---
379 
380   IF ( PRESENT(MY_COMM)) THEN
381 #ifdef HAVE_MPI
382     op%have_MPI = .TRUE.
383     op%MY_COMM = MY_COMM
384     CALL MPI_Comm_rank(op%MY_COMM, op%rank, ierr)
385     CALL MPI_Comm_size(op%MY_COMM, op%size, ierr)
386 #else
387     CALL WARN("MPI is not used                                    ")
388     op%have_MPI = .FALSE.
389     op%MY_COMM = -1
390     op%rank = 0
391     op%size = 1
392 #endif
393   ELSE
394     op%have_MPI = .FALSE.
395     op%MY_COMM = -1
396     op%rank = 0
397     op%size = 1
398   END IF
399 
400   !IF ( op%rank .EQ. 0 ) THEN
401 !  WRITE(ostream,'(A20)') 'Job reniced with +19'
402     !CALL FLUSH(ostream)
403   !END IF
404 
405   IF ( bFile .EQV. .TRUE. ) THEN
406     IF ( op%rank .EQ. 0 ) THEN
407 
408       READ(istream,*) buffer(1) !iseed
409       READ(istream,*) buffer(2) !op%sweeps
410       READ(istream,*) buffer(3) !op%thermalization
411       READ(istream,*) buffer(4) !op%measurements
412       READ(istream,*) buffer(5) !op%flavors
413       READ(istream,*) buffer(6) !op%samples
414       READ(istream,*) buffer(7) !op%beta
415       READ(istream,*) buffer(8) !U
416       READ(istream,*) buffer(9) !iTech
417       READ(istream,*) buffer(11)!op%nspinor
418       !READ(istream,*) buffer(9) !Wmax
419 !#ifdef CTCtqmcoffdiag_ANALYSIS
420       !READ(istream,*) buffer(10) !order
421 !#endif
422     END IF
423 
424 #ifdef HAVE_MPI
425     IF ( op%have_MPI .EQV. .TRUE. ) &
426       CALL MPI_Bcast(buffer, 11, MPI_DOUBLE_PRECISION, 0,    &
427                    op%MY_COMM, ierr)
428 #endif
429   ELSE IF ( PRESENT(iBuffer) ) THEN
430     buffer(1:11) = iBuffer(1:11)
431   ELSE
432     CALL ERROR("Ctqmcoffdiag_init : No input parameters                    ")
433   END IF
434 
435   CALL Ctqmcoffdiag_setParameters(op, buffer)
436 
437   CALL Ctqmcoffdiag_allocateAll(op)
438 
439   CALL GreenHyboffdiag_init(op%Greens,op%samples, op%beta,INT(buffer(5)), &
440   iTech=INT(buffer(9)),MY_COMM=op%MY_COMM)
441 
442 
443 !  op%seg_added    = 0.d0
444 !  op%anti_added   = 0.d0
445 !  op%seg_removed  = 0.d0
446 !  op%anti_removed = 0.d0
447 !  op%seg_sign     = 0.d0
448 !  op%anti_sign    = 0.d0
449   op%stats(:)     = 0.d0
450  ! write(std_out,*) "op%stats",op%stats
451   op%signvalue    = 1.d0
452 !  op%signvaluecurrent    = 0.d0
453 !  op%signvaluemeas = 0.d0
454   op%swap         = 0.d0
455   op%runTime      = 0.d0
456 
457   CALL Vector_init(op%measNoise(1),op%sweeps/op%modNoise1)
458   CALL Vector_init(op%measNoise(2),(op%sweeps/op%modNoise1+1)*CTQMC_SLICE2)
459   !CALL Vector_init(op%measNoise(3),101)
460   !CALL Vector_init(op%measNoise(4),101)
461 
462   op%set  = op%para .AND. op%inF
463   op%done = .FALSE.
464   op%init = .TRUE.
465 
466 !#ifdef CTCtqmcoffdiag_CHECK
467   op%errorImpurity = 0.d0
468   op%errorBath     = 0.d0
469 !#endif
470 END SUBROUTINE Ctqmcoffdiag_init

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_loop [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_loop

FUNCTION

  Definition the main loop of the CT-QMC

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  itotal=number of sweeps to perform : thermalization or sweeps
  ilatex=unit of file to write movie if so

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

1925 SUBROUTINE Ctqmcoffdiag_loop(op,itotal,ilatex)
1926 
1927 !Arguments ------------------------------------
1928   TYPE(Ctqmcoffdiag), INTENT(INOUT)         :: op
1929   INTEGER    , INTENT(IN   )         :: itotal
1930   INTEGER    , INTENT(IN   )         :: ilatex
1931 !Local variables ------------------------------
1932   LOGICAL                            :: updated 
1933   LOGICAL                            :: updated_seg
1934   LOGICAL, DIMENSION(:), ALLOCATABLE :: updated_swap
1935 
1936   INTEGER                            :: flavors
1937   INTEGER                            :: measurements
1938   INTEGER                            :: modNoise1
1939   INTEGER                            :: modNoise2
1940   INTEGER                            :: modGlobalMove
1941   INTEGER                            :: sp1
1942   INTEGER                            :: itau   
1943   INTEGER                            :: ind
1944   INTEGER                            :: endDensity
1945   INTEGER                            :: indDensity
1946   INTEGER                            :: swapUpdate1
1947   INTEGER                            :: swapUpdate2
1948   INTEGER                            :: old_percent
1949   INTEGER                            :: new_percent
1950   INTEGER                            :: ipercent !,ii
1951   INTEGER                            :: iflavor,ifl1,iflavor_d
1952   INTEGER                            :: isweep
1953 
1954   DOUBLE PRECISION                   :: cpu_time1
1955   DOUBLE PRECISION                   :: cpu_time2
1956   DOUBLE PRECISION                   :: NRJ_old1
1957   DOUBLE PRECISION                   :: NRJ_old2
1958   DOUBLE PRECISION                   :: NRJ_new
1959   DOUBLE PRECISION                   :: total
1960   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old1
1961   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old2
1962   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_new
1963 
1964   CALL CPU_TIME(cpu_time1)
1965 
1966   flavors        = op%flavors
1967   measurements   = op%measurements
1968   modNoise1      = op%modNoise1
1969   modNoise2      = op%modNoise2
1970   modGlobalMove  = op%modGlobalMove(1)
1971   sp1            = op%samples+1
1972   IF ( op%opt_histo .GT. 0 ) THEN
1973     op%occup_histo_time= 0.d0
1974     op%occupconfig= 0.d0
1975     op%suscep= 0.d0
1976     op%chi= 0.d0
1977     op%chicharge= 0.d0
1978     op%ntot= 0.d0
1979   END IF
1980 
1981   old_percent    = 0
1982 
1983   MALLOC(updated_swap,(1:flavors))
1984   updated_swap(:) = .FALSE.
1985 
1986   NRJ_old1  = 0.d0
1987   NRJ_old2  = 0.d0
1988   NRJ_new   = 0.d0
1989 
1990   MALLOC(gtmp_new,(1,1))
1991   gtmp_new  = 0.d0
1992   MALLOC(gtmp_old1,(1,1))
1993   gtmp_old1 = 0.d0
1994   MALLOC(gtmp_old2,(1,1))
1995   gtmp_old2 = 0.d0
1996 
1997   endDensity = SIZE(op%density,2)
1998 
1999   IF ( op%opt_noise .GT. 0 ) THEN
2000     FREEIF(gtmp_new)
2001     MALLOC(gtmp_new,(1:sp1,1:flavors))
2002     FREEIF(gtmp_old1)
2003     MALLOC(gtmp_old1,(1:sp1,1:flavors))
2004     FREEIF(gtmp_old2)
2005     MALLOC(gtmp_old2,(1:sp1,1:flavors))
2006   END IF
2007 
2008   IF ( op%rank .EQ. 0 ) THEN
2009     WRITE(op%ostream, '(1x,103A)') &
2010     "|----------------------------------------------------------------------------------------------------|"
2011     WRITE(op%ostream,'(1x,A)', ADVANCE="NO") "|"
2012   END IF
2013 
2014   total = DBLE(itotal)
2015   !write(std_out,*) "itotal",itotal
2016   indDensity = 1
2017   !write(std_out,*) "op%stats",op%stats
2018   DO isweep = 1, itotal
2019   !ii if(op%prtopt==1) write(std_out,*) "======== Isweep = ",isweep
2020     !updated_seg=.FALSE.
2021     DO iflavor = 1, flavors
2022      ! if(isweep==itotal) write(std_out,*) "    Iflavor = ",iflavor,op%Impurity%Particles(iflavor)%tail
2023    !ii if(op%prtopt==1)  write(std_out,*) "      ===Iflavor = ",iflavor
2024       op%Impurity%activeFlavor=iflavor
2025       op%Bath%activeFlavor=iflavor ; op%Bath%MAddFlag= .FALSE. ; op%Bath%MRemoveFlag = .FALSE.
2026 
2027       !write(std_out,*) "before tryaddremove"
2028 
2029       ! For iflavor, Try a move
2030       !==========================
2031       CALL Ctqmcoffdiag_tryAddRemove(op,updated_seg)
2032     !sui!write(std_out,*) "after tryaddremove",updated_seg
2033 
2034       updated = updated_seg .OR.  updated_swap(iflavor).OR.(isweep==1)
2035       updated_swap(iflavor) = .FALSE.
2036       if ( op%opt_nondiag >0 )  iflavor_d=0
2037       if ( op%opt_nondiag==0 )  iflavor_d=iflavor
2038       CALL GreenHyboffdiag_measHybrid(op%Greens, op%Bath%M, op%Impurity%Particles, updated,op%signvalue,iflavor_d) 
2039 
2040       CALL Ctqmcoffdiag_measN        (op, iflavor, updated)
2041       IF ( op%opt_analysis .EQ. 1 ) &
2042         CALL Ctqmcoffdiag_measCorrelation (op, iflavor)
2043       IF ( op%opt_order .GT. 0 ) &
2044         CALL Ctqmcoffdiag_measPerturbation(op, iflavor)
2045     END DO
2046     !CALL GreenHyboffdiag_measHybrid(op%Greens, op%Bath%M, op%Impurity%Particles, updated,op%signvalue,iflavor_d) 
2047     !DO iflavor = 1,flavors
2048     !  CALL Ctqmcoffdiag_measN        (op, iflavor, updated)
2049     !END DO
2050 
2051     IF ( MOD(isweep,modGlobalMove) .EQ. 0 ) THEN
2052   ! !sui!write(std_out,*) "isweep,modGlobalMove,inside",isweep,modGlobalMove
2053       CALL Ctqmcoffdiag_trySwap(op,swapUpdate1, swapUpdate2)
2054      ! !write(std_out,*) "no global move yet for non diag hybridization"
2055       IF ( swapUpdate1 .NE. 0 .AND. swapUpdate2 .NE. 0 ) THEN
2056         updated_swap(swapUpdate1) = .TRUE.
2057         updated_swap(swapUpdate2) = .TRUE.
2058       END IF
2059     END IF
2060     
2061     IF ( MOD(isweep,measurements) .EQ. 0 ) THEN ! default is always 
2062       CALL ImpurityOperator_measDE(op%Impurity,op%measDE)
2063       IF ( op%opt_spectra .GE. 1 .AND. MOD(isweep,measurements*op%opt_spectra) .EQ. 0 ) THEN
2064         op%density(1:flavors,indDensity) = op%measN(3,1:flavors)
2065         indDensity = indDensity+1
2066       END IF
2067     END IF
2068 
2069     IF ( MOD(isweep,measurements) .EQ. 0 ) THEN
2070       IF ( op%opt_histo .GT. 0 ) THEN
2071         CALL ImpurityOperator_occup_histo_time(op%Impurity,op%occup_histo_time,op%occupconfig,op%suscep,op%samples,op%chi,&
2072 & op%chicharge,op%ntot,op%opt_histo,op%nspinor)
2073       END IF
2074     ENDIF
2075 
2076     IF ( MOD(isweep, modNoise1) .EQ. 0 ) THEN
2077       !modNext = isweep + modNoise2
2078       NRJ_new = op%measDE(1,1)
2079       CALL Vector_pushBack(op%measNoise(1),NRJ_new - NRJ_old1)
2080       NRJ_old1 = NRJ_new
2081 
2082       !! Try to limit accumulation error
2083       CALL ImpurityOperator_cleanOverlaps(op%Impurity)
2084 
2085       IF ( op%opt_noise .EQ. 1 ) THEN
2086         DO ifl1 = 1, flavors
2087           DO ind = 1, op%Greens%map(ifl1,ifl1)%tail
2088             itau = op%Greens%map(ifl1,ifl1)%listINT(ind)
2089             gtmp_new(itau,ifl1) = op%Greens%oper(itau,ifl1,ifl1) & 
2090                      +op%Greens%map(ifl1,ifl1)%listDBLE(ind)*DBLE(op%Greens%factor)
2091           END DO
2092           DO itau = 1, sp1
2093            CALL Vector_pushBack(op%measNoiseG(itau,ifl1,1), gtmp_new(itau,ifl1) - gtmp_old1(itau,ifl1))
2094            gtmp_old1(itau,ifl1) = gtmp_new(itau,ifl1)
2095           END DO
2096         END DO
2097       END IF
2098     END IF
2099 
2100     IF ( MOD(isweep,modNoise2) .EQ. 0 ) THEN
2101       NRJ_new = op%measDE(1,1)
2102       CALL Vector_pushBack(op%measNoise(2),NRJ_new - NRJ_old2)
2103       NRJ_old2 = NRJ_new
2104       IF ( op%opt_noise .EQ. 1 ) THEN
2105         DO ifl1 = 1, flavors
2106           DO ind = 1, op%Greens%map(ifl1,ifl1)%tail
2107             itau = op%Greens%map(ifl1,ifl1)%listINT(ind)
2108             gtmp_new(itau,ifl1) = op%Greens%oper(itau,ifl1,ifl1) & 
2109                   +op%Greens%map(ifl1,ifl1)%listDBLE(ind)*op%Greens%factor
2110           END DO
2111           DO itau = 1, sp1
2112             CALL Vector_pushBack(op%measNoiseG(itau,ifl1,2), gtmp_new(itau,ifl1) - gtmp_old2(itau,ifl1))
2113             gtmp_old2(itau,ifl1) = gtmp_new(itau,ifl1)
2114           END DO
2115         END DO 
2116       END IF
2117 
2118       IF ( op%rank .EQ. 0 ) THEN 
2119         new_percent = CEILING(DBLE(isweep)*100.d0/DBLE(itotal))
2120         DO ipercent = old_percent+1, new_percent 
2121           WRITE(op%ostream,'(A)',ADVANCE="NO") "-"
2122         END DO
2123         old_percent = new_percent
2124       END IF
2125     END IF
2126 
2127     IF ( op%opt_movie .EQ. 1 ) THEN
2128       WRITE(ilatex,'(A11,I9)') "%iteration ", isweep
2129       CALL ImpurityOperator_printLatex(op%Impurity,ilatex,isweep)
2130     END IF
2131 
2132   END DO
2133 
2134   IF ( op%rank .EQ. 0 ) THEN
2135     DO ipercent = old_percent+1, 100
2136       WRITE(op%ostream,'(A)',ADVANCE="NO") "-"
2137     END DO
2138     WRITE(op%ostream,'(A)') "|"
2139   END IF
2140  
2141   FREE(gtmp_new)
2142   FREE(gtmp_old1)
2143   FREE(gtmp_old2)
2144   FREE(updated_swap)
2145 
2146   IF ( op%opt_spectra .GE. 1 .AND. itotal .EQ. op%sweeps ) THEN
2147     IF ( endDensity .NE. indDensity-1 ) THEN
2148       op%density(:,endDensity) = -1.d0
2149     END IF
2150   END IF
2151 
2152   CALL CPU_TIME(cpu_time2)
2153 
2154   op%runTime = (cpu_time2 - cpu_time1)*1.05d0 ! facteur arbitraire de correction
2155 END SUBROUTINE Ctqmcoffdiag_loop

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_measCorrelation [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_measCorrelation

FUNCTION

  measure all correlations in times for a flavor

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  iflavor=the flavor to measure

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

2855 SUBROUTINE Ctqmcoffdiag_measCorrelation(op, iflavor)
2856 
2857 !Arguments ------------------------------------
2858   TYPE(Ctqmcoffdiag)             , INTENT(INOUT)       :: op
2859   !TYPE(ImpurityOperator), INTENT(IN   )       :: impurity
2860   INTEGER               , INTENT(IN   )       :: iflavor
2861 !Local variables ------------------------------
2862   INTEGER                                     :: iCdag
2863   INTEGER                                     :: iCdagBeta
2864   INTEGER                                     :: iC
2865   INTEGER                                     :: index
2866   INTEGER                                     :: size
2867   DOUBLE PRECISION                            :: tC
2868   DOUBLE PRECISION                            :: tCdag
2869   !DOUBLE PRECISION                            :: time
2870   DOUBLE PRECISION                            :: inv_dt
2871   DOUBLE PRECISION                            :: beta
2872 
2873   IF ( .NOT. op%set ) &
2874     CALL ERROR("Ctqmcoffdiag_measCorrelation : QMC not set                 ")
2875     !write(6,*) "not available"
2876     stop
2877 
2878   size = op%impurity%particles(op%impurity%activeFlavor)%tail
2879   beta = op%beta
2880 
2881   IF ( size .EQ. 0 ) RETURN
2882   
2883   inv_dt = op%inv_dt
2884 
2885   DO iCdag = 1, size ! first segments
2886     tCdag  = op%impurity%particles(op%impurity%activeFlavor)%list(iCdag,Cdag_)
2887     tC     = op%impurity%particles(op%impurity%activeFlavor)%list(iCdag,C_   )
2888     index = INT( ( (tC - tCdag)  * inv_dt ) + .5d0 ) + 1
2889     op%measCorrelation(index,1,iflavor) = op%measCorrelation(index,1,iflavor) + 1.d0
2890     MODCYCLE(iCdag+1,size,iCdagBeta)
2891     index = INT( ( ( &
2892                     op%impurity%particles(op%impurity%activeFlavor)%list(iCdagBeta,Cdag_) - tC &
2893                     + AINT(DBLE(iCdag)/DBLE(size))*beta &
2894                    )  * inv_dt ) + .5d0 ) + 1
2895     IF ( index .LT. 1 .OR. index .GT. op%samples+1 ) THEN
2896       CALL WARN("Ctqmcoffdiag_measCorrelation : bad index line 1095         ")
2897     ELSE
2898       op%measCorrelation(index,2,iflavor) = op%measCorrelation(index,2,iflavor) + 1.d0
2899     END IF
2900 !    DO iC = 1, size
2901 !      tC = impurity%particles(impurity%activeFlavor)%list(C_,iC)
2902 !      time = tC - tCdag
2903 !      IF ( time .LT. 0.d0 ) time = time + beta
2904 !      index = INT( ( time * inv_dt ) + .5d0 ) + 1
2905 !      op%measCorrelation(index,3,iflavor) = op%measCorrelation(index,3,iflavor) + 1.d0
2906 !    END DO
2907     DO iC = 1, size!  op%Greens(iflavor)%index_old%tail 
2908 !todoba        op%measCorrelation(op%Greens(iflavor)%map%listINT(iC+(iCdag-1)*size),3,iflavor) = &
2909 !todoba        op%measCorrelation(op%Greens(iflavor)%map%listINT(iC+(iCdag-1)*size),3,iflavor) + 1.d0
2910     END DO
2911   END DO
2912 
2913 END SUBROUTINE Ctqmcoffdiag_measCorrelation

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_measN [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_measN

FUNCTION

  measures the number of electron
  by taking into account the value for the move before before this one
  with the correct weight.

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  iflavor=which flavor to measure
  updated=something has changed since last time

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

2791 SUBROUTINE Ctqmcoffdiag_measN(op, iflavor, updated)
2792 
2793 !Arguments ------------------------------------
2794   TYPE(Ctqmcoffdiag)             , INTENT(INOUT)     :: op
2795   !TYPE(ImpurityOperator), INTENT(IN   )     :: impurity
2796   INTEGER               , INTENT(IN   )     :: iflavor
2797   LOGICAL               , INTENT(IN   )     :: updated
2798 
2799 !  IF ( .NOT. op%set ) &
2800 !    CALL ERROR("Ctqmcoffdiag_measN : QMC not set                           ")
2801 
2802   
2803   IF ( updated .EQV. .TRUE. ) THEN
2804 !  --- accumulate occupations with values op%measN(3,iflavor) from the last measurements with the corresponding weight
2805 !  ---  op*measN(4,iflavor)
2806     op%measN(1,iflavor) = op%measN(1,iflavor) + op%measN(3,iflavor)*op%measN(4,iflavor)
2807    ! write(6,*) "Cllll42"
2808 
2809 !  --- Compute total number of new measurements 
2810     op%measN(2,iflavor) = op%measN(2,iflavor) + op%measN(4,iflavor)
2811 
2812    ! write(6,*) "Allll42"
2813 !  --- Compute the occupation for this configuration (will be put in
2814 !  --- op%measN(1,iflavor) at the next occurence of updated=.true.), with
2815 !  --- the corresponding weight  op%measN(4,iflavor) (we do not now it yet)
2816     op%measN(3,iflavor) = ImpurityOperator_measN(op%impurity)
2817 
2818 !  --- set weight: as update=true, it is a new measurement , so put it to one
2819     op%measN(4,iflavor) = 1.d0
2820 
2821   ELSE
2822 !  --- increased the count so that at new move, we will be able to update measN(1) correctly.
2823     op%measN(4,iflavor) = op%measN(4,iflavor) + 1.d0
2824    ! write(6,*) "Bllll42"
2825   END IF
2826 END SUBROUTINE Ctqmcoffdiag_measN

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_measPerturbation [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_measPerturbation

FUNCTION

  measure perturbation order

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  iflavor=the flavor to measure

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

2941 SUBROUTINE Ctqmcoffdiag_measPerturbation(op, iflavor)
2942 
2943 !Arguments ------------------------------------
2944   TYPE(Ctqmcoffdiag)             , INTENT(INOUT)     :: op
2945   !TYPE(ImpurityOperator), INTENT(IN   )     :: impurity
2946   INTEGER               , INTENT(IN   )     :: iflavor
2947 !Local variables ------------------------------
2948   INTEGER                                   :: index
2949 
2950   IF ( .NOT. op%set ) &
2951     CALL ERROR("Ctqmcoffdiag_measiPerturbation : QMC not set               ")
2952 
2953   index = op%impurity%particles(op%impurity%activeFlavor)%tail + 1
2954   IF ( index .LE. op%opt_order ) &
2955     op%measPerturbation(index,iflavor) = op%measPerturbation(index,iflavor) + 1.d0
2956   IF ( index == 1 ) THEN
2957     IF (op%impurity%particles(iflavor)%list(0,C_) < op%impurity%particles(iflavor)%list(0,Cdag_) ) THEN
2958       op%meas_fullemptylines(1,iflavor) = op%meas_fullemptylines(1,iflavor) + 1.d0
2959     ELSE
2960       op%meas_fullemptylines(2,iflavor) = op%meas_fullemptylines(2,iflavor) + 1.d0
2961     ENDIF
2962   ENDIF
2963 
2964 END SUBROUTINE Ctqmcoffdiag_measPerturbation

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printAll [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printAll

FUNCTION

  print different functions computed during the simulation

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4131 SUBROUTINE Ctqmcoffdiag_printAll(op)
4132 
4133 !Arguments ------------------------------------
4134   TYPE(Ctqmcoffdiag), INTENT(INOUT) :: op
4135 
4136   IF ( .NOT. op%done ) &
4137     CALL WARNALL("Ctqmcoffdiag_printAll : Simulation not run                 ")
4138 
4139 !sui!write(6,*) "op%stats",op%stats
4140   CALL Ctqmcoffdiag_printQMC(op)
4141 
4142   CALL Ctqmcoffdiag_printGreen(op)
4143 
4144   CALL Ctqmcoffdiag_printD(op)
4145 
4146 !  CALL Ctqmcoffdiag_printE(op)
4147 
4148 !#ifdef CTCtqmcoffdiag_ANALYSIS
4149   CALL Ctqmcoffdiag_printPerturbation(op)
4150 
4151   CALL Ctqmcoffdiag_printCorrelation(op)
4152 !#endif
4153 
4154   CALL Ctqmcoffdiag_printSpectra(op)
4155 
4156 END SUBROUTINE Ctqmcoffdiag_printAll

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printCorrelation [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printCorrelation

FUNCTION

  print correlation fonctions

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  oFileIn=file stream

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4634 SUBROUTINE Ctqmcoffdiag_printCorrelation(op, oFileIn)
4635 
4636 !Arguments ------------------------------------
4637   TYPE(Ctqmcoffdiag)          , INTENT(IN)             :: op
4638   INTEGER  , OPTIONAL, INTENT(IN)             :: oFileIn
4639 !Local variables ------------------------------
4640   INTEGER                                     :: oFile
4641   INTEGER                                     :: itime
4642   INTEGER                                     :: sp1
4643   INTEGER                                     :: iflavor
4644   INTEGER                                     :: i
4645   INTEGER                                     :: flavors
4646   CHARACTER(LEN=2)                            :: a
4647   CHARACTER(LEN=50)                           :: string
4648   DOUBLE PRECISION                            :: dt
4649 
4650   !IF ( op%rank .NE. MOD(5,op%size)) RETURN
4651   IF ( op%rank .NE. MOD(op%size+5,op%size)) RETURN
4652   IF ( op%opt_analysis .NE. 1 ) RETURN
4653 
4654   oFile = 44
4655   IF ( PRESENT(oFileIn) ) THEN
4656     oFile = oFileIn
4657   ELSE
4658     OPEN(UNIT=oFile, FILE="Correlation.dat")
4659   END IF
4660 
4661   sp1         =  op%samples
4662   dt          =  op%beta / sp1
4663   sp1         =  sp1 + 1
4664   flavors     =  op%flavors
4665 
4666   i = 3*flavors + 1
4667   WRITE(a,'(I2)') i
4668   WRITE(oFile,*) "# time  (/ (segement, antiseg, correl), i=1, flavor/)"
4669   string = '(1x,'//TRIM(ADJUSTL(a))//'F19.15)'
4670   DO itime = 1, sp1
4671     WRITE(oFile,string) DBLE(itime-1)*dt, &
4672                    (/ ( &
4673                    (/ ( op%measCorrelation(itime, i, iflavor), i=1,3) /) &
4674                    , iflavor=1, flavors) /)
4675   END DO
4676 
4677   IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4678 
4679 END SUBROUTINE Ctqmcoffdiag_printCorrelation

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printD [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printD

FUNCTION

  print individual double occupancy

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  oFileIn=file stream

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4458 SUBROUTINE Ctqmcoffdiag_printD(op,oFileIn)
4459 
4460 !Arguments ------------------------------------
4461   TYPE(Ctqmcoffdiag)          , INTENT(IN)    :: op
4462   INTEGER  , OPTIONAL, INTENT(IN)    :: oFileIn
4463 !Local variables ------------------------------
4464   INTEGER                            :: oFile
4465   INTEGER                            :: iflavor1
4466   INTEGER                            :: iflavor2
4467 
4468   !IF ( op%rank .NE. MOD(2,op%size)) RETURN
4469   IF ( op%rank .NE. MOD(op%size+2,op%size)) RETURN
4470 
4471   oFile = 41
4472   IF ( PRESENT(oFileIn) ) THEN
4473     oFile = oFileIn
4474   ELSE
4475     OPEN(UNIT=oFile, FILE="D.dat")
4476   END IF
4477 
4478   DO iflavor1 = 1, op%flavors
4479     DO iflavor2 = iflavor1+1, op%flavors
4480       WRITE(oFile,'(1x,A8,I4,A1,I4,A3,ES21.14)') "Orbitals", iflavor1, "-", iflavor2, " : ", op%measDE(iflavor2,iflavor1)
4481     END DO
4482   END DO
4483 
4484   IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4485 
4486 END SUBROUTINE Ctqmcoffdiag_printD

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printE [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printE

FUNCTION

  print energy and noise 

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  oFileIn=file stream

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4514 SUBROUTINE Ctqmcoffdiag_printE(op,oFileIn)
4515 
4516 !Arguments ------------------------------------
4517   TYPE(Ctqmcoffdiag)          , INTENT(IN)    :: op
4518   INTEGER  , OPTIONAL, INTENT(IN)    :: oFileIn
4519 !Local variables ------------------------------
4520   INTEGER                            :: oFile
4521   DOUBLE PRECISION                   :: E
4522   DOUBLE PRECISION                   :: Noise
4523 
4524   !IF ( op%rank .NE. MOD(3,op%size)) RETURN
4525   IF ( op%rank .NE. MOD(op%size+3,op%size)) RETURN
4526 
4527   oFile = 42
4528   IF ( PRESENT(oFileIn) ) THEN
4529     oFile = oFileIn
4530   ELSE
4531     OPEN(UNIT=oFile, FILE="BetaENoise.dat")
4532   END IF
4533 
4534   CALL Ctqmcoffdiag_getE(op,E,Noise)
4535 
4536   WRITE(oFile,'(1x,F3.2,A2,ES21.14,A2,ES21.14)') op%beta, "  ", E, "  ",  Noise
4537 
4538   IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4539 
4540 END SUBROUTINE Ctqmcoffdiag_printE

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printGreen [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printGreen

FUNCTION

  print green functions

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  oFileIn=file stream

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4339 SUBROUTINE Ctqmcoffdiag_printGreen(op, oFileIn)
4340 
4341 !Arguments ------------------------------------
4342   use m_io_tools, only : flush_unit
4343   TYPE(Ctqmcoffdiag)        , INTENT(IN)    :: op
4344   INTEGER  , OPTIONAL, INTENT(IN)    :: oFileIn
4345 !Local variables ------------------------------
4346   INTEGER                            :: oFile
4347   INTEGER                            :: itime
4348   INTEGER                            :: sp1
4349   INTEGER                            :: iflavor,iflavorb
4350   INTEGER                            :: flavors !, iflavor2 !,iflavor1,
4351   CHARACTER(LEN=4)                   :: cflavors
4352   CHARACTER(LEN=50)                  :: string
4353   DOUBLE PRECISION                   :: dt
4354   DOUBLE PRECISION                   :: sweeps
4355 
4356   !IF ( op%rank .NE. MOD(1,op%size)) RETURN
4357   IF ( op%rank .NE. MOD(op%size+1,op%size)) RETURN
4358 
4359   oFile = 40
4360   IF ( PRESENT(oFileIn) ) THEN
4361     oFile = oFileIn
4362   ELSE
4363     OPEN(UNIT=oFile, FILE="Gtau.dat")
4364   END IF
4365   OPEN(UNIT=43, FILE="Gtau_nd.dat")
4366   rewind(43)
4367   sp1     =  op%samples
4368   dt      =  op%beta / DBLE(sp1)
4369   sp1     =  sp1 + 1
4370   flavors =  op%flavors
4371   sweeps = DBLE(op%sweeps)*DBLE(op%size)
4372 
4373   IF ( op%opt_noise .EQ. 1) THEN
4374     WRITE(cflavors,'(I4)') (2*flavors+1)*2
4375     string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
4376     DO itime = 1, sp1
4377       WRITE(oFile,string) DBLE(itime-1)*dt, &
4378       (/ (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors) /), &
4379       (/ (op%abNoiseG(1,itime,iflavor)*(sweeps)**op%abNoiseG(2,itime,iflavor), iflavor=1, flavors) /)
4380     END DO
4381   ELSE
4382     WRITE(cflavors,'(I4)') (flavors+1)*2
4383     string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
4384     DO itime = 1, sp1
4385 !     WRITE(45,string) DBLE(itime-1)*dt, &
4386 !     (/ (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors) /)
4387       WRITE(oFile,string) DBLE(itime-1)*dt, &
4388       (/ (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors) /)
4389 !     WRITE(46,*) DBLE(itime-1)*dt, &
4390 !     & (/ ((op%Greens%oper(itime,iflavor,iflavorb), iflavor=1, flavors),iflavorb=1,flavors) /)
4391     END DO
4392 !   DO itime = 1, sp1
4393 !     WRITE(47,*) DBLE(itime-1)*dt, &
4394 !     & (/ ((op%Greens%oper(itime,iflavor,iflavorb), iflavor=1, flavors),iflavorb=1,flavors) /)
4395 !   END DO
4396 !  --- Print full non diagonal Gtau in Gtau_nd.dat
4397     WRITE(cflavors,'(I4)') (flavors*flavors+1)
4398 !   write(47,*) "cflavors",cflavors
4399     string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
4400 !   write(47,*) string
4401     DO itime = 1, sp1
4402       WRITE(43,string) DBLE(itime-1)*dt, &
4403       & (/ ((op%Greens%oper(itime,iflavor,iflavorb), iflavorb=1, flavors),iflavor=1,flavors) /)
4404 !     WRITE(44,*) DBLE(itime-1)*dt, &
4405 !     & (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors)
4406 !     WRITE(44,string) DBLE(itime-1)*dt, &
4407 !     & (op%Greens%oper(itime,iflavor,iflavor), iflavor=1, flavors)
4408     END DO
4409       WRITE(43,*) 
4410   END IF
4411 !    DO iflavor = 1, flavors
4412 !      DO iflavor2 = 1, flavors
4413 !          write(4436,*) "#",iflavor,iflavor2
4414 !        do  itime=1,sp1
4415 !          write(4436,*) DBLE(itime-1)*dt,real(op%Greens%oper(itime,iflavor,iflavor2))
4416 !        enddo
4417 !          write(4436,*) 
4418 !      END DO
4419 !    END DO
4420 !    close(4436)
4421 
4422   IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4423   CLOSE(43)
4424 ! CLOSE(44)
4425 ! CLOSE(45)
4426 ! CLOSE(46)
4427 ! CLOSE(47)
4428   !call flush_unit(43)
4429 
4430 END SUBROUTINE Ctqmcoffdiag_printGreen

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printPerturbation [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printPerturbation

FUNCTION

  print perturbation order

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  oFileIn=file stream

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

4570 SUBROUTINE Ctqmcoffdiag_printPerturbation(op, oFileIn)
4571 
4572 !Arguments ------------------------------------
4573   TYPE(Ctqmcoffdiag)          , INTENT(IN)           :: op
4574   INTEGER  , OPTIONAL,  INTENT(IN)          :: oFileIn
4575 !Local variables-------------------------------
4576   INTEGER                                   :: oFile
4577   INTEGER                                   :: iorder
4578   INTEGER                                   :: order
4579   INTEGER                                   :: iflavor
4580   INTEGER                                   :: flavors
4581   CHARACTER(LEN=2)                          :: a
4582   CHARACTER(LEN=50)                         :: string
4583 
4584   !IF ( op%rank .NE. MOD(4,op%size)) RETURN
4585   IF ( op%rank .NE. MOD(op%size+4,op%size)) RETURN
4586   IF ( op%opt_order .LE. 0 ) RETURN
4587 
4588   oFile = 43
4589   IF ( PRESENT(oFileIn) ) THEN
4590     oFile = oFileIn
4591   ELSE
4592     OPEN(UNIT=oFile, FILE="Perturbation.dat")
4593   END IF
4594     
4595   order        =  op%opt_order
4596   flavors      =  op%flavors
4597 
4598   WRITE(a,'(I2)') flavors
4599   string = '(I5,'//TRIM(ADJUSTL(a))//'F19.15)'
4600   DO iorder = 1, order
4601     WRITE(oFile,string) iorder-1, &
4602                 (/ (op%measPerturbation(iorder, iflavor), iflavor=1, flavors) /)
4603   END DO
4604 
4605   IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4606 END SUBROUTINE Ctqmcoffdiag_printPerturbation

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printQMC [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printQMC

FUNCTION

  print ctqmc statistics

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4183 SUBROUTINE Ctqmcoffdiag_printQMC(op)
4184 
4185 !Arguments ------------------------------------
4186   TYPE(Ctqmcoffdiag), INTENT(INOUT) :: op
4187 !Local variables ------------------------------
4188   INTEGER                  :: ostream
4189   INTEGER                  :: iflavor,iflavorbis,iorder
4190   DOUBLE PRECISION         :: sweeps
4191   DOUBLE PRECISION         :: invSweeps
4192   CHARACTER(LEN=2)         :: a
4193   CHARACTER(LEN=15)        :: string
4194 
4195   !IF ( op%rank .NE. 0) RETURN
4196   IF ( op%rank .NE. MOD(op%size,op%size)) RETURN
4197 
4198   ostream   = op%ostream
4199   sweeps    = DBLE(op%sweeps)
4200   invSweeps = 1.d0/sweeps
4201 
4202   WRITE(ostream,'(1x,F13.0,A11,F10.2,A12,I5,A5)') sweeps*DBLE(op%size), " sweeps in ", op%runTime, &
4203                  " seconds on ", op%size, " CPUs"
4204   WRITE(ostream,'(A28,F6.2)') "Segments added        [%] : ", op%stats(4)*invSweeps*100.d0
4205   WRITE(ostream,'(A28,F6.2)') "Segments removed      [%] : ", op%stats(5)*invSweeps*100.d0
4206   WRITE(ostream,'(A28,F6.2)') "Segments <0 sign      [%] : ", op%stats(6)*invSweeps*100.d0
4207   !WRITE(ostream,'(A28,F12.2)') "Number of meas        [%] : ", op%stats(6)
4208   WRITE(ostream,'(A28,F6.2)') "Anti-segments added   [%] : ", op%stats(1)*invSweeps*100.d0
4209   WRITE(ostream,'(A28,F6.2)') "Anti-segments removed [%] : ", op%stats(2)*invSweeps*100.d0
4210   WRITE(ostream,'(A28,F6.2)') "Anti-segments <0 sign [%] : ", op%stats(3)*invSweeps*100.d0
4211   !WRITE(ostream,'(A28,F12.2)') "Sum of sign       [%] : ", op%stats(3)
4212   WRITE(ostream,'(A28,F13.2)') "Signe value               : ", op%Greens%signvaluemeas
4213   IF ( op%modGlobalMove(1) .LT. op%sweeps + 1 ) THEN
4214     WRITE(ostream,'(A28,F6.2)') "Global Move           [%] : ", op%swap         *invSweeps*100.d0*op%modGlobalMove(1)
4215     WRITE(ostream,'(A28,F6.2)') "Global Move Reduced   [%] : ", op%swap         / DBLE(op%modGlobalMove(2))*100.d0
4216   END IF
4217 !#ifdef CTCtqmcoffdiag_CHECK
4218   IF ( op%opt_check .EQ. 1 .OR. op%opt_check .EQ. 3 ) &
4219     WRITE(ostream,'(A28,E22.14)') "Impurity test         [%] : ", op%errorImpurity*100.d0
4220   IF ( op%opt_check .GE. 2 ) &
4221       WRITE(ostream,'(A28,E22.14)') "Bath     test         [%] : ", op%errorBath    *100.d0
4222 !#endif
4223   WRITE(ostream,'(A28,ES22.14,A5,ES21.14)') "<Epot>                [U] : ", op%measDE(1,1), " +/- ",&
4224 !#ifdef HAVE_MPI
4225                                                               op%a_Noise*(sweeps*DBLE(op%size))**op%b_Noise
4226 !#else
4227 !                                                              op%a_Noise*(sweeps)**op%b_Noise
4228 !#endif
4229  !--------- Write double occupation between all pairs of orbitals --------------------------
4230   write(ostream,'(17x,a)') "Double occupation between pairs of orbitals"
4231   write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4232   write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4233   do iflavor=1, op%flavors
4234     write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(op%measDE(iflavor,iflavorbis),iflavorbis=1,op%flavors)
4235   enddo
4236   write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
4237  !------------------------------------------------------------------------------------------
4238 
4239  !--------- Write number of segments for each orbitals
4240  ! write(ostream,'(a)') "Number of segments for each orbitals"
4241  ! write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4242  ! write(ostream,'(17x,30a)') ("----------",iflavorbis=1,op%flavors)
4243  ! do iflavor=1, op%flavors
4244  !   write(ostream,'(i17,a,30f10.4)') iflavor,"|",(op%Impurity%particles(IT)%tail
4245  ! enddo
4246  ! write(ostream,'(17x,30a)') ("----------",iflavorbis=1,op%flavors)
4247  !------------------------------------------------------------------------------------------
4248  !--------- Write G(L)
4249   write(ostream,'(17x,a)') "G(L)"
4250   write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4251   write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4252   do iflavor=1, op%flavors
4253     write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(op%Greens%oper(op%samples,iflavor,iflavorbis),iflavorbis=1,op%flavors)
4254   enddo
4255   write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
4256  !------------------------------------------------------------------------------------------
4257  !--------- Write G(1)
4258   write(ostream,'(17x,a)') "G(1)"
4259   write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4260   write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4261   do iflavor=1, op%flavors
4262     write(ostream,'(7x,i10,a,30f10.4)') iflavor,"|",(op%Greens%oper(1,iflavor,iflavorbis),iflavorbis=1,op%flavors)
4263   enddo
4264   write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors),ch10
4265  !------------------------------------------------------------------------------------------
4266 
4267   WRITE(ostream,'(A28,F8.4,A3,F7.4)') "Noise                 [U] : ", op%a_Noise, " x^", op%b_Noise
4268   WRITE(ostream,'(A28,E10.2)')  "Niquist puls.     [/beta] : ", ACOS(-1.d0)*op%inv_dt
4269   WRITE(ostream,'(A28,E22.14)') "Max Acc. Epot Error   [U] : ", op%measDE(2,2)/(op%beta*op%modNoise1*2.d0)*sweeps
4270   
4271   !WRITE(ostream,'(A28,F7.4,A3,F7.4,A4,E20.14)') "Noise            [G(tau)] : ", op%a_Noise(2), "x^", op%b_Noise(2), " -> ", &
4272                                                               !op%a_Noise(2)*(sweeps*DBLE(op%size))**op%b_Noise(2)
4273  !----- PERTURBATION ORDER------------------------------------------------------------------
4274   IF ( op%opt_order .GT. 0 ) THEN 
4275     write(ostream,*) 
4276     WRITE(a,'(I2)') op%flavors
4277     string = '(A28,'//TRIM(ADJUSTL(a))//'(1x,I3))'
4278     WRITE(ostream,string) "Perturbation orders       : ",(/ (MAXLOC(op%measPerturbation(:, iflavor))-1, iflavor=1, op%flavors) /)
4279     write(ostream,'(17x,a)') "order of Perturbation for flavors"
4280     write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4281     write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4282     write(ostream,'(12x,a,30i10)') " max ",(/ (MAXLOC(op%measPerturbation(:, iflavor))-1, iflavor=1, op%flavors) /)
4283     write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4284     do iorder=0, op%opt_order-1
4285       write(ostream,'(7x,i10,a,30f10.4)') iorder,"|",(op%measPerturbation(iorder+1,iflavor),iflavor=1,op%flavors)
4286     enddo
4287   END IF
4288  !------------------------------------------------------------------------------------------
4289  !----- PERTURBATION ORDER------------------------------------------------------------------
4290   IF ( op%opt_order .GT. 0 ) THEN 
4291     write(ostream,*) 
4292     write(ostream,'(17x,a)') "Proportion of full and empty orbital for order 0"
4293     write(ostream,'(17x,30i10)') (iflavorbis,iflavorbis=1,op%flavors)
4294     write(ostream,'(18x,30a)') ("----------",iflavorbis=1,op%flavors)
4295     write(ostream,'(2x,a,30f10.4)') " full  orbital |",(op%meas_fullemptylines(1,iflavor),iflavor=1,op%flavors)
4296     write(ostream,'(2x,a,30f10.4)') " empty orbital |",(op%meas_fullemptylines(2,iflavor),iflavor=1,op%flavors)
4297   END IF
4298  !------------------------------------------------------------------------------------------
4299   !CALL FLUSH(op%ostream)
4300   IF ( ABS(((op%stats(4) *invSweeps*100.d0) / (op%stats(5) *invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 &
4301    .OR. ABS(((op%stats(1)*invSweeps*100.d0) / (op%stats(2)*invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 ) &
4302     THEN 
4303     CALL WARNALL("Ctqmcoffdiag_printQMC : bad statistic according to moves. Increase sweeps")
4304   END IF
4305   IF ( ABS(op%b_Noise+0.5)/0.5d0 .GE. 0.05d0 ) &
4306     CALL WARNALL("Ctqmcoffdiag_printQMC : bad statistic according to Noise. Increase sweeps")
4307 !  IF ( ISNAN(op%a_Noise) .OR. ISNAN(op%a_Noise) ) &
4308 !    CALL WARNALL("Ctqmcoffdiag_printQMC : NaN appeared. Increase sweeps    ")
4309 
4310 
4311 END SUBROUTINE Ctqmcoffdiag_printQMC

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_printSpectra [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_printSpectra

FUNCTION

  print fourier transform of time evolution of number of electrons

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  oFileIn=file stream

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

4708 SUBROUTINE Ctqmcoffdiag_printSpectra(op, oFileIn)
4709 
4710 !Arguments ------------------------------------
4711   TYPE(Ctqmcoffdiag)          , INTENT(IN)             :: op
4712   INTEGER  , OPTIONAL, INTENT(IN)             :: oFileIn
4713 !Local variables ------------------------------
4714   INTEGER                                     :: oFile
4715   INTEGER                                     :: flavors
4716   INTEGER                                     :: indDensity
4717   INTEGER                                     :: endDensity
4718   CHARACTER(LEN=4)                            :: a
4719   CHARACTER(LEN=16)                           :: formatSpectra
4720 
4721   !IF ( op%rank .NE. MOD(6,op%size)) RETURN
4722   IF ( op%opt_spectra .LT. 1 ) RETURN
4723 
4724   oFile = 45+op%rank
4725   a ="0000"
4726   WRITE(a,'(I4)') op%rank
4727   IF ( PRESENT(oFileIn) ) THEN
4728     oFile = oFileIn
4729   ELSE
4730     OPEN(UNIT=oFile, FILE="Markov_"//TRIM(ADJUSTL(a))//".dat")
4731   END IF
4732 
4733   flavors     =  op%flavors
4734   WRITE(a,'(I4)') flavors+1
4735   formatSpectra ='(1x,'//TRIM(ADJUSTL(a))//'ES22.14)'
4736   WRITE(oFile,*) "# freq[/hermalization] FFT"
4737 
4738   endDensity = SIZE(op%density,2)
4739   DO WHILE ( op%density(flavors+1,endDensity) .EQ. -1 )
4740     endDensity = endDensity -1
4741   END DO
4742 
4743   DO indDensity = 1, endDensity
4744     WRITE(oFile,formatSpectra) op%density(flavors+1,indDensity), op%density(1:flavors,indDensity)
4745   END DO
4746 
4747   IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
4748 
4749 END SUBROUTINE Ctqmcoffdiag_printSpectra

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_reset [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_reset

FUNCTION

  reset a ctqmc simulation

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

1109 SUBROUTINE Ctqmcoffdiag_reset(op)
1110 
1111 !Arguments ------------------------------------
1112   TYPE(Ctqmcoffdiag), INTENT(INOUT) :: op
1113 !Local variables ------------------------------
1114   !INTEGER                  :: iflavor
1115   DOUBLE PRECISION         :: sweeps
1116 
1117   CALL GreenHyboffdiag_reset(op%Greens)
1118   CALL Ctqmcoffdiag_clear(op)
1119   CALL ImpurityOperator_reset(op%Impurity)
1120   CALL BathOperatoroffdiag_reset    (op%Bath)
1121   op%measN(3,:) = 0.d0
1122   !complete restart -> measN=0
1123   op%done = .FALSE.
1124   op%set  = .FALSE.
1125   op%inF  = .FALSE.
1126   op%opt_movie = 0
1127   op%opt_analysis = 0
1128   op%opt_order = 0
1129   op%opt_check = 0
1130   op%opt_noise = 0
1131   op%opt_spectra = 0
1132   op%opt_levels = 0
1133   sweeps = DBLE(op%sweeps)*DBLE(op%size)
1134   CALL Ctqmcoffdiag_setSweeps(op, sweeps)
1135 !#ifdef HAVE_MPI
1136 !  CALL MPI_BARRIER(op%MY_COMM,iflavor)
1137 !  IF ( op%rank .EQ. 0 ) &
1138 !#endif
1139 !  WRITE(op%ostream,'(A9)') "QMC reset"
1140 !  CALL FLUSH(op%ostream)
1141 END SUBROUTINE Ctqmcoffdiag_reset

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_run [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_run

FUNCTION

  set all options and run a simulation

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  opt_order=maximal perturbation order to scope
  opt_movie=draw a movie of the simulation
  opt_analysis=compute correlation functions
  opt_check=check fast calculations
  opt_noise=compute noise for green function
  opt_spectra=fourier transform of the time evolution of the number of electrons
  opt_gMove=steps without global move

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

1756 SUBROUTINE Ctqmcoffdiag_run(op,opt_order,opt_histo,opt_movie,opt_analysis,opt_check,opt_noise,opt_spectra,opt_gMove)
1757 
1758 
1759 #ifdef HAVE_MPI1
1760 include 'mpif.h'
1761 #endif
1762 !Arguments ------------------------------------
1763   TYPE(Ctqmcoffdiag), INTENT(INOUT)           :: op
1764   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_order
1765   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_histo
1766   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_movie
1767   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_analysis
1768   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_check
1769   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_noise
1770   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_spectra
1771   INTEGER, OPTIONAL, INTENT(IN   )  :: opt_gMove
1772 !Local variables ------------------------------
1773 #ifdef HAVE_MPI
1774   INTEGER                            :: ierr
1775   DOUBLE PRECISION                   :: rtime(1)
1776 #endif
1777 !#ifdef CTCtqmcoffdiag_MOVIE
1778   INTEGER                            :: ilatex
1779   CHARACTER(LEN=4)                   :: Cchar
1780 !#endif
1781   DOUBLE PRECISION                   :: estimatedTime
1782 
1783   IF ( .NOT. op%set  ) &
1784     CALL ERROR("Ctqmcoffdiag_run : QMC not set up                          ")
1785   IF ( .NOT. op%setU ) &
1786     CALL ERROR("Ctqmcoffdiag_run : QMC does not have a U matrix            ")
1787 
1788 
1789 ! OPTIONS of the run
1790   IF ( PRESENT( opt_check ) ) THEN
1791     op%opt_check = opt_check
1792     CALL ImpurityOperator_doCheck(op%Impurity,opt_check)
1793     CALL BathOperatoroffdiag_doCheck(op%Bath,opt_check)
1794   END IF
1795   IF ( PRESENT( opt_movie ) ) &
1796     op%opt_movie = opt_movie
1797   IF ( PRESENT( opt_analysis ) ) &
1798     op%opt_analysis = opt_analysis
1799   IF ( PRESENT ( opt_order ) ) &
1800     op%opt_order = opt_order 
1801   IF ( PRESENT ( opt_histo ) ) &
1802     op%opt_histo = opt_histo 
1803   IF ( PRESENT ( opt_noise ) ) THEN
1804     op%opt_noise = opt_noise 
1805   END IF
1806   IF ( PRESENT ( opt_spectra ) ) &
1807     op%opt_spectra = opt_spectra
1808 
1809   op%modGlobalMove(1) = max(op%sweeps,op%thermalization)+1 ! No Global Move
1810 !!sui!write(std_out,*) "op%sweeps",op%thermalization,op%sweeps,opt_gMove
1811   op%modGlobalMove(2) = 0
1812   IF ( PRESENT ( opt_gMove ) ) THEN
1813     IF ( opt_gMove .LE. 0 .OR. opt_gMove .GT. op%sweeps ) THEN
1814      ! op%modGlobalMove(1) = op%sweeps+1
1815       op%modGlobalMove(1) = max(op%sweeps,op%thermalization)+1 ! No Global Move
1816       !write(std_out,*) "op%sweeps",op%sweeps, op%modGlobalMove(1)
1817       CALL WARNALL("Ctqmcoffdiag_run : global moves option is <= 0 or > sweeps/cpu -> No global Moves")
1818     ELSE 
1819       op%modGlobalMove(1) = opt_gMove 
1820     END IF
1821   END IF
1822 !sui!write(std_out,*) "op%sweeps",op%thermalization,op%sweeps
1823 
1824   CALL Ctqmcoffdiag_allocateOpt(op)
1825   
1826 !#ifdef CTCtqmcoffdiag_MOVIE  
1827   ilatex = 0
1828   IF ( op%opt_movie .EQ. 1 ) THEN
1829     Cchar ="0000"
1830     WRITE(Cchar,'(I4)') op%rank 
1831     ilatex = 87+op%rank
1832     OPEN(UNIT=ilatex, FILE="Movie_"//TRIM(ADJUSTL(Cchar))//".tex")
1833     WRITE(ilatex,'(A)') "\documentclass{beamer}"
1834     WRITE(ilatex,'(A)') "\usepackage{color}"
1835     WRITE(ilatex,'(A)') "\setbeamersize{sidebar width left=0pt}"
1836     WRITE(ilatex,'(A)') "\setbeamersize{sidebar width right=0pt}"
1837     WRITE(ilatex,'(A)') "\setbeamersize{text width left=0pt}"
1838     WRITE(ilatex,'(A)') "\setbeamersize{text width right=0pt}"
1839     WRITE(ilatex,*) 
1840     WRITE(ilatex,'(A)') "\begin{document}"
1841     WRITE(ilatex,*) 
1842   END IF
1843 !#endif
1844 
1845   IF ( op%rank .EQ. 0 ) THEN
1846     WRITE(op%ostream,'(A29)') "Starting QMC (Thermalization)"
1847   END IF
1848   
1849   !=================================
1850   ! STARTING THERMALIZATION 
1851   !=================================
1852   !write(std_out,*) "sweeps before thermalization",op%sweeps
1853   !write(std_out,*) "op%stats",op%stats
1854   CALL Ctqmcoffdiag_loop(op,op%thermalization,ilatex)
1855   !=================================
1856   ! ENDING   THERMALIZATION 
1857   !=================================
1858 
1859   estimatedTime = op%runTime
1860 #ifdef HAVE_MPI
1861   CALL MPI_REDUCE([op%runTime], rtime, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, op%MY_COMM, ierr)
1862   estimatedTime=rtime(1)
1863 #endif
1864 
1865   IF ( op%rank .EQ. 0 ) THEN
1866     WRITE(op%ostream,'(A26,I6,A11)') "Thermalization done in    ", CEILING(estimatedTime), "    seconds"
1867     WRITE(op%ostream,'(A25,I7,A15,I5,A5)') "The QMC should run in    ", &
1868            CEILING(estimatedTime*DBLE(op%sweeps)/DBLE(op%thermalization)),&
1869                         "    seconds on ", op%size, " CPUs"
1870   END IF
1871 
1872   !=================================
1873   ! CLEANING CTQMC          
1874   !=================================
1875   CALL Ctqmcoffdiag_clear(op)
1876 
1877   !=================================
1878   ! STARTING CTQMC          
1879   !=================================
1880   !write(std_out,*) "sweeps before loop",op%sweeps
1881   !write(std_out,*) "op%stats",op%stats
1882   CALL Ctqmcoffdiag_loop(op,op%sweeps,ilatex)
1883   !=================================
1884   ! ENDING   CTQMC          
1885   !=================================
1886 
1887   IF ( op%opt_movie .EQ. 1 ) THEN
1888     WRITE(ilatex,*) ""
1889     WRITE(ilatex,'(A14)') "\end{document}"
1890     CLOSE(ilatex)
1891   END IF
1892 
1893   op%done     = .TRUE.
1894 !sui!write(std_out,*) "op%stats en of ctqmc_run",op%stats
1895 
1896 END SUBROUTINE Ctqmcoffdiag_run

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_setG0wTab [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_setG0wTab

FUNCTION

  Set Gow from input array

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  Gomega=G0w
  opt_fk=F is already inversed with out iwn

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

908 SUBROUTINE Ctqmcoffdiag_setG0wTab(op,Gomega,opt_fk)
909 
910 !Arguments ------------------------------------
911   TYPE(Ctqmcoffdiag), INTENT(INOUT)                      :: op
912   COMPLEX(KIND=8), DIMENSION(:,:,:), INTENT(IN ) :: Gomega
913   INTEGER                         , INTENT(IN ) :: opt_fk
914 !Local variable -------------------------------
915   DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: F
916 
917   IF ( .NOT. op%para ) &
918     CALL ERROR("Ctqmcoffdiag_setG0wTab : Ctqmcoffdiag_setParameters never called    ") 
919 
920   MALLOC(F,(1:op%samples+1,1:op%flavors,1:op%flavors))
921   CALL Ctqmcoffdiag_computeF(op,Gomega, F, opt_fk)  ! mu is changed
922  !write(6,*) "eee111"
923   CALL BathOperatoroffdiag_setF(op%Bath, F)
924  ! CALL BathOperatoroffdiag_printF(op%Bath,333)
925  !write(6,*) "eee"
926   FREE(F)
927 
928   op%inF = .TRUE.
929   op%set = .TRUE. 
930 
931 END SUBROUTINE Ctqmcoffdiag_setG0wTab

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_sethybri_limit [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_sethybri_limit

FUNCTION

  use coefficient A such that F=-A/(iwn) given by DMFT code.

COPYRIGHT

  Copyright (C) 2013-2024 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

  hybri_limit(nflavor,nflavor)=contains the limit for each couple of flavors

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

  op(Ctqmcoffdiag_type) = is the ctqmc main variable
           op&limit is now filled

NOTES

SOURCE

1211 SUBROUTINE Ctqmcoffdiag_sethybri_limit(op, hybri_limit)
1212 
1213 !Arguments ------------------------------------
1214   TYPE(Ctqmcoffdiag)                     , INTENT(INOUT) :: op
1215   COMPLEX(KIND=8) , DIMENSION(:,:),  INTENT(IN ) :: hybri_limit
1216 
1217   IF ( op%flavors .NE. SIZE(hybri_limit,1) ) &
1218     CALL ERROR("Error in sethybri_limit")
1219 
1220   op%hybri_limit(:,:)=hybri_limit(:,:)  
1221   op%opt_hybri_limit = 1
1222 END SUBROUTINE Ctqmcoffdiag_sethybri_limit

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_setMagmom [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_setMagmom

FUNCTION

  set the Magnetic moment matrix for susceptibility

COPYRIGHT

  Copyright (C) 2013-2024 ABINIT group (F. Gendron)
  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

SIDE EFFECTS

NOTES

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

4878 SUBROUTINE Ctqmcoffdiag_setMagmom(op,Magmom_orb, Magmom_spin, Magmom_tot)
4879 
4880 !Arguments ------------------------------------
4881   TYPE(Ctqmcoffdiag), INTENT(INOUT) ::op
4882   DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_orb
4883   DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_spin
4884   DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: Magmom_tot
4885 !Local variables ------------------------------
4886 !  INTEGER :: iflavor1,iflavor2
4887 
4888  ! do iflavor1=1,10
4889  !   do iflavor2=1,10
4890  !      if(iflavor1==iflavor2) THEN
4891  !        write(6,*) iflavor1, iflavor2, Magmom(iflavor1,iflavor2)
4892  !      end if
4893  !   end do
4894  ! end do
4895 
4896   IF ( SIZE(Magmom_orb) .NE. op%flavors*op%flavors ) &
4897     CALL ERROR("Ctqmcoffdiag_setU : Wrong Magnetic Moment matrix (size)        ")
4898 
4899   CALL ImpurityOperator_setMagmommat(op%Impurity, Magmom_orb, Magmom_spin, Magmom_tot)
4900 
4901 END SUBROUTINE Ctqmcoffdiag_setMagmom

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_setMu [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_setMu

FUNCTION

  impose energy levels

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  levels=energy levels vector

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

1170 SUBROUTINE Ctqmcoffdiag_setMu(op, levels)
1171 
1172 !Arguments ------------------------------------
1173   TYPE(Ctqmcoffdiag)                     , INTENT(INOUT) :: op
1174   DOUBLE PRECISION, DIMENSION(:), INTENT(IN   ) :: levels
1175 
1176   IF ( op%flavors .NE. SIZE(levels,1) ) &
1177     CALL WARNALL("Ctqmcoffdiag_setMu : Taking energy levels from weiss G(iw)")
1178 
1179   op%mu(:)=-levels(:)  ! levels = \epsilon_j - \mu
1180   !op%mu =\tilde{\mu} = \mu -\epsilon_j
1181   op%opt_levels = 1
1182 END SUBROUTINE Ctqmcoffdiag_setMu

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_setParameters [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_setParameters

FUNCTION

  set all parameters and operators

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  buffer=input parameters

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

498 SUBROUTINE Ctqmcoffdiag_setParameters(op,buffer)
499 
500 !Arguments ------------------------------------
501   TYPE(Ctqmcoffdiag), INTENT(INOUT)                         :: op
502   DOUBLE PRECISION, DIMENSION(1:11), INTENT(IN   ) :: buffer
503 
504 
505   op%thermalization = INT(buffer(3)) !op%thermalization
506   CALL Ctqmcoffdiag_setSeed(op,INT(buffer(1)))
507   CALL Ctqmcoffdiag_setSweeps(op,buffer(2))
508 
509   op%measurements   = INT(buffer(4)) !op%measurements
510   op%flavors        = INT(buffer(5))
511   op%samples        = INT(buffer(6)) !op%samples
512   op%beta           = buffer(7)      !op%beta
513   op%U              = buffer(8)      !U
514   op%opt_nondiag    = INT(buffer(10))
515   op%nspinor        = INT(buffer(11)) !op%nspinor
516 !  op%mu             = buffer(9)      !op%mu
517   !op%Wmax           = INT(buffer(9)) !Freq
518 !#ifdef CTCtqmcoffdiag_ANALYSIS
519 !  op%order          = INT(buffer(10)) ! order
520   op%inv_dt         = op%samples / op%beta
521 !#endif
522 
523   !CALL ImpurityOperator_init(op%Impurity,op%flavors,op%beta, op%samples)
524   CALL ImpurityOperator_init(op%Impurity,op%flavors,op%beta)
525   IF ( op%U .GE. 0.d0 ) THEN
526     CALL ImpurityOperator_computeU(op%Impurity,op%U,0.d0)
527     op%setU = .TRUE.
528   END IF
529 !  op%mu = op%mu + op%Impurity%shift_mu
530 !sui!write(std_out,*) "op%opt_nondiag",op%opt_nondiag
531   CALL BathOperatoroffdiag_init(op%Bath, op%flavors, op%samples, op%beta, INT(buffer(9)), op%opt_nondiag)
532 
533   op%para = .TRUE.
534 
535 END SUBROUTINE Ctqmcoffdiag_setParameters

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_setSeed [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_setSeed

FUNCTION

  initialize random number generator

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  iseed=seed from imput

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

618 SUBROUTINE Ctqmcoffdiag_setSeed(op,iseed)
619 
620 !Arguments ------------------------------------
621   TYPE(Ctqmcoffdiag), INTENT(INOUT)           :: op
622   INTEGER  , INTENT(IN   )           :: iseed
623 !Local variables ------------------------------
624   !INTEGER                            :: n
625   !INTEGER                            :: i
626   !INTEGER, DIMENSION(:), ALLOCATABLE :: seed
627 
628 
629   !CALL RANDOM_SEED(size = n)
630   !MALLOC(seed,(n))
631   !seed =  iseed + (/ (i - 1, i = 1, n) /)
632 
633   !CALL RANDOM_SEED(PUT = seed+op%rank)
634 
635   !FREE(seed)
636 
637   op%seed=INT(iseed+op%rank,8)
638 
639 END SUBROUTINE Ctqmcoffdiag_setSeed

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_setSweeps [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_setSweeps

FUNCTION

  set the number of sweeps

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  sweeps=asked sweeps

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

563 SUBROUTINE Ctqmcoffdiag_setSweeps(op,sweeps)
564 
565 !Arguments ------------------------------------
566   TYPE(Ctqmcoffdiag)         , INTENT(INOUT) :: op
567   DOUBLE PRECISION  , INTENT(IN   ) :: sweeps
568 
569   op%sweeps = NINT(sweeps / DBLE(op%size))
570 !  !write(std_out,*) op%sweeps,NINT(sweeps / DBLE(op%size)),ANINT(sweeps/DBLE(op%size))
571   IF ( DBLE(op%sweeps) .NE. ANINT(sweeps/DBLE(op%size)) ) &
572     CALL ERROR("Ctqmcoffdiag_setSweeps : sweeps is negative or too big     ")
573   IF ( op%sweeps .LT. 2*CTQMC_SLICE1 ) THEN  !202
574     CALL WARNALL("Ctqmcoffdiag_setSweeps : # sweeps automtically changed     ")
575     op%sweeps = 2*CTQMC_SLICE1
576 !  ELSE IF ( op%sweeps .LT. op%thermalization ) THEN
577 !    CALL WARNALL("Ctqmcoffdiag_setSweeps : Thermalization > sweeps / cpu -> auto fix")
578 !    op%sweeps = op%thermalization
579   END IF
580   IF ( DBLE(NINT(DBLE(op%sweeps)*DBLE(op%size)/DBLE(CTQMC_SLICE1))) .NE.  &
581   ANINT(DBLE(op%sweeps)*DBLE(op%size)/DBLE(CTQMC_SLICE1)) ) THEN
582     op%modNoise1 = op%sweeps
583   ELSE
584     op%modNoise1    = MIN(op%sweeps,INT(DBLE(op%sweeps)*DBLE(op%size) / DBLE(CTQMC_SLICE1))) !101
585   END IF
586   op%modNoise2    = MAX(op%modNoise1 / CTQMC_SLICE2, 1)   ! 100
587 !  op%modGlobalMove(1) = op%thermalization / 10 + 1
588 !  op%modGlobalMove(2) = 0
589 
590 END SUBROUTINE Ctqmcoffdiag_setSweeps

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_setU [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_setU

FUNCTION

  set the interaction matrix

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  matU=interaction matrix

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

 990 SUBROUTINE Ctqmcoffdiag_setU(op,matU)
 991 
 992 !Arguments ------------------------------------
 993   TYPE(Ctqmcoffdiag), INTENT(INOUT) ::op
 994 !Local variables ------------------------------
 995   DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: matU
 996 
 997   IF ( SIZE(matU) .NE. op%flavors*op%flavors ) &
 998     CALL ERROR("Ctqmcoffdiag_setU : Wrong interaction matrix (size)        ")
 999 
1000   CALL ImpurityOperator_setUmat(op%Impurity, matU)
1001   op%setU = .TRUE.
1002 END SUBROUTINE Ctqmcoffdiag_setU

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_symmetrizeGreen [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_symmetrizeGreen

FUNCTION

  optionnaly symmetrize the green functions

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc
  syms=weight factors

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

SOURCE

3702 SUBROUTINE Ctqmcoffdiag_symmetrizeGreen(op, syms)
3703 
3704 !Arguments ------------------------------------
3705   TYPE(Ctqmcoffdiag)                     , INTENT(INOUT) :: op
3706   DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN   ) :: syms
3707 !Local variables ------------------------------
3708   !INTEGER :: iflavor1
3709   !INTEGER :: iflavor2
3710   !INTEGER :: flavors
3711   !DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: green_tmp
3712   !DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:  ) :: n_tmp
3713 
3714   ABI_UNUSED((/syms(1,1), op%swap/))
3715 
3716 !  flavors = op%flavors
3717 !  IF ( SIZE(syms,1) .NE. flavors .OR. SIZE(syms,2) .NE. flavors ) THEN
3718 !    CALL WARNALL("Ctqmcoffdiag_symmetrizeGreen : wrong opt_sym -> not symmetrizing")
3719 !    RETURN
3720 !  END IF
3721 ! 
3722 !  MALLOC(green_tmp,(1:op%samples+1,flavors))
3723 !  green_tmp(:,:) = 0.d0
3724 !  MALLOC(n_tmp,(1:flavors))
3725 !  n_tmp(:) = 0.d0
3726 !  DO iflavor1=1, flavors
3727 !    DO iflavor2=1,flavors
3728 !      green_tmp(:,iflavor1) = green_tmp(:,iflavor1) &
3729 !                             + syms(iflavor2,iflavor1) * op%Greens(iflavor2)%oper(:)
3730 !      n_tmp(iflavor1) = n_tmp(iflavor1) &
3731 !                             + syms(iflavor2,iflavor1) * op%measN(1,iflavor2)
3732 !    END DO
3733 !  END DO
3734 !  DO iflavor1=1, flavors
3735 !    op%Greens(iflavor1)%oper(:) = green_tmp(:,iflavor1)
3736 !    op%measN(1,iflavor1)          = n_tmp(iflavor1)
3737 !  END DO
3738 !  FREE(green_tmp)
3739 !  FREE(n_tmp)
3740 END SUBROUTINE Ctqmcoffdiag_symmetrizeGreen

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_tryAddRemove [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_tryAddRemove

FUNCTION

  Try to add or remove a segment and an anti-segment

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

  updated=something changed

SIDE EFFECTS

NOTES

SOURCE

2185 SUBROUTINE Ctqmcoffdiag_tryAddRemove(op,updated)
2186 
2187 !Arguments ------------------------------------
2188   TYPE(Ctqmcoffdiag)             , INTENT(INOUT) :: op
2189 !  TYPE(BathOperatoroffdiag)    , INTENT(INOUT) :: Bath 
2190 !  TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity 
2191   LOGICAL               , INTENT(  OUT) :: updated
2192 !Local variables ------------------------------
2193   INTEGER                               :: position
2194   INTEGER         , DIMENSION(1:2)     :: nature ! -2 for antiseg and 1 for seg
2195   INTEGER                               :: i! -2 for antiseg and 1 for seg
2196   !INTEGER                               :: it,it1 !ii,
2197   DOUBLE PRECISION                      :: action
2198   DOUBLE PRECISION                      :: beta
2199   DOUBLE PRECISION                      :: time1
2200   DOUBLE PRECISION                      :: time2
2201   DOUBLE PRECISION                      :: time_avail
2202   DOUBLE PRECISION                      :: det_ratio,sign_det_ratio
2203   DOUBLE PRECISION                      :: overlap
2204   DOUBLE PRECISION                      :: length
2205   DOUBLE PRECISION                      :: signe
2206   DOUBLE PRECISION                      :: tail
2207   INTEGER                      :: tailint
2208   DOUBLE PRECISION                      :: signdet, signdetprev
2209   DOUBLE PRECISION, DIMENSION(1:2)      :: CdagC_1
2210 
2211   IF ( .NOT. op%set ) &
2212     CALL ERROR("Ctqmcoffdiag_trySegment : QMC not set                       ")
2213 
2214         !write(std_out,*) "      TryAddRemove start"
2215   nature(1) = CTQMC_SEGME
2216   nature(2) = CTQMC_ANTIS
2217   beta      = op%beta
2218 
2219   updated = .FALSE.
2220   tailint  = (op%Impurity%particles(op%Impurity%activeFlavor)%tail)
2221   tail  = DBLE(tailint)
2222   !write(std_out,*) "op%Impurity%particles(op%Impurity%activeFlavor)%tail",op%Impurity%activeFlavor,tail
2223 
2224 
2225   !=====================================
2226   ! First choose segment or antisegment
2227   !=====================================
2228   DO i = 1, 2
2229     signe = SIGN(1.d0,DBLE(nature(i))) 
2230 !      -----  1: segment        signe= 1  ( CTQMC_SEGME =  1 )
2231 !      -----  2: antisegment    signe=-1  ( CTQMC_ANTIS = -2 )
2232 !    NB: Sign(a,b) = sign(b) * a
2233 
2234  !prt!if(op%prtopt==1) write(std_out,*) "       ==Starting configuration",i
2235  !prt!if(op%prtopt==1) write(std_out,*) "        = Segments:"
2236     tailint  = (op%Impurity%particles(op%Impurity%activeFlavor)%tail)
2237 !prt!    do ii=0, op%Impurity%Particles(op%Impurity%activeFlavor)%tail
2238  !prt!if(op%prtopt==1)  write(std_out,*) ii, op%Impurity%Particles(op%Impurity%activeFlavor)%list(ii,1), &
2239 !prt!&                    op%Impurity%Particles(op%Impurity%activeFlavor)%list(ii,2)
2240 !prt!    enddo
2241   !sui!write(std_out,*) "        = M Matrix",op%Bath%sumtails
2242 !prt!    do it=1,op%Bath%sumtails
2243     !sui!write(std_out,'(a,3x,500e10.3)') "        M start",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2244 !prt!    enddo
2245     CALL OurRng(op%seed,action)
2246 
2247     !==========================
2248     ! Add segment/antisegment
2249     !==========================
2250     IF ( action .LT. .5d0 ) THEN ! Add a segment or antisegment
2251      !ii  write(std_out,*) "        =try: Segment added of type",i,op%prtopt
2252 
2253       ! Select time1 (>0) in [0,beta]
2254       !==============================
2255       CALL OurRng(op%seed,time1)
2256       time1 = time1 * beta
2257 
2258       ! time_avail is the distance between between time1 and 
2259       !   - the next start of a segment for a segment addition
2260       !   - the next end of a segment for an antisegment addition
2261       ! ImpurityOperator_getAvailableTime > 0 for a segment      (signe>0) -> time_avail>0
2262       ! ImpurityOperator_getAvailableTime < 0 for an antisegment (signe<0) -> time_avail>0
2263       !====================================================================
2264       time_avail = ImpurityOperator_getAvailableTime(op%Impurity,time1,position) * signe
2265      !ii  write(std_out,*) "        =try: time_avail",time_avail,time1
2266       IF ( time_avail .GT. 0.d0 ) THEN
2267 
2268         ! Time2 is  the length of the proposed new (anti)segment
2269         !=======================================================
2270         CALL OurRng(op%seed,time2)
2271         IF ( time2 .EQ. 0.d0 ) CALL OurRng(op%seed,time2) ! Prevent null segment
2272 
2273         ! Now time2 is the time at the end of the proposed new (anti) segment
2274         ! time2 > time1 
2275         !====================================================================
2276         time2     = time1 + time2 * time_avail
2277       !sui!write(std_out,*) tailint+1,time1,time2,position
2278 !        CALL CdagC_init(CdagC_1,time1,time2)
2279 
2280         ! CdagC_1 gives the stard/end times for the proposed new segment/antisegment
2281         ! CdagC1(C_) can  be above beta.
2282         !  For a      segment CdagC_1(Cdag_) = time1 < CdagC_1(C_) = time2, l=time2-time1 > 0
2283         !  For a anti segment CdagC_1(Cdag_) = time2 > CdagC_1(C_) = time1, l=time1-time2 < 0
2284         !  time2 can be above beta and thus for a      segment CdagC_1(C_   ) > beta
2285         !  time2 can be above beta and thus for an antisegment CdagC_1(Cdag_) > beta
2286         !  length > 0 for     segment
2287         !  length < 0 for antisegment
2288         !====================================================================================
2289         CdagC_1(Cdag_) = ((1.d0+signe)*time1+(1.d0-signe)*time2)*0.5d0
2290         CdagC_1(C_   ) = ((1.d0+signe)*time2+(1.d0-signe)*time1)*0.5d0
2291 !        length    = CdagC_length(CdagC_1)
2292         length    = CdagC_1(C_   ) - CdagC_1(Cdag_)
2293         !write(std_out,*) "      try : times", CdagC_1(C_   ),CdagC_1(Cdag_)
2294         !write(std_out,*) "      length", length
2295 
2296 !      -----  Computes the determinant ratio
2297         det_ratio = BathOperatoroffdiag_getDetAdd(op%Bath,CdagC_1,position,op%Impurity%particles) 
2298 
2299 !      -----  Computes the overlap
2300         overlap   = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
2301         signdetprev  = ImpurityOperator_getsign(op%Impurity, time2, i, action, position)
2302 
2303         !write(std_out,*) "      overlap   ", overlap
2304         CALL OurRng(op%seed,time1)
2305         !write(std_out,*) "      Rnd", time1
2306         signdet=1.d0
2307         det_ratio=det_ratio*signdetprev
2308                  
2309         IF ( det_ratio .LT. 0.d0 ) THEN
2310         !sui!write(std_out,*) "         NEGATIVE DET",det_ratio,signdetprev
2311           det_ratio   = - det_ratio
2312           sign_det_ratio=-1
2313           op%stats(nature(i)+CTQMC_DETSI) = op%stats(nature(i)+CTQMC_DETSI) + 1.d0
2314        !  op%signvaluecurrent=-1.d0
2315         ELSE
2316           sign_det_ratio=1
2317        !  op%signvaluecurrent=+1.d0
2318          ! signdet=-1.d0
2319         !sui!write(std_out,*) "                  DET",det_ratio,signdetprev
2320         END IF
2321       !ii  write(std_out,*) "                  DET",det_ratio
2322        ! op%signvaluemeas=op%signvaluemeas+1.d0
2323         !write(std_out,*) "        .................",(time1 * (tail + 1.d0 )),beta * time_avail * det_ratio * DEXP(op%mu(op%Impurity%activeFlavor)*length + overlap)
2324         !write(std_out,*) "        .................",beta , time_avail , op%mu(op%Impurity%activeFlavor),op%Impurity%activeFlavor
2325 
2326         IF ( (time1 * (tail + 1.d0 )) &
2327              .LT. (beta * time_avail * det_ratio * DEXP(op%mu(op%Impurity%activeFlavor)*length + overlap) ) ) THEN
2328 !          write(*,*) "before"
2329 !          CALL ListCdagCoffdiag_print(op%Impurity%particles(op%Impurity%activeFlavor),6)
2330           CALL ImpurityOperator_add(op%Impurity,CdagC_1,position)
2331 !          write(*,*) "after "
2332 !          CALL ListCdagCoffdiag_print(op%Impurity%particles(op%Impurity%activeFlavor),6)
2333           CALL BathOperatoroffdiag_setMAdd(op%bath,op%Impurity%particles) 
2334           op%stats(nature(i)+CTQMC_ADDED) = op%stats(nature(i)+CTQMC_ADDED)  + 1.d0
2335           updated = .TRUE. .OR. updated
2336           tail = tail + 1.d0
2337           tailint = tailint + 1
2338 !          read(*,*) time1
2339           !ii  write(6,*) "        Accepted addition, new conf is",time1
2340          !prt! do ii=0, op%Impurity%Particles(op%Impurity%activeFlavor)%tail
2341           !prt!if(op%prtopt==1)  write(6,*) ii, op%Impurity%Particles(op%Impurity%activeFlavor)%list(ii,1),&
2342 !prt!&                          op%Impurity%Particles(op%Impurity%activeFlavor)%list(ii,2)
2343          !prt! enddo
2344         !sui!write(6,*) "        = M Matrix"
2345          !prt! do it=1,op%Bath%sumtails
2346           !sui!write(6,*) "        M new",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2347          !prt! enddo
2348 
2349           IF ( sign_det_ratio .LT. 0.d0 ) op%signvalue=-op%signvalue
2350         !sui!write(6,*) "                  signvalue",op%signvalue
2351         ELSE
2352      !ii  write(6,*) "        Refused      addition: proba",time1
2353         END IF 
2354       ELSE
2355     !sui!write(6,*) "        Refused      addition: time_avail <0"
2356       END IF 
2357 
2358     !========================================
2359     ! Remove segment/antisegment
2360     !========================================
2361     ELSE ! Remove a segment among the segment of the flavor activeflavor
2362       !ii if(op%prtopt==1)  write(6,*) "        =try: Segment removed of type",i
2363       IF ( tail .GT. 0.d0 ) THEN
2364         CALL OurRng(op%seed,time1)
2365         position = INT(((time1 * tail) + 1.d0) * signe )
2366         !prt!if(op%prtopt==1)  write(6,*) "         position",position 
2367         time_avail = ImpurityOperator_getAvailedTime(op%Impurity,position)
2368         det_ratio  = BathOperatoroffdiag_getDetRemove(op%Bath,position)
2369         !write(6,*) "        det_ratio", det_ratio
2370         CdagC_1    = ImpurityOperator_getSegment(op%Impurity,position)
2371 !        length     = CdagC_length(CdagC_1)
2372         length     = CdagC_1(C_) - CdagC_1(Cdag_)
2373         !write(6,*) "        length   ", length
2374         overlap    = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
2375         !write(6,*) "        overlap  ", overlap
2376         CALL OurRng(op%seed,time1)
2377         !write(6,*) "        Random   ",time1
2378         signdetprev = ImpurityOperator_getsign(op%Impurity, time2, i, action, position)
2379         det_ratio=det_ratio*signdetprev
2380         signdet=1.d0
2381         IF ( det_ratio .LT. 0.d0 ) THEN
2382         !sui!write(6,*) "         NEGATIVE DET",det_ratio,signdetprev
2383           det_ratio   = -det_ratio
2384           sign_det_ratio=-1
2385 !          op%seg_sign = op%seg_sign + 1.d0
2386           op%stats(nature(i)+CTQMC_DETSI) = op%stats(nature(i)+CTQMC_DETSI) + 1.d0
2387           signdet=-1.d0
2388         ELSE 
2389           sign_det_ratio=1
2390         !sui!write(6,*) "                  DET",det_ratio,signdetprev
2391         END IF
2392        !ii  write(6,*) "                  DET",det_ratio
2393         IF ( (time1 * beta * time_avail * DEXP(op%mu(op%Impurity%activeFlavor)*length+overlap)) &
2394              .LT. (tail * det_ratio ) ) THEN
2395           CALL ImpurityOperator_remove(op%Impurity,position)
2396           CALL BathOperatoroffdiag_setMRemove(op%Bath,op%Impurity%particles) 
2397           !op%seg_removed = op%seg_removed  + 1.d0
2398           op%stats(nature(i)+CTQMC_REMOV) = op%stats(nature(i)+CTQMC_REMOV)  + 1.d0
2399           updated = .TRUE. .OR. updated
2400           tail = tail -1.d0
2401           tailint = tailint -1
2402           !ii  write(6,*) "        Accepted removal, new conf is:",time1
2403       !prt!    do ii=0, op%Impurity%Particles(op%Impurity%activeFlavor)%tail
2404           !prt!if(op%prtopt==1)  write(6,*) ii, op%Impurity%Particles(op%Impurity%activeFlavor)%list(ii,1),&
2405 !prt!&                          op%Impurity%Particles(op%Impurity%activeFlavor)%list(ii,2)
2406       !prt!    enddo
2407         !sui!write(6,*) "        = M Matrix"
2408        !prt!   do it=1,op%Bath%sumtails
2409           !sui!write(6,*) "        M new",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2410         !prt!  enddo
2411           IF ( sign_det_ratio .LT. 0.d0 ) op%signvalue=-op%signvalue
2412         !sui!write(6,*) "                  signvalue",op%signvalue
2413         ELSE
2414      !ii  write(6,*) "        Refused      removal",time1
2415         END IF
2416       ELSE
2417       !sui!write(6,*) "        Refused      removal: no segment available"
2418       END IF
2419     END IF
2420     !========================================
2421     ! End Add/Remove Antisegment
2422     !========================================
2423   END DO
2424 END SUBROUTINE Ctqmcoffdiag_tryAddRemove

ABINIT/m_Ctqmcoffdiag/Ctqmcoffdiag_trySwap [ Functions ]

[ Top ] [ Functions ]

NAME

  Ctqmcoffdiag_trySwap

FUNCTION

  try a global move (swap to flavors)

COPYRIGHT

  Copyright (C) 2013-2024 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

  op=ctqmc

OUTPUT

  flav_i=first flavor swaped
  flav_j=second flavor swaped

SIDE EFFECTS

NOTES

SOURCE

2624 SUBROUTINE Ctqmcoffdiag_trySwap(op,flav_i,flav_j)
2625 
2626 !Arguments ------------------------------------
2627   TYPE(Ctqmcoffdiag)           , INTENT(INOUT) :: op
2628 !  TYPE(BathOperatoroffdiag)    , INTENT(INOUT) :: Bath 
2629 !  TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity 
2630   INTEGER               , INTENT(  OUT) :: flav_i
2631   INTEGER               , INTENT(  OUT) :: flav_j
2632 !Local variables ------------------------------
2633   INTEGER :: flavor_i
2634   INTEGER :: flavor_j !,ii,it,it1 !,iflavor
2635   DOUBLE PRECISION :: rnd
2636   DOUBLE PRECISION :: lengthi
2637   DOUBLE PRECISION :: lengthj
2638   DOUBLE PRECISION :: overlapic1
2639   DOUBLE PRECISION :: overlapjc1
2640   DOUBLE PRECISION :: overlapic2
2641   DOUBLE PRECISION :: overlapjc2
2642   !DOUBLE PRECISION :: detic1
2643   !DOUBLE PRECISION :: detjc1
2644   !DOUBLE PRECISION :: detic2
2645   !DOUBLE PRECISION :: detjc2
2646   DOUBLE PRECISION :: det_ratio,detnew,detold
2647   DOUBLE PRECISION :: local_ratio
2648  ! TYPE(BathOperatoroffdiag)  :: Bathnew
2649 
2650 
2651   !CALL RANDOM_NUMBER(rnd)
2652   CALL OurRng(op%seed,rnd)
2653   flavor_i = NINT(rnd*DBLE(op%flavors-1.d0))+1
2654   !CALL RANDOM_NUMBER(rnd)
2655   CALL OurRng(op%seed,rnd)
2656   flavor_j = NINT(rnd*DBLE(op%flavors-1.d0))+1
2657   !ii write(6,'(a,2i4)') "--------------- new swap --------------------------------",flavor_i,flavor_j
2658   
2659   flav_i = 0
2660   flav_j = 0
2661   !ii   do iflavor=1,op%flavors
2662   !ii     write(6,*) "BEFORE  GMOVE For flavor", iflavor,"size is",op%Impurity%particles(iflavor)%tail," and  Conf is :"
2663   !ii     do ii=1, op%Impurity%Particles(iflavor)%tail
2664   !ii       write(6,'(i4,100f12.3)') ii, op%Impurity%Particles(iflavor)%list(ii,1),&
2665   !ii  &                   op%Impurity%Particles(iflavor)%list(ii,2)
2666   !ii     enddo
2667   !ii   enddo
2668   !ii   write(6,*) "        = M Matrix"
2669   !ii   write(6,'(a,2x,100(i12))') "Flavor=",((iflavor,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2670   !ii   write(6,'(i21,100i12)') ((it,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2671   !ii   do it=1,op%Bath%sumtails
2672   !ii     write(6,'(a,100f12.3)') "        M before",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2673   !ii   enddo
2674 
2675   ! todoba this part
2676   IF ( flavor_i .NE. flavor_j ) THEN
2677     !CALL BathOperatoroffdiag_init(Bathnew, op%flavors, op%samples, op%beta, 0, op%opt_nondiag)
2678     ! On tente d'intervertir i et j
2679     ! Configuration actuelle :
2680 
2681     op%modGlobalMove(2) = op%modGlobalMove(2)+1
2682     ! ===========================================
2683     ! First use M matrix to compute determinant
2684     ! ===========================================
2685     detold     = BathOperatoroffdiag_getDetF(op%Bath) ! use op%Bath%M
2686 
2687     ! ===========================================
2688     ! Second build M_update matrix to compute determinant after.
2689     ! ===========================================
2690     !CALL ListCdagCoffdiag_print(particle)
2691     call BathOperatoroffdiag_recomputeM(op%Bath,op%impurity%particles,flavor_i,flavor_j) ! compute op%Bath%M_update
2692     detnew     = BathOperatoroffdiag_getDetF(op%Bath,option=1) ! use op%Bath%M_update
2693 
2694     lengthi    = ImpurityOperator_measN(op%Impurity,flavor_i)
2695     lengthj    = ImpurityOperator_measN(op%Impurity,flavor_j)
2696     overlapic1 = ImpurityOperator_overlapFlavor(op%Impurity,flavor_i)
2697     overlapjc1 = ImpurityOperator_overlapFlavor(op%Impurity,flavor_j)
2698     ! lengths unchanged
2699     overlapic2 = ImpurityOperator_overlapSwap(op%Impurity,flavor_i,flavor_j)
2700     overlapjc2 = ImpurityOperator_overlapSwap(op%Impurity,flavor_j,flavor_i)
2701 
2702 !    IF ( detic1*detjc1 .EQ. detic2*detjc2 ) THEN
2703 !      det_ratio = 1.d0
2704 !    ELSE IF ( detic1*detjc1 .EQ. 0.d0 ) THEN
2705 !      det_ratio = detic2*detjc2 ! evite de diviser par 0 si pas de segment
2706 !    ELSE
2707 
2708     det_ratio = detnew/detold ! because the determinant is the determinant of F
2709    !ii  write(6,*) "det_ratio, detold,detnew",det_ratio, detold,detnew, detold/detnew
2710 
2711 !    END IF
2712     local_ratio = DEXP(-overlapic2*overlapjc2+overlapic1*overlapjc1 &
2713                       +(lengthj-lengthi)*(op%mu(flavor_i)-op%mu(flavor_j)))
2714    !ii  write(6,*) "local_ratio",local_ratio
2715 
2716     ! Wloc = exp(muN-Uo)
2717     !CALL RANDOM_NUMBER(rnd)
2718     CALL OurRng(op%seed,rnd)
2719     IF ( rnd .LT. local_ratio*det_ratio ) THEN ! swap accepted
2720    !ii    write(6,*) "        = M Matrix before swap"
2721    !ii    write(6,'(a,2x,100(i12))') "Flavor=",((iflavor,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2722    !ii    write(6,'(i21,100i12)') ((it,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2723    !ii    do it=1,op%Bath%sumtails
2724    !ii      write(6,'(a,100f12.3)') "        M after ",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2725    !ii    enddo
2726    !ii    do it=1,op%Bath%sumtails
2727    !ii      write(6,'(a,100f12.3)') " update M after ",(op%Bath%M_update%mat(it,it1),it1=1,op%Bath%sumtails)
2728    !ii    enddo
2729    !ii    write(6,*) "Gmove accepted",rnd,local_ratio*det_ratio
2730       CALL ImpurityOperator_swap(op%Impurity, flavor_i,flavor_j)
2731       CALL BathOperatoroffdiag_swap    (op%Bath    , flavor_i,flavor_j) !  use op%Bath%M_update to built new op%Bath%M
2732       
2733       op%swap = op%swap + 1.d0
2734       flav_i = flavor_i
2735       flav_j = flavor_j
2736     ELSE
2737    !ii   write(6,*) "Gmove refused",rnd,local_ratio*det_ratio
2738 !      CALL WARN("Swap refused")
2739 !      WRITE(op%ostream,'(6E24.14)') local_ratio, det_ratio, detic1, detjc1, detic2, detjc2
2740     END IF
2741    ! CALL BathOperatoroffdiag_destroy(Bathnew)
2742   END IF
2743  !ii  do iflavor=1,op%flavors
2744  !ii    write(6,*) "AFTER   GMOVE For flavor", iflavor,"size is",op%Impurity%particles(iflavor)%tail," and  Conf is :"
2745  !ii    do ii=1, op%Impurity%Particles(iflavor)%tail
2746  !ii      write(6,'(15x,i4,100f12.3)') ii, op%Impurity%Particles(iflavor)%list(ii,1),&
2747  !ii &                   op%Impurity%Particles(iflavor)%list(ii,2)
2748  !ii    enddo
2749  !ii  enddo
2750  !ii  write(6,*) "        = M Matrix"
2751  !ii  write(6,'(a,2x,100(i12))') "Flavor=",((iflavor,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2752  !ii  write(6,'(i21,100i12)') ((it,it=1,op%Impurity%particles(iflavor)%tail),iflavor=1,op%flavors)
2753  !ii  do it=1,op%Bath%sumtails
2754  !ii    write(6,'(a,100f12.3)') "        M after ",(op%Bath%M%mat(it,it1),it1=1,op%Bath%sumtails)
2755  !ii  enddo
2756  !ii  do it=1,op%Bath%sumtails
2757  !ii    write(6,'(a,100f12.3)') " update M after ",(op%Bath%M_update%mat(it,it1),it1=1,op%Bath%sumtails)
2758  !ii  enddo
2759 
2760 END SUBROUTINE Ctqmcoffdiag_trySwap

m_Ctqmcoffdiag/Ctqmcoffdiag [ Types ]

[ Top ] [ m_Ctqmcoffdiag ] [ Types ]

NAME

  Ctqmcoffdiag

FUNCTION

  This structured datatype contains the necessary data

COPYRIGHT

  Copyright (C) 2013-2024 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

100 TYPE Ctqmcoffdiag
101 
102   LOGICAL :: init = .FALSE.
103 ! Flag: is MC initialized
104 
105   LOGICAL :: set  = .FALSE.
106 ! Flag: ??
107 
108   LOGICAL :: setU = .FALSE.
109 ! Flag: is U Set ?
110 
111   LOGICAL :: inF  = .FALSE.
112 ! Flag: is hybridization fct in input ?
113 
114   LOGICAL :: done = .FALSE.
115 ! Flag: is MC terminated ?
116 
117   LOGICAL :: para = .FALSE.
118 ! Flag:  do we have parameters in input
119 
120   LOGICAL :: have_MPI = .FALSE.
121 ! Flag: 
122 
123   INTEGER :: opt_movie = 0
124 !
125 
126   INTEGER :: opt_analysis = 0
127 ! correlations 
128 
129   INTEGER :: opt_check = 0
130 ! various check 0
131 ! various check 1 impurity
132 ! various check 2 bath
133 ! various check 3 both
134 
135   INTEGER :: opt_order = 0
136 ! nb of segments max for analysis
137 
138   INTEGER :: opt_histo = 0
139 ! Enable histo calc.
140 
141   INTEGER :: opt_noise = 0
142 ! compute noise
143 
144   INTEGER :: opt_spectra = 0
145 ! markov chain FT (correlation time)
146 
147   INTEGER :: opt_levels = 0
148 ! do we have energy levels
149 
150   INTEGER :: opt_hybri_limit = 0
151 ! do we have limit of hybridization (yes=1)
152 
153   INTEGER :: opt_nondiag = 0
154 ! if opt_nondiag = 1 F is non diagonal.
155 
156   INTEGER :: prtopt = 1
157 ! printing
158 
159   INTEGER :: flavors
160 ! number of flavors 
161 
162   INTEGER :: nspinor
163 ! number of spinor
164 
165   INTEGER :: measurements
166 !  The modulo used to measure the interaction energy and the number of electrons. Example : 2 means the measure is perform every two sweeps. 
167 
168   INTEGER :: samples
169 ! nb of L points (dmftqmc_l)
170 
171   INTEGER(8) :: seed
172 !
173 
174   INTEGER :: sweeps
175 !
176 
177   INTEGER :: thermalization
178 !
179 
180   INTEGER :: ostream
181 ! output file
182 
183   INTEGER :: istream
184 ! input file
185 
186   INTEGER :: modNoise1
187 ! measure the noise each modNoise1
188 
189   INTEGER :: modNoise2
190 ! measure the noise each modNoise2
191 
192   INTEGER :: activeFlavor
193 ! orbital on which one do sth now
194 
195   INTEGER, DIMENSION(1:2) :: modGlobalMove
196 ! 1: global move each modglobalmove(1)
197 ! 2: we have done modglobalmove(2) for two different orbitals.
198 
199   INTEGER :: Wmax
200 ! Max freq for FT
201 
202   DOUBLE PRECISION, DIMENSION(1:6) :: stats
203 ! to now how many negative determinant, antisegments,seeme.e.twfs...j
204 
205   DOUBLE PRECISION :: swap
206 ! nb of successfull GM
207 
208   DOUBLE PRECISION :: signvalue
209 
210   INTEGER :: MY_COMM
211 ! 
212 
213   INTEGER :: rank
214 !
215 
216   INTEGER :: size
217 ! size of MY_COMM
218 
219   DOUBLE PRECISION :: runTime ! time for the run routine
220 !  
221 
222   DOUBLE PRECISION :: beta
223 !
224 
225   DOUBLE PRECISION :: U
226 
227   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: mu
228 ! levels
229 
230   COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: hybri_limit
231 ! coeff A such that F=-A/(iwn)
232 
233   TYPE(GreenHyboffdiag)                        :: Greens 
234 ! Green's function
235 
236   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: measN 
237 ! measure of occupations (3or4,flavor) 
238 
239   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:  ) :: measDE
240 !  (flavor,flavor) double occupancies
241 !  (1,1): total energy of correlation.
242 
243   DOUBLE PRECISION :: a_Noise
244 ! Noise a exp (-bx) for the  noise
245 
246   DOUBLE PRECISION :: b_Noise
247 ! Noise a exp (-bx) for the  noise
248 
249   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: abNoiseG   !(ab,tau,flavor)
250 ! Noise but for G
251 
252   TYPE(Vector)             , DIMENSION(1:2) :: measNoise 
253   TYPE(Vector), ALLOCATABLE, DIMENSION(:,:,:) :: measNoiseG       !(tau,flavor,mod) 
254 ! accumulate each value relataed to measurenoise 1 2
255 
256 !#ifdef CTCtqmcoffdiag_ANALYSIS
257 !  INTEGER                                     :: order
258   DOUBLE PRECISION                            :: inv_dt
259 ! 1/(beta/L)
260 
261   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:  ) :: measPerturbation 
262 ! opt_order,nflavor
263 
264   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: occup_histo_time
265 ! nflavor
266 
267   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: occupconfig
268 ! 2**nflavor
269 
270   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: suscep
271 ! samples
272 
273   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: chi
274 ! samples
275 
276   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: chicharge
277 ! samples
278 
279   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ntot
280 ! occupation total, t2g, eg
281 
282   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:  ) :: meas_fullemptylines
283 ! opt_order,nflavor
284 
285   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: measCorrelation 
286 ! segment,antisegment,nflavor,nflavor
287 
288 !#endif
289 !#ifdef CTCtqmcoffdiag_CHECK
290   DOUBLE PRECISION :: errorImpurity
291 ! check 
292 
293   DOUBLE PRECISION :: errorBath
294 ! for check
295 
296 !#endif
297   TYPE(BathOperatoroffdiag)              :: Bath
298 
299 
300   TYPE(ImpurityOperator)          :: Impurity
301 
302   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: density
303 
304 END TYPE Ctqmcoffdiag