hdiff output

r29559/genrigid.f90 2015-12-03 17:30:13.783256919 +0000 r29558/genrigid.f90 2015-12-03 17:30:14.935272047 +0000
 30:       LOGICAL :: DUMMYLOGICAL = .FALSE. 30:       LOGICAL :: DUMMYLOGICAL = .FALSE.
 31:       INTEGER, ALLOCATABLE :: LRBNPERMGROUP(:), LRBNPERMSIZE(:,:), LRBPERMGROUP(:,:), LRBNSETS(:,:), LRBSETS(:,:,:) 31:       INTEGER, ALLOCATABLE :: LRBNPERMGROUP(:), LRBNPERMSIZE(:,:), LRBPERMGROUP(:,:), LRBNSETS(:,:), LRBSETS(:,:,:)
 32:  32: 
 33: !   vr274:  added lattice coordinates 33: !   vr274:  added lattice coordinates
 34: !           if HAS_LATTICE_COORDS is true, the last two atoms are treated 34: !           if HAS_LATTICE_COORDS is true, the last two atoms are treated
 35: !           as lattice coordintes and rigidcoords is in reduced lattice units 35: !           as lattice coordintes and rigidcoords is in reduced lattice units
 36:       LOGICAL HAS_LATTICE_COORDS 36:       LOGICAL HAS_LATTICE_COORDS
 37: !   csw34>  RIGIDISRIGID = logical array, size NATOMS, TRUE if atom is part of 37: !   csw34>  RIGIDISRIGID = logical array, size NATOMS, TRUE if atom is part of
 38: !           RB 38: !           RB
 39:       LOGICAL, ALLOCATABLE :: RIGIDISRIGID(:) 39:       LOGICAL, ALLOCATABLE :: RIGIDISRIGID(:)
 40:       INTEGER, ALLOCATABLE :: RB_BY_ATOM(:) ! sn402: records which RB an atom belongs to 
 41:  40: 
 42: !-----------------------------------------------------------------------------------! 41: !-----------------------------------------------------------------------------------!
 43: ! NRIGIDBODY  = number of rigid bodies 42: ! NRIGIDBODY  = number of rigid bodies
 44: ! DEGFREEDOMS = number of degrees of freedom = 6 * NRIGIDBODY + 3 * ADDITIONAL ATOMS 43: ! DEGFREEDOMS = number of degrees of freedom = 6 * NRIGIDBODY + 3 * ADDITIONAL ATOMS
 45: ! MAXSITE     = maximum number of sites in a rigid body 44: ! MAXSITE     = maximum number of sites in a rigid body
 46: ! NRELAXRIGIDR = rigid body minimisation for this number of steps 45: ! NRELAXRIGIDR = rigid body minimisation for this number of steps
 47: ! NRELAXRIGIDA = atom minimisation for this number of steps 46: ! NRELAXRIGIDA = atom minimisation for this number of steps
 48: ! NSITEPERBODY= number of rigid body sites, no need to be the same for all bodies 47: ! NSITEPERBODY= number of rigid body sites, no need to be the same for all bodies
 49: ! REFVECTOR   = reference vector for the atomistic to rigic coordinate transformation (used for checking only) 48: ! REFVECTOR   = reference vector for the atomistic to rigic coordinate transformation
 50: ! RIGIDSINGLES= list of atoms not in rigid bodies 49: ! RIGIDSINGLES= list of atoms not in rigid bodies
 51: ! RIGIDGROUPS = list of atoms in rigid bodies, need a file called rbodyconfig 50: ! RIGIDGROUPS = list of atoms in rigid bodies, need a file called rbodyconfig
 52: ! RIGIDCOORDS = 6 * NRIGIDBODY + 3 * ADDITIONAL ATOMS coordinates 51: ! RIGIDCOORDS = 6 * NRIGIDBODY + 3 * ADDITIONAL ATOMS coordinates
 53: ! SITESRIGIDBODY = coordinates of the rigid body sites 52: ! SITESRIGIDBODY = coordinates of the rigid body sites
 54: ! RIGIDINIT   = logical variable for generalised rigid body 53: ! RIGIDINIT   = logical variable for generalised rigid body
 55: ! ATOMRIGIDCOORDT, .TRUE. = atom coords active, .FALSE. = rigid coords active, used in mylbfgs & potential 54: ! ATOMRIGIDCOORDT, .TRUE. = atom coords active, .FALSE. = rigid coords active, used in mylbfgs & potential
 56: ! GENRIGIDT = generalised rigid body takestep taken if .TRUE. (GMIN only) 55: ! GENRIGIDT = generalised rigid body takestep taken if .TRUE.
 57: !-----------------------------------------------------------------------------------! 56: !-----------------------------------------------------------------------------------!
 58:  57: 
 59:  58: 
 60: CONTAINS 59: CONTAINS
 61: ! vr274> TODO: better initialization of SITESRIDIGBODY, why use maxsite? 60: ! vr274> TODO: better initialization of SITESRIDIGBODY, why use maxsite?
 62:  61: 
 63: !------------------------------------------- 62: !-------------------------------------------
 64: ! vr274> Initializes basic structures 63: ! vr274> Initializes basic structures
 65: ! hast to be the first call to GENRIGID function in order to setup basic structures. 64: ! hast to be the first call to GENRIGID function in order to setup basic structures.
 66: ! After that, the array which defines the sites can be filled. Then GENRIGID_INITIALIZE 65: ! After that, the array which defines the sites can be filled. Then GENRIGID_INITIALIZE
 81:  80: 
 82:   ! hk286 > Allocate NSITEPERBODY 81:   ! hk286 > Allocate NSITEPERBODY
 83:   ALLOCATE (NSITEPERBODY(NRIGIDBODY)) 82:   ALLOCATE (NSITEPERBODY(NRIGIDBODY))
 84:   ALLOCATE (SITESRIGIDBODY(MAXSITE,3,NRIGIDBODY)) 83:   ALLOCATE (SITESRIGIDBODY(MAXSITE,3,NRIGIDBODY))
 85:   ALLOCATE (RIGIDGROUPS(MAXSITE,NRIGIDBODY)) 84:   ALLOCATE (RIGIDGROUPS(MAXSITE,NRIGIDBODY))
 86:   ALLOCATE (REFVECTOR(NRIGIDBODY)) 85:   ALLOCATE (REFVECTOR(NRIGIDBODY))
 87:   ALLOCATE (GR_WEIGHTS(NATOMS)) 86:   ALLOCATE (GR_WEIGHTS(NATOMS))
 88:   ALLOCATE (IINVERSE(NRIGIDBODY,3,3)) 87:   ALLOCATE (IINVERSE(NRIGIDBODY,3,3))
 89: ! csw34>  88: ! csw34> 
 90:   ALLOCATE (RIGIDISRIGID(NATOMS)) 89:   ALLOCATE (RIGIDISRIGID(NATOMS))
 91:   ALLOCATE (RB_BY_ATOM(NATOMS)) 
 92:   RIGIDISRIGID=.FALSE. 90:   RIGIDISRIGID=.FALSE.
 93:   RB_BY_ATOM(:) = -1 
 94:  91: 
 95:   ! by default use center of geometry - check if ok for AMBER 9/8/12 92:   ! by default use center of geometry - check if ok for AMBER 9/8/12
 96:   IF ( ALLOCATED(ATMASS1) ) THEN 93:   IF ( ALLOCATED(ATMASS1) ) THEN
 97:      GR_WEIGHTS=ATMASS1 94:      GR_WEIGHTS=ATMASS1
 98:   ELSE IF ( ALLOCATED(AMBER12_MASSES) .AND. AMBER12T ) THEN 95:   ELSE IF ( ALLOCATED(AMBER12_MASSES) .AND. AMBER12T ) THEN
 99:      ! AMBER 12 atom masses 96:      ! AMBER 12 atom masses
