hdiff output

r32102/GMIN.tex 2017-03-16 16:30:12.096409735 +0000 r32101/GMIN.tex 2017-03-16 16:30:12.772418525 +0000
1106: The parameter {\it f} specifies what fraction of the Monte Carlo steps should be swaps between the1106: The parameter {\it f} specifies what fraction of the Monte Carlo steps should be swaps between the
1107: positions of a charged and a neutral particle, rather than a conventional step.1107: positions of a charged and a neutral particle, rather than a conventional step.
1108: $T_{\rm swap}$ is the temperature to be used in the acceptance criterion for swap moves, overriding1108: $T_{\rm swap}$ is the temperature to be used in the acceptance criterion for swap moves, overriding
1109: that specified using the {\it TEMPERATURE} keyword.  Generally, a lower temperature is more effective1109: that specified using the {\it TEMPERATURE} keyword.  Generally, a lower temperature is more effective
1110: at finding the lowest-energy permutation of charges.  The default value of $T_{\rm swap}$ is zero.1110: at finding the lowest-energy permutation of charges.  The default value of $T_{\rm swap}$ is zero.
1111: The reduced charge $q'$ is related to the actual charge $q$ by $q'=q/(4\pi\epsilon_0\epsilon\sigma)^{1/2}$,1111: The reduced charge $q'$ is related to the actual charge $q$ by $q'=q/(4\pi\epsilon_0\epsilon\sigma)^{1/2}$,
1112: where $\epsilon$ and $\sigma$ are the Lennard-Jones well depth and length parameter respectively.1112: where $\epsilon$ and $\sigma$ are the Lennard-Jones well depth and length parameter respectively.
1113: This way, the reduced energy of two charges is $E'=q'^2/r'$, where $r'=r/\sigma$ is the reduced distance1113: This way, the reduced energy of two charges is $E'=q'^2/r'$, where $r'=r/\sigma$ is the reduced distance
1114: bdetween the charges.1114: bdetween the charges.
1115: 1115: 
1116: \item {\it LJGAUSS mode rcut eps $r_0$} sepcifies that the Lennard-Jones Gauss double well potential should be used. 
1117: The Lennard-Jones potential has an energy minimum of 1 and a minimum at distance 1. 
1118: The relative weight of the Gaussian potential is {\it eps} and the Gaussian peak is at distance $r_0$. The Stoddard-Ford 
1119: procedure is used to make the potential and its first and second derivatives go smoothly to zero at {\it rcut}. 
1120: {\it mode} can be either 0, 1 or 2. 0 specifies the normal potential, and can be either a cluster or with periodic 
1121: boundary conditions with an orthogonal unit cell using the PERIODIC keyword. 1 or 2 requires periodic boundary conditions. 
1122: The PERIODIC keyword must be used, but the unit cell parameters from it will be ignored. Instead, the last line of the coords file will 
1123: be interpreted as the unit cell length parameters and will be optimised. Additionally, for {\it mode} 2, the second 
1124: last line of the coords file be interpreted as unit cell angle parameters for a triclinic unit cell, and optimised. Since there 
1125: are combinations of unit cell angles that are physically impossible, steps will be rejected if such combinations are detected. 
1126:  
1127: \item {\it LOCALSAMPLE abthresh acthresh\/}: Keyword currently under construction! For three groups of atoms defined in movableatoms1116: \item {\it LOCALSAMPLE abthresh acthresh\/}: Keyword currently under construction! For three groups of atoms defined in movableatoms
1128: (A,B,C), a step is quenched when the centre of coordinates of A->B is less than {\it abthresh} AND A->C is less than {\it acthresh}. 1117: (A,B,C), a step is quenched when the centre of coordinates of A->B is less than {\it abthresh} AND A->C is less than {\it acthresh}. 
1129: If this condition is broken AFTER the quench, it is automatically rejected.1118: If this condition is broken AFTER the quench, it is automatically rejected.
1130: 1119: 
1131: \item {\it LPERMDIST n thresh cut orbittol\/}: provides the preferred local1120: \item {\it LPERMDIST n thresh cut orbittol\/}: provides the preferred local
1132: permutational alignment procedure.1121: permutational alignment procedure.
1133: {\it n\/} is the maximum number of extra atoms to include in running1122: {\it n\/} is the maximum number of extra atoms to include in running
1134: {\bf minperm\/} for each permutable group (default 4), rather than1123: {\bf minperm\/} for each permutable group (default 4), rather than
1135: treating the whole system at the same time, as for {\it PERMDIST}.1124: treating the whole system at the same time, as for {\it PERMDIST}.
1136: The extra atoms that define the neighbourhood of the permutable group are1125: The extra atoms that define the neighbourhood of the permutable group are


r32102/keywords.f 2017-03-16 16:30:12.328412767 +0000 r32101/keywords.f 2017-03-16 16:30:13.000421481 +0000
6188:          LJCOULT=.TRUE.6188:          LJCOULT=.TRUE.
6189:          CALL READI(COULN)6189:          CALL READI(COULN)
6190:          CALL READF(COULQ)6190:          CALL READF(COULQ)
6191:          CALL READF(COULSWAP)6191:          CALL READF(COULSWAP)
6192:          CALL READF(COULTEMP)6192:          CALL READF(COULTEMP)
6193:          NRBSITES=16193:          NRBSITES=1
6194:          ALLOCATE(SITE(NRBSITES,3))6194:          ALLOCATE(SITE(NRBSITES,3))
6195: !        Maybe the above two lines are not necessary!6195: !        Maybe the above two lines are not necessary!
6196: 6196: 
6197: ! jwrm2>6197: ! jwrm2>
6198:       ELSE IF (WORD.EQ.'LJGAUSS') THEN6198:       ELSE IF (WORD.EQ.'LJ_GAUSS') THEN
6199:           LJ_GAUSST = .TRUE.6199:           LJ_GAUSST = .TRUE.
6200:           CALL READI(LJ_GAUSS_MODE)6200:           CALL READI(LJ_GAUSS_MODE)
6201:           CALL READF(LJ_GAUSS_RCUT)6201:           CALL READF(LJ_GAUSS_RCUT)
6202:           CALL READF(LJ_GAUSS_EPS)6202:           CALL READF(LJ_GAUSS_EPS)
6203:           CALL READF(LJ_GAUSS_R0)6203:           CALL READF(LJ_GAUSS_R0)
6204:           CALL LJ_GAUSS_INITIALISE()6204:           CALL LJ_GAUSS_INITIALISE()
6205: ! jwrm2> end6205: ! jwrm2> end
6206: 6206: 
6207:       ELSE IF (WORD.EQ.'STOCK') THEN6207:       ELSE IF (WORD.EQ.'STOCK') THEN
6208:          STOCKT=.TRUE.6208:          STOCKT=.TRUE.


