hdiff output

r31861/commons.f90 2017-01-30 19:30:24.958336832 +0000 r31860/commons.f90 2017-01-30 19:30:26.982364221 +0000
102:      &        MSGBT, MSTBINT, MMRSDPT, MSSTOCKT, LWOTPT, NCAPT, NPAHT, NTIPT, PAHW99T, ELLIPSOIDT, GAYBERNET,& 102:      &        MSGBT, MSTBINT, MMRSDPT, MSSTOCKT, LWOTPT, NCAPT, NPAHT, NTIPT, PAHW99T, ELLIPSOIDT, GAYBERNET,& 
103:      &        MULTPAHAT, PAPT, PAPBINT, PAPJANT, PTSTSTT, SHIFTED, SILANET, TDHDT, DDMT, WATERDCT, WATERKZT, CHECKDT, CHECKMARKOVT,& 103:      &        MULTPAHAT, PAPT, PAPBINT, PAPJANT, PTSTSTT, SHIFTED, SILANET, TDHDT, DDMT, WATERDCT, WATERKZT, CHECKDT, CHECKMARKOVT,& 
104:      &        TETHER, HISTSMOOTH, VISITPROP, ARMT, FixedEndMoveT, FIXCOM, RESTORET, QUADT, AMHT, MOVESHELLT, QDT, QD2T, &104:      &        TETHER, HISTSMOOTH, VISITPROP, ARMT, FixedEndMoveT, FIXCOM, RESTORET, QUADT, AMHT, MOVESHELLT, QDT, QD2T, &
105:      &        PARAMONOVPBCX, PARAMONOVPBCY, PARAMONOVPBCZ, PARAMONOVCUTOFF, GAYBERNEDCT, UNFREEZERES, FREEZEALL, &105:      &        PARAMONOVPBCX, PARAMONOVPBCY, PARAMONOVPBCZ, PARAMONOVCUTOFF, GAYBERNEDCT, UNFREEZERES, FREEZEALL, &
106:      &        PROJIT, PROJIHT, LEAPDIHE, DUMPQUT, DUMPBESTT, LJSITE, BLJSITE, LJSITEATTR, DUMPSTEPST, &106:      &        PROJIT, PROJIHT, LEAPDIHE, DUMPQUT, DUMPBESTT, LJSITE, BLJSITE, LJSITEATTR, DUMPSTEPST, &
107:      &        AMBERT, RANDOMSEEDT, PYGPERIODICT, LJCAPSIDT, PYBINARYT, SWAPMOVEST, MOVABLEATOMST, LIGMOVET,DUMPSTRUCTURES, &107:      &        AMBERT, RANDOMSEEDT, PYGPERIODICT, LJCAPSIDT, PYBINARYT, SWAPMOVEST, MOVABLEATOMST, LIGMOVET,DUMPSTRUCTURES, &
108:      &        LJSITECOORDST, VGW, ACKLANDT, G46, DF1T, PULLT, LOCALSAMPLET, CSMT, A9INTET, INTERESTORE, COLDFUSION, &108:      &        LJSITECOORDST, VGW, ACKLANDT, G46, DF1T, PULLT, LOCALSAMPLET, CSMT, A9INTET, INTERESTORE, COLDFUSION, &
109:      &        CSMGUIDET, MULTISITEPYT, CHAPERONINT, AVOIDRESEEDT, OHCELLT, UNFREEZEFINALQ, PERCOLATET, PERCT, PERCACCEPTED, PERCCOMPMARKOV, &109:      &        CSMGUIDET, MULTISITEPYT, CHAPERONINT, AVOIDRESEEDT, OHCELLT, UNFREEZEFINALQ, PERCOLATET, PERCT, PERCACCEPTED, PERCCOMPMARKOV, &
110:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF,&110:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF,&
111:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&111:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &
113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
121:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &121:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &
122:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS122:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS


r31861/convex_polyhedra.f90 2017-01-30 19:30:25.262340947 +0000 r31860/convex_polyhedra.f90 2017-01-30 19:30:27.306368605 +0000
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/convex_polyhedra.f90' in revision 31860
  2: !   Copyright (C) 1999-2017 David J. Wales 
  3: !   This file is part of GMIN. 
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: ! Potential for convex polyhedra. The polyhedra only interact when they are. In 
 20: ! this case, there is a harmonic repulsion based on the smallest distance 
 21: ! required to make them not overlap. Additionally, there is a harmonic 
 22: ! compression towards the origin on all particles, to promote dense packing. 
 23: ! Questions to jwrm2. 
 24:  
 25: MODULE CONVEX_POLYHEDRA_MODULE 
 26:  
 27: ! POLYHEDRON: type to store a polyhedron 
 28:     USE GJK_UTILS, ONLY: POLYHEDRON 
 29:  
 30:     IMPLICIT none 
 31:  
 32:     ! NPOLY: Number of polyhedra 
 33:     ! NVERTS: Number of vertices in each polyhedron 
 34:     ! X_SIZE: Size of the coordinates array 
 35:     INTEGER :: NPOLY, X_SIZE 
 36:  
 37:     ! The array of polyhedra 
 38:     TYPE(POLYHEDRON), ALLOCATABLE :: POLYHEDRA(:) 
 39:  
 40: CONTAINS 
 41: !------------------------------------------------------------------------------- 
 42: ! Initialise the system, including reading the vertex list from the file 
 43: ! "polyhedron_vertices.dat" 
 44: !------------------------------------------------------------------------------- 
 45:     SUBROUTINE INITIALISE_POLYHEDRA() 
 46:  
 47:         ! MYUNIT: File unit of GMIN_out 
 48:         ! NATOMS: the number of coordinate lines in the coords file, so twice 
 49:         !         the actual number of bodies 
 50:         ! DIM: Number of dimensions of the system 
 51:         ! NVERTS: the number of vertices per polyhedron 
 52:         ! VEC_CROSS: finds the cross product of two vectors 
 53:         USE COMMONS, ONLY: MYUNIT, NATOMS 
 54:         USE GJK_UTILS, ONLY: DIM, NVERTS 
 55:         USE VEC3, ONLY: VEC_CROSS 
 56:  
 57:         IMPLICIT NONE 
 58:  
 59:         ! GETUNIT: Function for finding an unused file unit, from utils.f 
 60:         ! VERT_UNIT: File unit for "polyhedron_vertices.dat" 
 61:         INTEGER :: GETUNIT, VERT_UNIT, I 
 62:  
 63:         ! TEMP_VERTS: Temporary vertex information 
 64:         ! NORMAL: Normal to the first face read 
 65:         ! TOLERANCE: Numerical tolerance for coplanar test 
 66:         DOUBLE PRECISION, ALLOCATABLE :: TEMP_VERTS(:,:) 
 67:         DOUBLE PRECISION :: NORMAL(3) 
 68:         DOUBLE PRECISION, PARAMETER :: TOL = 1.0D-6 
 69:  
 70:         ! COPLANAR: whether all the vertices are coplanar 
 71:         LOGICAL :: COPLANAR 
 72:  
 73:         ! Set global module parameters 
 74:         NPOLY = NATOMS / 2 
 75:         X_SIZE = NATOMS * 3 
 76:  
 77:         WRITE (MYUNIT, *) 'INITIALISE_POLYHEDRA> ', NPOLY, ' polyhedra' 
 78:  
 79:         ! Allocate array of polyhedra 
 80:         ! We don't need to worry about the positions and orientations here, 
 81:         ! they'll be dealt with the first time the potential is called. 
 82:         ALLOCATE(POLYHEDRA(NPOLY)) 
 83:  
 84:         ! Read the reference vertex positions from file. 
 85:         ! Format is 
 86:         !   NVERT 
 87:         !   VERT(X, 1) VERT(Y, 1) VERT(Z, 1) 
 88:         !   VERT(X, 2) VERT(Y, 2) VERT(Z, 2) 
 89:         !   ... 
 90:         VERT_UNIT = GETUNIT() 
 91:         OPEN(UNIT=VERT_UNIT, FILE="polyhedron_vertices.dat", STATUS="old") 
 92:         READ(VERT_UNIT, *) NVERTS 
 93:  
 94:         ! Check that we've got at least four vertices 
 95:         IF (NVERTS .LT. 4) THEN 
 96:             WRITE (MYUNIT, *) & 
 97:                 'INITIALISE_POLYHEDA> ERROR: require at least 4 vertices' 
 98:             STOP 
 99:         END IF 
