hdiff output

r31898/convex_polyhedra.f90 2017-02-09 15:30:17.023396444 +0000 r31897/convex_polyhedra.f90 2017-02-09 15:30:17.875407873 +0000
 29:  29: 
 30:     IMPLICIT none 30:     IMPLICIT none
 31:  31: 
 32:     ! NPOLY: Number of polyhedra 32:     ! NPOLY: Number of polyhedra
 33:     ! X_SIZE: Size of the coordinates array 33:     ! X_SIZE: Size of the coordinates array
 34:     INTEGER :: NPOLY, X_SIZE 34:     INTEGER :: NPOLY, X_SIZE
 35:  35: 
 36:     ! The array of polyhedra 36:     ! The array of polyhedra
 37:     TYPE(POLYHEDRON), ALLOCATABLE :: POLYHEDRA(:) 37:     TYPE(POLYHEDRON), ALLOCATABLE :: POLYHEDRA(:)
 38:  38: 
 39:     ! K_COMPRESS: Potential parameter, compression towards the centre 
 40:     ! K_OVERLAP: Potential parameter, cost of overlap 
 41:     ! MAX_VERT: distance squared of the furthest vertex from the origin 
 42:     DOUBLE PRECISION :: K_COMPRESS, K_OVERLAP, MAX_VERT_SQ 
 43:  
 44:     ! Array of the reference vertex positions 39:     ! Array of the reference vertex positions
 45:     DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: VERTS_0 40:     DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: VERTS_0
 46:  41: 
 47:     PRIVATE :: VERTS_0, NPOLY, X_SIZE, POLYHEDRA, UPDATE_POLYHEDRON 42:     PRIVATE :: VERTS_0, NPOLY, X_SIZE, POLYHEDRA, UPDATE_POLYHEDRON
 48:  43: 
 49: CONTAINS 44: CONTAINS
 50: !------------------------------------------------------------------------------- 45: !-------------------------------------------------------------------------------
 51: ! Initialise the system, including reading the vertex list from the file 46: ! Initialise the system, including reading the vertex list from the file
 52: ! "polyhedron_vertices.dat" 47: ! "polyhedron_vertices.dat"
 53: !------------------------------------------------------------------------------- 48: !-------------------------------------------------------------------------------
 63:         USE GJK_MODULE, ONLY: DIMS, NVERTS 58:         USE GJK_MODULE, ONLY: DIMS, NVERTS
 64:         USE VEC3, ONLY: VEC_CROSS 59:         USE VEC3, ONLY: VEC_CROSS
 65:  60: 
 66:         IMPLICIT NONE 61:         IMPLICIT NONE
 67:  62: 
 68:         ! GETUNIT: Function for finding an unused file unit, from utils.f 63:         ! GETUNIT: Function for finding an unused file unit, from utils.f
 69:         ! VERT_UNIT: File unit for "polyhedron_vertices.dat" 64:         ! VERT_UNIT: File unit for "polyhedron_vertices.dat"
 70:         INTEGER :: GETUNIT, VERT_UNIT, I 65:         INTEGER :: GETUNIT, VERT_UNIT, I
 71:  66: 
 72:         ! NORMAL: Normal to the first face read 67:         ! NORMAL: Normal to the first face read
 73:         ! TEMP_DIST: Distance squared of a vertex from the origin 
 74:         ! TOLERANCE: Numerical tolerance for coplanar test 68:         ! TOLERANCE: Numerical tolerance for coplanar test
 75:         DOUBLE PRECISION :: NORMAL(DIMS), TEMP_DIST 69:         DOUBLE PRECISION :: NORMAL(DIMS)
 76:         DOUBLE PRECISION, PARAMETER :: TOL = 1.0D-6 70:         DOUBLE PRECISION, PARAMETER :: TOL = 1.0D-6
 77:  71: 
 78:         ! COPLANAR: whether all the vertices are coplanar 72:         ! COPLANAR: whether all the vertices are coplanar
 79:         LOGICAL :: COPLANAR 73:         LOGICAL :: COPLANAR
 80:  74: 
 81:         ! Set global module parameters 75:         ! Set global module parameters
 82:         NPOLY = NATOMS / 2 76:         NPOLY = NATOMS / 2
 83:         X_SIZE = NATOMS * 3 77:         X_SIZE = NATOMS * 3
 84:  78: 
 85:         WRITE (MYUNIT, *) 'INITIALISE_POLYHEDRA> ', NPOLY, ' polyhedra' 79:         WRITE (MYUNIT, *) 'INITIALISE_POLYHEDRA> ', NPOLY, ' polyhedra'
 96:         !   VERT(X, 2) VERT(Y, 2) VERT(Z, 2) 90:         !   VERT(X, 2) VERT(Y, 2) VERT(Z, 2)
 97:         !   ... 91:         !   ...
 98:         VERT_UNIT = GETUNIT() 92:         VERT_UNIT = GETUNIT()
 99:         OPEN(UNIT=VERT_UNIT, FILE="polyhedron_vertices.dat", STATUS="old") 93:         OPEN(UNIT=VERT_UNIT, FILE="polyhedron_vertices.dat", STATUS="old")
