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 File: sirgp.F 20C Purpose: general purpose routines for the Sirius module 21C 22C /* Deck errwrk */ 23 SUBROUTINE ERRWRK (STRING,LNEED,LAVAIL) 24C 25C Version 6-Mar-1985 by hjaaj 26C 27 CHARACTER*(*) STRING 28C 29#include "priunit.h" 30C 31 CALL QENTER('ERRWRK') 32 IF (LNEED .GE. 0) THEN 33 WRITE (LUPRI,1010) STRING,LNEED,LAVAIL 34 ELSE 35 WRITE (LUPRI,1020) STRING,-LNEED,LAVAIL 36 END IF 37 CALL QTRACE(LUPRI) 38 CALL QUIT('ERROR, insufficient work space in memory') 39C 40 1010 FORMAT(/' FATAL ERROR, insufficient core for ',A, 41 * /T16,'( Need:',I10,', available:',I10,' )') 42 1020 FORMAT(/' FATAL ERROR, insufficient core for ',A, 43 * //T16,'Need :',I10,' + uncalculated amount', 44 * /T16,'Available :',I10) 45 END 46C /* Deck symort */ 47 SUBROUTINE SYMORT(CMO,SAO,NAO,NMO,WRK,LFREE,IRETUR) 48C 49C 28-Apr-1985 Hans Jorgen Aa. Jensen 50C Revised 9-Aug-1989 hjaaj (exit if diverging) 51C 1-Dec-1995 hjaaj (new exit for numerical problems, 52C see comments) 53C 54C Purpose: 55C Symmetrical orthonormalization of orbitals. 56C 57C Input: 58C CMO(*), non-orthogonal orbitals 59C SAO(*), ao overlap matrix 60C 61C Output: 62C CMO(*), normalized and symmetrically orthogonalized orbitals. 63C 64C Scratch: 65C WRK(LFREE) 66C 67#include "implicit.h" 68 DIMENSION CMO(NAO,*), SAO(*), WRK(*) 69 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, D3 = 3.0D0, DMP5 = -0.5D0 ) 70 PARAMETER ( ITSMAX = 20, TMIN = 1.D-6, CNVFAC = 0.1D0) 71 PARAMETER ( DIFMAX=1.0D-12, DIFMIN=1.0D-4) 72C 73#include "priunit.h" 74#include "infpri.h" 75C 76 CALL QENTER('SYMORT') 77C 78C Symmetric Orthogonalization with N-R iterative algorithm: 79C 80C C(p+1) = 1/2( 3C(p) - C(p) [Ct(p) Sao C(p)]) 81C = 3/2 C(p) - 1/2 C(p) Smo(p) 82C = C(p) [-1/2] [Smo(p) - 3 I] 83C 84 NCMOI = NAO*NMO 85 NNMO = NMO*(NMO+1)/2 86 KSMO = 1 87 KW = KSMO + NNMO 88 LNEED = KW + NAO*NMO - 1 89 IF (LNEED .GT. LFREE) CALL ERRWRK('SYMORT',LNEED,LFREE) 90C 91C Normalize initial vectors 92C (norm must be less than sqrt(3) for proper convergence, 93C if norm .gt. sqrt(5) the algorithm will diverge) 94C 951201: this loop is moved inside iteration loop. 95C 96C 951201/hjaaj: 97C Numerical round-off limit may be reached long before 98C DELMAX .lt. DIFMAX, for example, 99C DIFMAX used to be 1.0D-8, but that caused problems 100C (divergence from e.g. 1.4D-8 to 2.9D-8) for basis sets 101C with eigenvalues of overlap matrix down around 1.0D-6). 102C New exit without error print if DELMAX decreases 103C too slowly when DELMAX .lt. DIFMIN (=1.0D-4). 104C NOTE: because of these round-off problems you should 105C always follow SYMORT with Gram-Schmidt orthonormalizations 106C to come as close as possible to full orthonormalization. 107C By the way, test cases show that iterations (may) converge 108C even though DELMAX increases in initial iterations, but this 109C implementation will exit if DELMAX increases. 110C Test program also indicates that cubic convergence is 111C obtained if CMO vectors are renormalized in each iteration, 112C therefore 10-loop is moved inside iteration loop. 113C 114C Symmetric orthonormalization, using Newton-Raphson iterations 115C 116 PRVMAX = D0 117C ... to avoid compiler messages 118 ITS = 0 119 100 CONTINUE 120 ITS = ITS + 1 121 DO 10 I = 1,NMO 122 CALL MPAPV(NAO,SAO,CMO(1,I),WRK) 123 T = DDOT(NAO,WRK,1,CMO(1,I),1) 124 IF (T .LT. TMIN) THEN 125 IMO = I 126 GO TO 8820 127 END IF 128 T = D1 / SQRT(T) 129 CALL DSCAL(NAO,T,CMO(1,I),1) 130 10 CONTINUE 131C 132 CALL UTHU(SAO,WRK(KSMO),CMO,WRK(KW),NAO,NMO) 133C CALL UTHU(H,HT,U,WRK,NAO,NMO) 134 II = KSMO - 1 135 DO 200 I = 1,NMO 136 II = II + I 137 WRK(II) = WRK(II) - D3 138 200 CONTINUE 139 CALL DSCAL(NNMO,DMP5,WRK(KSMO),1) 140C 141C 142 CALL DZERO(WRK(KW),NCMOI) 143 IJ = KSMO - 1 144 IC1 = KW 145 DO 400 I = 1,NMO 146 JC1 = KW 147 DO 300 J = 1,I-1 148 IJ = IJ + 1 149 FAC = WRK(IJ) 150 CALL DAXPY(NAO,FAC,CMO(1,I),1,WRK(JC1),1) 151 CALL DAXPY(NAO,FAC,CMO(1,J),1,WRK(IC1),1) 152 JC1 = JC1 + NAO 153 300 CONTINUE 154 IJ = IJ + 1 155 FAC = WRK(IJ) 156 CALL DAXPY(NAO,FAC,CMO(1,I),1,WRK(IC1),1) 157 IC1 = IC1 + NAO 158 400 CONTINUE 159C 160 DELMAX = D0 161 IJ = KW - 1 162 DO 550 J = 1,NMO 163 DO 540 I = 1,NAO 164 IJ = IJ + 1 165 DEL = ABS(CMO(I,J)-WRK(IJ)) 166 DELMAX = MAX(DELMAX,DEL) 167 540 CONTINUE 168 550 CONTINUE 169 CALL DCOPY (NCMOI,WRK(KW),1,CMO,1) 170 IF (DELMAX.GT.DIFMAX) THEN 171 IF (ITS .GE. ITSMAX) GO TO 8810 172 IF (ITS .GT. 1) THEN 173 IF (DELMAX .LT. DIFMIN) THEN 174 IF (DELMAX .GT. CNVFAC*PRVMAX) GO TO 8840 175C ... convergence cannot be quadratic in this case, 176C exit because we probably have numeric lin.dep. 177 END IF 178 IF (DELMAX.GT.PRVMAX) GO TO 8830 179 END IF 180 PRVMAX = DELMAX 181 GO TO 100 182 END IF 183C 184 IF (P6FLAG(12)) WRITE (LUPRI,2211) ITS,DELMAX 185 2211 FORMAT(/' (SYMORT) Symmetric orthogonalization in', 186 * I3,' iterations.',/T11,'Max. error: ',F15.10) 187C 188 IRETUR = 0 189 GO TO 9999 190C 191 8810 IRETUR = ITS 192 WRITE (LUPRI,2211) ITS,DELMAX 193 GO TO 9999 194C 195 8820 IRETUR = -IMO 196 WRITE (LUPRI,2311) IMO, T 197 2311 FORMAT(/' (SYMORT) ***ERROR*** norm of input vector no.',I4, 198 * ' is',1P,D10.2) 199 GO TO 9999 200C 201 8830 WRITE (LUPRI,2411) ITS,DELMAX,PRVMAX 202 IRETUR = 9999 203 2411 FORMAT(/' (SYMORT) Symmetric orthogonalization diverging', 204 * /' Iteration no.',I5,' Max. error :',1P,T40,D16.8, 205 * /' Max. error previous iteration :',T40,D16.8) 206 GO TO 9999 207C 208 8840 IF (P6FLAG(12)) WRITE (LUPRI,2511) ITS,DELMAX,PRVMAX 209 IRETUR = 8888 210 2511 FORMAT(/' (SYMORT) Symmetric orthogonalization diverging', 211 & /' because of num. lin.dep. causing round-off errors', 212 * /' Iteration no.',I5,' Max. error :',1P,T40,D16.8, 213 * /' Max. error previous iteration :',T40,D16.8) 214C 215C End of "SYMORT" 216C 217 9999 CALL QEXIT('SYMORT') 218 RETURN 219 END 220C /* Deck vsymtr */ 221 SUBROUTINE VSYMTR(LEN,VEC,THREQL) 222C 223C 5-Jan-1989 Hans Joergen Aa. Jensen 224C Last revision 26-Nov-1992 hjaaj: do not compare zero elements 225C 226C Purpose: symmetrize vector so that numerical 227C inaccuracy won't break symmetry (e.g. in a CI vector) 228C 229C Quite slow because n ** 2 process 230C 231#include "implicit.h" 232 DIMENSION VEC(LEN) 233 PARAMETER (D0 = 0.0D0) 234C 235 CALL QENTER('VSYMTR') 236 IF (THREQL .LE. D0) GO TO 9999 237 DO 200 ICONF = 1,LEN-1 238 CTEST = VEC(ICONF) 239 IF (ABS(CTEST) .LE. THREQL) THEN 240 VEC(ICONF) = D0 241 GO TO 200 242 END IF 243 DO 100 JCONF = ICONF+1, LEN 244 IF (ABS(CTEST-VEC(JCONF)) .LE. THREQL) THEN 245 VEC(JCONF) = CTEST 246 ELSE IF (ABS(CTEST+VEC(JCONF)) .LE. THREQL) THEN 247 VEC(JCONF) = -CTEST 248 END IF 249 100 CONTINUE 250 200 CONTINUE 251C 252 9999 CALL QEXIT('VSYMTR') 253 RETURN 254 END 255C /* Deck upkwop */ 256 SUBROUTINE UPKWOP(NWOPT,JWOP,WOP,UPWOP) 257C 258C Written 1-Nov-1983 by Hans Jorgen Aa. Jensen in Uppsala 259C Revisions 260C 31-Jan-1984 hjaaj 261C 31-Oct-1989 hjaaj : N2ORBX instead of N2ORBT (thus symemtry 262C independent); zero UPWOP 263C 264C This subroutine unpacks ("UPK") an orbital rotation type vector 265C ("WOP") to matrix form. 266C This could for instance be the orbital gradient, or 267C the kappa vector for use in one-index transformation. 268C 269#include "implicit.h" 270 DIMENSION UPWOP(NORBT,NORBT),WOP(NWOPT) 271 INTEGER JWOP(2,NWOPT) 272C 273C Used from common blocks: 274C INFORB : NORBT,N2ORBX 275C 276#include "inforb.h" 277C 278 CALL DZERO(UPWOP,N2ORBX) 279 DO 200 IG = 1,NWOPT 280 I = JWOP(1,IG) 281 J = JWOP(2,IG) 282 UPWOP(I,J) = WOP(IG) 283 UPWOP(J,I) = -WOP(IG) 284 200 CONTINUE 285C 286 RETURN 287 END 288C /* Deck pkwop */ 289 SUBROUTINE PKWOP(NWOPT,JWOP,WOP,UPWOP) 290C 291C Written 30-Jul-1986 by Hans Jorgen Aa. Jensen 292C Revised 31-Oct-1989 hjaaj (removed symmetry blocking of UPWOP) 293C 294C This subroutine packs ("PK") an orbital rotation type vector 295C ("WOP") from matrix form (the reverse of UPKWOP) 296C 297#include "implicit.h" 298 DIMENSION UPWOP(NORBT,NORBT),WOP(NWOPT) 299 INTEGER JWOP(2,NWOPT) 300C 301C Used from common blocks: 302C INFORB : NORBT 303C 304#include "inforb.h" 305C 306 DO 200 IG = 1,NWOPT 307 I = JWOP(1,IG) 308 J = JWOP(2,IG) 309 WOP(IG) = UPWOP(I,J) 310 200 CONTINUE 311C 312 RETURN 313 END 314C /* Deck getwop */ 315 LOGICAL FUNCTION GETWOP(XLABEL) 316C 317C Purpose: read new orbital rotation specification into INFVAR. 318C 319C 9-Nov-1985 hjaaj 320C 1-May-2000 hjaaj: Removed reading of KLWOP (KLWOP not in infvar.h any more) 321C 322#include "implicit.h" 323 CHARACTER*8 XLABEL 324C 325C Used from common blocks: 326C INFORB: MULD2H() 327C INFVAR: MAXOCC,NWOPT,NVAR,NWOPH,NVARH,JWOP(2,*),JWOPSY 328C INFLIN: LSYMRF,LSYMPT,LSYMST,NCONST,NWOPPT,NVARPT 329C INFDIM: NWOPDI,NWOPMA,NCONMA,NVARMA 330C 331#include "maxorb.h" 332#include "priunit.h" 333#include "inforb.h" 334#include "infvar.h" 335#include "inflin.h" 336#include "infdim.h" 337C 338 LOGICAL FNDLAB 339C 340 LUNIT = -1 341 CALL GPOPEN(LUNIT,'LUINDF','UNKNOWN',' ','UNFORMATTED',IDUMMY, 342 & .FALSE.) 343 REWIND LUNIT 344 IF ( FNDLAB(XLABEL,LUNIT) ) THEN 345 GETWOP = .TRUE. 346 READ (LUNIT) NWOPT, NWOPH, JWOPSY 347 CALL READI (LUNIT,(2*NWOPH),JWOP) 348 NWOPDI = MAX(1,NWOPT) 349 NWOPMA = MAX(NWOPMA,NWOPT,NWOPH) 350 NVAR = NCONF + NWOPT 351 NVARH = NCONF + NWOPH 352 NVARMA = NCONMA+ NWOPMA 353 IF (LSYMST .NE. MULD2H(LSYMRF,JWOPSY)) THEN 354 NWARN = NWARN + 1 355 WRITE (LUPRI,'(//A/A/)') 356 & ' GETWOP WARNING: NVARPT not defined because NCONST', 357 & ' corresponds to wrong symmetry.' 358 WRITE (LUPRI,*) 'JWOPSY,LSYMPT',JWOPSY,LSYMPT 359 WRITE (LUPRI,*) 'LSYMRF,LSYMPT,LSYMST',LSYMRF,LSYMPT,LSYMST 360 WRITE (LUPRI,*) 'NCONRF,NCONST,NCONF ',NCONRF,NCONST,NCONF 361 WRITE (LUPRI,*) 'NWOPPT,NWOPT ,NVARPT',NWOPPT,NWOPT ,NVARPT 362 CALL QTRACE(LUPRI) 363 ELSE 364 LSYMPT = JWOPSY 365 NWOPPT = NWOPT 366 NVARPT = NCONST + NWOPPT 367 END IF 368 ELSE 369 GETWOP = .FALSE. 370 END IF 371 CALL GPCLOSE(LUNIT,'KEEP') 372C 373 RETURN 374 END 375C /* Deck make_klwop */ 376 SUBROUTINE MAKE_KLWOP(KLWOP) 377C 378C 1-May-2000 hjaaj 379C 380C Purpose: Generate KLWOP array from JWOP 381C 382#include "implicit.h" 383 DIMENSION KLWOP(NOCCT,NORBT) 384 CHARACTER*8 XLABEL 385C 386C Used from common blocks: 387C INFORB: NOCCT,NORBT 388C INFIND: ISW(*) 389C INFVAR: NWOPT,JWOP(2,*) 390C 391#include "priunit.h" 392#include "maxash.h" 393#include "maxorb.h" 394#include "inforb.h" 395#include "infind.h" 396#include "infvar.h" 397C 398 NKLWOP = NOCCT*NORBT 399 CALL IZERO(KLWOP,NKLWOP) 400C 401 DO IG = 1,NWOPT 402 KW = ISW(JWOP(1,IG)) 403 LW = ISW(JWOP(2,IG))-NISHT 404 IF (KLWOP(KW,LW) .EQ. 0) THEN 405 KLWOP(KW,LW) = IG 406 ELSE 407 WRITE (LUPRI,'(/A,I10,A,I10//A)') 408 & 'Fatal error, duplicate entry in JWOP found for elements', 409 & KLWOP(KW,LW),' and',IG,'Dump of JWOP:' 410 WRITE (LUPRI,'(3I10)') (I,JWOP(1,I),JWOP(2,I),I=1,NWOPT) 411 CALL QUIT 412 & ('MAKE_KLWOP: Fatal error, duplicate entry in JWOP.') 413 END IF 414 END DO 415C 416 RETURN 417 END 418C /* Deck prkap */ 419 SUBROUTINE PRKAP (NKAPT,XKAP,PRFAC,LUPRIN) 420C 421C Jan 90 hjaaj (descendant of PRWOP) 422C 423C Purpose: 424C Write orbital operator array XKAP to unit LUPRIN. 425C 426C if PRFAC .gt. 0, write only elements with absolute value 427C .ge. PRFAC * (max. absolute value) 428C if PRFAC .eq. 0, write all elements; 429C if PRFAC .lt. 0, write only elements with 430C value .le. PRFACX * min. pos value 431C where PRFACX = MAX(-PRFAC, -1.0/PRFAC) 432C 433#include "implicit.h" 434 DIMENSION XKAP(*) 435 PARAMETER (D0 = 0.0D0, D1 = 1.0D0) 436C 437C Used from common blocks: 438C INFVAR : NWOPT,NWOPH,JWOPSY,JWOP(2,*) 439C 440#include "maxorb.h" 441#include "infvar.h" 442C 443C Check if NKAPT legal 444C 445 IF (NKAPT .NE. NWOPT .AND. NKAPT .NE. NWOPH) THEN 446 WRITE (LUPRIN,'(//A/A,3I10)') 447 & ' ERROR in PRKAP, NKAPT .ne. NWOPT,NWOPH', 448 & ' NKAPT, NWOPT, NWOPH :',NKAPT,NWOPT,NWOPH 449 CALL QUIT(' ERROR in PRKAP, NKAPT .ne. NWOPT,NWOPH') 450 END IF 451C 452C Any elements? 453C 454 IF (NKAPT.LE.0) THEN 455 WRITE (LUPRIN,8000) JWOPSY 456 GO TO 9999 457 END IF 458 8000 FORMAT(/5X,'Orbital operator symmetry =',I2, 459 * /5X,'-- NO ELEMENTS --') 460C 461C Print XKAP as requested: 462C 463 IF (PRFAC .GT. D0) THEN 464 KMAX = IDAMAX(NKAPT,XKAP,1) 465 THPKAP = PRFAC*ABS(XKAP(KMAX)) 466 WRITE (LUPRIN,1010) JWOPSY,PRFAC 467 INPR = 0 468 DO 100 IG = 1,NKAPT 469 IF ( ABS(XKAP(IG)) .GT. THPKAP ) THEN 470 K = JWOP(1,IG) 471 L = JWOP(2,IG) 472 WRITE (LUPRIN,1100) IG,K,L,XKAP(IG) 473 ELSE 474 INPR = INPR + 1 475 END IF 476 100 CONTINUE 477 IF (INPR.GT.0) WRITE (LUPRIN,1210) INPR,THPKAP 478 ELSE IF (PRFAC .EQ. D0) THEN 479 WRITE (LUPRIN,1020) JWOPSY 480 INPR = 0 481 DO 200 IG = 1,NKAPT 482 IF ( XKAP(IG) .NE. D0 ) THEN 483 K = JWOP(1,IG) 484 L = JWOP(2,IG) 485 WRITE (LUPRIN,1100) IG,K,L,XKAP(IG) 486 ELSE 487 INPR = INPR + 1 488 END IF 489 200 CONTINUE 490 IF (INPR.GT.0) WRITE (LUPRIN,1220) INPR 491 ELSE 492C (PRFAC .LT. D0) 493 KMIN = IDAMIN(NKAPT,XKAP,1) 494 PRFACX = MAX( -PRFAC, -D1/PRFAC) 495 THPKAP = PRFACX*ABS(XKAP(KMIN)) 496 WRITE (LUPRIN,1030) JWOPSY,THPKAP 497 INPR = 0 498 DO 300 IG = 1,NKAPT 499 IF ( XKAP(IG) .LE. THPKAP ) THEN 500 K = JWOP(1,IG) 501 L = JWOP(2,IG) 502 WRITE (LUPRIN,1100) IG,K,L,XKAP(IG) 503 ELSE 504 INPR = INPR + 1 505 END IF 506 300 CONTINUE 507 IF (INPR.GT.0) WRITE (LUPRIN,1230) INPR,THPKAP 508 END IF 509C 510C 511 9999 CONTINUE 512 RETURN 513C 514 1010 FORMAT(/5X,'Orbital operator symmetry =',I2,/, 515 * 5X,'(only elements abs greater than',1P,D8.1, 516 * ' times max abs value)'//, 517 * 5X,' Index(r,s) r s operator value',/, 518 * 5X,' ---------- -- -- --------------') 519 1020 FORMAT(/5X,'Orbital operator symmetry =',I2,//, 520 * 5X,' Index(r,s) r s operator value',/, 521 * 5X,' ---------- -- -- --------------') 522 1030 FORMAT(/5X,'Orbital operator symmetry =',I2,/, 523 * 5X,'(only elements less than',1P,D10.2,')'// 524 * 5X,' Index(r,s) r s operator value',/, 525 * 5X,' ---------- -- -- --------------') 526 1100 FORMAT(I15,I7,I5,F20.10) 527 1210 FORMAT(/I8,' elements with absolute value less than', 528 * 1P,D9.2,' not printed.') 529 1220 FORMAT(/I8,' zero elements not printed.') 530 1230 FORMAT(/I8,' elements with value greater than', 531 * 1P,D9.2,' not printed.') 532C 533C -- end of subroutine PRKAP 534C 535 END 536C /* Deck prwop */ 537 SUBROUTINE PRWOP (WOP,IUNITX) 538C 539C Written 11-Apr-1984 -- hjaaj 540C Last revision 17-Oct-1984 hjaaj (IUNITX.lt.0 option) 541C 9-May-1984 hjaaj 542C 543C Purpose: 544C Write orbital operator array WOP to unit IUNIT. 545C 546C IUNIT = IABS(IUNITX) 547C 548C If IUNITX .lt. 0, write all elements; 549C if IUNITX .gt. 0, write only elements with absolute value 550C .ge. PRFAC * (max. absolute value) 551C 552#include "implicit.h" 553 DIMENSION WOP(*) 554 PARAMETER (PRFAC = 0.1D0, D0 = 0.0D0) 555C 556C Used from common blocks: 557C INFVAR : NWOPT,JWOPSY,JWOP(2,*) 558C 559#include "maxorb.h" 560#include "infvar.h" 561C 562 IUNIT = IABS(IUNITX) 563C 564C Any elements? 565C 566 IF (NWOPT.LE.0) GO TO 300 567C 568C 1. find maximum absolute value in WOP 569C 570 IF (IUNITX.GT.0) THEN 571 MAX = IDAMAX(NWOPT,WOP,1) 572 THPWOP = PRFAC*ABS(WOP(MAX)) 573 WRITE (IUNIT,1000) JWOPSY,PRFAC 574 ELSE 575 WRITE (IUNIT,1010) JWOPSY 576 THPWOP = D0 577 END IF 578 INPR = 0 579 DO 200 IG = 1,NWOPT 580 IF (ABS(WOP(IG)).GE.THPWOP) THEN 581 K = JWOP(1,IG) 582 L = JWOP(2,IG) 583 WRITE (IUNIT,1100) IG,K,L,WOP(IG) 584 ELSE 585 INPR = INPR + 1 586 END IF 587 200 CONTINUE 588 IF (INPR.GT.0) WRITE (IUNIT,1200) INPR,THPWOP 589C 590 RETURN 591C 592 300 CONTINUE 593 WRITE (IUNIT,3000) JWOPSY 594 RETURN 595C 596 1000 FORMAT(/5X,'Orbital operator symmetry =',I2,/, 597 * 5X,'(only elements abs greater than',1P,D8.1, 598 * ' times max abs value)'//, 599 * 5X,' Index(r,s) r s operator value',/, 600 * 5X,' ---------- -- -- --------------') 601 1010 FORMAT(/5X,'Orbital operator symmetry =',I2,//, 602 * 5X,' Index(r,s) r s operator value',/, 603 * 5X,' ---------- -- -- --------------') 604 1100 FORMAT(I15,I7,I5,F20.10) 605 1200 FORMAT(/I8,' elements with absolute value less than', 606 * 1P,D9.2,' not printed.') 607 3000 FORMAT(/5X,'Orbital operator symmetry =',I2, 608 * /5X,'-- NO ELEMENTS --') 609C 610C -- end of subroutine PRWOP 611C 612 END 613C /* Deck readmo */ 614 SUBROUTINE READMO(CMO,JRDMO) 615C 616C Last revision 7-Jan-1985 hjaaj -- ver. 2.6 617C 8-Jul-1992 hjaaj (removed H1VIRT and FCMO calls) 618C 6-May-1994 hjaaj (removed IRDMO=8 call H1MO) 619C 620C Purpose:Read MO coefficients from IUNIT or generate MO coefficients. 621C IRDMO = abs(JRDMO) determines how they are read. 622C 623C negative JRDMO means NO STOP but return in JRDMO number of errors 624C 625C IRDMO format 626C ----- ------ 627C 3 normal formatted input: (4F18.14) 628C 4 formatted input: (6F12.8) 629C 7 from LUIT1 with label "OLDORB" 630C 9 from LUIT1 with label "NEWORB" 631C 10 from SIRIFC (LUSIFC) 632C 11 from LUCIA program, formatted file LUCIA.CMO 633C 634C 635#include "implicit.h" 636#include "dummy.h" 637 DIMENSION CMO(*) 638 LOGICAL FOUND 639C 640#include "priunit.h" 641#include "inforb.h" 642#include "inftap.h" 643#include "infpri.h" 644C 645 LOGICAL NOSTOP 646 CHARACTER*8 MOLABL, RTNLBL(2) 647 CHARACTER*10 PFMT 648C 649 CALL QENTER('READMO') 650C 651 IRDMO = ABS(JRDMO) 652 IF (JRDMO .LT. 0) THEN 653 JRDMO = 0 654 NOSTOP = .TRUE. 655 ELSE 656 NOSTOP = .FALSE. 657 END IF 658 IF (IRDMO.LT.1 .OR. IRDMO.GT.11) THEN 659 WRITE (LUPRI,10000) IRDMO 660 WRITE (LUERR,10000) IRDMO 66110000 FORMAT(/' *** READMO-ERROR *** IRDMO =',I4,' is out of range') 662 CALL QTRACE(LUPRI) 663 IF (NOSTOP) THEN 664 JRDMO = JRDMO + 1 665 ELSE 666 CALL QUIT('*** ERROR *** illegal option in READMO') 667 END IF 668 END IF 669 GO TO (9000,9000,300,300,8000,8000,700,8000,700,1000, 670 & 1100),IRDMO 671C 672C-- IRDMO=3 or 4: 673C 674 300 IF (IRDMO .EQ. 3) THEN 675 PFMT = '(4F18.14)' 676 ELSE 677 PFMT = '(6F12.8)' 678 END IF 679C 680 JLUCMD = LUCMD 681 IF (JLUCMD .LT. 0) CALL GPOPEN(LUCMD,'DALTON.INP', 682 & 'OLD',' ','FORMATTED',IDUMMY,.FALSE.) 683 REWIND (LUCMD) 684 310 READ (LUCMD,'(A)',END=400,ERR=400) MOLABL 685 IF (MOLABL .NE. '**MOLORB' .AND. MOLABL .NE. '**NATORB') GO TO 310 686 IEND=0 687 DO 340 ISYM=1,NSYM 688 NBASI=NBAS(ISYM) 689 IF(NBASI.EQ.0) GO TO 340 690 NORBI=NORB(ISYM) 691 DO 335 NI=1,NBASI 692 IF (NI.LE.NORBI) THEN 693 IST=IEND+1 694 IEND=IEND+NBASI 695 READ(LUCMD,PFMT) (CMO(I),I=IST,IEND) 696 ELSE 697 READ(LUCMD,PFMT) (DUM,I=1,NBASI) 698 END IF 699 335 CONTINUE 700 340 CONTINUE 701 IF (JLUCMD .LT. 0) CALL GPCLOSE(LUCMD,'KEEP') 702 GO TO 9000 703 400 CONTINUE 704 WRITE (LUPRI,'(//A)') 705 * ' --- FATAL ERROR, label "**MOLORB" not found on SIRIUS input.' 706 CALL QTRACE(LUPRI) 707 IF (NOSTOP) THEN 708 JRDMO = JRDMO + 1 709 ELSE 710 CALL QUIT( 711 & 'ERROR (READMO) label "**MOLORB" not found in input.') 712 END IF 713C 714C-- IRDMO = 7 or 9: 715C 716 700 CONTINUE 717 JLUIT1 = LUIT1 718 IF (JLUIT1 .LE. 0) CALL GPOPEN(LUIT1,'SIRIUS.RST','OLD',' ', 719 & 'UNFORMATTED',IDUMMY,.FALSE.) 720 REWIND LUIT1 721 IF (IRDMO.EQ.7) THEN 722 MOLABL = 'OLDORB ' 723 ELSE 724 MOLABL = 'NEWORB ' 725 END IF 726 LBLERR = -1 727 CALL MOLLB2(MOLABL,RTNLBL,LUIT1,LBLERR) 728 IF (LBLERR .GE. 0) THEN 729 CALL READT(LUIT1,NCMOT,CMO) 730 IF (IPRSIR .GT. 0) WRITE(LUPRI,'(4A,1X,A)') 731 & '* Read MO coefficients from SIRIUS.RST with label "',MOLABL, 732 & '" and header info: ',RTNLBL(1:2) 733 ELSE IF (LBLERR .EQ. -1) THEN 734 WRITE (LUPRI,'(//3A)') 735 & ' --- FATAL ERROR, MO label "',MOLABL,'" not found on SIRIUS.RST' 736 CALL QTRACE(LUPRI) 737 IF (NOSTOP) THEN 738 JRDMO = JRDMO + 1 739 ELSE 740 CALL QUIT('ERROR (READMO) MO coefficients not found.') 741 END IF 742 ELSE 743 WRITE (LUPRI,'(//A/2A)') 744 & ' --- FATAL ERROR, could not read SIRIUS.RST', 745 & ' when searching for mo label :',MOLABL 746 CALL QTRACE(LUPRI) 747 IF (NOSTOP) THEN 748 JRDMO = JRDMO + 1 749 ELSE 750 CALL QUIT('ERROR (READMO) COULD NOT READ MO FILE.') 751 END IF 752 END IF 753 IF (JLUIT1 .LE. 0) CALL GPCLOSE(LUIT1,'KEEP') 754 GO TO 9000 755C 756C-- IRDMO = 10: read CMO from SIRIFC (LUSIFC) 757C 758 1000 CONTINUE 759 CALL RD_SIRIFC('CMO',FOUND,CMO) 760 IF (.NOT.FOUND) THEN 761 IF (NOSTOP) THEN 762 JRDMO = JRDMO + 1 763 ELSE 764 CALL QUIT('ERROR (READMO) MO''s not found on SIRIFC') 765 END IF 766 END IF 767 GO TO 9000 768C 769C-- IRDMO = 11: read CMO from LUCIA program, formatted file LUCIA.CMO 770C 771 1100 CONTINUE 772 LU_LUCIA = -1 773 JLUIT1 = LU_LUCIA 774 IF (JLUIT1 .LE. 0) CALL GPOPEN(LU_LUCIA,'LUCIA_91','OLD',' ', 775 & 'FORMATTED',IDUMMY,.FALSE.) 776 REWIND LU_LUCIA 777 778 call rd_lucia(cmo,nbas,norb,nsym,lu_lucia,found) 779 780 IF (JLUIT1 .LE. 0) CALL GPCLOSE(LU_LUCIA,'KEEP') 781 782 GO TO 9000 783C 784C-- IRDMO 5,6,8 are not defined any more. 785C 786 8000 CONTINUE 787C 788 WRITE (LUPRI,'(//A/A)') 789 * ' --- FATAL INTERNAL ERROR, IRDMO options 5, 6, 8', 790 * ' is not defined in READMO since May-95' 791 IF (NOSTOP) THEN 792 JRDMO = JRDMO + 1 793 ELSE 794 CALL QUIT('ERROR (READMO) options 5,6,8 is not defined') 795 END IF 796C 797C-- 798C 799 9000 CONTINUE 800C 801 CALL QEXIT('READMO') 802 RETURN 803 END 804C /* Deck neworb */ 805 SUBROUTINE NEWORB(LABEL3,CMO,REWIT1) 806C 807C 19-Oct-1989; 3-Aug-1990 (REWIT1; hjaaj) 808C 809C If (REWIT1) then 810C Write molecular orbitals with label 'NEWORB ' at 811C beginning of LUIT1 812C else 813C Overwrite 'NEWORB ' on LUIT1 814C end if 815C 816#include "implicit.h" 817 CHARACTER*8 LABEL3 818 DIMENSION CMO(NCMOT) 819 LOGICAL REWIT1 820C 821#include "dummy.h" 822C 823C Used from common blocks: 824C INFORB : NCMOT 825C INFTAP : LUIT1 826C exeinf.h : NEWCMO 827C 828#include "inforb.h" 829#include "inftap.h" 830#include "exeinf.h" 831C 832 LOGICAL FNDLAB 833 CHARACTER*8 MOLBL(3),LABTBL(3) 834 DATA MOLBL /'********',' ',' '/ 835 DATA LABTBL/'BASINFO ','NEWORB ','EODATA '/ 836C 837 CALL QENTER('NEWORB') 838C 839C Write orbitals to LUIT1 840C 841 CALL GETDAT(MOLBL(2),MOLBL(3)) 842 NCMOT4 = MAX(4,NCMOT) 843 IF (REWIT1) THEN 844 CALL NEWIT1 845 ELSE 846 REWIND LUIT1 847 IF (.NOT.FNDLAB(LABTBL(1),LUIT1)) THEN 848C no "BASINFO ", 849C must be a new SIRIUS.RST file, save BASINFO on it: 850 CALL NEWIT1 851 GO TO 100 852 END IF 853 REWIND LUIT1 854 IF (FNDLAB(LABTBL(2),LUIT1)) THEN 855C found "NEWORB ", 856C backspace before writing new "NEWORB " label 857 BACKSPACE LUIT1 858 GO TO 100 859 END IF 860 REWIND LUIT1 861 IF (FNDLAB(LABTBL(3),LUIT1)) THEN 862C no "NEWORB " but found "EODATA ", 863C backspace before writing "NEWORB " label 864 BACKSPACE LUIT1 865 GO TO 100 866 END IF 867C 868C no "NEWORB " and no "EODATA " 869C is positioned at end-of-file on LUIT1 after last FNDLAB 870C /hjaaj June 09 871 END IF 872 873 100 WRITE (LUIT1) MOLBL(1),MOLBL(2),LABEL3,LABTBL(2) 874C ...label 'NEWORB ' 875 CALL WRITT(LUIT1,NCMOT4,CMO) 876 WRITE (LUIT1) MOLBL,LABTBL(3) 877C ...label 'EODATA ' 878 REWIND LUIT1 879 NEWCMO = .TRUE. 880C 881 CALL QEXIT('NEWORB') 882 RETURN 883 END 884C /* Deck newit1 */ 885 SUBROUTINE NEWIT1 886C 887C 2-Dec-1992 Hans Joergen Aa. Jensen 888C 889C Write basis set information at beginning of LUIT1 890C 950825-hjaaj: also save NRHF(*) and IOPRHF for use with AUTOCC 891C 892#include "implicit.h" 893C 894C Used from common blocks: 895C INFORB: NSYM, NORB, NBAS, NRHF 896C SCBRHF: IOPRHF 897C INFTAP: LUIT1 898C 899#include "inforb.h" 900#include "scbrhf.h" 901#include "inftap.h" 902C 903C 904 INTEGER*4 NSYM_i4, NBAS_i4(8),NORB_i4(8), NRHF_i4(8), IOPRHF_i4 905 ! hjaaj Apr 2011: such that Dalton with 8 byte integers can read a 906 ! SIRIUS.RST from a Dalton with 4 byte integers, and vice versa. 907C 908 CHARACTER*8 BASINFO(4) 909 DATA BASINFO(1)/'********'/, BASINFO(4)/'BASINFO '/ 910C 911 CALL GETDAT(BASINFO(2),BASINFO(3)) 912 REWIND LUIT1 913 WRITE (LUIT1) BASINFO 914 NSYM_i4 = NSYM 915 NBAS_i4(:) = NBAS(:) 916 NORB_i4(:) = NORB(:) 917 NRHF_i4(:) = NRHF(:) 918 IOPRHF_i4 = IOPRHF 919 WRITE (LUIT1) NSYM_i4,NBAS_i4,NORB_i4,NRHF_i4,IOPRHF_i4 920 RETURN 921 END 922C /* Deck getac */ 923 SUBROUTINE GETAC(FC,FCAC) 924C 925C 7-Feb-1987 Hans Joergen Aa. Jensen 926C 927C Module : SIRGP 928C 929C Purpose: Extract block with active-active orbital indices out 930C of symmetry packed matrix (e.g. the Fock matrix). 931C 932#include "implicit.h" 933 DIMENSION FC(*), FCAC(*) 934C 935C INFORB : NASHT, NNASHX, NSYM, NASH(*), IIORB(*), 936C NISH(*), NOCC(*) 937C INFIND : IROW(*), ICH(*) 938C 939#include "maxash.h" 940#include "maxorb.h" 941#include "inforb.h" 942#include "infind.h" 943C 944 CALL QENTER('GETAC ') 945 IF (NASHT .EQ. 0) GO TO 9999 946C 947 CALL DZERO(FCAC,NNASHX) 948 DO 300 ISYM=1,NSYM 949 IF (NASH(ISYM).EQ.0) GO TO 300 950 NISHI = NISH(ISYM) 951 NOCI = NOCC(ISYM) 952 IJFC = IIORB(ISYM) + IROW(NISHI+1) 953 IORBI = IORB(ISYM) 954 DO 290 I = (NISHI+1),NOCI 955 NI = ICH(IORBI+I) 956 IROWNI = IROW(NI) 957 IJFC = IJFC + NISHI 958 DO 280 J = (NISHI+1),I 959 IJFC = IJFC + 1 960 NJ = ICH(IORBI+J) 961 IF (NI.GE.NJ) THEN 962 NIJ = IROWNI + NJ 963 ELSE 964 NIJ = IROW(NJ) + NI 965 END IF 966 FCAC(NIJ) = FC(IJFC) 967 280 CONTINUE 968 290 CONTINUE 969 300 CONTINUE 970C 971 9999 CALL QEXIT('GETAC ') 972 RETURN 973C 974C End of GETAC. 975C 976 END 977C /* Deck getac1 */ 978 SUBROUTINE GETAC1(HH,HHAC) 979C 980C 2-Nov-1989 Hans Joergen Aa. Jensen 981C 982C Module : SIRGP 983C 984C Purpose: Extract block with active-active orbital indices out 985C of full matrix. 986C 987#include "implicit.h" 988 DIMENSION HH(NORBT,NORBT), HHAC(NASHT,NASHT) 989C 990C INFORB : NORBT, NASHT, N2ASHX 991C INFIND : IOBTYP(*), ICH(*) 992C 993#include "maxash.h" 994#include "maxorb.h" 995#include "inforb.h" 996#include "infind.h" 997C 998 CALL QENTER('GETAC1') 999 IF (NASHT .EQ. 0) GO TO 9999 1000C 1001 DO 200 J = 1,NORBT 1002 IF (IOBTYP(J) .EQ. JTACT) THEN 1003 NJ = ICH(J) 1004 DO 100 I = 1,NORBT 1005 IF (IOBTYP(I) .EQ. JTACT) THEN 1006 NI = ICH(I) 1007 HHAC(NI,NJ) = HH(I,J) 1008 END IF 1009 100 CONTINUE 1010 END IF 1011 200 CONTINUE 1012C 1013 9999 CALL QEXIT('GETAC1') 1014 RETURN 1015C 1016C End of GETAC1. 1017C 1018 END 1019C /* Deck getac2 */ 1020 SUBROUTINE GETAC2(HH,HHAC) 1021C 1022C 6-May-1990 Hans Joergen Aa. Jensen 1023C 1024C Module : SIRGP 1025C 1026C Purpose: Extract NNASHX block with active-active orbital indices out 1027C of NNORBX matrix. 1028C 1029#include "implicit.h" 1030 DIMENSION HH(NNORBX), HHAC(NNASHX) 1031C 1032C INFORB : NORBT, NASHT, NNORBX,NNASHX 1033C INFIND : IOBTYP(*), ICH(*) 1034C 1035#include "maxash.h" 1036#include "maxorb.h" 1037#include "inforb.h" 1038#include "infind.h" 1039C 1040 CALL QENTER('GETAC2') 1041 IF (NASHT .EQ. 0) GO TO 9999 1042C 1043 DO 200 J = 1,NORBT 1044 IF (IOBTYP(J) .EQ. JTACT) THEN 1045 NJ = ICH(J) 1046 IRNJ= IROW(NJ) 1047 IRJ = IROW(J) 1048 DO 100 I = 1,J 1049 IF (IOBTYP(I) .EQ. JTACT) THEN 1050 NI = ICH(I) 1051 HHAC(IRNJ+NI) = HH(IRJ+I) 1052 END IF 1053 100 CONTINUE 1054 END IF 1055 200 CONTINUE 1056C 1057 9999 CALL QEXIT('GETAC2') 1058 RETURN 1059C 1060C End of GETAC2. 1061C 1062 END 1063C /* Deck getvir */ 1064 SUBROUTINE GETVIR(FC,FCVIR) 1065C 1066C 7-Nov-2017 Hans Joergen Aa. Jensen (based on GETAC) 1067C 1068C Module : SIRGP 1069C 1070C Purpose: Extract blocks with virtual-virtual orbital indices out 1071C of symmetry packed matrix (e.g. the Fock matrix). 1072C Both FC and FCVIR are symmetry packed. 1073C 1074#include "implicit.h" 1075 REAL*8 FC(*), FCVIR(*) 1076C 1077C INFORB : NSSHT, NNSSHT, NSYM, NSSH(*), IIORB(*), IISH(*), NOCC(*) 1078C 1079#include "inforb.h" 1080 1081 IROW(I) = (I*(I-1))/2 1082C 1083 CALL QENTER('GETVIR ') 1084C 1085 CALL DZERO(FCVIR,NNSSHT) 1086 IISSH_I = 0 1087 DO 300 ISYM=1,NSYM 1088 IF (NSSH(ISYM).EQ.0) GO TO 300 1089 NOCCI = NOCC(ISYM) 1090 NORBI = NORB(ISYM) 1091 IJFC = IIORB(ISYM) + IROW(NOCCI+1) 1092 IJFCVIR = IISSH_I 1093 DO 290 I = (NOCCI+1),NORBI 1094 IJFC = IJFC + NOCCI 1095 DO 280 J = (NOCCI+1),I 1096 IJFC = IJFC + 1 1097 IJFCVIR = IJFCVIR + 1 1098 FCVIR(IJFCVIR) = FC(IJFC) 1099 280 CONTINUE 1100 290 CONTINUE 1101 IISSH_I = IISSH_I + IROW(NSSH(ISYM)+1) 1102 300 CONTINUE 1103C 1104 CALL QEXIT('GETVIR ') 1105 RETURN 1106C 1107C End of GETVIR. 1108C 1109 END 1110C /* Deck pksym1 */ 1111 SUBROUTINE PKSYM1(ASP,APK,NDIM,NBLK,IWAY) 1112C 1113C Jan 90/May 96 Hans Joergen Aa. Jensen 1114C 1115C If IWAY .ge. 0 then 1116C (in this case ASP and APK may overlap as long as 1117C loc(APK) .le. loc(ASP)) 1118C If IWAY .eq. 2 then 1119C Pack NNDIMX matrix ASP of symmetry 1 to ASP in NNDIMT format 1120C i.e. call pksym1(asp,asp,ndim,nblk,+2) 1121C else 1122C Pack NNDIMX matrix ASP of symmetry 1 to APK in NNDIMT format 1123C end if 1124C else 1125C Unpack from APK to ASP 1126C 1127#include "implicit.h" 1128 DIMENSION ASP(*), APK(*), NDIM(NBLK) 1129 IF (IWAY .LT. 0) THEN 1130 NDIMT = ISUM(NBLK,NDIM,1) 1131 NNDIMX = (NDIMT*NDIMT+NDIMT)/2 1132 CALL DZERO(ASP,NNDIMX) 1133 END IF 1134 IF (IWAY .EQ. 2) THEN 1135C block 1 is not moved if equivalence (asp,apk) 1136 IBLK1 = 2 1137 IDIMI = NDIM(1) 1138 ISPOFF = IDIMI*(IDIMI+1)/2 1139 IPKOFF = ISPOFF 1140 ELSE 1141 IBLK1 = 1 1142 IDIMI = 0 1143 ISPOFF = 0 1144 IPKOFF = 0 1145 END IF 1146 DO 800 IBLK = IBLK1,NBLK 1147 NDIMI = NDIM(IBLK) 1148 DO 600 J = 1,NDIMI 1149 ISPOFF = ISPOFF + IDIMI 1150 IF (IWAY .GE. 0) THEN 1151C no synchronization problems as long as 1152C loc(apk) .le. loc(asp) (or no overlap) 1153 DO I = 1,J 1154 APK(IPKOFF+I) = ASP(ISPOFF+I) 1155 END DO 1156 ELSE 1157 DO I = 1,J 1158 ASP(ISPOFF+I) = APK(IPKOFF+I) 1159 END DO 1160 END IF 1161 ISPOFF = ISPOFF + J 1162 IPKOFF = IPKOFF + J 1163 600 CONTINUE 1164 IDIMI = IDIMI + NDIMI 1165 800 CONTINUE 1166 RETURN 1167 END 1168C /* Deck rdonel */ 1169 SUBROUTINE RDONEL(KEY,FOUND,A,NA) 1170C 1171C 31-Jul-1987 hjaaj 1172C 1173C Purpose : Read LUONEL 1174C 1175C Note: if FOUND true, it must not be updated. 1176C 1177C 1178#include "implicit.h" 1179#include "priunit.h" 1180C Global array 1181 DIMENSION A(NA,*) 1182 CHARACTER KEY*(*) 1183 LOGICAL FOUND, RDBUF 1184C 1185C Used from common blocks: 1186C INFINP : TITMOL(2), POTNUC, CENT(*), TYPE(*) 1187C INFORB : NBAS(8), NSYM, NBAST, NNBAST 1188C 1189#include "maxorb.h" 1190#include "infinp.h" 1191#include "inforb.h" 1192#include "inftap.h" 1193#include "r12int.h" 1194C Local arrays 1195 PARAMETER ( LBUF = 600 ) 1196 DIMENSION BUF(LBUF), DBUF(LBUF,3), QBUF(LBUF,6), IBUF(LBUF) 1197 EQUIVALENCE (BUF,DBUF,QBUF) 1198 LOGICAL OPEND, FEXIST, ERRSTP 1199 INTEGER, allocatable :: IJ(:), ITMPMO(:) 1200C 1201 CALL QENTER('RDONEL') 1202 ERRSTP = FOUND 1203 IF (.NOT.FOUND) FOUND = .TRUE. 1204C Note: if FOUND true, it must not be updated. 1205 RDBUF = .FALSE. 1206C 1207 CALL GPINQ(FNONEL,'OPENED',OPEND) 1208 IF (.NOT.OPEND) THEN 1209 CALL GPINQ(FNONEL,'EXIST',FEXIST) 1210 IF (FEXIST) THEN 1211 CALL GPOPEN(LUONEL,FNONEL,'OLD',' ','UNFORMATTED',IDUMMY, 1212 & .FALSE.) 1213 OPEND = .TRUE. 1214 IF (LBONEL .LT. 0) LBONEL = LBUF 1215 ELSE 1216 WRITE (LUPRI,'(/3A)') 1217 * ' RDONEL ERROR: LUONEL file ',FNONEL,' does not exist' 1218 WRITE (LUERR ,'(/3A)') 1219 * ' RDONEL ERROR: LUONEL file ',FNONEL,' does not exist' 1220 IF (ERRSTP) THEN 1221 CALL QUIT('RDONEL ERROR: AOONEINT file DOES NOT EXIST') 1222 ELSE 1223 FOUND = .FALSE. 1224 GO TO 9999 1225 END IF 1226 END IF 1227 END IF 1228C 1229 IF (NOAUXB) allocate( IJ(NNBASX) ) 1230C 1231 IF (KEY.EQ.'OPEN ') THEN 1232C ... file is opened automatically first time RDONEL is called. 1233 ELSE IF(KEY.EQ.'CLOSE ') THEN 1234 CALL GPCLOSE(LUONEL,'KEEP') 1235 ELSE IF(KEY.EQ.'MLCLINFO') THEN 1236 REWIND LUONEL 1237 READ (LUONEL) TITMOL 1238 READ (LUONEL) NSYM,(NBAS(ISYM),ISYM=1,NSYM),POTNUC 1239 NBAST = 0 1240 NNBAST = 0 1241 DO 110 I = 1,NSYM 1242 NBAST = NBAST + NBAS(I) 1243 NNBAST = NNBAST + NBAS(I)*(NBAS(I)+1)/2 1244 110 CONTINUE 1245 DO 120 I = NSYM+1,8 1246 NBAS(I) = 0 1247 120 CONTINUE 1248 allocate( ITMPMO(NBAST*2) ) 1249 CALL MOLLAB('SYMINPUT',LUONEL,LUPRI) 1250 READ (LUONEL) NBT,(ITMPMO(I),I=1,2*NBT) 1251C READ (LUONEL) NBT,(ICENT(I),I=1,NBT),(ITYPE(I),I=1,NBT) 1252 IF (NBT .NE. NBAST) THEN 1253 WRITE (LUPRI,'(//A,2(/A,I5))') 1254 * ' --- RDONEL ERROR', 1255 * ' NBAST from first record on AOONEINT :',NBAST, 1256 * ' NBAST from "SYMINPUT" on AOONEINT :',NBT 1257 CALL QUIT('RDONEL ERROR, two different NBAST on AOONEINT') 1258 END IF 1259 REWIND LUONEL 1260 DO I = 1,NBAST 1261 WRITE(CENT(I),'(A4)') ITMPMO(I) 1262 WRITE(TYPE(I),'(A4)') ITMPMO(NBAST+I) 1263 END DO 1264 deallocate( ITMPMO ) 1265 ELSE IF(KEY.EQ.'OVERLAP ') THEN 1266C Read overlap, in symmetry blocked triangular form 1267 REWIND LUONEL 1268 CALL MOLLAB('OVERLAP ',LUONEL,LUPRI) 1269 RDBUF = .TRUE. 1270 ELSE IF(KEY.EQ.'ONEHAMIL') THEN 1271C Read one electron part, in symmetry blocked triangular form 1272 REWIND LUONEL 1273 CALL MOLLAB('ONEHAMIL',LUONEL,LUPRI) 1274 RDBUF = .TRUE. 1275 ELSE IF(KEY.EQ.'KINETINT') THEN 1276 REWIND LUONEL 1277 CALL MOLLAB('KINETINT',LUONEL,LUPRI) 1278 RDBUF = .TRUE. 1279 ELSE IF(KEY.EQ.'DIPOLMOM') THEN 1280 REWIND LUONEL 1281 CALL MOLLAB('DIPOLMOM',LUONEL,LUPRI) 1282 IF (LBONEL .NE. LBUF) THEN 1283 WRITE(LUPRI,'(/A,2(/A,I8))') 1284 * ' ERROR RDONEL, LBONEL .ne. local LBUF for DIPOLMOM', 1285 * ' LBONEL =',LBONEL,' LBUF =',LBUF 1286 CALL QUIT('RDONEL ERROR, LBONEL .ne. LBUF for key DIPOLMOM') 1287 END IF 1288 IF (NOAUXB) THEN 1289 CALL SET_INOAUX(IJ) 1290 END IF 1291 501 READ (LUONEL) DBUF,IBUF,NUT 1292 IF (NUT.EQ.0) GO TO 501 1293 IF (NUT.LT.0) GO TO 502 1294 IF (NOAUXB) THEN 1295 DO I=1,NUT 1296 IR12 = IJ( IBUF(I) ) 1297 IF (IR12 .GT. 0) THEN 1298 A(IR12,1) = DBUF(I,1) 1299 A(IR12,2) = DBUF(I,2) 1300 A(IR12,3) = DBUF(I,3) 1301 END IF 1302 END DO 1303 ELSE 1304 DO I=1,NUT 1305 J = IBUF(I) 1306 A(J,1) = DBUF(I,1) 1307 A(J,2) = DBUF(I,2) 1308 A(J,3) = DBUF(I,3) 1309 END DO 1310 END IF 1311 GO TO 501 1312 502 CONTINUE 1313 ELSE IF(KEY.EQ.'QUADRP ') THEN 1314 REWIND LUONEL 1315 CALL MOLLAB('QUADRP ',LUONEL,LUPRI) 1316 IF (LBONEL .NE. LBUF) THEN 1317 WRITE(LUPRI,'(/A,2(/A,I8))') 1318 * ' ERROR RDONEL, LBONEL .ne. local LBUF for QUADRP ', 1319 * ' LBONEL =',LBONEL,' LBUF =',LBUF 1320 CALL QUIT('RDONEL ERROR, LBONEL .ne. LBUF for key QUADRP') 1321 END IF 1322 IF (NOAUXB) THEN 1323 CALL SET_INOAUX(IJ) 1324 END IF 1325 601 READ (LUONEL) QBUF,IBUF,NUT 1326 IF (NUT.EQ.0) GO TO 601 1327 IF (NUT.LT.0) GO TO 602 1328 IF (NOAUXB) THEN 1329 DO I=1,NUT 1330 IR12 = IJ( IBUF(I) ) 1331 IF (IR12 .GT. 0) THEN 1332 DO ICOMP = 1,6 1333 A(IR12,ICOMP) = QBUF(I,ICOMP) 1334 END DO 1335 END IF 1336 END DO 1337 ELSE 1338 DO I=1,NUT 1339 J = IBUF(I) 1340 DO ICOMP = 1,6 1341 A(J,ICOMP) = QBUF(I,ICOMP) 1342 END DO 1343 END DO 1344 END IF 1345 GO TO 601 1346 602 CONTINUE 1347 ELSE IF(KEY.EQ.'ISORDK ') THEN 1348 WRITE(LUPRI,'(2A)') KEY, ' IS NOT IMPLEMENTED YET IN RDONEL' 1349 IF (ERRSTP) THEN 1350 CALL QUIT('RDONEL ERROR, KEY ISORDK NOT IMPLEMENTED YET.') 1351 ELSE 1352 FOUND = .FALSE. 1353 END IF 1354 ELSE IF(KEY.EQ.'SCFINP ') THEN 1355 WRITE(LUPRI,'(2A)') KEY, ' IS NOT IMPLEMENTED YET IN RDONEL' 1356 IF (ERRSTP) THEN 1357 CALL QUIT('RDONEL ERROR, KEY SCFINP NOT IMPLEMENTED YET.') 1358 ELSE 1359 FOUND = .FALSE. 1360 END IF 1361 ELSE 1362 WRITE(LUPRI,'(2A)') KEY, ' IS AN INVALID KEY TO RDONEL' 1363 IF (ERRSTP) THEN 1364 CALL QUIT('ERROR RDONEL, INVALID KEY') 1365 ELSE 1366 FOUND = .FALSE. 1367 END IF 1368 END IF 1369 IF (RDBUF) THEN 1370 IF (LBONEL .GT. LBUF) THEN 1371 WRITE(LUPRI,'(/A,2(/A,I8))') 1372 * ' ERROR RDONEL, LBONEL .GT. local LBUF', 1373 * ' LBONEL =',LBONEL,' LBUF =',LBUF 1374 CALL QUIT('ERROR RDONEL, LBONEL .gt. LBUF') 1375 END IF 1376 IF (NOAUXB) THEN 1377 CALL SET_INOAUX(IJ) 1378 END IF 1379 CALL DZERO(A,NNBAST) 1380 2100 READ (LUONEL) (BUF(I),I=1,LBONEL),(IBUF(I),I=1,LBONEL),LENGTH 1381 IF (NOAUXB) THEN 1382 DO I = 1, LENGTH 1383 IR12 = IJ( IBUF(I) ) 1384 IF (IR12 .GT. 0) A( IR12 , 1 ) = BUF(I) 1385 END DO 1386 ELSE 1387 DO I = 1, LENGTH 1388 A( IBUF(I) , 1 ) = BUF(I) 1389 END DO 1390 END IF 1391 IF (LENGTH .GE. 0) GO TO 2100 1392 REWIND LUONEL 1393 END IF 1394 CALL GPCLOSE(LUONEL,'KEEP') 1395C 1396 if ( allocated(IJ) ) deallocate( IJ ) 1397 9999 CALL QEXIT('RDONEL') 1398 RETURN 1399 END 1400 subroutine rd_lucia(cmo,nbas,norb,nsym,lu,found) 1401! 1402! read LUCIA orbital file (fort.91) - v2012 - Stefan Knecht 1403! ------------------------------------------------------------------ 1404 implicit none 1405! ------------------------------------------------------------------ 1406 real(8), intent(inout) :: cmo(*) 1407 integer, intent(in) :: lu 1408 integer, intent(in) :: nsym 1409 integer, intent(in) :: nbas(nsym) 1410 integer, intent(in) :: norb(nsym) 1411 logical, intent(inout) :: found 1412! ------------------------------------------------------------------ 1413 integer :: nsym_lu 1414 integer :: nbas_lu(nsym) 1415 integer :: norb_lu(nsym) 1416 integer :: ncmoao 1417 integer :: ncmoao_lu 1418 integer :: nao_tot 1419 integer :: i 1420 character (len=4), allocatable :: ao_cent(:) 1421 character (len=4), allocatable :: ao_type(:) 1422! ------------------------------------------------------------------ 1423 1424 found = .false. 1425 nsym_lu = -1 1426 call izero(nbas_lu,nsym) 1427 call izero(norb_lu,nsym) 1428 1429! symmetry info 1430 read(lu,*) nsym_lu 1431 if(nsym_lu /= nsym) call quit('*** error in rd_lucia: nsym!') 1432 read(lu,*) (norb_lu(i),i=1, nsym_lu) 1433 read(lu,*) (nbas_lu(i),i=1, nsym_lu) 1434 read(lu,*) ncmoao_lu 1435 1436! MO coefficients 1437 ncmoao = 0 1438 do i = 1, nsym 1439 ncmoao = ncmoao + norb(i) * nbas(i) 1440 end do 1441 if(ncmoao /= ncmoao_lu)call quit('*** error in rd_lucia: ncmoao!') 1442 1443 read(lu,*) (cmo(i), i=1, ncmoao_lu) 1444 1445! AO labels 1446 nao_tot = 0 1447 do i = 1, nsym_lu 1448 nao_tot = nao_tot + nbas_lu(i) 1449 end do 1450 allocate(ao_cent(nao_tot*4)) 1451 allocate(ao_type(nao_tot*4)) 1452 read(lu,'(20A4)') (ao_cent(i),i = 1, nao_tot) 1453 read(lu,'(20A4)') (ao_type(i),i = 1, nao_tot) 1454 deallocate(ao_type) 1455 deallocate(ao_cent) 1456 1457 found = .true. 1458 1459 END 1460 1461 subroutine write_aodens(dv,nisht,nasht,nbast,nnashx,ncmot,nsym, 1462 & nbas,ibas,lupri,hsrohf,density) 1463! 1464! !> purpose: write MCSCF and HF densities in AO basis to disk (and print them to screen) 1465 1466! !> input: active 1p-density matrix in MO basis (from MCSCF) 1467 1468! !> notes: 1469! !> the active part of the density matrix is stored in dvao 1470! !> the in-active part of the density matrix is stored in dcao 1471! !> in case of open-shell HF (hsrohf == .true.) we add the active part of the density matrix to the inactive part and save the final result 1472 1473! \author> Stefan Knecht, January 2015 1474! ------------------------------------------------------------------ 1475 implicit none 1476! ------------------------------------------------------------------ 1477 real*8 , intent(inout) :: dv(nnashx) 1478 integer, intent(in) :: nisht 1479 integer, intent(in) :: nasht 1480 integer, intent(in) :: nbast 1481 integer, intent(in) :: nnashx 1482 integer, intent(in) :: ncmot 1483 integer, intent(in) :: nsym 1484 integer, intent(in) :: nbas(nsym) 1485 integer, intent(in) :: ibas(nsym) 1486 integer, intent(in) :: lupri 1487 logical, intent(in) :: hsrohf 1488 character(len=2), intent(in) :: density 1489! ------------------------------------------------------------------ 1490 integer :: i, j 1491 integer :: lscr 1492 real*8 :: vdummy(2) 1493 real*8, allocatable :: cmo(:) 1494 real*8, allocatable :: dvao(:,:,:) 1495 real*8, allocatable :: dcao(:,:,:) 1496 real*8, allocatable :: scratch(:) 1497 logical :: print_to_screen = .true. 1498 logical :: ex = .false. 1499! ------------------------------------------------------------------ 1500 1501 call qenter('write_aodens') 1502 1503 !> open the unformatted file AO-densities 1504 inquire(FILE='AO-densities',exist=ex) 1505 if(ex)then 1506 open(1234,file='AO-densities',status='old', 1507 & form='unformatted',access='sequential', 1508 & action='readwrite',position='append') 1509 else 1510 open(1234,file='AO-densities',status='replace', 1511 & form='unformatted',access='sequential', 1512 & action='readwrite',position='rewind') 1513 end if 1514 1515 !> write AO basis information to file 1516 write(1234) nsym, nbast 1517 write(1234) nbas(1:nsym) 1518 write(1234) ibas(1:nsym) 1519 1520 !> allocate scratch space 1521 lscr = nasht*nbast + nasht*nasht + 20 ! the 20 is for the mem-markers 1522 allocate(cmo(ncmot)); cmo = 0 1523 allocate(dcao(nbast,nbast,1)); dcao = 0 ! quadratic matrix of dimension nbast * nbast with nbast == total number of AO-basis functions 1524 allocate(dvao(nbast,nbast,1)); dvao = 0 ! quadratic matrix of dimension nbast * nbast with nbast == total number of AO-basis functions 1525 allocate(scratch(lscr)); scratch = 0 1526 1527 !> read MOs from SIRIUS.RST 1528 !> a. if HF these are the final HF orbitals 1529 !> b. if MC these are the final MCSCF orbitals 1530 call readmo(cmo,9) 1531 1532 select case(trim(density)) 1533 1534 case("HF") 1535 1536 !> construct active part (high-spin open-shell HF) of density matrix 1537 if(hsrohf)then 1538 call dzero(dv,nnashx) 1539 do i = 1, nasht 1540 j = i*(i+1)/2 1541 dv(j) = 1.0d0 1542 end do 1543 end if 1544 1545 call fckden((nisht>0),(nasht>0),dcao,dvao,cmo,dv,scratch,lscr) 1546 1547 if(hsrohf)then 1548 call daxpy(nbast**2,1.0d0,dvao,1,dcao,1) 1549 end if 1550 1551 !> write the combined (inactive + open-shell) AO-density to file 1552 write(1234) dcao 1553 1554 case("MC") 1555 1556 !> input: active part of 1-particle density matrix in dv (MO-basis) 1557 1558 !> note: change ".false." to ".true" below to get 1559 ! the inactive part of the density matrix in dcao 1560 call fckden(.false.,(nasht>0),dcao,dvao,cmo,dv,scratch,lscr) 1561 1562 !> write the AO-density (active part only) to file 1563 write(1234) dvao 1564 1565 if(.false.)then 1566 !> write the AO-density (inactive part only) to file 1567 write(1234) dcao 1568 end if 1569 1570 case default 1571 call quit("write_aodens: unknown density type") 1572 end select 1573 1574 !> print AO-density matrices to screen 1575 if(print_to_screen)then 1576 if(nisht > 0 .and. trim(density)== "HF")then 1577 write(lupri,'(/A)') 1578 & ' write_aodens: dcao = density matrix, inactive part (AO-basis)' 1579 call output(dcao,1,nbast,1,nbast,nbast,nbast,1,lupri) 1580 end if 1581 if(nasht > 0) then 1582 write(lupri,'(/a)') 1583 & ' write_aodens: dvao = density matrix, active part (AO-basis)' 1584 call output(dvao,1,nbast,1,nbast,nbast,nbast,1,lupri) 1585 end if 1586 end if 1587 1588 deallocate(cmo,dvao,dcao,scratch) 1589 1590 !> close the file "AO-densities" 1591 close(1234, status="keep") 1592 1593 call qexit('write_aodens') 1594 1595 end subroutine write_aodens 1596 1597C /* Deck ijkaux */ 1598 SUBROUTINE IJKAUX(II,JJ,KK,LL) 1599#include "implicit.h" 1600#include "priunit.h" 1601#include "ccorb.h" 1602C ... only called from cc/ routines, 1603C thus NSYM from ccorb.h and not from inforb.h /hjaaj checked aug 04 1604#include "maxorb.h" 1605#include "r12int.h" 1606 LOGICAL FIRST 1607 DATA FIRST /.TRUE./ 1608 DIMENSION IJ(MXCORB) 1609 SAVE IJ 1610 IF (FIRST) THEN 1611 KL = 0 1612 MN = 0 1613 DO ISYM = 1, NSYM 1614 DO I = 1, MBAS1(ISYM) + MBAS2(ISYM) 1615 KL = KL + 1 1616 IF (I .LE. MBAS1(ISYM)) THEN 1617 MN = MN + 1 1618 IJ(KL) = MN 1619 ELSE 1620 IJ(KL) = -1 1621 END IF 1622 END DO 1623 END DO 1624 FIRST = .FALSE. 1625 END IF 1626 II = IJ(II) 1627 JJ = IJ(JJ) 1628 KK = IJ(KK) 1629 LL = IJ(LL) 1630 RETURN 1631 END 1632C /* Deck set_inoaux */ 1633 SUBROUTINE SET_INOAUX(IJ) 1634#include "implicit.h" 1635#include "priunit.h" 1636#include "ccorb.h" 1637C ... only used in RDONEL when called from cc/ routines, 1638C thus NSYM from ccorb.h and not from inforb.h /hjaaj checked aug 04 1639#include "r12int.h" 1640 1641 DIMENSION IJ(*) 1642 1643 KL = 0 1644 MN = 0 1645 DO ISYM = 1, NSYM 1646 DO I = 1, MBAS1(ISYM) + MBAS2(ISYM) 1647 DO J = 1, I 1648 KL = KL + 1 1649 IF (I. LE. MBAS1(ISYM) .AND. J .LE. MBAS1(ISYM)) THEN 1650 MN = MN + 1 1651 IJ(KL) = MN 1652 ELSE 1653 IJ(KL) = -1 1654 END IF 1655 END DO 1656 END DO 1657 END DO 1658 1659 RETURN 1660 END 1661C /* Deck rdsupm */ 1662 SUBROUTINE RDSUPM(NSIM,FMAT,DMAT,ISYMDM,WORK,LFREE) 1663C 1664C Jan 90, Hans Joergen Aa. Jensen 1665C 2-Apr-1997 hjaaj (FMAT and DMAT now (NBAST,NBAST)) 1666C Oct-2014 hjaaj (ISYMDM in parameter list) 1667C 1668#include "implicit.h" 1669#include "priunit.h" 1670 DIMENSION FMAT(*), DMAT(*), ISYMDM(*), WORK(*) 1671C 1672C Used from common blocks: 1673C gnrinf.h: DOSRIN 1674C INFORB : NBAST,N2BASX 1675C INFTAP : LUSUPM 1676C CBIHRS : parameters for FORMSUP 1677C frame.h : POTNUC 1678C DFTCOM : HFXFAC 1679C 1680#include "gnrinf.h" 1681#include "inforb.h" 1682#include "inftap.h" 1683#include "cbihrs.h" 1684#include "frame.h" 1685#include "dftcom.h" 1686C 1687C Local variables: 1688C 1689 LOGICAL FNDLAB, FNSUPM_OK, F_EXIST, DOSRIN_save 1690 LOGICAL*4 NOSYM 1691 character*8 FNSUPM_local 1692 1693 CALL QENTER('RDSUPM') 1694C 1695 IF (LUSUPM .EQ. -1) THEN 1696C ... special code for AOSUPMAT not to be used 1697 WRITE (LUPRI,'(/A/A)') 1698 & 'PROGRAMMING ERROR: RDSUPM called with LUSUPM = -1', 1699 & 'Please report on Dalton forum.' 1700 CALL QTRACE(LUPRI) 1701 CALL QUIT('Programming ERROR: RDSUPM called with LUSUPM.eq.-1') 1702 END IF 1703 IF (ONLY_J) THEN 1704 HFXFACSUP = 0.0D0 1705 ELSE 1706 HFXFACSUP = HFXFAC 1707 END IF 1708 if (HFXFACSUP .eq. 0.0D0) then 1709 FNSUPM_local = 'AO2_JINT' 1710 else if (HFXFACSUP .eq. 1.0D0) then 1711 FNSUPM_local = FNSUPM 1712 else 1713 FNSUPM_local = 'AODFTINT' 1714 end if 1715 if (LUSUPM .gt. 0) call GPCLOSE(LUSUPM,'KEEP') 1716C 1717 I = 0 1718 100 IF (LUSUPM .LE. 0) CALL GPOPEN(LUSUPM,FNSUPM_local, 1719 & 'UNKNOWN',' ','UNFORMATTED',IDUMMY,.FALSE.) 1720 I = I + 1 1721 REWIND(LUSUPM) 1722 1723 FNSUPM_OK = FNDLAB('PXSUPMAT',LUSUPM) 1724 IF (FNSUPM_OK) THEN 1725 READ (LUSUPM) XFAC, XPOTNUC, NOSYM 1726 IF (abs(XFAC-HFXFACSUP) .gt. THRSUP) THEN 1727 WRITE(LUPRI,*) 1728 WRITE(LUPRI,*) 'Calculating ',FNSUPM_local, 1729 & ' because HFXFAC changed from',XFAC,' to',HFXFACSUP 1730 FNSUPM_OK = .FALSE. 1731 END IF 1732 IF (abs(XPOTNUC-POTNUC) .gt. THRSUP) THEN 1733 WRITE(LUPRI,'(/3A)') 'Calculating ',FNSUPM_local, 1734 & ' because geometry has changed.' 1735 FNSUPM_OK = .FALSE. 1736 END IF 1737 IF (NOSSUP .AND. .NOT.NOSYM) THEN 1738 FNSUPM_OK = .FALSE. 1739 WRITE(LUPRI,'(/3A)') 'Calculating ',FNSUPM_local, 1740 & ' because .NOSYM requested' 1741 FNSUPM_OK = .FALSE. 1742 END IF 1743 ELSE 1744 WRITE(LUPRI,'(/2A)') 'Calculating ',FNSUPM_local 1745 END IF 1746 1747 IF (.NOT. FNSUPM_OK) THEN 1748 IF (I.GT.1) CALL QUIT('CALL FORMSUP failed') 1749 inquire(file='AOTWOINT',EXIST=F_EXIST) 1750 IF (.NOT. F_EXIST) THEN 1751 DOSRIN_save = DOSRIN 1752 DOSRIN = .FALSE. ! do not generate AOSR2INT here when srDFT 1753 CALL newTIMER('START ') 1754 CALL MAKE_AOTWOINT(WORK,LFREE) 1755 CALL newTIMER('MAKE_AOTWOINT') 1756 CALL FLSHFO(LUPRI) 1757 DOSRIN = DOSRIN_save 1758 END IF 1759 CALL newTIMER('START ') 1760 CALL FORMSUP(WORK,LFREE,NOSSUP,HFXFACSUP,THRSUP,IPRSUP) 1761 CALL newTIMER('FORMSUP') 1762 CALL FLSHFO(LUPRI) 1763 GO TO 100 1764 END IF 1765C 1766C Fold density matrices 1767C 1768 IF (LFREE .LT. N2BASX) CALL ERRWRK('RDSUPM',-N2BASX,LFREE) 1769 DO I = 1,NSIM 1770 JDG = 1 + (I-1)*N2BASX 1771 CALL DCOPY(N2BASX,DMAT(JDG),1,WORK,1) 1772 CALL DGEFSP(NBAST,WORK,DMAT(JDG)) 1773 END DO 1774C 1775C note equivalence(WORK, IWORK) 1776C 1777 CALL RDSUPP(NSIM,FMAT,DMAT,ISYMDM,WORK,WORK,LFREE) 1778C 1779C square Fock matrices, restore DCAO 1780C 1781 DO I = 1,NSIM 1782 JDG = 1 + (I-1)*N2BASX 1783 CALL DCOPY(NNBASX,FMAT(JDG),1,WORK,1) 1784 CALL DSPTGE(NBAST,WORK,FMAT(JDG)) 1785 CALL DCOPY(NNBASX,DMAT(JDG),1,WORK,1) 1786 CALL DUNFLD(NBAST,WORK,DMAT(JDG)) 1787 END DO 1788C 1789 CALL GPCLOSE(LUSUPM,'KEEP') 1790 CALL QEXIT('RDSUPM') 1791 RETURN 1792 END 1793C /* Deck rdsupp */ 1794 SUBROUTINE RDSUPP(NSIM,FMAT,DMAT,ISYMDM,WORK,IWORK,LFREE) 1795C 1796C Jan 90 + Oct 16, Hans Joergen Aa. Jensen 1797C 1798C NOTE: WORK and IWORK must be equivalenced. 1799C 1800#include "implicit.h" 1801#include "thrzer.h" 1802 PARAMETER (D0 = 0.0D0) 1803C 1804 REAL*8 FMAT(N2BASX,NSIM), DMAT(N2BASX,NSIM), WORK(LFREE) 1805 INTEGER*4 IWORK(2*LFREE) 1806 INTEGER ISYMDM(NSIM) 1807C 1808C Used from common blocks: 1809C 1810C INFORB : N2BASX,MULD2H(8,8),... 1811C INFIND : IROW(:),ISAO(:) 1812C INFTAP : LUSUPM 1813C PRIUNIT: LUERR, IPRSTAT 1814C 1815#include "maxash.h" 1816#include "maxorb.h" 1817#include "priunit.h" 1818#include "inforb.h" 1819#include "infind.h" 1820#include "inftap.h" 1821C 1822! always *4 on LUSUPM 1823 LOGICAL*4 NOSYM 1824 INTEGER*4 MSYM, MBAS(8) 1825 INTEGER*4 ITYP,NP1,NQ1,NPL,NQL,NPQRS_4 1826C 1827 CALL QENTER('RDSUPP') 1828 REWIND(LUSUPM) 1829 CALL MOLLAB('PXSUPMAT',LUSUPM,LUPRI) 1830C 1831C *** Test if LUSUPM file is OK 1832C 1833 READ (LUSUPM) XFAC,XPOTNUC,NOSYM, MSYM, MBAS 1834 NERR = 0 1835 ISYMDM_ALL = ISYMDM(1) 1836 DO I = 1,NSIM 1837 IF (ISYMDM(I) .NE. ISYMDM_ALL) ISYMDM_ALL = -1 ! not all d.m. have same symmetry 1838 IF (ISYMDM(I) .NE. 1 .AND. .NOT.NOSYM) NERR = NERR + 1 1839 IF (ISYMDM(I) .LT. 1 .OR. ISYMDM(I) .GT. NSYM) NERR = NERR + 1 1840 END DO 1841 IF (MSYM .NE. NSYM) THEN 1842 NERR = NERR + 10 1843 ELSE 1844 DO 10 ISYM = 1,NSYM 1845 10 IF (MBAS(ISYM) .NE. NBAS(ISYM)) NERR = NERR + 10 1846 END IF 1847 IF (NERR .GT. 0) THEN 1848 WRITE (LUPRI,'(//5X,A)') 1849 & 'RDSUPP ERROR: AOSUPINT inconsistent with Sirius input' 1850 IF (MOD(NERR,10) .NE. 0) WRITE (LUPRI,'(/5X,2A)') 1851 & 'AOSUPINT must be sorted with NOSYM when ', 1852 & 'perturbation symmetry .ne. 1' 1853 IF (NERR .GE. 10) WRITE (LUPRI,'(/5X,A,2I3/,(5X,A,8I5))') 1854 & 'AOSUPINT vs. AOONEINT, NSYM =',MSYM,NSYM, 1855 & 'NBAS from AOSUPINT :',(MBAS(I),I=1,8), 1856 & 'NBAS from AOONEINT :',(NBAS(I),I=1,8) 1857 CALL QTRACE(LUPRI) 1858 CALL QUIT('RDSUPP ERROR: AOSUPINT inconsistent with AOONEINT') 1859 END IF 1860 IF (IPRSTAT .GT. 5) THEN 1861 WRITE (LUERR,*) 'Test output from RDSUPP, XFAC,NOSYM =', 1862 & XFAC,NOSYM 1863 END IF 1864C 1865C 1866C *** Start loop over psupmat file 1867C 1868 IREC = 0 1869 IREC1 = 0 1870 IREC2 = 0 1871 100 CONTINUE 1872 IREC = IREC + 1 1873 READ (LUSUPM) ITYP,NP1,NQ1,NPL,NQL,NPQRS_4 1874 NPQRS = NPQRS_4 ! in case default is integer*8 1875 IF (IPRSTAT .GT. 20) THEN 1876 WRITE (LUERR,*) 'ITYP,NP1,NQ1,NPL,NQL,NPQRS' 1877 WRITE (LUERR,'(7I10)') ITYP,NP1,NQ1,NPL,NQL,NPQRS 1878 END IF 1879C ITYP = -2 finished 1880C 1 P(1:NPQRS) from NP1,NQ1 to NPL,NQL 1881C 2 P(1:NPQRS), INDP(1:NPQRS) with (RS) addresses and np1=npl, nq1=nql 1882 IF (ITYP .EQ. -2) GO TO 900 1883C 1884 ISP1 = ISAO(NP1) 1885 ISPL = ISAO(NPL) 1886 ISQ1 = ISAO(NQ1) 1887 ISQL = ISAO(NQL) 1888 if (ISYMDM_ALL .gt. 0) then 1889 IF (ISP1.EQ.ISPL .AND. ISQ1.EQ.ISQL) THEN 1890 ISPQ = MULD2H(ISP1,ISQ1) 1891 IF (ISPQ .NE. ISYMDM_ALL) THEN 1892 READ (LUSUPM) 1893 GO TO 100 1894 END IF 1895 END IF 1896 end if 1897C 1898 IF (ITYP .EQ. 1) THEN 1899 IREC1 = IREC1 + 1 1900 IF (NPQRS .GT. LFREE) GO TO 810 1901 CALL READT(LUSUPM,NPQRS,WORK) 1902 DO 280 ISIM = 1,NSIM 1903 IOFF = 0 1904 DO 260 IP = NP1,NPL 1905 ISP = ISAO(IP) 1906 IF (IP .EQ. NPL) THEN 1907 NQEND = NQL 1908 ELSE 1909 NQEND = IP 1910 END IF 1911 IF (IP .EQ. NP1) THEN 1912 NQBEG = NQ1 1913 ELSE 1914 NQBEG = 1 1915 END IF 1916 DO 240 IQ = NQBEG, NQEND 1917 ISQ = ISAO(IQ) 1918 ISPQ = MULD2H(ISQ,ISP) 1919 IPQ = IROW(IP) + IQ 1920 IF (ISPQ .EQ. ISYMDM(ISIM)) THEN 1921 DPQ = DMAT(IPQ,ISIM) 1922C note NRS = IPQ 1923 IF (ABS(DPQ) .GT. THRZER) 1924 & CALL DAXPY(IPQ,DPQ,WORK(IOFF+1),1,FMAT(1,ISIM),1) 1925 FMAT( IPQ , ISIM ) = FMAT( IPQ , ISIM ) + 1926 & DDOT(IPQ,WORK(IOFF+1),1,DMAT(1,ISIM),1) 1927 END IF 1928 IOFF = IOFF + IPQ 1929 240 CONTINUE 1930 NQ1 = 1 1931 260 CONTINUE 1932 IF (IOFF .NE. NPQRS .OR. NQL .GT. NPL) THEN 1933 WRITE (LUPRI,'(//5X,A,3(/5X,A,I12)/5X,A/5X,5I5,2I12)') 1934 & 'RDSUPP ERROR, inconsistent type 1 specifications', 1935 & ' in record pair no.',IREC, 1936 & ' Number of integrals from NP1,NQ1,NPL,NQL:',IOFF, 1937 & ' Number of integrals from NPQRS :',NPQRS, 1938 & ' dump of ITYP,NP1,NQ1,NPL,NQL,NPQRS:', 1939 & ITYP,NP1,NQ1,NPL,NQL,NPQRS 1940 CALL QTRACE(LUPRI) 1941 CALL QUIT( 1942 & 'RDSUPP ERROR, inconsistent type 1 spec. on LUSUPM') 1943 END IF 1944 280 CONTINUE 1945 ELSE IF (ITYP .EQ. 2) THEN 1946 IREC2 = IREC2 + 1 1947 LENP = 3*NPQRS ! real*8 P + integer*4 INDP 1948 IF (LENP .GT. LFREE) GO TO 820 1949 CALL READI4(LUSUPM,LENP,IWORK) 1950 INDP = 2*NPQRS ! "2" as IWORK is always integer*4 1951 IPQ = IROW(NP1) + NQ1 1952 DO 380 ISIM = 1,NSIM 1953 DPQ = DMAT(IPQ,ISIM) 1954 IF (ABS(DPQ) .GT. THRZER) THEN 1955! no vector dependence in IRS loop 1956 DO 320 IRS = 1,NPQRS 1957 FMAT( IWORK(INDP+IRS) , ISIM ) = 1958 & FMAT( IWORK(INDP+IRS) , ISIM ) + DPQ * WORK(IRS) 1959 320 CONTINUE 1960 END IF 1961 FPQ = D0 1962 DO 340 IRS = 1,NPQRS 1963 FPQ = FPQ + DMAT( IWORK(INDP+IRS) ,ISIM ) * WORK(IRS) 1964 340 CONTINUE 1965C Use first line below if P(pq,pq) has been divided by two 1966C otherwise use second line. 1967 FMAT( IPQ , ISIM ) = FMAT( IPQ , ISIM ) + FPQ 1968C FMAT( IPQ , ISIM ) = FPQ 1969 380 CONTINUE 1970 ELSE 1971 WRITE (LUPRI,'(//5X,A,I5/5X,A,I12/5X,A/5X,5I5,2I12)') 1972 & 'RDSUPP ERROR, unrecognized type on LUSUPM:',ITYP, 1973 & ' record pair no.',IREC, 1974 & ' dump of ITYP,NP1,NQ1,NPL,NQL,NPQRS:', 1975 & ITYP,NP1,NQ1,NPL,NQL,NPQRS 1976 CALL QTRACE(LUPRI) 1977 CALL QUIT('RDSUPP ERROR, unrecognized type on AOSUPINT') 1978 END IF 1979 GO TO 100 1980C 1981 810 CONTINUE 1982 CALL ERRWRK('RDSUPP.TYP1',LENP,LFREE) 1983 820 CONTINUE 1984 CALL ERRWRK('RDSUPP.TYP2',LENP,LFREE) 1985 900 CONTINUE 1986 IF (IPRSTAT .GT. 5) THEN 1987 WRITE (LUERR,*) 'RDSUPP statistics.' 1988 WRITE (LUERR,*) 'Total number of records :',IREC 1989 WRITE (LUERR,*) 'No. of records of type 1:',IREC1 1990 WRITE (LUERR,*) 'No. of records of type 2:',IREC2 1991 END IF 1992C 1993 CALL QEXIT('RDSUPP') 1994 RETURN 1995 END 1996C /* Deck sirh1 */ 1997 SUBROUTINE SIRH1(H1AO,WRK,LWRK) 1998C 1999C Copyright 29_nov-1990 Hans Joergen Aa. Jensen 2000C 2001#include "implicit.h" 2002 DIMENSION H1AO(*), WRK(LWRK) 2003C 2004 PARAMETER ( D0 = 0.0D0 ) 2005C 2006C Used from common blocks: 2007C INFINP : NFIELD, EFIELD, LFIELD 2008C INFORB : NNBAST,... 2009C 2010#include "maxorb.h" 2011#include "priunit.h" 2012#include "thrzer.h" 2013#include "infinp.h" 2014#include "inforb.h" 2015#include "infpri.h" 2016C 2017 LOGICAL ANTSYM 2018C 2019 CALL QENTER('SIRH1') 2020 IF (P6FLAG(18)) WRITE (LUPRI,'(/A)') 2021 & ' ----- One-electron Hamiltonian matrix from SIRH1 :' 2022 CALL RDONEL('ONEHAMIL',.TRUE.,H1AO,NNBAST) 2023C 2024C Add field-dependent terms 2025C 2026 IF (NFIELD .LE. 0) GO TO 800 2027 KPRPAO = 1 2028 KH1AO = KPRPAO + NNBASX 2029 LNEED = KH1AO + NNBAST 2030 IF (LNEED .GT. LWRK) CALL ERRWRK('SIRH1.RDPROP',LNEED,LWRK) 2031 DO 100 IFIELD = 1,NFIELD 2032 IF (P6FLAG(18)) WRITE (LUPRI,'(A,1P,G15.6,A,A8)') 2033 & ' Field strength :',EFIELD(IFIELD), 2034 & ' Field type : ',LFIELD(IFIELD) 2035 IF (EFIELD(IFIELD) .NE. D0) THEN 2036 CALL RDPROP(LFIELD(IFIELD),WRK,ANTSYM) 2037 IF (.NOT. ANTSYM) THEN 2038 CALL PKSYM1(WRK(KPRPAO),WRK(KH1AO),NBAS,NSYM,1) 2039 DPX = DASUM(NNBASX,WRK(KPRPAO),1) 2040 DPT = DASUM(NNBAST,WRK( KH1AO),1) 2041 IF (ABS(DPX - DPT) .GT. THRZER*DPX) THEN 2042 WRITE (LUPRI,'(/3A,/2A,2(/A,1P,D25.15))') 2043 & ' FATAL ERROR in SIRH1: ',LFIELD(IFIELD), 2044 & ' is not totally symmetric,', 2045 & ' and finite field is only allowed with', 2046 & ' totally symmetric operators.', 2047 & ' Absolute sum of lower triangle of property matrix:', 2048 & DPX, 2049 & ' Absolute sum of the totally symmetric elements :', 2050 & DPT 2051 CALL QUIT( 2052 & 'ERROR: finite field with non tot.sym. operator') 2053 END IF 2054 CALL DAXPY(NNBAST,EFIELD(IFIELD),WRK(KH1AO),1,H1AO,1) 2055 ELSE 2056 WRITE (LUPRI,'(/3A/A)') ' FATAL ERROR in SIRH1: ', 2057 & LFIELD(IFIELD),' is antisymmetric (i.e. imaginary)', 2058 & ' and finite field with imaginary operators is not allowed.' 2059 CALL QUIT('ERROR: finite field with imaginary operator') 2060 END IF 2061 END IF 2062 100 CONTINUE 2063 800 IF (P6FLAG(18)) CALL OUTPKB(H1AO,NBAS,NSYM,1,LUPRI) 2064 CALL QEXIT('SIRH1') 2065 RETURN 2066 END 2067C /* Deck rdprop */ 2068 SUBROUTINE RDPROP(WORD,PRPAO,ANTSYM) 2069C 2070C Written 28-Nov-1990 HJAaJ 2071C 2072C Read AO property integrals with label WORD from LUPROP 2073C Output: 2074C PRPAO : AO property integrals 2075C ANTSYM : true if integral matrix is antisymmetric 2076C false if matrix is symmetric 2077C 2078#include "implicit.h" 2079#include "dummy.h" 2080 CHARACTER*8 WORD 2081 DIMENSION PRPAO(*) 2082 LOGICAL ANTSYM 2083C 2084C -- common blocks: 2085C INFORB : NNBASX, N2BASX 2086C 2087#include "priunit.h" 2088#include "inftap.h" 2089#include "inforb.h" 2090C 2091 CHARACTER*8 RTNLBL(2) 2092 LOGICAL OPEND, FNDLB2 2093C 2094 CALL QENTER('RDPROP') 2095 IF (LUPROP .GT. 0) CALL GPCLOSE(LUPROP,'KEEP') 2096 CALL GPOPEN(LUPROP,'AOPROPER', 2097 & 'OLD',' ','UNFORMATTED',IDUMMY,.FALSE.) 2098 REWIND LUPROP 2099 IF ( FNDLB2(WORD,RTNLBL,LUPROP)) THEN 2100 IF (RTNLBL(2).EQ.'SQUARE ') THEN 2101 CALL READT(LUPROP,N2BASX,PRPAO) 2102 ANTSYM = .FALSE. 2103 ELSE IF (RTNLBL(2).EQ.'SYMMETRI') THEN 2104 ANTSYM = .FALSE. 2105 CALL READT(LUPROP,NNBASX,PRPAO) 2106 ELSE IF (RTNLBL(2).EQ.'ANTISYMM') THEN 2107 ANTSYM = .TRUE. 2108 CALL READT(LUPROP,NNBASX,PRPAO) 2109 ELSE 2110 WRITE (LUPRI,'(//3A/3(A,1X))') 2111 & ' --- RDPROP error: property "',WORD, 2112 & '" on AOPROPER has no valid antisymmetry label.', 2113 & ' Labels read on file:',RTNLBL(1),RTNLBL(2) 2114 CALL QTRACE(LUPRI) 2115 CALL QUIT('RDPROP error: No antisymmetry label on AOPROPER') 2116 END IF 2117 ELSE 2118 WRITE (LUPRI,'(//3A)') ' --- RDPROP error: property "',WORD, 2119 * '" not found on AOPROPER.' 2120 CALL QTRACE(LUPRI) 2121 CALL QUIT('RDPROP error: property not found on AOPROPER.') 2122 END IF 2123 CALL GPCLOSE(LUPROP,'KEEP') 2124 CALL QEXIT('RDPROP') 2125 RETURN 2126 END 2127C /* Deck rd_sirifc */ 2128 SUBROUTINE RD_SIRIFC(KEY,FOUND,AMAT) 2129C 2130C Read specific information from Sirius interface file 2131C 2132C 10-Dec-1992 HJAaJ 2133C Revision 2000/05/24 HJAaJ: new options (read CMO and others) 2134C Revision 2011/05/02 HJAaJ: new keys FC, FCDIAG 2135C 2136C The following records have been written by WR_SIRIFC: 2137C 2138C 0) label LBSIFC ("SIR IPH ") 2139C 1) POTNUC,EMY,EACTIV,EMCSCF,ISTATE,ISPIN,NACTEL,LSYM,MS2 2140C 2) NISHT,NASHT,... 2141C 3) CMO 2142C 4) CREF 2143C 5) DV 2144C 6) FOCK 2145C 7) PV 2146C 8) FC 2147C 9) FV 2148C 10) FCAC 2149C 11) H2AC 2150C 12) GORB 2151C If (GUGA) then 2152C 13) label "CIDIAG1 " 2153C 14) CIDIAG1 2154C end if 2155C 15) label "CIDIAG2 " 2156C 16) CIDIAG2 2157C 17) label "ORBDIAG " 2158C 18) ORBDIAG 2159C 2160C *) label "SIRFLAGS" 2161C *) FLAG(1:NFLAG) 2162C *) GRDNRM,POTNUC,EMY,EACTIV,ESOLT,EMCSCF 2163C If (fields) then 2164C *) label "EXTFIELD" 2165C *) NFIELD 2166C *) (EFIELD(I),I=1,NFIEL4) 2167C *) (LFIELD(I),I=1,NFIEL4) 2168C end if 2169C If (solvent) then 2170C *) label "SOLVINFO" 2171C *) EPSOL,EPSTAT,EPPN,RSOL(1:3),LSOLMX,INERSI,INERSF 2172C *) GRDNRM,POTNUC,EMY,EACTIV,ESOLT,EMCSCF 2173C *) ERLM(LM,1), LM = 1,NLMSOL) 2174C *) (TRLM(LM), LM = 1,NLMSOL) where TRLM(i) = ERLM(i,2) 2175C *) NSYM, NBAS 2176C *) label "SOLVTMAT" 2177C *) TMAT(1:NNORBX) 2178C end if 2179C If (CSF's) then 2180C *) label "CREFDETS" 2181C *) CREF in dets 2182C end if 2183C *) label 'TRCCINT" 2184C *) NSYM, NORBT, ... 2185C *) FCDIA, ISMO 2186C *) CMO 2187C 2188#include "implicit.h" 2189#include "priunit.h" 2190C 2191 CHARACTER*(*) KEY 2192 LOGICAL FOUND 2193 DIMENSION AMAT(*) 2194C 2195C Used from common blocks: 2196C INFINP : NLMSOL 2197C INFORB : NSYM,NCMOT,NNORBT,NNASHX 2198C INFOPT : POTNUC, EMY 2199C INFTAP : LUSIFC,FNSIFC,LBSIFC 2200C 2201#include "maxorb.h" 2202#include "infinp.h" 2203#include "inforb.h" 2204#include "inftap.h" 2205C 2206 LOGICAL CLSIFC, FNDLAB 2207C 2208 CALL QENTER('RD_SIRIFC') 2209 IF (LUSIFC .GT. 0) THEN 2210 CLSIFC = .FALSE. 2211 ELSE 2212 CLSIFC = .TRUE. 2213 CALL GPOPEN(LUSIFC,FNSIFC, 2214 & 'UNKNOWN',' ','UNFORMATTED',IDUMMY,.FALSE.) 2215 END IF 2216 REWIND LUSIFC 2217 IF (KEY .EQ. 'CMO') THEN 2218 IF (.NOT. FNDLAB(LBSIFC,LUSIFC)) THEN 2219 REWIND LUSIFC 2220 IF (.NOT. FNDLAB('CIRESPON',LUSIFC)) GO TO 8888 2221 ENDIF 2222 READ (LUSIFC) 2223 READ (LUSIFC) 2224 CALL READT(LUSIFC,NCMOT,AMAT) 2225 ELSE IF (KEY .EQ. 'DV') THEN 2226 IF (.NOT. FNDLAB(LBSIFC,LUSIFC)) GO TO 8888 2227 READ (LUSIFC) 2228 READ (LUSIFC) 2229 READ (LUSIFC) 2230 READ (LUSIFC) 2231 CALL READT(LUSIFC,NNASHX,AMAT) 2232 ELSE IF (KEY .EQ. 'PV') THEN 2233 IF (.NOT. FNDLAB(LBSIFC,LUSIFC)) GO TO 8888 2234 READ (LUSIFC) 2235 READ (LUSIFC) 2236 READ (LUSIFC) 2237 READ (LUSIFC) 2238 READ (LUSIFC) 2239 READ (LUSIFC) 2240 CALL READT(LUSIFC,NNASHX*NNASHX,AMAT) 2241 ELSE IF (KEY .EQ. 'FOCK') THEN 2242 IF (.NOT. FNDLAB(LBSIFC,LUSIFC)) GO TO 8888 2243 READ (LUSIFC) 2244 READ (LUSIFC) 2245 READ (LUSIFC) 2246 READ (LUSIFC) 2247 READ (LUSIFC) 2248 CALL READT(LUSIFC,N2ORBT,AMAT) 2249 ELSE IF (KEY .EQ. 'FC for MP2') THEN ! 3-May-2011 hjaaj, specifically for use in MP2FCK 2250 IF (.NOT. FNDLAB(LBSIFC,LUSIFC)) GO TO 8888 2251 ! rec no. 1 2252! WRITE (LUSIFC) POTNUC,EMY,EACTIV,EMCSCF, 2253! & ISTATE,ISPIN,NACTEL,LSYM,MS2 2254 READ (LUSIFC) POTNUC_IFC, EMY,EACTIV,EMCSCF ! we need EMY for SIRMP2 2255 EMY = EMCSCF - POTNUC_IFC ! and we need EMY to include ESOLT 2256 IF (POTNUC .NE. POTNUC_IFC) 2257 & CALL QUIT('POTNUC on SIRIFC and input are different') 2258 ! rec no. 2 2259! WRITE (LUSIFC) NISHT,NASHT,NOCCT,NORBT,NBAST,NCONF,NWOPT,NWOPH, 2260! & NCDETS, NCMOT,NNASHX,NNASHY,NNORBT,N2ORBT, 2261! & NSYM,MULD2H, NRHF,NFRO, 2262! & NISH,NASH,NORB,NBAS, 2263! & NELMN1, NELMX1, NELMN3, NELMX3, MCTYPE, 2264! & NAS1, NAS2, NAS3 2265 READ (LUSIFC) NISHT_IFC,NASHT_IFC,idum,NORBT_IFC,NBAST_IFC, 2266 & (idum,i=1,9),NSYM_IFC 2267 IF (NISHT .ne. NISHT_IFC .or. NASHT .ne. NASHT_IFC .or. 2268 & NORBT .ne. NORBT_IFC .or. NBAST .ne. NBAST_IFC .or. 2269 & NSYM .ne. NSYM_IFC ) CALL 2270 & QUIT('SIRIFC info not consistent with this calculation.') 2271 2272 DO IREC = 3,7 ! skip records 3:7 2273 READ (LUSIFC) 2274 END DO 2275 CALL READT(LUSIFC,NNORBT,AMAT) ! FC matrix is in rec. no. 8 2276 ELSE IF (KEY .EQ. 'SOLTLM') THEN 2277 IF (.NOT. FNDLAB('SOLVINFO',LUSIFC)) GO TO 8888 2278 READ (LUSIFC) 2279 READ (LUSIFC) 2280 READ (LUSIFC) 2281 CALL READT(LUSIFC,NLMSOL,AMAT) 2282 ELSE IF (KEY .EQ. 'SOLVTMAT') THEN 2283 IF (.NOT. FNDLAB('SOLVTMAT',LUSIFC)) GO TO 8888 2284 CALL READT(LUSIFC,NNORBT,AMAT) 2285 ELSE IF (KEY .EQ. 'FCDIAG') THEN 2286 IF (.NOT. FNDLAB('TRCCINT ',LUSIFC)) GO TO 8888 2287 READ (LUSIFC) NSYM_SIRIFC 2288 IF (NSYM .NE. NSYM_SIRIFC) THEN 2289 CALL QUIT('RD_SIRIFC - key FCDIAG : NSYM inconsistent') 2290 END IF 2291 CALL READT(LUSIFC,NORBT,AMAT) 2292 ELSE 2293 WRITE(LUPRI,*) 'RD_SIRIFC error, invalid key word ',KEY 2294 CALL QUIT('RD_SIRIFC error, invalid key word ') 2295 END IF 2296 FOUND = .TRUE. 2297 GO TO 9999 2298C 2299 8888 FOUND = .FALSE. 2300 9999 IF (CLSIFC) CALL GPCLOSE (LUSIFC,'KEEP') 2301 CALL QEXIT('RD_SIRIFC') 2302 RETURN 2303 END 2304C /* Deck rdcref */ 2305 SUBROUTINE RDCREF(CREF,set_cref_xc_flag_lucita) 2306C 2307C 26-Sep-1990 Hans Joergen Aa. Jensen 2308C Read reference CI vector from LUIT1 2309C 2310 use lucita_mcscf_ci_cfg 2311#include "implicit.h" 2312 DIMENSION CREF(*) 2313 PARAMETER ( D1 = 1.0D0 ) 2314C 2315C Used from common blocks: 2316C INFINP : ISTATE 2317C INFVAR : NCONF 2318C INFTAP : LUIT1 2319C INFPRI : NWARN 2320C 2321#include "maxorb.h" 2322#include "priunit.h" 2323#include "infinp.h" 2324#include "infvar.h" 2325#include "inftap.h" 2326C 2327 logical, intent(in) :: set_cref_xc_flag_lucita 2328 CHARACTER*8 RTNLBL(2), STAR8 2329 DATA STAR8/'********'/ 2330C 2331 IF ( NCONF .GT. 1 ) THEN 2332 REWIND LUIT1 2333 CALL MOLLB2 ('STARTVEC',RTNLBL,LUIT1,LUPRI) 2334 IF (RTNLBL(1) .NE. STAR8) THEN 2335 READ (RTNLBL(1),'(2I4)') NRS,I 2336 IF (ABS(I) .NE. ISTATE) THEN 2337 IF (I .EQ. 0 .AND. RTNLBL(2).EQ.'CISAVE ') THEN 2338 IF (NRS .LT. ISTATE) THEN 2339 WRITE (LUPRI,'(//2(A,I4/))') 2340 & ' ERROR from RDCREF: ISTATE specified in input:',ISTATE, 2341 & ' is greater than # of CI vectors on LUIT1 (CISAVE) :',NRS 2342 CALL QUIT('RDCREF: ISTATE CI vector not available') 2343 END IF 2344 I = ISTATE 2345 ELSE 2346 NWARN = NWARN + 1 2347 WRITE (LUPRI,'(//2(A,I4/))') 2348 * ' WARNING from RDCREF: ISTATE specified in input:',ISTATE, 2349 * ' ISTATE read from LUIT1 :',ABS(I) 2350 END IF 2351 END IF 2352 ELSE 2353C this is an old LUIT1 file (before 5-Aug-1986) 2354 I = ISTATE 2355 END IF 2356 IF (I .GT. 0) THEN 2357 DO 110 I = 1,(ISTATE-1) 2358 READ (LUIT1) 2359 110 CONTINUE 2360C else only reference vector saved 2361 END IF 2362 CALL READT (LUIT1,NCONF,CREF) 2363 ELSE IF (NCONF .EQ. 1) THEN 2364 CREF(1) = D1 2365 END IF 2366 2367! set logical flag for new reference CI vector in WRK(KCREF) - 2368! in the interface to lucita we use this info to save/broadcast the new 2369! reference vector on lucita internal sequential/parallel MPI-I/O files 2370 vector_exchange_type1 = 1 2371 vector_update_mc2lu_lu2mc((2-1)*vector_exchange_types + 2372 & vector_exchange_type1) = 2373 & set_cref_xc_flag_lucita 2374C 2375 RETURN 2376 END 2377C /* Deck sirphp */ 2378 SUBROUTINE SIRPHP(DIAGL,EACTVN,XNDXCI,FCAC,H2AC,WRK,LWRK) 2379C 2380C Oct 1990 hjaaj 2381C 2382C Purpose: set up diagonal and PHP matrix for SIRIUS calculations. 2383C (Called by SIRNEO and SIRNR) 2384C 2385#include "implicit.h" 2386#include "priunit.h" 2387 DIMENSION DIAGL(*),XNDXCI(*),FCAC(*),H2AC(*),WRK(*) 2388C 2389 PARAMETER (D0 = 0.0D0, DP5 = 0.5D0, D1 = 1.0D0) 2390C 2391C Used from common blocks: 2392C INFINP : LSYM, FLAG() 2393C INFVAR : NCONF, NWOPT 2394C INFDIM : LPHPMX, MAXPHP 2395C INFTAP : LUIT2 2396C INFPRI : IPRI6 2397C 2398#include "maxorb.h" 2399#include "mxcent.h" 2400Clf#include "pcmdef.h" 2401Clf#include "pcm.h" 2402#include "pcmlog.h" 2403#include "infinp.h" 2404#include "infvar.h" 2405#include "infdim.h" 2406#include "inftap.h" 2407#include "infpri.h" 2408C 2409 CHARACTER*8 LABIT2(3) 2410 DATA LABIT2/'CIDIAG2 ','ORBDIAG ','SOLVDIAG'/ 2411C 2412 CALL QENTER('SIRPHP') 2413C 2414C Get diagonal of Hessian / L-matrix 2415C (divided by two for PHP routines) 2416C 2417 IF (NCONF .GT. 1) THEN 2418 REWIND LUIT2 2419C Find label CIDIAG2, for diagonal of CI Ham. matrix. 2420C Hessian diagonal is 2 * (CIDIAG2 - EACTVN). 2421 CALL MOLLAB(LABIT2(1),LUIT2,LUPRI) 2422 CALL READT(LUIT2,NCONF,DIAGL) 2423 IF (FLAG(16) .OR. PCM) THEN 2424C IF (SOLVNT) THEN 2425 MAXPHP = 0 2426CPHPMAERKE ... PHP not implemented for SOLVNT in this version. 2427 IF (LWRK .LT. NCONF) CALL ERRWRK('SIRPHP',NCONF,LWRK) 2428 CALL MOLLAB(LABIT2(3),LUIT2,LUPRI) 2429 CALL READT (LUIT2,NCONF,WRK) 2430 CALL DAXPY (NCONF,DP5,WRK,1,DIAGL,1) 2431 END IF 2432 DO I = 1,NCONF 2433 DIAGL(I) = DIAGL(I) - EACTVN 2434 END DO 2435 ELSE 2436 DIAGL(1) = D0 2437 END IF 2438C 2439 IF (NWOPT .GT. 0) THEN 2440 REWIND LUIT2 2441C Find label ORBDIAG, for diagonal of orbital Hessian. 2442 CALL MOLLAB(LABIT2(2),LUIT2,LUPRI) 2443 CALL READT(LUIT2,NWOPT,DIAGL(1+NCONF)) 2444 END IF 2445 IF (MAXPHP .GT. 0 .AND. NCONF .GT. 1) THEN 2446 ESHIFT = -EACTVN 2447 IPWAY = 1 2448 IPRPHP = MAX(0,IPRI6-4) 2449 CALL PHPGET(LSYM,NCONF,XNDXCI,FCAC,H2AC,DIAGL,ESHIFT, 2450 * IPWAY,IPRPHP,WRK,LWRK) 2451 END IF 2452 CALL SRWPHP(DIAGL,LPHPMX,.TRUE.) 2453C CALL SRWPHP(PHPVEC,LPHPMX,PHPWRT) 2454 CALL QEXIT('SIRPHP') 2455 RETURN 2456 END 2457C /* Deck srwphp */ 2458 SUBROUTINE SRWPHP(PHPVEC,LPHPMX,PHPWRT) 2459C 2460C 27-Oct-1990 Hans Joergen Aa. Jensen 2461C 2462C Sirius routine for Reading or Writing PHP information. 2463C 2464C Hessian DIAGONAL + PHP INFORMATION AS LAST RECORD ON LUIT2 2465C 2466C PHPWRT = .true. write information 2467C = .false. read information 2468C 2469#include "implicit.h" 2470#include "priunit.h" 2471 DIMENSION PHPVEC(*) 2472 LOGICAL PHPWRT 2473C 2474C Information used from common blocks; 2475C 2476C INFTAP : LUIT2 2477C 2478#include "inftap.h" 2479C 2480 LOGICAL FNDLAB 2481 CHARACTER*8 LAB123(3),LBPHP,LBEOD 2482 DATA LAB123(1),LBPHP,LBEOD/'********','PHPINFO ','EODATA '/ 2483C 2484 REWIND LUIT2 2485C 2486C read or write diagonal php information 2487C 2488 IF (PHPWRT) THEN 2489 CALL GETDAT(LAB123(2),LAB123(3)) 2490 IF( FNDLAB(LBPHP,LUIT2) ) THEN 2491 BACKSPACE LUIT2 2492 ELSE 2493 REWIND LUIT2 2494 IF( FNDLAB(LBEOD,LUIT2) ) BACKSPACE LUIT2 2495 END IF 2496 WRITE (LUIT2) LAB123,LBPHP 2497 LPHPM4 = MAX(4,LPHPMX) 2498 CALL WRITT(LUIT2,LPHPM4,PHPVEC) 2499 WRITE (LUIT2) LAB123,LBEOD 2500 ELSE 2501 IF (.NOT.FNDLAB(LBPHP,LUIT2)) THEN 2502 CALL QTRACE(LUPRI) 2503 CALL QUIT('ERROR SRWPHP: label PHPINFO not found on LUIT2') 2504 END IF 2505 CALL READT(LUIT2,LPHPMX,PHPVEC) 2506 END IF 2507C 2508C *** END OF SRWPHP 2509C 2510 RETURN 2511 END 2512C /* Deck upkcmo */ 2513 SUBROUTINE UPKCMO(CMO,UCMO) 2514C 2515C Copyright 6-May-1990 Hans Joergen Aa. Jensen 2516C 2517C Unpack CMO from symmetry packed form to UCMO 2518C 2519#include "implicit.h" 2520 DIMENSION CMO(*), UCMO(NBAST,NORBT) 2521C 2522C Used from common blocks: 2523C INFORB : NSYM,NBAST,NORBT,... 2524C 2525#include "inforb.h" 2526C 2527 CALL DZERO(UCMO,NBAST*NORBT) 2528 DO 400 ISYM = 1,NSYM 2529 IF (NORB(ISYM) .GT. 0) 2530 & CALL MCOPY(NBAS(ISYM),NORB(ISYM),CMO(ICMO(ISYM)+1),NBAS(ISYM), 2531 & UCMO(IBAS(ISYM)+1,IORB(ISYM)+1),NBAST) 2532C CALL MCOPY(NROWA,NCOLA,A,NRDIMA,B,NRDIMB) 2533 400 CONTINUE 2534 RETURN 2535 END 2536