hdiff output

r31191/genrigid.f90 2016-09-23 11:30:16.970809946 +0100 r31190/genrigid.f90 2016-09-23 11:30:17.974823326 +0100
672:   IF(HAS_LATTICE_COORDS) THEN672:   IF(HAS_LATTICE_COORDS) THEN
673:     XRIGIDCOORDS(DEGFREEDOMS - 5:DEGFREEDOMS) =  XCOORDS(3*NATOMS-5:3*NATOMS)673:     XRIGIDCOORDS(DEGFREEDOMS - 5:DEGFREEDOMS) =  XCOORDS(3*NATOMS-5:3*NATOMS)
674:   ENDIF674:   ENDIF
675: 675: 
676: ! sn402: the remainder of the subroutine is a safety check.676: ! sn402: the remainder of the subroutine is a safety check.
677:   IF(DEBUG) THEN677:   IF(DEBUG) THEN
678:     CALL TRANSFORMRIGIDTOC(1,NRIGIDBODY,TEMPXCOORDS,XRIGIDCOORDS)678:     CALL TRANSFORMRIGIDTOC(1,NRIGIDBODY,TEMPXCOORDS,XRIGIDCOORDS)
679:     DO J1 = 1, 3*NATOMS679:     DO J1 = 1, 3*NATOMS
680:         IF(ABS(TEMPXCOORDS(J1)-SAVEXCOORDS(J1)) .GT. 1.0E-7) THEN680:         IF(ABS(TEMPXCOORDS(J1)-SAVEXCOORDS(J1)) .GT. 1.0E-7) THEN
681:             IF(.NOT.(BULKT) .OR. (ABS(ABS(TEMPXCOORDS(J1)-SAVEXCOORDS(J1))-BULK_BOXVEC(MOD(J1-1,3)+1)) .GT. 1.0E-7)) THEN681:             IF(.NOT.(BULKT) .OR. (ABS(ABS(TEMPXCOORDS(J1)-SAVEXCOORDS(J1))-BULK_BOXVEC(MOD(J1-1,3)+1)) .GT. 1.0E-7)) THEN
682:                 WRITE(*,*) "Warning, coordinate transform may have failed for coordinate ", J1682:                 WRITE(*,*) "Warning, coordinate transform may have failed."
683:                 WRITE(*,*) "Original coordinate:", SAVEXCOORDS(J1)683:                 WRITE(*,*) "Original coordinate:", SAVEXCOORDS(J1)
684:                 WRITE(*,*) "New coordinate:", TEMPXCOORDS(J1)684:                 WRITE(*,*) "New coordinate:", TEMPXCOORDS(J1)
685:                 WRITE(*,*) "Difference:", TEMPXCOORDS(J1)-SAVEXCOORDS(J1)685:                 WRITE(*,*) "Difference:", TEMPXCOORDS(J1)-SAVEXCOORDS(J1)
686:                 FAILCOUNT = FAILCOUNT+1686:                 FAILCOUNT = FAILCOUNT+1
687:             ENDIF687:             ENDIF
688:         ENDIF688:         ENDIF
689: !        IF(BULKT) THEN689: !        IF(BULKT) THEN
690: !            IF(ABS(TEMPXCOORDS(J1)-SAVEXCOORDS(J1)-BULK_BOXVEC(MOD(J1-1,3)+1)) .GT. 1.0E-7) THEN690: !            IF(ABS(TEMPXCOORDS(J1)-SAVEXCOORDS(J1)-BULK_BOXVEC(MOD(J1-1,3)+1)) .GT. 1.0E-7) THEN
691: !                WRITE(*,*) "Warning, coordinate transform may have failed."691: !                WRITE(*,*) "Warning, coordinate transform may have failed."
692: !                WRITE(*,*) "Original coordinate:", SAVEXCOORDS(J1)692: !                WRITE(*,*) "Original coordinate:", SAVEXCOORDS(J1)
734:       NLATTICECOORDS=6734:       NLATTICECOORDS=6
735:   ENDIF735:   ENDIF
736: 736: 
737:   GTEST = .TRUE.737:   GTEST = .TRUE.
738:   GR(:) = 0.0D0738:   GR(:) = 0.0D0
739:   739:   
740:   DO J1 = 1, NRIGIDBODY740:   DO J1 = 1, NRIGIDBODY
741:      741:      
742:      PI = XR(3*NRIGIDBODY+3*J1-2 : 3*NRIGIDBODY+3*J1)742:      PI = XR(3*NRIGIDBODY+3*J1-2 : 3*NRIGIDBODY+3*J1)
743:      CALL RMDRVT(PI, RMI, DRMI1, DRMI2, DRMI3, GTEST)743:      CALL RMDRVT(PI, RMI, DRMI1, DRMI2, DRMI3, GTEST)
744: 744:      
745:      DO J2 = 1, NSITEPERBODY(J1)745:      DO J2 = 1, NSITEPERBODY(J1)
746:         J9 = RIGIDGROUPS(J2, J1)746:         J9 = RIGIDGROUPS(J2, J1)
747: 747: 
748: ! hk286 > translation748: ! hk286 > translation
749:         GR(3*J1-2:3*J1) = GR(3*J1-2:3*J1) + G(3*J9-2:3*J9)749:         GR(3*J1-2:3*J1) = GR(3*J1-2:3*J1) + G(3*J9-2:3*J9)
750:         750:         
751: ! hk286 > rotation751: ! hk286 > rotation
752:         DR1(:) = MATMUL(DRMI1,SITESRIGIDBODY(J2,:,J1))752:         DR1(:) = MATMUL(DRMI1,SITESRIGIDBODY(J2,:,J1))
753:         DR2(:) = MATMUL(DRMI2,SITESRIGIDBODY(J2,:,J1))753:         DR2(:) = MATMUL(DRMI2,SITESRIGIDBODY(J2,:,J1))
754:         DR3(:) = MATMUL(DRMI3,SITESRIGIDBODY(J2,:,J1))754:         DR3(:) = MATMUL(DRMI3,SITESRIGIDBODY(J2,:,J1))
1012:   NOPT = 3*NATOMS1012:   NOPT = 3*NATOMS
1013:   DO I=1,NIMAGE+21013:   DO I=1,NIMAGE+2
1014:      XRIGIDCOORDS(1:DEGFREEDOMS) = XYZ( NOPT*(I-1)+1:NOPT*(I-1)+DEGFREEDOMS )1014:      XRIGIDCOORDS(1:DEGFREEDOMS) = XYZ( NOPT*(I-1)+1:NOPT*(I-1)+DEGFREEDOMS )
1015:      CALL TRANSFORMRIGIDTOC (1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)1015:      CALL TRANSFORMRIGIDTOC (1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)
1016:      XYZ(NOPT*(I-1)+1:NOPT*I) = XCOORDS(1:NOPT)1016:      XYZ(NOPT*(I-1)+1:NOPT*I) = XCOORDS(1:NOPT)
1017:   ENDDO1017:   ENDDO
1018: 1018: 
1019: END SUBROUTINE1019: END SUBROUTINE
1020: 1020: 
1021: !-----------------------------------------------------------1021: !-----------------------------------------------------------
1022: ! Calculate the Hessian in rigid body coordinates from the Gradient and Hessian in Cartesians.1022: 
1023: ! Follows the procedure outlined in Mochizuki et al, PCCP 16 (2014) 
1024: SUBROUTINE TRANSFORMHESSIAN (H, G, XR, HR, RBAANORMALMODET)1023: SUBROUTINE TRANSFORMHESSIAN (H, G, XR, HR, RBAANORMALMODET)
1025:   1024:   
1026:   USE COMMONS, ONLY: NATOMS, DEBUG1025:   USE COMMONS, ONLY: NATOMS, DEBUG
1027:   IMPLICIT NONE1026:   IMPLICIT NONE
1028:   1027:   
1029:   INTEGER          :: J1, J2, J3, J4, J8, J9, K, L1028:   INTEGER          :: J1, J2, J3, J4, J8, J9, K, L
1030:   DOUBLE PRECISION, INTENT(IN) :: G(3*NATOMS), H(3*NATOMS,3*NATOMS), XR(DEGFREEDOMS)1029:   DOUBLE PRECISION, INTENT(IN) :: G(3*NATOMS), H(3*NATOMS,3*NATOMS), XR(DEGFREEDOMS)
1031:   DOUBLE PRECISION, INTENT(OUT) :: HR(DEGFREEDOMS,DEGFREEDOMS)1030:   DOUBLE PRECISION, INTENT(OUT) :: HR(DEGFREEDOMS,DEGFREEDOMS)
1032:   DOUBLE PRECISION :: PI(3)1031:   DOUBLE PRECISION :: PI(3)
1033:   DOUBLE PRECISION :: AD2R11(3),AD2R22(3),AD2R33(3),AD2R12(3),AD2R23(3),AD2R31(3) 1032:   DOUBLE PRECISION :: AD2R11(3),AD2R22(3),AD2R33(3),AD2R12(3),AD2R23(3),AD2R31(3) 
1040:   DOUBLE PRECISION :: BD2RMI11(3,3), BD2RMI22(3,3), BD2RMI33(3,3)1039:   DOUBLE PRECISION :: BD2RMI11(3,3), BD2RMI22(3,3), BD2RMI33(3,3)
1041:   DOUBLE PRECISION :: BD2RMI12(3,3), BD2RMI23(3,3), BD2RMI31(3,3)1040:   DOUBLE PRECISION :: BD2RMI12(3,3), BD2RMI23(3,3), BD2RMI31(3,3)
1042:   LOGICAL :: GTEST, STEST, RBAANORMALMODET1041:   LOGICAL :: GTEST, STEST, RBAANORMALMODET
1043:   DOUBLE PRECISION :: RMI0(3,3), DRMI10(3,3), DRMI20(3,3), DRMI30(3,3)1042:   DOUBLE PRECISION :: RMI0(3,3), DRMI10(3,3), DRMI20(3,3), DRMI30(3,3)
1044:   DOUBLE PRECISION :: D2RMI10(3,3), D2RMI20(3,3), D2RMI30(3,3), D2RMI120(3,3), D2RMI230(3,3), D2RMI310(3,3)1043:   DOUBLE PRECISION :: D2RMI10(3,3), D2RMI20(3,3), D2RMI30(3,3), D2RMI120(3,3), D2RMI230(3,3), D2RMI310(3,3)
1045:   1044:   
1046:   GTEST = .TRUE.1045:   GTEST = .TRUE.
1047:   STEST = .TRUE.1046:   STEST = .TRUE.
1048:   HR(:,:) = 0.0D01047:   HR(:,:) = 0.0D0
1049: 1048: 
1050:   IF(DEBUG) THEN 
1051:      IF (.NOT.(ANY(ABS(G(DEGFREEDOMS+1:3*NATOMS)) .GT. 1.0E-8))) THEN 
1052:            WRITE(*,*) "genrigid> Error: Gradient passed into TRANSFORM_HESSIAN seems to be in rigid body coordinates." 
1053:            WRITE(*,*) G(DEGFREEDOMS+1:) 
1054:            STOP 
1055:      ENDIF 
1056:   ENDIF 
1057:  
1058:   IF ( RBAANORMALMODET ) THEN1049:   IF ( RBAANORMALMODET ) THEN
1059:      PI = (/0.0D0, 0.0D0, 0.0D0/)1050:      PI = (/0.0D0, 0.0D0, 0.0D0/)
1060:      CALL RMDFAS(PI, RMI0, DRMI10, DRMI20, DRMI30, D2RMI10, D2RMI20, D2RMI30, D2RMI120, D2RMI230, D2RMI310, GTEST, STEST)1051:      CALL RMDFAS(PI, RMI0, DRMI10, DRMI20, DRMI30, D2RMI10, D2RMI20, D2RMI30, D2RMI120, D2RMI230, D2RMI310, GTEST, STEST)
1061:   ENDIF1052:   ENDIF
1062: 1053: 
1063:   DO J1 = 1, NRIGIDBODY1054:   DO J1 = 1, NRIGIDBODY
1064: 1055: 
1065:      PI = XR(3*NRIGIDBODY+3*J1-2 : 3*NRIGIDBODY+3*J1)1056:      PI = XR(3*NRIGIDBODY+3*J1-2 : 3*NRIGIDBODY+3*J1)
1066:      CALL RMDFAS(PI, ARMI, ADRMI1, ADRMI2, ADRMI3, AD2RMI11, AD2RMI22, AD2RMI33, AD2RMI12, AD2RMI23, AD2RMI31, GTEST, STEST)1057:      CALL RMDFAS(PI, ARMI, ADRMI1, ADRMI2, ADRMI3, AD2RMI11, AD2RMI22, AD2RMI33, AD2RMI12, AD2RMI23, AD2RMI31, GTEST, STEST)
1067:     1058:     
1068:      DO J2 = J1, NRIGIDBODY1059:      DO J2 = J1, NRIGIDBODY
1069: 1060:         
1070:         PI = XR(3*NRIGIDBODY+3*J2-2 : 3*NRIGIDBODY+3*J2)1061:         PI = XR(3*NRIGIDBODY+3*J2-2 : 3*NRIGIDBODY+3*J2)
1071:         CALL RMDFAS(PI, BRMI, BDRMI1, BDRMI2, BDRMI3, BD2RMI11, BD2RMI22, BD2RMI33, BD2RMI12, BD2RMI23, BD2RMI31, GTEST, STEST)1062:         CALL RMDFAS(PI, BRMI, BDRMI1, BDRMI2, BDRMI3, BD2RMI11, BD2RMI22, BD2RMI33, BD2RMI12, BD2RMI23, BD2RMI31, GTEST, STEST)
1072:              1063:         
1073:         DO J3 = 1, NSITEPERBODY(J1)1064:         DO J3 = 1, NSITEPERBODY(J1)
1074:            J8 = RIGIDGROUPS(J3, J1)1065:            J8 = RIGIDGROUPS(J3, J1)
1075: 1066: 
1076:            DO J4 = 1, NSITEPERBODY(J2)1067:            DO J4 = 1, NSITEPERBODY(J2)
1077:               J9 = RIGIDGROUPS(J4, J2)1068:               J9 = RIGIDGROUPS(J4, J2)
1078: 1069:                             
1079: ! hk286 > translation1070: ! hk286 > translation
1080:               HR(3*J1-2:3*J1, 3*J2-2:3*J2) = HR(3*J1-2:3*J1, 3*J2-2:3*J2) + H(3*J8-2:3*J8, 3*J9-2:3*J9)1071:               HR(3*J1-2:3*J1, 3*J2-2:3*J2) = HR(3*J1-2:3*J1, 3*J2-2:3*J2) + H(3*J8-2:3*J8, 3*J9-2:3*J9)
1081: 1072: 
1082: ! hk286 > rotations1073: ! hk286 > rotations
1083:               IF ( RBAANORMALMODET ) THEN1074:               IF ( RBAANORMALMODET ) THEN
1084:                  ADR1(:) = MATMUL(DRMI10,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1075:                  ADR1(:) = MATMUL(DRMI10,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1085:                  ADR2(:) = MATMUL(DRMI20,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1076:                  ADR2(:) = MATMUL(DRMI20,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1086:                  ADR3(:) = MATMUL(DRMI30,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1077:                  ADR3(:) = MATMUL(DRMI30,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1087:                  BDR1(:) = MATMUL(DRMI10,MATMUL(BRMI,SITESRIGIDBODY(J4,:,J2)))1078:                  BDR1(:) = MATMUL(DRMI10,MATMUL(BRMI,SITESRIGIDBODY(J4,:,J2)))
1088:                  BDR2(:) = MATMUL(DRMI20,MATMUL(BRMI,SITESRIGIDBODY(J4,:,J2)))1079:                  BDR2(:) = MATMUL(DRMI20,MATMUL(BRMI,SITESRIGIDBODY(J4,:,J2)))
1143:                     HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-1)=HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-1)+&1134:                     HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-1)=HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-1)+&
1144:                     & H(3*J8-3+K, 3*J9-3+L) * ADR3(K) * BDR2(L)1135:                     & H(3*J8-3+K, 3*J9-3+L) * ADR3(K) * BDR2(L)
1145:                     HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2  )=HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2  )+&1136:                     HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2  )=HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2  )+&
1146:                     & H(3*J8-3+K, 3*J9-3+L) * ADR1(K) * BDR3(L)1137:                     & H(3*J8-3+K, 3*J9-3+L) * ADR1(K) * BDR3(L)
1147:                     HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  )=HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  )+&1138:                     HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  )=HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  )+&
1148:                     & H(3*J8-3+K, 3*J9-3+L) * ADR2(K) * BDR3(L)1139:                     & H(3*J8-3+K, 3*J9-3+L) * ADR2(K) * BDR3(L)
1149:                     HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  )=HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  )+&1140:                     HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  )=HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  )+&
1150:                     & H(3*J8-3+K, 3*J9-3+L) * ADR3(K) * BDR3(L)   1141:                     & H(3*J8-3+K, 3*J9-3+L) * ADR3(K) * BDR3(L)   
1151:                  ENDDO1142:                  ENDDO
1152:               ENDDO1143:               ENDDO
1153:            ENDDO ! Loop over J41144:            ENDDO
1154:  
1155:            IF (J1 .EQ. J2) THEN1145:            IF (J1 .EQ. J2) THEN
1156:               IF ( RBAANORMALMODET ) THEN1146:               IF ( RBAANORMALMODET ) THEN
1157:                  AD2R11(:) = MATMUL(D2RMI10, MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1147:                  AD2R11(:) = MATMUL(D2RMI10, MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1158:                  AD2R22(:) = MATMUL(D2RMI20, MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1148:                  AD2R22(:) = MATMUL(D2RMI20, MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1159:                  AD2R33(:) = MATMUL(D2RMI30, MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1149:                  AD2R33(:) = MATMUL(D2RMI30, MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1160:                  AD2R12(:) = MATMUL(D2RMI120,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1150:                  AD2R12(:) = MATMUL(D2RMI120,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1161:                  AD2R23(:) = MATMUL(D2RMI230,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1151:                  AD2R23(:) = MATMUL(D2RMI230,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1162:                  AD2R31(:) = MATMUL(D2RMI310,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))1152:                  AD2R31(:) = MATMUL(D2RMI310,MATMUL(ARMI,SITESRIGIDBODY(J3,:,J1)))
1163:               ELSE1153:               ELSE
1164:                  AD2R11(:) = MATMUL(AD2RMI11,SITESRIGIDBODY(J3,:,J1))1154:                  AD2R11(:) = MATMUL(AD2RMI11,SITESRIGIDBODY(J3,:,J1))
1165:                  AD2R22(:) = MATMUL(AD2RMI22,SITESRIGIDBODY(J3,:,J1))1155:                  AD2R22(:) = MATMUL(AD2RMI22,SITESRIGIDBODY(J3,:,J1))
1166:                  AD2R33(:) = MATMUL(AD2RMI33,SITESRIGIDBODY(J3,:,J1))1156:                  AD2R33(:) = MATMUL(AD2RMI33,SITESRIGIDBODY(J3,:,J1))
1167:                  AD2R12(:) = MATMUL(AD2RMI12,SITESRIGIDBODY(J3,:,J1))1157:                  AD2R12(:) = MATMUL(AD2RMI12,SITESRIGIDBODY(J3,:,J1))
1168:                  AD2R23(:) = MATMUL(AD2RMI23,SITESRIGIDBODY(J3,:,J1))1158:                  AD2R23(:) = MATMUL(AD2RMI23,SITESRIGIDBODY(J3,:,J1))
1169:                  AD2R31(:) = MATMUL(AD2RMI31,SITESRIGIDBODY(J3,:,J1))1159:                  AD2R31(:) = MATMUL(AD2RMI31,SITESRIGIDBODY(J3,:,J1))
1170:               ENDIF1160:               ENDIF
1171:  
1172:               ! p_x, p_x 
1173:               HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-2) = HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-2) &1161:               HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-2) = HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-2) &
1174:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R11(:))1162:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R11(:))
1175:               ! p_y, p_y 
1176:               HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2-1) = HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2-1) &1163:               HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2-1) = HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2-1) &
1177:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R22(:))1164:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R22(:))
1178:               ! p_z, p_z 
1179:               HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  ) = HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  ) &1165:               HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  ) = HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2  ) &
1180:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R33(:))1166:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R33(:))
1181:               ! p_x, p_y 
1182:               HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-1) = HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-1) &1167:               HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-1) = HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J2-1) &
1183:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R12(:))1168:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R12(:))
1184:               ! p_y, p_z 
1185:               HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  ) = HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  ) &1169:               HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  ) = HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J2  ) &
1186:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R23(:))1170:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R23(:))
1187:               ! p_z, p_x 
1188:               HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-2) = HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-2) &1171:               HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-2) = HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J2-2) &
1189:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R31(:))1172:                    + DOT_PRODUCT(G(3*J8-2:3*J8),AD2R31(:))
1190:            ENDIF1173:            ENDIF
1191:         ENDDO ! Loop over J3 (i.e. J8)1174:         ENDDO
1192: 1175: 
1193:         IF (J1 .EQ. J2) THEN1176:         IF (J1 .EQ. J2) THEN
1194:            HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J1-2) = HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J1-1) 1177:            HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J1-2) = HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J1-1) 
1195:            HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J1-1) = HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J1  )1178:            HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J1-1) = HR(3*NRIGIDBODY+3*J1-1, 3*NRIGIDBODY+3*J1  )
1196:            HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J1  ) = HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J1-2)1179:            HR(3*NRIGIDBODY+3*J1-2, 3*NRIGIDBODY+3*J1  ) = HR(3*NRIGIDBODY+3*J1  , 3*NRIGIDBODY+3*J1-2)
1197:            HR(3*NRIGIDBODY+3*J1-2:3*NRIGIDBODY+3*J1, 3*J1-2:3*J1) = &1180:            HR(3*NRIGIDBODY+3*J1-2:3*NRIGIDBODY+3*J1, 3*J1-2:3*J1) = &
1198:                 TRANSPOSE(HR(3*J1-2:3*J1, 3*NRIGIDBODY+3*J1-2:3*NRIGIDBODY+3*J1))1181:                 TRANSPOSE(HR(3*J1-2:3*J1, 3*NRIGIDBODY+3*J1-2:3*NRIGIDBODY+3*J1))
1199:         ELSE1182:         ELSE
1200:            HR(3*J2-2:3*J2, 3*J1-2:3*J1) = TRANSPOSE(HR(3*J1-2:3*J1, 3*J2-2:3*J2))1183:            HR(3*J2-2:3*J2, 3*J1-2:3*J1) = TRANSPOSE(HR(3*J1-2:3*J1, 3*J2-2:3*J2))
1201:            HR(3*NRIGIDBODY+3*J2-2:3*NRIGIDBODY+3*J2, 3*J1-2:3*J1) = &1184:            HR(3*NRIGIDBODY+3*J2-2:3*NRIGIDBODY+3*J2, 3*J1-2:3*J1) = &
1381:   USE MODHESS     ! For the Hessian1364:   USE MODHESS     ! For the Hessian
1382:   USE MODCHARMM1365:   USE MODCHARMM
1383:   USE KEY1366:   USE KEY
1384:   IMPLICIT NONE1367:   IMPLICIT NONE
1385: 1368: 
1386:   ! Communicating with the rest of OPTIM1369:   ! Communicating with the rest of OPTIM
1387:   DOUBLE PRECISION, INTENT(IN)  :: X(3*NATOMS), ATOMMASS(3*NATOMS) ! Input coordinates and atom masses1370:   DOUBLE PRECISION, INTENT(IN)  :: X(3*NATOMS), ATOMMASS(3*NATOMS) ! Input coordinates and atom masses
1388:                                                 ! Note, ATOMMASS is currently not being used.1371:                                                 ! Note, ATOMMASS is currently not being used.
1389:   DOUBLE PRECISION, INTENT(OUT) :: DIAG(3*NATOMS) ! Output: will contain the normal mode frequencies1372:   DOUBLE PRECISION, INTENT(OUT) :: DIAG(3*NATOMS) ! Output: will contain the normal mode frequencies
1390:   LOGICAL                       :: ART  ! Temporary variable to save the value of ATOMRIGIDCOORDT1373:   LOGICAL                       :: ART  ! Temporary variable to save the value of ATOMRIGIDCOORDT
1391:   LOGICAL                       :: RBT  ! Temporary variable to save the value of RBAANORMALMODET 
1392: 1374: 
1393:   ! Constructing the metric tensor1375:   ! Constructing the metric tensor
1394:   DOUBLE PRECISION :: P(3) ! Angle-axis vector of a particular rigid body1376:   DOUBLE PRECISION :: P(3) ! Angle-axis vector of a particular rigid body
1395:   DOUBLE PRECISION :: W ! Total weight of a particular rigid body1377:   DOUBLE PRECISION :: W ! Total weight of a particular rigid body
1396:   DOUBLE PRECISION :: COM(3) ! Centre of mass for a particular rigid body1378:   DOUBLE PRECISION :: COM(3) ! Centre of mass for a particular rigid body
1397:   DOUBLE PRECISION :: R(3,3), DR1(3,3), DR2(3,3), DR3(3,3) ! Rotation matrix and its derivatives for an RB1379:   DOUBLE PRECISION :: R(3,3), DR1(3,3), DR2(3,3), DR3(3,3) ! Rotation matrix and its derivatives for an RB
1398:   DOUBLE PRECISION :: RB_IMT(6,6)  ! Inverse Metric Tensor (IMT) for a particular rigid body1380:   DOUBLE PRECISION :: RB_IMT(6,6)  ! Inverse Metric Tensor (IMT) for a particular rigid body
1399:   DOUBLE PRECISION :: LR(3,3) ! A temporary matrix used to invert the rotational component of the metric tensor1381:   DOUBLE PRECISION :: LR(3,3) ! A temporary matrix used to invert the rotational component of the metric tensor
1400:   DOUBLE PRECISION :: DET ! Determinant of the rotational component of the metric tensor1382:   DOUBLE PRECISION :: DET ! Determinant of the rotational component of the metric tensor
1401:   DOUBLE PRECISION :: INVERSE_METRIC_TENSOR(DEGFREEDOMS,DEGFREEDOMS) ! For the whole system1383:   DOUBLE PRECISION :: INVERSE_METRIC_TENSOR(DEGFREEDOMS,DEGFREEDOMS) ! For the whole system
1417: 1399: 
1418: ! The following required to call the LAPACK routines DGETRF, DGETRI and DGEEV1400: ! The following required to call the LAPACK routines DGETRF, DGETRI and DGEEV
1419:   INTEGER, INTENT(OUT) :: INFO  ! An exit code for the LAPACK functions1401:   INTEGER, INTENT(OUT) :: INFO  ! An exit code for the LAPACK functions
1420:   INTEGER, PARAMETER   :: LWORK = 1000000 ! Arbitrary dimension for WORK1402:   INTEGER, PARAMETER   :: LWORK = 1000000 ! Arbitrary dimension for WORK
1421:   DOUBLE PRECISION     :: WORK(LWORK) ! Internal arrays for the LAPACK functions1403:   DOUBLE PRECISION     :: WORK(LWORK) ! Internal arrays for the LAPACK functions
1422:   INTEGER              :: PIVOTS(3)1404:   INTEGER              :: PIVOTS(3)
1423: 1405: 
1424: ! Initialize variables1406: ! Initialize variables
1425:   DIAG(:) = 0.0D01407:   DIAG(:) = 0.0D0
1426:   ART = ATOMRIGIDCOORDT1408:   ART = ATOMRIGIDCOORDT
1427:   RBT = RBAANORMALMODET 
1428: 1409: 
1429:   INVERSE_METRIC_TENSOR(:,:) = 0.0D01410:   INVERSE_METRIC_TENSOR(:,:) = 0.0D0
1430:   EIG_REAL(:) = 0.0D01411:   EIG_REAL(:) = 0.0D0
1431:   EIG_IM(:) = 0.0D01412:   EIG_IM(:) = 0.0D0
1432:   DUMMY_EVECS(:,:) = 0.0D01413:   DUMMY_EVECS(:,:) = 0.0D0
1433: 1414: 
1434:   DUMMY = 01415:   DUMMY = 0
1435:   NFAILS = 01416:   NFAILS = 0
1436: 1417: 
1437: ! Ensure we are in rigid body coordinates1418: ! Ensure we are in rigid body coordinates
1568:     DO J3=1,31549:     DO J3=1,3
1569:         write(*,*) RB_IMT(J3,1:3)1550:         write(*,*) RB_IMT(J3,1:3)
1570:     ENDDO1551:     ENDDO
1571:     write(*,*) "LR:"1552:     write(*,*) "LR:"
1572:     DO J3=4,61553:     DO J3=4,6
1573:         write(*,*) RB_IMT(J3,4:6)1554:         write(*,*) RB_IMT(J3,4:6)
1574:     ENDDO1555:     ENDDO
1575: 1556: 
1576:     ! We finally have the IMT for this particular body!1557:     ! We finally have the IMT for this particular body!
1577:     ! Now put it into the whole system IMT (which is block-diagonal)1558:     ! Now put it into the whole system IMT (which is block-diagonal)
1578:     ! Translational block first1559:     INVERSE_METRIC_TENSOR(6*(J1-1)+1:6*J1,6*(J1-1)+1:6*J1) = RB_IMT(:,:)
1579:     INVERSE_METRIC_TENSOR(3*(J1-1)+1:3*J1,3*(J1-1)+1:3*J1) = RB_IMT(1:3,1:3) 
1580:     ! Then rotational block 
1581:     INVERSE_METRIC_TENSOR(3*NRIGIDBODY+3*(J1-1)+1:3*NRIGIDBODY+3*J1,3*NRIGIDBODY+3*(J1-1)+1:3*NRIGIDBODY+3*J1) = RB_IMT(4:6,4:6) 
1582: 1560: 
1583:     ! Update the number of atoms which have already been considered1561:     ! Update the number of atoms which have already been considered
1584:     DUMMY = DUMMY + NSITEPERBODY(J1)1562:     DUMMY = DUMMY + NSITEPERBODY(J1)
1585:   ENDDO ! End loop over rigid bodies1563:   ENDDO ! End loop over rigid bodies
1586: 1564: 
1587:   ! If we had any serious errors detected during the construction of the inverse metric tensor, stop now.1565:   ! If we had any serious errors detected during the construction of the inverse metric tensor, stop now.
1588:   IF(DEBUG .AND. NFAILS .GT. 0) THEN1566:   IF(DEBUG .AND. NFAILS .GT. 0) THEN
1589:     WRITE(*,*) "GENRIGID_NORMALMODES> Construction of the IMT failed ", NFAILS, "times"1567:     WRITE(*,*) "GENRIGID_NORMALMODES> Construction of the IMT failed ", NFAILS, "times"
1590:     STOP1568:     STOP
1591:   ENDIF1569:   ENDIF
1592: 1570: 
1593:   ! Now we need to compute the correct rigid body hessian1571:   ! Now we need to compute the correct rigid body hessian
1594:   RBAANORMALMODET = .FALSE.1572:   RBAANORMALMODET = .TRUE. ! This is used in TRANSFORMHESSIAN and similar. I need to work out exactly how.
1595:   XCOORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)1573:   XCOORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)
1596:   XCOORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D01574:   XCOORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D0
1597:   IF (ENDNUMHESS .OR. AMBERT .OR. AMBER12T) THEN1575:   IF (ENDNUMHESS .OR. AMBERT .OR. AMBER12T) THEN
1598:      CALL GENRIGID_MAKENUMHESS(XCOORDS,NATOMS,DEGFREEDOMS)1576:      CALL GENRIGID_MAKENUMHESS(XCOORDS,NATOMS,DEGFREEDOMS)
1599:   ELSE1577:   ELSE
 1578:     write(*,*) "Calling potential in GENRIGID_NORMALMODES"
