hdiff output

r32229/lj_gauss.f90 2017-03-29 19:30:08.842165184 +0100 r32228/lj_gauss.f90 2017-03-29 19:30:09.286171066 +0100
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 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 parameters 
 27: 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
 28: ! Public functions 
 29: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS 27: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS
 30: ! Private parameters 
 31: PRIVATE :: LJ_GAUSS_SIGMASQ, RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT, PI 28: PRIVATE :: LJ_GAUSS_SIGMASQ, RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT, PI
 32: PRIVATE :: SCL_FCT 29: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_GRAD
 33: ! Private functions 30: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES
 34: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_ENERGY_LJ 
 35: PRIVATE :: CALC_ENERGY_GAUSS, CALC_GRAD, CALC_FCTS, CONSTRAIN_PARAMETERS 
 36: PRIVATE :: BOX_DERIV, CALC_BOX_DERIV, CHECK_ANGLES, REJECT 
 37:  31: 
 38: ! EPS: strength of the Gaussian potential. The LJ strength is 1. 32: ! EPS: strength of the Gaussian potential. The LJ strength is 1.
 39: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ 33: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ
 40: !     rmin is 1. 34: !     rmin is 1.
 41: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of 35: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of
 42: !         both parts of the potential go smoothly to zero. 36: !         both parts of the potential go smoothly to zero.
 43: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT 37: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT
 44:  38: 
 45: ! LJ_GAUSS_MODE: 0 standard minimisation 39: ! LJ_GAUSS_MODE: 0 standard minimisation
 46: !                1 the final triplet of coordinates will be interpreted as the 40: !                1 the final triplet of coordinates will be interpreted as the
 47: !                  unit cell parameters of an orthogonal unit cell, and 41: !                  unit cell parameters of an orthogonal unit cell, and
 48: !                  optimised 42: !                  minimised
 49: !                2 the final triplet of coordinates are the unit cell length 
 50: !                  parameters and the second last triplet are the unit cell 
 51: !                  angles, giving a triclinic unit cell, which will be optimised 
 52: !                3 as for 2, but the third last triplet are the potential 
 53: !                  parameters LJ_GAUSS_EPS, LJ_GAUSS_R0, 0.D0 and the potential 
 54: !                  parameters will be optimised 
 55: INTEGER :: LJ_GAUSS_MODE 43: INTEGER :: LJ_GAUSS_MODE
 56:  44: 
 57: ! LJ_GAUSS_SIGMA: variation (width) of the Gaussian potential 45: ! LJ_GAUSS_SIGMA: variation (width) of the Gaussian potential
 58: DOUBLE PRECISION, PARAMETER :: LJ_GAUSS_SIGMASQ = 0.02D0 46: DOUBLE PRECISION, PARAMETER :: LJ_GAUSS_SIGMASQ = 0.02D0
 59:  47: 
 60: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance 48: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance
 61: DOUBLE PRECISION :: RCUTINV6, RCUTINV12 49: DOUBLE PRECISION :: RCUTINV6, RCUTINV12
 62:  50: 
 63: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle 51: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle
 64: !                      distance 52: !                      distance
 65: ! SCL_FCT: overall scale factor for the potential, to allow comparison between 53: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT
 66: !          different parameters 
 67: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT, SCL_FCT 
 68:  54: 
 69: ! PI: the usual meaning 55: ! PI: the usual meaning
 70: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0) 56: DOUBLE PRECISION, PARAMETER :: PI = 4*ATAN(1.D0)
 71:  57: 
 72: ! BOX_DERIV: stores information for calculating lattice derivatives 58: ! BOX_DERIV: stores information for calculating lattice derivatives
 73: TYPE :: BOX_DERIV 59: TYPE :: BOX_DERIV
 74:     ! ORTHOG: the orthogonal tansformation matrix, for changing from 60:     ! ORTHOG: the orthogonal tansformation matrix, for changing from
 75:     !         crystallographic space to orthogonal space first index is row, 61:     !         crystallographic space to orthogonal space first index is row,
 76:     !         second is column 62:     !         second is column
 77:     DOUBLE PRECISION :: ORTHOG(3,3) 63:     DOUBLE PRECISION :: ORTHOG(3,3)
 78:  64: 
 79:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters, 65:     ! DERIV: derivatives of the orthogonal matrix wrt lattice parameters,
 80:     !        first index indicates the row, second index the column and 66:     !        first index indicates the row, second index the column and
 81:     !        third index the parameter, with 1-3 being the angles and 4-6 being 67:     !        third index the parameter, with 1-3 being the angles and 4-6 being
 82:     !        the lengths 68:     !        the lengths
 83:     DOUBLE PRECISION :: DERIV(3, 3, 6) 69:     DOUBLE PRECISION :: DERIV(3, 3, 6)
 84:  70: 
 85:     ! RECIP: lengths of the reciprocal lattice vectors 71:     ! RECIP: lengths of the reciprocal lattice vectors
 86:     DOUBLE PRECISION :: RECIP(3) 72:     DOUBLE PRECISION :: RECIP(3)
 87:  
 88:     ! CELL_RANGE: number of cells in each direction that must be evaluated 
 89:     INTEGER :: CELL_RANGE(3) 
 90: END TYPE BOX_DERIV 73: END TYPE BOX_DERIV
 91:  74: 
 92: CONTAINS 75: CONTAINS
 93:  76: 
 94: !------------------------------------------------------------------------------- 77: !-------------------------------------------------------------------------------
 95: ! 78: !
 96: ! Main routine for calculating the energy and gradients. 79: ! Main routine for calculating the energy and gradients.
 97: ! X: coordinates array 80: ! X: coordinates array
 98: ! GRAD: gradients array 81: ! GRAD: gradients array
 99: ! ENERGY: calculated energy 82: ! ENERGY: calculated energy
