hdiff output

r30558/isotropic_potentials.f90 2016-07-06 15:36:57.124327452 +0100 r30557/isotropic_potentials.f90 2016-07-06 15:36:59.336357579 +0100
  1: ! This module provides a library of potential functions for isotropic potentials which return a local copy of the Hessian.  1: ! This module provides a library of potential functions for isotropic potentials which returning a local copy of the Hessian.
  2: ! Here, "isotropic" means that the energy depends only on the distance between particles. 
  3: ! It is designed to be used with the MULTIPOT module  2: ! It is designed to be used with the MULTIPOT module
  4:   3: 
  5: !All subroutines in this module must have the following signature to be used with MULTIPOT:  4: !All subroutines in this module must have the following signature to be used with MULTIPOT:
  6: !        SUBROUTINE POT(TMP_NATOMS, X, PARAMS, TMP_ENERGY, TMP_G, TMP_HESS, GTEST, STEST)  5: !        SUBROUTINE POT(TMP_NATOMS, X, PARAMS, TMP_ENERGY, TMP_G, TMP_HESS, GTEST, STEST)
  7: !            INTEGER, INTENT(IN)           :: TMP_NATOMS  6: !            INTEGER, INTENT(IN)           :: TMP_NATOMS
  8: !            DOUBLE PRECISION, INTENT(IN)  :: X(3*TMP_NATOMS)  7: !            DOUBLE PRECISION, INTENT(IN)  :: X(3*TMP_NATOMS)
  9: !            DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)      ! Maximum number of parameters is hardcoded here  8: !            DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)      ! Maximum number of parameters is hardcoded here
 10: !                                                             ! Each potential will use a subset of the elements of PARAMS  9: !                                                             ! Each potential will use a subset of the elements of PARAMS
 11: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY 10: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_ENERGY
 12: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*TMP_NATOMS), TMP_HESS(3*TMP_NATOMS,3*TMP_NATOMS) 11: !            DOUBLE PRECISION, INTENT(OUT) :: TMP_G(3*TMP_NATOMS), TMP_HESS(3*TMP_NATOMS,3*TMP_NATOMS)


r30558/multipot.f90 2016-07-06 15:36:57.552333240 +0100 r30557/multipot.f90 2016-07-06 15:36:59.692362207 +0100
  1: ! Module that allows different potential types to be specified for different atoms in a single system  1: ! Module that allows different potential types to be specified for different atoms in a single system
  2:  
  3: ! Checklist for adding a new potential: 
  4: !   a) Write a subroutine to calculate the energy (preferably saving it in isotropic_potentials.f90 or pairwise_potentials.f90) 
  5: !   b) Make sure the signature of your subroutine matches the interface in COMPUTE_PAIRWISE_POTENTIAL or COMPUTE_ISOTROPIC_POTENTIAL 
  6: !      (see below). Is neither of these is appropriate for your potential, you will need to write a new COMPUTE subroutine. 
  7: !   c) Add your potential to the SELECT CASE statement in MULTIPOT_CALL 
  8: !   d) If your potential requires atom indices to be input in any format other than a simple list, add your requirements to the 
  9: !      SELECT CASE statement in MULTIPOT_INITIALISE 
 10: !   e) Check that you've made equivalent changes in GMIN so that the two systems stay in sync. 
 11: !   f) Enjoy! 
 12:  
 13: MODULE MULTIPOT  2: MODULE MULTIPOT
 14:   3: 
 15:     USE PAIRWISE_POTENTIALS  4:     USE PAIRWISE_POTENTIALS
 16:     USE ISOTROPIC_POTENTIALS  5:     USE ISOTROPIC_POTENTIALS
 17:     USE COMMONS, ONLY: DEBUG  6:     USE COMMONS, ONLY: DEBUG
 18:   7: 
 19:     IMPLICIT NONE  8:     IMPLICIT NONE
 20:   9: 
 21:     ! We pre-calculate as much as possible, hence the need for lots of different storage arrays, documented here. 10:     ! We pre-calculate as much as possible, hence the need for lots of different storage arrays, documented here.
 22:     INTEGER :: NPOTTYPES = 0                        ! The number of different potential types being used 11:     INTEGER :: NPOTTYPES = 0                        ! The number of different potential types being used