100:      GR_WEIGHTS = AMBER12_MASSES 97:      GR_WEIGHTS = AMBER12_MASSES
101:   ELSE IF (MASSES) THEN 98:   ELSE IF (MASSES) THEN
102:      GR_WEIGHTS=ATMASS 99:      GR_WEIGHTS=ATMASS
103:   ELSE100:   ELSE
346:      READ(222,*) CHECK1, NSITEPERBODY(J1)343:      READ(222,*) CHECK1, NSITEPERBODY(J1)
347:      DO J2 = 1, NSITEPERBODY(J1)344:      DO J2 = 1, NSITEPERBODY(J1)
348:         READ(222,*) RIGIDGROUPS(J2,J1)345:         READ(222,*) RIGIDGROUPS(J2,J1)
349: ! csw34> check to make sure the current atom is not already in a rigid body346: ! csw34> check to make sure the current atom is not already in a rigid body
350:         IF (RIGIDISRIGID(RIGIDGROUPS(J2,J1))) THEN347:         IF (RIGIDISRIGID(RIGIDGROUPS(J2,J1))) THEN
351:             PRINT *," genrigid> ERROR: atom ",RIGIDGROUPS(J2,J1)," is in multiple rigid bodies! Stopping."348:             PRINT *," genrigid> ERROR: atom ",RIGIDGROUPS(J2,J1)," is in multiple rigid bodies! Stopping."
352:             STOP349:             STOP
353:         ELSE       350:         ELSE       
354: ! csw34> if not, flag the current atom351: ! csw34> if not, flag the current atom
355:             RIGIDISRIGID(RIGIDGROUPS(J2,J1))=.TRUE.352:             RIGIDISRIGID(RIGIDGROUPS(J2,J1))=.TRUE.
356:             RB_BY_ATOM(RIGIDGROUPS(J2,J1)) = J1 
357:         ENDIF353:         ENDIF
358: ! vr274> Moved initialization of coordinates to GENRIGID_INITIALISE, here only read the setup354: ! vr274> Moved initialization of coordinates to GENRIGID_INITIALISE, here only read the setup
359: !        SITESRIGIDBODY(J2,:,J1) = COORDS(3*DUMMY-2:3*DUMMY,1)355: !        SITESRIGIDBODY(J2,:,J1) = COORDS(3*DUMMY-2:3*DUMMY,1)
360:      ENDDO356:      ENDDO
361:   ENDDO357:   ENDDO
362:   CLOSE(222)358:   CLOSE(222)
363:   CALL GENRIGID_INITIALISE(INICOORDS)359:   CALL GENRIGID_INITIALISE(INICOORDS)
364: END SUBROUTINE GENRIGID_READ_FROM_FILE360: END SUBROUTINE GENRIGID_READ_FROM_FILE
365: 361: 
366: !-----------------------------------------------------------362: !-----------------------------------------------------------
418: ENDDO414: ENDDO
419: 415: 
420: RETURN416: RETURN
421: END SUBROUTINE SETUP_TENSORS417: END SUBROUTINE SETUP_TENSORS
422: !-----------------------------------------------------------418: !-----------------------------------------------------------
423: 419: 
424: SUBROUTINE TRANSFORMRIGIDTOC (CMIN, CMAX, XCOORDS, XRIGIDCOORDS)420: SUBROUTINE TRANSFORMRIGIDTOC (CMIN, CMAX, XCOORDS, XRIGIDCOORDS)
425:       421:       
426:   USE COMMONS, ONLY: NATOMS422:   USE COMMONS, ONLY: NATOMS
427:   IMPLICIT NONE423:   IMPLICIT NONE
 424:   
