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 anainp */ 20 SUBROUTINE ANAINP(WORD) 21C 22C 5-Jul-1985 Hans Jorgen Aa. Jensen 23C 24#include "implicit.h" 25#include "priunit.h" 26#include "mxcent.h" 27 PARAMETER (NTABLE = 4) 28 PARAMETER (MAXANG = 20) 29 LOGICAL NEWDEF 30 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 31#include "abainf.h" 32 LOGICAL SKIP 33 COMMON /CBIANA/ IANG(3,MAXANG),IDIHED(4,MAXANG),NANG,NDIHED, 34 * SKIP 35C 36 DATA TABLE /'.SKIP ', '.XXXXXX', '.ANGLES', '.DIHEDR'/ 37 DATA MANG/0/, MDIHED/0/ 38C 39 CALL ANAINI 40C 41 NEWDEF = (WORD .EQ. '*GEOANA') 42 ICHANG = 0 43 IF (NEWDEF) THEN 44 WORD1 = WORD 45 100 CONTINUE 46 READ (LUCMD, '(A7)') WORD 47 CALL UPCASE(WORD) 48 PROMPT = WORD(1:1) 49 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 50 GO TO 100 51 ELSE IF (PROMPT .EQ. '.') THEN 52 ICHANG = ICHANG + 1 53 DO 200 I = 1, NTABLE 54 IF (TABLE(I) .EQ. WORD) THEN 55 GO TO (1,2,3,4), I 56 END IF 57 200 CONTINUE 58 IF (WORD .EQ. '.OPTION') THEN 59 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 60 GO TO 100 61 END IF 62 WRITE (LUPRI,'(/,3A,/)') ' Keyword "',WORD, 63 * '" not recognized in ANAINP.' 64 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 65 CALL QUIT('Illegal keyword in ANAINP.') 66 1 CONTINUE 67 SKIP = .TRUE. 68 GO TO 100 69 2 CONTINUE 70 GO TO 100 71 3 CONTINUE 72 READ (LUCMD,*) NANG 73 MANG = MIN(MAXANG,NANG) 74 DO 310 I = 1,MANG 75 READ(LUCMD,*) (IANG(J,I),J=1,3) 76 310 CONTINUE 77 MANG = NANG - MANG 78 DO 320 I = 1,MANG 79 READ(LUCMD,'()') 80 320 CONTINUE 81 GO TO 100 82 4 CONTINUE 83 READ (LUCMD,*) NDIHED 84 MDIHED = MIN(MAXANG,NDIHED) 85 DO 410 I = 1,MDIHED 86 READ(LUCMD,*) (IDIHED(J,I),J=1,4) 87 410 CONTINUE 88 MDIHED = NDIHED - MDIHED 89 DO 420 I = 1,MDIHED 90 READ(LUCMD,'()') 91 420 CONTINUE 92 GO TO 100 93 ELSE IF (PROMPT .EQ. '*') THEN 94 GO TO 300 95 ELSE 96 WRITE (LUPRI,'(/,3A,/)') ' Prompt "',WORD, 97 * '" not recognized in ANAINP.' 98 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 99 CALL QUIT('Illegal prompt in ANAINP.') 100 END IF 101 END IF 102 300 CONTINUE 103 IF (ICHANG .GT. 0) THEN 104 CALL HEADER('Changes of defaults for GEOANA:',0) 105 IF (SKIP) THEN 106 WRITE (LUPRI,'(A)') ' GEOANA skipped in this run.' 107 ELSE 108 IF (NANG .GT. 0) THEN 109 WRITE (LUPRI,'(/A/)') 110 * ' Following angles will be calculated:' 111 DO 1310 I = 1,NANG 112 WRITE (LUPRI,'(I10,A,4I5)') I,' : ',(IANG(J,I),J=1,3) 113 1310 CONTINUE 114 IF (MANG .GT. 0) THEN 115 WRITE (LUPRI,'(/A,I3,A)') ' The last',MANG, 116 * ' angles specified go beyond current maximum', 117 * ' and will not be printed.' 118 END IF 119 END IF 120 IF (NDIHED .GT. 0) THEN 121 WRITE (LUPRI,'(/A/)') 122 * ' Following dihedral angles will be calculated:' 123 DO 1410 I = 1,NDIHED 124 WRITE (LUPRI,'(I10,A,4I5)')I,' : ',(IDIHED(J,I),J=1,4) 125 1410 CONTINUE 126 IF (MDIHED .GT. 0) THEN 127 WRITE (LUPRI,'(/A,I3,A)') ' The last',MDIHED, 128 * ' dihedral angles specified go beyond current', 129 * ' maximum and will not be printed.' 130 END IF 131 END IF 132 END IF 133 WRITE (LUPRI,'(/)') 134 END IF 135 RETURN 136 END 137C /* Deck anaini */ 138 SUBROUTINE ANAINI 139C 140C Initialize /CBIANA/ 141C 142#include "implicit.h" 143 PARAMETER (MAXANG = 20) 144 LOGICAL SKIP 145 COMMON /CBIANA/ IANG(3,MAXANG),IDIHED(4,MAXANG),NANG,NDIHED, 146 * SKIP 147C 148 NANG = 0 149 NDIHED = 0 150 SKIP = .FALSE. 151 RETURN 152 END 153C /* Deck geoana */ 154 SUBROUTINE GEOANA(COORD,PRINT,DIF,NBONDS,LUPUNCH,WORK,LWORK) 155C 156C 30-Jun-1985 Hans Jorgen Aa. Jensen 157C Modified for symmetry 25-Sep-1989 tuh 158C Modified for differential geomtries 18-Oct-1989 tuh 159C 160#include "implicit.h" 161#include "maxaqn.h" 162#include "mxcent.h" 163#include "maxorb.h" 164 LOGICAL PRINT, DIF 165 DIMENSION COORD(*),WORK(LWORK) 166#include "nuclei.h" 167C 168 CALL QENTER('GEOANA') 169 KFREE = 1 170 LFREE = LWORK 171 CALL MEMGET('REAL',KVEC ,3*NUCDEP*NUCDEP,WORK,KFREE,LFREE) 172 CALL MEMGET('LOGI',KBOND ,NUCDEP*NUCDEP,WORK,KFREE,LFREE) 173 CALL MEMGET('INTE',KCHRG ,NUCDEP, WORK,KFREE,LFREE) 174 CALL MEMGET('INTE',KPAIR,2*NUCDEP*NUCDEP,WORK,KFREE,LFREE) 175C 176 CALL GEOANA_1(COORD,PRINT,DIF,NBONDS,LUPUNCH,WORK(KVEC), 177 & WORK(KBOND),WORK(KCHRG),WORK(KPAIR)) 178C 179 CALL MEMREL('GEOANA',WORK,1,1,KFREE,LFREE) 180 CALL QEXIT('GEOANA') 181 RETURN 182 END 183C /* Deck GEOANA_1 */ 184 SUBROUTINE GEOANA_1(COORD,PRINT,DIF,NBONDS,LUPUNCH,VEC, 185 & BONDED,ICHARG,IPAIRS) 186C 187C Modified for more selective printing of bonded atoms, 188C Jan-1995 Hanne Heiberg 189C LUPUN .gt. 0: punching atom bonds for Gamess graphic output, K.Ruud-95 190C 191#include "implicit.h" 192#include "priunit.h" 193#include "maxaqn.h" 194#include "mxcent.h" 195#include "maxorb.h" 196#include "codata.h" 197#include "facang.h" 198 PARAMETER (MAXANG = 20) 199 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, D1P5 = 1.5D00) 200 LOGICAL SKIP, DIF, PRINT 201 COMMON /CBIANA/ IANG(3,MAXANG),IDIHED(4,MAXANG),NANG,NDIHED, 202 * SKIP 203C 204#ifdef PRG_DIRAC 205#include "dcbgen.h" 206#else 207#include "gnrinf.h" 208#endif 209#include "nuclei.h" 210#include "symmet.h" 211#include "qm3.h" 212C 213 DIMENSION COORD(3,*), DIST(MXCENT*(MXCENT+1)/2), 214 & ANGLE(MAXANG), DIHED(MAXANG), ICHARG(NUCDEP) 215 DIMENSION VEC(3,NUCDEP,NUCDEP), IPAIRS(2,NUCDEP*NUCDEP) 216 LOGICAL BONDED(NUCDEP,NUCDEP) 217 CHARACTER*6 NUCNAM(4) 218 SAVE DIST, ANGLE, DIHED 219C statement function: 220 ARCCOS(ARG) = FACANG*ACOS(ARG) 221 222C 223#ifndef PRG_DIRAC 224 IF (REDCNT) RETURN 225! ... REDCNT is a QM3 variable 226#endif 227 228 NBONDS = 0 229 IF (NUCDEP .EQ. 1) RETURN 230C 231C set up bond vectors in Angstrom in VEC 232C NUCIND is number of QM atoms 233C NCTOT .ge. NUCIND is number of QM atoms + number of MM atoms 234C 235 N_QMATOMS = 0 236 IATOMA = 0 237 DO 100 ICENTA = 1, NCTOT 238 DO 110 IA = 0, MAXOPR 239 IF (IAND(IA,ISTBNU(ICENTA)) .EQ. 0) THEN 240 IATOMA = IATOMA + 1 241 CXA = PT(IAND(ISYMAX(1,1),IA))*COORD(1,ICENTA) 242 CYA = PT(IAND(ISYMAX(2,1),IA))*COORD(2,ICENTA) 243 CZA = PT(IAND(ISYMAX(3,1),IA))*COORD(3,ICENTA) 244C 245 ICHARG(IATOMA) = IZATOM(ICENTA) 246C 247 IATOMB = 0 248 DO 200 ICENTB = 1, NCTOT 249 DO 210 IB = 0, MAXOPR 250 IF (IAND(IB,ISTBNU(ICENTB)) .EQ. 0) THEN 251 IATOMB = IATOMB + 1 252 IF (IATOMB .GT. IATOMA) GO TO 110 253C ... next IATOMA, only IATOMB .le. IATOMA needed 254 CXB=PT(IAND(ISYMAX(1,1),IB))*COORD(1,ICENTB) 255 CYB=PT(IAND(ISYMAX(2,1),IB))*COORD(2,ICENTB) 256 CZB=PT(IAND(ISYMAX(3,1),IB))*COORD(3,ICENTB) 257 VEC(1,IATOMB,IATOMA) = XTANG*(CXA - CXB) 258 VEC(2,IATOMB,IATOMA) = XTANG*(CYA - CYB) 259 VEC(3,IATOMB,IATOMA) = XTANG*(CZA - CZB) 260 VEC(1,IATOMA,IATOMB) = -VEC(1,IATOMB,IATOMA) 261 VEC(2,IATOMA,IATOMB) = -VEC(2,IATOMB,IATOMA) 262 VEC(3,IATOMA,IATOMB) = -VEC(3,IATOMB,IATOMA) 263 END IF 264 210 CONTINUE 265 200 CONTINUE 266 END IF 267 110 CONTINUE 268 IF (ICENTA .EQ. NUCIND) N_QMATOMS = IATOMA 269 100 CONTINUE 270 IF (IATOMA .NE. NUCDEP) THEN 271 WRITE(LUPRI,*) 'GEOANA error, IATOMA .ne. NUCDEP:',IATOMA,NUCDEP 272 CALL QUIT('NCTOT and NUCDEP inconsistent') 273 END IF 274C 275C 276C Set up distance matrix in Angstrom 277C 278 DIST_MAX = D0 279 ABSDIST_MAX= D0 280 IDIST_MAX = 0 281 JDIST_MAX = 0 282 ADIST_MAX = D0 283 IADIST_MAX = 0 284 JADIST_MAX = 0 285 286 N_SHORT_HX = 0 287 N_SHORT_YX = 0 288 ADISTHX_MIN= 1.D200 289 ADISTYX_MIN= 1.D200 290 291 IJ = 0 292 DO 400 I = 1,NUCDEP 293 DO 300 J = 1,I 294 IJ = IJ + 1 295 DISTAN = VEC(1,J,I)*VEC(1,J,I) + VEC(2,J,I)*VEC(2,J,I) 296 * + VEC(3,J,I)*VEC(3,J,I) 297 DISTAN = SQRT(DISTAN) 298 IF (DIF) THEN 299 DISTAN = DISTAN - DIST(IJ) 300 ELSE IF (I .LE. N_QMATOMS .AND. I.NE.J) THEN 301 IF (ICHARG(I) .EQ. 0 .OR. ICHARG(J) .EQ. 0) THEN 302 ! do nothing - we do not want to include floating orbitals in minimum atom distance 303 ELSE IF (ICHARG(I) .NE. 1 .AND. ICHARG(J) .NE. 1) THEN 304 ADISTYX_MIN = MIN(ADISTYX_MIN,DISTAN) 305 IF (DISTAN .LE. 1.0D0) N_SHORT_YX = N_SHORT_YX + 1 ! R(Y-X) .lt. 1.0 Angstrom is usually an error 306 ELSE 307 ADISTHX_MIN = MIN(ADISTHX_MIN,DISTAN) 308 IF (DISTAN .LE. 0.7D0) N_SHORT_HX = N_SHORT_HX + 1 ! R(H-X) .lt. 0.7 Angstrom is usually an error 309 END IF 310 END IF 311 DIST(IJ) = DISTAN 312 IF (ABS(DISTAN) .GT. ABSDIST_MAX) THEN 313 DIST_MAX = DISTAN 314 ABSDIST_MAX = ABS(DIST_MAX) 315 IDIST_MAX = I 316 JDIST_MAX = J 317 IF (I .LE. N_QMATOMS) THEN 318C ... these are QM atoms 319 ADIST_MAX = ABS(DISTAN) 320 IADIST_MAX = I 321 JADIST_MAX = J 322 END IF 323 END IF 324 300 CONTINUE 325 400 CONTINUE 326C 327 IF (PRINT) THEN 328 JPRIDIS = MAX(1,IPRUSR) 329 IF (DIST_MAX .NE. D0) THEN ! exclude atoms (DIST_MAX .eq. 0.D0) 330 331 IF (NUCDEP .LE. 16*JPRIDIS) THEN 332C hjaaj Oct 2003: this output is only useful for small molecules ... 333 IF (DIF) THEN 334 CALL HEADER 335 * ('Differential interatomic separations (in Angstrom):',2) 336 ELSE 337 CALL HEADER('Interatomic separations (in Angstrom):',2) 338 END IF 339 CALL PRIDIS(NAMDEP,DIST,NUCDEP) 340 END IF 341C 342 IF (DIF) THEN 343 WRITE(LUPRI,'(/A,F12.6,A/A,I5,A,I5,5A)') 344 & ' Max differential change in interatomic separation is', 345 & ADIST_MAX,' Angstrom', 346 & ' between atoms',IADIST_MAX,' and',JADIST_MAX, 347 & ', "',NAMDEP(IADIST_MAX),'" and "',NAMDEP(JADIST_MAX),'".' 348 ELSE 349 WRITE(LUPRI,'(/A,2(F10.4,A)/A,I5,A,I5,5A)') 350 & ' Max interatomic separation is', 351 & ADIST_MAX,' Angstrom (',ADIST_MAX/XTANG,' Bohr)', 352 & ' between atoms',IADIST_MAX,' and',JADIST_MAX, 353 & ', "',NAMDEP(IADIST_MAX),'" and "',NAMDEP(JADIST_MAX),'".' 354 IF (ADISTHX_MIN .LT. 1.D10) WRITE(LUPRI,'(/A,2(F10.4,A))') 355 & ' Min HX interatomic separation is', 356 & ADISTHX_MIN,' Angstrom (',ADISTHX_MIN/XTANG,' Bohr)' 357 IF (ADISTYX_MIN .LT. 1.D10) WRITE(LUPRI,'(/A,2(F10.4,A))') 358 & ' Min YX interatomic separation is', 359 & ADISTYX_MIN,' Angstrom (',ADISTYX_MIN/XTANG,' Bohr)' 360 361 IF (N_SHORT_YX .GT. 0 .OR. N_SHORT_HX .GT. 0) THEN 362 NWARN = NWARN + 1 363 WRITE(LUPRI,'(/A,2I5/A/A)') 364 & '@ WARNING: Number of short HX and YX bond lengths:', 365 & N_SHORT_HX,N_SHORT_YX, 366 & '@ WARNING: If not intentional, maybe your coordinates' 367 & //' were in Angstrom,', 368 & '@ WARNING: but "Angstrom" was not specified' 369 & //' in .mol file' 370 END IF 371 372 END IF 373 IF (N_QMATOMS .NE. NUCDEP) THEN 374 IF (DIF) THEN 375 WRITE(LUPRI,'(/A,F12.6,A/A,I5,A,I5,5A)') 376 & ' Max differential change in QM+MM interatomic '// 377 & 'separation is',DIST_MAX,' Angstrom', 378 & ' between the QM+MM centers',IDIST_MAX,' and',JDIST_MAX, 379 & ', "',NAMDEP(IDIST_MAX),'" and "',NAMDEP(JDIST_MAX),'".' 380 ELSE 381 WRITE(LUPRI,'(/A,2(F10.4,A)/A,I5,A,I5,5A)') 382 & ' Max QM+MM interatomic separation is', 383 & DIST_MAX,' Angstrom (',DIST_MAX/XTANG,' Bohr)', 384 & ' between the QM+MM centers',IDIST_MAX,' and',JDIST_MAX, 385 & ', "',NAMDEP(IDIST_MAX),'" and "',NAMDEP(JDIST_MAX),'".' 386 END IF 387 END IF 388 END IF 389 END IF 390C 391 IF (PRINT .AND. .NOT. DIF) THEN 392 393 IJ = 0 394 DO 10 J= 1,N_QMATOMS 395C ... only QM atoms, not MM atoms 396 RADJ = RADIUS(ICHARG(J)) 397 IF (RADJ .LT. 0.0D0 .AND. ICHARG(J) .GT. 0) THEN 398C do not print if cavity center or floating orbital /hjaaj 399 WRITE(LUPRI,*) 400 & 'INFO: RADIUS FOR ATOM WITH ATOMIC NUMBER ', 401 & ICHARG(J),' IS UNAVAILABLE, USING 2.0 AA' 402 RADJ = 2.0D0 403 END IF 404 DO 20 I= 1, J-1 405 IJ = IJ + 1 406 RADI = RADIUS(ICHARG(I)) 407 IF (RADI .LT. 0.0D0 .AND. ICHARG(I) .GT. 0) THEN 408 RADI = 2.0D0 409 END IF 410 IF (RADI .LT. 0 .OR. RADJ .LT. 0) THEN 411C not bonded if cavity center or floating orbital /hjaaj 412 BONDED(I,J) = .FALSE. 413 BONDED(J,I) = .FALSE. 414 ELSE IF (DIST(IJ) .LT. (1.2D0 * (RADI + RADJ))) THEN 415 NBONDS = NBONDS + 1 416 IPAIRS(1,NBONDS) = I 417 IPAIRS(2,NBONDS) = J 418 BONDED(I,J) = .TRUE. 419 BONDED(J,I) = .TRUE. 420 ELSE 421 BONDED(I,J) = .FALSE. 422 BONDED(J,I) = .FALSE. 423 END IF 424 20 CONTINUE 425 IJ = IJ + 1 426 BONDED(J,J) = .FALSE. 427 10 CONTINUE 428C 429 IF (.NOT. QM3) THEN 430 CALL HEADER('Bond distances (Angstrom):',1) 431 WRITE (LUPRI,'(18X,A/18X,A)') 432 & 'atom 1 atom 2 distance', 433 & '------ ------ --------' 434 IJ = 0 435 DO I = 1, N_QMATOMS 436 DO J = 1, I-1 437 IJ = IJ + 1 438 IF (BONDED(I,J)) THEN 439 NUCNAM(1) = NAMDEP(I) 440 NUCNAM(2) = NAMDEP(J) 441 WRITE(LUPRI,'(A,2X,A6,5X,A6,F15.6)') 442 & ' bond distance:', 443 & NUCNAM(1), NUCNAM(2), DIST(IJ) 444 END IF 445 END DO 446 IJ = IJ + 1 447 END DO 448C 449 IF (N_QMATOMS .GT. 2 .AND. NANG .LE. 0) THEN 450 CALL HEADER('Bond angles (degrees):',1) 451 WRITE (LUPRI,'(18X,A/18X,A)') 452 $ 'atom 1 atom 2 atom 3 angle', 453 $ '------ ------ ------ -----' 454C 455 IJK = 0 456 DO 40, I= 1,N_QMATOMS 457 DO 50, J= 1, N_QMATOMS - 1 458 DO 60, K= J + 1, N_QMATOMS 459 IF (BONDED(I,J) .AND. BONDED(I,K)) THEN 460 IJK = IJK + 1 461 NUCNAM(1) = NAMDEP(J) 462 NUCNAM(2) = NAMDEP(I) 463 NUCNAM(3) = NAMDEP(K) 464 ANG = FACANG*VECANG(VEC(1,I,J),VEC(1,I,K)) 465 WRITE(LUPRI,'(A,5X,A6,5X,A6,5X,A6,F14.3)') 466 & ' bond angle:',NUCNAM(1),NUCNAM(2),NUCNAM(3),ANG 467 END IF 468 60 CONTINUE 469 50 CONTINUE 470 40 CONTINUE 471 IF (IJK .EQ. 0) WRITE(LUPRI,'(5X,A)') 'No angles found' 472 END IF 473 END IF ! not QM3 474 END IF 475C 476C Punch bonding information in Gamess output format on unit LUPUNCH 477C 478 IF (LUPUNCH .GT. 0) THEN 479 IF(NBONDS.LE.6) THEN 480 WRITE(LUPUNCH,8010) (IPAIRS(1,I),IPAIRS(2,I),I=1,NBONDS) 481 ELSE 482 WRITE(LUPUNCH,8020) (IPAIRS(1,I),IPAIRS(2,I),I=1,6) 483 WRITE(LUPUNCH,8030) (IPAIRS(1,I),IPAIRS(2,I),I=7,NBONDS) 484 END IF 485 END IF 486C 487 IF (NANG .GT. 0) THEN 488 IF (PRINT) THEN 489 CALL HEADER('Angles according to input list:',2) 490 WRITE (LUPRI,'(A/A)') 491 * ' atom 1 atom 2 atom 3 angle (degrees)', 492 * ' ------ ------ ------ ---------------' 493 END IF 494 DO 1000 I = 1,NANG 495 I1 = IANG(1,I) 496 I2 = IANG(2,I) 497 I3 = IANG(3,I) 498 IMX = MAX(I1,I2,I3) 499 IF (IMX .GT. N_QMATOMS) THEN 500 IF (PRINT) WRITE (LUPRI,'(/A/)') 501 & ' *GEOANA input error for .ANGLES: non-existent atom' 502 GO TO 1000 503 END IF 504 NUCNAM(1) = NAMDEP(I1) 505 NUCNAM(2) = NAMDEP(I2) 506 NUCNAM(3) = NAMDEP(I3) 507 IF (I1 .NE. I2 .AND. I2 .NE. I3) THEN 508 ANG = FACANG*VECANG(VEC(1,I2,I1),VEC(1,I2,I3)) 509 IF (.NOT.DIF) THEN 510 ANGLE(I) = ANG 511 ELSE 512 ANGLE(I) = ANG - ANGLE(I) 513 END IF 514 IF (PRINT) WRITE (LUPRI,'(4X,A6,5X,A6,5X,A6,F20.3)') 515 * NUCNAM(1),NUCNAM(2),NUCNAM(3),ANGLE(I) 516 ELSE 517 IF (PRINT) WRITE (LUPRI,'(4X,A6,5X,A6,5X,A6,10X,A)') 518 * NUCNAM(1),NUCNAM(2),NUCNAM(3),'undefined' 519 END IF 520 1000 CONTINUE 521 END IF 522C 523 IF (NDIHED .GT. 0) THEN 524 IF (PRINT) WRITE (LUPRI,'(//A/A)') 525 * ' atom 1 atom 2 atom 3 atom 4'// 526 * ' dihedral angle (degrees)', 527 * ' ------ ------ ------ ------'// 528 * ' ------------------------' 529 DO 2000 I = 1,NDIHED 530 I1 = IDIHED(1,I) 531 I2 = IDIHED(2,I) 532 I3 = IDIHED(3,I) 533 I4 = IDIHED(4,I) 534 IMX = MAX(I1,I2,I3,I4) 535 IF (IMX .GT. N_QMATOMS) THEN 536 IF (PRINT) WRITE (LUPRI,'(/A/)') 537 & ' *GEOANA input error for .DIHEDR: non-existent atom' 538 GO TO 2000 539 END IF 540 NUCNAM(1) = NAMDEP(I1) 541 NUCNAM(2) = NAMDEP(I2) 542 NUCNAM(3) = NAMDEP(I3) 543 NUCNAM(4) = NAMDEP(I4) 544 X1 = VEC(2,I2,I1)*VEC(3,I2,I3) - VEC(2,I2,I3)*VEC(3,I2,I1) 545 X2 = VEC(3,I2,I1)*VEC(1,I2,I3) - VEC(3,I2,I3)*VEC(1,I2,I1) 546 X3 = VEC(1,I2,I1)*VEC(2,I2,I3) - VEC(1,I2,I3)*VEC(2,I2,I1) 547 Y1 = VEC(2,I3,I2)*VEC(3,I3,I4) - VEC(2,I3,I4)*VEC(3,I3,I2) 548 Y2 = VEC(3,I3,I2)*VEC(1,I3,I4) - VEC(3,I3,I4)*VEC(1,I3,I2) 549 Y3 = VEC(1,I3,I2)*VEC(2,I3,I4) - VEC(1,I3,I4)*VEC(2,I3,I2) 550 Z1 = X2*Y3 - X3*Y2 551 Z2 = X3*Y1 - X1*Y3 552 Z3 = X1*Y2 - X2*Y1 553 SENSE = Z1*VEC(1,I2,I3) + Z2*VEC(2,I2,I3) + Z3*VEC(3,I2,I3) 554 SENSE = SIGN(D1,SENSE) 555 ANG = X1*Y1 + X2*Y2 + X3*Y3 556 DDD = (X1*X1 + X2*X2 + X3*X3) * (Y1*Y1 + Y2*Y2 + Y3*Y3) 557 IF (DDD .GT. 1.D-10) THEN 558 ANG = ANG / SQRT(DDD) 559 IF (ABS(ANG) .GT. D1) ANG = SIGN(D1,ANG) 560 ANG = SENSE*ARCCOS(ANG) 561 IF (.NOT.DIF) THEN 562 DIHED(I) = ANG 563 ELSE 564 DIHED(I) = ANG - DIHED(I) 565 END IF 566 IF (PRINT) WRITE(LUPRI,'(4X,A6,5X,A6,5X,A6,5X,A6,F20.3)') 567 * NUCNAM(1),NUCNAM(2),NUCNAM(3),NUCNAM(4),DIHED(I) 568 ELSE 569 IF (PRINT) WRITE(LUPRI,'(4X,A6,5X,A6,5X,A6,5X,A6,10X,A)') 570 * NUCNAM(1),NUCNAM(2),NUCNAM(3),NUCNAM(4),'undefined' 571 END IF 572 2000 CONTINUE 573 END IF 574C 575 IF (PRINT) WRITE (LUPRI,'(/)') 576 RETURN 577 8010 FORMAT('BONDATOMS ',6(I4,I4,2X)) 578 8020 FORMAT('BONDATOMS ',6(I4,I4,2X),' >') 579 8030 FORMAT(7(I4,I4,2X),:,' >') 580 END 581C /* Deck pridis */ 582 SUBROUTINE PRIDIS (NAMDEP,DISMAT,NROW) 583C 584C 30-Jun-1985 Hans Jorgen Aa. Jensen 585C (based on OUTPAK by Nelson H.F. Beebe) 586C 587C Print bond distance matrix (or other matrix over atoms) 588C 589#include "implicit.h" 590#include "priunit.h" 591 PARAMETER (KCOL=6) 592 CHARACTER*6 NAMDEP(*) 593 DIMENSION DISMAT(*) 594 INTEGER BEGIN 595C 596 LAST = MIN(NROW,KCOL) 597 BEGIN = 1 598 1050 NCOL = 1 599 WRITE (LUPRI,1000) (NAMDEP(I),I = BEGIN,LAST) 600 WRITE (LUPRI,1000) ('------' ,I = BEGIN,LAST) 601 DO 40 K = BEGIN,NROW 602 KTOTAL = (K*(K-1))/2 + BEGIN - 1 603 WRITE (LUPRI,2000) NAMDEP(K), 604 * (DISMAT(KTOTAL+J),J = 1,NCOL) 605 IF (K .LT. (BEGIN+KCOL-1)) NCOL = NCOL + 1 606 40 CONTINUE 607 WRITE (LUPRI,'()') 608 LAST = MIN(LAST+KCOL,NROW) 609 BEGIN = BEGIN+NCOL 610 IF (BEGIN.LE.NROW) GO TO 1050 611 RETURN 612 1000 FORMAT (8X,6(4X,A6,2X)) 613 2000 FORMAT (1X,A6,':',6F12.6) 614 END 615C /* Deck radius */ 616 FUNCTION RADIUS(NCHARGE) 617#include "implicit.h" 618#include "priunit.h" 619C 620C Based on covalent radii and metallic radii in Angstrom. 621C Returns -1 where data is inavailable 622C Oct 2006 hjaaj: changed Hydrogen from 30 to 40 pm, 623C such that H2 is printed as bonded ;-) . 624C 625 DIMENSION RAD(100) 626 DATA (RAD(I), I = 1, 100)/ 627 & 40., 155., 160., 110., 628 & 90., 80., 70., 68., 65., 629 &154., 190., 160., 140., 110., 630 &110., 105., 105., 190., 238., 631 &200., 165., 145., 135., 130., 632 &125., 125., 125., 125., 125., 633 &140., 140., 130., 120., 120., 634 &120., 200., 255., 215., 180., 635 &160., 145., 140., 135., 130., 636 &130., 135., 140., 155., 160., 637 &160., 140., 140., 140., 220., 638 &270., 220., 185., 180., 180., 639 &180., 180., 180., 200., 180., 640 &175., 175., 175., 175., 170., 641 &170., 170., 155., 145., 140., 642 &135., 135., 135., 135., 145., 643 &155., 170., 175., 170., -100., 644 & -100., -100., -100., -100., -100., 645 & -100., -100., -100., -100., -100., 646 & -100., -100., -100., -100., -100., 647 & -100./ 648C 649 IF (NCHARGE .LE. 0) THEN 650Chj = 0: solvent cavity center or floating orbital 651Chj = -Z: multiple basis for r12 methods for nuclear charge Z 652Chj = -1234567890: point charges 653C 654 RADIUS = -1.0D0 655 ELSE IF (NCHARGE .LT. 1 .OR. NCHARGE .GT. 100) THEN 656 WRITE (LUPRI,*) 657 & 'ERROR, RADIUS called with CHARGE =',NCHARGE 658 CALL QUIT('RADIUS called with unvalid CHARGE') 659 ELSE 660 RADIUS = 0.01D0 * RAD(NCHARGE) 661 END IF 662 RETURN 663 END 664C /* Deck vdwrad */ 665 FUNCTION VDWRAD(NCHARGE) 666#include "implicit.h" 667#include "priunit.h" 668C Based on van der Waals radii in Angstrom. 669C Returns -1 where data is inavailable 670 DIMENSION RAD(100) 671 DATA (RAD(I), I = 1, 100)/ 672 & 110., 220., 122., 63., 673 &155., 155., 140., 135., 130., 674 &154., 190., 160., 140., 110., 675 &202., 220., 150., 150., 220., 676 &188., 181., 175., 277., 239., 677 & -100., -100., -100., -100., -100., 678 & -100., -100., -100., -100., -100., 679 & -100., -100., -100., -100., -100., 680 & -100., -100., -100., -100., -100., 681 & -100., -100., -100., -100., -100., 682 & -100., -100., -100., -100., -100., 683 & -100., -100., -100., -100., -100., 684 & -100., -100., -100., -100., -100., 685 & -100., -100., -100., -100., -100., 686 & -100., -100., -100., -100., -100., 687 & -100., -100., -100., -100., -100., 688 & -100., -100., -100., -100., -100., 689 & -100., -100., -100., -100., -100., 690 & -100., -100., -100., -100., -100., 691 & -100., -100., -100., -100., -100., 692 & -100./ 693C 694 IF (NCHARGE .LE. 0) THEN 695Chj = 0: solvent cavity center or floating orbital 696Chj = -Z: multiple basis for r12 methods for nuclear charge Z 697Chj = -1234567890: point charges 698C 699 VDWRAD = -1.0D0 700 ELSE IF (NCHARGE .GT. 100) THEN 701 WRITE (LUPRI,*) 'ERROR, VDWRAD called with CHARGE =',NCHARGE 702 CALL QUIT('VDWRAD called with illegal CHARGE') 703 ELSE 704 VDWRAD = 0.01D0 * RAD(NCHARGE) 705 IF (VDWRAD .LT. 0.0D0) THEN 706C if no table value, use covalent radius plus 0.65 AA 707C (as for B-F above) /Dec. 2006, hjaaj 708 WDWRAD = RADIUS(NCHARGE) 709 IF (VDWRAD .GT. 0.0D0) VDWRAD = VDWRAD + 0.65D0 710 END IF 711 END IF 712 RETURN 713 END 714