1 SUBROUTINE DCART (COORD,DXYZ) 2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 3 INCLUDE 'SIZES' 4 DIMENSION COORD(3,*), DXYZ(3,*) 5 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 6 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 7 2 NCLOSE,NOPEN,NDUMY,FRACT 8 COMMON /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK) 9C*********************************************************************** 10C 11C DCART CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE 12C CARTESIAN COORDINATES. THIS IS DONE BY FINITE DIFFERENCES. 13C 14C THE MAIN ARRAYS IN DCART ARE: 15C DXYZ ON EXIT CONTAINS THE CARTESIAN DERIVATIVES. 16C 17C*********************************************************************** 18 COMMON /KEYWRD/ KEYWRD 19 COMMON /EULER / TVEC(3,3), ID 20 COMMON /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE 21 COMMON /UCELL / L1L,L2L,L3L,L1U,L2U,L3U 22 COMMON /DCARTC/ K1L,K2L,K3L,K1U,K2U,K3U 23 COMMON /NUMCAL/ NUMCAL 24C COSMO change 25 LOGICAL ISEPS, USEPS , UPDA 26 COMMON /ISEPS/ ISEPS, USEPS, UPDA 27C end of COSMO change 28 CHARACTER*241 KEYWRD 29 DIMENSION PDI(171),PADI(171),PBDI(171), 30 1CDI(3,2),NDI(2),LSTOR1(6), LSTOR2(6), ENG(3) 31 LOGICAL DEBUG, FORCE, MAKEP, ANADER, LARGE 32 EQUIVALENCE (LSTOR1(1),L1L), (LSTOR2(1), K1L) 33 SAVE CHNGE, CHNGE2, ANADER, DEBUG, FORCE 34 DATA ICALCN/0/ 35 DATA CHNGE /1.D-4/ 36 CHNGE2=CHNGE*0.5D0 37* 38* CHNGE IS A MACHINE-PRECISION DEPENDENT CONSTANT 39* CHNGE2=CHNGE/2 40* 41 IF (ICALCN.NE.NUMCAL) THEN 42 ICALCN=NUMCAL 43 LARGE = (INDEX(KEYWRD,'LARGE') .NE. 0) 44 ANADER= (INDEX(KEYWRD,'ANALYT') .NE. 0) 45 DEBUG = (INDEX(KEYWRD,'DCART') .NE. 0) 46 FORCE = (INDEX(KEYWRD,'PREC')+INDEX(KEYWRD,'FORCE') .NE. 0) 47 ENDIF 48 NCELLS=(L1U-L1L+1)*(L2U-L2L+1)*(L3U-L3L+1) 49 DO 10 I=1,6 50 LSTOR2(I)=LSTOR1(I) 51 10 LSTOR1(I)=0 52 IOFSET=(NCELLS+1)/2 53 NUMTOT=NUMAT*NCELLS 54 DO 20 I=1,NUMTOT 55 DO 20 J=1,3 56 20 DXYZ(J,I)=0.D0 57 IF(ANADER) REWIND 2 58 DO 130 II=1,NUMAT 59 III=NCELLS*(II-1)+IOFSET 60 IM1=II 61 IF=NFIRST(II) 62 IM=NMIDLE(II) 63 IL=NLAST(II) 64 NDI(2)=NAT(II) 65 DO 30 I=1,3 66 30 CDI(I,2)=COORD(I,II) 67 DO 130 JJ=1,IM1 68 JJJ=NCELLS*(JJ-1) 69C FORM DIATOMIC MATRICES 70 JF=NFIRST(JJ) 71 JM=NMIDLE(JJ) 72 JL=NLAST(JJ) 73C GET FIRST ATOM 74 NDI(1)=NAT(JJ) 75 MAKEP=.TRUE. 76 DO 120 IK=K1L,K1U 77 DO 120 JK=K2L,K2U 78 DO 120 KL=K3L,K3U 79 JJJ=JJJ+1 80* KKK=KKK-1 81 DO 40 L=1,3 82 40 CDI(L,1)=COORD(L,JJ)+TVEC(L,1)*IK+TVEC(L,2)*JK+TVEC 83 1(L,3)*KL 84 IF(.NOT. MAKEP) GOTO 90 85 MAKEP=.FALSE. 86 IJ=0 87 DO 50 I=JF,JL 88 K=I*(I-1)/2+JF-1 89 DO 50 J=JF,I 90 IJ=IJ+1 91 K=K+1 92 PADI(IJ)=PA(K) 93 PBDI(IJ)=PB(K) 94 50 PDI(IJ)=P(K) 95C GET SECOND ATOM FIRST ATOM INTERSECTION 96 DO 80 I=IF,IL 97 L=I*(I-1)/2 98 K=L+JF-1 99 DO 60 J=JF,JL 100 IJ=IJ+1 101 K=K+1 102 PADI(IJ)=PA(K) 103 PBDI(IJ)=PB(K) 104 60 PDI(IJ)=P(K) 105 K=L+IF-1 106 DO 70 L=IF,I 107 K=K+1 108 IJ=IJ+1 109 PADI(IJ)=PA(K) 110 PBDI(IJ)=PB(K) 111 70 PDI(IJ)=P(K) 112 80 CONTINUE 113 90 CONTINUE 114 IF(II.EQ.JJ) GOTO 120 115 IF(ANADER)THEN 116 CALL ANALYT(PDI,PADI,PBDI,CDI,NDI,JF,JL,IF,IL 117 1, ENG) 118 DO 100 K=1,3 119 DXYZ(K,III)=DXYZ(K,III)-ENG(K) 120 100 DXYZ(K,JJJ)=DXYZ(K,JJJ)+ENG(K) 121 ELSE 122 IF( .NOT. FORCE) THEN 123 CDI(1,1)=CDI(1,1)+CHNGE2 124 CDI(2,1)=CDI(2,1)+CHNGE2 125 CDI(3,1)=CDI(3,1)+CHNGE2 126 CALL DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF,IM 127 1,IL, AA,1) 128 ENDIF 129 DO 110 K=1,3 130 IF( FORCE )THEN 131 CDI(K,2)=CDI(K,2)-CHNGE2 132 CALL DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF 133 1,IM,IL, AA,1) 134 ENDIF 135 CDI(K,2)=CDI(K,2)+CHNGE 136 CALL DHC(PDI,PADI,PBDI,CDI,NDI,JF,JM,JL,IF,IM 137 1,IL, EE,2) 138 CDI(K,2)=CDI(K,2)-CHNGE2 139 IF( .NOT. FORCE) CDI(K,2)=CDI(K,2)-CHNGE2 140 DERIV=(AA-EE)*23.061D0/CHNGE 141 DXYZ(K,III)=DXYZ(K,III)-DERIV 142 DXYZ(K,JJJ)=DXYZ(K,JJJ)+DERIV 143 110 CONTINUE 144 ENDIF 145 120 CONTINUE 146 130 CONTINUE 147 IF(NNHCO.NE.0)THEN 148C 149C NOW ADD IN MOLECULAR-MECHANICS CORRECTION TO THE H-N-C=O TORSION 150C 151 DEL=1.D-8 152 DO 160 I=1,NNHCO 153 DO 150 J=1,4 154 DO 140 K=1,3 155 COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))-DEL 156 CALL DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),NHCO(4, 157 1I),ANGLE) 158 REFH=HTYPE(ITYPE)*SIN(ANGLE)**2 159 COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))+DEL*2.D0 160 CALL DIHED(COORD,NHCO(1,I),NHCO(2,I),NHCO(3,I),NHCO(4, 161 1I),ANGLE) 162 COORD(K,NHCO(J,I))=COORD(K,NHCO(J,I))-DEL 163 HEAT=HTYPE(ITYPE)*SIN(ANGLE)**2 164 SUM=(REFH-HEAT)/(2.D0*DEL) 165 DXYZ(K,NHCO(J,I))=DXYZ(K,NHCO(J,I))-SUM 166 140 CONTINUE 167 150 CONTINUE 168 160 CONTINUE 169 ENDIF 170C COSMO change A. Klamt 171C analytic calculation of the gradient of the dielectric energy A.Klamt 172 IF (USEPS) CALL DIEGRD(COORD,DXYZ) 173C DO 170 I=1,6 174C 170 LSTOR1(I)=LSTOR2(I) 175 IF ( .NOT. DEBUG) RETURN 176 IW = 6 177 WRITE(IW,'(//10X,''CARTESIAN COORDINATE DERIVATIVES'',//3X, 178 1''NUMBER ATOM '',5X,''X'',12X,''Y'',12X,''Z'',/)') 179 IF(NCELLS.EQ.1)THEN 180 WRITE(IW,'(2I6,F13.6,2F13.6)') 181 1 (I,NAT(I),(DXYZ(J,I),J=1,3),I=1,NUMTOT) 182 ELSEIF(LARGE)THEN 183 WRITE(IW,'(2I6,F13.6,2F13.6)') 184 1 (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I),J=1,3),I=1,NUMTOT) 185 ELSE 186 WRITE(IW,'(2I6,F13.6,2F13.6)') 187 1 (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I)+DXYZ(J,I+1)+DXYZ(J,I+2) 188 2,J=1,3),I=1,NUMTOT,3) 189 ENDIF 190 IROT = 2 191 IF (ANADER) REWIND IROT 192C end of COSMO (A. Klamt) changes 193 IF ( .NOT. DEBUG) RETURN 194 WRITE(6,'(//10X,''CARTESIAN COORDINATE DERIVATIVES'',//3X, 195 1''NUMBER ATOM '',5X,''X'',12X,''Y'',12X,''Z'',/)') 196 IF(NCELLS.EQ.1)THEN 197 WRITE(6,'(2I6,F13.6,2F13.6)') 198 1 (I,NAT(I),(DXYZ(J,I),J=1,3),I=1,NUMTOT) 199 ELSEIF(LARGE)THEN 200 WRITE(6,'(2I6,F13.6,2F13.6)') 201 1 (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I),J=1,3),I=1,NUMTOT) 202 ELSE 203 WRITE(6,'(2I6,F13.6,2F13.6)') 204 1 (I,NAT((I-1)/NCELLS+1),(DXYZ(J,I)+DXYZ(J,I+1)+DXYZ(J,I+2) 205 2,J=1,3),I=1,NUMTOT,3) 206 ENDIF 207 IF (ANADER) REWIND 2 208 RETURN 209 END 210 SUBROUTINE DHC (P,PA,PB,XI,NAT,IF,IM,IL,JF,JM,JL,DENER,MODE) 211 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 212 DIMENSION P(*), PA(*), PB(*) 213 DIMENSION XI(3,*),NFIRST(2),NMIDLE(2),NLAST(2),NAT(*) 214C*********************************************************************** 215C 216C DHC CALCULATES THE ENERGY CONTRIBUTIONS FROM THOSE PAIRS OF ATOMS 217C THAT HAVE BEEN MOVED BY SUBROUTINE DERIV. 218C 219C*********************************************************************** 220 COMMON /KEYWRD/ KEYWRD 221 1 /ONELEC/ USS(107),UPP(107),UDD(107) 222 COMMON /EULER / TVEC(3,3), ID 223 COMMON /NUMCAL/ NUMCAL 224 SAVE ICALCN, WLIM, UHF 225 CHARACTER*241 KEYWRD 226 LOGICAL UHF, CUTOFF 227 DIMENSION H(171), SHMAT(9,9), F(171), 228 1 WJ(100), E1B(10), E2A(10), WK(100), W(100), 229 2 WJS(100), WKS(100) 230 DOUBLE PRECISION WJS, WKS 231 DATA ICALCN /0/ 232 IF( ICALCN.NE.NUMCAL) THEN 233 ICALCN=NUMCAL 234 WLIM=4.D0 235 IF(ID.EQ.0)WLIM=0.D0 236 UHF=(INDEX(KEYWRD,'UHF') .NE. 0) 237 ENDIF 238 NFIRST(1)=1 239 NMIDLE(1)=IM-IF+1 240 NLAST(1)=IL-IF+1 241 NFIRST(2)=NLAST(1)+1 242 NMIDLE(2)=NFIRST(2)+JM-JF 243 NLAST(2)=NFIRST(2)+JL-JF 244 LINEAR=(NLAST(2)*(NLAST(2)+1))/2 245 DO 10 I=1,LINEAR 246 F(I)=0.D0 247 10 H(I)=0.0D00 248 DO 20 I=1,LINEAR 249 20 F(I)=H(I) 250 JA=NFIRST(2) 251 JB=NLAST(2) 252 JC=NMIDLE(2) 253 IA=NFIRST(1) 254 IB=NLAST(1) 255 IC=NMIDLE(1) 256 J=2 257 I=1 258 NJ=NAT(2) 259 NI=NAT(1) 260 CALL H1ELEC(NI,NJ,XI(1,1),XI(1,2),SHMAT) 261 IF(NAT(1).EQ.102.OR.NAT(2).EQ.102) THEN 262 K=(JB*(JB+1))/2 263 DO 30 J=1,K 264 30 H(J)=0.D0 265 ELSE 266 J1=0 267 DO 40 J=JA,JB 268 JJ=J*(J-1)/2 269 J1=J1+1 270 I1=0 271 DO 40 I=IA,IB 272 JJ=JJ+1 273 I1=I1+1 274 H(JJ)=SHMAT(I1,J1) 275 F(JJ)=SHMAT(I1,J1) 276 40 CONTINUE 277 ENDIF 278 KR=1 279 IF(ID.EQ.0)THEN 280 CALL ROTATE (NJ,NI,XI(1,2),XI(1,1),W(KR),KR,E2A,E1B,ENUCLR,100. 281 1D0) 282 ELSE 283 CALL SOLROT (NJ,NI,XI(1,2),XI(1,1),WJ,WK,KR,E2A,E1B,ENUCLR,100. 284 1D0) 285 IF(MODE.EQ.1)CUTOFF=(WJ(1).LT.WLIM) 286 IF(CUTOFF)THEN 287 DO 50 I=1,KR-1 288 50 WK(I)=0.D0 289 ENDIF 290 DO 60 I=1,KR-1 291 WJS(I)=WJ(I) 292 WKS(I)=WK(I) 293 60 CONTINUE 294 ENDIF 295C 296C * ENUCLR IS SUMMED OVER CORE-CORE REPULSION INTEGRALS. 297C 298 I2=0 299 DO 70 I1=IA,IC 300 II=I1*(I1-1)/2+IA-1 301 DO 70 J1=IA,I1 302 II=II+1 303 I2=I2+1 304 H(II)=H(II)+E1B(I2) 305 70 F(II)=F(II)+E1B(I2) 306 DO 80 I1=IC+1,IB 307 II=(I1*(I1+1))/2 308 F(II)=F(II)+E1B(1) 309 80 H(II)=H(II)+E1B(1) 310 I2=0 311 DO 90 I1=JA,JC 312 II=I1*(I1-1)/2+JA-1 313 DO 90 J1=JA,I1 314 II=II+1 315 I2=I2+1 316 H(II)=H(II)+E2A(I2) 317 90 F(II)=F(II)+E2A(I2) 318 DO 100 I1=JC+1,JB 319 II=(I1*(I1+1))/2 320 F(II)=F(II)+E2A(1) 321 100 H(II)=H(II)+E2A(1) 322 CALL FOCK2(F,P,PA,W, WJS, WKS,2,NAT,NFIRST,NMIDLE,NLAST) 323 EE=HELECT(NLAST(2),PA,H,F) 324 IF( UHF ) THEN 325 DO 110 I=1,LINEAR 326 110 F(I)=H(I) 327 CALL FOCK2(F,P,PB,W, WJS, WKS,2,NAT,NFIRST,NMIDLE,NLAST) 328 EE=EE+HELECT(NLAST(2),PB,H,F) 329 ELSE 330 EE=EE*2.D0 331 ENDIF 332 DENER=EE+ENUCLR 333 RETURN 334C 335 END 336