hdiff output

r31604/Cv.BS.f90 2016-12-08 16:30:08.541320662 +0000 r31603/Cv.BS.f90 2016-12-08 16:30:08.793323906 +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(:), LNWEIGHT2(:)
 11: DOUBLE PRECISION, ALLOCATABLE :: MEANVQI(:), MEANQI(:), ANSHIFT(:), QENERGY(:) 11: DOUBLE PRECISION, ALLOCATABLE :: MEANVQI(:), MEANQI(:), ANSHIFT(:)
 12: INTEGER, ALLOCATABLE :: CLOSEST(:), QBIN(:), QTOI(:) 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, NQBINS, NCOUNT, NQBINSREAL, ND1 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
 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)
  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:    LNWEIGHT2(1:NBINS)=0.0D0
  47:    OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD')
  48:    DO 
  49:       READ(19,*,END=753) J1, J2, DUMMY2, DUMMY, DWEIGHT
  50:       MEANVQI(IBININDEX(J2))=MEANVQI(IBININDEX(J2))+EXP(DWEIGHT)*DUMMY2
  51:       MEANQI(IBININDEX(J2))=MEANQI(IBININDEX(J2))+EXP(DWEIGHT)
  52:       LNWEIGHT2(IBININDEX(J2))=LNWEIGHT2(IBININDEX(J2))+EXP(DWEIGHT)
  53:    ENDDO
  54: 753 CONTINUE
  55:    CLOSE(19)
  56:    PRINT '(A)','instantaneous pe bin,    energy,    weight,    weighted quench energy,    mean quench pe,       shift'
  57:    DO J1=1,NBINS
  58:       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)
  60:    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
 40:  68: 
 41: TINT=(TMAX-TMIN)/(NTEMP-1) 69: TINT=(TMAX-TMIN)/(NTEMP-1)
 42: DELTA=ENERGY(2)-ENERGY(1) 70: DELTA=ENERGY(2)-ENERGY(1)
 43: ! 71: !
 44: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa, 72: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa,
 45: !  which cancel. 73: !  which cancel.
 46: ! 74: !
 47: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN') 75: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN')
 48: DO J1=1,NTEMP 76: DO J1=1,NTEMP
 49:    Z0=0.0D0 77:    Z0=0.0D0
 73: OPEN(UNIT=1,FILE='prob.out.BS',STATUS='UNKNOWN')101: OPEN(UNIT=1,FILE='prob.out.BS',STATUS='UNKNOWN')
 74: DO J1=1,NTEMP102: DO J1=1,NTEMP
 75:    TEMPERATURE=TMIN+(J1-1)*TINT103:    TEMPERATURE=TMIN+(J1-1)*TINT
 76:    PSUM=0.0D0104:    PSUM=0.0D0
 77:    PSUM2=0.0D0105:    PSUM2=0.0D0
 78:    DO J2=1,NBINS106:    DO J2=1,NBINS
 79:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2))107:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2))
 80:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2)108:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2)
 81:       PSUM=PSUM+DP*(ENERGY(J2)-ENERGY(1)-MEANE(J1))109:       PSUM=PSUM+DP*(ENERGY(J2)-ENERGY(1)-MEANE(J1))
 82: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP110: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP
 83:       PSUM2=PSUM2+ZSAVE(J1)*(TEMPERATURE*DP)**2/DUMMY111:       PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2
 84: !     PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2 
 85:    ENDDO112:    ENDDO
 86:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN113:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN
 87:       ONEMEXP=-DELTA/TEMPERATURE114:       ONEMEXP=-DELTA/TEMPERATURE
 88:    ELSE115:    ELSE
 89:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE)116:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE)
 90:    ENDIF117:    ENDIF
 91:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM, &118:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM, &
 92:   &                                  KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM2119:   &                                  KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM2
 93: ENDDO120: ENDDO
 94: CLOSE(1)121: CLOSE(1)