100:  
101:         ! Read the reference vertices 
102:         ALLOCATE(TEMP_VERTS(DIM, NVERTS)) 
103:         DO I = 1, NVERTS 
104:             READ(VERT_UNIT, *) TEMP_VERTS(1,I), TEMP_VERTS(2,I), TEMP_VERTS(3,I) 
105:         END DO 
106:  
107:         ! Check that the vertices are not all coplanar 
108:         NORMAL = VEC_CROSS((TEMP_VERTS(:, 1) - TEMP_VERTS(:, 2)), & 
109:                                (TEMP_VERTS(:, 1) - TEMP_VERTS(:, 3))) 
110:         COPLANAR = .TRUE. 
111:         DO I = 4, NVERTS ! Loop over vertices 
112:             IF (DOT_PRODUCT((TEMP_VERTS(:, I) - TEMP_VERTS(:, 1)), NORMAL) & 
113:                 .GT. TOL) THEN 
114:                 COPLANAR = .FALSE. 
115:                 EXIT 
116:             END IF 
117:         END DO ! Loop over vertices 
118:  
119:         ! Error out if the vertices are all coplanar 
120:         IF (COPLANAR) THEN 
121:             WRITE (MYUNIT, *) & 
122:                 'INITIALISE_POLYHEDRA> ERROR: all vertices are coplanar' 
123:             STOP 
124:         END IF 
125:  
126:         ! Copy the reference vertices to the polyhedra 
127:         DO I = 1, NPOLY 
128:             ALLOCATE(POLYHEDRA(NPOLY)%VERTS(DIM, NVERTS)) 
129:             POLYHEDRA(NPOLY)%VERTS(:, :) = TEMP_VERTS(:, :) 
130:         END DO 
131:  
132:     END SUBROUTINE INITIALISE_POLYHEDRA 
133:  
134: !------------------------------------------------------------------------------- 
135: ! Updates a polyhedron with the rotation matrix and derivatives, 
136: ! and the position of the vertices. 
137: ! POLY: the polyhedron to operate on 
138: ! R: Coordinates of the centre of the particle 
139: ! P: Angle-axis style orientation of the particle 
140: !    INTENT(INOUT) because orientations are reranged 
141: ! GTEST: Whether to computer derivatives 
142: !------------------------------------------------------------------------------- 
143:     SUBROUTINE UPDATE_POLYHEDRON(POLY, R, P, GTEST) 
144:  
145:         ! DIM: Dimension of the space 
146:         ! POLYHEDRON: Type for storing a polyhedron 
147:         USE GJK_UTILS, ONLY: DIM, POLYHEDRON 
148:  
149:         IMPLICIT NONE 
150:  
151:         TYPE(POLYHEDRON), INTENT(OUT) :: POLY 
152:         DOUBLE PRECISION, INTENT(IN) :: R(DIM) 
153:         DOUBLE PRECISION, INTENT(INOUT) :: P(DIM) 
154:         LOGICAL, INTENT(IN) :: GTEST 
155:  
156:         DOUBLE PRECISION :: PI = 4.0D0 * ATAN(1.0D0) 
157:         ! PMOD: angle of the angle axis orientation 
158:         DOUBLE PRECISION :: PMOD 
159:  
160:         POLY%R(:) = R(:) 
161:  
162:         ! Make sure that 0 < |p| < 2*pi 
163:         PMOD = sqrt(dot_product(p, p)) 
164:         IF(PMOD > 2 * PI) P(:) = P(:) / (PMOD * MOD(PMOD, 2 * PI)) 
165:         POLY%P(:) = P(:) 
166:  
167:         ! Compute the rotation matrices and derivatives 
168:         CALL RMDRVT(POLY%P, POLY%RM, POLY%RMD(1,:,:), POLY%RMD(2,:,:), & 
169:             POLY%RMD(3,:,:), GTEST) 
170:  
171:     END SUBROUTINE UPDATE_POLYHEDRON 
172:  
173: !------------------------------------------------------------------------------- 
174: ! Calculates the potential, and optionally the gradients 
175: ! X: array of coordinates, all translations, then all orientations 
176: !    INTENT(INOUT) because orientation may be reranged 
177: ! GRAD: array of gradients, all translations, then all orientations 
178: ! ENERGY: the energy of the system 
179: ! GTEST: whether gradients are required 
180: !------------------------------------------------------------------------------- 
181:     SUBROUTINE CONVEX_POLYHEDRA(X, GRAD, ENERGY, GTEST) 
182:  
183:         ! GJK_INTERSECTION: Tests whether two shapes overlap 
184:         ! DIM: Dimension of the space 
185:         USE GJK_MODULE, ONLY: GJK_INTERSECTION 
186:         USE GJK_UTILS, ONLY: DIM 
187:  
188:         IMPLICIT NONE 
189:  
190:         DOUBLE PRECISION, INTENT(INOUT) :: X(X_SIZE) 
191:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(X_SIZE), ENERGY 
192:         LOGICAL, INTENT(IN) :: GTEST 
193:  
194:         ! MOLI_IND, MOLJ_IND: indices for looping over molecules 
195:         ! OFFSET: The start of the orientation coordinates in X 
196:         INTEGER :: MOLI_IND, MOLJ_IND, OFFSET 
197:  
198:         ! INIT_AXIS: Initial direction guess for GJK, just use the difference of 
199:         !            the centres 
200:         ! RESULT_SIMPLEX: a simplex enclosing the origin if GJK finds an 
201:         !                 intersection 
202:         DOUBLE PRECISION :: INIT_AXIS(DIM), RESULT_SIMPLEX(DIM, DIM+1) 
203:  
204:         ! RESULT: whether or not GJK finds an intersection 
205:         LOGICAL :: RESULT 
206:  
207:         ! Initialise ouput variables 
208:         ENERGY = 0.D0 
209:         IF (GTEST) GRAD(:) = 0.D0 
210:  
211:         ! Calculate offset from the number of particles 
212:         OFFSET = 3*NPOLY + 1 
213:  
214:         ! Update all the particle positions and rotation matrices 
215:         DO MOLI_IND = 1, NPOLY 
216:             CALL UPDATE_POLYHEDRON(POLYHEDRA(MOLI_IND), & 
217:                                    X(MOLI_IND*3-2:MOLI_IND*3), & 
218:                                    X(OFFSET+MOLI_IND*3-2:OFFSET+MOLI_IND*3), & 
219:                                    GTEST) 
220:         END DO 
221:  
222:         ! Just looking to test at the moment 
223:         ! Testing with two particles, just print whether they overlap or not 
224:         INIT_AXIS(:) = POLYHEDRA(1)%R(:) - POLYHEDRA(2)%R(:) 
225:         CALL GJK_INTERSECTION(POLYHEDRA(1), POLYHEDRA(2), INIT_AXIS(:), DIM, & 
226:                               RESULT, RESULT_SIMPLEX(:,:)) 
227:  
228:         WRITE(*,*) 'GJK returned ', RESULT 
229:         IF (RESULT) THEN 
230:             WRITE(*,*) 'Result simplex is:' 
231:             DO MOLJ_IND = 1, 4 
232:                 WRITE(*,*) RESULT_SIMPLEX(1,MOLJ_IND), & 
233:                            RESULT_SIMPLEX(2,MOLJ_IND), & 
234:                            RESULT_SIMPLEX(3,MOLJ_IND) 
235:             END DO 
236:         END IF 
237:         WRITE (*,*) 'Exiting now' 
238:         STOP 
239:  
240:     END SUBROUTINE CONVEX_POLYHEDRA 
241:  
242: !------------------------------------------------------------------------------ 
243: ! Prints out the vertex information to the file poly.xyz 
244: !------------------------------------------------------------------------------ 
245:     SUBROUTINE VIEW_POLYHEDRA() 
246:  
247:         IMPLICIT NONE 
248:  
249:         ! Doesn't actually do anything yet 
250:  
251:     END SUBROUTINE VIEW_POLYHEDRA 
252:  
253: END MODULE CONVEX_POLYHEDRA_MODULE 


