hdiff output

r31425/grad.f90 2016-11-02 17:30:14.126302221 +0000 r31424/grad.f90 2016-11-02 17:30:14.866312197 +0000
 25:      ! http://www-wales.ch.cam.ac.uk/~sat39/DNEBtests/  25:      ! http://www-wales.ch.cam.ac.uk/~sat39/DNEBtests/ 
 26:      SUBROUTINE NEBGRADIENT 26:      SUBROUTINE NEBGRADIENT
 27:           USE KEYGRAD 27:           USE KEYGRAD
 28:           USE KEYTAU 28:           USE KEYTAU
 29:           USE NEBDATA 29:           USE NEBDATA
 30:           USE KEYNEB,ONLY: NIMAGE 30:           USE KEYNEB,ONLY: NIMAGE
 31:           USE TANGENT 31:           USE TANGENT
 32:           USE SPFUNCTS 32:           USE SPFUNCTS
 33:           USE NEBUTILS 33:           USE NEBUTILS
 34:           USE MODUNRES, ONLY: C,NRES 34:           USE MODUNRES, ONLY: C,NRES
 35:           USE GENRIGID, ONLY: RIGIDINIT, DEGFREEDOMS, RB_DISTANCE 
 36:           USE KEY, ONLY: UNRST, FROZEN, FREEZE, NEBK, BULKT, NEBRESEEDT, NEBRESEEDEMAX, NEBRESEEDBMAX, & 35:           USE KEY, ONLY: UNRST, FROZEN, FREEZE, NEBK, BULKT, NEBRESEEDT, NEBRESEEDEMAX, NEBRESEEDBMAX, &
 37:   &                      RBAAT, NREPULSIVE, NEBRESEEDINT, INTCONFRAC, INTCONSTRAINTT, & 36:   &                      RBAAT, NREPULSIVE, NEBRESEEDINT, INTCONFRAC, INTCONSTRAINTT, &
 38:   &                      ORDERI,ORDERJ,EPSALPHA, REPPOW,PERMDIST, TWOD, RIGIDBODY, NEBRESEEDDEL1, & 37:   &                      ORDERI,ORDERJ,EPSALPHA, REPPOW,PERMDIST, TWOD, RIGIDBODY, NEBRESEEDDEL1, &
 39:   &                      DISTREF, ADDREPT, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, CHECKCONINT 38:   &                      DISTREF, ADDREPT, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, CHECKCONINT
 40:           USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3, DEBUG 39:           USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3, DEBUG
 41:           IMPLICIT NONE 40:           IMPLICIT NONE
 42:             41:            
 43:           INTEGER :: J1,J2,K,J3,J4,J5,JOUT,JIN,NSHRINK,NSHRINKATOM(NATOMS),J6,NATTACH,NATTACHS,NMAXINT,NMININT,LUNIT 42:           INTEGER :: J1,J2,K,J3,J4,J5,JOUT,JIN,NSHRINK,NSHRINKATOM(NATOMS),J6,NATTACH,NATTACHS,NMAXINT,NMININT,LUNIT
 44:           DOUBLE PRECISION :: PERPCOMP=0.1D0 43:           DOUBLE PRECISION :: PERPCOMP=0.1D0
 45:           DOUBLE PRECISION :: DIHEDIST,TMPRMS,QINT(NINTS*(NIMAGE+2)),DIFFM(NINTS),DIFFP(NINTS) ! JMC 44:           DOUBLE PRECISION :: DIHEDIST,TMPRMS,QINT(NINTS*(NIMAGE+2)),DIFFM(NINTS),DIFFP(NINTS) ! JMC
 46:           DOUBLE PRECISION DIST 45:           DOUBLE PRECISION DIST
 47:           DOUBLE PRECISION :: PI=3.141592653589793D0 46:           DOUBLE PRECISION :: PI=3.141592653589793D0
 48:           DOUBLE PRECISION TEMP1(NOPT), TEMP2(NOPT), GGGSAVE(NOPT*(NIMAGE+2)), DUMMY, DUMMY2, SPRING(NOPT*(NIMAGE+2)), & 47:           DOUBLE PRECISION TEMP1(NOPT), TEMP2(NOPT), GGGSAVE(NOPT*(NIMAGE+2)), DUMMY, DUMMY2, SPRING(NOPT*(NIMAGE+2)), &
 49:   &                        TRUEPERP(NOPT*(NIMAGE+2)), IMCOORDS(NOPT), RSITE(NOPT), RMAT(3,3), D, DIST2, LX(3), LV(3), & 48:   &                        TRUEPERP(NOPT*(NIMAGE+2)), IMCOORDS(NOPT), RSITE(NOPT), RMAT(3,3), D, DIST2, LX(3), LV(3), &
 50:   &                        REPGRAD(NOPT), MAGREP, MEAND, VS(3), VF(3), MP1S(3), MP2S(3), MP1F(3), MP2F(3), DS, DF, D1, D2, DI, DJ, & 49:   &                        REPGRAD(NOPT), MAGREP, MEAND, VS(3), VF(3), MP1S(3), MP2S(3), MP1F(3), MP2F(3), DS, DF, D1, D2, DI, DJ, &
 51:   &                        MP2SS(3), MP2FF(3), MP1SS(3), MP1FF(3), DSS, DFF, DISTA, DISTB 50:   &                        MP2SS(3), MP2FF(3), MP1SS(3), MP1FF(3), DSS, DFF, DISTA, DISTB
 52:  
 53:           ! sn402: rigid body variables 
 54:           DOUBLE PRECISION :: DISTL, DISTR, GRADL1(DEGFREEDOMS), GRADR1(DEGFREEDOMS), GRADL2(DEGFREEDOMS), GRADR2(DEGFREEDOMS) 
 55:  
 56:           INTEGER, ALLOCATABLE :: IREPTEMP(:) 51:           INTEGER, ALLOCATABLE :: IREPTEMP(:)
 57:           LOGICAL SHRINKT(NATOMS) 52:           LOGICAL SHRINKT(NATOMS)
 58:           DOUBLE PRECISION EEETMP(NIMAGE+2), MYGTMP(NOPT*(NIMAGE+2)), ETMP, RMSTMP, GETUNIT 53:           DOUBLE PRECISION EEETMP(NIMAGE+2), MYGTMP(NOPT*(NIMAGE+2)), ETMP, RMSTMP, GETUNIT
 59:  
 60:           GSPR(:) = 0.0D0  ! sn402: added 
 61:           GGG(:) = 0.0D0   ! sn402: added 
 62:           SSS(:) = 0.0D0   ! sn402: added 
 63:            54:           
 64:           CALL TRUEPOTEG 55:           CALL TRUEPOTEG
 65:  56: 
 66: !         IF (INTCONSTRAINTT.AND.(INTCONFRAC.LT.1.0D0)) THEN 57: !         IF (INTCONSTRAINTT.AND.(INTCONFRAC.LT.1.0D0)) THEN
 67: !            EEETMP(1:NIMAGE+2)=EEE(1:NIMAGE+2) 58: !            EEETMP(1:NIMAGE+2)=EEE(1:NIMAGE+2)
 68: !            MYGTMP(1:NOPT*(NIMAGE+2))=GGG(1:NOPT*(NIMAGE+2)) 59: !            MYGTMP(1:NOPT*(NIMAGE+2))=GGG(1:NOPT*(NIMAGE+2))
 69: !            IF (CHECKCONINT) THEN 60: !            IF (CHECKCONINT) THEN
 70: !               CALL CONGRAD2(NMAXINT,NMININT,ETMP,XYZ,GGG,EEE,IMGFREEZE,RMSTMP) 61: !               CALL CONGRAD2(NMAXINT,NMININT,ETMP,XYZ,GGG,EEE,IMGFREEZE,RMSTMP)
 71: !            ELSE 62: !            ELSE
 72: !               CALL CONGRAD(NMAXINT,NMININT,ETMP,XYZ,GGG,EEE,IMGFREEZE,RMSTMP) 63: !               CALL CONGRAD(NMAXINT,NMININT,ETMP,XYZ,GGG,EEE,IMGFREEZE,RMSTMP)
