hdiff output

r29884/keywords.f 2016-02-01 12:30:15.365472762 +0000 r29883/keywords.f 2016-02-01 12:30:15.953480585 +0000
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 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 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 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19:       SUBROUTINE KEYWORD 19:       SUBROUTINE KEYWORD
 20: !      USE NEW_INPUT_MOD, ONLY : NEW_INPUT 20: !      USE NEW_INPUT_MOD, ONLY : NEW_INPUT
 21:       use COMMONS 21:       use COMMONS
 22:       use VEC3 22:       use VEC3
 23:       use genrigid 23:       use genrigid
 24:       use multipot 
 25:       use MODMXATMS   ! NEEDED FOR charmm 24:       use MODMXATMS   ! NEEDED FOR charmm
 26:       USE modcharmm 25:       USE modcharmm
 27:       USE TWIST_MOD 26:       USE TWIST_MOD
 28: !       sf344> AMBER additions 27: !       sf344> AMBER additions
 29:       USE modamber9, only : coords1,amberstr,amberstr1,mdstept,inpcrd,amberenergiest, nocistransdna, nocistransrna, 28:       USE modamber9, only : coords1,amberstr,amberstr1,mdstept,inpcrd,amberenergiest, nocistransdna, nocistransrna,
 30:      &                      uachiral, ligrotscale, setchiral, STEEREDMINT, SMINATOMA, SMINATOMB, SMINK, SMINKINC, 29:      &                      uachiral, ligrotscale, setchiral, STEEREDMINT, SMINATOMA, SMINATOMB, SMINK, SMINKINC,
 31:      &                      SMINDISTSTART, SMINDISTFINISH, natomsina, natomsinb, natomsinc, atomsinalist, atomsinblist, 30:      &                      SMINDISTSTART, SMINDISTFINISH, natomsina, natomsinb, natomsinc, atomsinalist, atomsinblist,
 32:      &                      atomsinclist, atomsinalistlogical, atomsinblistlogical, atomsinclistlogical, ligcartstep, 31:      &                      atomsinclist, atomsinalistlogical, atomsinblistlogical, atomsinclistlogical, ligcartstep,
 33:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange, 32:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange,
 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT, 33:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT,
 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 40:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 42:       USE MBPOLMOD, ONLY: MBPOLINIT 41:       USE MBPOLMOD, ONLY: MBPOLINIT
 43:       USE SWMOD, ONLY: SWINIT, MWINIT 42:       USE SWMOD, ONLY: SWINIT, MWINIT
 44:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS, 43:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS,
 45:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA 44:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA
 46:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 45:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 47:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 46:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 48:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,  47:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 
 49:      &     PARSE_MLJ_PARAMS 48:      &     PARSE_MLJ_PARAMS
 50:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT 49:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT
 51:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE 
 52:       IMPLICIT NONE 50:       IMPLICIT NONE
 53:  51: 
 54:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, INDEX, J2, J3, J4 52:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, INDEX, J2, J3, J4
 55:       INTEGER DATA_UNIT 53:       INTEGER DATA_UNIT
 56:       INTEGER MOVABLEATOMINDEX 54:       INTEGER MOVABLEATOMINDEX
 57:       LOGICAL CAT, YESNO, PERMFILE, CONFILE 55:       LOGICAL CAT, YESNO, PERMFILE, CONFILE
 58:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR, 56:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR,
 59:      &                NERROR, ECHO, LAST, CAT 57:      &                NERROR, ECHO, LAST, CAT
 60:        DOUBLE PRECISION XX, ROH, ROM, WTHETA  58:        DOUBLE PRECISION XX, ROH, ROM, WTHETA 
 61:       LOGICAL END, SKIPBL, CLEAR, ECHO 59:       LOGICAL END, SKIPBL, CLEAR, ECHO