428:   INTEGER :: J1, J2, J5, J7, J9, NP        !NP = No of processors425:   INTEGER :: J1, J2, J5, J7, J9, NP        !NP = No of processors
429:   INTEGER :: CMIN, CMAX426:   INTEGER :: CMIN, CMAX
430:   DOUBLE PRECISION :: P(3), RMI(3,3), DRMI1(3,3), DRMI2(3,3), DRMI3(3,3)427:   DOUBLE PRECISION :: P(3), RMI(3,3), DRMI1(3,3), DRMI2(3,3), DRMI3(3,3)
431:   DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(3*NATOMS)428:   DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(3*NATOMS)
432:   DOUBLE PRECISION :: COM(3) ! center of mass429:   DOUBLE PRECISION :: COM(3) ! center of mass
433:   LOGICAL          :: GTEST !, ATOMTEST430:   LOGICAL          :: GTEST !, ATOMTEST
434:   DOUBLE PRECISION :: MLATTICE(3,3)431:   DOUBLE PRECISION :: MLATTICE(3,3)
435:   432:   
436:   GTEST = .FALSE.433:   GTEST = .FALSE.
437:   NP = 1434:   NP = 1
438: ! vr274 > are there additional lattice coordinates? If yes, setup transformation matrix435: ! vr274 > are there additional lattice coordinates? If yes, setup transformation matrix
439:   IF(HAS_LATTICE_COORDS) THEN436:   IF(HAS_LATTICE_COORDS) THEN
440:     CALL GET_LATTICE_MATRIX(XRIGIDCOORDS(DEGFREEDOMS-5:DEGFREEDOMS), MLATTICE)437:     CALL GET_LATTICE_MATRIX(XRIGIDCOORDS(DEGFREEDOMS-5:DEGFREEDOMS), MLATTICE)
441:   ELSE ! vr274 > otherwise identity matrix438:   ELSE ! vr274 > otherwise identity matrix
442:     MLATTICE = 0D0439:     MLATTICE = 0D0
443:     MLATTICE(1,1)=1d0440:     MLATTICE(1,1)=1d0
444:     MLATTICE(2,2)=1D0441:     MLATTICE(2,2)=1D0
445:     MLATTICE(3,3)=1D0442:     MLATTICE(3,3)=1D0
446:   ENDIF443:   ENDIF
447: 444: 
 445: 
