hdiff output

r30598/key.f90 2016-06-15 11:30:10.042957973 +0100 r30597/key.f90 2016-06-15 11:30:11.062971472 +0100
 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, & 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, &
 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, & 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, &
 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 55:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, SLERPT 55:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT
 56:  56: 
 57: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 57: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 58:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 58:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 59:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 59:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 60:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 60:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 61:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 61:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
 62:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads 62:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads
 63:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs 63:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs
 64:  64: 
 65: ! hk286 > generalised THOMSON problem 65: ! hk286 > generalised THOMSON problem


r30598/keywords.f 2016-06-15 11:30:10.242960620 +0100 r30597/keywords.f 2016-06-15 11:30:11.286974435 +0100
167:          NTYPEA=NATOMS167:          NTYPEA=NATOMS
168:          !168:          !
169:          ! hk286 - initialise to stationary frame169:          ! hk286 - initialise to stationary frame
170:          RBAANORMALMODET = .FALSE.170:          RBAANORMALMODET = .FALSE.
171:          ! generalised rigid body171:          ! generalised rigid body
172:          ATOMRIGIDCOORDT = .TRUE.172:          ATOMRIGIDCOORDT = .TRUE.
173:          RIGIDINIT = .FALSE.173:          RIGIDINIT = .FALSE.
174:          AACONVERGENCET = .FALSE.174:          AACONVERGENCET = .FALSE.
175: 175: 
176:          ALIGNRBST = .FALSE.176:          ALIGNRBST = .FALSE.
177:          SLERPT = .FALSE. 
178: 177: 
179:          AVOID_COLLISIONS = .FALSE.178:          AVOID_COLLISIONS = .FALSE.
180:          COLL_TOL = 0.0D0179:          COLL_TOL = 0.0D0
181: 180: 
182:          ! Multiple potential scheme181:          ! Multiple potential scheme
183:          MULTIPOTT= .FALSE.182:          MULTIPOTT= .FALSE.
184: 183: 
185:          ! Thomson problem184:          ! Thomson problem
186:          GTHOMSONT = .FALSE.185:          GTHOMSONT = .FALSE.
187:          GTHOMPOT = 1186:          GTHOMPOT = 1
5631:                CALL READF(PARAM3)5630:                CALL READF(PARAM3)
5632:             ENDIF5631:             ENDIF
5633:             IF (NITEMS.GT.4) THEN5632:             IF (NITEMS.GT.4) THEN
5634:                CALL READF(PARAM4)5633:                CALL READF(PARAM4)
5635:             ENDIF5634:             ENDIF
5636:          ELSE IF (WORD.EQ.'SIO2C6') THEN5635:          ELSE IF (WORD.EQ.'SIO2C6') THEN
5637:             SIO2C6T=.TRUE.5636:             SIO2C6T=.TRUE.
5638:             CALL READF(C6PP)5637:             CALL READF(C6PP)
5639:             CALL READF(C6MM)5638:             CALL READF(C6MM)
5640:             CALL READF(C6PM)5639:             CALL READF(C6PM)
5641:  
5642:          ! Activate the iSLERP interpolation scheme for rigid bodies. Only has any effect 
5643:          ! if RIGIDINIT is also set. 
5644:          ELSE IF (WORD.EQ.'SLERP') THEN 
5645:             SLERPT=.TRUE.            
5646:              
5647: ! 5640: ! 
5648: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using5641: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using
5649: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations5642: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations
5650: ! or if the RMS force goes below SQVVGuessRMSTol.5643: ! or if the RMS force goes below SQVVGuessRMSTol.
5651: ! 5644: ! 
5652:          ELSE IF (WORD == 'SQVV') THEN5645:          ELSE IF (WORD == 'SQVV') THEN
5653:             SQVVGUESS=.TRUE.5646:             SQVVGUESS=.TRUE.
5654:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)5647:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)
5655:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)5648:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)
5656:          ELSE IF (WORD.EQ.'SSH') THEN5649:          ELSE IF (WORD.EQ.'SSH') THEN


