hdiff output

r32288/key.f90 2017-04-19 12:30:14.682765596 +0100 r32287/key.f90 2017-04-19 12:30:15.578777359 +0100
 47:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, & 47:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, &
 48:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, & 48:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, &
 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 57:      &        MALONALDEHYDE, SIO2PT, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, & 57:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, &
 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, & 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, &
 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T
 60:  60: 
 61: ! sy349 > for testing the flatpath after dneb 61: ! sy349 > for testing the flatpath after dneb
 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:)
 63:       LOGICAL FLATPATHT 63:       LOGICAL FLATPATHT
 64:  64: 
 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)


r32288/keywords.f 2017-04-19 12:30:14.914768642 +0100 r32287/keywords.f 2017-04-19 12:30:15.806780348 +0100
196: 196: 
197:          DESMAXEJUMP = HUGE(1.0D0)197:          DESMAXEJUMP = HUGE(1.0D0)
198:          DESMAXAVGE = HUGE(1.0D0)198:          DESMAXAVGE = HUGE(1.0D0)
199:          UNSTRING='UNDEFINED'199:          UNSTRING='UNDEFINED'
200:          WELCH=.FALSE.200:          WELCH=.FALSE.
201:          TOSI=.FALSE.201:          TOSI=.FALSE.
202:          TOSIC6=.FALSE.202:          TOSIC6=.FALSE.
203:          SIO2T=.FALSE.203:          SIO2T=.FALSE.
204:          SIO2C6T=.FALSE.204:          SIO2C6T=.FALSE.
205:          TOSIPOL=.FALSE.205:          TOSIPOL=.FALSE.
206:          SIO2PT=.FALSE. 
207:          ALLOCATE(TAGFAC(NATOMS),TAGNUM(NATOMS))206:          ALLOCATE(TAGFAC(NATOMS),TAGNUM(NATOMS))
208:          TAGFAC(1:NATOMS)=1.0D0207:          TAGFAC(1:NATOMS)=1.0D0
209:          TAGNUM(1:NATOMS)=0208:          TAGNUM(1:NATOMS)=0
210:          NTAG=0209:          NTAG=0
211:          TAGT=.FALSE.210:          TAGT=.FALSE.
212: 211: 
213:          REPEL=.FALSE.212:          REPEL=.FALSE.
214: 213: 
215:          INR=-1214:          INR=-1
216: 215: 
6243:                CALL READF(PARAM3)6242:                CALL READF(PARAM3)
6244:             ENDIF6243:             ENDIF
6245:             IF (NITEMS.GT.4) THEN6244:             IF (NITEMS.GT.4) THEN
6246:                CALL READF(PARAM4)6245:                CALL READF(PARAM4)
6247:             ENDIF6246:             ENDIF
6248:          ELSE IF (WORD.EQ.'SIO2C6') THEN6247:          ELSE IF (WORD.EQ.'SIO2C6') THEN
6249:             SIO2C6T=.TRUE.6248:             SIO2C6T=.TRUE.
6250:             CALL READF(C6PP)6249:             CALL READF(C6PP)
6251:             CALL READF(C6MM)6250:             CALL READF(C6MM)
6252:             CALL READF(C6PM)6251:             CALL READF(C6PM)
6253:          ELSE IF (WORD.EQ.'SIO2P') THEN6252: 
6254:             SIO2PT=.TRUE. 
6255:          ! Activate the iSLERP interpolation scheme for rigid bodies. Only has any effect6253:          ! Activate the iSLERP interpolation scheme for rigid bodies. Only has any effect
6256:          ! if RIGIDINIT is also set.6254:          ! if RIGIDINIT is also set.
6257:          ELSE IF (WORD.EQ.'SLERP') THEN6255:          ELSE IF (WORD.EQ.'SLERP') THEN
6258:              SLERPT=.TRUE.           6256:             SLERPT=.TRUE.           
 6257:             
