hdiff output

r29245/congrad.f90 2015-11-17 23:32:39.604231212 +0000 r29244/congrad.f90 2015-11-17 23:32:39.796233787 +0000
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: SUBROUTINE CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 19: SUBROUTINE CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 20: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, & 20: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
 21:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, & 21:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
 22:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, INTMINFAC, & 22:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, &
 23:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, INTSPRINGACTIVET  23:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC
 24: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG 24: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
 25: USE PORFUNCS 25: USE PORFUNCS
 26: IMPLICIT NONE 26: IMPLICIT NONE
 27:             27:            
 28: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT 28: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT
 29: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS 29: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS
 30: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1 30: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
 31: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3) 31: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
 32: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI 32: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
 33: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2) 33: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)
131: !        PRINT '(A,I6)','image number ',J1131: !        PRINT '(A,I6)','image number ',J1
132: !        PRINT '(A,6F15.10)','R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ132: !        PRINT '(A,6F15.10)','R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
133: !        PRINT '(A,6F15.10)','R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ133: !        PRINT '(A,6F15.10)','R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
134: !     ENDIF134: !     ENDIF
135:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &135:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
136:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))136:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
137:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN137:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
138: !        D12=DSQI**6138: !        D12=DSQI**6
139:          D12=DSQI139:          D12=DSQI
140: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)140: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
141:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)141:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
142:          IF (J1.EQ.2) THEN142:          IF (J1.EQ.2) THEN
143:             EEE(J1)=EEE(J1)+DUMMY143:             EEE(J1)=EEE(J1)+DUMMY
144:             REPEINT(J1)=REPEINT(J1)+DUMMY144:             REPEINT(J1)=REPEINT(J1)+DUMMY
145:             NREPINT(J1)=NREPINT(J1)+1145:             NREPINT(J1)=NREPINT(J1)+1
146:          ELSE IF (J1.LT.INTIMAGE+2) THEN146:          ELSE IF (J1.LT.INTIMAGE+2) THEN
147:             EEE(J1)=EEE(J1)+DUMMY/2.0D0147:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
148:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0148:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
149:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0149:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0
150:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0150:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0
151:             NREPINT(J1)=NREPINT(J1)+1151:             NREPINT(J1)=NREPINT(J1)+1
152:             NREPINT(J1-1)=NREPINT(J1-1)+1152:             NREPINT(J1-1)=NREPINT(J1-1)+1
153:          ELSE IF (J1.EQ.INTIMAGE+2) THEN153:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
154:             EEE(J1-1)=EEE(J1-1)+DUMMY154:             EEE(J1-1)=EEE(J1-1)+DUMMY
155:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY155:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
156:             NREPINT(J1-1)=NREPINT(J1-1)+1156:             NREPINT(J1-1)=NREPINT(J1-1)+1
157:          ENDIF157:          ENDIF
158:          EREP=EREP+DUMMY158:          EREP=EREP+DUMMY
159: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)159: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
160:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)160:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
161:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)161:          REPGRAD(1:3)=DUMMY*G1INT(1:3)
162: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &162: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
163: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)163: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
164: !164: !
165: ! Gradient contributions for image J1-1165: ! Gradient contributions for image J1-1
166: !166: !
167:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)167:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
168:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)168:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
169:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)169:          REPGRAD(1:3)=DUMMY*G2INT(1:3)
170: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &170: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
171: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)171: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
172: !172: !
173: ! Gradient contributions for image J1173: ! Gradient contributions for image J1
174: !174: !
175:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)175:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
176:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)176:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
177:       ENDIF177:       ENDIF
178:    ENDDO178:    ENDDO
179: ENDDO179: ENDDO
183: !183: !
184: ESPRING=0.0D0184: ESPRING=0.0D0
185: IF (KINT.NE.0.0D0) THEN185: IF (KINT.NE.0.0D0) THEN
186:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1186:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1
187:       NI1=NOPT*(J1-1)187:       NI1=NOPT*(J1-1)
188:       NI2=NOPT*J1188:       NI2=NOPT*J1
189: !189: !
190: !  Edge between J1 and J1+1190: !  Edge between J1 and J1+1
191: !191: !
192:       DPLUS=0.0D0192:       DPLUS=0.0D0
193: ! 
194: !  Shouldn't we sum over active atoms only here? 
195: ! 
196:       DO J2=1,NATOMS193:       DO J2=1,NATOMS
197:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN194:          DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &
198:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &195:   &                 +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &
199:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &196:   &                 +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2
200:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2 
201:          ENDIF 
202:       ENDDO197:       ENDDO
203:       DPLUS=SQRT(DPLUS)198:       DPLUS=SQRT(DPLUS)
204:       IF (DPLUS.GT.IMSEPMAX) THEN199:       IF (DPLUS.GT.IMSEPMAX) THEN
205: !        DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2200: !        DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2
206:          DUMMY=KINT*0.5D0*DPLUS**2201:          DUMMY=KINT*0.5D0*DPLUS**2
207:          ESPRING=ESPRING+DUMMY202:          ESPRING=ESPRING+DUMMY
208: !        DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS203: !        DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS
209:          DUMMY=KINT204:          DUMMY=KINT
210:          DO J2=1,NATOMS205:          DO J2=1,NATOMS
211:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN206:             SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))
212:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))207:             GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)
213:                GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)208:             GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3)
214:                GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3) 
215:             ENDIF 
216:          ENDDO209:          ENDDO
217:       ENDIF210:       ENDIF
218:    ENDDO211:    ENDDO
219: ENDIF212: ENDIF
220: !213: !
221: ! Set gradients on frozen atoms to zero.214: ! Set gradients on frozen atoms to zero.
222: !215: !
223: IF (FREEZE) THEN216: IF (FREEZE) THEN
224:    DO J1=2,INTIMAGE+1  217:    DO J1=2,INTIMAGE+1  
225:       DO J2=1,NATOMS218:       DO J2=1,NATOMS
480: 473: 
481: !474: !
482: ! This version of congrad tests for an internal minimum in the475: ! This version of congrad tests for an internal minimum in the
483: ! constraint distances as well as the repulsions.476: ! constraint distances as well as the repulsions.
484: !477: !
485: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)478: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
486: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &479: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
487:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &480:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
488:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &481:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &
489:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &482:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &
490:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, INTMINFAC, INTSPRINGACTIVET 483:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT
491: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG484: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
492: IMPLICIT NONE485: IMPLICIT NONE
493:            486:            
494: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2)487: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2)
495: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS488: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS
496: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1489: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
497: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)490: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
498: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI491: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
499: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)492: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)
500: LOGICAL NOINT, LPRINT493: LOGICAL NOINT, LPRINT
567:       CALL MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &560:       CALL MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
568:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,.FALSE.)561:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,.FALSE.)
569: !562: !
570: ! Need to include both D2 and D1 contributions if they are both outside tolerance.563: ! Need to include both D2 and D1 contributions if they are both outside tolerance.
571: ! Otherwise we get discontinuities if they are very close and swap over.564: ! Otherwise we get discontinuities if they are very close and swap over.
572: !565: !
573: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)566: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)
574: !567: !
575: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.568: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.
576: !569: !
577:       IF (LPRINT) PRINT '(A,I6,5G15.5)', &570: !     IF (LPRINT) PRINT '(A,I6,5G15.5)', &
578:   &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL571: ! &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL
579:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 572:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 
580:          DUMMY=D2-CONDISTREFLOCAL(J2)  573:          DUMMY=D2-CONDISTREFLOCAL(J2)  
581:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)574:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
582:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)575:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
583:          EEE(J1)=EEE(J1)+DUMMY576:          EEE(J1)=EEE(J1)+DUMMY
584:          CONE(J1)=CONE(J1)+DUMMY577:          CONE(J1)=CONE(J1)+DUMMY
585:          ECON=ECON      +DUMMY578:          ECON=ECON      +DUMMY
586:          IF (LPRINT) PRINT '(A,4I6,G15.5)','min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &579: !        IF (LPRINT) PRINT '(A,4I6,G15.5)','min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
587:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)580: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
588:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)581:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
589:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)582:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
590:       ENDIF583:       ENDIF
591: !584: !
592: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.585: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.
593: !586: !
594: ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1.587: ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1.
595: !588: !
596:       IF (LPRINT) PRINT '(A,I6,5G15.5)', &589: !     IF (LPRINT) PRINT '(A,I6,5G15.5)', &
597:   &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL590: ! &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL
598:       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN  591:       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN  
599:          DUMMY=D1-CONDISTREFLOCAL(J2)  592:          DUMMY=D1-CONDISTREFLOCAL(J2)  
600:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3)593:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3)
601:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)594:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
602:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN595:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
603:             PRINT '(A,2I6,2L5,G20.10)','A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &596:             PRINT '(A,2I6,2L5,G20.10)','A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
604:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY597:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
605:          ENDIF598:          ENDIF
606:          EEE(J1-1)=EEE(J1-1)+DUMMY599:          EEE(J1-1)=EEE(J1-1)+DUMMY
607:          CONE(J1-1)=CONE(J1-1)+DUMMY600:          CONE(J1-1)=CONE(J1-1)+DUMMY
608:          ECON=ECON      +DUMMY601:          ECON=ECON      +DUMMY
609:          IF (LPRINT) PRINT '(A,4I6,G15.5)','max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &602: !        IF (LPRINT) PRINT '(A,4I6,G15.5)','max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
610:   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)603: ! &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
611:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)604:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
612:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)605:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
613:       ENDIF606:       ENDIF
614:       IF ((.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN607:       IF ((.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN
615:          DUMMY=DINT-CONDISTREFLOCAL(J2)  608:          DUMMY=DINT-CONDISTREFLOCAL(J2)  
616:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)609:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)
617:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)610:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
618:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)611:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
619:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)612:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)
620:          DUMMY=INTMINFAC*INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)613:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
621:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN614:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
622: !           PRINT '(A,2I6,2L5,G20.10)','B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &615:             PRINT '(A,2I6,2L5,G20.10)','B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
623: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY616:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
624:          ENDIF617:          ENDIF
625:          ECON=ECON+DUMMY618:          ECON=ECON+DUMMY
 619: !        IF ((J2.EQ.304).OR.(J2.EQ.309).OR.(J2.EQ.311)) THEN
 620: !           PRINT '(A,G20.10)','adding INT term DUMMY=',DUMMY
 621: !        ENDIF