448:   ! hk286 > coord transformations for rigid bodies CMIN to CMAX446:   ! hk286 > coord transformations for rigid bodies CMIN to CMAX
449:   DO J1 = CMIN, CMAX447:   DO J1 = CMIN, CMAX
450:      J5   = 3*J1448:      J5   = 3*J1
451:      J7   = 3*NRIGIDBODY + J5449:      J7   = 3*NRIGIDBODY + J5
452:      P(:) = XRIGIDCOORDS(J7-2:J7)450:      P(:) = XRIGIDCOORDS(J7-2:J7)
453:      CALL RMDRVT(P, RMI, DRMI1, DRMI2, DRMI3, GTEST)451:      CALL RMDRVT(P, RMI, DRMI1, DRMI2, DRMI3, GTEST)
454: 452: 
455: ! vr274 > MLATTICE can have lattice transformation or be identity matrix453: ! vr274 > MLATTICE can have lattice transformation or be identity matrix
456:      COM = matmul(MLATTICE, XRIGIDCOORDS(J5-2:J5))454:      COM = matmul(MLATTICE, XRIGIDCOORDS(J5-2:J5))
457:      DO J2 = 1,  NSITEPERBODY(J1)455:      DO J2 = 1,  NSITEPERBODY(J1)
458:         J9 = RIGIDGROUPS(J2, J1)456:         J9 = RIGIDGROUPS(J2, J1)
459:         XCOORDS(3*J9-2:3*J9) = COM + MATMUL(RMI(:,:),SITESRIGIDBODY(J2,:,J1))457:         XCOORDS(3*J9-2:3*J9) = COM + MATMUL(RMI(:,:),SITESRIGIDBODY(J2,:,J1))
460:      ENDDO458:      ENDDO
 459:      
461:   ENDDO460:   ENDDO
462:   461:   
463: ! hk286 > now the single atoms462: ! hk286 > now the single atoms
464: ! vr274 > this copies lattice coordinates as well which is stored in last 2 atoms463: ! vr274 > this copies lattice coordinates as well which is stored in last 2 atoms
465:   IF (DEGFREEDOMS > 6 * NRIGIDBODY) THEN464:   IF (DEGFREEDOMS > 6 * NRIGIDBODY) THEN
466:      DO J1 = 1, (DEGFREEDOMS - 6*NRIGIDBODY)/3465:      DO J1 = 1, (DEGFREEDOMS - 6*NRIGIDBODY)/3
467:         J9 = RIGIDSINGLES(J1)466:         J9 = RIGIDSINGLES(J1)
468:         XCOORDS(3*J9-2:3*J9) = XRIGIDCOORDS(6*NRIGIDBODY + 3*J1-2:6*NRIGIDBODY + 3*J1)467:         XCOORDS(3*J9-2:3*J9) = XRIGIDCOORDS(6*NRIGIDBODY + 3*J1-2:6*NRIGIDBODY + 3*J1)
469:      ENDDO468:      ENDDO
470:   ENDIF469:   ENDIF
 470: !      WRITE(*,*) "Finally"
 471: !      WRITE(*,*) XCOORDS
