hdiff output

r29883/commons.f90 2016-03-16 18:33:22.402946093 +0000 r29882/commons.f90 2016-03-16 18:33:23.602958435 +0000
 51:      &                 HISTFACMUL, HPERCENT, AVOIDDIST, MAXERISE, MAXEFALL, TSTART, MATDIFF, STICKYSIG, SDTOL, & 51:      &                 HISTFACMUL, HPERCENT, AVOIDDIST, MAXERISE, MAXEFALL, TSTART, MATDIFF, STICKYSIG, SDTOL, &
 52:      &                 MinimalTemperature, MaximalTemperature, SwapProb, hdistconstraint, COREFRAC, TSTAR, & 52:      &                 MinimalTemperature, MaximalTemperature, SwapProb, hdistconstraint, COREFRAC, TSTAR, &
 53:      &                 RK_R, RK_THETA,ARMA,ARMB, ExtrapolationPercent, lnHarmFreq, PTEMIN, PTEMAX, PTTMIN, PTTMAX, EXCHPROB, & 53:      &                 RK_R, RK_THETA,ARMA,ARMB, ExtrapolationPercent, lnHarmFreq, PTEMIN, PTEMAX, PTTMIN, PTTMAX, EXCHPROB, &
 54:      &                 PTSTEPS, NEQUIL, NQUENCH, COLDFUSIONLIMIT, NEWRES_TEMP, MINOMEGA, LJSIGMA, LJEPSILON, TAUMAX, & 54:      &                 PTSTEPS, NEQUIL, NQUENCH, COLDFUSIONLIMIT, NEWRES_TEMP, MINOMEGA, LJSIGMA, LJEPSILON, TAUMAX, &
 55:      &                 TAUMAXFULL, CPFACTORSG, CPFACTORFG, VGWTOL, ABTHRESH, ACTHRESH, CSMPMAT(3,3), & 55:      &                 TAUMAXFULL, CPFACTORSG, CPFACTORFG, VGWTOL, ABTHRESH, ACTHRESH, CSMPMAT(3,3), &
 56:      &                 RADIUS_CONTAINER, HYDROPHOBIC, RESTRICTREGIONX0, RESTRICTREGIONY0, RESTRICTREGIONZ0, & 56:      &                 RADIUS_CONTAINER, HYDROPHOBIC, RESTRICTREGIONX0, RESTRICTREGIONY0, RESTRICTREGIONZ0, &
 57:      &                 RESTRICTREGIONRADIUS, HARMONICSTR, DUMPUNIQUEEPREV, DUMPUNIQUEEMARKOV, FREEZESAVEE, & 57:      &                 RESTRICTREGIONRADIUS, HARMONICSTR, DUMPUNIQUEEPREV, DUMPUNIQUEEMARKOV, FREEZESAVEE, &
 58:      &                 TBPMIN, TBPSTEP, TBPHF, TBPCF, TBPINCR, SHIFTV, GEOMDIFFTOL, LJATTOC, GCMU, & 58:      &                 TBPMIN, TBPSTEP, TBPHF, TBPCF, TBPINCR, SHIFTV, GEOMDIFFTOL, LJATTOC, GCMU, &
 59:      &                 SRATIO, TRATIO, EXCHINT, DDMCUT, SUMTEMP, SUMSTEP, SUMOSTEP, EXPANDFACTOR, ROTATEFACTOR, EPSRIGID, & 59:      &                 SRATIO, TRATIO, EXCHINT, DDMCUT, SUMTEMP, SUMSTEP, SUMOSTEP, EXPANDFACTOR, ROTATEFACTOR, EPSRIGID, &
 60:      &                 CONTOURBOUNDS(3,2), KCOMP_RIGID, RIGIDCOMDIST, PALPHA, PBETA, PGAMMA, LAT(3,3), MFETPCTL, MFETTRGT, QUIPEQDIST, & 60:      &                 CONTOURBOUNDS(3,2), KCOMP_RIGID, RIGIDCOMDIST, PALPHA, PBETA, PGAMMA, LAT(3,3), MFETPCTL, MFETTRGT, QUIPEQDIST, &
 61:  61: !
 62: !   parameters for anisotropic potentials 62: !   parameters for anisotropic potentials
 63: ! 63: !
 64: !    DC430 > 64: !    DC430 >
 65:      &                 CAPEPS2, CAPRAD, CAPRHO, CAPHEIGHT1, CAPHEIGHT2, & 65:      &                 CAPEPS2, CAPRAD, CAPRHO, CAPHEIGHT1, CAPHEIGHT2, &
 66:      &                 EPSR, GBKAPPA, GBKAPPRM, GBMU,GBNU, GBSIGNOT, GBEPSNOT, GBCHI, GBCHIPRM, & 66:      &                 EPSR, GBKAPPA, GBKAPPRM, GBMU,GBNU, GBSIGNOT, GBEPSNOT, GBCHI, GBCHIPRM, &
 67:      &                 SIGNOT, EPSNOT, SIGMAF, INVKAP, ESA(3), LPRSQ, LSQDFR, GBDPMU, GBDPEPS, GBDPFCT, & 67:      &                 SIGNOT, EPSNOT, SIGMAF, INVKAP, ESA(3), LPRSQ, LSQDFR, GBDPMU, GBDPEPS, GBDPFCT, &
 68:      &                 PYSIGNOT, PYEPSNOT, PYA1(3), PYA2(3), PYDPMU, PYDPEPS, PYDPFCT, & 68:      &                 PYSIGNOT, PYEPSNOT, PYA1(3), PYA2(3), PYDPMU, PYDPEPS, PYDPFCT, &
 69:      &                 LWRCUT, LWCNSTA, LWCNSTB, LWRCUTSQ, LWRCUT2SQ, DELRC, PAPALP, PAPS, PAPCD, PAPEPS, PAPANG1, PAPANG2, &     69:      &                 LWRCUT, LWCNSTA, LWCNSTB, LWRCUTSQ, LWRCUT2SQ, DELRC, PAPALP, PAPS, PAPCD, PAPEPS, PAPANG1, PAPANG2, &    
 70:      &                 DBEPSBB, DBEPSAB, DBSIGBB, DBSIGAB, DBPMU, EFIELD, YKAPPA, YEPS, GEMRC, MREQ, HSEFF, BEPS, & 70:      &                 DBEPSBB, DBEPSAB, DBSIGBB, DBSIGAB, DBPMU, EFIELD, YKAPPA, YEPS, GEMRC, MREQ, HSEFF, BEPS, &
 71:      &                 EPS11, EPS22, EPS12, MRHO11, MRHO22, MRHO12, REQ11, REQ22, REQ12, & 71:      &                 EPS11, EPS22, EPS12, MRHO11, MRHO22, MRHO12, REQ11, REQ22, REQ12, &
108:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&108:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
109:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &109:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &
110:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &110:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
111:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &111:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
112:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &112:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
113:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &113:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
114:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 114:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
115:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &115:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
116:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &116:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
117:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &117:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
118:      &        MLP3T, MKTRAPT, MLPB3T, MULTIPOTT118:      &        MLP3T, MKTRAPT, MLPB3T
119: !119: !
120:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 120:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 
121:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)121:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)
122:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 122:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 
123: 123: 
124: ! csw34> FREEZEGROUP variables124: ! csw34> FREEZEGROUP variables
125: !125: !
126:       INTEGER :: GROUPCENTRE126:       INTEGER :: GROUPCENTRE
127:       DOUBLE PRECISION :: GROUPRADIUS127:       DOUBLE PRECISION :: GROUPRADIUS
128:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE128:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE
188:       LOGICAL :: EXPANDRIGIDT, NORMALISEEXPANDT, DOEXPANDRIGID=.FALSE.188:       LOGICAL :: EXPANDRIGIDT, NORMALISEEXPANDT, DOEXPANDRIGID=.FALSE.
189: 189: 
190: ! ROTATERIGID variables190: ! ROTATERIGID variables
191:       INTEGER :: ROTATERIGIDFREQ191:       INTEGER :: ROTATERIGIDFREQ
192:       LOGICAL :: ROTATERIGIDT, DOROTATERIGID=.FALSE.192:       LOGICAL :: ROTATERIGIDT, DOROTATERIGID=.FALSE.
193: ! TRANSLATERIGID variables193: ! TRANSLATERIGID variables
194:       DOUBLE PRECISION :: TRANSLATEFACTOR = 1D0194:       DOUBLE PRECISION :: TRANSLATEFACTOR = 1D0
195:       INTEGER :: TRANSLATERIGIDFREQ = 1195:       INTEGER :: TRANSLATERIGIDFREQ = 1
196:       LOGICAL :: TRANSLATERIGIDT = .FALSE.,DOTRANSLATERIGID=.FALSE.196:       LOGICAL :: TRANSLATERIGIDT = .FALSE.,DOTRANSLATERIGID=.FALSE.
197: 197: 
198: ! ROTATEHINGE variables 
199:       LOGICAL :: ROTATEHINGET = .FALSE., DOROTATEHINGE = .FALSE. 
200:       INTEGER :: ROTATEHINGEFREQ = 1 
201:       DOUBLE PRECISION :: HINGEROTATEFACTOR = 1.0D0 
202:  
203: ! csw34> HBONDMATRIX VARIABLES198: ! csw34> HBONDMATRIX VARIABLES
204:       LOGICAL :: HBONDMATRIX, HBONDMATCH, HBONDACCEPT, HBONDLIGAND199:       LOGICAL :: HBONDMATRIX, HBONDMATCH, HBONDACCEPT, HBONDLIGAND
205:       CHARACTER(LEN=80) :: HBONDDONORSACCEPTORS, HBONDRESIDUES, NHBONDGROUPSCHAR, HBONDGROUPNAME200:       CHARACTER(LEN=80) :: HBONDDONORSACCEPTORS, HBONDRESIDUES, NHBONDGROUPSCHAR, HBONDGROUPNAME
206:       CHARACTER(LEN=80) :: CHBONDDCUT, HBONDACUT, HBONDDCUT, HBONDACUTLOOSE, HBONDACUTTIGHT201:       CHARACTER(LEN=80) :: CHBONDDCUT, HBONDACUT, HBONDDCUT, HBONDACUTLOOSE, HBONDACUTTIGHT
207:       CHARACTER(LEN=80) :: HBONDDCUTLOOSE, HBONDDCUTTIGHT202:       CHARACTER(LEN=80) :: HBONDDCUTLOOSE, HBONDDCUTTIGHT
208:       CHARACTER(LEN=6) :: HBONDTYPE203:       CHARACTER(LEN=6) :: HBONDTYPE
209:       INTEGER :: NHBONDGROUPS, HBONDLIGANDN204:       INTEGER :: NHBONDGROUPS, HBONDLIGANDN
210:       INTEGER, ALLOCATABLE :: HBONDGROUPS(:,:,:), HBONDMAT(:,:), HBONDGROUPPOP(:),HBONDGROUPIDS(:)205:       INTEGER, ALLOCATABLE :: HBONDGROUPS(:,:,:), HBONDMAT(:,:), HBONDGROUPPOP(:),HBONDGROUPIDS(:)
211:       INTEGER, ALLOCATABLE :: HBONDSOFTMAT(:,:)206:       INTEGER, ALLOCATABLE :: HBONDSOFTMAT(:,:)
212:       DOUBLE PRECISION :: HBONDESAVE, HBONDQUZEROE207:       DOUBLE PRECISION :: HBONDESAVE, HBONDQUZEROE