r30598/newneb.f90 2016-06-15 11:30:09.658952891 +0100 r30597/newneb.f90 2016-06-15 11:30:10.666966231 +0100
 29:           USE MINIMISER3 29:           USE MINIMISER3
 30:           USE NEBOUTPUT 30:           USE NEBOUTPUT
 31:           USE NEBUTILS 31:           USE NEBUTILS
 32:           USE KEY, ONLY : UNRST, GROWSTRINGT, FREEZENODEST, DESMDEBUG, & 32:           USE KEY, ONLY : UNRST, GROWSTRINGT, FREEZENODEST, DESMDEBUG, &
 33:                & NEBMUPDATE, MUPDATE, BFGSSTEPS, NEBRESEEDT, & 33:                & NEBMUPDATE, MUPDATE, BFGSSTEPS, NEBRESEEDT, &
 34:                & INTCONMAX, ORDERI, ORDERJ, EPSALPHA, REDOBFGSSTEPS, &  34:                & INTCONMAX, ORDERI, ORDERJ, EPSALPHA, REDOBFGSSTEPS, & 
 35:                & NREPMAX, DISTREF, NEBKINITIAL, ADDREPT, REPPOW, REDOTSIM, MIN1REDO, MIN2REDO, PUSHOFF, & 35:                & NREPMAX, DISTREF, NEBKINITIAL, ADDREPT, REPPOW, REDOTSIM, MIN1REDO, MIN2REDO, PUSHOFF, &
 36:                & CONI, CONJ, AMHT, NUMGLY, REPI, REPJ, BULKT, D1INIT, D2INIT, & 36:                & CONI, CONJ, AMHT, NUMGLY, REPI, REPJ, BULKT, D1INIT, D2INIT, &
 37:                & REDOKADD, REDOPATH1, INTCONSTRAINTT, INTNEBIMAGES, & 37:                & REDOKADD, REDOPATH1, INTCONSTRAINTT, INTNEBIMAGES, &
 38:                & REDOPATH2, NREPI, NREPJ, REPCUT, NREPCUT, TWOD, RIGIDBODY, PERMDIST, WHOLEDNEB, & 38:                & REDOPATH2, NREPI, NREPJ, REPCUT, NREPCUT, TWOD, RIGIDBODY, PERMDIST, WHOLEDNEB, &
 39:                & CPPNEBT, VARIABLES, SLERPT 39:                & CPPNEBT, VARIABLES
 40:           USE GROWSTRINGUTILS, ONLY: GROWSTRING, TOTSTEPS 40:           USE GROWSTRINGUTILS, ONLY: GROWSTRING, TOTSTEPS
 41:           USE GSDATA, ONLY : KEYGSPRINT 41:           USE GSDATA, ONLY : KEYGSPRINT
 42:           USE MODGUESS,ONLY: GUESSPATHT,NINTERP 42:           USE MODGUESS,ONLY: GUESSPATHT,NINTERP
 43:           USE MODMEC,ONLY: MECCANOT           43:           USE MODMEC,ONLY: MECCANOT          
 44:           USE INTCOMMONS, ONLY : DESMINT, INTINTERPT, NINTIM, NDIH, DIHINFO, ALIGNDIR, PREVDIH, NINTC 44:           USE INTCOMMONS, ONLY : DESMINT, INTINTERPT, NINTIM, NDIH, DIHINFO, ALIGNDIR, PREVDIH, NINTC
 45:           USE INTCUTILS, ONLY : INTINTERPOLATE, CART2INT 45:           USE INTCUTILS, ONLY : INTINTERPOLATE, CART2INT
 46:           USE SPFUNCTS, ONLY : DUMPCOORDS 46:           USE SPFUNCTS, ONLY : DUMPCOORDS
 47:           USE NEBTOCONNECT 47:           USE NEBTOCONNECT
 48:           USE AMHGLOBALS, ONLY : NMRES 48:           USE AMHGLOBALS, ONLY : NMRES
 49:           USE COMMONS,ONLY: PARAM1,PARAM2,PARAM3,REDOPATHNEB,ZSYM,DEBUG 49:           USE COMMONS,ONLY: PARAM1,PARAM2,PARAM3,REDOPATHNEB,ZSYM,DEBUG
