hdiff output

r32888/commons.f90 2017-06-30 16:30:11.218631451 +0100 r32887/commons.f90 2017-06-30 16:30:13.158657344 +0100
114:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &114:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
115:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &115:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
116:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &116:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
117:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, LJ_VAR_RADT, READMASST, SPECMASST, NEWTSALLIST, &117:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, LJ_VAR_RADT, READMASST, SPECMASST, NEWTSALLIST, &
118:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, &118:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, &
119:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &119:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
120:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &120:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
121:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &121:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
122:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &122:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &
123:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS, PYADDT, PYADD2T, LJADD3T, REORDERADDT,  LJADD4T, &123:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS, PYADDT, PYADD2T, LJADD3T, REORDERADDT,  LJADD4T, &
124:      &        SQNMT, SQNM_DEBUGT, SQNM_BIOT, BENZRIGIDEWALDT, ORTHO, EWALDT, WATERMETHANET124:      &        SQNMT, SQNM_DEBUGT, SQNM_BIOT, BENZRIGIDEWALDT, ORTHO, EWALDT
125: !125: !
126:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:)126:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:)
127:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)127:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)
128:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:)128:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:)
129: 129: 
130: ! dj337: Ewald summation variables130: ! dj337: Ewald summation variables
131:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: RERHOARRAY, IMRHOARRAY131:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: RERHOARRAY, IMRHOARRAY
132: 132: 
133: ! csw34> FREEZEGROUP variables133: ! csw34> FREEZEGROUP variables
134: !134: !


r32888/genrigid.f90 2017-06-30 16:30:11.710638019 +0100 r32887/genrigid.f90 2017-06-30 16:30:13.438661080 +0100
379:   CLOSE(222)379:   CLOSE(222)
380: 380: 
381: !  vr274> Calling function for allocation to make more general (setup rigid bodies from code)381: !  vr274> Calling function for allocation to make more general (setup rigid bodies from code)
382:   CALL GENRIGID_ALLOCATE(NRIGIDBODY,MAXSITE)382:   CALL GENRIGID_ALLOCATE(NRIGIDBODY,MAXSITE)
383: 383: 
384: ! hk286 > initialise SITESRIGIDBODY, RIGIDGROUPS, RIGIDSINGLES 384: ! hk286 > initialise SITESRIGIDBODY, RIGIDGROUPS, RIGIDSINGLES 
385:   OPEN(UNIT=222,FILE='rbodyconfig',status='unknown')385:   OPEN(UNIT=222,FILE='rbodyconfig',status='unknown')
386:   DO J1 = 1, NRIGIDBODY386:   DO J1 = 1, NRIGIDBODY
387:      READ(222,*) CHECK1, NSITEPERBODY(J1)387:      READ(222,*) CHECK1, NSITEPERBODY(J1)
388:      DO J2 = 1, NSITEPERBODY(J1)388:      DO J2 = 1, NSITEPERBODY(J1)
389:         write(*,*) J1, J2, NSITEPERBODY(J1), CHECK1 
390:         READ(222,*) RIGIDGROUPS(J2,J1)389:         READ(222,*) RIGIDGROUPS(J2,J1)
391: ! csw34> check to make sure the current atom is not already in a rigid body390: ! csw34> check to make sure the current atom is not already in a rigid body
392:         IF (RIGIDISRIGID(RIGIDGROUPS(J2,J1))) THEN391:         IF (RIGIDISRIGID(RIGIDGROUPS(J2,J1))) THEN
393:             PRINT *," genrigid> ERROR: atom ",RIGIDGROUPS(J2,J1)," is in multiple rigid bodies! Stopping."392:             PRINT *," genrigid> ERROR: atom ",RIGIDGROUPS(J2,J1)," is in multiple rigid bodies! Stopping."
394:             STOP393:             STOP
395:         ELSE       394:         ELSE       
396: ! csw34> if not, flag the current atom395: ! csw34> if not, flag the current atom
397:             RIGIDISRIGID(RIGIDGROUPS(J2,J1))=.TRUE.396:             RIGIDISRIGID(RIGIDGROUPS(J2,J1))=.TRUE.
398:             RB_BY_ATOM(RIGIDGROUPS(J2,J1)) = J1397:             RB_BY_ATOM(RIGIDGROUPS(J2,J1)) = J1
399:         ENDIF398:         ENDIF
400: ! vr274> Moved initialization of coordinates to GENRIGID_INITIALISE, here only read the setup399: ! vr274> Moved initialization of coordinates to GENRIGID_INITIALISE, here only read the setup
401: !        SITESRIGIDBODY(J2,:,J1) = COORDS(3*DUMMY-2:3*DUMMY,1)400: !        SITESRIGIDBODY(J2,:,J1) = COORDS(3*DUMMY-2:3*DUMMY,1)
402:      ENDDO401:      ENDDO
403:   ENDDO402:   ENDDO
404:   CLOSE(222)403:   CLOSE(222)
 404: 
405:   CALL GENRIGID_INITIALISE(INICOORDS)405:   CALL GENRIGID_INITIALISE(INICOORDS)
406: END SUBROUTINE GENRIGID_READ_FROM_FILE406: END SUBROUTINE GENRIGID_READ_FROM_FILE
407: 407: 
408: !-----------------------------------------------------------408: !-----------------------------------------------------------
409: 409: 
410: !-----------------------------------------------------------410: !-----------------------------------------------------------
411: 411: 
412: SUBROUTINE TRANSFORMRIGIDTOC (CMIN, CMAX, XCOORDS, XRIGIDCOORDS)412: SUBROUTINE TRANSFORMRIGIDTOC (CMIN, CMAX, XCOORDS, XRIGIDCOORDS)
413:       413:       
414:   USE COMMONS, ONLY: NATOMS414:   USE COMMONS, ONLY: NATOMS