626:          IF (J1.EQ.2) THEN622:          IF (J1.EQ.2) THEN
627:             EEE(J1)=EEE(J1)+DUMMY623:             EEE(J1)=EEE(J1)+DUMMY
628:             CONEINT(J1)=CONEINT(J1)+DUMMY624:             CONEINT(J1)=CONEINT(J1)+DUMMY
629:             NCONINT(J1)=NCONINT(J1)+1625:             NCONINT(J1)=NCONINT(J1)+1
630:          ELSE IF (J1.LT.INTIMAGE+2) THEN626:          ELSE IF (J1.LT.INTIMAGE+2) THEN
631:             EEE(J1)=EEE(J1)+DUMMY/2.0D0627:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
632:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0628:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
633:             CONEINT(J1)=CONEINT(J1)+DUMMY/2.0D0629:             CONEINT(J1)=CONEINT(J1)+DUMMY/2.0D0
634:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY/2.0D0630:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY/2.0D0
635:             NCONINT(J1)=NCONINT(J1)+1631:             NCONINT(J1)=NCONINT(J1)+1
641:          ENDIF637:          ENDIF
642: !        PRINT '(A,4I6,G15.5)','in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &638: !        PRINT '(A,4I6,G15.5)','in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
643: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)639: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
644:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)640:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
645:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)641:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
646:       ENDIF642:       ENDIF
647:    ENDDO643:    ENDDO
648: ENDDO644: ENDDO
649: 645: 
650: ! INTCONST=INTCONSTRAINREPCUT**13646: ! INTCONST=INTCONSTRAINREPCUT**13
651: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes 
652: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes 
653: 647: 
654: DO J2=1,NNREPULSIVE648: DO J2=1,NNREPULSIVE
655: !  INTCONST=NREPCUT(J2)**13649: !  INTCONST=NREPCUT(J2)**13
656:    INTCONST=NREPCUT(J2)**3650:    INTCONST=NREPCUT(J2)**3
657:    DO J1=2,INTIMAGE+2651:    DO J1=2,INTIMAGE+2
658:       If (FREEZENODEST) THEN652:       If (FREEZENODEST) THEN
659:          IF (J1.EQ.2) THEN653:          IF (J1.EQ.2) THEN
660:             IF (IMGFREEZE(1)) CYCLE654:             IF (IMGFREEZE(1)) CYCLE
661:          ELSE IF (J1.EQ.INTIMAGE+2) THEN655:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
662:             IF (IMGFREEZE(INTIMAGE)) CYCLE656:             IF (IMGFREEZE(INTIMAGE)) CYCLE
670: !  ENDIF664: !  ENDIF
671:       NI1=NOPT*(J1-2)+3*(NREPI(J2)-1)665:       NI1=NOPT*(J1-2)+3*(NREPI(J2)-1)
672:       NI2=NOPT*(J1-1)+3*(NREPI(J2)-1)666:       NI2=NOPT*(J1-1)+3*(NREPI(J2)-1)
673:       NJ1=NOPT*(J1-2)+3*(NREPJ(J2)-1)667:       NJ1=NOPT*(J1-2)+3*(NREPJ(J2)-1)
674:       NJ2=NOPT*(J1-1)+3*(NREPJ(J2)-1)668:       NJ2=NOPT*(J1-1)+3*(NREPJ(J2)-1)
675:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)669:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
676:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)670:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
677:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)671:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
678:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)672:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
679:       IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN673:       IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
680: !        PRINT '(A,I6,A,2I6)','A repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)674:          PRINT '(A,I6,A,2I6)','A repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
681: !        PRINT '(A,I6)','image number ',J1675:          PRINT '(A,I6)','image number ',J1
682: !        PRINT '(A,6F15.10)','R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ676:          PRINT '(A,6F15.10)','R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
683: !        PRINT '(A,6F15.10)','R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ677:          PRINT '(A,6F15.10)','R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
684:       ENDIF678:       ENDIF
685:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &679:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
686:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))680:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
687:       IF ((NREPI(J2).EQ.34).AND.(NREPJ(J2).EQ.400)) THEN681: !     IF ((NREPI(J2).EQ.135).AND.(NREPJ(J2).EQ.192)) THEN
688: !        PRINT '(A,3G20.10)',' congrad2> R1AX,R1AY,R1AZ=',R1AX,R1AY,R1AZ682: !        PRINT '(A,3G20.10)',' congrad2> R1AX,R1AY,R1AZ=',R1AX,R1AY,R1AZ
689: !        PRINT '(A,3G20.10)',' congrad2> R1BX,R1BY,R1BZ=',R1BX,R1BY,R1BZ683: !        PRINT '(A,3G20.10)',' congrad2> R1BX,R1BY,R1BZ=',R1BX,R1BY,R1BZ
690: !        PRINT '(A,3G20.10)',' congrad2> R2AX,R2AY,R2AZ=',R2AX,R2AY,R2AZ684: !        PRINT '(A,3G20.10)',' congrad2> R2AX,R2AY,R2AZ=',R2AX,R2AY,R2AZ
691: !        PRINT '(A,3G20.10)',' congrad2> R2BX,R2BY,R2BZ=',R2BX,R2BY,R2BZ685: !        PRINT '(A,3G20.10)',' congrad2> R2BX,R2BY,R2BZ=',R2BX,R2BY,R2BZ
692: !        PRINT '(A,I6,A,2I6)',' congrad2> J1=',J1,' edge between images: ',J1-1,J1686: !        PRINT '(A,I6,A,2I6)',' congrad2> J1=',J1,' edge between images: ',J1-1,J1
693: !        PRINT '(A,L5,3G20.10)',' congrad2> NOINT,D2,D1,DINT=',NOINT,D2,D1,DINT687: !        PRINT '(A,L5,3G20.10)',' congrad2> NOINT,D2,D1,DINT=',NOINT,D2,D1,DINT
694:       ENDIF688: !     ENDIF
695:       DUMMY=0.0D0 689:       DUMMY=0.0D0 
696: !690: !
697: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.691: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
698: !692: !
699: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1693: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
700:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1694:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
701: !        D12=DSQ2**6695: !        D12=DSQ2**6
702:          D12=DSQ2696:          D12=DSQ2
703: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)697: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
704: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)698: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
705:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)699:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
706: !        WRITE(*, '(A,5G20.10)') 'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST 
707:          EEE(J1)=EEE(J1)+DUMMY700:          EEE(J1)=EEE(J1)+DUMMY
708:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN701:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
709:             PRINT '(A,2I6,2L5,G20.10)','R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &702:             PRINT '(A,2I6,2L5,G20.10)','R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
710:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY703:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
711:          ENDIF704:          ENDIF
712:          REPE(J1)=REPE(J1)+DUMMY705:          REPE(J1)=REPE(J1)+DUMMY
713:          EREP=EREP+DUMMY706:          EREP=EREP+DUMMY
714: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)707: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
715:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)708:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
716:          REPGRAD(1:3)=DUMMY*G2(1:3)709:          REPGRAD(1:3)=DUMMY*G2(1:3)
745:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)738:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
746:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)739:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
747:       ENDIF740:       ENDIF
748:       DUMMY=0.0D0741:       DUMMY=0.0D0
749: !     IF ((.NOT.NOINT).AND.(DINT.LT.INTCONSTRAINREPCUT)) THEN742: !     IF ((.NOT.NOINT).AND.(DINT.LT.INTCONSTRAINREPCUT)) THEN
750:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN743:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
751: !        D12=DSQI**6744: !        D12=DSQI**6
752:          D12=DSQI745:          D12=DSQI
753: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*INTCONSTRAINREPCUT)/INTCONST)746: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
754: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)747: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
755:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)748:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
756:          EREP=EREP+DUMMY749:          EREP=EREP+DUMMY
757: !        IF (DUMMY.GT.1.0D7) PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &750: !        IF (DUMMY.GT.1.0D7) PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &
758: ! &                                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY751: ! &                                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY
 752: !        IF (((NREPI(J2).EQ.143).AND.(NREPJ(J2).EQ.191)).OR.  &
 753: ! &          ((NREPJ(J2).EQ.143).AND.(NREPI(J2).EQ.191))) THEN
 754: !            PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &
 755: ! &                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY
759: !        ENDIF756: !        ENDIF
760:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN757:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
761:             PRINT '(A,2I6,2L5,G20.10)','R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &758:             PRINT '(A,2I6,2L5,G20.10)','R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
762:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY759:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
763:          ENDIF760:          ENDIF
764:          IF (J1.EQ.2) THEN761:          IF (J1.EQ.2) THEN
765:             EEE(J1)=EEE(J1)+DUMMY762:             EEE(J1)=EEE(J1)+DUMMY
766:             REPEINT(J1)=REPEINT(J1)+DUMMY763:             REPEINT(J1)=REPEINT(J1)+DUMMY
767:             NREPINT(J1)=NREPINT(J1)+1764:             NREPINT(J1)=NREPINT(J1)+1
768:          ELSE IF (J1.LT.INTIMAGE+2) THEN765:          ELSE IF (J1.LT.INTIMAGE+2) THEN
772:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0769:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0
773:             NREPINT(J1)=NREPINT(J1)+1770:             NREPINT(J1)=NREPINT(J1)+1
774:             NREPINT(J1-1)=NREPINT(J1-1)+1771:             NREPINT(J1-1)=NREPINT(J1-1)+1
775:          ELSE IF (J1.EQ.INTIMAGE+2) THEN772:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
776:             EEE(J1-1)=EEE(J1-1)+DUMMY773:             EEE(J1-1)=EEE(J1-1)+DUMMY
777:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY774:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
778:             NREPINT(J1-1)=NREPINT(J1-1)+1775:             NREPINT(J1-1)=NREPINT(J1-1)+1
779:          ENDIF776:          ENDIF
780: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)777: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
781:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)778:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
782:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)779:          REPGRAD(1:3)=DUMMY*G1INT(1:3)
783: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &780: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
784: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)781: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
785:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)782:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
786:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)783:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
787:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)784:          REPGRAD(1:3)=DUMMY*G2INT(1:3)
788: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &785: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
789: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)786: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
790:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)787:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
791:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)788:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
792:       ENDIF789:       ENDIF
793:    ENDDO790:    ENDDO
794: ENDDO791: ENDDO
795: !792: !
796: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise793: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise
797: ! energy terms between images except for the end points.794: ! energy terms between images except for the end points.
798: !795: !
799: ESPRING=0.0D0796: ESPRING=0.0D0
800: IF (KINT.NE.0.0D0) THEN797: IF (KINT.NE.0.0D0) THEN
801:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1798:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1
802:       NI1=NOPT*(J1-1)799:       NI1=NOPT*(J1-1)
803:       NI2=NOPT*J1800:       NI2=NOPT*J1
804: !801: !
805: !  Edge between J1 and J1+1802: !  Edge between J1 and J1+1
806: !  Shouldn't we sum over active atoms only here? 
807: !803: !
808:       DPLUS=0.0D0804:       DPLUS=0.0D0
809:       DO J2=1,NATOMS805:       DO J2=1,NATOMS
810:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN806:          DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &
811:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &807:   &                 +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &
812:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &808:   &                 +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2
813:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2 
814:          ENDIF 
815:       ENDDO809:       ENDDO
816:       DPLUS=SQRT(DPLUS)810:       DPLUS=SQRT(DPLUS)
817:       IF (DPLUS.GT.IMSEPMAX) THEN811:       IF (DPLUS.GT.IMSEPMAX) THEN
818:          DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2812:          DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2
819: !        DUMMY=KINT*0.5D0*DPLUS**2813: !        DUMMY=KINT*0.5D0*DPLUS**2
820:          ESPRING=ESPRING+DUMMY814:          ESPRING=ESPRING+DUMMY
821: !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING 
822:          DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS815:          DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS
823: !        DUMMY=KINT816: !        DUMMY=KINT
824:          DO J2=1,NATOMS817:          DO J2=1,NATOMS
825:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN818:             SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))
826:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))819:             GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)
827:                GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)820:             GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3)
828:                GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3) 
829:             ENDIF 
830:          ENDDO821:          ENDDO
831:       ENDIF822:       ENDIF
832:    ENDDO823:    ENDDO
833: ENDIF824: ENDIF
834:          IF (PRINTE) THEN825:          IF (PRINTE) THEN
835:             PRINT '(A,G20.10)','ESPRING=',ESPRING826:             PRINT '(A,G20.10)','ESPRING=',ESPRING
836:          ENDIF827:          ENDIF
837: !828: !
838: ! Set gradients on frozen atoms to zero.829: ! Set gradients on frozen atoms to zero.
839: !830: !
869:    GGG(1:NOPT)=0.0D0860:    GGG(1:NOPT)=0.0D0
870:    GGG((INTIMAGE+1)*NOPT+1:(INTIMAGE+2)*NOPT)=0.0D0861:    GGG((INTIMAGE+1)*NOPT+1:(INTIMAGE+2)*NOPT)=0.0D0
871: ENDIF862: ENDIF
872: RMS=0.0D0863: RMS=0.0D0
873: RMSIMAGE(1:INTIMAGE+2)=0.0D0864: RMSIMAGE(1:INTIMAGE+2)=0.0D0
874: DO J1=2,INTIMAGE+1865: DO J1=2,INTIMAGE+1
875:    DO J2=1,NOPT866:    DO J2=1,NOPT
876:       RMSIMAGE(J1)=RMSIMAGE(J1)+GGG(NOPT*(J1-1)+J2)**2867:       RMSIMAGE(J1)=RMSIMAGE(J1)+GGG(NOPT*(J1-1)+J2)**2
877:    ENDDO868:    ENDDO
878:    RMS=RMS+RMSIMAGE(J1)869:    RMS=RMS+RMSIMAGE(J1)
879:    IF (LPRINT) PRINT '(A,I6,2G20.10,L5)',' congrad2> J1,EEE,RMSIMAGE=',J1,EEE(J1),RMSIMAGE(J1)870:    IF (LPRINT) PRINT '(A,I6,2G20.10,L5)',' congrad2> J1,EEE,RMSIMAGE,freeze=', &
 871:   &                                                  J1,EEE(J1),RMSIMAGE(J1),IMGFREEZE(J1)
880: ENDDO872: ENDDO
881: IF (INTIMAGE.NE.0) THEN873: IF (INTIMAGE.NE.0) THEN
882:    RMS=SQRT(RMS/(NOPT*INTIMAGE))874:    RMS=SQRT(RMS/(NOPT*INTIMAGE))
883: ENDIF875: ENDIF
884: !876: !
885: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,877: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,
886: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)878: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)
887: !879: !
888: ETOTAL=0.0D0880: ETOTAL=0.0D0
889: MAXINT=-1.0D100881: MAXINT=-1.0D100


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0