101: !         ENDIF 92: !         ENDIF
102: !           93: !          
103: ! GGGSAVE contains the true gradient on images 94: ! GGGSAVE contains the true gradient on images
104: ! 95: !
105:           GGGSAVE(1:NOPT*(NIMAGE+2))=GGG(1:NOPT*(NIMAGE+2)) 96:           GGGSAVE(1:NOPT*(NIMAGE+2))=GGG(1:NOPT*(NIMAGE+2))
106:           ! GGG contains the perpendicular component of the true gradient after this block 97:           ! GGG contains the perpendicular component of the true gradient after this block
107:           DO J1=2,NIMAGE+1  !  MEP-perpendicular component of true gradient only: G = G - G|| = G - (G,TAU)*TAU 98:           DO J1=2,NIMAGE+1  !  MEP-perpendicular component of true gradient only: G = G - G|| = G - (G,TAU)*TAU
108:                IF (UNRST) THEN 99:                IF (UNRST) THEN
109:                   GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS) = GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS) - &100:                   GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS) = GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS) - &
110:                &    DOT_PRODUCT(GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS),TANVEC(:NINTS,J1-1))*TANVEC(:NINTS,J1-1)101:                &    DOT_PRODUCT(GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS),TANVEC(:NINTS,J1-1))*TANVEC(:NINTS,J1-1)
111:                ELSEIF (RIGIDINIT) THEN 
112:                   GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) = GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) - & 
113:                &    DOT_PRODUCT(GGG(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS),TANVEC(:DEGFREEDOMS,J1-1))*TANVEC(:DEGFREEDOMS,J1-1) 
114:                ELSE102:                ELSE
115:                   GGG(NOPT*(J1-1)+1:NOPT*J1) = GGG(NOPT*(J1-1)+1:NOPT*J1) - &103:                   GGG(NOPT*(J1-1)+1:NOPT*J1) = GGG(NOPT*(J1-1)+1:NOPT*J1) - &
116:                &    DOT_PRODUCT(GGG(NOPT*(J1-1)+1:NOPT*J1),TANVEC(:,J1-1))*TANVEC(:,J1-1)104:                &    DOT_PRODUCT(GGG(NOPT*(J1-1)+1:NOPT*J1),TANVEC(:,J1-1))*TANVEC(:,J1-1)
117:                ENDIF105:                ENDIF
118:                TRUEPERP(1:NOPT*(NIMAGE+2))=GGG(1:NOPT*(NIMAGE+2))106:                TRUEPERP(1:NOPT*(NIMAGE+2))=GGG(1:NOPT*(NIMAGE+2))
119:           ENDDO107:           ENDDO
120:           IF (UNRST) THEN108:           IF (UNRST) THEN
121:              TMPRMS=0.0D0109:              TMPRMS=0.0D0
122:              DO J1=1,NIMAGE110:              DO J1=1,NIMAGE
123:                 TMPRMS = TMPRMS + DOT_PRODUCT(G(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS),G(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS))111:                 TMPRMS = TMPRMS + DOT_PRODUCT(G(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS),G(NOPT*(J1-1)+1:NOPT*(J1-1)+NINTS))
124:              ENDDO112:              ENDDO
125:              RMS=SQRT( TMPRMS/(NIMAGE*NINTS) ) 113:              RMS=SQRT( TMPRMS/(NIMAGE*NINTS) ) 
126:           ELSE114:           ELSE
127:              RMS=SQRT( DOT_PRODUCT(G,G)/(NIMAGE*NOPT) ) 115:              RMS=SQRT( DOT_PRODUCT(G,G)/(NIMAGE*NOPT) ) 
128:           ENDIF116:           ENDIF
129: 117: 
130:           IF (BULKT.AND.(GRADTYPE.NE.'dneb')) THEN118:           IF (BULKT.AND.(GRADTYPE.NE.'dneb')) THEN
131:              PRINT '(A)',' nebgradient> ERROR - minimum image distances only coded for GRADTYPE dneb'119:              PRINT '(A)',' nebgradient> ERROR - minimum image distances only coded for GRADTYPE dneb'
132:              STOP120:              STOP
133:           ENDIF121:           ENDIF
134: 122:           
135:           ! spring gradient on images123:           ! spring gradient on images
136:           SELECT CASE (GRADTYPE) ! THIS CALCULATES SPRING GRADIENT124:           SELECT CASE (GRADTYPE) ! THIS CALCULATES SPRING GRADIENT
137:  
138:              CASE("jnew") ! springs in new implementation of neb by Jonsson125:              CASE("jnew") ! springs in new implementation of neb by Jonsson
139:                 IF(RIGIDINIT) THEN  ! Warning: not tested126:                DO J1=1,NIMAGE
140:                    WRITE(*,*) "Warning: using the new rigid body distance measure with NEB gradient type ", GRADTYPE127:                  GSPR(NOPT*(J1-1)+1:NOPT*J1) = &
141:                    WRITE(*,*) "This combination has not been tested"128:                  & - NEBK*( SQRT(SUM( ( XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) &
142:                    DO J1=1,NIMAGE129:                  &          - SQRT(SUM( ( XYZ(NOPT*(J1-1)+1:NOPT*J1)     - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) )*TANVEC(:,J1)
143:                       CALL RB_DISTANCE(DISTR, XYZ(NOPT*(J1+1)+1:NOPT*(J1+1)+DEGFREEDOMS), XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), &130:                ENDDO
144:                       &                   GRADL1, GRADR1, .FALSE.) 
145:                       CALL RB_DISTANCE(DISTL, XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), & 
146:                       &                   GRADL2, GRADR2, .FALSE.) 
147:                       GSPR(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) = - NEBK*(DISTR - DISTL)*TANVEC(:DEGFREEDOMS,J1)  ! Not sure what to do about TANVEC 
148:                    ENDDO 
149:                 ELSE 
150:                    DO J1=1,NIMAGE 
151:                       GSPR(NOPT*(J1-1)+1:NOPT*J1) = & 
152:                       & - NEBK*( SQRT(SUM( ( XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) & 
153:                       &         - SQRT(SUM( ( XYZ(NOPT*(J1-1)+1:NOPT*J1)     - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) )*TANVEC(:,J1) 
154:                    ENDDO 
155:                 ENDIF 
156:  
157:              CASE("spr") ! full springs:  Gspr = k [2x(j) -x(j-1) -x(j+1)]131:              CASE("spr") ! full springs:  Gspr = k [2x(j) -x(j-1) -x(j+1)]
158:                   IF(RIGIDINIT) THEN   ! Not tested132:                DO J1=1,NIMAGE
159:                    WRITE(*,*) "Warning: using the new rigid body distance measure with NEB gradient type ", GRADTYPE133:                  GSPR(NOPT*(J1-1)+1:NOPT*J1) = &
160:                    WRITE(*,*) "This combination has not been tested"134:                  & NEBK*( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) )
161:                      DO J1=1,NIMAGE135:                ENDDO
162:                         ! GRADL1 is the distance derivative with respect to the first coordinate argument (X_j1) 
163:                         CALL RB_DISTANCE(DISTR, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1+1)+1:NOPT*(J1+1)+DEGFREEDOMS), & 
164:                       &                   GRADL1, GRADR1, .TRUE.) 
165:                         CALL RB_DISTANCE(DISTL, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), & 
166:                       &                   GRADL2, GRADR2, .TRUE.) 
167:                         GSPR(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) = 0.5D0*NEBK*(GRADL1(:)+GRADL2(:)) 
168:                      ENDDO 
169:                   ELSE 
170:                      DO J1=1,NIMAGE 
171:                         GSPR(NOPT*(J1-1)+1:NOPT*J1) = & 
172:                         & NEBK*( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ) 
173:                      ENDDO 
174:                   ENDIF 
175:  
176:              CASE("jold") ! springs MEP-parallel component only (original formulation): G = (Gspr,Tau)*Tau136:              CASE("jold") ! springs MEP-parallel component only (original formulation): G = (Gspr,Tau)*Tau
177:                 IF(RIGIDINIT) THEN  ! Not tested137:                DO J1=1,NIMAGE
178:                    WRITE(*,*) "Warning: using the new rigid body distance measure with NEB gradient type ", GRADTYPE138:                  GSPR(NOPT*(J1-1)+1:NOPT*J1) = &
179:                    WRITE(*,*) "This combination has not been tested"139:                  & NEBK*DOT_PRODUCT( (2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2))), &
180:                    DO J1=1,NIMAGE140:                  & TANVEC(:,J1) )*TANVEC(:,J1)
181:                       ! GRADL1 is the distance derivative with respect to the first coordinate argument (X_j1)141:                ENDDO
182:                       CALL RB_DISTANCE(DISTR, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1+1)+1:NOPT*(J1+1)+DEGFREEDOMS), &142:               CASE("nper")! springs MEP-parallel component (original formulation) + n% of perpendicular component:
183:                     &                   GRADL1, GRADR1, .TRUE.)143:                DO J1=1,NIMAGE   !GJ = GJ + G|_*0.1 + GSPR|| == GJ + [ GSPR - (GSPR,TAU)*TAU ]*N + ...
184:                       CALL RB_DISTANCE(DISTL, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), &144:                  GSPR(NOPT*(J1-1)+1:NOPT*J1) = ( &
185:                     &                   GRADL2, GRADR2, .TRUE.)145:                  &   NEBK*           ( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) )& 
186:                       GSPR(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) = 0.5D0*NEBK*DOT_PRODUCT((GRADL1(:)+GRADL2(:)), &146:                  & - NEBK*DOT_PRODUCT( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ,&
187:                     &         TANVEC(:DEGFREEDOMS,J1))*TANVEC(:DEGFREEDOMS,J1)147:                  & TANVEC(:,J1) )*TANVEC(:,J1) &
188:                    ENDDO148:                                                & )*PERPCOMP
189:                 ELSE 
190:                    DO J1=1,NIMAGE 
191:                       GSPR(NOPT*(J1-1)+1:NOPT*J1) = & 
192:                    &  NEBK*DOT_PRODUCT( (2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2))), & 
193:                    &  TANVEC(:,J1) )*TANVEC(:,J1) 
194:                    ENDDO 
195:                 ENDIF 
196:  
197:              CASE("nper")! springs MEP-parallel component (original formulation) + n% of perpendicular component: 
198:                 IF(RIGIDINIT) THEN  ! Not tested 
199:                    WRITE(*,*) "Warning: using the new rigid body distance measure with NEB gradient type ", GRADTYPE 
200:                    WRITE(*,*) "This combination has not been tested" 
201:                    DO J1=1,NIMAGE 
202:                       ! GRADL1 is the distance derivative with respect to the first coordinate argument (X_j1) 
203:                       CALL RB_DISTANCE(DISTR, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1+1)+1:NOPT*(J1+1)+DEGFREEDOMS), & 
204:                     &                   GRADL1, GRADR1, .TRUE.) 
205:                       CALL RB_DISTANCE(DISTL, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), & 
206:                     &                   GRADL2, GRADR2, .TRUE.) 
207:                       GSPR(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) = 0.5D0*NEBK*(GRADL1(:)+GRADL2(:)) - & 
208:                     &      0.5D0*NEBK*DOT_PRODUCT((GRADL1(:)+GRADL2(:)),TANVEC(:DEGFREEDOMS,J1))*TANVEC(:DEGFREEDOMS,J1)*PERPCOMP 
209:  
210:                       GSPR(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) = GSPR(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) + & 
211:                     &      0.5*NEBK*DOT_PRODUCT((GRADL1(:)+GRADL2(:)),TANVEC(:DEGFREEDOMS,J1))*TANVEC(:DEGFREEDOMS,J1) 
212:                    ENDDO 
213:                 ELSE 
214:                    DO J1=1,NIMAGE   !GJ = GJ + G|_*0.1 + GSPR|| == GJ + [ GSPR - (GSPR,TAU)*TAU ]*N + ... 
215:                       GSPR(NOPT*(J1-1)+1:NOPT*J1) = ( & 
216:                    &   NEBK*           ( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) )& 
217:                    & - NEBK*DOT_PRODUCT( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ,& 
218:                    &   TANVEC(:,J1) )*TANVEC(:,J1) )*PERPCOMP 
219:                  149:                  
220:                       GSPR(NOPT*(J1-1)+1:NOPT*J1) = GSPR(NOPT*(J1-1)+1:NOPT*J1) &  !... (GSPR,TAU)*TAU150:                  GSPR(NOPT*(J1-1)+1:NOPT*J1) = GSPR(NOPT*(J1-1)+1:NOPT*J1) &  !... (GSPR,TAU)*TAU
221:                    &  + NEBK*DOT_PRODUCT( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ,&151:                  & + NEBK*DOT_PRODUCT( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ,&
222:                    &    TANVEC(:,J1) )*TANVEC(:,J1)152:                  & TANVEC(:,J1) )*TANVEC(:,J1)
223:                    ENDDO153:                ENDDO
224:                 ENDIF154:               CASE("nper2")! same as above but optimised for efficiency?
225: 155:                DO J1=1,NIMAGE
226:              CASE("nper2")! same as above but optimised for efficiency?156:                  ! calculate gs (spring gradient on image j1)
227:                 IF(RIGIDINIT) THEN  ! Not tested157:                  GSPR(NOPT*(J1-1)+1:NOPT*J1) = NEBK*(&
228:                    WRITE(*,*) "Warning: using the new rigid body distance measure with NEB gradient type ", GRADTYPE158:                  2*XYZ(NOPT*J1+1:NOPT*(J1+1))-XYZ(NOPT*(J1-1)+1:NOPT*J1)-XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) &
229:                    WRITE(*,*) "This combination has not been tested" 
230:                    DO J1=1,NIMAGE 
231:                       ! GRADL1 is the distance derivative with respect to the first coordinate argument (X_j1) 
232:                       CALL RB_DISTANCE(DISTR, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1+1)+1:NOPT*(J1+1)+DEGFREEDOMS), & 
233:                     &                   GRADL1, GRADR1, .TRUE.) 
234:                       CALL RB_DISTANCE(DISTL, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), & 
235:                     &                   GRADL2, GRADR2, .TRUE.) 
236:                       GSPR(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS) = 0.5D0*NEBK*(GRADL1(:)+GRADL2(:))*PERPCOMP + & 
237:                     &     (1-PERPCOMP)*DOT_PRODUCT(0.5D0*NEBK*(GRADL1(:)+GRADL2(:)),TANVEC(:DEGFREEDOMS,J1))*TANVEC(:DEGFREEDOMS,J1) 
238:                    ENDDO 
239:                 ELSE 
240:                    DO J1=1,NIMAGE 
241:                       ! calculate gs (spring gradient on image j1) 
242:                       GSPR(NOPT*(J1-1)+1:NOPT*J1) = NEBK*(& 
243:                             2*XYZ(NOPT*J1+1:NOPT*(J1+1))-XYZ(NOPT*(J1-1)+1:NOPT*J1)-XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) & 
244:                                                        &)159:                                                        &)
245:                       ! calculate [ f*gs + (1-f)(gs,tau)tau ]160:                  ! calculate [ f*gs + (1-f)(gs,tau)tau ]
246:                       GSPR(NOPT*(J1-1)+1:NOPT*J1) = &161:                  GSPR(NOPT*(J1-1)+1:NOPT*J1) = &
247:                  &    PERPCOMP*GSPR(NOPT*(J1-1)+1:NOPT*J1) &162:                  & PERPCOMP*GSPR(NOPT*(J1-1)+1:NOPT*J1) &
248:                  &        +(1-PERPCOMP)*DOT_PRODUCT(GSPR(NOPT*(J1-1)+1:NOPT*J1),TANVEC(:,J1))*TANVEC(:,J1)163:                  & +(1-PERPCOMP)*DOT_PRODUCT(GSPR(NOPT*(J1-1)+1:NOPT*J1),TANVEC(:,J1))*TANVEC(:,J1)
249:                    ENDDO164:                ENDDO
250:                 ENDIF 
251:               CASE("dneb")  ! this appears to be the default165:               CASE("dneb")  ! this appears to be the default
252:                  DO J1=1,NIMAGE166:                DO J1=1,NIMAGE
253:                  ! calculate SSS = ~g perp in DNEB paper, the perpendicular part of the spring gradient vector167:                  ! calculate SSS = ~g perp in DNEB paper, the perpendicular part of the spring gradient vector
254: !168: !
255: !  Here we are calculating spring gradient vector - (spring gradient vector . tangent vector) tangent vector169: !  Here we are calculating spring gradient vector - (spring gradient vector . tangent vector) tangent vector
256: !  This gives the spring gradient perpendicular component, which is put in SSS.170: !  This gives the spring gradient perpendicular component, which is put in SSS.
257: !171: !
258: !  Changed to use dynamically adjusted NEBK DJW 21/10/08172: !  Changed to use dynamically adjusted NEBK DJW 21/10/08
259: !  Only done for default CASE("dneb")173: !  Only done for default CASE("dneb")
260: !174: !
261:                     IF(RIGIDINIT) THEN ! sn402: added175:                   IF (BULKT) THEN ! minimum image convention for distances!
262:                        ! calculate SSS = ~g perp in DNEB paper, the perpendicular part of the spring gradient vector176:                      DO K=1,NATOMS
263:                        ! Note, the RB_DISTANCE call takes account of PBCs, if necessary.177:                         TEMP1(3*(K-1)+1)=XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1 ) &
264:                        CALL RB_DISTANCE(DISTL, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), &178:    &                       -PARAM1*NINT((XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1))/PARAM1)
265:                     &                   GRADL1, GRADR1, .TRUE.)179:                         TEMP1(3*(K-1)+2)=XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2 ) &
266:                        CALL RB_DISTANCE(DISTR, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1+1)+1:NOPT*(J1+1)+DEGFREEDOMS), &180:    &                       -PARAM2*NINT((XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2))/PARAM2)
267:                     &                   GRADL2, GRADR2, .TRUE.)181:                         IF (.NOT.TWOD) TEMP1(3*(K-1)+3)=XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3 ) &
268: 182:    &                       -PARAM3*NINT((XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3))/PARAM3)
269:                        SSS(NOPT*J1+1:NOPT*J1+DEGFREEDOMS) = 0.5D0*NEWNEBK(J1+1)*GRADL1(:)+0.5D0*NEWNEBK(J1)*GRADL2(:) - &183: 
270:                     &      0.5D0*NEWNEBK(J1+1)*DOT_PRODUCT( GRADL1(:),TANVEC(1:DEGFREEDOMS,J1) )*TANVEC(1:DEGFREEDOMS,J1) - &184:                         TEMP2(3*(K-1)+1)=XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1+1)+3*(K-1)+1 ) &
271:                     &      0.5D0*NEWNEBK(J1)*DOT_PRODUCT( GRADL2(:),TANVEC(1:DEGFREEDOMS,J1) )*TANVEC(1:DEGFREEDOMS,J1)185:    &                       -PARAM1*NINT((XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1+1)+3*(K-1)+1))/PARAM1)
272: 186:                         TEMP2(3*(K-1)+2)=XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1+1)+3*(K-1)+2 ) &
273:                     ELSEIF (BULKT) THEN ! minimum image convention for distances!187:    &                       -PARAM2*NINT((XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1+1)+3*(K-1)+2))/PARAM2)
274:                        DO K=1,NATOMS188:                         IF (.NOT.TWOD) TEMP2(3*(K-1)+3)=XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1+1)+3*(K-1)+3 ) &
275:                           ! TEMP1 holds the minimum-image displacement of image J1 from image J1-1189:    &                       -PARAM3*NINT((XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1+1)+3*(K-1)+3))/PARAM3)
276:                           TEMP1(3*(K-1)+1)=XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1 ) &190:                      ENDDO
277:    &                         -PARAM1*NINT((XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1))/PARAM1)191:                      SSS(NOPT*J1+1:NOPT*(J1+1)) = NEWNEBK(J1)*TEMP1(1:NOPT)+NEWNEBK(J1+1)*TEMP2(1:NOPT) &
278:                           TEMP1(3*(K-1)+2)=XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2 ) &192:    &                    -NEWNEBK(J1)  *DOT_PRODUCT(TEMP1(1:NOPT),TANVEC(1:NOPT,J1) )*TANVEC(1:NOPT,J1) &
279:    &                         -PARAM2*NINT((XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2))/PARAM2)193:    &                    -NEWNEBK(J1+1)*DOT_PRODUCT(TEMP2(1:NOPT),TANVEC(1:NOPT,J1) )*TANVEC(1:NOPT,J1)
280:                           IF (.NOT.TWOD) TEMP1(3*(K-1)+3)=XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3 ) &194:                   ELSE
281:    &                         -PARAM3*NINT((XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3))/PARAM3)195:                      SPRING(NOPT*J1+1:NOPT*(J1+1))=  &
282:                           ! TEMP2 holds the minimum-image displacement of image J1 from image J1+1196:       &   NEWNEBK(J1)  *( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) ) & 
283:                           TEMP2(3*(K-1)+1)=XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1+1)+3*(K-1)+1 ) &197:       & + NEWNEBK(J1+1)*( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ) 
284:    &                          -PARAM1*NINT((XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1+1)+3*(K-1)+1))/PARAM1)198: 
285:                           TEMP2(3*(K-1)+2)=XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1+1)+3*(K-1)+2 ) &199:                      SSS(NOPT*J1+1:NOPT*(J1+1)) = &
286:    &                         -PARAM2*NINT((XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1+1)+3*(K-1)+2))/PARAM2)200:       &   NEWNEBK(J1)  *( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) ) &
287:                           IF (.NOT.TWOD) TEMP2(3*(K-1)+3)=XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1+1)+3*(K-1)+3 ) &201:       & + NEWNEBK(J1+1)*( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ) &
288:    &                         -PARAM3*NINT((XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1+1)+3*(K-1)+3))/PARAM3)202:       & - NEWNEBK(J1)  *DOT_PRODUCT( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1)    , TANVEC(:,J1) )*TANVEC(:,J1) &
289:                        ENDDO203:       & - NEWNEBK(J1+1)*DOT_PRODUCT( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)), TANVEC(:,J1) )*TANVEC(:,J1)
290:                        SSS(NOPT*J1+1:NOPT*(J1+1)) = NEWNEBK(J1)*TEMP1(1:NOPT)+NEWNEBK(J1+1)*TEMP2(1:NOPT) &204:                   ENDIF
291:    &                      -NEWNEBK(J1)  *DOT_PRODUCT(TEMP1(1:NOPT),TANVEC(1:NOPT,J1) )*TANVEC(1:NOPT,J1) & 
292:    &                      -NEWNEBK(J1+1)*DOT_PRODUCT(TEMP2(1:NOPT),TANVEC(1:NOPT,J1) )*TANVEC(1:NOPT,J1) 
293:                     ELSE 
294:                        ! sn402: So far as I can tell, this variable (SPRING) is not used for anything. I will experiment with 
295:                        ! commenting it out. 
296:                        SPRING(NOPT*J1+1:NOPT*(J1+1))=  & 
297:    &                        NEWNEBK(J1)  *( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) ) & 
298:    &                        + NEWNEBK(J1+1)*( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ) 
299:  
300:                        SSS(NOPT*J1+1:NOPT*(J1+1)) = & 
301:    &                       NEWNEBK(J1)  *( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) ) & 
302:    &                       + NEWNEBK(J1+1)*( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ) & 
303:    &      - NEWNEBK(J1)  *DOT_PRODUCT( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1)    , TANVEC(:,J1) )*TANVEC(:,J1) & 
304:    &      - NEWNEBK(J1+1)*DOT_PRODUCT( XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)), TANVEC(:,J1) )*TANVEC(:,J1) 
305:                     ENDIF 
306: !                 PRINT '(A,I6)','grad> SSS spring gradient perpendicular for image ',J1205: !                 PRINT '(A,I6)','grad> SSS spring gradient perpendicular for image ',J1
307: !                 PRINT '(6G20.10)',SSS(1+NOPT*J1:(J1+1)*NOPT)206: !                 PRINT '(6G20.10)',SSS(1+NOPT*J1:(J1+1)*NOPT)
308: 207: 
309: !208: !
310: ! SSS now contains the spring gradient perpendicular part.209: ! SSS now contains the spring gradient perpendicular part.
311: ! SPRING contains the complete spring gradient.210: ! SPRING contains the complete spring gradient.
312: !211: !
313:                        CALL NUDGE(J1,J1)212:                   CALL NUDGE(J1,J1)
314: !                 PRINT '(A,I6)','grad> SSS  g~* for image ',J1213: !                 PRINT '(A,I6)','grad> SSS  g~* for image ',J1
315: !                 PRINT '(6G20.10)',SSS(1+NOPT*J1:(J1+1)*NOPT)214: !                 PRINT '(6G20.10)',SSS(1+NOPT*J1:(J1+1)*NOPT)
316: !215: !
317: ! On exit from NUDGE SSS contains g~* from equation (13) of the DNEB paper.216: ! On exit from NUDGE SSS contains g~* from equation (13) of the DNEB paper.
318: ! GGG still contains g perp, the perpendicular part of the true gradient.217: ! GGG still contains g perp, the perpendicular part of the true gradient.
319: ! These are two of the three terms in equation (12) of the DNEB paper. 
320: !218: !
321: 219:                ENDDO
322:                  ENDDO  ! End of loop over images (J1) 
323: !220: !
324: ! The call to NUDGE calculated a new SSS containing the nudged spring gradient g~*.221: ! The call to NUDGE calculated a new SSS containing the nudged spring gradient g~*.
325: ! We now add this to GGG. We are still missing ~g parallel.222: ! We now add this to GGG. We are still missing ~g parallel.
326: ! GGG now contains those bits of equation (12) that we have calculated thus far. 
327: !223: !
328:                  GGG(NOPT+1:NOPT*(NIMAGE+1)) = GGG(NOPT+1:NOPT*(NIMAGE+1)) + SSS(NOPT+1:NOPT*(NIMAGE+1))224:                GGG(NOPT+1:NOPT*(NIMAGE+1)) = GGG(NOPT+1:NOPT*(NIMAGE+1)) + SSS(NOPT+1:NOPT*(NIMAGE+1))
329: !                 PRINT '(A,I6)','grad> GGG g perp + g~* for image ',J1225: !                 PRINT '(A,I6)','grad> GGG g perp + g~* for image ',J1
330: !                 PRINT '(6G20.10)',GGG(1+NOPT*J1:(J1+1)*NOPT)226: !                 PRINT '(6G20.10)',GGG(1+NOPT*J1:(J1+1)*NOPT)
331: !227: !
332: ! Calculate ~g parallel and store in SSS. This is not added to GGG until after the whole SELECT block.228: ! Calculate ~g parallel and store in SSS. This is not added to GGG until after the whole SELECT block.
333: ! In fact ~g parallel is not used as the formula expected from projection, but uses an229: ! In fact ~g parallel is not used as the formula expected from projection, but uses an
334: ! alternative formulation from equation (5) of the DNEB paper.230: ! alternative formulation from equation (5) of the DNEB paper.
335: !231: !
336:                  IF(RIGIDINIT) THEN232:                DO J1=1,NIMAGE 
337:                     DO J1=1,NIMAGE233:                   IF (BULKT) THEN ! minimum image convention for distances!
338:                        CALL RB_DISTANCE(DISTR, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1+1)+1:NOPT*(J1+1)+DEGFREEDOMS), & 
339:                     &                   GRADL1, GRADR1, .FALSE.) 
340:                        CALL RB_DISTANCE(DISTL, XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), & 
341:                     &                   GRADL2, GRADR2, .FALSE.) 
342:                        SSS(NOPT*J1+1:NOPT*J1+DEGFREEDOMS)= -0.5D0*(NEWNEBK(J1+1)*DISTR - NEWNEBK(J1)*DISTL)*TANVEC(1:DEGFREEDOMS,J1) 
343:                     ENDDO 
344:                  ELSE 
345:                     DO J1=1,NIMAGE 
346:                        IF (BULKT) THEN ! minimum image convention for distances! 
347: 234: 
348:                           DISTA=0.0D0; DISTB=0.0D0235:                      DISTA=0.0D0; DISTB=0.0D0
349:                           DO J2=1,NATOMS236:                      DO J2=1,NATOMS
350:                              DISTA=DISTA + (XYZ(NOPT*(J1+1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1) - &237:                         DISTA=DISTA+ &
351:   &                            PARAM1*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &238:   &                           (XYZ(NOPT*(J1+1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1) - &
352:   &                                        (XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &239:   &               PARAM1*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &
353:   &                            PARAM2*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2240:   &                           (XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &
354:                 IF (.NOT.TWOD) DISTA=DISTA+(XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &241:   &               PARAM2*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2
355:   &                            PARAM3*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2242:    IF (.NOT.TWOD) DISTA=DISTA+(XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &
356: 243:   &               PARAM3*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2
357:                              DISTB=DISTB + (XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1) - &244:                         DISTB=DISTB+ &
358:   &                            PARAM1*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &245:   &                           (XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1) - &
359:   &                                        (XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &246:   &               PARAM1*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &
360:   &                            PARAM2*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2247:   &                           (XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &
361:                 IF (.NOT.TWOD) DISTB=DISTB+(XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &248:   &               PARAM2*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2
362:   &                            PARAM3*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2249:    IF (.NOT.TWOD) DISTB=DISTB+(XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &
363:                           ENDDO250:   &               PARAM3*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2
364:                           SSS(NOPT*J1+1:NOPT*(J1+1))= - ( NEWNEBK(J1+1)*SQRT(DISTA) - NEWNEBK(J1)*SQRT(DISTB) )*TANVEC(1:NOPT,J1)251:                      ENDDO
365:                        ELSE252:                      SSS(NOPT*J1+1:NOPT*(J1+1))= - ( NEWNEBK(J1+1)*SQRT(DISTA) - NEWNEBK(J1)*SQRT(DISTB) )*TANVEC(1:NOPT,J1)
366:                           SSS(NOPT*J1+1:NOPT*(J1+1))= &253:                   ELSE
 254:                      SSS(NOPT*J1+1:NOPT*(J1+1))= &
367:   &               - ( NEWNEBK(J1+1)*SQRT(SUM( ( XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) &255:   &               - ( NEWNEBK(J1+1)*SQRT(SUM( ( XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) &
368:   &                 - NEWNEBK(J1)  *SQRT(SUM( ( XYZ(NOPT*(J1-1)+1:NOPT*J1)     - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) )*TANVEC(:,J1)256:   &                 - NEWNEBK(J1)  *SQRT(SUM( ( XYZ(NOPT*(J1-1)+1:NOPT*J1)     - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) )*TANVEC(:,J1)
369:                        ENDIF257:                   ENDIF
370:                     ENDDO258:                ENDDO
371:                  ENDIF 
372:  
373:               CASE("dneb2") ! same as "dneb", except spring gradient usage is more consistent;259:               CASE("dneb2") ! same as "dneb", except spring gradient usage is more consistent;
374:                DO J1=1,NIMAGE ! THANKS TO DR.~DOMINIC R. ALFONSO FOR POINTING THAT OUT260:                DO J1=1,NIMAGE ! THANKS TO DR.~DOMINIC R. ALFONSO FOR POINTING THAT OUT
375:                  ! calculate Gspr perp261:                  ! calculate Gspr perp
376:                  SSS(NOPT*J1+1:NOPT*(J1+1)) = &262:                  SSS(NOPT*J1+1:NOPT*(J1+1)) = &
377:                  &   NEBK*           ( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) )&263:                  &   NEBK*           ( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) )&
378:                  & - NEBK*DOT_PRODUCT( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ,&264:                  & - NEBK*DOT_PRODUCT( 2*XYZ(NOPT*J1+1:NOPT*(J1+1)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) - XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) ,&
379:                  & TANVEC(:,J1) )*TANVEC(:,J1)265:                  & TANVEC(:,J1) )*TANVEC(:,J1)
380:                  CALL NUDGE(J1,J1)266:                  CALL NUDGE(J1,J1)
381:                ENDDO267:                ENDDO
382: !268: !
602:           ETOTAL=SUM(EIMAGE)488:           ETOTAL=SUM(EIMAGE)
603: 489: 
604:           IF (ORT) THEN ! REMOVE OVERALL ROTATION/TRANSLATION BY FREEZING 6 DEGREES OF FREEDOM490:           IF (ORT) THEN ! REMOVE OVERALL ROTATION/TRANSLATION BY FREEZING 6 DEGREES OF FREEDOM
605:                G(1::NOPT)=0.0D0491:                G(1::NOPT)=0.0D0
606:                G(2::NOPT)=0.0D0492:                G(2::NOPT)=0.0D0
607:                G(3::NOPT)=0.0D0493:                G(3::NOPT)=0.0D0
608:                G(4::NOPT)=0.0D0494:                G(4::NOPT)=0.0D0
609:                G(5::NOPT)=0.0D0495:                G(5::NOPT)=0.0D0
610:                G(7::NOPT)=0.0D0496:                G(7::NOPT)=0.0D0
611:           ENDIF497:           ENDIF
612:  
613:      END SUBROUTINE NEBGRADIENT498:      END SUBROUTINE NEBGRADIENT
614:           499:           
615:      SUBROUTINE NUDGE(JS,JT) ! SUBROUTINE DOES HIGHER-ORDER NUDGINGS; IT MODIFIES GSPR, EATS BOTH GSPR AND G; PREREQ.500:      SUBROUTINE NUDGE(JS,JT) ! SUBROUTINE DOES HIGHER-ORDER NUDGINGS; IT MODIFIES GSPR, EATS BOTH GSPR AND G; PREREQ.
616:           ! - Gspr and G are calculated; js - image on which gsper is nudged by gtper on image jt501:           ! - Gspr and G are calculated; js - image on which gsper is nudged by gtper on image jt
617:           USE NEBDATA502:           USE NEBDATA
618:           USE KEY, ONLY: UNRST503:           USE KEY, ONLY: UNRST
619:           USE GENRIGID, ONLY: RIGIDINIT, DEGFREEDOMS 
620:           IMPLICIT NONE504:           IMPLICIT NONE
621:           INTEGER,INTENT(IN) :: JS,JT505:           INTEGER,INTENT(IN) :: JS,JT
622:           DOUBLE PRECISION DUMMY506:           DOUBLE PRECISION DUMMY
623: !507: !
624: ! On entry, SSS contains the perpendicular spring gradient component508: ! On entry, SSS contains the perpendicular spring gradient component
625: ! GGG must contain a non-unit vector with the perpendicular component of the true potential509: ! GGG must contain a non-unit vector with the perpendicular component of the true potential
626: ! On exit, SSS contains g~* from equation (13) of the DNEB paper510: ! On exit, SSS contains g~* from equation (13) of the DNEB paper
627: !511: !
628:           IF (UNRST) THEN512:           IF (UNRST) THEN
629:              IF (DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*JS+NINTS),GGG(NOPT*JT+1:NOPT*JT+NINTS)) < 0.0D0) THEN513:              IF (DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*JS+NINTS),GGG(NOPT*JT+1:NOPT*JT+NINTS)) < 0.0D0) THEN
630:                   ! calculate the beast Gspr perp - (Gspr perp, Gt perp)* Gt perp/(Gt perp, Gt perp)514:                   ! calculate the beast Gspr perp - (Gspr perp, Gt perp)* Gt perp/(Gt perp, Gt perp)
631:                   SSS(NOPT*JS+1:NOPT*JS+NINTS) = SSS(NOPT*JS+1:NOPT*JS+NINTS) - &515:                   SSS(NOPT*JS+1:NOPT*JS+NINTS) = SSS(NOPT*JS+1:NOPT*JS+NINTS) - &
632:                   & DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*JS+NINTS),GGG(NOPT*JT+1:NOPT*JT+NINTS))*GGG(NOPT*JT+1:NOPT*JT+NINTS)/ &516:                   & DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*JS+NINTS),GGG(NOPT*JT+1:NOPT*JT+NINTS))*GGG(NOPT*JT+1:NOPT*JT+NINTS)/ &
633:                   & DOT_PRODUCT(GGG(NOPT*JT+1:NOPT*JT+NINTS),GGG(NOPT*JT+1:NOPT*JT+NINTS))517:                   & DOT_PRODUCT(GGG(NOPT*JT+1:NOPT*JT+NINTS),GGG(NOPT*JT+1:NOPT*JT+NINTS))
634:              ENDIF518:              ENDIF
635:           ELSE IF(RIGIDINIT) THEN 
636: !             sn402: I don't know why the next line is here. So I've commented it out 
637: !             IF (DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*JS+DEGFREEDOMS),GGG(NOPT*JT+1:NOPT*JT+DEGFREEDOMS)) < 0.0D0) THEN 
638:                   ! calculate the beast Gspr perp - (Gspr perp, Gt perp)* Gt perp/(Gt perp, Gt perp) 
639:                 DUMMY=DOT_PRODUCT(GGG(NOPT*JT+1:NOPT*JT+DEGFREEDOMS),GGG(NOPT*JT+1:NOPT*JT+DEGFREEDOMS)) 
640:                 IF (DUMMY.EQ.0.0D0) DUMMY=1.0D-100 
641:                   SSS(NOPT*JS+1:NOPT*JS+DEGFREEDOMS) = SSS(NOPT*JS+1:NOPT*JS+DEGFREEDOMS) - & 
642:                   & DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*JS+DEGFREEDOMS),GGG(NOPT*JT+1:NOPT*JT+DEGFREEDOMS))*  & 
643:                   & GGG(NOPT*JT+1:NOPT*JT+DEGFREEDOMS) / DUMMY 
644: !             ENDIF 
645:           ELSE519:           ELSE
646:              IF (DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*(JS+1)),GGG(NOPT*JT+1:NOPT*(JT+1))) < 0.0D0) THEN520:              IF (DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*(JS+1)),GGG(NOPT*JT+1:NOPT*(JT+1))) < 0.0D0) THEN
647:                   ! calculate the beast Gspr perp - (Gspr perp, Gt perp)* Gt perp/(Gt perp, Gt perp)521:                   ! calculate the beast Gspr perp - (Gspr perp, Gt perp)* Gt perp/(Gt perp, Gt perp)
648:                 DUMMY=DOT_PRODUCT(GGG(NOPT*JT+1:NOPT*(JT+1)),GGG(NOPT*JT+1:NOPT*(JT+1)))522:                 DUMMY=DOT_PRODUCT(GGG(NOPT*JT+1:NOPT*(JT+1)),GGG(NOPT*JT+1:NOPT*(JT+1)))
649:                 IF (DUMMY.EQ.0.0D0) DUMMY=1.0D-100523:                 IF (DUMMY.EQ.0.0D0) DUMMY=1.0D-100
650:                   SSS(NOPT*JS+1:NOPT*(JS+1)) = SSS(NOPT*JS+1:NOPT*(JS+1)) - &524:                   SSS(NOPT*JS+1:NOPT*(JS+1)) = SSS(NOPT*JS+1:NOPT*(JS+1)) - &
651:                   & DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*(JS+1)),GGG(NOPT*JT+1:NOPT*(JT+1)))*GGG(NOPT*JT+1:NOPT*(JT+1))/ &525:                   & DOT_PRODUCT(SSS(NOPT*JS+1:NOPT*(JS+1)),GGG(NOPT*JT+1:NOPT*(JT+1)))*GGG(NOPT*JT+1:NOPT*(JT+1))/ &
652:                   & DUMMY526:                   & DUMMY
653:              ENDIF527:              ENDIF
654:           ENDIF528:           ENDIF


