hdiff output

r33422/intlbfgs.f90 2017-10-27 12:32:07.812129667 +0100 r33421/intlbfgs.f90 2017-10-27 12:32:08.164133616 +0100
 11: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 11: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 12: !   GNU General Public License for more details. 12: !   GNU General Public License for more details.
 13: ! 13: !
 14: !   You should have received a copy of the GNU General Public License 14: !   You should have received a copy of the GNU General Public License
 15: !   along with this program; if not, write to the Free Software 15: !   along with this program; if not, write to the Free Software
 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 17: ! 17: !
 18: SUBROUTINE INTLBFGS(QSTART,QFINISH) 18: SUBROUTINE INTLBFGS(QSTART,QFINISH)
 19: USE PORFUNCS 19: USE PORFUNCS
 20: USE KEY, ONLY : FREEZENODEST, FREEZETOL, MAXBFGS, CONVR, ATOMSTORES, & 20: USE KEY, ONLY : FREEZENODEST, FREEZETOL, MAXBFGS, CONVR, ATOMSTORES, &
 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, INTMUPDATE, INTDGUESS, & 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, MUPDATE, INTDGUESS, &
 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, & 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, &
 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, & 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, &
 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, & 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, &
 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, & 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, &
 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, & 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, &
 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, & 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, &
 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, & 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, &
 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, & 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, &
 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, & 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, &
 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, & 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, &
 49: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY 49: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY
 50: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF 50: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF
 51: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2 51: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2
 52: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE, NBONDED(NATOMS), BONDEDLIST(NATOMS,6), NBOND 52: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE, NBONDED(NATOMS), BONDEDLIST(NATOMS,6), NBOND
 53: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS), ACID 53: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS), ACID
 54: LOGICAL CHIRALSR, CHIRALSRP  54: LOGICAL CHIRALSR, CHIRALSRP 
 55: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS) 55: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS)
 56: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, & 56: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, &
 57:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), & 57:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), &
 58:   &                 BESTWORST, WORST, COORDSA(3*NATOMS), COORDSB(3*NATOMS), COORDSC(3*NATOMS) 58:   &                 BESTWORST, WORST, COORDSA(3*NATOMS), COORDSB(3*NATOMS), COORDSC(3*NATOMS)
 59: DOUBLE PRECISION, DIMENSION(INTMUPDATE)     :: RHO1,ALPHA 59: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA
 60: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS) 60: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS)
 61: LOGICAL SWITCHED, AABACK(NATOMS), BACKDONE 61: LOGICAL SWITCHED, AABACK(NATOMS), BACKDONE
 62: DOUBLE PRECISION, POINTER :: X(:), G(:) 62: DOUBLE PRECISION, POINTER :: X(:), G(:)
 63: ! 63: !
 64: ! efk: for freezenodes 64: ! efk: for freezenodes
 65: ! 65: !
 66: DOUBLE PRECISION :: TESTG, TOTGNORM 66: DOUBLE PRECISION :: TESTG, TOTGNORM
 67: INTEGER :: IM 67: INTEGER :: IM
 68: ! 68: !
 69: ! Dimensions involving INTIMAGE 69: ! Dimensions involving INTIMAGE
