hdiff output

r29826/efol.f90 2016-01-25 12:30:06.988782506 +0000 r29825/efol.f90 2016-01-25 12:30:08.808806567 +0000
397:          ZT(5)=.FALSE.397:          ZT(5)=.FALSE.
398:          ZT(6)=.FALSE.398:          ZT(6)=.FALSE.
399:          NZERO=6399:          NZERO=6
400:       ENDIF400:       ENDIF
401:       DO J1=1,NOPT-NEV401:       DO J1=1,NOPT-NEV
402:          ZT(J1)=.FALSE.402:          ZT(J1)=.FALSE.
403:       ENDDO403:       ENDDO
404: 404: 
405:       IF (VALUEST.AND.(MOD(ITER-1,NVALUES).EQ.0)) THEN405:       IF (VALUEST.AND.(MOD(ITER-1,NVALUES).EQ.0)) THEN
406:          IF (PTEST) WRITE(*,'(A)') ' Eigenvalues of the Hessian matrix:'406:          IF (PTEST) WRITE(*,'(A)') ' Eigenvalues of the Hessian matrix:'
407:          IF (PTEST) WRITE (*,'(6(G12.5,1X))') (DIAG(I),I=MAX(NZERO+1,NOPT-NEV+1),NOPT)407:          IF (PTEST) WRITE (*,'(6(F12.5,1X))') (DIAG(I),I=MAX(NZERO+1,NOPT-NEV+1),NOPT)
408:       ENDIF408:       ENDIF
409:       IF (VECTORST.AND.(MOD(ITER-1,NVECTORS).EQ.0)) THEN409:       IF (VECTORST.AND.(MOD(ITER-1,NVECTORS).EQ.0)) THEN
410:          IF (PTEST) WRITE(*,'(A)') ' Eigenvectors of the Hessian matrix:'410:          IF (PTEST) WRITE(*,'(A)') ' Eigenvectors of the Hessian matrix:'
411:          CALL HESSOUT(1)411:          CALL HESSOUT(1)
412:       ENDIF412:       ENDIF
413: !     PRINT '(A)','Predicted zero eigenvectors for Thomson'413: !     PRINT '(A)','Predicted zero eigenvectors for Thomson'
414: !     VECX(1:3*NATOMS)=0.0D0; VECY(1:3*NATOMS)=0.0D0; VECZ(1:3*NATOMS)=0.0D0414: !     VECX(1:3*NATOMS)=0.0D0; VECY(1:3*NATOMS)=0.0D0; VECZ(1:3*NATOMS)=0.0D0
415: !     DO J1=1,3*NATOMS,2415: !     DO J1=1,3*NATOMS,2
416: !        VECX(J1)=SIN(QTS(J1+1))416: !        VECX(J1)=SIN(QTS(J1+1))
417: !        VECY(J1)=COS(QTS(J1+1))417: !        VECY(J1)=COS(QTS(J1+1))


