hdiff output

r31190/genrigidtip.f90 2016-09-23 11:30:15.878795387 +0100 r31189/genrigidtip.f90 2016-09-23 11:30:16.398802319 +0100
 32:  32: 
 33:       STCHRG(:) = (/0.D0, 0.52D0, 0.52D0, -1.04D0/) 33:       STCHRG(:) = (/0.D0, 0.52D0, 0.52D0, -1.04D0/)
 34:       C6        = 2552.24D0 ! LJ coefficients in kJ/mol Angstrom**6 or Angstrom**12 34:       C6        = 2552.24D0 ! LJ coefficients in kJ/mol Angstrom**6 or Angstrom**12
 35:       C12       = 2510.4D3 35:       C12       = 2510.4D3
 36:       CH2O      = 1389.354848D0 ! Conversion factor for coulomb energy 36:       CH2O      = 1389.354848D0 ! Conversion factor for coulomb energy
 37:  37: 
 38:       ENERGY  = 0.D0 38:       ENERGY  = 0.D0
 39:       IF (GTEST) G(:) = 0.D0 39:       IF (GTEST) G(:) = 0.D0
 40:       IF (STEST) HESS(:,:) = 0.D0 40:       IF (STEST) HESS(:,:) = 0.D0
 41:  41: 
 42:       ! Loop over pairs of different rigid bodies - exclude intramolecular interactions 42: 
 43:       DO J1 = 1, NRIGIDBODY 43:       DO J1 = 1, NRIGIDBODY
 44:  44: 
 45:          J3 = 3*(J1-1)*SITESPERBODY 45:          J3 = 3*(J1-1)*SITESPERBODY
 46:  46: 
 47:          DO J2 = J1 + 1, NRIGIDBODY 47:          DO J2 = J1 + 1, NRIGIDBODY
 48:  48: 
 49:             J4 = 3*(J2-1)*SITESPERBODY 49:             J4 = 3*(J2-1)*SITESPERBODY
 50:  50: 
 51: !     O-O LJ CONTRIBUTION 51: !     O-O LJ CONTRIBUTION
 52:  52: 
 66:  66: 
 67:             IF (GTEST .OR. STEST) THEN 67:             IF (GTEST .OR. STEST) THEN
 68: !     DVDR = DVDR/R 68: !     DVDR = DVDR/R
 69: !               DVDR(J1,J2) =-6.D0*(2.D0*C12*R12 - C6*R6)*R2 69: !               DVDR(J1,J2) =-6.D0*(2.D0*C12*R12 - C6*R6)*R2
 70: !               DVDR(J2,J1) = DVDR(J1,J2) 70: !               DVDR(J2,J1) = DVDR(J1,J2)
 71:                DVDR(J5,J6) =-6.D0*(2.D0*C12*R12 - C6*R6)*R2 71:                DVDR(J5,J6) =-6.D0*(2.D0*C12*R12 - C6*R6)*R2
 72:                DVDR(J6,J5) = DVDR(J5,J6) 72:                DVDR(J6,J5) = DVDR(J5,J6)
 73:  73: 
 74: !               G(J7:J7+2) = G(J7:J7+2) + DVDR(J1,J2)*RSS(:) 74: !               G(J7:J7+2) = G(J7:J7+2) + DVDR(J1,J2)*RSS(:)
 75: !               G(J8:J8+2) = G(J8:J8+2) - DVDR(J1,J2)*RSS(:) 75: !               G(J8:J8+2) = G(J8:J8+2) - DVDR(J1,J2)*RSS(:)
 76:                 G(J7:J7+2) = G(J7:J7+2) + DVDR(J5,J6)*RSS(:) 76:                G(J7:J7+2) = G(J7:J7+2) + DVDR(J5,J6)*RSS(:)
 77:                 G(J8:J8+2) = G(J8:J8+2) - DVDR(J5,J6)*RSS(:) 77:                G(J8:J8+2) = G(J8:J8+2) - DVDR(J5,J6)*RSS(:)
 78:  78: 
 79: !               D2VDR2(J1,J2) = (168.D0*C12*R12 - 48.D0*C6*R6)*R2*R2 79: !               D2VDR2(J1,J2) = (168.D0*C12*R12 - 48.D0*C6*R6)*R2*R2
 80:                D2VDR2(J5,J6) = (168.D0*C12*R12 - 48.D0*C6*R6)*R2*R2 80:                D2VDR2(J5,J6) = (168.D0*C12*R12 - 48.D0*C6*R6)*R2*R2
 81:             ENDIF 81:             ENDIF
 82:  82: 
 83:             IF (STEST) THEN 83:             IF (STEST) THEN
 84: !                D2VDR2(J2,J1) = D2VDR2(J1,J2) 84: !                D2VDR2(J2,J1) = D2VDR2(J1,J2)
 85:                 D2VDR2(J6,J5) = D2VDR2(J5,J6) 85:                 D2VDR2(J6,J5) = D2VDR2(J5,J6)
 86:             ENDIF 86:             ENDIF
 87:  87: 