115:    GOTO 653115:    GOTO 653
116: 654 CONTINUE116: 654 CONTINUE
117:    INTIMAGE=INTIMAGE-2117:    INTIMAGE=INTIMAGE-2
118:    WRITE(*,'(A,I10,A)') 'intlbfgs> Rereading ',INTIMAGE,' frames'118:    WRITE(*,'(A,I10,A)') 'intlbfgs> Rereading ',INTIMAGE,' frames'
119:    CLOSE(LUNIT)119:    CLOSE(LUNIT)
120: ENDIF120: ENDIF
121: 121: 
122: ALLOCATE(TRUEEE(INTIMAGE+2), &122: ALLOCATE(TRUEEE(INTIMAGE+2), &
123:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &123:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &
124:   &      GTMP(3*NATOMS*INTIMAGE), &124:   &      GTMP(3*NATOMS*INTIMAGE), &
125:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:INTMUPDATE,(3*NATOMS)*INTIMAGE), &125:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE), &
126:   &      GDIF(0:INTMUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &126:   &      GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &
127:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &127:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &
128:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))128:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))
129: 129: 
130: SWITCHED=.FALSE.130: SWITCHED=.FALSE.
131: INTIMAGESAVE=INTIMAGE131: INTIMAGESAVE=INTIMAGE
132: NBACKTRACK=1132: NBACKTRACK=1
133: CALL MYCPU_TIME(STIME,.FALSE.)133: CALL MYCPU_TIME(STIME,.FALSE.)
134: WRITE(*,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1134: WRITE(*,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1
135: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',INTMUPDATE,' maximum step size=',MAXBFGS135: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',MUPDATE,' maximum step size=',MAXBFGS
136: ADDATOM=.FALSE.136: ADDATOM=.FALSE.
137: NFAIL=0137: NFAIL=0
138: IMGFREEZE(1:INTIMAGE)=.FALSE.138: IMGFREEZE(1:INTIMAGE)=.FALSE.
139: D=(3*NATOMS)*INTIMAGE139: D=(3*NATOMS)*INTIMAGE
140: U=INTMUPDATE140: U=MUPDATE
141: NITERDONE=1141: NITERDONE=1
142: NITERUSE=1142: NITERUSE=1
143: NQDONE=0143: NQDONE=0
144: 144: 
145: !145: !
146: ! Read AMBER topology and count bonded neighbours for each atom.146: ! Read AMBER topology and count bonded neighbours for each atom.
147: ! Assume that we won't have more than 6 bonds!147: ! Assume that we won't have more than 6 bonds!
148: !148: !
149: IF (QCIAMBERT) THEN149: IF (QCIAMBERT) THEN
150:    CALL TOPOLOGY_READER(NBOND)150:    CALL TOPOLOGY_READER(NBOND)
1002:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))1002:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
1003:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))1003:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
1004:             DEALLOCATE(XYZ)1004:             DEALLOCATE(XYZ)
1005:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))1005:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))
1006:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)1006:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)
1007:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &1007:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &
1008:   &                                               + DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1)))/2.0D01008:   &                                               + DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1)))/2.0D0
1009:             XYZ(3*NATOMS*(JMAX+1)+1:3*NATOMS*(INTIMAGE+3))=DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+2))1009:             XYZ(3*NATOMS*(JMAX+1)+1:3*NATOMS*(INTIMAGE+3))=DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+2))
1010: !1010: !
1011: ! Save step-taking memories in SEARCHSTEP and GDIF.1011: ! Save step-taking memories in SEARCHSTEP and GDIF.
1012: ! These arrays run from 0 to INTMUPDATE over memories and1012: ! These arrays run from 0 to MUPDATE over memories and
1013: ! 1:(3*NATOMS)*INTIMAGE over only the variable images.1013: ! 1:(3*NATOMS)*INTIMAGE over only the variable images.
1014: !1014: !
1015:             DEALLOCATE(DPTMP)1015:             DEALLOCATE(DPTMP)
1016:             ALLOCATE(D2TMP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE))1016:             ALLOCATE(D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE))
1017:             D2TMP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)1017:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
1018:             DEALLOCATE(SEARCHSTEP)1018:             DEALLOCATE(SEARCHSTEP)
1019:             ALLOCATE(SEARCHSTEP(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))1019:             ALLOCATE(SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))
1020:             DO J1=0,INTMUPDATE1020:             DO J1=0,MUPDATE
1021:                IF (JMAX.GT.1) SEARCHSTEP(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))1021:                IF (JMAX.GT.1) SEARCHSTEP(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))
1022:                IF (JMAX.LT.INTIMAGE+1) SEARCHSTEP(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &1022:                IF (JMAX.LT.INTIMAGE+1) SEARCHSTEP(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &
1023:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)1023:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)
1024:                SEARCHSTEP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &1024:                SEARCHSTEP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &
1025:   &                             D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))1025:   &                             D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))
1026:             ENDDO1026:             ENDDO
1027: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1027: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1028:             SEARCHSTEP(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE+1))=0.0D01028:             SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1))=0.0D0
1029: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1029: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1030:             D2TMP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)1030:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
1031:             DEALLOCATE(GDIF)1031:             DEALLOCATE(GDIF)
1032:             ALLOCATE(GDIF(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))1032:             ALLOCATE(GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))
1033:             DO J1=0,INTMUPDATE1033:             DO J1=0,MUPDATE
1034:                IF (JMAX.GT.1) GDIF(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))1034:                IF (JMAX.GT.1) GDIF(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))
1035:                IF (JMAX.LT.INTIMAGE+1) GDIF(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &1035:                IF (JMAX.LT.INTIMAGE+1) GDIF(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &
1036:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)1036:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)
1037:                GDIF(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &1037:                GDIF(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &
1038:   &                       D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))1038:   &                       D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))
1039:             ENDDO1039:             ENDDO
1040: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1040: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1041:             GDIF(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE+1))=0.0D01041:             GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1))=0.0D0
1042: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1042: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1043:             DEALLOCATE(D2TMP)1043:             DEALLOCATE(D2TMP)
1044: 1044: 
1045:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &1045:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &
1046:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)1046:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)
1047:             ALLOCATE(TRUEEE(INTIMAGE+3), &1047:             ALLOCATE(TRUEEE(INTIMAGE+3), &
1048:   &                  EEETMP(INTIMAGE+3), MYGTMP(3*NATOMS*(INTIMAGE+1)), &1048:   &                  EEETMP(INTIMAGE+3), MYGTMP(3*NATOMS*(INTIMAGE+1)), &
1049:   &                  GTMP(3*NATOMS*(INTIMAGE+1)), &1049:   &                  GTMP(3*NATOMS*(INTIMAGE+1)), &
1050:   &                  DIAG(3*NATOMS*(INTIMAGE+1)), STP(3*NATOMS*(INTIMAGE+1)), &1050:   &                  DIAG(3*NATOMS*(INTIMAGE+1)), STP(3*NATOMS*(INTIMAGE+1)), &
1051:   &                  GLAST((3*NATOMS)*(INTIMAGE+1)), &1051:   &                  GLAST((3*NATOMS)*(INTIMAGE+1)), &
1087:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))1087:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
1088:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))1088:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
1089:             DEALLOCATE(XYZ)1089:             DEALLOCATE(XYZ)
1090:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+1)))1090:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+1)))
1091:             XYZ(1:3*NATOMS*(JMIN-1))=DPTMP(1:3*NATOMS*(JMIN-1))1091:             XYZ(1:3*NATOMS*(JMIN-1))=DPTMP(1:3*NATOMS*(JMIN-1))
1092:             XYZ(3*NATOMS*(JMIN-1)+1:3*NATOMS*(INTIMAGE+1))=DPTMP(3*NATOMS*JMIN+1:3*NATOMS*(INTIMAGE+2))1092:             XYZ(3*NATOMS*(JMIN-1)+1:3*NATOMS*(INTIMAGE+1))=DPTMP(3*NATOMS*JMIN+1:3*NATOMS*(INTIMAGE+2))
1093: 1093: 
1094:             DEALLOCATE(DPTMP)1094:             DEALLOCATE(DPTMP)
1095: !1095: !
1096: ! Save step-taking memories in SEARCHSTEP and GDIF.1096: ! Save step-taking memories in SEARCHSTEP and GDIF.
1097: ! These arrays run from 0 to INTMUPDATE over memories and1097: ! These arrays run from 0 to MUPDATE over memories and
1098: ! 1:(3*NATOMS)*INTIMAGE over only the variable images.1098: ! 1:(3*NATOMS)*INTIMAGE over only the variable images.
1099: !1099: !
1100:             ALLOCATE(D2TMP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE))1100:             ALLOCATE(D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE))
1101:             D2TMP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)1101:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
1102:             DEALLOCATE(SEARCHSTEP)1102:             DEALLOCATE(SEARCHSTEP)
1103:             ALLOCATE(SEARCHSTEP(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))1103:             ALLOCATE(SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))
1104:             DO J1=0,INTMUPDATE1104:             DO J1=0,MUPDATE
1105:                SEARCHSTEP(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))1105:                SEARCHSTEP(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))
1106:                SEARCHSTEP(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &1106:                SEARCHSTEP(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &
1107:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)1107:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)
1108:             ENDDO1108:             ENDDO
1109: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1109: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1110:             SEARCHSTEP(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE-1))=0.0D01110:             SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1))=0.0D0
1111: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1111: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1112:             D2TMP(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:INTMUPDATE,1:(3*NATOMS)*INTIMAGE)1112:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
1113:             DEALLOCATE(GDIF)1113:             DEALLOCATE(GDIF)
1114:             ALLOCATE(GDIF(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))1114:             ALLOCATE(GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))
1115:             DO J1=0,INTMUPDATE1115:             DO J1=0,MUPDATE
1116:                GDIF(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))1116:                GDIF(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))
1117:                GDIF(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &1117:                GDIF(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &
1118:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)1118:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)
1119:             ENDDO1119:             ENDDO
1120: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1120: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1121:             GDIF(0:INTMUPDATE,1:(3*NATOMS)*(INTIMAGE-1))=0.0D01121:             GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1))=0.0D0
1122: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1122: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1123:             DEALLOCATE(D2TMP)1123:             DEALLOCATE(D2TMP)
1124: 1124: 
1125:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &1125:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &
1126:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)1126:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)
1127:             ALLOCATE(TRUEEE(INTIMAGE+1),&1127:             ALLOCATE(TRUEEE(INTIMAGE+1),&
1128:   &                  EEETMP(INTIMAGE+1), MYGTMP(3*NATOMS*(INTIMAGE-1)), &1128:   &                  EEETMP(INTIMAGE+1), MYGTMP(3*NATOMS*(INTIMAGE-1)), &
1129:   &                  GTMP(3*NATOMS*(INTIMAGE-1)), &1129:   &                  GTMP(3*NATOMS*(INTIMAGE-1)), &
1130:   &                  DIAG(3*NATOMS*(INTIMAGE-1)), STP(3*NATOMS*(INTIMAGE-1)), &1130:   &                  DIAG(3*NATOMS*(INTIMAGE-1)), STP(3*NATOMS*(INTIMAGE-1)), &
1131:   &                  GLAST((3*NATOMS)*(INTIMAGE-1)), &1131:   &                  GLAST((3*NATOMS)*(INTIMAGE-1)), &
1506:       EXIT1506:       EXIT
1507:    ENDIF1507:    ENDIF
1508:    777 CONTINUE1508:    777 CONTINUE
1509: !1509: !
1510: ! Compute the new step and gradient change1510: ! Compute the new step and gradient change
1511: !1511: !
1512:    NPT=POINT*D1512:    NPT=POINT*D
1513:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)1513:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)
1514:    GDIF(POINT,:)=G-GTMP1514:    GDIF(POINT,:)=G-GTMP
1515:    1515:    
1516:    POINT=POINT+1; IF (POINT==INTMUPDATE) POINT=01516:    POINT=POINT+1; IF (POINT==MUPDATE) POINT=0
1517: 1517: 
1518:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)1518:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)
1519:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1519:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1520: 1520: 
1521:    NITERDONE=NITERDONE+11521:    NITERDONE=NITERDONE+1
1522:    NITERUSE=NITERUSE+11522:    NITERUSE=NITERUSE+1
1523: 1523: 
1524:    IF (NITERDONE.GT.NSTEPSMAX) EXIT1524:    IF (NITERDONE.GT.NSTEPSMAX) EXIT
1525:    IF (NACTIVE.EQ.NATOMS) THEN1525:    IF (NACTIVE.EQ.NATOMS) THEN
1526:       IF (.NOT.SWITCHED) THEN1526:       IF (.NOT.SWITCHED) THEN
2976: ! IF (DEBUG) WRITE(*,'(A,F15.5)') ' checkperc> Final constraint tolerance parameter ',LINTCONSTRAINTTOL2976: ! IF (DEBUG) WRITE(*,'(A,F15.5)') ' checkperc> Final constraint tolerance parameter ',LINTCONSTRAINTTOL
2977: 2977: 
2978: ! WRITE(*,'(A,I6,3(A,F15.5))') ' checkperc> Total distance constraints=',NCONSTRAINT, &2978: ! WRITE(*,'(A,I6,3(A,F15.5))') ' checkperc> Total distance constraints=',NCONSTRAINT, &
2979: !   &                    ' shortest=',MINCONDIST,' longest=',MAXCONDIST,' tolerance=',LINTCONSTRAINTTOL2979: !   &                    ' shortest=',MINCONDIST,' longest=',MAXCONDIST,' tolerance=',LINTCONSTRAINTTOL
2980: 2980: 
2981: CALLED=.TRUE.2981: CALLED=.TRUE.
2982: 2982: 
2983: END SUBROUTINE CHECKPERC2983: END SUBROUTINE CHECKPERC
2984: 2984: 
2985: SUBROUTINE MAKESTEP(NITERDONE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA)2985: SUBROUTINE MAKESTEP(NITERDONE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA)
2986: USE KEY, ONLY : INTMUPDATE, INTDGUESS2986: USE KEY, ONLY : MUPDATE, INTDGUESS
2987: USE COMMONS, ONLY: NATOMS2987: USE COMMONS, ONLY: NATOMS
2988: IMPLICIT NONE2988: IMPLICIT NONE
2989: INTEGER NITERDONE, POINT, BOUND, NPT, D, CP, INTIMAGE, I2989: INTEGER NITERDONE, POINT, BOUND, NPT, D, CP, INTIMAGE, I
2990: DOUBLE PRECISION DIAG(3*NATOMS*INTIMAGE),SEARCHSTEP(0:INTMUPDATE,(3*NATOMS)*INTIMAGE),G((3*NATOMS)*INTIMAGE), &2990: DOUBLE PRECISION DIAG(3*NATOMS*INTIMAGE),SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE),G((3*NATOMS)*INTIMAGE), &
2991:   &  GTMP(3*NATOMS*INTIMAGE), GNORM, STP(3*NATOMS*INTIMAGE), YS, GDIF(0:INTMUPDATE,(3*NATOMS)*INTIMAGE), YY, &2991:   &  GTMP(3*NATOMS*INTIMAGE), GNORM, STP(3*NATOMS*INTIMAGE), YS, GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE), YY, &
2992:   &  SQ, YR, BETA2992:   &  SQ, YR, BETA
2993: DOUBLE PRECISION, DIMENSION(INTMUPDATE)     :: RHO1,ALPHA2993: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA
2994: LOGICAL CHANGEIMAGE2994: LOGICAL CHANGEIMAGE
2995: SAVE2995: SAVE
2996: 2996: 
2997: MAIN: IF (NITERDONE==1) THEN2997: MAIN: IF (NITERDONE==1) THEN
2998:      POINT = 02998:      POINT = 0
2999:      DIAG(1:D)=INTDGUESS2999:      DIAG(1:D)=INTDGUESS
3000:      SEARCHSTEP(0,1:D)= -G(1:D)*INTDGUESS            ! NR STEP FOR DIAGONAL INVERSE HESSIAN3000:      SEARCHSTEP(0,1:D)= -G(1:D)*INTDGUESS            ! NR STEP FOR DIAGONAL INVERSE HESSIAN
3001:      GTMP(1:D)        = SEARCHSTEP(0,1:D)3001:      GTMP(1:D)        = SEARCHSTEP(0,1:D)
3002:      GNORM            = MAX(SQRT(DOT_PRODUCT(G(1:D),G(1:D))),1.0D-100)3002:      GNORM            = MAX(SQRT(DOT_PRODUCT(G(1:D),G(1:D))),1.0D-100)
3003:      STP(1:D)         = MIN(1.0D0/GNORM, GNORM) ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS3003:      STP(1:D)         = MIN(1.0D0/GNORM, GNORM) ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS
3004: ELSE MAIN3004: ELSE MAIN
3005:      BOUND=NITERDONE-13005:      BOUND=NITERDONE-1
3006:      IF (NITERDONE.GT.INTMUPDATE) BOUND=INTMUPDATE3006:      IF (NITERDONE.GT.MUPDATE) BOUND=MUPDATE
3007:      YS=DOT_PRODUCT( GDIF(NPT/D,:), SEARCHSTEP(NPT/D,:)  )3007:      YS=DOT_PRODUCT( GDIF(NPT/D,:), SEARCHSTEP(NPT/D,:)  )
3008:      IF (YS==0.0D0) YS=1.0D03008:      IF (YS==0.0D0) YS=1.0D0
3009:     3009:     
3010: ! Update estimate of diagonal inverse Hessian elements.3010: ! Update estimate of diagonal inverse Hessian elements.
3011: ! We divide by both YS and YY at different points, so they had better not be zero!3011: ! We divide by both YS and YY at different points, so they had better not be zero!
3012: 3012: 
3013:      YY=DOT_PRODUCT( GDIF(NPT/D,:) , GDIF(NPT/D,:) )3013:      YY=DOT_PRODUCT( GDIF(NPT/D,:) , GDIF(NPT/D,:) )
3014:      IF (YY==0.0D0) YY=1.0D03014:      IF (YY==0.0D0) YY=1.0D0
3015: !    DIAG = ABS(YS/YY)3015: !    DIAG = ABS(YS/YY)
3016:      DIAG(1) = YS/YY3016:      DIAG(1) = YS/YY
3017:       3017:       
3018: ! COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, 3018: ! COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, 
3019: ! "Updating quasi-Newton matrices with limited storage",3019: ! "Updating quasi-Newton matrices with limited storage",
3020: ! Mathematics of Computation, Vol.35, No.151, pp. 773-7823020: ! Mathematics of Computation, Vol.35, No.151, pp. 773-782
3021: 3021: 
3022:      CP= POINT; IF (POINT==0) CP = INTMUPDATE3022:      CP= POINT; IF (POINT==0) CP = MUPDATE
3023:      RHO1(CP)=1.0D0/YS3023:      RHO1(CP)=1.0D0/YS
3024:      GTMP(1:D) = -G(1:D)3024:      GTMP(1:D) = -G(1:D)
3025:      CP= POINT 3025:      CP= POINT 
3026:                    3026:                    
3027:      DO I= 1,BOUND 3027:      DO I= 1,BOUND 
3028:           CP = CP - 1; IF (CP == -1) CP = INTMUPDATE - 13028:           CP = CP - 1; IF (CP == -1) CP = MUPDATE - 1
3029:           SQ= DOT_PRODUCT( SEARCHSTEP(CP,1:D),GTMP(1:D) )3029:           SQ= DOT_PRODUCT( SEARCHSTEP(CP,1:D),GTMP(1:D) )
3030:           ALPHA(CP+1) = RHO1(CP+1) * SQ3030:           ALPHA(CP+1) = RHO1(CP+1) * SQ
3031:           GTMP(1:D)        = -ALPHA(CP+1)*GDIF(CP,1:D) + GTMP(1:D)3031:           GTMP(1:D)        = -ALPHA(CP+1)*GDIF(CP,1:D) + GTMP(1:D)
3032:      ENDDO3032:      ENDDO
3033:               3033:               
3034:      GTMP(1:D)=DIAG(1)*GTMP(1:D)3034:      GTMP(1:D)=DIAG(1)*GTMP(1:D)
3035: 3035: 
3036:      DO I=1,BOUND3036:      DO I=1,BOUND
3037:           YR= DOT_PRODUCT( GDIF(CP,1:D) , GTMP )3037:           YR= DOT_PRODUCT( GDIF(CP,1:D) , GTMP )
3038:           BETA= RHO1(CP+1)*YR3038:           BETA= RHO1(CP+1)*YR
3039:           BETA= ALPHA(CP+1)-BETA3039:           BETA= ALPHA(CP+1)-BETA
3040: !         WRITE(*,'(A,I8,4G20.10)') 'makestep> I,YR,BETA,RHO1,ALPHA=',I,YR,BETA,RHO1(CP+1),ALPHA(CP+1)3040: !         WRITE(*,'(A,I8,4G20.10)') 'makestep> I,YR,BETA,RHO1,ALPHA=',I,YR,BETA,RHO1(CP+1),ALPHA(CP+1)
3041:           GTMP(1:D) = BETA*SEARCHSTEP(CP,1:D) + GTMP(1:D)3041:           GTMP(1:D) = BETA*SEARCHSTEP(CP,1:D) + GTMP(1:D)
3042:           CP=CP+13042:           CP=CP+1
3043: !         IF (CP==M) CP=03043: !         IF (CP==M) CP=0
3044:           IF (CP==INTMUPDATE) CP=03044:           IF (CP==MUPDATE) CP=0
3045:      ENDDO3045:      ENDDO
3046:               3046:               
3047:      STP(1:D) = 1.0D03047:      STP(1:D) = 1.0D0
3048: ENDIF MAIN3048: ENDIF MAIN
3049: 3049: 
3050: !  Store the new search direction3050: !  Store the new search direction
3051: IF (NITERDONE.GT.1) SEARCHSTEP(POINT,1:D)=GTMP(1:D)3051: IF (NITERDONE.GT.1) SEARCHSTEP(POINT,1:D)=GTMP(1:D)
3052: 3052: 
3053: END SUBROUTINE MAKESTEP3053: END SUBROUTINE MAKESTEP
3054: 3054: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0