r31861/finalio.f90 2017-01-30 19:30:25.598345496 +0000 r31860/finalio.f90 2017-01-30 19:30:27.606372662 +0000
 23:   USE MODAMBER9, ONLY : COORDS1,IH,M04,AMBFINALIO_NODE 23:   USE MODAMBER9, ONLY : COORDS1,IH,M04,AMBFINALIO_NODE
 24:   USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, & 24:   USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, &
 25:        AMBER12_WRITE_XYZ 25:        AMBER12_WRITE_XYZ
 26:   USE OPEP_INTERFACE_MOD, ONLY: OPEP_FINISH, OPEP_WRITE_PDB 26:   USE OPEP_INTERFACE_MOD, ONLY: OPEP_FINISH, OPEP_WRITE_PDB
 27:   USE QMODULE 27:   USE QMODULE
 28:   USE MODCHARMM 28:   USE MODCHARMM
 29:   USE AMHGLOBALS, ONLY:NMRES,IRES 29:   USE AMHGLOBALS, ONLY:NMRES,IRES
 30:   USE BGUPMOD 30:   USE BGUPMOD
 31:   USE PERMU 31:   USE PERMU
 32:   USE MODHESS, ONLY : HESS, MASSWT 32:   USE MODHESS, ONLY : HESS, MASSWT
 33:   USE CONVEX_POLYHEDRA_MODULE, ONLY: VIEW_POLYHEDRA 
 34:    33:   
 35:   IMPLICIT NONE 34:   IMPLICIT NONE
 36:    35:   
 37:   !   MCP 36:   !   MCP
 38:   INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, AMHUNIT2, LUNIT, NSYMOPS 37:   INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, AMHUNIT2, LUNIT, NSYMOPS
 39:   INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NC, NRBS1, NRBS2, MYUNIT4 38:   INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NC, NRBS1, NRBS2, MYUNIT4
 40:   DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY, ITDET 39:   DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY, ITDET
 41:   ! 40:   !
 42:   !ds656> Symmetry detection for clusters on a substrate... 41:   !ds656> Symmetry detection for clusters on a substrate...
 43:   INTEGER :: ISITES(MIEF_NSITES), NSITES, COMBO_N 42:   INTEGER :: ISITES(MIEF_NSITES), NSITES, COMBO_N
739:            DCOORDS(3*(J2-1)+2)=QMINP(J1,3*(J2-1)+2)738:            DCOORDS(3*(J2-1)+2)=QMINP(J1,3*(J2-1)+2)
740:            DCOORDS(3*(J2-1)+3)=QMINP(J1,3*(J2-1)+3)739:            DCOORDS(3*(J2-1)+3)=QMINP(J1,3*(J2-1)+3)
741:         ENDDO740:         ENDDO
742:         CALL CHARMMDUMP(DCOORDS,DBNAME(J1))741:         CALL CHARMMDUMP(DCOORDS,DBNAME(J1))
743:         742:         
744:         !    DC430 >743:         !    DC430 >
745:         !    |gd351> added patchy744:         !    |gd351> added patchy
746:         745:         
747:      ELSE IF (CAPBINT.OR.DBPT.OR.DBPTDT.OR.DMBLMT.OR.DMBLPYT.OR.LINRODT.OR.LWOTPT.OR.MSTBINT.OR.MSSTOCKT.OR.NCAPT.OR.NPAHT &746:      ELSE IF (CAPBINT.OR.DBPT.OR.DBPTDT.OR.DMBLMT.OR.DMBLPYT.OR.LINRODT.OR.LWOTPT.OR.MSTBINT.OR.MSSTOCKT.OR.NCAPT.OR.NPAHT &
748:           .OR. NTIPT .OR. STOCKAAT .OR. PAHAT .OR. PAHW99T .OR. TDHDT .OR. WATERDCT .OR. WATERKZT .OR. PATCHY .OR. PAPT   &747:           .OR. NTIPT .OR. STOCKAAT .OR. PAHAT .OR. PAHW99T .OR. TDHDT .OR. WATERDCT .OR. WATERKZT .OR. PATCHY .OR. PAPT   &
749:           .OR. PAPBINT .OR. PAPJANT .OR. POLYT .OR. PTSTSTT .OR. MORSEDPT) THEN748:           .OR. PAPBINT .OR. PAPJANT .OR. PTSTSTT .OR. MORSEDPT) THEN
750:         DO J2 = 1, NATOMS/2749:         DO J2 = 1, NATOMS/2
751:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)750:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)
752:         ENDDO751:         ENDDO
753:         DO J2 = 1, NATOMS/2752:         DO J2 = 1, NATOMS/2
754:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)753:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)
755:         ENDDO754:         ENDDO
756:         755:         
757:      ELSE IF (GBT.OR.GBDT.OR.GBDPT.OR.MSGBT) THEN756:      ELSE IF (GBT.OR.GBDT.OR.GBDPT.OR.MSGBT) THEN
758:         DO J2 = 1, NATOMS/2757:         DO J2 = 1, NATOMS/2
759:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)758:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)
1547:      1546:      
1548:   ELSE IF (PAPBINT) THEN1547:   ELSE IF (PAPBINT) THEN
1549:      1548:      
1550:      CALL VIEWPAPBIN()1549:      CALL VIEWPAPBIN()
1551:      RETURN1550:      RETURN
1552:      1551:      
1553:   ELSE IF (PAPJANT) THEN1552:   ELSE IF (PAPJANT) THEN
1554:      1553:      
1555:      CALL VIEWPAPJANUS()1554:      CALL VIEWPAPJANUS()
1556:      RETURN1555:      RETURN
1557:  
1558:   ELSE IF (POLYT) THEN 
1559:  
1560:      CALL VIEW_POLYHEDRA() 
1561:      RETURN 
1562:      1556:      
1563:   ELSE IF (PTSTSTT) THEN1557:   ELSE IF (PTSTSTT) THEN
1564:      1558:      
1565:      CALL VIEWPTSTST()1559:      CALL VIEWPTSTST()
1566:      RETURN1560:      RETURN
1567:      1561:      
1568:      !|gd351>1562:      !|gd351>
1569:      1563:      
1570:   ELSE IF (PATCHY) THEN1564:   ELSE IF (PATCHY) THEN
1571:      1565:      


