hdiff output

r31574/Cv.BS.f90 2016-12-03 15:30:12.016926628 +0000 r31573/Cv.BS.f90 2016-12-03 15:30:12.256929705 +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.BS.f90' in revision 31573
  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 
 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:) 
 11: INTEGER, ALLOCATABLE :: MINP(:), MINM(:) 
 12: INTEGER J1, NTEMP, J2, KAPPA, NBINS, NPLUS, NMINUS, NDUMMY 
 13: LOGICAL YESNO 
 14:  
 15: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD') 
 16: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NBINS 
 17: CLOSE(10) 
 18: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP)) 
 19:  
 20: PRINT '(3(A,I10))','bins=',NBINS 
 21: PRINT '(A,2G15.5,A,I8)','minimum and maximum T for Cv calculation=',TMIN,TMAX,' number of Cv data points=',NTEMP 
 22: PRINT '(A,I8)','number of degrees of freedom=',KAPPA 
 23:  
 24: ALLOCATE(ENERGY(NBINS),LNWEIGHT(NBINS)) 
 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) 
 31: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD') 
 32: DO J1=1,NBINS 
 33:    READ(10,*) NDUMMY,DUMMY,ENERGY(J1),DUMMY,LNWEIGHT(J1) 
 34: ENDDO 
 35: CLOSE(10) 
 36:  
 37: TINT=(TMAX-TMIN)/(NTEMP-1) 
 38: DELTA=ENERGY(2)-ENERGY(1) 
 39: ! 
 40: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa, 
 41: !  which cancel. 
 42: ! 
 43: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN') 
 44: DO J1=1,NTEMP 
 45:    Z0=0.0D0 
 46:    Z1=0.0D0 
 47:    Z2=0.0D0 
 48:    TEMPERATURE=TMIN+(J1-1)*TINT 
 49:    DO J2=1,NBINS 
 50:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
 51:       Z0=Z0+DUMMY 
 52:       Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1)) 
 53:       Z2=Z2+DUMMY*(ENERGY(J2)-ENERGY(1))**2 
 54:    ENDDO 
 55:    MEANE(J1)=Z1/Z0 
 56:    ZSAVE(J1)=Z0 
 57:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN 
 58:       ONEMEXP=-DELTA/TEMPERATURE 
 59:    ELSE 
 60:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE) 
 61:    ENDIF 
 62:    WRITE(1,'(7G20.10)') TEMPERATURE, Z0, Z1, Z2, & 
 63:   &                     KAPPA - (Z1/(Z0*TEMPERATURE))**2 + Z2/(Z0*TEMPERATURE**2), MEANE(J1), ZSAVE(J1) 
 64: ENDDO 
 65: CLOSE(1) 
 66:  
 67: OPEN(UNIT=1,FILE='prob.out.BS',STATUS='UNKNOWN') 
 68: DO J1=1,NTEMP 
 69:    TEMPERATURE=TMIN+(J1-1)*TINT 
 70:    PSUM=0.0D0 
 71:    PSUM2=0.0D0 
 72:    DO J2=1,NBINS 
 73:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
 74:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2) 
 75:       PSUM=PSUM+DP*(ENERGY(J2)-ENERGY(1)-MEANE(J1)) 
 76: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP 
 77:       PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2 
 78:    ENDDO 
 79:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM, & 
 80:   &                                  KAPPA+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM2 
 81: ENDDO 
 82: CLOSE(1) 
 83:  
 84: INQUIRE(FILE='Tanal.BS',EXIST=YESNO) 
 85:  
 86: IF (YESNO) THEN 
 87:    NPLUS=0.0D0 
 88:    NMINUS=0.0D0 
 89:    DPLUS=0.0D0 
 90:    DMINUS=0.0D0 
 91:    ALLOCATE(PERMINP(NBINS),PERMINM(NBINS),MINP(NBINS),MINM(NBINS)) 
 92:    OPEN(UNIT=1,FILE='Tanal') 
 93:    READ(1,*) TANAL 
 94:    COLOURTHRESH=0.9D0 
 95:    READ(1,*,END=666) COLOURTHRESH 
 96: 666 CONTINUE 
 97:    PRINT '(2(A,G20.10))','Analysing heat capacity contributions at T=',TANAL,' colour threshold=',COLOURTHRESH 
 98:    CLOSE(1) 
 99:    TEMPERATURE=TANAL 