100:         READ(VERT_UNIT, *) NVERTS 94:         READ(VERT_UNIT, *) NVERTS
101:  95: 
102:         ! Check that we've got at least four vertices 96:         ! Check that we've got at least four vertices
103:         IF (NVERTS .LT. 4) THEN 97:         IF (NVERTS .LT. 4) THEN
104:             WRITE (MYUNIT, *) & 98:             WRITE (MYUNIT, *) &
105:                 'INITIALISE_POLYHEDA> ERROR: require at least 4 vertices' 99:                 'INITIALISE_POLYHEDA> ERROR: require at least 4 vertices'
106:             CLOSE(MYUNIT) 
107:             STOP100:             STOP
108:         END IF101:         END IF
109: 102: 
110:         ! Read the reference vertices103:         ! Read the reference vertices
111:         ALLOCATE(VERTS_0(DIMS, NVERTS))104:         ALLOCATE(VERTS_0(DIMS, NVERTS))
112:         DO I = 1, NVERTS105:         DO I = 1, NVERTS
113:             READ(VERT_UNIT, *) VERTS_0(1,I), VERTS_0(2,I), VERTS_0(3,I)106:             READ(VERT_UNIT, *) VERTS_0(1,I), VERTS_0(2,I), VERTS_0(3,I)
114:         END DO107:         END DO
115: 108: 
116:         ! Check that the vertices are not all coplanar109:         ! Check that the vertices are not all coplanar
122:                 .GT. TOL) THEN115:                 .GT. TOL) THEN
123:                 COPLANAR = .FALSE.116:                 COPLANAR = .FALSE.
124:                 EXIT117:                 EXIT
125:             END IF118:             END IF
126:         END DO ! Loop over vertices119:         END DO ! Loop over vertices
127: 120: 
128:         ! Error out if the vertices are all coplanar121:         ! Error out if the vertices are all coplanar
129:         IF (COPLANAR) THEN122:         IF (COPLANAR) THEN
130:             WRITE (MYUNIT, *) &123:             WRITE (MYUNIT, *) &
131:                 'INITIALISE_POLYHEDRA> ERROR: all vertices are coplanar'124:                 'INITIALISE_POLYHEDRA> ERROR: all vertices are coplanar'
132:             CLOSE(MYUNIT) 
133:             STOP125:             STOP
134:         END IF126:         END IF
135: 127: 
136:         ! Work out the maximum vertex distance 
137:         MAX_VERT_SQ = 0.D0 
138:         DO I = 1, NVERTS 
139:             TEMP_DIST = DOT_PRODUCT(VERTS_0(:, I), VERTS_0(:, I)) 
140:             MAX_VERT_SQ = MAX(MAX_VERT_SQ, TEMP_DIST) 
141:         END DO 
142:  
143:         ! Copy the reference vertices to the polyhedra128:         ! Copy the reference vertices to the polyhedra
144:         DO I = 1, NPOLY129:         DO I = 1, NPOLY
145:             ALLOCATE(POLYHEDRA(I)%VERTS(DIMS, NVERTS))130:             ALLOCATE(POLYHEDRA(I)%VERTS(DIMS, NVERTS))
146:             POLYHEDRA(I)%VERTS(:, :) = VERTS_0(:, :)131:             POLYHEDRA(I)%VERTS(:, :) = VERTS_0(:, :)
147:         END DO132:         END DO
148: 133: 
149:     END SUBROUTINE INITIALISE_POLYHEDRA134:     END SUBROUTINE INITIALISE_POLYHEDRA
150: 135: 
151: !-------------------------------------------------------------------------------136: !-------------------------------------------------------------------------------
152: ! Updates a polyhedron with the rotation matrix and derivatives,137: ! Updates a polyhedron with the rotation matrix and derivatives,
155: ! R: Coordinates of the centre of the particle140: ! R: Coordinates of the centre of the particle
156: ! P: Angle-axis style orientation of the particle141: ! P: Angle-axis style orientation of the particle
157: !    INTENT(INOUT) because orientations are reranged142: !    INTENT(INOUT) because orientations are reranged
158: ! GTEST: Whether to computer derivatives143: ! GTEST: Whether to computer derivatives
159: !-------------------------------------------------------------------------------144: !-------------------------------------------------------------------------------
160:     SUBROUTINE UPDATE_POLYHEDRON(POLY, R, P, GTEST)145:     SUBROUTINE UPDATE_POLYHEDRON(POLY, R, P, GTEST)
161: 146: 
162:         ! DIMS: Dimension of the space147:         ! DIMS: Dimension of the space
163:         ! NVERTS: Number of vertices per polyhedron148:         ! NVERTS: Number of vertices per polyhedron
164:         ! POLYHEDRON: Type for storing a polyhedron149:         ! POLYHEDRON: Type for storing a polyhedron
165:         ! INVERT3X3: Matrix inversion 
166:         USE GJK_MODULE, ONLY: DIMS, NVERTS, POLYHEDRON150:         USE GJK_MODULE, ONLY: DIMS, NVERTS, POLYHEDRON
167:         USE VEC3, ONLY: INVERT3X3 
168: 151: 
169:         IMPLICIT NONE152:         IMPLICIT NONE
170: 153: 
171:         TYPE(POLYHEDRON), INTENT(OUT) :: POLY154:         TYPE(POLYHEDRON), INTENT(OUT) :: POLY
172:         DOUBLE PRECISION, INTENT(IN) :: R(DIMS)155:         DOUBLE PRECISION, INTENT(IN) :: R(DIMS)
173:         DOUBLE PRECISION, INTENT(INOUT) :: P(DIMS)156:         DOUBLE PRECISION, INTENT(INOUT) :: P(DIMS)
174:         LOGICAL, INTENT(IN) :: GTEST157:         LOGICAL, INTENT(IN) :: GTEST
175: 158: 
176:         INTEGER :: I159:         INTEGER :: I
177:         DOUBLE PRECISION :: PI = 4.0D0 * ATAN(1.0D0)160:         DOUBLE PRECISION :: PI = 4.0D0 * ATAN(1.0D0)
179:         DOUBLE PRECISION :: PMOD162:         DOUBLE PRECISION :: PMOD
180: 163: 
181:         POLY%R(:) = R(:)164:         POLY%R(:) = R(:)
182: 165: 
183:         ! Make sure that 0 < |p| < 2*pi166:         ! Make sure that 0 < |p| < 2*pi
184:         PMOD = sqrt(dot_product(p, p))167:         PMOD = sqrt(dot_product(p, p))
185:         IF(PMOD > 2 * PI) P(:) = P(:) / (PMOD * MOD(PMOD, 2 * PI))168:         IF(PMOD > 2 * PI) P(:) = P(:) / (PMOD * MOD(PMOD, 2 * PI))
186:         POLY%P(:) = P(:)169:         POLY%P(:) = P(:)
187: 170: 
188:         ! Compute the rotation matrices and derivatives171:         ! Compute the rotation matrices and derivatives
189:         CALL RMDRVT(POLY%P, POLY%RM, POLY%RMD(:,:,1), POLY%RMD(:,:,2), &172:         CALL RMDRVT(POLY%P, POLY%RM, POLY%RMD(1,:,:), POLY%RMD(2,:,:), &
190:             POLY%RMD(:,:,3), GTEST)173:             POLY%RMD(3,:,:), GTEST)
191:  
192:         ! We need the inverse of the rotation matrix for derivatives 
193:         IF (GTEST) CALL INVERT3X3(POLY%RM(:,:), POLY%INVRM(:,:)) 
194: 174: 
195:         ! Apply rotation matrix and translational offset to all the vertices175:         ! Apply rotation matrix and translational offset to all the vertices
196:         ALLOCATE(POLY%VERTS(DIMS,NVERTS))176:         ALLOCATE(POLY%VERTS(DIMS,NVERTS))
197:         DO I = 1, NVERTS177:         DO I = 1, NVERTS
198:             POLY%VERTS(:, I) = POLY%R(:) + MATMUL(POLY%RM, VERTS_0(:, I))178:             POLY%VERTS(:, I) = POLY%R(:) + MATMUL(POLY%RM, VERTS_0(:, I))
199:         END DO179:         END DO
200: 180: 
201:     END SUBROUTINE UPDATE_POLYHEDRON181:     END SUBROUTINE UPDATE_POLYHEDRON
202: 182: 
203: !-------------------------------------------------------------------------------183: !-------------------------------------------------------------------------------
207: ! GRAD: array of gradients, all translations, then all orientations187: ! GRAD: array of gradients, all translations, then all orientations
208: ! ENERGY: the energy of the system188: ! ENERGY: the energy of the system
209: ! GTEST: whether gradients are required189: ! GTEST: whether gradients are required
210: !-------------------------------------------------------------------------------190: !-------------------------------------------------------------------------------
211:     SUBROUTINE CONVEX_POLYHEDRA(X, GRAD, ENERGY, GTEST)191:     SUBROUTINE CONVEX_POLYHEDRA(X, GRAD, ENERGY, GTEST)
212: 192: 
213:         ! GJK_INTERSECTION: Tests whether two shapes overlap193:         ! GJK_INTERSECTION: Tests whether two shapes overlap
214:         ! EXPANDING_POLYTOP: Finds the intersection information of two shapes194:         ! EXPANDING_POLYTOP: Finds the intersection information of two shapes
215:         ! DIMS: Dimension of the space195:         ! DIMS: Dimension of the space
216:         ! SUPPORT_POINT: Type for storing points in Minkowski difference space196:         ! SUPPORT_POINT: Type for storing points in Minkowski difference space
217:         ! NVERTS: testing only, the number of vertices 
218:         USE GJK_MODULE, ONLY: GJK_INTERSECTION, EXPANDING_POLYTOPE, DIMS, &197:         USE GJK_MODULE, ONLY: GJK_INTERSECTION, EXPANDING_POLYTOPE, DIMS, &
219:                               SUPPORT_POINT !, NVERTS198:                               SUPPORT_POINT
220: 199: 
221:         IMPLICIT NONE200:         IMPLICIT NONE
222: 201: 
223:         DOUBLE PRECISION, INTENT(INOUT) :: X(X_SIZE)202:         DOUBLE PRECISION, INTENT(INOUT) :: X(X_SIZE)
224:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(X_SIZE), ENERGY203:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(X_SIZE), ENERGY
225:         LOGICAL, INTENT(IN) :: GTEST204:         LOGICAL, INTENT(IN) :: GTEST
226: 205: 
227:         ! MOLI_IND, MOLJ_IND: indices for looping over molecules206:         ! MOLI_IND, MOLJ_IND: indices for looping over molecules
228:         ! OFFSET: The end of the position coordinates in X207:         ! OFFSET: The end of the position coordinates in X
229:         INTEGER :: MOLI_IND, MOLJ_IND, OFFSET208:         INTEGER :: MOLI_IND, MOLJ_IND, OFFSET
230: 209: 
231:         ! RESULT_SIMPLEX: A simplex enclosing the origin if GJK finds an210:         ! RESULT_SIMPLEX: A simplex enclosing the origin if GJK finds an
232:         !                 intersection211:         !                 intersection
233:         TYPE(SUPPORT_POINT) :: RESULT_SIMPLEX(DIMS+1)212:         TYPE(SUPPORT_POINT) :: RESULT_SIMPLEX(DIMS+1)
234: 213: 
235:         ! INIT_AXIS: Initial direction guess for GJK, just use the difference of214:         ! INIT_AXIS: Initial direction guess for GJK, just use the difference of
236:         !            the centres215:         !            the centres
237:         ! CENTRE_DIST: distance squared from the centre 
238:         ! OVERLAP_VECT: Normalised vector in the overlap direction, from EPA216:         ! OVERLAP_VECT: Normalised vector in the overlap direction, from EPA
239:         ! OVERLAP_DIST: Penetration distance, from EPA217:         ! OVERLAP_DIST: Penetration distance, from EPA
240:         ! WITNESS1, WITNESS2: The witness points, from EPA218:         ! WITNESS1, WITNESS2: The witness points, from EPA
241:         ! MOLI_REF, MOLJ_REF: the witness points in the reference geometries 
242:         DOUBLE PRECISION :: INIT_AXIS(DIMS)219:         DOUBLE PRECISION :: INIT_AXIS(DIMS)
243:         DOUBLE PRECISION :: CENTRE_DIST, OVERLAP_VECT(DIMS), OVERLAP_DIST220:         DOUBLE PRECISION :: OVERLAP_VECT(DIMS), OVERLAP_DIST
244:         DOUBLE PRECISION :: WITNESS1(DIMS), WITNESS2(DIMS)221:         DOUBLE PRECISION :: WITNESS1(DIMS), WITNESS2(DIMS)
245:         DOUBLE PRECISION :: MOLI_REF(DIMS), MOLJ_REF(DIMS) 
246: 222: 
247:         ! RESULT: whether or not GJK finds an intersection223:         ! RESULT: whether or not GJK finds an intersection
248:         LOGICAL :: RESULT224:         LOGICAL :: RESULT
249: 225: 
250:         ! Testing 
251: !        INTEGER :: GETUNIT, TEST_UNIT, I 
252: !        TEST_UNIT = GETUNIT() 
253: !        OPEN(UNIT=TEST_UNIT, FILE="test_vertices.xyz", STATUS="REPLACE") 
254:  
255:         ! Initialise ouput variables226:         ! Initialise ouput variables
256:         ENERGY = 0.D0227:         ENERGY = 0.D0
257:         IF (GTEST) GRAD(:) = 0.D0228:         IF (GTEST) GRAD(:) = 0.D0
258: 229: 
259:         ! Print out all coordinates, for testing 
260: !        DO MOLI_IND = 1, NPOLY 
261: !            WRITE (*,*) 'Particle ', MOLI_IND, ', trans ', X(MOLI_IND*3-2:MOLI_IND*3) 
262: !            WRITE (*,*) 'Particle ', MOLI_IND, ', rot ', X(MOLI_IND*3-2+OFFSET:MOLI_IND*3+OFFSET) 
263: !        END DO 
264:  
265:         ! Calculate offset from the number of particles230:         ! Calculate offset from the number of particles
266:         OFFSET = 3*NPOLY231:         OFFSET = 3*NPOLY
267: 232: 
268:         ! Update all the particle positions and rotation matrices233:         ! Update all the particle positions and rotation matrices
269:         DO MOLI_IND = 1, NPOLY234:         DO MOLI_IND = 1, NPOLY
270:             CALL UPDATE_POLYHEDRON(POLYHEDRA(MOLI_IND), &235:             CALL UPDATE_POLYHEDRON(POLYHEDRA(MOLI_IND), &
271:                                    X(MOLI_IND*3-2:MOLI_IND*3), &236:                                    X(MOLI_IND*3-2:MOLI_IND*3), &
272:                                    X(OFFSET+MOLI_IND*3-2:OFFSET+MOLI_IND*3), &237:                                    X(OFFSET+MOLI_IND*3-2:OFFSET+MOLI_IND*3), &
273:                                    GTEST)238:                                    GTEST)
274:         END DO239:         END DO
275: 240: 
276:         ! Just looking to test at the moment241:         ! Just looking to test at the moment
277:         ! Testing with two particles, just print whether they overlap or not242:         ! Testing with two particles, just print whether they overlap or not
278: !        INIT_AXIS(:) = POLYHEDRA(1)%R(:) - POLYHEDRA(2)%R(:)243:         INIT_AXIS(:) = POLYHEDRA(1)%R(:) - POLYHEDRA(2)%R(:)
279: !        CALL GJK_INTERSECTION(POLYHEDRA(1), POLYHEDRA(2), INIT_AXIS(:), &244:         CALL GJK_INTERSECTION(POLYHEDRA(1), POLYHEDRA(2), INIT_AXIS(:), &
280: !                              RESULT, RESULT_SIMPLEX(:))245:                               RESULT, RESULT_SIMPLEX(:))
281: 246: 
282: !        WRITE(*,*) 'GJK returned ', RESULT247:         WRITE(*,*) 'GJK returned ', RESULT
283: !        IF (RESULT) THEN248:         IF (RESULT) THEN
284: !            WRITE(*,*) 'Result simplex is:'249:             WRITE(*,*) 'Result simplex is:'
285: !            DO MOLJ_IND = 1, 4250:             DO MOLJ_IND = 1, 4
286: !                WRITE(*,*) RESULT_SIMPLEX(MOLJ_IND)%V(1), &251:                 WRITE(*,*) RESULT_SIMPLEX(MOLJ_IND)%V(1), &
287: !                           RESULT_SIMPLEX(MOLJ_IND)%V(2), &252:                            RESULT_SIMPLEX(MOLJ_IND)%V(2), &
288: !                           RESULT_SIMPLEX(MOLJ_IND)%V(3)253:                            RESULT_SIMPLEX(MOLJ_IND)%V(3)
289: !            END DO254:             END DO
290: !        END IF 
291:  
292: !        IF (RESULT) THEN 
293: !            WRITE(*,*) 'Calling EXPANDING_POLYTOPE' 
294: !            CALL EXPANDING_POLYTOPE(POLYHEDRA(1), POLYHEDRA(2), & 
295: !                                    RESULT_SIMPLEX(:), OVERLAP_DIST, & 
296: !                                    OVERLAP_VECT(:), WITNESS1(:), WITNESS2(:)) 
297: !            WRITE (*,*) 'Overlap distance is ', OVERLAP_DIST 
298: !            WRITE (*,*) 'Overlap vector is ', OVERLAP_VECT(:) 
299: !            WRITE (*,*) 'Witness point on shape 1 is ', WITNESS1(:) 
300: !            WRITE (*,*) 'Witness point on shape 2 is ', WITNESS2(:) 
301: !        END IF 
302:  
303: !        WRITE (*,*) 'Exiting now' 
304: !        STOP 
305:  
306:     ! Harmonic compression towards the origin 
307:     DO MOLI_IND = 1, NPOLY 
308:         CENTRE_DIST = DOT_PRODUCT(POLYHEDRA(MOLI_IND)%R, POLYHEDRA(MOLI_IND)%R) 
309:         ENERGY = ENERGY + 0.5 * K_COMPRESS * CENTRE_DIST 
310:  
311:         IF (GTEST) THEN 
312:             GRAD(MOLI_IND*3-2:MOLI_IND*3) = GRAD(MOLI_IND*3-2:MOLI_IND*3) + & 
313:                                            K_COMPRESS*X(MOLI_IND*3-2:MOLI_IND*3) 
314:         END IF255:         END IF
315:     END DO ! Harmonic compression loop 
316:  
317:     ! Loop over all pairs of particles 
318:     DO MOLI_IND = 1, NPOLY-1 
319:         DO MOLJ_IND = MOLI_IND+1, NPOLY 
320: 256: 
321:             ! Get the vector between the body centres257:         IF (RESULT) THEN
322:             INIT_AXIS(:) = POLYHEDRA(MOLI_IND)%R(:) - POLYHEDRA(MOLJ_IND)%R(:)258: !            WRITE(*,*) 'Calling EXPANDING_POLYTOPE'
323:  
324:             ! If the distance betweent the body centres is greater than twice 
325:             ! the max vertex distance, we can completely skip this pair 
326:             IF (DOT_PRODUCT(INIT_AXIS(:), INIT_AXIS(:)) .GT. 4.D0*MAX_VERT_SQ) & 
327:                 CYCLE 
328:  
329:             ! Otherwise we have to call GJK to see if the shapes intersect 
330:             CALL GJK_INTERSECTION(POLYHEDRA(MOLI_IND), POLYHEDRA(MOLJ_IND), & 
331:                                   INIT_AXIS(:), RESULT, RESULT_SIMPLEX(:)) 
332:  
333:             ! If they don't intersect, we can skip this pair 
334:             IF (.NOT. RESULT) CYCLE 
335:  
336:             ! If they do intersect, we run EPA to get the penetration distance 
337:             ! and witness points. 
338:             CALL EXPANDING_POLYTOPE(POLYHEDRA(1), POLYHEDRA(2), &259:             CALL EXPANDING_POLYTOPE(POLYHEDRA(1), POLYHEDRA(2), &
339:                                     RESULT_SIMPLEX(:), OVERLAP_DIST, &260:                                     RESULT_SIMPLEX(:), OVERLAP_DIST, &
340:                                     OVERLAP_VECT(:), WITNESS1(:), WITNESS2(:))261:                                     OVERLAP_VECT(:), WITNESS1(:), WITNESS2(:))
 262:             WRITE (*,*) 'Overlap distance is ', OVERLAP_DIST
 263:             WRITE (*,*) 'Overlap vector is ', OVERLAP_VECT(:)
 264:             WRITE (*,*) 'Witness point on shape 1 is ', WITNESS1(:)
 265:             WRITE (*,*) 'Witness point on shape 2 is ', WITNESS2(:)
 266:         END IF