r31861/gjk.f90 2017-01-30 19:30:25.910349714 +0000 r31860/gjk.f90 2017-01-30 19:30:27.994377913 +0000
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/gjk.f90' in revision 31860
  2: !   Copyright (C) 1999-2017 David J. Wales 
  3: !   This file is part of GMIN. 
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: ! An implementation of the Gilbert–Johnson–Keerthi algorithm to test whether 
 20: ! convex shapes overlap and the expanding polytope algorithm to find the minimum 
 21: ! distance move to get rid of overlap. See 
 22: ! http://www.dyn4j.org/2010/04/gjk-gilbert-johnson-keerthi/#gjk-support 
 23: ! (particularly the video tutorial) for an explanation of GJK and 
 24: ! http://www.dyn4j.org/2010/05/epa-expanding-polytope-algorithm/ for EPA. 
 25: ! Questions to jwrm2 
 26:  
 27: MODULE GJK_MODULE 
 28: CONTAINS 
 29: !------------------------------------------------------------------------------- 
 30: ! Performs the GJK intersection test 
 31: ! SHAPE1: The first shape of the test 
 32: ! SHAPE2: The second shape of the test 
 33: ! INIT_AXIS: The initial guess, just needs to be a point on the surface of the 
 34: !            Minkowski difference of the shapes 
 35: ! DIM: dimension of the space (probably 3) 
 36: ! RESULT: whether or not the shapes intersect 
 37: ! RESULT_SIMPLEX: if the shapes intersect, the simplex containing the origin 
 38: !------------------------------------------------------------------------------- 
 39:     SUBROUTINE GJK_INTERSECTION(SHAPE1, SHAPE2, INIT_AXIS, DIM, RESULT, & 
 40:                                 RESULT_SIMPLEX) 
 41:  
 42:         ! MYUNIT: File unit of GMIN_out 
 43:         ! SUPPORT_FUNC: Support function for convex polyhedra 
 44:         ! POLYHEDRON: Type to store a polyhedron 
 45:         USE COMMONS, ONLY: MYUNIT 
 46:         USE GJK_UTILS, ONLY: POLYHEDRON, SUPPORT_FUNC 
 47:  
 48:         IMPLICIT NONE 
 49:  
 50:         TYPE(POLYHEDRON), INTENT(IN) :: SHAPE1, SHAPE2 
 51:         INTEGER, INTENT(IN) :: DIM 
 52:         DOUBLE PRECISION, INTENT(IN) :: INIT_AXIS(DIM) 
 53:         DOUBLE PRECISION, INTENT(OUT) :: RESULT_SIMPLEX(DIM, DIM+1) 
 54:         LOGICAL, INTENT(OUT) :: RESULT 
 55:  
 56:         ! WORK_SIMPLEX: the working simplex, unused parts set to 1.D100 
 57:         ! SEARCH_DIRECTION: the current search direction 
 58:         ! NEW_POINT: the next point to be considered for addition to the simplex 
 59:         ! LEN_SIMPLEX: how many vertices the simplex currently contains 
 60:         DOUBLE PRECISION :: WORK_SIMPLEX(DIM, DIM+1), SEARCH_DIRECTION(DIM) 
 61:         DOUBLE PRECISION :: NEW_POINT(DIM) 
 62:         INTEGER :: LEN_SIMPLEX 
 63:  
 64:         ! Initialise the output variables 
 65:         RESULT = .FALSE. 
 66:         RESULT_SIMPLEX(:,:) = 1.D100 
 67:  
 68:         ! Initialise the working simplex 
 69:         WORK_SIMPLEX(:,:) = 1.D100 
 70:         WORK_SIMPLEX(:, 1) = SUPPORT_FUNC(SHAPE1,  INIT_AXIS(:)) & 
 71:                            - SUPPORT_FUNC(SHAPE2, -INIT_AXIS(:)) 
 72:         SEARCH_DIRECTION(:) = - WORK_SIMPLEX(:, 1) 
 73:         LEN_SIMPLEX = 1 
 74:  
 75:         ! Start the GJK refinement loop 
 76:         DO 
 77:             ! Go as far as we can within the Minkowski difference in the search 
 78:             ! direction. 
 79:             NEW_POINT(:) = SUPPORT_FUNC(SHAPE1,  SEARCH_DIRECTION(:)) & 
 80:                          - SUPPORT_FUNC(SHAPE2, -SEARCH_DIRECTION(:)) 
 81:  
 82:             ! If we haven't gone as far as the origin, the shapes definitely 
 83:             ! don't intersect 
 84:             IF (DOT_PRODUCT(SEARCH_DIRECTION(:), NEW_POINT(:)) < 0.D0) THEN 
 85:                 RESULT = .FALSE. 
 86:                 EXIT 
 87:             END IF 
 88:  
 89:             ! Add the new point to the simplex 
 90:             LEN_SIMPLEX = LEN_SIMPLEX + 1 
 91:             WORK_SIMPLEX(:, LEN_SIMPLEX) = NEW_POINT(:) 
 92:  
 93:             ! Check if the new simplex contains the origin 
 94:             ! Update the simplex and the search direction 
 95:             CALL NEAREST_SIMPLEX(WORK_SIMPLEX, LEN_SIMPLEX, DIM, & 
 96:                                  SEARCH_DIRECTION, RESULT) 
 97:  
 98:             ! If the simplex contains the origin, we're done 
 99:             IF (RESULT) THEN 