118:         INTEGER :: I, J, NMOL101:         INTEGER :: I, J, NMOL
119: 102: 
120:         ! RIJ: Interparticle vector103:         ! RIJ: Interparticle vector
121:         ! MODRIJ: distance between a pair of particles104:         ! MODRIJ: distance between a pair of particles
122:         ! DVDR: derivative of pair potential wrt MODRIJ105:         ! DVDR: derivative of pair potential wrt MODRIJ
123:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR106:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR
124: 107: 
125:         ! Factors for triclinic lattice108:         ! Factors for triclinic lattice
126:         TYPE(BOX_DERIV) :: BD109:         TYPE(BOX_DERIV) :: BD
127: 110: 
128:  !       DO I = 1, NATOMS 
129:  !           WRITE (MYUNIT, *) X(3*I-2:3*I) 
130:  !       END DO 
131:  
132:         ! Initialise output variables111:         ! Initialise output variables
133:         ENERGY = 0.D0112:         ENERGY = 0.D0
134:         GRAD(:) = 0.D0113:         GRAD(:) = 0.D0
135: 114: 
136:         ! Calculate factors that do depend on parameters, but not on RIJ 
137:         CALL CALC_FCTS(X(NATOMS*3-8:NATOMS*3-7)) 
138:  
139:         IF (.NOT. PERIODIC) THEN115:         IF (.NOT. PERIODIC) THEN
140:             IF (LJ_GAUSS_MODE .EQ. 0) THEN116:             IF (LJ_GAUSS_MODE .EQ. 0) THEN
141:                 ! Normal cluster117:                 ! Normal cluster
142:                 NMOL = NATOMS118:                 NMOL = NATOMS
143:                 DO I = 1, NMOL-1 ! Outer loop over particles119:                 DO I = 1, NMOL-1 ! Outer loop over particles
144:                     DO J = I+1, NMOL ! Inner loop over particles120:                     DO J = I+1, NMOL ! Inner loop over particles
145:                         ! Get the particle-particle distance121:                         ! Get the particle-particle distance
146:                         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)
147:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))123:                         MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
148: 124: 
201:                                           GRAD(3*NATOMS-2:3*NATOMS))177:                                           GRAD(3*NATOMS-2:3*NATOMS))
202:                     END DO ! Inner loop over particles178:                     END DO ! Inner loop over particles
203:                 END DO ! Outer loop over particles179:                 END DO ! Outer loop over particles
204: 180: 
205:             CASE(2)181:             CASE(2)
206:                 ! Triclinic varying lattice182:                 ! Triclinic varying lattice
207:                 ! Reset all particles to within the box183:                 ! Reset all particles to within the box
208:                 CALL PERIODIC_RESET(X(:))184:                 CALL PERIODIC_RESET(X(:))
209: 185: 
210:                 ! Check whether the unit cell angles are physically possible.186:                 ! Check whether the unit cell angles are physically possible.
211:                 ! Reject if not. 
212:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN 
213:                     CALL REJECT(ENERGY, GRAD) 
214:                     RETURN 
215:                 END IF 
216:  
217:                 ! Calculate box derivative factors 
218:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), & 
219:                                     X(3*NATOMS-2:3*NATOMS  )) 
220:  
221:                 ! Reject steps where the cell range has got bigger than 20. Such 
222:                 ! unit cells are highly distorted and there are probably 
223:                 ! equivalent ones of a more usual shape. Allowing them leads to 
224:                 ! very slow run times. 
225:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. 20)) THEN 
226:                     CALL REJECT(ENERGY, GRAD) 
227:                     RETURN 
228:                 END IF 
229:  
230: !                WRITE(MYUNIT, *) 'Box lengths: ', X(3*NATOMS-2:3*NATOMS) 
231: !                WRITE(MYUNIT, *) 'Box angles: ', X(3*NATOMS-5:3*NATOMS-3) 
232: !                WRITE(MYUNIT, *) 'CELL_RANGE(:) ', BD%CELL_RANGE(:) 
233:  
234:                 NMOL = NATOMS-2 
235:                 ! We need to include self-interactions of particles in different 
236:                 ! unit cells 
237:                 DO I = 1, NMOL! Outer loop over particles 
238:                     DO J = I, NMOL ! Inner loop over particles 
239:                         ! Add energy and gradients 
240:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), & 
241:                             GRAD(3*J-2:3*J), ENERGY, GTEST, & 
242:                             GRAD(3*NATOMS-5:3*NATOMS), BD) 
243:                     END DO ! Inner loop over particles 
244:                 END DO ! Outer loop over particles 
245:             CASE(3) 
246:                 ! Triclinic varying lattice, with varying parameters 
247:                 ! Reset all particles to within the box 
248:                 CALL PERIODIC_RESET(X(:)) 
249:  
250:                 ! Check whether the unit cell angles are physically possible. 
251:                 ! If we have a disallowed unit cell, we set the energy to very187:                 ! If we have a disallowed unit cell, we set the energy to very
252:                 ! high and the gradients to very small, so the step is188:                 ! high and the gradients to very small, so the step is
253:                 ! 'converged' but gets rejected.189:                 ! 'converged' but gets rejected.
254:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN190:                 IF (.NOT. CHECK_ANGLES(X(3*NATOMS-5:3*NATOMS-3))) THEN
255:                     ENERGY = 1.D20191:                     ENERGY = 1.D20
256:                     GRAD(:) = 1.D-20192:                     GRAD(:) = 1.D-20
257:                     RETURN193:                     RETURN
258:                 END IF194:                 END IF
259: 195: 
260:                 ! Calculate box derivative factors196:                 ! Calculate box derivative factors
261:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &197:                 BD = CALC_BOX_DERIV(X(3*NATOMS-5:3*NATOMS-3), &
262:                                     X(3*NATOMS-2:3*NATOMS  ))198:                                     X(3*NATOMS-2:3*NATOMS  ))
263: 199: 
264:                 ! Reject steps where the cell range has got bigger than 20. Such200:                 NMOL = NATOMS-2
265:                 ! unit cells are highly distorted and there are probably 
266:                 ! equivalent ones of a more usual shape. Allowing them leads to 
267:                 ! very slow run times. 
268:                 IF (.NOT. ALL(BD%CELL_RANGE .LE. 20)) THEN 
269:                     CALL REJECT(ENERGY, GRAD) 
270:                     RETURN 
271:                 END IF 
272:  
273:                 NMOL = NATOMS-3 
274:                 ! We need to include self-interactions of particles in different201:                 ! We need to include self-interactions of particles in different
275:                 ! unit cells202:                 ! unit cells
276:                 DO I = 1, NMOL! Outer loop over particles203:                 DO I = 1, NMOL! Outer loop over particles
277:                     DO J = I, NMOL ! Inner loop over particles204:                     DO J = I, NMOL ! Inner loop over particles
278:                         ! Add energy and gradients205:                         ! Add energy and gradients
279:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &206:                         CALL LJ_GAUSS_PER_TRI(I, J, X(:), GRAD(3*I-2:3*I), &
280:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &207:                             GRAD(3*J-2:3*J), ENERGY, GTEST, &
281:                             GRAD(3*NATOMS-5:3*NATOMS), BD, &208:                             GRAD(3*NATOMS-5:3*NATOMS), BD)
282:                             GRAD(3*NATOMS-8:3*NATOMS-7)) 
283:                     END DO ! Inner loop over particles209:                     END DO ! Inner loop over particles
284:                 END DO ! Outer loop over particles210:                 END DO ! Outer loop over particles
285:  
286:                 ! Constrain the potential parameters with a harmonic repulsion 
287:                 CALL CONSTRAIN_PARAMETERS(ENERGY, GRAD(3*NATOMS-8:3*NATOMS-7), & 
288:                                           GTEST) 
289:             END SELECT ! Mode selections211:             END SELECT ! Mode selections
290:         END IF ! Periodic212:         END IF ! Periodic
291: 213: 
292:         ! Apply the overall scaling 
293:         ENERGY = ENERGY * SCL_FCT 
294:         GRAD(:) = GRAD(:) * SCL_FCT 
295:  
296:         ! The gradients wrt the potential parameters already include SCL_FCT 
297:         IF (LJ_GAUSS_MODE .EQ. 3) THEN 
298:             GRAD(3*NATOMS-8:3*NATOMS-7) = GRAD(3*NATOMS-8:3*NATOMS-7) / SCL_FCT 
299:         END IF 
300:  
301:     END SUBROUTINE LJ_GAUSS214:     END SUBROUTINE LJ_GAUSS
302: 215: 
303: !-------------------------------------------------------------------------------216: !-------------------------------------------------------------------------------
304: !217: !
305: ! Rejects a step by setting the energy very high and the gradients very low. 
306: ! This tricks L-BFGS into thinking it's converged, but the high energy of the 
307: ! quench will ensure it is rejected. 
308: ! ENERGY: system energy 
309: ! GRAD: system gradients 
310: ! 
311: !------------------------------------------------------------------------------- 
312:  
313:     SUBROUTINE REJECT(ENERGY, GRAD) 
314:  
315:         ! NATOMS: number of atoms 
316:         USE COMMONS, ONLY: NATOMS 
317:  
318:         IMPLICIT NONE 
319:  
320:         DOUBLE PRECISION, INTENT(OUT) :: ENERGY, GRAD(3*NATOMS) 
321:  
322:         ENERGY = 1.D20 
323:         GRAD(:) = 1.D-20 
324:  
325:     END SUBROUTINE REJECT 
326:  
327: !------------------------------------------------------------------------------- 
328: ! 
329: ! Calculates the factors that depend on the potential parameters but not on the 
330: ! distance between particles. 
331: ! PARAMS: the potential parameters when those are to be adjusted, otherwise 
332: !         ignored 
333: ! 
334: !------------------------------------------------------------------------------- 
335:  
336:     SUBROUTINE CALC_FCTS(PARAMS) 
337:  
338:         IMPLICIT NONE 
339:  
340:         DOUBLE PRECISION :: PARAMS(2) 
341:  
342:         ! Adjust the potential parameters, if the mode requires it 
343:         IF (LJ_GAUSS_MODE .EQ. 3) THEN 
344:             LJ_GAUSS_EPS = PARAMS(1) 
345:             LJ_GAUSS_R0 = PARAMS(2) 
346:         END IF 
347:  
348:         ! Overall scale factor 
349:         SCL_FCT = 1.D0 / (SQRT(2.D0 * PI * LJ_GAUSS_SIGMASQ) * LJ_GAUSS_EPS + & 
350:                           2.D0**(5.D0/6.D0) * 12.D0/55.D0) 
351:  
352:         ! Factors for the Gaussian 
353:         EXP_RCUT = LJ_GAUSS_EPS * EXP(- (LJ_GAUSS_RCUT - LJ_GAUSS_R0)**2 / & 
354:                                         (2.D0*LJ_GAUSS_SIGMASQ)) 
355:         PREF_RCUT = (LJ_GAUSS_RCUT - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ 
356:  
357:     END SUBROUTINE CALC_FCTS 
358:  
359: !------------------------------------------------------------------------------- 
360: ! 
361: ! Calculates the energy contribution for a pair of particles.218: ! Calculates the energy contribution for a pair of particles.
362: ! MODRIJ: distance beteen the particle pair219: ! MODRIJ: distance beteen the particle pair
363: !220: !
364: !-------------------------------------------------------------------------------221: !-------------------------------------------------------------------------------
365:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)222:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)
366: 223: 
367:         IMPLICIT NONE224:         IMPLICIT NONE
368: 225: 
369:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ226:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ
370: 227: 
371:         CALC_ENERGY = CALC_ENERGY_LJ(MODRIJ) + CALC_ENERGY_GAUSS(MODRIJ) 
372:  
373:     END FUNCTION CALC_ENERGY 
374:  
375: !------------------------------------------------------------------------------- 
376: ! 
377: ! Calculates the LJ contribution to the energy 
378: ! MODRIJ: distance beteen the particle pair 
379: ! 
380: !------------------------------------------------------------------------------- 
381:  
382:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY_LJ (MODRIJ) 
383:  
384:         IMPLICIT NONE 
385:  
386:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
387:  
388:         ! LJ6, LJ12: factors for the LJ potential228:         ! LJ6, LJ12: factors for the LJ potential
 229:         DOUBLE PRECISION :: LJ6, LJ12
 230:         ! EXPRIJ: factor for the Gauss potential
 231:         DOUBLE PRECISION :: EXP_RIJ