1600:      CALL POTENTIAL(XCOORDS,ENERGY,G,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)1579:      CALL POTENTIAL(XCOORDS,ENERGY,G,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
1601:   ENDIF1580:   ENDIF
 1581:   RBAANORMALMODET = .FALSE.
1602: 1582: 
1603:   ! We need to left-multiply the Hessian by the inverse MT we calculated earlier, then find its eigenvalues and1583:   ! We need to left-multiply the Hessian by the inverse MT we calculated earlier, then find its eigenvalues and
1604:   ! eigenvectors. See eq.31 in the paper.1584:   ! eigenvectors. See eq.31 in the paper.
1605:   METRICHESS(:,:) = MATMUL(INVERSE_METRIC_TENSOR,HESS(1:DEGFREEDOMS,1:DEGFREEDOMS))1585:   METRICHESS(:,:) = MATMUL(INVERSE_METRIC_TENSOR,HESS(1:DEGFREEDOMS,1:DEGFREEDOMS))
1606: 1586: 
1607:   IF (DUMPV) THEN ! Save right-eigenvectors as well.1587:   IF (DUMPV) THEN ! Save right-eigenvectors as well.
1608:      CALL DGEEV('N','V', DEGFREEDOMS, METRICHESS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, EIG_REAL, EIG_IM, &1588:      CALL DGEEV('N','V', DEGFREEDOMS, METRICHESS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, EIG_REAL, EIG_IM, &
1609:                 & DUMMY_EVECS, DEGFREEDOMS, EVECS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, WORK, LWORK, INFO)1589:                 & DUMMY_EVECS, DEGFREEDOMS, EVECS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, WORK, LWORK, INFO)
1610:   ELSE1590:   ELSE
1611:      CALL DGEEV('N','N', DEGFREEDOMS, METRICHESS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, EIG_REAL, EIG_IM, &1591:      CALL DGEEV('N','N', DEGFREEDOMS, METRICHESS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, EIG_REAL, EIG_IM, &
1612:                 & DUMMY_EVECS, DEGFREEDOMS, EVECS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, WORK, LWORK, INFO)1592:                 & DUMMY_EVECS, DEGFREEDOMS, EVECS(1:DEGFREEDOMS,1:DEGFREEDOMS), DEGFREEDOMS, WORK, LWORK, INFO)
1613:   ENDIF1593:   ENDIF
1614: 1594: 
1615:   IF(DEBUG .AND. ANY(ABS(EIG_IM(:)).GT.1.0e-8)) THEN1595:   IF(DEBUG .AND. ANY(ABS(EIG_IM(:)).GT.1.0e-8)) THEN
1616:     WRITE(*,*) "GENRIGID_NORMALMODES> Error: complex eigenvalues returned from metric-tensor Hessian."1596:     WRITE(*,*) "GENRIGID_NORMALMODES> Error: complex eigenvalues returned from metric-tensor Hessian."
1617:     STOP1597:     STOP
1618:   ENDIF1598:   ENDIF
1619: 1599: 
1620:  CALL EIGENSORT_VAL_ASC(EIG_REAL,METRICHESS,DEGFREEDOMS,DEGFREEDOMS)1600:  CALL EIGENSORT_VAL_ASC(EIG_REAL,METRICHESS,DEGFREEDOMS,DEGFREEDOMS)
1621: 1601: 
1622:   IF (MWVECTORS .AND. (AMBERT .OR. AMBER12T)) THEN1602:   IF (MWVECTORS) THEN
1623:      DO J1 = 1, DEGFREEDOMS1603:      DO J1 = 1, DEGFREEDOMS
1624:         IF (EIG_REAL(J1) > 0.0D0) THEN1604:         IF (EIG_REAL(J1) > 0.0D0) THEN
1625:            EIG_REAL(J1) = SQRT((EIG_REAL(J1)))*108.521605:            EIG_REAL(J1) = SQRT((EIG_REAL(J1)))*108.52
1626:         ELSE1606:         ELSE
1627:            EIG_REAL(J1) = -SQRT((-EIG_REAL(J1)))*108.521607:            EIG_REAL(J1) = -SQRT((-EIG_REAL(J1)))*108.52
1628:         ENDIF1608:         ENDIF
1629:      ENDDO1609:      ENDDO
1630:   ENDIF1610:   ENDIF
1631: 1611: 
1632:   ! Copy the eigenvalues to DIAG1612:   ! Copy the eigenvalues to DIAG
1633:   DIAG(1:DEGFREEDOMS) = EIG_REAL(1:DEGFREEDOMS)1613:   DIAG(1:DEGFREEDOMS) = EIG_REAL(1:DEGFREEDOMS)
1634:   ! Fill in the remaining elements with a very large number so that after sorting, these get ignored.1614:   ! Fill in the remaining elements with a very large number so that after sorting, these get ignored.
1635:   DIAG(DEGFREEDOMS+1:3*NATOMS) = 1.0D101615:   DIAG(DEGFREEDOMS+1:3*NATOMS) = 1.0D10
1636: 1616: 
1637: ! Set the Hessian matrix to its eigenvector matrix1617: ! Set the Hessian matrix to its eigenvector matrix
1638:   HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = EVECS(1:DEGFREEDOMS,1:DEGFREEDOMS)1618:   HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = EVECS(1:DEGFREEDOMS,1:DEGFREEDOMS)
1639: 1619: 
1640:   ! Reset this variable to its original value1620:   ! Reset this variable to its original value
1641:   ATOMRIGIDCOORDT = ART1621:   ATOMRIGIDCOORDT = ART
1642:   RBAANORMALMODET = RBT 
1643: 1622: 
1644:   OPEN(UNIT = 28, FILE = 'LRBNORMALMODES')1623:   OPEN(UNIT = 28, FILE = 'LRBNORMALMODES')
1645:   WRITE(28, *) ENERGY, 3*NATOMS, DEGFREEDOMS1624:   WRITE(28, *) ENERGY, 3*NATOMS, DEGFREEDOMS
1646:   DO J1 = 3, DEGFREEDOMS/31625:   DO J1 = 3, DEGFREEDOMS/3
1647:      WRITE(28, *) DIAG(3*J1-2), DIAG(3*J1-1), DIAG(3*J1)1626:      WRITE(28, *) DIAG(3*J1-2), DIAG(3*J1-1), DIAG(3*J1)
1648:   ENDDO1627:   ENDDO
1649:   CLOSE(UNIT = 28)1628:   CLOSE(UNIT = 28)
1650: 1629: 
1651:   IF (DUMPV) THEN1630:   IF (DUMPV) THEN
1652:      CALL GENRIGID_EIGENVECTORS(METRICHESS, XRIGIDCOORDS, ATOMMASS)1631:      CALL GENRIGID_EIGENVECTORS(METRICHESS, XRIGIDCOORDS, ATOMMASS)
1653:   ENDIF1632:   ENDIF
1654: 1633: 
1655: END SUBROUTINE GENRIGID_NORMALMODES1634: END SUBROUTINE GENRIGID_NORMALMODES
1656: !----------------------------------------------------------------------------------------------1635: !----------------------------------------------------------------------------------------------
1657: ! See section 2.2.1 of Mochizuki et al, PCCP 16 (2014), 2842-53  DOI 10.1039/C3CP53537A1636: 
1658: ! However, the notation used is closer to Wales and Ohmine, JCP 98 (1993), 7257-7268 (KBLOCK, U) 
1659: ! See also rigidb.f90>NRMLMD, an obsolete version of the same procedure which has comments from the 
1660: ! original authors. The comments on this routine were mostly added by sn402. 
1661: SUBROUTINE GENRIGID_EIGENVALUES(X, ATOMMASS, DIAG, INFO)1637: SUBROUTINE GENRIGID_EIGENVALUES(X, ATOMMASS, DIAG, INFO)
1662:   1638:   
1663:   USE COMMONS1639:   USE COMMONS
1664:   USE MODHESS1640:   USE MODHESS
1665:   USE MODCHARMM1641:   USE MODCHARMM
1666:   USE KEY1642:   USE KEY
1667:   IMPLICIT NONE1643:   IMPLICIT NONE
1668: 1644: 
1669:   INTEGER          :: IR, IC, OFFSET, NDIM1645:   INTEGER          :: IR, IC, OFFSET, NDIM
1670:   INTEGER          :: I, J, J1, J2, J3, J5, J8, K1, K2, ISTART, IFINISH, JSTART, JFINISH 1646:   INTEGER          :: I, J, J1, J2, J3, J5, J8, K1, K2, ISTART, IFINISH, JSTART, JFINISH 
1672:   DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(3*NATOMS), G(3*NATOMS)1648:   DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(3*NATOMS), G(3*NATOMS)
1673:   DOUBLE PRECISION :: X(3*NATOMS), FRQN(DEGFREEDOMS), ATOMMASS(3*NATOMS), DIAG(3*NATOMS)1649:   DOUBLE PRECISION :: X(3*NATOMS), FRQN(DEGFREEDOMS), ATOMMASS(3*NATOMS), DIAG(3*NATOMS)
1674:   DOUBLE PRECISION :: P(3), RMI(3,3), DRMI(3,3), DR(3)1650:   DOUBLE PRECISION :: P(3), RMI(3,3), DRMI(3,3), DR(3)
1675:   DOUBLE PRECISION :: KD(DEGFREEDOMS), AP(DEGFREEDOMS,DEGFREEDOMS), RMS1651:   DOUBLE PRECISION :: KD(DEGFREEDOMS), AP(DEGFREEDOMS,DEGFREEDOMS), RMS
1676: ! the following required to call the LAPACK routine DSYEV1652: ! the following required to call the LAPACK routine DSYEV
1677:   INTEGER          :: INFO1653:   INTEGER          :: INFO
1678:   INTEGER, PARAMETER :: LWORK = 1000000 ! the dimension is set arbitrarily1654:   INTEGER, PARAMETER :: LWORK = 1000000 ! the dimension is set arbitrarily
1679:   DOUBLE PRECISION :: WORK(LWORK)1655:   DOUBLE PRECISION :: WORK(LWORK)
1680:   LOGICAL          :: ART1656:   LOGICAL          :: ART
1681: 1657: 
1682:   ! The METRICTENSOR keyword instructs us to override all calls to this subroutine with a call 
1683:   ! to the newer subroutine GENRIGID_NORMALMODES. They should have exactly the same behaviour. 
1684:   IF(METRICTENSOR) THEN 
1685:      CALL GENRIGID_NORMALMODES(X, ATOMMASS, DIAG, INFO) 
1686:      RETURN 
1687:   ENDIF 
1688:  
1689: ! Initialize1658: ! Initialize
1690: ! hk2861659: ! hk286
1691:   OFFSET = 3*NRIGIDBODY1660:   OFFSET = 3*NRIGIDBODY
1692:   U(:,:) = 0.D01661:   U(:,:) = 0.D0
1693:   IR     = 0  ! IR+1 is the index of the first coordinate for the current rigid body1662:   IR     = 0
1694:   IC     = 01663:   IC     = 0
1695:   ART = ATOMRIGIDCOORDT1664:   ART = ATOMRIGIDCOORDT
1696: 1665: 
1697: ! Transform coordinates1666: ! Transform coordinates
1698:   IF ( ATOMRIGIDCOORDT .EQV. .TRUE. ) THEN1667:   IF ( ATOMRIGIDCOORDT .EQV. .TRUE. ) THEN
 1668:      write(*,*) "Called with atomistic coordinates."
1699:      CALL TRANSFORMCTORIGID (X, XRIGIDCOORDS)1669:      CALL TRANSFORMCTORIGID (X, XRIGIDCOORDS)
1700:   ELSE1670:   ELSE
 1671:      write(*,*) "Called with RB coordinates"
1701:      XRIGIDCOORDS(1:DEGFREEDOMS) = X(1:DEGFREEDOMS)1672:      XRIGIDCOORDS(1:DEGFREEDOMS) = X(1:DEGFREEDOMS)
1702:   ENDIF1673:   ENDIF
 1674: 
 1675:   DO J1 = 1, NRIGIDBODY
1703:      1676:      
1704:      J3 = 3*J11677:      J3 = 3*J1
1705:      J5 = OFFSET + J31678:      J5 = OFFSET + J3
1706:      P  = XRIGIDCOORDS(J5-2:J5)1679:      P  = XRIGIDCOORDS(J5-2:J5)
1707:      ! KBLOCK is the moment of inertia matrix for this rigid body.1680:      ! KBLOCK is the weighted gyration tensor for this particular rigid body?
1708:      KBLOCK(:,:) = 0.D0     1681:      KBLOCK(:,:) = 0.D0     
1709:      ! Get the rotation matrix RMI that corresponds to P (doesn't calculate any derivatives)1682:      ! Get the rotation matrix RMI that corresponds to P (doesn't calculate any derivatives)
1710:      CALL RMDFAS(P, RMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, .FALSE., .FALSE.)1683:      CALL RMDFAS(P, RMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, DRMI, .FALSE., .FALSE.)
1711: 1684: 
1712:      TMASS = 0.0D0 ! The total mass of the rigid body1685:      TMASS = 0.0D0
1713:      DO J2 = 1, NSITEPERBODY(J1)1686:      DO J2 = 1, NSITEPERBODY(J1)
1714:         J8 = RIGIDGROUPS(J2, J1)1687:         J8 = RIGIDGROUPS(J2, J1)
1715:         ! DR is the coordinates of this atom relative to the rigid body CoM.1688:         ! DR is the coordinates of this atom relative to the molecular CoM.
1716:         DR(:)  = MATMUL(RMI(:,:),SITESRIGIDBODY(J2,:,J1))1689:         DR(:)  = MATMUL(RMI(:,:),SITESRIGIDBODY(J2,:,J1))
1717:         TMASS = TMASS + ATOMMASS(J8)1690:         TMASS = TMASS + ATOMMASS(J8)
1718:         DO I = 1, 31691:         DO I = 1, 3
1719:            ! Diagonal terms in the moment of inertia about the rigid-body CoM 
1720:            ! sn402: Is this right? Shouldn't Ixx = sum_i m_i * (y_i^2+z_i^2) ? Here we have Ixx = sum_i m_i * (x_i^2+y_i^2+z_i^2) 
1721:            KBLOCK(I,I) = KBLOCK(I,I) + ATOMMASS(J8)*(DR(1)*DR(1) + DR(2)*DR(2) + DR(3)*DR(3))1692:            KBLOCK(I,I) = KBLOCK(I,I) + ATOMMASS(J8)*(DR(1)*DR(1) + DR(2)*DR(2) + DR(3)*DR(3))
1722:            DO J = 1, 3    ! could have been J = 1, I; KBLOCK is a symmetric matrix1693:            DO J = 1, 3    ! could have been J = 1, I; KBLOCK is a symmetric matrix
1723:               KBLOCK(I,J) = KBLOCK(I,J) - ATOMMASS(J8)*DR(I)*DR(J)1694:               KBLOCK(I,J) = KBLOCK(I,J) - ATOMMASS(J8)*DR(I)*DR(J)
1724:            ENDDO1695:            ENDDO
1725:         ENDDO1696:         ENDDO
1726:      ENDDO1697:      ENDDO
1727:      ! Diagonalise KBLOCK1698:      ! Diagonalise KBLOCK
1728:      ! KBEGNV are the KBLOCK eigenvalues, KBLOCK now contains the eigenvectors1699:      ! KBEGNV are the KBLOCK eigenvalues, KBLOCK now contains the eigenvectors
1729:      CALL DSYEV('V','L',3,KBLOCK,3,KBEGNV,WORK,LWORK,INFO)1700:      CALL DSYEV('V','L',3,KBLOCK,3,KBEGNV,WORK,LWORK,INFO)
1730:      IF (INFO /= 0) THEN1701:      IF (INFO /= 0) THEN
1731:         WRITE(*,*) 'NRMLMD > Error in DSYEV with KBLOCK, INFO =', INFO1702:         WRITE(*,*) 'NRMLMD > Error in DSYEV with KBLOCK, INFO =', INFO
1732:         STOP1703:         STOP
1733:      ENDIF1704:      ENDIF
1734: 1705: 
1735: !    Construction of the matrix U, which diagonalises the coordinate vector x = (r,p) for this rigid body.1706: !     Construction of the matrix U
1736: !    The reason we want this is that making the transform w = matmul(U,p) gives us coordinates1707: !     First: translation coordinates
1737: !    which are diagonal in the kinetic energy. 
1738:  
1739: !    First: translation coordinates (recall IR+1 is the first coordinate for this rigid body) 
1740: !    No diagonalisation required. U is identity. 
1741:      U(IR+1,IC+1) = 1.D0; U(IR+2,IC+2) = 1.D0; U(IR+3,IC+3) = 1.D01708:      U(IR+1,IC+1) = 1.D0; U(IR+2,IC+2) = 1.D0; U(IR+3,IC+3) = 1.D0
1742:      ! KD is the diagonalised KBLOCK. Only diagonal elements are recorded. 
1743:      ! Scaling by sqrt(TMASS) is required to move to the CoM frame, as usual in a normal mode calculation. 
1744:      KD(IC+1:IC+3) = 1.D0/SQRT(TMASS)1709:      KD(IC+1:IC+3) = 1.D0/SQRT(TMASS)
 1710: !     Now rotational coordinates
 1711:      U(OFFSET+IR+1:OFFSET+IR+3,OFFSET+IC+1:OFFSET+IC+3) = KBLOCK(:,:) 
 1712:      KD(OFFSET+IC+1:OFFSET+IC+3) = 1.D0/SQRT(KBEGNV(:))
1745: 1713: 
1746: !    Now rotational coordinates.1714:      ! IR = IC = 3*(J1-1) ?
1747: !    Recall that KBLOCK now contains the eigenvectors of the original matrix. These are required to diagonalise 
1748: !    the rotational component of the coordinate vector. 
1749:      U(OFFSET+IR+1:OFFSET+IR+3,OFFSET+IC+1:OFFSET+IC+3) = KBLOCK(:,:) 
1750:      KD(OFFSET+IC+1:OFFSET+IC+3) = 1.D0/SQRT(KBEGNV(:))               ! See Mochizuki et al 
1751:  
1752:      ! IR = IC = 3*(J1-1) - these two variables indicate the first coordinate position of the current rigid body. 
1753:      IR = IR + 31715:      IR = IR + 3
1754:      IC = IC + 3 1716:      IC = IC + 3 
1755: 1717: 
1756:   ENDDO ! Loop over rigid bodies (J1)1718:   ENDDO ! Loop over rigid bodies (J1)
1757: 1719: 
1758:   ! Deal with the free atoms. No off-diagonal terms in the kinetic energy here, so they are automatically diagonalised and only 
1759:   ! require mass-rescaling. 
1760:   DO J1 = 1, (DEGFREEDOMS - 6*NRIGIDBODY)/31720:   DO J1 = 1, (DEGFREEDOMS - 6*NRIGIDBODY)/3
1761:      J8 = RIGIDSINGLES(J1)1721:      J8 = RIGIDSINGLES(J1)
1762:      U(6*NRIGIDBODY+3*J1-2,6*NRIGIDBODY+3*J1-2) = 1.D01722:      U(6*NRIGIDBODY+3*J1-2,6*NRIGIDBODY+3*J1-2) = 1.D0
1763:      U(6*NRIGIDBODY+3*J1-1,6*NRIGIDBODY+3*J1-1) = 1.D01723:      U(6*NRIGIDBODY+3*J1-1,6*NRIGIDBODY+3*J1-1) = 1.D0
1764:      U(6*NRIGIDBODY+3*J1  ,6*NRIGIDBODY+3*J1  ) = 1.D01724:      U(6*NRIGIDBODY+3*J1  ,6*NRIGIDBODY+3*J1  ) = 1.D0
1765:      KD(6*NRIGIDBODY+3*J1-2:6*NRIGIDBODY+3*J1 ) = 1.D0/SQRT(ATOMMASS(J8))     1725:      KD(6*NRIGIDBODY+3*J1-2:6*NRIGIDBODY+3*J1 ) = 1.D0/SQRT(ATOMMASS(J8))     
1766:   ENDDO1726:   ENDDO
1767: 1727: 
1768:   ! We have now obtained the required diagonal coordinates for the total kinetic energy. 
1769:   ! We next need to obtain the Hessian in terms of these coordinates. 
1770:  
1771:   RBAANORMALMODET = .TRUE.1728:   RBAANORMALMODET = .TRUE.
1772:   ATOMRIGIDCOORDT = .FALSE.1729:   ATOMRIGIDCOORDT = .FALSE.
1773:   XCOORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)1730:   XCOORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)
1774:   XCOORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D01731:   XCOORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D0
1775:   IF (ENDNUMHESS .OR. AMBERT .OR. AMBER12T) THEN1732:   IF (ENDNUMHESS .OR. AMBERT .OR. AMBER12T) THEN
1776:      CALL GENRIGID_MAKENUMHESS(XCOORDS,NATOMS,DEGFREEDOMS)1733:      CALL GENRIGID_MAKENUMHESS(XCOORDS,NATOMS,DEGFREEDOMS)
1777:   ELSE1734:   ELSE
1778:      ! STEST is set to TRUE so we get the Hessian calculated.1735:     write(*,*) "Calling potential here"
1779:      ! When TRANSFORM_HESSIAN is called from POTENTIAL, it will behave differently because of RBAANORMALMODET=TRUE. 
1780:      ! See TRANSFORM_HESSIAN for comments. 
1781:      CALL POTENTIAL(XCOORDS,ENERGY,G,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)     1736:      CALL POTENTIAL(XCOORDS,ENERGY,G,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)     
1782:   ENDIF1737:   ENDIF
1783:   RBAANORMALMODET = .FALSE.1738:   RBAANORMALMODET = .FALSE.
1784: !  ATOMRIGIDCOORDT = .TRUE.1739: !  ATOMRIGIDCOORDT = .TRUE.
1785: 1740: 
1786:   NDIM = DEGFREEDOMS  1741:   NDIM = DEGFREEDOMS  
1787:   AP(:,:) = 0.D0  ! This will be our Hessian in diagonalised general coordinates.1742:   AP(:,:) = 0.D0
1788:   ! Fill in the lower-left triangle of the transformed Hessian 
1789:   DO I = 1, NDIM1743:   DO I = 1, NDIM
1790:      DO J = 1, I1744:      DO J = 1, I
1791:         1745: 
1792:         ! sn402: I'm not quite sure what's going on here. 
1793:         ! This block seems to suggest that the block AP(1:3,1:3) gets written to more than everything else. 
1794:         ! Wouldn't it be easier to use the upper-right triangle and use 
1795:         !DO I=1,NDIM 
1796:         !   DO J=I,NDIM 
1797:         !      DO K1=1,NDIM 
1798:         !         DO K2=1,NDIM 
1799:         !            ? 
1800:         ! Possibly this way is more efficient because U is 0 away from the main diagonal? 
1801:         IF (I .LE. 2) THEN1746:         IF (I .LE. 2) THEN
1802:            ISTART = 11747:            ISTART = 1
1803:         ELSE1748:         ELSE
1804:            ISTART = I-21749:            ISTART = I-2
1805:         ENDIF1750:         ENDIF
1806:         IF (I .GE. NDIM -2) THEN1751:         IF (I .GE. NDIM -2) THEN
1807:            IFINISH = NDIM1752:            IFINISH = NDIM
1808:         ELSE1753:         ELSE
1809:            IFINISH = I+21754:            IFINISH = I+2
1810:         ENDIF1755:         ENDIF
1814:            JSTART = J-21759:            JSTART = J-2
1815:         ENDIF1760:         ENDIF
1816:         IF (J .GE. NDIM-2) THEN1761:         IF (J .GE. NDIM-2) THEN
1817:            JFINISH = NDIM1762:            JFINISH = NDIM
1818:         ELSE1763:         ELSE
1819:            JFINISH = J+21764:            JFINISH = J+2
1820:         ENDIF1765:         ENDIF
1821: 1766: 
1822:         DO K1 = ISTART, IFINISH1767:         DO K1 = ISTART, IFINISH
1823:            DO K2 = JSTART, JFINISH1768:            DO K2 = JSTART, JFINISH
1824:               ! The actual coordinate transform. Recall that U is diagonal for translational coordinates and single 
1825:               ! atoms (the elements are square roots of the corresponding masses) 
1826:               ! But U and HESS are both different for the rotational coordinates. 
1827:               AP(I,J) = AP(I,J) + U(K1,I)*HESS(K1,K2)*U(K2,J)1769:               AP(I,J) = AP(I,J) + U(K1,I)*HESS(K1,K2)*U(K2,J)
1828:            ENDDO1770:            ENDDO
1829:         ENDDO1771:         ENDDO
1830:         AP(I,J) = KD(I)*AP(I,J)*KD(J) ! Weight the coordinates by the eigenvalues of the inertia tensor (see above)1772:         AP(I,J) = KD(I)*AP(I,J)*KD(J)
1831:      ENDDO1773:      ENDDO
1832:   ENDDO1774:   ENDDO
1833: 1775: 
1834:   ! Now diagonalise the modified Hessian. The eigenvalues (returned as FRQN) are the normal mode frequencies. 
1835:   ! If DUMPV is set, then on return AP contains the eigenvectors of the transformed Hessian - i.e. the normal modes. 
1836:   IF (DUMPV) THEN1776:   IF (DUMPV) THEN
1837:      CALL DSYEV('V','L',NDIM,AP(1:NDIM,1:NDIM),NDIM,FRQN,WORK,LWORK,INFO)1777:      CALL DSYEV('V','L',NDIM,AP(1:NDIM,1:NDIM),NDIM,FRQN,WORK,LWORK,INFO)
1838:   ELSE1778:   ELSE
1839:      CALL DSYEV('N','L',NDIM,AP(1:NDIM,1:NDIM),NDIM,FRQN,WORK,LWORK,INFO)1779:      CALL DSYEV('N','L',NDIM,AP(1:NDIM,1:NDIM),NDIM,FRQN,WORK,LWORK,INFO)
1840:   ENDIF1780:   ENDIF
1841: 1781: 
1842: !  call eigensort_val_asc(FRQN,AP,NDIM,NDIM)1782: !  call eigensort_val_asc(FRQN,AP,NDIM,NDIM)
1843: ! Mass-weight the normal mode frequencies here?1783:   IF (MWVECTORS) THEN
1844: ! Looks like someone has hard-coded the units in. Need to change that. 
1845:   IF (MWVECTORS .AND. (AMBERT .OR. AMBER12T)) THEN ! .OR. CHARMMT)) THEN  ! I need to get this included at some point 
1846:      ! Hessian diagonalisation returns square frequencies in internal units. For AMBER/CHARMM, these units are 
1847:      ! (kCal mol^-1)/(amu Angstrom^2) = 4.184E26 s^-2. 
1848:      ! So to convert to rad s^-1, we want sqrt(FRQN(I)*4.184D26) = 2.045E13*sqrt(FRQN(I)) 
1849:      ! For Hz: sqrt(FRQN(I)*4.184D26)/2*pi = 3.255E12*sqrt(FRQN(I)) 
1850:      ! For cm^-1: sqrt(FRQN(I)*4.184D26)/(2*pi*c*100) = 108.52*sqrt(FRQN(I)) 
1851:      DO I = 1, NDIM1784:      DO I = 1, NDIM
1852:         IF (FRQN(I) > 0.0D0) THEN1785:         IF (FRQN(I) > 0.0D0) THEN
1853:            FRQN(I) = SQRT((FRQN(I)))*108.521786:            FRQN(I) = SQRT((FRQN(I)))*108.52
1854:         ELSE1787:         ELSE
1855:            FRQN(I) = -SQRT((-FRQN(I)))*108.521788:            FRQN(I) = -SQRT((-FRQN(I)))*108.52
1856:         ENDIF1789:         ENDIF
1857:      ENDDO1790:      ENDDO
1858:   ENDIF1791:   ENDIF
1859: 1792: 
1860:   DIAG(1:DEGFREEDOMS) = FRQN1793:   DIAG(1:DEGFREEDOMS) = FRQN
1861:   DIAG(DEGFREEDOMS+1:3*NATOMS) = 1.0D101794:   DIAG(DEGFREEDOMS+1:3*NATOMS) = 1.0D10
1862: 1795: 
1863: ! hk286 - set the Hessian matrix to its eigenvector matrix1796: ! hk286 - set the Hessian matrix to its eigenvector matrix
1864:   HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = AP1797:   HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = AP
1865: ! Restore this variable to its saved value. 
1866:   ATOMRIGIDCOORDT = ART1798:   ATOMRIGIDCOORDT = ART
1867: 1799: 
1868: !  OPEN(UNIT = 28, FILE = 'LRBNORMALMODES')1800: !  OPEN(UNIT = 28, FILE = 'LRBNORMALMODES')
1869: !  WRITE(28, *) ENERGY, 3*NATOMS, DEGFREEDOMS1801: !  WRITE(28, *) ENERGY, 3*NATOMS, DEGFREEDOMS
1870: !  DO J1 = 3, DEGFREEDOMS/31802: !  DO J1 = 3, DEGFREEDOMS/3
1871: !     WRITE(28, *) FRQN(3*J1-2), FRQN(3*J1-1), FRQN(3*J1)1803: !     WRITE(28, *) FRQN(3*J1-2), FRQN(3*J1-1), FRQN(3*J1)
1872: !  ENDDO1804: !  ENDDO
1873: !  CLOSE(UNIT = 28)1805: !  CLOSE(UNIT = 28)
1874: 1806: 
1875:   ! hk2861807:   ! hk286
1965: 1897: 
1966:   ENDDO1898:   ENDDO
1967: 1899: 
1968: END SUBROUTINE GENRIGID_EIGENVECTORS1900: END SUBROUTINE GENRIGID_EIGENVECTORS
1969: 1901: 
1970: !----------------------------------------------------------------------------------------------1902: !----------------------------------------------------------------------------------------------
1971: 1903: 
1972: SUBROUTINE GENRIGID_VDUMP(DIAG,ZT,N,M)1904: SUBROUTINE GENRIGID_VDUMP(DIAG,ZT,N,M)
1973:   USE KEY1905:   USE KEY
1974:   USE MODHESS1906:   USE MODHESS
1975:   USE MODCHARMM, ONLY: CHRMMT 
1976:   USE PORFUNCS1907:   USE PORFUNCS
1977:   IMPLICIT NONE1908:   IMPLICIT NONE
1978:   INTEGER M, N, J1, J2, ISTAT, MCOUNT1909:   INTEGER M, N, J1, J2, ISTAT, MCOUNT
1979:   DOUBLE PRECISION DIAG(M)1910:   DOUBLE PRECISION DIAG(M)
1980:   LOGICAL ZT(M)1911:   LOGICAL ZT(M)
1981: 1912: 
1982: !1913: !
1983: !  dump the eigenvectors which correspond to non-zero eigenvalues1914: !  dump the eigenvectors which correspond to non-zero eigenvalues
1984: !  in file vectors.dump1915: !  in file vectors.dump
1985: !1916: !
1990:         IF (ZT(J1)) MCOUNT=MCOUNT+11921:         IF (ZT(J1)) MCOUNT=MCOUNT+1
1991:      ENDDO1922:      ENDDO
1992:      OPEN(UNIT=499,FILE='nmodes.dat',STATUS='UNKNOWN')1923:      OPEN(UNIT=499,FILE='nmodes.dat',STATUS='UNKNOWN')
1993:      WRITE(499,'(I6)') MCOUNT1924:      WRITE(499,'(I6)') MCOUNT
1994:      CLOSE(499)1925:      CLOSE(499)
1995:      DO J1=1,M1926:      DO J1=1,M
1996:         IF (ZT(J1)) THEN1927:         IF (ZT(J1)) THEN
1997: ! If printing the mass weighted vectors (normal modes), convert omega^21928: ! If printing the mass weighted vectors (normal modes), convert omega^2
1998: ! into the vibrational frequency in wavenumbers (cm^(-1)). 108.52 is the1929: ! into the vibrational frequency in wavenumbers (cm^(-1)). 108.52 is the
1999: ! conversion factor from (kcal mol-1 kg-1 A-2)^2 to cm-1.1930: ! conversion factor from (kcal mol-1 kg-1 A-2)^2 to cm-1.
2000:            IF (MWVECTORS .AND. (AMBERT .OR. AMBER12T .OR. CHRMMT)) THEN1931:            IF (MWVECTORS) THEN
2001:               WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.521932:               WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.52
2002:            ELSE1933:            ELSE
2003:               WRITE(44,'(F20.10)') DIAG(J1)1934:               WRITE(44,'(F20.10)') DIAG(J1)
2004:            ENDIF1935:            ENDIF
2005:            WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N)1936:            WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N)
2006:         ENDIF1937:         ENDIF
2007:      ENDDO1938:      ENDDO
2008:   ELSE1939:   ELSE
2009:      DO J1=M,1,-11940:      DO J1=M,1,-1
2010:         IF (ZT(J1)) THEN1941:         IF (ZT(J1)) THEN
2011: ! As above1942: ! As above
2012:            IF (MWVECTORS .AND. (AMBERT .OR. AMBER12T .OR. CHRMMT)) THEN1943:            IF (MWVECTORS) THEN
2013:               WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.521944:               WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.52
2014:            ELSE1945:            ELSE
2015:               WRITE(44,'(F20.10)') DIAG(J1)1946:               WRITE(44,'(F20.10)') DIAG(J1)
2016:            ENDIF1947:            ENDIF
2017:            WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N)1948:            WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N)
2018:            CALL FLUSH(44)1949:            CALL FLUSH(44)
2019:            RETURN1950:            RETURN
2020:         ENDIF1951:         ENDIF
2021:      ENDDO1952:      ENDDO
2022:   ENDIF1953:   ENDIF


