hdiff output

r33313/centre.f 2017-09-15 17:30:11.746064290 +0100 r33312/centre.f 2017-09-15 17:30:12.430073343 +0100
 19: C 19: C
 20: C************************************************************************** 20: C**************************************************************************
 21: C 21: C
 22: C  Subroutine CENTRE moves the centre of mass to the origin. 22: C  Subroutine CENTRE moves the centre of mass to the origin.
 23: C 23: C
 24: C********************************************************************* 24: C*********************************************************************
 25: C 25: C
 26:       SUBROUTINE CENTRE2(X) 26:       SUBROUTINE CENTRE2(X)
 27:       USE commons 27:       USE commons
 28:       use genrigid 28:       use genrigid
 29:       USE OPP_MOD, ONLY: OPP_MODE 
 30:       IMPLICIT NONE 29:       IMPLICIT NONE
 31:       DOUBLE PRECISION XMASS, YMASS, ZMASS, DIST, DISTMAX, X(3*NATOMS) 30:       DOUBLE PRECISION XMASS, YMASS, ZMASS, DIST, DISTMAX, X(3*NATOMS)
 32: !hk286 - additional variables for generalised rigid body 31: !hk286 - additional variables for generalised rigid body
 33:       DOUBLE PRECISION XRIGIDCOORDS(DEGFREEDOMS) 32:       DOUBLE PRECISION XRIGIDCOORDS(DEGFREEDOMS)
 34:       INTEGER I, J1 33:       INTEGER I, J1
 35: !hk286 34: !hk286
 36:  
 37: ! jwrm2> If using OPP with alchemical moves, the last coordinates are the 
 38: ! potential parameters. Fiddle NATOMS so these aren't included. 
 39:       IF (OPPT .AND. OPP_MODE .EQ. 3) THEN 
 40:          NATOMS = NATOMS - 1 
 41:       END IF 
 42:           35:          
 43: ! hk286 > if generalised rigid body is used, be careful when averaging 36: ! hk286 > if generalised rigid body is used, be careful when averaging
 44: ! hk286 > no need to shift the rotational degrees of freedom! 37: ! hk286 > no need to shift the rotational degrees of freedom!
 45:       IF ( .NOT. ATOMRIGIDCOORDT ) THEN 38:       IF ( .NOT. ATOMRIGIDCOORDT ) THEN
 46:          XMASS=0.0D0 39:          XMASS=0.0D0
 47:          YMASS=0.0D0 40:          YMASS=0.0D0
 48:          ZMASS=0.0D0 41:          ZMASS=0.0D0
 49:          DO I=1,NRIGIDBODY 42:          DO I=1,NRIGIDBODY
 50:             XMASS=XMASS+X(3*(I-1)+1) 43:             XMASS=XMASS+X(3*(I-1)+1)
 51:             YMASS=YMASS+X(3*(I-1)+2) 44:             YMASS=YMASS+X(3*(I-1)+2)
136: C        DO J1=1,NATOMS129: C        DO J1=1,NATOMS
137: C           X(3*J1-2)=X(3*J1-2)*DISTMAX130: C           X(3*J1-2)=X(3*J1-2)*DISTMAX
138: C           X(3*J1-1)=X(3*J1-1)*DISTMAX131: C           X(3*J1-1)=X(3*J1-1)*DISTMAX
139: C           X(3*J1)  =X(3*J1)  *DISTMAX132: C           X(3*J1)  =X(3*J1)  *DISTMAX
140: C        ENDDO133: C        ENDDO
141: C     ENDIF134: C     ENDIF
142: 135: 
143:       ENDIF136:       ENDIF
144: !     hk286 - endif of generalised rigid body if statement137: !     hk286 - endif of generalised rigid body if statement
145: 138: 
146: ! jwrm2> Reset NATOMS 
147:       IF (OPPT .AND. OPP_MODE .EQ. 3) THEN 
148:          NATOMS = NATOMS + 1 
149:       END IF 
150:  
151:       RETURN139:       RETURN
152:       END SUBROUTINE CENTRE2140:       END
153: 141: 
154:       SUBROUTINE SETCENTRE(X)142:       SUBROUTINE SETCENTRE(X)
155:       USE commons143:       USE commons
156:       use genrigid144:       use genrigid
157: 145: 
158:       IMPLICIT NONE146:       IMPLICIT NONE
159:       DOUBLE PRECISION XMASS, YMASS, ZMASS, X(3*NATOMS)147:       DOUBLE PRECISION XMASS, YMASS, ZMASS, X(3*NATOMS)
160: !hk286 - additional variables for generalised rigid body148: !hk286 - additional variables for generalised rigid body
161:       DOUBLE PRECISION XRIGIDCOORDS(DEGFREEDOMS)149:       DOUBLE PRECISION XRIGIDCOORDS(DEGFREEDOMS)
162:       INTEGER I, J1150:       INTEGER I, J1
234:             X(3*(I-1)+2)=X(3*(I-1)+2)-YMASS+CENTY222:             X(3*(I-1)+2)=X(3*(I-1)+2)-YMASS+CENTY
235:             X(3*(I-1)+3)=X(3*(I-1)+3)-ZMASS+CENTZ223:             X(3*(I-1)+3)=X(3*(I-1)+3)-ZMASS+CENTZ
236:          ENDDO224:          ENDDO
237:       ENDIF225:       ENDIF
238:       IF (DEBUG) WRITE(MYUNIT,'(A,3F12.4)') 'centre of mass moved to ',CENTX,CENTY,CENTZ226:       IF (DEBUG) WRITE(MYUNIT,'(A,3F12.4)') 'centre of mass moved to ',CENTX,CENTY,CENTZ
239: 227: 
240:       ENDIF228:       ENDIF
241: !     hk286 - endif of generalised rigid body if statement229: !     hk286 - endif of generalised rigid body if statement
242: 230: 
243:       RETURN231:       RETURN
244:       END SUBROUTINE SETCENTRE232:       END
245: 233: 


r33313/keywords.f 2017-09-15 17:30:11.970067251 +0100 r33312/keywords.f 2017-09-15 17:30:12.662076426 +0100
 38:       USE PORFUNCS 38:       USE PORFUNCS
 39:       USE MYGA_PARAMS 39:       USE MYGA_PARAMS
 40:       USE BGUPMOD 40:       USE BGUPMOD
 41:       USE GLJYMOD 41:       USE GLJYMOD
 42:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 42:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 43:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP 43:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP
 44:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, 44:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS,
 45:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS, 45:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS,
 46:      &                        LJ_GAUSS_INITIALISE 46:      &                        LJ_GAUSS_INITIALISE
 47:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS, 47:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS,
 48:      &                   OPP_INITIALISE, OPP_K_LOWER, OPP_K_UPPER 48:      &                   OPP_INITIALISE
 49:       USE MBPOLMOD, ONLY: MBPOLINIT 49:       USE MBPOLMOD, ONLY: MBPOLINIT
 50:       USE SWMOD, ONLY: SWINIT, MWINIT 50:       USE SWMOD, ONLY: SWINIT, MWINIT
 51:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS 51:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS
 52: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP, 52: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP,
 53: !     &                                 AMBER12_RESIDUES, 53: !     &                                 AMBER12_RESIDUES,
 54: !     &                                 POPULATE_ATOM_DATA 54: !     &                                 POPULATE_ATOM_DATA
 55:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 55:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 56:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 56:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 57:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 57:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,
 58:      &     PARSE_MLJ_PARAMS 58:      &     PARSE_MLJ_PARAMS
6576:           CALL READF(LJ_GAUSS_SIGMASQ)6576:           CALL READF(LJ_GAUSS_SIGMASQ)
6577:           IF (LJ_GAUSS_MODE .EQ. 3) CALL READI(LJ_GAUSS_PARAMS)6577:           IF (LJ_GAUSS_MODE .EQ. 3) CALL READI(LJ_GAUSS_PARAMS)
6578:           CALL LJ_GAUSS_INITIALISE()6578:           CALL LJ_GAUSS_INITIALISE()
6579: 6579: 
6580:       ELSE IF (WORD.EQ.'OPP') THEN6580:       ELSE IF (WORD.EQ.'OPP') THEN
6581:           OPPT = .TRUE.6581:           OPPT = .TRUE.
6582:           CALL READI(OPP_MODE)6582:           CALL READI(OPP_MODE)
6583:           CALL READF(OPP_RCUT)6583:           CALL READF(OPP_RCUT)
6584:           CALL READF(OPP_K)6584:           CALL READF(OPP_K)
6585:           CALL READF(OPP_PHI)6585:           CALL READF(OPP_PHI)
6586:           IF (OPP_MODE .EQ. 2 .OR. OPP_MODE .EQ. 3) THEN6586:           IF (OPP_MODE .EQ. 2) CALL READI(OPP_PARAMS)
6587:               CALL READI(OPP_PARAMS) 
6588:               IF (BTEST(OPP_PARAMS, 0)) THEN 
6589:                 CALL READF(OPP_K_LOWER) 
6590:                 CALL READF(OPP_K_UPPER) 
6591:               END IF 
6592:           END IF 
6593:           CALL OPP_INITIALISE()6587:           CALL OPP_INITIALISE()
6594: ! jwrm2> end6588: ! jwrm2> end
6595: 6589: 
6596:       ELSE IF (WORD.EQ.'STOCK') THEN6590:       ELSE IF (WORD.EQ.'STOCK') THEN
6597:          STOCKT=.TRUE.6591:          STOCKT=.TRUE.
6598:          RIGID=.TRUE.6592:          RIGID=.TRUE.
6599:          NRBSITES=16593:          NRBSITES=1
6600:          CALL READF(STOCKMU)6594:          CALL READF(STOCKMU)
6601:          CALL READF(STOCKLAMBDA)6595:          CALL READF(STOCKLAMBDA)
6602:          ALLOCATE(SITE(NRBSITES,3))6596:          ALLOCATE(SITE(NRBSITES,3))


