hdiff output

r30137/nnutils.f90 2016-03-16 18:34:18.611524073 +0000 r30136/nnutils.f90 2016-03-16 18:34:18.807526092 +0000
946: 946: 
947:         ! Lists of parallel and perpendicular quaternions for body J1 in all the images.947:         ! Lists of parallel and perpendicular quaternions for body J1 in all the images.
948:         ALLOCATE(QTN(NIMAGE+2,4))948:         ALLOCATE(QTN(NIMAGE+2,4))
949:         ALLOCATE(PTN(NIMAGE+2,4))949:         ALLOCATE(PTN(NIMAGE+2,4))
950: 950: 
951: !       CONVERT FROM AA -> QUATERNION FOR THE INITIAL AND FINAL FRAMES951: !       CONVERT FROM AA -> QUATERNION FOR THE INITIAL AND FINAL FRAMES
952:         DO K = 1, (NIMAGE+2), (NIMAGE+1)952:         DO K = 1, (NIMAGE+2), (NIMAGE+1)
953:            ! P is the angle-axis vector for rigid body J1 in the initial or final frame953:            ! P is the angle-axis vector for rigid body J1 in the initial or final frame
954:            J = (K-1)*NOPT + RBOFFSET + 3*(J1-1)954:            J = (K-1)*NOPT + RBOFFSET + 3*(J1-1)
955:            P(:)   = XYZ(J+1:J+3)955:            P(:)   = XYZ(J+1:J+3)
956:  
957:            THETA  = DSQRT(DOT_PRODUCT(P(:),P(:)))956:            THETA  = DSQRT(DOT_PRODUCT(P(:),P(:)))
958:  
959:            IF (ABS(THETA) .LT. 1.0D-8) THEN957:            IF (ABS(THETA) .LT. 1.0D-8) THEN
960:               QTN(K,1)   = 1.D0958:               QTN(K,1)   = 1.D0
961:               QTN(K,2:4) = 0.D0959:               QTN(K,2:4) = 0.D0
962:            ELSE960:            ELSE
963:               THETAH     = 0.5D0*THETA961:               THETAH     = 0.5D0*THETA
964:               ST         = SIN(THETAH)962:               ST         = SIN(THETAH)
965:               QTN(K,1)   = COS(THETAH)963:               QTN(K,1)   = COS(THETAH)
966:               QTN(K,2:4) = P(:)*ST/THETA964:               QTN(K,2:4) = P(:)*ST/THETA
967:            ENDIF965:            ENDIF
968:         ENDDO966:         ENDDO
984:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)982:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)
985:          ENDIF983:          ENDIF
986: 984: 
987:         ! We may sometimes wish to interpolate around the longer of the two angles, for example to985:         ! We may sometimes wish to interpolate around the longer of the two angles, for example to
988:         ! avoid clash of two rigid bodies. Inverting the final quarternion causes iSLERP to986:         ! avoid clash of two rigid bodies. Inverting the final quarternion causes iSLERP to
989:         ! interpolate in the other direction, so this is what we do here.987:         ! interpolate in the other direction, so this is what we do here.
990:          IF (LONG_WAY) THEN988:          IF (LONG_WAY) THEN
991:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)989:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)
992:            CT =-CT  ! Inverting the quarternion obviously inverts the inner product990:            CT =-CT  ! Inverting the quarternion obviously inverts the inner product
993:          ENDIF991:          ENDIF
994: 992:  
995:          IF (ABS(CT-1.0D0) .LT. 1.0D-10) THEN 993:          THETA    = ACOS(CT)
996:             ! For some reason the ACOS function requires an argument with a magnitude < 1 
997:             ! So if this happens, we set THETA directly 
998:             THETA = 0.0D0 
999:          ELSEIF (ABS(-1*CT-1.0D0) .LT. 1.0D-10) THEN 
1000:             THETA = 4*ATAN(1.0D0) 
1001:          ELSE 
1002:             THETA    = ACOS(CT) 
1003:          ENDIF 
1004:          ST       = SIN(THETA)994:          ST       = SIN(THETA)
1005: 995: 
1006:          IF(ABS(THETA).LT.1.0D-10) THEN ! This happens with ALIGNRBST, for example996:          IF(ABS(THETA).LT.1.0D-10) THEN ! This happens with ALIGNRBST, for example
1007:             DO I = 2, NIMAGE+1997:             DO I = 2, NIMAGE+1
1008:                 J = NOPT*(I-1) + RBOFFSET + 3*(J1-1)998:                 J = NOPT*(I-1) + RBOFFSET + 3*(J1-1)
1009:                 XYZ(J+1:J+3) = P(:)999:                 XYZ(J+1:J+3) = P(:)
1010:             ENDDO1000:             ENDDO
1011:             RETURN1001:             RETURN
1012:          ENDIF1002:          ENDIF
1013: !     INCREMENTAL APPROACH: TANGENT QUATERNION1003: !     INCREMENTAL APPROACH: TANGENT QUATERNION
1031: !     CONVERT FROM QUATERNION -> AA1021: !     CONVERT FROM QUATERNION -> AA
1032:              THETA  = 2.D0*ACOS(QTN(I,1))1022:              THETA  = 2.D0*ACOS(QTN(I,1))
1033: 1023: 
1034:              J = NOPT*(I-1) + RBOFFSET + 3*(J1-1)1024:              J = NOPT*(I-1) + RBOFFSET + 3*(J1-1)
1035:              IF (ABS(THETA) .LT. 1.0D-8) THEN1025:              IF (ABS(THETA) .LT. 1.0D-8) THEN
1036:                  XYZ(J+1:J+3) = 0.D01026:                  XYZ(J+1:J+3) = 0.D0
1037:              ELSE1027:              ELSE
1038:                  FCT    = DSQRT(DOT_PRODUCT(QTN(I,2:4),QTN(I,2:4)))1028:                  FCT    = DSQRT(DOT_PRODUCT(QTN(I,2:4),QTN(I,2:4)))
1039:                  XYZ(J+1:J+3) = THETA*QTN(I,2:4)/FCT1029:                  XYZ(J+1:J+3) = THETA*QTN(I,2:4)/FCT
1040:              ENDIF1030:              ENDIF
 1031: 