r32888/keywords.f 2017-06-30 16:30:12.030642289 +0100 r32887/keywords.f 2017-06-30 16:30:13.730664977 +0100
1117:       POLIRT=.FALSE.1117:       POLIRT=.FALSE.
1118: ! jdf43> MBPOL1118: ! jdf43> MBPOL
1119:       MBPOLT=.FALSE.1119:       MBPOLT=.FALSE.
1120:       MOLECULART=.FALSE.1120:       MOLECULART=.FALSE.
1121: ! jdf43>1121: ! jdf43>
1122:       REPMATCHT=.FALSE.1122:       REPMATCHT=.FALSE.
1123: 1123: 
1124:       UNIFORMMOVE=.FALSE.1124:       UNIFORMMOVE=.FALSE.
1125:       ORBITTOL=1.0D-31125:       ORBITTOL=1.0D-3
1126:       NOINVERSION=.FALSE.1126:       NOINVERSION=.FALSE.
1127: !cv320> rigid body watermethane1127: !
1128:       WATERMETHANET=.FALSE. 
1129: ! ds656> Parallelised generalised basin-hopping1128: ! ds656> Parallelised generalised basin-hopping
1130:       GBHT=.FALSE.1129:       GBHT=.FALSE.
1131: 1130: 
1132: ! General mixed LJ systems1131: ! General mixed LJ systems
1133:       GLJT=.FALSE.1132:       GLJT=.FALSE.
1134: ! ds656> Multicomponent LJ system (different implementation to GLJ!)1133: ! ds656> Multicomponent LJ system (different implementation to GLJ!)
1135:       MLJT=.FALSE.1134:       MLJT=.FALSE.
1136: ! jwrm2> LJ with each atom having a different radius. Works with RIGIDINIT.1135: ! jwrm2> LJ with each atom having a different radius. Works with RIGIDINIT.
1137:       LJ_VAR_RADT = .FALSE.1136:       LJ_VAR_RADT = .FALSE.
1138: 1137: 
3672:          ALLOCATE(STCHRG(NRBSITES))3671:          ALLOCATE(STCHRG(NRBSITES))
3673: 3672: 
3674:          CALL DEFPAHARIGID()3673:          CALL DEFPAHARIGID()
3675:          NCARBON = 63674:          NCARBON = 6
3676:          CALL DEFBENZENERIGID()3675:          CALL DEFBENZENERIGID()
3677: 3676: 
3678: ! ----------------------------------------------------------------------------------------3677: ! ----------------------------------------------------------------------------------------
3679: 3678: 
3680:       ELSE IF (WORD.EQ.'FAL') THEN3679:       ELSE IF (WORD.EQ.'FAL') THEN
3681:          FAL=.TRUE.3680:          FAL=.TRUE.
3682:       ELSE IF (WORD.EQ.'WATERMETHANE') THEN3681: 
3683:          WATERMETHANET=.TRUE. 
3684:       ELSE IF (WORD == 'FEBH') THEN3682:       ELSE IF (WORD == 'FEBH') THEN
3685:          CALL READF(FETEMP)3683:          CALL READF(FETEMP)
3686:          FEBHT = .TRUE.3684:          FEBHT = .TRUE.
3687:          FE_FILE_UNIT = GETUNIT()3685:          FE_FILE_UNIT = GETUNIT()
3688:          OPEN(UNIT = FE_FILE_UNIT, FILE = 'free_energy', STATUS = 'REPLACE')3686:          OPEN(UNIT = FE_FILE_UNIT, FILE = 'free_energy', STATUS = 'REPLACE')
3689:          WRITE(FE_FILE_UNIT, '(6A20)') '       Quench       ', '  Potential energy  ',3687:          WRITE(FE_FILE_UNIT, '(6A20)') '       Quench       ', '  Potential energy  ',
3690:      &   '   Harmonic term    ', '    Free energy     ', '   Markov energy    ', '        Time        '3688:      &   '   Harmonic term    ', '    Free energy     ', '   Markov energy    ', '        Time        '
3691:          WRITE(FE_FILE_UNIT, '(6A20)') ' ------------------ ', ' ------------------ ',3689:          WRITE(FE_FILE_UNIT, '(6A20)') ' ------------------ ', ' ------------------ ',
3692:      &   ' ------------------ ', ' ------------------ ', ' ------------------ ', ' ------------------ '3690:      &   ' ------------------ ', ' ------------------ ', ' ------------------ ', ' ------------------ '
3693:          IF (NITEMS .GT. 2) THEN3691:          IF (NITEMS .GT. 2) THEN