r32102/lj_gauss.f90 2017-03-16 16:30:12.548415620 +0000 r32101/lj_gauss.f90 2017-03-16 16:30:13.220424336 +0000
 18: ! 18: !
 19: ! 19: !
 20: ! This module provides a routine for calculating the energy and gradients with 20: ! This module provides a routine for calculating the energy and gradients with
 21: ! a double well Lennard-Jones + Gaussian potential. 21: ! a double well Lennard-Jones + Gaussian potential.
 22: ! Questions to jwrm2. 22: ! Questions to jwrm2.
 23: ! 23: !
 24: MODULE LJ_GAUSS_MOD 24: MODULE LJ_GAUSS_MOD
 25:  25: 
 26: PUBLIC :: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, LJ_GAUSS_R0 26: PUBLIC :: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, LJ_GAUSS_R0
 27: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS 27: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS
 28: PRIVATE :: LJ_GAUSS_SIGMASQ, RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT, PI 28: PRIVATE :: LJ_GAUSS_SIGMASQ, RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT
 29: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_GRAD 29: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_GRAD
 30: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES 
 31:  30: 
 32: ! EPS: strength of the Gaussian potential. The LJ strength is 1. 31: ! EPS: strength of the Gaussian potential. The LJ strength is 1.
 33: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ 32: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ
 34: !     rmin is 1. 33: !     rmin is 1.
 35: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of 34: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of
 36: !         both parts of the potential go smoothly to zero. 35: !         both parts of the potential go smoothly to zero.
 37: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT 36: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT
 38:  37: 
 39: ! LJ_GAUSS_MODE: 0 standard minimisation 38: ! LJ_GAUSS_MODE: 0 standard minimisation
 40: !                1 the final triplet of coordinates will be interpreted as the 39: !                1 the final triplet of coordinates will be interpreted as the
 45: ! LJ_GAUSS_SIGMA: variation (width) of the Gaussian potential 44: ! LJ_GAUSS_SIGMA: variation (width) of the Gaussian potential
 46: DOUBLE PRECISION, PARAMETER :: LJ_GAUSS_SIGMASQ = 0.02D0 45: DOUBLE PRECISION, PARAMETER :: LJ_GAUSS_SIGMASQ = 0.02D0
 47:  46: 
 48: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance 47: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance
 49: DOUBLE PRECISION :: RCUTINV6, RCUTINV12 48: DOUBLE PRECISION :: RCUTINV6, RCUTINV12
 50:  49: 
 51: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle 50: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle
 52: !                      distance 51: !                      distance
 53: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT 52: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT
 54:  53: 
 55: ! PI: the usual meaning 
 56: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0) 
 57:  
 58: ! BOX_DERIV: stores information for calculating lattice derivatives 
 59: TYPE :: BOX_DERIV 
 60:     ! ORTHOG: the orthogonal tansformation matrix, for changing from 
 61:     !         crystallographic space to orthogonal space first index is row, 
 62:     !         second is column 
 63:     DOUBLE PRECISION :: ORTHOG(3,3) 
 64:  
 65:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters, 
 66:     !        first index indicates the row, second index the column and 
 67:     !        third index the parameter, with 1-3 being the angles and 4-6 being 
 68:     !        the lengths 
 69:     DOUBLE PRECISION :: DERIV(3, 3, 6) 
 70:  
 71:     ! RECIP: lengths of the reciprocal lattice vectors 
 72:     DOUBLE PRECISION :: RECIP(3) 
 73: END TYPE BOX_DERIV 
 74:  
 75: CONTAINS 54: CONTAINS
 76:  55: 
 77: !------------------------------------------------------------------------------- 56: !-------------------------------------------------------------------------------
 78: ! 57: !
 79: ! Main routine for calculating the energy and gradients. 58: ! Main routine for calculating the energy and gradients.
 80: ! X: coordinates array 59: ! X: coordinates array
 81: ! GRAD: gradients array 60: ! GRAD: gradients array
 82: ! ENERGY: calculated energy 61: ! ENERGY: calculated energy
 83: ! GTEST: whether gradients should be calculated 62: ! GTEST: whether gradients should be calculated
 84: ! 63: !
 85: !------------------------------------------------------------------------------- 64: !-------------------------------------------------------------------------------
 86:     SUBROUTINE LJ_GAUSS(X, GRAD, ENERGY, GTEST) 65:     SUBROUTINE LJ_GAUSS(X, GRAD, ENERGY, GTEST)
 87:  66: 
 88:         ! MYUNIT: file unit for main output file 
 89:         ! NATOMS: number of particles 67:         ! NATOMS: number of particles
 90:         ! PERIODIC: whether to use periodic boundary conditions 68:         ! PERIODIC: whether to use periodic boundary conditions
 91:         USE COMMONS, ONLY: MYUNIT, NATOMS, PERIODIC 69:         USE COMMONS, ONLY: NATOMS, PERIODIC
 92:  70: 
 93:         IMPLICIT NONE 71:         IMPLICIT NONE
 94:  72: 
 95:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS) 73:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS)
 96:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY 74:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY
 97:         LOGICAL, INTENT(IN) :: GTEST 75:         LOGICAL, INTENT(IN) :: GTEST
 98:  76: 
 99:         ! I, J: indices for looping over particles 77:         ! I, J: indices for looping over particles
