hdiff output

r31606/Cv.BS2.f90 2016-12-09 11:30:06.699661182 +0000 r31605/Cv.BS2.f90 2016-12-09 11:30:06.959664528 +0000
  1: !  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/UTILS/GMIN/Cv.BS2.f90' in revision 31605
  2: !  Calculate Cv curve from weights and visit statistics 
  3: !  for replicas at different temperatures REPT(J1). 
  4: !  Minimise chi^2 statistic to extract best probabilities. 
  5: ! 
  6: PROGRAM CVBS 
  7: IMPLICIT NONE 
  8: DOUBLE PRECISION TMIN, TMAX, TINT, Z0, Z1, Z2, TEMPERATURE, DUMMY, DELTA, ONEMEXP, PSUM, DP, PSUM2, TANAL, DPLUS, DMINUS 
  9: DOUBLE PRECISION COLOURTHRESH, HSHIFT, DUMMY2, DWEIGHT, DMIN, DUMMY3, DUMMY4, LN2DSUM 
 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:), EMIN(:), LNW2D(:,:) 
 11: DOUBLE PRECISION, ALLOCATABLE :: QENERGY(:), CQP(:), CQM(:) 
 12: INTEGER, ALLOCATABLE :: CLOSEST(:), QBIN(:) 
 13: INTEGER, ALLOCATABLE :: MINP(:), MINM(:) 
 14: INTEGER J1, NTEMP, J2, KAPPA, NPLUS, NMINUS, NDUMMY, NMIN, MAXBIN, NQBINS, NCOUNT 
 15: LOGICAL YESNO, DO2DT 
 16: LOGICAL, ALLOCATABLE :: NO1D(:),NO2D(:,:) 
 17:  
 18: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD') 
 19: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, MAXBIN 
 20: CLOSE(10) 
 21: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP)) 
 22:  
 23: PRINT '(3(A,I10))','bins=',MAXBIN 
 24: PRINT '(A,2G15.5,A,I8)','minimum and maximum T for Cv calculation=',TMIN,TMAX,' number of Cv data points=',NTEMP 
 25: PRINT '(A,I8)','number of degrees of freedom=',KAPPA 
 26:  
 27: ALLOCATE(ENERGY(MAXBIN),LNWEIGHT(MAXBIN),EMIN(NMIN),QBIN(NMIN),NO1D(MAXBIN)) 
 28:  
 29: OPEN(UNIT=10,FILE='min.data',STATUS='OLD') 
 30: DO J1=1,NMIN 
 31:    READ(10,*) EMIN(J1) 
 32: !  PRINT '(A,I8,G20.10)','J1,emin=',J1,EMIN(J1) 
 33: ENDDO 
 34: CLOSE(10) 
 35:  
 36: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD') 
 37: NO1D(1:MAXBIN)=.TRUE. 
 38: DO J1=1,MAXBIN 
 39:    READ(10,*,END=432) NDUMMY,DUMMY,DUMMY2,DUMMY3,DUMMY4 
 40:    ENERGY(NDUMMY)=DUMMY2 
 41:    LNWEIGHT(NDUMMY)=DUMMY4 
 42:    NO1D(NDUMMY)=.FALSE. 
 43: ENDDO 
 44: 432 CONTINUE 
 45: CLOSE(10) 
 46:  
 47: TINT=(TMAX-TMIN)/(NTEMP-1) 
 48: DELTA=ENERGY(2)-ENERGY(1) 
 49: ! 
 50: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa, 
 51: !  which cancel. 
 52: ! 
 53: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN') 
 54: DO J1=1,NTEMP 
 55:    Z0=0.0D0 
 56:    Z1=0.0D0 
 57:    Z2=0.0D0 
 58:    TEMPERATURE=TMIN+(J1-1)*TINT 
 59:    DO J2=1,MAXBIN 
 60:       IF (NO1D(J2)) CYCLE 
 61:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
 62:       Z0=Z0+DUMMY 
 63:       Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1)) 
 64:       Z2=Z2+DUMMY*(ENERGY(J2)-ENERGY(1))**2 
 65:    ENDDO 
 66:    MEANE(J1)=Z1/Z0 
 67:    ZSAVE(J1)=Z0 
 68:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN 
 69:       ONEMEXP=-DELTA/TEMPERATURE 
 70:    ELSE 
 71:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE) 
 72:    ENDIF 
 73:    WRITE(1,'(7G20.10)') TEMPERATURE, Z0, Z1, Z2, & 
 74:   &                     KAPPA/2.0D0 + & 
 75:   &                      1.0D0*(1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)) & 
 76:   &                     - (Z1/(Z0*TEMPERATURE))**2 + Z2/(Z0*TEMPERATURE**2), MEANE(J1), ZSAVE(J1) 
 77: ENDDO 
 78: CLOSE(1) 
 79:  
 80: OPEN(UNIT=1,FILE='prob.out.BS',STATUS='UNKNOWN') 
 81: DO J1=1,NTEMP 
 82:    TEMPERATURE=TMIN+(J1-1)*TINT 
 83:    PSUM=0.0D0 
 84:    PSUM2=0.0D0 
 85:    DO J2=1,MAXBIN 
 86:       IF (NO1D(J2)) CYCLE 
 87:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
 88:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2) 
 89:       PSUM=PSUM+DP*(ENERGY(J2)-ENERGY(1)-MEANE(J1)) 
 90: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP 
 91:       PSUM2=PSUM2+ZSAVE(J1)*(TEMPERATURE*DP)**2/DUMMY 
 92: !     PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2 
 93:    ENDDO 
 94:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN 
 95:       ONEMEXP=-DELTA/TEMPERATURE 
 96:    ELSE 
 97:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE) 
 98:    ENDIF 
 99:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM, & 