100:                 ! The simplex length should be DIM+1 (ie a tetrahedron in 3D), 
101:                 ! or something horrible has happened. 
102:                 IF (LEN_SIMPLEX .NE. DIM + 1) THEN 
103:                     WRITE(MYUNIT, *) 'GJK_INTERSECTION> ERROR, length of ', & 
104:                     ' simpled is ', LEN_SIMPLEX, ' but dimension is ', DIM 
105:                     STOP 
106:                 END IF 
107:  
108:                 RESULT_SIMPLEX(:,:) = WORK_SIMPLEX(:,:) 
109:                 EXIT 
110:             END IF 
111:  
112:         END DO ! GJK refinement loop 
113:  
114:     END SUBROUTINE GJK_INTERSECTION 
115:  
116: !------------------------------------------------------------------------------- 
117: ! Finds the nearest simplex to the origin of rank lower than the supplied 
118: ! simplex. If the simplex is at the maximum size, check if the origin lies 
119: ! within the simplex. 
120: ! Provides the new search direction based on the closest simplex. 
121: ! TEST_SIMPLEX: the simplex to check 
122: ! LEN_SIMPLEX: current length of the simplex 
123: ! DIM: dimension of the space (probably 3) 
124: ! SEARCH_DIRECTION: the new search direction 
125: ! RESULT: whether the origin lies within the simplex 
126: !------------------------------------------------------------------------------- 
127:     SUBROUTINE NEAREST_SIMPLEX(TEST_SIMPLEX, LEN_SIMPLEX, DIM, & 
128:                                SEARCH_DIRECTION, RESULT) 
129:  
130:         ! MYUNIT: File unit of GMIN_out 
131:         ! VEC_CROSS: vector cross product 
132:         USE COMMONS, ONLY: MYUNIT 
133:         USE VEC3, ONLY: VEC_CROSS 
134:  
135:         IMPLICIT NONE 
136:  
137:         INTEGER, INTENT(IN) :: DIM 
138:         DOUBLE PRECISION, INTENT(INOUT) :: TEST_SIMPLEX(DIM, DIM+1) 
139:         DOUBLE PRECISION, INTENT(OUT) :: SEARCH_DIRECTION(DIM) 
140:         INTEGER, INTENT(INOUT) :: LEN_SIMPLEX 
141:         LOGICAL, INTENT(OUT) :: RESULT 
142:  
143:         ! O: the origin 
144:         ! AB, AC, AD, AO, ABPERP, ACPERP, ADPERP, ABCPERP, ACDPERP, ADBPERP: 
145:         !     storage for plane tests and search direction calculation 
146:         DOUBLE PRECISION :: AB(DIM), AC(DIM), AD(DIM), AO(DIM) 
147:         DOUBLE PRECISION :: ABPERP(DIM), ACPERP(DIM), ADPERP(DIM) 
148:         DOUBLE PRECISION :: ABCPERP(DIM), ACDPERP(DIM), ADBPERP(DIM) 
149:         ! OVERABC, OVERACD, OVERADB: Results of plane tests 
150:         LOGICAL :: OVERABC, OVERACD, OVERADB 
151:  
152:         ! Initialise output 
153:         SEARCH_DIRECTION(:) = 1.D100 
154:         RESULT = .FALSE. 
155:  
156:         ! The point referred to as A is always the mose recent added to the 
157:         ! simplex, B the second most recent, etc. 
158:         ! No point can be the closest to the origin. In the GJK loop we've 
159:         ! checked to make sure the origin has been passed. Since it has, there's 
160:         ! no way the points can be closer than the lines connecting them to 
161:         ! other points. 
162:         ! Likewise, only lines and faces involving A are candidates. 
163:  
164:         ! Now assuming that DIM is 3. May need to rewrite. 
165:         IF (LEN_SIMPLEX .EQ. 2) THEN 
166:             ! Simplest case of a line. The nearest simplex can only be the line. 
167:             RESULT = .FALSE. 
168:             AO(:) = - TEST_SIMPLEX(:, 2) 
169:             AB(:) = TEST_SIMPLEX(:, 1) - TEST_SIMPLEX(:, 2) 
170:             SEARCH_DIRECTION(:) = CROSS_121(AB(:), AO(:), DIM) 
171:  
172:         ELSE IF (LEN_SIMPLEX .EQ. 3) THEN 
173:             ! Simplex is a triangle. Nearest simplex still cannot be a point. 
174:             ! It could be either of the lines with the latest point, but not the 
175:             ! third. 
176:             ! Origin could lie above or below the triangle. 
177:             RESULT = .FALSE. 
178:  
179:             ! Generate the necessary vectors 
180:             AO(:) = - TEST_SIMPLEX(:, 3) 
181:             AB(:) = TEST_SIMPLEX(:, 2) - TEST_SIMPLEX(:, 3) 
182:             AC(:) = TEST_SIMPLEX(:, 1) - TEST_SIMPLEX(:, 3) 
183:             ABCPERP(:) = VEC_CROSS(AB(:), AC(:)) 
184:             ABPERP(:) = VEC_CROSS(AB(:), ABCPERP(:)) 
185:             ACPERP(:) = VEC_CROSS(ABCPERP(:), AC(:)) 
186:  
187:             ! Check whether the origin is outside the triangle on the AB side 
188:             IF (SAME_DIRECTION(ABPERP(:), AO(:), DIM, MYUNIT)) THEN 
189:                 ! AB is the closest simplex 
190:                 ! Delete C from the simplex 
191:                 TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 2) 
192:                 TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 3) 
193:                 LEN_SIMPLEX = 2 
194:                 SEARCH_DIRECTION(:) = CROSS_121(AB(:), AO(:), DIM) 
195:  
196:             ! Check whether the origin is outside the triangle on the AC side 
197:             ELSE IF (SAME_DIRECTION(ACPERP(:), AO(:), DIM, MYUNIT)) THEN 
198:                 ! AC is the closest simplex 
199:                 ! Delete B from the simplex 
200:                 TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 3) 
201:                 LEN_SIMPLEX = 2 
202:                 SEARCH_DIRECTION(:) = CROSS_121(AC(:), AO(:), DIM) 
203:             ELSE 
204:                 ! ABC is the closest simplex 
205:                 ! Need to pick the correct direction (above or below face) 
206:                 IF (SAME_DIRECTION(ABCPERP(:), AO(:), DIM, MYUNIT)) THEN 
207:                     SEARCH_DIRECTION(:) = ABCPERP(:) 
208:                 ELSE 
209:                     SEARCH_DIRECTION(:) = -ABCPERP(:) 
210:                     ! Swap the simplex ordering, so the point is above. Makes it 
211:                     ! easier for tetrahedra if this is guaranteed. 
212:                     TEST_SIMPLEX(:, 4) = TEST_SIMPLEX(:, 1) 
213:                     TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 2) 
214:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
215:                 END IF 
216:             END IF 
217:  
218:         ELSE IF (LEN_SIMPLEX .EQ. 4) THEN 
219:             ! Simplex is a tetrahedron. It might enclose the origin. 
220:             ! Only lines and faces involving A need to be considered. 
221:             ! See http://vec3.ca/gjk/implementation/ for an explanation 
222:  
223:             ! Generate the necessary vectors 
224:             AO(:) = - TEST_SIMPLEX(:, 4) 
225:             ! Edge vectors 
226:             AB(:) = TEST_SIMPLEX(:, 3) - TEST_SIMPLEX(:, 4) 
227:             AC(:) = TEST_SIMPLEX(:, 2) - TEST_SIMPLEX(:, 4) 
228:             AD(:) = TEST_SIMPLEX(:, 1) - TEST_SIMPLEX(:, 4) 
229:             ! Perpendiculars to faces (pointing outwards) 
230:             ABCPERP(:) = VEC_CROSS(AB(:), AC(:)) 
231:             ACDPERP(:) = VEC_CROSS(AC(:), AD(:)) 
232:             ADBPERP(:) = VEC_CROSS(AD(:), AB(:)) 
233:  
234:             ! Store the results of plane tests, so we only have to do them once 
235:             OVERABC = SAME_DIRECTION(ABCPERP(:), AO(:), DIM, MYUNIT) 
236:             OVERACD = SAME_DIRECTION(ACDPERP(:), AO(:), DIM, MYUNIT) 
237:             OVERADB = SAME_DIRECTION(ADBPERP(:), AO(:), DIM, MYUNIT) 
238:  
239:             IF (.NOT. (OVERABC .OR. OVERACD .OR. OVERADB)) THEN 
240:                 ! Origin is inside the simplex 
241:                 RESULT = .TRUE. 
242:  
243:             ! Deal with the two faces case 
244:             ELSE IF (OVERABC .AND. OVERACD .AND. .NOT. OVERADB) THEN 
245:                 ! Check the common edge, relative to the first face 
246:                 ACPERP(:) = VEC_CROSS(ABCPERP(:), AC(:)) 
247:                 IF (SAME_DIRECTION(ACPERP(:), AO(:), DIM, MYUNIT)) THEN 
248:                     ! Now need to do the full check for the second face. 
249:                     ! Fiddle the logicals and let it fall through. 
250:                     OVERABC = .FALSE. 
251:                 ELSE 
252:                     ! Just need to decide between the first face and the other 
253:                     ! edge. However, it's simpler (if slightly less efficient) 
254:                     ! to do the full check for the first face. 
255:                     OVERACD = .FALSE. 
256:                 ENDIF 
257:             ELSE IF (OVERACD .AND. OVERADB .AND. .NOT. OVERABC) THEN 
258:                 ! Check the common edge, relative to the first face 
259:                 ADPERP(:) = VEC_CROSS(ACDPERP(:), AD(:)) 
260:                 IF (SAME_DIRECTION(ADPERP(:), AO(:), DIM, MYUNIT)) THEN 
261:                     ! Now need to do the full check for the second face. 
262:                     ! Fiddle the logicals and let it fall through. 
263:                     OVERACD = .FALSE. 
264:                 ELSE 
265:                     ! Just need to decide between the first face and the other 
266:                     ! edge. However, it's simpler (if slightly less efficient) 
267:                     ! to do the full check for the first face. 
268:                     OVERADB = .FALSE. 
269:                 ENDIF 
270:             ELSE IF (OVERADB .AND. OVERABC .AND. .NOT. OVERACD) THEN 
271:                 ! Check the common edge, relative to the first face 
272:                 ABPERP(:) = VEC_CROSS(ADBPERP(:), AB(:)) 
273:                 IF (SAME_DIRECTION(ABPERP(:), AO(:), DIM, MYUNIT)) THEN 
274:                     ! Now need to do the full check for the second face. 
275:                     ! Fiddle the logicals and let it fall through. 
276:                     OVERADB = .FALSE. 
277:                 ELSE 
278:                     ! Just need to decide between the first face and the other 
279:                     ! edge. However, it's simpler (if slightly less efficient) 
280:                     ! to do the full check for the first face. 
281:                     OVERABC = .FALSE. 
282:                 ENDIF 
283:             ENDIF ! Two faces case 
284:  
285:             ! Deal with the only one face cases. 
286:             ! Need to check the bordering lines. 
287:             IF (OVERABC .AND. .NOT. OVERACD .AND. .NOT. OVERADB) THEN 
288:                 ! Perpendiculars to the edges, in plane of this face 
289:                 ABPERP(:) = VEC_CROSS(AB(:), ABCPERP(:)) 
290:                 ACPERP(:) = VEC_CROSS(ABCPERP(:), AC(:)) 
291:                 IF (SAME_DIRECTION(ABPERP(:), AO(:), DIM, MYUNIT)) THEN 
292:                     ! In the AB line region 
293:                     TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 3) 
294:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
295:                     LEN_SIMPLEX = 2 
296:                     SEARCH_DIRECTION(:) = CROSS_121(AB(:), AO(:), DIM) 
297:                 ELSE IF (SAME_DIRECTION(ACPERP(:), AO(:), DIM, MYUNIT)) THEN 
298:                     ! In the AC line region 
299:                     TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 2) 
300:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
301:                     LEN_SIMPLEX = 2 
302:                     SEARCH_DIRECTION(:) = CROSS_121(AC(:), AO(:), DIM) 
303:                 ELSE 
304:                     ! In the ABC region 
305:                     TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 2) 
306:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 3) 
307:                     TEST_SIMPLEX(:, 3) = TEST_SIMPLEX(:, 4) 
308:                     LEN_SIMPLEX = 3 
309:                     SEARCH_DIRECTION(:) = ABCPERP(:) 
310:                 END IF 
311:             ELSE IF (OVERACD .AND. .NOT. OVERADB .AND. .NOT. OVERABC) THEN 
312:                 ! Perpendiculars to the edges, in plane of this face 
313:                 ACPERP(:) = VEC_CROSS(AC(:), ACDPERP(:)) 
314:                 ADPERP(:) = VEC_CROSS(ACDPERP(:), AD(:)) 
315:                 IF (SAME_DIRECTION(ACPERP(:), AO(:), DIM, MYUNIT)) THEN 
316:                     ! In the AC line region 
317:                     TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 2) 
318:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
319:                     LEN_SIMPLEX = 2 
320:                     SEARCH_DIRECTION(:) = CROSS_121(AC(:), AO(:), DIM) 
321:                 ELSE IF (SAME_DIRECTION(ADPERP(:), AO(:), DIM, MYUNIT)) THEN 
322:                     ! In the AD line region 
323:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
324:                     LEN_SIMPLEX = 2 
325:                     SEARCH_DIRECTION(:) = CROSS_121(AD(:), AO(:), DIM) 
326:                 ELSE 
327:                     ! In the ACD region 
328:                     TEST_SIMPLEX(:, 3) = TEST_SIMPLEX(:, 4) 
329:                     LEN_SIMPLEX = 3 
330:                     SEARCH_DIRECTION(:) = ACDPERP(:) 
331:                 END IF 
332:             ELSE IF (OVERADB .AND. .NOT. OVERABC .AND. .NOT. OVERACD) THEN 
333:                 ! Perpendiculars to the edges, in plane of this face 
334:                 ADPERP(:) = VEC_CROSS(AD(:), ADBPERP(:)) 
335:                 ABPERP(:) = VEC_CROSS(ADBPERP(:), AB(:)) 
336:                 IF (SAME_DIRECTION(ADPERP(:), AO(:), DIM, MYUNIT)) THEN 
337:                     ! In the AD line region 
338:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
339:                     LEN_SIMPLEX = 2 
340:                     SEARCH_DIRECTION(:) = CROSS_121(AD(:), AO(:), DIM) 
341:                 ELSE IF (SAME_DIRECTION(ABPERP(:), AO(:), DIM, MYUNIT)) THEN 
342:                     ! In the AB line region 
343:                     TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 3) 
344:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
345:                     LEN_SIMPLEX = 2 
346:                     SEARCH_DIRECTION(:) = CROSS_121(AB(:), AO(:), DIM) 
347:                 ELSE 
348:                     ! In the ADC region 
349:                     TEST_SIMPLEX(:, 3) = TEST_SIMPLEX(:, 4) 
350:                     TEST_SIMPLEX(:, 4) = TEST_SIMPLEX(:, 1) 
351:                     TEST_SIMPLEX(:, 1) = TEST_SIMPLEX(:, 2) 
352:                     TEST_SIMPLEX(:, 2) = TEST_SIMPLEX(:, 4) 
353:                     LEN_SIMPLEX = 3 
354:                     SEARCH_DIRECTION(:) = ABCPERP(:) 
355:                 END IF 
356:             END IF ! One face case 
357:         ELSE 
358:             ! Simplex is not the right size, something horrible has happened 
359:             WRITE (MYUNIT, *) 'NEAREST_SIMPLEX> ERROR, simplex size is ', & 
360:                 LEN_SIMPLEX, ' This should not happen.' 
361:             STOP 
362:         END IF ! Simplex size 
363:  
364:     END SUBROUTINE NEAREST_SIMPLEX 
365:  
366: !------------------------------------------------------------------------------- 
367: ! Utility function for dot products. 
368: ! All we care about is whether or not the dot is positive or negative 
369: ! Print a warning if it's zero 
370: ! VEC1, VEC2: the vectors to dot 
371: ! DIM: dimension of vectors 
372: ! OUT_UNIT: File unit for warnings 
373: !------------------------------------------------------------------------------- 
374:     LOGICAL FUNCTION SAME_DIRECTION(VEC1, VEC2, DIM, OUT_UNIT) 
375:  
376:         IMPLICIT NONE 
377:  
378:         INTEGER, INTENT(IN) :: DIM, OUT_UNIT 
379:         DOUBLE PRECISION, INTENT(IN) :: VEC1(DIM), VEC2(DIM) 
380:  
381:         ! DOT: store the dot product 
382:         DOUBLE PRECISION :: DOT 
383:  
384:         DOT = DOT_PRODUCT(VEC1(:), VEC2(:)) 
385:  
386:         IF (DOT .EQ. 0.D0) THEN 
387:             WRITE(OUT_UNIT, *) 'NEAREST_SIMPLEX> WARNING: exact overlap' 
388:             SAME_DIRECTION = .FALSE. 
389:         ELSE 
390:             SAME_DIRECTION = (DOT .GT. 0.D0) 
391:         END IF 
392:     END FUNCTION SAME_DIRECTION 
393:  
394: !------------------------------------------------------------------------------- 
395: ! Utility function for cross products. 
396: ! Carries out VEC1 x VEC2 x VEC1 
397: ! VEC1, VEC2: the vectors to cross 
398: ! DIM: Dimension of the space (which I suppose must be 3 for cross products) 
399: !------------------------------------------------------------------------------- 
400:     FUNCTION CROSS_121(VEC1, VEC2, DIM) 
401:  
402:         ! VEC_CROSS: cross producr 
403:         USE VEC3, ONLY: VEC_CROSS 
404:  
405:         IMPLICIT NONE 
406:  
407:         INTEGER, INTENT(IN) :: DIM 
408:         DOUBLE PRECISION, INTENT(IN) :: VEC1(DIM), VEC2(DIM) 
409:         DOUBLE PRECISION :: CROSS_121(DIM) 
410:  
411:         CROSS_121(:) = VEC_CROSS(VEC_CROSS(VEC1(:), VEC2(:)), VEC1(:)) 
412:     END FUNCTION CROSS_121 
413:  
414: END MODULE GJK_MODULE 