471: END SUBROUTINE TRANSFORMRIGIDTOC472: END SUBROUTINE TRANSFORMRIGIDTOC
472: 473: 
473: !----------------------------------------------------------474: !----------------------------------------------------------
474: 475: 
475: SUBROUTINE ROTATEINITIALREF ()476: SUBROUTINE ROTATEINITIALREF ()
476: IMPLICIT NONE477: IMPLICIT NONE
477: DOUBLE PRECISION :: P(3)478: DOUBLE PRECISION :: P(3)
478: INTEGER J1479: INTEGER J1
479: 480: 
480: ! hk286 - rotate the system - new481: ! hk286 - rotate the system - new


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


r29559/key.f90 2015-12-03 17:30:14.167261962 +0000 r29558/key.f90 2015-12-03 17:30:15.315277038 +0000
 44:      &        CHECKCONINT, ATOMMATCHFULL, NIH2LEPST,NIHEAM7T,NIHLEPST, NIHPAIRONLYT, & 44:      &        CHECKCONINT, ATOMMATCHFULL, NIH2LEPST,NIHEAM7T,NIHLEPST, NIHPAIRONLYT, &
 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, & 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, &
 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, & 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, &
 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST
 55:  55: 
 56: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 56: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 57:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 57:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 58:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 58:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 59:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 59:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 60:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 60:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
 61:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads 61:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads
 62:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs 62:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs
 63:  63: 
 64: ! hk286 > generalised THOMSON problem 64: ! hk286 > generalised THOMSON problem


r29559/keywords.f 2015-12-03 17:30:14.359264483 +0000 r29558/keywords.f 2015-12-03 17:30:15.511279611 +0000
 50:          USE INTCOMMONS, ONLY : NATINT, INTNEWT, BBCART, INTINTERPT, INTERPSIMPLE, 50:          USE INTCOMMONS, ONLY : NATINT, INTNEWT, BBCART, INTINTERPT, INTERPSIMPLE,
 51:      $   INTMINPERMT, INTERPCHOICE, NINTIM, CARTRESSTART, INTPARFILE, 51:      $   INTMINPERMT, INTERPCHOICE, NINTIM, CARTRESSTART, INTPARFILE,
 52:      $   MINBACKTCUT, INTERPBACKTCUT, PRINTCOORDS, DESMINT, NURINGS, URINGS, 52:      $   MINBACKTCUT, INTERPBACKTCUT, PRINTCOORDS, DESMINT, NURINGS, URINGS,
 53:      $   NUBONDS, UBONDS, USEPARFILE, CHICDNEB, OLDINTMINPERMT, 53:      $   NUBONDS, UBONDS, USEPARFILE, CHICDNEB, OLDINTMINPERMT,
 54:      $   GLYCART, INTDISTANCET, RIGIDBONDS, GMAXINT, BONDSFROMFILE 54:      $   GLYCART, INTDISTANCET, RIGIDBONDS, GMAXINT, BONDSFROMFILE
 55:          ! MCP 55:          ! MCP
 56:          USE AMHGLOBALS 56:          USE AMHGLOBALS
 57:          USE BGUPMOD 57:          USE BGUPMOD
 58:          ! hk286 58:          ! hk286
 59:          USE GENRIGID 59:          USE GENRIGID
 60:          USE MULTIPOT 
 61:          ! includes toggle for switching frames - required for RBAA 60:          ! includes toggle for switching frames - required for RBAA
 62:          USE MODHESS, ONLY : RBAANORMALMODET 61:          USE MODHESS, ONLY : RBAANORMALMODET
 63:          ! jdf43> for MMEINITWRAPPER 62:          ! jdf43> for MMEINITWRAPPER
 64:          USE ISO_C_BINDING, ONLY: C_NULL_CHAR 63:          USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 65: ! AMBER12 64: ! AMBER12
 66:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ATOM, AMBER12_RESIDUE, AMBER12_ATOMS, 65:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ATOM, AMBER12_RESIDUE, AMBER12_ATOMS,
 67:      &                                    AMBER12_RESIDUES, AMBER12_GET_COORDS 66:      &                                    AMBER12_RESIDUES, AMBER12_GET_COORDS
 68:          USE CHIRALITY 67:          USE CHIRALITY
 69:  68: 
 70:          IMPLICIT NONE 69:          IMPLICIT NONE
