hdiff output

r31583/Cv.BS.f90 2016-12-05 15:30:12.525346985 +0000 r31582/Cv.BS.f90 2016-12-05 15:30:12.777350235 +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
 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:), EMIN(:) 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:)
 11: DOUBLE PRECISION, ALLOCATABLE :: MEANVQI(:), MEANQI(:), ANSHIFT(:) 11: INTEGER, ALLOCATABLE :: MINP(:), MINM(:)
 12: INTEGER, ALLOCATABLE :: CLOSEST(:) 12: INTEGER J1, NTEMP, J2, KAPPA, NBINS, NPLUS, NMINUS, NDUMMY
 13: INTEGER, ALLOCATABLE :: MINP(:), MINM(:), IBININDEX(:) 13: LOGICAL YESNO
 14: INTEGER J1, NTEMP, J2, KAPPA, NBINS, NPLUS, NMINUS, NDUMMY, NMIN, MAXBIN 
 15: LOGICAL YESNO, DO2DT 
 16:  14: 
 17: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD') 15: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD')
 18: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, NBINS, MAXBIN 16: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NBINS
 19: CLOSE(10) 17: CLOSE(10)
 20: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP)) 18: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP))
 21:  19: 
 22: PRINT '(3(A,I10))','bins=',NBINS 20: 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 21: 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 22: PRINT '(A,I8)','number of degrees of freedom=',KAPPA
 25:  23: 
 26: ALLOCATE(ENERGY(NBINS),LNWEIGHT(NBINS),EMIN(NMIN),IBININDEX(MAXBIN)) 24: ALLOCATE(ENERGY(NBINS),LNWEIGHT(NBINS))
 27:  
 28: OPEN(UNIT=10,FILE='min.data',STATUS='OLD') 
 29: DO J1=1,NMIN 
 30:    READ(10,*) EMIN(J1) 
 31: ENDDO 
 32: CLOSE(10) 
 33:  25: 
  26: ! OPEN(UNIT=10,FILE='min.data',STATUS='OLD')
  27: ! DO J1=1,NMIN
  28: !    READ(10,*) ENERGY(J1),LNFRQS(J1),HORDER(J1)
  29: ! ENDDO
  30: ! CLOSE(10)
 34: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD') 31: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD')
 35: DO J1=1,NBINS 32: DO J1=1,NBINS
 36:    READ(10,*) NDUMMY,DUMMY,ENERGY(J1),DUMMY,LNWEIGHT(J1) 33:    READ(10,*) NDUMMY,DUMMY,ENERGY(J1),DUMMY,LNWEIGHT(J1)
 37:    IBININDEX(NDUMMY)=J1 
 38: ENDDO 34: ENDDO
 39: CLOSE(10) 35: CLOSE(10)
 40: INQUIRE(FILE='weights.2D.BS',EXIST=DO2DT) 
 41: IF (DO2DT) THEN 
 42:   PRINT '(A)','Reading 2D weights from weights.2D.BS' 
 43:   ALLOCATE(MEANVQI(NBINS),MEANQI(NBINS),ANSHIFT(NBINS),CLOSEST(NMIN)) 
 44:   MEANVQI(1:NBINS)=0.0D0 
 45:   MEANQI(1:NBINS)=0.0D0 
 46:   OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD') 
 47:   DO  
 48:      READ(19,*,END=753) J1, J2, DUMMY2, DUMMY, DWEIGHT 
 49:      MEANVQI(IBININDEX(J2))=MEANVQI(IBININDEX(J2))+EXP(DWEIGHT)*DUMMY2 
 50:      MEANQI(IBININDEX(J2))=MEANQI(IBININDEX(J2))+EXP(DWEIGHT) 
 51:   ENDDO 
 52: 753 CONTINUE 
 53:   CLOSE(19) 
 54:   PRINT '(A)','instantaneous pe bin,    energy,    weight,    weighted quench energy,    mean quench pe,       shift' 
 55:   DO J1=1,NBINS 
 56:      ANSHIFT(J1)=ENERGY(J1)-MEANVQI(J1)/MAX(MEANQI(J1),1.0D-300) 
 57:      PRINT '(I8,5G20.10)',J1,ENERGY(J1),MEANQI(J1),MEANVQI(J1),MEANVQI(J1)/MAX(MEANQI(J1),1.0D-300),ANSHIFT(J1) 
 58:   ENDDO 
 59: ENDIF 
 60:  36: 
 61: TINT=(TMAX-TMIN)/(NTEMP-1) 37: TINT=(TMAX-TMIN)/(NTEMP-1)
 62: DELTA=ENERGY(2)-ENERGY(1) 38: DELTA=ENERGY(2)-ENERGY(1)
 63: ! 39: !
 64: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa, 40: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa,
 65: !  which cancel. 41: !  which cancel.
 66: ! 42: !
 67: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN') 43: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN')
 68: DO J1=1,NTEMP 44: DO J1=1,NTEMP
 69:    Z0=0.0D0 45:    Z0=0.0D0