198:                DUMMY           = - D2VDR2(J1,J2)*RSS(3)*RSS(1)198:                DUMMY           = - D2VDR2(J1,J2)*RSS(3)*RSS(1)
199:                HESS(J3,J4-2)   = HESS(J3,J4-2) + DUMMY199:                HESS(J3,J4-2)   = HESS(J3,J4-2) + DUMMY
200:                HESS(J4-2,J3)   = HESS(J4-2,J3) + DUMMY200:                HESS(J4-2,J3)   = HESS(J4-2,J3) + DUMMY
201: 201: 
202:             ENDDO202:             ENDDO
203: 203: 
204:          ENDDO204:          ENDDO
205: 205: 
206:       ENDIF206:       ENDIF
207: 207: 
208: !write(*,*) "2nd derivative matrix" 
209: !DO J1=1,UBOUND(D2VDR2,1) 
210: !   DO J2=1,UBOUND(D2VDR2,2) 
211: !      write(*,*) J1, J2, D2VDR2(J1,J2) 
212: !   ENDDO 
213: !ENDDO 
214: !   STOP 
215:  
216: ! sn402: remove208: ! sn402: remove
217: !    IF(DEBUG .AND. STEST) THEN209:     IF(DEBUG .AND. STEST) THEN
218: !      DO J1 = 1,3*NATOMS210:       DO J1 = 1,3*NATOMS
219: !        DO J2 = 1,3*NATOMS211:         DO J2 = 1,3*NATOMS
220: !            IF(ABS(HESS(J1,J2)-HESS(J2,J1)) .GT. 1.0E-7) THEN212:             IF(ABS(HESS(J1,J2)-HESS(J2,J1)) .GT. 1.0E-7) THEN
221: !                write(*,*) "genrigidtip> Asymmetric Hessian, coords ", J1, J2213:                 write(*,*) "genrigidtip> Asymmetric Hessian, coords ", J1, J2
222: !            ENDIF214:             ENDIF
223: !        ENDDO215:         ENDDO
224: !      ENDDO216:       ENDDO
225: !      write(*,*) "Cartesian HESSIAN:"217: !      write(*,*) "Cartesian HESSIAN:"
226: !      DO J1=1,3*NATOMS218: !      DO J1=1,3*NATOMS
227: !          write(*,*) HESS(J1,:)219: !          write(*,*) HESS(J1,:)
228: !     ENDDO220: !     ENDDO
229: !    ENDIF221:     ENDIF
230: 222: 
231:       END SUBROUTINE GENRIGIDTIP223:       END SUBROUTINE GENRIGIDTIP
232: 224: 
233: !     ---------------------------------------------------------------------------------------------225: !     ---------------------------------------------------------------------------------------------