r31191/geopt.f 2016-09-23 11:30:17.222813304 +0100 r31190/geopt.f 2016-09-23 11:30:18.226826687 +0100
 85:  85: 
 86:       DOUBLE PRECISION ETS,EPLUS,EMINUS 86:       DOUBLE PRECISION ETS,EPLUS,EMINUS
 87:       COMMON /OEPATH/ ETS,EPLUS,EMINUS 87:       COMMON /OEPATH/ ETS,EPLUS,EMINUS
 88:  88: 
 89: ! hk286 89: ! hk286
 90:       DOUBLE PRECISION :: XCOORDS (3*NATOMS) 90:       DOUBLE PRECISION :: XCOORDS (3*NATOMS)
 91:  91: 
 92: ! cs778 92: ! cs778
 93:       DOUBLE PRECISION :: SOVER 93:       DOUBLE PRECISION :: SOVER
 94:  94: 
 95:       PROD = 0.0D0  ! Initialise here, so it can be used by every other IF block 
 96:  
 97:       INFO=0 95:       INFO=0
 98:       LWORK=33*3*NATOMS 96:       LWORK=33*3*NATOMS
 99:       ILWORK=33*3*NATOMS 97:       ILWORK=33*3*NATOMS
100:       IF (NENDHESS.LE.0) NENDHESS=NOPT 98:       IF (NENDHESS.LE.0) NENDHESS=NOPT
101: C 99: C
102: C  *************** two-ended pathways ********************100: C  *************** two-ended pathways ********************
103: C101: C
104:       TTDONE=.FALSE.102:       TTDONE=.FALSE.
105:       IF (TWOENDS) THEN103:       IF (TWOENDS) THEN
106:          CALL TWOEND(ENERGY,VNEW,VECS,Q)104:          CALL TWOEND(ENERGY,VNEW,VECS,Q)
301:             IF (BFGSMINT.OR.BSMIN.OR.RKMIN) CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)299:             IF (BFGSMINT.OR.BSMIN.OR.RKMIN) CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
302:             CALL CHECKIND(Q,MFLAG,INEG,ENERGY,EVALMIN,EVALMAX,.FALSE.)300:             CALL CHECKIND(Q,MFLAG,INEG,ENERGY,EVALMIN,EVALMAX,.FALSE.)
303:          ENDIF301:          ENDIF
304:       ENDIF302:       ENDIF
305: !303: !
306: ! If we have just finished the BHINTERPT or BISECTT procedure and geopt304: ! If we have just finished the BHINTERPT or BISECTT procedure and geopt
307: ! is being called to make min.data.info MFLAG is not set, and could be false.305: ! is being called to make min.data.info MFLAG is not set, and could be false.
308: !306: !
309:       IF ((BHINTERPT.OR.BISECTT).AND.(.NOT.NEWCONNECTT).AND.(.NOT.REOPTIMISEENDPOINTS)) MFLAG=.TRUE.307:       IF ((BHINTERPT.OR.BISECTT).AND.(.NOT.NEWCONNECTT).AND.(.NOT.REOPTIMISEENDPOINTS)) MFLAG=.TRUE.
310:       IF (.NOT.MFLAG) GOTO 11 ! skip the calculations for min.data.info construction308:       IF (.NOT.MFLAG) GOTO 11 ! skip the calculations for min.data.info construction
311:  
312: ! Set up the required information (i.e. frequencies) for min.data.info to be printed 
313: C309: C
314: C DAE310: C DAE
315: C may want to calculate rate from system which has no second derivatives in the potential311: C may want to calculate rate from system which has no second derivatives in the potential
316: C There are two possibilites for doing this312: C There are two possibilites for doing this
317: C One is constructing a numerical hessian and diagonalise - remember mass weighting313: C One is constructing a numerical hessian and diagonalise - remember mass weighting
318: C Or could make approximation that only the highest (lowest?!) frequency mode is important and use314: C Or could make approximation that only the highest (lowest?!) frequency mode is important and use
319: C the XMYLBFGS routine to give estimate of this for stationary points - again need mass weighting315: C the XMYLBFGS routine to give estimate of this for stationary points - again need mass weighting
320: C316: C
321: C Alternatively may want to do BFGS minimisation on a system where second derivatives317: C Alternatively may want to do BFGS minimisation on a system where second derivatives
322: C are available but only use the sec.der. at the end318: C are available but only use the sec.der. at the end
323: C319: C
324: C In all cases aim to to mimic output of efol of 'Log product ...', so that320: C In all cases aim to to mimic output of efol of 'Log product ...', so that
325: C grep in Filthy or pathsample can cope without modification.321: C grep in Filthy or pathsample can cope without modification.
326: C should remove need for paul;s gmfreq and crap (sic!) programs which he introduced to do322: C should remove need for paul;s gmfreq and crap (sic!) programs which he introduced to do
327: C mass weighting323: C mass weighting
328: C324: C
329: C May want to rewrite above CHECKINDEX code to avoid duplication of effort325: C May want to rewrite above CHECKINDEX code to avoid duplication of effort
330: C326: C
331: 327:       IF (NOFRQS) THEN
332: ! If we're not attempting to calculate frequencies, set the log product to 1 as a dummy, and do no more. 
333:       IF (NOFRQS) THEN   
334:          PROD=1.0D0328:          PROD=1.0D0
335: !329: !
336: ! REOPTIMISEENDPOINTS is set to false after a bhinterp run is finished330: ! REOPTIMISEENDPOINTS is set to false after a bhinterp run is finished
337: ! so that second derivatives are only calculated at the end.331: ! so that second derivatives are only calculated at the end.
338: !332: !
339: 333: 
340: ! ENDHESS should be set if you want to print frequencies to min.data.info 
341: !     ELSEIF (ENDHESS .AND. ((.NOT. (REOPTIMISEENDPOINTS.AND.(BHINTERPT.OR.BISECTT))).AND.DUMPDATAT)) THEN334: !     ELSEIF (ENDHESS .AND. ((.NOT. (REOPTIMISEENDPOINTS.AND.(BHINTERPT.OR.BISECTT))).AND.DUMPDATAT)) THEN
342:       ELSEIF (ENDHESS .AND. (.NOT. (REOPTIMISEENDPOINTS.AND.(BHINTERPT.OR.BISECTT)))) THEN335:       ELSEIF (ENDHESS .AND. (.NOT. (REOPTIMISEENDPOINTS.AND.(BHINTERPT.OR.BISECTT)))) THEN
343: 336: 
344: ! This first block calculates the Hessian (either analytical or numerical) to extract the frequencies. 
345: ! jmc49 skipping the calculation of the Hessian here if RIGIDINIT is true (unless LOWESTFRQT is also true)337: ! jmc49 skipping the calculation of the Hessian here if RIGIDINIT is true (unless LOWESTFRQT is also true)
346: ! to avoid duplication of effort because GENRIGID_EIGENVALUES will calculate everything later 338: ! to avoid duplication of effort because GENRIGID_EIGENVALUES will calculate everything later 
347:        IF ((.NOT. RIGIDINIT) .OR. (RIGIDINIT .AND. LOWESTFRQT)) THEN339:        IF ((.NOT. RIGIDINIT) .OR. (RIGIDINIT .AND. LOWESTFRQT)) THEN
348:          IF (ENDNUMHESS) THEN340:          IF (ENDNUMHESS) THEN
349:             IF (UNRST) THEN341:             IF (UNRST) THEN
350:                CALL MAKENUMINTHESS(NINTS,NATOMS)342:                CALL MAKENUMINTHESS(NINTS,NATOMS)
351:                WRITE (*,'(A)') ' geopt> Value below calculated from numerical hessian' 343:                WRITE (*,'(A)') ' geopt> Value below calculated from numerical hessian' 
352:                CALL GETSTUFF(KD,NNZ,NINTB)344:                CALL GETSTUFF(KD,NNZ,NINTB)
353:                CALL INTSECDET(Q,3*NATOMS,KD,NNZ,NINTB,EVALUES)345:                CALL INTSECDET(Q,3*NATOMS,KD,NNZ,NINTB,EVALUES)
354:                NEXMODES=0346:                NEXMODES=0
371: C363: C
372: ! hk286 - tells potential that we switch to moving frame - needed for normal modes calculations364: ! hk286 - tells potential that we switch to moving frame - needed for normal modes calculations
373:             RBAANORMALMODET = .TRUE.365:             RBAANORMALMODET = .TRUE.
374:             CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)  366:             CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)  
375:             RBAANORMALMODET = .FALSE.367:             RBAANORMALMODET = .FALSE.
376: ! hk286368: ! hk286
377:             WRITE (*,'(A)') ' geopt> Value below calculated from true hessian'369:             WRITE (*,'(A)') ' geopt> Value below calculated from true hessian'
378:          ENDIF370:          ENDIF
379:        ENDIF371:        ENDIF
380: 372: 
381: ! Next, we diagonalise the Hessian. We need to know how many zero eigenvalues there will be (NEXMODES),  
382: ! so we can ignore these when we come to print output. 
383:          IF (.NOT.UNRST) THEN373:          IF (.NOT.UNRST) THEN
384: C374: C
385: C dae DSEYV sorts eigenvalues so 6 zero ones should be at bottom375: C dae DSEYV sorts eigenvalues so 6 zero ones should be at bottom
386: C therefore do product up to 3*NATOMS-6376: C therefore do product up to 3*NATOMS-6
387: C this will only work for isolated non-linear systems - fine for charmm377: C this will only work for isolated non-linear systems - fine for charmm
388: C For general case need code involving ZT from efol378: C For general case need code involving ZT from efol
389: C379: C
390:  
391: ! sn402: possible todo: add an IF block here for RIGIDINIT? 
392:             NEXMODES=6380:             NEXMODES=6
393:             IF (BFGSMINT.AND.(.NOT.BFGSTST)) THEN381:             IF (BFGSMINT.AND.(.NOT.BFGSTST)) THEN
394:                NEXMODES=6382:                NEXMODES=6
395:             ELSEIF (BFGSTST) THEN383:             ELSEIF (BFGSTST) THEN
396:                ! If we have a transition state, there will be a negative mode as well as the 6 zeros. Exclude this. 
397:                NEXMODES=7384:                NEXMODES=7
398:             ELSEIF (INR.EQ.2) THEN ! this will detect higher order saddles385:             ELSEIF (INR.EQ.2) THEN ! this will detect higher order saddles
399:                NEXMODES=7386:                NEXMODES=7
400:             ENDIF387:             ENDIF
401:             ! No rotational modes in bulk systems 
402:             IF (BULKT) NEXMODES=NEXMODES-3388:             IF (BULKT) NEXMODES=NEXMODES-3
403:             IF (BULKT.AND.TWOD) NEXMODES=NATOMS+2389:             IF (BULKT.AND.TWOD) NEXMODES=NATOMS+2
404:             IF (PULLT.OR.EFIELDT) NEXMODES=4390:             IF (PULLT.OR.EFIELDT) NEXMODES=4
405: ! hk286391: ! hk286
406:             IF (GTHOMSONT .AND. (GTHOMMET .EQ. 5) ) NEXMODES = 3392:             IF (GTHOMSONT .AND. (GTHOMMET .EQ. 5) ) NEXMODES = 3
407:             IF (GTHOMSONT .AND. (GTHOMMET < 5) ) NEXMODES = 1393:             IF (GTHOMSONT .AND. (GTHOMMET < 5) ) NEXMODES = 1
408:             IF (TWOD) NEXMODES=NEXMODES+NATOMS394:             IF (TWOD) NEXMODES=NEXMODES+NATOMS
409:             IF (FREEZE) THEN395:             IF (FREEZE) THEN
410:                NEXMODES=3*NFREEZE396:                NEXMODES=3*NFREEZE
411:             ENDIF397:             ENDIF
429:             IF (STEALTHYT) THEN415:             IF (STEALTHYT) THEN
430:                IF (ENERGY.LT.1.0D-8) THEN416:                IF (ENERGY.LT.1.0D-8) THEN
431:                   NEXMODES=NATOMS*3 - STM*2417:                   NEXMODES=NATOMS*3 - STM*2
432:                ELSE418:                ELSE
433:                   NEXMODES=3419:                   NEXMODES=3
434:                ENDIF420:                ENDIF
435:             ENDIF421:             ENDIF
436:             IF (VARIABLES.OR.MIEFT .OR. NOTRANSROTT) NEXMODES=0422:             IF (VARIABLES.OR.MIEFT .OR. NOTRANSROTT) NEXMODES=0
437:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES423:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES
438: 424: 
439: ! Keyword-specific block 
440:             IF (LOWESTFRQT) THEN425:             IF (LOWESTFRQT) THEN
441: C426: C
442: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file427: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file
443: C  No mass-weighting here!428: C  No mass-weighting here!
444: C  'U' specifies that the upper triangle contains the Hessian.429: C  'U' specifies that the upper triangle contains the Hessian.
445: C  Need to save HESS before this call and restore afterwards.430: C  Need to save HESS before this call and restore afterwards.
446: C431: C
447:     432:     
448:                ABSTOL=DLAMCH('Safe  minimum')433:                ABSTOL=DLAMCH('Safe  minimum')
449:                IF (ALLOCATED(ZWK)) DEALLOCATE(ZWK)434:                IF (ALLOCATED(ZWK)) DEALLOCATE(ZWK)
450: !              ALLOCATE(ZWK(NOPT,NEXMODES+1))435: !              ALLOCATE(ZWK(NOPT,NEXMODES+1))
451:                ALLOCATE(ZWK(1,1))436:                ALLOCATE(ZWK(1,1))
452:                ALLOCATE(ZSAVE(NOPT,NOPT))437:                ALLOCATE(ZSAVE(NOPT,NOPT))
453:                ZSAVE(1:NOPT,1:NOPT)=HESS(1:NOPT,1:NOPT)438:                ZSAVE(1:NOPT,1:NOPT)=HESS(1:NOPT,1:NOPT)
454: 439:                IF (RBAAT) THEN  ! hk286
455:                IF (RIGIDINIT) THEN 
456:                   RBAANORMALMODET = .TRUE. 
457:                   CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO) 
458:                   RBAANORMALMODET = .FALSE. 
459:                   WRITE(*,*) "geopt> This combination of keywords (LOWESTFRQ) has not been checked - sn402" 
460:                ELSEIF (RBAAT) THEN  ! hk286 
461:                   RBAANORMALMODET = .TRUE.440:                   RBAANORMALMODET = .TRUE.
462:                   CALL NRMLMD (Q, EVALUES, .TRUE.)441:                   CALL NRMLMD (Q, EVALUES, .TRUE.)
463:                   RBAANORMALMODET = .FALSE.442:                   RBAANORMALMODET = .FALSE.
464:                   PRINT *, "geopt> This combination of keywords has not been checked carefully - hk286"443:                   PRINT *, "geopt> This combination of keywords has not been checked carefully - hk286"
465:                ELSE444:                ELSE
466:                   CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NEXMODES+1,ABSTOL,445:                   CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NEXMODES+1,ABSTOL,
467:      &                        NFOUND,EVALUES,446:      &                        NFOUND,EVALUES,
468:      &                        ZWK,NOPT,ISUPPZ,WORK,447:      &                        ZWK,NOPT,ISUPPZ,WORK,
469:      &                        LWORK, IWORK, ILWORK, INFO )448:      &                        LWORK, IWORK, ILWORK, INFO )
470:                ENDIF449:                ENDIF
471:                MINCURVE=EVALUES(NEXMODES+1)450:                MINCURVE=EVALUES(NEXMODES+1)
472:                DEALLOCATE(ZWK)451:                DEALLOCATE(ZWK)
473:                HESS(1:NOPT,1:NOPT)=ZSAVE(1:NOPT,1:NOPT)452:                HESS(1:NOPT,1:NOPT)=ZSAVE(1:NOPT,1:NOPT)
474:                DEALLOCATE(ZSAVE)453:                DEALLOCATE(ZSAVE)
475:                PRINT '(A,G20.10)',' geopt> lowest non-zero positive eigenvalue=',MINCURVE454:                PRINT '(A,G20.10)',' geopt> lowest non-zero positive eigenvalue=',MINCURVE
476:                DO J1=1,NEXMODES+1455:                DO J1=1,NEXMODES+1
477:                   PRINT '(I8,G20.10)',J1,EVALUES(J1)456:                   PRINT '(I8,G20.10)',J1,EVALUES(J1)
478:                ENDDO457:                ENDDO
479:             ENDIF  ! End of IF(LOWESTFRQT)458:             ENDIF
480:  
481: !           PRINT '(A,I6,A)',' geopt> Ignoring the ',NEXMODES,' lowest Hessian eigenvalues'459: !           PRINT '(A,I6,A)',' geopt> Ignoring the ',NEXMODES,' lowest Hessian eigenvalues'
482: !           WRITE(*,*) 'ATMASS=', ATMASS(:)460: !           WRITE(*,*) 'ATMASS=', ATMASS(:)
483:  
484: ! hk286461: ! hk286
485: ! Mass-weight the Hessian 
486: ! We don't mass-weight rigid body systems here. It's done in the GENRIGID_EIGENVALUES routine (NRMLMD for the old RBAA scheme) 
487:             IF (.NOT. RIGIDINIT) THEN462:             IF (.NOT. RIGIDINIT) THEN
488:                IF (CHRMMT) THEN463:                IF (CHRMMT) THEN
489:                   CALL MASSWT2(NATOMS,ATMASS,Q,VNEW,.TRUE.) ! just does the Hessian464:                   CALL MASSWT2(NATOMS,ATMASS,Q,VNEW,.TRUE.) ! just does the Hessian
490:                ELSEIF (RINGPOLYMERT) THEN465:                ELSEIF (RINGPOLYMERT) THEN
491:                   CALL MASSWTRP(NOPT,RPMASSES,RPDOF) ! just does the Hessian466:                   CALL MASSWTRP(NOPT,RPMASSES,RPDOF) ! just does the Hessian
492:                ELSEIF (RBAAT.OR.VARIABLES) THEN ! hk286467:                ELSEIF (RBAAT.OR.VARIABLES) THEN ! hk286
493:                ELSE468:                ELSE
494:                   PRINT '(A)',' geopt> calling MASSWT' !! DJW469:                   PRINT '(A)',' geopt> calling MASSWT' !! DJW
495:                   CALL MASSWT(NATOMS,ATMASS,Q,VNEW,.TRUE.) ! does the Hessian, coordinates and gradient vector470:                   CALL MASSWT(NATOMS,ATMASS,Q,VNEW,.TRUE.) ! does the Hessian, coordinates and gradient vector
496:                ENDIF471:                ENDIF
497:             ENDIF472:             ENDIF
498:  
499: ! Keyword-specific block 
500: C473: C
501: C calling DSEYV with 'N' as we;re not interested in the eigenvectors474: C calling DSEYV with 'N' as we;re not interested in the eigenvectors
502: C Hard-wired calculation of eigenvalues and eigenvectors for diala system475: C Hard-wired calculation of eigenvalues and eigenvectors for diala system
503: C in order to get phi and psi numerical derivatives.476: C in order to get phi and psi numerical derivatives.
504: C477: C
505: C In the new scheme we only include the derivative that is largest in magnitude to the478: C In the new scheme we only include the derivative that is largest in magnitude to the
506: C corresponding order parameter.479: C corresponding order parameter.
507: C480: C
508:             IF (ORDERPARAMT) THEN481:             IF (ORDERPARAMT) THEN
509:  
510:                IF (RBAAT .OR. RIGIDINIT) THEN 
511:                   WRITE(*,*) "geopt> Warning- ORDERPARAM probably gives the wrong answers when used with rigid body systems" 
512:                ENDIF 
513:  
514: C If AMBER or NAB, need to undo the mass weighting  482: C If AMBER or NAB, need to undo the mass weighting  
515:                IF(.NOT.CHRMMT) THEN483:                IF(.NOT.CHRMMT) THEN
516:                  DO J1=1,NATOMS484:                  DO J1=1,NATOMS
517:                     AMASS=1/SQRT(ATMASS(J1))485:                     AMASS=1/SQRT(ATMASS(J1))
518:                     J3=3*J1486:                     J3=3*J1
519:                     Q(J3-2)=AMASS*Q(J3-2)487:                     Q(J3-2)=AMASS*Q(J3-2)
520:                     Q(J3-1)=AMASS*Q(J3-1)488:                     Q(J3-1)=AMASS*Q(J3-1)
521:                     Q(J3)=AMASS*Q(J3)489:                     Q(J3)=AMASS*Q(J3)
522:                  ENDDO490:                  ENDDO
523:                ENDIF491:                ENDIF
872: 444            CONTINUE840: 444            CONTINUE
873:                CLOSE(UNIT=9123)841:                CLOSE(UNIT=9123)
874:                CLOSE(UNIT=9124)842:                CLOSE(UNIT=9124)
875:                ENERGY=ESAVE843:                ENERGY=ESAVE
876:                EVALUES(1:3*NATOMS)=EVSAVE(1:3*NATOMS)844:                EVALUES(1:3*NATOMS)=EVSAVE(1:3*NATOMS)
877: !!!!!!!!!!!!!!!!!!!!!!!845: !!!!!!!!!!!!!!!!!!!!!!!
878: !846: !
879: ! End of order parameter block847: ! End of order parameter block
880: !848: !
881: !!!!!!!!!!!!!!!!!!!!!!!849: !!!!!!!!!!!!!!!!!!!!!!!
882:  
883: ! Back into the main branch of the subroutine. 
884: ! This is where we actually calculate the eigenvectors of the (mass-weighted) Hessian. 
885:             ELSE850:             ELSE
886:                IF (DUMPV) THEN ! use job type 'V' to get eigenvectors851:                IF (DUMPV) THEN ! use job type 'V' to get eigenvectors
887: ! hk286852: ! hk286
888:                   IF (RIGIDINIT) THEN853:                   IF (RIGIDINIT) THEN
889:                      RBAANORMALMODET = .TRUE.854:                      RBAANORMALMODET = .TRUE.
890:                      ! Eigenvectors will be calculated because DUMPV = .TRUE. 
891:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)855:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)
892:                      RBAANORMALMODET = .FALSE.856:                      RBAANORMALMODET = .FALSE.
893:                   ELSEIF (RBAAT) THEN857:                   ELSEIF (RBAAT) THEN
894:                      RBAANORMALMODET = .TRUE.858:                      RBAANORMALMODET = .TRUE.
895:                      ! Last argument to NRMLMD is EIGENVECTORT. Set this to .TRUE. so eigenvectors are calculated 
896:                      CALL NRMLMD (Q, EVALUES, .TRUE.)859:                      CALL NRMLMD (Q, EVALUES, .TRUE.)
897:                      RBAANORMALMODET = .FALSE.860:                      RBAANORMALMODET = .FALSE.
898:                   ELSE861:                   ELSE
899:                      CALL DSYEV('V','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)862:                      CALL DSYEV('V','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)
900:                   ENDIF863:                   ENDIF
901:                ELSE   ! Then we don't need to calculate the eigenvectors. Eigenvalues only.864:                ELSE 
902: ! hk286 - optional number of ENDHESS not yet implemented865: ! hk286 - optional number of ENDHESS not yet implemented
903:                   IF (RIGIDINIT) THEN866:                   IF (RIGIDINIT) THEN
904:                      RBAANORMALMODET = .TRUE.867:                      RBAANORMALMODET = .TRUE.
905:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)868:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)
906:                      RBAANORMALMODET = .FALSE.869:                      RBAANORMALMODET = .FALSE.
907:                   ELSE870:                   ELSE
908:                      IF ((NENDHESS.GE.NOPT).OR.(VARIABLES.AND.(NENDHESS.GE.NATOMS))) THEN871:                      IF ((NENDHESS.GE.NOPT).OR.(VARIABLES.AND.(NENDHESS.GE.NATOMS))) THEN
909: C                       CALL DSYEV('N','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)872: C                       CALL DSYEV('N','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)
910: C  csw34> changed the call used here to be the same as if NENDHESS < NOPT below873: C  csw34> changed the call used here to be the same as if NENDHESS < NOPT below
911: C         as this leads to a 20% speed increase!874: C         as this leads to a 20% speed increase!
930:                         ELSE893:                         ELSE
931:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NENDHESS,ABSTOL,894:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NENDHESS,ABSTOL,
932:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )895:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )
933:                         ENDIF896:                         ENDIF
934:                         DEALLOCATE(ZWK)897:                         DEALLOCATE(ZWK)
935:                      ENDIF898:                      ENDIF
936:                   ENDIF899:                   ENDIF
937:                ENDIF900:                ENDIF
938:             ENDIF901:             ENDIF
939: 902: 
940: ! By this point we have an array EVALUES containing the eigenvalues (probably ordered from smallest to largest) and 
941: ! if VDUMP is set, the matrix HESS contains the corresponding eigenvectors as its columns. 
942:  
943: C903: C
944: C  MASSWT2 and MASSWTRP do not mass weight Q and VNEW, but MASSWT does. Need to undo this904: C  MASSWT2 and MASSWTRP do not mass weight Q and VNEW, but MASSWT does. Need to undo this
945: C  if DUMPDATAT is .TRUE. for comparison with pathsample (which uses unit masses).905: C  if DUMPDATAT is .TRUE. for comparison with pathsample (which uses unit masses).
946: C  Probably best to always undo it, in case we need non-mass-weighted Q somewhere906: C  Probably best to always undo it, in case we need non-mass-weighted Q somewhere
947: C  further down.907: C  further down.
948: C908: C
949: C           IF (DUMPDATAT.AND.(.NOT.(CHRMMT.OR.RINGPOLYMERT))) THEN909: C           IF (DUMPDATAT.AND.(.NOT.(CHRMMT.OR.RINGPOLYMERT))) THEN
950:             IF (.NOT. (CHRMMT .OR. RBAAT)) THEN ! hk286910:             IF (.NOT. (CHRMMT .OR. RBAAT)) THEN ! hk286
951:                IF (.NOT.(RIGIDINIT.OR.VARIABLES)) THEN ! jmc49 bugfix, because the mass weighting is not done if rigidinit 911:                IF (.NOT.(RIGIDINIT.OR.VARIABLES)) THEN ! jmc49 bugfix, because the mass weighting is not done if rigidinit 
952:                   DO J1=1,NATOMS912:                   DO J1=1,NATOMS
953:                      AMASS=1/SQRT(ATMASS(J1))913:                      AMASS=1/SQRT(ATMASS(J1))
954:                      J3=3*J1914:                      J3=3*J1
955:                      Q(J3-2)=AMASS*Q(J3-2)915:                      Q(J3-2)=AMASS*Q(J3-2)
956:                      Q(J3-1)=AMASS*Q(J3-1)916:                      Q(J3-1)=AMASS*Q(J3-1)
957:                      Q(J3)=AMASS*Q(J3)917:                      Q(J3)=AMASS*Q(J3)
958:                   ENDDO918:                   ENDDO
959:                ENDIF919:                ENDIF
960:             ENDIF920:             ENDIF
961: 921: 
962:             ! This orders the eigenvalues from largest to smallest (despite the name of the subroutine...) 
963:             if (evalues(1).lt.evalues(NENDHESS)) call eigensort_val_asc(evalues,hess,NENDHESS,3*natoms)922:             if (evalues(1).lt.evalues(NENDHESS)) call eigensort_val_asc(evalues,hess,NENDHESS,3*natoms)
964:             IF (INFO.NE.0) PRINT*,'WARNING - INFO=',INFO,' in DSYEV'923:             IF (INFO.NE.0) PRINT*,'WARNING - INFO=',INFO,' in DSYEV'
965: 924:          ENDIF
966:          ENDIF   ! End of IF(.NOT. UNRST) (That was a very long block!) 
967: !925: !
968: ! Shift zero eigenvalues for Thomson problem. DJW926: ! Shift zero eigenvalues for Thomson problem. DJW
969: !927: !
970:          IF (GTHOMSONT) THEN928:          IF (GTHOMSONT) THEN
971:             EVALUES(2*NATOMS+1:3*NATOMS)=1.0D10929:             EVALUES(2*NATOMS+1:3*NATOMS)=1.0D10
972:             CALL EIGENSORT_VAL_ASC(EVALUES,HESS,NENDHESS,3*NATOMS)930:             CALL EIGENSORT_VAL_ASC(EVALUES,HESS,NENDHESS,3*NATOMS)
973:          ENDIF931:          ENDIF
974: !932: !
975: ! The test below will not necessarily spot a stationary point of the wrong index.933: ! The test below will not necessarily spot a stationary point of the wrong index.
976: ! We could read a zero eigenvalue that is > 0 and no error message will result.934: ! We could read a zero eigenvalue that is > 0 and no error message will result.
977: !935: !
978: ! sn402: Could this block be moved further down so that we have all the unit conversion stuff in the same place? 
979:          IF (DEBUG.OR.AMHT.OR.CASTEP.OR.VASP.OR.QCHEM.OR.RINGPOLYMERT.OR.ONETEP) THEN936:          IF (DEBUG.OR.AMHT.OR.CASTEP.OR.VASP.OR.QCHEM.OR.RINGPOLYMERT.OR.ONETEP) THEN
980: ! hk286937: ! hk286
981:             IF (RIGIDINIT) THEN938:             IF (RIGIDINIT) THEN
982:                PRINT '(A)', ' geopt> RIGIDINIT IS ON'939:                PRINT '(A)', ' geopt> RIGIDINIT IS ON'
983:                PRINT '(A,I6,A)',' geopt> ',DEGFREEDOMS,' Hessian eigenvalues:'940:                PRINT '(A,I6,A)',' geopt> ',DEGFREEDOMS,' Hessian eigenvalues:'
984:                PRINT '(6G20.10)',EVALUES(NENDHESS-DEGFREEDOMS+1:NENDHESS)941:                PRINT '(6G20.10)',EVALUES(NENDHESS-DEGFREEDOMS+1:NENDHESS)
985:             ELSEIF (RBAAT) THEN942:             ELSEIF (RBAAT) THEN
986:                ! sn402: need to change this here, when I've decided what I'm doing about unit conversion. 
987:                PRINT '(A,I6,A)',' geopt> ',NENDHESS,' Hessian eigenvalues in cm^-1:'943:                PRINT '(A,I6,A)',' geopt> ',NENDHESS,' Hessian eigenvalues in cm^-1:'
988:                DO J1 = 1, NENDHESS944:                DO J1 = 1, NENDHESS
989:                   ! sn402: this should absolutely not happen! We use EVALUES later on, and we do NOT want the behaviour to depend on whether DEBUG is set!945:                   IF (EVALUES(J1) < 0.0D0) EVALUES(J1) = - EVALUES(J1)
990: !                  IF (EVALUES(J1) < 0.0D0) EVALUES(J1) = - EVALUES(J1) 
991:                ENDDO946:                ENDDO
992:                ! sn402: need to sort this bit out. Use a new array to store ABS(EVALUES). Also need to decide what I'm doing with the unit conversions 
993:                PRINT '(6G20.10)', DSQRT(EVALUES(1:NENDHESS)) / 8.0D10 / ATAN(1.0D0) / 2.998D0947:                PRINT '(6G20.10)', DSQRT(EVALUES(1:NENDHESS)) / 8.0D10 / ATAN(1.0D0) / 2.998D0
994:             ELSE948:             ELSE
995:                PRINT '(A,I6,A)',' geopt> ',NENDHESS,' Hessian eigenvalues:'949:                PRINT '(A,I6,A)',' geopt> ',NENDHESS,' Hessian eigenvalues:'
996:                PRINT '(6G20.10)',EVALUES(1:NENDHESS)950:                PRINT '(6G20.10)',EVALUES(1:NENDHESS)
997:             ENDIF951:             ENDIF
998:             IF (CASTEP.OR.ONETEP) THEN952:             IF (CASTEP.OR.ONETEP) THEN
999:                PRINT '(A,I6,A)',' geopt> ',NENDHESS, 953:                PRINT '(A,I6,A)',' geopt> ',NENDHESS, 
1000:      &                          ' normal mode frequencies in Hz and wavenumbers, assuming eV and Angstrom units for input:'954:      &                          ' normal mode frequencies in Hz and wavenumbers, assuming eV and Angstrom units for input:'
1001:                IF (ONETEP) THEN955:                IF (ONETEP) THEN
1002:                   PRINT '(A,I6,A)',' geopt> ',NENDHESS, 956:                   PRINT '(A,I6,A)',' geopt> ',NENDHESS, 
1003:      &             ' normal mode frequencies in Hz and wavenumbers, assuming hartree and bohr units for input:'957:      &             ' normal mode frequencies in Hz and wavenumbers, assuming hartree and bohr units for input:'
1004:                   DO J1=1,NENDHESS958:                   DO J1=1,NENDHESS
1005:                      IF (EVALUES(J1).GT.0.0D0) THEN959:                      IF (EVALUES(J1).GT.0.0D0) THEN
1006:                         PRINT '(I6,2G20.10)',J1,SQRT(EVALUES(J1)*9.3757D29)/(2*3.141592654D0), 960:                         PRINT '(I6,2G20.10)',J1,SQRT(EVALUES(J1)*9.3757D29)/(2*3.141592654D0), 
1007:      &                                         SQRT(EVALUES(J1)*9.3757D29)/(2*3.141592654D0*2.998D10)961:      &                                         SQRT(EVALUES(J1)*9.3757D29)/(2*3.141592654D0*2.998D10)
1008:                      ELSE962:                      ELSE
1009:                         PRINT '(I6,2(G20.10,A2))',J1,SQRT(-EVALUES(J1)*9.3757D29)/(2*3.141592654D0),' i',963:                         PRINT '(I6,2(G20.10,A2))',J1,SQRT(-EVALUES(J1)*9.3757D29)/(2*3.141592654D0),' i',
1010:      &                                      SQRT(-EVALUES(J1)*9.3757D29)/(2*3.141592654D0*2.998D10),' i'964:      &                                      SQRT(-EVALUES(J1)*9.3757D29)/(2*3.141592654D0*2.998D10),' i'
1011:                      ENDIF965:                         ENDIF
1012:                   ENDDO966:                   ENDDO
1013:                ELSE967:                ELSE
1014:                   PRINT '(A,I6,A)',' geopt> ',NENDHESS, 968:                   PRINT '(A,I6,A)',' geopt> ',NENDHESS, 
1015:      &             ' normal mode frequencies in Hz and wavenumbers, assuming eV and Angstrom units for input:'969:      &             ' normal mode frequencies in Hz and wavenumbers, assuming eV and Angstrom units for input:'
1016:                   DO J1=1,NENDHESS970:                   DO J1=1,NENDHESS
1017:                      IF (EVALUES(J1).GT.0.0D0) THEN971:                      IF (EVALUES(J1).GT.0.0D0) THEN
1018:                         PRINT '(I6,2G20.10)',J1,SQRT(EVALUES(J1)*9.75586D27)/(2*3.141592654D0), 972:                         PRINT '(I6,2G20.10)',J1,SQRT(EVALUES(J1)*9.75586D27)/(2*3.141592654D0), 
1019:      &                                         SQRT(EVALUES(J1)*9.75586D27)/(2*3.141592654D0*2.998D10)973:      &                                         SQRT(EVALUES(J1)*9.75586D27)/(2*3.141592654D0*2.998D10)
1020:                      ELSE974:                      ELSE
1021:                         PRINT '(I6,2(G20.10,A2))',J1,SQRT(-EVALUES(J1)*9.75586D27)/(2*3.141592654D0),' i',975:                         PRINT '(I6,2(G20.10,A2))',J1,SQRT(-EVALUES(J1)*9.75586D27)/(2*3.141592654D0),' i',
1022:      &                                      SQRT(-EVALUES(J1)*9.75586D27)/(2*3.141592654D0*2.998D10),' i'976:      &                                      SQRT(-EVALUES(J1)*9.75586D27)/(2*3.141592654D0*2.998D10),' i'
1023:                      ENDIF977:                         ENDIF
1024:                   ENDDO978:                   ENDDO
1025:                ENDIF979:                ENDIF
1026:            ! (Still in CASTEP/ONETEP block) 
1027: C980: C
1028: C  Added transformation back to Cartesian basis for Hessian eigenvectors,981: C  Added transformation back to Cartesian basis for Hessian eigenvectors,
1029: C  as in "Energy Landscapes" equation (2.51). Otherwise the eigenvector982: C  as in "Energy Landscapes" equation (2.51). Otherwise the eigenvector
1030: C  components refer to mass-weighted coordinates, not Cartesians. DJW 7/11/09983: C  components refer to mass-weighted coordinates, not Cartesians. DJW 7/11/09
1031: C  The eigenvectors of the mass-weighted Hessian correspond to the A matrix984: C  The eigenvectors of the mass-weighted Hessian correspond to the A matrix
1032: C  components A_{alpha gamma} for eigenvector gamma, and these vectors985: C  components A_{alpha gamma} for eigenvector gamma, and these vectors
1033: C  are orthonormal.986: C  are orthonormal.
1034: C  Second index of HESS labels the eigenvector, first index runs over components.987: C  Second index of HESS labels the eigenvector, first index runs over components.
1035: C  The transformed eigenvectors in the Cartesian basis are not orthogonal.988: C  The transformed eigenvectors in the Cartesian basis are not orthogonal.
1036: C989: C
1081: !  This isn't right - the formula in section V has reciprocal factors of g1034: !  This isn't right - the formula in section V has reciprocal factors of g
1082: !  in the frequencies as well. Need to check further.1035: !  in the frequencies as well. Need to check further.
1083: !1036: !
1084:                   QFAC=QFAC+(0.5D0)*DUMMY+0.5D0*LOG(RPBN*RPIMAGES/(6.283185307D0*RPBETA))1037:                   QFAC=QFAC+(0.5D0)*DUMMY+0.5D0*LOG(RPBN*RPIMAGES/(6.283185307D0*RPBETA))
1085:      &                              -0.5D0*PROD-(RPIMAGES-2)*LOG(RPBETA/RPIMAGES)1038:      &                              -0.5D0*PROD-(RPIMAGES-2)*LOG(RPBETA/RPIMAGES)
1086:                ELSE1039:                ELSE
1087:                ENDIF1040:                ENDIF
1088:                PRINT '(2(A,G20.10))',' geopt> ln(k_instanton^+ * Q^+)=',QFAC-(ETS-EPLUS)*RPBETA/RPIMAGES,' E+=',EPLUS1041:                PRINT '(2(A,G20.10))',' geopt> ln(k_instanton^+ * Q^+)=',QFAC-(ETS-EPLUS)*RPBETA/RPIMAGES,' E+=',EPLUS
1089:                PRINT '(2(A,G20.10))',' geopt> ln(k_instanton^- * Q^-)=',QFAC-(ETS-EMINUS)*RPBETA/RPIMAGES,' E-=',EMINUS1042:                PRINT '(2(A,G20.10))',' geopt> ln(k_instanton^- * Q^-)=',QFAC-(ETS-EMINUS)*RPBETA/RPIMAGES,' E-=',EMINUS
1090:             ENDIF1043:             ENDIF
1091:          ENDIF ! End of IF(DEBUG.OR. systems)1044:          ENDIF
1092:  
1093:          ! ZT and LZT determines which normal mode frequencies will be printed out. Generally, this means identifying the 
1094:          ! zero eigenvalues, which will be at the front (or perhaps the end) of the eigenvalue vector. 
1095:          ! In this case, all the normal modes are set to print - but this conflicts with the documentation which suggests 
1096:          ! that only non-zero eigenvalues should be printed. 
1097:          LZT(1:NOPT)=.TRUE.1045:          LZT(1:NOPT)=.TRUE.
1098:          ! sn402: Adding this bit to make the behaviour match that described in the documentation. 
1099:          DO J1=1,NOPT 
1100:             IF (ABS(EVALUES(J1)).LT.EVCUT) LZT(J1)=.FALSE. ! EVCUT determines the minimum magnitude for a non-zero eigenvalue 
1101:          ENDDO 
1102:          ! If we've sorted eigenvalues in descending order, the first set will be dummies for RIGIDINIT. 
1103:          IF(RIGIDINIT) LZT(:NENDHESS-DEGFREEDOMS) = .FALSE. 
1104: 1046: 
1105: ! sn402: This is a problem, because it has the effect of undoing the sorting from the earlier subroutine call. Since we've already1047: 
1106: ! got the data we would obtain from this call, I'm going to comment it out. This has required editing some earlier code, but that's 
1107: ! better than doing another hessian diagonalisation unnecessarily. 
1108: ! hk286 - compute normal modes if rigid body angle-axis framework is used, in preparation for VDUMP subroutine1048: ! hk286 - compute normal modes if rigid body angle-axis framework is used, in preparation for VDUMP subroutine
1109: !         IF (RBAAT) THEN1049:          IF (RBAAT) THEN
1110: !            RBAANORMALMODET = .TRUE.1050:             RBAANORMALMODET = .TRUE.
1111: !            CALL NRMLMD (Q, EVALUES, .TRUE.)1051:             CALL NRMLMD (Q, EVALUES, .TRUE.)
1112: !            RBAANORMALMODET = .FALSE.1052:             RBAANORMALMODET = .FALSE.
1113: !         ENDIF1053:          ENDIF
1114: ! hk2861054: ! hk286
1115:          IF (DUMPV) THEN1055:          IF (DUMPV) THEN
1116:            IF (RIGIDINIT) THEN1056:             IF (RIGIDINIT) THEN
1117:               ! sn402: should this be LZT here? Looks to me as though ZT is undefined.1057:                CALL GENRIGID_VDUMP(DIAG(1:DEGFREEDOMS),ZT(1:DEGFREEDOMS),NOPT,DEGFREEDOMS)
1118:               ! GENRIGID_EIGENVALUES is set up so that all dummy entries in DIAG (from rigidified degrees of freedom) are at 
1119:               ! the end of the vector,so will never be printed. 
1120: ! sn402: changing the next line on a trial basis 
1121: !               CALL GENRIGID_VDUMP(DIAG(1:DEGFREEDOMS),ZT(1:DEGFREEDOMS),NOPT,DEGFREEDOMS) 
1122:                CALL GENRIGID_VDUMP(EVALUES,LZT,NOPT,NOPT) 
1123:             ELSE1058:             ELSE
1124:                CALL VDUMP(EVALUES,LZT,NOPT,3*NATOMS)1059:                CALL VDUMP(EVALUES,LZT,NOPT,3*NATOMS)
1125:             ENDIF1060:             ENDIF
1126:          ENDIF1061:          ENDIF
1127:          PROD=0.0D01062:          PROD=0.0D0
1128:  
1129: ! Now we calculate the log product of positive frequencies, which will be written to min.data.info. This must be done regardless of whether we are dumping eigenvectors. 
1130: ! Unit conversions are added to the product in a separate block - see below. At the moment, the eigenvalues should be in internal units. 
1131: ! hk2861063: ! hk286
1132:          IF (RIGIDINIT) THEN1064:          IF (RIGIDINIT) THEN
1133:             ! Note, there is no separate IF(RBAAT) block here, because there are no dummy eigenvalues in that case. 
1134:             ! MINFRQ2 is the (log of the) smallest positive hessian eigenvalue (i.e. smallest square frequency) 
1135:             IF (NENDHESS-NEXMODES.GT.0) THEN1065:             IF (NENDHESS-NEXMODES.GT.0) THEN
1136:                MINFRQ2=LOG(EVALUES(NENDHESS-NEXMODES))1066:                MINFRQ2=LOG(EVALUES(NENDHESS-NEXMODES))
1137:             ELSE1067:             ELSE
1138:                ! Dummy value, if there are no positive eigenvalues. 
1139:                MINFRQ2=1.0D01068:                MINFRQ2=1.0D0
1140:             ENDIF1069:             ENDIF
1141:             EWARN=.FALSE.1070:             EWARN=.FALSE.
1142:             ! This slightly awkward index choice arises because the first block of EVALUES now contains dummy values set to 10^10 
1143:             ! and the real eigenvalues are then listed in descending order. 
1144:             DO I1=NENDHESS-DEGFREEDOMS+1,NENDHESS-NEXMODES1071:             DO I1=NENDHESS-DEGFREEDOMS+1,NENDHESS-NEXMODES
1145:                IF (I1.GT.NENDHESS-DEGFREEDOMS+1) THEN1072:                IF (I1.GT.NENDHESS-DEGFREEDOMS+1) THEN
1146:                   IF (EVALUES(I1-1).NE.0.0D0) THEN1073:                   IF (EVALUES(I1-1).NE.0.0D0) THEN
1147:                      IF (ABS(EVALUES(I1)/EVALUES(I1-1)).LT.1.0D-2) THEN1074:                      IF (ABS(EVALUES(I1)/EVALUES(I1-1)).LT.1.0D-2) THEN
1148:                         PRINT '(A,G20.10,A,G20.10)',' geopt> WARNING - decrease in magnitude of eigenvalues from ',EVALUES(I1-1),1075:                         PRINT '(A,G20.10,A,G20.10)',' geopt> WARNING - decrease in magnitude of eigenvalues from ',EVALUES(I1-1),
1149:      &                       ' to ',EVALUES(I1)1076:      &                       ' to ',EVALUES(I1)
1150:                         PRINT '(A)',' geopt> WARNING - this could indicate a stationary point of the wrong index'1077:                         PRINT '(A)',' geopt> WARNING - this could indicate a stationary point of the wrong index'
1151:                         EWARN=.TRUE.1078:                         EWARN=.TRUE.
1152:                      ENDIF1079:                      ENDIF
1153:                   ENDIF1080:                   ENDIF
1154:                ENDIF1081:                ENDIF
1155:                ! Calculate the actual log product here! 
1156:                IF (EVALUES(I1).GT.0.0D0) THEN1082:                IF (EVALUES(I1).GT.0.0D0) THEN
1157:                   PROD=PROD+DLOG(EVALUES(I1))1083:                   PROD=PROD+DLOG(EVALUES(I1))
1158:                ELSE1084:                ELSE
1159:                   IF (I1.LT.(NENDHESS-NEXMODES)) PRINT *,'Higher order saddle detected: eigenvalue ',EVALUES(I1)1085:                   IF (I1.LT.(NENDHESS-NEXMODES)) PRINT *,'Higher order saddle detected: eigenvalue ',EVALUES(I1)
1160:                   ! jmc put in this test mainly for pathsample purposes...1086:                   ! jmc put in this test mainly for pathsample purposes...
1161:                ENDIF 1087:                ENDIF 
1162:             ENDDO1088:             ENDDO
1163:          ELSE1089:          ELSE
1164:             IF (NENDHESS-NEXMODES.GT.0) THEN1090:             IF (NENDHESS-NEXMODES.GT.0) THEN
1165:                MINFRQ2=LOG(EVALUES(NENDHESS-NEXMODES))1091:                MINFRQ2=LOG(EVALUES(NENDHESS-NEXMODES))
1192:                      PRINT '(I6,2G20.10E3)',J1,SQRT(EVALUES(J1)*4.184D26)/(2*3.141592654D0),1118:                      PRINT '(I6,2G20.10E3)',J1,SQRT(EVALUES(J1)*4.184D26)/(2*3.141592654D0),
1193:      &                 SQRT(EVALUES(J1)*4.184D26)/(2*3.141592654D0*2.998D10)1119:      &                 SQRT(EVALUES(J1)*4.184D26)/(2*3.141592654D0*2.998D10)
1194:                   ELSE1120:                   ELSE
1195:                      PRINT '(I6,2(G20.10,A2))',J1,SQRT(-EVALUES(J1)*4.184D26)/(2*3.141592654D0),' i',1121:                      PRINT '(I6,2(G20.10,A2))',J1,SQRT(-EVALUES(J1)*4.184D26)/(2*3.141592654D0),' i',
1196:      &                 SQRT(-EVALUES(J1)*4.184D26)/(2*3.141592654D0*2.998D10),' i'1122:      &                 SQRT(-EVALUES(J1)*4.184D26)/(2*3.141592654D0*2.998D10),' i'
1197:                   ENDIF1123:                   ENDIF
1198:                ENDDO1124:                ENDDO
1199:             ENDIF1125:             ENDIF
1200:          ENDIF1126:          ENDIF
1201: 1127: 
1202: ! Applying a unit conversion to PROD for the potentials which need it. Remember that PROD is a sum of logs, so multiplying 
1203: ! the frequencies by a conversion factor is the same as adding the conversion factor to the sum (which is what we do here). 
1204:          IF (CHRMMT.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.SDT.OR.TTM3T.OR.QTIP4PFT) THEN1128:          IF (CHRMMT.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.SDT.OR.TTM3T.OR.QTIP4PFT) THEN
1205: C1129: C
1206: C if charmm need to convert this to (radian/s)^2, rather than charmm units1130: C if charmm need to convert this to (radian/s)^2, rather than charmm units
1207: C conversion factor for this is 4.184 x 10^261131: C conversion factor for this is 4.184 x 10^26
1208: C same for AMBER1132: C same for AMBER
1209: C1133: C
1210:             MINFRQ2=MINFRQ2+LOG(4.184D26)1134:             MINFRQ2=MINFRQ2+LOG(4.184D26)
1211: ! hk2861135: ! hk286
1212:             IF (RIGIDINIT) THEN1136:             IF (RIGIDINIT) THEN
1213:                PROD=PROD+(DEGFREEDOMS-NEXMODES)*DLOG(4.184D26)1137:                PROD=PROD+(DEGFREEDOMS-NEXMODES)*DLOG(4.184D26)
1287:             IF (DUMPFRQST) THEN1211:             IF (DUMPFRQST) THEN
1288:                 IF (FILTH.EQ.0) THEN1212:                 IF (FILTH.EQ.0) THEN
1289:                    OPEN(91220,FILE='frqs.dump',POSITION='APPEND',ACTION='WRITE',STATUS='UNKNOWN')1213:                    OPEN(91220,FILE='frqs.dump',POSITION='APPEND',ACTION='WRITE',STATUS='UNKNOWN')
1290:                 ELSE1214:                 ELSE
1291:                    OPEN(91220,FILE='frqs.dump.'//TRIM(ADJUSTL(FILTHSTR)),POSITION='APPEND',ACTION='WRITE',STATUS='UNKNOWN')1215:                    OPEN(91220,FILE='frqs.dump.'//TRIM(ADJUSTL(FILTHSTR)),POSITION='APPEND',ACTION='WRITE',STATUS='UNKNOWN')
1292:                 END IF1216:                 END IF
1293:                 WRITE(91220,*)  PROD1217:                 WRITE(91220,*)  PROD
1294:                 CLOSE(91220)1218:                 CLOSE(91220)
1295:             ENDIF1219:             ENDIF
1296:          ENDIF1220:          ENDIF
1297: 1221:       ENDIF
1298:       ENDIF  ! End of IF(ENDHESS .AND. ...) (This block contains most of the subroutine, at the time of writing) 
1299: 1222: 
1300:       IF (CHRMMT.AND.CALCDIHE) THEN1223:       IF (CHRMMT.AND.CALCDIHE) THEN
1301:          STOP 'Necessary CHARMM routines not implemented yet for NSEG>1'1224:          STOP 'Necessary CHARMM routines not implemented yet for NSEG>1'
1302: C         LSELECT=.FALSE.1225: C         LSELECT=.FALSE.
1303: C         CALL CHCALCRGYR(RGYR,Q,LSELECT)1226: C         CALL CHCALCRGYR(RGYR,Q,LSELECT)
1304: C         LNATIVE=.FALSE.1227: C         LNATIVE=.FALSE.
1305: C         CALL CHCALCNUMHB(NUMHB,Q,LNATIVE)1228: C         CALL CHCALCNUMHB(NUMHB,Q,LNATIVE)
1306: C         CALL CHCALCRMSD(RMSD,Q)1229: C         CALL CHCALCRMSD(RMSD,Q)
1307: C         WRITE(*,'(A,4X,F20.10)') 'Final Rmsd from reference structure=',RMSD1230: C         WRITE(*,'(A,4X,F20.10)') 'Final Rmsd from reference structure=',RMSD
1308: C         WRITE(*,'(A,4X,F20.10)') 'Final Radius of gyration=',RGYR1231: C         WRITE(*,'(A,4X,F20.10)') 'Final Radius of gyration=',RGYR
1316:          CALL UNRESCALCRGYR(RGYR,Q)1239:          CALL UNRESCALCRGYR(RGYR,Q)
1317:          WRITE(*,'(A,4X,F20.10)') ' Dihedral angle order parameter=',DIHE1240:          WRITE(*,'(A,4X,F20.10)') ' Dihedral angle order parameter=',DIHE
1318:          WRITE(*,'(A,4X,F20.10)') ' All angle order parameter=',ALLANG1241:          WRITE(*,'(A,4X,F20.10)') ' All angle order parameter=',ALLANG
1319:          WRITE(*,'(A,4X,F20.10)') ' Radius of gyration=',RGYR1242:          WRITE(*,'(A,4X,F20.10)') ' Radius of gyration=',RGYR
1320:       ENDIF1243:       ENDIF
1321: 1244: 
1322: ! jmc49 Unconditionally executing the lines below, to make sure the Cartesian coords 1245: ! jmc49 Unconditionally executing the lines below, to make sure the Cartesian coords 
1323: ! are written correctly to points.final etc.1246: ! are written correctly to points.final etc.
1324: ! hk2861247: ! hk286
1325:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN1248:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN
1326:          ! DUMQ will hold the cartesian coordinates. We don't change ATOMRIGIDCOORDT because we're only changing 
1327:          ! coordinates of a dummy array, not Q itself. 
1328:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, DUMQ, Q(1:DEGFREEDOMS))1249:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, DUMQ, Q(1:DEGFREEDOMS))
1329:       ELSE IF (GTHOMSONT) THEN1250:       ELSE IF (GTHOMSONT) THEN
1330: !jwrm2> With GTHOMSON, want Cartesian coordinates for moment of inertia calculation1251: !jwrm2> With GTHOMSON, want Cartesian coordinates for moment of inertia calculation
1331:          CALL GTHOMSONANGTOC(DUMQ(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)1252:          CALL GTHOMSONANGTOC(DUMQ(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)
1332:       ELSE1253:       ELSE
1333:          DO I1=1,3*NATOMS1254:          DO I1=1,3*NATOMS
1334:             DUMQ(I1)=Q(I1)1255:             DUMQ(I1)=Q(I1)
1335:          ENDDO1256:          ENDDO
1336:       ENDIF1257:       ENDIF
1337: 1258: 
1338:       INERTIAT=.TRUE.1259:       INERTIAT=.TRUE.
1339: C     IF ((UNRST.OR.CHRMMT).AND.INERTIAT) THEN1260: C     IF ((UNRST.OR.CHRMMT).AND.INERTIAT) THEN
1340:  
1341:       ! The next block is used for dumping to a min.data.info file. The frequencies calculated 
1342:       ! previously will be used, unless NOFRQS is set. 
1343:       IF ((DUMPDATAT).AND.INERTIAT) THEN1261:       IF ((DUMPDATAT).AND.INERTIAT) THEN
1344: C1262: C
1345: C  Pathsample uses unit masses, so we need to change to unit mass temporarily here,1263: C  Pathsample uses unit masses, so we need to change to unit mass temporarily here,
1346: C  and then put the right values back, otherwise the frequencies will be wrong1264: C  and then put the right values back, otherwise the frequencies will be wrong
1347: C  when we call geopt again from newconnect in a bhinterp run!!!!1265: C  when we call geopt again from newconnect in a bhinterp run!!!!
1348: C1266: C
1349:          DO J1=1,NATOMS1267:          DO J1=1,NATOMS
1350:             ATMASSSAVE(J1)=ATMASS(J1)1268:             ATMASSSAVE(J1)=ATMASS(J1)
1351:             ATMASS(J1)=1.0D0 1269:             ATMASS(J1)=1.0D0 
1352:          ENDDO1270:          ENDDO
1353: 1271: 
1354:          ! Get the inertia tensor 
1355:          CALL INERTIA2(DUMQ,ITX,ITY,ITZ)1272:          CALL INERTIA2(DUMQ,ITX,ITY,ITZ)
1356:          WRITE(*,'(A,4X,F20.10)') ' geopt> X component of inertia tensor=',ITX1273:          WRITE(*,'(A,4X,F20.10)') ' geopt> X component of inertia tensor=',ITX
1357:          WRITE(*,'(A,4X,F20.10)') ' geopt> Y component of inertia tensor=',ITY1274:          WRITE(*,'(A,4X,F20.10)') ' geopt> Y component of inertia tensor=',ITY
1358:          WRITE(*,'(A,4X,F20.10)') ' geopt> Z component of inertia tensor=',ITZ1275:          WRITE(*,'(A,4X,F20.10)') ' geopt> Z component of inertia tensor=',ITZ
1359:          HORDER=11276:          HORDER=1
1360:          DO J1=1,NATOMS1277:          DO J1=1,NATOMS
1361:             ATMASS(J1)=ATMASSSAVE(J1)1278:             ATMASS(J1)=ATMASSSAVE(J1)
1362:          ENDDO1279:          ENDDO
1363:       ENDIF1280:       ENDIF
1364: 1281: 
1387: !1304: !
1388:          IF (.NOT.(VARIABLES.OR.RINGPOLYMERT).AND.(.NOT.(ZSYM(NATOMS).EQ.'TH'))) CALL SYMMETRY(HORDER,.TRUE.,DUMQ,IT)1305:          IF (.NOT.(VARIABLES.OR.RINGPOLYMERT).AND.(.NOT.(ZSYM(NATOMS).EQ.'TH'))) CALL SYMMETRY(HORDER,.TRUE.,DUMQ,IT)
1389:          if (MACHINE) then1306:          if (MACHINE) then
1390:              call WriteOutFile(DUMQ,FNAMEF)1307:              call WriteOutFile(DUMQ,FNAMEF)
1391:          else1308:          else
1392:              CALL DUMPIT(DUMQ,FNAMEF)1309:              CALL DUMPIT(DUMQ,FNAMEF)
1393:          endif1310:          endif
1394:       ENDIF1311:       ENDIF
1395: 11    CONTINUE1312: 11    CONTINUE
1396:       IF ((BHINTERPT.OR.BISECTT).AND.(.NOT.NEWCONNECTT).AND.(.NOT.REOPTIMISEENDPOINTS)) THEN ! do nothing1313:       IF ((BHINTERPT.OR.BISECTT).AND.(.NOT.NEWCONNECTT).AND.(.NOT.REOPTIMISEENDPOINTS)) THEN ! do nothing
1397:       ! sn402: I think MFLAG is used when we have converged to a minimum and want to write it to min.data.info. 
1398:       ELSE IF (MFLAG) THEN1314:       ELSE IF (MFLAG) THEN
1399:         ! csw34> Mass weighted normal mode dumping 1315: ! csw34> Mass weighted normal mode dumping 
1400:          IF (ENDHESS .AND. (.NOT. (REOPTIMISEENDPOINTS.AND.(BHINTERPT.OR.BISECTT)))) THEN1316:          IF (ENDHESS .AND. (.NOT. (REOPTIMISEENDPOINTS.AND.(BHINTERPT.OR.BISECTT)))) THEN
1401: ! in this case we have already called vdump. We don;t want to do it again!1317: ! in this case we have already called vdump. We don;t want to do it again!
1402: ! This IF statement (ENDHESS AND...) is equivalent to the one near the start of the subroutine. 
1403: ! So quite a lot of the code here duplicates the earlier behaviour. 
1404:          ELSEIF (DUMPV.AND.ALLVECTORS) THEN1318:          ELSEIF (DUMPV.AND.ALLVECTORS) THEN
1405:                 DO J1=1,3*NATOMS1319:                 DO J1=1,3*NATOMS
1406:                    ZT(J1)=.TRUE.1320:                    ZT(J1)=.TRUE.
1407:                 ENDDO1321:                 ENDDO
1408: ! Set the first six normal modes to not be printed. For a non-linear1322: ! Set the first six normal modes to not be printed. For a non-linear
1409: ! molecule, these correspond to the pure rotations and translations1323: ! molecule, these correspond to the pure rotations and translations
1410:                 DO J1=1,MIN(6,NOPT)1324:                 DO J1=1,MIN(6,NOPT)
1411:                    ZT(J1)=.FALSE.1325:                    ZT(J1)=.FALSE.
1412:                 ENDDO1326:                 ENDDO
1413: 1327:                                                                
1414: ! Here we compute the frequencies/eigenvalues themselves.                             
1415: ! Can't mass weight an already diagonalised Hessian so call potential again1328: ! Can't mass weight an already diagonalised Hessian so call potential again
1416: ! hk286 - tells potential that we switch to moving frame - needed for normal modes calculations1329: ! hk286 - tells potential that we switch to moving frame - needed for normal modes calculations
1417:                 RBAANORMALMODET = .TRUE.1330:                 RBAANORMALMODET = .TRUE.
1418:                 CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)1331:                 CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
1419:                 RBAANORMALMODET = .FALSE.1332:                 RBAANORMALMODET = .FALSE.
1420: ! hk286 - do better than putting this here?1333: ! hk286 - do better than putting this here?
1421:                 IF (RIGIDINIT) THEN1334:                 IF (RIGIDINIT) THEN
1422:                    CALL GENRIGID_EIGENVALUES(Q, ATMASS, DIAG, INFO)1335:                    CALL GENRIGID_EIGENVALUES(Q, ATMASS, DIAG, INFO)
1423:                 ELSE1336:                 ELSE
1424: ! Mass weight hessian - using MASSWT would also change the coordinates1337: ! Mass weight hessian - using MASSWT would also change the coordinates
1425:                    IF (MWVECTORS) THEN1338:                    IF (MWVECTORS) THEN
1426:                       IF (RINGPOLYMERT) THEN1339:                       IF (RINGPOLYMERT) THEN
1427:                          CALL MASSWTRP(NOPT,RPMASSES,RPDOF)1340:                          CALL MASSWTRP(NOPT,RPMASSES,RPDOF)
1428:                       ELSE1341:                       ELSE
1429:                          ! sn402: we're wasting a bit of time here when RBAAT is TRUE - we call this function 
1430:                          ! but we then overwrite the mass-weighted hessian in NRMLMD 
1431:                          CALL MASSWT2(NATOMS,ATMASS,Q,VNEW,.TRUE.)1342:                          CALL MASSWT2(NATOMS,ATMASS,Q,VNEW,.TRUE.)
1432:                       ENDIF1343:                       ENDIF
1433:                    ENDIF1344:                    ENDIF
1434: ! Diagonalise1345: ! Diagonalise
1435: ! hk286 - normal modes calculation, for rigid body angle-axis the diagonalization is rather complex1346: ! hk286 - normal modes calculation, for rigid body angle-axis the diagonalization is rather complex
1436:                    IF (RBAAT) THEN1347:                    IF (RBAAT) THEN
1437:                       RBAANORMALMODET = .TRUE.1348:                       RBAANORMALMODET = .TRUE.
1438:                       CALL NRMLMD (Q, DIAG, .TRUE.)1349:                       CALL NRMLMD (Q, DIAG, .TRUE.)
1439:                       RBAANORMALMODET = .FALSE.1350:                       RBAANORMALMODET = .FALSE.
1440:                    ELSE1351:                    ELSE
1443: ! hk2861354: ! hk286
1444: ! If we're freezing atoms, some zero eiganvalue modes can creep in that1355: ! If we're freezing atoms, some zero eiganvalue modes can creep in that
1445: ! are NOT real i.e. one atoms moving and all others stationary. Here, we1356: ! are NOT real i.e. one atoms moving and all others stationary. Here, we
1446: ! check and remove these using ZT1357: ! check and remove these using ZT
1447:                    IF (FREEZE) THEN1358:                    IF (FREEZE) THEN
1448:                       DO J1=1,3*NATOMS1359:                       DO J1=1,3*NATOMS
1449:                          IF (ABS(DIAG(J1)).LT.0.000001) ZT(J1)=.FALSE.1360:                          IF (ABS(DIAG(J1)).LT.0.000001) ZT(J1)=.FALSE.
1450:                       ENDDO1361:                       ENDDO
1451:                    ENDIF1362:                    ENDIF
1452:                 ENDIF1363:                 ENDIF
1453: ! Dump vectors - we will only get here if ENDHESS was set, and only get to the earlier block if it wasn't.1364: ! Dump vectors
1454: !hk2861365: !hk286
1455:                 IF (RIGIDINIT) THEN1366:                 IF (RIGIDINIT) THEN
1456:                    CALL GENRIGID_VDUMP(DIAG(1:DEGFREEDOMS),ZT(1:DEGFREEDOMS),NOPT,DEGFREEDOMS)1367:                    CALL GENRIGID_VDUMP(DIAG(1:DEGFREEDOMS),ZT(1:DEGFREEDOMS),NOPT,DEGFREEDOMS)
1457:                 ELSE1368:                 ELSE
1458:                    CALL VDUMP(DIAG,ZT,NOPT,3*NATOMS)1369:                    CALL VDUMP(DIAG,ZT,NOPT,3*NATOMS)
1459:                 ENDIF1370:                 ENDIF
1460: ! hk286 - do I need to change the following for CHARMM and AMBER?1371: ! hk286 - do I need to change the following for CHARMM and AMBER?
1461: 1372: 
1462: ! If using CHARMM, call the CHARMMDUMPMODES subroutine to output a PDB1373: ! If using CHARMM, call the CHARMMDUMPMODES subroutine to output a PDB
1463: ! containing one frame per mode and scaled atomic and residue1374: ! containing one frame per mode and scaled atomic and residue
1464: ! displacements - coords are in Q here!1375: ! displacements - coords are in Q here!
1465:                 IF (CHRMMT.AND.MWVECTORS) CALL CHARMMDUMPMODES(Q,DIAG,ZT,NOPT,3*NATOMS)1376:                 IF (CHRMMT.AND.MWVECTORS) CALL CHARMMDUMPMODES(Q,DIAG,ZT,NOPT,3*NATOMS)
1466: ! AMBER call to routine in amberinterface.f1377: ! AMBER call to routine in amberinterface.f
1467:                 IF ((AMBERT.OR.NABT).AND.MWVECTORS) THEN1378:                 IF ((AMBERT.OR.NABT).AND.MWVECTORS) THEN
1468:                         CALL A9DUMPMODES(DIAG,ZT,NOPT,3*NATOMS)1379:                         CALL A9DUMPMODES(DIAG,ZT,NOPT,3*NATOMS)
1469:                 ENDIF1380:                 ENDIF
1470:               IF(NORMALMODET) CALL VISNORMODES(Q)1381:               IF(NORMALMODET) CALL VISNORMODES(Q)
1471:          ENDIF ! End of IF(ENDHESS AND...)1382:          ENDIF
1472: 1383: 
1473:          GOODSTRUC1 = .TRUE.1384:          GOODSTRUC1 = .TRUE.
1474:          GOODSTRUC2 = .TRUE.1385:          GOODSTRUC2 = .TRUE.
1475:          IF (AMBER12T) THEN1386:          IF (AMBER12T) THEN
1476:             CALL CIS_TRANS_CHECK(Q, GOODSTRUC1)1387:             CALL CIS_TRANS_CHECK(Q, GOODSTRUC1)
1477:             CALL CHIRALITY_CHECK(Q, GOODSTRUC2)1388:             CALL CHIRALITY_CHECK(Q, GOODSTRUC2)
1478:          END IF1389:          END IF
1479:          IF (.NOT. GOODSTRUC1) THEN1390:          IF (.NOT. GOODSTRUC1) THEN
1480:             WRITE(*,'(A)') ' geopt> cistrans problem' 1391:             WRITE(*,'(A)') ' geopt> cistrans problem' 
1481:          END IF1392:          END IF
1499:             CALL H2OMODES(NATOMS/2,IPOT,Q,DIAG)1410:             CALL H2OMODES(NATOMS/2,IPOT,Q,DIAG)
1500:             PRINT '(A,I6,A)',' geopt> TIP normal mode frequencies in wavenumbers'1411:             PRINT '(A,I6,A)',' geopt> TIP normal mode frequencies in wavenumbers'
1501:             DO J1=1,3*NATOMS1412:             DO J1=1,3*NATOMS
1502:                IF (DIAG(J1).GT.0.0D0) THEN1413:                IF (DIAG(J1).GT.0.0D0) THEN
1503:                   PRINT '(I6,2G20.10)',J1,DIAG(J1)1414:                   PRINT '(I6,2G20.10)',J1,DIAG(J1)
1504:                ELSE1415:                ELSE
1505:                   PRINT '(I6,2(G20.10,A2))',J1,-DIAG(J1),' i'1416:                   PRINT '(I6,2(G20.10,A2))',J1,-DIAG(J1),' i'
1506:                ENDIF1417:                ENDIF
1507:             ENDDO1418:             ENDDO
1508:          ENDIF1419:          ENDIF
1509:       ELSE  ! End of IF(MFLAG)1420:       ELSE 
1510:          ! If we haven't converged to a minimum, print a bit and exit. 
1511:          IF (GRADSQ) WRITE(*,'(A,4F20.10)') ' g^2, RMS force and real energy and RMS=',ENERGY,RMS2,EREAL,RMS1421:          IF (GRADSQ) WRITE(*,'(A,4F20.10)') ' g^2, RMS force and real energy and RMS=',ENERGY,RMS2,EREAL,RMS
1512:          PRINT*1422:          PRINT*
1513:          IF (BFGSMINT.OR.(BFGSTST.AND.(HINDEX.EQ.0))) THEN1423:          IF (BFGSMINT.OR.(BFGSTST.AND.(HINDEX.EQ.0))) THEN
1514:               PRINT*,ITDONE,' steps completed without convergence to required tolerance'1424:               PRINT*,ITDONE,' steps completed without convergence to required tolerance'
1515:          ELSE1425:          ELSE
1516:               PRINT*,NSTEPS,' steps completed without convergence to required tolerance'1426:               PRINT*,NSTEPS,' steps completed without convergence to required tolerance'
1517:          ENDIF1427:          ENDIF
1518:          PRINT*1428:          PRINT*
1519:          CALL FLUSH(6)1429:          CALL FLUSH(6)
1520: !        STOP1430: !        STOP
1521:          IF (.NOT.MULTIJOBT) CALL EXIT(0)1431:          IF (.NOT.MULTIJOBT) CALL EXIT(0)
1522:       ENDIF1432:       ENDIF
1523:       CALL FLUSH(6)1433:       CALL FLUSH(6)
1524:  
1525: ! If we are making a min.data.info file, do it here. 
1526:       IF (DUMPDATAT) THEN 1434:       IF (DUMPDATAT) THEN 
1527: !        IF (EWARN.AND.MULTIJOBT) MFLAG=.FALSE.1435: !        IF (EWARN.AND.MULTIJOBT) MFLAG=.FALSE.
1528:          IF (MFLAG) THEN1436:          IF (MFLAG) THEN
1529: !1437: !
1530: ! min.data.info file is now opened on unit 881 in keyword.f1438: ! min.data.info file is now opened on unit 881 in keyword.f
1531: !1439: !
1532: !           OPEN(UNIT=100,FILE='min.data.info',STATUS='UNKNOWN')1440: !           OPEN(UNIT=100,FILE='min.data.info',STATUS='UNKNOWN')
1533:             ! sn402: need Cartesian coords to calculate the inertia tensor1441:             ! sn402: need Cartesian coords to calculate the inertia tensor
1534:             IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN1442:             IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN
1535:               CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, DUMQ, Q(1:DEGFREEDOMS))1443:                 CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, DUMQ, Q(1:DEGFREEDOMS))
1536:             ENDIF1444:             ENDIF
 1445: 