r29883/genrigid.f90 2016-03-16 18:33:22.622948355 +0000 r29882/genrigid.f90 2016-03-16 18:33:23.826960738 +0000
 20:       INTEGER, ALLOCATABLE :: LRBNPERMGROUP(:), LRBNPERMSIZE(:,:), LRBPERMGROUP(:,:), LRBNSETS(:,:), LRBSETS(:,:,:) 20:       INTEGER, ALLOCATABLE :: LRBNPERMGROUP(:), LRBNPERMSIZE(:,:), LRBPERMGROUP(:,:), LRBNSETS(:,:), LRBSETS(:,:,:)
 21:  21: 
 22: !   vr274:  added lattice coordinates 22: !   vr274:  added lattice coordinates
 23: !           if HAS_LATTICE_COORDS is true, the last two atoms are treated 23: !           if HAS_LATTICE_COORDS is true, the last two atoms are treated
 24: !           as lattice coordintes and rigidcoords is in reduced lattice units 24: !           as lattice coordintes and rigidcoords is in reduced lattice units
 25:       LOGICAL HAS_LATTICE_COORDS 25:       LOGICAL HAS_LATTICE_COORDS
 26:  26: 
 27: !   jdf43:  RIGIDISRIGID = logical array, size NATOMS, TRUE if atom is part of 27: !   jdf43:  RIGIDISRIGID = logical array, size NATOMS, TRUE if atom is part of
 28: !           RB 28: !           RB
 29:       LOGICAL, ALLOCATABLE :: RIGIDISRIGID(:) 29:       LOGICAL, ALLOCATABLE :: RIGIDISRIGID(:)
 30:       INTEGER, ALLOCATABLE :: RB_BY_ATOM(:) ! sn402: records which RB an atom belongs to 
 31: !   jdf43:  FROZENRIGIDBODY = logical array, size NATOMS, TRUE if RB is frozen 30: !   jdf43:  FROZENRIGIDBODY = logical array, size NATOMS, TRUE if RB is frozen
 32:       LOGICAL, ALLOCATABLE :: FROZENRIGIDBODY(:) 31:       LOGICAL, ALLOCATABLE :: FROZENRIGIDBODY(:)
 33:  32: 
 34: !-----------------------------------------------------------------------------------! 33: !-----------------------------------------------------------------------------------!
 35: ! NRIGIDBODY  = number of rigid bodies 34: ! NRIGIDBODY  = number of rigid bodies
 36: ! DEGFREEDOMS = number of degrees of freedom = 6 * NRIGIDBODY + 3 * ADDITIONAL ATOMS 35: ! DEGFREEDOMS = number of degrees of freedom = 6 * NRIGIDBODY + 3 * ADDITIONAL ATOMS
 37: ! MAXSITE     = maximum number of sites in a rigid body 36: ! MAXSITE     = maximum number of sites in a rigid body
 38: ! NRELAXRIGIDR = rigid body minimisation for this number of steps 37: ! NRELAXRIGIDR = rigid body minimisation for this number of steps
 39: ! NRELAXRIGIDA = atom minimisation for this number of steps 38: ! NRELAXRIGIDA = atom minimisation for this number of steps
 40: ! NSITEPERBODY= number of rigid body sites, no need to be the same for all bodies 39: ! NSITEPERBODY= number of rigid body sites, no need to be the same for all bodies
 72:  71: 
 73:   ! hk286 > Allocate NSITEPERBODY 72:   ! hk286 > Allocate NSITEPERBODY
 74:   ALLOCATE (NSITEPERBODY(NRIGIDBODY)) 73:   ALLOCATE (NSITEPERBODY(NRIGIDBODY))
 75:   ALLOCATE (SITESRIGIDBODY(MAXSITE,3,NRIGIDBODY)) 74:   ALLOCATE (SITESRIGIDBODY(MAXSITE,3,NRIGIDBODY))
 76:   ALLOCATE (RIGIDGROUPS(MAXSITE,NRIGIDBODY)) 75:   ALLOCATE (RIGIDGROUPS(MAXSITE,NRIGIDBODY))
 77:   ALLOCATE (REFVECTOR(NRIGIDBODY)) 76:   ALLOCATE (REFVECTOR(NRIGIDBODY))
 78:   ALLOCATE (GR_WEIGHTS(NATOMS)) 77:   ALLOCATE (GR_WEIGHTS(NATOMS))
 79:   ALLOCATE (IINVERSE(NRIGIDBODY,3,3)) 78:   ALLOCATE (IINVERSE(NRIGIDBODY,3,3))
 80: !jdf43> 79: !jdf43>
 81:   ALLOCATE (RIGIDISRIGID(NATOMS)) 80:   ALLOCATE (RIGIDISRIGID(NATOMS))
 82:   ALLOCATE (RB_BY_ATOM(NATOMS)) 
 83:   RIGIDISRIGID=.FALSE. 81:   RIGIDISRIGID=.FALSE.
 84:   RB_BY_ATOM(:) = -1 
 85:  82: 
 86:    ! If ATMASS1 is allocated and using AMBER 9, use ATMASS1 to scale. 83:    ! If ATMASS1 is allocated and using AMBER 9, use ATMASS1 to scale.
 87:    ! If AMBER12_MASSES is allocated and using AMBER 12, use AMBER12_MASSES. 84:    ! If AMBER12_MASSES is allocated and using AMBER 12, use AMBER12_MASSES.
 88:    ! Otherwise, use centre of coordinates. 85:    ! Otherwise, use centre of coordinates.
 89:    IF ( ALLOCATED(ATMASS1) .AND. AMBERT ) THEN 86:    IF ( ALLOCATED(ATMASS1) .AND. AMBERT ) THEN
 90:       ! AMBER 9 atom masses 87:       ! AMBER 9 atom masses
 91:       GR_WEIGHTS = ATMASS1 88:       GR_WEIGHTS = ATMASS1
 92:    ELSE IF ( ALLOCATED(AMBER12_MASSES) .AND. AMBER12T ) THEN 89:    ELSE IF ( ALLOCATED(AMBER12_MASSES) .AND. AMBER12T ) THEN
 93:       ! AMBER 12 atom masses 90:       ! AMBER 12 atom masses
 94:       GR_WEIGHTS = AMBER12_MASSES 91:       GR_WEIGHTS = AMBER12_MASSES
386:      READ(222,*) CHECK1, NSITEPERBODY(J1)383:      READ(222,*) CHECK1, NSITEPERBODY(J1)
387:      DO J2 = 1, NSITEPERBODY(J1)384:      DO J2 = 1, NSITEPERBODY(J1)
388:         READ(222,*) RIGIDGROUPS(J2,J1)385:         READ(222,*) RIGIDGROUPS(J2,J1)
389: ! csw34> check to make sure the current atom is not already in a rigid body386: ! csw34> check to make sure the current atom is not already in a rigid body
390:         IF (RIGIDISRIGID(RIGIDGROUPS(J2,J1))) THEN387:         IF (RIGIDISRIGID(RIGIDGROUPS(J2,J1))) THEN
391:             PRINT *," genrigid> ERROR: atom ",RIGIDGROUPS(J2,J1)," is in multiple rigid bodies! Stopping."388:             PRINT *," genrigid> ERROR: atom ",RIGIDGROUPS(J2,J1)," is in multiple rigid bodies! Stopping."
392:             STOP389:             STOP
393:         ELSE       390:         ELSE       
394: ! csw34> if not, flag the current atom391: ! csw34> if not, flag the current atom
395:             RIGIDISRIGID(RIGIDGROUPS(J2,J1))=.TRUE.392:             RIGIDISRIGID(RIGIDGROUPS(J2,J1))=.TRUE.
396:             RB_BY_ATOM(RIGIDGROUPS(J2,J1)) = J1 
397:         ENDIF393:         ENDIF
398: ! vr274> Moved initialization of coordinates to GENRIGID_INITIALISE, here only read the setup394: ! vr274> Moved initialization of coordinates to GENRIGID_INITIALISE, here only read the setup
399: !        SITESRIGIDBODY(J2,:,J1) = COORDS(3*DUMMY-2:3*DUMMY,1)395: !        SITESRIGIDBODY(J2,:,J1) = COORDS(3*DUMMY-2:3*DUMMY,1)
400:      ENDDO396:      ENDDO
401:   ENDDO397:   ENDDO
402:   CLOSE(222)398:   CLOSE(222)
403: 399: 
404:   CALL GENRIGID_INITIALISE(INICOORDS)400:   CALL GENRIGID_INITIALISE(INICOORDS)
405: END SUBROUTINE GENRIGID_READ_FROM_FILE401: END SUBROUTINE GENRIGID_READ_FROM_FILE
406: 402: 