r31425/newneb.f90 2016-11-02 17:30:14.374305564 +0000 r31424/newneb.f90 2016-11-02 17:30:15.114315538 +0000
220:              ENDIF220:              ENDIF
221:              IF (DESMINT) THEN221:              IF (DESMINT) THEN
222:                 CALL GROWSTRING(QQ, FINFIN, NIMAGE, XCART, EIMAGE, TANPTR,RMS,GSMFLAG)222:                 CALL GROWSTRING(QQ, FINFIN, NIMAGE, XCART, EIMAGE, TANPTR,RMS,GSMFLAG)
223:                 XYZ(:) = 0.0D0223:                 XYZ(:) = 0.0D0
224:              ELSE224:              ELSE
225:                 CALL GROWSTRING(QQ, FINFIN, NIMAGE, X, EIMAGE, TANPTR,RMS,GSMFLAG)225:                 CALL GROWSTRING(QQ, FINFIN, NIMAGE, X, EIMAGE, TANPTR,RMS,GSMFLAG)
226:                 CALL IMAGEDISTRIBUTION226:                 CALL IMAGEDISTRIBUTION
227:              ENDIF227:              ENDIF
228:              NITERDONE = TOTSTEPS228:              NITERDONE = TOTSTEPS
229:              229:              
230:           ELSE  ! construct the band230:           ELSE
231: 231:              ! construct the band
232:              ! Read an initial band in from a file, if specified 
233:              IF (READGUESS.OR.(GUESSPATHT.AND.UNRST.AND.(NINTERP.GT.1)).OR.(MECCANOT)) THEN232:              IF (READGUESS.OR.(GUESSPATHT.AND.UNRST.AND.(NINTERP.GT.1)).OR.(MECCANOT)) THEN
234:                 CALL RWG("r",.True.,1)233:                 CALL RWG("r",.True.,1)
235:                 READGUESS = .FALSE.234:                 READGUESS = .FALSE.
236: 235: 
237:                 IF (RIGIDINIT) THEN236:                 IF (RIGIDINIT) THEN
238:                     ! Read coordinates in as atomistic xyz, but we need rigid body coordinates.237:                     ! Read coordinates in as atomistic xyz, but we need rigid body coordinates.
239:                     ! It's possible that the xyz input file doesn't respect the rigid body constraints238:                     ! It's possible that the xyz input file doesn't respect the rigid body constraints
240:                     ! If so, converting to rigid body coordinates will restore the rigid body constraints in the 239:                     ! If so, converting to rigid body coordinates will restore the rigid body constraints in the 
241:                     ! best alignment MINPERMDIST can find to the cartesian input coords.240:                     ! best alignment MINPERMDIST can find to the cartesian input coords.
242:                     ! This will not work very well if the 'rigid bodies' in the guessed path are very different241:                     ! This will not work very well if the 'rigid bodies' in the guessed path are very different
247:                                     ! read in don't exactly match the rigid body constraints. But this is expected246:                                     ! read in don't exactly match the rigid body constraints. But this is expected
248:                                     ! behaviour in this case.247:                                     ! behaviour in this case.
249:                     CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)248:                     CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)
250:                     DEBUG = LDEBUG249:                     DEBUG = LDEBUG
251:                 ENDIF 250:                 ENDIF 
252:              ELSE251:              ELSE
253:                 IF (UNRST) THEN ! JMC252:                 IF (UNRST) THEN ! JMC
254:                    CALL UNRESDIHENEB(QQ,FINFIN,MYPTS)253:                    CALL UNRESDIHENEB(QQ,FINFIN,MYPTS)
255:                    XYZ(NOPT+1:NOPT*(NIMAGE+1))=MYPTS(1:NOPT*NIMAGE)254:                    XYZ(NOPT+1:NOPT*(NIMAGE+1))=MYPTS(1:NOPT*NIMAGE)
256:                 ELSEIF (REDOPATHNEB) THEN255:                 ELSEIF (REDOPATHNEB) THEN
257:                    ! sn402: There's probably no rigid body support with this option. 
258:                    REDOKADD=.TRUE.256:                    REDOKADD=.TRUE.
259:                    REDOPATH1=.TRUE.257:                    REDOPATH1=.TRUE.
260:                    REDOTSIM=NIMAGE*D1INIT/(D1INIT+D2INIT)+1258:                    REDOTSIM=NIMAGE*D1INIT/(D1INIT+D2INIT)+1
261:                    XYZ(NOPT*REDOTSIM+1:NOPT*(REDOTSIM+1))=TSREDO(1:NOPT)259:                    XYZ(NOPT*REDOTSIM+1:NOPT*(REDOTSIM+1))=TSREDO(1:NOPT)
262:                    ALLOCATE(DELTAX(NOPT))260:                    ALLOCATE(DELTAX(NOPT))
263:                    IF (BULKT) THEN261:                    IF (BULKT) THEN
264:                       DO K=1,NATOMS262:                       DO K=1,NATOMS
265:                          DELTAX(3*(K-1)+1)=XYZ(NOPT*REDOTSIM+3*(K-1)+1) - XYZ(3*(K-1)+1) &263:                          DELTAX(3*(K-1)+1)=XYZ(NOPT*REDOTSIM+3*(K-1)+1) - XYZ(3*(K-1)+1) &
266:   &                          -PARAM1*NINT((XYZ(NOPT*REDOTSIM+3*(K-1)+1) - XYZ(3*(K-1)+1))/PARAM1)264:   &                          -PARAM1*NINT((XYZ(NOPT*REDOTSIM+3*(K-1)+1) - XYZ(3*(K-1)+1))/PARAM1)
267:                          DELTAX(3*(K-1)+2)=XYZ(NOPT*REDOTSIM+3*(K-1)+2) - XYZ(3*(K-1)+2) &265:                          DELTAX(3*(K-1)+2)=XYZ(NOPT*REDOTSIM+3*(K-1)+2) - XYZ(3*(K-1)+2) &
324:                    KNOWE=.FALSE.; KNOWG=.FALSE.; KNOWH=.FALSE.322:                    KNOWE=.FALSE.; KNOWG=.FALSE.; KNOWH=.FALSE.
325:                    REDOKADD=.FALSE.323:                    REDOKADD=.FALSE.
326:                    REDOPATH2=.FALSE.324:                    REDOPATH2=.FALSE.
327:                 ELSEIF (REDOPATH) THEN325:                 ELSEIF (REDOPATH) THEN
328:                    XYZ(NOPT+1:NOPT*2) = TSREDO(1:NOPT)326:                    XYZ(NOPT+1:NOPT*2) = TSREDO(1:NOPT)
329:                    PRINT '(A)','newneb> Setting tangent vector and hence initial eigenvector guess to difference between minima'327:                    PRINT '(A)','newneb> Setting tangent vector and hence initial eigenvector guess to difference between minima'
330:                    TANVEC(1:NOPT,1)=XYZ(1:NOPT)-XYZ(NOPT*2+1:NOPT*3)328:                    TANVEC(1:NOPT,1)=XYZ(1:NOPT)-XYZ(NOPT*2+1:NOPT*3)
331:                    EEE(2)=1.0D100329:                    EEE(2)=1.0D100
332:                 ELSE330:                 ELSE
333:                    IF (INTCONSTRAINTT) THEN331:                    IF (INTCONSTRAINTT) THEN
334:                       ! QCI method to construct initial band. 
335:  
336: !                     XYZ(1:NOPT)=QQ(1:NOPT)332: !                     XYZ(1:NOPT)=QQ(1:NOPT)
337: !                     XYZ(NOPT+1:NOPT*(NIMAGE+1))=INTNEBIMAGES(1:NOPT*NIMAGE)333: !                     XYZ(NOPT+1:NOPT*(NIMAGE+1))=INTNEBIMAGES(1:NOPT*NIMAGE)
338:                       IF (.NOT.ALLOCATED(QCIXYZ)) THEN334:                       IF (.NOT.ALLOCATED(QCIXYZ)) THEN
339:                          PRINT *,'newneb> ERROR *** QCIXYZ is not allocated'335:                          PRINT *,'newneb> ERROR *** QCIXYZ is not allocated'
340:                          STOP336:                          STOP
341:                       ENDIF337:                       ENDIF
342:                       IF (MOREPRINTING) CALL KEYNEBPRINT338:                       IF (MOREPRINTING) CALL KEYNEBPRINT
343:                       XYZ(1:NOPT*(NIMAGE+2))=QCIXYZ(1:NOPT*(NIMAGE+2))339:                       XYZ(1:NOPT*(NIMAGE+2))=QCIXYZ(1:NOPT*(NIMAGE+2))
344:                       DEALLOCATE(QCIXYZ)340:                       DEALLOCATE(QCIXYZ)
345: 341: 
421:                          XYZ(NOPT+1:NOPT*(NIMAGE+1))=INTNEBIMAGES(1:NOPT*NIMAGE)417:                          XYZ(NOPT+1:NOPT*(NIMAGE+1))=INTNEBIMAGES(1:NOPT*NIMAGE)
422:                          DEALLOCATE(INTNEBIMAGES)418:                          DEALLOCATE(INTNEBIMAGES)
423:                          LDEBUG=DEBUG419:                          LDEBUG=DEBUG
424:                          DEBUG=.FALSE. 420:                          DEBUG=.FALSE. 
425:                          CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)421:                          CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)
426:                          DEBUG=LDEBUG422:                          DEBUG=LDEBUG
427:                       ENDIF423:                       ENDIF
428:                    ELSEIF (INTINTERPT) THEN424:                    ELSEIF (INTINTERPT) THEN
429:                       CALL INTINTERPOLATE(QQ,FINFIN,NINTIM,NIMAGE,X,DESMDEBUG,FAILED)425:                       CALL INTINTERPOLATE(QQ,FINFIN,NINTIM,NIMAGE,X,DESMDEBUG,FAILED)
430:                    ELSE426:                    ELSE
431:                       ! This is where we usually end up! Construct a band completely from scratch using 
432:                       ! (spherical) linear interpolation. 
433:  
434:                       IF(RIGIDINIT .AND. (RIGIDMOLECULEST .OR. SLERPT)) THEN427:                       IF(RIGIDINIT .AND. (RIGIDMOLECULEST .OR. SLERPT)) THEN
435:                       ! MAKEIMAGE doesn't actually use the initial/final coordinates passed into it unless428:                       ! MAKEIMAGE doesn't actually use the initial/final coordinates passed into it unless
436:                       ! we have MORPHT set (it just uses the endpoints in XYZ instead) so the following429:                       ! we have MORPHT set (it just uses the endpoints in XYZ instead) so the following
437:                       ! line is usually unnecessary.430:                       ! line is usually unnecessary.
438:                         CALL MAKEIMAGE(EINITIAL,EFINAL,RIGIDQ,RIGIDFIN)431:                         CALL MAKEIMAGE(EINITIAL,EFINAL,RIGIDQ,RIGIDFIN)
439:                       ELSE432:                       ELSE
440:                         CALL MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN)433:                         CALL MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN)
441:                       ENDIF434:                       ENDIF
442:                    ENDIF435:                    ENDIF
443:                 ENDIF436:                 ENDIF
444:              ENDIF437:              ENDIF
445:              IF (UNRST) DEALLOCATE(MYPTS) ! JMC438:              IF (UNRST) DEALLOCATE(MYPTS) ! JMC
446: 439: 
447: ! hk286440: ! hk286
448:             ! Currently, the default with RIGIDINIT is still to do unrigidified Cartesian interpolation, which allows441:              IF (RIGIDINIT .AND. .NOT. (SLERPT .OR. RIGIDMOLECULEST)) THEN  ! sn402: This whole block could possibly
449:             ! the rigid bodies to deform.442:              ! be removed if the earlier block is changed as suggested
450:             ! Then when we go through this block and convert to rigid body coordinates the atoms get forced 
451:             ! back into their correct rigid body structures. 
452:             ! If the keywords RIGIDMOLECULES (intended for a system entirely composed of identical rigid bodies) 
453:             ! or SLERP are set, then instead we convert to rigid body coordinates before interpolation and do the interpolation 
454:             ! using the iSLERP procedure. In that case, the following block is unnecessary. 
455:             ! In the future, we may change the default behaviour so that rigid bodies are always interpolated correctly. 
456:              IF (RIGIDINIT .AND. .NOT. (SLERPT .OR. RIGIDMOLECULEST)) THEN 
457:                 CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)443:                 CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)
458:                 ATOMRIGIDCOORDT = .FALSE.444:                 ATOMRIGIDCOORDT = .FALSE.
459:              ENDIF445:              ENDIF
460: ! hk286446: ! hk286
461: 447: 
462:              ! preoptimise if requested        448:              ! preoptimise if requested        
463:              IF (SQVVGUESS) THEN449:              IF (SQVVGUESS) THEN
464:                 STEPTOT = 5.0D0 ! TO AVOID GRADIENT SCALING DURING SQVV450:                 STEPTOT = 5.0D0 ! TO AVOID GRADIENT SCALING DURING SQVV
465:                 CALL NEBSQVV(NOPT*NIMAGE)451:                 CALL NEBSQVV(NOPT*NIMAGE)
466:              ENDIF452:              ENDIF
467: 453: 
468: !!!!!!!!!!!!! Perform the actual NEB job here !!!!!!!!!!!!! 
469:  
470:              NPERSIST=0454:              NPERSIST=0
471:              ! use the c++ neb routines455:              ! use the c++ neb routines
472:              IF (CPPNEBT) THEN456:              IF (CPPNEBT) THEN
473:                 IF(RIGIDINIT) THEN457:                 IF(RIGIDINIT) THEN
474:                     ! Seems strange to call this with the number of coords to458:                     ! Seems strange to call this with the number of coords to
475:                     ! optimise set to NOPT (which is the atomistic value)459:                     ! optimise set to NOPT (which is the atomistic value)
476:                     ! rather than DEGFREEDOMS. But in fact this is used to460:                     ! rather than DEGFREEDOMS. But in fact this is used to
477:                     ! initialise arrays, so if we use DEGFREEDOMS then we end up461:                     ! initialise arrays, so if we use DEGFREEDOMS then we end up
478:                     ! running over the ends of arrays.462:                     ! running over the ends of arrays.
479:                     call neb_example(NIMAGE,NOPT,XYZ,EEE,NITERMAX)463:                     call neb_example(NIMAGE,NOPT,XYZ,EEE,NITERMAX)
 464: !                    call neb_example(RIGIDQ,RIGIDFIN,NIMAGE,NOPT,XYZ,EEE,NITERMAX)