1041:          ENDDO1032:          ENDDO
1042: 1033: 
1043:          DEALLOCATE(QTN)1034:          DEALLOCATE(QTN)
1044:          DEALLOCATE(PTN)1035:          DEALLOCATE(PTN)
1045: 1036: 
1046:     END SUBROUTINE iSLERP1037:     END SUBROUTINE iSLERP
1047: ! ------------------------------------------------------------------------------------------------------1038: ! ------------------------------------------------------------------------------------------------------
1048:     ! Extract the index and energy of the highest-energy image in an array of images (BAND)1039:     ! Extract the index and energy of the highest-energy image in an array of images (BAND)
1049:     ! Normally, BAND would be XYZ from NEBDATA, but we sometimes want to find the highest image in a reduced band1040:     ! Normally, BAND would be XYZ from NEBDATA, but we sometimes want to find the highest image in a reduced band
1050:     ! (either using a subset of the images or a subset of the system in each image)1041:     ! (either using a subset of the images or a subset of the system in each image)
1056:     DOUBLE PRECISION, INTENT(OUT) :: EWORST1047:     DOUBLE PRECISION, INTENT(OUT) :: EWORST
1057:     INTEGER, INTENT(OUT) :: JWORST ! Energy and index of the highest-energy index, respectively1048:     INTEGER, INTENT(OUT) :: JWORST ! Energy and index of the highest-energy index, respectively
1058:     INTEGER :: J11049:     INTEGER :: J1
1059:     DOUBLE PRECISION :: E1, RMSDUMMY, LGDUMMY(NOPT), COORDS(NOPT)1050:     DOUBLE PRECISION :: E1, RMSDUMMY, LGDUMMY(NOPT), COORDS(NOPT)
1060: 1051: 
1061:     EWORST=-HUGE(1.0D0)1052:     EWORST=-HUGE(1.0D0)
1062:     1053:     
1063:     DO J1=1,NIMAGE1054:     DO J1=1,NIMAGE
1064: 1055: 
1065:         COORDS = BAND(NOPT*J1+1:NOPT*(J1+1))1056:         COORDS = BAND(NOPT*J1+1:NOPT*(J1+1))
1066:         write(*,*) "Image coords"1057: 
1067:         write(*,*) COORDS 
1068:         ! Determine the energy of this image1058:         ! Determine the energy of this image
1069:         CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)1059:         CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)
1070: 1060: 
1071:         IF (E1.GT.EWORST) THEN1061:         IF (E1.GT.EWORST) THEN
1072:             JWORST=J11062:             JWORST=J1
1073:             EWORST=E11063:             EWORST=E1
1074:         ENDIF1064:         ENDIF
1075:         IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> image ',J1,' energy=',E11065:         IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> image ',J1,' energy=',E1
1076:     ENDDO1066:     ENDDO
1077:     IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> highest image is ',JWORST,' with energy=',EWORST 1067:     IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> highest image is ',JWORST,' with energy=',EWORST 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0