147:                 POTLISTS(J1,2*J2-1,1) = ATOM1  ! Make an entry for each of these atoms in POTLISTS136:                 POTLISTS(J1,2*J2-1,1) = ATOM1  ! Make an entry for each of these atoms in POTLISTS
148:                 POTLISTS(J1,2*J2,1) = ATOM2137:                 POTLISTS(J1,2*J2,1) = ATOM2
149: 138: 
150:                 N_ATOM_PARTNERS(J1,2*J2-1)=1 ! ATOM1 partners with ATOM2...139:                 N_ATOM_PARTNERS(J1,2*J2-1)=1 ! ATOM1 partners with ATOM2...
151:                 POTLISTS(J1,2*J2-1,2) = ATOM2140:                 POTLISTS(J1,2*J2-1,2) = ATOM2
152:                 N_ATOM_PARTNERS(J1,2*J2)=0   ! ...but atom B doesn't partner with any (we've already included its interaction with A)141:                 N_ATOM_PARTNERS(J1,2*J2)=0   ! ...but atom B doesn't partner with any (we've already included its interaction with A)
153:             ENDDO142:             ENDDO
154: 143: 
155:             IF(DEBUG) THEN144:             IF(DEBUG) THEN
156:                 WRITE(*,*) "Potential:", POTTYPES(J1)145:                 WRITE(*,*) "Potential:", POTTYPES(J1)
157:                 WRITE(*,*) "Bonds:"146:                 WRITE(*,*) "N_ATOM_PARTNERS:"
 147:                 WRITE(*,*) N_ATOM_PARTNERS(J1,:NATOM_BY_POT(J1))
 148:                 WRITE(*,*) "POTLISTS:"
158:                 DO J2 = 1,NATOM_BY_POT(J1)149:                 DO J2 = 1,NATOM_BY_POT(J1)
159:                     IF (N_ATOM_PARTNERS(J1,J2).EQ.1) WRITE(*,*) POTLISTS(J1,J2,1:2)150:                     WRITE(*,*) "Atom number", J2, "in this potential"
 151:                     WRITE(*,*) POTLISTS(J1,J2,:N_ATOM_PARTNERS(J1,J2))
160:                 ENDDO152:                 ENDDO
161:             ENDIF153:             ENDIF
162: 154: 
163:         ! For "iso_" potentials. Every specified atom interacts with all the others, but this is done by means of a single155:         ! For "iso_" potentials. Every specified atom interacts with all the others, but this is done by means of a single
164:         ! 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 of156:         ! 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
165:         ! interacting atoms. This is less flexible, but significantly more efficient for large systems.157:         ! interacting atoms. This is less flexible, but significantly more efficient for large systems.
166:         ! The input format is straightforward: simply list all the atoms using this potential on a single line, separated by spaces.158:         ! The input format is straightforward: simply list all the atoms using this potential on a single line, separated by spaces.
167:         ! They don't have to be in index order, but everything will make more sense if they are!159:         ! They don't have to be in index order, but everything will make more sense if they are!
168:         CASE('ILJ','IWCA')160:         CASE('ILJ','IWCA')
169:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))161:             READ(22,*) ATOMLIST(:NATOM_BY_POT(J1))


