1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck ecktrn */ 20 SUBROUTINE ECKTRN(GRDCAR,HESCAR,DIP0,SUSTOT,GTRAN,QUADT,TMAT, 21 & GTRANT,POLARS,POLDD,DIANQC,SPNTOT, 22 & MAGSUS,MOLGFA,QUADRU,SHIELD,SPINRO,POLAR,ALFA,NQCC, 23 & SPNSPN,NCART,NCRTOT,MXFR,NFRVAL,FRVAL,GEOM, 24 & AMASS,NATTYP,ECK,GEOM2,IPRINT,HESSIA,WORK,LWORK) 25C 26C Define an Eckart frame, and transform properties to this frame. 27C Based on the program ECKART by Juha Vaara: 28C Juha Vaara 130398 29C last change 310398 30C 31C Modified for use in Dalton by Kenneth Ruud, San Diego April 1999 32C 33#include "implicit.h" 34#include "priunit.h" 35#include "mxcent.h" 36#include "maxaqn.h" 37#include "maxorb.h" 38#include "nuclei.h" 39C 40 PARAMETER (NUCMAX = 15,NSTMAX = 10,NTRIAL = 1500, 41 & TOLX = 1.0D-12,TOLF = 1.0D-12) 42 PARAMETER (D0 = 0.0D0, D1 = 1.0D0) 43 LOGICAL MAGSUS, MOLGFA, QUADRU, SHIELD, SPINRO, POLAR, ALFA, 44 & NQCC, SPNSPN, PLANAR, LINEAR, HESSIA 45C 46 DIMENSION GEOM(3,MXCENT), AMASS(MXCENT), NATTYP(MXCENT), 47 & ECK(3,MXCENT), GEOM2(NUCDEP,3), EQCM(3), XN(3) 48 DIMENSION CMAT(3,3), EIGVAL(3), ANGMOM(3), TINERT(3,3), 49 & OMEGA(3,3), BVEC(3), XJMAT(3,3), ECKP(3) 50 DIMENSION GRDCAR(NCRTOT), HESCAR(NCART,NCART), DIP0(3), 51 & SUSTOT(3,3), GTRAN(3,3), QUADT(3,3), TMAT(3,3,MXCENT), 52 & GTRANT(3,3,MXCENT), POLARS(3,3), POLDD(3,3,MXFR), 53 & DIANQC(3,3,MXCENT), SPNTOT(MXCOOR,MXCOOR), 54 & FRVAL(NFRVAL), WORK(LWORK) 55C 56#include "orgcom.h" 57#include "cbiwlk.h" 58C 59C Define instantaneous geometry 60C 61 CALL CMMASS(GEOM,AMASS,NATTYP,WORK,IPRINT) 62 IF (IPRINT .GE. 4) THEN 63 CALL HEADER('Instantaneous geometry:',0) 64 DO I = 1, NUCDEP 65 WRITE(LUPRI,'(I2,3F15.10)') I,GEOM(1,I),GEOM(2,I), 66 & GEOM(3,I) 67 END DO 68 END IF 69C 70C Center of mass for equilibrium geometry 71C 72 IF (IPRINT .GE. 4) THEN 73 CALL HEADER('Equilibrium geometry',0) 74 DO I = 1, NUCDEP 75 WRITE(LUPRI,'(I2,3F15.8)') I,ECKGEO(1,I),ECKGEO(2,I), 76 & ECKGEO(3,I) 77 END DO 78 END IF 79 TOTMAS = D0 80 CALL DZERO(EQCM,3) 81 DO I = 1, NUCDEP 82 RMASS = AMASS(I) 83 TOTMAS = TOTMAS + RMASS 84 DO J = 1, 3 85 EQCM(J) = EQCM(J) + RMASS*ECKGEO(J,I) 86 END DO 87 END DO 88 DO J = 1, 3 89 EQCM(J) = EQCM(J)/TOTMAS 90 END DO 91C 92C We now start to build the Eckart frame 93C 94 IF (IPRINT .GE. 2) THEN 95 WRITE(LUPRI,'(A,3F15.6)') ' equilibrium CM : ',EQCM(1), 96 & EQCM(2),EQCM(3) 97 WRITE(LUPRI,'(A,3F15.6)') ' instantaneous CM: ',CMXYZ(1), 98 & CMXYZ(2),CMXYZ(3) 99 END IF 100C 101C construct B vector and J matrix for this isotopomer 102C 103 CALL DZERO(BVEC,3) 104 CALL DZERO(XJMAT,9) 105 DO 16 J = 1, NUCDEP 106 GEOM(1,J) = GEOM(1,J) - CMXYZ(1) 107 GEOM(2,J) = GEOM(2,J) - CMXYZ(2) 108 GEOM(3,J) = GEOM(3,J) - CMXYZ(3) 109 ECK(1,J) = ECKGEO(1,J) - EQCM(1) 110 ECK(2,J) = ECKGEO(2,J) - EQCM(2) 111 ECK(3,J) = ECKGEO(3,J) - EQCM(3) 112 CALL CROSSP(ECK(1,J),GEOM(1,J),ANGMOM) 113 BVEC(1) = BVEC(1) + AMASS(J)*ANGMOM(1) 114 BVEC(2) = BVEC(2) + AMASS(J)*ANGMOM(2) 115 BVEC(3) = BVEC(3) + AMASS(J)*ANGMOM(3) 116 TMP = DDOT(3,ECK(1,J),1,GEOM(1,J),1) 117 XJMAT(1,1) = XJMAT(1,1) + AMASS(J)*(TMP -ECK(1,J)*GEOM(1,J)) 118 XJMAT(2,2) = XJMAT(2,2) + AMASS(J)*(TMP -ECK(2,J)*GEOM(2,J)) 119 XJMAT(3,3) = XJMAT(3,3) + AMASS(J)*(TMP -ECK(3,J)*GEOM(3,J)) 120 XJMAT(1,2) = XJMAT(1,2) - AMASS(J)*ECK(1,J)*GEOM(2,J) 121 XJMAT(1,3) = XJMAT(1,3) - AMASS(J)*ECK(1,J)*GEOM(3,J) 122 XJMAT(2,1) = XJMAT(2,1) - AMASS(J)*ECK(2,J)*GEOM(1,J) 123 XJMAT(2,3) = XJMAT(2,3) - AMASS(J)*ECK(2,J)*GEOM(3,J) 124 XJMAT(3,1) = XJMAT(3,1) - AMASS(J)*ECK(3,J)*GEOM(1,J) 125 XJMAT(3,2) = XJMAT(3,2) - AMASS(J)*ECK(3,J)*GEOM(2,J) 126 16 CONTINUE 127C 128C Mismatch in the ordering of atom and coordinate for the geometry in 129C CMMASS and WLKDIN 130C 131 DO I = 1, NUCDEP 132 DO J = 1, 3 133 GEOM2(I,J) = GEOM(J,I) 134 END DO 135 END DO 136C 137C Moment of inertia for instantaneous geometry 138C 139 CALL WLKDIN(GEOM2,AMASS,NUCDEP,ANGMOM,TINERT,OMEGA, 140 & EIGVAL,CMAT,.TRUE.,PLANAR,LINEAR) 141C 142C We will need rotational g tensors and spin-rotation constants in 143C space-fixed xyz axis system 144C 145 IF (MOLGFA .AND. .NOT. LINEAR) THEN 146 CALL DGEMM('N','N',3,3,3,1.D0, 147 & CMAT,3, 148 & GTRAN,3,0.D0, 149 & OMEGA,3) 150 CALL DGEMM('N','T',3,3,3,1.D0, 151 & OMEGA,3, 152 & CMAT,3,0.D0, 153 & GTRAN,3) 154 IF (IPRINT .GE. 4) THEN 155 CALL HEADER('Rotational g tensor in Cartesian '// 156 & 'input frame',0) 157 CALL OUTPUT(GTRAN,1,3,1,3,3,3,1,LUPRI) 158 END IF 159 END IF 160C 161 IF (SPINRO) THEN 162 DO I = 1, NUCDEP 163 CALL DGEMM('N','N',3,3,3,1.D0, 164 & CMAT,3, 165 & GTRANT(1,1,I),3,0.D0, 166 & OMEGA,3) 167 CALL DGEMM('N','T',3,3,3,1.D0, 168 & OMEGA,3, 169 & CMAT,3,0.D0, 170 & GTRANT(1,1,I),3) 171 END DO 172 IF (IPRINT .GE. 4) THEN 173 DO I = 1, NUCDEP 174 CALL HEADER('Spin-rotation tensor in Cartesian '// 175 & 'input frame',0) 176 CALL OUTPUT(GTRANT(1,1,I),1,3,1,3,3,3,1,LUPRI) 177 END DO 178 END IF 179 END IF 180C 181 IF (IPRINT .GE. 3) THEN 182 CALL HEADER('Equilibrium geometry in the CM frame:',0) 183 DO 17 J = 1, NUCDEP 184 WRITE(LUPRI,'(I2,3F15.10)') J,ECK(1,J),ECK(2,J), 185 & ECK(3,J) 186 17 CONTINUE 187 CALL HEADER('Instantaneous geometry in the CM frame:',0) 188 DO 18 J = 1, NUCDEP 189 WRITE(LUPRI,'(I2,3F15.10)') J,GEOM(1,J),GEOM(2,J),GEOM(3,J) 190 18 CONTINUE 191 END IF 192 IF (IPRINT .GE. 2) THEN 193 CALL HEADER('B-vector: ',0) 194 WRITE(LUPRI,'(A,F14.6,2F15.6)') 'B',(BVEC(I), I=1,3) 195 CALL HEADER('J-matrix:',0) 196 CALL OUTPUT(XJMAT,1,3,1,3,3,3,1,LUPRI) 197 END IF 198C 199 ECKP(1) = 1.0D0 200 ECKP(2) = 1.0D0 201 ECKP(3) = 1.0D0 202C 203 CALL ECK_MNEWT(NTRIAL,ECKP,BVEC,XJMAT,TOLX,TOLF) 204 XN(1) = SIN(ECKP(1))*COS(ECKP(2)) 205 XN(2) = SIN(ECKP(1))*SIN(ECKP(2)) 206 XN(3) = COS(ECKP(1)) 207 AN = ECKP(3) 208C 209 IF (IPRINT .GE. 2) THEN 210 CALL HEADER('Parameters of the Eckart frame:',0) 211 WRITE(LUPRI,'(A,F13.6,2F15.6)') 'NV',(XN(I), I = 1, 3) 212 WRITE(LUPRI,'(A,F13.6)') 'AN',AN 213 END IF 214C 215C 216 IF (ABS(XN(1)).LT.1.0D-7 .AND. ABS(XN(2)).LT.1.0D-7 .AND. 217 & ABS(D1-XN(3)).LT.1.0D-7 ) THEN 218 WRITE(LUPRI,'(A)') ' (rotation in the xy-plane only;' 219 & //' zeroing a redundant parameter)' 220 ECKP(2) = D0 221 ECKP(1) = D0 222 ENDIF 223C 224 CALL CLCGLD(ECKP,CMAT) 225 IF (IPRINT .GE. 2) THEN 226 CALL HEADER('C-matrix:',0) 227 CALL OUTPUT(CMAT,1,3,1,3,3,3,1,LUPRI) 228 END IF 229C 230 IF (IPRINT .GE. 4) 231 & CALL HEADER('Instantaneous geometry in the Eckart frame:',0) 232 DO 20 J = 1, NUCDEP 233 CALL DGEMM('N','N',3,1,3,1.D0, 234 & CMAT,3, 235 & GEOM(1,J),3,0.D0, 236 & ANGMOM,3) 237 CALL DCOPY(3,ANGMOM,1,GEOM(1,J),1) 238 IF (IPRINT .GE. 4) WRITE(LUPRI,'(A,I2,3F15.10)') 239 & 'IN',J, GEOM(1,J), GEOM(2,J), GEOM(3,J) 240 20 CONTINUE 241C 242C Transform all properties to Eckart frame 243C 244 CALL DGEMM('N','N',3,1,3,1.D0, 245 & CMAT,3, 246 & DIP0,3,0.D0, 247 & ANGMOM,3) 248 CALL DCOPY(3,ANGMOM,1,DIP0,1) 249 IF (IPRINT .GE. 2) THEN 250 CALL HEADER('Dipole moment in Eckart frame',0) 251 CALL OUTPUT(DIP0,1,3,1,1,3,1,1,LUPRI) 252 END IF 253C 254 IF (MAGSUS) THEN 255 CALL DGEMM('N','N',3,3,3,1.D0, 256 & CMAT,3, 257 & SUSTOT,3,0.D0, 258 & OMEGA,3) 259 CALL DGEMM('N','T',3,3,3,1.D0, 260 & OMEGA,3, 261 & CMAT,3,0.D0, 262 & SUSTOT,3) 263 IF (IPRINT .GE. 2) THEN 264 CALL HEADER('Magnetizability tensor ' 265 & //'in the Eckart frame:',0) 266 CALL OUTPUT(SUSTOT,1,3,1,3,3,3,1,LUPRI) 267 ENDIF 268 END IF 269C 270 IF (MOLGFA) THEN 271 CALL DGEMM('N','N',3,3,3,1.D0, 272 & CMAT,3, 273 & GTRAN,3,0.D0, 274 & OMEGA,3) 275 CALL DGEMM('N','T',3,3,3,1.D0, 276 & OMEGA,3, 277 & CMAT,3,0.D0, 278 & GTRAN,3) 279 IF (IPRINT .GE. 2) THEN 280 CALL HEADER('Molecular rotational g-tensor ' 281 & //'in the Eckart frame:',0) 282 CALL OUTPUT(GTRAN,1,3,1,3,3,3,1,LUPRI) 283 ENDIF 284 END IF 285C 286 IF (QUADRU) THEN 287 CALL DGEMM('N','N',3,3,3,1.D0, 288 & CMAT,3, 289 & QUADT,3,0.D0, 290 & OMEGA,3) 291 CALL DGEMM('N','T',3,3,3,1.D0, 292 & OMEGA,3, 293 & CMAT,3,0.D0, 294 & QUADT,3) 295 IF (IPRINT .GE. 2) THEN 296 CALL HEADER('Quadrupole moment tensor ' 297 & //'in the Eckart frame:',0) 298 CALL OUTPUT(QUADT,1,3,1,3,3,3,1,LUPRI) 299 ENDIF 300 END IF 301C 302C 303 IF (HESSIA) THEN 304C 305C Gradient 306C 307 DO I = 1, 3*NUCDEP, 3 308 CALL DGEMM('N','N',3,3,3,1.0D0, 309 & CMAT,3, 310 & GRDCAR(I),3,0.0D0, 311 & OMEGA,3) 312 CALL DGEMM('N','T',3,3,3,1.0D0, 313 & OMEGA,3, 314 & CMAT,3,0.0D0, 315 & GRDCAR(I),3) 316 END DO 317C 318C Hessian 319C 320 DO I = 1, NCART, 3 321 DO J = 1, NCART, 3 322 CALL DGEMM('N','N',3,3,3,1.0D0, 323 & CMAT,3, 324 & HESCAR(I,J),3,0.0D0, 325 & OMEGA,3) 326 CALL DGEMM('N','T',3,3,3,1.0D0, 327 & OMEGA,3, 328 & CMAT,3,0.0D0, 329 & HESCAR(I,J),3) 330 END DO 331 END DO 332 IF (IPRINT .GE. 2) THEN 333 CALL HEADER('Molecular Hessian in the Eckart frame:',0) 334 CALL OUTPAK(HESCAR,NCART,1,LUPRI) 335 END IF 336 END IF 337C 338 IF (SHIELD) THEN 339 DO I =1 , NUCDEP 340 CALL DGEMM('N','N',3,3,3,1.D0, 341 & CMAT,3, 342 & TMAT(1,1,I),3,0.D0, 343 & OMEGA,3) 344 CALL DGEMM('N','T',3,3,3,1.D0, 345 & OMEGA,3, 346 & CMAT,3,0.D0, 347 & TMAT(1,1,I),3) 348 END DO 349 IF (IPRINT .GE. 2) THEN 350 CALL HEADER('Shielding tensors ' 351 & //'in the Eckart frame:',0) 352 DO I = 1, NUCDEP 353 CALL OUTPUT(TMAT(1,1,I),1,3,1,3,3,3,1,LUPRI) 354 END DO 355 ENDIF 356 END IF 357C 358 IF (SPINRO) THEN 359 DO I =1 , NUCDEP 360 CALL DGEMM('N','N',3,3,3,1.D0, 361 & CMAT,3, 362 & GTRANT(1,1,I),3,0.D0, 363 & OMEGA,3) 364 CALL DGEMM('N','T',3,3,3,1.D0, 365 & OMEGA,3, 366 & CMAT,3,0.D0, 367 & GTRANT(1,1,I),3) 368 END DO 369 IF (IPRINT .GE. 2) THEN 370 CALL HEADER('Nuclear spin-rotation tensors ' 371 & //'in the Eckart frame:',0) 372 DO I = 1, NUCDEP 373 CALL OUTPUT(GTRANT(1,1,I),1,3,1,3,3,3,1,LUPRI) 374 END DO 375 ENDIF 376 END IF 377C 378 IF (POLAR) THEN 379 CALL DGEMM('N','N',3,3,3,1.D0, 380 & CMAT,3, 381 & POLARS,3,0.D0, 382 & OMEGA,3) 383 CALL DGEMM('N','T',3,3,3,1.D0, 384 & OMEGA,3, 385 & CMAT,3,0.D0, 386 & POLARS,3) 387 IF (IPRINT .GE. 2) THEN 388 CALL HEADER('Polarzability tensor ' 389 & //'in the Eckart frame:',0) 390 CALL OUTPUT(POLARS,1,3,1,3,3,3,1,LUPRI) 391 ENDIF 392 END IF 393C 394 IF (ALFA) THEN 395 DO IFR = 1, NFRVAL 396 CALL DGEMM('N','N',3,3,3,1.D0, 397 & CMAT,3, 398 & POLDD(1,1,IFR),3,0.D0, 399 & OMEGA,3) 400 CALL DGEMM('N','T',3,3,3,1.D0, 401 & OMEGA,3, 402 & CMAT,3,0.D0, 403 & POLDD(1,1,IFR),3) 404 END DO 405 IF (IPRINT .GE. 2) THEN 406 CALL HEADER('Polarizability tensors ' 407 & //'in the Eckart frame:',0) 408 DO IFR = 1, NFRVAL 409 CALL OUTPUT(POLDD(1,1,IFR),1,3,1,3,3,3,1,LUPRI) 410 END DO 411 ENDIF 412 END IF 413C 414 IF (NQCC) THEN 415 DO I = 1, NUCDEP 416 CALL DGEMM('N','N',3,3,3,1.D0, 417 & CMAT,3, 418 & DIANQC(1,1,I),3,0.D0, 419 & OMEGA,3) 420 CALL DGEMM('N','T',3,3,3,1.D0, 421 & OMEGA,3, 422 & CMAT,3,0.D0, 423 & DIANQC(1,1,I),3) 424 END DO 425 IF (IPRINT .GE. 2) THEN 426 CALL HEADER('Nuclear quadrupole moment tensors ' 427 & //'in the Eckart frame:',0) 428 DO I = 1, NUCDEP 429 CALL OUTPUT(DIANQC(1,1,I),1,3,1,3,3,3,1,LUPRI) 430 END DO 431 ENDIF 432 END IF 433C 434 RETURN 435 END 436C 437 SUBROUTINE CROSSP(A,B,C) 438C 439C calculates (AX,AY,AZ) x (BX,BY,BZ) = CZ,CY,CZ) 440C 441#include "implicit.h" 442 DIMENSION A(3), B(3), C(3) 443C 444 C(1) = A(2)*B(3) - A(3)*B(2) 445 C(2) = A(3)*B(1) - A(1)*B(3) 446 C(3) = A(1)*B(2) - A(2)*B(1) 447 RETURN 448 END 449 SUBROUTINE ECK_USRFUN(ECKP,FVEC,FJAC,BX,BY,BZ, 450 & JXX,JXY,JXZ,JYX,JYY,JYZ,JZX,JZY,JZZ) 451C 452C the system of equations to be solved for obtaining the Eckart frame 453C is specified with this subroutine 454C 455 IMPLICIT NONE 456C 457 DOUBLE PRECISION ECKP(3),FVEC(3),FJAC(3,3),SINK,COSK, 458 & SINL,COSL,SINM,COSM,JXX,JXY,JXZ,JYX,JYY,JYZ,JZX,JZY,JZZ, 459 & BX,BY,BZ 460C 461 SINK = SIN(ECKP(1)) 462 COSK = COS(ECKP(1)) 463 SINL = SIN(ECKP(2)) 464 COSL = COS(ECKP(2)) 465 SINM = SIN(ECKP(3)) 466 COSM = COS(ECKP(3)) 467C 468 FVEC(1) = BX*COSM + (1 - COSM)*(JZY*COSK**2 + 469 & JXY*COSK*COSL*SINK + JYY*COSK*SINK*SINL 470 & -JZZ*COSK*SINK*SINL - JXZ*COSL*SINK**2*SINL 471 & -JYZ*SINK**2*SINL**2) 472 & +(JXZ*COSK + JXX*COSL*SINK + JXY*SINK*SINL)*SINM 473 FVEC(2) = BY*COSM + (1 - COSM)*(-(JZX*COSK**2) - 474 & JXX*COSK*COSL*SINK + JZZ*COSK*COSL*SINK 475 & +JXZ*COSL**2*SINK**2 - JYX*COSK*SINK*SINL + 476 & JYZ*COSL*SINK**2*SINL) + (JYZ*COSK + 477 & JYX*COSL*SINK + JYY*SINK*SINL)*SINM 478 FVEC(3) = BZ*COSM + (1 - COSM)*(-(JZY*COSK*COSL*SINK) 479 & - JXY*COSL**2*SINK**2 + JZX*COSK*SINK*SINL + 480 & JXX*COSL*SINK**2*SINL - JYY*COSL*SINK**2*SINL + 481 & JYX*SINK**2*SINL**2) + (JZZ*COSK + JZX*COSL*SINK 482 & +JZY*SINK*SINL)*SINM 483C 484 FJAC(1,1) = (1 - COSM)*(JXY*COSK**2*COSL - 485 & 2*JZY*COSK*SINK - JXY*COSL*SINK**2 + 486 & JYY*COSK**2*SINL - JZZ*COSK**2*SINL - 487 & 2*JXZ*COSK*COSL*SINK*SINL - 488 & JYY*SINK**2*SINL + JZZ*SINK**2*SINL - 489 & 2*JYZ*COSK*SINK*SINL**2) + 490 & (JXX*COSK*COSL - JXZ*SINK + JXY*COSK*SINL)*SINM 491 FJAC(2,1) = (1 - COSM)*(-(JXX*COSK**2*COSL) + 492 & JZZ*COSK**2*COSL + 2*JZX*COSK*SINK + 493 & 2*JXZ*COSK*COSL**2*SINK + JXX*COSL*SINK**2 494 & - JZZ*COSL*SINK**2 - 495 & JYX*COSK**2*SINL + 2*JYZ*COSK*COSL*SINK*SINL 496 & + JYX*SINK**2*SINL) + 497 & (JYX*COSK*COSL - JYZ*SINK + JYY*COSK*SINL)*SINM 498 FJAC(3,1) = (1 - COSM)*(-(JZY*COSK**2*COSL) - 499 & 2*JXY*COSK*COSL**2*SINK + JZY*COSL*SINK**2 + 500 & JZX*COSK**2*SINL + 501 & 2*JXX*COSK*COSL*SINK*SINL - 502 & 2*JYY*COSK*COSL*SINK*SINL - 503 & JZX*SINK**2*SINL + 2*JYX*COSK*SINK*SINL**2) + 504 & (JZX*COSK*COSL - JZZ*SINK + JZY*COSK*SINL)*SINM 505 FJAC(1,2) = (1 - COSM)*(JYY*COSK*COSL*SINK - 506 & JZZ*COSK*COSL*SINK - JXZ*COSL**2*SINK**2 - 507 & JXY*COSK*SINK*SINL - 508 & 2*JYZ*COSL*SINK**2*SINL + JXZ*SINK**2*SINL**2) + 509 & (JXY*COSL*SINK - JXX*SINK*SINL)*SINM 510 FJAC(2,2) = (1 - COSM)*(-(JYX*COSK*COSL*SINK) 511 & + JYZ*COSL**2*SINK**2 + JXX*COSK*SINK*SINL - 512 & JZZ*COSK*SINK*SINL - 2*JXZ*COSL*SINK**2*SINL 513 & - JYZ*SINK**2*SINL**2) + 514 & (JYY*COSL*SINK - JYX*SINK*SINL)*SINM 515 FJAC(3,2) = (1 - COSM)*(JZX*COSK*COSL*SINK + 516 & JXX*COSL**2*SINK**2 - JYY*COSL**2*SINK**2 + 517 & JZY*COSK*SINK*SINL + 518 & 2*JXY*COSL*SINK**2*SINL + 2*JYX*COSL*SINK**2*SINL - 519 & JXX*SINK**2*SINL**2 + JYY*SINK**2*SINL**2) + 520 & (JZY*COSL*SINK - JZX*SINK*SINL)*SINM 521 FJAC(1,3) = COSM*(JXZ*COSK + JXX*COSL*SINK + 522 & JXY*SINK*SINL) - BX*SINM + 523 & (JZY*COSK**2 + JXY*COSK*COSL*SINK + 524 & JYY*COSK*SINK*SINL - JZZ*COSK*SINK*SINL - 525 & JXZ*COSL*SINK**2*SINL - JYZ*SINK**2*SINL**2)*SINM 526 FJAC(2,3) = COSM*(JYZ*COSK + JYX*COSL*SINK + 527 & JYY*SINK*SINL) - BY*SINM + 528 & (-(JZX*COSK**2) - JXX*COSK*COSL*SINK + 529 & JZZ*COSK*COSL*SINK + 530 & JXZ*COSL**2*SINK**2 - JYX*COSK*SINK*SINL 531 & + JYZ*COSL*SINK**2*SINL)*SINM 532 FJAC(3,3) = COSM*(JZZ*COSK + JZX*COSL*SINK + 533 & JZY*SINK*SINL) - BZ*SINM + 534 & (-(JZY*COSK*COSL*SINK) - JXY*COSL**2*SINK**2 + 535 & JZX*COSK*SINK*SINL + 536 & JXX*COSL*SINK**2*SINL - 537 & JYY*COSL*SINK**2*SINL + JYX*SINK**2*SINL**2)*SINM 538C 539 RETURN 540 END 541C 542C 543 SUBROUTINE CLCGLD(ECKP,C) 544C 545C calculates the matrix required for transformation 546C from the input to the Eckart frame 547C The 'Goldstein' way (Classical Mechanics, p.165) and 548C thus avoiding unnecessary interchange of axes! 549C 550 IMPLICIT NONE 551C 552 DOUBLE PRECISION ECKP(3),C(3,3),NX,NY,NZ,SINPH,M1COS,COSPH 553C 554 NX = SIN( ECKP(1) ) * COS( ECKP(2) ) 555 NY = SIN( ECKP(1) ) * SIN( ECKP(2) ) 556 NZ = COS( ECKP(1) ) 557 COSPH = COS( ECKP(3) ) 558 M1COS = 1.D0 - COSPH 559 SINPH = SIN( ECKP(3) ) 560C 561 C(1,1) = COSPH + NX * NX * M1COS 562 C(2,2) = COSPH + NY * NY * M1COS 563 C(3,3) = COSPH + NZ * NZ * M1COS 564C 565 C(1,2) = NX * NY * M1COS - NZ * SINPH 566 C(1,3) = NX * NZ * M1COS + NY * SINPH 567 C(2,1) = NY * NX * M1COS + NZ * SINPH 568 C(2,3) = NY * NZ * M1COS - NX * SINPH 569 C(3,1) = NZ * NX * M1COS - NY * SINPH 570 C(3,2) = NZ * NY * M1COS + NX * SINPH 571C 572 RETURN 573 END 574C 575C 576 SUBROUTINE ECK_mnewt(ntrial,x,bvec,xjmat,tolx,tolf) 577C 578 INTEGER ntrial 579 DOUBLE PRECISION tolf,tolx,x(3) 580 INTEGER i,k,indx(3) 581 DOUBLE PRECISION errf,errx,fjac(3,3),fvec(3),p(3), 582 & bvec(3), xjmat(3,3) 583 do k=1,ntrial 584 call ECK_usrfun(x,fvec,fjac,bvec(1),bvec(2),bvec(3), 585 & xjmat(1,1),xjmat(1,2),xjmat(1,3),xjmat(2,1), 586 & xjmat(2,2),xjmat(2,3),xjmat(3,1),xjmat(3,2), 587 & xjmat(3,3)) 588 errf =0.0d0 589 do i=1,3 590 errf=errf+abs(fvec(i)) 591 enddo 592 if(errf.le.tolf)return 593 do i=1,3 594 p(i)=-fvec(i) 595 enddo 596 call ECK_lucdmpHJ(fjac,indx) 597 call ECK_lubksb(fjac,indx,p) 598 errx=0.0d0 599 do i=1,3 600 errx=errx+abs(p(i)) 601 x(i)=x(i)+p(i) 602 enddo 603 if(errx.le.tolx)return 604 enddo 605 return 606 END 607C 608C 609 SUBROUTINE ECK_lucdmpHJ(a,indx) 610C 611C Rev. Mar. 2006 hjaaj so it can treat singular 612C matrices (which occur for linear molecules) 613C 614 INTEGER indx(3) 615 DOUBLE PRECISION a(3,3),TINY 616 PARAMETER (TINY=1.0d-20) 617 INTEGER i,imax,j,k 618 DOUBLE PRECISION aamax,dum,sum,vv(3) 619#include "priunit.h" 620 do i=1,3 621 aamax=0.0d0 622 do j=1,3 623 if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 624 enddo 625CHJ if (aamax.eq.0.0d0) call quit('singular matrix in ludcmp') 626CHJ vv(i)=1.d0/aamax 627 if (aamax.le.TINY) then 628 a(i,i)=TINY 629 vv(i)=1.d0/TINY 630 else 631 vv(i)=1.d0/aamax 632 end if 633CHJ end 634 enddo 635 do j=1,3 636 do i=1,j-1 637 sum=a(i,j) 638 do k=1,i-1 639 sum=sum-a(i,k)*a(k,j) 640 enddo 641 a(i,j)=sum 642 enddo 643 aamax=0.0d0 644 imax=j 645 do i=j,3 646 sum=a(i,j) 647 do k=1,j-1 648 sum=sum-a(i,k)*a(k,j) 649 enddo 650 a(i,j)=sum 651 dum=vv(i)*abs(sum) 652 if (dum.ge.aamax) then 653 imax=i 654 aamax=dum 655 endif 656 enddo 657 if (j.ne.imax)then 658 do k=1,3 659 dum=a(imax,k) 660 a(imax,k)=a(j,k) 661 a(j,k)=dum 662 enddo 663 vv(imax)=vv(j) 664 endif 665 indx(j)=imax 666 if(a(j,j).eq.0.0d0)a(j,j)=TINY 667 if(j.ne.3)then 668 dum=1.d0/a(j,j) 669 do i=j+1,3 670 a(i,j)=a(i,j)*dum 671 enddo 672 endif 673 enddo 674 return 675 END 676C 677C 678 SUBROUTINE ECK_lubksb(a,indx,b) 679C 680C 681 INTEGER indx(3) 682 DOUBLE PRECISION a(3,3),b(3) 683 INTEGER i,ii,j,ll 684 DOUBLE PRECISION sum 685 ii=0 686 do i=1,3 687 ll=indx(i) 688 sum=b(ll) 689 b(ll)=b(i) 690 if (ii.ne.0)then 691 do j=ii,i-1 692 sum=sum-a(i,j)*b(j) 693 enddo 694 else if (sum .ne. 0.0d0) then 695 ii=i 696 endif 697 b(i)=sum 698 enddo 699 do i=3,1,-1 700 sum=b(i) 701 do j=i+1,3 702 sum=sum-a(i,j)*b(j) 703 enddo 704 b(i)=sum/a(i,i) 705 enddo 706 return 707 END 708