100:   &                                  KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM2 
101: ENDDO 
102: CLOSE(1) 
103:  
104: INQUIRE(FILE='Tanal.BS',EXIST=YESNO) 
105:  
106: IF (YESNO) THEN 
107:    NPLUS=0.0D0 
108:    NMINUS=0.0D0 
109:    DPLUS=0.0D0 
110:    DMINUS=0.0D0 
111:    ALLOCATE(PERMINP(MAXBIN),PERMINM(MAXBIN),MINP(MAXBIN),MINM(MAXBIN)) 
112:    OPEN(UNIT=1,FILE='Tanal.BS') 
113:    READ(1,*) TANAL 
114:    COLOURTHRESH=0.9D0 
115:    READ(1,*,END=666) COLOURTHRESH 
116: 666 CONTINUE 
117:    PRINT '(2(A,G20.10))','Analysing heat capacity contributions at T=',TANAL,' colour threshold=',COLOURTHRESH 
118:    CLOSE(1) 
119:    TEMPERATURE=TANAL 
120:    Z0=0.0D0 
121:    Z1=0.0D0 
122:    DO J2=1,MAXBIN 
123:       IF (NO1D(J2)) CYCLE 
124:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
125:       Z0=Z0+DUMMY 
126:       Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1)) 
127:    ENDDO 
128:    MEANE(1)=Z1/Z0 
129:    ZSAVE(1)=Z0 
130:    DO J2=1,MAXBIN 
131:       IF (NO1D(J2)) CYCLE 
132:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
133:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(1))/(ZSAVE(1)*TEMPERATURE**2) 
134:       DUMMY=DP*(ENERGY(J2)-ENERGY(1)-MEANE(1)) 
135:       IF (DP.LT.0.0D0) THEN 
136:          NMINUS=NMINUS+1 
137:          MINM(NMINUS)=J2 
138:          PERMINM(NMINUS)=DUMMY 
139:          DMINUS=DMINUS+DUMMY 
140:       ELSE 
141:          NPLUS=NPLUS+1 
142:          MINP(NPLUS)=J2 
143:          PERMINP(NPLUS)=DUMMY 
144:          DPLUS=DPLUS+DUMMY 
145:       ENDIF 
146: !     PRINT '(A,I6,3G20.10,2I6)','J2,ENERGY(J2)-ENERGY(1),MEANE,DP,NPLUS,NMINUS=',J2,ENERGY(J2)-ENERGY(1),MEANE(1),DP,NPLUS,NMINUS 
147:    ENDDO 
148:    DO J2=1,NMINUS 
149:       PERMINM(J2)=PERMINM(J2)/DMINUS 
150:    ENDDO 
151:    DO J2=1,NPLUS 
152:       PERMINP(J2)=PERMINP(J2)/DPLUS 
153:    ENDDO 
154:    CALL SORT(NMINUS,MAXBIN,PERMINM,MINM) 
155:    CALL SORT(NPLUS,MAXBIN,PERMINP,MINP) 
156:    DUMMY=0.0D0 
157:    PRINT '(A)','i bins with negative dp/dT below threshold:' 
158:    DO J2=1,NMINUS 
159:       DUMMY=DUMMY+PERMINM(J2) 
160:       PRINT '(2I6,4G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY,ENERGY(MINM(J2)),ENERGY(MINM(J2))-ENERGY(1) 
161:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
162:    ENDDO 
163:    DUMMY=0.0D0 
164:    PRINT '(A)','i bins with positive dp/dT:' 
165:    DO J2=1,NPLUS 
166:       DUMMY=DUMMY+PERMINP(J2) 
167:       PRINT '(2I6,4G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY,ENERGY(MINP(J2)),ENERGY(MINP(J2))-ENERGY(1) 
168:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
169:    ENDDO 
170:  
171:    INQUIRE(FILE='weights.2D.BS',EXIST=DO2DT) 
172:    IF (DO2DT) THEN 
173:       PRINT '(A)','Reading 2D weights from weights.2D.BS' 
174:       OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD') 
175: ! 
176: ! This itme, look for NQBINS parameter at the end for allocations 
177: ! 
178:       READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, MAXBIN, NQBINS 
179:       CLOSE(10) 
180:       ALLOCATE(LNW2D(NQBINS,MAXBIN),NO2D(NQBINS,MAXBIN),CQP(NQBINS),CQM(NQBINS),QENERGY(NQBINS)) 
181:       LNW2D(1:NQBINS,1:MAXBIN)=-HUGE(1.0D0) 
182:       NO2D(1:NQBINS,1:MAXBIN)=.TRUE. 
183:       OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD') 
184:       NCOUNT=0 
185:       DUMMY=0.0D0 
186:       DO  
187:          READ(19,*,END=753) J1, J2, DUMMY2, DUMMY, DWEIGHT 
188:          QENERGY(J1)=DUMMY2 
189:          NO2D(J1,J2)=.FALSE. 
190:          LNW2D(J1,J2)=DWEIGHT 
191:          DUMMY=DUMMY+EXP(DWEIGHT) 
192:       ENDDO 
193: 753   CONTINUE 
194:       LN2DSUM=LOG(DUMMY) 
195:       CLOSE(19) 
196:  
197:       DO J2=1,NQBINS 
198:          CQP(J2)=0.0D0 
199:          CQM(J2)=0.0D0 
200:          DO J1=1,NPLUS 
201:             IF (NO2D(J2,MINP(J1))) CYCLE 
202:             IF (NO1D(MINP(J1))) CYCLE 
203:             CQP(J2)=CQP(J2)+PERMINP(J1)*DPLUS*EXP(LNW2D(J2,MINP(J1))-LN2DSUM-LNWEIGHT(MINP(J1))) 
204:          ENDDO 
205:          DO J1=1,NMINUS 
206:             IF (NO2D(J2,MINM(J1))) CYCLE 
207:             IF (NO1D(MINM(J1))) CYCLE 
208:             CQM(J2)=CQM(J2)+PERMINM(J1)*DMINUS*EXP(LNW2D(NCOUNT,MINM(J1))-LN2DSUM-LNWEIGHT(MINM(J1))) 
209:          ENDDO 
210:       ENDDO 
211: ! 
212: ! Breakdown of Cv over q bins at TANAL 
213: ! 
214:       PRINT '(A,G20.10,A)','Breakdown of Cv over q bins at T=',TEMPERATURE,' q bin, plus, minus, sum, cumulative Cv' 
215:       DUMMY2=0.0D0 
216:       DO J2=1,NQBINS 
217:          DUMMY2=DUMMY2+CQP(J2)+CQM(J2) 
218:          PRINT '(I8,4G20.10)',J2,CQP(J2),CQM(J2),CQP(J2)+CQM(J2),DUMMY2 
219:       ENDDO 
220:     
221:       STOP 
222:              
223:    ENDIF 
224:  
225:    IF (DO2DT) THEN 
226:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.anhshift.BS min.plus.anhshift.BS and ' 
227:       PRINT '(A)','CHOOSECOLOURS in the dinfo file' 
228:       PRINT '(A)','Assignment of minima to quench bins: minimum, pe, q bin, q energy, diff' 
229:       DO J1=1,NMIN 
230:          DMIN=HUGE(1.0D0) 
231:          DO J2=1,NQBINS 
232:             IF (ABS(EMIN(J1)-QENERGY(J2)).LT.DMIN) THEN 
233:                QBIN(J1)=J2 
234:                DMIN=ABS(EMIN(J1)-QENERGY(J2)) 
235:             ENDIF 
236:          ENDDO 
237:          PRINT '(I8,G20.10,I8,2G20.10)',J1,EMIN(J1),QBIN(J1),QENERGY(QBIN(J1)),DMIN 
238:       ENDDO 
239:  
240: !     OPEN(UNIT=1,FILE='min.minus.anhshift.BS',STATUS='UNKNOWN') 
241: !     DUMMY=0.0D0 
242: !     DO J2=1,NMINUS 
243: !        DUMMY=DUMMY+PERMINM(J2) 
244: !        DO J1=1,NMIN 
245: !           IF (QTOI(QBIN(J1)).EQ.MINM(J2)) THEN 
246: !              WRITE(1,'(I6)') J1 
247: !              WRITE(*,'(A,3I8,3G20.10)') 'J1, J2, MINM(J2), EMIN, diff=', & 
248: ! &                                   J1, J2, MINM(J2), EMIN(J1), & 
249: ! &                                   ABS(EMIN(J1)-ENERGY(MINM(J2))+ANSHIFT(QBIN(J1))) 
250: !           ENDIF 
251: !        ENDDO 
252: !        IF (DUMMY.GT.COLOURTHRESH) EXIT 
253: !     ENDDO 
254: !     CLOSE(1) 
255:  
256: !     OPEN(UNIT=1,FILE='min.plus.anhshift.BS',STATUS='UNKNOWN') 
257: !     DUMMY=0.0D0 
258: !     DO J2=1,NPLUS 
259: !        DUMMY=DUMMY+PERMINP(J2) 
260: !        DO J1=1,NMIN 
261: !           IF (QTOI(QBIN(J1)).EQ.MINP(J2)) THEN 
262: !              WRITE(1,'(I6)') J1 
263: !              WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, diff=', & 
264: ! &                                   J1, J2, MINP(J2), EMIN(J1), & 
265: ! &                                   ABS(EMIN(J1)-ENERGY(MINP(J2))+ANSHIFT(QBIN(J1))) 
266: !           ENDIF 
267: !        ENDDO 
268: !        IF (DUMMY.GT.COLOURTHRESH) EXIT 
269: !     ENDDO  
270: !     CLOSE(1) 
271:  
272:    ENDIF 
273: ENDIF 
274:  
275: END PROGRAM CVBS 
276:  
277: SUBROUTINE SORT(N,J3,A,NA) 
278: IMPLICIT NONE 
279: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP 
280: DOUBLE PRECISION TEMP, A(J3) 
281: ! 
282: DO 20 J1=1,N-1 
283:    L=J1 
284:    DO 10 J2=J1+1,N 
285:       IF (A(L).LT.A(J2)) L=J2 
286: 10 CONTINUE 
287:    TEMP=A(L) 
288:    A(L)=A(J1) 
289:    A(J1)=TEMP 
290:    NTEMP=NA(L) 
291:    NA(L)=NA(J1) 
292:    NA(J1)=NTEMP 
293: 20 CONTINUE 
294: RETURN 
295: END SUBROUTINE SORT 
296:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0