r31190/newtip.f90 2016-09-23 11:30:16.126798691 +0100 r31189/newtip.f90 2016-09-23 11:30:16.658805783 +0100
 30:       D2VDR2(:,:) = 0.D0; DVDR(:,:) = 0.D0 30:       D2VDR2(:,:) = 0.D0; DVDR(:,:) = 0.D0
 31:       DOTI1(:,:) = 0.D0; DOTI2(:,:) = 0.D0; DOTI3(:,:) = 0.D0 31:       DOTI1(:,:) = 0.D0; DOTI2(:,:) = 0.D0; DOTI3(:,:) = 0.D0
 32:       DOTJ1(:,:) = 0.D0; DOTJ2(:,:) = 0.D0; DOTJ3(:,:) = 0.D0 32:       DOTJ1(:,:) = 0.D0; DOTJ2(:,:) = 0.D0; DOTJ3(:,:) = 0.D0
 33:  33: 
 34:       CALL DEFTIP4(C12, C6, CH2O) 34:       CALL DEFTIP4(C12, C6, CH2O)
 35:   35:  
 36:       ENERGY  = 0.D0 36:       ENERGY  = 0.D0
 37:       IF (GTEST) G(:) = 0.D0 37:       IF (GTEST) G(:) = 0.D0
 38:       IF (STEST) HESS(:,:) = 0.D0 38:       IF (STEST) HESS(:,:) = 0.D0
 39:  39: 
 40:       ! REALNATOMS is the number of TIP4P molecules (=NOPT/6 because 3xCoM+3xAA degrees of freedom per molecule) 
 41:       ! NOPT = 3*NATOMS (number of degrees of freedom for an atomic system), so REALNATOMS=NATOMS/2 
 42:       REALNATOMS = NATOMS/2 40:       REALNATOMS = NATOMS/2
 43:       ! OFFSET is the number of CoM coords. 
 44:       OFFSET     = 3*REALNATOMS 41:       OFFSET     = 3*REALNATOMS
 45:   42:  
 46:       ! The potential is defined by pairwise site-site interactions. 
 47:       ! In order to compute it, we need to obtain cartesian coordinates for the system. 
 48:       ! Do this here. 
 49:       DO J1 = 1, REALNATOMS 43:       DO J1 = 1, REALNATOMS
 50:  44: 
 51:          J3 = 3*J1 45:          J3 = 3*J1
 52:          J5 = OFFSET + J3 46:          J5 = OFFSET + J3
 53:          RI = X(J3-2:J3) ! CoM coords for molecule J1 47:          RI = X(J3-2:J3)
 54:          P  = X(J5-2:J5) ! AA coords for molecule J1 48:          P  = X(J5-2:J5)
 55:  49: 
 56:          CALL RMDFAS(P, RMI, DRMI1, DRMI2, DRMI3, D2RMI1, D2RMI2, D2RMI3, D2RMI12, D2RMI23, D2RMI31, GTEST, STEST) 50:          CALL RMDFAS(P, RMI, DRMI1, DRMI2, DRMI3, D2RMI1, D2RMI2, D2RMI3, D2RMI12, D2RMI23, D2RMI31, GTEST, STEST)
 57:  51: 
 58: ! hk286 52: ! hk286
 59:          IF ( RBAANORMALMODET ) THEN 53:          IF ( RBAANORMALMODET ) THEN
 60:             P = (/0.0D0, 0.0D0, 0.0D0/) 54:             P = (/0.0D0, 0.0D0, 0.0D0/)
 61:             CALL RMDFAS(P, RMI0, DRMI10, DRMI20, DRMI30, D2RMI10, D2RMI20, D2RMI30, D2RMI120, D2RMI230, D2RMI310, GTEST, STEST) 55:             CALL RMDFAS(P, RMI0, DRMI10, DRMI20, DRMI30, D2RMI10, D2RMI20, D2RMI30, D2RMI120, D2RMI230, D2RMI310, GTEST, STEST)
 62:          ENDIF 56:          ENDIF
 63:  57: 
 64:          ! Loop over sites (i.e. atoms) in this molecule 
 65:          DO J2 = 1, NRBSITES 58:          DO J2 = 1, NRBSITES
 66:  59: 
 67:             ! J4 is an index for this site relative to a complete list of sites in all molecules 
 68:             J4        = NRBSITES*(J1-1) + J2 60:             J4        = NRBSITES*(J1-1) + J2
 69:             ! R(J4,:) contains the cartesian coordinates for atom J4 
 70:             R(J4,:)   = RI(:) + MATMUL(RMI,RBSITE(J2,:)) 61:             R(J4,:)   = RI(:) + MATMUL(RMI,RBSITE(J2,:))
 71:  62: 
 72:             IF (GTEST .OR. STEST) THEN 63:             IF (GTEST .OR. STEST) THEN
 73:  64: 
 74: ! hk286 65: ! hk286
 75:                IF ( RBAANORMALMODET ) THEN 66:                IF ( RBAANORMALMODET ) THEN
 76:                   DR1(J4,:) = MATMUL(DRMI10,MATMUL(RMI,RBSITE(J2,:))) 67:                   DR1(J4,:) = MATMUL(DRMI10,MATMUL(RMI,RBSITE(J2,:)))
 77:                   DR2(J4,:) = MATMUL(DRMI20,MATMUL(RMI,RBSITE(J2,:))) 68:                   DR2(J4,:) = MATMUL(DRMI20,MATMUL(RMI,RBSITE(J2,:)))
 78:                   DR3(J4,:) = MATMUL(DRMI30,MATMUL(RMI,RBSITE(J2,:))) 69:                   DR3(J4,:) = MATMUL(DRMI30,MATMUL(RMI,RBSITE(J2,:)))
 79:                ELSE 70:                ELSE