1537:             IF (BHINTERPT.AND.(.NOT.REOPTIMISEENDPOINTS)) ENERGY=BHENERGY1446:             IF (BHINTERPT.AND.(.NOT.REOPTIMISEENDPOINTS)) ENERGY=BHENERGY
1538:             IF (BISECTT.AND.(.NOT.REOPTIMISEENDPOINTS)) ENERGY=BISECTENERGY1447:             IF (BISECTT.AND.(.NOT.REOPTIMISEENDPOINTS)) ENERGY=BISECTENERGY
1539:             IF (LOWESTFRQT) THEN1448:             IF (LOWESTFRQT) THEN
1540:                IF (CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T) THEN1449:                IF (CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T) THEN
1541:                   WRITE(881,'(2F20.10,I6,5F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ,MINCURVE,MINFRQ21450:                   WRITE(881,'(2F20.10,I6,5F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ,MINCURVE,MINFRQ2
1542:                ELSE1451:                ELSE
1543:                   CALL INERTIA2(Q,ITX,ITY,ITZ)1452:                   CALL INERTIA2(Q,ITX,ITY,ITZ)
1544:                   WRITE(881,'(2F20.10,I6,5F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ,MINCURVE,MINFRQ21453:                   WRITE(881,'(2F20.10,I6,5F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ,MINCURVE,MINFRQ2
1545:                ENDIF1454:                ENDIF
1546:             ELSE   1455:             ELSE   
1547:                IF (CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T) THEN1456:                IF (CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T) THEN
1548:                   WRITE(881,'(2F20.10,I6,4F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ1457:                   WRITE(881,'(2F20.10,I6,4F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ
1549:                ELSE1458:                ELSE
1550:                   IF(RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN1459:                   IF(RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN
1551:                       CALL INERTIA2(DUMQ,ITX,ITY,ITZ)1460:                       CALL INERTIA2(DUMQ,ITX,ITY,ITZ)
1552:                   ELSE IF (GTHOMSONT) THEN1461:                   ELSE IF (GTHOMSONT) THEN
1553:                      CALL GTHOMSONANGTOC(TMPC(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)1462:                      CALL GTHOMSONANGTOC(TMPC(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)
1554:                      CALL INERTIA2(TMPC,ITX,ITY,ITZ)1463:                      CALL INERTIA2(TMPC,ITX,ITY,ITZ)
1555:                   ENDIF1464:                   ENDIF
1556:  
1557:                   ! Write the header line to min.data.info. PROD should be 1 if NOFRQS is set, or a positive real value 
1558:                   ! if ENDHESS is set. If neither of these is set, PROD is probably undefined (for some reason we don't 
1559:                   ! calculate it in the MFLAG section - I'll look into adding it) -sn402 
1560:                   WRITE(881,'(2F20.10,I6,4F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ1465:                   WRITE(881,'(2F20.10,I6,4F20.10)') ENERGY,PROD,HORDER,ITX,ITY,ITZ
1561:                ENDIF1466:                ENDIF
1562:             ENDIF1467:             ENDIF
1563:             NRES=NMRES1468:             NRES=NMRES
1564: 1469: 
1565:             IF (AMHT) THEN1470:             IF (AMHT) THEN
1566:                GLY_COUNT = 01471:                GLY_COUNT = 0
1567:                SDUMMY='AM'1472:                SDUMMY='AM'
1568:                DO J2=1,NRES1473:                DO J2=1,NRES
1569:                  IF (SEQ(J2).EQ.8) THEN1474:                  IF (SEQ(J2).EQ.8) THEN
1574:                    GLY_COUNT = GLY_COUNT +11479:                    GLY_COUNT = GLY_COUNT +1
1575:                  ELSE1480:                  ELSE
1576:                      WRITE(881,*)Q(9*(J2-1)+1-GLY_COUNT*3),Q(9*(J2-1)+2-GLY_COUNT*3),Q(9*(J2-1)+3-GLY_COUNT*3)   1481:                      WRITE(881,*)Q(9*(J2-1)+1-GLY_COUNT*3),Q(9*(J2-1)+2-GLY_COUNT*3),Q(9*(J2-1)+3-GLY_COUNT*3)   
1577:                      WRITE(881,*)Q(9*(J2-1)+4-GLY_COUNT*3),Q(9*(J2-1)+5-GLY_COUNT*3),Q(9*(J2-1)+6-GLY_COUNT*3)1482:                      WRITE(881,*)Q(9*(J2-1)+4-GLY_COUNT*3),Q(9*(J2-1)+5-GLY_COUNT*3),Q(9*(J2-1)+6-GLY_COUNT*3)
1578:                      WRITE(881,*)Q(9*(J2-1)+7-GLY_COUNT*3),Q(9*(J2-1)+8-GLY_COUNT*3),Q(9*(J2-1)+9-GLY_COUNT*3)1483:                      WRITE(881,*)Q(9*(J2-1)+7-GLY_COUNT*3),Q(9*(J2-1)+8-GLY_COUNT*3),Q(9*(J2-1)+9-GLY_COUNT*3)
1579:                  ENDIF1484:                  ENDIF
1580: 1485: 
1581:                ENDDO1486:                ENDDO
1582:                CALL FLUSH(881)1487:                CALL FLUSH(881)
1583:             ELSE1488:             ELSE
1584:  
1585:                ! Now we write the coordinates, which is the main part of the min.data.info file.  
1586: ! hk2861489: ! hk286
1587:                IF (RIGIDINIT) THEN1490:                IF (RIGIDINIT) THEN
1588:                   CALL TRANSFORMRIGIDTOC (1, NRIGIDBODY, XCOORDS, Q(1:DEGFREEDOMS))1491:                   CALL TRANSFORMRIGIDTOC (1, NRIGIDBODY, XCOORDS, Q(1:DEGFREEDOMS))
1589:                   WRITE(881,'(3F25.15)') XCOORDS(1:3*NATOMS)1492:                   WRITE(881,'(3F25.15)') XCOORDS(1:3*NATOMS)
1590: 1493: 
1591:                ELSE IF (GTHOMSONT) THEN1494:                ELSE IF (GTHOMSONT) THEN
1592:                   CALL GTHOMSONANGTOC(TMPC(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)1495:                   CALL GTHOMSONANGTOC(TMPC(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)
1593:                   WRITE(881,'(3F25.15)') TMPC(1:3*NATOMS)1496:                   WRITE(881,'(3F25.15)') TMPC(1:3*NATOMS)
1594:                ELSE IF (VARIABLES) THEN1497:                ELSE IF (VARIABLES) THEN
1595:                   WRITE(881,'(F25.15)') Q(1:NATOMS)1498:                   WRITE(881,'(F25.15)') Q(1:NATOMS)
1598:                ENDIF1501:                ENDIF
1599:                CALL FLUSH(881)1502:                CALL FLUSH(881)
1600:             ENDIF1503:             ENDIF
1601:          ELSE1504:          ELSE
1602:             PRINT '(A)',' geopt> WARNING - DUMPDATA is set, but MFLAG is false; not creating entry in min.data.info'1505:             PRINT '(A)',' geopt> WARNING - DUMPDATA is set, but MFLAG is false; not creating entry in min.data.info'
1603:          ENDIF1506:          ENDIF
1604:       ENDIF1507:       ENDIF
1605: 1508: 
1606:       RETURN1509:       RETURN
1607: 1510: 
1608:       END  ! END GEOPT!1511:       END
1609: 1512: 
1610:       SUBROUTINE MASSWT(NATOMS,ATMASS,Q,VNEW,STEST)1513:       SUBROUTINE MASSWT(NATOMS,ATMASS,Q,VNEW,STEST)
1611:       USE MODHESS1514:       USE MODHESS
1612:       IMPLICIT NONE1515:       IMPLICIT NONE
1613:       INTEGER NATOMS, J1, J2, J3, J41516:       INTEGER NATOMS, J1, J2, J3, J4
1614:       LOGICAL STEST1517:       LOGICAL STEST
1615:       DOUBLE PRECISION ATMASS(NATOMS),Q(3*NATOMS),VNEW(3*NATOMS),AMASS,BMASS,PMASS1518:       DOUBLE PRECISION ATMASS(NATOMS),Q(3*NATOMS),VNEW(3*NATOMS),AMASS,BMASS,PMASS
1616: 1519: 
1617:       IF (SIZE(HESS,2).LT.3*NATOMS) RETURN ! for example, for 2D XY model1520:       IF (SIZE(HESS,2).LT.3*NATOMS) RETURN ! for example, for 2D XY model
1618:       IF (.NOT.ALLOCATED(HESS)) THEN1521:       IF (.NOT.ALLOCATED(HESS)) THEN


r31191/rigidb.f90 2016-09-23 11:30:17.474816663 +0100 r31190/rigidb.f90 2016-09-23 11:30:18.482830098 +0100
2124:             AP(I,J) = KD(I)*AP(I,J)*KD(J)2124:             AP(I,J) = KD(I)*AP(I,J)*KD(J)
2125:          ENDDO2125:          ENDDO
2126:       ENDDO2126:       ENDDO
2127: 2127: 
2128:       IF (EIGENVECTORT) THEN2128:       IF (EIGENVECTORT) THEN
2129:          CALL DSYEV('V','L',NDIM,AP,NDIM,FRQN,WORK,LWORK,INFO)2129:          CALL DSYEV('V','L',NDIM,AP,NDIM,FRQN,WORK,LWORK,INFO)
2130:       ELSE2130:       ELSE
2131:          CALL DSYEV('N','L',NDIM,AP,NDIM,FRQN,WORK,LWORK,INFO)2131:          CALL DSYEV('N','L',NDIM,AP,NDIM,FRQN,WORK,LWORK,INFO)
2132:       ENDIF2132:       ENDIF
2133: 2133: 
2134:       ! sn402: commented this out: it doesn't play well with geopt, because this subroutine actually sorts2134:       call eigensort_val_asc(FRQN,AP,NDIM,3*NATOMS)
2135:       ! eigenvalues into descending order but geopt expects ascending order. 
2136: !      call eigensort_val_asc(FRQN,AP,NDIM,3*NATOMS) 
2137:  
2138: !      DO I = 1, NDIM2135: !      DO I = 1, NDIM
2139: !         IF (FRQN(I) > 0.0D0) THEN2136: !         IF (FRQN(I) > 0.0D0) THEN
2140: !            FRQN(I) = FRQCNV * SQRT((FRQN(I)))2137: !            FRQN(I) = FRQCNV * SQRT((FRQN(I)))
2141: !         ELSE2138: !         ELSE
2142: !            FRQN(I) = -FRQCNV * SQRT((-FRQN(I)))2139: !            FRQN(I) = -FRQCNV * SQRT((-FRQN(I)))
2143: !         ENDIF2140: !         ENDIF
2144: !      ENDDO2141: !      ENDDO
2145: !      FRQN(:) = FRQCNV*SQRT(ABS(FRQN(:)))2142: !      FRQN(:) = FRQCNV*SQRT(ABS(FRQN(:)))
2146:  
2147:       ! sn402: This conversion factor is not general and unexplained, so I'm getting rid of it for now. 
2148:       ! It is presumably converting squared frequencies in internal units (kJmol^-1 amu^-1 Angs^-2) to s^-2. 
2149:       ! But this is not the right place for unit conversions. I'm going to implement a general solution to the 
2150:       ! unit problem in the near future. 
2151:       IF(.FALSE.) THEN 
2152:       FRQN(:) = FRQN(:) * 1.0D262143:       FRQN(:) = FRQN(:) * 1.0D26
2153:       ENDIF 
2154: !      print *, 'FRQN'2144: !      print *, 'FRQN'
2155: !      print *, FRQN2145: !      print *, FRQN
2156: 2146: 
2157: !      PRINT *, "CALLED FREQ"2147: !      PRINT *, "CALLED FREQ"
2158: !      PRINT *, FRQN2148: !      PRINT *, FRQN
2159: 2149: 
2160:       IF (EIGENVECTORT) THEN      2150:       IF (EIGENVECTORT) THEN      
2161:          HESS = AP2151:          HESS = AP
2162:       ENDIF2152:       ENDIF
2163: 2153: 


r31191/vdump.f 2016-09-23 11:30:17.722819971 +0100 r31190/vdump.f 2016-09-23 11:30:18.738833511 +0100
 12: C   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: C   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: C   GNU General Public License for more details. 13: C   GNU General Public License for more details.
 14: C 14: C
 15: C   You should have received a copy of the GNU General Public License 15: C   You should have received a copy of the GNU General Public License
 16: C   along with this program; if not, write to the Free Software 16: C   along with this program; if not, write to the Free Software
 17: C   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: C   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: C 18: C
 19:       SUBROUTINE VDUMP(DIAG,ZT,N,M) 19:       SUBROUTINE VDUMP(DIAG,ZT,N,M)
 20:       USE KEY 20:       USE KEY
 21:       USE MODHESS 21:       USE MODHESS
 22:       USE MODCHARMM, ONLY: CHRMMT 
 23:       USE PORFUNCS 22:       USE PORFUNCS
 24:       IMPLICIT NONE 23:       IMPLICIT NONE
 25:       INTEGER M, N, J1, J2, ISTAT, MCOUNT 24:       INTEGER M, N, J1, J2, ISTAT, MCOUNT
 26:       DOUBLE PRECISION DIAG(M) 25:       DOUBLE PRECISION DIAG(M)
 27:       LOGICAL ZT(M) 26:       LOGICAL ZT(M)
 28:    
 29: C 27: C
 30: C  dump the eigenvectors which correspond to non-zero eigenvalues 28: C  dump the eigenvectors which correspond to non-zero eigenvalues
 31: C  in file vectors.dump 29: C  in file vectors.dump
 32: C 30: C
 33:       IF (.NOT.ALLSTEPS) REWIND(44) 31:       IF (.NOT.ALLSTEPS) REWIND(44)
 34:       IF (ALLVECTORS) THEN 32:       IF (ALLVECTORS) THEN
 35: !       csw34> Find how many modes are going to be printed and echo this 33: !       csw34> Find how many modes are going to be printed and echo this
 36: !       to a file for the analysis script (this could be removed when all the 34: !       to a file for the analysis script (this could be removed when all the
 37: !       dumping is done in OPTIM). This is a problem with FREEZE where 35: !       dumping is done in OPTIM). This is a problem with FREEZE where
 38: !       some non-real modes have zero eiganvalue. 36: !       some non-real modes have zero eiganvalue.
 41:                 IF (ZT(J1)) MCOUNT=MCOUNT+1 39:                 IF (ZT(J1)) MCOUNT=MCOUNT+1
 42:             ENDDO 40:             ENDDO
 43:             OPEN(UNIT=499,FILE='nmodes.dat',STATUS='UNKNOWN') 41:             OPEN(UNIT=499,FILE='nmodes.dat',STATUS='UNKNOWN')
 44:             WRITE(499,'(I6)') MCOUNT 42:             WRITE(499,'(I6)') MCOUNT
 45:             CLOSE(499) 43:             CLOSE(499)
 46:          DO J1=1,N 44:          DO J1=1,N
 47:             IF (ZT(J1)) THEN 45:             IF (ZT(J1)) THEN
 48: ! If printing the mass weighted vectors (normal modes), convert omega^2 46: ! If printing the mass weighted vectors (normal modes), convert omega^2
 49: ! into the vibrational frequency in wavenumbers (cm^(-1)). 108.52 is the 47: ! into the vibrational frequency in wavenumbers (cm^(-1)). 108.52 is the
 50: ! conversion factor from (kcal mol-1 kg-1 A-2)^2 to cm-1. 48: ! conversion factor from (kcal mol-1 kg-1 A-2)^2 to cm-1.
 51:               IF (MWVECTORS .AND. (AMBERT .OR. AMBER12T .OR. CHRMMT)) THEN 49:               IF (MWVECTORS) THEN
 52:                         WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.52 50:                         WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.52
 53:               ELSE 51:               ELSE
 54:                         WRITE(44,'(F20.10)') DIAG(J1) 52:                         WRITE(44,'(F20.10)') DIAG(J1)
 55:               ENDIF 53:               ENDIF
 56: !              IF (ANGLEAXIS2) N=N/2 54: !              IF (ANGLEAXIS2) N=N/2
 57:                WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N) 55:                WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N)
 58: !              IF (ANGLEAXIS2) N=N*2 56: !              IF (ANGLEAXIS2) N=N*2
 59:             ENDIF 57:             ENDIF
 60:          ENDDO 58:          ENDDO
 61:       ELSE 59:       ELSE
 62:          DO J1=N,1,-1 60:          DO J1=N,1,-1
 63:             IF (ZT(J1)) THEN 61:             IF (ZT(J1)) THEN
 64: ! As above 62: ! As above
 65:                IF (MWVECTORS .AND. (AMBERT .OR. AMBER12T .OR. CHRMMT)) THEN 63:                IF (MWVECTORS) THEN
 66:                         WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.52 64:                         WRITE(44,'(F20.10)') DSQRT(DIAG(J1))*108.52
 67:                ELSE 65:                ELSE
 68:                         WRITE(44,'(F20.10)') DIAG(J1) 66:                         WRITE(44,'(F20.10)') DIAG(J1)
 69:                ENDIF 67:                ENDIF
 70:                WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N) 68:                WRITE(44,'(3F20.10)') (HESS(J2,J1),J2=1,N)
 71:                CALL FLUSH(44) 69:                CALL FLUSH(44)
 72:                RETURN 70:                RETURN
 73:             ENDIF 71:             ENDIF
 74:          ENDDO 72:          ENDDO
 75:       ENDIF 73:       ENDIF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0