100:         ! NMOL: actual number of particles 78:         ! NMOL: actual number of particles
101:         INTEGER :: I, J, NMOL 79:         INTEGER :: I, J, NMOL
102:  80: 
103:         ! RIJ: Interparticle vector 81:         ! RIJ: Interparticle vector
104:         ! MODRIJ: distance between a pair of particles 82:         ! MODRIJ: distance between a pair of particles
105:         ! DVDR: derivative of pair potential wrt MODRIJ 83:         ! DVDR: derivative of pair potential wrt MODRIJ
106:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR 84:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR
107:  85: 
108:         ! Factors for triclinic lattice 
109:         TYPE(BOX_DERIV) :: BD 
110:  
111:         ! Initialise output variables 86:         ! Initialise output variables
112:         ENERGY = 0.D0 87:         ENERGY = 0.D0
113:         GRAD(:) = 0.D0 88:         GRAD(:) = 0.D0
114:  89: 
115:         IF (.NOT. PERIODIC) THEN 90:         IF (.NOT. PERIODIC .AND. LJ_GAUSS_MODE .EQ. 0) THEN
116:             IF (LJ_GAUSS_MODE .EQ. 0) THEN 91:             ! Normal cluster
117:                 ! Normal cluster 92:             NMOL = NATOMS
118:                 NMOL = NATOMS 93:             DO I = 1, NMOL-1 ! Outer loop over particles
119:                 DO I = 1, NMOL-1 ! Outer loop over particles 94:                 DO J = I+1, NMOL ! Inner loop over particles
120:                     DO J = I+1, NMOL ! Inner loop over particles 95:                     ! Get the particle-particle distance
121:                         ! Get the particle-particle distance 96:                     RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)
122:                         RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J) 97:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
123:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
124:  98: 
125:                         ! Add energy and gradients 99:                     ! Add energy and gradients
126:                         IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN100:                     IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN
127:                             ENERGY = ENERGY + CALC_ENERGY(MODRIJ)101:                         ENERGY = ENERGY + CALC_ENERGY(MODRIJ)
128: 102: 
129:                             IF (GTEST) THEN103:                         IF (GTEST) THEN
130:                                 DVDR = CALC_GRAD(MODRIJ)104:                             DVDR = CALC_GRAD(MODRIJ)
131:                                 GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + &105:                             GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + &
132:                                                   DVDR*RIJ(:)/MODRIJ106:                                               DVDR*RIJ(:)/MODRIJ
133:                                 GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - &107:                             GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - &
134:                                                   DVDR*RIJ(:)/MODRIJ108:                                               DVDR*RIJ(:)/MODRIJ
135:                             END IF ! GTEST109:                         END IF ! GTEST
136:                         END IF ! If less than cutoff distance110:                     END IF ! If less than cutoff distance
137:                     END DO ! Inner loop over particles111:                 END DO ! Inner loop over particles
138:                 END DO ! Outer loop over particles112:             END DO ! Outer loop over particles
139:             ELSE ! mode not equal to zero113: 
140:                 WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR, PERIODIC must be ', &114:         ELSE IF (PERIODIC .AND. LJ_GAUSS_MODE .EQ. 0) THEN
141:                                   'specified if mode is not 0'115:             ! Periodic, fixed box
142:                 STOP116:             ! Reset all particles to within the box
143:             END IF117:             CALL PERIODIC_RESET(X(:))
144:         ELSE ! PERIODIC 
145:             SELECT CASE (LJ_GAUSS_MODE) 
146:             CASE(0) 
147:                 ! Periodic, fixed box 
148:                 ! Reset all particles to within the box 
149:                 CALL PERIODIC_RESET(X(:)) 
150:  
151:                 NMOL = NATOMS 
152:                 ! We need to include self-interactions of particles in different 
153:                 ! unit cells 
154:                 DO I = 1, NMOL ! Outer loop over particles 
155:                     DO J = I, NMOL ! Inner loop over particles 
156:                         ! Add energy and gradients 
157:                         CALL LJ_GAUSS_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), & 
158:                                           GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), & 
159:                                           ENERGY, GTEST) 
160:                     END DO ! Inner loop over particles 
161:                 END DO ! Outer loop over particles 
162:  
163:             CASE(1) 
164:                 ! Periodic varying lattice 
165:                 ! Reset all particles to within the box 
166:                 CALL PERIODIC_RESET(X(:)) 
167:  
168:                 NMOL = NATOMS-1 
169:                 ! We need to include self-interactions of particles in different 
170:                 ! unit cells 
171:                 DO I = 1, NMOL! Outer loop over particles 
172:                     DO J = I, NMOL ! Inner loop over particles 
173:                         ! Add energy and gradients 
174:                         CALL LJ_GAUSS_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), & 
175:                                           GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), & 
176:                                           ENERGY, GTEST, & 
177:                                           GRAD(3*NATOMS-2:3*NATOMS)) 
178:                     END DO ! Inner loop over particles 
179:                 END DO ! Outer loop over particles 
180:  
181:             CASE(2) 
182:                 ! Triclinic varying lattice 
183:                 ! Reset all particles to within the box 
184:                 CALL PERIODIC_RESET(X(:)) 
185:  
186:                 ! Check whether the unit cell angles are physically possible. 
187:                 ! If we have a disallowed unit cell, we set the energy to very 
188:                 ! high and the gradients to very small, so the step is 
189:                 ! 'converged' but gets rejected. 
190:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN 
191:                     ENERGY = 1.D20 
192:                     GRAD(:) = 1.D-20 
193:                     RETURN 
194:                 END IF 
195: 118: 
196:                 ! Calculate box derivative factors119:             NMOL = NATOMS
197:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &120:             ! We need to include self-interactions of particles in different
198:                                     X(3*NATOMS-2:3*NATOMS  ))121:             ! unit cells
199: 122:             DO I = 1, NMOL ! Outer loop over particles
200:                 NMOL = NATOMS-2123:                 DO J = I, NMOL ! Inner loop over particles
201:                 ! We need to include self-interactions of particles in different124:                     ! Add energy and gradients
202:                 ! unit cells125:                     CALL LJ_GAUSS_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &
203:                 DO I = 1, NMOL! Outer loop over particles126:                                       GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), &
204:                     DO J = I, NMOL ! Inner loop over particles127:                                       ENERGY, GTEST)
205:                         ! Add energy and gradients128:                 END DO ! Inner loop over particles
206:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &129:             END DO ! Outer loop over particles
207:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &130: 
208:                             GRAD(3*NATOMS-5:3*NATOMS), BD)131:         ELSE IF (PERIODIC .AND. LJ_GAUSS_MODE .EQ. 1) THEN
209:                     END DO ! Inner loop over particles132:             ! Periodic varying lattice
210:                 END DO ! Outer loop over particles133:             ! Reset all particles to within the box
211:             END SELECT ! Mode selections134:             CALL PERIODIC_RESET(X(:))
212:         END IF ! Periodic135: 
 136:             NMOL = NATOMS-1
 137:             ! We need to include self-interactions of particles in different
 138:             ! unit cells
 139:             DO I = 1, NMOL! Outer loop over particles
 140:                 DO J = I, NMOL ! Inner loop over particles
 141:                     ! Add energy and gradients
 142:                     CALL LJ_GAUSS_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &
 143:                                       GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), &
 144:                                       ENERGY, GTEST, GRAD(3*NATOMS-2:3*NATOMS))
 145:                 END DO ! Inner loop over particles
 146:             END DO ! Outer loop over particles
 147:         END IF ! Mode selection
