hdiff output

r32610/fetchz.f 2017-05-25 16:30:09.482211063 +0100 r32609/fetchz.f 2017-05-25 16:30:10.382222880 +0100
266:             WRITE(*,'(A,G20.10,A)') '               1/kT (mol/kcal)=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))266:             WRITE(*,'(A,G20.10,A)') '               1/kT (mol/kcal)=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))
267:             WRITE(*,'(A,2I6)')     '               Number of O and H atoms=',NOXYGEN,(RPDOF/3)-NOXYGEN267:             WRITE(*,'(A,2I6)')     '               Number of O and H atoms=',NOXYGEN,(RPDOF/3)-NOXYGEN
268:             CALL SDINIT(NOXYGEN,NCHARGE)268:             CALL SDINIT(NOXYGEN,NCHARGE)
269:             NATOMS=NOPT/3269:             NATOMS=NOPT/3
270:          ELSEIF (RPSYSTEM(1:2).EQ.'TT') THEN270:          ELSEIF (RPSYSTEM(1:2).EQ.'TT') THEN
271:             WRITE(*,'(A,G20.10,A)')  '              1/kT a.u.=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))271:             WRITE(*,'(A,G20.10,A)')  '              1/kT a.u.=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))
272:             RPBETA=RPBETA*0.00159360144367 ! convert to kcal/mol units272:             RPBETA=RPBETA*0.00159360144367 ! convert to kcal/mol units
273:             WRITE(*,'(A,G20.10,A)')  '              1/kT (mol/kcal)=',RPBETA273:             WRITE(*,'(A,G20.10,A)')  '              1/kT (mol/kcal)=',RPBETA
274:             WRITE(*,'(A,I0,A)') '              TTM3-F potential for a cluster of ', RPDOF/9, ' water molecules'274:             WRITE(*,'(A,I0,A)') '              TTM3-F potential for a cluster of ', RPDOF/9, ' water molecules'
275:             NATOMS=NOPT/3275:             NATOMS=NOPT/3
276:          ELSEIF (RPSYSTEM(1:2).EQ.'MB') THEN 
277:             WRITE(*,'(A,G20.10,A)')  '              1/kT a.u.=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM)) 
278:             RPBETA=RPBETA*0.00159360144367 ! convert to kcal/mol units 
279:             WRITE(*,'(A,G20.10,A)')  '              1/kT (mol/kcal)=',RPBETA 
280:             WRITE(*,'(A,I0,A)') '              MBPOL potential for a cluster of ', RPDOF/9, ' water molecules' 
281:             NATOMS=NOPT/3 
282:          ELSEIF (RPSYSTEM(1:3).EQ.'MAL') THEN276:          ELSEIF (RPSYSTEM(1:3).EQ.'MAL') THEN
283:             WRITE(*,'(A,G20.10,A)')  '              1/kT a.u.=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))277:             WRITE(*,'(A,G20.10,A)')  '              1/kT a.u.=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))
284:             WRITE(*,'(A,I0,A)') '              malonaldehyde potential '278:             WRITE(*,'(A,I0,A)') '              malonaldehyde potential '
285:             NATOMS=NOPT/3279:             NATOMS=NOPT/3
286:          ELSEIF (RPSYSTEM(1:3).EQ.'MCY') THEN280:          ELSEIF (RPSYSTEM(1:3).EQ.'MCY') THEN
287:             WRITE(*,'(A,G20.10,A)')  '              1/kT (a.u.)=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))281:             WRITE(*,'(A,G20.10,A)')  '              1/kT (a.u.)=',RPBETA,' potential type=' // TRIM(ADJUSTL(RPSYSTEM))
288:             WRITE(*,'(A,I1,A,I0,A)') '              VRT(MCY-5f) flexible water dimer potential'282:             WRITE(*,'(A,I1,A,I0,A)') '              VRT(MCY-5f) flexible water dimer potential'
289:             NOXYGEN=RPDOF/9283:             NOXYGEN=RPDOF/9
290:             IF (NOXYGEN /= 2) THEN284:             IF (NOXYGEN /= 2) THEN
291:                PRINT *, 'ERROR: VRT(MCY-5f) only for water dimer; NOXYGEN = ', NOXYGEN285:                PRINT *, 'ERROR: VRT(MCY-5f) only for water dimer; NOXYGEN = ', NOXYGEN


r32610/keywords.f 2017-05-25 16:30:09.710214057 +0100 r32609/keywords.f 2017-05-25 16:30:10.606225822 +0100
 83:          ! DOUBLE PRECISION :: ZBB,R0AA,R0AB,R0BB 83:          ! DOUBLE PRECISION :: ZBB,R0AA,R0AB,R0BB
 84:          DOUBLE PRECISION :: XX, EPSAB, EPSBB, SIGAB, SIGBB, RANDOM, DPRAND 84:          DOUBLE PRECISION :: XX, EPSAB, EPSBB, SIGAB, SIGBB, RANDOM, DPRAND
 85:          LOGICAL END, SKIPBL, CLEAR, ECHO, CAT, CISTRANS, RBSYMTEST, YESNO 85:          LOGICAL END, SKIPBL, CLEAR, ECHO, CAT, CISTRANS, RBSYMTEST, YESNO
 86:          CHARACTER WORD*25, WW*20, PBC*3 86:          CHARACTER WORD*25, WW*20, PBC*3
 87:          CHARACTER(LEN=80) :: FILENAME,DUMMYS 87:          CHARACTER(LEN=80) :: FILENAME,DUMMYS
 88:          CHARACTER WORD2*25 88:          CHARACTER WORD2*25
 89:          ! COMMON /BIN/ NTYPEA,AAA,AAB,ABB,PAA,PAB,PBB,QAA,QAB,QBB,ZAA,ZAB 89:          ! COMMON /BIN/ NTYPEA,AAA,AAB,ABB,PAA,PAB,PBB,QAA,QAB,QBB,ZAA,ZAB
 90:          ! COMMON /BIN/ ZBB,R0AA,R0AB,R0BB 90:          ! COMMON /BIN/ ZBB,R0AA,R0AB,R0BB
 91:          COMMON /BIN/ EPSAB, EPSBB, SIGAB, SIGBB,NTYPEA 91:          COMMON /BIN/ EPSAB, EPSBB, SIGAB, SIGBB,NTYPEA
 92:  92: 
 93:          INTEGER NATOM, DMODE, NATOMSSAVE 93:          INTEGER NATOM, DMODE
 94:          DOUBLE PRECISION CHX(MXATMS), CHY(MXATMS), CHZ(MXATMS), CHMASS(MXATMS) 94:          DOUBLE PRECISION CHX(MXATMS), CHY(MXATMS), CHZ(MXATMS), CHMASS(MXATMS)
 95:          DOUBLE PRECISION DPERT 95:          DOUBLE PRECISION DPERT
 96:          DOUBLE PRECISION CHPMIN, CHPMAX, CHNMIN, CHNMAX 96:          DOUBLE PRECISION CHPMIN, CHPMAX, CHNMIN, CHNMAX
 97:          DOUBLE PRECISION, ALLOCATABLE :: LCONGEOM(:,:) 97:          DOUBLE PRECISION, ALLOCATABLE :: LCONGEOM(:,:)
 98:          DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:) 98:          DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:)
 99:          INTEGER ISEED 99:          INTEGER ISEED