r29883/isotropic_potentials.f90 2016-03-16 18:33:22.818950372 +0000 r29882/isotropic_potentials.f90 2016-03-16 18:33:24.018962713 +0000
  1: ! This module provides a library of potential functions for isotropic potentials which returning a local copy of the Hessian.  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/isotropic_potentials.f90' in revision 29882
  2: ! It is designed to be used with the MULTIPOT module 
  3:  
  4: !All subroutines in this module must have the following signature to be used with MULTIPOT: 
  5: !        SUBROUTINE POT(TMP_NATOMS, X, PARAMS, TMP_ENERGY, TMP_G, TMP_HESS, GTEST, STEST) 
  6: !            INTEGER, INTENT(IN)           :: TMP_NATOMS 
  7: !            DOUBLE PRECISION, INTENT(IN)  :: X(3*TMP_NATOMS) 
  8: !            DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)      ! Maximum number of parameters is hardcoded here 
  9: !                                                             ! Each potential will use a subset of the elements of PARAMS 
 10: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY 
 11: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*TMP_NATOMS), TMP_HESS(3*TMP_NATOMS,3*TMP_NATOMS) 
 12: !            LOGICAL, INTENT(IN)           :: GTEST, STEST    ! GTEST is true when calculating the gradient as well as energy. 
 13: !                                                             ! STEST is true when calculating the Hessian as well as the other two. 
 14: !        END SUBROUTINE POT 
 15:  
 16: MODULE ISOTROPIC_POTENTIALS 
 17: USE COMMONS, ONLY: MYUNIT 
 18: IMPLICIT NONE 
 19:  
 20: DOUBLE PRECISION, PARAMETER :: WCA_CUT=1.2599210498948731647672106072782283505702514647015079 ! = 2^(1/3) 
 21:  
 22: CONTAINS 
 23:  
 24: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 25: ! Simple LJ potential, no cutoff. Atom radius sigma is given as a variable, but will often be set to 1. 
 26: ! This is almost an exact copy of ljdiff.f, but with sigma added in and the Hessian stored locally to be returned. 
 27: SUBROUTINE ISO_LJ(N, X, PARAMS, TMP_ENERGY, TMP_G, TMP_HESS, GTEST, STEST) 
 28:     IMPLICIT NONE 
 29:  
 30:     INTEGER, INTENT(IN)           :: N ! Number of atoms 
 31:     DOUBLE PRECISION, INTENT(IN)  :: X(3*N) 
 32:     DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)  ! Maximum number of parameters is hardcoded here 
 33:     DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY 
 34:     DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*N), TMP_HESS(3*N,3*N) 
 35:     LOGICAL, INTENT(IN)           :: GTEST, STEST 
 36:  
 37:     ! Various powers of the distance between the atoms, and the atom radius 
 38:     DOUBLE PRECISION :: DIST, R2(N,N), R6, R8(N,N), R14(N,N), SIG, SIG6, SIG12 
 39:     DOUBLE PRECISION :: G(N,N), F(N,N)  ! G tensor and F tensor (see ISOTROPIC_GRAD and ISOTROPIC_HESSIAN, below) 
 40:     INTEGER :: J1, J2, J3, J4 
 41:  
 42:     SIG = PARAMS(1) 
 43:     SIG6 = SIG**6 
 44:     SIG12 = SIG6**2 
 45:  
 46:     TMP_ENERGY=0.0D0 
 47:  
 48:     ! The arrangement of these IF statements seems slightly odd, but it's been done to make sure we only need to set up one pair 
 49:     ! of DO loops for each call. 
 50:     IF (GTEST) THEN 
 51:         IF (STEST) THEN  ! Gradient + Hessian 
 52:             DO J1=1,N 
 53:                 J3=3*J1 
 54:                 R2(J1,J1)=0.0D0 
 55:                 R8(J1,J1)=0.0D0 
 56:                 R14(J1,J1)=0.0D0 
 57:                 G(J1,J1)=0.0D0 
 58:                 F(J1,J1)=0.0D0 
 59:                 DO J2=J1+1,N 
 60:                     J4=3*J2 
 61:  
 62:                     R2(J2,J1)=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
 63:                     R2(J2,J1)=1.0D0/R2(J2,J1) 
 64:                     R6=R2(J2,J1)**3 
 65:                     TMP_ENERGY=TMP_ENERGY+SIG6*R6*(SIG6*R6-1.0D0) 
 66:  
 67:                     ! Set up storage arrays to use in the gradient- and hessian-calculating routines 
 68:                     R8(J2,J1)=R2(J2,J1)**4 
 69:                     R14(J2,J1)=R8(J2,J1)*R8(J2,J1)/R2(J2,J1) 
 70:                     R2(J1,J2)=R2(J2,J1) 
 71:                     G(J2,J1)=-24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2(J1,J2)*SIG6*R6 
 72:                     G(J1,J2)=G(J2,J1) 
 73:                     F(J2,J1)=672.0D0*R14(J2,J1)*SIG12-192.0D0*R8(J2,J1)*SIG6 
 74:                     F(J1,J2)=F(J2,J1) 
 75:                 ENDDO 
 76:             ENDDO 
 77:         ELSE             ! Energy + Gradient only 
 78:             DO J1=1,N 
 79:                 J3=3*J1 
 80:                 G(J1,J1)=0.0D0 
 81:                 DO J2=J1+1,N 
 82:                     J4=3*J2 
 83:  
 84:                     DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
 85:                     DIST=1.0D0/DIST 
 86:                     R6=DIST**3 
 87:  
 88:                     TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) 
 89:  
 90:                     G(J2,J1)=-24.0D0*(2.0D0*R6*SIG6-1.0D0)*DIST*R6*SIG6 
 91:                     G(J1,J2)=G(J2,J1) 
 92:                 ENDDO 
 93:             ENDDO 
 94:         ENDIF 
 95:     ELSE                ! Energy only 
 96:         DO J1=1,N 
 97:             J3=3*J1 
 98:             G(J1,J1)=0.0D0 
 99:             DO J2=J1+1,N 