r32888/minpermdist.f90 2017-06-30 16:30:12.318646133 +0100 r32887/minpermdist.f90 2017-06-30 16:30:14.002668607 +0100
 77: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES) 77: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES)
 78: SAVE NORBIT1, NORBIT2 78: SAVE NORBIT1, NORBIT2
 79: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS) 79: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS)
 80: CHARACTER(LEN=5) ZSYMSAVE 80: CHARACTER(LEN=5) ZSYMSAVE
 81: COMMON /SYS/ ZSYMSAVE 81: COMMON /SYS/ ZSYMSAVE
 82:  82: 
 83: DOUBLE PRECISION :: DINV 83: DOUBLE PRECISION :: DINV
 84:  84: 
 85: DISTANCE = 0.0D0 85: DISTANCE = 0.0D0
 86: RMAT=0.0D0 86: RMAT=0.0D0
 87: !cv320: commented out this write statement, looks like it was just for debugging. 87: WRITE(*,*) 'RMAT', RMAT
 88: ! WRITE(*,*) 'RMAT', RMAT 
 89: ! sf344 temporary debug 88: ! sf344 temporary debug
 90: !      CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGID,DEBUG,RMAT) 89: !      CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGID,DEBUG,RMAT)
 91: !      RMATBEST = RMAT 90: !      RMATBEST = RMAT
 92: !   RETURN 91: !   RETURN
 93: ! end sf344 92: ! end sf344
 94:  93: 
 95: ! jwrm2 94: ! jwrm2
 96: IF (GTHOMSONT) THEN 95: IF (GTHOMSONT) THEN
 97:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,DISTANCE,DIST2,RMATBEST) 96:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,DISTANCE,DIST2,RMATBEST)
 98:    RETURN 97:    RETURN


r32888/potential.f90 2017-06-30 16:30:12.618650136 +0100 r32887/potential.f90 2017-06-30 16:30:14.298672558 +0100
 42:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS 42:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS
 43:    USE GAUSS_MOD, ONLY: GFIELD 43:    USE GAUSS_MOD, ONLY: GFIELD
 44:    USE RAD_MOD, ONLY: RAD, RADR 44:    USE RAD_MOD, ONLY: RAD, RADR
 45:    USE PREC, ONLY: INT32, REAL64 45:    USE PREC, ONLY: INT32, REAL64
 46:    USE TWIST_MOD, ONLY: TWISTT, TWIST 46:    USE TWIST_MOD, ONLY: TWISTT, TWIST
 47:    USE CONVEX_POLYHEDRA_MODULE, ONLY: CONVEX_POLYHEDRA 47:    USE CONVEX_POLYHEDRA_MODULE, ONLY: CONVEX_POLYHEDRA
 48:    USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS 48:    USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS
 49:    USE OPP_MOD, ONLY: OPP 49:    USE OPP_MOD, ONLY: OPP
 50:    USE LJ_VAR_RAD_MOD, ONLY: LJ_VAR_RAD 50:    USE LJ_VAR_RAD_MOD, ONLY: LJ_VAR_RAD
 51:    USE EWALD 51:    USE EWALD
 52:    USE watermethane_mod 
 53:    IMPLICIT NONE 52:    IMPLICIT NONE
 54: ! Arguments 53: ! Arguments
 55: ! TODO: Work out intents 54: ! TODO: Work out intents
 56: ! TODO: Fix array dimensions? 55: ! TODO: Fix array dimensions?
 57:    REAL(REAL64)               :: X(*) 56:    REAL(REAL64)               :: X(*)
 58:    REAL(REAL64)               :: GRAD(*) 57:    REAL(REAL64)               :: GRAD(*)
 59:    REAL(REAL64)               :: EREAL 58:    REAL(REAL64)               :: EREAL
 60:    LOGICAL, INTENT(IN)        :: GRADT 59:    LOGICAL, INTENT(IN)        :: GRADT
 61:    LOGICAL, INTENT(IN)        :: SECT 60:    LOGICAL, INTENT(IN)        :: SECT
 62: ! Variables 61: ! Variables
1282: 1281: 
1283:    ELSE IF (MBPOLT) THEN1282:    ELSE IF (MBPOLT) THEN
1284:       CALL MBPOL (X, EREAL, GRAD, GRADT)1283:       CALL MBPOL (X, EREAL, GRAD, GRADT)
1285: 1284: 
1286:    ELSE IF (WATERDCT) THEN1285:    ELSE IF (WATERDCT) THEN
1287:       CALL WATERPDC (X, GRAD, EREAL, GRADT)1286:       CALL WATERPDC (X, GRAD, EREAL, GRADT)
1288: 1287: 
1289:    ELSE IF (WATERKZT) THEN1288:    ELSE IF (WATERKZT) THEN
1290:       CALL WATERPKZ (X, GRAD, EREAL, GRADT)1289:       CALL WATERPKZ (X, GRAD, EREAL, GRADT)
1291: 1290: 
1292:    ELSE IF (WATERMETHANET) THEN 
1293:       IF (.NOT. RIGIDINIT) THEN 
1294:          WRITE(*,*) "WATER-METHANE POTENTIAL ONLY IMPLEMENTED FOR GENRIGID FRAMEWORK" 
1295:          STOP 
1296:       END IF 
1297:       IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN 
1298:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS) 
1299:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS) 
1300:       END IF 
1301:       CALL wmrb(X, GRAD, EREAL, GRADT) 
1302:       IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN 
1303:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
1304:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS) 
1305:          IF (GRADT) THEN 
1306:             GRADATOMS(1:3*NATOMS) = GRAD(1:3*NATOMS) ! Needed for AACONVERGENCE 
1307:             CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD) 
1308:             GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
1309:             GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS) 
1310:          END IF 
1311:      END IF 
1312:  
1313:  
1314:    ELSE IF (GAYBERNET) THEN1291:    ELSE IF (GAYBERNET) THEN
1315:       CALL GAYBERNE(X, GRAD, EREAL, GRADT,.FALSE.)1292:       CALL GAYBERNE(X, GRAD, EREAL, GRADT,.FALSE.)
1316: 1293: 
1317:    ELSE IF (PYGPERIODICT .OR. PYBINARYT) THEN1294:    ELSE IF (PYGPERIODICT .OR. PYBINARYT) THEN
1318:       IF (PYADDT) THEN1295:       IF (PYADDT) THEN
1319:          CALL PYGPERIODICADD(X, GRAD, EREAL, GRADT)1296:          CALL PYGPERIODICADD(X, GRAD, EREAL, GRADT)
1320:       ELSEIF (PYADD2T) THEN1297:       ELSEIF (PYADD2T) THEN
1321:          CALL PYGPERIODICADD2(X, GRAD, EREAL, GRADT)1298:          CALL PYGPERIODICADD2(X, GRAD, EREAL, GRADT)
1322:       ELSE1299:       ELSE
1323:          CALL PYGPERIODIC(X, GRAD, EREAL, GRADT)1300:          CALL PYGPERIODIC(X, GRAD, EREAL, GRADT)