684:       CHSTANDARDT=.FALSE.682:       CHSTANDARDT=.FALSE.
685:       CHLOOPMODT=.FALSE.683:       CHLOOPMODT=.FALSE.
686:       CHCARTMOVET=.FALSE.684:       CHCARTMOVET=.FALSE.
687:       CHCLUSTERT=.FALSE.685:       CHCLUSTERT=.FALSE.
688:       CHNEIGHBOURT=.FALSE.686:       CHNEIGHBOURT=.FALSE.
689:       CHBBT=.FALSE.687:       CHBBT=.FALSE.
690:       CHSCT=.FALSE.688:       CHSCT=.FALSE.
691:       SECPREDT=.FALSE.689:       SECPREDT=.FALSE.
692:       SECPREDFILE='UNDEFINED'690:       SECPREDFILE='UNDEFINED'
693: 691: 
694: ! sn402: hinge moves for the plate potential 
695:       ROTATEHINGET=.FALSE. 
696:       HINGEROTATEFACTOR=1.0D0 
697:       ROTATEHINGEFREQ=1 
698:  
699:       OSASAT=.FALSE.692:       OSASAT=.FALSE.
700:       RPRO=1.4D0693:       RPRO=1.4D0
701:       ODIHET=.FALSE.694:       ODIHET=.FALSE.
702:       ORGYT=.FALSE.695:       ORGYT=.FALSE.
703:       OEINTT=.FALSE.696:       OEINTT=.FALSE.
704:       MON1(1:2)=1697:       MON1(1:2)=1
705:       MON2(1:2)=1698:       MON2(1:2)=1
706: 699: 
707:       BSMIN=.FALSE.700:       BSMIN=.FALSE.
708:       RKMIN=.FALSE.701:       RKMIN=.FALSE.
991:       RELAXRIGIDT = .FALSE.984:       RELAXRIGIDT = .FALSE.
992:       NRELAXRIGIDR = 1000000000985:       NRELAXRIGIDR = 1000000000
993:       NRELAXRIGIDA = 1000000000986:       NRELAXRIGIDA = 1000000000
994: 987: 
995:       RIGIDOPTIMROTAT = .FALSE.988:       RIGIDOPTIMROTAT = .FALSE.
996:       OPTIMROTAVALUES(:) = 0.0D0989:       OPTIMROTAVALUES(:) = 0.0D0
997:       FREEZERIGIDBODYT = .FALSE.      990:       FREEZERIGIDBODYT = .FALSE.      
998: 991: 
999:       AACONVERGENCET = .FALSE.992:       AACONVERGENCET = .FALSE.
1000: 993: 
1001: ! sn402 > Multiple potential scheme 
1002:       MULTIPOTT = .FALSE. 
1003:  
1004: !--------------------------------!994: !--------------------------------!
1005: ! hk286 > Generalised Thomson    !995: ! hk286 > Generalised Thomson    !
1006: !--------------------------------!996: !--------------------------------!
1007:       GTHOMSONT = .FALSE.997:       GTHOMSONT = .FALSE.
1008:       GTHOMPOT = 1998:       GTHOMPOT = 1
1009: 999: 
1010: ! hk286 > Damped group moves1000: ! hk286 > Damped group moves
1011:       DAMPEDGMOVET = .FALSE.1001:       DAMPEDGMOVET = .FALSE.
1012:       DMOVEFREQ = 11002:       DMOVEFREQ = 1
1013: 1003: 
4422:          MSORIGT=.TRUE.4412:          MSORIGT=.TRUE.
4423: 4413: 
4424:       ELSE IF (WORD.EQ.'MSTRANS') THEN4414:       ELSE IF (WORD.EQ.'MSTRANS') THEN
4425:          MSTRANST=.TRUE.4415:          MSTRANST=.TRUE.
4426: 4416: 
4427:       ELSE IF (WORD.EQ.'MULLERBROWN') THEN4417:       ELSE IF (WORD.EQ.'MULLERBROWN') THEN
4428:          MULLERBROWNT=.TRUE.4418:          MULLERBROWNT=.TRUE.
4429: 4419: 
4430:       ELSE IF (WORD.EQ.'MULTIPLICITY') THEN4420:       ELSE IF (WORD.EQ.'MULTIPLICITY') THEN
4431:          CALL READI(XMUL)4421:          CALL READI(XMUL)
4432:  
4433:       ELSE IF (WORD .EQ. 'MULTIPOT') THEN 
4434:          MULTIPOTT = .TRUE. 
4435:          CALL MULTIPOT_INITIALISE         
4436:           
4437: !4422: !
4438: ! NOT DOCUMENTED4423: ! NOT DOCUMENTED
4439: !4424: !
4440:       ELSE IF (WORD.EQ.'MYSD') THEN4425:       ELSE IF (WORD.EQ.'MYSD') THEN
4441:          MYSDT=.TRUE.4426:          MYSDT=.TRUE.
4442:          CALL READF(SDTOL)4427:          CALL READF(SDTOL)
4443: 4428: 
4444:       ELSE IF (WORD.EQ.'NATB') THEN4429:       ELSE IF (WORD.EQ.'NATB') THEN
4445:          NATBT=.TRUE.4430:          NATBT=.TRUE.
4446: 4431: 
5480:          CALL READI(RMSSAVE)5465:          CALL READI(RMSSAVE)
5481:          CALL READI(J1)5466:          CALL READI(J1)
5482:          CALL READI(J2)5467:          CALL READI(J2)
5483:          IF(J1.EQ.1) THEN5468:          IF(J1.EQ.1) THEN
5484:            SELECTT=.TRUE.5469:            SELECTT=.TRUE.
5485:          ELSE5470:          ELSE
5486:            SELECTT=.FALSE.5471:            SELECTT=.FALSE.
5487:          ENDIF5472:          ENDIF
5488:          IF (J2.EQ.1) PROGRESS=.TRUE.5473:          IF (J2.EQ.1) PROGRESS=.TRUE.
5489:          WRITE(MYUNIT,'(A)') 'RMST set'5474:          WRITE(MYUNIT,'(A)') 'RMST set'
5490:  
5491:        
5492:        ELSE IF (WORD.EQ.'ROTATEHINGE') THEN 
5493:           ROTATEHINGET=.TRUE. 
5494: !         Read ROTATEHINGEFREQ  
5495:           IF (NITEMS.GT.1) CALL READI(ROTATEHINGEFREQ) 
5496: !         Read in ROTATEFACTOR 
5497:           IF (NITEMS.GT.2) CALL READF(HINGEROTATEFACTOR) 
5498:           WRITE(MYUNIT,'(A)') ' keywords> Enabled hinge rotation moves for the plate potential' 
5499:           WRITE(MYUNIT,'(A,I2,A)') ' keywords> Hinges will be rotated every ',ROTATEHINGEFREQ,' steps' 
5500:           WRITE(MYUNIT,'(A,F20.10)') ' keywords> Hinge angle rotatefactor =',HINGEROTATEFACTOR 
5501:           CALL HINGE_INITIALISE 
5502:           
5503: !5475: !
5504: ! csw34> Rigid body rotation moves. Each rigid body is randomly rotated about its COM every ROTATERIGIDFREQ steps.5476: ! csw34> Rigid body rotation moves. Each rigid body is randomly rotated about its COM every ROTATERIGIDFREQ steps.
5505: !        ROTATEFACTOR scales the maximum rotation with 1.0 being complete freedom to rotate.5477: !        ROTATEFACTOR scales the maximum rotation with 1.0 being complete freedom to rotate.
5506: !5478: !
5507:       ELSE IF (WORD.EQ.'ROTATERIGID') THEN5479:       ELSE IF (WORD.EQ.'ROTATERIGID') THEN
5508:          ROTATERIGIDT=.TRUE.5480:          ROTATERIGIDT=.TRUE.
5509: ! Read ROTATERIGIDFREQ 5481: ! Read ROTATERIGIDFREQ 
5510:          IF (NITEMS.GT.1) CALL READI(ROTATERIGIDFREQ)5482:          IF (NITEMS.GT.1) CALL READI(ROTATERIGIDFREQ)
5511: ! Read in ROTATEFACTOR5483: ! Read in ROTATEFACTOR
5512:          IF (NITEMS.GT.2) CALL READF(ROTATEFACTOR)5484:          IF (NITEMS.GT.2) CALL READF(ROTATEFACTOR)