480:                 ELSE465:                 ELSE
 466: !                    call neb_example(QQ,FINFIN,NIMAGE,NOPT,XYZ,EEE,NITERMAX)
481:                     CALL NEB_EXAMPLE(NIMAGE,NOPT,XYZ,EEE,NITERMAX)467:                     CALL NEB_EXAMPLE(NIMAGE,NOPT,XYZ,EEE,NITERMAX)
482:                 ENDIF468:                 ENDIF
483:              ELSEIF ((.NOT.REDOPATH).OR.REDOPATHNEB) THEN469:              ELSEIF ((.NOT.REDOPATH).OR.REDOPATHNEB) THEN
484:                 SELECT CASE(MINTYPE)470:                 SELECT CASE(MINTYPE)
485:                 CASE("lbfgs")471:                 CASE("lbfgs")
486:                    IF (UNRST) THEN472:                    IF (UNRST) THEN
487:                       CALL NEBBFGSINT(NINTS*NIMAGE,NEBMUPDATE)473:                       CALL NEBBFGSINT(NINTS*NIMAGE,NEBMUPDATE)
488:                    ELSE474:                    ELSE
489:                       CALL NEBBFGS(NOPT*NIMAGE,NEBMUPDATE,NPERSIST,PERSISTENT)475:                       CALL NEBBFGS(NOPT*NIMAGE,NEBMUPDATE,NPERSIST,PERSISTENT)
490:                    END IF476:                    END IF
491:                 CASE("sqvv")477:                 CASE("sqvv")
492:                    CALL NEBSQVV(NOPT*NIMAGE)478:                    CALL NEBSQVV(NOPT*NIMAGE)
493:                 END SELECT479:                 END SELECT
494:              ENDIF480:              ENDIF
495:  
496: !!!!!!!!!!!!! End of main NEB job !!!!!!!!!!!!! 
497:  
498:           ENDIF481:           ENDIF
499:           ! save final NEB coordinates and energy profile482:           ! save final NEB coordinates and energy profile
500:           IF (DEBUG) THEN483:           IF (DEBUG) THEN
501:              PRINT '(A,F12.4)',' newneb> mean image separation is ',SEPARATION/(NIMAGE+1)484:              PRINT '(A,F12.4)',' newneb> mean image separation is ',SEPARATION/(NIMAGE+1)
502: !            DO J1=1,NIMAGE+1485: !            DO J1=1,NIMAGE+1
503: !               PRINT '(A,F12.4,A,I8,A,F12.4)',' newneb> NEB k is ',NEWNEBK(J1),' for gap ',J1,' value=',DVEC(J1)486: !               PRINT '(A,F12.4,A,I8,A,F12.4)',' newneb> NEB k is ',NEWNEBK(J1),' for gap ',J1,' value=',DVEC(J1)
504: !            ENDDO487: !            ENDDO
505:           ENDIF488:           ENDIF
506: 489: 
507:           IF (DUMPNEBEOS) CALL WRITEPROFILE(0)490:           IF (DUMPNEBEOS) CALL WRITEPROFILE(0)
524:                EMAX=MAXVAL(EIMAGE)507:                EMAX=MAXVAL(EIMAGE)
525:                DO J1=1,NIMAGE508:                DO J1=1,NIMAGE
526:                     IF (EMAX == EIMAGE(J1)) JMAX=J1509:                     IF (EMAX == EIMAGE(J1)) JMAX=J1
527:                ENDDO510:                ENDDO
528:                QQ = X(NOPT*(JMAX-1)+1:NOPT*JMAX)511:                QQ = X(NOPT*(JMAX-1)+1:NOPT*JMAX)
529:           ENDIF512:           ENDIF
530:           513:           
531:           CALL MYCPU_TIME(ENDTIME,.FALSE.)514:           CALL MYCPU_TIME(ENDTIME,.FALSE.)
532: 515: 
533:           WRITE(*,*) "Time to go through NEB: ", ENDTIME-STARTTIME516:           WRITE(*,*) "Time to go through NEB: ", ENDTIME-STARTTIME
534:           STOP ! sn402: remove!517: 
535:           NMINFOUND=0518:           NMINFOUND=0
536:           IF (TSRESET) NTSFOUND=0519:           IF (TSRESET) NTSFOUND=0
537: !520: !
538: !  Current structure precludes searching the NEB profile for521: !  Current structure precludes searching the NEB profile for
539: !  both minima and ts. Could perhaps change this. Current philosophy is522: !  both minima and ts. Could perhaps change this. Current philosophy is
540: !  that if we have persistent minima we should start the whole DNEB again523: !  that if we have persistent minima we should start the whole DNEB again
541: !  without ts searches.524: !  without ts searches.
542: !525: !
543:           IF (NPERSIST.GT.0) THEN526:           IF (NPERSIST.GT.0) THEN
544:              PTEST=.FALSE.527:              PTEST=.FALSE.