r33313/opp.f90 2017-09-15 17:30:12.190070168 +0100 r33312/opp.f90 2017-09-15 17:30:12.886079393 +0100
 21: ! an oscillating pair potential (OPP). The possibility of varying 21: ! an oscillating pair potential (OPP). The possibility of varying
 22: ! and optimising the potential parameters is included. 22: ! and optimising the potential parameters is included.
 23: ! Also included is a framework for a varying triclinic unit cell. This could be 23: ! Also included is a framework for a varying triclinic unit cell. This could be
 24: ! separated out and generalised to all isotropic pair potentials. 24: ! separated out and generalised to all isotropic pair potentials.
 25: ! Questions to jwrm2. 25: ! Questions to jwrm2.
 26: ! 26: !
 27: MODULE OPP_MOD 27: MODULE OPP_MOD
 28:  28: 
 29: ! Public parameters 29: ! Public parameters
 30: PUBLIC :: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS 30: PUBLIC :: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS
 31: PUBLIC :: OPP_K_LOWER, OPP_K_UPPER 
 32: ! Public functions 31: ! Public functions
 33: PUBLIC :: OPP, OPP_INITIALISE, VIEW_OPP, OPP_TAKESTEP 32: PUBLIC :: OPP, OPP_INITIALISE, VIEW_OPP, OPP_TAKESTEP
 34: ! Private parameters 33: ! Private parameters
 35: PRIVATE :: V_EPS, V_SIGMA, FD_RCUT, IMAGE_CUTOFF!, SD_RCUT 34: PRIVATE :: V_EPS, V_SIGMA, FD_RCUT, IMAGE_CUTOFF!, SD_RCUT
 36: PRIVATE :: MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_K_STEP, MAX_PHI_STEP 35: PRIVATE :: MAX_LENGTH_STEP, MAX_ANGLE_STEP, MAX_K_STEP, MAX_PHI_STEP
 37: PRIVATE :: K_REPULSION 36: PRIVATE :: K_REPULSION, K_LOWER, K_UPPER
 38: ! Private functions 37: ! Private functions
 39: PRIVATE :: PERIODIC_RESET, OPP_PER, OPP_PER_TRI, CALC_ENERGY 38: PRIVATE :: PERIODIC_RESET, OPP_PER, OPP_PER_TRI, CALC_ENERGY
 40: PRIVATE :: CALC_GRAD, CALC_SEC, CALC_FCTS, CONSTRAIN_PARAMETERS 39: PRIVATE :: CALC_GRAD, CALC_SEC, CALC_FCTS, CONSTRAIN_PARAMETERS
 41: PRIVATE :: CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME 40: PRIVATE :: CHECK_ANGLES, REJECT, CONSTRAIN_VOLUME
 42:  41: 
 43: ! K: frequency parameter, controls the number of wells 42: ! K: frequency parameter, controls the number of wells
 44: ! PHI: phase parameter, controls the position of the first minimum 43: ! PHI: phase parameter, controls the position of the first minimum
 45: ! RCUT: cutoff distance at which the energy, gradients and Hessian of 44: ! RCUT: cutoff distance at which the energy, gradients and Hessian of
 46: !       the potential go smoothly to zero. 45: !       the potential go smoothly to zero.
 47: DOUBLE PRECISION :: OPP_K, OPP_PHI, OPP_RCUT 46: DOUBLE PRECISION :: OPP_K, OPP_PHI, OPP_RCUT
 48:  47: 
 49: ! OPP_MODE: 0 standard minimisation 48: ! OPP_MODE: 0 standard minimisation
 50: !           1 the final triplet of coordinates are the unit cell length 49: !           1 the final triplet of coordinates are the unit cell length
 51: !             parameters and the second last triplet are the unit cell 50: !             parameters and the second last triplet are the unit cell
 52: !             angles, giving a triclinic unit cell, which will be optimised 51: !             angles, giving a triclinic unit cell, which will be optimised
 53: !           2 as for 1, but the third last triplet are the potential 52: !           2 as for 1, but the third last triplet are the potential
 54: !             parameters OPP_K, OPP_PHI and 0 and 53: !             parameters OPP_K, OPP_PHI and 0 and
 55: !             the potential parameters will be optimised, depending on the 54: !             the potential parameters will be optimised, depending on the
 56: !             value of OPP_PARAMS 55: !             value of OPP_PARAMS
 57: !           3 for clusters, but with the last set of coordinates as the 
 58: !             potential parameters 
 59: ! OPP_PARAMS: Specifies which of the three potential parameters will be 56: ! OPP_PARAMS: Specifies which of the three potential parameters will be
 60: !             optimised. An integer between 1 and 3. 57: !             optimised. An integer between 1 and 3.
 61: !           1 OPP_K will be optimised 58: !           1 OPP_K will be optimised
 62: !           2 OPP_PHI will be optimised 59: !           2 OPP_PHI will be optimised
 63: !           3 both parameters will be optimised 60: !           3 both parameters will be optimised
 64: INTEGER :: OPP_MODE, OPP_PARAMS 61: INTEGER :: OPP_MODE, OPP_PARAMS
 65:  62: 
 66: ! V_RCUT: the unmodified potential evaluated at the cutoff 63: ! V_RCUT: the unmodified potential evaluated at the cutoff
 67: ! FD_RCUT: the first derivative evaluated at the cutoff 64: ! FD_RCUT: the first derivative evaluated at the cutoff
 68: ! SD_RCUT: the second derivative evaluated at the cutoff 65: ! SD_RCUT: the second derivative evaluated at the cutoff
 77: DOUBLE PRECISION, PARAMETER :: IMAGE_CUTOFF = 30 74: DOUBLE PRECISION, PARAMETER :: IMAGE_CUTOFF = 30
 78:  75: 
 79: ! Parameters specifying the largest steps in OPP_TAKESTEP 76: ! Parameters specifying the largest steps in OPP_TAKESTEP
 80: ! MAX_LENGTH_STEP: largest allowed step in the length of a lattice 77: ! MAX_LENGTH_STEP: largest allowed step in the length of a lattice
 81: !                  vector 78: !                  vector
 82: ! MAX_ANGLE_STEP: largest allowed step in the angle of a lattice vector 79: ! MAX_ANGLE_STEP: largest allowed step in the angle of a lattice vector
 83: ! MAX_K_STEP: largest allowed step in OPP_K 80: ! MAX_K_STEP: largest allowed step in OPP_K
 84: ! MAX_PHI_STEP: largest allowed step in OPP_PHI 81: ! MAX_PHI_STEP: largest allowed step in OPP_PHI
 85: DOUBLE PRECISION, PARAMETER :: MAX_LENGTH_STEP = 0.3D0 82: DOUBLE PRECISION, PARAMETER :: MAX_LENGTH_STEP = 0.3D0
 86: DOUBLE PRECISION, PARAMETER :: MAX_ANGLE_STEP = 0.1D0 83: DOUBLE PRECISION, PARAMETER :: MAX_ANGLE_STEP = 0.1D0
 87: DOUBLE PRECISION, PARAMETER :: MAX_K_STEP = 3.D0 84: DOUBLE PRECISION, PARAMETER :: MAX_K_STEP = 1.D0
 88: DOUBLE PRECISION, PARAMETER :: MAX_PHI_STEP = 0.5D0 85: DOUBLE PRECISION, PARAMETER :: MAX_PHI_STEP = 0.5D0
 89:  86: 
 90: ! Parameters used for constraining the unit cell volume and the potential 87: ! Parameters used for constraining the unit cell volume and the potential
 91: ! parameters 88: ! parameters
 92: ! V_EPS: scaling factor for the unit cell volume contraint energy 89: ! V_EPS: scaling factor for the unit cell volume contraint energy
 93: ! V_SIGMA: WCA cutoff distance 90: ! V_SIGMA: WCA cutoff distance
 94: ! OPP_K_LOWER, OPP_K_UPPER: bounds on OPP_K 91: ! K_LOWER, K_UPPER: bounds on OPP_K
 95: ! K_REPULSION: harmonic force constant for OPP_K 92: ! K_REPULSION: harmonic force constant for OPP_K
 96: ! OPP_PHI is periodic over the range 0 to 2*PI, so doesn't need constraints 93: ! OPP_PHI is periodic over the range 0 to 2*PI, so doesn't need constraints
 97: DOUBLE PRECISION, PARAMETER :: V_EPS = 1.D-3 94: DOUBLE PRECISION, PARAMETER :: V_EPS = 1.D-3
 98: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1 95: DOUBLE PRECISION, PARAMETER :: V_SIGMA = 3.D-1
 99: DOUBLE PRECISION :: OPP_K_LOWER, OPP_K_UPPER 96: DOUBLE PRECISION, PARAMETER :: K_LOWER = 5.D0, K_UPPER = 20.D0
100: DOUBLE PRECISION, PARAMETER :: K_REPULSION = 1.D4 97: DOUBLE PRECISION, PARAMETER :: K_REPULSION = 1.D4
101:  98: 
102: CONTAINS 99: CONTAINS
103: 100: 
104: !-------------------------------------------------------------------------------101: !-------------------------------------------------------------------------------
105: !102: !
106: ! Main routine for calculating the energy and gradients.103: ! Main routine for calculating the energy and gradients.
107: ! X: coordinates array104: ! X: coordinates array
108: ! GRAD: gradients array105: ! GRAD: gradients array
109: ! ENERGY: calculated energy106: ! ENERGY: calculated energy
137: 134: 
138:         ! RIJ: Interparticle vector135:         ! RIJ: Interparticle vector
139:         ! MODRIJ: distance between a pair of particles136:         ! MODRIJ: distance between a pair of particles
140:         ! DVDR: derivative of pair potential wrt MODRIJ137:         ! DVDR: derivative of pair potential wrt MODRIJ
141:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ138:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ
142:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR, D2VDR2139:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR, D2VDR2
143: 140: 
144:         ! Factors for triclinic lattice141:         ! Factors for triclinic lattice
145:         TYPE(BOX_DERIV) :: BD142:         TYPE(BOX_DERIV) :: BD
146: 143: 
147: !        WRITE (MYUNIT, *) 'OPP> beginning coords' 
148: !        DO I = 1, NATOMS 
149: !            WRITE (MYUNIT, *) X(3*I-2:3*I) 
150: !        END DO 
151:  
152:         ! Initialise output variables144:         ! Initialise output variables
153:         ENERGY = 0.D0145:         ENERGY = 0.D0
154:         GRAD(:) = 0.D0146:         GRAD(:) = 0.D0
155:         IF (STEST) HESS(:,:) = 0.D0147:         IF (STEST) HESS(:,:) = 0.D0
156: 148: 
 149:         ! Calculate factors that do depend on parameters, but not on RIJ
 150:         IF (OPP_MODE .EQ. 2) THEN
 151:             CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-6))
 152:         ELSE
 153:             CALL CALC_FCTS()
 154:         END IF
 155: 
157:         IF (.NOT. PERIODIC) THEN156:         IF (.NOT. PERIODIC) THEN
158:             IF (OPP_MODE .EQ. 0) THEN157:             IF (OPP_MODE .EQ. 0) THEN
159:                 ! Normal cluster158:                 ! Normal cluster
160:                 CALL CALC_FCTS() 
161:                 NMOL = NATOMS159:                 NMOL = NATOMS
162:                 DO I = 1, NMOL-1 ! Outer loop over particles160:                 DO I = 1, NMOL ! Outer loop over particles
163:                     DO J = I+1, NMOL ! Inner loop over particles161:                     DO J = 1, NMOL ! Inner loop over particles
164:  
165:                         ! Get the particle-particle distance 
166:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J) 
167:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
168:  
169:                         ! Check against cutoff 
170:                         IF (MODRIJ .LT. OPP_RCUT) THEN 
171:                             ! Add energy and gradients 
172:                             ENERGY = ENERGY + CALC_ENERGY(MODRIJ) 
173: 162: 
174:                             IF (GTEST) THEN163:                         ! No self interaction
175:                                 DVDR = CALC_GRAD(MODRIJ) / MODRIJ164:                         IF (I .EQ. J) CYCLE
176:                                 GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + DVDR*RIJ(:) 
177:                                 GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - DVDR*RIJ(:) 
178:                             END IF ! GTEST 
179:  
180:                             IF (STEST) THEN 
181:                                 CALL HESS_CONTRIB(RIJ(:), I, J) 
182:                             END IF ! End second derivative 
183:                         END IF ! End less than cutoff 
184:                     END DO ! Inner loop over particles 
185:                 END DO ! Outer loop over particles 
186:             ELSE IF (OPP_MODE .EQ. 3) THEN 
187:                 ! Cluster with potential parameter minimisation 
188:                 CALL CALC_FCTS(X(NATOMS*3-2:NATOMS*3)) 
189:                 ! Reset phi parameter 
190:                 CALL PERIODIC_RESET(X(:)) 
191:                 NMOL = NATOMS - 1 
192:                 DO I = 1, NMOL-1 ! Outer loop over particles 
193:                     DO J = I+1, NMOL ! Inner loop over particles 
194: 165: 
195:                         ! Get the particle-particle distance166:                         ! Get the particle-particle distance
196:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)167:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)
197:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))168:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
198: 169: 
199:                         ! Check against cutoff170:                         ! Add energy and gradients
200:                         IF (MODRIJ .LT. OPP_RCUT) THEN171:                         IF (MODRIJ .LT. OPP_RCUT) THEN
201:                             ! Add energy and gradients172: 
202:                             ENERGY = ENERGY + CALC_ENERGY(MODRIJ)173:                             ! Don't double count energy
 174:                             IF (J .GT. I) THEN
 175:                                 ENERGY = ENERGY + CALC_ENERGY(MODRIJ)
 176:                             END IF