137:    ENDDO164:    ENDDO
138:    DO J2=1,NMINUS165:    DO J2=1,NMINUS
139:       PERMINM(J2)=PERMINM(J2)/DMINUS166:       PERMINM(J2)=PERMINM(J2)/DMINUS
140:    ENDDO167:    ENDDO
141:    DO J2=1,NPLUS168:    DO J2=1,NPLUS
142:       PERMINP(J2)=PERMINP(J2)/DPLUS169:       PERMINP(J2)=PERMINP(J2)/DPLUS
143:    ENDDO170:    ENDDO
144:    CALL SORT(NMINUS,NBINS,PERMINM,MINM)171:    CALL SORT(NMINUS,NBINS,PERMINM,MINM)
145:    CALL SORT(NPLUS,NBINS,PERMINP,MINP)172:    CALL SORT(NPLUS,NBINS,PERMINP,MINP)
146:    DUMMY=0.0D0173:    DUMMY=0.0D0
147:    PRINT '(A)','i bins with negative dp/dT below threshold:'174:    PRINT '(A)','bins with negative dp/dT:'
148:    DO J2=1,NMINUS175:    DO J2=1,NMINUS
149:       DUMMY=DUMMY+PERMINM(J2)176:       DUMMY=DUMMY+PERMINM(J2)
150:       PRINT '(2I6,4G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY,ENERGY(MINM(J2)),ENERGY(MINM(J2))-ENERGY(1)177:       PRINT '(2I6,2G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY
151:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
152:    ENDDO178:    ENDDO
153:    DUMMY=0.0D0179:    DUMMY=0.0D0
154:    PRINT '(A)','i bins with positive dp/dT:'180:    PRINT '(A)','bins with positive dp/dT:'
155:    DO J2=1,NPLUS181:    DO J2=1,NPLUS
156:       DUMMY=DUMMY+PERMINP(J2)182:       DUMMY=DUMMY+PERMINP(J2)
157:       PRINT '(2I6,4G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY,ENERGY(MINP(J2)),ENERGY(MINP(J2))-ENERGY(1)183:       PRINT '(2I6,2G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY
158:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
159:    ENDDO184:    ENDDO
160: !185: !
161: ! Harmonic shift to map i bins to q bins. V_i = V_q + kappa*k*T/2186: ! Harmonic shift to map i bins to q bins. V_i = V_q + kappa*k*T/2
162: !187: !
163:    HSHIFT=KAPPA*TANAL/2.0D0188:    HSHIFT=KAPPA*TANAL/2.0D0
164:    PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.hshift.BS min.plus.hshift.BS and '189:    PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.hshift.BS min.plus.hshift.BS and '
165:    PRINT '(A)','CHOOSECOLOURS in the dinfo file'190:    PRINT '(A)','CHOOSECOLOURS in the dinfo file'
166:    OPEN(UNIT=1,FILE='min.minus.hshift.BS',STATUS='UNKNOWN')191:    OPEN(UNIT=1,FILE='min.minus.hshift.BS',STATUS='UNKNOWN')
167: 192: 
168:    DUMMY=0.0D0193:    DUMMY=0.0D0
169:    DO J2=1,NMINUS194:    DO J2=1,NMINUS
170:       DUMMY=DUMMY+PERMINM(J2)195:       DUMMY=DUMMY+PERMINM(J2)
171:       DO J1=1,NMIN196:       DO J1=1,NMIN
172:          IF (ABS(EMIN(J1)-ENERGY(MINM(J2))+HSHIFT).LT.DELTA/2.0D0) THEN197:          IF (ABS(EMIN(J1)-ENERGY(MINM(J2))+HSHIFT).LT.DELTA/2.0D0) THEN
173:             WRITE(1,'(I6)') J1198:             WRITE(1,'(I6)') J1
174: !           WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINM(J2), EMIN, ENERGY(MINM(J2))-HSHIFT, abs=',J1, J2, MINM(J2), EMIN(J1), &199:             WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINM(J2), EMIN, ENERGY(MINM(J2))-HSHIFT, abs=',J1, J2, MINM(J2), EMIN(J1), &
175: ! &                                     ENERGY(MINM(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINM(J2))+HSHIFT)200:   &                                     ENERGY(MINM(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINM(J2))+HSHIFT)
176:          ENDIF201:          ENDIF
177:       ENDDO202:       ENDDO
178:       IF (DUMMY.GT.COLOURTHRESH) EXIT203:       IF (DUMMY.GT.COLOURTHRESH) EXIT
179:    ENDDO204:    ENDDO
180:    CLOSE(1)205:    CLOSE(1)
181:    OPEN(UNIT=1,FILE='min.plus.hshift.BS',STATUS='UNKNOWN')206:    OPEN(UNIT=1,FILE='min.plus.hshift.BS',STATUS='UNKNOWN')
182:    DMIN=HUGE(1.0D0)207:    DMIN=HUGE(1.0D0)
183:    DO J1=1,NMIN208:    DO J1=1,NMIN
184:    ENDDO209:    ENDDO
185: 210: 
186:    DUMMY=0.0D0211:    DUMMY=0.0D0
187:    DO J2=1,NPLUS212:    DO J2=1,NPLUS
188:       DUMMY=DUMMY+PERMINP(J2)213:       DUMMY=DUMMY+PERMINP(J2)
189:       DO J1=1,NMIN214:       DO J1=1,NMIN
190:          IF (ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT).LT.DELTA/2.0D0) THEN215:          IF (ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT).LT.DELTA/2.0D0) THEN
191:             WRITE(1,'(I6)') J1216:             WRITE(1,'(I6)') J1
192: !           WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, ENERGY(MINP(J2))-HSHIFT, abs=',J1, J2, MINP(J2), EMIN(J1), &217:             WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, ENERGY(MINP(J2))-HSHIFT, abs=',J1, J2, MINP(J2), EMIN(J1), &
193: ! &                                     ENERGY(MINP(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT)218:   &                                     ENERGY(MINP(J2))-HSHIFT,  ABS(EMIN(J1)-ENERGY(MINP(J2))+HSHIFT)
194:          ENDIF219:          ENDIF
195:       ENDDO220:       ENDDO
196:       IF (DUMMY.GT.COLOURTHRESH) EXIT221:       IF (DUMMY.GT.COLOURTHRESH) EXIT
197:    ENDDO222:    ENDDO
198:    CLOSE(1)223:    CLOSE(1)
199: 224: 
200:    INQUIRE(FILE='weights.2D.BS',EXIST=DO2DT) 
201:    IF (DO2DT) THEN225:    IF (DO2DT) THEN
202:       PRINT '(A)','Reading 2D weights from weights.2D.BS'226:       Z0=0.0D0
203:       OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD')227:       Z1=0.0D0
204: !228:       DO J2=1,NBINS
205: ! This itme, look for NQBINS parameter at the end for allocations229:          DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT2(J2))
206: !230:          Z0=Z0+DUMMY
207:       READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, NBINS, MAXBIN, NQBINS231:          Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1))
208:       CLOSE(10) 
209:       ALLOCATE(MEANVQI(NQBINS),MEANQI(NQBINS),ANSHIFT(NQBINS),CLOSEST(NQBINS),QBIN(NMIN),QENERGY(NQBINS),QTOI(NQBINS)) 
210:       MEANVQI(1:NQBINS)=0.0D0 
211:       MEANQI(1:NQBINS)=0.0D0 
212:       OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD') 
213:       NCOUNT=0 
214:       ND1=HUGE(1) 
215:       DO  
216:          READ(19,*,END=753) J1, J2, DUMMY2, DUMMY, DWEIGHT 
217:          IF (J2.LT.ND1) THEN 
218:             NCOUNT=NCOUNT+1 
219:             QENERGY(NCOUNT)=DUMMY2 
220:          ENDIF 
221:          ND1=J2 
222:          IF (DUMMY-DUMMY2.GT.0.0D0) THEN 
223:             MEANVQI(NCOUNT)=MEANVQI(NCOUNT)+EXP(DWEIGHT-(DUMMY-ENERGY(1))/TEMPERATURE)*(DUMMY-DUMMY2) 
224:             MEANQI(NCOUNT)=MEANQI(NCOUNT)  +EXP(DWEIGHT-(DUMMY-ENERGY(1))/TEMPERATURE) 
225:          ENDIF 
226:       ENDDO232:       ENDDO
227: 753   CONTINUE233:       MEANE(1)=Z1/Z0
228: 234:       ZSAVE(1)=Z0
229:       NQBINSREAL=NCOUNT235:       NMINUS=0
230:       PRINT '(A,I8,A,I8)','Mean thermal PE evaluated for ',NQBINSREAL,' quanch bins, maximum bins was ',NQBINS236:       NPLUS=0
231:       CLOSE(19)237:       DPLUS=0.0D0
232:       PRINT '(A)','    q bin,    energy,    weight,    weighted thermal energy,    mean quench pe,       shift,  shift/harmonic'238:       DMINUS=0.0D0
233:       DO J1=1,NQBINSREAL239:       DO J2=1,NBINS
234:          ANSHIFT(J1)=MEANVQI(J1)/MAX(MEANQI(J1),1.0D-300)240:          DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT2(J2))
235:          PRINT '(I8,5G20.10)',J1,QENERGY(J1),MEANQI(J1),ANSHIFT(J1),ANSHIFT(J1)/HSHIFT241:          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
236:       ENDDO255:       ENDDO
237:    ENDIF256:       DO J2=1,NMINUS
238: 257:          PERMINM(J2)=PERMINM(J2)/DMINUS
239:    IF (DO2DT) THEN 
240:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.anhshift.BS min.plus.anhshift.BS and ' 
241:       PRINT '(A)','CHOOSECOLOURS in the dinfo file' 
242: ! 
243: ! which q bin is this minimum in?  
244: ! which i bin does this q bin correspond to in terms of mean thermal pe? 
245: ! is this one of the dominant probability flux i bins? 
246: ! 
247:       PRINT '(A)','Assignment of minima to quench bins: minimum, pe, q bin, q energy, diff' 
248:       DO J1=1,NMIN 
249:          DMIN=HUGE(1.0D0) 
250:          DO J2=1,NQBINSREAL 
251:             IF (ABS(EMIN(J1)-QENERGY(J2)).LT.DMIN) THEN 
252:                QBIN(J1)=J2 
253:                DMIN=ABS(EMIN(J1)-QENERGY(J2)) 
254:             ENDIF 
255:          ENDDO 
256:          PRINT '(I8,G20.10,I8,2G20.10)',J1,EMIN(J1),QBIN(J1),QENERGY(QBIN(J1)),DMIN 
257:       ENDDO258:       ENDDO
258:       PRINT '(A)','Assignment of i bins to q bins for calculated anharmonic thermal pe: q bin, q energy, i bin, i energy, diff'259:       DO J2=1,NPLUS
259:       DO J2=1,NQBINSREAL260:          PERMINP(J2)=PERMINP(J2)/DPLUS
260:          DMIN=HUGE(1.0D0) 
261:          DO J1=1,NBINS 
262:             IF (ABS(ENERGY(J1)-QENERGY(J2)-ANSHIFT(J2)).LT.DMIN) THEN 
263:                QTOI(J2)=J1 
264:                DMIN=ABS(ENERGY(J1)-QENERGY(J2)-ANSHIFT(J2)) 
265:             ENDIF 
266:          ENDDO 
267:          PRINT '(I8,G20.10,I8,2G20.10)',J2,QENERGY(J2),QTOI(J2),ENERGY(QTOI(J2)),DMIN 
268:       ENDDO261:       ENDDO
269: 262:       CALL SORT(NMINUS,NBINS,PERMINM,MINM)
270:       OPEN(UNIT=1,FILE='min.minus.anhshift.BS',STATUS='UNKNOWN')263:       CALL SORT(NPLUS,NBINS,PERMINP,MINP)
271:       DUMMY=0.0D0264:       DUMMY=0.0D0
 265:       PRINT '(A)','bins with negative dp/dT:'