162:          161:          
163:          !ds656> Initialise NTYPEA here162:          !ds656> Initialise NTYPEA here
164:          NTYPEA=NATOMS163:          NTYPEA=NATOMS
165:          !164:          !
166:          ! hk286 - initialise to stationary frame165:          ! hk286 - initialise to stationary frame
167:          RBAANORMALMODET = .FALSE.166:          RBAANORMALMODET = .FALSE.
168:          ! generalised rigid body167:          ! generalised rigid body
169:          ATOMRIGIDCOORDT = .TRUE.168:          ATOMRIGIDCOORDT = .TRUE.
170:          RIGIDINIT = .FALSE.169:          RIGIDINIT = .FALSE.
171:          AACONVERGENCET = .FALSE.170:          AACONVERGENCET = .FALSE.
172:  
173:          ! Multiple potential scheme 
174:          MULTIPOTT= .FALSE. 
175:  
176:          ! Thomson problem171:          ! Thomson problem
177:          GTHOMSONT = .FALSE.172:          GTHOMSONT = .FALSE.
178:          GTHOMPOT = 1173:          GTHOMPOT = 1
179: 174: 
180:          DESMAXEJUMP = HUGE(1.0D0)175:          DESMAXEJUMP = HUGE(1.0D0)
181:          DESMAXAVGE = HUGE(1.0D0)176:          DESMAXAVGE = HUGE(1.0D0)
182:          UNSTRING='UNDEFINED'177:          UNSTRING='UNDEFINED'
183:          WELCH=.FALSE.178:          WELCH=.FALSE.
184:          TOSI=.FALSE.179:          TOSI=.FALSE.
185:          TOSIC6=.FALSE.180:          TOSIC6=.FALSE.
3739:             PRINT '(A)',' keywords> Multiple jobs will be run, reading additional starting coordinates from file '3734:             PRINT '(A)',' keywords> Multiple jobs will be run, reading additional starting coordinates from file '
3740:      &      // TRIM(ADJUSTL(MULTISTART))3735:      &      // TRIM(ADJUSTL(MULTISTART))
3741:             MULTIFUNIT=-13736:             MULTIFUNIT=-1
3742:             IF (NITEMS.GT.2) THEN3737:             IF (NITEMS.GT.2) THEN
3743:                CALL READA(MULTIFINISH)3738:                CALL READA(MULTIFINISH)
3744:                MULTIFUNIT=GETUNIT()3739:                MULTIFUNIT=GETUNIT()
3745:                OPEN(MULTIFUNIT,FILE=TRIM(ADJUSTL(MULTIFINISH)),STATUS='OLD')3740:                OPEN(MULTIFUNIT,FILE=TRIM(ADJUSTL(MULTIFINISH)),STATUS='OLD')
3746:                PRINT '(A)',' keywords>                            reading additional finish coordinates from file '3741:                PRINT '(A)',' keywords>                            reading additional finish coordinates from file '
3747:      &         // TRIM(ADJUSTL(MULTIFINISH))3742:      &         // TRIM(ADJUSTL(MULTIFINISH))
3748:             ENDIF3743:             ENDIF
3749:  
3750:              ELSE IF (WORD .EQ. 'MULTIPOT') THEN ! Activate the multiple-potential scheme 
3751:              
3752:             MULTIPOTT = .TRUE. 
3753:             CALL MULTIPOT_INITIALISE 
3754:  
3755:          ELSE IF (WORD .EQ. 'MULTISITEPY') THEN3744:          ELSE IF (WORD .EQ. 'MULTISITEPY') THEN
3756: ! Syntax: MULTISITEPY sig_0 eps_0 [cut] [XYZ boxx boxy boxz]3745: ! Syntax: MULTISITEPY sig_0 eps_0 [cut] [XYZ boxx boxy boxz]
3757: ! Notes: The cutoff length is the raw length. It is not scaled3746: ! Notes: The cutoff length is the raw length. It is not scaled
3758: ! by the PY sigma_0 since it is also used for the LJ potential.3747: ! by the PY sigma_0 since it is also used for the LJ potential.
3759: ! The box length is in units of cutoff distance, so it should be3748: ! The box length is in units of cutoff distance, so it should be
3760: ! >= 2.3749: ! >= 2.
3761: 3750: 
3762:             MULTISITEPYT = .TRUE.3751:             MULTISITEPYT = .TRUE.
3763: ! ANGLEAXIS2 = .TRUE.3752: ! ANGLEAXIS2 = .TRUE.
3764:             RBAAT = .TRUE.3753:             RBAAT = .TRUE.


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


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


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0