100: 100: 
101:          DOUBLE PRECISION UNRX(NATOMS), UNRY(NATOMS), UNRZ(NATOMS) ! UNRES101:          DOUBLE PRECISION UNRX(NATOMS), UNRY(NATOMS), UNRZ(NATOMS) ! UNRES
102:          DOUBLE PRECISION DUMMY1(NATOMS)102:          DOUBLE PRECISION DUMMY1(NATOMS)
103: 103: 
727:          ! 727:          ! 
728:          GLJT=.FALSE.728:          GLJT=.FALSE.
729:          NGLJ=1 ! number of atom types729:          NGLJ=1 ! number of atom types
730:          ! 730:          ! 
731:          ! ds656> substrate field(s)731:          ! ds656> substrate field(s)
732:          MIEFT=.FALSE.732:          MIEFT=.FALSE.
733:          MIEF_PBCT=.FALSE.733:          MIEF_PBCT=.FALSE.
734:          MIEF_CUTT=.FALSE.734:          MIEF_CUTT=.FALSE.
735:          MIEF_BOX(1:3) = 1.0D9735:          MIEF_BOX(1:3) = 1.0D9
736:          MIEF_RCUT= 1.0D9736:          MIEF_RCUT= 1.0D9
737:          MAXIMFACTOR=10.0D0 
738:          !737:          !
739:          ! UNDOCUMENTED keywords/parameters738:          ! UNDOCUMENTED keywords/parameters
740:          ! 739:          ! 
741:          TWISTT=.FALSE.740:          TWISTT=.FALSE.
742:          LJADDT=.FALSE.741:          LJADDT=.FALSE.
743:          LJADD2T=.FALSE.742:          LJADD2T=.FALSE.
744:          LJADD3T=.FALSE.743:          LJADD3T=.FALSE.
745:          LJADD4T=.FALSE.744:          LJADD4T=.FALSE.
746:          NADDTARGET=1745:          NADDTARGET=1
747:          PYADDT=.FALSE.746:          PYADDT=.FALSE.
3852:             ELSE IF (WW .EQ. 'OFF') THEN3851:             ELSE IF (WW .EQ. 'OFF') THEN
3853:                MASST=.FALSE.3852:                MASST=.FALSE.
3854:             ENDIF3853:             ENDIF
3855: !3854: !
3856: !3855: !
3857: !3856: !
3858:             ELSE IF (WORD.EQ.'MALONALDEHYDE') THEN3857:             ELSE IF (WORD.EQ.'MALONALDEHYDE') THEN
3859:                MALONALDEHYDE=.TRUE.3858:                MALONALDEHYDE=.TRUE.
3860: 3859: 
3861:             ELSE IF (WORD.EQ.'MBPOL') THEN3860:             ELSE IF (WORD.EQ.'MBPOL') THEN
3862:                FRQCONV=5.123934D14 
3863:                CALL MBPOLINIT3861:                CALL MBPOLINIT
3864:                MBPOLT=.TRUE.3862:                MBPOLT=.TRUE.
3865: 3863: 
3866: ! 3864: ! 
3867: ! Maximum value for the smaller barrier height that is allowed to constitute a connection during the3865: ! Maximum value for the smaller barrier height that is allowed to constitute a connection during the
3868: ! Dijkstra connection procedure.3866: ! Dijkstra connection procedure.
3869: ! MAXMAXBARRIER specifies a maximum for the maximum barrier.3867: ! MAXMAXBARRIER specifies a maximum for the maximum barrier.
3870: ! MAXBARRIER requires both sides to be greater than MAXBARRIER to discard.3868: ! MAXBARRIER requires both sides to be greater than MAXBARRIER to discard.
3871: ! 3869: ! 
3872:          ELSE IF (WORD.EQ.'MAXBARRIER') THEN3870:          ELSE IF (WORD.EQ.'MAXBARRIER') THEN
3890: ! The deafult is 4.3888: ! The deafult is 4.
3891: ! 3889: ! 
3892:          ELSE IF (WORD.EQ.'MAXCON') THEN3890:          ELSE IF (WORD.EQ.'MAXCON') THEN
3893:             CALL READI(MAXCONUSE)3891:             CALL READI(MAXCONUSE)
3894: ! 3892: ! 
3895: ! The maximum energy increase above which mylbfgs will reject a proposed step.3893: ! The maximum energy increase above which mylbfgs will reject a proposed step.
3896: ! 3894: ! 
3897:          ELSE IF (WORD.EQ.'MAXERISE') THEN3895:          ELSE IF (WORD.EQ.'MAXERISE') THEN
3898:             CALL READF(MAXERISE)3896:             CALL READF(MAXERISE)
3899:             IF (NITEMS.GT.1) CALL READF(XMAXERISE)3897:             IF (NITEMS.GT.1) CALL READF(XMAXERISE)
3900: ! 
3901: ! This is the factor used in NEB/output.f90 to identify local maxima for hybrid EF ts searches. 
3902: ! It multiplies EDIFFTOL. 
3903: ! 
3904:          ELSE IF (WORD.EQ.'MAXIMFACTOR') THEN 
3905:             CALL READF(MAXIMFACTOR) 
3906: 3898: 
3907: ! Maximum number of failures allowed in a minimisation before giving up.3899: ! Maximum number of failures allowed in a minimisation before giving up.
3908: ! 3900: ! 
3909:          ELSE IF (WORD.EQ.'MAXFAIL') THEN3901:          ELSE IF (WORD.EQ.'MAXFAIL') THEN
3910:             IF (NITEMS.GT.1) CALL READI(NFAILMAX)3902:             IF (NITEMS.GT.1) CALL READI(NFAILMAX)
3911: !3903: !
3912: ! Specify single connection attempt for largest gap in current Dijkstra selected chain.3904: ! Specify single connection attempt for largest gap in current Dijkstra selected chain.
3913: ! 3905: ! 
3914:          ELSE IF (WORD.EQ.'MAXGAP') THEN3906:          ELSE IF (WORD.EQ.'MAXGAP') THEN
3915:             MAXGAPT=.TRUE.3907:             MAXGAPT=.TRUE.
6150: ! 6142: ! 
6151: ! Sanity checks.6143: ! Sanity checks.
6152: ! 6144: ! 
6153:             TEMPSTRING=TRIM(ADJUSTL(RPSYSTEM))6145:             TEMPSTRING=TRIM(ADJUSTL(RPSYSTEM))
6154:             IF (TEMPSTRING(1:2).EQ.' ') THEN6146:             IF (TEMPSTRING(1:2).EQ.' ') THEN
6155:                PRINT '(A)','keyword> ERROR *** Ring polymer potential type is not set'6147:                PRINT '(A)','keyword> ERROR *** Ring polymer potential type is not set'
6156:             ENDIF6148:             ENDIF
6157:             IF (RPIMAGES.LT.1) THEN6149:             IF (RPIMAGES.LT.1) THEN
6158:                PRINT '(A)','keyword> ERROR *** Ring polymer images too small, value is ',RPIMAGES6150:                PRINT '(A)','keyword> ERROR *** Ring polymer images too small, value is ',RPIMAGES
6159:             ENDIF6151:             ENDIF
6160:             IF (MBPOLT) THEN 
6161:                NATOMSSAVE=NATOMS 
6162:                NATOMS=NATOMS/(3*RPIMAGES) 
6163:                CALL MBPOLINIT 
6164:                DO J1=1,RPIMAGES    ! number of copies of water dimer, trimer, etc. 
6165:                   DO J2=1,NATOMS/3 ! number of water molecules in the cluster 
6166:                      ZSYM(NATOMS*(J1-1)+3*(J2-1)+1)='O ' 
6167:                      ZSYM(NATOMS*(J1-1)+3*(J2-1)+2)='H ' 
6168:                      ZSYM(NATOMS*(J1-1)+3*(J2-1)+3)='H ' 
6169:                   ENDDO 
6170:                ENDDO 
6171:                NATOMS=NATOMSSAVE 
6172:             ENDIF 
6173:  
6174:             RETURN6152:             RETURN
6175: 6153: 
6176: 6154: 
6177: ! 6155: ! 
6178: ! RKMIN calculates a steepest-descent path using gradient only information6156: ! RKMIN calculates a steepest-descent path using gradient only information
6179: ! with convergence criterion GMAX for the RMS force and initial precision6157: ! with convergence criterion GMAX for the RMS force and initial precision
6180: ! EPS. A fifth order Runga-Kutta algorithm is used.6158: ! EPS. A fifth order Runga-Kutta algorithm is used.
6181: ! 6159: ! 
6182:          ELSE IF (WORD.EQ.'RKMIN') THEN6160:          ELSE IF (WORD.EQ.'RKMIN') THEN
6183:             RKMIN=.TRUE.6161:             RKMIN=.TRUE.