r30558/pairwise_potentials.f90 2016-07-06 15:36:57.892337834 +0100 r30557/pairwise_potentials.f90 2016-07-06 15:37:00.060367178 +0100
  1: ! This module provides a library of potential functions which operate only on a pair of atoms.  1: ! This module provides a library of potential functions which operate only on a pair of atoms
  2: ! It is designed to be used with the MULTIPOT module  2: ! It is designed to be used with the MULTIPOT module
  3:   3: 
  4: !All subroutines in this module must have the following signature to be used with MULTIPOT:  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)  5: !        SUBROUTINE POT(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST)
  6: !            DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3)  6: !            DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3)
  7: !            DOUBLE PRECISION, INTENT(IN)  :: PARAMS(10)       ! Maximum number of parameters is hardcoded here  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  8: !                                                             ! Each potential will use a subset of the elements of PARAMS
  9: !            DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY  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 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. 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. 12: !                                                             ! STEST is true when calculating the Hessian as well as the other two.
 13: !        END SUBROUTINE POT 13: !        END SUBROUTINE POT
 14:  14: 
 15: MODULE PAIRWISE_POTENTIALS 15: MODULE PAIRWISE_POTENTIALS
 16:  16: 
 17: IMPLICIT NONE 17: IMPLICIT NONE
 18:  18: 
 19: DOUBLE PRECISION, PARAMETER :: WCA_CUT=1.2599210498948731647672106072782283505702514647015079 ! = 2^(1/3) 19: DOUBLE PRECISION, PARAMETER :: WCA_CUT=1.2599210498948731647672106072782283505702514647015079 ! = 2^(1/3)
 20:  20: 
 21: CONTAINS 21: CONTAINS
 22:  22: 
 26:     IMPLICIT NONE 26:     IMPLICIT NONE
 27:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10) 27:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10)
 28:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY 28:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY
 29:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6) 29:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6)
 30:     LOGICAL, INTENT(IN)           :: GTEST, STEST 30:     LOGICAL, INTENT(IN)           :: GTEST, STEST
 31:  31: 
 32:     DOUBLE PRECISION :: R2, R6, R8, R14, SIG, SIG6, SIG12  ! Various powers of the distance between the atoms, and the atom radius 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) 33:         DOUBLE PRECISION :: G, F  ! G tensor and F tensor (see PAIRWISE_GRAD and PAIRWISE_HESSIAN, below)
 34:  34: 
 35:     SIG = PARAMS(1) 35:     SIG = PARAMS(1)
 36:     SIG6 = SIG**6 36:         SIG6 = SIG**6
 37:     SIG12 = SIG6**2 37:         SIG12 = SIG6**2
 38:  38: 
 39:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2 39:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2
 40:     R2 = 1.0D0/R2 40:     R2 = 1.0D0/R2
 41:     R6 = R2**3 41:     R6 = R2**3
 42:  42: 
 43:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0) 43:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0)
 44:  44: 
 45:     IF (.NOT.GTEST) RETURN 45:     IF (.NOT.GTEST) RETURN
 46:  46: 
 47:     G = -24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2*SIG6*R6 47:     G = -24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2*SIG6*R6
 65: ! WCA potential, which is simply the LJ potential but truncated at the first minimum of the function and shifted so that 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. 66: ! the potential is entirely repulsive.
 67: ! The energy and gradient are both continuous at the cutoff. 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. 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) 69: SUBROUTINE PAIRWISE_WCA(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST)
 70:     IMPLICIT NONE 70:     IMPLICIT NONE
 71:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10) ! Maximum number of params is hardcoded here 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 72:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY
 73:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6) 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. 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. 75:                                                                            ! STEST is true when calculating the Hessian as well as the other two.
 76:  76: 
 77:     DOUBLE PRECISION :: R2, R6, R8, R14, SIG, SIG6, SIG12  ! Various powers of the distance between the atoms, and the atom radius 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) 78:     DOUBLE PRECISION :: G, F  ! G tensor and F tensor (see PAIRWISE_GRAD and PAIRWISE_HESSIAN, below)
 79:  79: 
 80:     SIG = PARAMS(1) 80:     SIG = PARAMS(1)
 81:     SIG6 = SIG**6 81:         SIG6 = SIG**6
 82:     SIG12 = SIG6**2 82:         SIG12 = SIG6**2
 83:  83: 
 84:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2 84:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2
 85:  85: 
 86:     IF(R2.GT.WCA_CUT*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file. 86:         IF(R2.GT.WCA_CUT*SIG) THEN   ! WCA_CUT is a PARAMETER - see top of file.
 87:     PAIR_ENERGY = 0.0D0 87:             PAIR_ENERGY = 0.0D0
 88:     PG(:) = 0.0D0 88:             PG(:) = 0.0D0
 89:     P_HESS(:,:) = 0.0D0 89:             P_HESS(:,:) = 0.0D0
 90:     RETURN 90:             RETURN
 91:     ENDIF 91:         ENDIF
 92:  92: 
 93:     R2 = 1.0D0/R2 93:     R2 = 1.0D0/R2
 94:     R6 = R2**3 94:     R6 = R2**3
 95:  95: 
 96:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0) + 1 96:     PAIR_ENERGY=4.0D0*SIG6*R6*(SIG6*R6-1.0D0) + 1
 97:  97: 
 98:     IF (.NOT.GTEST) RETURN 98:     IF (.NOT.GTEST) RETURN
 99:  99: 