r31425/tau.f90 2016-11-02 17:30:14.618308853 +0000 r31424/tau.f90 2016-11-02 17:30:15.362318880 +0000
 20:      IMPLICIT NONE 20:      IMPLICIT NONE
 21:      CONTAINS 21:      CONTAINS
 22:  22: 
 23:      SUBROUTINE MEPTANGENT 23:      SUBROUTINE MEPTANGENT
 24:           USE NEBDATA 24:           USE NEBDATA
 25:           USE KEYNEB,ONLY: NIMAGE 25:           USE KEYNEB,ONLY: NIMAGE
 26:           USE MODUNRES 26:           USE MODUNRES
 27:           USE KEYTAU 27:           USE KEYTAU
 28:           USE KEY, ONLY: BULKT, TWOD 28:           USE KEY, ONLY: BULKT, TWOD
 29:           USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3 29:           USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3
 30:           USE GENRIGID, ONLY: RIGIDINIT, RB_DISTANCE, DEGFREEDOMS 
 31:           IMPLICIT NONE 30:           IMPLICIT NONE
 32:  31: 
 33:           INTEGER :: J1,J2,K 32:           INTEGER :: J1,J2,K
 34:           DOUBLE PRECISION :: WPLUS,WMINUS,DUMMY,TEMP1(NOPT),TEMP2(NOPT) 33:           DOUBLE PRECISION :: WPLUS,WMINUS,DUMMY,TEMP1(NOPT),TEMP2(NOPT)
 35:           DOUBLE PRECISION :: QINT(NINTS*(NIMAGE+2)) 34:           DOUBLE PRECISION :: QINT(NINTS*(NIMAGE+2))
 36:           DOUBLE PRECISION :: DIFFM(NINTS),DIFFP(NINTS) 35:           DOUBLE PRECISION :: DIFFM(NINTS),DIFFP(NINTS)
 37:           DOUBLE PRECISION :: PI=3.141592653589793D0 36:           DOUBLE PRECISION :: PI=3.141592653589793D0
 38:           DOUBLE PRECISION :: DIST, GRADL(DEGFREEDOMS), GRADR(DEGFREEDOMS), SAVEGRADR(DEGFREEDOMS) 
 39:  37: 
 40:           IF ((TANTYPE.NE.1).AND.(BULKT)) THEN 38:           IF ((TANTYPE.NE.1).AND.(BULKT)) THEN
 41:              PRINT '(A)','meptangent> ERROR - minimum image distances only coded for TANTYPE 1' 39:              PRINT '(A)','meptangent> ERROR - minimum image distances only coded for TANTYPE 1'
 42:              STOP 40:              STOP
 43:           ENDIF 41:           ENDIF
 44:           IF ((TANTYPE.NE.1).AND.(RIGIDINIT)) THEN 
 45:              PRINT '(A)','meptangent> ERROR - Rigid body distance measure only coded for TANTYPE 1' 
 46:              STOP 
 47:           ENDIF 
 48:  42: 
 49:           SELECT CASE (TANTYPE) 43:           SELECT CASE (TANTYPE)
 50:              CASE(1) ! Henkelman and Jonsson's improved tangent from JCP, 113, 9978, 2000 44:              CASE(1) ! Henkelman and Jonsson's improved tangent from JCP, 113, 9978, 2000
 51:                 IF(RIGIDINIT) THEN 45:                DO J1=2,NIMAGE+1
 52:                    ! Get the gradient of the leftmost spring's distance with respect to the first image's coordinates (as SAVEGRADR) 46:                   CALL WS(EEE(J1-1),EEE(J1),EEE(J1+1),WMINUS,WPLUS)
 53:                    CALL RB_DISTANCE(DIST, XYZ(1:DEGFREEDOMS), XYZ(NOPT+1:NOPT+DEGFREEDOMS), & 47:                   TEMP1(1:NOPT)=0.0D0
 54:                     &                   GRADL, SAVEGRADR, .TRUE.) 48:                   TEMP2(1:NOPT)=0.0D0
 55:                    DO J1=2,NIMAGE+1 49:                   IF (BULKT) THEN ! minimum image convention for distances!
 56:                       ! Get the weights for interpolating the tangent 50:                      DO K=1,NATOMS
 57:                       CALL WS(EEE(J1-1),EEE(J1),EEE(J1+1),WMINUS,WPLUS) 
 58:                       ! Get the next pair of distance gradients  (image 1 to 2, then 2 to 3, and so on) 
 59:                       CALL RB_DISTANCE(DIST, XYZ(NOPT*(J1-1)+1:NOPT*(J1-1)+DEGFREEDOMS), XYZ(NOPT*J1+1:NOPT*J1+DEGFREEDOMS), & 
 60:                     &                   GRADL, GRADR, .TRUE.) 
 61:  
 62:                       TANVEC(:DEGFREEDOMS,J1-1) = WMINUS*SAVEGRADR(:) - WPLUS*GRADL(:) 
 63:                       TANVEC(DEGFREEDOMS+1:,J1-1)=0.0D0 
 64:  
 65:                       SAVEGRADR(:)=GRADR(:) 
 66:  
 67:                    ENDDO 
 68:                 ELSE 
 69:                    DO J1=2,NIMAGE+1 
 70:                       CALL WS(EEE(J1-1),EEE(J1),EEE(J1+1),WMINUS,WPLUS) 
 71:                       TEMP1(1:NOPT)=0.0D0 
 72:                       TEMP2(1:NOPT)=0.0D0 
 73:                       IF (BULKT) THEN ! minimum image convention for distances! 
 74:                          DO K=1,NATOMS 
 75:                         TEMP1(3*(K-1)+1)=XYZ(NOPT*(J1-1)+3*(K-1)+1) - XYZ(NOPT*(J1-2)+3*(K-1)+1 ) & 51:                         TEMP1(3*(K-1)+1)=XYZ(NOPT*(J1-1)+3*(K-1)+1) - XYZ(NOPT*(J1-2)+3*(K-1)+1 ) &
 76:    &                       -PARAM1*NINT((XYZ(NOPT*(J1-1)+3*(K-1)+1) - XYZ(NOPT*(J1-2)+3*(K-1)+1))/PARAM1) 52:    &                       -PARAM1*NINT((XYZ(NOPT*(J1-1)+3*(K-1)+1) - XYZ(NOPT*(J1-2)+3*(K-1)+1))/PARAM1)
 77:                         TEMP1(3*(K-1)+2)=XYZ(NOPT*(J1-1)+3*(K-1)+2) - XYZ(NOPT*(J1-2)+3*(K-1)+2 ) & 53:                         TEMP1(3*(K-1)+2)=XYZ(NOPT*(J1-1)+3*(K-1)+2) - XYZ(NOPT*(J1-2)+3*(K-1)+2 ) &
 78:    &                       -PARAM2*NINT((XYZ(NOPT*(J1-1)+3*(K-1)+2) - XYZ(NOPT*(J1-2)+3*(K-1)+2))/PARAM2) 54:    &                       -PARAM2*NINT((XYZ(NOPT*(J1-1)+3*(K-1)+2) - XYZ(NOPT*(J1-2)+3*(K-1)+2))/PARAM2)
 79:                         IF (.NOT.TWOD) TEMP1(3*(K-1)+3)=XYZ(NOPT*(J1-1)+3*(K-1)+3) - XYZ(NOPT*(J1-2)+3*(K-1)+3 ) & 55:                         IF (.NOT.TWOD) TEMP1(3*(K-1)+3)=XYZ(NOPT*(J1-1)+3*(K-1)+3) - XYZ(NOPT*(J1-2)+3*(K-1)+3 ) &
 80:    &                       -PARAM3*NINT((XYZ(NOPT*(J1-1)+3*(K-1)+3) - XYZ(NOPT*(J1-2)+3*(K-1)+3))/PARAM3) 56:    &                       -PARAM3*NINT((XYZ(NOPT*(J1-1)+3*(K-1)+3) - XYZ(NOPT*(J1-2)+3*(K-1)+3))/PARAM3)
 81:  57: 
 82:                         TEMP2(3*(K-1)+1)=XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1 ) & 58:                         TEMP2(3*(K-1)+1)=XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1 ) &
 83:    &                       -PARAM1*NINT((XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1))/PARAM1) 59:    &                       -PARAM1*NINT((XYZ(NOPT*J1+3*(K-1)+1) - XYZ(NOPT*(J1-1)+3*(K-1)+1))/PARAM1)
 84:                         TEMP2(3*(K-1)+2)=XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2 ) & 60:                         TEMP2(3*(K-1)+2)=XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2 ) &
 85:    &                       -PARAM2*NINT((XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2))/PARAM2) 61:    &                       -PARAM2*NINT((XYZ(NOPT*J1+3*(K-1)+2) - XYZ(NOPT*(J1-1)+3*(K-1)+2))/PARAM2)
 86:                         IF (.NOT.TWOD) TEMP2(3*(K-1)+3)=XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3 ) & 62:                         IF (.NOT.TWOD) TEMP2(3*(K-1)+3)=XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3 ) &
 87:    &                       -PARAM3*NINT((XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3))/PARAM3) 63:    &                       -PARAM3*NINT((XYZ(NOPT*J1+3*(K-1)+3) - XYZ(NOPT*(J1-1)+3*(K-1)+3))/PARAM3)
 88:                          ENDDO 64:                      ENDDO
 89:                          TANVEC(1:NOPT,J1-1)=WMINUS*TEMP1(1:NOPT) + WPLUS*TEMP2(1:NOPT) 65:                      TANVEC(1:NOPT,J1-1)=WMINUS*TEMP1(1:NOPT) + WPLUS*TEMP2(1:NOPT)
 90:                       ELSE 66:                   ELSE
 91:                         TANVEC(:,J1-1) = WMINUS*( XYZ(NOPT*(J1-1)+1:NOPT*J1)         - XYZ(NOPT*(J1-2)+1:NOPT*(J1-1))         ) & 67:                      TANVEC(:,J1-1) = WMINUS*( XYZ(NOPT*(J1-1)+1:NOPT*J1)         - XYZ(NOPT*(J1-2)+1:NOPT*(J1-1))         ) &
 92:                                  &  + WPLUS *( XYZ(NOPT*J1+1:NOPT*(J1+1))         - XYZ(NOPT*(J1-1)+1:NOPT*J1    )         ) 68:                                  &  + WPLUS *( XYZ(NOPT*J1+1:NOPT*(J1+1))         - XYZ(NOPT*(J1-1)+1:NOPT*J1    )         )
 93:                       ENDIF 69:                   ENDIF
 94:  70:                ENDDO
 95:                    ENDDO 
 96:                 ENDIF 
 97:              CASE(2) !NORMALIZED LINE SEGMENT (TANGENT ESTIMATED FROM TWO ADJACENT IMAGES) 71:              CASE(2) !NORMALIZED LINE SEGMENT (TANGENT ESTIMATED FROM TWO ADJACENT IMAGES)
 98:                DO J1=1,NIMAGE 72:                DO J1=1,NIMAGE
 99:                     TANVEC(:,J1) = XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*(J1-1)+1:NOPT*J1) 73:                     TANVEC(:,J1) = XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*(J1-1)+1:NOPT*J1)