r31861/gjk_utils.f90 2017-01-30 19:30:26.178353341 +0000 r31860/gjk_utils.f90 2017-01-30 19:30:28.350382731 +0000
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/gjk_utils.f90' in revision 31860
  2: !   Copyright (C) 1999-2017 David J. Wales 
  3: !   This file is part of GMIN. 
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: ! Utilities for the implementation of the GJK and EPA algorithms. 
 20: ! Put here to avoid circular dependencies 
 21:  
 22: MODULE GJK_UTILS 
 23:  
 24:     ! DIM: Working in 3 dimensions 
 25:     INTEGER, PARAMETER :: DIM = 3 
 26:  
 27:     ! NVERTS: number of vertices 
 28:     INTEGER :: NVERTS 
 29:  
 30:     ! Stores a single polyhedron 
 31:     TYPE POLYHEDRON 
 32:         ! R: Coordinates of the centre of the particle 
 33:         ! P: Angle-axis style orientation of the particle 
 34:         ! RM: Rotation matrix, derived from P 
 35:         ! RMD: Derivatives of the rotation matrix 
 36:         DOUBLE PRECISION :: R(DIM), P(DIM), RM(DIM, DIM), RMD(DIM, DIM, DIM) 
 37:         ! VERT: Vertex positions 
 38:         DOUBLE PRECISION, ALLOCATABLE :: VERTS(:,:) 
 39:     END TYPE POLYHEDRON 
 40:  
 41: CONTAINS 
 42: !------------------------------------------------------------------------------- 
 43: ! Support function for the GJK algorithm. Returns the vertex furthest along 
 44: ! the direction DIR 
 45: ! POLY: the polyhedron of interest 
 46: ! DIR: the search direction 
 47: !------------------------------------------------------------------------------- 
 48:     FUNCTION SUPPORT_FUNC(POLY, DIR) 
 49:  
 50:         IMPLICIT NONE 
 51:  
 52:         DOUBLE PRECISION :: SUPPORT_FUNC(DIM) 
 53:         DOUBLE PRECISION :: DIR(:) 
 54:         TYPE(POLYHEDRON), INTENT(IN) :: POLY 
 55:  
 56:         ! BEST_INDEX: index of the best vertex 
 57:         INTEGER :: I, BEST_INDEX 
 58:         ! DOT: temporary dot product 
 59:         ! MAX_DOT: largest dot product between DIR and a vertex 
 60:         DOUBLE PRECISION :: DOT, MAX_DOT 
 61:  
 62:         BEST_INDEX = 1 
 63:         MAX_DOT = -1.0D-100 
 64:  
 65:         ! Loop over the vertices 
 66:         DO I = 1, NVERTS 
 67:             DOT = DOT_PRODUCT(POLY%VERTS(:,I), DIR(:)) 
 68:             IF (DOT .GT. MAX_DOT) THEN 
 69:                 MAX_DOT = DOT 
 70:                 BEST_INDEX = I 
 71:             END IF 
 72:         END DO ! Loop over vertices 
 73:  
 74:         SUPPORT_FUNC(:) = POLY%VERTS(:,BEST_INDEX) 
 75:     END FUNCTION SUPPORT_FUNC 
 76:  
 77: END MODULE GJK_UTILS 