341: 267: 
342: !            WRITE (*,*) 'Overlap distance is ', OVERLAP_DIST268:         WRITE (*,*) 'Exiting now'
343: !            WRITE (*,*) 'Overlap vector is ', OVERLAP_VECT(:)269:         STOP
344: !            WRITE (*,*) 'Witness point on shape 1 is ', WITNESS1(:) 
345: !            WRITE (*,*) 'Witness point on shape 2 is ', WITNESS2(:) 
346:  
347: !            WRITE (TEST_UNIT, *) (2*NVERTS + 2) 
348: !            WRITE (TEST_UNIT, *) 'Whatever' 
349: !            DO I = 1, NVERTS 
350: !                WRITE (TEST_UNIT, *) 'C ', POLYHEDRA(MOLI_IND)%VERTS(:, I) 
351: !            END DO 
352: !            DO I = 1, NVERTS 
353: !                WRITE (TEST_UNIT, *) 'N ', POLYHEDRA(MOLJ_IND)%VERTS(:, I) 
354: !            END DO 
355: !            WRITE (TEST_UNIT, *) 'O ', WITNESS1(:) 
356: !            WRITE (TEST_UNIT, *) 'F ', WITNESS2(:) 
357:  
358:             ! Unnormalised overlap vector is more useful 
359:             OVERLAP_VECT(:) = OVERLAP_VECT(:) * OVERLAP_DIST 
360:  
361:             ! Add the harmonic overlap repulsion 
362:             ENERGY = ENERGY + 0.5D0 * K_OVERLAP * OVERLAP_DIST**2 
363:  
364:             IF (GTEST) THEN 
365:                 ! Translational derivatives 
366:                 GRAD(3*MOLI_IND-2:3*MOLI_IND) = GRAD(3*MOLI_IND-2:3*MOLI_IND) & 
367:                                                 + K_OVERLAP * OVERLAP_VECT(:) 
368:                 GRAD(3*MOLJ_IND-2:3*MOLJ_IND) = GRAD(3*MOLJ_IND-2:3*MOLJ_IND) & 
369:                                                 - K_OVERLAP * OVERLAP_VECT(:) 
370:  
371:                 ! Rotational derivatives. We need to get the witness points in 
372:                 ! the reference geometries 
373:                 MOLI_REF(:) = MATMUL(POLYHEDRA(MOLI_IND)%INVRM, WITNESS1(:) - & 
374:                                      POLYHEDRA(MOLI_IND)%R(:)) 
375:                 MOLJ_REF(:) = MATMUL(POLYHEDRA(MOLJ_IND)%INVRM, WITNESS1(:) - & 
376:                                      POLYHEDRA(MOLJ_IND)%R(:)) 
377:  
378: !                WRITE (*,*) 'MOLI_REF = ', MOLI_REF(:) 
379: !                WRITE (*,*) 'MOLJ_REF = ', MOLJ_REF(:) 
380:  
381:                 ! Get the rotational derivatives by taking the dot product of 
382:                 ! the overlap vector with the derivative of the rotation matrix 
383:                 ! acting on the witness point in the reference geometry 
384:                 GRAD(3*MOLI_IND-2+OFFSET) = GRAD(3*MOLI_IND-2+OFFSET) + & 
385:                     K_OVERLAP*DOT_PRODUCT(MATMUL(POLYHEDRA(MOLI_IND)%RMD(:,:,1)& 
386:                     , MOLI_REF(:)), OVERLAP_VECT(:)) 
387:                 GRAD(3*MOLI_IND-1+OFFSET) = GRAD(3*MOLI_IND-1+OFFSET) + & 
388:                     K_OVERLAP*DOT_PRODUCT(MATMUL(POLYHEDRA(MOLI_IND)%RMD(:,:,2)& 
389:                     , MOLI_REF(:)), OVERLAP_VECT(:)) 
390:                 GRAD(3*MOLI_IND + OFFSET) = GRAD(3*MOLI_IND + OFFSET) + & 
391:                     K_OVERLAP*DOT_PRODUCT(MATMUL(POLYHEDRA(MOLI_IND)%RMD(:,:,3)& 
392:                     , MOLI_REF(:)), OVERLAP_VECT(:)) 
393:                 GRAD(3*MOLJ_IND-2+OFFSET) = GRAD(3*MOLJ_IND-2+OFFSET) - & 
394:                     K_OVERLAP*DOT_PRODUCT(MATMUL(POLYHEDRA(MOLJ_IND)%RMD(:,:,1)& 
395:                     , MOLJ_REF(:)), OVERLAP_VECT(:)) 
396:                 GRAD(3*MOLJ_IND-1+OFFSET) = GRAD(3*MOLJ_IND-1+OFFSET) - & 
397:                     K_OVERLAP*DOT_PRODUCT(MATMUL(POLYHEDRA(MOLJ_IND)%RMD(:,:,2)& 
398:                     , MOLJ_REF(:)), OVERLAP_VECT(:)) 
399:                 GRAD(3*MOLJ_IND + OFFSET) = GRAD(3*MOLJ_IND + OFFSET) - & 
400:                     K_OVERLAP*DOT_PRODUCT(MATMUL(POLYHEDRA(MOLJ_IND)%RMD(:,:,3)& 
401:                     , MOLJ_REF(:)), OVERLAP_VECT(:)) 
402:             END IF 
403:  
404:         END DO ! Inner particle loop 
405:     END DO ! Outer particle loop 
406:  
407:     ! Testing 
408: !    CLOSE(TEST_UNIT) 
409: 270: 
410:     END SUBROUTINE CONVEX_POLYHEDRA271:     END SUBROUTINE CONVEX_POLYHEDRA
411: 272: 
412: !------------------------------------------------------------------------------273: !------------------------------------------------------------------------------
413: ! Prints out the vertex information to the file poly.xyz274: ! Prints out the vertex information to the file poly.xyz
414: !------------------------------------------------------------------------------275: !------------------------------------------------------------------------------
415:     SUBROUTINE VIEW_POLYHEDRA()276:     SUBROUTINE VIEW_POLYHEDRA()
416: 277: 
417:         IMPLICIT NONE278:         IMPLICIT NONE
418: 279: 