203: 177: 
204:                             IF (GTEST) THEN178:                             IF (GTEST) THEN
205:                                 DVDR = CALC_GRAD(MODRIJ) / MODRIJ179:                                 DVDR = CALC_GRAD(MODRIJ) / MODRIJ
206:                                 GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + DVDR*RIJ(:)180:                                 GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + DVDR*RIJ(:)
207:                                 GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - DVDR*RIJ(:)181: !                                GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - DVDR*RIJ(:)
208:                                 GRAD(3*NATOMS-2:3*NATOMS) = & 
209:                                     GRAD(3*NATOMS-2:3*NATOMS) + & 
210:                                     CALC_GRAD_PARAMS(MODRIJ) 
211:                             END IF ! GTEST182:                             END IF ! GTEST
212: 183: 
213:                             IF (STEST) THEN184:                             IF (STEST) THEN
214:                                 CALL HESS_CONTRIB(RIJ(:), I, J)185:                                 ! Hessian terms
 186:                                 DVDR = DVDR / MODRIJ
 187:                                 D2VDR2 = (CALC_SEC(MODRIJ) - DVDR) / MODRIJ**2
 188:                                 ! Same atom, same coordinate
 189:                                 HESS(3*I-2,3*I-2) = HESS(3*I-2,3*I-2) + &
 190:                                     D2VDR2 * RIJ(1)**2 + DVDR
 191:                                 HESS(3*I-1,3*I-1) = HESS(3*I-1,3*I-1) + &
 192:                                     D2VDR2 * RIJ(2)**2 + DVDR
 193:                                 HESS(3*I  ,3*I  ) = HESS(3*I  ,3*I  ) + &
 194:                                     D2VDR2 * RIJ(3)**2 + DVDR
 195:                                 ! Same atom, different coordinate
 196:                                 HESS(3*I-2,3*I-1) = HESS(3*I-2,3*I-1) + &
 197:                                     D2VDR2 * RIJ(1) * RIJ(2)
 198:                                 HESS(3*I-1,3*I  ) = HESS(3*I-1,3*I  ) + &
 199:                                     D2VDR2 * RIJ(2) * RIJ(3)
 200:                                 HESS(3*I  ,3*I-2) = HESS(3*I  ,3*I-2) + &
 201:                                     D2VDR2 * RIJ(3) * RIJ(1)
 202:                                 HESS(3*I-1,3*I-2) = HESS(3*I-1,3*I-2) + &
 203:                                     D2VDR2 * RIJ(2) * RIJ(1)
 204:                                 HESS(3*I  ,3*I-1) = HESS(3*I  ,3*I-1) + &
 205:                                     D2VDR2 * RIJ(3) * RIJ(2)
 206:                                 HESS(3*I-2,3*I  ) = HESS(3*I-2,3*I) + &
 207:                                     D2VDR2 * RIJ(1) * RIJ(3)
 208:                                 ! Different atom, same coordinate
 209:                                 HESS(3*I-2,3*J-2) = - D2VDR2 * RIJ(1)**2 - DVDR
 210:                                 HESS(3*I-1,3*J-1) = - D2VDR2 * RIJ(2)**2 - DVDR
 211:                                 HESS(3*I  ,3*J  ) = - D2VDR2 * RIJ(3)**2 - DVDR
 212:                                 HESS(3*J-2,3*I-2) = - D2VDR2 * RIJ(1)**2 - DVDR
 213:                                 HESS(3*J-1,3*I-1) = - D2VDR2 * RIJ(2)**2 - DVDR
 214:                                 HESS(3*J  ,3*I  ) = - D2VDR2 * RIJ(3)**2 - DVDR
 215:                                 ! Different atom, different coordinate
 216:                                 HESS(3*I-2,3*J-1) = - D2VDR2 * RIJ(1) * RIJ(2)
 217:                                 HESS(3*I-1,3*J  ) = - D2VDR2 * RIJ(2) * RIJ(3)
 218:                                 HESS(3*I-2,3*J  ) = - D2VDR2 * RIJ(3) * RIJ(1)
 219:                                 HESS(3*J-1,3*I-2) = - D2VDR2 * RIJ(1) * RIJ(2)
 220:                                 HESS(3*J  ,3*I-1) = - D2VDR2 * RIJ(2) * RIJ(3)
 221:                                 HESS(3*J  ,3*I-2) = - D2VDR2 * RIJ(3) * RIJ(1)
215:                             END IF ! End second derivative222:                             END IF ! End second derivative
216:                         END IF ! End less than cutoff223:                         END IF ! If less than cutoff distance
217:                     END DO ! Inner loop over particles224:                     END DO ! Inner loop over particles
218:                 END DO ! Outer loop over particles225:                 END DO ! Outer loop over particles
219:  
220:                 ! Constrain the potential parameters with a harmonic repulsion 
221:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-2:3*NATOMS), & 
222:                                           GTEST, STEST) 
223:             ELSE ! mode not equal to zero226:             ELSE ! mode not equal to zero
224:                 WRITE (MYUNIT, *) 'OPP> ERROR, PERIODIC must be ', &227:                 WRITE (MYUNIT, *) 'OPP> ERROR, PERIODIC must be ', &
225:                                   'specified if mode is not 0'228:                                   'specified if mode is not 0'
226:                 STOP229:                 STOP
227:             END IF230:             END IF
228:         ELSE ! PERIODIC231:         ELSE ! PERIODIC
229:             SELECT CASE (OPP_MODE)232:             SELECT CASE (OPP_MODE)
230:             CASE(0)233:             CASE(0)
231:                 ! Periodic, fixed box234:                 ! Periodic, fixed box
232:                 CALL CALC_FCTS() 
233:                 ! Reset all particles to within the box235:                 ! Reset all particles to within the box
234:                 CALL PERIODIC_RESET(X(:))236:                 CALL PERIODIC_RESET(X(:))
235: 237: 
236:                 NMOL = NATOMS238:                 NMOL = NATOMS
237:                 ! We need to include self-interactions of particles in different239:                 ! We need to include self-interactions of particles in different
238:                 ! unit cells240:                 ! unit cells
239:                 DO I = 1, NMOL ! Outer loop over particles241:                 DO I = 1, NMOL ! Outer loop over particles
240:                     DO J = 1, NMOL ! Inner loop over particles242:                     DO J = 1, NMOL ! Inner loop over particles
241:                         ! Add energy and gradients243:                         ! Add energy and gradients
242:                         CALL OPP_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &244:                         CALL OPP_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &
243:                                      GRAD(3*I-2:3*I), ENERGY, GTEST, STEST)245:                                      GRAD(3*I-2:3*I), ENERGY, GTEST, STEST)
244:                     END DO ! Inner loop over particles246:                     END DO ! Inner loop over particles
245:                 END DO ! Outer loop over particles247:                 END DO ! Outer loop over particles
246: 248: 
247:             CASE(1)249:             CASE(1)
248:                 ! Triclinic varying lattice250:                 ! Triclinic varying lattice
249:                 CALL CALC_FCTS() 
250:                 ! Reset all particles to within the box251:                 ! Reset all particles to within the box
251:                 CALL PERIODIC_RESET(X(:))252:                 CALL PERIODIC_RESET(X(:))
252: 253: 
253: !                DO I = 1, NATOMS254: !                DO I = 1, NATOMS
254: !                    WRITE (MYUNIT, *) X(3*I-2:3*I)255: !                    WRITE (MYUNIT, *) X(3*I-2:3*I)
255: !               END DO256: !               END DO
256: 257: 
257:                 ! Calculate box derivative factors258:                 ! Calculate box derivative factors
258:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &259:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
259:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)260:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)
286:                     END DO ! Inner loop over particles287:                     END DO ! Inner loop over particles
287:                 END DO ! Outer loop over particles288:                 END DO ! Outer loop over particles
288: 289: 
289:                 ! Constrain the box volume to be greater than 0 with a WCA-style290:                 ! Constrain the box volume to be greater than 0 with a WCA-style
290:                 ! repulsion291:                 ! repulsion
291:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &292:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &
292:                                       GTEST, STEST)293:                                       GTEST, STEST)
293: 294: 
294:             CASE(2)295:             CASE(2)
295:                 ! Triclinic varying lattice, with varying parameters296:                 ! Triclinic varying lattice, with varying parameters
296:                 CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-6)) 
297:                 ! Reset all particles to within the box297:                 ! Reset all particles to within the box
298:                 CALL PERIODIC_RESET(X(:))298:                 CALL PERIODIC_RESET(X(:))
299: 299: 
300:                 ! Calculate box derivative factors300:                 ! Calculate box derivative factors
301:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &301:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
302:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)302:                                     X(3*NATOMS-2:3*NATOMS  ), OPP_RCUT)
303: 303: 
304:                 ! Check whether the unit cell angles are physically possible.304:                 ! Check whether the unit cell angles are physically possible.
305:                 ! If we have a disallowed unit cell, we set the energy to very305:                 ! If we have a disallowed unit cell, we set the energy to very
306:                 ! high and the gradients to very small, so the step is306:                 ! high and the gradients to very small, so the step is
332:                 END DO ! Outer loop over particles332:                 END DO ! Outer loop over particles
333: 333: 
334:                 ! Constrain the box volume to be greater than 0 with a WCA-style334:                 ! Constrain the box volume to be greater than 0 with a WCA-style
335:                 ! repulsion335:                 ! repulsion
336:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &336:                 CALL CONSTRAIN_VOLUME(ENERGY, GRAD(3*NATOMS-5:3*NATOMS-3), BD, &
337:                                       GTEST, STEST)337:                                       GTEST, STEST)
338: 338: 
339:                 ! Constrain the potential parameters with a harmonic repulsion339:                 ! Constrain the potential parameters with a harmonic repulsion
340:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-6), &340:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-6), &
341:                                           GTEST, STEST)341:                                           GTEST, STEST)
342:  
343:                 ! Fiddle to lock in cubic 
344: !                X(3*NATOMS-1) = X(3*NATOMS-2) 
345: !                X(3*NATOMS  ) = X(3*NATOMS-2) 
346: !                GRAD(3*NATOMS-1) = GRAD(3*NATOMS-2) 
347: !                GRAD(3*NATOMS  ) = GRAD(3*NATOMS-2) 
348:             END SELECT ! Mode selections342:             END SELECT ! Mode selections
349:         END IF ! Periodic343:         END IF ! Periodic
350: 344: 
351:     END SUBROUTINE OPP345:     END SUBROUTINE OPP
352: 346: 
353: !-------------------------------------------------------------------------------347: !-------------------------------------------------------------------------------
354: !348: !
355: ! Rejects a step by setting the energy very high and the gradients very low.349: ! Rejects a step by setting the energy very high and the gradients very low.
356: ! This tricks L-BFGS into thinking it's converged, but the high energy of the350: ! This tricks L-BFGS into thinking it's converged, but the high energy of the
357: ! quench will ensure it is rejected.351: ! quench will ensure it is rejected.
386:     SUBROUTINE CALC_FCTS(PARAMS)380:     SUBROUTINE CALC_FCTS(PARAMS)
387: 381: 
388:         IMPLICIT NONE382:         IMPLICIT NONE
389: 383: 
390:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: PARAMS(3)384:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: PARAMS(3)
391: 385: 
392:         ! S, C: factors which occur repeatedly386:         ! S, C: factors which occur repeatedly
393:         DOUBLE PRECISION :: S, C387:         DOUBLE PRECISION :: S, C
394: 388: 
395:         ! Adjust the potential parameters, if the mode requires it389:         ! Adjust the potential parameters, if the mode requires it
396:         IF (PRESENT(PARAMS)) THEN390:         IF (OPP_MODE .EQ. 2 .AND. PRESENT(PARAMS)) THEN
397:             ! For each parameter, if it is varying, set it from the coordinate.391:             ! For each parameter, if it is varying, set it from the coordinate.
398:             ! If it's not varying, reset the coordinate, which TAKESTEP will392:             ! If it's not varying, reset the coordinate, which TAKESTEP will
399:             ! have changed.393:             ! have changed.
400:             IF (BTEST(OPP_PARAMS, 0)) THEN394:             IF (BTEST(OPP_PARAMS, 0)) THEN
401:                 OPP_K = PARAMS(1)395:                 OPP_K = PARAMS(1)
402:             ELSE396:             ELSE
403:                 PARAMS(1) = OPP_K397:                 PARAMS(1) = OPP_K
404:             END IF398:             END IF
405:             IF (BTEST(OPP_PARAMS, 1)) THEN399:             IF (BTEST(OPP_PARAMS, 1)) THEN
406:                 OPP_PHI = PARAMS(2)400:                 OPP_PHI = PARAMS(2)
426:     END SUBROUTINE CALC_FCTS420:     END SUBROUTINE CALC_FCTS
427: 421: 
428: !-------------------------------------------------------------------------------422: !-------------------------------------------------------------------------------
429: !423: !
430: ! Calculates the energy contribution for a pair of particles.424: ! Calculates the energy contribution for a pair of particles.
431: ! MODRIJ: distance beteen the particle pair425: ! MODRIJ: distance beteen the particle pair
432: !426: !
433: !-------------------------------------------------------------------------------427: !-------------------------------------------------------------------------------
434:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)428:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)
435: 429: 
436:         ! PERIODIC: whether peiodic boundary conditions are used 
437:         USE COMMONS, ONLY: PERIODIC 
438:  
439:         IMPLICIT NONE430:         IMPLICIT NONE
440: 431: 
441:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ432:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ
442: 433: 
443:         ! ENERGY: working variable for energy434:         ! ENERGY: working variable for energy
444:         DOUBLE PRECISION :: ENERGY435:         DOUBLE PRECISION :: ENERGY
445: 436: 
446:         ! Unmodified potential437:         ! Unmodified potential
447:         ENERGY = MODRIJ**(-15) + MODRIJ**(-3) * &438:         ENERGY = MODRIJ**(-15) + MODRIJ**(-3) * &
448:                  COS(OPP_K * (MODRIJ - 1) - OPP_PHI)439:                  COS(OPP_K * (MODRIJ - 1) - OPP_PHI)
462:     END FUNCTION CALC_ENERGY453:     END FUNCTION CALC_ENERGY
463: 454: 
464: !-------------------------------------------------------------------------------455: !-------------------------------------------------------------------------------
465: !456: !
466: ! Calculates the pair potential gradient for a pair of particles.457: ! Calculates the pair potential gradient for a pair of particles.
467: ! MODRIJ: distance beteen the particle pair458: ! MODRIJ: distance beteen the particle pair
468: !459: !
469: !-------------------------------------------------------------------------------460: !-------------------------------------------------------------------------------
470:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ)461:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ)
471: 462: 
472:         ! PERIODIC: whether peiodic boundary conditions are used 
473:         USE COMMONS, ONLY: PERIODIC 
474:  
475:         IMPLICIT NONE463:         IMPLICIT NONE
476: 464: 
477:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ465:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ
478: 466: 
479:         ! DVDR: temporary store for the result467:         ! DVDR: temporary store for the result
480:         DOUBLE PRECISION :: DVDR468:         DOUBLE PRECISION :: DVDR
481: 469: 
482:         ! Unmodified gradient470:         ! Unmodified gradient
483:         DVDR = - 15.D0 * MODRIJ**(-16) - &471:         DVDR = - 15.D0 * MODRIJ**(-16) - &
484:                OPP_K * MODRIJ**(-3) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) - & 472:                OPP_K * MODRIJ**(-3) * SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) - & 
520: !508: !
521: !-------------------------------------------------------------------------------509: !-------------------------------------------------------------------------------
522:     SUBROUTINE OPP_INITIALISE()510:     SUBROUTINE OPP_INITIALISE()
523: 511: 
524:         ! MYUNIT: file unit for main GMIN output512:         ! MYUNIT: file unit for main GMIN output
525:         USE COMMONS, ONLY: MYUNIT513:         USE COMMONS, ONLY: MYUNIT
526: 514: 
527:         IMPLICIT NONE515:         IMPLICIT NONE
528: 516: 
529:         ! Check we have sensible value for the mode517:         ! Check we have sensible value for the mode
530:         IF (OPP_MODE .LT. 0 .OR. OPP_MODE .GT. 3 ) THEN518:         IF (OPP_MODE .NE. 0 .AND. OPP_MODE .NE. 1 .AND. &
531:             WRITE (MYUNIT, *) 'OPP> ERROR: mode must be between 0 and 3'519:             OPP_MODE .NE. 2 ) THEN
 520:             WRITE (MYUNIT, *) 'OPP> ERROR: mode must be 0, 1, 2'