r32610/mbpol.f90 2017-05-25 16:30:09.926216892 +0100 r32609/mbpol.f90 2017-05-25 16:30:10.826228711 +0100
 15: implicit none 15: implicit none
 16: nnode=natoms/3 16: nnode=natoms/3
 17: end subroutine mbpolinit 17: end subroutine mbpolinit
 18:  18: 
 19: subroutine mbpol(coordinates, energy, gradient, gradt) 19: subroutine mbpol(coordinates, energy, gradient, gradt)
 20: implicit none 20: implicit none
 21:          logical, intent(in)  :: gradt 21:          logical, intent(in)  :: gradt
 22: double precision, intent(in)  :: coordinates(nnode*9) 22: double precision, intent(in)  :: coordinates(nnode*9)
 23: double precision, intent(out) :: gradient(nnode*9), energy 23: double precision, intent(out) :: gradient(nnode*9), energy
 24:  24: 
 25: ! print '(A,I10,L5)','nnode, gradt=',nnode, gradt 
 26: ! print '(A)','coordinates:' 
 27: ! print '(3F20.10)',coordinates(1:nnode*9) 
 28: if (gradt) then 25: if (gradt) then
 29:   call mbpolenergygradient(nnode,energy,coordinates,gradient) 26:   call mbpolenergygradient(nnode,energy,coordinates,gradient)
 30: else 27: else
 31:   call mbpolenergy(nnode,energy,coordinates) 28:   call mbpolenergy(nnode,energy,coordinates)
 32: endif 29: endif
 33:  30: 
 34: end subroutine mbpol 31: end subroutine mbpol
 35:  32: 
 36: subroutine mbpolstep(coordinates,displacement,angle,translate,rotate) 33: subroutine mbpolstep(coordinates,displacement,angle,translate,rotate)
 37: implicit none 34: implicit none