100:                 J4=3*J2 
101:  
102:                 DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
103:                 DIST=1.0D0/DIST 
104:                 R6=DIST**3 
105:  
106:                 TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) 
107:             ENDDO 
108:         ENDDO 
109:  
110:     ENDIF 
111:     TMP_ENERGY=4.0D0*TMP_ENERGY 
112:  
113:     IF (.NOT.GTEST) RETURN 
114:  
115:     ! G should already be set 
116:     CALL ISOTROPIC_GRAD(N, X, G, TMP_G) 
117:  
118:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient. 
119:  
120:     CALL ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS) 
121:  
122:     RETURN 
123:  
124: END SUBROUTINE ISO_LJ 
125:  
126:  
127: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
128: ! WCA potential, which is simply the LJ potential but truncated at the first minimum of the function and shifted so that 
129: ! the potential is entirely repulsive. 
130: ! The energy and gradient are both continuous at the cutoff. 
131: ! Atom radius sigma is given as a variable (SIG), but will often be set to 1. 
132: SUBROUTINE ISO_WCA(N, X, PARAMS, TMP_ENERGY, TMP_G, TMP_HESS, GTEST, STEST) 
133:     IMPLICIT NONE 
134:  
135:     INTEGER, INTENT(IN)           :: N ! Number of atoms 
136:     DOUBLE PRECISION, INTENT(IN)  :: X(3*N) 
137:     DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)  ! Maximum number of parameters is hardcoded here 
138:     DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY 
139:     DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*N), TMP_HESS(3*N,3*N) 
140:     LOGICAL, INTENT(IN)           :: GTEST, STEST 
141:  
142:     ! Various powers of the distance between the atoms, and the atom radius 
143:     DOUBLE PRECISION :: DIST, R2(N,N), R6, R8(N,N), R14(N,N), SIG, SIG6, SIG12 
144:     DOUBLE PRECISION :: G(N,N), F(N,N)  ! G tensor and F tensor (see ISOTROPIC_GRAD and ISOTROPIC_HESSIAN, below) 
145:     INTEGER :: J1, J2, J3, J4 
146:  
147:     SIG = PARAMS(1) 
148:     SIG6 = SIG**6 
149:     SIG12 = SIG6**2 
150:  
151:     TMP_ENERGY=0.0D0 
152:  
153:     ! The arrangement of these IF statements seems slightly odd, but it's been done to make sure we only need to set up one pair 
154:     ! of DO loops for each call. 
155:     IF (GTEST) THEN 
156:         IF (STEST) THEN  ! Gradient + Hessian 
157:             DO J1=1,N 
158:                 J3=3*J1 
159:                 R2(J1,J1)=0.0D0 
160:                 R8(J1,J1)=0.0D0 
161:                 R14(J1,J1)=0.0D0 
162:                 G(J1,J1)=0.0D0 
163:                 F(J1,J1)=0.0D0 
164:                 DO J2=J1+1,N 
165:                     J4=3*J2 
166:  
167:                     R2(J2,J1)=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
168:  
169:                     IF(R2(J2,J1).GT.WCA_CUT*SIG*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 
170:                     ! We don't compute the energy for this pair of atoms. Set G and F to 0 so that the gradient and 
171:                     ! Hessian terms will go to 0 also. 
172: !                        write(*,*) "E+G+H. Separation>cutoff:", J1, J2, R2(J2,J1) 
173:                         G(J1,J2) = 0.0D0 
174:                         G(J2,J1) = 0.0D0 
175:                         F(J1,J2) = 0.0D0 
176:                         F(J2,J1) = 0.0D0 
177:                         CYCLE 
178:                     ENDIF 
179:  
180:                     R2(J2,J1)=1.0D0/R2(J2,J1) 
181:                     R6=R2(J2,J1)**3 
182:  
183:                     TMP_ENERGY=TMP_ENERGY+SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
184: !                    write(*,*) "E+G+H. Pair, Dist, Energy:", J1, J2, 1.0D0/R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
185:                     ! Set up storage arrays to use in the gradient- and hessian-calculating routines 
186:                     R8(J2,J1)=R2(J2,J1)**4 
187:                     R14(J2,J1)=R8(J2,J1)*R8(J2,J1)/R2(J2,J1) 
188:                     R2(J1,J2)=R2(J2,J1) 
189:                     G(J2,J1)=-24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2(J1,J2)*SIG6*R6 
190:                     G(J1,J2)=G(J2,J1) 
191:                     F(J2,J1)=672.0D0*R14(J2,J1)*SIG12-192.0D0*R8(J2,J1)*SIG6 
192:                     F(J1,J2)=F(J2,J1) 
193:                 ENDDO 
194:             ENDDO 
195:         ELSE             ! Energy + Gradient only 
196:             DO J1=1,N 
197:                 J3=3*J1 
198:                 G(J1,J1)=0.0D0 
199:                 DO J2=J1+1,N 
200:                     J4=3*J2 
201:  
202:                     DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
203:  
204:                     IF(DIST.GT.WCA_CUT*SIG*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 
205: !                        write(*,*) "E+G. Separation>cutoff:", J1, J2, DIST 
206:                         G(J1,J2) = 0.0D0 
207:                         G(J2,J1) = 0.0D0 
208:                         CYCLE 
209:                     ENDIF 
210:  
211:                     DIST=1.0D0/DIST 
212:                     R6=DIST**3 
213:  
214: !                    IF ((J2.GT.165).AND.(J2.LT.176)) THEN 
215: !                        IF (J1.LT.56) THEN 
216: !                            write(*,*) "WCA pair", J1, J2 
217: !                            write(*,*) X(J3-2:J3) 
218: !                            write(*,*) X(J4-2:J4) 
219: !                            write(*,*) "DISTANCE, cutoff", DIST, WCA_CUT*SIG*SIG 
220: !                            write(*,*) "Energy, G", 4*R6*SIG6*(R6*SIG6-1.0D0) + 1.0D0, -24.0D0*(2.0D0*R6*SIG6-1.0D0)*DIST*R6*SIG6 
221: !                        ENDIF 
222: !                    ENDIF 
223:  
224:                     TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) + 0.25D0 
225: !                    write(*,*) "E+G. Pair, Dist, Energy:", J1, J2, (X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2, SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
226:                     G(J2,J1)=-24.0D0*(2.0D0*R6*SIG6-1.0D0)*DIST*R6*SIG6 
227:                     G(J1,J2)=G(J2,J1) 
228:                 ENDDO 
229:             ENDDO 
230:         ENDIF 
231:     ELSE                ! Energy only 
232:         DO J1=1,N 
233:             J3=3*J1 
234:             G(J1,J1)=0.0D0 
235:             DO J2=J1+1,N 
236:                 J4=3*J2 
237:  
238:                 DIST=(X(J3-2)-X(J4-2))**2+(X(J3-1)-X(J4-1))**2+(X(J3)-X(J4))**2 
239:  
240:                 IF(DIST.GT.WCA_CUT*SIG*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 
241: !                    write(*,*) "E. Separation>cutoff:", J1, J2, R2(J2,J1) 
242:                     CYCLE 
243:                 ENDIF 
244:  
245:                 DIST=1.0D0/DIST 
246:                 R6=DIST**3 
247:  
248:                 TMP_ENERGY=TMP_ENERGY+R6*SIG6*(R6*SIG6-1.0D0) + 0.25D0 
249: !                write(*,*) "E. Pair, Dist, Energy:", J1, J2, R2(J2,J1), SIG6*R6*(SIG6*R6-1.0D0) + 0.25D0 
250:             ENDDO 
251:         ENDDO 
252:  
253:     ENDIF 
254:     TMP_ENERGY=4.0D0*TMP_ENERGY 
255:  
256:     IF (.NOT.GTEST) RETURN 
257:  
258:     ! G should already be set 
259:     CALL ISOTROPIC_GRAD(N, X, G, TMP_G) 
260:  
261:     IF (.NOT.STEST) RETURN  ! It is assumed we will never need the Hessian without also needing the gradient. 
262:  
263:     CALL ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS) 
264:  
265:     RETURN 
266:  
267: END SUBROUTINE ISO_WCA 
268:  
269: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
270: ! END OF ISOTROPIC POTENTIALS 
271: ! Next, two functions that are general to all isotropic potentials. To be clear, by "isotropic", I mean "depends only 
272: ! on the distance between the two particles" 
273: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
274:  
275: ! For all isotropic potentials the gradient is calculated in the same way. We simply need to know the positions of the 
276: ! two particles, and the G tensor: G = (1/r)dU/dr for an isotropic potential U(r) 
277: SUBROUTINE ISOTROPIC_GRAD(N, X, G, TMP_G) 
278:     IMPLICIT NONE 
279:     INTEGER, INTENT(IN)           :: N  ! Number of atoms 
280:     DOUBLE PRECISION, INTENT(IN)  :: X(3*N), G(N,N) 
281:     DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*N) 
282:     INTEGER :: J1, J2, J3, J4 
283:     DOUBLE PRECISION :: DUMMYX, DUMMYY, DUMMYZ, XMUL 
284:  
285:     DO J1=1,N  ! The atom for which we are calculating the gradient 
286:         J3=3*J1 
287:         DUMMYX = 0.0D0 
288:         DUMMYY = 0.0D0 
289:         DUMMYZ = 0.0D0 
290:         DO J4=1,N  ! Consider the interaction with each other atom in turn. 
291:             ! This inner loop is the slow bit. Minimise the number of operations required here. 
292:             J2=3*J4 
293:             XMUL=G(J4,J1)  ! Only do the array-access once to save time 
294:             DUMMYX=DUMMYX+XMUL*(X(J3-2)-X(J2-2)) 
295:             DUMMYY=DUMMYY+XMUL*(X(J3-1)-X(J2-1)) 
296:             DUMMYZ=DUMMYZ+XMUL*(X(J3)-X(J2)) 
297:         ENDDO 
298:         TMP_G(J3-2) = DUMMYX 
299:         TMP_G(J3-1) = DUMMYY 
300:         TMP_G(J3) = DUMMYZ 
301:     ENDDO 
302:  
303:     RETURN 
304: END SUBROUTINE ISOTROPIC_GRAD 
305:  
306: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
307: ! For all isotropic potentials the Hessian is calculated in the same way. We simply need to know the positions of the 
308: ! two particles, the G tensor (see above) and the F tensor: F = r*d[(1/r)dU/dr]/dr 
309: SUBROUTINE ISOTROPIC_HESSIAN(N, X, G, F, R2, TMP_HESS) 
310:     IMPLICIT NONE 
311:     INTEGER, INTENT(IN)           :: N  ! Number of atoms 
312:     DOUBLE PRECISION, INTENT(IN)  :: X(3*N), G(N,N), F(N,N), R2(N,N) 
313:     DOUBLE PRECISION, INTENT(OUT) :: TMP_HESS(3*N,3*N) 
314:     INTEGER :: J1, J2, J3, J4, J5, J6 
315:     DOUBLE PRECISION :: DUMMY 
316:  
317: !       Compute the hessian. First are the entirely diagonal terms (3N of them) 
318: !       These correspond to 2nd derivatives wrt the same coordinate twice 
319:     DO J1=1,N 
320:         DO J2=1,3 
321:             J3=3*(J1-1)+J2 
322:             DUMMY=0.0D0 
323:             DO J4=1,N 
324:                 DUMMY=DUMMY+F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2 + G(J4,J1) 
325:             ENDDO 
326:             TMP_HESS(J3,J3)=DUMMY 
327:         ENDDO 
328:     ENDDO 
329:  
330: !  Next are the terms where x_i and x_j are on the same atom 
331: !  but are different, e.g. y and z. 
332:  
333:     DO J1=1,N 
334:         DO J2=1,3 
335:             J3=3*(J1-1)+J2 
336:             DO J4=J2+1,3 
337:                 DUMMY=0.0D0 
338:                 DO J5=1,N 
339:                     DUMMY=DUMMY + F(J5,J1)*R2(J5,J1)*(X(J3)-X(3*(J5-1)+J2))*(X(3*(J1-1)+J4)-X(3*(J5-1)+J4)) 
340:                 ENDDO 
341:                 TMP_HESS(3*(J1-1)+J4,J3)=DUMMY 
342:             ENDDO 
343:         ENDDO 
344:     ENDDO 
345:  
346: !  Case III, different atoms, same cartesian coordinate. 
347:  
348:     DO J1=1,N 
349:         DO J2=1,3 
350:             J3=3*(J1-1)+J2 
351:             DO J4=J1+1,N 
352:                 TMP_HESS(3*(J4-1)+J2,J3)=-F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))**2-G(J4,J1) 
353:            ENDDO 
354:         ENDDO 
355:     ENDDO 
356:  
357: !  Case IV: different atoms and different cartesian coordinates. 
358:  
359:     DO J1=1,N 
360:         DO J2=1,3 
361:             J3=3*(J1-1)+J2 
362:             DO J4=J1+1,N 
363:                 DO J5=1,J2-1 
364:                     J6=3*(J4-1)+J5 
365:                     TMP_HESS(J6,J3)=-F(J4,J1)*R2(J4,J1)*(X(J3)-X(3*(J4-1)+J2))*(X(3*(J1-1)+J5)-X(J6)) 
366:                     TMP_HESS(3*(J4-1)+J2,3*(J1-1)+J5)=TMP_HESS(J6,J3) 
367:                 ENDDO 
368:             ENDDO 
369:         ENDDO 
370:     ENDDO 
371:  
372: !  Symmetrise Hessian 
373:  
374:     DO J1=1,3*N 
375:         DO J2=J1+1,3*N 
376:             TMP_HESS(J1,J2)=TMP_HESS(J2,J1) 
377:         ENDDO 
378:     ENDDO 
379:  
380:     RETURN 
381: END SUBROUTINE ISOTROPIC_HESSIAN 
382:  
383: END MODULE ISOTROPIC_POTENTIALS 