104:                   D2R23(J4,:) = MATMUL(D2RMI23,RBSITE(J2,:)) 95:                   D2R23(J4,:) = MATMUL(D2RMI23,RBSITE(J2,:))
105:                   D2R31(J4,:) = MATMUL(D2RMI31,RBSITE(J2,:)) 96:                   D2R31(J4,:) = MATMUL(D2RMI31,RBSITE(J2,:))
106:                ENDIF 97:                ENDIF
107:  98: 
108:             ENDIF 99:             ENDIF
109: 100: 
110:          ENDDO101:          ENDDO
111: 102: 
112:       ENDDO103:       ENDDO
113: 104: 
114: !            write(*,*) "Cartesian Coordinates of the system" 
115: !            do j6=1,UBOUND(R,1) 
116: !               write(*,*) R(j6,:) 
117: !            enddo 
118: !            stop 
119:  
120:       ! Now compute the actual potential and derivatives (analogous to the G and F tensors in ljpshift, for example) 
121:       ! Loop over molecules 
122:       DO J1 = 1, REALNATOMS 105:       DO J1 = 1, REALNATOMS 
123: 106: 
124:          ! Base index for the CoM coords 
125:          J3 = 3*J1107:          J3 = 3*J1
126:          ! Base index for the AA coords 
127:          J5 = OFFSET + J3108:          J5 = OFFSET + J3
128: 109: 
129:          ! Loop over other molecules - we only compute interactions between different molecules. 
130:          DO J2 = J1 + 1, REALNATOMS110:          DO J2 = J1 + 1, REALNATOMS
131: 111: 
132:             J4 = 3*J2112:             J4 = 3*J2
133:             J6 = OFFSET + J4113:             J6 = OFFSET + J4
134: 114: 
135: !     O-O LJ CONTRIBUTION115: !     O-O LJ CONTRIBUTION
136: 116: 
137:             ! J7 and J8 are site indices for the O atoms on the two different molecules 
138:             J7 = NRBSITES*(J1-1) + 1117:             J7 = NRBSITES*(J1-1) + 1
139:             J8 = NRBSITES*(J2-1) + 1118:             J8 = NRBSITES*(J2-1) + 1
140:             ! Compute the displacement between the O atoms 
141:             RSS(:) = R(J7,:) - R(J8,:)119:             RSS(:) = R(J7,:) - R(J8,:)
142:             ! Construct the LJ energy term 
143:             RSS2   = DOT_PRODUCT(RSS(:),RSS(:))120:             RSS2   = DOT_PRODUCT(RSS(:),RSS(:))
144:             ABSR   = DSQRT(RSS2)121:             ABSR   = DSQRT(RSS2)
145:             R2     = 1.D0/RSS2 ! Note, from now on the R2, R6, R12 terms are all inverse lengths122:             R2     = 1.D0/RSS2
146:             R6     = R2*R2*R2123:             R6     = R2*R2*R2
147:             R12    = R6*R6124:             R12    = R6*R6
148:             ENERGY = ENERGY + C12*R12 - C6*R6125:             ENERGY = ENERGY + C12*R12 - C6*R6
149: 126: 
150:             IF (GTEST .OR. STEST) THEN127:             IF (GTEST .OR. STEST) THEN
151: !     DVDR = DVDR/R  (referred to as the G tensor in some other potential functions)128: !     DVDR = DVDR/R
152:                DVDR(J7,J8) =-6.D0*(2.D0*C12*R12 - C6*R6)*R2129:                DVDR(J7,J8) =-6.D0*(2.D0*C12*R12 - C6*R6)*R2
153:                DVDR(J8,J7) = DVDR(J7,J8)130:                DVDR(J8,J7) = DVDR(J7,J8)
154: !               write(*,*) "Atoms ", J7, J8131: !               write(*,*) "Atoms ", J7, J8
155: !               write(*,*) "Derivative ", DVDR(J7, J8)132: !               write(*,*) "Derivative ", DVDR(J7, J8)
156: 133: 
157:                DOTI1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J7,:))134:                DOTI1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J7,:))
158:                DOTI2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J7,:))135:                DOTI2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J7,:))
159:                DOTI3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J7,:))136:                DOTI3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J7,:))
160: 137: 
161:                DOTJ1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J8,:))138:                DOTJ1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J8,:))
162:                DOTJ2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J8,:))139:                DOTJ2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J8,:))
163:                DOTJ3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J8,:))140:                DOTJ3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J8,:))
164: 141: !write(*,*) "Gradient for atoms ", J7, J8
165:                ! Calculate the gradient here.142: !write(*,*) DVDR(J7,J8)*RSS(:)
166:                ! CoM terms are quite simple, just add up the site terms (J7,J8) 
167:                ! that belong to each molecule (J3,J4) 
168:                G(J3-2:J3) = G(J3-2:J3) + DVDR(J7,J8)*RSS(:)143:                G(J3-2:J3) = G(J3-2:J3) + DVDR(J7,J8)*RSS(:)
169:                G(J4-2:J4) = G(J4-2:J4) - DVDR(J7,J8)*RSS(:)144:                G(J4-2:J4) = G(J4-2:J4) - DVDR(J7,J8)*RSS(:)
170: 145: 
171:                ! AA terms are more complicated to obtain from cartesian terms. 
172:                G(J5-2) = G(J5-2) + DVDR(J7,J8)*DOTI1(J7,J8)146:                G(J5-2) = G(J5-2) + DVDR(J7,J8)*DOTI1(J7,J8)
173:                G(J5-1) = G(J5-1) + DVDR(J7,J8)*DOTI2(J7,J8)147:                G(J5-1) = G(J5-1) + DVDR(J7,J8)*DOTI2(J7,J8)
174:                G(J5)   = G(J5)   + DVDR(J7,J8)*DOTI3(J7,J8)148:                G(J5)   = G(J5)   + DVDR(J7,J8)*DOTI3(J7,J8)
175: 149: 
176:                G(J6-2) = G(J6-2) - DVDR(J7,J8)*DOTJ1(J7,J8)150:                G(J6-2) = G(J6-2) - DVDR(J7,J8)*DOTJ1(J7,J8)
177:                G(J6-1) = G(J6-1) - DVDR(J7,J8)*DOTJ2(J7,J8)151:                G(J6-1) = G(J6-1) - DVDR(J7,J8)*DOTJ2(J7,J8)
178:                G(J6)   = G(J6)   - DVDR(J7,J8)*DOTJ3(J7,J8)152:                G(J6)   = G(J6)   - DVDR(J7,J8)*DOTJ3(J7,J8)
179: 153: 
180:                ! D2VDR2 = 1/r d/dr(1/r dV/dr) = 1/r d/dr(DVDR) is different from the  
181:                ! F tensor referred to in the LJ potential, for instance. 
182:                D2VDR2(J7,J8) = (168.D0*C12*R12 - 48.D0*C6*R6)*R2*R2154:                D2VDR2(J7,J8) = (168.D0*C12*R12 - 48.D0*C6*R6)*R2*R2
 155: !               write(*,*) "D2VDR2", J7, J8
 156: !               write(*,*) D2VDR2(J7,J8)