389:         ! ENERGY: temporary energy store232:         ! ENERGY: temporary energy store
390:         DOUBLE PRECISION :: LJ6, LJ12, ENERGY233:         DOUBLE PRECISION :: ENERGY
 234: 
 235:         ! Initialise output variables
 236:         ENERGY = 0.D0
391: 237: 
392:         ! Evaluate LJ energy238:         ! Evaluate LJ energy
393:         LJ6 = (1.D0/MODRIJ)**6239:         LJ6 = (1.D0/MODRIJ)**6
394:         LJ12 = LJ6**2240:         LJ12 = LJ6**2
395:         ENERGY = LJ12 - 2.D0*LJ6241:         ENERGY = ENERGY + LJ12 - 2.D0*LJ6
396:         ! Correction to LJ to make potential go to zero at RCUT242:         ! Correction to LJ to make potential go to zero at RCUT
397:         ENERGY = ENERGY - (RCUTINV12 - 2*RCUTINV6)243:         ENERGY = ENERGY - (RCUTINV12 - 2*RCUTINV6)
398:         ! Correction to LJ to make derivative go to zero at RCUT244:         ! Correction to LJ to make derivative go to zero at RCUT
399:         ENERGY = ENERGY + (12.D0/LJ_GAUSS_RCUT) * &245:         ENERGY = ENERGY + (12.D0/LJ_GAUSS_RCUT) * &
400:                           (RCUTINV12 - RCUTINV6)* &246:                           (RCUTINV12 - RCUTINV6)* &
401:                           (MODRIJ - LJ_GAUSS_RCUT)247:                           (MODRIJ - LJ_GAUSS_RCUT)
402:         ! Correction to LJ to make second derivatives go to zero at248:         ! Correction to LJ to make second derivatives go to zero at
403:         ! RCUT249:         ! RCUT
404:         ENERGY = ENERGY + (1.D0/LJ_GAUSS_RCUT**2) * &250:         ENERGY = ENERGY + (1.D0/LJ_GAUSS_RCUT**2) * &
405:                           (42.D0*RCUTINV6 - 78.D0*RCUTINV12) * &251:                           (42.D0*RCUTINV6 - 78.D0*RCUTINV12) * &
406:                           (MODRIJ - LJ_GAUSS_RCUT)**2252:                           (MODRIJ - LJ_GAUSS_RCUT)**2
407:          
408:         ! Set output 
409:         CALC_ENERGY_LJ = ENERGY 
410:  
411:     END FUNCTION CALC_ENERGY_LJ 
412:  
413: !------------------------------------------------------------------------------- 
414: ! 
415: ! Calculates the Gaussian contribution to the energy 
416: ! MODRIJ: distance beteen the particle pair 
417: ! 
418: !------------------------------------------------------------------------------- 
419:  
420:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY_GAUSS (MODRIJ) 
421:  
422:         IMPLICIT NONE 
423:  
424:         DOUBLE PRECISION, INTENT(IN) :: MODRIJ 
425:  
426:         ! EXP_RIJ: factor for the Gauss potential 
427:         ! ENERGY: temporary energy store 
428:         DOUBLE PRECISION :: EXP_RIJ, ENERGY 
429:          
430:         ! Set output 
431:         CALC_ENERGY_GAUSS = ENERGY 
432:  
433:         ! Evaluate Gaussian energy253:         ! Evaluate Gaussian energy
434:         EXP_RIJ = LJ_GAUSS_EPS * &254:         EXP_RIJ = LJ_GAUSS_EPS * &
435:                   EXP(- (MODRIJ - LJ_GAUSS_R0)**2 / &255:                   EXP(- (MODRIJ - LJ_GAUSS_R0)**2 / &
436:                         (2.D0*LJ_GAUSS_SIGMASQ))256:                         (2.D0*LJ_GAUSS_SIGMASQ))
437:         ENERGY = - EXP_RIJ257:         ENERGY = ENERGY - EXP_RIJ
438:         ! Correction to Gauss to make potential go to zero at RCUT258:         ! Correction to Gauss to make potential go to zero at RCUT
439:         ENERGY = ENERGY + EXP_RCUT259:         ENERGY = ENERGY + EXP_RCUT
440:         ! Correction to Gauss to make derivatives go to zero at RCUT260:         ! Correction to Gauss to make derivatives go to zero at RCUT
441:         ENERGY = ENERGY - PREF_RCUT*EXP_RCUT*(MODRIJ-LJ_GAUSS_RCUT)261:         ENERGY = ENERGY - PREF_RCUT*EXP_RCUT*(MODRIJ-LJ_GAUSS_RCUT)
442:         ! Correction to Gauss to make second derivatives go to zero at262:         ! Correct to Gauss to make second derivatives go to zero at
443:         ! RCUT263:         ! RCUT
444:         ENERGY = ENERGY + EXP_RCUT * (PREF_RCUT**2 - 1.D0/LJ_GAUSS_SIGMASQ) * &264:         ENERGY = ENERGY - (EXP_RCUT / LJ_GAUSS_SIGMASQ - &
 265:                            (PREF_RCUT**2) * EXP_RCUT) * &