r32610/potential.f 2017-05-25 16:30:10.158219939 +0100 r32609/potential.f 2017-05-25 16:30:11.050231652 +0100
 43:          IMPLICIT NONE 43:          IMPLICIT NONE
 44:  44: 
 45:          DOUBLE PRECISION COORDS(NOPT) 45:          DOUBLE PRECISION COORDS(NOPT)
 46:          DOUBLE PRECISION ENERGY 46:          DOUBLE PRECISION ENERGY
 47:          DOUBLE PRECISION VNEW(NOPT) 47:          DOUBLE PRECISION VNEW(NOPT)
 48:          LOGICAL GTEST, STEST 48:          LOGICAL GTEST, STEST
 49:          DOUBLE PRECISION RMS 49:          DOUBLE PRECISION RMS
 50:          LOGICAL PTEST, BOXTEST, file_exists 50:          LOGICAL PTEST, BOXTEST, file_exists
 51:  51: 
 52:          INTEGER J1, J2, J3, NN, MM, IPOT, NELEMENTS, NTYPE(105), NCOUNT, ISTART, NDUM, NPCALL, ECALL, FCALL, SCALL, 52:          INTEGER J1, J2, J3, NN, MM, IPOT, NELEMENTS, NTYPE(105), NCOUNT, ISTART, NDUM, NPCALL, ECALL, FCALL, SCALL,
 53:      1   J4, J5, J6, JSTART, ISTAT, NDUMMY,IDUM1,IDUM2, NATOMSSAVE 53:      1   J4, J5, J6, JSTART, ISTAT, NDUMMY,IDUM1,IDUM2
 54:          DOUBLE PRECISION P2, P3, BA, XLAMBDA, AMAT(3,3), AINV(3,3), TEMPX, BOXLX, 54:          DOUBLE PRECISION P2, P3, BA, XLAMBDA, AMAT(3,3), AINV(3,3), TEMPX, BOXLX,
 55:      1   TEMPXX, TEMPH(3,3), TEMPV(3), ENERGY1, ENERGY2,DOUBLEDUM1, 55:      1   TEMPXX, TEMPH(3,3), TEMPV(3), ENERGY1, ENERGY2,DOUBLEDUM1,
 56:      1   EDOUBLE, XTEMP, YTEMP, ZTEMP, GAMESR(3,3), GAMEST(3), 56:      1   EDOUBLE, XTEMP, YTEMP, ZTEMP, GAMESR(3,3), GAMEST(3),
 57:      2   TEMPA, TEMP(6), COORDSO(3*NATOMS), GRADO(3*NATOMS), GEMAX, 57:      2   TEMPA, TEMP(6), COORDSO(3*NATOMS), GRADO(3*NATOMS), GEMAX,
 58:      4   ETIME, FTIME, STIME, TIME, TIME0, DUMMY1, DUMMY2, DUMMY3, EPLUS, EMINUS, DUMMY, DOTOPT, DIST, DIST1, DIST2, RMAT(3,3) 58:      4   ETIME, FTIME, STIME, TIME, TIME0, DUMMY1, DUMMY2, DUMMY3, EPLUS, EMINUS, DUMMY, DOTOPT, DIST, DIST1, DIST2, RMAT(3,3)
 59:          INTEGER NSTART, NFINISH, NSTEP 59:          INTEGER NSTART, NFINISH, NSTEP
 60:          INTEGER LUNIT, GETUNIT 60:          INTEGER LUNIT, GETUNIT
 61:          DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:,:) ::  HGAUSS 61:          DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:,:) ::  HGAUSS
 62:          CHARACTER(LEN=87) ESTRING 62:          CHARACTER(LEN=87) ESTRING
 63:          CHARACTER(LEN=80) GPSTRING, NSTRING, FSTRING, FNAME, FNAME2, GSTRING,FROMFILE 63:          CHARACTER(LEN=80) GPSTRING, NSTRING, FSTRING, FNAME, FNAME2, GSTRING,FROMFILE
 64:          CHARACTER(LEN=132) STRING 64:          CHARACTER(LEN=132) STRING
 65:          COMMON /STRINGS/ ESTRING, GPSTRING, NSTRING, FSTRING 65:          COMMON /STRINGS/ ESTRING, GPSTRING, NSTRING, FSTRING
 66:          DOUBLE PRECISION C1, C2, C3, IZ, ROT, EDISP, EIND 66:          DOUBLE PRECISION C1, C2, C3, IZ, ROT, EDISP, EIND
 67:          LOGICAL ETEST, SSTEST, YESNO, AMIDEFAIL,GAUSSIANTEST 67:          LOGICAL ETEST, SSTEST, YESNO, AMIDEFAIL,GAUSSIANTEST
 68:          COMMON /CAS/ AMAT, AINV, NELEMENTS, NTYPE 68:          COMMON /CAS/ AMAT, AINV, NELEMENTS, NTYPE
 69:          COMMON /PCALL/ NPCALL, ECALL, FCALL, SCALL, ETIME, FTIME, STIME 69:          COMMON /PCALL/ NPCALL, ECALL, FCALL, SCALL, ETIME, FTIME, STIME
 70:          LOGICAL KNOWE, KNOWG, KNOWH 70:          LOGICAL KNOWE, KNOWG, KNOWH
 71:          COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 71:          COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 72:  72: 
 73:          DOUBLE PRECISION VPLUS(3*NATOMS), VMINUS(3*NATOMS), DIFF, EPLUS1, EPLUS2, EMINUS1, EMINUS2, EDUMMY 73:          DOUBLE PRECISION VPLUS(3*NATOMS), VMINUS(3*NATOMS), DIFF, EPLUS1, EPLUS2, EMINUS1, EMINUS2
 74:  74: 
 75:          ! double precision upperE, lowerE, deltaCoord, numericalGrad(3*NATOMS), RMSdiff 75:          ! double precision upperE, lowerE, deltaCoord, numericalGrad(3*NATOMS), RMSdiff
 76:          ! double precision dummyGrad(3*NATOMS), upperGrad(3*NATOMS), lowerGrad(3*NATOMS) 76:          ! double precision dummyGrad(3*NATOMS), upperGrad(3*NATOMS), lowerGrad(3*NATOMS)
 77:          ! double precision numericalSD, tempHess(3*NATOMS,3*NATOMS) 77:          ! double precision numericalSD, tempHess(3*NATOMS,3*NATOMS)
 78:          double precision HESSDUM(3*NATOMS,3*NATOMS) 78:          double precision HESSDUM(3*NATOMS,3*NATOMS)
 79:  79: 
 80:          ! sf344> NAB & AMBER additions 80:          ! sf344> NAB & AMBER additions
 81:          DOUBLE PRECISION,dimension(:),allocatable  ::  temphess 81:          DOUBLE PRECISION,dimension(:),allocatable  ::  temphess
 82:          DOUBLE PRECISION  :: GRAD1(3*NATOMS) 82:          DOUBLE PRECISION  :: GRAD1(3*NATOMS)
 83:          integer i,j,k 83:          integer i,j,k
 84:          ! hk286 - local rigid body 84:          ! hk286 - local rigid body
 85:          DOUBLE PRECISION :: XCOORDS(3*NATOMS), XRIGIDCOORDS(DEGFREEDOMS), XRIGIDGRAD(DEGFREEDOMS) 85:          DOUBLE PRECISION :: XCOORDS(3*NATOMS), XRIGIDCOORDS(DEGFREEDOMS), XRIGIDGRAD(DEGFREEDOMS)
 86:          DOUBLE PRECISION :: XRIGIDHESS(DEGFREEDOMS, DEGFREEDOMS) 86:          DOUBLE PRECISION :: XRIGIDHESS(DEGFREEDOMS, DEGFREEDOMS)
 87:          DOUBLE PRECISION,ALLOCATABLE :: FOCK(:,:), MOCOEFF(:,:) 87:          DOUBLE PRECISION,ALLOCATABLE :: FOCK(:,:), MOCOEFF(:,:)
 88:          DOUBLE PRECISION :: GRADATOMS(3*NATOMS) 88:          DOUBLE PRECISION :: GRADATOMS(3*NATOMS)
 89:          TYPE(POT_ENE_REC_C) :: ENERGY_DECOMP 89:          TYPE(POT_ENE_REC_C) :: ENERGY_DECOMP
 90:          SAVE 90:          SAVE
 91:  91: 
 92: ! #ifndef _OPTIMLIBRARY 
 93:         IF (RIGIDMOLECULEST) THEN   ! Should this be changed to IF(RIGIDINIT)? 92:         IF (RIGIDMOLECULEST) THEN   ! Should this be changed to IF(RIGIDINIT)?
 94:             ! sn402: addition to allow generalised rigid bodies to be used for arbitrary potential. 93:             ! sn402: addition to allow generalised rigid bodies to be used for arbitrary potential.
 95:             IF (ATOMRIGIDCOORDT .EQV. .FALSE.) THEN 94:             IF (ATOMRIGIDCOORDT .EQV. .FALSE.) THEN
 96:                 IF (DEBUG) THEN 95:                 IF (DEBUG) THEN
 97:                     IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .GT. 1.0E-12)) THEN 96:                     IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .GT. 1.0E-12)) THEN
 98:                         WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = FALSE but coords appear to be atomistic" 97:                         WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = FALSE but coords appear to be atomistic"
 99: !                        WRITE(*,*) "Largest absolute value:", MAXVAL(ABS(COORDS(DEGFREEDOMS+1:))) 98: !                        WRITE(*,*) "Largest absolute value:", MAXVAL(ABS(COORDS(DEGFREEDOMS+1:)))