183:             ENDIF157:             ENDIF
184: 158: 
185:             IF (STEST) THEN159:             IF (STEST) THEN
186: 160: 
187:                 D2VDR2(J8,J7) = D2VDR2(J7,J8)161:                 D2VDR2(J8,J7) = D2VDR2(J7,J8)
188:                 DOTI1(J8,J7)  =-DOTJ1(J7,J8)162:                 DOTI1(J8,J7)  =-DOTJ1(J7,J8)
189:                 DOTI2(J8,J7)  =-DOTJ2(J7,J8)163:                 DOTI2(J8,J7)  =-DOTJ2(J7,J8)
190:                 DOTI3(J8,J7)  =-DOTJ3(J7,J8)164:                 DOTI3(J8,J7)  =-DOTJ3(J7,J8)
191:                 DOTJ1(J8,J7)  =-DOTI1(J7,J8)165:                 DOTJ1(J8,J7)  =-DOTI1(J7,J8)
192:                 DOTJ2(J8,J7)  =-DOTI2(J7,J8)166:                 DOTJ2(J8,J7)  =-DOTI2(J7,J8)
193:                 DOTJ3(J8,J7)  =-DOTI3(J7,J8)167:                 DOTJ3(J8,J7)  =-DOTI3(J7,J8)
194: 168: 
195:             ENDIF169:             ENDIF
196: 170: 
197:             ! Now compute the coulomb part of the interaction and derivative tensors. 
198:             ! Because there are multiple coulomb sites on each molecule, we need a loop over these sites. 
199:             DO I = 2, NRBSITES 171:             DO I = 2, NRBSITES 
200:                ! Index for this particular interaction site (either an H atom or the dummy charge)172: 
201:                J7 = NRBSITES*(J1-1) + I173:                J7 = NRBSITES*(J1-1) + I
202: 174: 
203:                DO J = 2, NRBSITES 175:                DO J = 2, NRBSITES 
204: 176: 
205:                   J8     = NRBSITES*(J2-1) + J177:                   J8     = NRBSITES*(J2-1) + J
206:                   ! Site-site displacement vector 
207:                   RSS(:) = R(J7,:) - R(J8,:)178:                   RSS(:) = R(J7,:) - R(J8,:)
208:                   ! Compute the energy 
209:                   RSS2   = DOT_PRODUCT(RSS(:),RSS(:))179:                   RSS2   = DOT_PRODUCT(RSS(:),RSS(:))
210:                   R2     = 1.D0/RSS2  ! This is an inverse (square) length180:                   R2     = 1.D0/RSS2
211:                   ABSR   = DSQRT(RSS2) ! This is a real length181:                   ABSR   = DSQRT(RSS2)
212:                   ENERGY = ENERGY + CH2O*STCHRG(I)*STCHRG(J)/ABSR182:                   ENERGY = ENERGY + CH2O*STCHRG(I)*STCHRG(J)/ABSR
213: 183: 
214:                   ! Calculate the gradient and second derivatives as before. 
215:                   IF (GTEST .OR. STEST) THEN184:                   IF (GTEST .OR. STEST) THEN
216: !     DVDR = DVDR/R = -CH20 * (q_i * q_j / (r_ij)^3)185: !     DVDR = DVDR/R
217:                      DVDR(J7,J8) =-CH2O*STCHRG(I)*STCHRG(J)*R2/ABSR186:                      DVDR(J7,J8) =-CH2O*STCHRG(I)*STCHRG(J)*R2/ABSR
218:                      DVDR(J8,J7) = DVDR(J7,J8)187:                      DVDR(J8,J7) = DVDR(J7,J8)
219: 188: !               write(*,*) "Atoms ", J7, J8   ! sn402
 189: !               write(*,*) "Derivative ", DVDR(J7, J8)
