hdiff output

r33449/rotate_hinge_takestep.f90 2017-11-08 21:30:16.020035807 +0000 r33448/rotate_hinge_takestep.f90 2017-11-08 21:30:16.320039786 +0000
 10:     INTEGER, ALLOCATABLE :: SET_SIZES(:)     ! A list of the number of plates which need to move for each hinge. 10:     INTEGER, ALLOCATABLE :: SET_SIZES(:)     ! A list of the number of plates which need to move for each hinge.
 11:     INTEGER, ALLOCATABLE :: PLATE_SETS(:,:)  ! For each hinge, this will hold the list of plates which  11:     INTEGER, ALLOCATABLE :: PLATE_SETS(:,:)  ! For each hinge, this will hold the list of plates which 
 12:                                              ! need to be moved when the hinge angle is changed.  12:                                              ! need to be moved when the hinge angle is changed. 
 13:  13: 
 14: CONTAINS 14: CONTAINS
 15:  15: 
 16: SUBROUTINE HINGE_INITIALISE 16: SUBROUTINE HINGE_INITIALISE
 17:  17: 
 18:     IMPLICIT NONE 18:     IMPLICIT NONE
 19:  19: 
 20:     INTEGER :: J1, HINGEUNIT, JUNK 20:     INTEGER :: J1
 21:     INTEGER :: GETUNIT 
 22:  21: 
 23:     ! Input file: hingeconfig 22:     ! Input file: hingeconfig
 24:     ! Format as follows. 23:     ! Format as follows.
 25:     ! 24:     !
 26:     ! 25:     !
 27:     ! The first line contains the number of (flexible) hinges in the system. 26:     ! The first line contains the number of (flexible) hinges in the system.
 28:     ! Each hinge is specified by a set of 3 lines. Sets may be separated by a  27:     ! Each hinge is specified by a set of 3 lines. Sets may be separated by a blank line, for clarity, but this is not compulsory.
 29:     ! blank line, for clarity, but this is not compulsory. 
 30:     ! 28:     !
 31:     ! We need to know the following: 29:     ! We need to know the following:
 32:     ! - where the hinge is located (we get this by specifying 4 atom indices, 2 30:     ! - where the hinge is located (we get this by specifying 4 atom indices, 2 of which lie on either side of the hinge)
 33:     !   of which lie on either side of the hinge) 31:     ! - which plates need to move (one side of the hinge is the "moving" side; the plate on that side needs to move, along with all other plates hinged to it)
 34:     ! - which plates need to move (one side of the hinge is the "moving" side; 32:     !
 35:     !   the plate on that side needs to move, along with all other plates hinged 33:     ! The first line of each set contains a number, corresponding to the number of plates on the "moving" side of the hinge.
 36:     !   to it) 34:     ! It will be more efficient to use the side of the hinge with the smallest number of connected plates as the "moving side"
 37:     ! 35:     !
 38:     ! The first line of each set contains a number, corresponding to the number 36:     ! The second line of each set contains a list of indices specifying which plates belong to the moving side. These indices must match the order in which the plates are specified in rbodyconfig.
 39:     ! of plates on the "moving" side of the hinge. It will be more efficient to 37:     !
 40:     ! use the side of the hinge with the smallest number of connected plates as 38:     ! The third line contains a list of 4 atom indices. These are the reference atoms for the hinge. They should fall into 2 pairs
 41:     ! the "moving side" 39:     ! (listed one after the other), with each pair straddling the hinge and joined by a stiff spring.
 42:     ! 40:     ! The two atoms on each side of the hinge should be connected by a vector parallel to the hinge.
 43:     ! The second line of each set contains a list of indices specifying which 41: 
 44:     ! plates belong to the moving side. These indices must match the order in 42: 
 45:     ! which the plates are specified in rbodyconfig. 43:     OPEN(UNIT=22,FILE='hingeconfig',status='old')
 46:     ! 44:     READ(22,*) NHINGES
 47:     ! The third line contains a list of 4 atom indices. These are the reference 
 48:     ! atoms for the hinge. They should fall into 2 pairs (listed one after the 
 49:     ! other), with each pair straddling the hinge and joined by a stiff spring. 
 50:     ! The two atoms on each side of the hinge should be connected by a vector 
 51:     ! parallel to the hinge. 
 52:  
 53:     HINGEUNIT = GETUNIT() 
 54:     OPEN(UNIT=HINGEUNIT,FILE='hingeconfig',status='old') 
 55:     READ(HINGEUNIT,*) NHINGES 
 56:     IF(DEBUG) WRITE(MYUNIT,*) "Reading hingeconfig. NHINGES:", NHINGES 45:     IF(DEBUG) WRITE(MYUNIT,*) "Reading hingeconfig. NHINGES:", NHINGES
 57:     ! We may eventually want to move to 6 reference atoms so as to compute the 46:     ALLOCATE(REFATOMS(NHINGES,4))  ! We may eventually want to move to 6 reference atoms so as to compute the angle between the two hinges as well.
 58:     ! angle between the two hinges as well. 47:     ALLOCATE(PLATE_SETS(NHINGES,1))
 59:     ALLOCATE(REFATOMS(NHINGES,4)) 
 60:     ALLOCATE(SET_SIZES(NHINGES)) 48:     ALLOCATE(SET_SIZES(NHINGES))
 61:  49: 
 62: ! jwrm2> First go through the whole file to find the maximum number of plates in 50:     DO J1=1,NHINGES
 63: !        a set 51:        READ(22,*) SET_SIZES(J1)
 64:     DO J1 = 1, NHINGES 
 65:         ! Number of plates this hinge moves 
 66:         READ(HINGEUNIT, *) SET_SIZES(J1) 
 67:         ! Identity of plates to move, gobble this time sice the array to store 
 68:         ! them isn't allocated 
 69:         READ(HINGEUNIT, *) JUNK 
 70:         ! Identity of reference atoms 
 71:         READ(HINGEUNIT, *) REFATOMS(J1, :) 
 72:         IF(DEBUG) THEN 
 73:             WRITE(MYUNIT,*) "Reference atoms:" 
 74:             WRITE(MYUNIT,*) REFATOMS(J1,:) 
 75:         ENDIF 
 76:     END DO 
 77:  52: 
 78:     IF (DEBUG) THEN 53:        IF(DEBUG) WRITE(MYUNIT,'(A,I2,A)') "Hinge ", J1,". Moving side:"
 79:         WRITE (MYUNIT, *) "SET_SIZES = ", SET_SIZES(:) 
 80:     END IF 
 81:  
 82: ! jwrm2> Now we can allocate enough space for the plate sets 
 83:     ALLOCATE(PLATE_SETS(NHINGES, MAXVAL(SET_SIZES))) 
 84:  
 85: ! jwrm2> Read the file again to get the plate sets 
 86:     REWIND (UNIT=HINGEUNIT) 
 87:  
 88: ! jwrm2> Gobble number of hinges 
 89:     READ (HINGEUNIT, *) JUNK 
 90:  
 91:     DO J1 = 1, NHINGES 
 92:         ! Gobble set size 
 93:         READ(HINGEUNIT,*) JUNK 
 94:  
 95:         ! Get plate set 
 96:         IF(DEBUG) WRITE(MYUNIT,'(A,I2,A)') "Hinge ", J1, ". Moving side:" 
 97:         READ(HINGEUNIT,*) PLATE_SETS(J1,1:SET_SIZES(J1)) 
 98:         IF(DEBUG) WRITE(MYUNIT,*) PLATE_SETS(J1,1:SET_SIZES(J1)) 
 99:  54: 