272:       DO J2=1,NMINUS266:       DO J2=1,NMINUS
273:          DUMMY=DUMMY+PERMINM(J2)267:          DUMMY=DUMMY+PERMINM(J2)
274:          DO J1=1,NMIN268:          PRINT '(2I6,2G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY
275:             IF (QTOI(QBIN(J1)).EQ.MINM(J2)) THEN 
276:                WRITE(1,'(I6)') J1 
277:                WRITE(*,'(A,3I8,3G20.10)') 'J1, J2, MINM(J2), EMIN, diff=', & 
278:   &                                   J1, J2, MINM(J2), EMIN(J1), & 
279:   &                                   ABS(EMIN(J1)-ENERGY(MINM(J2))+ANSHIFT(QBIN(J1))) 
280:             ENDIF 
281:          ENDDO 
282:          IF (DUMMY.GT.COLOURTHRESH) EXIT 
283:       ENDDO269:       ENDDO
284:       CLOSE(1) 
285:  
286:       OPEN(UNIT=1,FILE='min.plus.anhshift.BS',STATUS='UNKNOWN') 
287:       DUMMY=0.0D0270:       DUMMY=0.0D0
 271:       PRINT '(A)','bins with positive dp/dT:'
288:       DO J2=1,NPLUS272:       DO J2=1,NPLUS
289:          DUMMY=DUMMY+PERMINP(J2)273:          DUMMY=DUMMY+PERMINP(J2)
290:          DO J1=1,NMIN274:          PRINT '(2I6,2G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY
291:             IF (QTOI(QBIN(J1)).EQ.MINP(J2)) THEN275:       ENDDO
292:                WRITE(1,'(I6)') J1276: 
293:                WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, diff=', &277:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.anhshift.BS min.plus.anhshift.BS and '
294:   &                                   J1, J2, MINP(J2), EMIN(J1), &278:       PRINT '(A)','CHOOSECOLOURS in the dinfo file'
295:   &                                   ABS(EMIN(J1)-ENERGY(MINP(J2))+ANSHIFT(QBIN(J1)))279:       OPEN(UNIT=1,FILE='min.minus.anhshift.BS',STATUS='UNKNOWN')
296:             ENDIF280:       DO J1=1,NMIN
297:          ENDDO281:            DMIN=HUGE(1.0D0)
298:          IF (DUMMY.GT.COLOURTHRESH) EXIT282:            DO J2=1,NBINS
299:       ENDDO 283:               IF (ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2)).LT.DMIN) THEN
300:       CLOSE(1)284:                  DMIN=ABS(EMIN(J1)-ENERGY(J2)+ANSHIFT(J2))
 285:                  CLOSEST(J1)=J2
 286: !                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))
 288:               ENDIF
 289: !             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))
 291:            ENDDO
 292:            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)))
 294:         ENDDO
 295: 
 296:         DUMMY=0.0D0
 297:         DO J2=1,NMINUS
 298:            DUMMY=DUMMY+PERMINM(J2)
 299:            DO J1=1,NMIN
 300:               IF (CLOSEST(J1).EQ.MINM(J2)) THEN
 301:                  WRITE(1,'(I6)') J1
 302:                  WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINM(J2), EMIN, ENERGY(MINM(J2))-ANSHIFT(MINM(J2)), abs=', &
 303:   &                                     J1, J2, MINM(J2), EMIN(J1), &
 304:   &                                     ENERGY(MINM(J2))-ANSHIFT(MINM(J2)), ABS(EMIN(J1)-ENERGY(MINM(J2))+ANSHIFT(MINM(J2)))
 305:               ENDIF
 306:            ENDDO
 307:            IF (DUMMY.GT.COLOURTHRESH) EXIT
 308:         ENDDO
 309:         CLOSE(1)
 310:         OPEN(UNIT=1,FILE='min.plus.anhshift.BS',STATUS='UNKNOWN')
 311: 
 312:         DUMMY=0.0D0
 313:         DO J2=1,NPLUS
 314:            DUMMY=DUMMY+PERMINP(J2)
 315:            DO J1=1,NMIN
 316:               IF (CLOSEST(J1).EQ.MINP(J2)) THEN
 317:                  WRITE(1,'(I6)') J1
 318:                  WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, ENERGY(MINP(J2))-ANSHIFT(MINP(J2)), abs=', &
 319:   &                                     J1, J2, MINP(J2), EMIN(J1), &
 320:   &                                     ENERGY(MINP(J2))-ANSHIFT(MINP(J2)), ABS(EMIN(J1)-ENERGY(MINP(J2))+ANSHIFT(MINP(J2)))
 321:               ENDIF
 322:            ENDDO
 323:            IF (DUMMY.GT.COLOURTHRESH) EXIT
 324:         ENDDO 
 325:         CLOSE(1)
301: 326: 
302:    ENDIF327:    ENDIF
303: ENDIF328: ENDIF
304: 329: 
305: END PROGRAM CVBS330: END PROGRAM CVBS
306: 331: 
307: SUBROUTINE SORT(N,J3,A,NA)332: SUBROUTINE SORT(N,J3,A,NA)
308: IMPLICIT NONE333: IMPLICIT NONE
309: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP334: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP
310: DOUBLE PRECISION TEMP, A(J3)335: DOUBLE PRECISION TEMP, A(J3)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0