r29883/multipot.f90 2016-03-16 18:33:23.014952388 +0000 r29882/multipot.f90 2016-03-16 18:33:24.210964687 +0000
  1: ! Module that allows different potential types to be specified for different atoms in a single system  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/multipot.f90' in revision 29882
  2: MODULE MULTIPOT 
  3:  
  4:     USE PAIRWISE_POTENTIALS 
  5:     USE ISOTROPIC_POTENTIALS 
  6:     USE COMMONS, ONLY: DEBUG, MYUNIT 
  7:  
  8:     IMPLICIT NONE 
  9:  
 10:     ! We pre-calculate as much as possible, hence the need for lots of different storage arrays, documented here. 
 11:     INTEGER :: NPOTTYPES = 0                        ! The number of different potential types being used 
 12:     CHARACTER(LEN=10), ALLOCATABLE :: POTTYPES(:)   ! An array of identifiers for the different types of potential being 
 13:                                                     ! used for this particular system 
 14:     DOUBLE PRECISION, ALLOCATABLE :: POTSCALES(:)   ! A list of the energy scale factors required to match the energy units of the 
 15:                                                     ! different potentials 
 16:     DOUBLE PRECISION, ALLOCATABLE :: POTPARAMS(:,:) ! Holds the parameters required for each potential type 
 17:     INTEGER, ALLOCATABLE :: NATOM_BY_POT(:)         ! Holds the number of atoms using each potential 
 18:  
 19:  
 20:                                                     ! The next two variables are for use by pairwise potentials only (the entries in 
 21:                                                     ! these arrays for non-pairwise potentials must be filled by dummies) 
 22:     INTEGER, ALLOCATABLE :: N_ATOM_PARTNERS(:,:)    ! Holds the number of interaction partners for each atom, sorted by potential type 
 23:     INTEGER, ALLOCATABLE :: POTLISTS(:,:,:)         ! Holds the indices of all atoms belonging to each potential type, and 
 24:                                                     ! the indices of the partners which interact with each atom 
 25:                                                     ! POTLISTS(m,n,1) contains the index of the nth atom which has an attraction of 
 26:                                                     ! type POTTYPES(m), and POTLISTS(m,n,2:N_ATOM_PARTNERS(m,n)) are the indices of 
 27:                                                     ! its partners. 
 28:  
 29:     INTEGER, ALLOCATABLE :: DEGS_BY_POT(:,:)        ! For use by isotropic potentials, where all atoms interact with all other atoms. 
 30:                                                     ! Contains a list of the degrees of freedom belonging to all the atoms using this 
 31:                                                     ! potential. 
 32:  
 33: CONTAINS 
 34: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 35:     ! Parse the input file which is required to specify the different types of interaction we have in the system 
 36: SUBROUTINE MULTIPOT_INITIALISE 
 37:  
 38:     USE COMMONS, ONLY: NATOMS 
 39:     USE GENRIGID, ONLY: RIGIDINIT, RB_BY_ATOM 
 40:  
 41:     IMPLICIT NONE 
 42:  
 43:     INTEGER NPARAMS, ATOM1, ATOM2 
 44:     INTEGER :: ATOMLIST(NATOMS) 
 45:     CHARACTER(LEN=1000) :: DUMMYCHAR 
 46:     LOGICAL :: END 
 47:     INTEGER :: J1, J2, J3, iostatus, COUNTER 
 48:  
 49:     ! Input file: multipotconfig 
 50:     ! Format as follows. 
 51:     ! 
 52:     ! Introduce each type of interaction with a line: 'POT' 
 53:     ! To comment out a potential type, simply insert a '#' at the start of the 'POT' line. 
 54:     ! The next line has the form: POTTYPE n scale nparams 
 55:     !   POTTYPE is a string specifying the type of potential 
 56:     !   n is the number of atoms which will use this type of potential 
 57:     !   scale is the factor by which the energy of this interaction must be scaled to match the others 
 58:     !   nparams is the number of arguments to the pairwise potential, which will be potential specific. nparams<=10, currently. 
 59:     ! The next line is the list of input parameters (up to 10 of them), separated by spaces.  
 60:     !   Note that nparams needs to be at least 1. If your potential has no extra parameters, set nparams to 1 
 61:     !   and have this line contain a single "0", which you then don't need to use. 
 62:     ! The next line(s) contain lists of the atoms which will interact according to the current potential. 
 63:     !   See below for example formats which may be used. 
 64:  
 65: !   Read multipotconfig once to count the number of types of potential 
 66:     NPOTTYPES=0 
 67:     OPEN(UNIT=22,FILE='multipotconfig',status='old') 
 68:     DO 
 69:        READ(22,*,IOSTAT=iostatus) DUMMYCHAR 
 70:        IF (iostatus<0) THEN 
 71:           CLOSE(22) 
 72:           EXIT 
 73:        ELSE IF (DUMMYCHAR(1:3).EQ.'POT') THEN 
 74:           NPOTTYPES = NPOTTYPES + 1 
 75:        ENDIF 
 76:     END DO 
 77:     CLOSE(22) 
 78:  
 79:     ALLOCATE(POTTYPES(NPOTTYPES)) 
 80:     ALLOCATE(NATOM_BY_POT(NPOTTYPES)) 
 81:     ALLOCATE(DEGS_BY_POT(NPOTTYPES,3*NATOMS)) 
 82:     ALLOCATE(POTSCALES(NPOTTYPES)) 
 83:     ALLOCATE(POTPARAMS(NPOTTYPES,10))  ! Currently each potential is allowed to have up to 10 parameters 
 84:  
 85: !   Determine the names of the potentials and the number of atoms using each potential 
 86:     COUNTER = 1 ! Counts the number of types of potential that have been read thus far 
 87:     OPEN(UNIT=22,FILE='multipotconfig',status='old') 
 88:     DO 
 89:         READ(22,*,IOSTAT=iostatus) DUMMYCHAR 
 90:         IF (iostatus<0) THEN 
 91:             CLOSE(22) 
 92:             EXIT 
 93:         ELSE IF (DUMMYCHAR(1:3).EQ.'POT') THEN ! Read in information and parameters associated with this potential 
 94:             ! NPARAMS is the number of parameters to be read from the next line 
 95:             READ(22,*) POTTYPES(COUNTER), NATOM_BY_POT(COUNTER), POTSCALES(COUNTER), NPARAMS 
 96:             READ(22,*) POTPARAMS(COUNTER,:NPARAMS) 
 97:             IF(NATOM_BY_POT(COUNTER).GT.NATOMS) WRITE(MYUNIT,*) "WARNING: NATOM_BY_POT for potential ", POTTYPES(COUNTER), & 
 98:                 "is larger than NATOMS" 
 99:  