r31861/keywords.f 2017-01-30 19:30:26.454357075 +0000 r31860/keywords.f 2017-01-30 19:30:29.310395722 +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 
 43:       USE MBPOLMOD, ONLY: MBPOLINIT 42:       USE MBPOLMOD, ONLY: MBPOLINIT
 44:       USE SWMOD, ONLY: SWINIT, MWINIT 43:       USE SWMOD, ONLY: SWINIT, MWINIT
 45:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS, 44:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS,
 46:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA 45:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA
 47:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 46:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 48:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 47:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 49:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,  48:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 
 50:      &     PARSE_MLJ_PARAMS 49:      &     PARSE_MLJ_PARAMS
 51:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT 50:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT
 52:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE 51:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE
846:       MSGBT       = .FALSE. 845:       MSGBT       = .FALSE. 
847:       MSTBINT     = .FALSE.846:       MSTBINT     = .FALSE.
848:       MSSTOCKT    = .FALSE.847:       MSSTOCKT    = .FALSE.
849:       MULTPAHAT   = .FALSE.848:       MULTPAHAT   = .FALSE.
850:       NCAPT       = .FALSE.849:       NCAPT       = .FALSE.
851:       NPAHT       = .FALSE.850:       NPAHT       = .FALSE.
852:       NTIPT       = .FALSE.851:       NTIPT       = .FALSE.
853:       NZERO=0                   ! jdf43>852:       NZERO=0                   ! jdf43>
854:       PAHAT       = .FALSE.853:       PAHAT       = .FALSE.
855:       PAPT        = .FALSE.854:       PAPT        = .FALSE.
856:       POLYT       = .FALSE. 
857:       PTSTSTT     = .FALSE.855:       PTSTSTT     = .FALSE.
858:       RBSYMT      = .FALSE.     ! jdf43>856:       RBSYMT      = .FALSE.     ! jdf43>
859:       SHIFTV=1.0D6              ! jdf43>857:       SHIFTV=1.0D6              ! jdf43>
860:       SANDBOXT    = .FALSE.858:       SANDBOXT    = .FALSE.
861:       SILANET     = .FALSE.859:       SILANET     = .FALSE.
862:       STOCKAAT    = .FALSE.860:       STOCKAAT    = .FALSE.
863:       TDHDT       = .FALSE.861:       TDHDT       = .FALSE.
864:       WATERDCT    = .FALSE.862:       WATERDCT    = .FALSE.
865:       WATERKZT    = .FALSE.863:       WATERKZT    = .FALSE.
866: !|gd351>864: !|gd351>
6590:          CALL READF(YKAPPA)6588:          CALL READF(YKAPPA)
6591: 6589: 
6592:          NRBSITES = 26590:          NRBSITES = 2
6593:          ALLOCATE(RBSTLA(NRBSITES,3))6591:          ALLOCATE(RBSTLA(NRBSITES,3))
6594: 6592: 
6595:          CALL DEFPAPJANUS()6593:          CALL DEFPAPJANUS()
6596: 6594: 
6597:          PAPJANT = .TRUE.6595:          PAPJANT = .TRUE.
6598:          RIGID   = .TRUE.6596:          RIGID   = .TRUE.
6599: 6597: 
6600:       ELSE IF (WORD .EQ. 'POLYHEDRA') THEN 
6601:  
6602:         POLYT = .TRUE. 
6603:         RIGID = .TRUE. 
6604:         CALL INITIALISE_POLYHEDRA() 
6605:  
6606:       ELSE IF (WORD .EQ. 'PTSTST') THEN6598:       ELSE IF (WORD .EQ. 'PTSTST') THEN
6607: 6599: 
6608:          CALL READI(NPATCH)6600:          CALL READI(NPATCH)
6609:          CALL READF(PAPEPS)6601:          CALL READF(PAPEPS)
6610:          CALL READF(PAPCD)6602:          CALL READF(PAPCD)
6611:          CALL READF(YKAPPA)6603:          CALL READF(YKAPPA)
6612: 6604: 
6613:          NRBSITES = NPATCH6605:          NRBSITES = NPATCH
6614:          ALLOCATE(SITE(NRBSITES,3))6606:          ALLOCATE(SITE(NRBSITES,3))
6615: 6607: 


