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 fcsini*/ 20 SUBROUTINE FCSINI 21C ************************************************************************ 22C *** This routine initializes the variables needed for a more general *** 23C *** symmetry treatment. *** 24C ************************************************************************ 25#include "implicit.h" 26C 27#include "fcsym.h" 28C 29 IF (FCLASS .EQ. 'C1 ') THEN 30C 31C *** Not supposed to do anything symmetry related. *** 32C 33 ELSE IF ((FCLASS(2:2).EQ.'i') .OR. (FCLASS(2:2).EQ.'s')) THEN 34 IF (FCLASS(2:2).EQ.'s') VPLANE = .TRUE. 35 IF (FCLASS(2:2).EQ.'i') ICNTR = .TRUE. 36 NCORDR = 1 37 NGORDR = 2 38 NGVERT = 2 39 NCVERT = 2 40 NUMELM = 1 41 NMTRX = 1 42 NONEDI = 1 43 NDEGIR = 0 44 N1DIME = 2 45 N2DIME = 0 46 ELSE 47 IF ((FCLASS(1:1).EQ.'S').OR.((FCLASS(1:1).EQ.'D').AND. 48 & (FCLASS(3:3).EQ.'d'))) THEN 49 ROTARE = .TRUE. 50 READ (FCLASS(2:2),'(I1)') NCORDR 51 NONEDI = 2 52 IF (FCLASS(3:3).EQ.'d') THEN 53 NCORDR = 2*NCORDR 54 END IF 55 NDEGIR = (NCORDR-NONEDI)/2 56 SEPDEG = .TRUE. 57 ELSE 58 MROTAX = .TRUE. 59 READ (FCLASS(2:2),'(I1)') NCORDR 60 NONEDI = 1 61 IF (MOD(NCORDR,2).EQ.0) NONEDI = 2 62 NDEGIR = (NCORDR-NONEDI)/2 63 END IF 64C 65 IF ((FCLASS(1:1).NE.'D').AND.(FCLASS(3:3) .EQ. ' ')) THEN 66 NGORDR = NCORDR 67 NGVERT = NONEDI + 4*NDEGIR 68 NCVERT = NONEDI + NDEGIR 69 NUMELM = 0 70 N1DIME = NONEDI 71 N2DIME = NDEGIR 72 ELSE IF ((FCLASS(1:1).EQ.'C').AND.(FCLASS(3:3).EQ.'h')) THEN 73 NGORDR = 2*NCORDR 74 NGVERT = 2*(NONEDI + 4*NDEGIR) 75 NCVERT = 2*(NONEDI + NDEGIR) 76 NUMELM = 1 77 N1DIME = 2*NONEDI 78 N2DIME = 2*NDEGIR 79 SEPDEG = .TRUE. 80 ELSE 81 NGORDR = 2*NCORDR 82 NUMELM = 1 83 N1DIME = 2*NONEDI 84 N2DIME = NDEGIR 85 IF ((FCLASS(1:1).EQ.'D').AND.(FCLASS(3:3).NE.'d')) THEN 86 NGORDR = 2*NGORDR 87 NUMELM = 2 88 N1DIME = 2*N1DIME 89 N2DIME = 2*N2DIME 90 END IF 91 NGVERT = NGORDR 92 NCVERT = N1DIME + N2DIME 93 END IF 94C 95 NMTRX = NUMELM + 1 96 END IF 97C 98 NIREP = N1DIME + N2DIME 99C 100 IF ((FCLASS(1:1).EQ.'D').AND.(FCLASS(3:3).NE.'d')) ROTAX2 = .TRUE. 101 IF (FCLASS(3:3).EQ.'v') VPLANE = .TRUE. 102 IF (FCLASS(3:3).EQ.'h') HPLANE = .TRUE. 103C 104 RETURN 105 END 106C 107C 108C /*Deck grpchr*/ 109 SUBROUTINE GRPCHR(COOR,SYMCOR,GRIREP,CHRCTR,WORK,ICRIRP,LWORK, 110 & IPRINT) 111C 112#include "implicit.h" 113#include "mxcent.h" 114#include "priunit.h" 115 PARAMETER (D1 = 1.0D0) 116#include "nuclei.h" 117#include "trkoor.h" 118#include "fcsym.h" 119C 120 DIMENSION COOR(NCOOR), SYMCOR(NCOOR,NCOOR), GRIREP(NGORDR,NGVERT), 121 & CHRCTR(NGORDR,NCVERT), WORK(LWORK), ICRIRP(NCOOR,2) 122C 123 IF (FCLASS .NE. 'C1 ') THEN 124 CALL HEADER('Making symmetry adapted coor in '//FCLASS(1:3),-1) 125C 126 KSIRRP = 1 127 KLAST = KSIRRP + NCORDR**2 128 LWRK = LWORK - KLAST 129 CALL CYCLIC(WORK(KSIRRP),GRIREP,CHRCTR,WORK(KLAST),LWRK,IPRINT) 130C 131 KATMTR = 1 132 KCORTR = KATMTR + NUCDEP**2*NMTRX 133 KTRMTX = KCORTR + 9 *NMTRX 134 KTMPSM = KTRMTX + 9*NUCDEP**2*NGORDR 135 KTMPTR = KTMPSM + NCOOR 136 KTMPVC = KTMPTR + 2*NCOOR**2 137 KNSTBC = KTMPVC + NCOOR 138 KLAST = KNSTBC + NUCDEP 139 LWRK = LWORK - KLAST 140 CALL MKSYMC(COOR,SYMCOR,WORK(KATMTR),WORK(KCORTR),WORK(KTRMTX), 141 & WORK(KTMPSM),GRIREP,WORK(KTMPTR),WORK(KTMPVC), 142 & WORK(KLAST),ICRIRP,WORK(KNSTBC),LWRK,IPRINT) 143 ELSE 144 CALL HEADER('There is only C1-symmetry, using cart. coor',-1) 145C 146 KDIM = NCOOR**2 147 CALL DZERO(SYMCOR,KDIM) 148 DO 100 ICOOR = 1, NCOOR 149 SYMCOR(ICOOR,ICOOR) = D1 150 ICRIRP(ICOOR,1) = 1 151 ICRIRP(ICOOR,2) = 0 152 100 CONTINUE 153 N1DIME = 1 154 END IF 155C 156 RETURN 157 END 158C 159C /*Deck cyclic*/ 160 SUBROUTINE CYCLIC(SIRREP,GRIREP,CHRCTR,WORK,LWORK,IPRINT) 161C 162#include "implicit.h" 163#include "pi.h" 164#include "priunit.h" 165 PARAMETER (NDEG = 2) 166 PARAMETER (D1 = 1.0D0) 167 PARAMETER (D0 = 0.0D0, D2 = 2.0D0) 168C 169 DIMENSION SIRREP(NCORDR,NCORDR), GRIREP(NGORDR,NGVERT), 170 & CHRCTR(NGORDR,NCVERT), WORK(LWORK) 171#include "fcsym.h" 172C 173 DO 100 J = 0, NONEDI-1 174 DO 100 I = 1, NCORDR 175 SIRREP(I,J+1) = DBLE((-1)**(J*(I-1))) 176 100 CONTINUE 177C 178 IORDR = NONEDI 179 DO 200 N = 1, NDEGIR 180 DO 200 M = 1, NDEG 181 IORDR = IORDR + 1 182 DO 300 K = 1, NCORDR 183 PHASE = D1 184 IF (M .EQ. 2) PHASE = -D1 185C 186 SIRREP(K,IORDR) = PHASE*DBLE(N*(K-1)) 187 300 CONTINUE 188 200 CONTINUE 189C 190 CONST = D2*PI/DBLE(NCORDR) 191C 192 IF ((FCLASS(1:1).EQ.'S').OR.((FCLASS(1:1).EQ.'D').AND. 193 & (FCLASS(3:3).EQ.'d'))) THEN 194 KRDREP = 1 195 CALL REDRST(WORK(KRDREP),SIRREP,NDEG,IPRINT) 196 ELSE 197 KRDREP = 1 198 CALL REDTRA(WORK(KRDREP),SIRREP,NDEG,IPRINT) 199 END IF 200C 201 IF ((FCLASS(1:1).EQ.'C').AND.(FCLASS(3:3).EQ.'h')) THEN 202 KRDREP = 1 203 CALL CHREP(GRIREP,WORK(KRDREP),SIRREP,NDEG,IPRINT) 204 ELSE 205 KRDREP = 1 206 KTMPRD = KRDREP + NCORDR*NDEGIR*NDEG**2 207 CALL GENREP(GRIREP,WORK(KRDREP),SIRREP,WORK(KTMPRD),NDEG, 208 & IPRINT) 209 END IF 210C 211C *** Calculate character table. *** 212C 213 CALL CALCHR(GRIREP,CHRCTR) 214C 215C *** Print section. *** 216C 217 IF (IPRINT .GT. 7) THEN 218 WRITE (LUPRI,'(A)') ' ' 219 WRITE (LUPRI,'(A)') 'The char. for the cyclic "basis" '//FCLASS 220 DO J = 1, NONEDI 221 WRITE (LUPRI,'(I4,A,10F6.1)') 222 & J, '. irep', (SIRREP(I,J),I=1,NCORDR) 223 WRITE (LUPRI,'(A)') ' ' 224 END DO 225C 226 IORDR = NONEDI 227 DO K = 1, NDEGIR 228 DO J = 1, NDEG 229 IORDR = IORDR + 1 230 WRITE (LUPRI,'(I4,A,10F6.1)') NONEDI+K, 231 & '. irep', (SIRREP(I,IORDR),I=1,NCORDR) 232 END DO 233 WRITE (LUPRI,'(A)') ' ' 234 END DO 235C 236 END IF 237C 238 CALL HEADER('Group matrix elements for the group '//FCLASS(1:3), 239 & -1) 240C 241 DO ID1 = 1, N1DIME 242 WRITE (LUPRI,'(I2,A)') ID1, '. irep:' 243 WRITE (LUPRI,'(24F6.2)') (GRIREP(I,ID1), I = 1, NGORDR) 244 END DO 245C 246 IORDR = N1DIME 247 DO ID2 = 1, N2DIME 248 JEXT = N1DIME + ID2 249 WRITE (LUPRI,'(I2,A)') JEXT, '. irep:' 250 DO I = 1, NDEG**2 251 IORDR = IORDR + 1 252 WRITE (LUPRI,'(24F6.2)') 253 & (GRIREP(J,IORDR), J = 1, NGORDR) 254 END DO 255 END DO 256 WRITE (LUPRI,'(A)') ' ' 257C 258 CALL HEADER('Characters for the group '// FCLASS(1:3),0) 259 DO IREP = 1, N1DIME+N2DIME 260 WRITE (LUPRI,'(I2,A)') IREP, '. irep:' 261 WRITE (LUPRI,'(24F5.1)') (CHRCTR(I,IREP), I = 1, NGORDR) 262 END DO 263 WRITE (LUPRI,'(A)') ' ' 264C 265 RETURN 266 END 267C 268C 269C /*Deck chrep*/ 270 SUBROUTINE CHREP(GRIREP,REDREP,SIRREP,NDEG,IPRINT) 271#include "implicit.h" 272 PARAMETER (D1 = 1.0D0) 273C 274#include "fcsym.h" 275 DIMENSION GRIREP(NGORDR,NGVERT), SIRREP(NCORDR,NCORDR), 276 & REDREP(NDEG,NDEG,NCORDR,NDEGIR) 277C 278 DO 100 J = 1, NONEDI 279 DO 100 I = 1, NCORDR 280 JNM = 2*(J-1) + 1 281 GRIREP(I ,JNM ) = SIRREP(I,J) 282 GRIREP(I+NCORDR,JNM ) = SIRREP(I,J) 283 GRIREP(I ,JNM+1) = SIRREP(I,J) 284 GRIREP(I+NCORDR,JNM+1) = -SIRREP(I,J) 285 100 CONTINUE 286C 287 DO 200 M = 1, 2 288 DO 200 L = 1, NDEGIR 289 DO 200 IORDR = 1, NCORDR 290 PHASE = D1 291 IF (M.EQ.2) PHASE = -D1 292 IPLACE = 2*NONEDI + 4*((M-1)*NDEGIR + (L-1)) 293C 294 DO 300 J = 1, NDEG 295 DO 300 I = 1, NDEG 296 IPLACE = IPLACE + 1 297 GRIREP(IORDR ,IPLACE) = REDREP(I,J,IORDR,L) 298 GRIREP(IORDR+NCORDR,IPLACE) = PHASE*REDREP(I,J,IORDR,L) 299 300 CONTINUE 300 200 CONTINUE 301C 302 RETURN 303 END 304C 305C 306C /*Deck redtra*/ 307 SUBROUTINE REDTRA(REDREP,SIRREP,NDEG,IPRINT) 308C 309#include "implicit.h" 310#include "pi.h" 311#include "priunit.h" 312C 313 PARAMETER (D2 = 2.0D0) 314#include "fcsym.h" 315 DIMENSION REDREP(NDEG,NDEG,NCORDR,NDEGIR), SIRREP(NCORDR,NCORDR) 316C 317 CONST = D2*PI/DBLE(NCORDR) 318C 319 DO 100 J = 1, NDEGIR 320 DO 100 I = 1, NCORDR 321 JEXT = 2*J-1 + NONEDI 322 ANGLE = CONST*SIRREP(I,JEXT) 323C 324 IF (IPRINT .GT. 20) THEN 325 WRITE (LUPRI,'(A)') 'Main rotation angle:' 326 DEGANG = 180.0D0*ANGLE/PI 327 WRITE (LUPRI,'(F10.6,F12.6)') ANGLE, DEGANG 328 END IF 329C 330 REDREP(1,1,I,J) = COS(ANGLE) 331 REDREP(1,2,I,J) = -SIN(ANGLE) 332 REDREP(2,1,I,J) = SIN(ANGLE) 333 REDREP(2,2,I,J) = COS(ANGLE) 334 100 CONTINUE 335C 336 IF (IPRINT .GT. 7) THEN 337 WRITE (LUPRI,'(/A)') 338 & 'The transformed cyclic group degenerate representations' 339 DO 200 L = 1, NDEGIR 340 DO 200 J = 1, NDEG 341 DO 200 I = 1, NDEG 342 WRITE (LUPRI,'(10F8.4)') (REDREP(I,J,K,L),K=1,NCORDR) 343 200 CONTINUE 344 END IF 345C 346 RETURN 347 END 348C 349C 350C /*Deck redrst*/ 351 SUBROUTINE REDRST(REDREP,SIRREP,NDEG,IPRINT) 352C 353#include "implicit.h" 354#include "pi.h" 355#include "priunit.h" 356C 357 PARAMETER (D2 = 2.0D0, D4 = 4.0D0) 358#include "fcsym.h" 359 DIMENSION REDREP(NDEG,NDEG,NCORDR,NDEGIR), SIRREP(NCORDR,NCORDR) 360C 361 CONST = D2*PI/DBLE(NCORDR) 362C 363 DO 100 J = 1, NDEGIR 364 DO 100 I = 1, NCORDR 365 JEXT = 2*J-1 + NONEDI 366 ANGLE = CONST*SIRREP(I,JEXT) 367C 368 IF (IPRINT .GT. 20) THEN 369 WRITE (LUPRI,'(A)') 'Main rotation angle:' 370 DEGANG = 180.0D0*ANGLE/PI 371 WRITE (LUPRI,'(F10.6,F12.6)') ANGLE, DEGANG 372 END IF 373C 374 REDREP(1,1,I,J) = COS(ANGLE) 375 REDREP(1,2,I,J) = -SIN(ANGLE) 376 REDREP(2,1,I,J) = SIN(ANGLE) 377 REDREP(2,2,I,J) = COS(ANGLE) 378 100 CONTINUE 379C 380 IF (IPRINT .GT. 7) THEN 381 WRITE (LUPRI,'(/A)') 'The transformed improper rotational' // 382 & ' group degenerate representations' 383 DO 200 L = 1, NDEGIR 384 WRITE (LUPRI,'(A,i4)') 'Degenerate irep', L 385 DO 300 J = 1, NDEG 386 DO 300 I = 1, NDEG 387 WRITE (LUPRI,'(10F8.4)') (REDREP(I,J,K,L),K=1,NCORDR) 388 300 CONTINUE 389 200 CONTINUE 390 END IF 391C 392 RETURN 393 END 394C 395C 396C 397 SUBROUTINE GENREP(GRIREP,REDREP,SIRREP,TMPRED,NDEG,IPRINT) 398C 399#include "implicit.h" 400 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, D2 = 2.0D0) 401#include "fcsym.h" 402C 403 DIMENSION GRIREP(NGORDR,NGVERT), SIRREP(NCORDR,NCORDR), 404 & REDREP(NDEG,NDEG,NCORDR,NDEGIR), 405 & TMPRED(2*NCORDR,2*NCORDR) 406C 407 IF (NUMELM .EQ. 0) THEN 408 DO 10 J = 1, NONEDI 409 DO 10 I = 1, NCORDR 410 GRIREP(I,J) = SIRREP(I,J) 411 10 CONTINUE 412C 413 IPLACE = NONEDI 414 DO 20 L = 1, NDEGIR 415 DO 20 K = 1, NDEG 416 DO 20 J = 1, NDEG 417 IPLACE = IPLACE + 1 418 DO 30 I = 1, NCORDR 419 GRIREP(I,IPLACE) = REDREP(J,K,I,L) 420 30 CONTINUE 421 20 CONTINUE 422 ELSE 423 DO 100 L = 1, 2 424 DO 100 K = 1, 2 425 ISTART = (K-1)*NCORDR 426 DO 200 J = 1, NONEDI 427 DO 200 I = 1, NCORDR 428 PHASE = D1 429 IF ((K.EQ.2).AND.(L.EQ.2)) PHASE = -D1 430C 431 INM = ISTART + I 432 JNM = NONEDI*(L-1) + J 433 GRIREP(INM,JNM) = PHASE*SIRREP(I,J) 434 200 CONTINUE 435 100 CONTINUE 436C 437 DO 300 M = 1, 2 438 ISTART = (M-1)*NCORDR 439 IPLACE = 2*NONEDI 440 DO 400 L = 1, NDEGIR 441 DO 400 K = 1, NDEG 442 DO 400 J = 1, NDEG 443 IPLACE = IPLACE + 1 444 PHASE = D1 445 IF ((M.EQ.2).AND.(K.EQ.2)) PHASE = -D1 446C 447 DO 500 I = 1, NCORDR 448 INM = ISTART + I 449 GRIREP(INM,IPLACE) = PHASE*REDREP(J,K,I,L) 450 500 CONTINUE 451 400 CONTINUE 452 300 CONTINUE 453 END IF 454C 455 IF (NUMELM .EQ. 2) THEN 456 457 DO 600 J = 1, 2*NCORDR 458 DO 600 I = 1, 2*NCORDR 459 TMPRED(I,J) = GRIREP(I,J) 460 600 CONTINUE 461C 462 DO 700 L = 1, 2 463 DO 700 K = 1, 2 464C 465 PHASE = D1 466 IF ((K.EQ.2).AND.(L.EQ.2)) PHASE = -D1 467C 468 DO 800 J = 1, 2*NONEDI 469 DO 800 I = 1, 2*NCORDR 470 INM = (L-1)*2*NCORDR+I 471 JNM = (K-1)*2*NONEDI+J 472 GRIREP(INM,JNM) = PHASE*TMPRED(I,J) 473 800 CONTINUE 474 700 CONTINUE 475C 476 ISTART = 2*NONEDI 477 DO 900 M = 1, 2 478 DO 900 L = 1, 2 479 JVAL = ISTART 480 DO 1100 K = 1, NDEGIR 481 DO 1100 J = 1, NDEG**2 482 JVAL = JVAL + 1 483 PHASE = D1 484 IF ((L.EQ.2).AND.(M.EQ.2)) PHASE = -D1 485C 486 DO 1100 I = 1, 2*NCORDR 487 INM = (M-1)*2*NCORDR + I 488 JNM = (L-1)*4*NDEGIR + 2*NONEDI + JVAL 489 GRIREP(INM,JNM) = PHASE*TMPRED(I,JVAL) 490 1100 CONTINUE 491C 492 900 CONTINUE 493C 494 END IF 495 RETURN 496 END 497C 498C 499C /*Deck calchr*/ 500 SUBROUTINE CALCHR(GRIREP,CHRCTR) 501#include "implicit.h" 502#include "priunit.h" 503C 504#include "fcsym.h" 505C 506 DIMENSION GRIREP(NGORDR,NGVERT), CHRCTR(NGORDR,NCVERT) 507C 508 DO 100 IREP = 1, N1DIME 509 DO 100 IORDR = 1, NGORDR 510 CHRCTR(IORDR,IREP) = GRIREP(IORDR,IREP) 511 100 CONTINUE 512C 513 IGPLAC = N1DIME + 1 514 DO 200 I = 1, N2DIME 515 IREP = N1DIME + I 516 DO 300 IORDR = 1, NGORDR 517 CHRCTR(IORDR,IREP) = GRIREP(IORDR,IGPLAC ) 518 & + GRIREP(IORDR,IGPLAC+3) 519 300 CONTINUE 520 IGPLAC = IGPLAC + 4 521 200 CONTINUE 522C 523 RETURN 524 END 525C 526C 527C /*Deck atmtra*/ 528 SUBROUTINE TRMTRX(COOR,ATMTRA,CORTRA,IPRINT) 529C 530#include "implicit.h" 531#include "pi.h" 532#include "mxcent.h" 533#include "priunit.h" 534 PARAMETER (D1 = 1.0D0, D2 = 2.0D0) 535 PARAMETER (DTHR = 1.0D-6) 536#include "nuclei.h" 537#include "trkoor.h" 538#include "fcsym.h" 539C 540 DIMENSION COOR(NCOOR), ATMTRA(NUCDEP,NUCDEP,NMTRX), 541 & CORTRA(3,3,NMTRX), TMPCOR(3) 542C 543 INMTRX = 0 544C 545C *** Check how atoms are connected by the main rotation *** 546C 547 IF (MROTAX) THEN 548 INMTRX = INMTRX + 1 549 CALL DZERO(CORTRA(1,1,INMTRX),9) 550C 551 ROTM = D2*PI/DBLE(NCORDR) 552C 553 CORTRA(1,1,INMTRX) = COS(ROTM) 554 CORTRA(2,1,INMTRX) = -SIN(ROTM) 555 CORTRA(1,2,INMTRX) = SIN(ROTM) 556 CORTRA(2,2,INMTRX) = COS(ROTM) 557 CORTRA(3,3,INMTRX) = D1 558C 559 CALL FNDATR(ATMTRA(1,1,INMTRX),COOR,CORTRA(1,1,INMTRX),TMPCOR, 560 & 'main rotational axis.') 561 ELSE IF (ROTARE) THEN 562 INMTRX = INMTRX + 1 563 ROTM = D2*PI/DBLE(NCORDR) 564C 565 CORTRA(1,1,INMTRX) = COS(ROTM) 566 CORTRA(2,1,INMTRX) = -SIN(ROTM) 567 CORTRA(1,2,INMTRX) = SIN(ROTM) 568 CORTRA(2,2,INMTRX) = COS(ROTM) 569 CORTRA(3,3,INMTRX) = -D1 570C 571 CALL FNDATR(ATMTRA(1,1,INMTRX),COOR,CORTRA(1,1,INMTRX),TMPCOR, 572 & 'improper rotational axis.') 573C 574 END IF 575C 576 IF (ROTAX2) THEN 577 INMTRX = INMTRX + 1 578 CALL DZERO(CORTRA(1,1,INMTRX),9) 579C 580 CORTRA(1,1,INMTRX) = D1 581 CORTRA(2,2,INMTRX) = -D1 582 CORTRA(3,3,INMTRX) = -D1 583C 584 CALL FNDATR(ATMTRA(1,1,INMTRX),COOR,CORTRA(1,1,INMTRX),TMPCOR, 585 & 'rotational axis.') 586 END IF 587C 588 IF (HPLANE) THEN 589 INMTRX = INMTRX + 1 590 CALL DZERO(CORTRA(1,1,INMTRX),9) 591C 592 CORTRA(1,1,INMTRX) = D1 593 CORTRA(2,2,INMTRX) = D1 594 CORTRA(3,3,INMTRX) = -D1 595C 596 CALL FNDATR(ATMTRA(1,1,INMTRX),COOR,CORTRA(1,1,INMTRX),TMPCOR, 597 & 'horizontal plane.') 598 END IF 599C 600 IF (VPLANE) THEN 601 INMTRX = INMTRX + 1 602 CALL DZERO(CORTRA(1,1,INMTRX),9) 603C 604 CORTRA(1,1,INMTRX) = D1 605 CORTRA(2,2,INMTRX) = -D1 606 CORTRA(3,3,INMTRX) = D1 607C 608 CALL FNDATR(ATMTRA(1,1,INMTRX),COOR,CORTRA(1,1,INMTRX),TMPCOR, 609 & 'vertical mirror plane.') 610 END IF 611C 612 IF (ICNTR) THEN 613 INMTRX = INMTRX + 1 614 CALL DZERO(CORTRA(1,1,INMTRX),9) 615C 616 CORTRA(1,1,INMTRX) = -D1 617 CORTRA(2,2,INMTRX) = -D1 618 CORTRA(3,3,INMTRX) = -D1 619C 620 CALL FNDATR(ATMTRA(1,1,INMTRX),COOR,CORTRA(1,1,INMTRX),TMPCOR, 621 & 'inversion center.') 622 END IF 623C 624 IF (IPRINT .GT. 20) THEN 625C 626 INMTRX = 0 627C 628 IF (MROTAX) THEN 629 INMTRX = INMTRX + 1 630 WRITE (LUPRI,'(A)') 'Atom transformation matrix' // 631 & ' for main rotation.' 632 DO J = 1, NUCDEP 633 WRITE (LUPRI,'(20F5.1)') (ATMTRA(I,J,INMTRX),I=1, NUCDEP) 634 END DO 635 WRITE (LUPRI,'(A)') ' ' 636 WRITE (LUPRI,'(A)') 'Coordinate transformation matrix' // 637 & ' for main rotation.' 638 DO J = 1, 3 639 WRITE (LUPRI,'(20F8.4)') (CORTRA(I,J,INMTRX),I=1,3) 640 END DO 641 WRITE (LUPRI,'(A)') ' ' 642 END IF 643C 644 IF (ROTARE) THEN 645 INMTRX = INMTRX + 1 646 WRITE (LUPRI,'(A)') 'Atom transformation matrix' // 647 & ' for improper rotation.' 648 DO J = 1, NUCDEP 649 WRITE (LUPRI,'(20F5.1)') (ATMTRA(I,J,INMTRX),I=1, NUCDEP) 650 END DO 651 WRITE (LUPRI,'(A)') ' ' 652 WRITE (LUPRI,'(A)') 'Coordinate transformation matrix' // 653 & ' for improper rotation.' 654 DO J = 1, 3 655 WRITE (LUPRI,'(20F8.4)') (CORTRA(I,J,INMTRX),I=1,3) 656 END DO 657 WRITE (LUPRI,'(A)') ' ' 658 END IF 659C 660 IF (ROTAX2) THEN 661 INMTRX = INMTRX + 1 662 WRITE (LUPRI,'(A)') 'Atom transformation matrix for' // 663 & ' 2. rotation axis' 664 DO J = 1, NUCDEP 665 WRITE (LUPRI,'(20F5.1)') (ATMTRA(I,J,INMTRX),I=1, NUCDEP) 666 END DO 667 WRITE (LUPRI,'(A)') ' ' 668C 669 WRITE (LUPRI,'(A)') 'Coordinate transformation matrix' // 670 & ' for 2. rotational axis.' 671 DO J = 1, 3 672 WRITE (LUPRI,'(20F8.4)') (CORTRA(I,J,INMTRX),I=1,3) 673 END DO 674 WRITE (LUPRI,'(A)') ' ' 675 END IF 676C 677 IF (HPLANE) THEN 678 INMTRX = INMTRX + 1 679 WRITE (LUPRI,'(A)') 'Atom transformation matrix for' // 680 & ' horizontal mirror plane' 681 DO J = 1, NUCDEP 682 WRITE (LUPRI,'(20F5.1)') (ATMTRA(I,J,INMTRX),I=1, NUCDEP) 683 END DO 684 WRITE (LUPRI,'(A)') ' ' 685C 686 WRITE (LUPRI,'(A)') 'Coordinate transformation matrix' // 687 & ' for horizontal mirror plane.' 688 DO J = 1, 3 689 WRITE (LUPRI,'(20F8.4)') (CORTRA(I,J,INMTRX),I=1,3) 690 END DO 691 WRITE (LUPRI,'(A)') ' ' 692 END IF 693C 694 IF (VPLANE) THEN 695 INMTRX = INMTRX + 1 696 WRITE (LUPRI,'(A)') 'Atom transformation matrix for' // 697 & ' vertical mirror plane.' 698 DO J = 1, NUCDEP 699 WRITE (LUPRI,'(20F5.1)') (ATMTRA(I,J,INMTRX),I=1, NUCDEP) 700 END DO 701 WRITE (LUPRI,'(A)') ' ' 702C 703 WRITE (LUPRI,'(A)') 'Coordinate transformation matrix' // 704 & ' for vertical mirror plane.' 705 DO J = 1, 3 706 WRITE (LUPRI,'(20F8.4)') (CORTRA(I,J,INMTRX),I=1,3) 707 END DO 708 WRITE (LUPRI,'(A)') ' ' 709 END IF 710C 711 IF (ICNTR) THEN 712 INMTRX = INMTRX + 1 713 WRITE (LUPRI,'(A)') 'Atom transformation matrix for' // 714 & ' inversion center' 715 DO J = 1, NUCDEP 716 WRITE (LUPRI,'(20F5.1)') (ATMTRA(I,J,INMTRX),I=1, NUCDEP) 717 END DO 718 WRITE (LUPRI,'(A)') ' ' 719C 720 WRITE (LUPRI,'(A)') 'Coordinate transformation matrix' // 721 & ' for inversion center.' 722 DO J = 1, 3 723 WRITE (LUPRI,'(20F8.4)') (CORTRA(I,J,INMTRX),I=1,3) 724 END DO 725 WRITE (LUPRI,'(A)') ' ' 726 END IF 727 728 END IF 729C 730 RETURN 731 END 732C 733C 734C /*Deck fndatr*/ 735 SUBROUTINE FNDATR(AMAT,ATMCOR,TRAMAT,TMPCOR,TEXT) 736C 737#include "implicit.h" 738#include "mxcent.h" 739#include "priunit.h" 740C 741 PARAMETER (D1 = 1.0D0, DTHR = 1.0D-4) 742C 743#include "nuclei.h" 744#include "trkoor.h" 745#include "fcsym.h" 746 CHARACTER*(*) TEXT 747 LOGICAL FOUND 748 DIMENSION AMAT(NUCDEP,NUCDEP), ATMCOR(NCOOR), TRAMAT(3,3), 749 & TMPTRA(3,3), TMPCOR(3) 750C 751 CALL DZERO(AMAT,NUCDEP**2) 752C 753 DO 200 IATOM2 = 1, NUCDEP 754 CALL DZERO(TMPCOR,3) 755C 756 ICOOR2 = 3*(IATOM2-1) 757 DO 300 IX2 = 1, 3 758 DO 300 IX1 = 1, 3 759 TMPCOR(IX2) = TMPCOR(IX2) 760 & + TRAMAT(IX1,IX2)*ATMCOR(ICOOR2+IX1) 761 300 CONTINUE 762C 763 FOUND = .FALSE. 764 DO 400 IATOM1 = 1, NUCDEP 765 ICOOR1 = 3*(IATOM1-1) 766 IF ((ABS(TMPCOR(1)-ATMCOR(ICOOR1+1)).LT.DTHR) .AND. 767 & (ABS(TMPCOR(2)-ATMCOR(ICOOR1+2)).LT.DTHR) .AND. 768 & (ABS(TMPCOR(3)-ATMCOR(ICOOR1+3)).LT.DTHR)) THEN 769 FOUND = .TRUE. 770 AMAT(IATOM2,IATOM1) = D1 771 END IF 772 400 CONTINUE 773 IF (.NOT. FOUND) THEN 774 775 WRITE (LUPRI,'(/A/A,I5)') 776 & 'Error in transformation matrix for ' // TEXT, 777 & 'Atom ',IATOM2 778 CALL OUTPUT(TRAMAT,1,3,1,3,3,3,1,LUPRI) 779 WRITE (LUPRI,'(/A,3F20.10//A)') 780 & 'Transformed coordinates:',(TMPCOR(IX2),IX2=1,3), 781 & 'All atom coordinates:' 782 CALL OUTPUT(ATMCOR,1,3,1,NUCDEP,3,NUCDEP,1,LUPRI) 783 784 CALL QUIT('You may have entered wrong group information.') 785 END IF 786 200 CONTINUE 787C 788 RETURN 789 END 790C 791C 792C /*Deck mksymc*/ 793 SUBROUTINE MKSYMC(COOR,SYMCOR,ATMTRA,CORTRA,TRAMTX,TMPSYM,GRIREP, 794 & TMPTRA,TMPVEC,WORK,ICRIRP,NSTBCT,LWORK,IPRINT) 795C 796#include "implicit.h" 797#include "mxcent.h" 798#include "priunit.h" 799 PARAMETER (D1 = 1.0D0, DTHR = 1.0D-6) 800C 801#include "nuclei.h" 802#include "trkoor.h" 803#include "fcsym.h" 804 LOGICAL EXIST 805 DIMENSION COOR(NCOOR), SYMCOR(NCOOR,NCOOR), 806 & ATMTRA(NUCDEP,NUCDEP,NMTRX), CORTRA(3,3,NMTRX), 807 & TRAMTX(NCOOR,NCOOR,NGORDR), TMPSYM(NCOOR), 808 & GRIREP(NGORDR,NGVERT), TMPTRA(NCOOR,NCOOR,2), 809 & TMPVEC(NCOOR), WORK(LWORK) 810 DIMENSION NSTBCT(NUCDEP), ICRIRP(NCOOR,2) 811C 812 CALL DZERO(SYMCOR,NCOOR**2 ) 813 CALL DZERO(TRAMTX,NCOOR**2*NGORDR) 814 CALL DZERO(ATMTRA,NUCDEP*NUCDEP*NMTRX) 815 CALL DZERO(CORTRA,9*NMTRX) 816C 817 CALL TRMTRX(COOR,ATMTRA,CORTRA,IPRINT) 818C 819 CALL FNDSTB(COOR,ATMTRA,TRAMTX,TMPTRA,TMPVEC,NSTBCT,KINDCT,IPRINT) 820C 821 CALL DZERO(TRAMTX,NCOOR**2*NGORDR) 822C 823 DO 100 IORDR = 1, NCOOR 824 TRAMTX(IORDR,IORDR,1) = D1 825 100 CONTINUE 826C 827 IF (.NOT. MROTAX .AND. .NOT. ROTARE) THEN 828C 829 DO 200 INMTRX = 1, NMTRX 830 DO 300 IATOM2 = 1, NUCDEP 831 DO 300 IATOM1 = 1, NUCDEP 832C 833 DO 400 IX2 = 1, 3 834 DO 400 IX1 = 1, 3 835 IORDR1 = 3*(IATOM1-1) + IX1 836 IORDR2 = 3*(IATOM2-1) + IX2 837C 838 TRAMTX(IORDR1,IORDR2,INMTRX+1) 839 & = CORTRA(IX1,IX2,INMTRX)*ATMTRA(IATOM1,IATOM2,INMTRX) 840 400 CONTINUE 841 300 CONTINUE 842 200 CONTINUE 843C 844 ELSE 845 DO 500 IATOM2 = 1, NUCDEP 846 DO 500 IATOM1 = 1, NUCDEP 847C 848 DO 600 IX2 = 1, 3 849 DO 600 IX1 = 1, 3 850 IORDR1 = 3*(IATOM1-1) + IX1 851 IORDR2 = 3*(IATOM2-1) + IX2 852C 853 TRAMTX(IORDR1,IORDR2,2) 854 & = CORTRA(IX1,IX2,1)*ATMTRA(IATOM1,IATOM2,1) 855 600 CONTINUE 856 500 CONTINUE 857C 858 DO 700 IROT = 2, NCORDR-1 859C 860 DO 701 ICOOR2 = 1, NCOOR 861 DO 701 ICOOR1 = 1, NCOOR 862 TMPTRA(ICOOR1,ICOOR2,2) = TRAMTX(ICOOR1,ICOOR2,2) 863 701 CONTINUE 864C 865 DO 800 ITMP = 2, IROT 866 IF (MOD(ITMP,2).EQ.0) THEN 867 IM1 = 1 868 IM2 = 2 869 ELSE 870 IM1 = 2 871 IM2 = 1 872 END IF 873C 874 CALL DZERO(TMPTRA(1,1,IM1),NCOOR**2) 875C 876 DO 900 ICOOR3 = 1, NCOOR 877 DO 900 ICOOR2 = 1, NCOOR 878 DO 900 ICOOR1 = 1, NCOOR 879 TMPTRA(ICOOR1,ICOOR3,IM1) = TMPTRA(ICOOR1,ICOOR3,IM1) 880 & + TMPTRA(ICOOR1,ICOOR2,IM2) 881 & *TRAMTX(ICOOR2,ICOOR3,2 ) 882 900 CONTINUE 883 800 CONTINUE 884C 885 DO 1100 ICOOR2 = 1, NCOOR 886 DO 1100 ICOOR1 = 1, NCOOR 887 TRAMTX(ICOOR1,ICOOR2,IROT+1) = TMPTRA(ICOOR1,ICOOR2,IM1) 888 1100 CONTINUE 889C 890 700 CONTINUE 891C 892 DO 1200 IELM = 1, NUMELM 893 MLTMAX = IELM*NCORDR 894C 895 DO 1300 IATOM2 = 1, NUCDEP 896 DO 1300 IATOM1 = 1, NUCDEP 897C 898 DO 1400 IX2 = 1, 3 899 DO 1400 IX1 = 1, 3 900 ICOOR1 = 3*(IATOM1-1) + IX1 901 ICOOR2 = 3*(IATOM2-1) + IX2 902C 903 TRAMTX(ICOOR1,ICOOR2,MLTMAX+1) 904 & = CORTRA(IX1 ,IX2 ,IELM+1) 905 & *ATMTRA(IATOM1,IATOM2,IELM+1) 906 1400 CONTINUE 907 1300 CONTINUE 908C 909 DO 1500 IOPR = 2, MLTMAX 910 CALL DZERO(TRAMTX(1,1,MLTMAX+IOPR),NCOOR**2) 911C 912 DO 1600 ICOOR3 = 1, NCOOR 913 DO 1600 ICOOR2 = 1, NCOOR 914 DO 1600 ICOOR1 = 1, NCOOR 915 TRAMTX(ICOOR1,ICOOR3,MLTMAX+IOPR) 916 & = TRAMTX(ICOOR1,ICOOR3,MLTMAX+IOPR) 917 & + TRAMTX(ICOOR1,ICOOR2,MLTMAX+1 ) 918 & *TRAMTX(ICOOR2,ICOOR3, IOPR) 919 1600 CONTINUE 920 1500 CONTINUE 921 1200 CONTINUE 922 END IF 923C 924 ISYMCO = 0 925 IPLACE = 1 926 DO 2100 IREP = 1, NIREP 927 NSDIM = 1 928 IF (IREP.GT.N1DIME) NSDIM = 2 929 DO 2101 IPL = 1, NSDIM 930 DO 2101 ICENT = 1, KINDCT 931 DO 2101 ICOOR2 = 3*(NSTBCT(ICENT)-1)+1, 3*(NSTBCT(ICENT)-1)+3 932 CALL DZERO(TMPSYM,NCOOR) 933 IF (IREP .LE. N1DIME) THEN 934 IPLACE1 = IREP 935 ELSE 936 IF (IPL.EQ.1) THEN 937 IPLACE1 = (N1DIME + 1) + 4*(IREP - N1DIME - 1) 938 IPLACE2 = IPLACE1 + 1 939 ELSE 940 IPLACE1 = N1DIME + 4*(IREP - N1DIME) 941 IPLACE2 = IPLACE1 - 1 942 END IF 943 END IF 944C 945 DO 2200 ICOOR1 = 1, NCOOR 946 DO 2200 IGORDR = 1, NGORDR 947 TMPSYM(ICOOR1) = TMPSYM(ICOOR1) 948 & + GRIREP(IGORDR,IPLACE1) 949 & *TRAMTX(ICOOR1,ICOOR2,IGORDR) 950 2200 CONTINUE 951C 952 RNORM2 = 0.0D0 953 DO 2300 ICOOR1 = 1, NCOOR 954 RNORM2 = RNORM2 + (TMPSYM(ICOOR1))**2 955 2300 CONTINUE 956 IF (RNORM2 .LT. DTHR) GOTO 2700 957C 958 DO 2350 ICOOR1 = 1, NCOOR 959 TMPSYM(ICOOR1) = TMPSYM(ICOOR1)/SQRT(RNORM2) 960 2350 CONTINUE 961C 962C *** Orthogonalize it. *** 963C 964 EXIST = .FALSE. 965 CALL ORTSCP(TMPSYM,SYMCOR,ICRIRP,ISYMCO,IREP,IPL,IPRINT,EXIST) 966C 967 DO 2400 IPHASE = 1, 2 968 DO 2400 ICOO2 = 1, ISYMCO 969 IF (.NOT.EXIST) THEN 970 EXIST = .TRUE. 971 PHASE = D1 972 IF (IPHASE.EQ.2) PHASE = -D1 973C 974 DO 2500 ICOO1 = 1, NCOOR 975 DIFF2 = (PHASE*TMPSYM(ICOO1) 976 & - SYMCOR(ICOO1,ICOO2))**2 977C 978 IF (DIFF2 .GT. DTHR) EXIST = .FALSE. 979 2500 CONTINUE 980C 981 END IF 982 2400 CONTINUE 983 IF (EXIST) GOTO 2700 984C 985C *** We have now found a new symmetry coordinate, and find *** 986C *** an address for it. *** 987C 988 ISYMCO = ISYMCO + 1 989C 990C *** If the projection operator for the second row is used *** 991C *** then we need to put the symmetry coordinate 1 index on *** 992C 993 ICRIRP(ISYMCO,2) = IPL-1 994 ICRIRP(ISYMCO,1) = IREP 995 DO ICOOR1 = 1, NCOOR 996 SYMCOR(ICOOR1,ISYMCO) = TMPSYM(ICOOR1) 997 END DO 998C 999C *** Using the shift function to create a partner function. *** 1000C 1001 IF (IREP .GT. N1DIME) THEN 1002 CALL DZERO(TMPSYM,NCOOR) 1003C 1004 DO 2800 IGORDR = 1, NGORDR 1005 DO 2800 ICOO2 = 1, NCOOR 1006 CALL DZERO(TMPVEC,NCOOR) 1007 DO ICOO1 = 1, NCOOR 1008 TMPVEC(ICOO2) = TMPVEC(ICOO2) 1009 & + TRAMTX(ICOO1,ICOO2,IGORDR) 1010 & *SYMCOR(ICOO1,ISYMCO) 1011 END DO 1012C 1013 TMPSYM(ICOO2) = TMPSYM(ICOO2) 1014 & + GRIREP(IGORDR,IPLACE2) 1015 & *TMPVEC(ICOO2) 1016 2800 CONTINUE 1017C 1018 RNORM2 = 0.0D0 1019 DO ICOOR1 = 1, NCOOR 1020 RNORM2 = RNORM2 + TMPSYM(ICOOR1)**2 1021 END DO 1022C 1023C *** Normalizing partner geometry. *** 1024C 1025 DO ICOOR1 = 1, NCOOR 1026 TMPSYM(ICOOR1) = TMPSYM(ICOOR1)/SQRT(RNORM2) 1027 END DO 1028C 1029C *** Orthogonalizing partner coordinate. *** 1030C 1031 CALL ORTSCP(TMPSYM,SYMCOR,ICRIRP,ISYMCO,IREP,IPL,IPRINT, 1032 & EXIST) 1033C 1034 ISYMCO = ISYMCO + 1 1035C 1036 ICRIRP(ISYMCO,1) = IREP 1037 ICRIRP(ISYMCO,2) = 1 1038 IF (IPL.EQ.2) ICRIRP(ISYMCO,2) = 0 1039C 1040 DO ICOOR1 = 1, NCOOR 1041 SYMCOR(ICOOR1,ISYMCO) = TMPSYM(ICOOR1) 1042 END DO 1043C 1044 END IF 1045C 1046 2700 CONTINUE 1047 2101 CONTINUE 1048 2100 CONTINUE 1049C 1050C 1051 IF (N1DIME.LT.NIREP) THEN 1052 KTMPCR = 1 1053 KICTMP = KTMPCR + NCOOR**2 1054 CALL SRTSCR(SYMCOR,WORK(KTMPCR),ICRIRP,WORK(KICTMP),IPRINT) 1055 END IF 1056C 1057 NUMTIM = (NCOOR-1)/10 + 1 1058 CALL HEADER ('Symmetry coordinates',-1) 1059 DO I = 1, NUMTIM 1060 ILEFT = NCOOR - 10*(I-1) 1061 ISTRT = 10*(I-1) + 1 1062 IEND =(ISTRT-1) + MIN(ILEFT,10) 1063 WRITE (LUPRI,'(A,I4,10I8)') 'Symmetry', 1064 & (ICRIRP(ICOOR2,1),ICOOR2=ISTRT,IEND) 1065 DO ICOOR1 = 1, NCOOR 1066 WRITE (LUPRI,'(A,10F8.4)') ' ', 1067 & (SYMCOR(ICOOR1,ICOOR2),ICOOR2=ISTRT,IEND) 1068 END DO 1069 WRITE (LUPRI,'(/)') 1070 END DO 1071 WRITE (LUPRI,'()') 1072C 1073 IF (IPRINT .GT. 20) THEN 1074 WRITE (LUPRI,'(/A/)') 'Transformation matrices for the ' // 1075 & 'different symmetry operations.' 1076 DO IGORDR = 1, NGORDR 1077 WRITE (LUPRI,'(A,I4)') ' Matrix number: ', IGORDR 1078 DO I = 1, NUMTIM 1079 ILEFT = NCOOR - 18*(I-1) 1080 ISTRT = 18*(I-1) + 1 1081 IEND = (ISTRT-1) + MIN(ILEFT,18) 1082 DO IORDR2 = 1, NCOOR 1083 WRITE (LUPRI,'(18F8.4)') 1084 & (TRAMTX(IORDR1,IORDR2,IGORDR),IORDR1=ISTRT,IEND) 1085 END DO 1086 WRITE (LUPRI,'(/)') 1087 END DO 1088 WRITE (LUPRI,'()') 1089 END DO 1090 END IF 1091C 1092 IF (ISYMCO.LT.NCOOR) CALL QUIT('Unable to make enough sym. coor') 1093C 1094C *** Orthogonality test for coordinates. *** 1095C 1096 ITEST = 0 1097 DO ICOOR1 = 1, NCOOR 1098 DO ICOOR2 = 1, ICOOR1-1 1099 GCCNST = 0.0D0 1100 DO ICOOR3 = 1, NCOOR 1101 GCCNST = GCCNST + SYMCOR(ICOOR3,ICOOR1) 1102 & *SYMCOR(ICOOR3,ICOOR2) 1103 END DO 1104C 1105 IF (ABS(GCCNST) .GT. DTHR) THEN 1106 ITEST = ITEST + 1 1107 WRITE (LUPRI,'(5X,A,2I4)') 1108 & 'ERROR, Non orthogonal symmetry coordinates:', ICOOR1, ICOOR2 1109 END IF 1110 END DO 1111 END DO 1112 IF (ITEST.GT.0) CALL QUIT('Problems with symmetry coordinates,' // 1113 & 'numerical differentiation. Please contact dalton admin.') 1114C 1115 RETURN 1116 END 1117C 1118C 1119C /*Deck fndstb*/ 1120 SUBROUTINE FNDSTB(COOR,ATMTRA,TRAMTX,TMPTRA,TMPVEC,NSTBCT,KINDCT, 1121 & IPRINT) 1122C 1123#include "implicit.h" 1124#include "mxcent.h" 1125#include "priunit.h" 1126 PARAMETER (D1 = 1.0D0, DTHR = 1.0D-6) 1127#include "nuclei.h" 1128#include "trkoor.h" 1129#include "fcsym.h" 1130C 1131 LOGICAL EXIST 1132 DIMENSION ATMTRA(NUCDEP,NUCDEP,NMTRX), TMPTRA(NCOOR,NCOOR,2), 1133 & TRAMTX(NCOOR,NCOOR,NGORDR), NSTBCT(NUCDEP), 1134 & TMPVEC(NCOOR), COOR(NCOOR) 1135C 1136 DO 100 IATOM = 1, NUCDEP 1137 TRAMTX(IATOM,IATOM,1) = D1 1138 100 CONTINUE 1139C 1140 IF (.NOT. MROTAX .AND. .NOT. ROTARE) THEN 1141C 1142 DO 200 INMTRX = 1, NMTRX 1143 DO 200 IATOM2 = 1, NUCDEP 1144 DO 200 IATOM1 = 1, NUCDEP 1145 TRAMTX(IATOM1,IATOM2,INMTRX+1)= ATMTRA(IATOM1,IATOM2,INMTRX) 1146 200 CONTINUE 1147C 1148 ELSE 1149 DO 300 IATOM2 = 1, NUCDEP 1150 DO 300 IATOM1 = 1, NUCDEP 1151 TRAMTX(IATOM1,IATOM2,2) = ATMTRA(IATOM1,IATOM2,1) 1152 300 CONTINUE 1153C 1154 DO 400 IROT = 2, NCORDR-1 1155C 1156 DO 500 IATOM2 = 1, NCOOR 1157 DO 500 IATOM1 = 1, NCOOR 1158 TMPTRA(IATOM1,IATOM2,2) = TRAMTX(IATOM1,IATOM2,2) 1159 500 CONTINUE 1160C 1161 DO 600 ITMP = 2, IROT 1162 IF (MOD(ITMP,2).EQ.0) THEN 1163 IM1 = 1 1164 IM2 = 2 1165 ELSE 1166 IM1 = 2 1167 IM2 = 1 1168 END IF 1169C 1170 CALL DZERO(TMPTRA(1,1,IM1),NCOOR**2) 1171C 1172 DO 700 IATOM3 = 1, NUCDEP 1173 DO 700 IATOM2 = 1, NUCDEP 1174 DO 700 IATOM1 = 1, NUCDEP 1175 TMPTRA(IATOM1,IATOM3,IM1) = TMPTRA(IATOM1,IATOM3,IM1) 1176 & + TMPTRA(IATOM1,IATOM2,IM2) 1177 & *TRAMTX(IATOM2,IATOM3,2 ) 1178 700 CONTINUE 1179 600 CONTINUE 1180C 1181 DO 800 IATOM2 = 1, NCOOR 1182 DO 800 IATOM1 = 1, NCOOR 1183 TRAMTX(IATOM1,IATOM2,IROT+1) = TMPTRA(IATOM1,IATOM2,IM1) 1184 800 CONTINUE 1185C 1186 400 CONTINUE 1187C 1188 DO 900 IELM = 1, NUMELM 1189 MLTMAX = IELM*NCORDR 1190C 1191 DO 1100 IATOM2 = 1, NUCDEP 1192 DO 1100 IATOM1 = 1, NUCDEP 1193C 1194 TRAMTX(IATOM1,IATOM2,MLTMAX+1) 1195 & = ATMTRA(IATOM1,IATOM2,IELM+1) 1196 1100 CONTINUE 1197C 1198 DO 1200 IOPR = 2, MLTMAX 1199 CALL DZERO(TRAMTX(1,1,MLTMAX+IOPR),NCOOR**2) 1200C 1201 DO 1300 IATOM3 = 1, NUCDEP 1202 DO 1300 IATOM2 = 1, NUCDEP 1203 DO 1300 IATOM1 = 1, NUCDEP 1204 TRAMTX(IATOM1,IATOM3,MLTMAX+IOPR) 1205 & = TRAMTX(IATOM1,IATOM3,MLTMAX+IOPR) 1206 & + TRAMTX(IATOM1,IATOM2,MLTMAX+1 ) 1207 & *TRAMTX(IATOM2,IATOM3, IOPR) 1208 1300 CONTINUE 1209 1200 CONTINUE 1210 900 CONTINUE 1211 END IF 1212C 1213 KINDCT = 0 1214 CALL IZERO(NSTBCT,NUCDEP) 1215 DO 1400 IATOM = 1, NUCDEP 1216C 1217 EXIST = .FALSE. 1218 DO 1500 IORDR = 1, NGORDR 1219 DO 1500 IINDCT = 1, KINDCT 1220 IF (ABS(TRAMTX(NSTBCT(IINDCT),IATOM,IORDR)).GT.DTHR) THEN 1221 GOTO 1600 1222 END IF 1223 1500 CONTINUE 1224C 1225 KINDCT = KINDCT + 1 1226 NSTBCT(KINDCT) = IATOM 1227C 1228 1600 CONTINUE 1229 1400 CONTINUE 1230C 1231c DO 1700 IINDCT = 1, KINDCT 1232c IF ((FCLASS(1:1).EQ.'S').OR.(FCLASS(3:3).EQ.'d')) THEN 1233c IYCOOR = 3*(NSTBCT(INDCT)-1) + 1 1234c ELSE 1235c IYCOOR = 3*(NSTBCT(IINDCT)-1) + 2 1236c END IF 1237cC 1238c IF (ABS(COOR(IYCOOR)) .LT. DTHR) THEN 1239cC 1240c DO 1800 IORDR = 2, NGORDR 1241cC 1242c DO 1900 IATOM = 1, NUCDEP 1243c TMPVEC(IATOM) = TRAMTX(NSTBCT(IINDCT),IATOM,IORDR) 1244c IF (ABS(TMPVEC(IATOM)).GT.DTHR) 1245c & NSTBCT(IINDCT) = IATOM 1246c 1900 CONTINUE 1247cC 1248c IYCOOR = 3*(NSTBCT(IINDCT)-1) + 2 1249c IF (ABS(COOR(IYCOOR)).GT.DTHR) GOTO 2100 1250c 1800 CONTINUE 1251c END IF 1252cC 1253c 2100 CONTINUE 1254cC 1255c 1700 CONTINUE 1256C 1257 IF (IPRINT .GT. 20) THEN 1258 WRITE (LUPRI,'(A)') 'Atom transformation matrices:' 1259C 1260 DO IMTX = 1, NGORDR 1261 WRITE (LUPRI,'(A,I2)') 'Operation ', IMTX 1262 DO I = 1, NUCDEP 1263 WRITE (LUPRI,'(12F5.1)') (TRAMTX(I,J,IMTX),J=1, NUCDEP) 1264 END DO 1265 WRITE (LUPRI,'(A)') ' ' 1266 END DO 1267C 1268 WRITE (LUPRI,'(A)') 'Symmetry independent centers:' 1269 DO IINDCT = 1, KINDCT 1270 WRITE (LUPRI,'(A,I4)') ' ', NSTBCT(IINDCT) 1271 END DO 1272 END IF 1273C 1274 RETURN 1275 END 1276C 1277C 1278C /*Deck fcscrn*/ 1279 SUBROUTINE FCSCRN(GRIREP,WORK,KDPMTX,INDSTP,ICRIRP,IRPIND,IDXTMP, 1280 & IDDBTP,IRPDEG,IREPST,NPRTNR,LWORK,NLDPMX,LDPMTX, 1281 & IFRSTD,IORDR,IRSRDR,MAXINR,IINNER,NMPRTN,IPRINT, 1282 & CLNRGY,PRTNR,ALRCAL,SCND) 1283#include "implicit.h" 1284#include "mxcent.h" 1285#include "priunit.h" 1286C 1287#include "trkoor.h" 1288#include "numder.h" 1289#include "fcsym.h" 1290#include "pvibav.h" 1291 LOGICAL CLNRGY, PRTNR, SCND, C2NRGY, DEPFC, FOUND, ALRCAL 1292 DIMENSION GRIREP(NGORDR,NGVERT), WORK(LWORK), 1293 & INDSTP(NTORDR), KDPMTX(LDPMTX,NSTRDR,IFRSTD), 1294 & IRPDEG(NMORDR), IREPST(NMORDR), ICRIRP(NCOOR,2), 1295 & IRPIND(NMORDR), IDXTMP(NMORDR), IDDBTP(NMORDR), 1296 & NPRTNR(MAXINR) 1297C 1298 IF ((FCLASS .NE. 'C1 ').AND.(NAORDR.LT.1).AND..NOT.CNMPRP) THEN 1299C 1300C ***Check symmetry for the original component *** 1301C 1302 KKIRPD = 1 1303 KIDDEG = KKIRPD + NMORDR 1304 KIDEGI = KIDDEG + NMORDR 1305C 1306 C2NRGY = .FALSE. 1307 CALL SYMMLT(GRIREP,INDSTP,ICRIRP,IREPST,IRPDEG,WORK(KKIRPD), 1308 & WORK(KIDDEG),WORK(KIDEGI),NMORDR,NTORDR,IRSRDR, 1309 & IORDR,IPRINT,C2NRGY) 1310C 1311C *** Checking whether this is a dependent force constant *** 1312C *** (as a result of using groups with degenerate ireps) *** 1313C 1314 DEPFC = .FALSE. 1315 IF (C2NRGY) THEN 1316 CALL FFCDEP(KDPMTX,INDSTP,IDXTMP,NLDPMX,LDPMTX,IFRSTD, 1317 & IRSRDR,IORDR,0,IPRINT,DEPFC) 1318 END IF 1319C 1320C *** Final test for this component *** 1321C 1322 CLNRGY = C2NRGY.AND..NOT.DEPFC 1323C 1324C *** Check if it is possible to construct a force-constant *** 1325C *** with two additional equal and arbitary indices. This *** 1326C *** force constant will contain this energy value. *** 1327C 1328 IF (((NMORDR-IORDR).GE.2).AND.(.NOT.CLNRGY) 1329 & .AND.(.NOT.SCND )) THEN 1330 ISTRT = IRSRDR + 1 1331 INUM = (NMORDR-IORDR)/2 1332C 1333 NSTP = NDCOOR 1334 IDDBTP(1) = 0 1335 DO 100 II = 2, INUM 1336 NSTP = NSTP*(NDCOOR+II-1) 1337 IDDBTP(II) = 1 1338 100 CONTINUE 1339C 1340 DO 200 ISTP = 1, NSTP 1341 IF (.NOT.CLNRGY) THEN 1342 FOUND = .FALSE. 1343 DO 300 II = INUM-1, 1, -1 1344 IF ((IDDBTP(II).GT.IDDBTP(II+1)) 1345 & .AND.(.NOT.FOUND)) THEN 1346 FOUND = .TRUE. 1347 IDDBTP(II+1) = IDDBTP(II+1) + 1 1348 END IF 1349 300 CONTINUE 1350 IF (.NOT.FOUND) IDDBTP(1) = IDDBTP(1) + 1 1351C 1352 IDX = 0 1353 DO 400 IJ = 1, INUM 1354 DO 400 II = 1, 2 1355 IDX = IDX + 1 1356 INDSTP(ISTRT+IDX) = IDDBTP(IJ) 1357 400 CONTINUE 1358C 1359 C2NRGY = .FALSE. 1360 CALL SYMMLT(GRIREP,INDSTP,ICRIRP,IREPST,IRPDEG, 1361 & WORK(KKIRPD),WORK(KIDDEG),WORK(KIDEGI), 1362 & NMORDR,NTORDR,IRSRDR+2*INUM,IORDR+2*INUM, 1363 & IPRINT,C2NRGY) 1364C 1365C *** Checking if this is a dependent component. *** 1366C 1367 DEPFC = .FALSE. 1368 IF (C2NRGY) THEN 1369 CALL FFCDEP(KDPMTX,INDSTP,IDXTMP,NLDPMX,LDPMTX, 1370 & IFRSTD,IRSRDR,IORDR,2*INUM,IPRINT, 1371 & DEPFC) 1372 END IF 1373C 1374C *** Final test for this component. *** 1375C 1376 CLNRGY = C2NRGY.AND..NOT.DEPFC 1377 END IF 1378 200 CONTINUE 1379C 1380 END IF 1381C 1382C *** Check symmetry for use of component in higher derivatives *** 1383C *** MIN(NMORDR-IORDR,IORDR) -> if calculating derivatives of *** 1384C *** order > 2*iordr then clnrgy = .true. *** 1385C 1386 IF ((IORDR.LT.NMORDR) .AND. .NOT.CLNRGY .AND. .NOT.SCND) THEN 1387 DO 500 IN = 1, NMORDR-IORDR 1388C 1389 NRP = 1 1390 DO 600 IRP = 1, IN 1391 NRP = NRP*(IRSRDR + IRP)/IRP 1392 600 CONTINUE 1393C 1394C *** Which components do we have to add *** 1395C 1396 CALL IZERO(IRPIND,NMORDR) 1397 IF (IORDR.EQ.1) THEN 1398 IRPIND(1) = 1 1399 ELSE 1400 IRPIND(1) = 0 1401 END IF 1402 DO 700 I = 2, IN 1403 IRPIND(I) = 1 1404 700 CONTINUE 1405C 1406 DO 800 IRP = 1, NRP 1407C 1408 IF (IORDR.GT.1) THEN 1409 IF (IN.GT.1) THEN 1410 FOUND = .FALSE. 1411 DO 900 IIN = 1, IN 1412 IF ((IRPIND(IIN).LE.IRPIND(IIN+1)).AND. 1413 & (.NOT.FOUND)) THEN 1414 FOUND = .TRUE. 1415 IRPIND(IIN) = IRPIND(IIN) + 1 1416 END IF 1417 900 CONTINUE 1418 ELSE 1419 IRPIND(1) = IRPIND(1) + 1 1420 END IF 1421 END IF 1422C 1423C *** Assign the new components already in the array, *** 1424C *** and screen again. *** 1425C 1426 IF (.NOT. CLNRGY) THEN 1427 ISTRT = IRSRDR + 1 1428 DO 1100 IIND = 1, IN 1429 INDSTP(ISTRT+IIND) = INDSTP(IRPIND(IIND)) 1430 1100 CONTINUE 1431C 1432 C2NRGY = .FALSE. 1433 CALL SYMMLT(GRIREP,INDSTP,ICRIRP,IREPST,IRPDEG, 1434 & WORK(KKIRPD),WORK(KIDDEG),WORK(KIDEGI), 1435 & NMORDR,NTORDR,IRSRDR+IN,IORDR+IN, 1436 & IPRINT,C2NRGY) 1437C 1438C *** Checking if this is a dependent component. *** 1439C 1440 DEPFC = .FALSE. 1441 IF (C2NRGY) THEN 1442 CALL FFCDEP(KDPMTX,INDSTP,IDXTMP,NLDPMX,LDPMTX, 1443 & IFRSTD,IRSRDR,IORDR,IN,IPRINT,DEPFC) 1444 END IF 1445C 1446C *** Final test for this component. *** 1447C 1448 CLNRGY = C2NRGY.AND..NOT.DEPFC 1449 END IF 1450C 1451C *** Check if it is possible to construct a force-constant *** 1452C *** with two additional equal and arbitary indices. This *** 1453C *** force constant will contain this energy value. *** 1454C 1455 IF (((NMORDR-(IORDR+IN)).GE.2).AND.(.NOT.CLNRGY)) THEN 1456C 1457 ISTRT = IRSRDR + 1 + IN 1458 INUM = (NMORDR-(IORDR+IN))/2 1459C 1460 NSTP = NDCOOR 1461 IDDBTP(1) = 0 1462 DO 1200 II = 2, INUM 1463 NSTP = NSTP*(NDCOOR+II-1) 1464 IDDBTP(II) = 1 1465 1200 CONTINUE 1466C 1467 DO 1300 ISTP = 1, NSTP 1468 IF (.NOT.CLNRGY) THEN 1469 FOUND = .FALSE. 1470 DO 1400 II = INUM-1, 1, -1 1471 IF ((IDDBTP(II).GT.IDDBTP(II+1)) 1472 & .AND.(.NOT.FOUND)) THEN 1473 FOUND = .TRUE. 1474 IDDBTP(II+1) = IDDBTP(II+1) + 1 1475 END IF 1476 1400 CONTINUE 1477 IF (.NOT.FOUND) IDDBTP(1) = IDDBTP(1) + 1 1478C 1479 IDX = 0 1480 DO 1500 IJ = 1, INUM 1481 DO 1500 II = 1, 2 1482 IDX = IDX + 1 1483 INDSTP(ISTRT+IDX) = IDDBTP(IJ) 1484 1500 CONTINUE 1485C 1486 C2NRGY = .FALSE. 1487 CALL SYMMLT(GRIREP,INDSTP,ICRIRP,IREPST, 1488 & IRPDEG,WORK(KKIRPD),WORK(KIDDEG), 1489 & WORK(KIDEGI),NMORDR,NTORDR, 1490 & IRSRDR+IN+2*INUM,IORDR+IN+2*INUM, 1491 & IPRINT,C2NRGY) 1492C 1493C *** Checking if this is a dependent component. *** 1494C 1495 DEPFC = .FALSE. 1496 IF (C2NRGY) THEN 1497 CALL FFCDEP(KDPMTX,INDSTP,IDXTMP,NLDPMX, 1498 & LDPMTX,IFRSTD,IRSRDR,IORDR, 1499 & IN+2*INUM,IPRINT,DEPFC) 1500 END IF 1501C 1502C *** Final test for this component. *** 1503C 1504 CLNRGY = C2NRGY.AND..NOT.DEPFC 1505 END IF 1506 1300 CONTINUE 1507C 1508 END IF 1509C 1510 800 CONTINUE 1511 500 CONTINUE 1512 END IF 1513C 1514C *** Check if further simplyfication is possible through *** 1515C *** the use of partner geometry's *** 1516C 1517 IF (CLNRGY.AND.(NAORDR.EQ.0).AND.(.NOT.SCND)) THEN 1518 KIDCTP = 1 1519 KINDTP = KIDCTP + NMORDR 1520 KLAST = KINDTP + 5 1521 LWRK = LWORK - KLAST + 1 1522 CALL PRTNRP(GRIREP,WORK(KLAST),ICRIRP,INDSTP,NPRTNR,MAXINR, 1523 & IINNER,IORDR,IRSRDR,NMPRTN,WORK(KIDCTP), 1524 & WORK(KINDTP),LWRK,CLNRGY,PRTNR,ALRCAL) 1525 END IF 1526C 1527 ELSE 1528 CLNRGY = .TRUE. 1529 END IF 1530C 1531c if (iordr.eq.2) then 1532c if ((icrirp(indstp(1),1).eq.5).and. 1533c & (icrirp(indstp(1),2).eq.1).and. 1534c & (icrirp(indstp(2),1).eq.5).and. 1535c & (icrirp(indstp(2),2).eq.0)) 1536c & stop ' ' 1537c end if 1538 RETURN 1539 END 1540C 1541C 1542C /*Deck symmlt*/ 1543 SUBROUTINE SYMMLT(GRIREP,INDSTP,ICRIRP,IREPST,IRPDEG,KIRPDG, 1544 & IDDEGI,IDEGID,NMORDR,NTORDR,IRSRDR,IORDR,IPRINT, 1545 & CLNRGY) 1546#include "implicit.h" 1547#include "mxcent.h" 1548#include "priunit.h" 1549 PARAMETER (D0 = 0.0D0, DTHR = 1.0D-12) 1550#include "trkoor.h" 1551#include "fcsym.h" 1552 LOGICAL CLNRGY 1553 DIMENSION GRIREP(NGORDR,NGVERT) 1554 DIMENSION INDSTP(NTORDR), KIRPDG(NMORDR), IRPDEG(NMORDR), 1555 & IREPST(NMORDR), IDEGID(NMORDR), IDDEGI(NMORDR), 1556 & ICRIRP(NCOOR,2) 1557C 1558C *** Starting points for the relevant rows in the relevant *** 1559C *** irreps For instance, the starting point for the 1. row *** 1560C *** of E in C3v is 2, while the second row is 4. *** 1561C *** IREPST -> Starting point of irrep *** 1562C *** KIRPDG -> dimension of the irrep *** 1563C *** IDEGGI -> coordinates that belong to 2 dim. irreps *** 1564C 1565 NMAX = 0 1566 IC = 0 1567 DO 100 IC2 = 1, IRSRDR+1 1568 NC1 = 1 1569 IF (IC2.EQ.1) NC1 = IORDR-IRSRDR 1570 DO 200 IC1 = 1, NC1 1571 IC = IC + 1 1572 IREPST(IC) = IRPSTR(GRIREP,ICRIRP,NCOOR,INDSTP(IC2),IPRINT) 1573 IF (ICRIRP(INDSTP(IC2),1) .GT. N1DIME) THEN 1574 NMAX = NMAX + 1 1575 KIRPDG(IC) = 2 1576 IDDEGI(NMAX) = IC 1577 ELSE 1578 KIRPDG(IC) = 1 1579 END IF 1580 200 CONTINUE 1581 100 CONTINUE 1582C 1583C *** IRPDEG -> The row number we want use in our product - *** 1584C *** IREPST (irep-start). For two 2-dimensional irreps this *** 1585C *** is 11, 12, 21 and 22 *** 1586C 1587 DO 300 IC = 1, IORDR 1588 IRPDEG(IC) = 1 1589 300 CONTINUE 1590C 1591C *** Number of 2 dimensional irreps. *** 1592C 1593 DO 400 IMAX = 0, NMAX 1594 INUM = 1 1595 IDEN = 1 1596 DO 500 I = 1, IMAX 1597 INUM = INUM*(NMAX-I+1) 1598 IDEN = IDEN*I 1599 500 CONTINUE 1600 NPOSBL = INUM/IDEN 1601C 1602C *** Number of possible permutations for IMAX values of 2 *** 1603C 1604 DO 600 IPOSBL = 1, NPOSBL 1605C 1606 DO 700 IC = 1, IORDR 1607 IRPDEG(IC) = 1 1608 700 CONTINUE 1609C 1610C *** IMAX .gt. 0 -> at least 1 two dimensional irrep. *** 1611C *** IDEGID -> which values of IDDEGID is 2 *** 1612C 1613 IF (IMAX .GT. 0) THEN 1614 IF (IPOSBL.EQ.1) THEN 1615 DO 800 IRDR = 1, IMAX 1616 IDEGID(IRDR) = IMAX - IRDR + 1 1617 800 CONTINUE 1618 ELSE 1619 DO 900 IRDR = 1, IMAX 1620 IF (IDEGID(IRDR).LE.NMAX-IRDR) THEN 1621 IDEGID(IRDR) = IDEGID(IRDR) + 1 1622 GOTO 1100 1623 END IF 1624 900 CONTINUE 1625 END IF 1626 END IF 1627C 1628 1100 CONTINUE 1629C 1630C *** Assigning the nesscecary values of 2 *** 1631C 1632 DO 1200 IRDR = 1, IMAX 1633 IRPDEG(IDDEGI(IDEGID(IRDR))) = 2 1634 1200 CONTINUE 1635C 1636C *** Calculating the screening condition -> P. Taylor, *** 1637C *** Lecture notes in quantum chemistry 56, Springer *** 1638C 1639 SMTXEL = D0 1640 DO 2100 IORDER = 1, NGORDR 1641C 1642C *** Finding the starting point for the *** 1643C *** screening needed for the analytical *** 1644C *** derivative. Output: RTMPCO *** 1645C 1646 CALL ANLSYM(GRIREP,ICRIRP,INDSTP,IREPST,RTMPCO,NCOOR, 1647 & IORDER,IORDR,IPRINT) 1648C 1649 DO 2200 IC = 1, IORDR 1650C 1651C *** Special care for groups based on separable *** 1652C *** degenerate groups, Sn, Cnh Dnd. GOT doesn't *** 1653C *** apply, only LOT. *** 1654C 1655 IF (SEPDEG) THEN 1656 IFAC = 2 1657 IF (ICRIRP(INDSTP(IC),2).EQ.1) IFAC = -2 1658 FMULT = GRIREP(IORDER,IREPST(IC)+IRPDEG(IC) ) 1659 & + GRIREP(IORDER,IREPST(IC)+IRPDEG(IC)+IFAC) 1660 ELSE 1661 FMULT = GRIREP(IORDER,IREPST(IC)+IRPDEG(IC)) 1662 END IF 1663 RTMPCO = RTMPCO*FMULT 1664 2200 CONTINUE 1665C 1666 SMTXEL = SMTXEL + RTMPCO 1667 2100 CONTINUE 1668C 1669 IF (ABS(SMTXEL) .GT. DTHR) CLNRGY = .TRUE. 1670 600 CONTINUE 1671 400 CONTINUE 1672C 1673 RETURN 1674 END 1675C 1676C 1677C /* Deck irpstr */ 1678 FUNCTION IRPSTR(GRIREP,ICRIRP,NCOOR,ICOOR,IPRINT) 1679C ****************************************************** 1680C *** Function that takes a coordinate as input, and *** 1681C *** returns the starting point of the irrep that *** 1682C *** span that irrep in GRIREP. *** 1683C ****************************************************** 1684#include "implicit.h" 1685#include "priunit.h" 1686C 1687#include "fcsym.h" 1688 DIMENSION GRIREP(NGORDR,NGVERT) 1689 DIMENSION ICRIRP(NCOOR,2) 1690C 1691C *** Finding the starting point. *** 1692C 1693 IF (ICRIRP(ICOOR,1) .GT. N1DIME) THEN 1694 ISTRTT = N1DIME + 4*(ICRIRP(ICOOR,1)-N1DIME-1) 1695 IF (ICRIRP(ICOOR,2).EQ.1) THEN 1696 ISTRTT = ISTRTT + 2 1697 END IF 1698 ELSE 1699 ISTRTT = ICRIRP(ICOOR,1)-1 1700 END IF 1701C 1702C *** Assigning the value. *** 1703C 1704 IRPSTR = ISTRTT 1705C 1706 RETURN 1707 END 1708C 1709C 1710C /* Deck anasym */ 1711 SUBROUTINE ANLSYM(GRIREP,ICRIRP,INDSTP,IREPST,RTMPCO,NCOOR,IORDER, 1712 & IORDR,IPRINT) 1713C ******************************************************* 1714C *** Subroutine that decides how the symmetry is for *** 1715C *** analytical derivatives for this numerical *** 1716C *** distortion. Output parameter is RTMPCO, which is*** 1717C *** assigned a value according to the symmetry prop.*** 1718C *** of the analytical derivative. *** 1719C ******************************************************* 1720#include "implicit.h" 1721#include "priunit.h" 1722 PARAMETER (D1 = 1.0D0, D0 = 0.0D0) 1723#include "numder.h" 1724#include "fcsym.h" 1725 LOGICAL SBDELM, SMEIRP 1726 DIMENSION GRIREP(NGORDR,NGVERT) 1727 DIMENSION INDSTP(NTORDR), IREPST(NMORDR), ICRIRP(NCOOR,2) 1728C 1729 IF (NAORDR.EQ.0) THEN 1730 RTMPCO=D1 1731 ELSE IF (NAORDR.EQ.1) THEN 1732 RTMPCO = D0 1733 DO ICOOR = 1, NCOOR 1734 IPSTRT = IRPSTR(GRIREP,ICRIRP,NCOOR,ICOOR,IPRINT) 1735 SBDELM = SMEIRP(GRIREP,ICRIRP,INDSTP,IREPST,ICOOR,IPSTRT, 1736 & IORDR,NCOOR,IPRINT) 1737 IF (SBDELM) THEN 1738 NSTEP = 1 1739 IF (ICRIRP(ICOOR,1).GT.N1DIME) NSTEP = 2 1740 DO ISTP = 1, NSTEP 1741 RTMPCO = RTMPCO + GRIREP(IORDER,IPSTRT+ISTP) 1742 END DO 1743 END IF 1744 END DO 1745 END IF 1746C 1747 RETURN 1748 END 1749C 1750C 1751C /* Deck smeirp */ 1752 LOGICAL FUNCTION SMEIRP(GRIREP,ICRIRP,INDSTP,IREPST,ICOOR,IPSTRT, 1753 & IORDR,NCOOR,IPRINT) 1754C ********************************************************* 1755C *** Subroutine that checks if the coordinate ICOOR is *** 1756C *** totally symmetric in the distorted geometry. The *** 1757C *** distortions is done along the INDSTP coordinates. *** 1758C ********************************************************* 1759#include "implicit.h" 1760#include "priunit.h" 1761 PARAMETER (D1=1.0D0, DTHRS=1.0D-12) 1762#include "fcsym.h" 1763#include "numder.h" 1764 LOGICAL EXSTEL, TSMEIR 1765 DIMENSION GRIREP(NGORDR,NGVERT) 1766 DIMENSION INDSTP(NTORDR), ICRIRP(NCOOR,2), IREPST(NMORDR) 1767C 1768 TSMEIR = .TRUE. 1769 DO IEL = 1, NGORDR 1770C 1771C *** Checkin if this element is in the subgroup. *** 1772C 1773 EXSTEL = .TRUE. 1774 DO IC = 1, IORDR 1775 KSTP = ICRIRP(INDSTP(IC),2) + 1 1776 IF ((GRIREP(IEL,IREPST(IC)+KSTP)-D1)**2 .GT. DTHRS) 1777 & EXSTEL = .FALSE. 1778 END DO 1779C 1780C *** If this element exist, check if irep is totally *** 1781C *** symmetric. *** 1782C 1783 KSTP = ICRIRP(ICOOR,2) + 1 1784 IF (EXSTEL.AND.((GRIREP(IEL,IPSTRT+KSTP)-D1)**2.GT.DTHRS)) 1785 & TSMEIR = .FALSE. 1786C 1787 END DO 1788C 1789C *** Assigning the result to the function. *** 1790C 1791 SMEIRP = TSMEIR 1792C 1793 RETURN 1794 END 1795C 1796C 1797C /* Deck ffcdep */ 1798 SUBROUTINE FFCDEP(KDPMTX,INDSTP,IDXTMP,NLDPMX,LDPMTX,IFRSTD, 1799 & IRSRDR,IORDR,IEIIN,IPRINT,DEPFC) 1800#include "implicit.h" 1801Chjaaj DEBUG priunit.h 1802#include "priunit.h" 1803#include "mxcent.h" 1804 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, DTHR = 1.0D-12) 1805#include "fcsym.h" 1806#include "numder.h" 1807 LOGICAL DONE, DEPFC 1808 DIMENSION KDPMTX(LDPMTX,NSTRDR,IFRSTD), INDSTP(NTORDR), 1809 & IDXTMP(NMORDR) 1810C 1811 IDXTMP(1) = INDSTP(1) 1812C 1813 DO 100 II = 2, IORDR-IRSRDR 1814 IDXTMP(II) = INDSTP(1) 1815 100 CONTINUE 1816C 1817 ISTRT = IORDR-IRSRDR 1818 DO 200 II = 1, IRSRDR 1819 IDXTMP(ISTRT+II) = INDSTP(II+1) 1820 200 CONTINUE 1821C 1822 DO 300 II = ISTRT+IRSRDR+1, ISTRT+IRSRDR+IEIIN 1823 DONE = .FALSE. 1824 DO 400 IJ = ISTRT, II-1 1825 IF ((IDXTMP(IJ).LT.INDSTP(II)).AND.(.NOT.DONE)) THEN 1826 DO 500 IK = II, IJ+1, -1 1827 IDXTMP(IK) = IDXTMP(IK-1) 1828 500 CONTINUE 1829 IDXTMP(IJ) = INDSTP(II) 1830 DONE = .TRUE. 1831 END IF 1832 400 CONTINUE 1833 IF (.NOT.DONE) IDXTMP(II) = INDSTP(II) 1834 300 CONTINUE 1835C 1836 ITHDIM = ISTRT+IRSRDR+IEIIN 1837 if (ITHDIM+1 .gt. NSTRDR) then 1838 write(lupri,*) 'FFCDEP error: ITHDIM+1 .gt. NSTRDR' 1839 write(lupri,*) 1840 & 'FFCDEP error: ITHDIM,ISTRT,IRSRDR,IEIIN,NSTRDR', 1841 & ITHDIM,ISTRT,IRSRDR,IEIIN,NSTRDR 1842 call quit('error in FFCDEP') 1843 end if 1844 MINTST = MIN(ITHDIM,NMORDR) 1845 DO 600 II = 1, NLDPMX 1846 IF (((KDPMTX(II,ITHDIM+1,1).EQ.0).OR.(MINTST.EQ.NMORDR)).AND. 1847 & (.NOT.DEPFC)) THEN 1848 DEPFC = .TRUE. 1849 DO 700 IRDR = 1, ITHDIM 1850 DEPFC = DEPFC.AND.(KDPMTX(II,IRDR,1).EQ.IDXTMP(IRDR)) 1851 700 CONTINUE 1852 END IF 1853 600 CONTINUE 1854C 1855 RETURN 1856 END 1857C 1858C 1859C /*Deck ortscp*/ 1860 SUBROUTINE ORTSCP(TRIAL,SYMCOR,ICRIRP,ISYMCO,KIREP,KIPL,IPRINT, 1861 & EXIST) 1862#include "implicit.h" 1863#include "priunit.h" 1864#include "mxcent.h" 1865 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, DTHR = 1.0D-10) 1866#include "trkoor.h" 1867 LOGICAL EXIST 1868 DIMENSION TRIAL(NCOOR), SYMCOR(NCOOR,NCOOR), ICRIRP(NCOOR,2) 1869C 1870 DO ICOOR = 1, ISYMCO 1871 IREP = ICRIRP(ICOOR,1) 1872C 1873 IF ((IREP .EQ. KIREP) .AND. (ICRIRP(ICOOR,2).EQ.(KIPL-1))) THEN 1874C 1875 RLENGS = D0 1876 DO ICOOR1 = 1, NCOOR 1877 RLENGS = RLENGS + SYMCOR(ICOOR1,ICOOR)**2 1878 END DO 1879 RINVS = D1/RLENGS 1880C 1881 GCCNST = D0 1882 DO ICOOR1 = 1, NCOOR 1883 GCCNST = GCCNST + SYMCOR(ICOOR1,ICOOR)*TRIAL(ICOOR1) 1884 END DO 1885C 1886 IF (ABS(GCCNST) .GT. DTHR) THEN 1887C 1888 DO ICOOR1 = 1, NCOOR 1889 TRIAL(ICOOR1) = TRIAL(ICOOR1) 1890 & - RINVS*GCCNST*SYMCOR(ICOOR1,ICOOR) 1891 END DO 1892C 1893 END IF 1894C 1895C *** Normalizing. *** 1896C 1897 RLENGT = D0 1898 DO ICOOR1 = 1, NCOOR 1899 RLENGT = RLENGT+ TRIAL(ICOOR1)**2 1900 END DO 1901 RLENGT = SQRT(RLENGT) 1902 IF (RLENGT.GT.DTHR) THEN 1903 RINV = D1/RLENGT 1904 DO ICOOR1 = 1, NCOOR 1905 TRIAL(ICOOR1) = TRIAL(ICOOR1)*RINV 1906 END DO 1907 END IF 1908 END IF 1909 END DO 1910C 1911C *** Has it become the zero vector?*** 1912C 1913 RLENGT = D0 1914 DO ICOOR1 = 1, NCOOR 1915 RLENGT = RLENGT+ TRIAL(ICOOR1)**2 1916 END DO 1917 IF (RLENGT.LT.DTHR) EXIST = .TRUE. 1918C 1919C *** Test printing. *** 1920C 1921 IF ((IPRINT .GT. 20) .AND. .NOT.EXIST) THEN 1922 WRITE (LUPRI,'(A)') 'Orthoganlized symmetry coordinate:' 1923C 1924 WRITE (LUPRI,'(A,2I4)') 'Symmetry', KIREP, KIPL-1 1925 DO I = 1, NCOOR 1926 WRITE (LUPRI,'(10X,F8.4)') TRIAL(I) 1927 END DO 1928 END IF 1929C 1930 RETURN 1931 END 1932C 1933C 1934C /* Deck srtscr */ 1935 SUBROUTINE SRTSCR(SYMCOR,TMPCOR,ICRIRP,ICTMP,IPRINT) 1936C ****************************************************************** 1937C ***** This subroutine sorts the symmetry adapted coordinates ***** 1938C ***** such that all the coordinates that transforms as 1. ***** 1939C ***** of a 2 dimensional irrep (for a given irrep), lies ***** 1940C ***** before the coordinates that transforms as 2. row of ***** 1941C ***** that irrep. ***** 1942C ****************************************************************** 1943C 1944#include "implicit.h" 1945#include "priunit.h" 1946#include "mxcent.h" 1947C 1948#include "fcsym.h" 1949#include "trkoor.h" 1950 DIMENSION SYMCOR(NCOOR,NCOOR), TMPCOR(NCOOR,NCOOR) 1951 DIMENSION ICRIRP(NCOOR,2), ICTMP(NCOOR,2) 1952 LOGICAL FND2D 1953C 1954C *** Copying all the coordinates belonging to 1. dim ireps *** 1955C *** over to TMPCOR *** 1956C 1957 IC = 0 1958 DO 200 IC2 = 1, NCOOR 1959 IF (ICRIRP(IC2,1) .LE. N1DIME) THEN 1960 IC = IC + 1 1961 ICTMP(IC,1) = ICRIRP(IC2,1) 1962 ICTMP(IC,2) = 0 1963 DO 300 IC1 = 1, NCOOR 1964 TMPCOR(IC1,IC2) = SYMCOR(IC1,IC2) 1965 300 CONTINUE 1966 END IF 1967 200 CONTINUE 1968 N1DCOR = IC 1969C 1970C *** Sorting the rest of the coordinates according *** 1971C *** to the principle above. *** 1972C 1973 IC = N1DCOR 1974 ITMPC = N1DCOR 1975 DO 400 IIREP = N1DIME+1, NIREP 1976C 1977 DO ITIM = 0, 1 1978 DO ICOOR1 = N1DCOR, NCOOR 1979 IF ((ICRIRP(ICOOR1,1).EQ.IIREP).AND. 1980 & (ICRIRP(ICOOR1,2).EQ.ITIM )) THEN 1981 IC = IC + 1 1982 ICTMP(IC,1) = IIREP 1983 ICTMP(IC,2) = ITIM 1984 DO ICOOR2 = 1, NCOOR 1985 TMPCOR(ICOOR2,IC) = SYMCOR(ICOOR2,ICOOR1) 1986 END DO 1987 END IF 1988 END DO 1989 END DO 1990 400 CONTINUE 1991C 1992C *** Test print *** 1993C 1994 IF (IPRINT .GT. 20) THEN 1995 CALL HEADER('Symmetry coordinates before sorting',-1) 1996 DO 700 IC1 = 1, NCOOR 1997 WRITE (LUPRI,'(12F12.8)') (SYMCOR(IC1,IC2),IC2=1,NCOOR) 1998 700 CONTINUE 1999 END IF 2000C 2001C *** Assigning the new coordinate order to SYMCOR *** 2002C 2003 DO IC2 = 1, NCOOR 2004 ICRIRP(IC2,1) = ICTMP(IC2,1) 2005 ICRIRP(IC2,2) = ICTMP(IC2,2) 2006 DO IC1 = 1, NCOOR 2007 SYMCOR(IC1,IC2) = TMPCOR(IC1,IC2) 2008 END DO 2009 END DO 2010C 2011C *** Test print *** 2012C 2013 IF (IPRINT .GT. 20) THEN 2014 CALL HEADER('The rows are placed as follows',-1) 2015 WRITE (LUPRI,'(15I6)') (ICRIRP(IC,1), IC = 1, NCOOR) 2016 WRITE (LUPRI,'(15I6)') (ICRIRP(IC,2), IC = 1, NCOOR) 2017 CALL HEADER('Symmetry coordinates after sorting',-1) 2018 DO 900 IC1 = 1, NCOOR 2019 WRITE (LUPRI,'(12F15.8)') (SYMCOR(IC1,IC2),IC2=1,NCOOR) 2020 900 CONTINUE 2021 END IF 2022C 2023 RETURN 2024 END 2025C 2026C 2027C /*Deck fsdcst*/ 2028 SUBROUTINE FSDCST(SYMCOR,GRIREP,DCOEFF,WORK,KDPMTX,NMIDPC,ICRIRP, 2029 & LDPMTX,IFRSTD,NLDPMX,LWORK,IPRINT) 2030#include "implicit.h" 2031#include "priunit.h" 2032#include "mxcent.h" 2033C 2034#include "trkoor.h" 2035#include "numder.h" 2036#include "fcsym.h" 2037 DIMENSION SYMCOR(NCOOR,NCOOR), GRIREP(NGORDR,NGVERT), 2038 & DCOEFF(LDPMTX,IFRSTD), WORK(LWORK) 2039 DIMENSION ICRIRP(NCOOR,2), KDPMTX(LDPMTX,NSTRDR,IFRSTD), 2040 & NMIDPC(LDPMTX) 2041C 2042 call flshfo(5) 2043 DO 100 IORDR = 1, NMORDR 2044 IF (IORDR .EQ. 2) THEN 2045 CALL SCNFCS(DCOEFF,KDPMTX,ICRIRP,LDPMTX,IFRSTD,NLDPMX, 2046 & IPRINT) 2047C 2048 DO 200 II = 1, NLDPMX 2049 NMIDPC(II) = 1 2050 200 CONTINUE 2051 END IF 2052C 2053 IF (IORDR .EQ. 3) THEN 2054 KEQUMT = 1 2055 KIDXTS = KEQUMT + 2**(2*IORDR) 2056 KLAST = KIDXTS + 8*NMORDR 2057 LWRK = LWORK - KLAST 2058 CALL THRDFC(DCOEFF,GRIREP,WORK(KEQUMT),WORK(KLAST),KDPMTX, 2059 & NMIDPC,ICRIRP,WORK(KIDXTS),LDPMTX,IFRSTD,NLDPMX, 2060 & LWRK,IPRINT) 2061 END IF 2062C 2063 IF (IORDR .EQ. 4) THEN 2064 KEQUMT = 1 2065 KIDXTS = KEQUMT + 2**(2*IORDR) 2066 KLAST = KIDXTS + 16*NMORDR 2067 LWRK = LWORK - KLAST 2068 CALL FRTHFC(DCOEFF,GRIREP,WORK(KEQUMT),WORK(KLAST),KDPMTX, 2069 & NMIDPC,ICRIRP,WORK(KIDXTS),LDPMTX,IFRSTD,NLDPMX, 2070 & LWRK,IPRINT) 2071 END IF 2072 100 CONTINUE 2073C 2074 RETURN 2075 END 2076C 2077C 2078C /* Deck thrdfc*/ 2079 SUBROUTINE THRDFC(DCOEFF,GRIREP,EQUMTX,WORK,KDPMTX,NMIDPC,ICRIRP, 2080 & IDXTST,LDPMTX,IFRSTD,NLDPMX,LWORK,IPRINT) 2081#include "implicit.h" 2082#include "priunit.h" 2083#include "mxcent.h" 2084 PARAMETER (D0=0.0D0, D1=1.0D0, DMTHR = 1.0D-6) 2085#include "trkoor.h" 2086#include "numder.h" 2087#include "fcsym.h" 2088 LOGICAL DPNDCY, EXIST, IEXIST, TRIVIA 2089 DIMENSION GRIREP(NGORDR,NGVERT), DCOEFF(LDPMTX,IFRSTD), 2090 & TCOEFF(8,8), EQUMTX(8,8),WORK(LWORK) 2091 DIMENSION KDPMTX(LDPMTX,NSTRDR,IFRSTD), NMIDPC(LDPMTX), 2092 & ICRIRP(NCOOR,2),IDEP(8,8), IDXTST(NMORDR,8), 2093 & ICRWMX(8,3) 2094C 2095 NLDPST = NLDPMX 2096C 2097 DO 100 IC3 = 1, NDCOOR 2098 DO 100 IC2 = 1, IC3 2099 DO 100 IC1 = 1, IC2 2100C 2101 IF (ICRIRP(IC1,1).GT.N1DIME) THEN 2102 ISTII = N1DIME + 4*(ICRIRP(IC1,1)-N1DIME-1) 2103 IF (ICRIRP(IC1,2).EQ.0) THEN 2104 N1STR = 1 2105 N1END = 2 2106 N1STP = 1 2107 ELSE 2108 N1STR = 2 2109 N1END = 1 2110 N1STP = -1 2111 END IF 2112 ELSE 2113 ISTII = ICRIRP(IC1,1)-1 2114 N1STR = 1 2115 N1END = 1 2116 N1STP = 1 2117 END IF 2118C 2119 IF (ICRIRP(IC2,1).GT.N1DIME) THEN 2120 ISTIJ = N1DIME + 4*(ICRIRP(IC2,1)-N1DIME-1) 2121 IF (ICRIRP(IC2,2).EQ.0) THEN 2122 N2STR = 1 2123 N2END = 2 2124 N2STP = 1 2125 ELSE 2126 N2STR = 2 2127 N2END = 1 2128 N2STP = -1 2129 END IF 2130 ELSE 2131 ISTIJ = ICRIRP(IC2,1)-1 2132 N2STR = 1 2133 N2END = 1 2134 N2STP = 1 2135 END IF 2136C 2137 IF (ICRIRP(IC3,1).GT.N1DIME) THEN 2138 ISTIK = N1DIME + 4*(ICRIRP(IC3,1)-N1DIME-1) 2139 IF (ICRIRP(IC3,2).EQ.0) THEN 2140 N3STR = 1 2141 N3END = 2 2142 N3STP = 1 2143 ELSE 2144 N3STR = 2 2145 N3END = 1 2146 N3STP = -1 2147 END IF 2148 ELSE 2149 ISTIK = ICRIRP(IC3,1)-1 2150 N3STR = 1 2151 N3END = 1 2152 N3STP = 1 2153 END IF 2154C 2155C *** Number of coordinates that transform as the *** 2156C *** same row of the same irep (which is equal to *** 2157C *** the spacing between function 1 and 2 of the *** 2158C *** same mode in a 2-dimensional irep) *** 2159C 2160 N1ICOR = 0 2161 N2ICOR = 0 2162 N3ICOR = 0 2163 DO 200 IC = 1, NDCOOR 2164 IF ((ICRIRP(IC,1).EQ.ICRIRP(IC1,1)).AND. 2165 & (ICRIRP(IC,2).EQ. 0)) THEN 2166 N1ICOR = N1ICOR + 1 2167 END IF 2168C 2169 IF ((ICRIRP(IC,1).EQ.ICRIRP(IC2,1)).AND. 2170 & (ICRIRP(IC,2).EQ. 0)) THEN 2171 N2ICOR = N2ICOR + 1 2172 END IF 2173C 2174 IF ((ICRIRP(IC,1).EQ.ICRIRP(IC3,1)).AND. 2175 & (ICRIRP(IC,2).EQ. 0)) THEN 2176 N3ICOR = N3ICOR + 1 2177 END IF 2178 200 CONTINUE 2179C 2180 KDIM = 64 2181 CALL DZERO(EQUMTX,KDIM) 2182C 2183 ICIJK = 0 2184 DO 300 IK = N3STR, N3END, N3STP 2185 DO 300 IJ = N2STR, N2END, N2STP 2186 DO 300 II = N1STR, N1END, N1STP 2187 ICIJK = ICIJK + 1 2188C 2189 ICRWMX(ICIJK,1) = IC3 - ICRIRP(IC3,2)*N3ICOR + (IK-1)*N3ICOR 2190 ICRWMX(ICIJK,2) = IC2 - ICRIRP(IC2,2)*N2ICOR + (IJ-1)*N2ICOR 2191 ICRWMX(ICIJK,3) = IC1 - ICRIRP(IC1,2)*N1ICOR + (II-1)*N1ICOR 2192C 2193 ICTUV = 0 2194 DO 400 IV = N3STR, N3END, N3STP 2195 DO 400 IU = N2STR, N2END, N2STP 2196 DO 400 IT = N1STR, N1END, N1STP 2197 ICTUV = ICTUV + 1 2198C 2199 IV1 = ISTII + 2*(II-1) + IT 2200 IV2 = ISTIJ + 2*(IJ-1) + IU 2201 IV3 = ISTIK + 2*(IK-1) + IV 2202 DO 500 IGORDR = 1, NGORDR 2203 EQUMTX(ICTUV,ICIJK) = EQUMTX(ICTUV,ICIJK) 2204 & +(GRIREP(IGORDR,IV1) 2205 & *GRIREP(IGORDR,IV2) 2206 & *GRIREP(IGORDR,IV3))/DBLE(NGORDR) 2207 500 CONTINUE 2208 400 CONTINUE 2209 300 CONTINUE 2210C 2211 DPNDCY = .FALSE. 2212 IF (ICIJK.GT.1) THEN 2213 SUMMTX = D0 2214 DO 600 IIIJK = 1, ICIJK 2215 DO 600 IITUV = 1, ICTUV 2216 SUMMTX = SUMMTX + ABS(EQUMTX(IITUV,IIIJK)) 2217 600 CONTINUE 2218 IF (SUMMTX .GT. DMTHR) DPNDCY = .TRUE. 2219 END IF 2220C 2221 IF (DPNDCY) THEN 2222 DO 700 II = 1, ICIJK 2223 EQUMTX(II,II) = EQUMTX(II,II)-D1 2224 700 CONTINUE 2225C 2226 IF (IPRINT .GT. 20) THEN 2227 WRITE (LUPRI,'(A)') 'Matrix to determine dependent' // 2228 & ' force constants' 2229 WRITE (LUPRI,'(2X,A,I4)') 'Number of components', 2230 & (N1STR+N1END-1)*(N2STR+N2END-1)*(N3STR+N3END-1) 2231 WRITE (LUPRI,'(2X,A,3I4)') 'Component', IC1, IC2, IC3 2232 DO IIIJK = 1, ICIJK 2233 WRITE (LUPRI,'(8F8.4)') 2234 & (EQUMTX(IITUV,IIIJK),IITUV=1,ICTUV) 2235 END DO 2236 WRITE (LUPRI,'(A)') ' ' 2237 WRITE (LUPRI,'(A)') ' ' 2238 END IF 2239C 2240C *** Sorting the components, such that the force-constants *** 2241C *** calculated is situated last. *** 2242C 2243 KITMPE = 1 2244 KITMPR = KITMPE + 8 2245 CALL SRTNCF(EQUMTX,WORK(KITMPR),ICRWMX,WORK(KITMPE),8,3, 2246 & ICIJK,IPRINT) 2247 ICTUV = ICIJK 2248C 2249C *** Solving the homogeneous linear equation system. *** 2250C 2251 KIRNDX = 1 2252 KZINDX = KIRNDX + 8 2253 KIDXTP = KZINDX + 8 2254 CALL DIAGUD(EQUMTX,TCOEFF,IDEP,WORK(KIRNDX),WORK(KZINDX), 2255 & WORK(KIDXTP),8,ICIJK,NROW,NIDEP,IPRINT) 2256C 2257C *** Independent component. *** 2258C *** Checking if it is av valid component, and assigning the *** 2259C *** "coordinate-numbers" *** 2260C 2261 DO 800 IDP = 1, NIDEP 2262 DO 900 ILNGTH = 1, ICIJK 2263 IF (ILNGTH.EQ.IDEP(1,IDP+1)) THEN 2264 IDXTST(1,IDP) = ICRWMX(IDEP(1,IDP+1),1) 2265 IDXTST(2,IDP) = ICRWMX(IDEP(1,IDP+1),2) 2266 IDXTST(3,IDP) = ICRWMX(IDEP(1,IDP+1),3) 2267C 2268C *** The other reason for discarding the batch *** 2269C *** based on the independent component, is that *** 2270C *** the batch already EXIST *** 2271C 2272 CALL CHKCMP(KDPMTX,IDXTST(1,IDP),3,LDPMTX,IFRSTD, 2273 & NLDPST,NLDPMX,NIDEP,IPRINT,EXIST) 2274 END IF 2275 900 CONTINUE 2276 800 CONTINUE 2277C 2278C *** Finding the dependent components *** 2279C 2280 NSTART = NLDPMX 2281 IF (.NOT.EXIST) THEN 2282 DO 1100 IROW = 1, NROW-1 2283 DO 1200 ILNGTH = 1, ICIJK 2284 IF (ILNGTH.EQ.IDEP(IROW,1)) THEN 2285C 2286C *** Checking whether there are several equal solutions for *** 2287C *** one dependent force constant. *** 2288C 2289 IDP3 = ICRWMX(IDEP(IROW,1),1) 2290 IDP2 = ICRWMX(IDEP(IROW,1),2) 2291 IDP1 = ICRWMX(IDEP(IROW,1),3) 2292 IEXIST = .FALSE. 2293 DO 1300 ICOUNT = NSTART+1, NLDPMX 2294 IF (.NOT.IEXIST) THEN 2295 IEXIST = (IDP3.EQ.KDPMTX(ICOUNT,1,1)).AND. 2296 & (IDP2.EQ.KDPMTX(ICOUNT,2,1)).AND. 2297 & (IDP1.EQ.KDPMTX(ICOUNT,3,1)) 2298 END IF 2299 1300 CONTINUE 2300C 2301C *** Checking if this is a trivial solution of the dependent *** 2302C *** force constant of the type K_{aaa} = K_{aaa} *** 2303C 2304 IF (NIDEP.EQ.1) THEN 2305 TRIVIA = (IDP3.EQ.IDXTST(1,1)).AND. 2306 & (IDP2.EQ.IDXTST(2,1)).AND. 2307 & (IDP1.EQ.IDXTST(3,1)) 2308 END IF 2309C 2310 IF ((.NOT.IEXIST).AND.(.NOT.TRIVIA)) THEN 2311 NLDPMX = NLDPMX + 1 2312 KDPMTX(NLDPMX,1,1) = IDP3 2313 KDPMTX(NLDPMX,2,1) = IDP2 2314 KDPMTX(NLDPMX,3,1) = IDP1 2315 NMIDPC(NLDPMX ) = NIDEP 2316 DO 1400 IDP = 1, NIDEP 2317 KDPMTX(NLDPMX,1,IDP+1) = IDXTST(1,IDP) 2318 KDPMTX(NLDPMX,2,IDP+1) = IDXTST(2,IDP) 2319 KDPMTX(NLDPMX,3,IDP+1) = IDXTST(3,IDP) 2320 DCOEFF(NLDPMX ,IDP ) = TCOEFF(IROW,IDP) 2321 1400 CONTINUE 2322 ELSE 2323 NROW = NROW - 1 2324 END IF 2325 END IF 2326 1200 CONTINUE 2327 1100 CONTINUE 2328C 2329C 2330C 2331 KIDXTP = 1 2332 KNMIDP = KIDXTP + 3 2333 KMVTMP = KNMIDP + 3 2334 KICMPI = KMVTMP + 3 2335 CALL CHKLCP(DCOEFF,KDPMTX,ICRIRP,WORK(KIDXTP), 2336 & WORK(KNMIDP),WORK(KMVTMP),WORK(KICMPI), 2337 & LDPMTX,IFRSTD,3,NROW-1,NLDPMX-NROW+NIDEP+1, 2338 & NIDEP,IPRINT) 2339C 2340C *** Test print *** 2341C 2342 CALL HEADER('Coordinate dependency in cartesians.',-1) 2343 WRITE (LUPRI,'(5X,A,I4)') 'Number of independent' // 2344 & 'force constants.', NMIDPC(NLDPMX) 2345 WRITE (LUPRI,'(5X,A)') 'Dependent components ' // 2346 & 'Independent components Coeffisient' 2347 DO II = NLDPMX-NROW+2, NLDPMX 2348 WRITE (LUPRI,'(A,3I4,A,3I4,A,F8.4)') ' ', 2349 & (KDPMTX(II,IX,1),IX=1,3), ' ', 2350 & (KDPMTX(II,IX,2),IX=1,3), ' ', 2351 & DCOEFF(II,1) 2352 DO IDP = 2, NIDEP 2353 WRITE (LUPRI,'(A,3I4,A,F8.4)') 2354 & ' ', 2355 & (KDPMTX(II,IX,IDP+1),IX=1,3), ' ', 2356 & DCOEFF(II,IDP) 2357 END DO 2358 END DO 2359 WRITE (LUPRI,'(A)')' ' 2360 WRITE (LUPRI,'(A)')' ' 2361 ELSE 2362 IF (IPRINT.GT.20) THEN 2363 WRITE (LUPRI,'(A)') ' ' 2364 WRITE (LUPRI,'(A)') ' ' 2365 WRITE (LUPRI,'(A)') 2366 & 'Component not going through due to:' 2367 WRITE (LUPRI,'(A,L1)')'Exist : ', EXIST 2368 WRITE (LUPRI,'(A)') ' ' 2369 WRITE (LUPRI,'(A)') ' ' 2370 END IF 2371 END IF 2372 END IF 2373C 2374 100 CONTINUE 2375C 2376 RETURN 2377 END 2378C 2379C 2380C /* Deck scnfcs */ 2381 SUBROUTINE SCNFCS(DCOEFF,KDPMTX,ICRIRP,LDPMTX,IFRSTD,NLDPMX, 2382 & IPRINT) 2383C *************************************************** 2384C **** Figures out which second derivatives are **** 2385C **** symmetrical dependent. The dependent **** 2386C **** components are put in first row of KDPMTX **** 2387C **** The independent are put in the second row **** 2388C **** The dependency coeffecients are put in **** 2389C **** DCOEFF. **** 2390C *************************************************** 2391#include "implicit.h" 2392#include "priunit.h" 2393#include "mxcent.h" 2394 PARAMETER (D1 = 1.0D0) 2395#include "trkoor.h" 2396#include "numder.h" 2397#include "fcsym.h" 2398 DIMENSION DCOEFF(LDPMTX,IFRSTD), KDPMTX(LDPMTX,NSTRDR,IFRSTD), 2399 & ICRIRP(NCOOR,2) 2400C 2401C *** Loop over all the 2 dimensional ireps *** 2402C 2403 KDPSTR = NLDPMX+1 2404 DO 100 IIRP = N1DIME+1, NCVERT 2405C 2406C *** Find how many functions there are of this irrep *** 2407C 2408 INUM = 0 2409 DO 200 IC = 1, NDCOOR 2410 IF ((ICRIRP(IC,1).EQ.IIRP).AND.(ICRIRP(IC,2).EQ.0)) THEN 2411 IF (INUM.EQ.0) ISTART = IC 2412 INUM = INUM + 1 2413 END IF 2414 200 CONTINUE 2415C 2416C *** Assign dependencies *** 2417C 2418 DO 300 IC2 = ISTART, ISTART+INUM-1 2419 DO 300 IC1 = ISTART, IC2 2420 NLDPMX = NLDPMX + 1 2421 KDPMTX(NLDPMX,1,1) = IC2 + INUM 2422 KDPMTX(NLDPMX,2,1) = IC1 + INUM 2423 KDPMTX(NLDPMX,1,2) = IC2 2424 KDPMTX(NLDPMX,2,2) = IC1 2425 DCOEFF(NLDPMX,1 ) = D1 2426 300 CONTINUE 2427C 2428 100 CONTINUE 2429C 2430C *** Test print *** 2431C 2432 IF (NLDPMX .GE. KDPSTR) THEN 2433 CALL HEADER('Dependent and indepenent force coefficients',-1) 2434 WRITE (LUPRI,'(5X,A)') 'Dependent components Independent ' // 2435 & 'components coefficients' 2436 DO II = KDPSTR, NLDPMX 2437 WRITE (LUPRI,'(I15,I4,I19,I4,15X,F6.2)') 2438 & (KDPMTX(II,J,1),J=1,2), (KDPMTX(II,J,2),J=1,2), 2439 & DCOEFF(II,1) 2440 END DO 2441 END IF 2442C 2443 RETURN 2444 END 2445C 2446C 2447C /* Deck frthfc*/ 2448 SUBROUTINE FRTHFC(DCOEFF,GRIREP,EQUMTX,WORK,KDPMTX,NMIDPC,ICRIRP, 2449 & IDXTST,LDPMTX,IFRSTD,NLDPMX,LWORK,IPRINT) 2450#include "implicit.h" 2451#include "priunit.h" 2452#include "mxcent.h" 2453 PARAMETER (D0=0.0D0, D1=1.0D0, DMTHR = 1.0D-6) 2454 PARAMETER (ITHDDM = 16) 2455#include "trkoor.h" 2456#include "numder.h" 2457#include "fcsym.h" 2458 LOGICAL DPNDCY, EXIST, IEXIST, TRIVIA 2459 DIMENSION GRIREP(NGORDR,NGVERT), DCOEFF(LDPMTX,IFRSTD), 2460 & TCOEFF(ITHDDM,ITHDDM), EQUMTX(ITHDDM,ITHDDM),WORK(10000) 2461 DIMENSION KDPMTX(LDPMTX,NSTRDR,IFRSTD), NMIDPC(LDPMTX), 2462 & ICRIRP(NCOOR,2),IDEP(ITHDDM,ITHDDM), 2463 & IDXTST(NMORDR,ITHDDM), ICRWMX(ITHDDM,4), ITLRGR(ITHDDM) 2464C 2465c iprint = 55 2466 NLDPST = NLDPMX 2467 DO 100 IC4 = 1, NDCOOR 2468 DO 100 IC3 = 1, IC4 2469 DO 100 IC2 = 1, IC3 2470 DO 100 IC1 = 1, IC2 2471C 2472 IF (ICRIRP(IC1,1).GT.N1DIME) THEN 2473 ISTII = N1DIME + 4*(ICRIRP(IC1,1)-N1DIME-1) 2474 IF (ICRIRP(IC1,2).EQ.0) THEN 2475 N1STR = 1 2476 N1END = 2 2477 N1STP = 1 2478 ELSE 2479 N1STR = 2 2480 N1END = 1 2481 N1STP = -1 2482 END IF 2483 ELSE 2484 ISTII = ICRIRP(IC1,1)-1 2485 N1STR = 1 2486 N1END = 1 2487 N1STP = 1 2488 END IF 2489C 2490 IF (ICRIRP(IC2,1).GT.N1DIME) THEN 2491 ISTIJ = N1DIME + 4*(ICRIRP(IC2,1)-N1DIME-1) 2492 IF (ICRIRP(IC2,2).EQ.0) THEN 2493 N2STR = 1 2494 N2END = 2 2495 N2STP = 1 2496 ELSE 2497 N2STR = 2 2498 N2END = 1 2499 N2STP = -1 2500 END IF 2501 ELSE 2502 ISTIJ = ICRIRP(IC2,1)-1 2503 N2STR = 1 2504 N2END = 1 2505 N2STP = 1 2506 END IF 2507C 2508 IF (ICRIRP(IC3,1).GT.N1DIME) THEN 2509 ISTIK = N1DIME + 4*(ICRIRP(IC3,1)-N1DIME-1) 2510 IF (ICRIRP(IC3,2).EQ.0) THEN 2511 N3STR = 1 2512 N3END = 2 2513 N3STP = 1 2514 ELSE 2515 N3STR = 2 2516 N3END = 1 2517 N3STP = -1 2518 END IF 2519 ELSE 2520 ISTIK = ICRIRP(IC3,1)-1 2521 N3STR = 1 2522 N3END = 1 2523 N3STP = 1 2524 END IF 2525C 2526 IF (ICRIRP(IC4,1).GT.N1DIME) THEN 2527 ISTIL = N1DIME + 4*(ICRIRP(IC4,1)-N1DIME-1) 2528 IF (ICRIRP(IC4,2).EQ.0) THEN 2529 N4STR = 1 2530 N4END = 2 2531 N4STP = 1 2532 ELSE 2533 N4STR = 2 2534 N4END = 1 2535 N4STP = -1 2536 END IF 2537 ELSE 2538 ISTIL = ICRIRP(IC4,1)-1 2539 N4STR = 1 2540 N4END = 1 2541 N4STP = 1 2542 END IF 2543C 2544C *** Number of coordinates that transform as the *** 2545C *** same row of the same irep (which is equal to *** 2546C *** the spacing between function 1 and 2 of the *** 2547C *** same mode in a 2-dimensional irep) *** 2548C 2549 N1ICOR = 0 2550 N2ICOR = 0 2551 N3ICOR = 0 2552 N4ICOR = 0 2553 DO 200 IC = 1, NDCOOR 2554 IF ((ICRIRP(IC,1).EQ.ICRIRP(IC1,1)).AND. 2555 & (ICRIRP(IC,2).EQ. 0)) THEN 2556 N1ICOR = N1ICOR + 1 2557 END IF 2558C 2559 IF ((ICRIRP(IC,1).EQ.ICRIRP(IC2,1)).AND. 2560 & (ICRIRP(IC,2).EQ. 0)) THEN 2561 N2ICOR = N2ICOR + 1 2562 END IF 2563C 2564 IF ((ICRIRP(IC,1).EQ.ICRIRP(IC3,1)).AND. 2565 & (ICRIRP(IC,2).EQ. 0)) THEN 2566 N3ICOR = N3ICOR + 1 2567 END IF 2568C 2569 IF ((ICRIRP(IC,1).EQ.ICRIRP(IC4,1)).AND. 2570 & (ICRIRP(IC,2).EQ. 0)) THEN 2571 N4ICOR = N4ICOR + 1 2572 END IF 2573 200 CONTINUE 2574C 2575 KDIM = ITHDDM**2 2576 CALL DZERO(EQUMTX,KDIM) 2577C 2578 ICIJKL = 0 2579 DO 300 IL = N4STR, N4END, N4STP 2580 DO 300 IK = N3STR, N3END, N3STP 2581 DO 300 IJ = N2STR, N2END, N2STP 2582 DO 300 II = N1STR, N1END, N1STP 2583 ICIJKL = ICIJKL + 1 2584C 2585 ICRWMX(ICIJKL,1)= IC4 - ICRIRP(IC4,2)*N4ICOR + (IL-1)*N4ICOR 2586 ICRWMX(ICIJKL,2)= IC3 - ICRIRP(IC3,2)*N3ICOR + (IK-1)*N3ICOR 2587 ICRWMX(ICIJKL,3)= IC2 - ICRIRP(IC2,2)*N2ICOR + (IJ-1)*N2ICOR 2588 ICRWMX(ICIJKL,4)= IC1 - ICRIRP(IC1,2)*N1ICOR + (II-1)*N1ICOR 2589 2590C 2591 ICTUVX = 0 2592 DO 400 IX = N4STR, N4END, N4STP 2593 DO 400 IV = N3STR, N3END, N3STP 2594 DO 400 IU = N2STR, N2END, N2STP 2595 DO 400 IT = N1STR, N1END, N1STP 2596 ICTUVX = ICTUVX + 1 2597C 2598 IV1 = ISTII + 2*(II-1) + IT 2599 IV2 = ISTIJ + 2*(IJ-1) + IU 2600 IV3 = ISTIK + 2*(IK-1) + IV 2601 IV4 = ISTIL + 2*(IL-1) + IX 2602 DO 500 IGORDR = 1, NGORDR 2603 EQUMTX(ICTUVX,ICIJKL) = EQUMTX(ICTUVX,ICIJKL) 2604 & +(GRIREP(IGORDR,IV1) 2605 & *GRIREP(IGORDR,IV2) 2606 & *GRIREP(IGORDR,IV3) 2607 & *GRIREP(IGORDR,IV4))/DBLE(NGORDR) 2608 500 CONTINUE 2609 400 CONTINUE 2610 300 CONTINUE 2611C 2612 DPNDCY = .FALSE. 2613 IF (ICIJKL.GT.1) THEN 2614 SUMMTX = D0 2615 DO 600 IIIJK = 1, ICIJKL 2616 DO 600 IITUV = 1, ICTUVX 2617 SUMMTX = SUMMTX + ABS(EQUMTX(IITUV,IIIJK)) 2618 600 CONTINUE 2619 IF (SUMMTX .GT. DMTHR) DPNDCY = .TRUE. 2620 END IF 2621C 2622 IF (DPNDCY) THEN 2623 DO 700 II = 1, ICIJKL 2624 EQUMTX(II,II) = EQUMTX(II,II)-D1 2625 700 CONTINUE 2626C 2627 IF (IPRINT .GT. 20) THEN 2628 WRITE (LUPRI,'(A)') 'Matrix to determine dependent' // 2629 & ' force constants' 2630 WRITE (LUPRI,'(2X,A,I4)') 'Number of components', 2631 & (N1STR+N1END-1)*(N2STR+N2END-1) 2632 & *(N3STR+N3END-1)*(N4STR+N4END-1) 2633 WRITE (LUPRI,'(2X,A,4I4)') 'Component', IC1,IC2,IC3,IC4 2634 DO IIIJKL = 1, ICIJKL 2635 WRITE (LUPRI,'(16F8.4)') 2636 & (EQUMTX(IITUVX,IIIJKL),IITUVX=1,ICTUVX) 2637 END DO 2638 WRITE (LUPRI,'(A)') ' ' 2639 WRITE (LUPRI,'(A)') ' ' 2640 END IF 2641C 2642C *** Sorting the components, such thet the highest number *** 2643C *** component is first. In addition matrix elements for *** 2644C *** equal force-constants are added together *** 2645C 2646 KITMPE = 1 2647 KITMPR = KITMPE + ITHDDM 2648 CALL SRTNCF(EQUMTX,WORK(KITMPR),ICRWMX,WORK(KITMPE),ITHDDM, 2649 & 4,ICIJKL,IPRINT) 2650C 2651C *** Solving the homogeneous linear equation system. *** 2652C 2653 KIRNDX = 1 2654 KZINDX = KIRNDX + ICIJKL 2655 KIDXTP = KZINDX + ICIJKL 2656 CALL DIAGUD(EQUMTX,TCOEFF,IDEP,WORK(KIRNDX),WORK(KZINDX), 2657 & WORK(KIDXTP),ITHDDM,ICIJKL,NROW,NIDEP,IPRINT) 2658C 2659C *** Independent component. *** 2660C 2661 DO 800 IDP = 1, NIDEP 2662 DO 900 ILNGTH = 1, ICIJKL 2663 IF (ILNGTH.EQ.IDEP(1,IDP+1)) THEN 2664 IDXTST(1,IDP) = ICRWMX(IDEP(1,IDP+1),1) 2665 IDXTST(2,IDP) = ICRWMX(IDEP(1,IDP+1),2) 2666 IDXTST(3,IDP) = ICRWMX(IDEP(1,IDP+1),3) 2667 IDXTST(4,IDP) = ICRWMX(IDEP(1,IDP+1),4) 2668C 2669C *** One reason of reason for discarding the batch *** 2670C *** based on the independent component, is that *** 2671C *** the batch already EXIST *** 2672C 2673 CALL CHKCMP(KDPMTX,IDXTST(1,IDP),4,LDPMTX,IFRSTD, 2674 & NLDPST,NLDPMX,NIDEP,IPRINT,EXIST) 2675 END IF 2676 900 CONTINUE 2677 800 CONTINUE 2678C 2679C *** Finding the dependent components *** 2680C 2681 NSTART = NLDPMX 2682 IF (.NOT.EXIST) THEN 2683 DO 1100 IROW = 1, NROW-NIDEP 2684 DO 1200 ILNGTH = 1, ICIJKL 2685 IF (ILNGTH.EQ.IDEP(IROW,1)) THEN 2686C 2687C *** Checking whether there are several equal *** 2688C *** solutions for one dependent force constant. *** 2689C 2690C 2691 IDP4 = ICRWMX(IDEP(IROW,1),1) 2692 IDP3 = ICRWMX(IDEP(IROW,1),2) 2693 IDP2 = ICRWMX(IDEP(IROW,1),3) 2694 IDP1 = ICRWMX(IDEP(IROW,1),4) 2695 IEXIST = .FALSE. 2696 DO 1300 ICOUNT = NSTART+1, NLDPMX 2697 IF (.NOT.IEXIST) THEN 2698 IEXIST = (IDP4.EQ.KDPMTX(ICOUNT,1,1)).AND. 2699 & (IDP3.EQ.KDPMTX(ICOUNT,2,1)).AND. 2700 & (IDP2.EQ.KDPMTX(ICOUNT,3,1)).AND. 2701 & (IDP1.EQ.KDPMTX(ICOUNT,4,1)) 2702 END IF 2703 1300 CONTINUE 2704C 2705C *** Checking if this is a trivial solution of the dependent *** 2706C *** force constant of the type K_{abcd} = K_{abcd} *** 2707C 2708 IF (NIDEP.EQ.1) THEN 2709 TRIVIA = (IDP4.EQ.IDXTST(1,1)).AND. 2710 & (IDP3.EQ.IDXTST(2,1)).AND. 2711 & (IDP2.EQ.IDXTST(3,1)).AND. 2712 & (IDP1.EQ.IDXTST(4,1)) 2713 END IF 2714C 2715 IF ((.NOT.IEXIST).AND.(.NOT.TRIVIA)) THEN 2716 NLDPMX = NLDPMX + 1 2717 KDPMTX(NLDPMX,1,1) = IDP4 2718 KDPMTX(NLDPMX,2,1) = IDP3 2719 KDPMTX(NLDPMX,3,1) = IDP2 2720 KDPMTX(NLDPMX,4,1) = IDP1 2721 NMIDPC(NLDPMX ) = NIDEP 2722 write (lupri,*) 'hallo', NLDPMX 2723 DO 1400 IDP = 1, NIDEP 2724 KDPMTX(NLDPMX,1,IDP+1) = IDXTST(1,IDP) 2725 KDPMTX(NLDPMX,2,IDP+1) = IDXTST(2,IDP) 2726 KDPMTX(NLDPMX,3,IDP+1) = IDXTST(3,IDP) 2727 KDPMTX(NLDPMX,4,IDP+1) = IDXTST(4,IDP) 2728 DCOEFF(NLDPMX ,IDP ) = TCOEFF(IROW,IDP) 2729 1400 CONTINUE 2730 ELSE 2731 NROW = NROW - 1 2732 END IF 2733 END IF 2734 1200 CONTINUE 2735 1100 CONTINUE 2736C 2737C *** Checking whether the calulated independent components *** 2738C *** are the the components with least energy calculations. *** 2739C 2740 KIDXTP = 1 2741 KNMIDP = KIDXTP + 4 2742 KMVTMP = KNMIDP + 4 2743 KICMPI = KMVTMP + 4 2744 CALL CHKLCP(DCOEFF,KDPMTX,ICRIRP,WORK(KIDXTP), 2745 & WORK(KNMIDP),WORK(KMVTMP),WORK(KICMPI), 2746 & LDPMTX,IFRSTD,4,NROW-1,NLDPMX-NROW+NIDEP+1, 2747 & NIDEP,IPRINT) 2748C 2749C *** Test print *** 2750C 2751 IF (NROW.GT.1) THEN 2752 CALL HEADER('Coordinate dependency in cartesians.',-1) 2753 WRITE (LUPRI,'(5X,A,I4)') 'Number of independent' // 2754 & 'force constants.', NMIDPC(NLDPMX) 2755 WRITE (LUPRI,'(5X,A)') 'Dependent components ' // 2756 & 'Independent components Coeffisient' 2757 DO II = NLDPMX-NROW+NIDEP+1, NLDPMX 2758 write (lupri,*) 'component', ii 2759 WRITE (LUPRI,'(A,4I4,A,4I4,A,F8.4)') 2760 & ' ', (KDPMTX(II,IX,1),IX=1,4), ' ', 2761 & (KDPMTX(II,IX,2),IX=1,4), ' ', 2762 & DCOEFF(II,1) 2763 DO IDP = 2, NIDEP 2764 WRITE (LUPRI,'(A,4I4,A,F8.4)') 2765 & ' ', 2766 & (KDPMTX(II,IX,IDP+1),IX=1,4), ' ', 2767 & DCOEFF(II,IDP) 2768 END DO 2769 END DO 2770 WRITE (LUPRI,'(A)')' ' 2771 WRITE (LUPRI,'(A)')' ' 2772 ELSE 2773 WRITE (LUPRI,'(A)') ' ' 2774 WRITE (LUPRI,'(A)') ' ' 2775 WRITE (LUPRI,'(5X,A)') 'No valid component of the' // 2776 & ' independent force constant' 2777 WRITE (LUPRI,'(A)') ' ' 2778 WRITE (LUPRI,'(A)') ' ' 2779 END IF 2780 ELSE 2781 IF (IPRINT.GT.20) THEN 2782 WRITE (LUPRI,'(A)') ' ' 2783 WRITE (LUPRI,'(A)') ' ' 2784 WRITE (LUPRI,'(A)') 2785 & 'Component not going through due to:' 2786 WRITE (LUPRI,'(A,L1)')'Exist : ', EXIST 2787 WRITE (LUPRI,'(A)') ' ' 2788 WRITE (LUPRI,'(A)') ' ' 2789 END IF 2790 END IF 2791 END IF 2792C 2793 100 CONTINUE 2794 iprint = 0 2795C 2796 RETURN 2797 END 2798C 2799C 2800C /*Deck addpfc*/ 2801 SUBROUTINE ADDPFC(DERIV,DCOEFF,KDPMTX,NMIDPC,LDPMTX,IFRSTD,NDERIV, 2802 & NLDPMX,IPRINT) 2803C ******************************************************************* 2804C **** This assigns values to the dependent force-constants from **** 2805C **** the independent. **** 2806C ******************************************************************* 2807#include "implicit.h" 2808#include "priunit.h" 2809#include "mxcent.h" 2810 PARAMETER (D0 = 0.0D0) 2811#include "numder.h" 2812#include "fcsym.h" 2813 DIMENSION DERIV(NDERIV), DCOEFF(LDPMTX,IFRSTD) 2814 DIMENSION KDPMTX(LDPMTX,NSTRDR,IFRSTD), NMIDPC(LDPMTX) 2815 2816#include "trkoor.h" 2817 REAL*8 ERGMOL, GRDMOL(NCOOR), HESMOL(NCOOR,NCOOR) ! automatic arrays 2818C 2819 CALL ABAREAD_TAYMOL(ERGMOL,GRDMOL,HESMOL,NCOOR) 2820 2821 NSTART = 1 2822 DO 100 IRDR = 1, NMORDR 2823 IF (IRDR .EQ. 2) THEN 2824 ISTART = 0 2825 MINTST = MIN(NMORDR,3) 2826 DO 200 ILDP = 1, NLDPMX 2827 IF ((KDPMTX(ILDP,MINTST,1).EQ.0).OR. 2828 & (MINTST.EQ.2)) THEN 2829 HESMOL(KDPMTX(ILDP,1,1),KDPMTX(ILDP,2,1)) = 2830 & HESMOL(KDPMTX(ILDP,1,2),KDPMTX(ILDP,2,2)) 2831 ISTART = ISTART + 1 2832 END IF 2833 200 CONTINUE 2834 NSTART = NSTART + ISTART 2835 ELSE IF (IRDR.EQ.3) THEN 2836 ISTART = 0 2837 MINTST = MIN(NMORDR,4) 2838 DO 300 ILDP = NSTART, NLDPMX 2839 IF (((KDPMTX(ILDP,MINTST,1).EQ.0).AND.(MINTST.GT.3)) 2840 & .OR.((MINTST.EQ.3).AND.(KDPMTX(ILDP,3,1).NE.0))) THEN 2841C 2842 ISTART = ISTART + 1 2843C 2844 IDPDRV = ((KDPMTX(ILDP,1,1)-1)*KDPMTX(ILDP,1,1) 2845 & * (KDPMTX(ILDP,1,1)+1))/6 2846 & +((KDPMTX(ILDP,2,1)-1)*KDPMTX(ILDP,2,1))/2 2847 & + KDPMTX(ILDP,3,1) 2848 DERIV(IDPDRV) = D0 2849 DO 400 II = 1, NMIDPC(ILDP) 2850 IIDPDR = 2851 & ((KDPMTX(ILDP,1,II+1)-1)*KDPMTX(ILDP,1,II+1) 2852 & * (KDPMTX(ILDP,1,II+1)+1))/6 2853 & +((KDPMTX(ILDP,2,II+1)-1)*KDPMTX(ILDP,2,II+1))/2 2854 & + KDPMTX(ILDP,3,II+1) 2855 DERIV(IDPDRV) = DERIV(IDPDRV) 2856 & + DCOEFF(ILDP,II)*DERIV(IIDPDR) 2857 400 CONTINUE 2858 END IF 2859 300 CONTINUE 2860 NSTART = NSTART + ISTART 2861 ELSE IF (IRDR.EQ.4) THEN 2862 MINTST = MIN(NMORDR,5) 2863 IOFFST =(NDCOOR*(NDCOOR+1)*(NDCOOR+2))/6 2864C 2865 DO 500 ILDP = NSTART, NLDPMX 2866 IF (((KDPMTX(ILDP,MINTST,1).EQ.0).AND.(MINTST.GT.4)) 2867 & .OR.((MINTST.EQ.4).AND.(KDPMTX(ILDP,4,1).NE.0))) THEN 2868 IDPDRV = IOFFST + 2869 & ((KDPMTX(ILDP,1,1)-1)* KDPMTX(ILDP,1,1) 2870 & *(KDPMTX(ILDP,1,1)+1)*(KDPMTX(ILDP,1,1)+2))/24 2871 & +((KDPMTX(ILDP,2,1)-1)*KDPMTX(ILDP,2,1) 2872 & * (KDPMTX(ILDP,2,1)+1))/6 2873 & +((KDPMTX(ILDP,3,1)-1)*KDPMTX(ILDP,3,1))/2 2874 & + KDPMTX(ILDP,4,1) 2875 DERIV(IDPDRV) = D0 2876 DO 600 II = 1, NMIDPC(ILDP) 2877 IJ = II+1 2878 IIDPDR = IOFFST + 2879 & ((KDPMTX(ILDP,1,IJ)-1)* KDPMTX(ILDP,1,IJ) 2880 & *(KDPMTX(ILDP,1,IJ)+1)*(KDPMTX(ILDP,1,IJ)+2))/24 2881 & +((KDPMTX(ILDP,2,IJ)-1)* KDPMTX(ILDP,2,IJ) 2882 & *(KDPMTX(ILDP,2,IJ)+1))/6 2883 & +((KDPMTX(ILDP,3,IJ)-1)*KDPMTX(ILDP,3,IJ))/2 2884 & + KDPMTX(ILDP,4,IJ) 2885 DERIV(IDPDRV) = DERIV(IDPDRV) 2886 & + DCOEFF(ILDP,II)*DERIV(IIDPDR) 2887 600 CONTINUE 2888 END IF 2889 500 CONTINUE 2890 END IF 2891 100 CONTINUE 2892C 2893 CALL ABAWRIT_TAYMOL(ERGMOL,GRDMOL,HESMOL,NCOOR) 2894 RETURN 2895 END 2896C 2897C 2898C /* Deck diagud */ 2899 SUBROUTINE DIAGUD(RMTRIX,COEFF,IDEP,IRINDX,IZINDX,IDXTMP,NRDIM, 2900 & NDIM,NROW,NIDEP,IPRINT) 2901C ****************************************************************** 2902C *** Subroutine that solves a underdetermined set of homogenous *** 2903C *** linear equations. *** 2904C *** In : RMTRIX - coeffisient matrix to the equation set. *** 2905C *** Out: NIDEP -> Number of independent variables *** 2906C *** IDEP -> First row: Dependendent variables *** 2907C *** Scn. row: Dependend on this variable *** 2908C *** COEFF-> Dependence coeffisient. *** 2909C ****************************************************************** 2910#include "implicit.h" 2911#include "priunit.h" 2912 LOGICAL EXIT 2913 PARAMETER (D0 = 0.0D0, D1=1.0D0, DMTHR = 1.0D-6) 2914C 2915 DIMENSION RMTRIX(NRDIM,NRDIM), COEFF(NRDIM,NRDIM) 2916 DIMENSION IDEP(NRDIM,NRDIM), IRINDX(NRDIM), IZINDX(NRDIM), 2917 & IDXTMP(NRDIM,2) 2918C 2919 CALL RMZROR(RMTRIX,IRINDX,IZINDX,NRDIM,NDIM,NROW,IPRINT) 2920C 2921C *** The matrix is then "Gaussian-eliminated" to get it *** 2922C *** almost on diagonal form. *** 2923C 2924C *** Forward substitution. *** 2925C 2926 DO 100 ICOL = 1 , NROW 2927 IF (ABS(RMTRIX(ICOL,ICOL)) .GT. DMTHR) THEN 2928 BCOEFF = D1/RMTRIX(ICOL,ICOL) 2929C 2930 DO 200 IROW = ICOL+1, NROW 2931 TCOEFF = RMTRIX(IROW,ICOL)*BCOEFF 2932 DO 300 IICOL = ICOL, NROW 2933 RMTRIX(IROW,IICOL) = RMTRIX(IROW,IICOL) 2934 & - TCOEFF*RMTRIX(ICOL,IICOL) 2935 300 CONTINUE 2936 200 CONTINUE 2937C 2938 END IF 2939 100 CONTINUE 2940C 2941C *** Back substitution. *** 2942C 2943 DO 400 ICOL = NROW, 1, -1 2944 IF (ABS(RMTRIX(ICOL,ICOL)) .GT. DMTHR) THEN 2945 BCOEFF = D1/RMTRIX(ICOL,ICOL) 2946 DO 500 IROW = 1, ICOL-1 2947 TCOEFF = RMTRIX(IROW,ICOL)*BCOEFF 2948 DO 600 IICOL = ICOL, NROW 2949 RMTRIX(IROW,IICOL) = RMTRIX(IROW,IICOL) 2950 & - TCOEFF*RMTRIX(ICOL,IICOL) 2951 600 CONTINUE 2952 500 CONTINUE 2953 END IF 2954 400 CONTINUE 2955C 2956C *** Each row is then normalized ** 2957C 2958 DO 700 IROW = 1, NROW 2959 IF (ABS(RMTRIX(IROW,IROW)).GT.DMTHR) THEN 2960 BCOEFF = D1/RMTRIX(IROW,IROW) 2961 DO 800 ICOL = IROW, NROW 2962 RMTRIX(IROW,ICOL) = BCOEFF*RMTRIX(IROW,ICOL) 2963 800 CONTINUE 2964 END IF 2965 700 CONTINUE 2966C 2967C *** Finding rows in the matrix that are zero *** 2968C *** These are the independent variables. *** 2969C 2970 IZRO = 0 2971 IDP = 0 2972 EXIT = .FALSE. 2973 DO 900 IROW = 1, NROW 2974 SUMROW = D0 2975 DO 1100 ICOL = 1, NROW 2976 SUMROW = SUMROW + ABS(RMTRIX(IROW,ICOL)) 2977 1100 CONTINUE 2978 IF (SUMROW.LT.DMTHR) THEN 2979 IZRO = IZRO + 1 2980 IDXTMP(IZRO,2) = IROW 2981 ELSE 2982 IDP = IDP + 1 2983 IDXTMP(IDP,1) = IROW 2984 IDEP(IDP,1) = IRINDX(IROW) 2985 END IF 2986 900 CONTINUE 2987 NIDEP = IZRO 2988 NDDEP = IDP 2989C 2990C *** Assigning dependencies. *** 2991C 2992 DO 1200 IDP = 1, NDDEP 2993 DO 1200 IZRO = 1, NIDEP 2994 IDEP (IDP,IZRO+1) = IRINDX( IDXTMP(IZRO,2)) 2995 COEFF(IDP,IZRO ) =-RMTRIX(IDXTMP(IDP,1),IDXTMP(IZRO,2)) 2996 1200 CONTINUE 2997C 2998C *** Test print *** 2999C 3000 IF (IPRINT .GT. 20) THEN 3001 WRITE (LUPRI,'(A,I4)') 'Number of nonzero rows', NROW-NIDEP 3002 WRITE (LUPRI,'(A,16I4)') 'The non-zero rows are:', 3003 & (IRINDX(II),II=1,NROW) 3004 DO II = 1, NDDEP 3005 WRITE (LUPRI,'(16F8.4)') (RMTRIX(II,IJ),IJ=1,NROW) 3006 END DO 3007 WRITE (LUPRI,'(A)') ' ' 3008 WRITE (LUPRI,'(A)') ' ' 3009C 3010 WRITE (LUPRI,'(2X,A)') 'Dependent component ' // 3011 & 'independent component coeffisients' 3012 DO II = 1, NDDEP 3013 WRITE (LUPRI,'(10X,I4,20X,I3,15X,F8.4)') IDEP(II,1), 3014 & IDEP(II,2), COEFF(II,1) 3015 DO IDP = 2, NIDEP 3016 WRITE (LUPRI,'(34X,I3,15X,F8.4)') IDEP(II,IDP+1), 3017 & COEFF(II,IDP ) 3018 END DO 3019 END DO 3020 WRITE (LUPRI,'(A)') ' ' 3021 WRITE (LUPRI,'(A)') ' ' 3022 END IF 3023C 3024 RETURN 3025 END 3026C 3027C 3028C /*Deck rmzror*/ 3029 SUBROUTINE RMZROR(RMTRIX,IRINDX,IZINDX,NRDIM,NDIM,NROW,IPRINT) 3030C **************************************************************** 3031C *** Subroutine that removes the rows of a set of homogenous *** 3032C *** linear equations, that does not contribute -> type: *** 3033C *** a_i*x_i = 0 *** 3034C **************************************************************** 3035#include "implicit.h" 3036#include "priunit.h" 3037 PARAMETER (DMTHR = 1.0D-6) 3038 DIMENSION RMTRIX(NRDIM,NRDIM) 3039 DIMENSION IRINDX(NRDIM), IZINDX(NRDIM) 3040C 3041C *** Test print before removal of rows *** 3042C 3043 IF (IPRINT .GT. 20) THEN 3044 WRITE (LUPRI,'(A)') 'Matrix before removal of rows' 3045 DO II = 1, NDIM 3046 WRITE (LUPRI,'(16F8.4)') (RMTRIX(II,IJ),IJ=1,NDIM) 3047 END DO 3048 END IF 3049C 3050C *** Indexing for removal of rows *** 3051C *** IRINDX -> which rows are left after the removal. *** 3052C *** IZINDX -> which rows are supposed to be removed. *** 3053C 3054 IRROW = 0 3055 IZROW = 0 3056 MXDIM = NDIM 3057 DO 100 IDIM1 = 1, NDIM 3058 INOZRO = 0 3059 DO 200 IDIM2= 1, NDIM 3060 IF (ABS(RMTRIX(IDIM1,IDIM2)).GT.DMTHR) INOZRO = INOZRO + 1 3061 200 CONTINUE 3062C 3063 IF (INOZRO.GT.1) THEN 3064 IRROW = IRROW + 1 3065 IRINDX(IRROW) = IDIM1 3066 ELSE 3067 IZROW = IZROW + 1 3068 IZINDX(IZROW) = IDIM1 3069 END IF 3070 100 CONTINUE 3071 NROW = IRROW 3072 NZROW = IZROW 3073C 3074C *** removal of the rows. *** 3075C 3076 MXDIM = NDIM 3077 DO 300 IROW = 1, NZROW 3078 DO 400 IDIM2 = 1, NDIM 3079 DO 400 IDIM1 = IZINDX(IROW)-IROW+1, MXDIM 3080 RMTRIX(IDIM1,IDIM2) = RMTRIX(IDIM1+1,IDIM2) 3081 400 CONTINUE 3082 MXDIM = MXDIM - 1 3083 300 CONTINUE 3084C 3085C *** removal of the corresponding coloumns *** 3086C 3087 MXDIM = NDIM 3088 DO 500 ICOL = 1, NZROW 3089 DO 600 IDIM2 = IZINDX(ICOL)-ICOL+1, MXDIM 3090 DO 600 IDIM1 = 1, NDIM 3091 RMTRIX(IDIM1,IDIM2) = RMTRIX(IDIM1,IDIM2+1) 3092 600 CONTINUE 3093 MXDIM = MXDIM - 1 3094 500 CONTINUE 3095C 3096C *** Test print after removal of rows *** 3097C 3098 IF (IPRINT .GT. 20) THEN 3099 WRITE (LUPRI,'(A)') 'Matrix after removal of rows' 3100 DO II = 1, NROW 3101 WRITE (LUPRI,'(16F8.4)') (RMTRIX(II,IJ),IJ=1,NROW) 3102 END DO 3103 END IF 3104C 3105 RETURN 3106 END 3107C 3108C 3109C /* Deck chkcmp */ 3110 SUBROUTINE CHKCMP(KDPMTX,IDXTST,KORDR,LDPMTX,IFRSTD,NSTART,NLDPMX, 3111 & NIDEP,IPRINT,EXIST) 3112C ************************************************************************* 3113C *** This subroutine checks whether the force constant IDXTST exist in *** 3114C *** KDPMTX. It returns exist = .true. if this is the case. *** 3115C ************************************************************************* 3116#include "implicit.h" 3117#include "priunit.h" 3118C 3119#include "numder.h" 3120#include "fcsym.h" 3121 LOGICAL EXIST 3122 DIMENSION KDPMTX(LDPMTX,NSTRDR,IFRSTD), IDXTST(NMORDR) 3123C 3124 DO 100 IDP = 1, NIDEP+1 3125 DO 100 II = NSTART, NLDPMX 3126 EXIST = .TRUE. 3127 DO 200 IORDR = 1, KORDR 3128 EXIST = EXIST.AND.(KDPMTX(II,IORDR,IDP).EQ.IDXTST(IORDR)) 3129 200 CONTINUE 3130 IF (EXIST) GOTO 300 3131 100 CONTINUE 3132C 3133C *** if it gets here with exist = .true. (Through the goto), *** 3134C *** the component exist. *** 3135C 3136 300 CONTINUE 3137C 3138 IF (IPRINT.GT.20) THEN 3139 WRITE (LUPRI,'(A,I4,A,I4)') 'Search startst at',NSTART, 3140 & ' and ends at', NLDPMX 3141 WRITE (LUPRI,'(A,4I4)') 'The component', 3142 & (IDXTST(II),II=1,KORDR) 3143 WRITE (LUPRI,'(L1)') EXIST 3144 END IF 3145C 3146 RETURN 3147 END 3148C 3149C 3150C /* Deck chklcp*/ 3151 SUBROUTINE CHKLCP(DCOEFF,KDPMTX,ICRIRP,IDXTMP,INMIDP,IMVTMP, 3152 & ICMPID,LDPMTX,IFRSTD,NDMCMP,NUMCMP,NSTART,NIDEP, 3153 & IPRINT) 3154C ********************************************************************** 3155C *** This subroutine checks if the independent component chosen *** 3156C *** is the component with the least energy calculations. If this *** 3157C *** is not the case, the components are rearranged such that it is *** 3158C ********************************************************************** 3159#include "implicit.h" 3160#include "priunit.h" 3161#include "mxcent.h" 3162 PARAMETER (D1 = 1.0D0, DTHR = 1.0D-6) 3163#include "trkoor.h" 3164#include "numder.h" 3165#include "fcsym.h" 3166 LOGICAL INCMNT, FOUND 3167 DIMENSION DCOEFF(LDPMTX,IFRSTD) 3168 DIMENSION KDPMTX(LDPMTX,NSTRDR,IFRSTD), ICRIRP(NCOOR,2), 3169 & IDXTMP(NDMCMP), ICMPID(NDMCMP), INMIDP(NIDEP), 3170 & IMVTMP(NDMCMP) 3171C 3172C *** Test print *** 3173C 3174 IF (IPRINT .GT. 20) THEN 3175 CALL HEADER('Coordinate dependency before change.',-1) 3176 WRITE (LUPRI,'(5X,A)') 'Dependent components ' // 3177 & 'Independent components Coeffisient' 3178 DO II = NSTART, NSTART+NUMCMP-NIDEP 3179 IF (NDMCMP.EQ.3) THEN 3180 WRITE (LUPRI,FMT=1003) (KDPMTX(II,IX,1),IX=1,NDMCMP), 3181 & (KDPMTX(II,IX,2),IX=1,NDMCMP), DCOEFF(II,1) 3182 DO IDP = 2, NIDEP 3183 WRITE (LUPRI,FMT=1006) 3184 & (KDPMTX(II,IX,IDP+1),IX=1,NDMCMP), DCOEFF(II,IDP) 3185 END DO 3186 ELSE IF (NDMCMP.EQ.4) THEN 3187 write (lupri,'(a,I4)') 'testkomp', II 3188 WRITE (LUPRI,FMT=1004) (KDPMTX(II,IX,1),IX=1,NDMCMP), 3189 & (KDPMTX(II,IX,2),IX=1,NDMCMP), DCOEFF(II,1) 3190 DO IDP = 2, NIDEP 3191 WRITE (LUPRI,FMT=1008) 3192 & (KDPMTX(II,IX,IDP+1),IX=1,NDMCMP), DCOEFF(II,IDP) 3193 END DO 3194 END IF 3195 END DO 3196 WRITE (LUPRI,'(A)')' ' 3197 WRITE (LUPRI,'(A)')' ' 3198 END IF 3199C 3200C *** Finding the number of alike variables in the *** 3201C *** independent variables. *** 3202C 3203 DO 100 IDP = 1, NIDEP 3204C 3205 CALL IZERO(IDXTMP,NDMCMP) 3206C 3207 INUM = 1 3208 IDXTMP(1) = 1 3209 ICMPID(1) = KDPMTX(NSTART,1,IDP+1) 3210 DO 200 II = 2, NDMCMP 3211 INCMNT = .TRUE. 3212 DO 300 IJ = 1, INUM 3213 IF (KDPMTX(NSTART,II,IDP+1).EQ.ICMPID(IJ)) THEN 3214 IDXTMP(IJ) = IDXTMP(IJ) + 1 3215 INCMNT = .FALSE. 3216 END IF 3217 300 CONTINUE 3218C 3219 IF (INCMNT) THEN 3220 INUM = INUM + 1 3221 IDXTMP(INUM) = 1 3222 ICMPID(INUM) = KDPMTX(NSTART,II,IDP+1) 3223 END IF 3224 200 CONTINUE 3225 NUMB = INUM 3226C 3227 INMIDP(IDP) = 0 3228 DO 400 II = 1, NUMB 3229 INMIDP(IDP) = MAX(INMIDP(IDP),IDXTMP(II)) 3230 400 CONTINUE 3231C 3232 100 CONTINUE 3233C 3234C *** Comparing this with the dependent variables. *** 3235C 3236 DO 500 ICMP = NSTART, NSTART+NUMCMP-NIDEP 3237 CALL IZERO(IDXTMP,NDMCMP) 3238C 3239 INUM = 1 3240 IDXTMP(1) = 1 3241 ICMPID(1) = KDPMTX(ICMP,1,1) 3242 DO 600 II = 2, NDMCMP 3243 INCMNT = .TRUE. 3244 DO 700 IJ = 1, INUM 3245 IF (KDPMTX(ICMP,II,1).EQ.ICMPID(IJ)) THEN 3246 IDXTMP(IJ) = IDXTMP(IJ) + 1 3247 INCMNT = .FALSE. 3248 END IF 3249 700 CONTINUE 3250C 3251 IF (INCMNT) THEN 3252 INUM = INUM + 1 3253 IDXTMP(INUM) = 1 3254 ICMPID(INUM) = KDPMTX(ICMP,II,1) 3255 END IF 3256 600 CONTINUE 3257 NUMB = INUM 3258C 3259 MAXDP = 0 3260 DO 800 II = 1, NUMB 3261 MAXDP = MAX(MAXDP,IDXTMP(II)) 3262 800 CONTINUE 3263C 3264 FOUND = .FALSE. 3265 DO 900 IDP = 1, NIDEP 3266 IF ((MAXDP.GT.INMIDP(IDP)).AND.(.NOT.FOUND).AND. 3267 & (ABS(DCOEFF(ICMP,IDP)).GT.1.0D-4)) THEN 3268 FOUND = .TRUE. 3269 IMIN = 1 3270 NMIN = INMIDP(1) 3271 DO 1100 IDP2 = 2, NIDEP 3272 IF (INMIDP(IDP2).LT.NMIN) THEN 3273 IMIN = IDP2 3274 NMIN = INMIDP(IDP2) 3275 END IF 3276 1100 CONTINUE 3277C 3278C *** Swapping coordinates with respect to this criteria *** 3279C 3280 DO 1200 II = 1, NDMCMP 3281 IMVTMP(II) = KDPMTX(ICMP,II,1) 3282 KDPMTX(ICMP,II,1) = KDPMTX(NSTART,II,IMIN+1) 3283 DO 1300 IST = NSTART, NSTART+NUMCMP-NIDEP 3284 KDPMTX(IST,II,IMIN+1) = IMVTMP(II) 3285 1300 CONTINUE 3286 1200 CONTINUE 3287 INMIDP(IDP) = MAXDP 3288C 3289C *** Arranging the coeffisients for the new and *** 3290C *** better point. And placing the old as *** 3291C *** independent. *** 3292C 3293 BCOEFF = D1/DCOEFF(ICMP,IMIN) 3294 DO 1400 IDP2 = 1, NIDEP 3295 IF (IDP2 .EQ. IMIN) THEN 3296 DCOEFF(ICMP,IDP2) = BCOEFF 3297 ELSE 3298 DCOEFF(ICMP,IDP2) = -DCOEFF(ICMP,IDP2)*BCOEFF 3299 END IF 3300 1400 CONTINUE 3301C 3302C *** Changing the other coeffisients to fit the *** 3303C *** new and better point *** 3304C 3305 DO 1500 ICMP2 = NSTART, NSTART+NUMCMP-NIDEP 3306 IF (ICMP2.NE.ICMP) THEN 3307 DO 1600 IDP2 = 1, NIDEP 3308 IF (IDP2 .NE. IMIN) THEN 3309 DCOEFF(ICMP2,IDP2) = DCOEFF(ICMP2,IDP2) 3310 & + DCOEFF(ICMP2,IMIN) 3311 & *DCOEFF(ICMP ,IDP2) 3312 END IF 3313 1600 CONTINUE 3314 DCOEFF(ICMP2,IMIN) = DCOEFF(ICMP2,IMIN)*BCOEFF 3315 END IF 3316 1500 CONTINUE 3317 3318 IF (IPRINT .GT. 20) THEN 3319 CALL HEADER('Coordinate dependency after 1. change.', 3320 & -1) 3321 WRITE (LUPRI,'(5X,A)') 'Dependent components ' // 3322 & 'Independent components Coeffisient' 3323 DO II = NSTART, NSTART+NUMCMP-NIDEP 3324 IF (NDMCMP.EQ.3) THEN 3325 WRITE (LUPRI,FMT=1003) 3326 & (KDPMTX(II,IX,1),IX=1,NDMCMP), 3327 & (KDPMTX(II,IX,2),IX=1,NDMCMP), DCOEFF(II,1) 3328 DO IDP2 = 2, NIDEP 3329 WRITE (LUPRI,FMT=1006) 3330 & (KDPMTX(II,IX,IDP2+1),IX=1,NDMCMP), 3331 & DCOEFF(II,IDP2) 3332 END DO 3333 ELSE IF (NDMCMP.EQ.4) THEN 3334 WRITE (LUPRI,FMT=1004) 3335 & (KDPMTX(II,IX,1),IX=1,NDMCMP), 3336 & (KDPMTX(II,IX,2),IX=1,NDMCMP), DCOEFF(II,1) 3337 DO IDP2 = 2, NIDEP 3338 WRITE (LUPRI,FMT=1008) 3339 & (KDPMTX(II,IX,IDP2+1),IX=1,NDMCMP), 3340 & DCOEFF(II,IDP2) 3341 END DO 3342 END IF 3343 END DO 3344 WRITE (LUPRI,'(A)')' ' 3345 WRITE (LUPRI,'(A)')' ' 3346 END IF 3347 END IF 3348 900 CONTINUE 3349C 3350C *** Checking wether the component is better with respect to *** 3351C *** which row it belongs to (row 1 of a 2 dimensional irrep *** 3352C *** is preffered. *** 3353C 3354 IF (.NOT.FOUND) THEN 3355C 3356C *** Finds the worst independent variable available *** 3357C 3358 MXIDEP = 0 3359 MXIDX = 0 3360 DO 1700 IDP = 1, NIDEP 3361 IIDEP = 0 3362 DO 1800 II = 1, NDMCMP 3363 IIDEP = IIDEP + ICRIRP(KDPMTX(ICMP,II,IDP+1),2) 3364 1800 CONTINUE 3365 IF ((IIDEP.GT.MXIDEP) .AND. 3366 & (ABS(DCOEFF(ICMP,IDP)).GT.DTHR)) THEN 3367 MXIDX = IDP 3368 MXIDEP = IIDEP 3369 END IF 3370 1700 CONTINUE 3371C 3372 IDEP = 0 3373 DO 1900 II = 1, NDMCMP 3374 IDEP = IDEP + ICRIRP(KDPMTX(ICMP,II,1),2) 3375 1900 CONTINUE 3376C 3377C *** If mxidep .gt. idep then the number coordinates that *** 3378C *** transforms as row 2 i larger in iidep then in idep. *** 3379C *** We then change places. *** 3380C 3381 IF (MXIDEP .GT. IDEP) THEN 3382 FOUND = .TRUE. 3383C 3384 DO 2100 II = 1, NDMCMP 3385 IMVTMP(II) = KDPMTX(ICMP,II,1) 3386 KDPMTX(ICMP,II,1) = KDPMTX(NSTART,II,MXIDX+1) 3387 DO 2200 IST = NSTART, NSTART+NUMCMP-NIDEP 3388 KDPMTX(IST,II,MXIDX+1) = IMVTMP(II) 3389 2200 CONTINUE 3390 2100 CONTINUE 3391C 3392C *** Arranging the coeffisients for the new and *** 3393C *** better point. And placing the old as *** 3394C *** independent. *** 3395C 3396 BCOEFF = D1/DCOEFF(ICMP,MXIDX) 3397 DO 2300 IDP2 = 1, NIDEP 3398 IF (IDP2 .EQ. MXIDX) THEN 3399 DCOEFF(ICMP,IDP2) = BCOEFF 3400 ELSE 3401 DCOEFF(ICMP,IDP2) = -DCOEFF(ICMP,IDP2)*BCOEFF 3402 END IF 3403 2300 CONTINUE 3404C 3405C *** Changing the other coeffisients to fit the *** 3406C *** new and better point *** 3407C 3408 DO 2400 ICMP2 = NSTART, NSTART+NUMCMP-NIDEP 3409 IF (ICMP2.NE.ICMP) THEN 3410 DO 2500 IDP2 = 1, NIDEP 3411 IF (IDP2 .NE. MXIDX) THEN 3412 DCOEFF(ICMP2,IDP2) = DCOEFF(ICMP2,IDP2 ) 3413 & + DCOEFF(ICMP2,MXIDX) 3414 & *DCOEFF(ICMP ,IDP2 ) 3415 END IF 3416 2500 CONTINUE 3417 DCOEFF(ICMP2,MXIDX)=DCOEFF(ICMP2,MXIDX)*BCOEFF 3418 END IF 3419 2400 CONTINUE 3420 END IF 3421 END IF 3422 IF (IPRINT .GT. 20) THEN 3423 CALL HEADER('Coordinate dependency after 2. change. ',-1) 3424 WRITE (LUPRI,'(5X,A)') 'Dependent components ' // 3425 & 'Independent components Coeffisient' 3426 DO II = NSTART, NSTART+NUMCMP-NIDEP 3427 IF (NDMCMP.EQ.3) THEN 3428 WRITE (LUPRI,FMT=1003) (KDPMTX(II,IX,1),IX=1,NDMCMP), 3429 & (KDPMTX(II,IX,2),IX=1,NDMCMP), DCOEFF(II,1) 3430 DO IDP = 2, NIDEP 3431 WRITE (LUPRI,FMT=1006) 3432 & (KDPMTX(II,IX,IDP+1),IX=1,NDMCMP), DCOEFF(II,IDP) 3433 END DO 3434 ELSE IF (NDMCMP.EQ.4) THEN 3435 WRITE (LUPRI,FMT=1004) (KDPMTX(II,IX,1),IX=1,NDMCMP), 3436 & (KDPMTX(II,IX,2),IX=1,NDMCMP), DCOEFF(II,1) 3437 DO IDP = 2, NIDEP 3438 WRITE (LUPRI,FMT=1008) 3439 & (KDPMTX(II,IX,IDP+1),IX=1,NDMCMP), DCOEFF(II,IDP) 3440 END DO 3441 END IF 3442 END DO 3443 WRITE (LUPRI,'(A)')' ' 3444 WRITE (LUPRI,'(A)')' ' 3445 END IF 3446 500 CONTINUE 3447C 3448 1003 FORMAT (10X,3I4,10X,3I4,15X,F8.4) 3449 1006 FORMAT (32X,3I4,15X,F8.4) 3450 1004 FORMAT (6X,4I4,10X,4I4,11X,F8.4) 3451 1008 FORMAT (28X,4I4,15X,F8.4) 3452C 3453 RETURN 3454 END 3455C 3456C 3457C /* Deck srtncf */ 3458 SUBROUTINE SRTNCF(EQUMTX,ITMPRM,ICRWMX,ITMPEQ,NRDIM,KORDR,LDIM, 3459 & IPRINT) 3460C **************************************************** 3461C *** Sorting the components such that they become *** 3462C *** calculated ones (X1 .GE. X2 .GE. X3 ....) *** 3463C **************************************************** 3464#include "implicit.h" 3465#include "priunit.h" 3466C 3467 LOGICAL FOUND, EXIST 3468 DIMENSION EQUMTX(NRDIM,NRDIM) 3469 DIMENSION ICRWMX(NRDIM,KORDR), ITMPRM(KORDR), ITMPEQ(NRDIM) 3470C 3471C *** Test print before for comparison *** 3472C 3473 IF (IPRINT .GE. 20) THEN 3474 WRITE (LUPRI,'(5X,A)') 'Components before sorting' 3475C 3476 DO II =1 , LDIM 3477 WRITE (LUPRI,'(10X,20I4)') (ICRWMX(II,IJ),IJ=1,KORDR) 3478 END DO 3479 END IF 3480C 3481 DO 100 II = 1, LDIM 3482C 3483C *** Sorting such that the indeces are aligned *** 3484C *** in descending order *** 3485C 3486 IMAX = 0 3487 ITMPRM(1) = ICRWMX(II,1) 3488 DO 200 IJ = 2, KORDR 3489C 3490 FOUND = .FALSE. 3491 IMAX = IMAX + 1 3492 DO 300 IK = 1, IMAX 3493 IF ((ITMPRM(IK).LT.ICRWMX(II,IJ)).AND.(.NOT.FOUND)) THEN 3494 FOUND = .TRUE. 3495 DO 400 IL = IMAX, IK, -1 3496 ITMPRM(IL+1) = ITMPRM(IL) 3497 400 CONTINUE 3498 ITMPRM(IK) = ICRWMX(II,IJ) 3499 END IF 3500 300 CONTINUE 3501 IF (.NOT.FOUND) THEN 3502 ITMPRM(IMAX+1) = ICRWMX(II,IJ) 3503 END IF 3504 200 CONTINUE 3505C 3506C *** Copying it back to the index array *** 3507C 3508 DO 500 IJ = 1, KORDR 3509 ICRWMX(II,IJ) = ITMPRM(IJ) 3510 500 CONTINUE 3511C 3512 100 CONTINUE 3513C 3514C *** Removing redundant rows (with respect to equal force *** 3515C *** constants) in the equation matrix *** 3516C 3517 IREM = 0 3518 LSTART = LDIM 3519 DO 600 II = 1, LSTART 3520C 3521C *** Finding out if this force constant really already *** 3522C *** exist, and which force constant this is. *** 3523C 3524 EXIST = .FALSE. 3525 DO 700 IJ = 1, II-1 3526 IF (.NOT.EXIST) THEN 3527 EXIST = .TRUE. 3528 DO 800 IK = 1, KORDR 3529 IF (ICRWMX(II,IK).NE.ICRWMX(IJ,IK)) EXIST = .FALSE. 3530 800 CONTINUE 3531 IF (EXIST) IEXIST = IJ 3532 END IF 3533 700 CONTINUE 3534C 3535C *** If this component exist the matrix elements of *** 3536C *** component will be added to the old. *** 3537C 3538 IF (EXIST) THEN 3539 DO 900 IK = 1, LSTART 3540 EQUMTX(IK,IEXIST) = EQUMTX(IK,IEXIST) + EQUMTX(IK,II) 3541 900 CONTINUE 3542 IREM = IREM + 1 3543 ITMPEQ(IREM) = II 3544 END IF 3545 600 CONTINUE 3546C 3547C *** Removing the appropriate rows and coloumns. *** 3548C 3549 DO 1100 IK = 1, IREM 3550 DO 1200 IJ = ITMPEQ(IK), LDIM-1 3551 DO 1200 II = 1 , LDIM 3552 EQUMTX(II,IJ) = EQUMTX(II,IJ+1) 3553 1200 CONTINUE 3554C 3555 DO 1300 IJ = 1 , LDIM 3556 DO 1300 II = ITMPEQ(IK), LDIM-1 3557 EQUMTX(II,IJ) = EQUMTX(II+1,IJ) 3558 1300 CONTINUE 3559C 3560 DO 1400 IJ = 1, KORDR 3561 DO 1400 II = ITMPEQ(IK), LDIM 3562 ICRWMX(II,IJ) = ICRWMX(II+1,IJ) 3563 1400 CONTINUE 3564C 3565 DO 1500 II = IK+1, IREM 3566 ITMPEQ(II) = ITMPEQ(II) - 1 3567 1500 CONTINUE 3568C 3569 LDIM = LDIM - 1 3570 1100 CONTINUE 3571C 3572C *** Test print *** 3573C 3574 IF (IPRINT .GE. 20) THEN 3575C 3576 WRITE (LUPRI,'(5X,A,I4)') 'Size of matrix:', LDIM 3577 WRITE (LUPRI,'(5X,A)') 'Components after sorting' 3578C 3579 DO II =1 , LDIM 3580 WRITE (LUPRI,'(10X,20I4)') (ICRWMX(II,IJ),IJ=1,KORDR) 3581 END DO 3582C 3583 WRITE (LUPRI,'(A)') 'Matrix after sorting:' 3584 DO II = 1, LDIM 3585 WRITE (LUPRI,'(10X,20F8.4)') (EQUMTX(II,IJ),IJ=1,LDIM) 3586 END DO 3587 END IF 3588C 3589 RETURN 3590 END 3591C 3592C 3593C /* Deck prtnrp */ 3594 SUBROUTINE PRTNRP(GRIREP,WORK,ICRIRP,INDSTP,NPRTNR,MAXINR,IINNER, 3595 & IORDR,IRSRDR,NMPRTN,IDXTMP,INDTMP,LWORK,CLNRGY, 3596 & PRTNR,ALRCAL) 3597#include "implicit.h" 3598#include "priunit.h" 3599 PARAMETER (D0 = 0.0D0, DTHRSH = 1.0D-6, MXIDCR=5) 3600#include "fcsym.h" 3601#include "numder.h" 3602 LOGICAL PRTNR, CLNRGY, C2NRGY, TSM, ALRCAL 3603 DIMENSION GRIREP(NGORDR,NGVERT), WORK(LWORK) 3604 DIMENSION ICRIRP(NDCOOR,2), INDSTP(NTORDR), NPRTNR(MAXINR), 3605 & INDTMP(MXIDCR), IDXTMP(NMORDR) 3606C 3607 PRTNR = .TRUE. 3608C 3609C *** Checking wether this energy is already calculated *** 3610C *** using another geometry as "partner-geometry. *** 3611C 3612 ALRCAL = .FALSE. 3613 DO 100 INM = 1, NMPRTN 3614 IF ((NPRTNR(INM).EQ.IINNER).AND.(.NOT.ALRCAL)) THEN 3615 CLNRGY = .FALSE. 3616 PRTNR = .FALSE. 3617 ALRCAL = .TRUE. 3618 END IF 3619 100 CONTINUE 3620C 3621C *** Checking if all stepping coordinates are tot. sym. *** 3622C 3623 ITST = 1 3624 TSM = .FALSE. 3625 DO 200 I = 1, IRSRDR+1 3626 ITST = ITST*ICRIRP(INDSTP(I),1) 3627 200 CONTINUE 3628 IF (ITST.EQ.1) TSM = .TRUE. 3629C 3630C *** Checking if this geometry has a partner geometry *** 3631C 3632 IF ((IORDR.LE.2) .AND. (NMORDR.EQ.2) .AND.PRTNR 3633 & .AND..NOT.TSM) THEN 3634C 3635C *** Stepping in two directions. *** 3636C 3637 CALL SCNPRT(GRIREP,WORK,INDSTP,ICRIRP,INDTMP,IORDR,LWORK,PRTNR) 3638c ELSE ((IORDR.LE.3) .AND. (PRTNR).AND.(.NOT.TSM)) THEN 3639c 3640 ELSE 3641 PRTNR = .FALSE. 3642 END IF 3643C 3644C *** If there is a partner geometry, it's about *** 3645C *** time to find it. *** 3646C 3647 IF (PRTNR) THEN 3648 NMPRTN = NMPRTN + 1 3649 ITINNR = IINNER - 1 3650 DO 400 I = IORDR, 1, -1 3651 IS = 2**(I-1) 3652 IDXTMP(I) = (ITINNR-MOD(ITINNR,IS))/IS 3653 ITINNR = ITINNR - IDXTMP(I)*IS 3654 400 CONTINUE 3655C 3656 DO 500 I = 1, IORDR 3657 IF (ICRIRP(INDSTP(I),1).NE.1) THEN 3658 NPRTNR(NMPRTN) = NPRTNR(NMPRTN) 3659 & + ABS(IDXTMP(I)-1)*2**(I-1) 3660 END IF 3661 500 CONTINUE 3662 NPRTNR(NMPRTN) = NPRTNR(NMPRTN) + 1 3663C 3664 END IF 3665C 3666 RETURN 3667 END 3668C 3669C 3670 SUBROUTINE SCNPRT(GRIREP,WORK,INDSTP,ICRIRP,INDTMP,IORDR,LWORK, 3671 & PRTNR) 3672C ******************************************************** 3673C *** Checking if this geometry has a partner geometry *** 3674C *** Sufficient demand for stepping in 1-2 directions *** 3675C *** is that all possible F_{x_1,x_2,x_3} in taylor *** 3676C *** expansion is zero. *** 3677C ******************************************************** 3678#include "implicit.h" 3679#include "priunit.h" 3680 PARAMETER (MXIDCR=5) 3681#include "fcsym.h" 3682#include "numder.h" 3683 LOGICAL PRTNR, C2NRGY 3684 DIMENSION GRIREP(NGORDR,NGVERT), WORK(LWORK) 3685 DIMENSION ICRIRP(NDCOOR,2), INDSTP(NTORDR), INDTMP(MXIDCR) 3686C 3687 DO 100 IIRDR3 = 1, IORDR 3688 DO 100 IIRDR2 = 1, IIRDR3 3689 DO 100 IIRDR1 = 1, IIRDR2 3690 IIREP3 = ICRIRP(INDSTP(IIRDR3),1) 3691 IIREP2 = ICRIRP(INDSTP(IIRDR2),1) 3692 IIREP1 = ICRIRP(INDSTP(IIRDR1),1) 3693C 3694 IF ((IIREP3.NE.1).OR.(IIREP2.NE.1).OR.(IIREP1.NE.1)) THEN 3695C 3696 INDTMP(1) = INDSTP(IIRDR1) 3697 INDTMP(2) = INDSTP(IIRDR2) 3698 INDTMP(3) = INDSTP(IIRDR3) 3699C 3700 KIRPST = 1 3701 KIRPDG = KIRPST + MXIDCR 3702 KKIRPD = KIRPDG + MXIDCR 3703 KIDDEG = KKIRPD + MXIDCR 3704 KIDEGI = KIDDEG + MXIDCR 3705 KLAST = KIDEGI + MXIDCR 3706 IF (KLAST.GT.LWORK) CALL QUIT('Memory exceeded in' // 3707 & ' PRTPNR.') 3708C 3709 C2NRGY = .FALSE. 3710 CALL SYMMLT(GRIREP,INDTMP,ICRIRP,WORK(KIRPST),WORK(KIRPDG), 3711 & WORK(KKIRPD),WORK(KIDDEG),WORK(KIDEGI),MXIDCR, 3712 & NTORDR,2,3,IPRINT,C2NRGY) 3713 PRTNR = PRTNR .AND. .NOT. C2NRGY 3714 END IF 3715 100 CONTINUE 3716C 3717 RETURN 3718 END 3719