532:             STOP521:             STOP
533:         END IF522:         END IF
534: 523: 
535:         ! Check we have sensible values for the params524:         ! Check we have sensible values for the params
536:         IF (OPP_MODE .EQ. 2 .OR. OPP_MODE .EQ. 3) THEN525:         IF (OPP_MODE .EQ. 2) THEN
537:             IF (OPP_PARAMS .LT. 1 .OR. OPP_PARAMS .GT. 3) THEN526:             IF (OPP_PARAMS .LT. 1 .OR. OPP_PARAMS .GT. 3) THEN
538:                 WRITE (MYUNIT, *) 'OPP> ERROR: params must be between 1 and 3'527:                 WRITE (MYUNIT, *) 'OPP> ERROR: params must be between 1 and 3'
539:                 STOP528:                 STOP
540:             END IF529:             END IF
541:         END IF          530:         END IF          
542: 531: 
543:     END SUBROUTINE OPP_INITIALISE532:     END SUBROUTINE OPP_INITIALISE
544: 533: 
545: !-------------------------------------------------------------------------------534: !-------------------------------------------------------------------------------
546: !535: !
547: ! Takes a step, perturbing the Cartesian coordinates, the unit cell parameters536: ! Takes a step, perturbing the Cartesian coordinates, the unit cell parameters
548: ! and the potential parameters, according to the specified mode.537: ! and the potential parameters, according to the specified mode.
549: ! NP: the index in the main COORDS array to take the step from538: ! NP: the index in the main COORDS array to take the step from
550: !539: !
551: !-------------------------------------------------------------------------------540: !-------------------------------------------------------------------------------
552: 541: 
553:     SUBROUTINE OPP_TAKESTEP(NP)542:     SUBROUTINE OPP_TAKESTEP(NP)
554: 543: 
555:         ! COORDS: full set of Markov chain coordinates544:         ! COORDS: full set of Markov chain coordinates
556:         ! FROZEN: array of atom indices not to move 
557:         ! MYUNIT: file unit of the main GMIN output file545:         ! MYUNIT: file unit of the main GMIN output file
558:         ! NATOMS: number of coordinates546:         ! NATOMS: number of coordinates
559:         ! PERIODIC: whetehr periodic boundary conditions are used547:         ! PERIODIC: whetehr periodic boundary conditions are used
560:         ! PERCOLATET: whether percolation is used to constrain particles548:         ! PERCOLATET: whether percolation is used to constrain particles
561:         ! RADIUS: container radius549:         ! RADIUS: container radius
562:         ! STEP: maximum step size for this move550:         ! STEP: maximum step size for this move
563:         ! TMOVE: specifies which step to do a translational move on551:         ! TMOVE: specifies which step to do a translational move on
564:         ! TWOD: whetehr the system is two dimensional552:         ! TWOD: whetehr the system is two dimensional
565:         USE COMMONS, ONLY: COORDS, FROZEN, MYUNIT, NATOMS, PERIODIC, &553:         USE COMMONS, ONLY: COORDS, MYUNIT, NATOMS, PERIODIC, PERCOLATET, & 
566:                            PERCOLATET, RADIUS, STEP, TMOVE, TWOD554:                            RADIUS, STEP, TMOVE, TWOD
567:         555:         
568:         ! VEC_LEN: function, returns length of a 3D vector556:         ! VEC_LEN: function, returns length of a 3D vector
569:         ! VEC_RANDOM: function, returns a randomly orientated 3D unit vector557:         ! VEC_RANDOM: function, returns a randomly orientated 3D unit vector
570:         USE VEC3, ONLY: VEC_LEN, VEC_RANDOM558:         USE VEC3, ONLY: VEC_LEN, VEC_RANDOM
571: 559: 
572:         IMPLICIT NONE560:         IMPLICIT NONE
573: 561: 
574:         INTEGER, INTENT(IN) :: NP562:         INTEGER, INTENT(IN) :: NP
575: 563: 
576:         ! X: local copy of the coordinates564:         ! X: local copy of the coordinates
577:         ! RAND_STEP: unit vector in a random step direction565:         ! RAND_STEP: unit vector in a random step direction
578:         ! STEP_SIZE: random fraction of the max step size to move566:         ! STEP_SIZE: random fraction of the max step size to move
579:         ! NE_ANGLES: new set of unit cell angles, to be checked567:         ! NE_ANGLES: new set of unit cell angles, to be checked
580:         DOUBLE PRECISION :: X(3*NATOMS), RAND_STEP(3), STEP_SIZE, NEW_ANGLES(3)568:         DOUBLE PRECISION :: X(3*NATOMS), RAND_STEP(3), STEP_SIZE, NEW_ANGLES(3)
581: 569: 
582:         ! I: iteration index570:         ! I: iteration index
583:         ! NMOL: actual number of atoms571:         ! NMOL: actual number of atoms
584:         INTEGER :: I, NMOL572:         INTEGER :: I, NMOL
585: 573: 
586: !        WRITE (MYUNIT, *) 'OPP_TAKESTEP> beginning' 
587: !        FLUSH(MYUNIT) 
588:  
589: !        WRITE (MYUNIT, *) 'Beginning takestep, coordinates are:' 
590: !        DO I = 1, NATOMS 
591: !            WRITE (MYUNIT, *) COORDS(3*I-2:3*I, NP) 
592: !        END DO 
593:  
594:         ! Set the local copy of the coordinates574:         ! Set the local copy of the coordinates
595:         X(:) = COORDS(:, NP)575:         X(:) = COORDS(:, NP)
596: 576: 
597:         ! Set NMOL577:         ! Set NMOL
598:         SELECT CASE(OPP_MODE)578:         NMOL = NATOMS - OPP_MODE - 1
599:         CASE(0) 
600:             NMOL = NATOMS 
601:         CASE(1) 
602:             NMOL = NATOMS - 2 
603:         CASE(2) 
604:             NMOL = NATOMS - 3 
605:         CASE(3) 
606:             NMOL = NATOMS - 1 
607:         END SELECT 
608: 579: 
609:         ! Random translational steps for each of the particles580:         ! Random translational steps for each of the particles
610:         IF (TMOVE(NP)) THEN581:         IF (TMOVE(NP)) THEN
611:             DO I = 1, NMOL582:             DO I = 1, NMOL
612:                 ! Skip if frozen 
613:                 IF (FROZEN(I)) CYCLE 
614:  
615:                 ! Generate either a 2D or 3D random vector583:                 ! Generate either a 2D or 3D random vector
616:                 ! The magnitude will be uniformly distributed in a circle or 584:                 ! The magnitude will be uniformly distributed in a circle or 
617:                 ! sphere with radius given by the value of STEP for this step585:                 ! sphere with radius given by the value of STEP for this step
618:                 IF (TWOD) THEN586:                 IF (TWOD) THEN
619:                     DO587:                     DO
620:                         CALL RANDOM_NUMBER(RAND_STEP(1))588:                         CALL RANDOM_NUMBER(RAND_STEP(1))
621:                         CALL RANDOM_NUMBER(RAND_STEP(2))589:                         CALL RANDOM_NUMBER(RAND_STEP(2))
622:                         RAND_STEP(3) = 0.D0590:                         RAND_STEP(3) = 0.D0
623:                         IF(VEC_LEN(RAND_STEP) .LE. 1.D0) EXIT591:                         IF(VEC_LEN(RAND_STEP) .LE. 1.D0) EXIT
624:                     END DO592:                     END DO
625:                     RAND_STEP(:) = RAND_STEP(:) * STEP(NP)593:                     RAND_STEP(:) = RAND_STEP(:) * STEP(NP)
626:                 ELSE ! 3D594:                 ELSE ! 3D
627:                     RAND_STEP(:) = VEC_RANDOM()595:                     RAND_STEP(:) = VEC_RANDOM()
628:                     CALL RANDOM_NUMBER(STEP_SIZE)596:                     CALL RANDOM_NUMBER(STEP_SIZE)
629:                     ! Sqrt for random volume distribution597:                     ! Sqrt for random volume distribution
630:                     RAND_STEP(:) = RAND_STEP(:) * SQRT(STEP_SIZE * STEP(NP))598:                     RAND_STEP(:) = RAND_STEP(:) * SQRT(STEP_SIZE * STEP(NP))
631:                 END IF599:                 END IF
632:                 ! Add the geenrated step onto the coordinates600:                 ! Add the geenrated ste onto the coordinates
633:                 X(I*3-2:I*3) = X(I*3-2:I*3) + RAND_STEP(:)601:                 X(I*3-2:I*3) = X(I*3-2:I*3) + RAND_STEP(:)
634:             END DO602:             END DO
635:         END IF603:         END IF
636: 604: 
637:         ! If not periodic and we're using a radius container, bring particles in605:         ! If not periodic and we're using a radius container, bring particles in
638:         IF (.NOT. PERIODIC .AND. .NOT. PERCOLATET) THEN606:         IF (.NOT. PERIODIC .AND. .NOT. PERCOLATET) THEN
639:             DO I = 1, NMOL607:             DO I = 1, NMOL
640:                 IF (VEC_LEN(X(I*3-2:I*3)) .GT. RADIUS) THEN608:                 IF (VEC_LEN(X(I*3-2:I*3)) .GT. RADIUS) THEN
641:                     WRITE(MYUNIT, *) 'OPP_TAKESTEP> ', &609:                     WRITE(MYUNIT, *) 'OPP_TAKESTEP> ', &
642:                         'coord outside container, bringing in'610:                         'coord outside container, bringing in'
643:  
644:                     ! Skip if frozen 
645:                     IF (FROZEN(I)) CYCLE 
646:                     X(I*3-2:I*3) = X(I*3-2:I*3) - &611:                     X(I*3-2:I*3) = X(I*3-2:I*3) - &
647:                         SQRT(RADIUS) * NINT(X(I*3-2:I*3) / SQRT(RADIUS))612:                         SQRT(RADIUS) * NINT(X(I*3-2:I*3) / SQRT(RADIUS))
648:                 END IF613:                 END IF
649:             END DO614:             END DO
650:         END IF615:         END IF
651: 616: 
652:         ! Box length steps617:         ! Box length steps
653:         IF ((OPP_MODE .EQ. 1 .OR. OPP_MODE .EQ. 2) .AND. &618:         IF (OPP_MODE .GE. 1) THEN
654:             .NOT. FROZEN(NATOMS)) THEN 
655:             X(3*NATOMS-2:3*NATOMS) = X(3*NATOMS-2:3*NATOMS) + &619:             X(3*NATOMS-2:3*NATOMS) = X(3*NATOMS-2:3*NATOMS) + &
656:                                      VEC_RANDOM() * MAX_LENGTH_STEP620:                                      VEC_RANDOM() * MAX_LENGTH_STEP
657:             ! Fiddle to lock in cubic 
658: !            X(3*NATOMS-1) = X(3*NATOMS-2) 
659: !            X(3*NATOMS  ) = X(3*NATOMS-2) 
660:         END IF621:         END IF
661: 622: 
662:         ! Box angle steps623:         ! Box angle steps
663:         IF ((OPP_MODE .EQ. 1 .OR. OPP_MODE .EQ. 2) .AND. &624:         IF (OPP_MODE .GE. 1) THEN
664:             .NOT. FROZEN(NATOMS-1)) THEN 
665:             DO625:             DO
666:                 NEW_ANGLES(:) = X(3*NATOMS-5:3*NATOMS-3) + &626:                 NEW_ANGLES(:) = X(3*NATOMS-5:3*NATOMS-3) + &
667:                                 VEC_RANDOM() * MAX_ANGLE_STEP627:                                 VEC_RANDOM() * MAX_ANGLE_STEP
668: 628: 
669:                 ! See whether the unit cell we've generated is valid629:                 ! See whether the unit cell we've generated is valid
670:                 IF (CHECK_ANGLES(NEW_ANGLES(:))) THEN630:                 IF (CHECK_ANGLES(NEW_ANGLES(:))) THEN
671:                     X(3*NATOMS-5:3*NATOMS-3) =  NEW_ANGLES(:)631:                     X(3*NATOMS-5:3*NATOMS-3) =  NEW_ANGLES(:)
672:                     EXIT632:                     EXIT
673:                 END IF633:                 END IF
674:             END DO634:             END DO
675:         END IF635:         END IF
676: 636: 
677:         ! Parameter steps637:         ! Parameter steps
678:         IF (OPP_MODE .EQ. 2 .AND. .NOT. FROZEN(NATOMS-2)) THEN638:         IF (OPP_MODE .EQ. 2) THEN
679:             IF (BTEST(OPP_PARAMS, 0)) THEN639:             IF (BTEST(OPP_PARAMS, 0)) THEN
680:                 CALL RANDOM_NUMBER(STEP_SIZE)640:                 CALL RANDOM_NUMBER(STEP_SIZE)
681:                 X(3*NATOMS-8) = X(3*NATOMS-8) + &641:                 X(3*NATOMS-8) = X(3*NATOMS-8) + &
682:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_K_STEP642:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_K_STEP
683:             END IF643:             END IF
684:             IF (BTEST(OPP_PARAMS, 1)) THEN644:             IF (BTEST(OPP_PARAMS, 1)) THEN
685:                 CALL RANDOM_NUMBER(STEP_SIZE)645:                 CALL RANDOM_NUMBER(STEP_SIZE)
686:                 X(3*NATOMS-7) = X(3*NATOMS-7) + &646:                 X(3*NATOMS-7) = X(3*NATOMS-7) + &
687:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_PHI_STEP647:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_PHI_STEP
688:             END IF 
689:         ELSE IF (OPP_MODE .EQ. 3 .AND. .NOT. FROZEN(NATOMS)) THEN 
690:             IF (BTEST(OPP_PARAMS, 0)) THEN 
691:                 CALL RANDOM_NUMBER(STEP_SIZE) 
692:                 X(3*NATOMS-2) = X(3*NATOMS-2) + & 
693:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_K_STEP 
694:             END IF 
695:             IF (BTEST(OPP_PARAMS, 1)) THEN 
696:                 CALL RANDOM_NUMBER(STEP_SIZE) 
697:                 X(3*NATOMS-1) = X(3*NATOMS-1) + & 
698:                                 (STEP_SIZE - 0.5D0) * 2.D0 * MAX_PHI_STEP 
699:             END IF            648:             END IF            
700:         END IF649:         END IF
701: 650: 
702:         COORDS(:, NP) = X(:)651:         COORDS(:, NP) = X(:)
703: 652: 
704: !        WRITE (MYUNIT, *) 'OPP_TAKESTEP> ending' 
705: !        FLUSH(MYUNIT) 
706: !        WRITE (MYUNIT, *) 'Ending takestep, coordinates are:' 
707: !        DO I = 1, NATOMS 
708: !            WRITE (MYUNIT, *) COORDS(3*I-2:3*I, NP) 
709: !        END DO 
710:  
711:     END SUBROUTINE OPP_TAKESTEP653:     END SUBROUTINE OPP_TAKESTEP
712: 654: 
713: !-------------------------------------------------------------------------------655: !-------------------------------------------------------------------------------
714: !656: !
715: ! Resets all coordinates to within the unit cell657: ! Resets all coordinates to within the unit cell
716: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled658: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled
717: ! by the box lengths.659: ! by the box lengths.
718: ! COORDS: coordinates of the particles660: ! COORDS: coordinates of the particles
719: !661: !
720: !-------------------------------------------------------------------------------662: !-------------------------------------------------------------------------------
722: 664: 
723:         ! BOXLX, BOXLY, BOXLZ: dimensions of box665:         ! BOXLX, BOXLY, BOXLZ: dimensions of box
724:         ! NATOMS: number of particles666:         ! NATOMS: number of particles
725:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ, NATOMS667:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ, NATOMS
726: 668: 
727:         IMPLICIT NONE669:         IMPLICIT NONE
728: 670: 
729:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)671:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)
730: 672: 
731:         ! Sort out the unit cell sizes, if they are varying673:         ! Sort out the unit cell sizes, if they are varying
732:         IF (OPP_MODE .EQ. 1 .OR. OPP_MODE .EQ. 2) THEN674:         IF (OPP_MODE .GE. 1) THEN
733:             BOXLX = COORDS(3*NATOMS - 2)675:             BOXLX = COORDS(3*NATOMS - 2)
734:             BOXLY = COORDS(3*NATOMS - 1)676:             BOXLY = COORDS(3*NATOMS - 1)
735:             BOXLZ = COORDS(3*NATOMS)677:             BOXLZ = COORDS(3*NATOMS)
736:         END IF678:         END IF
737: 679: 
738:         ! Reset coordinates to within the box680:         ! Reset coordinates to within the box
739:         SELECT CASE (OPP_MODE)681:         SELECT CASE (OPP_MODE)
740:         CASE(0)682:         CASE(0)
741:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:))683:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:))
742:         CASE(1)684:         CASE(1)
749:         CASE(2)691:         CASE(2)
750:            COORDS(1:3*NATOMS-9) = COORDS(1:3*NATOMS-9) - &692:            COORDS(1:3*NATOMS-9) = COORDS(1:3*NATOMS-9) - &
751:                                    FLOOR(COORDS(1:3*NATOMS-9))693:                                    FLOOR(COORDS(1:3*NATOMS-9))
752:             ! Make the lattice angles modulo 2*pi694:             ! Make the lattice angles modulo 2*pi
753:             ! They will be checked later695:             ! They will be checked later
754:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &696:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &
755:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))697:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))
756:             ! Make OPP_PHI modulo 2*pi698:             ! Make OPP_PHI modulo 2*pi
757:             COORDS(3*NATOMS-7) = COORDS(3*NATOMS-7) - &699:             COORDS(3*NATOMS-7) = COORDS(3*NATOMS-7) - &
758:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-7)/(2.D0*PI))700:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-7)/(2.D0*PI))
759:             ! Reset the unused coordinate to zero 
760:             COORDS(3*NATOMS-6) = 0.D0 
761:         CASE(3) 
762:             ! Make OPP_PHI modulo 2*pi 
763:             COORDS(3*NATOMS-1) = COORDS(3*NATOMS-1) - & 
764:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-1)/(2.D0*PI)) 
765:             ! Reset the unused coordinate to zero 
766:             COORDS(3*NATOMS) = 0.D0 
767:         END SELECT701:         END SELECT
768: 702: 
769:     END SUBROUTINE PERIODIC_RESET703:     END SUBROUTINE PERIODIC_RESET
770: 704: 
771: !-------------------------------------------------------------------------------705: !-------------------------------------------------------------------------------
772: !706: !
773: ! Finds the energy and gradients for the othogonal fixed periodic system707: ! Finds the energy and gradients for the othogonal fixed periodic system
774: ! IDI, IDJ: id numbers of the particles708: ! IDI, IDJ: id numbers of the particles
775: ! COORDS_I: coordinates of particle I709: ! COORDS_I: coordinates of particle I
776: ! COORDS_J: coordinates of particle J710: ! COORDS_J: coordinates of particle J
1066:                             ! Divide gradient by 2 for images of the same atom,1000:                             ! Divide gradient by 2 for images of the same atom,
1067:                             ! to avoid double counting1001:                             ! to avoid double counting
1068:                             IF(PRESENT(GRAD_PARAMS) .AND. OPP_MODE .EQ. 2)&1002:                             IF(PRESENT(GRAD_PARAMS) .AND. OPP_MODE .EQ. 2)&
1069:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + &1003:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + &
1070:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &1004:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &
1071:                                     CALC_GRAD_PARAMS(MODYIJ)1005:                                     CALC_GRAD_PARAMS(MODYIJ)
1072: 1006: 
1073:                         END IF ! GTEST1007:                         END IF ! GTEST
1074: 1008: 
1075:                         IF (STEST) THEN1009:                         IF (STEST) THEN
1076:                             CALL HESS_CONTRIB_PER(RIJ(:), BD, IDI, IDJ)1010:                             CALL HESS_CONTRIB(RIJ(:), BD, IDI, IDJ)
1077:                         END IF ! End second derivative1011:                         END IF ! End second derivative
1078:                     END IF ! End if less than cutoff1012:                     END IF ! End if less than cutoff
1079:                 END DO ! z loop1013:                 END DO ! z loop
1080:             END DO ! y loop1014:             END DO ! y loop
1081:         END DO ! x loop1015:         END DO ! x loop
1082: 1016: 
1083:     END SUBROUTINE OPP_PER_TRI1017:     END SUBROUTINE OPP_PER_TRI
1084: 1018: 
1085: !-------------------------------------------------------------------------------1019: !-------------------------------------------------------------------------------
1086: !1020: !
1087: ! Calculates the second derivative contributions for a cluster1021: ! Calculates the second derivative contributions
1088: ! RIJ: vector between particles in fractional coordinates 
1089: ! BD: box derivative information 
1090: ! IDI, IDJ: indices of the atom pair 
1091: ! 
1092: !------------------------------------------------------------------------------- 
1093:     SUBROUTINE HESS_CONTRIB(RIJ, IDI, IDJ) 
1094:  
1095:         ! Number of atoms, including unit cell parameters 
1096:         USE COMMONS, ONLY: NATOMS 
1097:  
1098:         ! HESS: the Hessian matrix 
1099:         USE MODHESS, ONLY: HESS 
1100:  
1101:         IMPLICIT NONE 
1102:  
1103:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3) 
1104:         INTEGER, INTENT(IN) :: IDI, IDJ 
1105:  
1106:         ! MODRIJ: particle-particle distance in Cartesian space 
1107:         ! DVDR: derivative of pair potential wrt MODYIJ 
1108:         ! D2VDR2: D(DVDR/MODRIJ)/DMODRIJ * 1/MODRIJ 
1109:         ! TEMP, TEMP_GRAD: temporary calculation store 
1110:         DOUBLE PRECISION :: MODRIJ, DVDR, D2VDR2, TEMP, TEMP_GRAD(3) 
1111:  
1112:         ! I, J: loop indices 
1113:         INTEGER :: I, J 
1114:  
1115:         ! Calculate some factors 
1116:         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
1117:         DVDR = CALC_GRAD(MODRIJ) / MODRIJ 
1118:         D2VDR2 = (CALC_SEC(MODRIJ) - DVDR) / MODRIJ**2 
1119:  
1120:         DO I = 1, 3 ! outer loop over x, y, z 
1121:             DO J = 1, 3 ! inner loop over x, y, z 
1122:                 ! Same atom, same or different coordinate 
1123:                 HESS(3*IDI-3+I,3*IDI-3+J) = HESS(3*IDI-3+I,3*IDI-3+J) + & 
1124:                     D2VDR2 * RIJ(I) * RIJ(J) 
1125:                 HESS(3*IDJ-3+I,3*IDJ-3+J) = HESS(3*IDJ-3+I,3*IDJ-3+J) + & 
1126:                     D2VDR2 * RIJ(I) * RIJ(J) 
1127:                 ! Same atom, same coordinate 
1128:                 IF (I .EQ. J) THEN 
1129:                     HESS(3*IDI-3+I,3*IDI-3+J) = HESS(3*IDI-3+I,3*IDI-3+J) + DVDR 
1130:                     HESS(3*IDJ-3+I,3*IDJ-3+J) = HESS(3*IDJ-3+I,3*IDJ-3+J) + DVDR 
1131:                 END IF 
1132:                 ! Different atom, same or different coordinate 
1133:                 HESS(3*IDI-3+I,3*IDJ-3+J) = HESS(3*IDI-3+I,3*IDJ-3+J) - & 
1134:                     D2VDR2 * RIJ(I) * RIJ(J) 
1135:                 HESS(3*IDJ-3+I,3*IDI-3+J) = HESS(3*IDJ-3+I,3*IDI-3+J) - & 
1136:                     D2VDR2 * RIJ(I) * RIJ(J) 
1137:                 ! Same atom, same coordinate 
1138:                 IF (I .EQ. J) THEN 
1139:                     HESS(3*IDI-3+I,3*IDJ-3+J) = HESS(3*IDI-3+I,3*IDJ-3+J) - DVDR 
1140:                     HESS(3*IDJ-3+I,3*IDI-3+J) = HESS(3*IDJ-3+I,3*IDI-3+J) - DVDR 
1141:                 END IF 
1142:             END DO ! end inner loop over x, y, z 
1143:         END DO ! end outer loop over x, y, z 
1144:  
1145:         IF (OPP_MODE .EQ. 3) THEN 
1146:             ! Potential parameter contributions 
1147:             CALL CALC_HESS_PARAMS(RIJ, IDI, IDJ) 
1148:         END IF 
1149:  
1150:     END SUBROUTINE HESS_CONTRIB 
1151:  
1152: !------------------------------------------------------------------------------- 
1153: ! 
1154: ! Calculates the second derivative contributions for a periodic system 
1155: ! RIJ: vector between particles in fractional coordinates1022: ! RIJ: vector between particles in fractional coordinates
1156: ! BD: box derivative information1023: ! BD: box derivative information
1157: ! IDI, IDJ: indices of the atom pair1024: ! IDI, IDJ: indices of the atom pair
1158: !1025: !
1159: !-------------------------------------------------------------------------------1026: !-------------------------------------------------------------------------------
1160:     SUBROUTINE HESS_CONTRIB_PER(RIJ, BD, IDI, IDJ)1027:     SUBROUTINE HESS_CONTRIB(RIJ, BD, IDI, IDJ)
1161: 1028: 
1162:         ! Number of atoms, including unit cell parameters1029:         ! Number of atoms, including unit cell parameters
1163:         USE COMMONS, ONLY: NATOMS1030:         USE COMMONS, ONLY: NATOMS
1164: 1031: 
1165:         ! HESS: the Hessian matrix1032:         ! HESS: the Hessian matrix
1166:         USE MODHESS, ONLY: HESS1033:         USE MODHESS, ONLY: HESS
1167: 1034: 
1168:         ! BOX_DERIV: type for storing box derivative information1035:         ! BOX_DERIV: type for storing box derivative information
1169:         USE TRICLINIC_MOD, ONLY: BOX_DERIV1036:         USE TRICLINIC_MOD, ONLY: BOX_DERIV
1170: 1037: 
1226: 1093: 
1227:         ! Pure unit cell parameters1094:         ! Pure unit cell parameters
1228:         DO I = 1, 6 ! outer loop over cell parameters1095:         DO I = 1, 6 ! outer loop over cell parameters
1229:             ! Diagonals1096:             ! Diagonals
1230:             TEMP = D2VDY2 * &1097:             TEMP = D2VDY2 * &
1231:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:)))**2 + &1098:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:)))**2 + &
1232:                 DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &1099:                 DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &
1233:                 MATMUL(BD%DERIV(:,:,I), RIJ(:))) + &1100:                 MATMUL(BD%DERIV(:,:,I), RIJ(:))) + &
1234:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,I), RIJ(:))))1101:                 DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,I), RIJ(:))))
1235:             HESS(3*NATOMS-6+I,3*NATOMS-6+I) = HESS(3*NATOMS-6+I,3*NATOMS-6+I) + TEMP1102:             HESS(3*NATOMS-6+I,3*NATOMS-6+I) = HESS(3*NATOMS-6+I,3*NATOMS-6+I) + TEMP
 1103:             IF (I .EQ. 1) THEN
 1104: !                WRITE (*, *) 'IDI = ', IDI, 'IDJ = ', IDJ, 'Hessian contribution = ', TEMP
 1105: !                WRITE (*, *) 'MODYIJ = ', MODYIJ
 1106: !                WRITE (*, *) 'RIJ = ', RIJ
 1107: !                WRITE (*, *) 'YIJ = ', YIJ
 1108: !                WRITE (*, *) 'CALC_GRAD(MODYIJ) = ', CALC_GRAD(MODYIJ)
 1109: !                WRITE (*, *) 'DVDY =', DVDY
 1110: !                WRITE (*, *) 'D2VDY2 = ', D2VDY2
 1111:             END IF
 1112:             ! Off diagonals