183:              DO J1 = 2,NIMAGE+2183:              DO J1 = 2,NIMAGE+2
184:                 ! align all other dihedrals to start184:                 ! align all other dihedrals to start
185:                 DIHINFO(J1,:) = DIHINFO(1,:)185:                 DIHINFO(J1,:) = DIHINFO(1,:)
186:              ENDDO186:              ENDDO
187: 187: 
188:              ALIGNDIR = .TRUE.188:              ALIGNDIR = .TRUE.
189:              PREVDIH => DIHINFO(NIMAGE+2,:)189:              PREVDIH => DIHINFO(NIMAGE+2,:)
190:              CALL CART2INT(FINFIN,XYZ(NOPT*(NIMAGE+1)+1:))190:              CALL CART2INT(FINFIN,XYZ(NOPT*(NIMAGE+1)+1:))
191:              ALIGNDIR = .FALSE.191:              ALIGNDIR = .FALSE.
192:           ELSE192:           ELSE
193:              IF(RIGIDINIT .AND. (RIGIDMOLECULEST .OR. SLERPT)) THEN  ! Should we replace this by IF(RIGIDINIT)?193:              IF(RIGIDMOLECULEST) THEN  ! Should we replace this by IF(RIGIDINIT)?
194:              ! XYZ contains NIMAGE+2 sets of NOPT coordinates, one for each image194:              ! XYZ contains NIMAGE+2 sets of NOPT coordinates, one for each image
195:              ! and each endpoint.195:              ! and each endpoint.
196:              ! Fill in the first DEGFREEDOMS coordinates of each image with the196:              ! Fill in the first DEGFREEDOMS coordinates of each image with the
197:              ! rigid-body coordinates, then the remainder with 0's (they get ignored197:              ! rigid-body coordinates, then the remainder with 0's (they get ignored
198:              ! by calls to the potential anyway)198:              ! by calls to the potential anyway)
199:                 CALL TRANSFORMCTORIGID(QQ,RIGIDQ)199:                 CALL TRANSFORMCTORIGID(QQ,RIGIDQ)
200:                 CALL TRANSFORMCTORIGID(FINFIN,RIGIDFIN)200:                 CALL TRANSFORMCTORIGID(FINFIN,RIGIDFIN)
201:                 XYZ(:DEGFREEDOMS) = RIGIDQ(:)201:                 XYZ(:DEGFREEDOMS) = RIGIDQ(:)
202:                 XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+DEGFREEDOMS) = RIGIDFIN(:)202:                 XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+DEGFREEDOMS) = RIGIDFIN(:)
203:                 ATOMRIGIDCOORDT = .FALSE.203:                 ATOMRIGIDCOORDT = .FALSE.
381:                                   IF (NDONE.GT.NIMAGE) EXIT imageloop1381:                                   IF (NDONE.GT.NIMAGE) EXIT imageloop1
382:                                ENDIF382:                                ENDIF
383:                             ENDDO383:                             ENDDO
384:                             LDTOTAL=LDTOTAL+XDUMMY384:                             LDTOTAL=LDTOTAL+XDUMMY
385:                          ENDDO imageloop1385:                          ENDDO imageloop1
386:                          XYZ(NOPT+1:NOPT*(NIMAGE+1))=INTNEBIMAGES(1:NOPT*NIMAGE)386:                          XYZ(NOPT+1:NOPT*(NIMAGE+1))=INTNEBIMAGES(1:NOPT*NIMAGE)
387:                       ENDIF387:                       ENDIF
388:                    ELSEIF (INTINTERPT) THEN388:                    ELSEIF (INTINTERPT) THEN
389:                       CALL INTINTERPOLATE(QQ,FINFIN,NINTIM,NIMAGE,X,DESMDEBUG,FAILED)389:                       CALL INTINTERPOLATE(QQ,FINFIN,NINTIM,NIMAGE,X,DESMDEBUG,FAILED)
390:                    ELSE390:                    ELSE
391:                       IF(RIGIDINIT .AND. (RIGIDMOLECULEST .OR. SLERPT)) THEN391:                       IF(RIGIDMOLECULEST) THEN
392:                       ! MAKEIMAGE doesn't actually use the initial/final coordinates passed into it unless392:                       ! MAKEIMAGE doesn't actually use the initial/final coordinates passed into it unless
393:                       ! we have MORPHT set (it just uses the endpoints in XYZ instead) so the following393:                       ! we have MORPHT set (it just uses the endpoints in XYZ instead) so the following
394:                       ! line is usually unnecessary.394:                       ! line is usually unnecessary.
395:                         CALL MAKEIMAGE(EINITIAL,EFINAL,RIGIDQ,RIGIDFIN)395:                         CALL MAKEIMAGE(EINITIAL,EFINAL,RIGIDQ,RIGIDFIN)
396:                       ELSE396:                       ELSE
397:                         CALL MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN)397:                         CALL MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN)
398:                       ENDIF398:                       ENDIF
399:                    ENDIF399:                    ENDIF
400:                 ENDIF400:                 ENDIF
401:              ENDIF401:              ENDIF
402:              IF (UNRST) DEALLOCATE(MYPTS) ! JMC402:              IF (UNRST) DEALLOCATE(MYPTS) ! JMC
403: 403: 
404: ! hk286404: ! hk286
405:              IF (RIGIDINIT .AND. .NOT. (SLERPT .OR. RIGIDMOLECULEST)) THEN  ! sn402: This whole block could possibly405:              IF (RIGIDINIT .AND. .NOT. RIGIDMOLECULEST) THEN  ! sn402: This whole block could possibly
406:              ! be removed if the earlier block is changed as suggested406:              ! be removed if the earlier block is changed as suggested
407:                 CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)407:                 CALL GENRIGID_IMAGE_CTORIGID(NIMAGE, XYZ)
408:                 ATOMRIGIDCOORDT = .FALSE.408:                 ATOMRIGIDCOORDT = .FALSE.
409:              ENDIF409:              ENDIF
410: ! hk286410: ! hk286
411: 411: 
412:              ! preoptimise if requested        412:              ! preoptimise if requested        
413:              IF (SQVVGUESS) THEN413:              IF (SQVVGUESS) THEN
414:                 STEPTOT = 5.0D0 ! TO AVOID GRADIENT SCALING DURING SQVV414:                 STEPTOT = 5.0D0 ! TO AVOID GRADIENT SCALING DURING SQVV
415:                 CALL NEBSQVV(NOPT*NIMAGE)415:                 CALL NEBSQVV(NOPT*NIMAGE)


