hdiff output

r31592/Cv.BS.f90 2016-12-06 13:30:10.258520851 +0000 r31591/Cv.BS.f90 2016-12-06 13:30:10.506524056 +0000
  1: !  1: !
  2: !  Calculate Cv curve from weights and visit statistics  2: !  Calculate Cv curve from weights and visit statistics
  3: !  for replicas at different temperatures REPT(J1).  3: !  for replicas at different temperatures REPT(J1).
  4: !  Minimise chi^2 statistic to extract best probabilities.  4: !  Minimise chi^2 statistic to extract best probabilities.
  5: !  5: !
  6: PROGRAM CVBS  6: PROGRAM CVBS
  7: IMPLICIT NONE  7: IMPLICIT NONE
  8: DOUBLE PRECISION TMIN, TMAX, TINT, Z0, Z1, Z2, TEMPERATURE, DUMMY, DELTA, ONEMEXP, PSUM, DP, PSUM2, TANAL, DPLUS, DMINUS  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  9: DOUBLE PRECISION COLOURTHRESH, HSHIFT, DUMMY2, DWEIGHT, DMIN
 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:), EMIN(:), LNWEIGHT2(:) 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:), EMIN(:)
 11: DOUBLE PRECISION, ALLOCATABLE :: MEANVQI(:), MEANQI(:), ANSHIFT(:) 11: DOUBLE PRECISION, ALLOCATABLE :: MEANVQI(:), MEANQI(:), ANSHIFT(:)
 12: INTEGER, ALLOCATABLE :: CLOSEST(:) 12: INTEGER, ALLOCATABLE :: CLOSEST(:)
 13: INTEGER, ALLOCATABLE :: MINP(:), MINM(:), IBININDEX(:) 13: INTEGER, ALLOCATABLE :: MINP(:), MINM(:), IBININDEX(:)
 14: INTEGER J1, NTEMP, J2, KAPPA, NBINS, NPLUS, NMINUS, NDUMMY, NMIN, MAXBIN 14: INTEGER J1, NTEMP, J2, KAPPA, NBINS, NPLUS, NMINUS, NDUMMY, NMIN, MAXBIN
 15: LOGICAL YESNO, DO2DT 15: LOGICAL YESNO, DO2DT
 16:  16: 
 17: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD') 17: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD')
 18: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, NBINS, MAXBIN 18: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, NBINS, MAXBIN
 19: CLOSE(10) 19: CLOSE(10)
 20: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP)) 20: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP))
 21:  21: 
 22: PRINT '(3(A,I10))','bins=',NBINS 22: PRINT '(3(A,I10))','bins=',NBINS
 23: PRINT '(A,2G15.5,A,I8)','minimum and maximum T for Cv calculation=',TMIN,TMAX,' number of Cv data points=',NTEMP 23: PRINT '(A,2G15.5,A,I8)','minimum and maximum T for Cv calculation=',TMIN,TMAX,' number of Cv data points=',NTEMP
 24: PRINT '(A,I8)','number of degrees of freedom=',KAPPA 24: PRINT '(A,I8)','number of degrees of freedom=',KAPPA
 25:  25: 
 26: ALLOCATE(ENERGY(NBINS),LNWEIGHT2(NBINS),LNWEIGHT(NBINS),EMIN(NMIN),IBININDEX(MAXBIN)) 26: ALLOCATE(ENERGY(NBINS),LNWEIGHT(NBINS),EMIN(NMIN),IBININDEX(MAXBIN))
 27:  27: 
 28: OPEN(UNIT=10,FILE='min.data',STATUS='OLD') 28: OPEN(UNIT=10,FILE='min.data',STATUS='OLD')
 29: DO J1=1,NMIN 29: DO J1=1,NMIN
 30:    READ(10,*) EMIN(J1) 30:    READ(10,*) EMIN(J1)
 31: ENDDO 31: ENDDO
 32: CLOSE(10) 32: CLOSE(10)
 33:  33: 
 34: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD') 34: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD')
 35: DO J1=1,NBINS 35: DO J1=1,NBINS
 36:    READ(10,*) NDUMMY,DUMMY,ENERGY(J1),DUMMY,LNWEIGHT(J1) 36:    READ(10,*) NDUMMY,DUMMY,ENERGY(J1),DUMMY,LNWEIGHT(J1)
 37:    IBININDEX(NDUMMY)=J1 37:    IBININDEX(NDUMMY)=J1
 38: ENDDO 38: ENDDO
 39: CLOSE(10) 39: CLOSE(10)
 40: INQUIRE(FILE='weights.2D.BS',EXIST=DO2DT) 40: INQUIRE(FILE='weights.2D.BS',EXIST=DO2DT)
 41: IF (DO2DT) THEN 41: IF (DO2DT) THEN
 42:    PRINT '(A)','Reading 2D weights from weights.2D.BS' 42:   PRINT '(A)','Reading 2D weights from weights.2D.BS'
 43:    ALLOCATE(MEANVQI(NBINS),MEANQI(NBINS),ANSHIFT(NBINS),CLOSEST(NMIN)) 43:   ALLOCATE(MEANVQI(NBINS),MEANQI(NBINS),ANSHIFT(NBINS),CLOSEST(NMIN))
 44:    MEANVQI(1:NBINS)=0.0D0 44:   MEANVQI(1:NBINS)=0.0D0
 45:    MEANQI(1:NBINS)=0.0D0 45:   MEANQI(1:NBINS)=0.0D0
 46:    LNWEIGHT2(1:NBINS)=0.0D0 46:   OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD')
 47:    OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD') 47:   DO 
 48:    DO  48:      READ(19,*,END=753) J1, J2, DUMMY2, DUMMY, DWEIGHT
 49:       READ(19,*,END=753) J1, J2, DUMMY2, DUMMY, DWEIGHT 49:      MEANVQI(IBININDEX(J2))=MEANVQI(IBININDEX(J2))+EXP(DWEIGHT)*DUMMY2
 50:       MEANVQI(IBININDEX(J2))=MEANVQI(IBININDEX(J2))+EXP(DWEIGHT)*DUMMY2 50:      MEANQI(IBININDEX(J2))=MEANQI(IBININDEX(J2))+EXP(DWEIGHT)
 51:       MEANQI(IBININDEX(J2))=MEANQI(IBININDEX(J2))+EXP(DWEIGHT) 51:   ENDDO
 52:       LNWEIGHT2(IBININDEX(J2))=LNWEIGHT2(IBININDEX(J2))+EXP(DWEIGHT) 
 53:    ENDDO 
 54: 753 CONTINUE 52: 753 CONTINUE
 55:    CLOSE(19) 53:   CLOSE(19)
 56:    PRINT '(A)','instantaneous pe bin,    energy,    weight,    weighted quench energy,    mean quench pe,       shift' 54:   PRINT '(A)','instantaneous pe bin,    energy,    weight,    weighted quench energy,    mean quench pe,       shift'
 57:    DO J1=1,NBINS 55:   DO J1=1,NBINS
 58:       ANSHIFT(J1)=ENERGY(J1)-MEANVQI(J1)/MAX(MEANQI(J1),1.0D-300) 56:      ANSHIFT(J1)=ENERGY(J1)-MEANVQI(J1)/MAX(MEANQI(J1),1.0D-300)
 59:       PRINT '(I8,5G20.10)',J1,ENERGY(J1),MEANQI(J1),MEANVQI(J1),MEANVQI(J1)/MAX(MEANQI(J1),1.0D-300),ANSHIFT(J1) 57:      PRINT '(I8,5G20.10)',J1,ENERGY(J1),MEANQI(J1),MEANVQI(J1),MEANVQI(J1)/MAX(MEANQI(J1),1.0D-300),ANSHIFT(J1)
 60:    ENDDO 58:   ENDDO
 61:    PRINT '(A)','ln i bin weights from weights.BS and weights.2D.BS and difference' 
 62:    DUMMY=LOG(LNWEIGHT2(1)) 
 63:    DO J1=1,NBINS 
 64:       LNWEIGHT2(J1)=LOG(LNWEIGHT2(J1))-DUMMY 
 65:       PRINT '(I8,3G20.10)',J1,LNWEIGHT(J1),LNWEIGHT2(J1),LNWEIGHT(J1)-LNWEIGHT2(J1) 
 66:    ENDDO 
 67: ENDIF 59: ENDIF
 68:  60: 
 69: TINT=(TMAX-TMIN)/(NTEMP-1) 61: TINT=(TMAX-TMIN)/(NTEMP-1)
 70: DELTA=ENERGY(2)-ENERGY(1) 62: DELTA=ENERGY(2)-ENERGY(1)
 71: ! 63: !
 72: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa, 64: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa,
 73: !  which cancel. 65: !  which cancel.
 74: ! 66: !
 75: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN') 67: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN')
 76: DO J1=1,NTEMP 68: DO J1=1,NTEMP