100:         ! Gobble reference atoms 55:        READ(22,*) PLATE_SETS(J1,1:SET_SIZES(J1))
101:         READ(HINGEUNIT,*) JUNK 56: 
102:     END DO 57:        IF(DEBUG) WRITE(MYUNIT,*) PLATE_SETS(J1,:)
  58: 
  59:        READ(22,*) REFATOMS(J1,:)
103:  60: 
104:     CLOSE(HINGEUNIT) 61:        IF(DEBUG) THEN
  62:            WRITE(MYUNIT,*) "Reference atoms:"
  63:            WRITE(MYUNIT,*) REFATOMS(J1,:)
  64:        ENDIF
  65: 
  66:     END DO
  67:     CLOSE(22)
105:  68: 
106: END SUBROUTINE HINGE_INITIALISE 69: END SUBROUTINE HINGE_INITIALISE
107:  70: 
108: SUBROUTINE HINGE_ROTATE(XCOORDS, ROTATEFACTOR) 71: SUBROUTINE HINGE_ROTATE(XCOORDS, ROTATEFACTOR)
109:  72: 
110:     USE COMMONS, ONLY: NATOMS 73:     USE COMMONS, ONLY: NATOMS
111:     USE ROTATIONS, ONLY: rot_aa2mx 74:     USE ROTATIONS, ONLY: rot_aa2mx
112:     USE GENRIGID, ONLY: NSITEPERBODY, RIGIDGROUPS 75:     USE GENRIGID, ONLY: NSITEPERBODY, RIGIDGROUPS
113:  76: 
114:     IMPLICIT NONE 77:     IMPLICIT NONE
115:  78: 
116:     INTEGER :: J1, J2, J3, ATOM, PLATE_INDEX 79:     INTEGER :: J1, J2, J3, ATOM, PLATE_INDEX
117:     DOUBLE PRECISION :: REFPOINT(3), PI, TWOPI, DPRAND 80:     DOUBLE PRECISION :: REFPOINT(3), PI, TWOPI, DPRAND
118:     DOUBLE PRECISION :: AXIS(3), MATRIX(3,3), TOROTATE(3), ATOMROTATED(3) 81:     DOUBLE PRECISION :: AXIS(3), MATRIX(3,3), TOROTATE(3), ATOMROTATED(3)
119:     DOUBLE PRECISION :: ANGLE, NORM 82:     DOUBLE PRECISION :: ANGLE, NORM
120:     DOUBLE PRECISION, INTENT(INOUT) :: XCOORDS(3*NATOMS) 83:     DOUBLE PRECISION, INTENT(INOUT) :: XCOORDS(3*NATOMS)
121:     DOUBLE PRECISION, INTENT(IN) :: ROTATEFACTOR 84:     DOUBLE PRECISION, INTENT(IN) :: ROTATEFACTOR
122:  85: 
123:     ! Current implementation - somewhat crude - is to change every hinge every 86:     ! Current implementation - somewhat crude - is to change every hinge every time the subroutine is called. We pick a random rotation angle between -pi*rotatefactor and pi*rotatefactor
124:     ! time the subroutine is called. We pick a random rotation angle between 
125:     ! -pi*rotatefactor and pi*rotatefactor 
126:  87: 
127: !    write(*,*) "Before hinge takestep" 88: !    write(*,*) "Before hinge takestep"
128: !    write(*,*) XCOORDS 89: !    write(*,*) XCOORDS
129:  90: 
130:     MATRIX(:,:) = 0.0D0 91:     MATRIX(:,:) = 0.0D0
131:     TOROTATE(:) = 0.0D0 92:     TOROTATE(:) = 0.0D0
132:     ! Define some constants 93:     ! Define some constants
133:     PI=ATAN(1.0D0)*4 94:     PI=ATAN(1.0D0)*4
134:     TWOPI=2.0D0*PI 95:     TWOPI=2.0D0*PI
135:  96: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0