r31898/gjk.f90 2017-02-09 15:30:17.323400468 +0000 r31897/gjk.f90 2017-02-09 15:30:18.111411039 +0000
 37:         ! V: the point in difference space 37:         ! V: the point in difference space
 38:         ! S1, S2: the points of the individual support functions 38:         ! S1, S2: the points of the individual support functions
 39:         DOUBLE PRECISION :: V(DIMS), S1(DIMS), S2(DIMS) 39:         DOUBLE PRECISION :: V(DIMS), S1(DIMS), S2(DIMS)
 40:     END TYPE SUPPORT_POINT 40:     END TYPE SUPPORT_POINT
 41:  41: 
 42:     ! Stores a single polyhedron 42:     ! Stores a single polyhedron
 43:     TYPE :: POLYHEDRON 43:     TYPE :: POLYHEDRON
 44:         ! R: Coordinates of the centre of the particle 44:         ! R: Coordinates of the centre of the particle
 45:         ! P: Angle-axis style orientation of the particle 45:         ! P: Angle-axis style orientation of the particle
 46:         ! RM: Rotation matrix, derived from P 46:         ! RM: Rotation matrix, derived from P
 47:         ! INVRM: Inverse of RM 
 48:         ! RMD: Derivatives of the rotation matrix 47:         ! RMD: Derivatives of the rotation matrix
 49:         DOUBLE PRECISION :: R(DIMS), P(DIMS), RM(DIMS, DIMS), INVRM(DIMS, DIMS) 48:         DOUBLE PRECISION :: R(DIMS), P(DIMS), RM(DIMS, DIMS)
 50:         DOUBLE PRECISION :: RMD(DIMS, DIMS, DIMS) 49:         DOUBLE PRECISION :: RMD(DIMS, DIMS, DIMS)
 51:         ! VERTS: Vertex positions 50:         ! VERTS: Vertex positions
 52:         DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: VERTS 51:         DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: VERTS
 53:     END TYPE POLYHEDRON 52:     END TYPE POLYHEDRON
 54:  53: 
 55:     ! Stores a triangular face for the EPA algorithm 54:     ! Stores a triangular face for the EPA algorithm
 56:     TYPE :: FACE 55:     TYPE :: FACE
 57:         ! A, B, C: vertex positions of the faces 56:         ! A, B, C: vertex positions of the faces
 58:         ! NORMAL: normal to the face (AB x AC) 57:         ! NORMAL: normal to the face (AB x AC)
 59:         ! ORIG_DIST: distance to the origin 58:         ! ORIG_DIST: distance to the origin