220:                      DOTI1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J7,:))190:                      DOTI1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J7,:))
221:                      DOTI2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J7,:))191:                      DOTI2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J7,:))
222:                      DOTI3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J7,:))192:                      DOTI3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J7,:))
223: 193: 
224:                      DOTJ1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J8,:))194:                      DOTJ1(J7,J8) = DOT_PRODUCT(RSS(:),DR1(J8,:))
225:                      DOTJ2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J8,:))195:                      DOTJ2(J7,J8) = DOT_PRODUCT(RSS(:),DR2(J8,:))
226:                      DOTJ3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J8,:))196:                      DOTJ3(J7,J8) = DOT_PRODUCT(RSS(:),DR3(J8,:))
227: 197: 
228:                      ! Calculate the gradient here. 
229:                      ! CoM terms are quite simple, just add up the site terms (J7,J8) 
230:                      ! that belong to each molecule (J3,J4) 
231:                      G(J3-2:J3)  = G(J3-2:J3) + DVDR(J7,J8)*RSS(:)198:                      G(J3-2:J3)  = G(J3-2:J3) + DVDR(J7,J8)*RSS(:)
232:                      G(J4-2:J4)  = G(J4-2:J4) - DVDR(J7,J8)*RSS(:)199:                      G(J4-2:J4)  = G(J4-2:J4) - DVDR(J7,J8)*RSS(:)
233: 200: !write(*,*) "Gradient for atoms ", J7, J8
234:                      ! AA terms are more complicated to obtain from cartesian terms.201: !write(*,*) DVDR(J7,J8)*RSS(:)
235:                      G(J5-2)     = G(J5-2) + DVDR(J7,J8)*DOTI1(J7,J8)202:                      G(J5-2)     = G(J5-2) + DVDR(J7,J8)*DOTI1(J7,J8)
236:                      G(J5-1)     = G(J5-1) + DVDR(J7,J8)*DOTI2(J7,J8)203:                      G(J5-1)     = G(J5-1) + DVDR(J7,J8)*DOTI2(J7,J8)
237:                      G(J5)       = G(J5)   + DVDR(J7,J8)*DOTI3(J7,J8)204:                      G(J5)       = G(J5)   + DVDR(J7,J8)*DOTI3(J7,J8)
238: 205: 
239:                      G(J6-2)     = G(J6-2) - DVDR(J7,J8)*DOTJ1(J7,J8)206:                      G(J6-2)     = G(J6-2) - DVDR(J7,J8)*DOTJ1(J7,J8)
240:                      G(J6-1)     = G(J6-1) - DVDR(J7,J8)*DOTJ2(J7,J8)207:                      G(J6-1)     = G(J6-1) - DVDR(J7,J8)*DOTJ2(J7,J8)
241:                      G(J6)       = G(J6)   - DVDR(J7,J8)*DOTJ3(J7,J8)208:                      G(J6)       = G(J6)   - DVDR(J7,J8)*DOTJ3(J7,J8)
242: 209: 
243:                      ! D2VDR2 = 1/r d/dr(1/r dv/dr) 
244:                      D2VDR2(J7,J8) = 3.D0*CH2O*STCHRG(I)*STCHRG(J)*R2*R2/ABSR210:                      D2VDR2(J7,J8) = 3.D0*CH2O*STCHRG(I)*STCHRG(J)*R2*R2/ABSR
 211: !                write(*,*) "D2VDR2", J7, J8
 212: !               write(*,*) D2VDR2(J7,J8)
