1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! Copyright 2010. Los Alamos National Security, LLC. This material was ! 3! produced under U.S. Government contract DE-AC52-06NA25396 for Los Alamos ! 4! National Laboratory (LANL), which is operated by Los Alamos National ! 5! Security, LLC for the U.S. Department of Energy. The U.S. Government has ! 6! rights to use, reproduce, and distribute this software. NEITHER THE ! 7! GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, ! 8! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS ! 9! SOFTWARE. If software is modified to produce derivative works, such ! 10! modified software should be clearly marked, so as not to confuse it ! 11! with the version available from LANL. ! 12! ! 13! Additionally, this program is free software; you can redistribute it ! 14! and/or modify it under the terms of the GNU General Public License as ! 15! published by the Free Software Foundation; version 2.0 of the License. ! 16! Accordingly, this program is distributed in the hope that it will be ! 17! useful, but WITHOUT ANY WARRANTY; without even the implied warranty of ! 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General ! 19! Public License for more details. ! 20!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 21 22SUBROUTINE FERMIEXPANS 23 24 USE CONSTANTS_MOD 25 USE SETUPARRAY 26 USE FERMICOMMON 27 USE SPINARRAY 28 USE MYPRECISION 29 30 IMPLICIT NONE 31 32 INTEGER :: I, J, ITER 33 INTEGER :: BREAKLOOP 34 REAL(LATTEPREC) :: OCC, OCCERROR 35 REAL(LATTEPREC) :: PREVERROR, PREVERROR2, PREVERROR3 36 REAL(LATTEPREC) :: BOVER2M, TRX, TRXOMX 37 REAL(LATTEPREC) :: SHIFTCP, HFACT 38 REAL(LATTEPREC), PARAMETER :: MAXSHIFT = ONE 39 IF (EXISTERROR) RETURN 40 41 OCC = BNDFIL*FLOAT(HDIM) 42 43 ITER = 0 44 45 BREAKLOOP = 0 46 47 BOVER2M = ONE/(KBT*(TWO**(2 + FERMIM))) 48 49 SCFS = 0 50 51 PREVERROR = ZERO 52 PREVERROR2 = ZERO 53 PREVERROR3 = ZERO 54 55 DO WHILE (BREAKLOOP .EQ. 0) 56 57 SCFS = SCFS + 1 58 59 ITER = ITER + 1 60 61 IF (ITER .EQ. 100) THEN 62 CALL PANIC 63 CALL ERRORS("fermiexpans","Fermi expansion is not converging: STOP!") 64 ENDIF 65 66 IF (SPINON .EQ. 0) THEN 67 68 BO = -BOVER2M*H 69 70 HFACT = HALF + BOVER2M*CHEMPOT 71 72 DO I = 1, HDIM 73 BO(I,I) = HFACT + BO(I,I) 74 ENDDO 75 76 ELSE 77 78 RHOUP = -BOVER2M*HUP 79 RHODOWN = -BOVER2M*HDOWN 80 81 HFACT = HALF + BOVER2M*CHEMPOT 82 83 DO I = 1, HDIM 84 RHOUP(I,I) = HFACT + RHOUP(I,I) 85 RHODOWN(I,I) = HFACT + RHODOWN(I,I) 86 ENDDO 87 88 ENDIF 89#ifdef GPUOFF 90 91 DO I = 1, FERMIM 92 93 IF (CGORLIB .EQ. 0) THEN 94 95 ! Call GESV-based solver 96 97 CALL SOLVEMATLAPACK 98 99 ELSE 100 101 ! Call the conjugate gradient solver on the CPU 102 103 CALL SOLVEMATCG 104 105 ENDIF 106 107 ENDDO 108 109#elif defined(GPUON) 110 111 ! 112 ! This calls Sanville's CUDA routines on the GPU 113 ! Now modified by MJC to send down FERMIM too 114 ! 115 116 CALL SOLVE_MATRIX_CG(BO, RHOUP, RHODOWN, HDIM, CGTOL2, & 117 SPINON, LATTEPREC, FERMIM) 118 119#endif 120 121 ! Modifying chemical potential 122 123 TRX = ZERO 124 TRXOMX = ZERO 125 126 IF (SPINON .EQ. 0) THEN 127 128 DO I = 1, HDIM 129 DO J = I, HDIM 130 131 IF (I .EQ. J) THEN 132 133 TRXOMX = TRXOMX + BO(I,I)*(ONE - BO(I,I)) 134 135 ELSE 136 137 TRXOMX = TRXOMX - TWO*BO(J,I)*BO(J,I) 138 139 ENDIF 140 141 ENDDO 142 143 TRX = TRX + BO(I,I) 144 145 ENDDO 146 147 SHIFTCP = KBT*(OCC - TRX)/TRXOMX 148 149 PREVERROR3 = PREVERROR2 150 PREVERROR2 = PREVERROR 151 PREVERROR = OCCERROR 152 OCCERROR = ABS(OCC - TRX) 153 ! PRINT*, ITER, OCCERROR 154 155 ELSE 156 157 DO I = 1, HDIM 158 DO J = I, HDIM 159 160 IF (I .EQ. J) THEN 161 162 TRXOMX = TRXOMX + RHOUP(I,I)*(ONE - RHOUP(I,I)) + & 163 RHODOWN(I,I)*(ONE - RHODOWN(I,I)) 164 165 ELSE 166 167 TRXOMX = TRXOMX - TWO*(RHOUP(J,I)*RHOUP(J,I) + & 168 RHODOWN(J,I)*RHODOWN(J,I)) 169 170 ENDIF 171 172 ENDDO 173 174 TRX = TRX + RHOUP(I,I) + RHODOWN(I,I) 175 176 ENDDO 177 178 SHIFTCP = KBT*(TOTNE - TRX)/TRXOMX 179 180 PREVERROR3 = PREVERROR2 181 PREVERROR2 = PREVERROR 182 PREVERROR = OCCERROR 183 OCCERROR = ABS(TOTNE - TRX) 184 185 ENDIF 186 187 ! PRINT*, ITER, OCCERROR 188 189 IF (ABS(SHIFTCP) .GT. MAXSHIFT) THEN 190 SHIFTCP = SIGN(MAXSHIFT, SHIFTCP) 191 ENDIF 192 193 CHEMPOT = CHEMPOT + SHIFTCP 194 195#ifdef DOUBLEPREC 196 197 IF (ITER .GE. 3 .AND. OCCERROR .LT. BREAKTOL) THEN 198 BREAKLOOP = 1 199 ENDIF 200 201#elif defined(SINGLEPREC) 202 203 IF ( ITER .GE. 3 .AND. ((OCCERROR .LT. BREAKTOL) .AND. & 204 (OCCERROR .EQ. PREVERROR .OR. & 205 OCCERROR .EQ. PREVERROR2 .OR. & 206 OCCERROR .EQ. PREVERROR3)) .OR. & 207 ITER .EQ. 25) THEN 208 209 BREAKLOOP = 1 210 211 ENDIF 212 213#endif 214 215 ENDDO 216 217 IF (SPINON .EQ. 0) THEN 218 BO = TWO*BO 219 ENDIF 220 221 RETURN 222 223END SUBROUTINE FERMIEXPANS 224