1236:             DO J = I + 1, 6 ! inner loop over cell parameters1113:             DO J = I + 1, 6 ! inner loop over cell parameters
1237:                 TEMP = D2VDY2 * &1114:                 TEMP = D2VDY2 * &
1238:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) * &1115:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,I), RIJ(:))) * &
1239:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &1116:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &
1240:                     DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &1117:                     DVDY * (DOT_PRODUCT(MATMUL(BD%DERIV(:,:,I), RIJ(:)), &
1241:                     MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &1118:                     MATMUL(BD%DERIV(:,:,J), RIJ(:))) + &
1242:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,J), RIJ(:))))1119:                     DOT_PRODUCT(YIJ(:), MATMUL(BD%DERIV2(:,:,I,J), RIJ(:))))
1243:                 HESS(3*NATOMS-6+I,3*NATOMS-6+J) = &1120:                 HESS(3*NATOMS-6+I,3*NATOMS-6+J) = &
1244:                     HESS(3*NATOMS-6+I,3*NATOMS-6+J) + TEMP1121:                     HESS(3*NATOMS-6+I,3*NATOMS-6+J) + TEMP
1245:                 HESS(3*NATOMS-6+J,3*NATOMS-6+I) = &1122:                 HESS(3*NATOMS-6+J,3*NATOMS-6+I) = &
1246:                     HESS(3*NATOMS-6+J,3*NATOMS-6+I) + TEMP1123:                     HESS(3*NATOMS-6+J,3*NATOMS-6+I) + TEMP
1247:             END DO ! end inner loop over cell parameters1124:             END DO ! end inner loop over cell parameters
1248:         END DO ! end outer loop over cell parameters1125:         END DO ! end outer loop over cell parameters
1249: 1126: 
1250:         IF (OPP_MODE .EQ. 2) THEN1127:         IF (OPP_MODE .EQ. 2) THEN
1251:             ! Potential parameter contributions1128:             ! Potential parameter contributions
1252:             CALL CALC_HESS_PARAMS_PER(RIJ, BD, IDI, IDJ)1129:             CALL CALC_HESS_PARAMS(RIJ, BD, IDI, IDJ)
1253:         END IF1130:         END IF
1254: 1131: 
1255:     END SUBROUTINE HESS_CONTRIB_PER1132:     END SUBROUTINE HESS_CONTRIB
1256: 1133: 
1257: !-------------------------------------------------------------------------------1134: !-------------------------------------------------------------------------------
1258: !1135: !
1259: ! Calculates the contribution to the derivatives of the potential wrt the1136: ! Calculates the contribution to the derivatives of the potential wrt the
1260: ! potential parameters OPP_K and OPP_PHI1137: ! potential parameters OPP_K and OPP_PHI
1261: ! MODRIJ: distance between the particles for this contribution1138: ! MODYIJ: distance between the particles for this contribution
1262: !1139: !
1263: !-------------------------------------------------------------------------------1140: !-------------------------------------------------------------------------------
1264: 1141: 
1265:     PURE FUNCTION CALC_GRAD_PARAMS(MODRIJ)1142:     PURE FUNCTION CALC_GRAD_PARAMS(MODYIJ)
1266:  
1267:         USE COMMONS, ONLY: PERIODIC 
1268: 1143: 
1269:         IMPLICIT NONE1144:         IMPLICIT NONE
1270: 1145: 
1271:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(3)1146:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(3)
1272:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ1147:         DOUBLE PRECISION, INTENT(IN) :: MODYIJ
1273: 1148: 
1274:         ! S, C: sin and cos factor that appears repeatedly1149:         ! S, C: sin and cos factor that appears repeatedly
1275:         ! T: working store for the output1150:         ! T: working store for the output
1276:         DOUBLE PRECISION :: S, C1151:         DOUBLE PRECISION :: S, C
1277:         DOUBLE PRECISION :: T(3)1152:         DOUBLE PRECISION :: T(3)
1278: 1153: 
1279:         ! Initialise output1154:         ! Initialise output
1280:         T(:) = 0.D01155:         T(:) = 0.D0
1281: 1156: 
1282:         ! Calculate factors1157:         ! Calculate factors
1283:         S = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI)1158:         S = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
1284:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)1159:         C = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI)
1285: 1160: 
1286:         IF (BTEST(OPP_PARAMS, 0)) THEN1161:         IF (BTEST(OPP_PARAMS, 0)) THEN
1287:             ! OPP_K derivative1162:             ! OPP_K derivative
1288:             ! Unmodified potential term1163:             ! Unmodified potential term
1289:             T(1) = - SIN(OPP_K * (MODRIJ-1) - OPP_PHI) * (MODRIJ-1) / MODRIJ**31164:             T(1) = - SIN(OPP_K * (MODYIJ-1) - OPP_PHI) * (MODYIJ-1) / MODYIJ**3
1290:             ! Term to make the potential go to zero at OPP_RCUT1165:             ! Term to make the potential go to zero at OPP_RCUT
1291:             T(1) = T(1) + S * (OPP_RCUT - 1) / OPP_RCUT **31166:             T(1) = T(1) + S * (OPP_RCUT - 1) / OPP_RCUT **3
1292:             ! Terms to make the first derivative go to zero at OPP_RCUT1167:             ! Terms to make the first derivative go to zero at OPP_RCUT
1293:             T(1) = T(1) + (MODRIJ - OPP_RCUT) * &1168:             T(1) = T(1) + (MODYIJ - OPP_RCUT) * &
1294:                 (C * OPP_K * (OPP_RCUT - 1) + S - &1169:                 (C * OPP_K * (OPP_RCUT - 1) + S - &
1295:                 S * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT ) / OPP_RCUT**31170:                  S * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT ) / OPP_RCUT**3
1296:             ! Terms to make the second derivative go to zero at OPP_RCUT1171:             ! Terms to make the second derivative go to zero at OPP_RCUT
1297: !            T(1) = T(1) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &1172: !            T(1) = T(1) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &
1298: !                (S * OPP_K**2 * (OPP_RCUT - 1) / OPP_RCUT**3 - &1173: !                   (S * OPP_K**2 * (OPP_RCUT - 1) / OPP_RCUT**3 - &
1299: !                C * 2.D0 * OPP_K / OPP_RCUT**3 + &1174: !                    C * 2.D0 * OPP_K / OPP_RCUT**3 + &
1300: !                C * 6.D0 * OPP_K * (OPP_RCUT - 1) / OPP_RCUT**4 + &1175: !                    C * 6.D0 * OPP_K * (OPP_RCUT - 1) / OPP_RCUT**4 + &
1301: !                S * 6.D0 / OPP_RCUT**4 - &1176: !                    S * 6.D0 / OPP_RCUT**4 - &
1302: !                S * 12.D0 * (OPP_RCUT - 1) / OPP_RCUT**5)1177: !                    S * 12.D0 * (OPP_RCUT - 1) / OPP_RCUT**5)
1303:         END IF1178:         END IF
1304:         IF (BTEST(OPP_PARAMS, 1)) THEN1179:         IF (BTEST(OPP_PARAMS, 1)) THEN
1305:             ! OPP_PHI derivative1180:             ! OPP_PHI derivative
1306:             ! Unmodified potential term1181:             ! Unmodified potential term
1307:             T(2) = SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) / MODRIJ**31182:             T(2) = SIN(OPP_K * (MODYIJ - 1) - OPP_PHI) / MODYIJ**3
1308:             ! Term to make the potential go to zero at OPP_RCUT1183:             ! Term to make the potential go to zero at OPP_RCUT
1309:             T(2) = T(2) - S / OPP_RCUT **31184:             T(2) = T(2) - S / OPP_RCUT **3
1310:             ! Terms to make the first derivative go to zero at OPP_RCUT1185:             ! Terms to make the first derivative go to zero at OPP_RCUT
1311:             T(2) = T(2) + (MODRIJ - OPP_RCUT) * &1186:             T(2) = T(2) + (MODYIJ - OPP_RCUT) * &
1312:                 ( - C * OPP_K + S * 3.D0 / OPP_RCUT ) / OPP_RCUT**31187:                    ( - C * OPP_K + S * 3.D0 / OPP_RCUT ) / OPP_RCUT**3
1313:             ! Terms to make the second derivative go to zero at OPP_RCUT1188:             ! Terms to make the second derivative go to zero at OPP_RCUT
1314: !            T(2) = T(2) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &1189: !            T(2) = T(2) - 0.5D0 * (MODYIJ - OPP_RCUT)**2 * &
1315: !                ( - S * OPP_K**2 / OPP_RCUT**3 - &1190: !                   ( - S * OPP_K**2 / OPP_RCUT**3 - &
1316: !                C * 6.D0 * OPP_K / OPP_RCUT**4 + S * 12.D0 / OPP_RCUT**5)1191: !                    C * 6.D0 * OPP_K / OPP_RCUT**4 + S * 12.D0 / OPP_RCUT**5)           
1317:         END IF1192:         END IF
1318: 1193: 
1319:         CALC_GRAD_PARAMS(:) = T(:)1194:         CALC_GRAD_PARAMS(:) = T(:)
1320: 1195: 
1321:     END FUNCTION CALC_GRAD_PARAMS1196:     END FUNCTION CALC_GRAD_PARAMS
1322: 1197: 
1323: !-------------------------------------------------------------------------------1198: !-------------------------------------------------------------------------------
1324: !1199: !
1325: ! Calculates the second derivative contributions for the potential parameters1200: ! Calculates the second derivative contributions for the potential parameters
1326: ! for a cluster 
1327: ! RIJ: vector between particles 
1328: ! IDI, IDJ: indices of the atom pair 
1329: ! 
1330: !------------------------------------------------------------------------------- 
1331:  
1332:     SUBROUTINE CALC_HESS_PARAMS(RIJ, IDI, IDJ) 
1333:  
1334:         ! Number of atoms, including unit cell parameters 
1335:         USE COMMONS, ONLY: NATOMS 
1336:  
1337:         ! HESS: the Hessian matrix 
1338:         USE MODHESS, ONLY: HESS 
1339:  
1340:         IMPLICIT NONE 
1341:  
1342:         DOUBLE PRECISION, INTENT(IN) :: RIJ(3) 
1343:         INTEGER, INTENT(IN) :: IDI, IDJ 
1344:  
1345:         ! MODRIJ: particle-particle distance 
1346:         ! S, C, SC, CC: recurring sin and cos terms 
1347:         ! TEMP, TEMP_GRAD: temporary calculation stores 
1348:         ! D2VDK2, D2VDPHI2, D2VDKDPHI: pure parameter second derivatives 
1349:         ! D2VDRDK, D2VDRDPHI; mixed second derivatives of pair potential 
1350:         DOUBLE PRECISION :: MODRIJ, S, C, SC, CC, TEMP, TEMP_GRAD(3) 
1351:         DOUBLE PRECISION :: D2VDK2, D2VDPHI2, D2VDKDPHI, D2VDRDK, D2VDRDPHI 
1352:  
1353:         ! I: loop index 
1354:         INTEGER :: I 
1355:  
1356:         ! Calculate some factors 
1357:         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
1358:         S = SIN(OPP_K * (MODRIJ - 1) - OPP_PHI) 
1359:         C = COS(OPP_K * (MODRIJ - 1) - OPP_PHI) 
1360:         SC = SIN(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
1361:         CC = COS(OPP_K * (OPP_RCUT - 1) - OPP_PHI) 
1362:  
1363:         IF (BTEST(OPP_PARAMS, 0)) THEN 
1364:             ! OPP_K pure derivative 
1365:             ! Unmodified potential term 
1366:             D2VDK2 = - C * (MODRIJ - 1)**2 / MODRIJ**3 
1367:             HESS(3*NATOMS-2,3*NATOMS-2) = HESS(3*NATOMS-2,3*NATOMS-2) + D2VDK2 
1368:  
1369:             ! OPP_K and position mixed derivative 
1370:             ! Unmodified potential term 
1371:             D2VDRDK = (- C * OPP_K * (MODRIJ - 1) - S + & 
1372:                 S * 3.D0 * (MODRIJ - 1) / MODRIJ) / MODRIJ**3 
1373:             D2VDRDK = D2VDRDK / MODRIJ 
1374:             TEMP_GRAD(:) = D2VDRDK * RIJ(:) 
1375:             HESS(3*IDI-2:3*IDI,3*NATOMS-2) = & 
1376:                 HESS(3*IDI-2:3*IDI,3*NATOMS-2) + TEMP_GRAD(:) 
1377:             HESS(3*NATOMS-2,3*IDI-2:3*IDI) = & 
1378:                 HESS(3*NATOMS-2,3*IDI-2:3*IDI) + TEMP_GRAD(:) 
1379:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-2) = & 
1380:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-2) - TEMP_GRAD(:) 
1381:             HESS(3*NATOMS-2,3*IDJ-2:3*IDJ) = & 
1382:                 HESS(3*NATOMS-2,3*IDJ-2:3*IDJ) - TEMP_GRAD(:) 
1383:         END IF         
1384:  
1385:         IF (BTEST(OPP_PARAMS, 1)) THEN 
1386:             ! OPP_PHI pure derivative 
1387:             ! Unmodified potential term 
1388:             D2VDPHI2 = - C / MODRIJ**3 
1389:             HESS(3*NATOMS-1,3*NATOMS-1) = HESS(3*NATOMS-1,3*NATOMS-1) + D2VDPHI2 
1390:  
1391:             ! OPP_PHI and position mixed derivative 
1392:             ! Unmodified potential term 
1393:             D2VDRDPHI = (C * OPP_K - S * 3.D0 / MODRIJ) / MODRIJ**3 
1394:             D2VDRDPHI = D2VDRDPHI  / MODRIJ 
1395:             TEMP_GRAD(:) = D2VDRDPHI * RIJ(:) 
1396:             HESS(3*IDI-2:3*IDI,3*NATOMS-1) = & 
1397:                 HESS(3*IDI-2:3*IDI,3*NATOMS-1) + TEMP_GRAD(:) 
1398:             HESS(3*NATOMS-1,3*IDI-2:3*IDI) = & 
1399:                 HESS(3*NATOMS-1,3*IDI-2:3*IDI) + TEMP_GRAD(:) 
1400:             HESS(3*IDJ-2:3*IDJ,3*NATOMS-1) = & 
1401:                 HESS(3*IDJ-2:3*IDJ,3*NATOMS-1) - TEMP_GRAD(:) 
1402:             HESS(3*NATOMS-1,3*IDJ-2:3*IDJ) = & 
1403:                 HESS(3*NATOMS-1,3*IDJ-2:3*IDJ) - TEMP_GRAD(:) 
1404:         END IF 
1405:  
1406:         IF (BTEST(OPP_PARAMS, 0) .AND. BTEST(OPP_PARAMS, 1)) THEN 
1407:             ! OPP_K and OPP_PHI derivative 
1408:             ! Unmodified potential term 
1409:             D2VDKDPHI = C * (MODRIJ - 1) / MODRIJ**3 
1410:             HESS(3*NATOMS-2,3*NATOMS-1) = HESS(3*NATOMS-2,3*NATOMS-1)+D2VDKDPHI 
1411:             HESS(3*NATOMS-1,3*NATOMS-2) = HESS(3*NATOMS-1,3*NATOMS-2)+D2VDKDPHI 
1412:         END IF 
1413:  
1414:     END SUBROUTINE CALC_HESS_PARAMS 
1415:  
1416: !------------------------------------------------------------------------------- 
1417: ! 
1418: ! Calculates the second derivative contributions for the potential parameters in 
1419: ! a periodic system. 
1420: ! RIJ: vector between particles in fractional coordinates1201: ! RIJ: vector between particles in fractional coordinates
1421: ! BD: box derivative information1202: ! BD: box derivative information
1422: ! IDI, IDJ: indices of the atom pair1203: ! IDI, IDJ: indices of the atom pair
1423: !1204: !
1424: !-------------------------------------------------------------------------------1205: !-------------------------------------------------------------------------------
1425: 1206: 
1426:     SUBROUTINE CALC_HESS_PARAMS_PER(RIJ, BD, IDI, IDJ)1207:     SUBROUTINE CALC_HESS_PARAMS(RIJ, BD, IDI, IDJ)
1427: 1208: 
1428:         ! Number of atoms, including unit cell parameters1209:         ! Number of atoms, including unit cell parameters
1429:         USE COMMONS, ONLY: NATOMS1210:         USE COMMONS, ONLY: NATOMS
1430: 1211: 
1431:         ! HESS: the Hessian matrix1212:         ! HESS: the Hessian matrix
1432:         USE MODHESS, ONLY: HESS1213:         USE MODHESS, ONLY: HESS
1433: 1214: 
1434:         ! BOX_DERIV: type for storing box derivative information1215:         ! BOX_DERIV: type for storing box derivative information
1435:         USE TRICLINIC_MOD, ONLY: BOX_DERIV1216:         USE TRICLINIC_MOD, ONLY: BOX_DERIV
1436: 1217: 
1550:             D2VDKDPHI = D2VDKDPHI - CC * (OPP_RCUT - 1) / OPP_RCUT**31331:             D2VDKDPHI = D2VDKDPHI - CC * (OPP_RCUT - 1) / OPP_RCUT**3
1551:             ! Terms to make the first derivative go to zero at OPP_RCUT1332:             ! Terms to make the first derivative go to zero at OPP_RCUT
1552:             D2VDKDPHI = D2VDKDPHI + (MODYIJ - OPP_RCUT) * &1333:             D2VDKDPHI = D2VDKDPHI + (MODYIJ - OPP_RCUT) * &
1553:                 (SC * OPP_K * (OPP_RCUT - 1) -  CC + &1334:                 (SC * OPP_K * (OPP_RCUT - 1) -  CC + &
1554:                 CC * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT) / OPP_RCUT**31335:                 CC * 3.D0 * (OPP_RCUT - 1) / OPP_RCUT) / OPP_RCUT**3
1555:             D2VDKDPHI = D2VDKDPHI * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)1336:             D2VDKDPHI = D2VDKDPHI * MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ)
1556:             HESS(3*NATOMS-8,3*NATOMS-7) = HESS(3*NATOMS-8,3*NATOMS-7)+D2VDKDPHI1337:             HESS(3*NATOMS-8,3*NATOMS-7) = HESS(3*NATOMS-8,3*NATOMS-7)+D2VDKDPHI
1557:             HESS(3*NATOMS-7,3*NATOMS-8) = HESS(3*NATOMS-7,3*NATOMS-8)+D2VDKDPHI1338:             HESS(3*NATOMS-7,3*NATOMS-8) = HESS(3*NATOMS-7,3*NATOMS-8)+D2VDKDPHI
1558:         END IF1339:         END IF
1559: 1340: 
1560:     END SUBROUTINE CALC_HESS_PARAMS_PER1341:     END SUBROUTINE CALC_HESS_PARAMS
1561: 1342: 
1562: !-------------------------------------------------------------------------------1343: !-------------------------------------------------------------------------------
1563: !1344: !
1564: ! Repel the unit cell volume from approaching 0, which repels the unit cell from1345: ! Repel the unit cell volume from approaching 0, which repels the unit cell from
1565: ! adopting physically impossible angle combinations. Uses a WCA potential1346: ! adopting physically impossible angle combinations. Uses a WCA potential
1566: ! ENERGY: energy of the system1347: ! ENERGY: energy of the system
1567: ! GRAD_ANGLES: gradients for the box angles1348: ! GRAD_ANGLES: gradients for the box angles
1568: ! BD: information on the box and its derivatives1349: ! BD: information on the box and its derivatives
1569: ! GTEST: whether to calculate gradients1350: ! GTEST: whether to calculate gradients
1570: ! STEST: whether to calculate second derivatives1351: ! STEST: whether to calculate second derivatives
1644:         USE MODHESS, ONLY: HESS        1425:         USE MODHESS, ONLY: HESS        
1645: 1426: 
1646:         IMPLICIT NONE1427:         IMPLICIT NONE
1647: 1428: 
1648:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(3)1429:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(3)
1649:         LOGICAL, INTENT(IN) :: GTEST, STEST1430:         LOGICAL, INTENT(IN) :: GTEST, STEST
1650: 1431: 
1651:         ! EN_CON: constraint contribution to the energy1432:         ! EN_CON: constraint contribution to the energy
1652:         DOUBLE PRECISION :: EN_CON1433:         DOUBLE PRECISION :: EN_CON
1653: 1434: 
1654:         ! OFFSET: position of the end of the atom indices of the potential 
1655:         !         paramters 
1656:         INTEGER :: OFFSET 
1657:  
1658:         ! Initialise output variable1435:         ! Initialise output variable
1659:         EN_CON = 0.D01436:         EN_CON = 0.D0
1660: 1437: 
1661:         IF (OPP_MODE .EQ. 2) THEN 
1662:             OFFSET = - 6 
1663:         ELSE IF (OPP_MODE .EQ. 3) THEN 
1664:             OFFSET = 0 
1665:         END IF 
1666:  
1667:         ! OPP_K energy contribution1438:         ! OPP_K energy contribution
1668:         IF (BTEST(OPP_PARAMS, 0)) THEN1439:         IF (BTEST(OPP_PARAMS, 0)) THEN
1669:             EN_CON = EN_CON + 0.5D0 * K_REPULSION * &1440:             EN_CON = EN_CON + 0.5D0 * K_REPULSION * &
1670:                 (MERGE(OPP_K - OPP_K_UPPER, 0.D0, &1441:                 (MERGE(OPP_K - K_UPPER, 0.D0, &
1671:                        OPP_K .GT. OPP_K_UPPER)**2 + &1442:                        OPP_K .GT. K_UPPER)**2 + &
1672:                  MERGE(OPP_K - OPP_K_LOWER, 0.D0, &1443:                  MERGE(OPP_K - K_LOWER, 0.D0, &
1673:                        OPP_K .LT. OPP_K_LOWER)**2)1444:                        OPP_K .LT. K_LOWER)**2)
1674: 1445: 
1675:             IF (GTEST) THEN1446:             IF (GTEST) THEN
1676:                 ! OPP_K gradient contribution1447:                 ! OPP_K gradient contribution
1677:                 GRAD_PARAMS(1) = GRAD_PARAMS(1) + K_REPULSION * &1448:                 GRAD_PARAMS(1) = GRAD_PARAMS(1) + K_REPULSION * &
1678:                     (MERGE(OPP_K - OPP_K_UPPER, 0.D0, &1449:                     (MERGE(OPP_K - K_UPPER, 0.D0, &
1679:                            OPP_K .GT. OPP_K_UPPER) + &1450:                            OPP_K .GT. K_UPPER) + &
1680:                      MERGE(OPP_K - OPP_K_LOWER, 0.D0, &1451:                      MERGE(OPP_K - K_LOWER, 0.D0, &
1681:                            OPP_K .LT. OPP_K_LOWER))1452:                            OPP_K .LT. K_LOWER))
1682:             END IF1453:             END IF
1683: 1454: 
1684:             IF (STEST) THEN1455:             IF (STEST) THEN
1685:                 ! OPP_K Hessian contribution1456:                 ! OPP_K Hessian contribution
1686:                 HESS(3*NATOMS+OFFSET-2,3*NATOMS+OFFSET-2) = &1457:                 HESS(3*NATOMS-8,3*NATOMS-8) = HESS(3*NATOMS-8,3*NATOMS-8) + &
1687:                     HESS(3*NATOMS+OFFSET-2,3*NATOMS+OFFSET-2) + &1458:                     MERGE(K_REPULSION, 0.D0, OPP_K .GT. K_UPPER) + &
1688:                     MERGE(K_REPULSION, 0.D0, OPP_K .GT. OPP_K_UPPER) + &1459:                     MERGE(K_REPULSION, 0.D0, OPP_K .LT. K_LOWER)
1689:                     MERGE(K_REPULSION, 0.D0, OPP_K .LT. OPP_K_LOWER) 
1690:             END IF 1460:             END IF 
1691:         END IF1461:         END IF
1692: 1462: 
1693:         ! Add the energy contribution1463:         ! Add the energy contribution
1694:         ENERGY = ENERGY + EN_CON1464:         ENERGY = ENERGY + EN_CON
1695: 1465: 
1696:     END SUBROUTINE CONSTRAIN_PARAMETERS1466:     END SUBROUTINE CONSTRAIN_PARAMETERS
1697: 1467: 
1698: !-------------------------------------------------------------------------------1468: !-------------------------------------------------------------------------------
1699: !1469: !


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0