245:                   ENDIF213:                   ENDIF
246: 214: 
247:                   IF (STEST) THEN215:                   IF (STEST) THEN
248: 216: 
249:                      D2VDR2(J8,J7) = D2VDR2(J7,J8)217:                      D2VDR2(J8,J7) = D2VDR2(J7,J8)
250:                      DOTI1(J8,J7)  =-DOTJ1(J7,J8)218:                      DOTI1(J8,J7)  =-DOTJ1(J7,J8)
251:                      DOTI2(J8,J7)  =-DOTJ2(J7,J8)219:                      DOTI2(J8,J7)  =-DOTJ2(J7,J8)
252:                      DOTI3(J8,J7)  =-DOTJ3(J7,J8)220:                      DOTI3(J8,J7)  =-DOTJ3(J7,J8)
253:                      DOTJ1(J8,J7)  =-DOTI1(J7,J8)221:                      DOTJ1(J8,J7)  =-DOTI1(J7,J8)
254:                      DOTJ2(J8,J7)  =-DOTI2(J7,J8)222:                      DOTJ2(J8,J7)  =-DOTI2(J7,J8)
257:                   ENDIF225:                   ENDIF
258: 226: 
259:                ENDDO227:                ENDDO
260: 228: 
261:             ENDDO229:             ENDDO
262: 230: 
263:          ENDDO231:          ENDDO
264: 232: 
265:       ENDDO233:       ENDDO
266: 234: 
267: !      IF (.FALSE.) THEN 
268: !         write(*,*) "Calculating Cartesian Hessian for debugging." 
269: !         CALL Cart_hess(R,DVDR,D2VDR2) 
270: !         STOP 
271: !      ENDIF 
272:  
273:       ! Now use the derivatives saved in the previous section to compute the Hessian. 
274:       ! Note, the form of the expressions here is slightly different to that used for ljpshift etc.,  
275:       ! because the 2nd derivative tensor has been defined differently. 
276:       IF (STEST) THEN235:       IF (STEST) THEN
277: 236: 
278:          DO J1 = 1, REALNATOMS237:          DO J1 = 1, REALNATOMS
279: 238: 
280:             J3 = 3*J1  ! Index of the last CoM coord for this molecule239:             J3 = 3*J1
281:             J5 = OFFSET + J3  ! Index of the last AA coord for this molecule240:             J5 = OFFSET + J3
282: 241: 
283:             DO J2 = 1, REALNATOMS242:             DO J2 = 1, REALNATOMS
284: 243: 
285:                IF (J1 == J2) CYCLE244:                IF (J1 == J2) CYCLE
286: 245: 
287:                J4 = 3*J2246:                J4 = 3*J2
288:                J6 = OFFSET + J4247:                J6 = OFFSET + J4
289: 248: 
290:                DO I = 1, NRBSITES249:                DO I = 1, NRBSITES
291: 250: 
376: !     pi2,pi3335: !     pi2,pi3
377:                      DUMMY           = D2VDR2(J7,J8)*DOTI2(J7,J8)*DOTI3(J7,J8) + DVDR(J7,J8)*DOT_PRODUCT(DR3(J7,:),DR2(J7,:)) &336:                      DUMMY           = D2VDR2(J7,J8)*DOTI2(J7,J8)*DOTI3(J7,J8) + DVDR(J7,J8)*DOT_PRODUCT(DR3(J7,:),DR2(J7,:)) &
378:                                      + DVDR(J7,J8)*DOT_PRODUCT(RSS,D2R23(J7,:))337:                                      + DVDR(J7,J8)*DOT_PRODUCT(RSS,D2R23(J7,:))
379:                      HESS(J5-1,J5)   = HESS(J5-1,J5) + DUMMY338:                      HESS(J5-1,J5)   = HESS(J5-1,J5) + DUMMY
380:                      HESS(J5,J5-1)   = HESS(J5,J5-1) + DUMMY339:                      HESS(J5,J5-1)   = HESS(J5,J5-1) + DUMMY
381: !     pi3,pi1340: !     pi3,pi1
382:                      DUMMY           = D2VDR2(J7,J8)*DOTI3(J7,J8)*DOTI1(J7,J8) + DVDR(J7,J8)*DOT_PRODUCT(DR1(J7,:),DR3(J7,:)) &341:                      DUMMY           = D2VDR2(J7,J8)*DOTI3(J7,J8)*DOTI1(J7,J8) + DVDR(J7,J8)*DOT_PRODUCT(DR1(J7,:),DR3(J7,:)) &
383:                                      + DVDR(J7,J8)*DOT_PRODUCT(RSS,D2R31(J7,:))342:                                      + DVDR(J7,J8)*DOT_PRODUCT(RSS,D2R31(J7,:))
384:                      HESS(J5,J5-2)   = HESS(J5,J5-2) + DUMMY343:                      HESS(J5,J5-2)   = HESS(J5,J5-2) + DUMMY
385:                      HESS(J5-2,J5)   = HESS(J5-2,J5) + DUMMY344:                      HESS(J5-2,J5)   = HESS(J5-2,J5) + DUMMY
 345: 