r31861/potential.f90 2017-01-30 19:30:26.718360648 +0000 r31860/potential.f90 2017-01-30 19:30:29.818402596 +0000
 37:    USE MODAMBER9, ONLY: ATMASS1 37:    USE MODAMBER9, ONLY: ATMASS1
 38:    USE OPEP_INTERFACE_MOD, ONLY: OPEP_ENERGY_AND_GRADIENT, OPEP_NUM_HESS 38:    USE OPEP_INTERFACE_MOD, ONLY: OPEP_ENERGY_AND_GRADIENT, OPEP_NUM_HESS
 39:    USE POLIRMOD, ONLY: POLIR 39:    USE POLIRMOD, ONLY: POLIR
 40:    USE MBPOLMOD, ONLY: MBPOL 40:    USE MBPOLMOD, ONLY: MBPOL
 41:    USE SWMOD, ONLY: SWTYPE 41:    USE SWMOD, ONLY: SWTYPE
 42:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS 42:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS
 43:    USE GAUSS_MOD, ONLY: GFIELD 43:    USE GAUSS_MOD, ONLY: GFIELD
 44:    USE RAD_MOD, ONLY: RAD, RADR 44:    USE RAD_MOD, ONLY: RAD, RADR
 45:    USE PREC, ONLY: INT32, REAL64 45:    USE PREC, ONLY: INT32, REAL64
 46:    USE TWIST_MOD, ONLY: TWISTT, TWIST 46:    USE TWIST_MOD, ONLY: TWISTT, TWIST
 47:    USE CONVEX_POLYHEDRA_MODULE, ONLY: CONVEX_POLYHEDRA 
 48:    IMPLICIT NONE 47:    IMPLICIT NONE
 49: ! Arguments 48: ! Arguments
 50: ! TODO: Work out intents 49: ! TODO: Work out intents
 51: ! TODO: Fix array dimensions? 50: ! TODO: Fix array dimensions?
 52:    REAL(REAL64)               :: X(*) 51:    REAL(REAL64)               :: X(*)
 53:    REAL(REAL64)               :: GRAD(*) 52:    REAL(REAL64)               :: GRAD(*)
 54:    REAL(REAL64)               :: EREAL 53:    REAL(REAL64)               :: EREAL
 55:    LOGICAL, INTENT(IN)        :: GRADT 54:    LOGICAL, INTENT(IN)        :: GRADT
 56:    LOGICAL, INTENT(IN)        :: SECT 55:    LOGICAL, INTENT(IN)        :: SECT
 57: ! Variables       56: ! Variables      
1066: 1065: 
1067:    ELSE IF (PAPT) THEN1066:    ELSE IF (PAPT) THEN
1068:       CALL PAP (X, GRAD, EREAL, GRADT)1067:       CALL PAP (X, GRAD, EREAL, GRADT)
1069: 1068: 
1070:    ELSE IF (PAPBINT) THEN1069:    ELSE IF (PAPBINT) THEN
1071:       CALL PAPBIN (X, GRAD, EREAL, GRADT)1070:       CALL PAPBIN (X, GRAD, EREAL, GRADT)
1072: 1071: 
1073:    ELSE IF (PAPJANT) THEN1072:    ELSE IF (PAPJANT) THEN
1074:       CALL PAPJANUS (X, GRAD, EREAL, GRADT)1073:       CALL PAPJANUS (X, GRAD, EREAL, GRADT)
1075: 1074: 
1076:    ELSE IF (POLYT) THEN 
1077:       CALL CONVEX_POLYHEDRA(X, GRAD, EREAL, GRADT) 
1078:  
1079:    ELSE IF (PTSTSTT) THEN1075:    ELSE IF (PTSTSTT) THEN
1080:       CALL PTSTST (X, GRAD, EREAL, GRADT)1076:       CALL PTSTST (X, GRAD, EREAL, GRADT)
1081: 1077: 
1082:    ELSE IF (MULTPAHAT) THEN1078:    ELSE IF (MULTPAHAT) THEN
1083:       CALL MULTPAHA (X, GRAD, EREAL, GRADT)1079:       CALL MULTPAHA (X, GRAD, EREAL, GRADT)
1084: 1080: 
1085:    ELSE IF (PYT .AND. (.NOT.RIGIDCONTOURT)) THEN1081:    ELSE IF (PYT .AND. (.NOT.RIGIDCONTOURT)) THEN
1086:       CALL PY(X, GRAD, EREAL, GRADT)1082:       CALL PY(X, GRAD, EREAL, GRADT)
1087: 1083: 
1088:    ELSE IF (MSTBINT) THEN1084:    ELSE IF (MSTBINT) THEN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0