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! 18 SUBROUTINE D0FCDIAG(D0DIAG,FCDIAG,FC) 19#include "implicit.h" 20#include "priunit.h" 21#include "infrsp.h" 22#include "inforb.h" 23 24 DIMENSION D0DIAG(*),FCDIAG(*),FC(*) 25 26 DO 300 ISYM=1, NSYM 27 IOFF=IIORB(ISYM)+1 28 29 DO 100 I=1, NISH(ISYM) 30 D0DIAG(IORB(ISYM)+I) =2.0D0 31 100 CONTINUE 32 33 34 DO 200 I=1,NORB(ISYM) 35 FCDIAG(IORB(ISYM)+I)=FC(IOFF) 36 IOFF = IOFF + I + 1 37 200 CONTINUE 38 300 CONTINUE 39 40 IF (IPRRSP.GE.10) THEN 41 CALL PRINTVEC('FCDIAG: ',1,NORBT,FCDIAG) 42 CALL PRINTVEC('D0DIAG: ',1,NORBT,D0DIAG) 43 44 END IF 45 RETURN 46 END 47C 48C END OF D0FCDIAG 49C 50 51 SUBROUTINE GETH1MO(H1MO,CMO,WRK,LWRK) 52#include "implicit.h" 53#include "priunit.h" 54#include "infrsp.h" 55#include "inforb.h" 56 57 DIMENSION H1MO(NORBT,NORBT),CMO(*),WRK(*) 58 59C 60C Calls sirh1 to get the one electron integrals 61C in AO basis (h1AO) and transforms the triangular 62C output to the full symmetrical matrix NORBT*NORBT 63C 64 65 KH1AO = 1 66 KH1AOT = KH1AO + N2BASX 67 68 KWRK1 = KH1AOT + NNORBX+1 69 LWRK1 = LWRK - KWRK1 70 71 CALL SIRH1(WRK(KH1AOT),WRK(KWRK1),LWRK1) 72 73 K = KH1AOT 74 KOFF = KH1AO 75 76 DO 300 ISYM=1,NSYM 77 DO 200 I=1,NBAS(ISYM) 78 DO 100 J=1,I 79 WRK(KOFF+(I-1)*NBAST+J-1)=WRK(K) 80 WRK(KOFF+(J-1)*NBAST+I-1)=WRK(K) 81 K=K+1 82 100 CONTINUE 83 200 CONTINUE 84 KOFF = KOFF + NBAS(ISYM)*NBAST + NBAS(ISYM) 85 300 CONTINUE 86 87 IF (IPRRSP .GE. 10) THEN 88 CALL PRINTMAT('H1AO ',1,NBAST,WRK(KH1AO)) 89 END IF 90 91C Transform the integrals to MO basis 92 93 CALL AO2MO(WRK(KH1AO),H1MO,CMO,WRK(KWRK1),LWRK1) 94 95 RETURN 96 END 97 98 SUBROUTINE GETESG_DENMAT_FOCMAT(D,F,NBAST,NNBAST) 99#include "implicit.h" 100#include "rspprp.h" 101#include "esg.h" 102#include "dummy.h" 103#include "priunit.h" 104 105 DIMENSION D(NNBAST), F(NNBAST) 106 107 WRITE (LUPRI,'(A,/)') 108 & 'Reading density and fock matrices in 1-el. part' 109 110 CALL FLSHFO(LUPRI) 111 112 READ (LUESG) D 113 READ (LUESG) F 114 115 CALL GPCLOSE(LUESG,'KEEP') 116 117 RETURN 118 END 119 120 SUBROUTINE GETESG_DENMAT(DMAT,NDMAT,WRK,LWRK) 121#include "implicit.h" 122#include "mxcent.h" 123#include "rspprp.h" 124#include "esg.h" 125#include "dummy.h" 126#include "inforb.h" 127#include "abainf.h" 128#include "priunit.h" 129 130 DIMENSION DMAT(NBAST,NBAST,NDMAT),WRK(*) 131 PARAMETER (MXDMAT=50) 132 DIMENSION ISYMDM(MXDMAT) 133 134 WRITE(LUPRI,'(A,/)') 'Reading density matrices in 2-el. part' 135 136 KDSO = 1 137 KWRK = KDSO + N2BASX 138 139 DO 100 I=1,NDMAT 140 ISYMDM(I) = 0 141 100 CONTINUE 142 ISYMDM(3) = ISYME - 1 143 ISYMDM(4) = ISYME - 1 144 145 IF (KWRK.GT.LWRK) CALL STOPIT('GETESG_DENMAT',' ',KWRK,LWRK) 146 147 CALL GPOPEN(LUESG2,'ESG_DMAT',' ',' ',' ',IDUMMY,.FALSE.) 148 REWIND (LUESG2) 149 150 DO IDMAT=2,NDMAT 151 CALL READT(LUESG2,N2BASX,WRK) 152 CALL DSOTAO(WRK,DMAT(1,1,IDMAT),NBAST,ISYMDM(IDMAT),IPRESG) 153 END DO 154 155 CALL GPCLOSE(LUESG2,'KEEP') 156 157 IF ( IPRESG .GE. 10 ) THEN 158 WRITE(LUPRI,*) 'Reading density matrices in 2-el. part' 159 CALL PRINTMAT('D0AO: ',1,NBAST,DMAT(1,1,1)) 160 CALL PRINTMAT('D2AO: ',1,NBAST,DMAT(1,1,2)) 161 CALL PRINTMAT('DXAO: ',1,NBAST,DMAT(1,1,3)) 162 CALL PRINTMAT('DXSAO: ',1,NBAST,DMAT(1,1,4)) 163 END IF 164 165 RETURN 166 END 167 168 SUBROUTINE INVE2VEC(CMO,UDV,PV,FC,FV,FCAC, 169 & H2AC,NSIM,INPVECS,OUTVECS,XINDX,WRK,LWRK) 170C 171C CALCULATES INV( E2 ) * VECTOR 172C 173#include "implicit.h" 174#include "dummy.h" 175#include "codata.h" 176#include "priunit.h" 177#include "infopt.h" 178#include "infrsp.h" 179#include "wrkrsp.h" 180#include "rspprp.h" 181#include "infpp.h" 182#include "inflr.h" 183#include "inforb.h" 184#include "infdim.h" 185#include "infpri.h" 186#include "inftap.h" 187#include "infsop.h" 188 189 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) 190 DIMENSION XINDX(*),WRK(*) 191 DOUBLE PRECISION INPVECS(*),OUTVECS(*) 192 CHARACTER*8 BLANK 193C 194 PARAMETER ( MAXSIM = 15, BLANK = ' ' ) 195 196 PARAMETER ( D2 = 2.0D0, D100 = 100.0D0, D0 = 0.0D0) 197 PARAMETER ( D2R3 = (2.0D0/3.0D0)) 198C 199C space allocation for reduced E(2) and reduced S(2) 200C 201 KREDE = 1 202 KREDS = KREDE + MAXRM*MAXRM 203 KIBTYP = KREDS + MAXRM*MAXRM 204 KEIVAL = KIBTYP + MAXRM 205 KRESID = KEIVAL + MAXRM 206 KEIVEC = KRESID + MAXRM 207 KREDGD = KEIVEC + MAXRM*MAXRM 208 KWRK1 = KREDGD + MAXRM*MAXRM 209 LWRK1 = LWRK + 1 - KWRK1 210 211 CALL DZERO(WRK,KWRK1) 212 213 IF (IPRPP .GT. 2 .OR. LWRK1 .LT. 2*KZYVAR) THEN 214 WRITE(LUPRI,'(/A)') ' --- IN INVE2VEC:' 215 WRITE(LUPRI,*)' THCPP,MAXRM ',THCPP,MAXRM 216 WRITE(LUPRI,*)' KSYMOP,NGPPP(KSYMOP) ',KSYMOP,NGPPP(KSYMOP) 217 WRITE(LUPRI,*)' LWRK ,LWRK1 ',LWRK,LWRK1 218 END IF 219C 220C ALLOCATE WORK SPACE FOR EIGENVECTORS AND TRANSITION MOMENTS 221C 222 LNEED = 100 + 2*KZYVAR 223C 224C MAXIMUM NUMBER OF SIMULTANEOUS SOLUTION VECTORS 225C 226C NSIM = MIN(KEXCNV, MAXSIM, (LWRK1-LNEED)/KZYVAR ) 227 IF (IPRPP .GT. 2 .OR. NSIM .LE. 0) THEN 228 LWRK2 = KWRK1 + LNEED + KZYVAR 229C ... need at least space for one KZYVAR (NSIM = 1) 230 WRITE (LUPRI,*) ' KEXCNV,NSIM,LWRK2 ',KEXCNV,NSIM,LWRK2 231 IF (NSIM.LE.0) CALL ERRWRK('RSPPP work space',-LWRK2,LWRK) 232 END IF 233 234 KWRK2 = KWRK1 + NSIM*KZYVAR 235 236 LWRK2 = LWRK - KWRK2 237 238 THCRSP = THCPP 239 IPRRSP = IPRPP 240 MAXIT = MAXITP 241 242 CALL DZERO(WRK(KEIVAL),MAXRM) 243 CALL DZERO(WRK(KEIVEC),MAXRM) 244 245 246 DO 1100 ISIM=1,NSIM 247 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC,.TRUE.,BLANK, 248 * BLANK,INPVECS((ISIM-1)*KZYVAR+1),WRK(KREDGD), 249 * WRK(KREDE),WRK(KREDS), 250 * WRK(KIBTYP),WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC), 251 * XINDX,WRK(KWRK1),LWRK1) 252 WRITE(LUPRI,*) 'Finished RSPCTL for ISIM=', ISIM 253 CALL RSPEVE(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC), 254 & OUTVECS((ISIM-1)*KZYVAR+1), 255 & WRK(KWRK1),NSIM,ISIM-1) 256 1100 CONTINUE 257 258 RETURN 259 END 260C 261C *** END OF INVE2VEC 262C 263 264 265 SUBROUTINE FOLD(AIN,AOUT,N) 266C 267C This subroutine folds the square matrix AIN into the triangular 268C matrix AOUT 269C 270#include "implicit.h" 271#include "priunit.h" 272 DIMENSION AIN(N,N), AOUT(*) 273 274 IJ = 0 275 DO I=1,N 276 DO J=1,I-1 277 IJ = IJ + 1 278 AOUT(IJ) = AIN(I,J) + AIN(J,I) 279 END DO 280 IJ = IJ + 1 281 AOUT(IJ) = AIN(I,I) 282 END DO 283 284 RETURN 285 END 286 287 SUBROUTINE FMAT2VEC(NSIM,ZYVEC,ZYMAT) 288#include "implicit.h" 289#include "priunit.h" 290 DIMENSION ZYVEC(KZYWOP,*), ZYMAT(NORBT,NORBT,*) 291C 292#include "maxorb.h" 293#include "maxash.h" 294#include "infvar.h" 295#include "inforb.h" 296#include "infind.h" 297#include "infrsp.h" 298#include "wrkrsp.h" 299C 300 PARAMETER ( D1 = 1.0D0 ) 301 CALL DZERO(ZYVEC,NSIM*KZYVAR) 302 DO 100 ISIM=1,NSIM 303 DO 200 IG=1,KZWOPT 304 I=JWOP(1,IG) 305 J=JWOP(2,IG) 306 ZYVEC(IG,ISIM)=ZYMAT(J,I,ISIM)-ZYMAT(I,J,ISIM) 307 ZYVEC(IG+KZWOPT,ISIM)=-ZYVEC(IG,ISIM) 308 200 CONTINUE 309 100 CONTINUE 310 311 RETURN 312 END 313 314 SUBROUTINE SYMMETRIZE(MINP,MOUT,NSIM,NORBT) 315#include "implicit.h" 316#include "priunit.h" 317 318 DIMENSION MINP(NORBT,NORBT,*),MOUT(NORBT,NORBT,*) 319 DOUBLE PRECISION MINP, MOUT 320 321 DO 300 ISIM=1,NSIM 322 DO 200 I=1,NORBT 323 DO 100 J=1,I 324 MOUT(I,J,ISIM) = 0.5D0*(MINP(I,J,ISIM) 325 & +MINP(J,I,ISIM)) 326 MOUT(J,I,ISIM) = 0.5D0*(MINP(I,J,ISIM) 327 & +MINP(J,I,ISIM)) 328 100 CONTINUE 329 200 CONTINUE 330 300 CONTINUE 331 332 RETURN 333 END 334 335C 336C Old AO2MO for inputs without the symmetry 337C 338C SUBROUTINE AO2MO(MAO,MMO,CMO,NORBT,WRK,LWRK) 339C#include "implicit.h" 340C#include "priunit.h" 341C 342C 343C Purpose: transforms a two-index matrix from 344C AO to MO basis (using WRK as scratch ) 345C 346C MMO = CMO^T * MAO * CMO 347C 348C 349C DIMENSION MAO(*),MMO(*),CMO(*),WRK(*) 350C 351C CALL DGEMM('T','N',NORBT,NORBT,NORBT,1.D0, 352C & CMO,NORBT,MAO,NORBT,0.D0, 353C & WRK,NORBT) 354C 355C CALL DGEMM('N','N',NORBT,NORBT,NORBT,1.D0, 356C & WRK,NORBT,CMO,NORBT,0.D0, 357C & MMO,NORBT) 358C 359C RETURN 360C END 361C 362C 363C End of AO2MO_OLD 364C 365 366 SUBROUTINE AO2MO(AAO,BMO,CMO,WRK,LWRK) 367#include "implicit.h" 368#include "priunit.h" 369#include "inforb.h" 370 371C 372C Purpose: transforms a two-index matrix from 373C AO to MO basis (using WRK as scratch ) 374C 375C BMO = CMO^T * AAO * CMO 376C 377C CMO is in rectangular symmetry blocks, 378C AAO, BMO are full without any symmetry reduction 379C 380 381 DIMENSION AAO(*),BMO(*),CMO(*),WRK(*) 382 383 CALL DZERO(WRK,N2ORBX) 384 CALL DZERO(BMO,N2ORBX) 385 386C CALL PRINTMAT('AO2MO: beg. ',1,NORBT,AAO) 387 388 ICOFF = 1 389 DO 100 ISYM=1,NSYM 390 IF ( NBAS(ISYM) .GT. 0 ) THEN 391C IOFF = IORB(ISYM)*NORBT+IORB(ISYM) +1 392 IOFF = IORB(ISYM) + 1 393 CALL DGEMM('T','N', 394C & NORB(ISYM),NORB(ISYM),NBAS(ISYM),1.D0, 395 & NORB(ISYM),NORBT,NBAS(ISYM),1.D0, 396 & CMO(ICOFF),NBAS(ISYM), 397 & AAO(IOFF),NORBT,1.0D0, 398 & WRK(IOFF),NORBT) 399 ICOFF = ICOFF + N2ORB(ISYM) 400 END IF 401 100 CONTINUE 402 403 ICOFF = 1 404 DO 200 ISYM=1,NSYM 405 IF ( NBAS(ISYM) .GT. 0 ) THEN 406C IOFF = IORB(ISYM)*NORBT+IORB(ISYM)+1 407 IOFF = IORB(ISYM)*NORBT+1 408 CALL DGEMM('N','N', 409C & NORB(ISYM),NORB(ISYM),NBAS(ISYM),1.0D0, 410 & NORBT,NORB(ISYM),NBAS(ISYM),1.0D0, 411 & WRK(IOFF),NORBT, 412 & CMO(ICOFF),NBAS(ISYM),1.0D0, 413 & BMO(IOFF),NORBT) 414 ICOFF = ICOFF + N2ORB(ISYM) 415 END IF 416 200 CONTINUE 417 418C CALL PRINTMAT('AO2MO: final. ',1,NORBT,BMO) 419 420 RETURN 421 END 422 423C 424C End of AO2MO 425C 426 SUBROUTINE MO2AO(AMO,AAO,CMO,WRK,LWRK) 427#include "implicit.h" 428#include "priunit.h" 429#include "inforb.h" 430 431C 432C Purpose: transforms a two-index matrix from 433C AO to MO basis (using WRK as scratch ) 434C 435C AAO = CMO * AMO * CMO^T 436C 437C AMO, AAO are full matrix without a symmetry reduction 438C 439C We make a loop in which we multiply CMO block 440C of certain symmetry with the MO matrix 441C 442 DIMENSION AAO(*),AMO(*),CMO(*),WRK(*) 443 INTEGER COFF 444 445 CALL DZERO(WRK,NORBT**2) 446 CALL DZERO(AAO,NORBT**2) 447 448 DO ISYM=1,NSYM 449 IF ( NBAS(ISYM) .GT. 0 ) THEN 450 CALL DGEMM('N','N', 451 & NBAS(ISYM),NORBT,NORB(ISYM),1.D0, 452 & CMO(ICMO(ISYM) + 1),NBAS(ISYM), 453 & AMO(IORB(ISYM) + 1),NORBT,0.D0, 454 & WRK(IBAS(ISYM)+1),NBAST) 455 END IF 456 END DO 457 458 DO ISYM=1,NSYM 459 IF ( NBAS(ISYM) .GT. 0 ) THEN 460 CALL DGEMM('N','T', 461 & NORBT,NBAS(ISYM),NORB(ISYM),1.D0, 462 & WRK(IORB(ISYM)*NBAST+1),NBAST, 463 & CMO(ICMO(ISYM) + 1),NBAS(ISYM),0.D0, 464 & AAO(IBAS(ISYM)*NBAST+1),NBAST) 465 END IF 466 END DO 467 468 RETURN 469 END 470 471C 472C This version was working for non-symmetry calc 473C 474C 475C SUBROUTINE MO2AO_OLD(MMO,MAO,CMO,NORBT,WRK,LWRK) 476C#include "implicit.h" 477C#include "priunit.h" 478C 479C 480C Purpose: transforms a two-index matrix from 481C AO to MO basis (using WRK as scratch ) 482C 483C MAO = CMO * MMO * CMO^T 484C 485C 486C DIMENSION MAO(*),MMO(*),CMO(*),WRK(*) 487C DOUBLE PRECISION MAO,MMO,CMO,WRK 488C 489C CALL DGEMM('N','N',NORBT,NORBT,NORBT,1.D0, 490C & CMO,NORBT,MMO,NORBT,0.D0, 491C & WRK,NORBT) 492C 493C CALL DGEMM('N','T',NORBT,NORBT,NORBT,1.D0, 494C & WRK,NORBT,CMO,NORBT,0.D0, 495C & MAO,NORBT) 496C 497C RETURN 498C END 499 500C 501C End of MO2AO 502C 503 504 SUBROUTINE ESG_PVAL(DODFT,HFXFAC,PVAL,DMAT,A,B,C,D,F) 505#include "implicit.h" 506#include "inforb.h" 507 INTEGER A, B, C, D 508 PARAMETER (DP25 = 0.25 D00, DP5 = 0.5 D00, 509 & D1 = 1.0 D00, D2 = 2.0 D00, ZERADD = 1.D-15, 510 & D0 = 0.0 D00, DP125 = 0.125D00, D3=3.0D0, D4=4.0D0) 511 DIMENSION PVAL(*),DMAT(NBAST,NBAST,*) 512 DOUBLE PRECISION PVAL,DMAT 513 LOGICAL DODFT 514 515 IF (.NOT. DODFT) THEN 516 FACTOR = 1.0D0 517 ELSE 518 FACTOR = HFXFAC 519 END IF 520 PVAL(1) = D2*F*(DMAT(A,B,2)*DMAT(D,C,1)+DMAT(A,B,1)*DMAT(D,C,2) 521 & - DP25*FACTOR*(DMAT(C,A,2)*DMAT(D,B,1)+DMAT(D,A,2)*DMAT(C,B,1)+ 522 & DMAT(C,A,1)*DMAT(D,B,2)+DMAT(D,A,1)*DMAT(C,B,2)) 523 & + ( DMAT(A,B,4)*DMAT(D,C,4) 524 & - DP125*FACTOR*(DMAT(A,D,3)*DMAT(B,C,3)+DMAT(D,A,3)*DMAT(C,B,3)+ 525 & DMAT(A,C,3)*DMAT(B,D,3)+DMAT(C,A,3)*DMAT(D,B,3)))) 526 527 RETURN 528 END 529 530C---------------------------------------- 531C TEMPORARY SUBROUTINES FOR DEBUGING: 532C---------------------------------------- 533 SUBROUTINE PRINTVEC( TEXT, NSIM, KZYVAR, V ) 534#include "implicit.h" 535#include "priunit.h" 536 537 DIMENSION V(*) 538 CHARACTER*15 TEXT 539 540 WRITE (LUPRI,'(A)') TEXT 541 542 DO 200 I=1,NSIM 543 WRITE( LUPRI, * ) 'VECTOR NUMBER : ', I 544 DO 100 J=1,KZYVAR 545 WRITE( LUPRI, 300 ) I, J, 546 & V( (I-1)*KZYVAR + J ) 547 100 CONTINUE 548 200 CONTINUE 549 550 300 FORMAT(3X,' V',I1,'(',I2,')=',F12.6) 551 552 RETURN 553 END 554C 555C *** END OF PRINTVEC 556C 557 SUBROUTINE PRINTVEC2( TEXT, NSIM, KZVAR, V ) 558#include "implicit.h" 559#include "priunit.h" 560 561 DIMENSION V(*) 562 CHARACTER*15 TEXT 563 564 WRITE (LUPRI,'(A)') TEXT 565 566 DO ISIM=1,NSIM 567 IF (NSIM.GT.1) WRITE( LUPRI, * ) 'MATRIX NUMBER : ', ISIM 568 CALL OUTPUT(V,1,KZVAR,1,2,KZVAR,2,1,LUPRI) 569 END DO 570 571c DO 200 I=1,NSIM 572c WRITE( LUPRI, * ) 'VECTOR NUMBER : ', I 573c DO 100 J=1,KZYVAR 574c WRITE( LUPRI, 300 ) I, J, 575c & V( (I-1)*KZYVAR + J ) 576c 100 CONTINUE 577c 200 CONTINUE 578 579c 300 FORMAT(3X,' V',I1,'(',I2,')=',F12.6) 580 581 RETURN 582 END 583C 584C *** END OF PRINTVEC 585C 586 587 588 SUBROUTINE PRINTMAT( TEXT, NSIM, NORBT, A ) 589#include "implicit.h" 590#include "priunit.h" 591 592 593 DIMENSION A(NORBT,NORBT,*) 594 CHARACTER*15 TEXT 595 596 WRITE (LUPRI,'(A)') TEXT 597 DO ISIM=1,NSIM 598 IF (NSIM.GT.1) WRITE( LUPRI, * ) 'MATRIX NUMBER : ', ISIM 599 CALL OUTPUT(A,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 600 END DO 601c DO 200 ISIM=1,NSIM 602c WRITE( LUPRI, * ) 'MATRIX NUMBER : ', ISIM 603c DO 100 I=1,NORBT 604c WRITE( LUPRI, 400 ) 605c & (A(I,J,ISIM),J=1,NORBT) 606c 100 CONTINUE 607c 200 CONTINUE 608c 609c 400 FORMAT(6(3X,F12.6)) 610 611 RETURN 612 END 613C 614C *** END OF PRINTVEC 615C 616