100:                         WRITE(*,*) COORDS(DEGFREEDOMS+1:) 99:                         WRITE(*,*) COORDS(DEGFREEDOMS+1:)
101: !                         WRITE(*,*) COORDS100: !                         WRITE(*,*) COORDS
102:                     ENDIF101:                     ENDIF
104:                 XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)103:                 XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)
105:                 CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)104:                 CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)
106:             ELSE IF (DEBUG) THEN105:             ELSE IF (DEBUG) THEN
107:                 IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .LT. 1.0D-12)) THEN106:                 IF(ANY(ABS(COORDS(DEGFREEDOMS+1:)) .LT. 1.0D-12)) THEN
108:                     WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = TRUE but coords appear to be angleaxis"107:                     WRITE(*,*) "Warning: called POTENTIAL with ATOMRIGIDCOORDT = TRUE but coords appear to be angleaxis"
109: !                    WRITE(*,*) "Smallest absolute value:", MINVAL(ABS(COORDS(DEGFREEDOMS+1:)))108: !                    WRITE(*,*) "Smallest absolute value:", MINVAL(ABS(COORDS(DEGFREEDOMS+1:)))
110:                     WRITE(*,*) COORDS(DEGFREEDOMS+1:)109:                     WRITE(*,*) COORDS(DEGFREEDOMS+1:)
111:                 ENDIF110:                 ENDIF
112:             ENDIF111:             ENDIF
113:         ENDIF112:         ENDIF
114: ! #endif 
115:          ! 113:          ! 
116:          ! SSTEST needs to be true if we really want an analytic Hessian this step.114:          ! SSTEST needs to be true if we really want an analytic Hessian this step.
117:          ! 115:          ! 
118:          CALL MYCPU_TIME(TIME0,.FALSE.)116:          CALL MYCPU_TIME(TIME0,.FALSE.)
119:          NPCALL=NPCALL+1117:          NPCALL=NPCALL+1
120:          CALL CHANGEP118:          CALL CHANGEP
121:          SSTEST=.FALSE.119:          SSTEST=.FALSE.
122:          IF (STEST.AND.(.NOT.HUPDATE)) SSTEST=.TRUE.120:          IF (STEST.AND.(.NOT.HUPDATE)) SSTEST=.TRUE.
123:          IF (STEST.AND.(NHUP+1.EQ.NSTHUP)) SSTEST=.TRUE.121:          IF (STEST.AND.(NHUP+1.EQ.NSTHUP)) SSTEST=.TRUE.
124:          IF (STEST.AND.((INTHUP.GT.0).AND.(NHUP+1.GT.NSTHUP))) THEN122:          IF (STEST.AND.((INTHUP.GT.0).AND.(NHUP+1.GT.NSTHUP))) THEN
238:                WRITE(*,'(A)') ' potential> READHESS currently only works for CADPAC, GAMESS-US and GAMESS-UK run types'236:                WRITE(*,'(A)') ' potential> READHESS currently only works for CADPAC, GAMESS-US and GAMESS-UK run types'
239:             ENDIF237:             ENDIF
240:             CLOSE(15)238:             CLOSE(15)
241:          ENDIF239:          ENDIF
242:          ! 240:          ! 
243:          ! Spherical container241:          ! Spherical container
244:          ! 242:          ! 
245:          IF (CONTAINER) CALL RAD(COORDS)243:          IF (CONTAINER) CALL RAD(COORDS)
246:          ! IF (CONTAINER) CALL RAD(COORDS,ENERGY,VNEW,GTEST)244:          ! IF (CONTAINER) CALL RAD(COORDS,ENERGY,VNEW,GTEST)
247: 245: 
248: ! #ifndef _OPTIMLIBRARY 
249:          IF (RINGPOLYMERT) THEN246:          IF (RINGPOLYMERT) THEN
250:             ! 247:             ! 
251:             ! Get the energy and derivatives corresponding to the true potential for each bead and248:             ! Get the energy and derivatives corresponding to the true potential for each bead and
252:             ! construct the total energy, gradient and (if necessary) Hessian.249:             ! construct the total energy, gradient and (if necessary) Hessian.
253:             ! 250:             ! 
254:             IF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'AECK') THEN ! asymmetric Eckart barrier251:             IF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'AECK') THEN ! asymmetric Eckart barrier
255:                NZERO=0252:                NZERO=0
256:                ENERGY=0.0D0253:                ENERGY=0.0D0
257:                IF (GTEST) VNEW(1:NOPT)=0.0D0254:                IF (GTEST) VNEW(1:NOPT)=0.0D0
258:                IF (STEST) HESS(1:NOPT,1:NOPT)=0.0D0255:                IF (STEST) HESS(1:NOPT,1:NOPT)=0.0D0
354:      &               FINDIFHESS_POT(COORDS(RPDOF*(J1-1)+1:RPDOF*J1),MCYPOT,1.0D-3)351:      &               FINDIFHESS_POT(COORDS(RPDOF*(J1-1)+1:RPDOF*J1),MCYPOT,1.0D-3)
355:                   ENDDO352:                   ENDDO
356:                ELSEIF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'JB') THEN ! Joel Bowman's water potential353:                ELSEIF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'JB') THEN ! Joel Bowman's water potential
357:                   DUMMY2=1.0D0/(RPBETA/RPIMAGES)**2 ! atomic units (hbar = 1)354:                   DUMMY2=1.0D0/(RPBETA/RPIMAGES)**2 ! atomic units (hbar = 1)
358:                   DO J1=1,RPIMAGES355:                   DO J1=1,RPIMAGES
359:                      ENERGY=ENERGY+BOWMANPOT(COORDS(RPDOF*(J1-1)+1:RPDOF*J1)) ! Hartrees356:                      ENERGY=ENERGY+BOWMANPOT(COORDS(RPDOF*(J1-1)+1:RPDOF*J1)) ! Hartrees
360:                      IF (GTEST) VNEW(RPDOF*(J1-1)+1:RPDOF*J1)=FINDIFGRAD(COORDS(RPDOF*(J1-1)+1:RPDOF*J1), BOWMANPOT, 1.0D-3, GRAD4T)357:                      IF (GTEST) VNEW(RPDOF*(J1-1)+1:RPDOF*J1)=FINDIFGRAD(COORDS(RPDOF*(J1-1)+1:RPDOF*J1), BOWMANPOT, 1.0D-3, GRAD4T)
361:                      IF (STEST) HESS(RPDOF*(J1-1)+1:RPDOF*J1,RPDOF*(J1-1)+1:RPDOF*J1) =358:                      IF (STEST) HESS(RPDOF*(J1-1)+1:RPDOF*J1,RPDOF*(J1-1)+1:RPDOF*J1) =
362:      &               FINDIFHESS_POT(COORDS(RPDOF*(J1-1)+1:RPDOF*J1),BOWMANPOT,1.0D-3)359:      &               FINDIFHESS_POT(COORDS(RPDOF*(J1-1)+1:RPDOF*J1),BOWMANPOT,1.0D-3)
363:                   ENDDO360:                   ENDDO
364:                ELSEIF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'MB') THEN ! MBPOL 
365:                   DUMMY2=1.0D0/(RPBETA/RPIMAGES)**2 ! atomic units (hbar = 1) 
366:                   DO J1=1,RPIMAGES 
367:                      CALL MBPOL(COORDS(RPDOF*(J1-1)+1:RPDOF*J1),EDUMMY,VNEW(RPDOF*(J1-1)+1:RPDOF*J1),GTEST) 
368:                      ENERGY=ENERGY+EDUMMY 
369:                      IF (STEST) THEN 
370:                         PRINT *, 'no hessian for MBPOL' 
371:                         STOP 
372:                      ENDIF 
373:                   ENDDO 
374:  
375:                ELSEIF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'MAL') THEN ! malonaldehyde361:                ELSEIF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'MAL') THEN ! malonaldehyde
376:                   DUMMY2=1.0D0/(RPBETA/RPIMAGES)**2 ! atomic units (hbar = 1)362:                   DUMMY2=1.0D0/(RPBETA/RPIMAGES)**2 ! atomic units (hbar = 1)
377:                   DO J1=1,RPIMAGES363:                   DO J1=1,RPIMAGES
378:                      CALL MALPES(9, COORDS(RPDOF*(J1-1)+1:RPDOF*J1),  VNEW(RPDOF*(J1-1)+1:RPDOF*J1), DUMMY, GTEST, SSTEST)364:                      CALL MALPES(9, COORDS(RPDOF*(J1-1)+1:RPDOF*J1),  VNEW(RPDOF*(J1-1)+1:RPDOF*J1), DUMMY, GTEST, SSTEST)
379:                      ENERGY=ENERGY+DUMMY ! Hartrees365:                      ENERGY=ENERGY+DUMMY ! Hartrees
380:                      IF (SSTEST) THEN366:                      IF (SSTEST) THEN
381:                         PRINT *, 'no hessian for TTM3-F'367:                         PRINT *, 'no hessian for TTM3-F'
382:                         STOP368:                         STOP
383:                      ENDIF369:                      ENDIF
384:                   ENDDO370:                   ENDDO
3712:             ENDIF3698:             ENDIF
3713:          ELSE IF (ZSYM(NATOMS).EQ.'SS') THEN3699:          ELSE IF (ZSYM(NATOMS).EQ.'SS') THEN
3714:             ! param1 is the soft sphere power3700:             ! param1 is the soft sphere power
3715:             CALL SOFT_SPHERE_POTENTIAL(NATOMS,COORDS,VNEW,HESS,ENERGY, 3701:             CALL SOFT_SPHERE_POTENTIAL(NATOMS,COORDS,VNEW,HESS,ENERGY, 
3716:      &              BULK_BOXVEC, SOFT_SPHERE_RADII, PARAM1, GTEST,SSTEST)3702:      &              BULK_BOXVEC, SOFT_SPHERE_RADII, PARAM1, GTEST,SSTEST)
3717:             IF (PTEST) THEN3703:             IF (PTEST) THEN
3718:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY3704:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY
3719:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY3705:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY
3720:             ENDIF3706:             ENDIF
3721:             ! davidg: introduced userpot here3707:             ! davidg: introduced userpot here
3722: ! 
3723: ! USERPOT is the only potential required for the Biovia build. 
3724: ! 
3725:          ELSE IF (USERPOTT) THEN3708:          ELSE IF (USERPOTT) THEN
3726: ! #endif 
3727: ! #ifdef _OPTIMLIBRARY 
3728: !        IF (USERPOTT) THEN 
3729: ! #endif 
3730:             CALL USERPOT_POTENTIALHESS(3*NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST,HESS)3709:             CALL USERPOT_POTENTIALHESS(3*NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST,HESS)
3731:             IF (PTEST) THEN3710:             IF (PTEST) THEN
3732:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' meV'3711:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' meV'
3733:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' meV'3712:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' meV'
3734:             ENDIF3713:             ENDIF
3735: ! #ifdef _OPTIMLIBRARY3714: 
3736: !        ENDIF 
3737: ! #else 
3738:             ELSE IF (STEALTHYT) THEN3715:             ELSE IF (STEALTHYT) THEN
3739:                CALL STEALTHY(COORDS, VNEW, ENERGY, GTEST, SSTEST)3716:                CALL STEALTHY(COORDS, VNEW, ENERGY, GTEST, SSTEST)
3740:  3717:  
3741:                !DIFF=1.0D-63718:                !DIFF=1.0D-6
3742:                !OPEN(UNIT=1500, FILE="grad-test")3719:                !OPEN(UNIT=1500, FILE="grad-test")
3743:                !WRITE(1500,*) 'ANALYTIC AND NUMERICAL GRADIENTS:NATOMS=',NATOMS3720:                !WRITE(1500,*) 'ANALYTIC AND NUMERICAL GRADIENTS:NATOMS=',NATOMS
3744:                !DO J1=1, 3*NATOMS3721:                !DO J1=1, 3*NATOMS
3745:                !   COORDS(J1)=COORDS(J1)+DIFF3722:                !   COORDS(J1)=COORDS(J1)+DIFF
3746:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)3723:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)
3747:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF3724:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF
3801:             ! IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND.3778:             ! IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND.
3802:             ! 1             (ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D0)) THEN3779:             ! 1             (ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D0)) THEN
3803:             ! WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'3780:             ! WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'
3804:             ! ELSE3781:             ! ELSE
3805:             ! WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)3782:             ! WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)
3806:             ! ENDIF3783:             ! ENDIF
3807:             ! ENDDO3784:             ! ENDDO
3808:             ! ENDDO3785:             ! ENDDO
3809:             ! DC430 >3786:             ! DC430 >
3810:          ENDIF3787:          ENDIF
3811: !   #endif 
3812:          ! 3788:          ! 
3813:          ! End of possible potentials - now add extra terms required, if any ------------------------------3789:          ! End of possible potentials - now add extra terms required, if any ------------------------------
3814:          !3790:          !
3815:          !ds656> Apply a substrate field for keyword MIE_FIELD3791:          !ds656> Apply a substrate field for keyword MIE_FIELD
3816:           
3817: ! #ifndef _OPTIMLIBRARY          
3818: !ds656> Test gradient and Hessian with MIEF 
3819:          IF(MIEFT) CALL MIEF(NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST)3792:          IF(MIEFT) CALL MIEF(NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST)
 3793:          
 3794:          !ds656> Test gradient and Hessian with MIEF