r29826/MKtrap.f 2016-01-25 12:30:06.604777429 +0000 r29825/MKtrap.f 2016-01-25 12:30:08.424801490 +0000
  1:       SUBROUTINE MKTRAP(N,X,VNEW,ENERGY,SSTEST)  1:       SUBROUTINE MKTRAP(N,X,VNEW,ENERGY,SSTEST)
  2:       USE MODHESS  2: C     SUBROUTINE TO COMPUTE THE POTENTIAL ENERGY OF THE INDIVIDUAL ATOMS.     
   3: C     (COMPARE NOTES EQS. (3) - (18)).
  3:       IMPLICIT NONE  4:       IMPLICIT NONE
  4:       INTEGER N,I,J,J1,J2,J3,J4,J5,J6  5:       INTEGER N,I,J,J1
  5:       DOUBLE PRECISION V0,DELTA,OMEGA,M,B,A,VW,C1,SUM,DIST,  6:       DOUBLE PRECISION V0,DELTA,OMEGA,M,B,A,VW,C1,SUM,SUM11,SUM12,SUM13,SUM21,SUM22,SUM23
  6:      & SUM2_11,SUM2_12,SUM2_13,SUM2_22,SUM2_23,SUM2_33,VNEW(3*N),ENERGY,SUM1_1,SUM1_2,SUM1_3  7:      &SUM1_3,SUM2_11,SUM2_12,SUM2_13,SUM2_22,SUM2_23,SUM2_33,VNEW(3*N),ENERGY,SUM1_1,SUM1_2,SUM1_3
  7:   8: !     PARAMETER(N=100)
  8:       PARAMETER(OMEGA=0.37D-6)  9:       PARAMETER(V0=0.33*10**(-6.0D0),DELTA=0.001575) ! original
  9:       PARAMETER(M=1.6D4) 10: !     PARAMETER(V0=0.33*10**(-6.0D0),DELTA=0.01575)  ! testing
 10:       PARAMETER(V0=0.33D-6) 11:       PARAMETER(OMEGA=0.37*10**(-6.0D0))
 11:       PARAMETER(DELTA=0.001575D0) ! weak wall 12:       PARAMETER(M=1.6*10**(4.0D0),B=0.0139*10**(-4.0D0))
 12: !     PARAMETER(DELTA=0.004725D0) ! strong wall 13:       PARAMETER(A=5.29*10**(-6.0D0))
 13:       PARAMETER(B=0.0139D-4) 
 14:       PARAMETER(A=5.29D-6) ! 5.29 micrometres 
 15:  
 16:       DOUBLE PRECISION X(3*N) 14:       DOUBLE PRECISION X(3*N)
 17:       DOUBLE PRECISION R(N,N),R3(N,N),R5(N,N),V(N) 15:       DOUBLE PRECISION R(N,N),R3(N,N),R5(N,N),V(N)
 18:       DOUBLE PRECISION V1(N,3) 16:       DOUBLE PRECISION V1(N,3)
 19:       DOUBLE PRECISION V2(N,9) 17:       DOUBLE PRECISION V2(N,9)
 20:       LOGICAL SSTEST 18:       LOGICAL SSTEST
 21:  19: 
 22:       VW=V0*DELTA 20:       VW=V0*DELTA
 23:       IF (SSTEST) HESS(1:3*N,1:3*N)=0.0D0 21: c     DELTA=0.001575 (WEAK WALL) AND 0.004725 (STRONG WALL).      
 24:  22: 
 25: C     N IS THE TOTAL NUMBER OF ATOMS. 23: C     N IS THE TOTAL NUMBER OF ATOMS.
 26: C     V0,VW,OMEGA,M,B ARE PARAMETERS CHARACTERIZING THE POTENTIAL. 24: C     V0,VW,OMEGA,M,B ARE PARAMETERS CHARACTERIZING THE POTENTIAL.
 27: C     A IS THE LENGTH UNIT. 25: C     A IS THE LENGTH UNIT.
 28: C     X(3*N) CONTAIN THE COORDINATES OF THE ATOMS IN THE 3 SPATIAL DIRECTIONS X,Y,Z. 26: C     X(3*N) CONTAIN THE COORDINATES OF THE ATOMS IN THE 3 SPATIAL DIRECTIONS X,Y,Z.
 29: C     R(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)  27: C     R(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS) 
 30: C     R3(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**3  28: C     R3(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**3 
 31: C     R5(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**5  29: C     R5(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**5 
 32: C     V(I) IS THE POTENTIAL ENERGY OF THE ITH ATOM. 30: C     V(I) IS THE POTENTIAL ENERGY OF THE ITH ATOM.
 33: C     V1(I,1) IS \PARTIAL_X V(I). 31: C     V1(I,1) IS \PARTIAL_X V(I).
 38: C     V2(I,3) IS \PARTIAL_XZ V(I). 36: C     V2(I,3) IS \PARTIAL_XZ V(I).
 39: C     V2(I,4) IS \PARTIAL_YX V(I). 37: C     V2(I,4) IS \PARTIAL_YX V(I).
 40: C     V2(I,5) IS \PARTIAL_YY V(I). 38: C     V2(I,5) IS \PARTIAL_YY V(I).
 41: C     V2(I,6) IS \PARTIAL_YZ V(I). 39: C     V2(I,6) IS \PARTIAL_YZ V(I).
 42: C     V2(I,7) IS \PARTIAL_ZX V(I). 40: C     V2(I,7) IS \PARTIAL_ZX V(I).
 43: C     V2(I,8) IS \PARTIAL_ZY V(I). 41: C     V2(I,8) IS \PARTIAL_ZY V(I).
 44: C     V2(I,9) IS \PARTIAL_ZZ V(I). 42: C     V2(I,9) IS \PARTIAL_ZZ V(I).
 45:  43: 
 46:       C1=0.5*(OMEGA-M*OMEGA**2-V0) 44:       C1=0.5*(OMEGA-M*OMEGA**2-V0)
 47:  45: 
  46: C     GIVEN COORDINATES OF THE ATOMS.
  47: !     DO I=1,3*N
  48: !      READ(*,*)X(I)
  49: !     ENDDO      
  50: 
 48: C     COMPUTING THE DISTANCE MATRIX. 51: C     COMPUTING THE DISTANCE MATRIX.
 49:       CALL RMATMK(N,X,R,R3,R5) 52:       CALL RMATMK(N,X,R,R3,R5)
 50:  53: 
 51:       ENERGY=0.0D0 54:       ENERGY=0.0D0
 52:       DO I=1,N 55:       DO I=1,N
 53:  56: 
 54: c     THE POTENTIAL.  57: c     THE POTENTIAL. 
 55:        SUM=0.0D0 58:        SUM=0.0D0
 56:        DO J=I+1,I+N-1 59:        DO J=I+1,I+N-1
 57:         IF(J.GT.N)THEN 60:         IF(J.GT.N)THEN
 59:         ELSE 62:         ELSE
 60:          J1=J 63:          J1=J
 61:         ENDIF 64:         ENDIF
 62:         SUM=SUM+R(I,J1) 65:         SUM=SUM+R(I,J1)
 63:        ENDDO  66:        ENDDO 
 64:        V(I)=V0*X(3*(I-1)+3)**2 67:        V(I)=V0*X(3*(I-1)+3)**2
 65:      & +C1*(X(3*(I-1)+1)**2+X(3*(I-1)+2)**2) 68:      & +C1*(X(3*(I-1)+1)**2+X(3*(I-1)+2)**2)
 66:      & +VW*(X(3*(I-1)+1)**2-X(3*(I-1)+2)**2) 69:      & +VW*(X(3*(I-1)+1)**2-X(3*(I-1)+2)**2)
 67:      & +B*SUM 70:      & +B*SUM
 68:         ENERGY=ENERGY+V(I) 71:         ENERGY=ENERGY+V(I)
  72: !     PRINT *,'mktrap> I,SUM,V(I)=',I,SUM,V(I)
  73: !     PRINT *,'mktrap> I,energy=',I,ENERGY
 69:  74: 
 70: c     THE FIRST DERIVATIVES OF THE POTENTIAL.  75: c     THE FIRST DERIVATIVES OF THE POTENTIAL. 
 71:        SUM1_1=0.0D0 76:        SUM1_1=0.0D0
 72:        SUM1_2=0.0D0 77:        SUM1_2=0.0D0
 73:        SUM1_3=0.0D0 78:        SUM1_3=0.0D0
 74:        DO J=I+1,I+N-1 79:        DO J=I+1,I+N-1
 75:         IF(J.GT.N)THEN 80:         IF(J.GT.N)THEN
 76:          J1=MOD(J,N) 81:          J1=MOD(J,N)
 77:         ELSE 82:         ELSE
 78:          J1=J 83:          J1=J
 83:        ENDDO  88:        ENDDO 
 84:        V1(I,1)=2*(C1+VW)*X(3*(I-1)+1)-B*SUM1_1 89:        V1(I,1)=2*(C1+VW)*X(3*(I-1)+1)-B*SUM1_1
 85:        V1(I,2)=2*(C1-VW)*X(3*(I-1)+2)-B*SUM1_2 90:        V1(I,2)=2*(C1-VW)*X(3*(I-1)+2)-B*SUM1_2
 86:        V1(I,3)=2*V0*X(3*(I-1)+3)-B*SUM1_3 91:        V1(I,3)=2*V0*X(3*(I-1)+3)-B*SUM1_3
 87:        VNEW(3*(I-1)+1)=V1(I,1) 92:        VNEW(3*(I-1)+1)=V1(I,1)
 88:        VNEW(3*(I-1)+2)=V1(I,2) 93:        VNEW(3*(I-1)+2)=V1(I,2)
 89:        VNEW(3*(I-1)+3)=V1(I,3) 94:        VNEW(3*(I-1)+3)=V1(I,3)
 90:  95: 
 91: c     THE SECOND DERIVATIVES OF THE POTENTIAL. 96: c     THE SECOND DERIVATIVES OF THE POTENTIAL.
 92:       IF (SSTEST) THEN 97:       IF (SSTEST) THEN
  98:           PRINT *,'MKtrap> need to set Hessian!'
  99:           STOP
 93:        SUM=0.D0100:        SUM=0.D0
 94:        SUM2_11=0.0D0101:        SUM2_11=0.0D0
 95:        SUM2_12=0.0D0102:        SUM2_12=0.0D0
 96:        SUM2_13=0.0D0103:        SUM2_13=0.0D0
 97:        SUM2_22=0.0D0104:        SUM2_22=0.0D0
 98:        SUM2_23=0.0D0105:        SUM2_23=0.0D0
 99:        SUM2_33=0.0D0106:        SUM2_33=0.0D0
100:        DO J=I+1,I+N-1107:        DO J=I+1,I+N-1
101:         IF(J.GT.N)THEN108:         IF(J.GT.N)THEN
102:          J1=MOD(J,N)109:          J1=MOD(J,N)
113:      &  *R5(I,J1)120:      &  *R5(I,J1)
114:         SUM2_22=SUM2_22+((X(3*(I-1)+2)-X(3*(J1-1)+2))**2)*R5(I,J1)121:         SUM2_22=SUM2_22+((X(3*(I-1)+2)-X(3*(J1-1)+2))**2)*R5(I,J1)
115:         SUM2_23=SUM2_23+122:         SUM2_23=SUM2_23+
116:      &  ((X(3*(I-1)+2)-X(3*(J1-1)+2))*(X(3*(I-1)+3)-X(3*(J1-1)+3)))123:      &  ((X(3*(I-1)+2)-X(3*(J1-1)+2))*(X(3*(I-1)+3)-X(3*(J1-1)+3)))
117:      &  *R5(I,J1)124:      &  *R5(I,J1)
118:         SUM2_33=SUM2_33+((X(3*(I-1)+3)-X(3*(J1-1)+3))**2)*R5(I,J1)125:         SUM2_33=SUM2_33+((X(3*(I-1)+3)-X(3*(J1-1)+3))**2)*R5(I,J1)
119:        ENDDO 126:        ENDDO 
120:        V2(I,1)=2*(C1+VW)-B*SUM+3*B*SUM2_11127:        V2(I,1)=2*(C1+VW)-B*SUM+3*B*SUM2_11
121:        V2(I,2)=3*B*SUM2_12128:        V2(I,2)=3*B*SUM2_12
122:        V2(I,3)=3*B*SUM2_13129:        V2(I,3)=3*B*SUM2_13
 130:        V2(I,4)=V2(I,2)
123:        V2(I,5)=2*(C1-VW)-B*SUM+3*B*SUM2_22131:        V2(I,5)=2*(C1-VW)-B*SUM+3*B*SUM2_22
 132:        V2(I,6)=3*B*SUM2_23
 133:        V2(I,7)=V2(I,3)
124:        V2(I,8)=V2(I,6)134:        V2(I,8)=V2(I,6)
125:        V2(I,9)=2*V0-B*SUM+3*B*SUM2_33135:        V2(I,9)=2*V0-B*SUM+3*B*SUM2_33
126:        HESS(3*(I-1)+1,3*(I-1)+1)=V2(I,1) 
127:        HESS(3*(I-1)+1,3*(I-1)+2)=V2(I,2) 
128:        HESS(3*(I-1)+1,3*(I-1)+3)=V2(I,3) 
129:        HESS(3*(I-1)+2,3*(I-1)+1)=V2(I,2) 
130:        HESS(3*(I-1)+2,3*(I-1)+2)=V2(I,5) 
131:        HESS(3*(I-1)+2,3*(I-1)+3)=V2(I,6) 
132:        HESS(3*(I-1)+3,3*(I-1)+1)=V2(I,3) 
133:        HESS(3*(I-1)+3,3*(I-1)+2)=V2(I,6) 
134:        HESS(3*(I-1)+3,3*(I-1)+3)=V2(I,9) 
135:       ENDIF136:       ENDIF
136:       ENDDO 137:       ENDDO 
137: 138: 
138:       ENERGY=ENERGY/2.0D0 ! energy was out by a factor of two DJW 25/1/16 
139:  
140:       IF (.NOT.SSTEST) RETURN 
141:  
142: ! 
143: !  Include missing second derivatives. DJW 25/1/16 
144: ! 
145:  
146: C 
147: C  Case III, different atoms, same cartesian coordinate. 
148: C 
149:       DO J1=1,N 
150:          DO J2=1,3 
151:             J3=3*(J1-1)+J2 
152:             DO J4=J1+1,N 
153:                J5=3*(J4-1)+J2 
154:                HESS(J3,3*(J4-1)+J2)=-B*3.0D0*R5(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2+B*R3(J4,J1) 
155:                HESS(3*(J4-1)+J2,J3)=HESS(J3,3*(J4-1)+J2) 
156:            ENDDO 
157:          ENDDO 
158:       ENDDO 
159: C 
160: C  Case IV: different atoms and different cartesian coordinates. 
161: C 
162:       DO J1=1,N 
163:          DO J2=1,3 
164:             J3=3*(J1-1)+J2 
165:             DO J4=J1+1,N 
166:                DO J5=J2+1,3 
167:                   J6=3*(J4-1)+J5 
168:                   HESS(J3,J6)=-B*R5(J4,J1)*3.0D0 
169:      1                    *(X(J3)-X(3*(J4-1)+J2)) 
170:      2                    *(X(3*(J1-1)+J5)-X(J6)) 
171:                   HESS(J6,J3)=HESS(J3,J6) 
172:                   HESS(3*(J4-1)+J2,3*(J1-1)+J5)=HESS(J3,J6) 
173:                   HESS(3*(J1-1)+J5,3*(J4-1)+J2)=HESS(J3,J6) 
174:                ENDDO 
175:             ENDDO 
176:          ENDDO 
177:       ENDDO 
178:  
179:       RETURN139:       RETURN
180:       END140:       END
181:       141:       
182:       SUBROUTINE RMATMK(N,X,R,R3,R5)142:       SUBROUTINE RMATMK(N,X,R,R3,R5)
183:       IMPLICIT NONE143:       IMPLICIT NONE
184:       INTEGER N,J1,J2144:       INTEGER N,J1,J2
185:       DOUBLE PRECISION X(3*N),DISTSQ145:       DOUBLE PRECISION X(3*N),DISTSQ
186:       DOUBLE PRECISION R(N,N),R3(N,N),R5(N,N) 146:       DOUBLE PRECISION R(N,N),R3(N,N),R5(N,N) 
187:       DO J1=1,N147:       DO 10 J1=1,N
188:          R(J1,J1)=0.0D0148:          R(J1,J1)=0.0D0
189:          DO J2=J1+1,N149:          DO 20 J2=J1+1,N
190:           DISTSQ=(X(3*(J1-1)+1)-X(3*(J2-1)+1))**2150:           DISTSQ=(X(3*(J1-1)+1)-X(3*(J2-1)+1))**2
191:      &        +(X(3*(J1-1)+2)-X(3*(J2-1)+2))**2151:      &        +(X(3*(J1-1)+2)-X(3*(J2-1)+2))**2
192:      &        +(X(3*(J1-1)+3)-X(3*(J2-1)+3))**2152:      &        +(X(3*(J1-1)+3)-X(3*(J2-1)+3))**2
193:           R(J1,J2)=1.0D0/DISTSQ**(0.5D0)153:           R(J1,J2)=1.0D0/DISTSQ**(0.5D0)
194:           R3(J1,J2)=1.0D0/DISTSQ**(1.5D0)154:           R3(J1,J2)=1.0D0/DISTSQ**(1.5D0)
195:           R5(J1,J2)=1.0D0/DISTSQ**(2.5D0)155:           R5(J1,J2)=1.0D0/DISTSQ**(2.5D0)
196:           R(J2,J1)=R(J1,J2)156:           R(J2,J1)=R(J1,J2)
197:           R3(J2,J1)=R3(J1,J2)157:           R3(J2,J1)=R3(J1,J2)
198:           R5(J2,J1)=R5(J1,J2)158:           R5(J2,J1)=R5(J1,J2)
199:          ENDDO159: 20       CONTINUE
200:       ENDDO160: 10    CONTINUE
201:       RETURN161:       RETURN
202:       END162:       END
203: 163: 


r29826/OPTIM.F 2016-01-25 12:30:06.796779969 +0000 r29825/OPTIM.F 2016-01-25 12:30:08.616804028 +0000
413:          ELSE IF (.NOT.BULKT) THEN413:          ELSE IF (.NOT.BULKT) THEN
414:             SHIFTL(1)=SHIFTV414:             SHIFTL(1)=SHIFTV
415:             SHIFTL(2)=SHIFTV415:             SHIFTL(2)=SHIFTV
416:             SHIFTL(6)=SHIFTV416:             SHIFTL(6)=SHIFTV
417:             WRITE(*,'(A)') ' OPTIM> x,y translational and z rotational ev shifting'417:             WRITE(*,'(A)') ' OPTIM> x,y translational and z rotational ev shifting'
418:          ELSE418:          ELSE
419:             SHIFTL(1)=SHIFTV419:             SHIFTL(1)=SHIFTV
420:             SHIFTL(2)=SHIFTV420:             SHIFTL(2)=SHIFTV
421:             WRITE(*,'(A)') ' OPTIM> x,y translational ev shifting'421:             WRITE(*,'(A)') ' OPTIM> x,y translational ev shifting'
422:          ENDIF422:          ENDIF
423:       ELSE IF (MKTRAPT) THEN 
424:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted' 
425:       ELSE IF (EYTRAPT.OR.(ZSYMSAVE.EQ.'BE')) THEN423:       ELSE IF (EYTRAPT.OR.(ZSYMSAVE.EQ.'BE')) THEN
426:          SHIFTL(4)=SHIFTV424:          SHIFTL(4)=SHIFTV
427:          SHIFTL(5)=SHIFTV425:          SHIFTL(5)=SHIFTV
428:          SHIFTL(6)=SHIFTV426:          SHIFTL(6)=SHIFTV
429:          WRITE(*,'(A,G20.10)') ' OPTIM> Using rotational ev shift=',SHIFTV427:          WRITE(*,'(A,G20.10)') ' OPTIM> Using rotational ev shift=',SHIFTV
430:       ELSE IF (BULKT) THEN428:       ELSE IF (BULKT) THEN
431:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational ev shift=',SHIFTV429:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational ev shift=',SHIFTV
432:          DO J1=1,3430:          DO J1=1,3
433:             SHIFTL(J1)=SHIFTV431:             SHIFTL(J1)=SHIFTV
434:          ENDDO432:          ENDDO
445:             WRITE(*,'(A)') ' OPTIM> z rotational ev shifting'443:             WRITE(*,'(A)') ' OPTIM> z rotational ev shifting'
446:          ENDIF444:          ENDIF
447:          SHIFTL(1)=SHIFTV445:          SHIFTL(1)=SHIFTV
448:          SHIFTL(2)=SHIFTV446:          SHIFTL(2)=SHIFTV
449:          SHIFTL(3)=SHIFTV447:          SHIFTL(3)=SHIFTV
450:          SHIFTL(4)=SHIFTV448:          SHIFTL(4)=SHIFTV
451:          SHIFTL(5)=SHIFTV449:          SHIFTL(5)=SHIFTV
452:          SHIFTL(6)=SHIFTV450:          SHIFTL(6)=SHIFTV
453:       ELSE IF (MIEFT) THEN451:       ELSE IF (MIEFT) THEN
454:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted'452:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted'
 453:       ELSE IF (MKTRAPT) THEN
 454:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted'
455:       ELSE455:       ELSE
456:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV456:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV
457:          DO J1=1,6457:          DO J1=1,6
458:             SHIFTL(J1)=SHIFTV458:             SHIFTL(J1)=SHIFTV
459:          ENDDO459:          ENDDO
460:       ENDIF460:       ENDIF
461: 461: 
462:       !js850> 462:       !js850> 
463:       IF (ZSYM(NATOMS).EQ.'SS') THEN463:       IF (ZSYM(NATOMS).EQ.'SS') THEN
464:          ALLOCATE(SOFT_SPHERE_RADII(NATOMS))464:          ALLOCATE(SOFT_SPHERE_RADII(NATOMS))


r29826/potential.f 2016-01-25 12:30:08.232801264 +0000 r29825/potential.f 2016-01-25 12:30:09.004809158 +0000
1348:             ENDIF1348:             ENDIF
1349:          ELSE IF (EYTRAPT) THEN1349:          ELSE IF (EYTRAPT) THEN
1350:             CALL EYETRAP(NATOMS,COORDS,VNEW,ENERGY,C1,C2,C3)1350:             CALL EYETRAP(NATOMS,COORDS,VNEW,ENERGY,C1,C2,C3)
1351:             IF (SSTEST) CALL EYDTRAP(NATOMS,COORDS,VNEW,C1,C2,C3,GTEST,SSTEST)1351:             IF (SSTEST) CALL EYDTRAP(NATOMS,COORDS,VNEW,C1,C2,C3,GTEST,SSTEST)
1352:             IF (PTEST) THEN1352:             IF (PTEST) THEN
1353:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' hartree'1353:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' hartree'
1354:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' hartree'1354:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' hartree'
1355:             ENDIF1355:             ENDIF
1356:          ELSE IF (MKTRAPT) THEN1356:          ELSE IF (MKTRAPT) THEN
1357:             CALL MKTRAP(NATOMS,COORDS,VNEW,ENERGY,SSTEST)1357:             CALL MKTRAP(NATOMS,COORDS,VNEW,ENERGY,SSTEST)
1358:  
1359: !              DIFF=1.0D-4 
1360: !              PRINT*,'analytic and numerical gradients:' 
1361: !              IF (.NOT.(ALLOCATED(HESS))) ALLOCATE(HESS(NOPT,NOPT)) 
1362: !              CALL MKTRAP(NATOMS, COORDS, VNEW, ENERGY,.TRUE.) 
1363: !              PRINT '(A,I8)','SIZE(HESS)=',SIZE(HESS) 
1364: !              HESSDUM(1:NOPT,1:NOPT)=HESS(1:NOPT,1:NOPT) 
1365: !              DO J1=1,NATOMS 
1366: !                 COORDS(J1)=COORDS(J1)+DIFF 
1367: !                 CALL MKTRAP(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.) 
1368: !                 COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
1369: !                 CALL MKTRAP(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.) 
1370: !                 COORDS(J1)=COORDS(J1)+DIFF 
1371: !                 IF ((ABS(VNEW(J1)).NE.0.0D0).AND.ABS(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.-1.0D0) THEN 
1372: !                    WRITE(*,'(I5,2F20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF) 
1373: !                 ENDIF 
1374: !              ENDDO 
1375: !              PRINT*,'analytic and numerical second derivatives:' 
1376: !              DO J1=1,NATOMS 
1377: !                 COORDS(J1)=COORDS(J1)+DIFF 
1378: !                 CALL MKTRAP(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.) 
1379: !                 COORDS(J1)=COORDS(J1)-2.0D0*DIFF 
1380: !                 CALL MKTRAP(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.) 
1381: !                 COORDS(J1)=COORDS(J1)+DIFF 
1382: !                 DO J2=1,NATOMS 
1383: !                    IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND. 
1384: !    &                   (ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D-3)) THEN 
1385: !                    WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X' 
1386: !                    ELSE 
1387: !                       WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF) 
1388: !                    ENDIF 
1389: !                 ENDDO 
1390: !              ENDDO 
1391: !              STOP 
1392:  
1393:             IF (PTEST) THEN1358:             IF (PTEST) THEN
1394:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY ! ,' hartree'1359:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY ! ,' hartree'
1395:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY ! ,' hartree'1360:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY ! ,' hartree'
1396:             ENDIF1361:             ENDIF
1397:          ELSE IF (ZSYM(NATOMS).EQ.'C1') THEN1362:          ELSE IF (ZSYM(NATOMS).EQ.'C1') THEN
1398:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'1363:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'
1399:             CALL C10(COORDS,NATOMS,VNEW,ENERGY,GTEST,SSTEST)1364:             CALL C10(COORDS,NATOMS,VNEW,ENERGY,GTEST,SSTEST)
1400:             IF (PTEST) THEN1365:             IF (PTEST) THEN
1401:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' arbs'1366:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' arbs'
1402:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' arbs'1367:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' arbs'


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0