r30598/nnutils.f90 2016-06-15 11:30:09.850955432 +0100 r30597/nnutils.f90 2016-06-15 11:30:10.862968824 +0100
120:           USE KEY,ONLY: MORPHT, MSTEPS, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, &120:           USE KEY,ONLY: MORPHT, MSTEPS, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, &
121:                         BULKT, RBAAT, NTSITES, AMBERT, LOCALPERMDIST, NRBGROUP, RBGROUP, RBNINGROUP, NRBTRIES, NABT,TWOD, &121:                         BULKT, RBAAT, NTSITES, AMBERT, LOCALPERMDIST, NRBGROUP, RBGROUP, RBNINGROUP, NRBTRIES, NABT,TWOD, &
122:                         RIGIDBODY, AVOID_COLLISIONS, COLL_TOL, BULK_BOXVEC122:                         RIGIDBODY, AVOID_COLLISIONS, COLL_TOL, BULK_BOXVEC
123:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &123:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &
124:                                  CHICDNEB124:                                  CHICDNEB
125:           USE COMMONS, ONLY : ZSYM, NRBSITES, DEBUG125:           USE COMMONS, ONLY : ZSYM, NRBSITES, DEBUG
126:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT126:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT
127:           USE MODAMBER9, ONLY: AMBERICT, AMBICDNEBT, NICTOT127:           USE MODAMBER9, ONLY: AMBERICT, AMBICDNEBT, NICTOT
128:           USE PORFUNCS 128:           USE PORFUNCS 
129:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE129:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE
130:           USE KEY, ONLY: FILTH,FILTHSTR, MULTISITEPYT, ALIGNRBST, N_TO_ALIGN, TO_ALIGN, SLERPT 130:           USE KEY, ONLY: FILTH,FILTHSTR, MULTISITEPYT, ALIGNRBST, N_TO_ALIGN, TO_ALIGN !sn402
131:           USE GENRIGID, ONLY: RIGIDINIT, RIGIDMOLECULEST, DEGFREEDOMS, NRIGIDBODY, GENRIGID_IMAGE_CTORIGID, GENRIGID_IMAGE_RIGIDTOC !sn402131:           USE GENRIGID, ONLY: RIGIDMOLECULEST, DEGFREEDOMS, NRIGIDBODY, GENRIGID_IMAGE_CTORIGID, GENRIGID_IMAGE_RIGIDTOC !sn402
132:           IMPLICIT NONE132:           IMPLICIT NONE
133:           INTEGER :: I, J1, ITDONE, INTERVAL, NDONE, J2, J, K, ISTAT, J3, J4, J5, NDUMMY, J6, J7, JWORST2, JWORST3133:           INTEGER :: I, J1, ITDONE, INTERVAL, NDONE, J2, J, K, ISTAT, J3, J4, J5, NDUMMY, J6, J7, JWORST2, JWORST3
134:           DOUBLE PRECISION,ALLOCATABLE :: DELTAX(:)134:           DOUBLE PRECISION,ALLOCATABLE :: DELTAX(:)
135: !          DOUBLE PRECISION, ALLOCATABLE :: QTN(:,:), PTN(:,:)135: !          DOUBLE PRECISION, ALLOCATABLE :: QTN(:,:), PTN(:,:)
136:           DOUBLE PRECISION DPRAND, SHIFT(NOPT), DUMMY, DLENGTH, DUMMY2, ENERGY, VNEW(NOPT), LRMS, &136:           DOUBLE PRECISION DPRAND, SHIFT(NOPT), DUMMY, DLENGTH, DUMMY2, ENERGY, VNEW(NOPT), LRMS, &
137:    &                       QS(NOPT), QF(NOPT), EINITIAL, EFINAL, COORDS(NOPT), SFRAC137:    &                       QS(NOPT), QF(NOPT), EINITIAL, EFINAL, COORDS(NOPT), SFRAC
138: !          DOUBLE PRECISION THETAH, ST, CT, P(3), FCT138: !          DOUBLE PRECISION THETAH, ST, CT, P(3), FCT
139:           DOUBLE PRECISION THETA, XSHIFT, YSHIFT, ZSHIFT, EWORST, EWORST_NEW139:           DOUBLE PRECISION THETA, XSHIFT, YSHIFT, ZSHIFT, EWORST, EWORST_NEW
140:           DOUBLE PRECISION, ALLOCATABLE :: XINITIAL(:), XIMAGE(:,:), SAVE_BAND(:)140:           DOUBLE PRECISION, ALLOCATABLE :: XINITIAL(:), XIMAGE(:,:), SAVE_BAND(:)
141:                     DOUBLE PRECISION,DIMENSION(:)         :: QQ,FINFIN141:                     DOUBLE PRECISION,DIMENSION(:)         :: QQ,FINFIN
246:              GCIMAGE=NIMAGE246:              GCIMAGE=NIMAGE
247:              QS(1:3*NATOMS)=XYZ(1:NOPT)247:              QS(1:3*NATOMS)=XYZ(1:NOPT)
248:              QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)248:              QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)
249:              CALL GCLBFGS(QS,QF,XIMAGE,3*NATOMS+1,GCUPDATE,XINITIAL,.FALSE.,GCCONV,MFLAG,ENERGY,RMS,GCSTEPS,.TRUE.,ITDONE,.TRUE.)249:              CALL GCLBFGS(QS,QF,XIMAGE,3*NATOMS+1,GCUPDATE,XINITIAL,.FALSE.,GCCONV,MFLAG,ENERGY,RMS,GCSTEPS,.TRUE.,ITDONE,.TRUE.)
250:              DO I=2,NIMAGE+1250:              DO I=2,NIMAGE+1
251:                   XYZ(NOPT*(I-1)+1:NOPT*I) = XIMAGE(I-1,1:NOPT)251:                   XYZ(NOPT*(I-1)+1:NOPT*I) = XIMAGE(I-1,1:NOPT)
252:              ENDDO252:              ENDDO
253:              DEALLOCATE(XINITIAL,XIMAGE)253:              DEALLOCATE(XINITIAL,XIMAGE)
254: 254: 
255: !     DC430 >255: !     DC430 >
256:           ELSEIF (RBAAT .OR. RIGIDMOLECULEST .OR. (RIGIDINIT .AND. SLERPT)) THEN  ! sn402: should figure out whether this bit is256:           ELSEIF (RBAAT .OR. RIGIDMOLECULEST) THEN  ! sn402: should figure out whether this bit is
257:           ! applicable to RIGIDINIT in general, or just to RIGIDMOLECULEST257:           ! applicable to RIGIDINIT in general, or just to RIGIDMOLECULEST
258:              258: 
259:             ALLOCATE(DELTAX(NOPT))259:              ALLOCATE(DELTAX(NOPT))
260: 260: 
261:             ! DELTAX is the change in each coordinate between each consecutive pair of images261:             ! DELTAX is the change in each coordinate between each consecutive pair of images
262:             ! Note: we are performing some unnecessary calculations here because we262:             ! Note: we are performing some unnecessary calculations here because we
263:             ! compute DELTAX for the angle-axis coordinates as well as the CoM coordinates, but we never263:             ! compute DELTAX for the angle-axis coordinates as well as the CoM coordinates, but we never
264:             ! use the former. I hope this won't be a big problem for performance so I've left it in to keep264:             ! use the former. I hope this won't be a big problem for performance so I've left it in to keep
265:             ! the code tidy.265:             ! the code tidy.
266:             IF (BULKT) THEN266:             IF (BULKT) THEN
267:                 DO K=1,NATOMS267:                 DO K=1,NATOMS
268:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &268:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &
269:   &                    -BULK_BOXVEC(1)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1))/BULK_BOXVEC(1))269:   &                    -BULK_BOXVEC(1)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1))/BULK_BOXVEC(1))
273:   &                    -BULK_BOXVEC(3)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3))/BULK_BOXVEC(3))273:   &                    -BULK_BOXVEC(3)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3))/BULK_BOXVEC(3))
274:                 ENDDO274:                 ENDDO
275:                 DELTAX(1:NOPT)=DELTAX(1:NOPT)/(NIMAGE+1)275:                 DELTAX(1:NOPT)=DELTAX(1:NOPT)/(NIMAGE+1)
276:             ELSE276:             ELSE
277:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)277:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)
278:             ENDIF278:             ENDIF
279:             279:             
280:             ! Linear interpolation of the centre of mass coordinates for all images280:             ! Linear interpolation of the centre of mass coordinates for all images
281:              DO I = 2, NIMAGE+1281:              DO I = 2, NIMAGE+1
282:                 J = NOPT*(I-1)282:                 J = NOPT*(I-1)
283:                 IF(RIGIDINIT) THEN283:                 IF(RIGIDMOLECULEST) THEN
284:                     XYZ(J+1:J+3*NRIGIDBODY) = XYZ(1:3*NRIGIDBODY) + DELTAX(1:3*NRIGIDBODY)*(I-1)284:                     XYZ(J+1:J+3*NRIGIDBODY) = XYZ(1:3*NRIGIDBODY) + DELTAX(1:3*NRIGIDBODY)*(I-1)
285:                 ELSE285:                 ELSE
286:                     XYZ(J+1:J+NOPT/2) = XYZ(1:NOPT/2) + DELTAX(1:NOPT/2)*(I-1)286:                     XYZ(J+1:J+NOPT/2) = XYZ(1:NOPT/2) + DELTAX(1:NOPT/2)*(I-1)
287: !                   PRINT *, I, XYZ(J+9), XYZ(J+9) - XYZ(J+9-NOPT)287: !                   PRINT *, I, XYZ(J+9), XYZ(J+9) - XYZ(J+9-NOPT)
288:                 ENDIF288:                 ENDIF
289:              ENDDO289:              ENDDO
290:             ! Interpolate the angle-axis vector components290:             ! Interpolate the angle-axis vector components
291:             IF(RIGIDINIT) THEN291:             IF(RIGIDMOLECULEST) THEN
292:                 DO J1 = 1, NRIGIDBODY292:                 DO J1 = 1, NRIGIDBODY
293:                     ! Interpolate AA vector for body J1 along the shortest path between the endpoints293:                     ! Interpolate AA vector for body J1 along the shortest path between the endpoints
294:                     CALL iSLERP(J1, .FALSE., 3*NRIGIDBODY)294:                     CALL iSLERP(J1, .FALSE., 3*NRIGIDBODY)
295:                 ENDDO295:                 ENDDO
296:             ELSE  ! i.e. RBAAT296:             ELSE  ! i.e. RBAAT
297:             ! Recall, NATOMS is 1/3*(length of the coords array), i.e. NATOMS = 2*(no. of rigid bodies).297:             ! Recall, NATOMS is 1/3*(length of the coords array), i.e. NATOMS = 2*(no. of rigid bodies).
298:                 DO J1 = 1, NATOMS/2298:                 DO J1 = 1, NATOMS/2
299:                     CALL iSLERP(J1, .FALSE., NOPT/2)299:                     CALL iSLERP(J1, .FALSE., NOPT/2)
300:                 ENDDO300:                 ENDDO
301:             ENDIF301:             ENDIF
302: 302: 
303: 303: 
304:             IF(RIGIDINIT) THEN304:             IF(RIGIDMOLECULEST) THEN
305:             ! Deal with any loose atoms (normal linear interpolation)305:             ! Deal with any loose atoms (normal linear interpolation)
306:                 DO I=2, NIMAGE+1306:                 DO I=2, NIMAGE+1
307:                     J = NOPT*(I-1)307:                     J = NOPT*(I-1)
308:                     XYZ(J+6*NRIGIDBODY+1:J+DEGFREEDOMS) = XYZ(6*NRIGIDBODY+1:DEGFREEDOMS) + &308:                     XYZ(J+6*NRIGIDBODY+1:J+DEGFREEDOMS) = XYZ(6*NRIGIDBODY+1:DEGFREEDOMS) + &
309:                                                             & DELTAX(6*NRIGIDBODY+1:DEGFREEDOMS)*(I-1)309:                                                             & DELTAX(6*NRIGIDBODY+1:DEGFREEDOMS)*(I-1)
310:                 ENDDO310:                 ENDDO
311:             ENDIF311:             ENDIF
312: 312: 
313:             IF (RIGIDINIT .AND. AVOID_COLLISIONS) THEN313:             IF (RIGIDMOLECULEST .AND. AVOID_COLLISIONS) THEN
314:             ! Test to see if the interpolation has given us a bad path in which two rigid bodies collide,314:             ! Test to see if the interpolation has given us a bad path in which two rigid bodies collide,
315:             ! and if so we reinterpolate this rigid body with the rotation in the opposite sense.315:             ! and if so we reinterpolate this rigid body with the rotation in the opposite sense.
316:             ! This section is currently coded only for RIGIDMOLECULEST (i.e. GENRIGID)316:             ! This section is currently coded only for RIGIDMOLECULEST (i.e. GENRIGID)
317: 317: 
318:                 ! Compute EWORST, the energy of the highest-energy image318:                 ! Compute EWORST, the energy of the highest-energy image
319:                 CALL COMPUTE_HIGHEST_IMAGE(XYZ, EWORST, NDUMMY)319:                 CALL COMPUTE_HIGHEST_IMAGE(XYZ, EWORST, NDUMMY)
320: 320: 
321:                 IF (EWORST .GT. COLL_TOL) THEN321:                 IF (EWORST .GT. COLL_TOL) THEN
322: 322: 
323:                     IF(DEBUG) WRITE(*,*) "makeimage> Bad interpolation; suspect two rigid bodies have collided"323:                     IF(DEBUG) WRITE(*,*) "makeimage> Bad interpolation; suspect two rigid bodies have collided"
940: 940: 
941:         LOGICAL, INTENT(IN) :: LONG_WAY941:         LOGICAL, INTENT(IN) :: LONG_WAY
942:         DOUBLE PRECISION,ALLOCATABLE :: QTN(:,:), PTN(:,:)942:         DOUBLE PRECISION,ALLOCATABLE :: QTN(:,:), PTN(:,:)
943:         DOUBLE PRECISION THETA, THETAH, ST, CT, P(3), FCT943:         DOUBLE PRECISION THETA, THETAH, ST, CT, P(3), FCT
944:         ! J1 is the index of the rigid body we are interested in. RBOFFSET is the point within the list944:         ! J1 is the index of the rigid body we are interested in. RBOFFSET is the point within the list
945:         ! of coordinates for this image at which the angle-axis components begin to be listed.945:         ! of coordinates for this image at which the angle-axis components begin to be listed.
946:         ! RBOFFSET is the index of the COORDS array at which the angle-axis coordinates start. The946:         ! RBOFFSET is the index of the COORDS array at which the angle-axis coordinates start. The
947:         ! expression for this value is different depending on which rigid-body representation is being used947:         ! expression for this value is different depending on which rigid-body representation is being used
948:         ! (GENRIGID or RBAA) so it must be passed as a parameter.948:         ! (GENRIGID or RBAA) so it must be passed as a parameter.
949:         INTEGER J1, K, J, I, RBOFFSET949:         INTEGER J1, K, J, I, RBOFFSET
950:         950: 
951:         ! Lists of parallel and perpendicular quaternions for body J1 in all the images.951:         ! Lists of parallel and perpendicular quaternions for body J1 in all the images.
952:         ALLOCATE(QTN(NIMAGE+2,4))952:         ALLOCATE(QTN(NIMAGE+2,4))
953:         ALLOCATE(PTN(NIMAGE+2,4))953:         ALLOCATE(PTN(NIMAGE+2,4))
954: 954: 
955: !       CONVERT FROM AA -> QUATERNION FOR THE INITIAL AND FINAL FRAMES955: !       CONVERT FROM AA -> QUATERNION FOR THE INITIAL AND FINAL FRAMES
956:         DO K = 1, (NIMAGE+2), (NIMAGE+1)956:         DO K = 1, (NIMAGE+2), (NIMAGE+1)
957:            ! P is the angle-axis vector for rigid body J1 in the initial or final frame957:            ! P is the angle-axis vector for rigid body J1 in the initial or final frame
958:            J = (K-1)*NOPT + RBOFFSET + 3*(J1-1)958:            J = (K-1)*NOPT + RBOFFSET + 3*(J1-1)
959:            P(:)   = XYZ(J+1:J+3)959:            P(:)   = XYZ(J+1:J+3)
960: 960: 