100:    Z0=0.0D0 
101:    Z1=0.0D0 
102:    DO J2=1,NBINS 
103:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
104:       Z0=Z0+DUMMY 
105:       Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1)) 
106:    ENDDO 
107:    MEANE(1)=Z1/Z0 
108:    ZSAVE(1)=Z0 
109:    DO J2=1,NBINS 
110:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2)) 
111:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(1))/(ZSAVE(1)*TEMPERATURE**2) 
112:       DUMMY=DP*(ENERGY(J2)-ENERGY(1)-MEANE(1)) 
113:       IF (DP.LT.0.0D0) THEN 
114:          NMINUS=NMINUS+1 
115:          MINM(NMINUS)=J2 
116:          PERMINM(NMINUS)=DUMMY 
117:          DMINUS=DMINUS+DUMMY 
118:       ELSE 
119:          NPLUS=NPLUS+1 
120:          MINP(NPLUS)=J2 
121:          PERMINP(NPLUS)=DUMMY 
122:          DPLUS=DPLUS+DUMMY 
123:       ENDIF 
124: !     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 
125:    ENDDO 
126:    DO J2=1,NMINUS 
127:       PERMINM(J2)=PERMINM(J2)/DMINUS 
128:    ENDDO 
129:    DO J2=1,NPLUS 
130:       PERMINP(J2)=PERMINP(J2)/DPLUS 
131:    ENDDO 
132:    CALL SORT(NMINUS,NBINS,PERMINM,MINM) 
133:    CALL SORT(NPLUS,NBINS,PERMINP,MINP) 
134:    DUMMY=0.0D0 
135:    PRINT '(A)','bins with negative dp/dT:' 
136:    DO J2=1,NMINUS 
137:       DUMMY=DUMMY+PERMINM(J2) 
138:       PRINT '(2I6,2G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY 
139:    ENDDO 
140:    DUMMY=0.0D0 
141:    PRINT '(A)','bins with positive dp/dT:' 
142:    DO J2=1,NPLUS 
143:       DUMMY=DUMMY+PERMINP(J2) 
144:       PRINT '(2I6,2G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY 
145:    ENDDO 
146:    PRINT *,'need to translate +/- bins list into minima indices and write min.minus.BS min.plus.BS' 
147:    PRINT *,'need to count the NMIN minima' 
148: !  PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.BS min.plus.BS and ' 
149:    PRINT '(A)','CHOOSECOLOURS in the dinfo file' 
150:    OPEN(UNIT=1,FILE='min.minus',STATUS='UNKNOWN') 
151:    DUMMY=0.0D0 
152:    DO J2=1,NMINUS 
153:       DUMMY=DUMMY+PERMINM(J2) 
154:       WRITE(1,'(I6)') MINM(J2) 
155:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
156:    ENDDO 
157:    CLOSE(1) 
158:    OPEN(UNIT=1,FILE='min.plus',STATUS='UNKNOWN') 
159:    DUMMY=0.0D0 
160:    DO J2=1,NPLUS 
161:       DUMMY=DUMMY+PERMINP(J2) 
162:       WRITE(1,'(I6)') MINP(J2) 
163:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
164:    ENDDO 
165:    CLOSE(1) 
166: ENDIF 
167:  
168: END PROGRAM CVBS 
169:  
170: SUBROUTINE SORT(N,J3,A,NA) 
171: IMPLICIT NONE 
172: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP 
173: DOUBLE PRECISION TEMP, A(J3) 
174: ! 
175: DO 20 J1=1,N-1 
176:    L=J1 
177:    DO 10 J2=J1+1,N 
178:       IF (A(L).LT.A(J2)) L=J2 
179: 10 CONTINUE 
180:    TEMP=A(L) 
181:    A(L)=A(J1) 
182:    A(J1)=TEMP 
183:    NTEMP=NA(L) 
184:    NA(L)=NA(J1) 
185:    NA(J1)=NTEMP 
186: 20 CONTINUE 
187: RETURN 
188: END SUBROUTINE SORT 
189:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0