100:     G = -24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2*SIG6*R6100:     G = -24.0D0*(2.0D0*SIG6*R6-1.0D0)*R2*SIG6*R6
101: 101: 
102:     CALL PAIRWISE_GRAD(X1, X2, G, PG)102:     CALL PAIRWISE_GRAD(X1, X2, G, PG)
103: 103: 
104:     IF (.NOT.STEST) RETURN104:     IF (.NOT.STEST) RETURN
105: 105: 
106:     R8 = R2**4106:     R8 = R2**4
107:     R14 = R8*R8/R2107:     R14 = R8*R8/R2
108:     F = 672.0D0*SIG12*R14-192.0D0*SIG6*R8108:     F = 672.0D0*SIG12*R14-192.0D0*SIG6*R8
109: 109: 
110:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)110:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)
111: 111: 
112:     RETURN        112:         RETURN        
113: 113: 
114: END SUBROUTINE PAIRWISE_WCA114: END SUBROUTINE PAIRWISE_WCA
115: 115: 
116: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!116: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117: ! Hooke's Law spring potential117: ! 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.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)119: SUBROUTINE HARMONIC_SPRINGS(X1, X2, PARAMS, PG, PAIR_ENERGY, P_HESS, GTEST, STEST)
120:     IMPLICIT NONE120:     IMPLICIT NONE
121:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10)121:     DOUBLE PRECISION, INTENT(IN)  :: X1(3), X2(3), PARAMS(10)
122:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY122:     DOUBLE PRECISION, INTENT(OUT) :: PAIR_ENERGY
123:     DOUBLE PRECISION, INTENT(OUT) :: PG(6), P_HESS(6,6)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.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.125:                                                        ! STEST is true when calculating the Hessian as well as the other two.
126: 126: 
127:         DOUBLE PRECISION :: R0, R, R2, G, F127:         DOUBLE PRECISION :: R0, R, R2, G, F
128: 128: 
129:     R0 = PARAMS(1)129:     R0 = PARAMS(1)
130:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2130:     R2 = (X1(1)-X2(1))**2+(X1(2)-X2(2))**2+(X1(3)-X2(3))**2
131:     R = SQRT(R2)131:         R = SQRT(R2)
132: 132: 
133:     PAIR_ENERGY = 0.5D0*(R-R0)**2133:         PAIR_ENERGY = 0.5D0*(R-R0)**2
134: 134: 
135:     IF(.NOT.GTEST) RETURN135:         IF(.NOT.GTEST) RETURN
136: 136: 
137:     F = R0/R137:         F = R0/R
138:     G = 1.0D0 - F138:         G = 1.0D0 - F
139: 139: 
140:     CALL PAIRWISE_GRAD(X1,X2,G,PG)140:         CALL PAIRWISE_GRAD(X1,X2,G,PG)
141: 141: 
142:     IF (.NOT.STEST) RETURN142:         IF (.NOT.STEST) RETURN
143: 143: 
144:     CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)144:         CALL PAIRWISE_HESSIAN(X1,X2,G,F,R2,P_HESS)
145: 145: 
146:     RETURN146:         RETURN
147: 147: 
148: END SUBROUTINE HARMONIC_SPRINGS148: END SUBROUTINE HARMONIC_SPRINGS
149: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!149: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150: ! END OF PAIRWISE POTENTIALS150: ! END OF PAIRWISE POTENTIALS
151: ! Next, two functions that are general to all isotropic pairwise potentials. To be clear, by "isotropic", I mean151: ! 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"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 the153: ! 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 and154: ! 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 the155: ! 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.156: ! unique value in the matrix.