139:             ! Update the simplex and the search direction138:             ! Update the simplex and the search direction
140:             CALL NEAREST_SIMPLEX(WORK_SIMPLEX, LEN_SIMPLEX, &139:             CALL NEAREST_SIMPLEX(WORK_SIMPLEX, LEN_SIMPLEX, &
141:                                  SEARCH_DIRECTION, RESULT)140:                                  SEARCH_DIRECTION, RESULT)
142: 141: 
143:             ! If the simplex contains the origin, we're done142:             ! If the simplex contains the origin, we're done
144:             IF (RESULT) THEN143:             IF (RESULT) THEN
145:                 ! The simplex length should be DIMS+1 (ie a tetrahedron in 3D),144:                 ! The simplex length should be DIMS+1 (ie a tetrahedron in 3D),
146:                 ! or something horrible has happened.145:                 ! or something horrible has happened.
147:                 IF (LEN_SIMPLEX .NE. DIMS + 1) THEN146:                 IF (LEN_SIMPLEX .NE. DIMS + 1) THEN
148:                     WRITE(MYUNIT, *) 'GJK_INTERSECTION> ERROR, length of ', &147:                     WRITE(MYUNIT, *) 'GJK_INTERSECTION> ERROR, length of ', &
149:                     ' simpled is ', LEN_SIMPLEX, ' but dimension is ', DIMS148:                     ' simpled is ', LEN_SIMPLEX, ' but DIMSension is ', DIMS
150:                     CLOSE(MYUNIT) 
151:                     STOP149:                     STOP
152:                 END IF150:                 END IF
153: 151: 
154:                 RESULT_SIMPLEX(:) = WORK_SIMPLEX(:)152:                 RESULT_SIMPLEX(:) = WORK_SIMPLEX(:)
155:                 EXIT153:                 EXIT
156:             END IF154:             END IF
157: 155: 
158:         END DO ! GJK refinement loop156:         END DO ! GJK refinement loop
159: 157: 
160:     END SUBROUTINE GJK_INTERSECTION158:     END SUBROUTINE GJK_INTERSECTION
223:         INTEGER, INTENT(INOUT) :: LEN_SIMPLEX221:         INTEGER, INTENT(INOUT) :: LEN_SIMPLEX
224:         LOGICAL, INTENT(OUT) :: RESULT222:         LOGICAL, INTENT(OUT) :: RESULT
225: 223: 
226:         ! O: the origin224:         ! O: the origin
227:         ! AB, AC, AD, AO, ABPERP, ACPERP, ADPERP, ABCPERP, ACDPERP, ADBPERP:225:         ! AB, AC, AD, AO, ABPERP, ACPERP, ADPERP, ABCPERP, ACDPERP, ADBPERP:
228:         !     storage for plane tests and search direction calculation226:         !     storage for plane tests and search direction calculation
229:         ! TOL: For checking whether things are zero227:         ! TOL: For checking whether things are zero
230:         DOUBLE PRECISION :: AB(DIMS), AC(DIMS), AD(DIMS), AO(DIMS)228:         DOUBLE PRECISION :: AB(DIMS), AC(DIMS), AD(DIMS), AO(DIMS)
231:         DOUBLE PRECISION :: ABPERP(DIMS), ACPERP(DIMS), ADPERP(DIMS)229:         DOUBLE PRECISION :: ABPERP(DIMS), ACPERP(DIMS), ADPERP(DIMS)
232:         DOUBLE PRECISION :: ABCPERP(DIMS), ACDPERP(DIMS), ADBPERP(DIMS)230:         DOUBLE PRECISION :: ABCPERP(DIMS), ACDPERP(DIMS), ADBPERP(DIMS)
233:         DOUBLE PRECISION, PARAMETER :: TOL = 1.D-010231:         DOUBLE PRECISION, PARAMETER :: TOL = 1.D-06
234:         ! OVERABC, OVERACD, OVERADB: Results of plane tests232:         ! OVERABC, OVERACD, OVERADB: Results of plane tests
235:         LOGICAL :: OVERABC, OVERACD, OVERADB233:         LOGICAL :: OVERABC, OVERACD, OVERADB
236: 234: 
237:         ! Initialise output235:         ! Initialise output
238:         SEARCH_DIRECTION(:) = 1.D100236:         SEARCH_DIRECTION(:) = 1.D100
239:         RESULT = .FALSE.237:         RESULT = .FALSE.
240: 238: 
241:         ! The point referred to as A is always the mose recent added to the239:         ! The point referred to as A is always the mose recent added to the
242:         ! simplex, B the second most recent, etc.240:         ! simplex, B the second most recent, etc.
243:         ! No point can be the closest to the origin. In the GJK loop we've241:         ! No point can be the closest to the origin. In the GJK loop we've
382:                 TEST_SIMPLEX(1) = TEST_SIMPLEX(3)380:                 TEST_SIMPLEX(1) = TEST_SIMPLEX(3)
383:                 TEST_SIMPLEX(3) = TEST_SIMPLEX(4)381:                 TEST_SIMPLEX(3) = TEST_SIMPLEX(4)
384:                 ! Check the face and it's edges382:                 ! Check the face and it's edges
385:                 CALL CHECK_FACE(TEST_SIMPLEX, LEN_SIMPLEX, SEARCH_DIRECTION, &383:                 CALL CHECK_FACE(TEST_SIMPLEX, LEN_SIMPLEX, SEARCH_DIRECTION, &
386:                                 AD, AB, ADBPERP)384:                                 AD, AB, ADBPERP)
387:             END IF ! One face case385:             END IF ! One face case
388:         ELSE386:         ELSE
389:             ! Simplex is not the right size, something horrible has happened387:             ! Simplex is not the right size, something horrible has happened
390:             WRITE (MYUNIT, *) 'NEAREST_SIMPLEX> ERROR, simplex size is ', &388:             WRITE (MYUNIT, *) 'NEAREST_SIMPLEX> ERROR, simplex size is ', &
391:                 LEN_SIMPLEX, ' This should not happen.'389:                 LEN_SIMPLEX, ' This should not happen.'
392:             CLOSE(MYUNIT) 
393:             STOP390:             STOP
394:         END IF ! Simplex size391:         END IF ! Simplex size
395: 392: 
396:     END SUBROUTINE NEAREST_SIMPLEX393:     END SUBROUTINE NEAREST_SIMPLEX
397: 394: 
398: !-------------------------------------------------------------------------------395: !-------------------------------------------------------------------------------
399: ! Assuming we know a given face or it's edges are the nearest simplex, find396: ! Assuming we know a given face or it's edges are the nearest simplex, find
400: ! which of the two edges or face it is, and generate the result simplex and397: ! which of the two edges or face it is, and generate the result simplex and
401: ! search direction.398: ! search direction.
402: ! TEST_SIMPLEX: the simplex we're checking. The third point is 'A'399: ! TEST_SIMPLEX: the simplex we're checking. The third point is 'A'
488:         TYPE(EDGE), ALLOCATABLE :: CULL_EDGES(:)485:         TYPE(EDGE), ALLOCATABLE :: CULL_EDGES(:)
489: 486: 
490:         ! NEW_POINT: the new point on the surface of the Minkowski difference487:         ! NEW_POINT: the new point on the surface of the Minkowski difference
491:         TYPE(SUPPORT_POINT) :: NEW_POINT488:         TYPE(SUPPORT_POINT) :: NEW_POINT
492: 489: 
493:         ! CLOSEST_BARY: Barycentric coordinates of the origin projected onto the490:         ! CLOSEST_BARY: Barycentric coordinates of the origin projected onto the
494:         !               closest triangle491:         !               closest triangle
495:         ! MIN_DISTANCE: Keep track of the closest distance to the origin492:         ! MIN_DISTANCE: Keep track of the closest distance to the origin
496:         ! TOL: Tolerance used for deciding if we've moved closer to the origin493:         ! TOL: Tolerance used for deciding if we've moved closer to the origin
497:         DOUBLE PRECISION :: CLOSEST_BARY(DIMS), MIN_DISTANCE494:         DOUBLE PRECISION :: CLOSEST_BARY(DIMS), MIN_DISTANCE
498:         DOUBLE PRECISION, PARAMETER :: TOL = 1.D-10495:         DOUBLE PRECISION, PARAMETER :: TOL = 1.D-06
499: 496: 
500:         ! CLOSEST_INDEX: Index of the closest face to the origin497:         ! CLOSEST_INDEX: Index of the closest face to the origin
501:         ! LEN_CULL: The number of faces to cull498:         ! LEN_CULL: The number of faces to cull
502:         ! MAX_FACES: Current maximum number of faces in array499:         ! MAX_FACES: Current maximum number of faces in array
503:         ! NUM_FACES: Current number of faces in the expanding simplex500:         ! NUM_FACES: Current number of faces in the expanding simplex
504:         ! NUM_NEW: The number of new faces to add501:         ! NUM_NEW: The number of new faces to add
505:         ! CULL_FACES: Which faces need to be removed at each addition502:         ! CULL_FACES: Which faces need to be removed at each addition
506:         ! TEST_UNIT: file unit for testing purposes503:         ! TEST_UNIT: file unit for testing purposes
507:         INTEGER :: CLOSEST_INDEX, LEN_CULL, MAX_FACES, NUM_FACES, NUM_NEW504:         INTEGER :: CLOSEST_INDEX, LEN_CULL, MAX_FACES, NUM_FACES, NUM_NEW
508:         INTEGER, ALLOCATABLE :: CULL_FACES(:)505:         INTEGER, ALLOCATABLE :: CULL_FACES(:)
523: 520: 
524:         ! Start by allocating 20 faces. Will be doubled if required521:         ! Start by allocating 20 faces. Will be doubled if required
525:         MAX_FACES = 20522:         MAX_FACES = 20
526:         ALLOCATE(FACES(MAX_FACES))523:         ALLOCATE(FACES(MAX_FACES))
527: 524: 
528:         ! Add the initial simplex to the array525:         ! Add the initial simplex to the array
529:         NUM_FACES = SIZE(START_SIMPLEX)526:         NUM_FACES = SIZE(START_SIMPLEX)
530:         IF (NUM_FACES .NE. 4) THEN527:         IF (NUM_FACES .NE. 4) THEN
531:             WRITE (MYUNIT, *) 'EXPANDING_POLYTOPE> ERROR, START_SIMPLEX size ',&528:             WRITE (MYUNIT, *) 'EXPANDING_POLYTOPE> ERROR, START_SIMPLEX size ',&
532:                               'is ', NUM_FACES, ' but should be 4'529:                               'is ', NUM_FACES, ' but should be 4'
533:             CLOSE(MYUNIT) 
534:             STOP530:             STOP
535:         END IF531:         END IF
536: 532: 
537:         ! ABC533:         ! ABC
538:         FACES(1) = NEW_FACE(START_SIMPLEX(4),START_SIMPLEX(3),START_SIMPLEX(2))534:         FACES(1) = NEW_FACE(START_SIMPLEX(4),START_SIMPLEX(3),START_SIMPLEX(2))
539:         ! ACD535:         ! ACD
540:         FACES(2) = NEW_FACE(START_SIMPLEX(4),START_SIMPLEX(2),START_SIMPLEX(1))536:         FACES(2) = NEW_FACE(START_SIMPLEX(4),START_SIMPLEX(2),START_SIMPLEX(1))
541:         ! ADB537:         ! ADB
542:         FACES(3) = NEW_FACE(START_SIMPLEX(4),START_SIMPLEX(1),START_SIMPLEX(3))538:         FACES(3) = NEW_FACE(START_SIMPLEX(4),START_SIMPLEX(1),START_SIMPLEX(3))
543:         ! CBD539:         ! CBD
569:                 IF (FACES(I)%ORIG_DIST .LT. MIN_DISTANCE) THEN565:                 IF (FACES(I)%ORIG_DIST .LT. MIN_DISTANCE) THEN
570:                     MIN_DISTANCE = FACES(I)%ORIG_DIST566:                     MIN_DISTANCE = FACES(I)%ORIG_DIST
571:                     CLOSEST_INDEX = I567:                     CLOSEST_INDEX = I
572:                 END IF568:                 END IF
573:             END DO ! Loop over faces569:             END DO ! Loop over faces
574: 570: 
575:             ! Check that we've actually found something571:             ! Check that we've actually found something
576:             IF (CLOSEST_INDEX .EQ. -1) THEN572:             IF (CLOSEST_INDEX .EQ. -1) THEN
577:                 WRITE (MYUNIT, *) 'EXPANDING_POLYTOPE> ERROR, failed to find ',&573:                 WRITE (MYUNIT, *) 'EXPANDING_POLYTOPE> ERROR, failed to find ',&
578:                                   'closest distance'574:                                   'closest distance'
579:                 CLOSE(MYUNIT) 
580:                 STOP575:                 STOP
581:             END IF576:             END IF
582: 577: 
583: !            WRITE (*,*) 'Found minimum distance ', MIN_DISTANCE578: !            WRITE (*,*) 'Found minimum distance ', MIN_DISTANCE
584:             ! Generate the new point on the surface of the Minkowski difference579:             ! Generate the new point on the surface of the Minkowski difference
585:             NEW_POINT%S1 = SUPPORT_FUNC(SHAPE1,  FACES(CLOSEST_INDEX)%NORMAL(:))580:             NEW_POINT%S1 = SUPPORT_FUNC(SHAPE1,  FACES(CLOSEST_INDEX)%NORMAL(:))
586:             NEW_POINT%S2 = SUPPORT_FUNC(SHAPE2, -FACES(CLOSEST_INDEX)%NORMAL(:))581:             NEW_POINT%S2 = SUPPORT_FUNC(SHAPE2, -FACES(CLOSEST_INDEX)%NORMAL(:))
587:             NEW_POINT%V(:) = NEW_POINT%S1(:) - NEW_POINT%S2(:)582:             NEW_POINT%V(:) = NEW_POINT%S1(:) - NEW_POINT%S2(:)
588: !            WRITE (*,*) 'NEW_POINT = ', NEW_POINT%V(1), NEW_POINT%V(2),583: !            WRITE (*,*) 'NEW_POINT = ', NEW_POINT%V(1), NEW_POINT%V(2),
589: !                                        NEW_POINT%V(3)584: !                                        NEW_POINT%V(3)
597:             ! the distance from the origin to the face592:             ! the distance from the origin to the face
598:             IF (DOT_PRODUCT(NEW_POINT%V(:), FACES(CLOSEST_INDEX)%NORMAL(:)) &593:             IF (DOT_PRODUCT(NEW_POINT%V(:), FACES(CLOSEST_INDEX)%NORMAL(:)) &
599:                 - FACES(CLOSEST_INDEX)%ORIG_DIST .LT. TOL) THEN594:                 - FACES(CLOSEST_INDEX)%ORIG_DIST .LT. TOL) THEN
600:                 ! We've found the closest point595:                 ! We've found the closest point
601:                 OVERLAP_DIST = FACES(CLOSEST_INDEX)%ORIG_DIST596:                 OVERLAP_DIST = FACES(CLOSEST_INDEX)%ORIG_DIST
602:                 OVERLAP_VECT = FACES(CLOSEST_INDEX)%NORMAL(:)597:                 OVERLAP_VECT = FACES(CLOSEST_INDEX)%NORMAL(:)
603: 598: 
604:                 ! Find the barycentric coordinates of the origin projected onto599:                 ! Find the barycentric coordinates of the origin projected onto
605:                 ! the triangle600:                 ! the triangle
606:                 CLOSEST_BARY(:) = BARYCENTRIC(FACES(CLOSEST_INDEX))601:                 CLOSEST_BARY(:) = BARYCENTRIC(FACES(CLOSEST_INDEX))
607: !                WRITE(*,*) 'CLOSEST_BARY = ', CLOSEST_BARY(:)602:                 WRITE(*,*) 'CLOSEST_BARY = ', CLOSEST_BARY(:)
608: !                WRITE(*,*) 'A%S1 = ', FACES(CLOSEST_INDEX)%A%S1(:)603:                 WRITE(*,*) 'A%S1 = ', FACES(CLOSEST_INDEX)%A%S1(:)
609: !                WRITE(*,*) 'A%S1 = ', FACES(CLOSEST_INDEX)%B%S1(:)604:                 WRITE(*,*) 'A%S1 = ', FACES(CLOSEST_INDEX)%B%S1(:)
610: !                WRITE(*,*) 'A%S1 = ', FACES(CLOSEST_INDEX)%C%S1(:)605:                 WRITE(*,*) 'A%S1 = ', FACES(CLOSEST_INDEX)%C%S1(:)
611: 606: 
612:                 ! Calculate the witness points607:                 ! Calculate the witness points
613:                 WITNESS1(:) = CLOSEST_BARY(1)*FACES(CLOSEST_INDEX)%A%S1(:) + &608:                 WITNESS1(:) = CLOSEST_BARY(1)*FACES(CLOSEST_INDEX)%A%S1(:) + &
614:                               CLOSEST_BARY(2)*FACES(CLOSEST_INDEX)%B%S1(:) + &609:                               CLOSEST_BARY(2)*FACES(CLOSEST_INDEX)%B%S1(:) + &
615:                               CLOSEST_BARY(3)*FACES(CLOSEST_INDEX)%C%S1(:)610:                               CLOSEST_BARY(3)*FACES(CLOSEST_INDEX)%C%S1(:)
616: 611: 
617:                 WITNESS2(:) = CLOSEST_BARY(1)*FACES(CLOSEST_INDEX)%A%S2(:) + &612:                 WITNESS2(:) = CLOSEST_BARY(1)*FACES(CLOSEST_INDEX)%A%S2(:) + &
618:                               CLOSEST_BARY(2)*FACES(CLOSEST_INDEX)%B%S2(:) + &613:                               CLOSEST_BARY(2)*FACES(CLOSEST_INDEX)%B%S2(:) + &
619:                               CLOSEST_BARY(3)*FACES(CLOSEST_INDEX)%C%S2(:)614:                               CLOSEST_BARY(3)*FACES(CLOSEST_INDEX)%C%S2(:)
620:                 EXIT615:                 EXIT
687:             IF (NUM_FACES - LEN_CULL + NUM_NEW .GT. MAX_FACES) THEN682:             IF (NUM_FACES - LEN_CULL + NUM_NEW .GT. MAX_FACES) THEN
688:                 ! Reallocate the array683:                 ! Reallocate the array
689:                 ALLOCATE(TEMP_FACES(NUM_FACES))684:                 ALLOCATE(TEMP_FACES(NUM_FACES))
690:                 TEMP_FACES(:) = FACES(1:NUM_FACES)685:                 TEMP_FACES(:) = FACES(1:NUM_FACES)
691:                 DEALLOCATE(FACES)686:                 DEALLOCATE(FACES)
692:                 MAX_FACES = MAX_FACES*2687:                 MAX_FACES = MAX_FACES*2
693:                 ALLOCATE(FACES(MAX_FACES))688:                 ALLOCATE(FACES(MAX_FACES))
694:                 FACES(1:NUM_FACES) = TEMP_FACES(:)689:                 FACES(1:NUM_FACES) = TEMP_FACES(:)
695:                 DEALLOCATE(TEMP_FACES)690:                 DEALLOCATE(TEMP_FACES)
696:                 ! Stop here for testing691:                 ! Stop here for testing
697:                 !STOP692:                 STOP
698:             END IF693:             END IF
699: 694: 
700:             ! Add all the new faces695:             ! Add all the new faces
701:             ! Keep track of which edge we're on with J696:             ! Keep track of which edge we're on with J
702:             J = 1697:             J = 1
703:             DO I = 1, LEN_CULL698:             DO I = 1, LEN_CULL
704:                 ! Start by replacing deleted faces699:                 ! Start by replacing deleted faces
705:                 DO WHILE (.NOT. INCLUDE_EDGES(J))700:                 DO WHILE (.NOT. INCLUDE_EDGES(J))
706:                     J = J + 1701:                     J = J + 1
707:                 END DO702:                 END DO
817: ! Checks whether two edges are the same as each other (regardless of direction)812: ! Checks whether two edges are the same as each other (regardless of direction)
818: ! EDGE1, EDGE2: The two edges to check813: ! EDGE1, EDGE2: The two edges to check
819: !-------------------------------------------------------------------------------814: !-------------------------------------------------------------------------------
820:     LOGICAL FUNCTION CHECK_EDGE(EDGE1, EDGE2)815:     LOGICAL FUNCTION CHECK_EDGE(EDGE1, EDGE2)
821: 816: 
822:         IMPLICIT NONE817:         IMPLICIT NONE
823: 818: 
824:         TYPE(EDGE), INTENT(IN) :: EDGE1, EDGE2819:         TYPE(EDGE), INTENT(IN) :: EDGE1, EDGE2
825: 820: 
826:         ! TOL: numerical tolerance for testing equality821:         ! TOL: numerical tolerance for testing equality
827:         DOUBLE PRECISION, PARAMETER :: TOL = 1.D-10822:         DOUBLE PRECISION, PARAMETER :: TOL = 1.D-06
828: 823: 
829:         ! DIFFAA, DIFFBB, DIFFAB, DIFFBA: differences between the vertices824:         ! DIFFAA, DIFFBB, DIFFAB, DIFFBA: differences between the vertices
830:         DOUBLE PRECISION :: DIFFAA(DIMS), DIFFBB(DIMS)825:         DOUBLE PRECISION :: DIFFAA(DIMS), DIFFBB(DIMS)
831:         DOUBLE PRECISION :: DIFFAB(DIMS), DIFFBA(DIMS)826:         DOUBLE PRECISION :: DIFFAB(DIMS), DIFFBA(DIMS)
832: 827: 
833:         DIFFAA(:) = ABS(EDGE1%A%V(:) - EDGE2%A%V(:))828:         DIFFAA(:) = ABS(EDGE1%A%V(:) - EDGE2%A%V(:))
834:         DIFFBB(:) = ABS(EDGE1%B%V(:) - EDGE2%B%V(:))829:         DIFFBB(:) = ABS(EDGE1%B%V(:) - EDGE2%B%V(:))
835:         DIFFAB(:) = ABS(EDGE1%A%V(:) - EDGE2%B%V(:))830:         DIFFAB(:) = ABS(EDGE1%A%V(:) - EDGE2%B%V(:))
836:         DIFFBA(:) = ABS(EDGE1%B%V(:) - EDGE2%A%V(:))831:         DIFFBA(:) = ABS(EDGE1%B%V(:) - EDGE2%A%V(:))
837: 832: 