213: 148: 
214:     END SUBROUTINE LJ_GAUSS149:     END SUBROUTINE LJ_GAUSS
215: 150: 
216: !-------------------------------------------------------------------------------151: !-------------------------------------------------------------------------------
217: !152: !
218: ! Calculates the energy contribution for a pair of particles.153: ! Calculates the energy contribution for a pair of particles.
219: ! MODRIJ: distance beteen the particle pair154: ! MODRIJ: distance beteen the particle pair
220: !155: !
221: !-------------------------------------------------------------------------------156: !-------------------------------------------------------------------------------
222:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)157:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)
332: !267: !
333: !-------------------------------------------------------------------------------268: !-------------------------------------------------------------------------------
334:     SUBROUTINE LJ_GAUSS_INITIALISE ()269:     SUBROUTINE LJ_GAUSS_INITIALISE ()
335: 270: 
336:         ! MYUNIT: file unit for main GMIN output271:         ! MYUNIT: file unit for main GMIN output
337:         USE COMMONS, ONLY: MYUNIT272:         USE COMMONS, ONLY: MYUNIT
338: 273: 
339:         IMPLICIT NONE274:         IMPLICIT NONE
340: 275: 
341:         ! Check we have sensible value for the mode276:         ! Check we have sensible value for the mode
342:         IF (LJ_GAUSS_MODE .NE. 1 .AND. LJ_GAUSS_MODE .NE. 0 .AND. &277:         IF (LJ_GAUSS_MODE .NE. 1 .AND. LJ_GAUSS_MODE .NE. 0) THEN
343:             LJ_GAUSS_MODE .NE. 2) THEN278:             WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: mode must be 0 or 1'
344:             WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: mode must be 0, 1 or 2' 
345:             STOP279:             STOP
346:         END IF280:         END IF
347: 281: 
348:         ! Calculate some factors for LJ that don't depend on RIJ282:         ! Calculate some factors for LJ that don't depend on RIJ
349:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6283:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6
350:         RCUTINV12 = RCUTINV6**2284:         RCUTINV12 = RCUTINV6**2
351: 285: 
352:         ! Calculate some factors for Gauss that don't depend on RIJ286:         ! Calculate some factors for Gauss that don't depend on RIJ
353:         EXP_RCUT = LJ_GAUSS_EPS * EXP(- (LJ_GAUSS_RCUT - LJ_GAUSS_R0)**2 / &287:         EXP_RCUT = LJ_GAUSS_EPS * EXP(- (LJ_GAUSS_RCUT - LJ_GAUSS_R0)**2 / &
354:                                         (2.D0*LJ_GAUSS_SIGMASQ))288:                                         (2.D0*LJ_GAUSS_SIGMASQ))
368: 302: 
369:         ! BOXLX, BOXLY, BOXLZ: dimensions of box303:         ! BOXLX, BOXLY, BOXLZ: dimensions of box
370:         ! NATOMS: number of particles304:         ! NATOMS: number of particles
371:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ, NATOMS305:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ, NATOMS
372: 306: 
373:         IMPLICIT NONE307:         IMPLICIT NONE
374: 308: 
375:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)309:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)
376: 310: 
377:         ! Sort out the unit cell sizes, if they are varying311:         ! Sort out the unit cell sizes, if they are varying
378:         IF (LJ_GAUSS_MODE .EQ. 1 .OR. LJ_GAUSS_MODE .EQ. 2) THEN312:         IF (LJ_GAUSS_MODE .EQ. 1) THEN
379:             BOXLX = COORDS(3*NATOMS - 2)313:             BOXLX = COORDS(3*NATOMS - 2)
380:             BOXLY = COORDS(3*NATOMS - 1)314:             BOXLY = COORDS(3*NATOMS - 1)
381:             BOXLZ = COORDS(3*NATOMS)315:             BOXLZ = COORDS(3*NATOMS)
382:         END IF316:         END IF
383: 317: 
384:         ! Reset coordinates to within the box318:         ! Reset coordinates to within the box
385:         IF (LJ_GAUSS_MODE .EQ. 0) THEN319:         IF (LJ_GAUSS_MODE .EQ. 0) THEN
386:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:))320:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:))
387:         ELSE IF (LJ_GAUSS_MODE .EQ. 1) THEN321:         ELSE IF (LJ_GAUSS_MODE .EQ. 1) THEN
388:             COORDS(1:3*NATOMS-3) = COORDS(1:3*NATOMS-3) - &322:             COORDS(1:3*NATOMS-3) = COORDS(1:3*NATOMS-3) - &
389:                                    FLOOR(COORDS(1:3*NATOMS-3))323:                                    FLOOR(COORDS(1:3*NATOMS-3))
390:         ELSE IF (LJ_GAUSS_MODE .EQ. 2) THEN 
391:             COORDS(1:3*NATOMS-6) = COORDS(1:3*NATOMS-6) - & 
392:                                    FLOOR(COORDS(1:3*NATOMS-6)) 
393:             ! Make the lattice angles modulo 2*pi 
394:             ! They will be checked later 
395:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - & 
396:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI)) 
397:         END IF324:         END IF
398: 325: 
399:     END SUBROUTINE PERIODIC_RESET326:     END SUBROUTINE PERIODIC_RESET
400: 327: 
401: !-------------------------------------------------------------------------------328: !-------------------------------------------------------------------------------
402: !329: !
403: ! Finds the energy and gradients for the periodic system330: ! Finds the energy and gradients for the periodic system
404: ! IDI, IDJ: id numbers of the particles331: ! ID1, ID2: id numbers of the particles
405: ! COORDS_I: coordinates of particle I332: ! COORDS_I: coordinates of particle I
406: ! COORDS_J: coordinates of particle J333: ! COORDS_J: coordinates of particle J
407: ! GRAD_I: gradients for particle I334: ! GRAD_I: gradients for particle I
408: ! GRAD_J: gradients for particle J335: ! GRAD_J: gradients for particle J
409: ! ENERGY: energy contribution for this pair336: ! ENERGY: energy contribution for this pair
410: ! RIJ: mininum length vector between the particles337: ! RIJ: mininum length vector between the particles
411: ! GTEST: whether to calculate gradients338: ! GTEST: whether to calculate gradients
412: ! GRAD_BOX: gradients for the box339: ! GRAD_BOX: gradients for the box
413: !340: !
414: !-------------------------------------------------------------------------------341: !-------------------------------------------------------------------------------
415: 342: 
416:     SUBROUTINE LJ_GAUSS_PER(IDI, IDJ, COORDS_I,COORDS_J,GRAD_I,GRAD_J,ENERGY, &343:     SUBROUTINE LJ_GAUSS_PER(ID1, ID2, COORDS_I,COORDS_J,GRAD_I,GRAD_J,ENERGY, &
417:                             GTEST, GRAD_BOX)344:                             GTEST, GRAD_BOX)
418: 345: 
419:         ! BOXLX, BOXLY, BOXLZ: dimensions of the box346:         ! BOXLX, BOXLY, BOXLZ: dimensions of the box
420:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ347:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ
421: 348: 
422:         IMPLICIT NONE349:         INTEGER, INTENT(IN) :: ID1, ID2
423:  
424:         INTEGER, INTENT(IN) :: IDI, IDJ 
425:         DOUBLE PRECISION, INTENT(IN) :: COORDS_I(3), COORDS_J(3)350:         DOUBLE PRECISION, INTENT(IN) :: COORDS_I(3), COORDS_J(3)
426:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY351:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY
427:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_BOX(3)352:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_BOX(3)
428:         LOGICAL, INTENT(IN) :: GTEST353:         LOGICAL, INTENT(IN) :: GTEST
429: 354: 
430:         ! XJ: copy of j coordinates355:         ! XJ: copy of j coordinates
431:         ! RIJ: particle-particle vector356:         ! RIJ(3): particle-particle vector
432:         ! MODRIJ: particle-particle distance357:         ! MODRIJ: particle-particle distance
433:         ! DVDR: derivative of pair potential wrt MODRIJ358:         ! DVDR: derivative of pair potential wrt MODRIJ
434:         ! BOX_SIZE: convenience copy of the box size359:         ! BOX_SIZE: convenience copy of the box size
435:         DOUBLE PRECISION :: XJ(3), RIJ(3), MODRIJ, DVDR, BOX_SIZE(3)360:         DOUBLE PRECISION :: XJ(3), RIJ(3), MODRIJ, DVDR, BOX_SIZE(3)
436: 361: 
437:         ! PER_k: integers for looping over cell possibilities362:         ! PER_k: integers for looping over cell possibilities
438:         ! CELL_RANGE: number of cells in each direction required363:         ! CELL_RANGE: number of cells in each direction required
439:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3)364:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3)
440: 365: 
441: !        DOUBLE PRECISION :: START_ENERGY366: !        DOUBLE PRECISION :: START_ENERGY
451: 376: 
452:         ! Loop over different coordinate possibilities.377:         ! Loop over different coordinate possibilities.
453:         DO PER_X = - CELL_RANGE(1), CELL_RANGE(1)378:         DO PER_X = - CELL_RANGE(1), CELL_RANGE(1)
454:             DO PER_Y = - CELL_RANGE(2), CELL_RANGE(2)379:             DO PER_Y = - CELL_RANGE(2), CELL_RANGE(2)
455:                 DO PER_Z = - CELL_RANGE(3), CELL_RANGE(3)380:                 DO PER_Z = - CELL_RANGE(3), CELL_RANGE(3)
456: 381: 
457:                     ! Skip the interaction of a particle with itself in the same382:                     ! Skip the interaction of a particle with itself in the same
458:                     ! cell and don't double count interactions within the unit383:                     ! cell and don't double count interactions within the unit
459:                     ! cell.384:                     ! cell.
460:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 &385:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 &
461:                         .AND. IDI .EQ. IDJ) CYCLE386:                         .AND. ID1 .EQ. ID2) CYCLE
462: 387: 
463:                     ! Adjust the j coordinates for the periodic image388:                     ! Adjust the j coordinates for the periodic image
464:                     XJ(1) = COORDS_J(1) + PER_X389:                     XJ(1) = COORDS_J(1) + PER_X
465:                     XJ(2) = COORDS_J(2) + PER_Y390:                     XJ(2) = COORDS_J(2) + PER_Y
466:                     XJ(3) = COORDS_J(3) + PER_Z391:                     XJ(3) = COORDS_J(3) + PER_Z
467: 392: 
468:                     ! Find the distance and scale by the box size393:                     ! Find the distance and scale by the box size
469:                     RIJ(:) = (COORDS_I(:) - XJ(:))*BOX_SIZE(:)394:                     RIJ(:) = (COORDS_I(:) - XJ(:))*BOX_SIZE(:)
470:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))395:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
471: 396: 
472:                     IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN397:                     IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN
473:                         ! Add energy and gradients398:                         ! Add energy and gradients
474:                         ! Divide by 2 for images of the same atom, to avoid399:                         ! Divide by 2 for images of the same atom, to avoid
475:                         ! double counting400:                         ! double counting
476:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &401:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, ID1 .EQ. ID2) * &
477:                                           CALC_ENERGY(MODRIJ)402:                                           CALC_ENERGY(MODRIJ)
478: 403: 
479:                         IF (GTEST) THEN404:                         IF (GTEST) THEN
480:                             ! Divide gradient by 2 for images of the same atom,405:                             ! Divide gradient by 2 for images of the same atom,
481:                             ! to avoid double counting406:                             ! to avoid double counting
482:                             DVDR = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * &407:                             DVDR = MERGE(0.5D0, 1.0D0, ID1 .EQ. ID2) * &
483:                                    CALC_GRAD(MODRIJ)408:                                    CALC_GRAD(MODRIJ)
484:                             ! We must undo the box size coordinate scaling409:                             ! We must undo the box size coordinate scaling
485:                             GRAD_I(:) = GRAD_I(:) + &410:                             GRAD_I(:) = GRAD_I(:) + &
486:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ411:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ
487:                             GRAD_J(:) = GRAD_J(:) - &412:                             GRAD_J(:) = GRAD_J(:) - &
488:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ413:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ
489: 414: 
490:                             IF (PRESENT(GRAD_BOX)) THEN415:                             IF (PRESENT(GRAD_BOX)) THEN
491:                                 ! Lattice derivatives416:                                 ! Lattice derivatives
492:                                 GRAD_BOX(:) = GRAD_BOX(:) + &417:                                 GRAD_BOX(:) = GRAD_BOX(:) + &
495:                         END IF ! GTEST420:                         END IF ! GTEST
496:                     END IF ! End if less than cutoff421:                     END IF ! End if less than cutoff
497:                 END DO ! z loop422:                 END DO ! z loop
498:             END DO ! y loop423:             END DO ! y loop
499:         END DO ! x loop424:         END DO ! x loop
500: 425: 
501:     END SUBROUTINE LJ_GAUSS_PER426:     END SUBROUTINE LJ_GAUSS_PER
502: 427: 
503: !-------------------------------------------------------------------------------428: !-------------------------------------------------------------------------------
504: !429: !
505: ! Calculates factors for the lattice derivatives which only depend on the cell 
506: ! parameters. 
507: ! ANGLES: alpha, beta and gamma lattice parameters 
508: ! SIZES: a, b and c lattice parameters 
509: ! 
510: !------------------------------------------------------------------------------- 
511:  
512:     PURE TYPE(BOX_DERIV) FUNCTION CALC_BOX_DERIV(ANGLES, SIZES) 
513:  
514:         IMPLICIT NONE 
515:  
516:         DOUBLE PRECISION, INTENT(IN) :: ANGLES(3), SIZES(3) 
517:  
518:         ! V: factor, related to (but not quite equal to) the unit cell volume 
519:         ! C, S: cos and sin of the lattice angles 
520:         ! ORTHOG, DERIV: parts of the output 
521:         DOUBLE PRECISION :: V, C(3), S(3) 
522:         DOUBLE PRECISION :: ORTHOG(3, 3), DERIV(3, 3, 6) 
523:  
524:         ! Initialise output variables 
525:         ORTHOG(:,:) = 0.D0 
526:         DERIV(:,:,:) = 0.D0 
527:  
528:         ! Calculate factors 
529:         C(:) = COS(ANGLES(:)) 
530:         S(:) = SIN(ANGLES(:)) 
531:         V = SQRT(1 - C(1)**2 - C(2)**2 - C(3)**2 + 2.D0 * C(1) * C(2) * C(3)) 
532:  
533:         ! Calculate the orthogonal transformation matrix 
534:         ORTHOG(1, 1) = SIZES(1) 
535:         ORTHOG(1, 2) = SIZES(2) * C(3) 
536:         ORTHOG(1, 3) = SIZES(3) * C(2) 
537:         ORTHOG(2, 2) = SIZES(2) * S(3) 
538:         ORTHOG(2, 3) = SIZES(3) * (C(1) - C(2) * C(3)) / S(3) 
539:         ORTHOG(3, 3) = SIZES(3) * V / S(3) 
540:  
541:         ! Calculate the derivatives of the othogonal matrix 
542:         ! Derivatives of the top row wrt angles 
543:         DERIV(1, 3, 2) = - SIZES(3) * S(2) 
544:         DERIV(1, 2, 3) = - SIZES(2) * S(3) 
545:         ! Derivatives of the top row wrt to lengths 
546:         DERIV(1, 1, 4) = 1.D0 
547:         DERIV(1, 2, 5) = C(3) 
548:         DERIV(1, 3, 6) = C(2) 
549:         ! Derivatives of the middle row wrt angles 
550:         DERIV(2, 3, 1) = - SIZES(3) * S(1) / S(3) 
551:         DERIV(2, 3, 2) = SIZES(3) * S(2) * C(3) / S(3) 
552:         DERIV(2, 2, 3) = SIZES(2) * C(3) 
553:         DERIV(2, 3, 3) = SIZES(3) * (C(2) - C(1) * C(3)) / S(3)**2 
554:         ! Derivatives of the middle row wrt to lengths 
555:         DERIV(2, 2, 5) = S(3) 
556:         DERIV(2, 3, 6) = (C(1) - C(2) * C(3)) / S(3) 
557:         ! Derivatives of the bottom row wrt angles 
558:         DERIV(3, 3, 1) = SIZES(3) * S(1) * (C(1) - C(2) * C(3)) / (S(3) * V) 
559:         DERIV(3, 3, 2) = SIZES(3) * S(2) * (C(2) - C(3) * C(1)) / (S(3) * V) 
560:         DERIV(3, 3, 3) = SIZES(3) * ((C(3) - C(1) * C(2))/V - C(3) * V/S(3)**2) 
561:         ! Derivatives of the bottom row wrt lengths 
562:         DERIV(3, 3, 6) = V / S(3) 
563:  
564:         ! Calculate the lengths of the reciprocal lattice vectors 
565:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V) 
566:  
567:         ! Set the output 
568:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:) 
569:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:) 
570:  
571:     END FUNCTION CALC_BOX_DERIV 
572:  
573: !------------------------------------------------------------------------------- 
574: ! 
575: ! Checks whether the unit cell angles are allowed. Itis possible to generate 
576: ! sets of angles which do not correspond to a physically possible unit cell. 
577: ! This would manifest by an imaginary unit cell volume. See 
578: ! https://journals.iucr.org/a/issues/2011/01/00/au5114/au5114.pdf 
579: ! 'On the allowed values for the triclinic unit-cell angles', 
580: ! J. Foadi and G. Evans, Acta Cryst. (2011) A67, 93–95 
581: ! ANGLES: the three unit cell angles 
582: ! 
583: !------------------------------------------------------------------------------- 
584:  
585:     PURE LOGICAL FUNCTION CHECK_ANGLES(ANGLES) 
586:  
587:         IMPLICIT NONE 
588:  
589:         DOUBLE PRECISION, INTENT(IN) :: ANGLES(3) 
590:  
591:         ! SUMS: sums and differences of the angles 
592:         DOUBLE PRECISION :: SUMS(4) 
593:  
594:         ! Calculate the necessary sums 
595:         SUMS(1) =   ANGLES(1) + ANGLES(2) + ANGLES(3) 
596:         SUMS(2) = - ANGLES(1) + ANGLES(2) + ANGLES(3) 
597:         SUMS(3) =   ANGLES(1) - ANGLES(2) + ANGLES(3) 
598:         SUMS(4) =   ANGLES(1) + ANGLES(2) - ANGLES(3) 
599:  
600:         ! Check all sums 0 < SUM < 2 pi and all angles 0 < ANGLE < pi 
601:         CHECK_ANGLES = ALL(SUMS .GT. 0.D0) .AND. ALL(SUMS .LT. 2*PI) .AND. & 
602:                        ALL(ANGLES .GT. 0.D0) .AND. ALL(ANGLES .LT. PI) 
603:  
604:     END FUNCTION CHECK_ANGLES 
605:  
606: !------------------------------------------------------------------------------- 
607: ! 
608: ! Finds the energy and gradients for the periodic system with a triclinic 
609: ! varying unit cell 
610: ! IDI, IDJ: id numbers of the particles 
611: ! COORDS: coordinates of particles (in crystallographic space) 
612: !         and lattice parameters 
613: ! GRAD_I: gradients for particle I 
614: ! GRAD_J: gradients for particle J 
615: ! ENERGY: energy contribution for this pair 
616: ! RIJ: mininum length vector between the particles 
617: ! GTEST: whether to calculate gradients 
618: ! GRAD_BOX: gradients for the box (angles, then lengths) 
619: ! BD: factors for the lattice derivatives 
620: ! 
621: !------------------------------------------------------------------------------- 
622:  
623:     SUBROUTINE LJ_GAUSS_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, & 
624:                                 GTEST, GRAD_BOX, BD) 
625:  
626:         ! NATOMS: number of coordinates (including lattice parameters) 
627:         USE COMMONS, ONLY: NATOMS 
628:  
629:         IMPLICIT NONE 
630:  
631:         INTEGER, INTENT(IN) :: IDI, IDJ 
632:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS) 
633:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY 
634:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6) 
635:         LOGICAL, INTENT(IN) :: GTEST 
636:         TYPE(BOX_DERIV), INTENT(IN) :: BD 
637:  
638:         ! RJ: particle coordinates in crystallographic space 
639:         ! RIJ: particle-particle vector in crystallographic space 
640:         ! YIJ: particle-particle vector in orthogonal space 
641:         ! MODYIJ: particle-particle distance in orthogonal space 
642:         ! DVDY: derivative of pair potential wrt MODYIJ 
643:         ! TEMP_GRAD: temporary gradient store 
644:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ 
645:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(6) 
646:  
647:         ! PER_k: integers for looping over cell possibilities 
648:         ! CELL_RANGE: number of cells in each direction required 
649:         ! I, J: loop indices 
650:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3) 
651:         INTEGER :: I, J 
652:  
653:         ! Set the number of cells we need to check, based on the cut off 
654:         ! distance and the distance between lattice planes in the directions of 
655:         ! the reciprocal lattice vectors. 
656:         CELL_RANGE(:) = CEILING(LJ_GAUSS_RCUT*BD%RECIP(:)) 
657:  
658:         ! Loop over different coordinate possibilities. 
659:         DO PER_X = - CELL_RANGE(1), CELL_RANGE(1) 
660:             DO PER_Y = - CELL_RANGE(2), CELL_RANGE(2) 
661:                 DO PER_Z = - CELL_RANGE(3), CELL_RANGE(3) 
662:  
663:                     ! Skip the interaction of a particle with itself in the same 
664:                     ! cell and don't double count interactions within the unit 
665:                     ! cell. 
666:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 & 
667:                         .AND. IDI .EQ. IDJ) CYCLE 
668:  
669:                     ! Adjust the j coordinates for the periodic image 
670:                     RJ(1) = COORDS(3*IDJ-2) + PER_X 
671:                     RJ(2) = COORDS(3*IDJ-1) + PER_Y 
672:                     RJ(3) = COORDS(3*IDJ  ) + PER_Z 
673:  
674:                     ! Find the distance in crystallographic space 
675:                     RIJ(:) = (COORDS(3*IDI-2:3*IDI) - RJ(:)) 
676:  
677:                     ! Find the distance in orthogonal space 
678:                     YIJ(:) = MATMUL(BD%ORTHOG(:,:), RIJ(:)) 
679:                     MODYIJ = SQRT(DOT_PRODUCT(YIJ(:), YIJ(:))) 
680:  
681:                     IF (MODYIJ .LT. LJ_GAUSS_RCUT) THEN 
682:                         ! Add energy and gradients 
683:                         ! Divide by 2 for images of the same atom, to avoid 
684:                         ! double counting 
685:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
686:                                           CALC_ENERGY(MODYIJ) 
687:  
688:                         IF (GTEST) THEN 
689:                             ! Divide gradient by 2 for images of the same atom, 
690:                             ! to avoid double counting 
691:                             DVDY = MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
692:                                    CALC_GRAD(MODYIJ) 
693:  
694:                             ! Atom position derivatives 
695:                             DO I = 1, 3 
696:                                 TEMP_GRAD(I) = DVDY / MODYIJ * & 
697:                                              DOT_PRODUCT(YIJ(:), BD%ORTHOG(:,I)) 
698:                             END DO 
699:                             GRAD_I(:) = GRAD_I(:) + TEMP_GRAD(1:3) 
700:                             GRAD_J(:) = GRAD_J(:) - TEMP_GRAD(1:3) 
701:  
702:                             ! Lattice parameter derivatives 
703:                             TEMP_GRAD(:) = 0.D0 
704:                             DO I = 1, 6 ! Loop over the different box parameters 
705:                                 DO J = 1, 3 ! Loop over x, y, z contributions 
706:                                     TEMP_GRAD(I) = TEMP_GRAD(I) + & 
707:                                     YIJ(J)*DOT_PRODUCT(BD%DERIV(J,:,I), RIJ(:)) 
708:                                 END DO 
709:                             END DO 
710:                             TEMP_GRAD(:) = TEMP_GRAD(:) * DVDY / MODYIJ 
711:                             GRAD_BOX(:) = GRAD_BOX(:) + TEMP_GRAD(:) 
712:                         END IF ! GTEST 
713:                     END IF ! End if less than cutoff 
714:                 END DO ! z loop 
715:             END DO ! y loop 
716:         END DO ! x loop 
717:  
718:     END SUBROUTINE LJ_GAUSS_PER_TRI 
719:  
720: !------------------------------------------------------------------------------- 
721: ! 
722: ! Outputs the coordinates to the file ljgauss.xyz430: ! Outputs the coordinates to the file ljgauss.xyz
723: !431: !
724: !-------------------------------------------------------------------------------432: !-------------------------------------------------------------------------------
725:     SUBROUTINE VIEW_LJ_GAUSS()433:     SUBROUTINE VIEW_LJ_GAUSS()
726: 434: 
727:         ! NATOMS: number of coordinates435:         ! NATOMS: number of coordinates
728:         ! NSAVE: number of minima to write436:         ! NSAVE: number of minima to write
729:         ! BOXx: Box dimensions437:         ! BOXx: Box dimensions
730:         ! PERIODIC: whether periodic boundary conditions are in use438:         ! PERIODIC: whether periodic boundary conditions are in use
731:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ439:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ
737: 445: 
738:         IMPLICIT NONE446:         IMPLICIT NONE
739: 447: 
740:         ! GETUNIT: function for fiding a free file unit448:         ! GETUNIT: function for fiding a free file unit
741:         ! FUNIT: output file unit449:         ! FUNIT: output file unit
742:         ! NMOL: actual number of particles450:         ! NMOL: actual number of particles
743:         ! I, J: loop iterators451:         ! I, J: loop iterators
744:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J452:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J
745: 453: 
746:         ! BOX_SIZE: convenience array for the box sizes454:         ! BOX_SIZE: convenience array for the box sizes
747:         ! BOX_ANGLES: convenience array for the box angles455:         DOUBLE PRECISION :: BOX_SIZE(3)
748:         DOUBLE PRECISION :: BOX_SIZE(3), BOX_ANGLES(3) 
749: 456: 
750:         ! Set NMOL and box parameters457:         ! Set NMOL and box parameters
751:         SELECT CASE (LJ_GAUSS_MODE)458:         IF (LJ_GAUSS_MODE .EQ. 1) THEN
752:         CASE (0) 
753:             NMOL = NATOMS 
754:             IF (PERIODIC) THEN 
755:                 BOX_SIZE(1) = BOXLX 
756:                 BOX_SIZE(2) = BOXLY 
757:                 BOX_SIZE(3) = BOXLZ 
758:             ELSE 
759:                 BOX_SIZE(:) = 1.D0 
760:             END IF 
761:             BOX_ANGLES(:) = PI/2.D0 
762:         CASE (1) 
763:             NMOL = NATOMS - 1459:             NMOL = NATOMS - 1
764:             BOX_ANGLES(:) = PI/2.D0460:         ELSE IF (LJ_GAUSS_MODE .EQ. 0 .AND. PERIODIC) THEN
765:         CASE (2)461:             NMOL = NATOMS
766:             NMOL = NATOMS - 2462:             BOX_SIZE(1) = BOXLX
767:         END SELECT463:             BOX_SIZE(2) = BOXLY
 464:             BOX_SIZE(3) = BOXLZ
 465:         ELSE
 466:             NMOL = NATOMS
 467:         END IF