r29884/mc.F 2016-02-01 12:30:15.561475370 +0000 r29883/mc.F 2016-02-01 12:30:16.153483245 +0000
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 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 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 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19:       SUBROUTINE MC(NSTEPS,SCALEFAC,SCREENC) 19:       SUBROUTINE MC(NSTEPS,SCALEFAC,SCREENC)
 20:       USE COMMONS 20:       USE COMMONS
 21:       use genrigid 21:       use genrigid
 22:       USE HINGE_MOVES 
 23:       use moves 22:       use moves
 24:       use grouprotmod 23:       use grouprotmod
 25:       use mbpolmod, only: mbpolstep 24:       use mbpolmod, only: mbpolstep
 26:  25: 
 27:       USE QMODULE , ONLY : QMIN, QMINP, INTEQMIN 26:       USE QMODULE , ONLY : QMIN, QMINP, INTEQMIN
 28:       USE modcharmm 27:       USE modcharmm
 29:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA, 28:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA,
 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT, 29:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT,
 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL, 30:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL,
 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB 31:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB
760:                        DOEXPANDRIGID=.TRUE.759:                        DOEXPANDRIGID=.TRUE.
761:                ENDIF760:                ENDIF
762: ! csw34> rigid body rotation move switch761: ! csw34> rigid body rotation move switch
763:                IF(ROTATERIGIDT.AND.MOD(J1,ROTATERIGIDFREQ).EQ.0) THEN762:                IF(ROTATERIGIDT.AND.MOD(J1,ROTATERIGIDFREQ).EQ.0) THEN
764:                        DOROTATERIGID=.TRUE.763:                        DOROTATERIGID=.TRUE.
765:                ENDIF764:                ENDIF
766: ! mo361> rigid body rotation move switch765: ! mo361> rigid body rotation move switch
767:                IF(TRANSLATERIGIDT.AND.MOD(J1,TRANSLATERIGIDFREQ).EQ.0) THEN766:                IF(TRANSLATERIGIDT.AND.MOD(J1,TRANSLATERIGIDFREQ).EQ.0) THEN
768:                        DOTRANSLATERIGID=.TRUE.767:                        DOTRANSLATERIGID=.TRUE.
769:                ENDIF768:                ENDIF
770:  
771:                IF(ROTATEHINGET .AND. MOD(J1,ROTATEHINGEFREQ).EQ.0) THEN 
772:                        DOROTATEHINGE=.TRUE. 
773:                ENDIF 
774:  
775: !769: !
776: ! csw34> Coordinates are saved so that moves can be undone770: ! csw34> Coordinates are saved so that moves can be undone
777: !771: !
778:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)772:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)
779: ! csw34> If you want to look at the effect of moves, you can dump out773: ! csw34> If you want to look at the effect of moves, you can dump out
780: ! the structure BEFORE the move here.774: ! the structure BEFORE the move here.
781: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")775: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")
782: !                 CALL CHARMMDUMP(COORDS(:,JP),'beforemove')776: !                 CALL CHARMMDUMP(COORDS(:,JP),'beforemove')
783: 777: 
784: !778: !
790:                         CALL GENRIGID_EXPAND(COORDS(:,JP), EXPANDFACTOR) 784:                         CALL GENRIGID_EXPAND(COORDS(:,JP), EXPANDFACTOR) 
791:                      ENDIF785:                      ENDIF
792: ! csw34> Rigid body rotation moves786: ! csw34> Rigid body rotation moves
793:                      IF (RIGIDINIT.AND.ROTATERIGIDT.AND.DOROTATERIGID) THEN787:                      IF (RIGIDINIT.AND.ROTATERIGIDT.AND.DOROTATERIGID) THEN
794:                         CALL GENRIGID_ROTATE(COORDS(:,JP), ROTATEFACTOR)788:                         CALL GENRIGID_ROTATE(COORDS(:,JP), ROTATEFACTOR)
795:                      ENDIF789:                      ENDIF
796: ! mo361> Rigid body translation moves790: ! mo361> Rigid body translation moves
797:                      IF (RIGIDINIT.AND.TRANSLATERIGIDT.AND.DOTRANSLATERIGID) THEN791:                      IF (RIGIDINIT.AND.TRANSLATERIGIDT.AND.DOTRANSLATERIGID) THEN
798:                         CALL GENRIGID_TRANSLATE(COORDS(:,JP),TRANSLATEFACTOR)792:                         CALL GENRIGID_TRANSLATE(COORDS(:,JP),TRANSLATEFACTOR)
799:                      ENDIF793:                      ENDIF
800:                      IF (RIGIDINIT.AND.ROTATEHINGET.AND.DOROTATEHINGE) THEN 
801:                         CALL HINGE_ROTATE(COORDS(:,JP), HINGEROTATEFACTOR) 
802:                      ENDIF 
803: !794: !
804: ! CHARMM STEP TAKING795: ! CHARMM STEP TAKING
805: !796: !
806:                IF (CHRMMT) THEN797:                IF (CHRMMT) THEN
807:                   IF (CHMDT.AND.MOD(J1,CHMDFREQ).EQ.0) THEN798:                   IF (CHMDT.AND.MOD(J1,CHMDFREQ).EQ.0) THEN
808:                      CALL CHMD(JP)799:                      CALL CHMD(JP)
809:                   ELSE800:                   ELSE
810:                      IF (CHRIGIDTRANST.AND.MOD(J1,FTRANS).EQ.0) CALL MKRIGIDTRANS(JP)801:                      IF (CHRIGIDTRANST.AND.MOD(J1,FTRANS).EQ.0) CALL MKRIGIDTRANS(JP)
811:                      IF (CHRIGIDROTT.AND.MOD(J1,FROT).EQ.0) CALL MKRIGIDROT(JP)802:                      IF (CHRIGIDROTT.AND.MOD(J1,FROT).EQ.0) CALL MKRIGIDROT(JP)
812:                      IF (MOD(J1,CHFREQ).EQ.0) THEN803:                      IF (MOD(J1,CHFREQ).EQ.0) THEN