r30598/potential.f 2016-06-15 11:30:10.438963213 +0100 r30597/potential.f 2016-06-15 11:30:11.482977029 +0100
 94:                 IF (DEBUG) THEN 94:                 IF (DEBUG) THEN
 95:                     IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .GT. 1.0E-12)) THEN 95:                     IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .GT. 1.0E-12)) THEN
 96:                         WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = FALSE but coords appear to be atomistic" 96:                         WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = FALSE but coords appear to be atomistic"
 97: !                        WRITE(*,*) "Largest absolute value:", MAXVAL(ABS(COORDS(DEGFREEDOMS+1:))) 97: !                        WRITE(*,*) "Largest absolute value:", MAXVAL(ABS(COORDS(DEGFREEDOMS+1:)))
 98:                         WRITE(*,*) COORDS(DEGFREEDOMS+1:) 98:                         WRITE(*,*) COORDS(DEGFREEDOMS+1:)
 99: !                         WRITE(*,*) COORDS 99: !                         WRITE(*,*) COORDS
100:                     ENDIF100:                     ENDIF
101:                 ENDIF101:                 ENDIF
102:                 XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)102:                 XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)
103:                 CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)103:                 CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)
104:             ELSE IF (DEBUG) THEN104:             ELSE IF (ATOMRIGIDCOORDT .AND. DEBUG) THEN
105:                 IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .LT. 1.0D-12)) THEN105:                 IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .LT. 1.0D-12)) THEN
106:                     WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = TRUE but coords appear to be angleaxis"106:                     WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = TRUE but coords appear to be angleaxis"
107: !                    WRITE(*,*) "Smallest absolute value:", MINVAL(ABS(COORDS(DEGFREEDOMS+1:)))107: !                    WRITE(*,*) "Smallest absolute value:", MINVAL(ABS(COORDS(DEGFREEDOMS+1:)))
108:                     WRITE(*,*) COORDS(DEGFREEDOMS+1:)108:                     WRITE(*,*) COORDS(DEGFREEDOMS+1:)
109:                 ENDIF109:                 ENDIF
110:             ENDIF110:             ENDIF
111:         ENDIF111:         ENDIF
112:          ! 112:          ! 
113:          ! SSTEST needs to be true if we really want an analytic Hessian this step.113:          ! SSTEST needs to be true if we really want an analytic Hessian this step.
114:          ! 114:          ! 
3879:                DO J1=1,NATOMS3879:                DO J1=1,NATOMS
3880:                   IF (FROZEN(J1)) THEN3880:                   IF (FROZEN(J1)) THEN
3881:                      VNEW(3*(J1-1)+1)=0.0D03881:                      VNEW(3*(J1-1)+1)=0.0D0
3882:                      VNEW(3*(J1-1)+2)=0.0D03882:                      VNEW(3*(J1-1)+2)=0.0D0
3883:                      VNEW(3*(J1-1)+3)=0.0D03883:                      VNEW(3*(J1-1)+3)=0.0D0
3884:                   ENDIF3884:                   ENDIF
3885:                ENDDO3885:                ENDDO
3886:             ENDIF3886:             ENDIF
3887: 3887: 
3888:             ! sn402: other half of hack to allow genrigid to work with arbitrary potential3888:             ! sn402: other half of hack to allow genrigid to work with arbitrary potential
3889:             IF (RIGIDMOLECULEST .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN3889:             IF (RIGIDMOLECULEST .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN
3890:                 IF (STEST) THEN3890:                 IF (STEST) THEN
3891:                      CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)3891:                      CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)
3892:                          HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D03892:                          HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0
3893:                          HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D03893:                          HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0
3894:                          HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)3894:                          HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)
3895:                 ENDIF3895:                 ENDIF
3896:                 CALL TRANSFORMGRAD(VNEW, XRIGIDCOORDS, XRIGIDGRAD)3896:                 CALL TRANSFORMGRAD(VNEW, XRIGIDCOORDS, XRIGIDGRAD)
3897:                 VNEW(1:DEGFREEDOMS) = XRIGIDGRAD(1:DEGFREEDOMS)3897:                 VNEW(1:DEGFREEDOMS) = XRIGIDGRAD(1:DEGFREEDOMS)
3898:                 VNEW(DEGFREEDOMS+1:3*NATOMS) = 0.0D03898:                 VNEW(DEGFREEDOMS+1:3*NATOMS) = 0.0D0
3899:             ENDIF3899:             ENDIF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0