r30558/rotate_hinge_takestep.f90 2016-07-06 15:36:56.768322637 +0100 r30557/rotate_hinge_takestep.f90 2016-07-06 15:36:58.972352447 +0100
  1: ! A takestep type for the plate-folding potential. 
  2: MODULE HINGE_MOVES  1: MODULE HINGE_MOVES
  3:   2: 
  4:     USE COMMONS, ONLY: DEBUG, MYUNIT  3:     USE COMMONS, ONLY: DEBUG, MYUNIT
  5:   4: 
  6:     IMPLICIT NONE  5:     IMPLICIT NONE
  7:   6: 
  8:     INTEGER :: NHINGES                       ! The number of hinges in the system  7:     INTEGER :: NHINGES                       ! The number of hinges in the system
  9:     INTEGER, ALLOCATABLE :: REFATOMS(:,:)    ! Will hold a set of 4 reference atoms for each hinge  8:     INTEGER, ALLOCATABLE :: REFATOMS(:,:)    ! Will hold a set of 4 reference atoms for each hinge
 10:     INTEGER, ALLOCATABLE :: SET_SIZES(:)     ! A list of the number of plates which need to move for each hinge.  9:     INTEGER, ALLOCATABLE :: SET_SIZES(:)     ! A list of the number of plates which need to move for each hinge.
 11:     INTEGER, ALLOCATABLE :: PLATE_SETS(:,:)  ! For each hinge, this will hold the list of plates which  10:     INTEGER, ALLOCATABLE :: PLATE_SETS(:,:)  ! For each hinge, this will hold the list of plates which 
 12:                                              ! need to be moved when the hinge angle is changed.  11:                                              ! need to be moved when the hinge angle is changed. 
 13:  12: 
 14: CONTAINS 13: CONTAINS
 15:  14: 
  15: 
  16: ! Eventually this will contain code for parsing an input file setting up the hinges.
  17: ! For now, I'm hard-coding in the net for the symmetrical tetrahedron.
 16: SUBROUTINE HINGE_INITIALISE 18: SUBROUTINE HINGE_INITIALISE
 17:  19: 
 18:     IMPLICIT NONE 20:     IMPLICIT NONE
 19:  21: 
 20:     INTEGER :: J1 22:     INTEGER :: J1
 21:  23: 
 22:     ! Input file: hingeconfig 24:     ! Input file: hingeconfig
 23:     ! Format as follows. 25:     ! Format as follows.
 24:     ! 26:     !
 25:     ! 27:     !
 26:     ! The first line contains the number of (flexible) hinges in the system. 28:     ! The first line contains the number of (flexible) hinges in the system.
 27:     ! Each hinge is specified by a set of 3 lines. Sets may be separated by a blank line, for clarity, but this is not compulsory. 29:     ! Each hinge is specified by a set of 3 lines. Sets may be separated by a blank line, for clarity, but this is not compulsory.
 28:     ! 30:     !
 29:     ! We need to know the following: 31:     ! We need to know the following:
 30:     ! - where the hinge is located (we get this by specifying 4 atom indices, 2 of which lie on either side of the hinge) 32:     ! - where the hinge is located (we get this by specifying 4 atom indices, 2 of which lie on either side of the hinge)
 31:     ! - which plates need to move (one side of the hinge is the "moving" side; the plate on that side needs to move, along with all other plates hinged to it) 33:     ! - which plates need to move (one side of the hinge is the "moving" side; the plate on that side needs to move, along with all other plates hinged to it)
 32:     ! 34:     !
 33:     ! The first line of each set contains a number, corresponding to the number of plates on the "moving" side of the hinge. 35:     ! The first line of each set contains a number, corresponding to the number of plates on the "moving" side of the hinge.
 34:     ! It will be more efficient to use the side of the hinge with the smallest number of connected plates as the "moving side" 36:     ! It will be more efficient to use the side of the hinge with the smallest number of connected plates as the "moving side"
 35:     ! 37:     !
 36:     ! The second line of each set contains a list of indices specifying which plates belong to the moving side. These indices must match the order in which the plates are specified in rbodyconfig. 38:     ! The second line of each set contains a list of indices specifying which hinges will move. These indices must match the order in which the plates are specified in rbodyconfig.
 37:     ! 39:     !
 38:     ! The third line contains a list of 4 atom indices. These are the reference atoms for the hinge. They should fall into 2 pairs 40:     ! The third line contains a list of 4 atom indices. These are the reference atoms for the hinge. They should fall into 2 pairs
 39:     ! (listed one after the other), with each pair straddling the hinge and joined by a stiff spring. 41:     ! (listed one after the other), with each pair straddling the hinge and joined by a stiff spring.
 40:     ! The two atoms on each side of the hinge should be connected by a vector parallel to the hinge. 42:     ! The two atoms on one side of the hinge should be connected by a vector parallel to the hinge.
 41:  43: 
 42:  44: 
 43:     OPEN(UNIT=22,FILE='hingeconfig',status='old') 45:     OPEN(UNIT=22,FILE='hingeconfig',status='old')
 44:     READ(22,*) NHINGES 46:     READ(22,*) NHINGES
 45:     IF(DEBUG) WRITE(MYUNIT,*) "Reading hingeconfig. NHINGES:", NHINGES 47:     IF(DEBUG) WRITE(MYUNIT,*) "Reading hingeconfig. NHINGES:", NHINGES
 46:     ALLOCATE(REFATOMS(NHINGES,4))  ! We may eventually want to move to 6 reference atoms so as to compute the angle between the two hinges as well. 48:     ALLOCATE(REFATOMS(NHINGES,4))  ! We may eventually want to move to 6 reference atoms so as to compute the angle between the two hinges as well.
 47:     ALLOCATE(PLATE_SETS(NHINGES,1)) 49:     ALLOCATE(PLATE_SETS(NHINGES,1))
 48:     ALLOCATE(SET_SIZES(NHINGES)) 50:     ALLOCATE(SET_SIZES(NHINGES))
 49:  51: 
 50:     DO J1=1,NHINGES 52:     DO J1=1,NHINGES


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0