100:             COUNTER = COUNTER + 1 
101:         ENDIF 
102:     END DO 
103:     CLOSE(22) 
104:  
105:     ALLOCATE(N_ATOM_PARTNERS(NPOTTYPES,MAXVAL(NATOM_BY_POT))) 
106:     ALLOCATE(POTLISTS(NPOTTYPES,MAXVAL(NATOM_BY_POT),NATOMS)) 
107:  
108:     ! Now the important bit: read in the actual atom lists 
109:     OPEN(UNIT=22,FILE='multipotconfig',status='unknown') 
110:     DO J1 = 1, NPOTTYPES 
111:         ! Read the three header lines 
112:         READ(22,*) DUMMYCHAR 
113:  
114:         IF (DUMMYCHAR(1:1).EQ.'#') THEN  ! This potential is commented out: skip it. 
115:             IF(DEBUG) WRITE(MYUNIT,*) "Skipping commented block in multipotconfig" 
116:             DO 
117:                 READ(22,*) DUMMYCHAR 
118:                 IF (DUMMYCHAR(1:3).EQ.'POT') EXIT  ! We have found the start of the next uncommented header 
119:             ENDDO 
120:         ENDIF 
121:  
122:         READ(22,*) DUMMYCHAR 
123:         READ(22,*) DUMMYCHAR 
124:         ! Now read in the atom numbers from the following line(s). 
125:         ! The required input format will vary between potentials, although most will use the format described under CASE DEFAULT 
126:  
127:         SELECT CASE(POTTYPES(J1)) 
128:  
129:         ! For potentials which only operate between specific pairs of atoms, such as harmonic springs, each line should contain an interacting pair. 
130:         CASE('HSPR') 
131:             DO J2 = 1, NATOM_BY_POT(J1)/2  ! J2 ranges over the number of springs (Note, NATOM_BY_POT contains the actual no of atoms) 
132:                 ! We only want to compute the potential once for each spring. 
133:  
134:                 READ(22,*) ATOM1, ATOM2 
135:  
136:                 POTLISTS(J1,2*J2-1,1) = ATOM1  ! Make an entry for each of these atoms in POTLISTS 
137:                 POTLISTS(J1,2*J2,1) = ATOM2 
138:  
139:                 N_ATOM_PARTNERS(J1,2*J2-1)=1 ! ATOM1 partners with ATOM2... 
140:                 POTLISTS(J1,2*J2-1,2) = ATOM2 
141:                 N_ATOM_PARTNERS(J1,2*J2)=0   ! ...but atom B doesn't partner with any (we've already included its interaction with A) 
142:             ENDDO 
143:  
144:             IF(DEBUG) THEN 
145:                 WRITE(MYUNIT,*) "Potential:", POTTYPES(J1) 
146:                 WRITE(MYUNIT,*) "N_ATOM_PARTNERS:" 
147:                 WRITE(MYUNIT,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1)) 
148:                 WRITE(MYUNIT,*) "POTLISTS:" 
149:                 DO J2 = 1,NATOM_BY_POT(J1) 
150:                     WRITE(MYUNIT,*) "Atom number", J2, "in this potential" 
151:                     WRITE(MYUNIT,*) POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2)) 
152:                 ENDDO 
153:             ENDIF 
154:  
155:         ! For "iso_" potentials. Every specified atom interacts with all the others, but this is done by means of a single 
156:         ! call to a hard-coded version of the potential (e.g. ISO_LJ) rather than a series of calls to PAIRWISE_LJ for every pair of 
157:         ! interacting atoms. This is less flexible, but significantly more efficient for large systems. 
158:         ! The input format is straightforward: simply list all the atoms using this potential on a single line, separated by spaces. 
159:         ! They don't have to be in index order, but everything will make more sense if they are! 
160:         CASE('ILJ','IWCA') 
161:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1)) 
162:  
163:             ! All we need to save for this type of potential is the list of whole-system degrees of freedom on which 
164:             ! the potential will depend. 
165:             DO J2=1,NATOM_BY_POT(J1) 
166:                 J3 = ATOMLIST(J2) 
167:                 DEGS_BY_POT(J1,3*J2-2) = 3*J3-2 
168:                 DEGS_BY_POT(J1,3*J2-1) = 3*J3-1 
169:                 DEGS_BY_POT(J1,3*J2)   = 3*J3 
170:             ENDDO 
171:  
172:             IF(DEBUG) THEN 
173:                 write(MYUNIT,*) "Potential:", POTTYPES(J1) 
174:                 write(MYUNIT,*) "DEGS_BY_POT:" 
175:                 write(MYUNIT,*) DEGS_BY_POT(J1,:NATOM_BY_POT(J1)*3) 
176:             ENDIF 
177:  
178:         CASE DEFAULT 
179:         ! Default: an isotropic potential handled pair-by-pair. NOTE: likely to be very inefficient for large systems. Use the "iso_" 
180:         ! version of the relevant potential instead (which has the same input format). 
181:         ! All atoms with this potential type interact with all others - unless they are in the same rigid body. 
182:         ! All atoms interacting according to that potential should be specified on a single line, separated by spaces 
183:         ! They don't have to be in index order, but everything will make more sense if they are! 
184:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1)) 
185:  
186:             IF(RIGIDINIT .AND. .FALSE.) THEN   ! We can avoid calculate interactions within a rigid body by only including atom partner pairs 
187:                                                ! that involve atoms from two different bodies 
188:                 DO J2=1,NATOM_BY_POT(J1) 
189:                     POTLISTS(J1,J2,1) = ATOMLIST(J2)   ! Remember, the first element in POTLISTS(x,y,:) is the index of atom y in the list 
190:                                                        ! for potential type x. The remaining elements are the partners of atom y. 
191:                     N_ATOM_PARTNERS(J1,J2) = 0 
192:                     DO J3=J2,NATOM_BY_POT(J1) 
193:                     ! So that each pair only gets calculated once, each partner pair is only included once. So the 2nd atom 
194:                     ! in each pair to appear in the list gets listed as a partner to the 1st atom, but not the other way round. 
195:                         IF (RB_BY_ATOM(ATOMLIST(J3)) .NE. RB_BY_ATOM(ATOMLIST(J2))) THEN 
196:                             ! Only add partners that belong to different rigid bodies 
197:                             N_ATOM_PARTNERS(J1,J2) = N_ATOM_PARTNERS(J1,J2) + 1 
198:                             POTLISTS(J1,J2,1+N_ATOM_PARTNERS(J1,J2)) = ATOMLIST(J3) 
199:                         ENDIF 
200:                     ENDDO 
201:                 ENDDO 
202:             ELSE   ! No rigid bodies so all atoms in ATOMLIST should partner with all others. 
203:                 DO J2 = 1,NATOM_BY_POT(J1)  ! We still only add each pair once. 
204:                     N_ATOM_PARTNERS(J1,J2) = NATOM_BY_POT(J1)-J2  ! So each atom partners with every atom that comes later in the list 
205:                     POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2)+1)=ATOMLIST(J2:NATOM_BY_POT(J1)) 
206:                 ENDDO 
207:             ENDIF 
208:  
209:             IF(DEBUG) THEN 
210:                 WRITE(MYUNIT,*) "Potential:", POTTYPES(J1) 
211:                 WRITE(MYUNIT,*) "N_ATOM_PARTNERS:" 
212:                 WRITE(MYUNIT,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1)) 
213:                 WRITE(MYUNIT,*) "POTLISTS:" 
214:                 DO J2 = 1,NATOM_BY_POT(J1) 
215:                     WRITE(MYUNIT,*) "Atom number", J2, "in this potential" 
216:                     WRITE(MYUNIT,*) POTLISTS(J1,J2,N_ATOM_PARTNERS(J1,J2)) 
217:                 ENDDO 
218:             ENDIF 
219:  
220:         END SELECT 
221:     ENDDO 
222:     CLOSE(22) 
223:  
224: END SUBROUTINE MULTIPOT_INITIALISE 
225:  
226: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
227: ! Calculate the potential energy and the gradient (if GTEST is set) and hessian (if STEST) 
228: SUBROUTINE MULTIPOT_CALL (X, G, ENERGY, GTEST, STEST) 
229:     USE COMMONS, ONLY: NATOMS 
230:     USE MODHESS 
231:  
232:     IMPLICIT NONE 
233:  
234:     DOUBLE PRECISION, INTENT(IN)  :: X(3*NATOMS) 
235:     DOUBLE PRECISION, INTENT(OUT) :: ENERGY, G(3*NATOMS) 
236:     LOGICAL, INTENT(IN)           :: GTEST, STEST 
237:     INTEGER                       :: J1 
238:  
239:     ENERGY  = 0.D0 
240:     IF (GTEST) G(:) = 0.D0 
241:     IF (STEST) THEN 
242:         HESS(:,:) = 0.D0 
243:     ENDIF 
244:  
245:     ! Cycle through the different types of interaction and perform the potential call for all the corresponding atoms 
246:     DO J1 = 1, NPOTTYPES 
247:         SELECT CASE(POTTYPES(J1)) 
248:             CASE('ILJ') 
249:                 ! Only one parameter for this potential: the particle radius, sigma 
250:                 CALL COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, ISO_LJ, J1, NATOM_BY_POT(J1)) 
251:             CASE('IWCA') 
252:                 ! Only one parameter for this potential: the particle radius, sigma 
253:                 CALL COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, ISO_WCA, J1, NATOM_BY_POT(J1)) 
254:             CASE('PLJ') 
255:                 ! Only one parameter for this potential: the particle radius, sigma 
256:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_LJ, J1) 
257:             CASE('PWCA') 
258:                 ! Only one parameter for this potential: the particle radius, sigma 
259:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, PAIRWISE_WCA, J1) 
260:             CASE('HSPR') 
261:                 ! Parameter is the equilibrium bond length, R0. Energy is returned in units of k_spr. 
262:                 CALL COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, HARMONIC_SPRINGS, J1) 
263:             CASE DEFAULT 
264:                 WRITE(MYUNIT,*) "multipot> Error: unspecified potential" 
265:                 STOP 
266:         END SELECT 
267:     ENDDO 
268:  
269: END SUBROUTINE MULTIPOT_CALL 
270:  
271: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
272: !   Evaluate the energy(+gradient(+hessian)) for a set of atoms interacting according to this particular potential 
273: SUBROUTINE COMPUTE_ISOTROPIC_POTENTIAL(X, G, ENERGY, GTEST, STEST, POT, POTID, TMP_NATOMS) 
274:     USE COMMONS, ONLY: NATOMS 
275:     USE MODHESS 
276:     IMPLICIT NONE 
277:     DOUBLE PRECISION, INTENT(IN)    :: X(3*NATOMS) 
278:     DOUBLE PRECISION, INTENT(INOUT) :: G(3*NATOMS) 
279:     DOUBLE PRECISION, INTENT(INOUT) :: ENERGY 
280:     LOGICAL, INTENT(IN)             :: GTEST, STEST 
281:     INTEGER, INTENT(IN)             :: POTID, TMP_NATOMS 
282:     DOUBLE PRECISION                :: TMP_X(3*TMP_NATOMS), TMP_G(3*TMP_NATOMS), TMP_HESS(3*TMP_NATOMS,3*TMP_NATOMS), TMP_E 
283:     INTEGER                         :: J1, J2, NDEGS 
284:  
285:     ! Interface to the potential type which is passed in 
286:     INTERFACE 
287:         SUBROUTINE POT(TMP_NATOMS, X, PARAMS, TMP_ENERGY, TMP_G, TMP_HESS, GTEST, STEST) 
288:             INTEGER, INTENT(IN)           :: TMP_NATOMS 
289:             DOUBLE PRECISION, INTENT(IN)  :: X(3*TMP_NATOMS) 
290:             DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)  ! Maximum number of parameters is hardcoded here 
291:             DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY 
292:             DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*TMP_NATOMS), TMP_HESS(3*TMP_NATOMS,3*TMP_NATOMS) 
293:             LOGICAL, INTENT(IN)           :: GTEST, STEST 
294:         END SUBROUTINE POT 
295:     END INTERFACE 
296:  
297:     NDEGS = 3*TMP_NATOMS 
298:  
299:     DO J1 = 1, NDEGS  ! Loop over all the atoms with this kind of potential 
300:         TMP_X(J1) = X(DEGS_BY_POT(POTID,J1)) 
301:     ENDDO 
302:  
303: !    IF(DEBUG) THEN 
304: !        WRITE(MYUNIT,*) "Calling potential", POTTYPES(J1) 
305: !        WRITE(MYUNIT,*) "Degrees of freedom being used:" 
306: !        WRITE(MYUNIT,*) DEGS_BY_POT(POTID,:NDEGS) 
307: !    ENDIF 
308:  
309:     ! The advantage of this slightly tortuous method is that we only need to call POT once. Compare with the pairwise 
310:     ! routine, where we need to make a function call for every pair of atoms in the system. 
311:     CALL POT(TMP_NATOMS, TMP_X, POTPARAMS(POTID,:), TMP_E, TMP_G, TMP_HESS, GTEST, STEST) 
312:  
313:     ENERGY = ENERGY + TMP_E*POTSCALES(POTID) 
314:  
315:     ! Unpack the gradient, if it is being calculated 
316:     IF(GTEST) THEN 
317:         DO J1 = 1,NDEGS 
318:             G(DEGS_BY_POT(POTID,J1)) = G(DEGS_BY_POT(POTID,J1)) + TMP_G(J1)*POTSCALES(POTID) 
319:         ENDDO 
320:     ENDIF 
321:  
322:     ! Unpack the Hessian, if it is being calculated 
323:     IF(STEST) THEN 
324:         DO J1 = 1,NDEGS 
325:             DO J2 = 1,NDEGS 
326:                 HESS(DEGS_BY_POT(POTID,J1), DEGS_BY_POT(POTID,J2)) = HESS(DEGS_BY_POT(POTID,J1), DEGS_BY_POT(POTID,J2)) & 
327:                                                                                + TMP_HESS(J1,J2)*POTSCALES(POTID) 
328:             ENDDO 
329:         ENDDO 
330:     ENDIF 
331:  
332: END SUBROUTINE COMPUTE_ISOTROPIC_POTENTIAL 
333:  
334: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
335: !   Evaluate the energy(+gradient(+hessian)) for each pair of atoms listed as interacting according to this particular potential 
336: SUBROUTINE COMPUTE_PAIRWISE_POTENTIAL(X, G, ENERGY, GTEST, STEST, POT, POTID) 
337:     USE COMMONS, ONLY: NATOMS 
338:     USE MODHESS 
339:     IMPLICIT NONE 
340:     DOUBLE PRECISION, INTENT(IN)    :: X(3*NATOMS) 
341:     DOUBLE PRECISION, INTENT(INOUT) :: G(3*NATOMS) 
342:     DOUBLE PRECISION, INTENT(INOUT) :: ENERGY 
343:     LOGICAL, INTENT(IN)             :: GTEST, STEST 
344:     INTEGER, INTENT(IN)             :: POTID 
345:     DOUBLE PRECISION                :: TMP_E, TMP_G(6), TMP_HESS(6,6) 
346:     INTEGER                         :: J1, J2, J3, ATOM1, ATOM2 
347:  
348:     ! Interface to the potential type which is passed in 
349:     INTERFACE 
350:         SUBROUTINE POT(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST) 
351:             DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3) 
352:             DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)  ! Maximum number of parameters is hardcoded here 
353:             DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY 
354:             DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6) 
355:             LOGICAL, INTENT(IN)           :: GTEST, STEST 
356:         END SUBROUTINE POT 
357:     END INTERFACE 
358:  
359:     DO J1 = 1, NATOM_BY_POT(POTID)  ! Loop over all the atoms with this kind of potential 
360:         ATOM1 = POTLISTS(POTID,J1,1) 
361:         DO J2 = 1, N_ATOM_PARTNERS(POTID,J1)  ! Loop over every partner with which this atom must interact 
362:             ATOM2 = POTLISTS(POTID,J1,1+J2) 
363:  
364:             !!!!!! The actual potential call! !!!!! 
365:             CALL POT(X(3*ATOM1-2:3*ATOM1),X(3*ATOM2-2:3*ATOM2),POTPARAMS(POTID,:),TMP_G,TMP_E,TMP_HESS,GTEST,STEST) 
366:             !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
367:  
368:             ENERGY = ENERGY + TMP_E*POTSCALES(POTID) 
369:  
370:             ! Unpack the pairwise gradient, if required 
371:             IF(GTEST) THEN 
372:                 G(3*ATOM1-2:3*ATOM1) = G(3*ATOM1-2:3*ATOM1) + TMP_G(:3)*POTSCALES(POTID) 
373:                 G(3*ATOM2-2:3*ATOM2) = G(3*ATOM2-2:3*ATOM2) + TMP_G(4:)*POTSCALES(POTID) 
374:  
375:                 ! Unpack the pairwise Hessian, if required 
376:                 IF(STEST) THEN 
377:                     DO J3 = 1, 3 
378:                         ! On-diagonal block for atom 1 
379:                         HESS(3*(ATOM1-1)+J3,3*ATOM1-2:3*ATOM1) = HESS(3*(ATOM1-1)+J3,3*ATOM1-2:3*ATOM1) + & 
380:                                                                         TMP_HESS(J3,1:3)*POTSCALES(POTID) 
381:                         ! On-diagonal block for atom 2 
382:                         HESS(3*(ATOM2-1)+J3,3*ATOM2-2:3*ATOM2) = HESS(3*(ATOM2-1)+J3,3*ATOM2-2:3*ATOM2) + & 
383:                                                                         TMP_HESS(J3+3,4:6)*POTSCALES(POTID) 
384:                         ! Top-right off-diagonal block 
385:                         HESS(3*(ATOM1-1)+J3,3*ATOM2-2:3*ATOM2) = HESS(3*(ATOM1-1)+J3,3*ATOM2-2:3*ATOM2) + & 
386:                                                                         TMP_HESS(J3,4:6)*POTSCALES(POTID) 
387:                         ! Bottom-left off-diagonal block 
388:                         HESS(3*(ATOM2-1)+J3,3*ATOM1-2:3*ATOM1) = HESS(3*(ATOM2-1)+J3,3*ATOM1-2:3*ATOM1) + & 
389:                                                                         TMP_HESS(J3+3,1:3)*POTSCALES(POTID) 
390:                     ENDDO 
391:                 ENDIF 
392:  
393:             ENDIF 
394:  
395:         ENDDO 
396:     ENDDO 
397:  
398: END SUBROUTINE COMPUTE_PAIRWISE_POTENTIAL 
399:  
400: END MODULE MULTIPOT 