r31898/keywords.f 2017-02-09 15:30:17.619404439 +0000 r31897/keywords.f 2017-02-09 15:30:18.403414957 +0000
 32:      &                      atomsinclist, atomsinalistlogical, atomsinblistlogical, atomsinclistlogical, ligcartstep, 32:      &                      atomsinclist, atomsinalistlogical, atomsinblistlogical, atomsinclistlogical, ligcartstep,
 33:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange, 33:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange,
 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT, 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT,
 35:      &                      SALTCON, macroiont, nmacroions, macroiondist 35:      &                      SALTCON, macroiont, nmacroions, macroiondist
 36:       USE modamber 36:       USE modamber
 37:       USE PORFUNCS 37:       USE PORFUNCS
 38:       USE MYGA_PARAMS 38:       USE MYGA_PARAMS
 39:       USE BGUPMOD 39:       USE BGUPMOD
 40:       USE GLJYMOD 40:       USE GLJYMOD
 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA
 43:       USE MBPOLMOD, ONLY: MBPOLINIT 43:       USE MBPOLMOD, ONLY: MBPOLINIT
 44:       USE SWMOD, ONLY: SWINIT, MWINIT 44:       USE SWMOD, ONLY: SWINIT, MWINIT
 45:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS, 45:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS,
 46:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA 46:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA
 47:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 47:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 48:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 48:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 49:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,  49:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 
 50:      &     PARSE_MLJ_PARAMS 50:      &     PARSE_MLJ_PARAMS
 51:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT 51:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT
 52:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE 52:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE
