hdiff output

r31545/Cv.HSA.f90 2016-11-29 17:30:06.912732748 +0000 r31544/Cv.HSA.f90 2016-11-29 17:30:07.220736833 +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 CALCCV  6: PROGRAM CALCCV
  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
  9: DOUBLE PRECISION COLOURTHRESH  9: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNFRQS(:), MEANE(:), ZSAVE(:)
 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNFRQS(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:) 10: INTEGER, ALLOCATABLE :: HORDER(:)
 11: INTEGER, ALLOCATABLE :: HORDER(:), MINP(:), MINM(:) 11: INTEGER J1, NTEMP, J2, KAPPA, NMIN
 12: INTEGER J1, NTEMP, J2, KAPPA, NMIN, NPLUS, NMINUS 
 13: LOGICAL YESNO 
 14:  12: 
 15: OPEN(UNIT=10,FILE='Cv.data',STATUS='OLD') 13: OPEN(UNIT=10,FILE='Cv.data',STATUS='OLD')
 16: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN 14: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN
 17: CLOSE(10) 15: CLOSE(10)
 18: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP)) 16: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP))
 19:  17: 
 20: PRINT '(3(A,I4))','minima=',NMIN 18: PRINT '(3(A,I4))','minima=',NMIN
 21: PRINT '(A,2G15.5,A,I8)','minimum and maximum T for Cv calculation=',TMIN,TMAX,' number of Cv data points=',NTEMP 19: 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 20: PRINT '(A,I8)','number of degrees of freedom=',KAPPA
 23:  21: 
 62:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+(LNFRQS(1)-LNFRQS(J2))/2.0D0)/HORDER(J2) 60:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+(LNFRQS(1)-LNFRQS(J2))/2.0D0)/HORDER(J2)
 63:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2) 61:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2)
 64:       PSUM=PSUM+DP*(ENERGY(J2)-ENERGY(1)-MEANE(J1)) 62:       PSUM=PSUM+DP*(ENERGY(J2)-ENERGY(1)-MEANE(J1))
 65: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP 63: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP
 66:       PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2 64:       PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2
 67:    ENDDO 65:    ENDDO
 68:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA+PSUM, KAPPA+PSUM2 66:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA+PSUM, KAPPA+PSUM2
 69: ENDDO 67: ENDDO
 70: CLOSE(1) 68: CLOSE(1)
 71:  69: 
 72: INQUIRE(FILE='Tanal',EXIST=YESNO) 
 73:  
 74: IF (YESNO) THEN 
 75:    NPLUS=0.0D0 
 76:    NMINUS=0.0D0 
 77:    DPLUS=0.0D0 
 78:    DMINUS=0.0D0 
 79:    ALLOCATE(PERMINP(NMIN),PERMINM(NMIN),MINP(NMIN),MINM(NMIN)) 
 80:    OPEN(UNIT=1,FILE='Tanal') 
 81:    READ(1,*) TANAL 
 82:    COLOURTHRESH=0.9D0 
 83:    READ(1,*,END=666) COLOURTHRESH 
 84: 666 CONTINUE 
 85:    PRINT '(2(A,G20.10))','Analysing heat capacity contributions at T=',TANAL,' colour threshold=',COLOURTHRESH 
 86:    CLOSE(1) 
 87:    TEMPERATURE=TANAL 
 88:    Z0=0.0D0 
 89:    Z1=0.0D0 
 90:    DO J2=1,NMIN 
 91:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+(LNFRQS(1)-LNFRQS(J2))/2.0D0)/HORDER(J2) 
 92:       Z0=Z0+DUMMY 
 93:       Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1)) 
 94:    ENDDO 
 95:    MEANE(1)=Z1/Z0 
 96:    ZSAVE(1)=Z0 
 97:    DO J2=1,NMIN 
 98:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+(LNFRQS(1)-LNFRQS(J2))/2.0D0)/HORDER(J2) 
 99:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(1))/(ZSAVE(1)*TEMPERATURE**2) 
100:       DUMMY=DP*(ENERGY(J2)-ENERGY(1)-MEANE(1)) 
101:       IF (DP.LT.0.0D0) THEN 
102:          NMINUS=NMINUS+1 
103:          MINM(NMINUS)=J2 
104:          PERMINM(NMINUS)=DUMMY 
105:          DMINUS=DMINUS+DUMMY 
106:       ELSE 
107:          NPLUS=NPLUS+1 
108:          MINP(NPLUS)=J2 
109:          PERMINP(NPLUS)=DUMMY 
110:          DPLUS=DPLUS+DUMMY 
111:       ENDIF 
112: !     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 
113:    ENDDO 
114:    DO J2=1,NMINUS 
115:       PERMINM(J2)=PERMINM(J2)/DMINUS 
116:    ENDDO 
117:    DO J2=1,NPLUS 
118:       PERMINP(J2)=PERMINP(J2)/DPLUS 
119:    ENDDO 
120:    CALL SORT(NMINUS,NMIN,PERMINM,MINM) 
121:    CALL SORT(NPLUS,NMIN,PERMINP,MINP) 
122:    DUMMY=0.0D0 
123:    PRINT '(A)','minima with negative dp/dT:' 
124:    DO J2=1,NMINUS 
125:       DUMMY=DUMMY+PERMINM(J2) 
126:       PRINT '(2I6,2G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY 
127:    ENDDO 
128:    DUMMY=0.0D0 
129:    PRINT '(A)','minima with positive dp/dT:' 
130:    DO J2=1,NPLUS 
131:       DUMMY=DUMMY+PERMINP(J2) 
132:       PRINT '(2I6,2G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY 
133:    ENDDO 
134:    PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus min.plus and ' 
135:    PRINT '(A)','CHOOSECOLOURS in the dinfo file' 
136:    OPEN(UNIT=1,FILE='min.minus',STATUS='UNKNOWN') 
137:    DUMMY=0.0D0 
138:    DO J2=1,NMINUS 
139:       DUMMY=DUMMY+PERMINM(J2) 
140:       WRITE(1,'(I6)') MINM(J2) 
141:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
142:    ENDDO 
143:    CLOSE(1) 
144:    OPEN(UNIT=1,FILE='min.plus',STATUS='UNKNOWN') 
145:    DUMMY=0.0D0 
146:    DO J2=1,NPLUS 
147:       DUMMY=DUMMY+PERMINP(J2) 
148:       WRITE(1,'(I6)') MINP(J2) 
149:       IF (DUMMY.GT.COLOURTHRESH) EXIT 
150:    ENDDO 
151:    CLOSE(1) 
152: ENDIF 
153:  
154: END PROGRAM CALCCV 70: END PROGRAM CALCCV
155:  71: 
156: SUBROUTINE SORT(N,J3,A,NA) 
157: IMPLICIT NONE 
158: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP 
159: DOUBLE PRECISION TEMP, A(J3) 
160: ! 
161: DO 20 J1=1,N-1 
162:    L=J1 
163:    DO 10 J2=J1+1,N 
164:       IF (A(L).LT.A(J2)) L=J2 
165: 10 CONTINUE 
166:    TEMP=A(L) 
167:    A(L)=A(J1) 
168:    A(J1)=TEMP 
169:    NTEMP=NA(L) 
170:    NA(L)=NA(J1) 
171:    NA(J1)=NTEMP 
172: 20 CONTINUE 
173: RETURN 
174: END SUBROUTINE SORT 
175:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0