216:             WRITE(1,'(I6)') J1208:             WRITE(1,'(I6)') J1
217:             WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, ENERGY(MINP(J2))-HSHIFT, abs=',J1, J2, MINP(J2), EMIN(J1), &209:             WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, ENERGY(MINP(J2))-HSHIFT, abs=',J1, J2, MINP(J2), EMIN(J1), &
218:   &                                     ENERGY(MINP(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT)210:   &                                     ENERGY(MINP(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT)
219:          ENDIF211:          ENDIF
220:       ENDDO212:       ENDDO
221:       IF (DUMMY.GT.COLOURTHRESH) EXIT213:       IF (DUMMY.GT.COLOURTHRESH) EXIT
222:    ENDDO214:    ENDDO
223:    CLOSE(1)215:    CLOSE(1)
224: 216: 
225:    IF (DO2DT) THEN217:    IF (DO2DT) THEN
226:       Z0=0.0D0 
227:       Z1=0.0D0 
228:       DO J2=1,NBINS 
229:          DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT2(J2)) 
230:          Z0=Z0+DUMMY 
231:          Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1)) 
232:       ENDDO 
233:       MEANE(1)=Z1/Z0 
234:       ZSAVE(1)=Z0 
235:       NMINUS=0 
236:       NPLUS=0 
237:       DPLUS=0.0D0 
238:       DMINUS=0.0D0 
239:       DO J2=1,NBINS 
240:          DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT2(J2)) 
241:          DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(1))/(ZSAVE(1)*TEMPERATURE**2) 
242:          DUMMY=DP*(ENERGY(J2)-ENERGY(1)-MEANE(1)) 
243:          IF (DP.LT.0.0D0) THEN 
244:             NMINUS=NMINUS+1 
245:             MINM(NMINUS)=J2 
246:             PERMINM(NMINUS)=DUMMY 
247:             DMINUS=DMINUS+DUMMY 
248:          ELSE 
249:             NPLUS=NPLUS+1 
250:             MINP(NPLUS)=J2 
251:             PERMINP(NPLUS)=DUMMY 
252:             DPLUS=DPLUS+DUMMY 
253:          ENDIF 
254: !     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 
255:       ENDDO 
256:       DO J2=1,NMINUS 
257:          PERMINM(J2)=PERMINM(J2)/DMINUS 
258:       ENDDO 
259:       DO J2=1,NPLUS 
260:          PERMINP(J2)=PERMINP(J2)/DPLUS 
261:       ENDDO 
262:       CALL SORT(NMINUS,NBINS,PERMINM,MINM) 
263:       CALL SORT(NPLUS,NBINS,PERMINP,MINP) 
264:       DUMMY=0.0D0 
265:       PRINT '(A)','bins with negative dp/dT:' 
266:       DO J2=1,NMINUS 
267:          DUMMY=DUMMY+PERMINM(J2) 
268:          PRINT '(2I6,2G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY 
269:       ENDDO 
270:       DUMMY=0.0D0 
271:       PRINT '(A)','bins with positive dp/dT:' 
272:       DO J2=1,NPLUS 
273:          DUMMY=DUMMY+PERMINP(J2) 
274:          PRINT '(2I6,2G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY 
275:       ENDDO 
276:  
277:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.anhshift.BS min.plus.anhshift.BS and '218:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.anhshift.BS min.plus.anhshift.BS and '
278:       PRINT '(A)','CHOOSECOLOURS in the dinfo file'219:       PRINT '(A)','CHOOSECOLOURS in the dinfo file'
279:       OPEN(UNIT=1,FILE='min.minus.anhshift.BS',STATUS='UNKNOWN')220:       OPEN(UNIT=1,FILE='min.minus.anhshift.BS',STATUS='UNKNOWN')
280:       DO J1=1,NMIN221:       DO J1=1,NMIN
281:            DMIN=HUGE(1.0D0)222:            DMIN=HUGE(1.0D0)
282:            DO J2=1,NBINS223:            DO J2=1,NBINS
283:               IF (ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2)).LT.DMIN) THEN224:               IF (ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2)).LT.DMIN) THEN
284:                  DMIN=ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2))225:                  DMIN=ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2))
285:                  CLOSEST(J1)=J2226:                  CLOSEST(J1)=J2
286: !                PRINT '(A,I6,A,I6,A,3G20.10)','changing closest i bin for min ',J1,' to ',J2,' emin, ei, diff=', &227: !                PRINT '(A,I6,A,I6,A,3G20.10)','changing closest i bin for min ',J1,' to ',J2,' emin, ei, diff=', &
287: ! &                     EMIN(J1),ENERGY(J2),ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2))228: ! &                     EMIN(J1),ENERGY(J2),ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2))
288:               ENDIF229:               ENDIF
289: !             PRINT '(A,I6,A,I6,A,3G20.10)','min ',J1,' ibin ',J2,' emin, ei, diff=', &230:               PRINT '(A,I6,A,I6,A,3G20.10)','min ',J1,' ibin ',J2,' emin, ei, diff=', &
290: ! &                  EMIN(J1),ENERGY(J2),ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2))231:   &                  EMIN(J1),ENERGY(J2),ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2))
291:            ENDDO232:            ENDDO
292:            PRINT '(A,I6,A,I6,A,3G20.10)','closest i bin for min ',J1,' is ',CLOSEST(J1),' emin, ei, diff=', &233:            PRINT '(A,I6,A,I6,A,3G20.10)','closest i bin for min ',J1,' is ',CLOSEST(J1),' emin, ei, diff=', &
293:   &               EMIN(J1),ENERGY(CLOSEST(J1)),ABS(EMIN(J1)-ENERGY(CLOSEST(J1))+ANSHIFT(CLOSEST(J1)))234:   &               EMIN(J1),ENERGY(CLOSEST(J1)),ABS(EMIN(J1)-ENERGY(CLOSEST(J1))+ANSHIFT(CLOSEST(J1)))
294:         ENDDO235:         ENDDO
295: 236: 
296:         DUMMY=0.0D0237:         DUMMY=0.0D0
297:         DO J2=1,NMINUS238:         DO J2=1,NMINUS
298:            DUMMY=DUMMY+PERMINM(J2)239:            DUMMY=DUMMY+PERMINM(J2)
299:            DO J1=1,NMIN240:            DO J1=1,NMIN
300:               IF (CLOSEST(J1).EQ.MINM(J2)) THEN241:               IF (CLOSEST(J1).EQ.MINM(J2)) THEN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0