r29884/rotate_hinge_takestep.f90 2016-02-01 12:30:15.753477924 +0000 r29883/rotate_hinge_takestep.f90 2016-02-01 12:30:16.341485747 +0000
  1: MODULE HINGE_MOVES  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/rotate_hinge_takestep.f90' in revision 29883
  2:  
  3:     USE COMMONS, ONLY: DEBUG, MYUNIT 
  4:  
  5:     IMPLICIT NONE 
  6:  
  7:     INTEGER :: NHINGES                       ! The number of hinges in the system 
  8:     INTEGER, ALLOCATABLE :: REFATOMS(:,:)    ! Will hold a set of 4 reference atoms for each hinge 
  9:     INTEGER, ALLOCATABLE :: SET_SIZES(:)     ! A list of the number of plates which need to move for each hinge. 
 10:     INTEGER, ALLOCATABLE :: PLATE_SETS(:,:)  ! For each hinge, this will hold the list of plates which  
 11:                                              ! need to be moved when the hinge angle is changed.  
 12:  
 13: CONTAINS 
 14:  
 15:  
 16: ! Eventually this will contain code for parsing an input file setting up the hinges. 
 17: ! For now, I'm hard-coding in the net for the symmetrical tetrahedron. 
 18: SUBROUTINE HINGE_INITIALISE 
 19:  
 20:     IMPLICIT NONE 
 21:  
 22:     INTEGER :: J1 
 23:  
 24:     ! Input file: hingeconfig 
 25:     ! Format as follows. 
 26:     ! 
 27:     ! 
 28:     ! The first line contains the number of (flexible) hinges in the system. 
 29:     ! Each hinge is specified by a set of 3 lines. Sets may be separated by a blank line, for clarity, but this is not compulsory. 
 30:     ! 
 31:     ! We need to know the following: 
 32:     ! - where the hinge is located (we get this by specifying 4 atom indices, 2 of which lie on either side of the hinge) 
 33:     ! - which plates need to move (one side of the hinge is the "moving" side; the plate on that side needs to move, along with all other plates hinged to it) 
 34:     ! 
 35:     ! The first line of each set contains a number, corresponding to the number of plates on the "moving" side of the hinge. 
 36:     ! It will be more efficient to use the side of the hinge with the smallest number of connected plates as the "moving side" 
 37:     ! 
 38:     ! The second line of each set contains a list of indices specifying which hinges will move. These indices must match the order in which the plates are specified in rbodyconfig. 
 39:     ! 
 40:     ! The third line contains a list of 4 atom indices. These are the reference atoms for the hinge. They should fall into 2 pairs 
 41:     ! (listed one after the other), with each pair straddling the hinge and joined by a stiff spring. 
 42:     ! The two atoms on one side of the hinge should be connected by a vector parallel to the hinge. 
 43:  
 44:  
 45:     OPEN(UNIT=22,FILE='hingeconfig',status='old') 
 46:     READ(22,*) NHINGES 
 47:  
 48:     ALLOCATE(REFATOMS(NHINGES,4))  ! We may eventually want to move to 6 reference atoms so as to compute the angle between the two hinges as well. 
 49:     ALLOCATE(PLATE_SETS(NHINGES,1)) 
 50:     ALLOCATE(SET_SIZES(NHINGES)) 
 51:  
 52:     DO J1=1,NHINGES 
 53:        READ(22,*) SET_SIZES(J1) 
 54:        READ(22,*) PLATE_SETS(J1,1:SET_SIZES(J1)) 
 55:        READ(22,*) REFATOMS(J1,:)  
 56:     END DO 
 57:     CLOSE(22) 
 58:  
 59:     IF(DEBUG) THEN 
 60:         WRITE(MYUNIT,*) "Read hingeconfig. NHINGES:", NHINGES 
 61:         DO J1=1,NHINGES 
 62:             WRITE(MYUNIT,'(A,I2,A)') "Hinge ", J1,". Moving side:" 
 63:             WRITE(MYUNIT,*) PLATE_SETS(J1,:) 
 64:             WRITE(MYUNIT,*) "Reference atoms:" 
 65:             WRITE(MYUNIT,*) REFATOMS(J1,:) 
 66:         ENDDO 
 67:     ENDIF 
 68:  
 69: END SUBROUTINE HINGE_INITIALISE 
 70:  
 71: SUBROUTINE HINGE_ROTATE(XCOORDS, ROTATEFACTOR) 
 72:  
 73:     USE COMMONS, ONLY: NATOMS 
 74:     USE ROTATIONS, ONLY: rot_aa2mx 
 75:     USE GENRIGID, ONLY: NSITEPERBODY, RIGIDGROUPS 
 76:  
 77:     IMPLICIT NONE 
 78:  
 79:     INTEGER :: J1, J2, J3, ATOM, PLATE_INDEX 
 80:     DOUBLE PRECISION :: REFPOINT(3), PI, TWOPI, DPRAND 
 81:     DOUBLE PRECISION :: AXIS(3), MATRIX(3,3), TOROTATE(3), ATOMROTATED(3) 
 82:     DOUBLE PRECISION :: ANGLE, NORM 
 83:     DOUBLE PRECISION, INTENT(INOUT) :: XCOORDS(3*NATOMS) 
 84:     DOUBLE PRECISION, INTENT(IN) :: ROTATEFACTOR 
 85:  
 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 
 87:  
 88:     write(*,*) "Before hinge takestep" 
 89:     write(*,*) XCOORDS 
 90:  
 91:     MATRIX(:,:) = 0.0D0 
 92:     TOROTATE(:) = 0.0D0 
 93:     ! Define some constants 
 94:     PI=ATAN(1.0D0)*4 
 95:     TWOPI=2.0D0*PI 
 96:  
 97:     ! Loop over all hinges 
 98:     DO J1 = 1, NHINGES 
 99:  