113: CLOSE(1) 89: CLOSE(1)
114:  90: 
115: INQUIRE(FILE='Tanal.BS',EXIST=YESNO) 91: INQUIRE(FILE='Tanal.BS',EXIST=YESNO)
116:  92: 
117: IF (YESNO) THEN 93: IF (YESNO) THEN
118:    NPLUS=0.0D0 94:    NPLUS=0.0D0
119:    NMINUS=0.0D0 95:    NMINUS=0.0D0
120:    DPLUS=0.0D0 96:    DPLUS=0.0D0
121:    DMINUS=0.0D0 97:    DMINUS=0.0D0
122:    ALLOCATE(PERMINP(NBINS),PERMINM(NBINS),MINP(NBINS),MINM(NBINS)) 98:    ALLOCATE(PERMINP(NBINS),PERMINM(NBINS),MINP(NBINS),MINM(NBINS))
123:    OPEN(UNIT=1,FILE='Tanal.BS') 99:    OPEN(UNIT=1,FILE='Tanal')
124:    READ(1,*) TANAL100:    READ(1,*) TANAL
125:    COLOURTHRESH=0.9D0101:    COLOURTHRESH=0.9D0
126:    READ(1,*,END=666) COLOURTHRESH102:    READ(1,*,END=666) COLOURTHRESH
127: 666 CONTINUE103: 666 CONTINUE
128:    PRINT '(2(A,G20.10))','Analysing heat capacity contributions at T=',TANAL,' colour threshold=',COLOURTHRESH104:    PRINT '(2(A,G20.10))','Analysing heat capacity contributions at T=',TANAL,' colour threshold=',COLOURTHRESH
129:    CLOSE(1)105:    CLOSE(1)
130:    TEMPERATURE=TANAL106:    TEMPERATURE=TANAL
131:    Z0=0.0D0107:    Z0=0.0D0
132:    Z1=0.0D0108:    Z1=0.0D0
133:    DO J2=1,NBINS109:    DO J2=1,NBINS
167:    DO J2=1,NMINUS143:    DO J2=1,NMINUS
168:       DUMMY=DUMMY+PERMINM(J2)144:       DUMMY=DUMMY+PERMINM(J2)
169:       PRINT '(2I6,2G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY145:       PRINT '(2I6,2G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY
170:    ENDDO146:    ENDDO
171:    DUMMY=0.0D0147:    DUMMY=0.0D0
172:    PRINT '(A)','bins with positive dp/dT:'148:    PRINT '(A)','bins with positive dp/dT:'
173:    DO J2=1,NPLUS149:    DO J2=1,NPLUS
174:       DUMMY=DUMMY+PERMINP(J2)150:       DUMMY=DUMMY+PERMINP(J2)
175:       PRINT '(2I6,2G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY151:       PRINT '(2I6,2G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY
176:    ENDDO152:    ENDDO
177: !153:    PRINT *,'need to translate +/- bins list into minima indices and write min.minus.BS min.plus.BS'
178: ! Harmonic shift to map i bins to q bins. V_i = V_q + kappa*k*T/2154:    PRINT *,'need to count the NMIN minima'
179: !155: !  PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.BS min.plus.BS and '
180:    HSHIFT=KAPPA*TANAL/2.0D0 
181:    PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.hshift.BS min.plus.hshift.BS and ' 
182:    PRINT '(A)','CHOOSECOLOURS in the dinfo file'156:    PRINT '(A)','CHOOSECOLOURS in the dinfo file'
183:    OPEN(UNIT=1,FILE='min.minus.hshift.BS',STATUS='UNKNOWN')157:    OPEN(UNIT=1,FILE='min.minus',STATUS='UNKNOWN')
184:  
185:    DUMMY=0.0D0158:    DUMMY=0.0D0
186:    DO J2=1,NMINUS159:    DO J2=1,NMINUS
187:       DUMMY=DUMMY+PERMINM(J2)160:       DUMMY=DUMMY+PERMINM(J2)
188:       DO J1=1,NMIN161:       WRITE(1,'(I6)') MINM(J2)
189:          IF (ABS(EMIN(J1)-ENERGY(MINM(J2))+HSHIFT).LT.DELTA/2.0D0) THEN 
190:             WRITE(1,'(I6)') J1 
191:             WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINM(J2), EMIN, ENERGY(MINM(J2))-HSHIFT, abs=',J1, J2, MINM(J2), EMIN(J1), & 
192:   &                                     ENERGY(MINM(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINM(J2))+HSHIFT) 
193:          ENDIF 
194:       ENDDO 
195:       IF (DUMMY.GT.COLOURTHRESH) EXIT162:       IF (DUMMY.GT.COLOURTHRESH) EXIT
196:    ENDDO163:    ENDDO
197:    CLOSE(1)164:    CLOSE(1)
198:    OPEN(UNIT=1,FILE='min.plus.hshift.BS',STATUS='UNKNOWN')165:    OPEN(UNIT=1,FILE='min.plus',STATUS='UNKNOWN')
199:    DMIN=HUGE(1.0D0) 
200:    DO J1=1,NMIN 
201:    ENDDO 
202:  
203:    DUMMY=0.0D0166:    DUMMY=0.0D0
204:    DO J2=1,NPLUS167:    DO J2=1,NPLUS
205:       DUMMY=DUMMY+PERMINP(J2)168:       DUMMY=DUMMY+PERMINP(J2)
206:       DO J1=1,NMIN169:       WRITE(1,'(I6)') MINP(J2)
207:          IF (ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT).LT.DELTA/2.0D0) THEN 
208:             WRITE(1,'(I6)') J1 
209:             WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, ENERGY(MINP(J2))-HSHIFT, abs=',J1, J2, MINP(J2), EMIN(J1), & 
210:   &                                     ENERGY(MINP(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT) 
211:          ENDIF 
212:       ENDDO 
213:       IF (DUMMY.GT.COLOURTHRESH) EXIT170:       IF (DUMMY.GT.COLOURTHRESH) EXIT
214:    ENDDO171:    ENDDO
215:    CLOSE(1)172:    CLOSE(1)
216:  
217:    IF (DO2DT) THEN 
218:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.anhshift.BS min.plus.anhshift.BS and ' 
219:       PRINT '(A)','CHOOSECOLOURS in the dinfo file' 
220:       OPEN(UNIT=1,FILE='min.minus.anhshift.BS',STATUS='UNKNOWN') 
221:       DO J1=1,NMIN 
222:            DMIN=HUGE(1.0D0) 
223:            DO J2=1,NBINS 
224:               IF (ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2)).LT.DMIN) THEN 
225:                  DMIN=ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2)) 
226:                  CLOSEST(J1)=J2 
227: !                PRINT '(A,I6,A,I6,A,3G20.10)','changing closest i bin for min ',J1,' to ',J2,' emin, ei, diff=', & 
228: ! &                     EMIN(J1),ENERGY(J2),ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2)) 
229:               ENDIF 
230:               PRINT '(A,I6,A,I6,A,3G20.10)','min ',J1,' ibin ',J2,' emin, ei, diff=', & 
231:   &                  EMIN(J1),ENERGY(J2),ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2)) 
232:            ENDDO 
233:            PRINT '(A,I6,A,I6,A,3G20.10)','closest i bin for min ',J1,' is ',CLOSEST(J1),' emin, ei, diff=', & 
234:   &               EMIN(J1),ENERGY(CLOSEST(J1)),ABS(EMIN(J1)-ENERGY(CLOSEST(J1))+ANSHIFT(CLOSEST(J1))) 
235:         ENDDO 
236:  
237:         DUMMY=0.0D0 
238:         DO J2=1,NMINUS 
239:            DUMMY=DUMMY+PERMINM(J2) 
240:            DO J1=1,NMIN 
241:               IF (CLOSEST(J1).EQ.MINM(J2)) THEN 
242:                  WRITE(1,'(I6)') J1 
243:                  WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINM(J2), EMIN, ENERGY(MINM(J2))-ANSHIFT(MINM(J2)), abs=', & 
244:   &                                     J1, J2, MINM(J2), EMIN(J1), & 
245:   &                                     ENERGY(MINM(J2))-ANSHIFT(MINM(J2)), ABS(EMIN(J1)-ENERGY(MINM(J2))+ANSHIFT(MINM(J2))) 
246:               ENDIF 
247:            ENDDO 
248:            IF (DUMMY.GT.COLOURTHRESH) EXIT 
249:         ENDDO 
250:         CLOSE(1) 
251:         OPEN(UNIT=1,FILE='min.plus.anhshift.BS',STATUS='UNKNOWN') 
252:  
253:         DUMMY=0.0D0 
254:         DO J2=1,NPLUS 
255:            DUMMY=DUMMY+PERMINP(J2) 
256:            DO J1=1,NMIN 
257:               IF (CLOSEST(J1).EQ.MINP(J2)) THEN 
258:                  WRITE(1,'(I6)') J1 
259:                  WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, ENERGY(MINP(J2))-ANSHIFT(MINP(J2)), abs=', & 
260:   &                                     J1, J2, MINP(J2), EMIN(J1), & 
261:   &                                     ENERGY(MINP(J2))-ANSHIFT(MINP(J2)), ABS(EMIN(J1)-ENERGY(MINP(J2))+ANSHIFT(MINP(J2))) 
262:               ENDIF 
263:            ENDDO 
264:            IF (DUMMY.GT.COLOURTHRESH) EXIT 
265:         ENDDO  
266:         CLOSE(1) 
267:  
268:    ENDIF 
269: ENDIF173: ENDIF
270: 174: 
271: END PROGRAM CVBS175: END PROGRAM CVBS
272: 176: 
273: SUBROUTINE SORT(N,J3,A,NA)177: SUBROUTINE SORT(N,J3,A,NA)
274: IMPLICIT NONE178: IMPLICIT NONE
275: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP179: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP
276: DOUBLE PRECISION TEMP, A(J3)180: DOUBLE PRECISION TEMP, A(J3)
277: !181: !
278: DO 20 J1=1,N-1182: DO 20 J1=1,N-1


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0