6259: ! 6258: ! 
6260: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using6259: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using
6261: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations6260: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations
6262: ! or if the RMS force goes below SQVVGuessRMSTol.6261: ! or if the RMS force goes below SQVVGuessRMSTol.
6263: ! 6262: ! 
6264:          ELSE IF (WORD == 'SQVV') THEN6263:          ELSE IF (WORD == 'SQVV') THEN
6265:             SQVVGUESS=.TRUE.6264:             SQVVGUESS=.TRUE.
6266:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)6265:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)
6267:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)6266:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)
6268:          ELSE IF (WORD.EQ.'SSH') THEN6267:          ELSE IF (WORD.EQ.'SSH') THEN
6269:             SSHT=.TRUE.6268:             SSHT=.TRUE.
6270: !6269: !
6271: ! KLIM sets the radius in then wave vector space.6270: ! KLIM sets the radius in then wave vector space.
6272: ! SCA sets the multiplier on the stealthy potential and gradients.6271: ! SCA sets the multiplier on the stealthy potential and gradients.
6273: !6272: !
6274:          ELSE IF (WORD.EQ.'STEALTHY') THEN6273:          ELSE IF (WORD.EQ.'STEALTHY') THEN
6275:              STEALTHYT=.TRUE.6274:             STEALTHYT=.TRUE.
6276:              CALL READF(KLIM)6275:             CALL READF(KLIM)
6277:              IF (NITEMS.GT.2) THEN6276:             IF (NITEMS.GT.2) THEN
6278:                 CALL READF(SCA)6277:                CALL READF(SCA)
6279:              ELSE6278:             ELSE
6280:                 SCA=16279:                SCA=1
6281:              END IF6280:             END IF
6282: 6281: 
6283:          ELSE IF (WORD.EQ.'STEALTHYTEST') THEN6282:          ELSE IF (WORD.EQ.'STEALTHYTEST') THEN
6284:              STEALTV=.TRUE.6283:             STEALTV=.TRUE.
6285: ! 6284: ! 
6286: ! NSTEPMIN sets the minimum number of steps allowed before convergence.6285: ! NSTEPMIN sets the minimum number of steps allowed before convergence.
6287: ! 6286: ! 
6288:          ELSE IF (WORD .EQ. 'STEPMIN') THEN6287:          ELSE IF (WORD .EQ. 'STEPMIN') THEN
6289:              CALL READI(NSTEPMIN)6288:             CALL READI(NSTEPMIN)
6290: ! 6289: ! 
6291: ! STEPS n sets the number of optimisation steps to perform6290: ! STEPS n sets the number of optimisation steps to perform
6292: ! per call to OPTIM                                    - default n=16291: ! per call to OPTIM                                    - default n=1
6293: ! If BFGSSTEPS is not specified then it is set to the same value as NSTEPS6292: ! If BFGSSTEPS is not specified then it is set to the same value as NSTEPS
6294: ! 6293: ! 
6295:          ELSE IF (WORD .EQ. 'STEPS') THEN6294:          ELSE IF (WORD .EQ. 'STEPS') THEN
6296:             CALL READI(NSTEPS)6295:             CALL READI(NSTEPS)
6297:             IF (BFGSSTEPS.EQ.1) BFGSSTEPS=NSTEPS6296:             IF (BFGSSTEPS.EQ.1) BFGSSTEPS=NSTEPS
6298: ! 6297: ! 
6299: ! Stillinger-David water potential - coded by Jeremy Richardson6298: ! Stillinger-David water potential - coded by Jeremy Richardson


r32288/potential.f 2017-04-19 12:30:15.134771529 +0100 r32287/potential.f 2017-04-19 12:30:16.034783340 +0100
970:             IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' hartree'970:             IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' hartree'
971:             IF (PTEST) WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' hartree'971:             IF (PTEST) WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' hartree'
972:             ! ELSE IF (ZSYM(NATOMS).EQ.'CL') THEN972:             ! ELSE IF (ZSYM(NATOMS).EQ.'CL') THEN
973:             ! PRINT*,' WARNING - GTEST and SSTEST ignored'973:             ! PRINT*,' WARNING - GTEST and SSTEST ignored'
974:             ! ELSE IF (ZSYM(NATOMS).EQ.'CL') THEN974:             ! ELSE IF (ZSYM(NATOMS).EQ.'CL') THEN
975:             ! PRINT*,' WARNING - GTEST and SSTEST ignored'975:             ! PRINT*,' WARNING - GTEST and SSTEST ignored'
976:             ! CALL KDIFF(NATOMS, COORDS, VNEW, ENERGY)976:             ! CALL KDIFF(NATOMS, COORDS, VNEW, ENERGY)
977:             ! CALL KPAIRS(NATOMS, 3*NATOMS, COORDS, ENERGY)977:             ! CALL KPAIRS(NATOMS, 3*NATOMS, COORDS, ENERGY)
978:             ! WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' hartree'978:             ! WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' hartree'
979:             ! WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' hartree'979:             ! WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' hartree'
980:          ELSE IF (SIO2PT) THEN 
981:              CALL SIO2PSHIFT(NATOMS, COORDS, VNEW, ENERGY, PARAM1,PARAM2, PARAM3, PARAM4, GTEST, SSTEST, PTEST, BOXTEST) 
982:          ELSE IF (ZSYM(NATOMS).EQ.'AZ') THEN980:          ELSE IF (ZSYM(NATOMS).EQ.'AZ') THEN
983:             PRINT*,' WARNING - GTEST and SSTEST ignored'981:             PRINT*,' WARNING - GTEST and SSTEST ignored'
984:             CALL AZIZ(NATOMS,COORDS,VNEW,ENERGY,1)982:             CALL AZIZ(NATOMS,COORDS,VNEW,ENERGY,1)
985:             IF (PARAM1.NE.0.0D0) THEN983:             IF (PARAM1.NE.0.0D0) THEN
986:                ZSTAR=PARAM1984:                ZSTAR=PARAM1
987:                CALL AXDIFF(NATOMS, COORDS, VNEW, ZSTAR, GTEST, SSTEST)985:                CALL AXDIFF(NATOMS, COORDS, VNEW, ZSTAR, GTEST, SSTEST)
988:                CALL AXPAIRS (NATOMS, COORDS, P2, P3, TEMPA, ZSTAR)986:                CALL AXPAIRS (NATOMS, COORDS, P2, P3, TEMPA, ZSTAR)
989:                P2=ENERGY987:                P2=ENERGY
990:                ENERGY=ENERGY+P3988:                ENERGY=ENERGY+P3
991:             ENDIF989:             ENDIF