r29883/pairwise_potentials.f90 2016-03-16 18:33:23.210954403 +0000 r29882/pairwise_potentials.f90 2016-03-16 18:33:24.402966662 +0000
  1: ! This module provides a library of potential functions which operate only on a pair of atoms  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/pairwise_potentials.f90' in revision 29882
  2: ! It is designed to be used with the MULTIPOT module 
  3:  
  4: !All subroutines in this module must have the following signature to be used with MULTIPOT: 
  5: !        SUBROUTINE POT(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST) 
  6: !            DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3) 
  7: !            DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)      ! Maximum number of parameters is hardcoded here 
  8: !                                                             ! Each potential will use a subset of the elements of PARAMS 
  9: !            DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY 
 10: !            DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(9)! Gradient and energy for the subsystem composed of this pair 
 11: !            LOGICAL, INTENT(IN)           :: GTEST, STEST    ! GTEST is true when calculating the gradient as well as energy. 
 12: !                                                             ! STEST is true when calculating the Hessian as well as the other two. 
 13: !        END SUBROUTINE POT 
 14:  
 15: MODULE PAIRWISE_POTENTIALS 
 16:  
 17: IMPLICIT NONE 
 18:  
 19: DOUBLE PRECISION, PARAMETER :: WCA_CUT=1.2599210498948731647672106072782283505702514647015079 ! = 2^(1/3) 
 20:  
 21: CONTAINS 
 22:  
 23: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 24: ! Simple LJ potential, no cutoff. Atom radius sigma is given as a variable, but will often be set to 1. 
 25: SUBROUTINE PAIRWISE_LJ(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST) 
 26:     IMPLICIT NONE 
 27:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10) 
 28:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY 
 29:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6) 
 30:     LOGICAL, INTENT(IN)           :: GTEST, STEST 
 31:  
 32:     DOUBLE PRECISION :: R2, R6, R8, R14, SIG, SIG6, SIG12  ! Various powers of the distance between the atoms, and the atom radius 
 33:     DOUBLE PRECISION :: G, F  ! G tensor and F tensor (see PAIRWISE_GRAD and PAIRWISE_HESSIAN, below) 
 34:  
 35:     SIG = PARAMS(1) 
 36:     SIG6 = SIG**6 
 37:     SIG12 = SIG6**2 
 38:  
 39:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2 
 40:     R2 = 1.0D0/R2 
 41:     R6 = R2**3 
 42:  
 43:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0) 
 44:  
 45:     IF (.NOT.GTEST) RETURN 
 46:  
 47:     G = -24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2*SIG6*R6 
 48:  
 49:     CALL PAIRWISE_GRAD(X1, X2, G, PG) 
 50:  
 51:     IF (.NOT.STEST) RETURN 
 52:  
 53:     R8 = R2**4 
 54:     R14 = R8*R8/R2 
 55:     F = 672.0D0*SIG12*R14-192.0D0*SIG6*R8 
 56:  
 57:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS) 
 58:  
 59:     RETURN 
 60:  
 61: END SUBROUTINE PAIRWISE_LJ 
 62:  
 63:  
 64: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 65: ! WCA potential, which is simply the LJ potential but truncated at the first minimum of the function and shifted so that 
 66: ! the potential is entirely repulsive. 
 67: ! The energy and gradient are both continuous at the cutoff. 
 68: ! Atom radius sigma is given as a variable (SIG), but will often be set to 1. 
 69: SUBROUTINE PAIRWISE_WCA(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST) 
 70:     IMPLICIT NONE 
 71:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10) ! Maximum number of params is hardcoded here 
 72:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY 
 73:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6) 
 74:     LOGICAL, INTENT(IN)           :: GTEST, STEST  ! GTEST is true when calculating the gradient as well as energy. 
 75:                                                    ! STEST is true when calculating the Hessian as well as the other two. 
 76:  
 77:     DOUBLE PRECISION :: R2, R6, R8, R14, SIG, SIG6, SIG12  ! Various powers of the distance between the atoms, and the atom radius 
 78:     DOUBLE PRECISION :: G, F  ! G tensor and F tensor (see PAIRWISE_GRAD and PAIRWISE_HESSIAN, below) 
 79:  
 80:     SIG = PARAMS(1) 
 81:     SIG6 = SIG**6 
 82:     SIG12 = SIG6**2 
 83:  
 84:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2 
 85:  
 86:     IF(R2.GT.WCA_CUT*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 
 87:         PAIR_ENERGY = 0.0D0 
 88:         PG(:) = 0.0D0 
 89:         P_HESS(:,:) = 0.0D0 
 90:         RETURN 
 91:     ENDIF 
 92:  
 93:     R2 = 1.0D0/R2 
 94:     R6 = R2**3 
 95:  
 96:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0) + 1 
 97:  
 98:     IF (.NOT.GTEST) RETURN 
 99:  
