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