r32288/SiO2pshift.f 2017-04-19 12:30:14.466762761 +0100 r32287/SiO2pshift.f 2017-04-19 12:30:15.366774572 +0100
  1: C   OPTIM: A program for optimizing geometries and calculating reaction pathways  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/SiO2pshift.f' in revision 32287
  2: C   Copyright (C) 1999-2006 David J. Wales 
  3: C   This file is part of OPTIM. 
  4: C 
  5: C   OPTIM is free software; you can redistribute it and/or modify 
  6: C   it under the terms of the GNU General Public License as published by 
  7: C   the Free Software Foundation; either version 2 of the License, or 
  8: C   (at your option) any later version. 
  9: C 
 10: C   OPTIM is distributed in the hope that it will be useful, 
 11: C   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: C   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: C   GNU General Public License for more details. 
 14: C 
 15: C   You should have received a copy of the GNU General Public License 
 16: C   along with this program; if not, write to the Free Software 
 17: C   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: C 
 19: C 
 20: C************************************************************************* 
 21: C 
 22: C  Subroutine Si02PSHIFT calculates the energy, cartesian gradient and second 
 23: C  derivative matrix analytically for Si02 using a shifted, truncated potential. 
 24: C 
 25: C  Adapted for the Si02 glass described by van Beest, Kramer, van Santen, Physical Review Letters, 64, 1990.  
 26: C   
 27: C  Atom types are A(Si) and B(O). The first 
 28: C  NTYPEA are A, the next NTYPEB=NATOMS-NTYPEA are B.   
 29: C  (input must be in this order!) 
 30: C 
 31: C 
 32: C  The shifting and truncation for the short-range part of the potential is  
 33: C  as described by Stoddard and Ford, Phys. Rev. A, 8, 1504, 1973.  
 34: C  The Coulomb-part is shifted and truncated via the Wolf-method,  
 35: C  (see Carre et al., Jour. Chem. Phys., 127, 114512, 2007 and  
 36: C  Fennell and Gezelter, Jour. Chem. Phys., 124, 234104, 2006)  
 37: C   
 38: C 
 39: C************************************************************************* 
 40: C 
 41:       SUBROUTINE SIO2PSHIFT(N, X, V, POTEL, BOXLX, BOXLY, BOXLZ, CUTOFF, GTEST, STEST, PTEST, BOXTEST) 
 42:       USE KEY 
 43:       USE MODHESS  
 44:       IMPLICIT NONE 
 45:       INTEGER N, J1, J2, J3, J4, NTYPEA, ANV(N,N,3) 
 46:       DOUBLE PRECISION X(3*N), VEC1, VEC2, VEC3,  
 47:      1                 V(3*N), R6, R2DUM, 
 48:      2                 Q_A, Q_B, QQ_AA, QQ_AB, QQ_BB, A_AA, A_AB, A_BB, 
 49:      3                 B_AA, B_AB, B_BB, C_AA, C_AB, C_BB, EPS_AA, SIG_AA,  
 50:      4                 EPS_AB, EPS_BB, SIG_AB, SIG_BB, POTEL,  
 51:      5                 BOXLX, BOXLY, BOXLZ, IRCUT2AA, IRCUT2AB, IRCUT2BB,    
 52:      6                 CUTOFF, CONSTAA, CONSTBB, CONSTAB, 
 53:      7                 RCONSTAA, RCONSTAB, RCONSTBB, CUTAA, CUTAB, CUTBB, 
 54:      8                 SIG6_AB, SIG30_AB, SIG6_BB, SIG30_BB 
 55:       LOGICAL GTEST, STEST, PTEST, BOXTEST 
 56:       COMMON /BIN/ NTYPEA 
 57:       COMMON /RCONSTBIN/ RCONSTAA, RCONSTAB, RCONSTBB, IRCUT2AA, IRCUT2AB, IRCUT2BB 
 58:       parameter(QQ_AA=82.9428D0,QQ_AB=-41.4714D0, QQ_BB=20.7357D0) 
 59: C      parameter(A_AA=0.0D0,A_AB=661.62554D0,A_BB=51.03644D0) ! [A] = Eh   
 60: C      parameter(B_AA=0.0D0,B_AB=2.57877477D0,B_BB=1.4605285D0) ! [B] = 1/a 
 61: C      parameter(C_AA=0.0D0,C_AB=223.485086D0,C_BB=292.8744D0) ! [C] = Eh a^6  
 62:       parameter(A_AA=0.0D0,A_AB=18003.7572D0,A_BB=1388.7730D0) ! [A] = eV   
 63:       parameter(B_AA=0.0D0,B_AB=4.87318D0,B_BB=2.76000D0) ! [B] = 1/Angstrom 
 64:       parameter(C_AA=0.0D0,C_AB=133.5381D0,C_BB=175.0000D0) ! [C] = eV Angstrom^6 
 65:       parameter(EPS_AA=0.0D0, EPS_AB=0.003097948D0, EPS_BB=0.00105105048D0) 
 66:       parameter(SIG_AA=0.0D0, SIG_AB=1.313635D0, SIG_BB=1.779239D0) 
 67:  
 68:       IF(MOD(N,3).EQ.0) THEN 
 69:         NTYPEA = N/3 
 70:       ELSE 
 71:         WRITE(*,*) "SiO2pshift.f> warning: number of atoms must correspond to stoichiometry of SiO2!" 
 72:         STOP 
 73:       ENDIF 
 74:     
 75: C      IF (PTEST) THEN 
 76: C         WRITE(102,*) "SiO2pshift.f> warning: PTEST not implemented" 
 77: C         STOP 
 78: C      ENDIF 
 79:       IF (BOXTEST) THEN 
 80:          WRITE(102,*) "SiO2pshift.f> warning: box derivatives not implemented" 
 81:          STOP 
 82:       ENDIF 
 83:  
 84:       SIG30_AB = SIG_AB**30 
 85:       SIG6_AB = SIG_AB**6 
 86:       SIG30_BB = SIG_BB**30 
 87:       SIG6_BB = SIG_BB**6       
 88:  
 89:       CUTAA=CUTOFF  
 90:       CUTAB=CUTOFF ! change? and set CUTOFF! 
 91:       CUTBB=CUTOFF ! 
 92:       IRCUT2AA = 1.D0/CUTAA**2 
 93:       IRCUT2AB = 1.D0/CUTAB**2 
 94:       IRCUT2BB = 1.D0/CUTBB**2 
 95:        
 96:       CONSTAA = 0.0D0 
 97:       CONSTAB = -A_AB*EXP(-B_AB*CUTAB)*(1.0D0+B_AB*CUTAB/2.0D0) 
 98:      &          +4.0D0*C_AB*IRCUT2AB**3 
 99:      &          +4.0D0*EPS_AB*(-16.D0*SIG30_AB/CUTAB**30+4.0D0*SIG6_AB/CUTAB**6) 