100:     G = -24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2*SIG6*R6 
101:  
102:     CALL PAIRWISE_GRAD(X1, X2, G, PG) 
103:  
104:     IF (.NOT.STEST) RETURN 
105:  
106:     R8 = R2**4 
107:     R14 = R8*R8/R2 
108:     F = 672.0D0*SIG12*R14-192.0D0*SIG6*R8 
109:  
110:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS) 
111:  
112:     RETURN 
113:  
114: END SUBROUTINE PAIRWISE_WCA 
115:  
116: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
117: ! Hooke's Law spring potential 
118: ! The equilibrium bond length is given as PARAMS(1) (saved as R0). Energy is given in units of the spring constant. 
119: SUBROUTINE HARMONIC_SPRINGS(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST) 
120:     IMPLICIT NONE 
121:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10) 
122:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY 
123:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6) 
124:     LOGICAL, INTENT(IN)           :: GTEST, STEST  ! GTEST is true when calculating the gradient as well as energy. 
125:                                                    ! STEST is true when calculating the Hessian as well as the other two. 
126:     DOUBLE PRECISION :: R0, R, R2, G, F 
127:  
128:     R0 = PARAMS(1) 
129:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2 
130:     R = SQRT(R2) 
131:  
132:     PAIR_ENERGY = 0.5D0*(R-R0)**2 
133:  
134:     IF(.NOT.GTEST) RETURN 
135:  
136:     F = R0/R 
137:     G = 1.0D0 - F 
138:  
139:     CALL PAIRWISE_GRAD(X1,X2,G,PG) 
140:  
141:         IF (.NOT.STEST) RETURN 
142:  
143:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS) 
144:  
145:     RETURN 
146:  
147: END SUBROUTINE HARMONIC_SPRINGS 
148: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
149: ! END OF PAIRWISE POTENTIALS 
150: ! Next, two functions that are general to all isotropic pairwise potentials. To be clear, by "isotropic", I mean 
151: ! "depends only on the distance between the two particles" 
152: ! Another note: I refer to the "G and F tensors", to be in line with comments in other source files. In the 
153: ! present case, where there are only 2 atoms, the tensor is represented by a 2x2 matrix in which two elements are 0 and 
154: ! the other two are equal. So I haven't bothered to represent it as a real tensor, but simply as a single scalar for the 
155: ! unique value in the matrix. 
156: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
157: ! For all isotropic potentials the gradient is calculated in the same way. We simply need to know the positions of the 
158: ! two particles, and the G tensor: G = (1/r)dU/dr for an isotropic potential U(r) 
159: SUBROUTINE PAIRWISE_GRAD(X1, X2, G, PG) 
160:     IMPLICIT NONE 
161:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), G 
162:     DOUBLE PRECISION, INTENT(OUT) :: PG(6) 
163:  
164:     ! Gradient with respect to atom 1 coordinates 
165:     PG(1) = G*(X1(1)-X2(1)) 
166:     PG(2) = G*(X1(2)-X2(2)) 
167:     PG(3) = G*(X1(3)-X2(3)) 
168:     ! Gradient with respect to atom 2 coordinates 
169:     PG(4) = G*(X2(1)-X1(1)) 
170:     PG(5) = G*(X2(2)-X1(2)) 
171:     PG(6) = G*(X2(3)-X1(3)) 
172:  
173:     RETURN 
174: END SUBROUTINE PAIRWISE_GRAD 
175:  
176: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
177: ! For all isotropic potentials the Hessian is calculated in the same way. We simply need to know the positions of the 
178: ! two particles, the G tensor (see above) and the F tensor: F = r*d[(1/r)dU/dr]/dr 
179: SUBROUTINE PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS) 
180:     IMPLICIT NONE 
181:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), G, F, R2 
182:     DOUBLE PRECISION, INTENT(OUT) :: P_HESS(6,6) 
183:     INTEGER :: J1,J2 
184:  
185: !       Compute the hessian. First are the entirely diagonal terms (6 of them) 
186: !       These correspond to 2nd derivatives wrt the same coordinate twice 
187:     DO J1=1,3 
188:         P_HESS(J1,J1) = F*R2*(X1(J1)-X2(J1))**2 + G 
189:         P_HESS(3+J1,3+J1) = F*R2*(X2(J1)-X1(J1))**2 + G 
190:     ENDDO 
191:  
192: !       Next the terms where x_i and x_j are on the same atom but different cartesian directions 
193: !       (12 of these but we only calculate the upper RH block, then symmetrise later) 
194:     DO J1=1,3 
195:         DO J2=J1+1,3 
196:             P_HESS(J1,J2) = F*R2*(X1(J1)-X2(J1))*(X1(J2)-X2(J2)) 
197:             P_HESS(3+J1,3+J2) = F*R2*(X2(J1)-X1(J1))*(X2(J2)-X1(J2)) 
198:         ENDDO 
199:     ENDDO 
200:  
201: !       Different atoms, same cartesian direction (6 of these but we only calculate the upper RH block, then symmetrise later) 
202:     DO J1=1,3 
203:         P_HESS(J1,3+J1) = -F*R2*(X1(J1)-X2(J1))**2 - G 
204:     ENDDO 
205:  
206: !       Finally, different atoms and different coordinates (12 of these, we only calculate 6) 
207:     DO J1=1,3 
208:         DO J2=J1+1,3 
209:             P_HESS(J1,3+J2) = -F*R2*(X1(J1)-X2(J1))*(X1(J2)-X2(J2)) 
210:             P_HESS(J2,3+J1) = P_HESS(J1,3+J2) 
211:         ENDDO 
212:     ENDDO 
213:  
214: !     Symmetrise Hessian 
215:     P_HESS(2,1)=P_HESS(1,2) 
216:     P_HESS(3,1)=P_HESS(1,3) 
217:     P_HESS(3,2)=P_HESS(2,3) 
218:     P_HESS(4,1)=P_HESS(1,4) 
219:     P_HESS(4,2)=P_HESS(2,4) 
220:     P_HESS(4,3)=P_HESS(3,4) 
221:     P_HESS(5,1)=P_HESS(1,5) 
222:     P_HESS(5,2)=P_HESS(2,5) 
223:     P_HESS(5,3)=P_HESS(3,5) 
224:     P_HESS(5,4)=P_HESS(4,5) 
225:     P_HESS(6,1)=P_HESS(1,6) 
226:     P_HESS(6,2)=P_HESS(2,6) 
227:     P_HESS(6,3)=P_HESS(3,6) 
228:     P_HESS(6,4)=P_HESS(4,6) 
229:     P_HESS(6,5)=P_HESS(5,6) 
230:  
231: !    DO J1=1,6 
232: !       DO J2=J1+1,6 
233: !          P_HESS(J2,J1)=P_HESS(J1,J2) 
234: !       ENDDO 
235: !    ENDDO 
236:  
237:     RETURN 
238: END SUBROUTINE PAIRWISE_HESSIAN 
239:  
240: END MODULE PAIRWISE_POTENTIALS 


r29883/potential.f90 2016-03-16 18:33:23.410956460 +0000 r29882/potential.f90 2016-03-16 18:33:24.602968717 +0000
 20: SUBROUTINE POTENTIAL(X, GRAD, EREAL, GRADT, SECT) 20: SUBROUTINE POTENTIAL(X, GRAD, EREAL, GRADT, SECT)
 21: ! This subroutine acts as a wrapper for all of the potentials implemented in GMIN. 21: ! This subroutine acts as a wrapper for all of the potentials implemented in GMIN.
 22: ! The potential is chosen in the 'data' file (default is Lennard-Jones particles). 22: ! The potential is chosen in the 'data' file (default is Lennard-Jones particles).
 23: !  23: ! 
 24: ! In addition, this routine also applies fields, if appropriate.  24: ! In addition, this routine also applies fields, if appropriate. 
 25:    USE COMMONS 25:    USE COMMONS
 26:    USE GENRIGID, ONLY: RIGIDINIT, NRIGIDBODY, DEGFREEDOMS, ATOMRIGIDCOORDT, & 26:    USE GENRIGID, ONLY: RIGIDINIT, NRIGIDBODY, DEGFREEDOMS, ATOMRIGIDCOORDT, &
 27:                        AACONVERGENCET, AACONVERGENCE, TRANSFORMRIGIDTOC, & 27:                        AACONVERGENCET, AACONVERGENCE, TRANSFORMRIGIDTOC, &
 28:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, & 28:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, &
 29:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS 29:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS
 30:    USE MULTIPOT, ONLY: MULTIPOT_CALL 
 31:    USE QMODULE, ONLY: QMIN, QMINP  30:    USE QMODULE, ONLY: QMIN, QMINP 
 32:    USE PERMU, ONLY: FIN 31:    USE PERMU, ONLY: FIN
 33:    USE PORFUNCS 32:    USE PORFUNCS
 34:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C 33:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C
 35:    USE MODHESS, ONLY: HESS, RBAANORMALMODET  34:    USE MODHESS, ONLY: HESS, RBAANORMALMODET 
 36:    USE MODAMBER9, ONLY: ATMASS1 35:    USE MODAMBER9, ONLY: ATMASS1
 37:    USE POLIRMOD, ONLY: POLIR 36:    USE POLIRMOD, ONLY: POLIR
 38:    USE MBPOLMOD, ONLY: MBPOL 37:    USE MBPOLMOD, ONLY: MBPOL
 39:    USE SWMOD, ONLY: SWTYPE 38:    USE SWMOD, ONLY: SWTYPE
 40:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS 39:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS
117: !      WRITE(88,'(G20.10)') (OSTEP(J1), J1=1, NPAR)116: !      WRITE(88,'(G20.10)') (OSTEP(J1), J1=1, NPAR)
118: !      CALL SYSTEM('rm ssave')117: !      CALL SYSTEM('rm ssave')
119: !      STOP118: !      STOP
120: !   END IF119: !   END IF
121: 120: 
122: ! Count the number of potential calls.121: ! Count the number of potential calls.
123:    NPCALL=NPCALL+1122:    NPCALL=NPCALL+1
124: 123: 
125: 10 CONTINUE124: 10 CONTINUE
126: 125: 
127:    IF (MULTIPOTT) THEN126:    IF (MSORIGT) THEN
128:        IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN 
129:            XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS) 
130:            CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)                
131:        END IF 
132:  
133: !       CALL RAD(X, GRAD, EREAL, GRADT) 
134:        CALL MULTIPOT_CALL(X, GRAD, EREAL, GRADT, SECT) 
135:  
136:        IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN 
137:  
138:            X(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
139:            X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS) 
140:  
141:            IF(GRADT) THEN 
142:               CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD) 
143:               GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
144:               GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS) 
145:            ENDIF 
146:  
147:            IF(SECT) THEN 
148:               CALL TRANSFORMHESSIAN(HESS, GRAD, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET) 
149:               HESS(DEGFREEDOMS+1:3*NATOMS,:)=0.0D0 
150:               HESS(:, DEGFREEDOMS+1:3*NATOMS)=0.0D0 
151:               HESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)=XRIGIDHESS(1:DEGFREEDOMS, 1:DEGFREEDOMS) 
152:            ENDIF 
153:  
154:        END IF 
155:  
156:    ELSE IF (MSORIGT) THEN 
157:       CALL RAD(X, GRAD, EREAL, GRADT)127:       CALL RAD(X, GRAD, EREAL, GRADT)
158:       IF (CUTT) THEN128:       IF (CUTT) THEN
159:          CALL MSORIGC(NATOMS, X, GRAD, EREAL, GRADT)129:          CALL MSORIGC(NATOMS, X, GRAD, EREAL, GRADT)
160:       ELSE130:       ELSE
161:          CALL MSORIG(NATOMS, X, GRAD, EREAL, GRADT)131:          CALL MSORIG(NATOMS, X, GRAD, EREAL, GRADT)
162:       END IF132:       END IF
163:       IF (FTEST) THEN133:       IF (FTEST) THEN
164:          RETURN134:          RETURN
165:       END IF135:       END IF
166:    ELSE IF (MSTRANST) THEN136:    ELSE IF (MSTRANST) THEN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0