100:         ! Calculate the centrepoint of the hinge 
101:         REFPOINT(:) = 0.0D0 
102:         DO J2 = 1,4 
103:             ATOM = REFATOMS(J1,J2) 
104:             REFPOINT(:) = REFPOINT(:) + XCOORDS(3*ATOM-2:3*ATOM) 
105:         ENDDO 
106:         REFPOINT(:) = REFPOINT(:)/4.0D0 
107:  
108:         ! Calculate the hinge axis (from the centrepoint of reference atoms 1 and 2 to the centrepoint of 3 and 4) 
109:         AXIS(:) = XCOORDS(3*REFATOMS(J1,3)-2:3*REFATOMS(J1,3)) + XCOORDS(3*REFATOMS(J1,4)-2:3*REFATOMS(J1,4)) - & 
110:                   XCOORDS(3*REFATOMS(J1,1)-2:3*REFATOMS(J1,1)) - XCOORDS(3*REFATOMS(J1,2)-2:3*REFATOMS(J1,2)) 
111: !       AXIS(:) = AXIS(:) * 0.5D0   ! This line is needed for the definition of the centrepoints, but since 
112:                                     ! we are about to normalise anyway, we don't really need it. 
113:         ! Normalise the axis 
114:         NORM = SQRT(DOT_PRODUCT(AXIS,AXIS)) 
115:         AXIS(:) = AXIS(:)/NORM 
116:  
117:         ! Centre the system about the hinge centrepoint 
118:         DO J2 = 1,NATOMS 
119:             XCOORDS(3*J2-2:3*J2) = XCOORDS(3*J2-2:3*J2) - REFPOINT(:) 
120:         ENDDO 
121:  
122: !        write(*,*) "After first translation" 
123: !        write(*,*) XCOORDS 
124:  
125:         ! Choose a random angle by which to rotate 
126:         ANGLE = (DPRAND()-0.5)*TWOPI*ROTATEFACTOR 
127:  
128:         ! Create an angle-axis vector corresponding to the desired rotation, and convert it to a rotation matrix 
129:         AXIS(:) = AXIS(:)*ANGLE 
130:         !write(*,*) "Angle Axis:", AXIS(:) 
131:         MATRIX = rot_aa2mx(AXIS) 
132:         !write(*,*) "Matrix" 
133:         !DO J2=1,3 
134:         !    write(*,*) MATRIX(J2,:) 
135:         !ENDDO 
136:  
137:         ! Apply the rotation matrix to the atoms in the rigid bodies indicated by PLATE_SETS(J1) 
138:         DO J2 = 1, SET_SIZES(J1) 
139:             PLATE_INDEX = PLATE_SETS(J1,J2) 
140:             DO J3 = 1, NSITEPERBODY(PLATE_INDEX) 
141:                 ATOM = RIGIDGROUPS(J3, PLATE_INDEX) 
142:                 !write(*,'(A,I3,A,I1,A,I1)') "Rotating atom ", ATOM, " in RB ", PLATE_INDEX, " for hinge ", J1 
143:                 TOROTATE = XCOORDS(3*ATOM-2:3*ATOM) 
144:                 !write(*,'(A,3F20.10)') "Before rotation: ", TOROTATE(1), TOROTATE(2), TOROTATE(3) 
145:                 ATOMROTATED=MATMUL(MATRIX,TOROTATE) 
146:                 !write(*,'(A,3F20.10)') "After rotation: ", ATOMROTATED(1), ATOMROTATED(2), ATOMROTATED(3) 
147:                 XCOORDS(3*ATOM-2:3*ATOM) = ATOMROTATED 
148:             ENDDO 
149:         ENDDO 
150:  
151: !        write(*,*) "After rotation, before back-translation" 
152: !        write(*,*) XCOORDS 
153:  
154:         ! Translate the system back to the old centre of mass 
155:         DO J2 = 1, NATOMS 
156:             XCOORDS(3*J2-2:3*J2) = XCOORDS(3*J2-2:3*J2) + REFPOINT(:) 
157:         ENDDO 
158:  
159: !        write(*,*) "After hinge", J1 
160: !        write(*,*) XCOORDS 
161:  
162:  
163:     ENDDO  ! End of loop over hinges 
164:  
165:     write(*,*) "After hinge takestep" 
166:     write(*,*) XCOORDS 
167: END SUBROUTINE HINGE_ROTATE 
168:  
169: END MODULE 
170:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0