386:                  346:                  
387: !     [3] DIAGONAL ELEMENTS ON OFF-DIAGONAL BLOCKS: DIFFERENT MOLECULES, SAME COORDINATE347: !     [3] DIAGONAL ELEMENTS ON OFF-DIAGONAL BLOCKS: DIFFERENT MOLECULES, SAME COORDINATE
388: 348: 
389: !     xi,xj349: !     xi,xj
390:                      HESS(J3-2,J4-2) = HESS(J3-2,J4-2) - D2VDR2(J7,J8)*RSS(1)*RSS(1) - DVDR(J7,J8)350:                      HESS(J3-2,J4-2) = HESS(J3-2,J4-2) - D2VDR2(J7,J8)*RSS(1)*RSS(1) - DVDR(J7,J8)
391: !     yi,yj351: !     yi,yj
392:                      HESS(J3-1,J4-1) = HESS(J3-1,J4-1) - D2VDR2(J7,J8)*RSS(2)*RSS(2) - DVDR(J7,J8)352:                      HESS(J3-1,J4-1) = HESS(J3-1,J4-1) - D2VDR2(J7,J8)*RSS(2)*RSS(2) - DVDR(J7,J8)
393: !     zi,zj353: !     zi,zj
394:                      HESS(J3,J4)     = HESS(J3,J4)     - D2VDR2(J7,J8)*RSS(3)*RSS(3) - DVDR(J7,J8)354:                      HESS(J3,J4)     = HESS(J3,J4)     - D2VDR2(J7,J8)*RSS(3)*RSS(3) - DVDR(J7,J8)
395: !     pi1,pj1355: !     pi1,pj1
469: 429: 
470:                   ENDDO430:                   ENDDO
471: 431: 
472:                ENDDO432:                ENDDO
473: 433: 
474:             ENDDO434:             ENDDO
475: 435: 
476:          ENDDO436:          ENDDO
477: 437: 
478:       ENDIF438:       ENDIF
479:  
480: !write(*,*) "2nd derivative matrix" 
481: !DO J1=1,UBOUND(D2VDR2,1) 
482: !   DO J2=1,UBOUND(D2VDR2,2) 
483: !      write(*,*) J1, J2, D2VDR2(J1,J2) 
484: !   ENDDO 
485: !ENDDO 
486: !   STOP 
487:  
488: !write(*,*) "Final gradient"   ! sn402439: !write(*,*) "Final gradient"   ! sn402
489: !write(*,*) G440: !write(*,*) G
490: !write(*,*) "Derivative matrix"441: !write(*,*) "Derivative matrix"
491: !write(*,*) DVDR442: !write(*,*) DVDR
492: !write(*,*) "Hessian"443: !write(*,*) "Hessian"
493: !!write(*,*) HESS444: !!write(*,*) HESS
494: !DO J1 = 1, 3*NATOMS, 3445: !DO J1 = 1, 3*NATOMS, 3
495: !    WRITE(*,*) HESS(J1:J1+2,:)446: !    WRITE(*,*) HESS(J1:J1+2,:)
496: !    write(*,*) " "447: !    write(*,*) " "
497: !ENDDO448: !ENDDO
696:                  WRITE(36,'(A15,F6.3,A1,F6.3,A1,F6.3,A3,F6.3,A1,F6.3,A1,F6.3,A15)') "draw cone {",    TV4(3*J9-2)," &647:                  WRITE(36,'(A15,F6.3,A1,F6.3,A1,F6.3,A3,F6.3,A1,F6.3,A1,F6.3,A15)') "draw cone {",    TV4(3*J9-2)," &
697:   &              ",TV4(3*J9-1)," ",TV4(3*J9),"} {",TV6(3*J9-2)," ",TV6(3*J9-1)," ",TV6(3*J9),"} radius 0.15"648:   &              ",TV4(3*J9-1)," ",TV4(3*J9),"} {",TV6(3*J9-2)," ",TV6(3*J9-1)," ",TV6(3*J9),"} radius 0.15"
698:               ENDDO649:               ENDDO
699:            ENDDO650:            ENDDO
700:            CLOSE (UNIT = 36)651:            CLOSE (UNIT = 36)
701:            652:            
702:         ENDDO653:         ENDDO
703:         654:         
704:       END SUBROUTINE VISUALISEMODESTIP4655:       END SUBROUTINE VISUALISEMODESTIP4
705: 656: 
706:  
707:       subroutine cart_hess(R, G, F) 
708:         ! Calculate (and print) the Hessian in cartesian coordinates using the tools developed for MULTIPOT.  
709:         ! This should only be used for debugging! 
710:  
711:         USE ISOTROPIC_POTENTIALS, ONLY: ISOTROPIC_HESSIAN 
712:         USE COMMONS, ONLY: NATOMS, NRBSITES 
713:  
714:         IMPLICIT NONE 
715:  
716:         DOUBLE PRECISION, INTENT(IN) :: R(NRBSITES*NATOMS/2,3) ! Cartesian coordinates of the system 
717:         ! Cartesian G and F tensors (defined in the body of the code, given above) 
718:         DOUBLE PRECISION, INTENT(IN) :: G(NRBSITES*NATOMS/2,NRBSITES*NATOMS/2), F(NRBSITES*NATOMS/2,NRBSITES*NATOMS/2) 
719:         DOUBLE PRECISION :: X(3*NRBSITES*NATOMS/2) 
720:         DOUBLE PRECISION :: R2DUMMY(NRBSITES*NATOMS/2,NRBSITES*NATOMS/2) 
721:         DOUBLE PRECISION :: TMP_HESS(3*NRBSITES*NATOMS/2,3*NRBSITES*NATOMS/2) 
722:         INTEGER :: J1, J2 
723:         write(*,*) "In cart_hess" 
724:  
725:         DO J1=1,NRBSITES*NATOMS/2 
726:            X(3*J1-2:3*J1)=R(J1,:) 
727:         ENDDO 
728:         R2DUMMY(:,:)=1.0D0 
729:  
730:         CALL ISOTROPIC_HESSIAN(NRBSITES*NATOMS/2,X,G,F,R2DUMMY,TMP_HESS) 
731:  
732:         DO J1=1,3*NRBSITES*NATOMS/2 
733:            DO J2=1,3*NRBSITES*NATOMS/2 
734:               write(*,*) J1, J2, TMP_HESS(J1,J2) 
735:            ENDDO 
736:         ENDDO 
737:  
738:      end subroutine cart_hess 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0