768: 468: 
769:         ! Open the file469:         ! Open the file
770:         FUNIT = GETUNIT()470:         FUNIT = GETUNIT()
771:         OPEN(UNIT=FUNIT, FILE='ljgauss.xyz', STATUS='UNKNOWN')471:         OPEN(UNIT=FUNIT, FILE='ljgauss.xyz', STATUS='UNKNOWN')
772: 472: 
773:         ! Loop over the configurations473:         ! Loop over the configurations
774:         DO I = 1, NSAVE474:         DO I = 1, NSAVE
775: 475: 
776:             WRITE(FUNIT,'(I6)') NMOL476:             WRITE(FUNIT,'(I6)') NMOL
777: 477: 
778:             WRITE(FUNIT, '(A,I6,A,F20.10,A,I8,A,F20.10)') 'Energy of minimum ',&478:             WRITE(FUNIT, '(A,I6,A,F20.10,A,I8,A,F20.10)') 'Energy of minimum ',&
779:                 I, '=', QMIN(I), ' first found at step ', FF(I), &479:                 I, '=', QMIN(I), ' first found at step ', FF(I), &
780:                 ' Energy per particle = ', QMIN(I)/NMOL480:                 ' Energy per particle = ', QMIN(I)/NMOL
781: 481: 
782:             ! Set box parameters for varying box482:             ! Set box parameters for varying box
783:             IF (LJ_GAUSS_MODE .EQ. 1 .OR. LJ_GAUSS_MODE .EQ. 2) THEN483:             IF (LJ_GAUSS_MODE .EQ. 1) THEN
784:                 BOX_SIZE(:) = QMINP(I, 3*NATOMS-2:3*NATOMS)484:                 BOX_SIZE(:) = QMINP(I, 3*NATOMS-2:3*NATOMS)
785:             END IF485:             END IF
786:             IF (LJ_GAUSS_MODE .EQ. 2) THEN 
787:                 BOX_ANGLES(:) = QMINP(I, 3*NATOMS-5:3*NATOMS-3) 
788:             END IF 
789: 486: 
790:             DO J = 1, NMOL487:             DO J = 1, NMOL
791: 488: 
792:                 IF (PERIODIC .AND. J .EQ. 1) THEN489:                 IF (PERIODIC) THEN
793:                     WRITE(FUNIT, '(A,3F20.10,A,3(A,F20.10))') 'O ', &490:                     IF (J .EQ. 1) THEN
794:                         QMINP(I, 3*J-2:3*J)*BOX_SIZE(:), ' bbox_xyz', &491:                         WRITE(FUNIT, '(A,3F20.10,A,3(A,F20.10))') 'O ', &
795:                         ' 0.0', BOX_SIZE(1), ' 0.0', BOX_SIZE(2), &492:                               QMINP(I, 3*J-2:3*J)*BOX_SIZE(:), ' bbox_xyz', &
796:                         ' 0.0', BOX_SIZE(3)493:                               ' 0.0', BOX_SIZE(1), ' 0.0', BOX_SIZE(2), &
797:                 ELSE IF (PERIODIC .AND. J .EQ. 2) THEN494:                               ' 0.0', BOX_SIZE(3)
798:                     WRITE(FUNIT, '(A,3F20.10,3(A,F20.10))') 'O ', &495:                     ELSE
799:                         QMINP(I, 3*J-2:3*J)*BOX_SIZE(:), ' alpha = ', &496:                         WRITE(FUNIT, '(A,3F20.10)') 'O ', &
800:                         BOX_ANGLES(1), ' beta = ', BOX_ANGLES(2), &497:                             QMINP(I, 3*J-2:3*J)*BOX_SIZE(:)
801:                         ' gamma = ', BOX_ANGLES(3)498:                     END IF
802:                 ELSE499:                 ELSE
803:                     WRITE(FUNIT, '(A,3F20.10)') 'O ', &500:                     WRITE(FUNIT, '(A,3F20.10)') 'O ', QMINP(I, 3*J-2:3*J)
804:                           QMINP(I, 3*J-2:3*J)*BOX_SIZE(:)*SIN(BOX_ANGLES(:)) 
805:                 END IF501:                 END IF
806: 502: 
807:             END DO ! Loop over particles503:             END DO ! Loop over particles
808: 504: 
809:         END DO ! Loop over the configurations505:         END DO ! Loop over the configurations
810: 506: 
811:         CLOSE(FUNIT)507:         CLOSE(FUNIT)
812: 508: 
813:     END SUBROUTINE VIEW_LJ_GAUSS509:     END SUBROUTINE VIEW_LJ_GAUSS
814: 510: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0