r32888/watermethane.f90 2017-06-30 16:30:12.906653980 +0100 r32887/watermethane.f90 2017-06-30 16:30:14.590676455 +0100
  1: module watermethane_mod  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/watermethane.f90' in revision 32887
  2: !Module providing a potential for water-methane dimer, written by C L Vaillant. Water-Methane interaction 
  3: !is taken from J Chem Phys 123, 134311 (2005). 
  4: !Energies are in kcal/mol, distances are in Angstroms. 
  5:  
  6: implicit none 
  7:  
  8: integer, parameter::             waterdof=7, methanedof=9 
  9: !Water charges, in a.u., in the order: H1 H2 Q D1 D2 T1 T2 
 10: double precision, parameter::    watercharge(waterdof)= (/0.494714d0, 0.494714d0, -1.830627d0, 0.420599d0, & 
 11:      0.420599d0, 0.0d0, 0.0d0/) 
 12: !Methane charges, in a.u., in the order: H1, H2, H3, H4, C, M1, M2, M3, M4 
 13: double precision, parameter::    methanecharge(methanedof)= (/0.279901d0,0.279901d0,0.279901d0,0.279901d0, 3.590472d0, & 
 14:      -1.177519d0,-1.177519d0,-1.177519d0,-1.177519d0/) 
 15: !Parameters in Angstroms and kcal/mol 
 16: double precision::               betaang(waterdof, methanedof), Aang0(waterdof, methanedof), Aang1(waterdof, methanedof), AangM(waterdof, methanedof) 
 17: double precision::               Cang6(waterdof, methanedof), Cang8(waterdof, methanedof), Cang10(waterdof, methanedof), deltaang6(waterdof, methanedof) 
 18: double precision::               deltaang8(waterdof, methanedof) 
 19: !Parameters in atomic units 
 20: double precision::               beta(waterdof, methanedof), A0(waterdof, methanedof), A1(waterdof, methanedof), AM(waterdof, methanedof) 
 21: double precision::               C6(waterdof, methanedof), C8(waterdof, methanedof), C10(waterdof, methanedof), delta6(waterdof, methanedof) 
 22: double precision::               delta8(waterdof, methanedof) 
 23: !-------------------------------------------------------------------------------------------------------------------- 
 24: !Input of data 
 25: !H1 and H2: 
 26: data betaang(1,:)/2.84808454d0,2.84808454d0,2.84808454d0,2.84808454d0,2.7971225d0,2.75581866d0,2.75581866d0,2.75581866d0,2.75581866d0/ 
 27: data betaang(2,:)/2.84808454d0,2.84808454d0,2.84808454d0,2.84808454d0,2.7971225d0,2.75581866d0,2.75581866d0,2.75581866d0,2.75581866d0/ 
 28: !Q: 
 29: data betaang(3,:)/2.86928398d0,2.86928398d0,2.86928398d0,2.86928398d0,2.3463075d0,2.31474866d0,2.31474866d0,2.31474866d0,2.31474866d0/ 
 30: !D1 and D2: 
 31: data betaang(4,:)/5.71995231d0,5.71995231d0,5.71995231d0,5.71995231d0,3.16999754d0,2.35594058d0,2.35594058d0,2.35594058d0,2.35594058d0/ 
 32: data betaang(5,:)/5.71995231d0,5.71995231d0,5.71995231d0,5.71995231d0,3.16999754d0,2.35594058d0,2.35594058d0,2.35594058d0,2.35594058d0/ 
 33: !T1 and T2: 
 34: data betaang(6,:)/6.24776382d0,6.24776382d0,6.24776382d0,6.24776382d0,2.31915671d0,2.28762859d0,2.28762859d0,2.28762859d0,2.28762859d0/ 
 35: data betaang(7,:)/6.24776382d0,6.24776382d0,6.24776382d0,6.24776382d0,2.31915671d0,2.28762859d0,2.28762859d0,2.28762859d0,2.28762859d0/ 
 36:  
 37: !H1 and H2: 
 38: data Aang0(1,:)/-752.765963d0,-752.765963d0,-752.765963d0,-752.765963d0,-40504.8858d0,5933.09667d0,5933.09667d0,5933.09667d0,5933.09667d0/ 
 39: data Aang0(2,:)/-752.765963d0,-752.765963d0,-752.765963d0,-752.765963d0,-40504.8858d0,5933.09667d0,5933.09667d0,5933.09667d0,5933.09667d0/ 
 40: !Q: 
 41: data Aang0(3,:)/4592.62807d0,4592.62807d0,4592.62807d0,4592.62807d0,43408.9282d0,-5121.6292d0,-5121.6292d0,-5121.6292d0,-5121.6292d0/ 
 42: !D1 and D2: 
 43: data Aang0(4,:)/5367.76805d0,5367.76805d0,5367.76805d0,5367.76805d0,55943.4633d0,-2584.27027d0,-2584.27027d0,-2584.27027d0,-2584.27027d0/ 
 44: data Aang0(5,:)/5367.76805d0,5367.76805d0,5367.76805d0,5367.76805d0,55943.4633d0,-2584.27027d0,-2584.27027d0,-2584.27027d0,-2584.27027d0/ 
 45: !T1 and T2: 
 46: data Aang0(6,:)/1258.12101d0,1258.12101d0,1258.12101d0,1258.12101d0,-19777.5292d0,2979.79274d0,2979.79274d0,2979.79274d0,2979.79274d0/ 
 47: data Aang0(7,:)/1258.12101d0,1258.12101d0,1258.12101d0,1258.12101d0,-19777.5292d0,2979.79274d0,2979.79274d0,2979.79274d0,2979.79274d0/ 
 48:  
 49: !H1 and H2: 
 50: data AangM(1,:)/908.685355d0,908.685355d0,908.685355d0,908.685355d0,26577.2947d0,-2622.41721d0,-2622.41721d0,-2622.41721d0,-2622.41721d0/ 
 51: data AangM(2,:)/908.685355d0,908.685355d0,908.685355d0,908.685355d0,26577.2947d0,-2622.41721d0,-2622.41721d0,-2622.41721d0,-2622.41721d0/ 
 52: !Q: 
 53: data AangM(3,:)/1252.30889d0,1252.30889d0,1252.30889d0,1252.30889d0,-23339.3574d0,8389.34399d0,8389.34399d0,8389.34399d0,8389.34399d0/ 
 54: !D1 and D2: 
 55: data AangM(4,:)/-1430.20075d0,-1430.20075d0,-1430.20075d0,-1430.20075d0,-84638.8183d0,1242.92288d0,1242.92288d0,1242.92288d0,1242.92288d0/ 
 56: data AangM(5,:)/-1430.20075d0,-1430.20075d0,-1430.20075d0,-1430.20075d0,-84638.8183d0,1242.92288d0,1242.92288d0,1242.92288d0,1242.92288d0/ 
 57: !T1 and T2: 
 58: data AangM(6,:)/-162.796205d0,-162.796205d0,-162.796205d0,-162.796205d0,8863.16664d0,-668.609725d0,-668.609725d0,-668.609725d0,-668.609725d0/ 
 59: data AangM(7,:)/-162.796205d0,-162.796205d0,-162.796205d0,-162.796205d0,8863.16664d0,-668.609725d0,-668.609725d0,-668.609725d0,-668.609725d0/ 
 60:  
 61: !H1 and H2: 
 62: data Aang1(1,:)/417.797177d0,417.797177d0,417.797177d0,417.797177d0,19352.3942d0,-3719.99887d0,-3719.99887d0,-3719.99887d0,-3719.99887d0/ 
 63: data Aang1(2,:)/417.797177d0,417.797177d0,417.797177d0,417.797177d0,19352.3942d0,-3719.99887d0,-3719.99887d0,-3719.99887d0,-3719.99887d0/ 
 64: !Q: 
 65: data Aang1(3,:)/-1789.69987d0,-1789.69987d0,-1789.69987d0,-1789.69987d0,-29232.8793d0,5058.10255d0,5058.10255d0,5058.10255d0,5058.10255d0/ 
 66: !D1 and D2: 
 67: data Aang1(4,:)/-14542.3959d0,-14542.3959d0,-14542.3959d0,-14542.3959d0,-5391.49926d0,651.930524d0,651.930524d0,651.930524d0,651.930524d0/ 
 68: data Aang1(5,:)/-14542.3959d0,-14542.3959d0,-14542.3959d0,-14542.3959d0,-5391.49926d0,651.930524d0,651.930524d0,651.930524d0,651.930524d0/ 
 69: !T1 and T2: 
 70: data Aang1(6,:)/-9354.24387d0,-9354.24387d0,-9354.24387d0,-9354.24387d0,10463.6792d0,-2062.86173d0,-2062.86173d0,-2062.86173d0,-2062.86173d0/ 
 71: data Aang1(7,:)/-9354.24387d0,-9354.24387d0,-9354.24387d0,-9354.24387d0,10463.6792d0,-2062.86173d0,-2062.86173d0,-2062.86173d0,-2062.86173d0/ 
 72:  
 73: !H1 and H2: 
 74: data Cang6(1,:)/-31.1396325d0,-31.1396325d0,-31.1396325d0,-31.1396325d0,-176.385261d0,26.1819133d0,26.1819133d0,26.1819133d0,26.1819133d0/ 
 75: data Cang6(2,:)/-31.1396325d0,-31.1396325d0,-31.1396325d0,-31.1396325d0,-176.385261d0,26.1819133d0,26.1819133d0,26.1819133d0,26.1819133d0/ 
 76: !Q: 
 77: data Cang6(3,:)/-1291.50705d0,-1291.50705d0,-1291.50705d0,-1291.50705d0,-23944.8325d0,7240.4699d0,7240.4699d0,7240.4699d0,7240.4699d0/ 
 78: !D1 and D2: 
 79: data Cang6(4,:)/172.100547d0,172.100547d0,172.100547d0,172.100547d0,4688.06162d0,-1425.39796,-1425.39796,-1425.39796,-1425.39796/ 
 80: data Cang6(5,:)/172.100547d0,172.100547d0,172.100547d0,172.100547d0,4688.06162d0,-1425.39796,-1425.39796,-1425.39796,-1425.39796/ 
 81: !T1 and T2: 
 82: data Cang6(6,:)/methanedof*0.0d0/ 
 83: data Cang6(7,:)/methanedof*0.0d0/ 
 84:  
 85: !H1 and H2: 
 86: data Cang8(1,:)/40.6973228d0,40.6973228d0,40.6973228d0,40.6973228d0,-470.183908d0,212.925753d0,212.925753d0,212.925753d0,212.925753d0/ 
 87: data Cang8(2,:)/40.6973228d0,40.6973228d0,40.6973228d0,40.6973228d0,-470.183908d0,212.925753d0,212.925753d0,212.925753d0,212.925753d0/ 
 88: !Q: 
 89: data Cang8(3,:)/7345.62345d0,7345.62345d0,7345.62345d0,7345.62345d0,132928.009d0,-46577.4854d0,-46577.4854d0,-46577.4854d0,-46577.4854d0/ 
 90: !D1 and D2: 
 91: data Cang8(4,:)/-1195.7863d0,-1195.7863d0,-1195.7863d0,-1195.7863d0,-32880.3016d0,11355.4445d0,11355.4445d0,11355.4445d0,11355.4445d0/ 
 92: data Cang8(5,:)/-1195.7863d0,-1195.7863d0,-1195.7863d0,-1195.7863d0,-32880.3016d0,11355.4445d0,11355.4445d0,11355.4445d0,11355.4445d0/ 
 93: !T1 and T2: 
 94: data Cang8(6,:)/methanedof*0.0d0/ 
 95: data Cang8(7,:)/methanedof*0.0d0/ 
 96:  
 97: !H1 and H2: 
 98: data Cang10(1,:)/-13.9555905d0,-13.9555905d0,-13.9555905d0,-13.9555905d0,334.00843d0,-620.561765d0,-620.561765d0,-620.561765d0,-620.561765d0/ 
 99: data Cang10(2,:)/-13.9555905d0,-13.9555905d0,-13.9555905d0,-13.9555905d0,334.00843d0,-620.561765d0,-620.561765d0,-620.561765d0,-620.561765d0/ 