445:                           0.5D0*((MODRIJ - LJ_GAUSS_RCUT)**2)266:                           0.5D0*((MODRIJ - LJ_GAUSS_RCUT)**2)
446: 267: 
447:         ! Set output268:         CALC_ENERGY = ENERGY
448:         CALC_ENERGY_GAUSS = ENERGY 
449: 269: 
450:     END FUNCTION CALC_ENERGY_GAUSS270:     END FUNCTION CALC_ENERGY
451: 271: 
452: !-------------------------------------------------------------------------------272: !-------------------------------------------------------------------------------
453: !273: !
454: ! Calculates the pair potential gradient for a pair of particles.274: ! Calculates the pair potential gradient for a pair of particles.
455: ! MODRIJ: distance beteen the particle pair275: ! MODRIJ: distance beteen the particle pair
456: !276: !
457: !-------------------------------------------------------------------------------277: !-------------------------------------------------------------------------------
458:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ)278:     PURE DOUBLE PRECISION FUNCTION CALC_GRAD (MODRIJ)
459: 279: 
460:         IMPLICIT NONE280:         IMPLICIT NONE
513: !-------------------------------------------------------------------------------333: !-------------------------------------------------------------------------------
514:     SUBROUTINE LJ_GAUSS_INITIALISE ()334:     SUBROUTINE LJ_GAUSS_INITIALISE ()
515: 335: 
516:         ! MYUNIT: file unit for main GMIN output336:         ! MYUNIT: file unit for main GMIN output
517:         USE COMMONS, ONLY: MYUNIT337:         USE COMMONS, ONLY: MYUNIT
518: 338: 
519:         IMPLICIT NONE339:         IMPLICIT NONE
520: 340: 
521:         ! Check we have sensible value for the mode341:         ! Check we have sensible value for the mode
522:         IF (LJ_GAUSS_MODE .NE. 1 .AND. LJ_GAUSS_MODE .NE. 0 .AND. &342:         IF (LJ_GAUSS_MODE .NE. 1 .AND. LJ_GAUSS_MODE .NE. 0 .AND. &
523:             LJ_GAUSS_MODE .NE. 2 .AND. LJ_GAUSS_MODE .NE. 3) THEN343:             LJ_GAUSS_MODE .NE. 2) THEN
524:             WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: mode must be 0, 1, 2 or 3'344:             WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: mode must be 0, 1 or 2'
525:             STOP345:             STOP
526:         END IF346:         END IF
527: 347: 
528:         ! Calculate some factors for LJ that don't depend on RIJ or the348:         ! Calculate some factors for LJ that don't depend on RIJ
529:         ! potential parameters 
530:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6349:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6
531:         RCUTINV12 = RCUTINV6**2350:         RCUTINV12 = RCUTINV6**2
532: 351: 
 352:         ! 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 / &
 354:                                         (2.D0*LJ_GAUSS_SIGMASQ))
 355:         PREF_RCUT = (LJ_GAUSS_RCUT - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ
 356: 
533:     END SUBROUTINE LJ_GAUSS_INITIALISE357:     END SUBROUTINE LJ_GAUSS_INITIALISE
534: 358: 
535: !-------------------------------------------------------------------------------359: !-------------------------------------------------------------------------------
536: !360: !
537: ! Resets all coordinates to within the unit cell361: ! Resets all coordinates to within the unit cell
538: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled362: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled
539: ! by the box lengths.363: ! by the box lengths.
540: ! COORDS: coordinates of the particles364: ! COORDS: coordinates of the particles
541: !365: !
542: !-------------------------------------------------------------------------------366: !-------------------------------------------------------------------------------
551:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)375:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)
552: 376: 
553:         ! Sort out the unit cell sizes, if they are varying377:         ! Sort out the unit cell sizes, if they are varying
554:         IF (LJ_GAUSS_MODE .EQ. 1 .OR. LJ_GAUSS_MODE .EQ. 2) THEN378:         IF (LJ_GAUSS_MODE .EQ. 1 .OR. LJ_GAUSS_MODE .EQ. 2) THEN
555:             BOXLX = COORDS(3*NATOMS - 2)379:             BOXLX = COORDS(3*NATOMS - 2)
556:             BOXLY = COORDS(3*NATOMS - 1)380:             BOXLY = COORDS(3*NATOMS - 1)
557:             BOXLZ = COORDS(3*NATOMS)381:             BOXLZ = COORDS(3*NATOMS)
558:         END IF382:         END IF
559: 383: 
560:         ! Reset coordinates to within the box384:         ! Reset coordinates to within the box
561:         SELECT CASE (LJ_GAUSS_MODE)385:         IF (LJ_GAUSS_MODE .EQ. 0) THEN
562:         CASE(0) 
563:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:))386:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:))
564:         CASE(1)387:         ELSE IF (LJ_GAUSS_MODE .EQ. 1) THEN
565:             COORDS(1:3*NATOMS-3) = COORDS(1:3*NATOMS-3) - &388:             COORDS(1:3*NATOMS-3) = COORDS(1:3*NATOMS-3) - &
566:                                    FLOOR(COORDS(1:3*NATOMS-3))389:                                    FLOOR(COORDS(1:3*NATOMS-3))
567:         CASE(2)390:         ELSE IF (LJ_GAUSS_MODE .EQ. 2) THEN
568:             COORDS(1:3*NATOMS-6) = COORDS(1:3*NATOMS-6) - &391:             COORDS(1:3*NATOMS-6) = COORDS(1:3*NATOMS-6) - &
569:                                    FLOOR(COORDS(1:3*NATOMS-6))392:                                    FLOOR(COORDS(1:3*NATOMS-6))
570:             ! Make the lattice angles modulo 2*pi393:             ! Make the lattice angles modulo 2*pi
571:             ! They will be checked later394:             ! They will be checked later
572:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &395:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - &
573:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))396:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI))
574:         CASE(3)397:         END IF
575:            COORDS(1:3*NATOMS-9) = COORDS(1:3*NATOMS-9) - & 
576:                                    FLOOR(COORDS(1:3*NATOMS-9)) 
577:             ! Make the lattice angles modulo 2*pi 
578:             ! They will be checked later 
579:             COORDS(3*NATOMS-5:3*NATOMS-3) = COORDS(3*NATOMS-5:3*NATOMS-3) - & 
580:                 2.D0*PI*FLOOR(COORDS(3*NATOMS-5:3*NATOMS-3)/(2.D0*PI)) 
581:             ! Reset the unused parameter corodinate to zero 
582:             COORDS(3*NATOMS-6) = 0.D0 
583:         END SELECT 
584: 398: 
585:     END SUBROUTINE PERIODIC_RESET399:     END SUBROUTINE PERIODIC_RESET
586: 400: 
587: !-------------------------------------------------------------------------------401: !-------------------------------------------------------------------------------
588: !402: !
589: ! Finds the energy and gradients for the periodic system403: ! Finds the energy and gradients for the periodic system
590: ! IDI, IDJ: id numbers of the particles404: ! IDI, IDJ: id numbers of the particles
591: ! COORDS_I: coordinates of particle I405: ! COORDS_I: coordinates of particle I
592: ! COORDS_J: coordinates of particle J406: ! COORDS_J: coordinates of particle J
593: ! GRAD_I: gradients for particle I407: ! GRAD_I: gradients for particle I
743:         ! Derivatives of the bottom row wrt angles557:         ! Derivatives of the bottom row wrt angles
744:         DERIV(3, 3, 1) = SIZES(3) * S(1) * (C(1) - C(2) * C(3)) / (S(3) * V)558:         DERIV(3, 3, 1) = SIZES(3) * S(1) * (C(1) - C(2) * C(3)) / (S(3) * V)
745:         DERIV(3, 3, 2) = SIZES(3) * S(2) * (C(2) - C(3) * C(1)) / (S(3) * V)559:         DERIV(3, 3, 2) = SIZES(3) * S(2) * (C(2) - C(3) * C(1)) / (S(3) * V)
746:         DERIV(3, 3, 3) = SIZES(3) * ((C(3) - C(1) * C(2))/V - C(3) * V/S(3)**2)560:         DERIV(3, 3, 3) = SIZES(3) * ((C(3) - C(1) * C(2))/V - C(3) * V/S(3)**2)
747:         ! Derivatives of the bottom row wrt lengths561:         ! Derivatives of the bottom row wrt lengths
748:         DERIV(3, 3, 6) = V / S(3)562:         DERIV(3, 3, 6) = V / S(3)
749: 563: 
750:         ! Calculate the lengths of the reciprocal lattice vectors564:         ! Calculate the lengths of the reciprocal lattice vectors
751:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V)565:         CALC_BOX_DERIV%RECIP(:) = S(:) / (SIZES(:) * V)
752: 566: 
753:         ! Calculate the number of cells to check in each direction 
754:         CALC_BOX_DERIV%CELL_RANGE(:) = CEILING(LJ_GAUSS_RCUT * & 
755:             CALC_BOX_DERIV%RECIP(:)) 
756:  
757:         ! Set the output567:         ! Set the output
758:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:)568:         CALC_BOX_DERIV%ORTHOG(:,:) = ORTHOG(:,:)
759:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:)569:         CALC_BOX_DERIV%DERIV(:,:,:) = DERIV(:,:,:)
760: 570: 
761:     END FUNCTION CALC_BOX_DERIV571:     END FUNCTION CALC_BOX_DERIV
762: 572: 
763: !-------------------------------------------------------------------------------573: !-------------------------------------------------------------------------------
764: !574: !
765: ! Checks whether the unit cell angles are allowed. Itis possible to generate575: ! Checks whether the unit cell angles are allowed. Itis possible to generate
766: ! sets of angles which do not correspond to a physically possible unit cell.576: ! sets of angles which do not correspond to a physically possible unit cell.
789: 599: 
790:         ! Check all sums 0 < SUM < 2 pi and all angles 0 < ANGLE < pi600:         ! Check all sums 0 < SUM < 2 pi and all angles 0 < ANGLE < pi
791:         CHECK_ANGLES = ALL(SUMS .GT. 0.D0) .AND. ALL(SUMS .LT. 2*PI) .AND. &601:         CHECK_ANGLES = ALL(SUMS .GT. 0.D0) .AND. ALL(SUMS .LT. 2*PI) .AND. &
792:                        ALL(ANGLES .GT. 0.D0) .AND. ALL(ANGLES .LT. PI)602:                        ALL(ANGLES .GT. 0.D0) .AND. ALL(ANGLES .LT. PI)
793: 603: 
794:     END FUNCTION CHECK_ANGLES604:     END FUNCTION CHECK_ANGLES
795: 605: 
796: !-------------------------------------------------------------------------------606: !-------------------------------------------------------------------------------
797: !607: !
798: ! Finds the energy and gradients for the periodic system with a triclinic608: ! Finds the energy and gradients for the periodic system with a triclinic
799: ! varying unit cell. Optionally with varying potential parameters609: ! varying unit cell
800: ! IDI, IDJ: id numbers of the particles610: ! IDI, IDJ: id numbers of the particles
801: ! COORDS: coordinates of particles (in crystallographic space)611: ! COORDS: coordinates of particles (in crystallographic space)
802: !         and lattice parameters612: !         and lattice parameters
803: ! GRAD_I: gradients for particle I613: ! GRAD_I: gradients for particle I
804: ! GRAD_J: gradients for particle J614: ! GRAD_J: gradients for particle J
805: ! ENERGY: energy contribution for this pair615: ! ENERGY: energy contribution for this pair
806: ! RIJ: mininum length vector between the particles616: ! RIJ: mininum length vector between the particles
807: ! GTEST: whether to calculate gradients617: ! GTEST: whether to calculate gradients
808: ! GRAD_BOX: gradients for the box (angles, then lengths)618: ! GRAD_BOX: gradients for the box (angles, then lengths)
809: ! BD: factors for the lattice derivatives619: ! BD: factors for the lattice derivatives
810: ! GRAD_PARAMS: optional, derivatives of potential parameters 
811: !620: !
812: !-------------------------------------------------------------------------------621: !-------------------------------------------------------------------------------
813: 622: 
814:     SUBROUTINE LJ_GAUSS_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, &623:     SUBROUTINE LJ_GAUSS_PER_TRI(IDI, IDJ, COORDS, GRAD_I, GRAD_J, ENERGY, &
815:                                 GTEST, GRAD_BOX, BD, GRAD_PARAMS)624:                                 GTEST, GRAD_BOX, BD)
816: 625: 
817:         ! NATOMS: number of coordinates (including lattice parameters)626:         ! NATOMS: number of coordinates (including lattice parameters)
818:         USE COMMONS, ONLY: NATOMS627:         USE COMMONS, ONLY: NATOMS
819: 628: 
820:         IMPLICIT NONE629:         IMPLICIT NONE
821: 630: 
822:         INTEGER, INTENT(IN) :: IDI, IDJ631:         INTEGER, INTENT(IN) :: IDI, IDJ
823:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)632:         DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)
824:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY633:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY
825:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)634:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_BOX(6)
826:         LOGICAL, INTENT(IN) :: GTEST635:         LOGICAL, INTENT(IN) :: GTEST
827:         TYPE(BOX_DERIV), INTENT(IN) :: BD636:         TYPE(BOX_DERIV), INTENT(IN) :: BD
828:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_PARAMS(2) 
829: 637: 
830:         ! RJ: particle coordinates in crystallographic space638:         ! RJ: particle coordinates in crystallographic space
831:         ! RIJ: particle-particle vector in crystallographic space639:         ! RIJ: particle-particle vector in crystallographic space
832:         ! YIJ: particle-particle vector in orthogonal space640:         ! YIJ: particle-particle vector in orthogonal space
833:         ! MODYIJ: particle-particle distance in orthogonal space641:         ! MODYIJ: particle-particle distance in orthogonal space
834:         ! DVDY: derivative of pair potential wrt MODYIJ642:         ! DVDY: derivative of pair potential wrt MODYIJ
835:         ! TEMP_GRAD: temporary gradient store643:         ! TEMP_GRAD: temporary gradient store
836:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ644:         DOUBLE PRECISION :: RJ(3), RIJ(3), YIJ(3), MODYIJ
837:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(6)645:         DOUBLE PRECISION :: DVDY, TEMP_GRAD(6)
838: 646: 
839:         ! PER_k: integers for looping over cell possibilities647:         ! PER_k: integers for looping over cell possibilities
840:         ! CELL_RANGE: number of cells in each direction required648:         ! CELL_RANGE: number of cells in each direction required
841:         ! I, J: loop indices649:         ! I, J: loop indices
842:         INTEGER :: PER_X, PER_Y, PER_Z650:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3)
843:         INTEGER :: I, J651:         INTEGER :: I, J
844: 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: 
845:         ! Loop over different coordinate possibilities.658:         ! Loop over different coordinate possibilities.
846:         DO PER_X = - BD%CELL_RANGE(1), BD%CELL_RANGE(1)659:         DO PER_X = - CELL_RANGE(1), CELL_RANGE(1)
847:             DO PER_Y = - BD%CELL_RANGE(2), BD%CELL_RANGE(2)660:             DO PER_Y = - CELL_RANGE(2), CELL_RANGE(2)
848:                 DO PER_Z = - BD%CELL_RANGE(3), BD%CELL_RANGE(3)661:                 DO PER_Z = - CELL_RANGE(3), CELL_RANGE(3)
849: 662: 
850:                     ! Skip the interaction of a particle with itself in the same663:                     ! Skip the interaction of a particle with itself in the same
851:                     ! cell and don't double count interactions within the unit664:                     ! cell and don't double count interactions within the unit
852:                     ! cell.665:                     ! cell.
853:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 &666:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 &
854:                         .AND. IDI .EQ. IDJ) CYCLE667:                         .AND. IDI .EQ. IDJ) CYCLE
855: 668: 
856:                     ! Adjust the j coordinates for the periodic image669:                     ! Adjust the j coordinates for the periodic image
857:                     RJ(1) = COORDS(3*IDJ-2) + PER_X670:                     RJ(1) = COORDS(3*IDJ-2) + PER_X
858:                     RJ(2) = COORDS(3*IDJ-1) + PER_Y671:                     RJ(2) = COORDS(3*IDJ-1) + PER_Y
889:                             ! Lattice parameter derivatives702:                             ! Lattice parameter derivatives
890:                             TEMP_GRAD(:) = 0.D0703:                             TEMP_GRAD(:) = 0.D0
891:                             DO I = 1, 6 ! Loop over the different box parameters704:                             DO I = 1, 6 ! Loop over the different box parameters
892:                                 DO J = 1, 3 ! Loop over x, y, z contributions705:                                 DO J = 1, 3 ! Loop over x, y, z contributions
893:                                     TEMP_GRAD(I) = TEMP_GRAD(I) + &706:                                     TEMP_GRAD(I) = TEMP_GRAD(I) + &
894:                                     YIJ(J)*DOT_PRODUCT(BD%DERIV(J,:,I), RIJ(:))707:                                     YIJ(J)*DOT_PRODUCT(BD%DERIV(J,:,I), RIJ(:))
895:                                 END DO708:                                 END DO
896:                             END DO709:                             END DO
897:                             TEMP_GRAD(:) = TEMP_GRAD(:) * DVDY / MODYIJ710:                             TEMP_GRAD(:) = TEMP_GRAD(:) * DVDY / MODYIJ
898:                             GRAD_BOX(:) = GRAD_BOX(:) + TEMP_GRAD(:)711:                             GRAD_BOX(:) = GRAD_BOX(:) + TEMP_GRAD(:)
899:  
900:                             ! If necessary, calculate the parameter derivatives 
901:                             ! Divide gradient by 2 for images of the same atom, 
902:                             ! to avoid double counting 
903:                             IF(PRESENT(GRAD_PARAMS) .AND. LJ_GAUSS_MODE .EQ. 3)& 
904:                                 GRAD_PARAMS(:) = GRAD_PARAMS(:) + & 
905:                                     MERGE(0.5D0, 1.0D0, IDI .EQ. IDJ) * & 
906:                                     CALC_GRAD_PARAMS(MODYIJ) 
907:  
908:                         END IF ! GTEST712:                         END IF ! GTEST
909:                     END IF ! End if less than cutoff713:                     END IF ! End if less than cutoff
910:                 END DO ! z loop714:                 END DO ! z loop
911:             END DO ! y loop715:             END DO ! y loop
912:         END DO ! x loop716:         END DO ! x loop
913: 717: 
914:     END SUBROUTINE LJ_GAUSS_PER_TRI718:     END SUBROUTINE LJ_GAUSS_PER_TRI
915: 719: 
916: !-------------------------------------------------------------------------------720: !-------------------------------------------------------------------------------
917: !721: !
918: ! Calculates the contribution to the derivatives of the potential wrt the 
919: ! potential parameters LJ_GAUSS_EPSILON and LJ_GAUSS_R0 
920: ! MODYIJ: distance between the particles for this contribution 
921: ! 
922: !------------------------------------------------------------------------------- 
923:  
924:     PURE FUNCTION CALC_GRAD_PARAMS(MODYIJ) 
925:  
926:         IMPLICIT NONE 
927:  
928:         DOUBLE PRECISION :: CALC_GRAD_PARAMS(2) 
929:         DOUBLE PRECISION, INTENT(IN) :: MODYIJ 
930:  
931:         ! DVDE: derivative of pair potential wrt LJ_GAUSS_EPS 
932:         ! DSFDE: derivative of SCL_FCT wrt LJ_GAUSS_EPS 
933:         ! DVDR0: derivative of pair potential wrt LJ_GAUSS_R0 
934:         ! PREF_RIJ, EXP_RIJ: factors for the Guass potential 
935:         DOUBLE PRECISION :: DVDE, DSFDE, DVDR0, PREF_RIJ, EXP_RIJ 
936:  
937:         ! Gauss factors 
938:         PREF_RIJ = (MODYIJ - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ 
939:         EXP_RIJ = LJ_GAUSS_EPS * EXP(- (MODYIJ - LJ_GAUSS_R0)**2 / & 
940:                                        (2.D0*LJ_GAUSS_SIGMASQ)) 
941:  
942:         ! Derivative of the pair potential wrt LJ_GAUSS_EPS just requires 
943:         ! dividing by LJ_GAUSS_EPS 
944:         DVDE = CALC_ENERGY_GAUSS(MODYIJ) / LJ_GAUSS_EPS 
945:  
946:         ! Derivative of the scale factor wrt LJ_GAUSS_EPS 
947:         DSFDE = - SQRT(2 * PI * LJ_GAUSS_SIGMASQ) * SCL_FCT**2 
948:  
949:         ! Now we have everything we need for the LJ_GAUSS_EPS derivative 
950:         CALC_GRAD_PARAMS(1) = SCL_FCT * DVDE + DSFDE * CALC_ENERGY(MODYIJ) 
951:  
952:         ! Next we look at the derivative wrt LJ_GAUSS_R0. This is simpler 
953:         ! because SCL_FCT does not depend on LJ_GAUSS_R0, but DVDR0 is more 
954:         ! complicated than DVDE 
955:         ! Gaussian derivative 
956:         DVDR0 = - PREF_RIJ * EXP_RIJ 
957:         ! Correction to make potential got to zero at RCUT 
958:         DVDR0 = DVDR0 + PREF_RCUT * EXP_RCUT 
959:         ! Correction to make the derivate go to zero at RCUT 
960:         DVDR0 = DVDR0 + (1.D0 / LJ_GAUSS_SIGMASQ - PREF_RCUT**2) * EXP_RCUT * & 
961:                         (MODYIJ - LJ_GAUSS_RCUT) 
962:         ! Correction to make the second derivative go to zero at RCUT 
963:         DVDR0 = DVDR0 - (3.D0 * PREF_RCUT/LJ_GAUSS_SIGMASQ - PREF_RCUT**3) * & 
964:                         EXP_RCUT * 0.5D0 * (MODYIJ - LJ_GAUSS_RCUT)**2 
965:  
966:         ! Now we can do the LJ_GAUSS_R0 derivative 
967:         CALC_GRAD_PARAMS(2) = SCL_FCT * DVDR0 
968:  
969:     END FUNCTION CALC_GRAD_PARAMS 
970:  
971: !------------------------------------------------------------------------------- 
972: ! 
973: ! Impose limits on the allowed range of the potential parameters LJ_GAUSS_R0 and 
974: ! LJ_GAUSS_EPSILON, by adding a harmonic repulsion outside the allowed range. 
975: ! We divide this energy contribution by SCL_FCT (it will be multiplied later) so 
976: ! the constraint is indepenednt of SCL_FCT. 
977: ! ENERGY: energy of the system 
978: ! GRAD_PARAMS: gradients for the potential parameters 
979: ! GTEST: whether to calculate gradients 
980: ! 
981: !------------------------------------------------------------------------------- 
982:     SUBROUTINE CONSTRAIN_PARAMETERS(ENERGY, GRAD_PARAMS, GTEST) 
983:  
984:         IMPLICIT NONE 
985:  
986:         DOUBLE PRECISION, INTENT(INOUT) :: ENERGY, GRAD_PARAMS(2) 
987:         LOGICAL, INTENT(IN) :: GTEST 
988:  
989:         ! EPS_LOWER, EPS_UPPER: bounds on LJ_GAUSS_EPSILON 
990:         ! R0_LOWER, R0_UPPER: bounds on LJ_GAUSS_R0 
991:         ! REPULSION: harmonic force constant 
992:         DOUBLE PRECISION, PARAMETER :: EPS_LOWER = 0.D0, EPS_UPPER = 5.D0 
993:         DOUBLE PRECISION, PARAMETER :: R0_LOWER = 1.D0, R0_UPPER = 2.1D0 
994:         DOUBLE PRECISION, PARAMETER :: REPULSION = 1.D4 
995:  
996:         ! EN_CON: constraint contribution to the energy 
997:         DOUBLE PRECISION :: EN_CON 
998:  
999:         ! Initialise output variable 
1000:         EN_CON = 0.D0 
1001:  
1002:         ! Epsilon energy contribution 
1003:         EN_CON = EN_CON + 0.5D0 * REPULSION * & 
1004:         (MERGE(LJ_GAUSS_EPS-EPS_UPPER, 0.D0, LJ_GAUSS_EPS .GT. EPS_UPPER)**2 + & 
1005:          MERGE(LJ_GAUSS_EPS-EPS_LOWER, 0.D0, LJ_GAUSS_EPS .LT. EPS_LOWER)**2) 
1006:  
1007:         ! R0 energy contribution 
1008:         EN_CON = EN_CON + 0.5D0 * REPULSION * & 
1009:             (MERGE(LJ_GAUSS_R0-R0_UPPER, 0.D0, LJ_GAUSS_R0 .GT. R0_UPPER)**2 + & 
1010:              MERGE(LJ_GAUSS_R0-R0_LOWER, 0.D0, LJ_GAUSS_R0 .LT. R0_LOWER)**2) 
1011:  
1012:         IF (GTEST) THEN 
1013:             ! Epsilon gradient contribution 
1014:             GRAD_PARAMS(1) = GRAD_PARAMS(1) + REPULSION * & 
1015:              (MERGE(LJ_GAUSS_EPS-EPS_UPPER, 0.D0, LJ_GAUSS_EPS .GT. EPS_UPPER)+& 
1016:               MERGE(LJ_GAUSS_EPS-EPS_LOWER, 0.D0, LJ_GAUSS_EPS .LT. EPS_LOWER)) 
1017:             ! R0 gradient contribution 
1018:             GRAD_PARAMS(2) = GRAD_PARAMS(2) + REPULSION * & 
1019:                (MERGE(LJ_GAUSS_R0-R0_UPPER, 0.D0, LJ_GAUSS_R0 .GT. R0_UPPER) + & 
1020:                 MERGE(LJ_GAUSS_R0-R0_LOWER, 0.D0, LJ_GAUSS_R0 .LT. R0_LOWER)) 
1021:         END IF 
1022:  
1023:         ! Scale the contribtion and add it to the energy 
1024:         ENERGY = ENERGY + EN_CON / SCL_FCT 
1025:  
1026:     END SUBROUTINE CONSTRAIN_PARAMETERS 
1027:  
1028: !------------------------------------------------------------------------------- 
1029: ! 
1030: ! Outputs the coordinates to the file ljgauss.xyz722: ! Outputs the coordinates to the file ljgauss.xyz
1031: !723: !
1032: !-------------------------------------------------------------------------------724: !-------------------------------------------------------------------------------
1033:     SUBROUTINE VIEW_LJ_GAUSS()725:     SUBROUTINE VIEW_LJ_GAUSS()
1034: 726: 
1035:         ! NATOMS: number of coordinates727:         ! NATOMS: number of coordinates
1036:         ! NSAVE: number of minima to write728:         ! NSAVE: number of minima to write
1037:         ! BOXx: Box dimensions729:         ! BOXx: Box dimensions
1038:         ! PERIODIC: whether periodic boundary conditions are in use730:         ! PERIODIC: whether periodic boundary conditions are in use
1039:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ731:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ
1046:         IMPLICIT NONE738:         IMPLICIT NONE
1047: 739: 
1048:         ! GETUNIT: function for fiding a free file unit740:         ! GETUNIT: function for fiding a free file unit
1049:         ! FUNIT: output file unit741:         ! FUNIT: output file unit
1050:         ! NMOL: actual number of particles742:         ! NMOL: actual number of particles
1051:         ! I, J: loop iterators743:         ! I, J: loop iterators
1052:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J744:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J
1053: 745: 
1054:         ! BOX_SIZE: convenience array for the box sizes746:         ! BOX_SIZE: convenience array for the box sizes
1055:         ! BOX_ANGLES: convenience array for the box angles747:         ! BOX_ANGLES: convenience array for the box angles
1056:         ! COORD: coordinate in orthogonal space748:         DOUBLE PRECISION :: BOX_SIZE(3), BOX_ANGLES(3)
1057:         DOUBLE PRECISION :: BOX_SIZE(3), BOX_ANGLES(3), COORD(3) 
1058:  
1059:         ! Information to transform from crystallographic to orthogonal space 
1060:         TYPE(BOX_DERIV) :: BD 
1061: 749: 
1062:         ! Set NMOL and box parameters750:         ! Set NMOL and box parameters
1063:         SELECT CASE (LJ_GAUSS_MODE)751:         SELECT CASE (LJ_GAUSS_MODE)
1064:         CASE (0)752:         CASE (0)
1065:             NMOL = NATOMS753:             NMOL = NATOMS
1066:             IF (PERIODIC) THEN754:             IF (PERIODIC) THEN
1067:                 BOX_SIZE(1) = BOXLX755:                 BOX_SIZE(1) = BOXLX
1068:                 BOX_SIZE(2) = BOXLY756:                 BOX_SIZE(2) = BOXLY
1069:                 BOX_SIZE(3) = BOXLZ757:                 BOX_SIZE(3) = BOXLZ
1070:             ELSE758:             ELSE
1071:                 BOX_SIZE(:) = 1.D0759:                 BOX_SIZE(:) = 1.D0
1072:             END IF760:             END IF
1073:             BOX_ANGLES(:) = PI/2.D0761:             BOX_ANGLES(:) = PI/2.D0
1074:         CASE (1)762:         CASE (1)
1075:             NMOL = NATOMS - 1763:             NMOL = NATOMS - 1
1076:             BOX_ANGLES(:) = PI/2.D0764:             BOX_ANGLES(:) = PI/2.D0
1077:         CASE (2)765:         CASE (2)
1078:             NMOL = NATOMS - 2766:             NMOL = NATOMS - 2
1079:         CASE (3) 
1080:             NMOL = NATOMS - 3 
1081:         END SELECT767:         END SELECT
1082: 768: 
1083:         ! Open the file769:         ! Open the file
1084:         FUNIT = GETUNIT()770:         FUNIT = GETUNIT()
1085:         OPEN(UNIT=FUNIT, FILE='ljgauss.xyz', STATUS='UNKNOWN')771:         OPEN(UNIT=FUNIT, FILE='ljgauss.xyz', STATUS='UNKNOWN')
1086: 772: 
1087:         ! Loop over the configurations773:         ! Loop over the configurations
1088:         DO I = 1, NSAVE774:         DO I = 1, NSAVE
1089: 775: 
1090:             WRITE(FUNIT,'(I6)') NMOL776:             WRITE(FUNIT,'(I6)') NMOL
1091: 777: 
1092:             WRITE(FUNIT, '(A,I6,A,F20.10,A,I8,A,F20.10)') 'Energy of minimum ',&778:             WRITE(FUNIT, '(A,I6,A,F20.10,A,I8,A,F20.10)') 'Energy of minimum ',&
1093:                 I, '=', QMIN(I), ' first found at step ', FF(I), &779:                 I, '=', QMIN(I), ' first found at step ', FF(I), &
1094:                 ' Energy per particle = ', QMIN(I)/NMOL780:                 ' Energy per particle = ', QMIN(I)/NMOL
1095: 781: 
1096:             ! Set box parameters for varying box782:             ! Set box parameters for varying box
1097:             IF (LJ_GAUSS_MODE .NE. 0) THEN783:             IF (LJ_GAUSS_MODE .EQ. 1 .OR. LJ_GAUSS_MODE .EQ. 2) THEN
1098:                 BOX_SIZE(:) = QMINP(I, 3*NATOMS-2:3*NATOMS)784:                 BOX_SIZE(:) = QMINP(I, 3*NATOMS-2:3*NATOMS)
1099:             END IF785:             END IF
1100:             IF (LJ_GAUSS_MODE .EQ. 2 .OR. LJ_GAUSS_MODE .EQ. 3) THEN786:             IF (LJ_GAUSS_MODE .EQ. 2) THEN
1101:                 BOX_ANGLES(:) = QMINP(I, 3*NATOMS-5:3*NATOMS-3)787:                 BOX_ANGLES(:) = QMINP(I, 3*NATOMS-5:3*NATOMS-3)
1102:                 BD = CALC_BOX_DERIV(BOX_ANGLES(:), BOX_SIZE(:)) 
1103:             END IF788:             END IF
1104: 789: 
1105:             DO J = 1, NMOL790:             DO J = 1, NMOL
1106: 791: 
1107:                 ! If necessary, transform the coordinate from crytallographic 
1108:                 ! space to orthogonal space 
1109:                 SELECT CASE (LJ_GAUSS_MODE) 
1110:                 CASE (0) 
1111:                     COORD(:) = QMINP(I, 3*J-2:3*J) 
1112:                 CASE (1) 
1113:                     COORD(:) = QMINP(I, 3*J-2:3*J) * BOX_SIZE(:) 
1114:                 CASE (2, 3) 
1115:                     COORD(:) = MATMUL(BD%ORTHOG(:,:), QMINP(I, 3*J-2:3*J)) 
1116:                 END SELECT 
1117:  
1118:                 IF (PERIODIC .AND. J .EQ. 1) THEN792:                 IF (PERIODIC .AND. J .EQ. 1) THEN
1119:                     WRITE(FUNIT, '(A,3F20.10,A,3(A,F20.10))') 'O ', &793:                     WRITE(FUNIT, '(A,3F20.10,A,3(A,F20.10))') 'O ', &
1120:                         COORD(:), ' bbox_xyz', &794:                         QMINP(I, 3*J-2:3*J)*BOX_SIZE(:), ' bbox_xyz', &
1121:                         ' 0.0', BOX_SIZE(1), ' 0.0', BOX_SIZE(2), &795:                         ' 0.0', BOX_SIZE(1), ' 0.0', BOX_SIZE(2), &
1122:                         ' 0.0', BOX_SIZE(3)796:                         ' 0.0', BOX_SIZE(3)
 797:                 ELSE IF (PERIODIC .AND. J .EQ. 2) THEN
 798:                     WRITE(FUNIT, '(A,3F20.10,3(A,F20.10))') 'O ', &
 799:                         QMINP(I, 3*J-2:3*J)*BOX_SIZE(:), ' alpha = ', &
 800:                         BOX_ANGLES(1), ' beta = ', BOX_ANGLES(2), &
 801:                         ' gamma = ', BOX_ANGLES(3)
1123:                 ELSE802:                 ELSE
1124:                     WRITE(FUNIT, '(A,3F20.10)') 'O ', COORD(:)803:                     WRITE(FUNIT, '(A,3F20.10)') 'O ', &
 804:                           QMINP(I, 3*J-2:3*J)*BOX_SIZE(:)*SIN(BOX_ANGLES(:))
1125:                 END IF805:                 END IF
1126: 806: 
1127:             END DO ! Loop over particles807:             END DO ! Loop over particles
1128: 808: 
1129:         END DO ! Loop over the configurations809:         END DO ! Loop over the configurations
1130: 810: 
1131:         CLOSE(FUNIT)811:         CLOSE(FUNIT)
1132: 812: 
1133:     END SUBROUTINE VIEW_LJ_GAUSS813:     END SUBROUTINE VIEW_LJ_GAUSS
1134: 814: 


r32229/modcudalbfgs.F90 2017-03-29 19:30:09.066168151 +0100 r32228/modcudalbfgs.F90 2017-03-29 19:30:09.506173979 +0100
276:     END SUBROUTINE CUDA_LBFGS_WRAPPER276:     END SUBROUTINE CUDA_LBFGS_WRAPPER
277: 277: 
278:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)278:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)
279: 279: 
280:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET280:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET
281: 281: 
282:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY282:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY
283:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS283:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS
284: #ifndef DUMMY_CUDA284: #ifndef DUMMY_CUDA
285:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT285:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT
286: #else /* ifdef DUMMY_CUDA */286: #else ifdef DUMMY_CUDA
287:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT287:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT
288: #endif288: #endif
289:         INTEGER :: X, I, J289:         INTEGER :: X, I, J
290: 290: 
291:         DOUBLE PRECISION :: TOTENERGY291:         DOUBLE PRECISION :: TOTENERGY
292:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS292:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS
293: 293: 
294: #ifndef DUMMY_CUDA294: #ifndef DUMMY_CUDA
295: 295: 
296:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN296:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN
319:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER319:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER
320: 320: 
321:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)321:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)
322:         IMPLICIT NONE322:         IMPLICIT NONE
323: 323: 
324:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET324:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET
325:         REAL(KIND=C_DOUBLE) :: C_ENERGY325:         REAL(KIND=C_DOUBLE) :: C_ENERGY
326:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)326:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)
327: #ifndef DUMMY_CUDA327: #ifndef DUMMY_CUDA
328:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT328:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT
329: #else /* ifdef DUMMY_CUDA */329: #else ifdef DUMMY_CUDA
330:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT330:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT
331: #endif331: #endif
332:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)332:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)
333:         DOUBLE PRECISION    :: DELTA333:         DOUBLE PRECISION    :: DELTA
334:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)334:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)
335:         INTEGER             :: I, J335:         INTEGER             :: I, J
336: 336: 
337: #ifndef DUMMY_CUDA337: #ifndef DUMMY_CUDA
338:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN338:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN
339:             DO J = 1,NADDTARGET339:             DO J = 1,NADDTARGET


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0