3820:          IF (ALLOCATED(HESS).AND.MIEFT.AND..FALSE.) THEN3795:          IF (ALLOCATED(HESS).AND.MIEFT.AND..FALSE.) THEN
3821:          DIFF=1.0D-43796:          DIFF=1.0D-4
3822:          PRINT*,'analytic and numerical gradients:'3797:          PRINT*,'analytic and numerical gradients:'
3823:          DO J1=1,3*NATOMS3798:          DO J1=1,3*NATOMS
3824:             COORDS(J1)=COORDS(J1)+DIFF3799:             COORDS(J1)=COORDS(J1)+DIFF
3825:             CALL BGUPTA(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)3800:             CALL BGUPTA(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)
3826:             CALL MIEF(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)3801:             CALL MIEF(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)
3827:             COORDS(J1)=COORDS(J1)-2.0D0*DIFF3802:             COORDS(J1)=COORDS(J1)-2.0D0*DIFF
3828:             CALL BGUPTA(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)3803:             CALL BGUPTA(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)
3829:             CALL MIEF(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)3804:             CALL MIEF(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)
3946:             ! ENDDO3921:             ! ENDDO
3947:             ! ENDDO3922:             ! ENDDO
3948:             ! STOP3923:             ! STOP
3949:          ENDIF3924:          ENDIF
3950:          ! 3925:          ! 
3951:          ! Add on attractive term for "closest" minimum if this is a REDOPATH run3926:          ! Add on attractive term for "closest" minimum if this is a REDOPATH run
3952:          ! Here we have to define closest in terms of which minimum we have moved3927:          ! Here we have to define closest in terms of which minimum we have moved
3953:          ! towards, not in terms of absolute distance, because of asymmetric pathways.3928:          ! towards, not in terms of absolute distance, because of asymmetric pathways.
3954:          ! Avoid division by zero for D1INIT and D2INIT!3929:          ! Avoid division by zero for D1INIT and D2INIT!
3955:          ! 3930:          ! 
3956: ! #endif 
3957:          IF (REDOKADD.AND.REDOPATH.AND.(.NOT.REDOPATHXYZ).AND.(REDOK.NE.0.0D0).AND.3931:          IF (REDOKADD.AND.REDOPATH.AND.(.NOT.REDOPATHXYZ).AND.(REDOK.NE.0.0D0).AND.
3958:      &   (ALLOCATED(MIN1REDO)).AND.(ALLOCATED(MIN2REDO)).AND.(D1INIT*D2INIT.NE.0.0D0)) THEN3932:      &   (ALLOCATED(MIN1REDO)).AND.(ALLOCATED(MIN2REDO)).AND.(D1INIT*D2INIT.NE.0.0D0)) THEN
3959: 3933: 
3960:             CALL NEWMINDIST(COORDS,MIN1REDO,NATOMS,DIST1,BULKT,TWOD,'AX   ',.FALSE.,RIGIDBODY,DEBUG,RMAT)3934:             CALL NEWMINDIST(COORDS,MIN1REDO,NATOMS,DIST1,BULKT,TWOD,'AX   ',.FALSE.,RIGIDBODY,DEBUG,RMAT)
3961:             CALL NEWMINDIST(COORDS,MIN2REDO,NATOMS,DIST2,BULKT,TWOD,'AX   ',.FALSE.,RIGIDBODY,DEBUG,RMAT)3935:             CALL NEWMINDIST(COORDS,MIN2REDO,NATOMS,DIST2,BULKT,TWOD,'AX   ',.FALSE.,RIGIDBODY,DEBUG,RMAT)
3962: 3936: 
3963:             IF     ((DIST1/D1INIT.LT.REDOFRAC).AND.REDOPATH1) THEN3937:             IF     ((DIST1/D1INIT.LT.REDOFRAC).AND.REDOPATH1) THEN
3964:                DUMMY1=0.0D03938:                DUMMY1=0.0D0
3965:                DUMMY2=0.0D03939:                DUMMY2=0.0D0
3966:                DUMMY3=0.0D03940:                DUMMY3=0.0D0
3991:                DUMMY1=DUMMY1/SQRT(DUMMY3)3965:                DUMMY1=DUMMY1/SQRT(DUMMY3)
3992:                IF (DEBUG) PRINT '(A,2G15.5,A,G20.10)',' potential> distances/initial distance are ',DIST1/D1INIT,DIST2/D2INIT,3966:                IF (DEBUG) PRINT '(A,2G15.5,A,G20.10)',' potential> distances/initial distance are ',DIST1/D1INIT,DIST2/D2INIT,
3993:      &         ' grad % for second minimum=',DUMMY1*100/SQRT(DUMMY2)3967:      &         ' grad % for second minimum=',DUMMY1*100/SQRT(DUMMY2)
3994:                ! IF (DUMMY1.GT.0.0D0) DUMMY1=-DUMMY13968:                ! IF (DUMMY1.GT.0.0D0) DUMMY1=-DUMMY1
3995:                DO J1=1,NOPT3969:                DO J1=1,NOPT
3996:                   VNEW(J1)=VNEW(J1)+REDOK*(COORDS(J1)-MIN2REDO(J1))3970:                   VNEW(J1)=VNEW(J1)+REDOK*(COORDS(J1)-MIN2REDO(J1))
3997:                   ! VNEW(J1)=VNEW(J1)+REDOK*(COORDS(J1)-MIN2REDO(J1))*SQRT(DUMMY2/DUMMY3)3971:                   ! VNEW(J1)=VNEW(J1)+REDOK*(COORDS(J1)-MIN2REDO(J1))*SQRT(DUMMY2/DUMMY3)
3998:                ENDDO3972:                ENDDO
3999:             ENDIF3973:             ENDIF
4000:          ENDIF3974:          ENDIF
4001: ! #ifndef _OPTIMLIBRARY 
4002:          ! 3975:          ! 
4003:          ! Add on terms for rotation about the z axis3976:          ! Add on terms for rotation about the z axis
4004:          ! 3977:          ! 
4005:          IF (RTEST) THEN3978:          IF (RTEST) THEN
4006:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'3979:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'
4007:             PRINT*,' WARNING - GTEST and SSTEST ignored'3980:             PRINT*,' WARNING - GTEST and SSTEST ignored'
4008:             IF (JZ.NE.0.0D0) THEN3981:             IF (JZ.NE.0.0D0) THEN
4009:                CALL ROTD(NATOMS, COORDS, VNEW, 1.0D0, JZ, .FALSE., ROT)3982:                CALL ROTD(NATOMS, COORDS, VNEW, 1.0D0, JZ, .FALSE., ROT)
4010:             ELSE3983:             ELSE
4011:                CALL ROTENERGY(NATOMS, COORDS, OMEGA, 1.0D0, IZ, ROT)3984:                CALL ROTENERGY(NATOMS, COORDS, OMEGA, 1.0D0, IZ, ROT)
4068:             IF (ETEST) THEN4041:             IF (ETEST) THEN
4069:                OPEN(UNIT=91,FILE='abenergy',STATUS='OLD')4042:                OPEN(UNIT=91,FILE='abenergy',STATUS='OLD')
4070:                READ(91,*) ENERGY4043:                READ(91,*) ENERGY
4071:                IF (PTEST) WRITE(*,'(A,27X,F20.10,A)') ' potential> Energy for last cycle=',ENERGY,' hartree'4044:                IF (PTEST) WRITE(*,'(A,27X,F20.10,A)') ' potential> Energy for last cycle=',ENERGY,' hartree'
4072:                CLOSE(91)4045:                CLOSE(91)
4073:             ELSE4046:             ELSE
4074:                WRITE(*,'(A)') ' potential> Error - abenergy file not found'4047:                WRITE(*,'(A)') ' potential> Error - abenergy file not found'
4075:                STOP4048:                STOP
4076:             ENDIF4049:             ENDIF
4077:          ENDIF4050:          ENDIF
4078: ! #endif 
4079: 4051: 
4080:          IF (GTEST) THEN4052:          IF (GTEST) THEN
4081:             ! PRINT '(A,L5)',' potential> FREEZE,VNEW=',FREEZE4053:             ! PRINT '(A,L5)',' potential> FREEZE,VNEW=',FREEZE
4082:             ! PRINT '(3G20.10)',VNEW(1:3*NATOMS)4054:             ! PRINT '(3G20.10)',VNEW(1:3*NATOMS)
4083:             IF (FREEZE) THEN4055:             IF (FREEZE) THEN
4084:                DO J1=1,NATOMS4056:                DO J1=1,NATOMS
4085:                   IF (FROZEN(J1)) THEN4057:                   IF (FROZEN(J1)) THEN
4086:                      IF (VARIABLES) THEN4058:                      IF (VARIABLES) THEN
4087:                         VNEW(J1)=0.0D04059:                         VNEW(J1)=0.0D0
4088:                      ELSE4060:                      ELSE
4089:                         VNEW(3*(J1-1)+1)=0.0D04061:                         VNEW(3*(J1-1)+1)=0.0D0
4090:                         VNEW(3*(J1-1)+2)=0.0D04062:                         VNEW(3*(J1-1)+2)=0.0D0
4091:                         VNEW(3*(J1-1)+3)=0.0D04063:                         VNEW(3*(J1-1)+3)=0.0D0
4092:                      ENDIF4064:                      ENDIF
4093:                   ENDIF4065:                   ENDIF
4094:                ENDDO4066:                ENDDO
4095:             ENDIF4067:             ENDIF
4096: ! #ifndef _OPTIMLIBRARY4068: 
4097:             ! sn402: other half of hack to allow genrigid to work with arbitrary potential4069:             ! sn402: other half of hack to allow genrigid to work with arbitrary potential
4098:             IF (RIGIDMOLECULEST .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN4070:             IF (RIGIDMOLECULEST .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN
4099:                 IF (STEST) THEN4071:                 IF (STEST) THEN
4100:                      CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)4072:                      CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)
4101:                          HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D04073:                          HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0
4102:                          HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D04074:                          HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0
4103:                          HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)4075:                          HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)
4104:                 ENDIF4076:                 ENDIF
4105:                 CALL TRANSFORMGRAD(VNEW, XRIGIDCOORDS, XRIGIDGRAD)4077:                 CALL TRANSFORMGRAD(VNEW, XRIGIDCOORDS, XRIGIDGRAD)
4106:                 VNEW(1:DEGFREEDOMS) = XRIGIDGRAD(1:DEGFREEDOMS)4078:                 VNEW(1:DEGFREEDOMS) = XRIGIDGRAD(1:DEGFREEDOMS)
4107:                 VNEW(DEGFREEDOMS+1:3*NATOMS) = 0.0D04079:                 VNEW(DEGFREEDOMS+1:3*NATOMS) = 0.0D0
4108:             ENDIF4080:             ENDIF
4109: ! #endif 
4110: 4081: 
4111:             IF (UNRST) THEN4082:             IF (UNRST) THEN
4112:                CALL VSTAT(VNEW,TEMP,NINTS,NOPT)4083:                CALL VSTAT(VNEW,TEMP,NINTS,NOPT)
4113:                IF (PTEST) WRITE(*,'(A,43X,F15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',4084:                IF (PTEST) WRITE(*,'(A,43X,F15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',
4114:      &         TEMP(5)*SQRT(1.0D0*(NINTS))4085:      &         TEMP(5)*SQRT(1.0D0*(NINTS))
4115:             ELSE IF ( RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN !hk2864086:             ELSE IF ( RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN !hk286
4116:                CALL VSTAT(VNEW(1:DEGFREEDOMS),TEMP,DEGFREEDOMS,DEGFREEDOMS)4087:                CALL VSTAT(VNEW(1:DEGFREEDOMS),TEMP,DEGFREEDOMS,DEGFREEDOMS)
4117:                IF (PTEST) WRITE(*,'(A,43X,F15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',4088:                IF (PTEST) WRITE(*,'(A,43X,F15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',
4118:      &         TEMP(5)*SQRT(1.0D0*(DEGFREEDOMS))4089:      &         TEMP(5)*SQRT(1.0D0*(DEGFREEDOMS))
4119:                IF (AACONVERGENCET .AND. (TEMP(5) < 5.0D0 * GMAX)) THEN4090:                IF (AACONVERGENCET .AND. (TEMP(5) < 5.0D0 * GMAX)) THEN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0