100:       CONSTBB = -A_BB*EXP(-B_BB*CUTBB)*(1.0D0+B_BB*CUTBB/2.0D0) 
101:      &          +4.0D0*C_BB*IRCUT2BB**3 
102:      &          +4.0D0*EPS_BB*(-16.D0*SIG30_BB/CUTBB**30+4.0D0*SIG6_BB/CUTBB**6) 
103:       RCONSTAA = 0.0D0 
104:       RCONSTAB = A_AB*B_AB/(2.0D0*CUTAB)*EXP(-B_AB*CUTAB) 
105:      &          -3.0D0*C_AB/CUTAB**8 
106:      &          +4.0D0*EPS_AB*(15.0D0*SIG30_AB/CUTAB**32-3.0D0*SIG6_AB/CUTAB**8) 
107:       RCONSTBB = A_BB*B_BB/(2.0D0*CUTBB)*EXP(-B_BB*CUTBB) 
108:      &          -3.0D0*C_BB/CUTBB**8 
109:      &          +4.0D0*EPS_BB*(15.0D0*SIG30_BB/CUTBB**32-3.0D0*SIG6_BB/CUTBB**8) 
110:  
111:  
112: C  Work out cutoff for potential. Two particles interact if r<c, but 
113: C  we will use the equivalent condition 1/r^2 > 1/c^2. 
114: C 
115: C  Deal with any atoms that have left the box. 
116: C 
117:       IF ((.NOT.FIXIMAGE).AND.(.NOT.NORESET)) THEN 
118:          DO J1=1,N 
119:             J2 = 3*(J1-1) 
120: C           IF (ANINT(X(J2+1)/BOXLX).NE.0) PRINT*,'resetting X for J1,X,Xnew=',J1,X(J2+1),X(J2+1) - BOXLX*ANINT(X(J2+1)/BOXLX) 
121: C           IF (ANINT(X(J2+2)/BOXLX).NE.0) PRINT*,'resetting Y for J1,X,Xnew=',J1,X(J2+2),X(J2+2) - BOXLY*ANINT(X(J2+2)/BOXLY) 
122: C           IF (ANINT(X(J2+3)/BOXLX).NE.0) PRINT*,'resetting Z for J1,X,Xnew=',J1,X(J2+3),X(J2+3) - BOXLZ*ANINT(X(J2+3)/BOXLZ) 
123:             X(J2+1)=X(J2+1) - BOXLX*ANINT(X(J2+1)/BOXLX) 
124:             X(J2+2)=X(J2+2) - BOXLY*ANINT(X(J2+2)/BOXLY) 
125:             X(J2+3)=X(J2+3) - BOXLZ*ANINT(X(J2+3)/BOXLZ) 
126:          ENDDO 
127:       ENDIF 
128: C 
129: C  Calculate interatomic vectors using the minimum image convention. 
130: C  VEC(i,j,alpha) is the alpha (x, y or z) component of the vector between 
131: C  atoms i and j. 
132: C 
133:  
134:       POTEL=0.0D0 
135:  
136:       ! sn402: We will profile this at some point and decide whether it needs to be rewritten 
137:       ! in a more efficient manner. 
138:       IF (STEST) THEN 
139:          ! calculate first and second derivatives of potential (hessian) 
140:          V(1:3*N)=0.0D0 ! not necessary if .not.GTEST .and. .not.STEST 
141:          HESS(:,:)=0.0D0  
142:  
143:          DO J1=1,NTYPEA 
144:             DO J2=J1+1,NTYPEA 
145:                CALL SIO2PSHIFT_UPDATE_PAIRS(N, X, QQ_AA, A_AA, B_AA, C_AA, 
146:      &           EPS_AA, SIG_AA, RCONSTAA, CONSTAA, CUTAA, IRCUT2AA, 
147:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
148:             ENDDO 
149:          ENDDO 
150:          DO J1=1,NTYPEA 
151:             DO J2=NTYPEA+1,N 
152:                CALL SIO2PSHIFT_UPDATE_PAIRS(N, X, QQ_AB, A_AB, B_AB, C_AB, 
153:      &           EPS_AB, SIG_AB, RCONSTAB, CONSTAB, CUTAB, IRCUT2AB, 
154:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
155:             ENDDO 
156:          ENDDO 
157:          DO J1=NTYPEA+1,N 
158:             DO J2=J1+1,N 
159:               CALL SIO2PSHIFT_UPDATE_PAIRS(N, X, QQ_BB, A_BB, B_BB, C_BB, 
160:      &           EPS_BB, SIG_BB, RCONSTBB, CONSTBB, CUTBB, IRCUT2BB, 
161:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
162:             ENDDO 
163:          ENDDO 
164:  
165:        ELSEIF (GTEST) THEN  
166:          ! calculate first derivatives 
167:          V(1:3*N)=0.0D0 ! not necessary if .not.GTEST .and. .not.STEST 
168:  
169:          DO J1=1,NTYPEA 
170:             DO J2=J1+1,NTYPEA 
171:                CALL SIO2PSHIFT_UPDATE_PAIRG(N, X, QQ_AA, A_AA, B_AA, C_AA, 
172:      &           EPS_AA, SIG_AA, RCONSTAA, CONSTAA, CUTAA, IRCUT2AA, 
173:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
174:             ENDDO 
175:          ENDDO 
176:          DO J1=1,NTYPEA 
177:             DO J2=NTYPEA+1,N 
178:                CALL SIO2PSHIFT_UPDATE_PAIRG(N, X, QQ_AB, A_AB, B_AB, C_AB, 
179:      &           EPS_AB, SIG_AB, RCONSTAB, CONSTAB, CUTAB, IRCUT2AB, 
180:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
181:             ENDDO 
182:          ENDDO 
183:          DO J1=NTYPEA+1,N 
184:             DO J2=J1+1,N 
185:                CALL SIO2PSHIFT_UPDATE_PAIRG(N, X, QQ_BB, A_BB, B_BB, C_BB, 
186:      &           EPS_BB, SIG_BB, RCONSTBB, CONSTBB, CUTBB, IRCUT2BB, 
187:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
188:             ENDDO 
189:          ENDDO 
190:  
191:        ELSE  
192:          ! calculate potential only 
193:  
194:          DO J1=1,NTYPEA 
195:             DO J2=J1+1,NTYPEA 
196:                ! A-A interaction 
197:                CALL SIO2PSHIFT_UPDATE_PAIR(N, X, QQ_AA, A_AA, B_AA, C_AA, 
198:      &           EPS_AA, SIG_AA, RCONSTAA, CONSTAA, CUTAA, IRCUT2AA, 
199:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE) 
200:             ENDDO 
201:          ENDDO 
202:          DO J1=1,NTYPEA 
203:             DO J2=NTYPEA+1,N 
204:                ! A-B interaction 
205:                CALL SIO2PSHIFT_UPDATE_PAIR(N, X, QQ_AB, A_AB, B_AB, C_AB, 
206:      &           EPS_AB, SIG_AB, RCONSTAB, CONSTAB, CUTAB, IRCUT2AB, 
207:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE) 
208:             ENDDO 
209:          ENDDO 
210:          DO J1=NTYPEA+1,N 
211:             DO J2=J1+1,N 
212:                ! B-B interaction 
213:                CALL SIO2PSHIFT_UPDATE_PAIR(N, X, QQ_BB, A_BB, B_BB, C_BB, 
214:      &           EPS_BB, SIG_BB, RCONSTBB, CONSTBB, CUTBB, IRCUT2BB, 
215:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE) 
216:             ENDDO 
217:          ENDDO 
218:  
219:       ENDIF 
220:  
221: C      PRINT *,'mb2098: POTEL=',POTEL 
222:       RETURN 
223:       END 
224:  
225:  
226:  
227: C***************************************************************************** 
228:  
229:       SUBROUTINE SIO2PSHIFT_UPDATE_PAIRS(N, X, QQ, A, B, C, 
230:      &           EPS, SIG, RCONST, CONST, RCUT, IRCUT2, 
231:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
232:         !Calculate the potential energy and gradient and hessian between atoms 
233:         !j1 and j2. 
234:         !Update also ANV 
235:         USE MODHESS 
236:  
237:         IMPLICIT NONE 
238:  
239:         INTEGER, INTENT(IN) :: J1, J2, N 
240:         DOUBLE PRECISION, INTENT(IN) :: QQ, A, B, C, RCONST, CONST, RCUT, IRCUT2 
241:         DOUBLE PRECISION, INTENT(IN) :: EPS, SIG 
242:         DOUBLE PRECISION, INTENT(IN) :: X(3*N) 
243:         DOUBLE PRECISION, INTENT(IN) :: BOXLX, BOXLY, BOXLZ 
244:         LOGICAL, INTENT(IN) :: FIXIMAGE 
245:         DOUBLE PRECISION, INTENT(OUT) :: POTEL 
246:         INTEGER, INTENT(INOUT) :: ANV(N,N,3) 
247:         DOUBLE PRECISION, INTENT(INOUT) :: V(3*N) 
248:  
249:         DOUBLE PRECISION :: R, R2, R6, VEC(3), G, R8, R30, R32, FR2, TEMP, COULPOT 
250:         DOUBLE PRECISION :: SIG6, SIG30 
251:         INTEGER :: J3, J4, J5, J6 
252:  
253:         SIG6 = SIG**6 
254:         SIG30 = SIG**30 
255:  
256:         J3=3*(J1-1) 
257:         J4=3*(J2-1) 
258:         !find particle separation with suitable boundary conditions 
259:         VEC(1)=X(J3+1)-X(J4+1) 
260:         VEC(2)=X(J3+2)-X(J4+2) 
261:         VEC(3)=X(J3+3)-X(J4+3) 
262:         IF (.NOT.FIXIMAGE) THEN 
263:           ANV(J2,J1,1)=NINT(VEC(1)/BOXLX) 
264:           ANV(J2,J1,2)=NINT(VEC(2)/BOXLY) 
265:           ANV(J2,J1,3)=NINT(VEC(3)/BOXLZ) 
266:           ANV(J1,J2,1)=-ANV(J2,J1,1) 
267:           ANV(J1,J2,2)=-ANV(J2,J1,2) 
268:           ANV(J1,J2,3)=-ANV(J2,J1,3) 
269:         ENDIF 
270:         VEC(1)=VEC(1)-BOXLX*ANV(J2,J1,1) 
271:         VEC(2)=VEC(2)-BOXLY*ANV(J2,J1,2) 
272:         VEC(3)=VEC(3)-BOXLZ*ANV(J2,J1,3) 
273:  
274:  
275:         R2=VEC(1)**2+VEC(2)**2+VEC(3)**2 
276:         R = SQRT(R2) 
277:         R2=1.0D0/R2 
278:         IF (R2.GT.IRCUT2) THEN 
279:           !update potential energy 
280:           R6=R2**3 
281:           R30=R2**15 
282:           COULPOT = QQ*(1.0D0/R - 1.0D0/RCUT + IRCUT2*(R-RCUT)) 
283:           POTEL = POTEL + COULPOT+A*EXP(-B*R)-C*R6+RCONST/R2+CONST 
284:      &            +4.0D0*EPS*(SIG30*R30-SIG6*R6)  
285:           !update derivative of potential energy 
286:           R8=R6*R2 
287:           R32=R8**4 
288:           G = QQ*(-1.0D0/R**3 + 1.0D0/R*IRCUT2) - A*B/R*EXP(-B*R) + 6.0D0*C*R8 + 2.0D0*RCONST 
289:      &        +4.0D0*EPS*(6.0D0*SIG6*R8-30.0D0*SIG30*R32) 
290:           DO J5=1,3 
291:             !careful with signs: both VEC and V are antisymmetric in J1<->J2 
292:             V(J3+J5)=V(J3+J5)+G*VEC(J5) 
293:             V(J4+J5)=V(J4+J5)-G*VEC(J5) 
294:           END DO 
295:  
296:           !calculate the second derivative, and update the components of the 
297:           !hessian 
298:           FR2 = (QQ*(3.0D0/R**3 - 1.0D0/R*IRCUT2) +A*B*EXP(-B*R)*(1.0D0/R+B)-48.0D0*C*R8 
299:      &          +4.0D0*EPS*(960.0D0*SIG30*R32-48.0D0*SIG6*R8))*R2 
300: C  js850> same cartesian coordinate 
301:           DO J5=1,3 
302:             TEMP = FR2 * VEC(J5)**2 + G 
303: C  Now do the hessian. First are the entirely diagonal terms. 
304:             HESS(J3+J5,J3+J5) = HESS(J3+J5,J3+J5) + TEMP 
305:             HESS(J4+J5,J4+J5) = HESS(J4+J5,J4+J5) + TEMP 
306: C  Case III, different atoms, same cartesian coordinate. 
307:             HESS(J4+J5,J3+J5) = - TEMP 
308:             HESS(J3+J5,J4+J5) = - TEMP 
309:           END DO 
310: C  js850> different cartesian coordinates 
311:           DO J5=1,3 
312:             DO J6=J5+1,3 
313:               TEMP = FR2 * VEC(J5) * VEC(J6)  
314: C  Next are the terms where x_i and x_j are on the same atom 
315: C  but are different, e.g. y and z. 
316:               HESS(J3+J5,J3+J6) = HESS(J3+J5,J3+J6) + TEMP 
317:               HESS(J4+J5,J4+J6) = HESS(J4+J5,J4+J6) + TEMP 
318:               ! J5 <-> J6 
319:               HESS(J3+J6,J3+J5) = HESS(J3+J6,J3+J5) + TEMP 
320:               HESS(J4+J6,J4+J5) = HESS(J4+J6,J4+J5) + TEMP 
321: C  Case IV: different atoms and different cartesian coordinates. 
322:               HESS(J4+J5,J3+J6) = -TEMP 
323:               HESS(J3+J5,J4+J6) = -TEMP 
324:               ! J5 <-> J6 
325:               HESS(J4+J6,J3+J5) = -TEMP 
326:               HESS(J3+J6,J4+J5) = -TEMP 
327:             END DO 
328:           END DO 
329:  
330:         ELSE 
331:           !set hessian to zero 
332:  
333:           DO J5=1,3 
334:             DO J6=J5,3 
335:               HESS(J4+J5,J3+J6) = 0.D0 
336:               HESS(J3+J5,J4+J6) = 0.D0 
337:               ! J5 <-> J6 
338:               HESS(J4+J6,J3+J5) = 0.D0 
339:               HESS(J3+J6,J4+J5) = 0.D0 
340:             END DO 
341:           END DO 
342:  
343:         ENDIF 
344:  
345:       END SUBROUTINE SIO2PSHIFT_UPDATE_PAIRS 
346:  
347:  
348:  
349: C***************************************************************************** 
350:  
351:       SUBROUTINE SIO2PSHIFT_UPDATE_PAIRG(N, X, QQ, A, B, C,  
352:      &           EPS, SIG, RCONST, CONST, RCUT, IRCUT2, 
353:      &           POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE, V) 
354:         !Calculate the potential energy and gradient between atoms j1 and j2. 
355:         !Update also ANV 
356:  
357:         IMPLICIT NONE 
358:  
359:         INTEGER, INTENT(IN) :: J1, J2, N 
360:         DOUBLE PRECISION, INTENT(IN) :: QQ, A, B, C, RCONST, CONST, RCUT, IRCUT2 
361:         DOUBLE PRECISION, INTENT(IN) :: EPS, SIG 
362:         DOUBLE PRECISION, INTENT(IN) :: X(3*N) 
363:         DOUBLE PRECISION, INTENT(IN) :: BOXLX, BOXLY, BOXLZ 
364:         LOGICAL, INTENT(IN) :: FIXIMAGE 
365:         DOUBLE PRECISION, INTENT(OUT) :: POTEL 
366:         INTEGER, INTENT(INOUT) :: ANV(N,N,3) 
367:         DOUBLE PRECISION, INTENT(INOUT) :: V(3*N) 
368:  
369:         DOUBLE PRECISION :: R, R2, R6, VEC(3), G, R8, R30, R32, COULPOT 
370:         DOUBLE PRECISION :: SIG6, SIG30 
371:         INTEGER :: J3, J4, J5 
372:  
373:         SIG6 = SIG**6 
374:         SIG30 = SIG**30 
375:  
376:         J3=3*(J1-1) 
377:         J4=3*(J2-1) 
378:         !find particle separation with suitable boundary conditions 
379:         VEC(1)=X(J3+1)-X(J4+1) 
380:         VEC(2)=X(J3+2)-X(J4+2) 
381:         VEC(3)=X(J3+3)-X(J4+3) 
382:         IF (.NOT.FIXIMAGE) THEN 
383:           ANV(J2,J1,1)=NINT(VEC(1)/BOXLX) 
384:           ANV(J2,J1,2)=NINT(VEC(2)/BOXLY) 
385:           ANV(J2,J1,3)=NINT(VEC(3)/BOXLZ) 
386:           ANV(J1,J2,1)=-ANV(J2,J1,1)  
387:           ANV(J1,J2,2)=-ANV(J2,J1,2) 
388:           ANV(J1,J2,3)=-ANV(J2,J1,3) 
389:         ENDIF 
390:         VEC(1)=VEC(1)-BOXLX*ANV(J2,J1,1) 
391:         VEC(2)=VEC(2)-BOXLY*ANV(J2,J1,2) 
392:         VEC(3)=VEC(3)-BOXLZ*ANV(J2,J1,3) 
393:  
394:         !calculate the potential 
395:         R2=VEC(1)**2+VEC(2)**2+VEC(3)**2 
396:         R=SQRT(R2) 
397:         R2=1.0D0/R2 
398:         IF (R2.GT.IRCUT2) THEN 
399:           R30=R2**15 
400:           R6=R2**3 
401:           COULPOT = QQ*(1.0D0/R - 1.0D0/RCUT + IRCUT2*(R-RCUT)) 
402:           POTEL = POTEL + COULPOT+A*EXP(-B*R)-C*R6+RCONST/R2+CONST 
403:      &            +4.0D0*EPS*(SIG30*R30-SIG6*R6) 
404:           !update derivative of potential energy 
405:           R8=R6*R2 
406:           R32=R8**4 
407:           G = QQ*(-1.0D0/R**3 + 1.0D0/R*IRCUT2) - A*B/R*EXP(-B*R) + 6.0D0*C*R8 + 2.0D0*RCONST 
408:      &        +4.0D0*EPS*(6.0D0*SIG6*R8-30.0D0*SIG30*R32) 
409:           DO J5=1,3 
410:             V(J3+J5)=V(J3+J5)+G*VEC(J5) 
411:             V(J4+J5)=V(J4+J5)-G*VEC(J5) 
412:           END DO   
413:         ENDIF 
414:  
415:       END SUBROUTINE SIO2PSHIFT_UPDATE_PAIRG 
416:  
417:  
418: C***************************************************************************** 
419:  
420:       SUBROUTINE SIO2PSHIFT_UPDATE_PAIR(N, X, QQ, A, B, C, EPS, SIG, RCONST,  
421:      &           CONST, RCUT, IRCUT2, POTEL, J1, J2, BOXLX, BOXLY, BOXLZ, ANV, FIXIMAGE) 
422:         !calculate the potential energy between atoms j1 and j2  
423:         !update ANV 
424:  
425:         IMPLICIT NONE 
426:  
427:         INTEGER, INTENT(IN) :: J1, J2, N 
428:         DOUBLE PRECISION, INTENT(IN) :: QQ, A, B, C, RCONST, CONST, RCUT, IRCUT2 
429:         DOUBLE PRECISION, INTENT(IN) :: EPS, SIG 
430:         DOUBLE PRECISION, INTENT(IN) :: X(3*N) 
431:         DOUBLE PRECISION, INTENT(IN) :: BOXLX, BOXLY, BOXLZ 
432:         LOGICAL, INTENT(IN) :: FIXIMAGE 
433:         DOUBLE PRECISION, INTENT(OUT) :: POTEL 
434:         INTEGER, INTENT(INOUT) :: ANV(N,N,3) 
435:  
436:         DOUBLE PRECISION :: R, R2, R6, R30, VEC1, VEC2, VEC3, COULPOT 
437:         DOUBLE PRECISION :: SIG6, SIG30 
438:         INTEGER :: J3, J4 
439:  
440:         SIG6 = SIG**6 
441:         SIG30 = SIG**30 
442:  
443:         J3=3*(J1-1) 
444:         J4=3*(J2-1) 
445:         !find particle separation with suitable boundary conditions 
446:         VEC1=X(J3+1)-X(J4+1) 
447:         VEC2=X(J3+2)-X(J4+2) 
448:         VEC3=X(J3+3)-X(J4+3) 
449:         IF (.NOT.FIXIMAGE) THEN 
450:           ANV(J2,J1,1)=NINT(VEC1/BOXLX) 
451:           ANV(J2,J1,2)=NINT(VEC2/BOXLY) 
452:           ANV(J2,J1,3)=NINT(VEC3/BOXLZ) 
453:           ANV(J1,J2,1)=-ANV(J2,J1,1)  
454:           ANV(J1,J2,2)=-ANV(J2,J1,2) 
455:           ANV(J1,J2,3)=-ANV(J2,J1,3) 
456:         ENDIF 
457:         VEC1=VEC1-BOXLX*ANV(J2,J1,1) 
458:         VEC2=VEC2-BOXLY*ANV(J2,J1,2) 
459:         VEC3=VEC3-BOXLZ*ANV(J2,J1,3) 
460:         !calculate the potential 
461:         R2=VEC1**2+VEC2**2+VEC3**2 
462:  
463:         R=SQRT(R2) 
464:         R2=1.0D0/R2 
465:         R30=R2**15 
466:  
467:         IF (R2.GT.IRCUT2) THEN 
468:           R6=R2**3 
469:           COULPOT = QQ*(1.0D0/R - 1.0D0/RCUT + IRCUT2*(R-RCUT)) 
470:           POTEL = POTEL + COULPOT+A*EXP(-B*R)-C*R6+RCONST/R2+CONST 
471:      &            +4.0D0*EPS*(SIG30*R30-SIG6*R6) 
472:         ENDIF 
473:  
474:       END SUBROUTINE SIO2PSHIFT_UPDATE_PAIR 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0