100: !Q: 
101: data Cang10(3,:)/-12518.903d0,-12518.903d0,-12518.903d0,-12518.903d0,-119240.978d0,62124.298d0,62124.298d0,62124.298d0,62124.298d0/ 
102: !D1 and D2: 
103: data Cang10(4,:)/1655.75062d0,1655.75062d0,1655.75062d0,1655.75062d0,-8388.62553,-13763.1195d0,-13763.1195d0,-13763.1195d0,-13763.1195d0/ 
104: data Cang10(5,:)/1655.75062d0,1655.75062d0,1655.75062d0,1655.75062d0,-8388.62553,-13763.1195d0,-13763.1195d0,-13763.1195d0,-13763.1195d0/ 
105: !T1 and T2: 
106: data Cang10(6,:)/methanedof*0.0d0/ 
107: data Cang10(7,:)/methanedof*0.0d0/ 
108:  
109: !H1 and H2: 
110: data deltaang6(1,:)/7.335799d0,7.335799d0,7.335799d0,7.335799d0,2.825277d0,1.943410d0,1.943410d0,1.943410d0,1.943410d0/ 
111: data deltaang6(2,:)/7.335799d0,7.335799d0,7.335799d0,7.335799d0,2.825277d0,1.943410d0,1.943410d0,1.943410d0,1.943410d0/ 
112: !Q: 
113: data deltaang6(3,:)/4.341591d0,4.341591d0,4.341591d0,4.341591d0,4.288189d0,4.259787d0,4.259787d0,4.259787d0,4.259787d0/ 
114: !D1 and D2: 
115: data deltaang6(4,:)/5.759895d0,5.759895d0,5.759895d0,5.759895d0,6.129260d0,3.737571d0,3.737571d0,3.737571d0,3.737571d0/ 
116: data deltaang6(5,:)/5.759895d0,5.759895d0,5.759895d0,5.759895d0,6.129260d0,3.737571d0,3.737571d0,3.737571d0,3.737571d0/ 
117: !T1 and T2: 
118: data deltaang6(6,:)/methanedof*0.0d0/ 
119: data deltaang6(7,:)/methanedof*0.0d0/ 
120:  
121: !H1 and H2: 
122: !H1, H2, H3, H4, C, M1, M2, M3, M4 
123: data deltaang8(1,:)/1.2192d-2,1.2192d-2,1.2192d-2,1.2192d-2,31.106042,5.747d-3,5.747d-3,5.747d-3,5.747d-3/ 
124: data deltaang8(2,:)/1.2192d-2,1.2192d-2,1.2192d-2,1.2192d-2,31.106042,5.747d-3,5.747d-3,5.747d-3,5.747d-3/ 
125: !Q: 
126: data deltaang8(3,:)/3.643903d0,3.643903d0,3.643903d0,3.643903d0,4.138380d0,4.368576d0,4.368576d0,4.368576d0,4.368576d0/ 
127: !D1 and D2: 
128: data deltaang8(4,:)/4.415080d0,4.415080d0,4.415080d0,4.415080d0,3.962102d0,3.741763d0,3.741763d0,3.741763d0,3.741763d0/ 
129: data deltaang8(5,:)/4.415080d0,4.415080d0,4.415080d0,4.415080d0,3.962102d0,3.741763d0,3.741763d0,3.741763d0,3.741763d0/ 
130: !T1 and T2: 
131: data deltaang8(6,:)/methanedof*0.0d0/ 
132: data deltaang8(7,:)/methanedof*0.0d0/ 
133:  
134: contains 
135:  
136: !------------------------------------------------------------------------------------------- 
137: subroutine wmrb(x, grad, ereal, gradt) 
138: !function that calculates the potential from the water-methane interaction 
139: !for rigid body 
140: !atoms/sites need to be: H H Q D D T T / H H H H C M M M M 
141: implicit none 
142: double precision::    x(48), grad(48), ereal 
143: double precision::    r12, rk 
144: logical::             gradt 
145: integer::             i,j,k, totdof 
146:  
147: !Convert these from crappy units to atomic units 
148: beta(:,:)= betaang(:,:)*0.529177d0 
149: delta6(:,:)= deltaang6(:,:)*0.529177d0 
150: delta8(:,:)= deltaang8(:,:)*0.529177d0 
151: A0(:,:)= Aang0(:,:)*1.59362d-3 
152: AM(:,:)= AangM(:,:)*1.59362d-3/0.529177d0 
153: A1(:,:)= Aang1(:,:)*1.59362d-3*0.529177d0 
154: C6(:,:)= Cang6(:,:)*1.59362d-3/(0.529177d0**6) 
155: C8(:,:)= Cang8(:,:)*1.59362d-3/(0.529177d0**8) 
156: C10(:,:)= Cang10(:,:)*1.59362d-3/(0.529177d0**10) 
157:  
158: totdof=waterdof+methanedof 
159: ! do i=1, totdof 
160: !    write(*,*) (x(3*(i-1)+j), j=1,3) 
161: ! end do 
162: ereal=0.0d0 
163: do i=1,waterdof 
164:    do j=1, methanedof 
165:       r12= calcr(x(3*(i-1)+1:3*(i-1)+3),x(3*(j-1)+22:3*(j-1)+24)) 
166:       ereal= ereal+ tangtoennies(r12, i,j) 
167:    end do 
168: end do 
169:  
170: grad(:)=0.0d0 
171: if (gradt) then 
172:    do i=1,waterdof 
173:       do j=1, methanedof 
174:          r12= calcr(x(3*(i-1)+1:3*(i-1)+3),x(3*(j-1)+22:3*(j-1)+24)) 
175:          do k=1,3 
176:             rk= x(3*(i-1) + k) - x(3*(j-1) + k+ 21) 
177:             grad(3*(i-1) + k) = grad(3*(i-1) + k) + (rk*gradtangtoennies(r12, i, j)/r12) 
178:             grad(3*(j-1) + k + 21) = grad(3*(j-1) + k+21) - (rk*gradtangtoennies(r12, i, j)/r12) 
179:          end do 
180:       end do 
181:    end do 
182: end if 
183:  
184: return 
185: end subroutine wmrb 
186: !------------------------------------------------------------------------------------------- 
187: function calcr(x1, x2) 
188: implicit none 
189: double precision::   x1(:), x2(:), calcr 
190: integer::            i 
191:  
192: calcr= 0.0d0 
193: do i=1,3 
194:    calcr= calcr+(x1(i)-x2(i))**2 
195: end do 
196: calcr= sqrt(calcr) 
197: return 
198: end function calcr 
199: !------------------------------------------------------------------------------------------- 
200: function tangtoennies(r, a, b) 
201: !Function returning the diatomic potential of Tang-Toennies, taking indeces 
202: !a and b as a label to the fitted parameters defined at the start. 
203: implicit none 
204: double precision::    r, eint, tangtoennies 
205: integer::             a, b 
206:  
207: eint= exp(-beta(a,b)*r)*(A0(a,b) + A1(a,b)*r + AM(a,b)/r) 
208: eint= eint + watercharge(a)*methanecharge(b)/(r) 
209: eint= eint + (C6(a,b)/r**6.0d0)*gammp(7.0d0,delta6(a,b)*r) 
210: eint= eint + (C8(a,b)/r**8.0d0)*gammp(9.0d0,delta8(a,b)*r) 
211: eint= eint + (C10(a,b)/r**10.0d0)*gammp(11.0d0,delta8(a,b)*r) 
212:  
213: tangtoennies= eint 
214: return 
215: end function tangtoennies 
216: !------------------------------------------------------------------------------------------- 
217: function gradtangtoennies(r, a, b) 
218: !Function returning the diatomic potential of Tang-Toennies, taking indeces 
219: !a and b as a label to the fitted parameters defined at the start. 
220: implicit none 
221: double precision::    r, grad,gradtangtoennies 
222: integer::             a, b, k 
223:  
224: grad= -beta(a,b)*exp(-beta(a,b)*r)*(A0(a,b) + A1(a,b)*r + AM(a,b)/r) 
225: grad= grad + exp(-beta(a,b)*r)*(A1(a,b) - AM(a,b)/r**2) 
226: grad= grad - watercharge(a)*methanecharge(b)/(r**2) 
227: grad= grad - 6.0d0*(C6(a,b)/r**7.0d0)*gammp(7.0d0,delta6(a,b)*r) 
228: grad= grad - 8.0d0*(C8(a,b)/r**9.0d0)*gammp(9.0d0,delta8(a,b)*r) 
229: grad= grad - 10.0d0*(C10(a,b)/r**11.0d0)*gammp(11.0d0,delta8(a,b)*r) 
230: grad= grad + C6(a,b)*(delta6(a,b)**7)*exp(-delta6(a,b)*r)/exp(gammln(7.0d0)) 
231: grad= grad + C8(a,b)*(delta8(a,b)**9)*exp(-delta8(a,b)*r)/exp(gammln(9.0d0)) 
232: grad= grad + C10(a,b)*(delta8(a,b)**11)*exp(-delta8(a,b)*r)/exp(gammln(11.0d0)) 
233:  
234: gradtangtoennies= grad 
235: return 
236: end function gradtangtoennies 
237:  
238: !----------------------------------------------------------------------------------------- 
239: !Numerical recipes functions to calculate normalized incomplete gamma function P 
240: ! (the Tang-Toennies damping function) 
241:  
242: function gammp(a,x) 
243:   double precision:: a,gammp,x 
244:   !Returns the incomplete gamma function P (a, x). 
245:   double precision:: gammcf,gamser,gln 
246:   if(x.lt.0.d0.or.a.le.0.d0) then 
247:      write(*,*) 'gammp' 
248:      stop 
249:   end if 
250:   if(x.lt.a+1.)then 
251:      ! Use the series representation. 
252:      call gser(gamser,a,x,gln) 
253:      gammp=gamser 
254:   else 
255:      ! Use the continued fraction representation 
256:      call gcf(gammcf,a,x,gln) 
257:      gammp=1.d0-gammcf 
258:      ! and take its complement. 
259:   endif 
260:   return 
261: end function gammp 
262:  
263: subroutine gser(gamser,a,x,gln) 
264:   integer:: ITMAX 
265:   double precision:: a,gamser,gln,x,EPS 
266:   PARAMETER (ITMAX=100,EPS=3.e-7) 
267:   ! Returns the incomplete gamma function P (a, x) evaluated by its series representation as 
268:   ! gamser. Also returns ln Γ(a) as gln. 
269:   integer:: n 
270:   double precision:: ap,del,sum 
271:   gln=gammln(a) 
272:   if(x.le.0.0d0) then 
273:      if(x.lt.0.0d0) stop('gser') 
274:      gamser=0.0d0 
275:      return 
276:   endif 
277:   ap=a 
278:   sum=1.0d0/a 
279:   del=sum 
280:   do n=1,ITMAX 
281:      ap=ap+1. 
282:      del=del*x/ap 
283:      sum=sum+del 
284:      if(abs(del).lt.abs(sum)*EPS) goto 1 
285:   enddo 
286:   write(*,*) 'gser' 
287:   stop 
288: 1  gamser=sum*exp(-x+a*log(x)-gln) 
289:   return 
290: END subroutine gser 
291:  
292: subroutine gcf(gammcf,a,x,gln) 
293:   integer:: itmax 
294:   double precision:: a,gammcf,gln,x,eps,fpmin 
295:   parameter (itmax=100,eps=3.e-7,fpmin=1.e-30) 
296:   ! Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction repre- 
297:   ! sentation as gammcf. Also returns ln Γ(a) as gln . 
298:   ! Parameters: ITMAX is the maximum allowed number of iterations; EPS is the relative accu- 
299:   ! racy; FPMIN is a number near the smallest representable floating-point number. 
300:   integer:: i 
301:   double precision:: an,b,c,d,del,h 
302:   gln=gammln(a) 
303:   b=x+1.-a 
304:   ! Set up for evaluating continued fraction by modified Lentz’s method (§5.2) with b0 = 0. 
305:   c=1.0d0/fpmin 
306:   d=1.0d0/b 
307:   h=d 
308:   do i=1,itmax 
309:      ! Iterate to convergence. 
310:      an=-i*(i-a) 
311:      b=b+2. 
312:      d=an*d+b 
313:      if(abs(d).lt.fpmin) d=fpmin 
314:      c=b+an/c 
315:      if(abs(c).lt.fpmin) c=fpmin 
316:      d=1.0d0/d 
317:      del=d*c 
318:      h=h*del 
319:      if(abs(del-1.d0).lt.EPS) goto 1 
320:   enddo 
321:   stop('gcf') 
322: 1 gammcf=exp(-x+a*log(x)-gln)*h 
323:   ! Put factors in front. 
324:   return 
325: end subroutine gcf 
326:  
327: function gammln(xx) 
328:   double precision:: gammln,xx 
329:   ! Returns the value ln[Gamma(xx)] for xx > 0. 
330:   integer j 
331:   double precision:: ser,stp,tmp,x,y,cof(6) 
332:   ! Internal arithmetic will be done in double precision, a nicety that you can omit if five-figure 
333:   ! accuracy is good enough. 
334:   data cof,stp/76.18009172947146d0,-86.50532032941677d0,& 
335:        24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,& 
336:        -.5395239384953d-5,2.5066282746310005d0/ 
337:   x=xx 
338:   y=x 
339:   tmp=x+5.5d0 
340:   tmp=(x+0.5d0)*log(tmp)-tmp 
341:   ser=1.000000000190015d0 
342:   do j=1,6 
343:      y=y+1.d0 
344:      ser=ser+cof(j)/y 
345:   enddo 
346:   gammln=tmp+log(stp*ser/x) 
347:   return 
348: end function gammln 
349:  
350:  
351: end module watermethane_mod 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0