100:                ENDDO 74:                ENDDO
101:              CASE(4) ! JMC FOR INTERNALS... 75:              CASE(4) ! JMC FOR INTERNALS...
102:                TANVEC=0.0D0 ! INITIALISE TANVEC 76:                TANVEC=0.0D0 ! INITIALISE TANVEC
103:                DO J1=1,NIMAGE+2 77:                DO J1=1,NIMAGE+2
104:                   DO J2=1,NRES 78:                   DO J2=1,NRES
105:                      C(1,J2)=XYZ(6*(J2-1)+1+NOPT*(J1-1)) 79:                      C(1,J2)=XYZ(6*(J2-1)+1+NOPT*(J1-1))
106:                      C(2,J2)=XYZ(6*(J2-1)+2+NOPT*(J1-1)) 80:                      C(2,J2)=XYZ(6*(J2-1)+2+NOPT*(J1-1))
135: 109: 
136:                IF (DUMMY==0.0D0) THEN110:                IF (DUMMY==0.0D0) THEN
137:                     PRINT '(1x,a)', 'Tau norm is zero! Probably image density is too big...'111:                     PRINT '(1x,a)', 'Tau norm is zero! Probably image density is too big...'
138:                     PRINT '(1x,a,i5)', ' This is image # ',j1112:                     PRINT '(1x,a,i5)', ' This is image # ',j1
139:                     BADTAU=.TRUE.113:                     BADTAU=.TRUE.
140:                     RETURN114:                     RETURN
141:                ENDIF115:                ENDIF
142: 116: 
143:                TANVEC(:,J1) = TANVEC(:,J1) / DUMMY117:                TANVEC(:,J1) = TANVEC(:,J1) / DUMMY
144:           ENDDO118:           ENDDO
145:  
146:      END SUBROUTINE MEPTANGENT119:      END SUBROUTINE MEPTANGENT
147: 120: 
148:      SUBROUTINE WS(EL,EM,ER,WM,WP) ! ASSUMES THAT EL,EM AND ER ARE CHECKED FOR NAN ELSEWHERE ...121:      SUBROUTINE WS(EL,EM,ER,WM,WP) ! ASSUMES THAT EL,EM AND ER ARE CHECKED FOR NAN ELSEWHERE ...
149:           IMPLICIT NONE122:           IMPLICIT NONE
150:           ! EL, EM, ER are energies of the Left, Middle and Right images 
151:           DOUBLE PRECISION,INTENT(IN) :: EL,EM,ER123:           DOUBLE PRECISION,INTENT(IN) :: EL,EM,ER
152:           ! WM, WP are weights for the left and right tangents, interpolated based on the energies of the corresponding minima 
153:           ! c.f. Henkelman and Jonsson, JCP 113 (2000), equations 8-11 
154:           DOUBLE PRECISION,INTENT(OUT) :: WM,WP124:           DOUBLE PRECISION,INTENT(OUT) :: WM,WP
155: 125: 
156:           ! El -->><<-- Em -->><<-- Er126:           ! El -->><<-- Em -->><<-- Er
157:           IF (EL < EM .AND. EM < ER) THEN127:           IF (EL < EM .AND. EM < ER) THEN
158:                WP=1.0D0;     WM=0.0D0128:                WP=1.0D0;     WM=0.0D0
159:           ELSE IF (EL > EM .AND. EM > ER) THEN129:           ELSE IF (EL > EM .AND. EM > ER) THEN
160:                WP=0.0D0;     WM=1.0D0130:                WP=0.0D0;     WM=1.0D0
161:           ELSE131:           ELSE
162:                IF (EL < ER) THEN132:                IF (EL < ER) THEN
163:                     WP=MAX(ABS(EL-EM),ABS(ER-EM));  WM=MIN(ABS(EL-EM),ABS(ER-EM))133:                     WP=MAX(ABS(EL-EM),ABS(ER-EM));  WM=MIN(ABS(EL-EM),ABS(ER-EM))


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0