6595: 6595: 
6596:          CALL DEFPAPJANUS()6596:          CALL DEFPAPJANUS()
6597: 6597: 
6598:          PAPJANT = .TRUE.6598:          PAPJANT = .TRUE.
6599:          RIGID   = .TRUE.6599:          RIGID   = .TRUE.
6600: 6600: 
6601:       ELSE IF (WORD .EQ. 'POLYHEDRA') THEN6601:       ELSE IF (WORD .EQ. 'POLYHEDRA') THEN
6602: 6602: 
6603:         POLYT = .TRUE.6603:         POLYT = .TRUE.
6604:         RIGID = .TRUE.6604:         RIGID = .TRUE.
6605:         CALL READF(K_OVERLAP) 
6606:         CALL READF(K_COMPRESS) 
6607:         CALL INITIALISE_POLYHEDRA()6605:         CALL INITIALISE_POLYHEDRA()
6608: 6606: 
6609:       ELSE IF (WORD .EQ. 'PTSTST') THEN6607:       ELSE IF (WORD .EQ. 'PTSTST') THEN
6610: 6608: 
6611:          CALL READI(NPATCH)6609:          CALL READI(NPATCH)
6612:          CALL READF(PAPEPS)6610:          CALL READF(PAPEPS)
6613:          CALL READF(PAPCD)6611:          CALL READF(PAPCD)
6614:          CALL READF(YKAPPA)6612:          CALL READF(YKAPPA)
6615: 6613: 
6616:          NRBSITES = NPATCH6614:          NRBSITES = NPATCH


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0