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 opgctl */ 20 SUBROUTINE OPGCTL(ORTOPG,RHSOPG,CMO,DV,SOINT,PRPACT,PRPICT,WORK, 21 & LWORK,IREP,NODC,NODV,IPRINT) 22C 23C tuh Dec 1989 24C 25C Control routine for calculation of reorthonormalization terms 26C (ORTOPG) and right-hand sides of response equation (RHSOPG) for 27C the molecular gradient of an arbitrary one-electron property. 28C Input are SO integrals SOINT of symmetry IREP and molecular 29C orbitals (CMO) and one-electron active density (DV). 30C 31#include "implicit.h" 32#include "mxcent.h" 33#include "nuclei.h" 34 DIMENSION ORTOPG(3*NUCDEP), RHSOPG(NVARPT), CMO(NCMOT), 35 * DV(NASHT,NASHT), SOINT(NBAST,NBAST), WORK(LWORK) 36 LOGICAL NODC, NODV 37#include "inforb.h" 38#include "inflin.h" 39C 40 CALL QENTER('OPGCTL') 41 KMOINT = 1 42 KFCKMO = KMOINT + NORBT*NORBT 43 KWRK = KFCKMO + NORBT*NORBT 44 LWRK = LWORK - KWRK + 1 45 IF (KWRK .GT. LWORK) CALL STOPIT('OPGCTL',' ',KWRK,LWORK) 46 CALL OPGCT1(ORTOPG,RHSOPG,CMO,DV,SOINT,WORK(KMOINT),WORK(KFCKMO), 47 & PRPACT,PRPICT,WORK(KWRK),LWRK,IREP,NODC,NODV, 48 & IPRINT) 49 CALL QEXIT('OPGCTL') 50 RETURN 51 END 52C /* Deck opgct1 */ 53 SUBROUTINE OPGCT1(ORTOPG,RHSOPG,CMO,DV,SOINT,DMOINT,FOCKMO,PRPACT, 54 & PRPICT,WORK,LWORK,IREP,NODC,NODV,IPRINT) 55C 56C December 1989, tuh 57C 58#include "implicit.h" 59#include "priunit.h" 60#include "mxcent.h" 61#include "nuclei.h" 62C 63 LOGICAL NODC, NODV 64 DIMENSION CMO(NCMOT), DV(NASHT,NASHT), 65 * SOINT(NBAST,NBAST), DMOINT(NORBT,NORBT), 66 * ORTOPG(3*NUCDEP), RHSOPG(NVARPT), 67 * FOCKMO(NORBT,NORBT), WORK(LWORK) 68C 69#include "abainf.h" 70#include "inforb.h" 71#include "inflin.h" 72C 73C *********************************** 74C ***** Fock matrix in MO basis ***** 75C *********************************** 76C 77 CALL MO1TRA(CMO,DMOINT,SOINT,WORK,LWORK,IREP,IPRINT) 78 CALL OPGFCK(DMOINT,FOCKMO,DV,IREP,NODC,NODV,IPRINT) 79C 80C ********************************************* 81C ***** Reorthonormalization contribution ***** 82C ********************************************* 83C 84 IF (DIPDER .OR. QPGRAD) THEN 85 CALL OPGORT(CMO,FOCKMO,ORTOPG,WORK,LWORK, 86 & IREP,IPRINT) 87 END IF 88C 89C *************************** 90C ***** Right-hand side ***** 91C *************************** 92C 93 CALL OPGRHS(RHSOPG,FOCKMO,DMOINT,PRPACT,PRPICT,WORK,LWORK, 94 & IPRINT) 95C 96 RETURN 97 END 98C /* Deck opgfck */ 99 SUBROUTINE OPGFCK(DMOINT,FOCK,DV,IREP,NODC,NODV,IPRINT) 100C 101C Purpose: Construction of property Fock matrix in MO basis 102C 103C tuh 131289 104C 105#include "implicit.h" 106#include "priunit.h" 107#include "maxorb.h" 108 PARAMETER(D2=2.0D0) 109C 110 INTEGER Q, U, U1, U2 111 LOGICAL NODC, NODV 112 DIMENSION DMOINT(NORBT,NORBT), FOCK(NORBT,NORBT), DV(NASHT,NASHT) 113C 114#include "inforb.h" 115C 116 ISYMPT = IREP + 1 117 IF (IPRINT .GE. 5) THEN 118 CALL TITLER('Output from OPGFCK','*',103) 119 WRITE (LUPRI,'(/A,I5)') ' ISYMPT ', ISYMPT 120 IF (IPRINT .GE. 10) THEN 121 CALL AROUND('Integrals in OPGFCK') 122 CALL OUTPUT(DMOINT,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 123 CALL AROUND('DV matrix in OPGFCK') 124 CALL OUTPUT(DV,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 125 END IF 126 END IF 127C 128 CALL DZERO(FOCK,N2ORBX) 129C 130C Inactive part 131C 132 IF (.NOT.NODC) THEN 133 DO 100 ISYMI = 1, NSYM 134 ISYMQ = MULD2H(ISYMI,ISYMPT) 135 NISHI = NISH(ISYMI) 136 NORBQ = NORB(ISYMQ) 137 IF (NISHI.GT.0 .AND. NORBQ.GT.0) THEN 138 ISTRI = IORB(ISYMI) + 1 139 IOFFQ = IORB(ISYMQ) 140 DO 110 Q = IOFFQ + 1, IOFFQ + NORBQ 141 CALL DCOPY(NISHI,DMOINT(ISTRI,Q),1,FOCK(ISTRI,Q),1) 142 CALL DSCAL(NISHI,D2,FOCK(ISTRI,Q),1) 143 110 CONTINUE 144 END IF 145 100 CONTINUE 146 END IF 147C 148C Active part 149C 150 IF (.NOT.NODV) THEN 151 DO 200 ISYMU = 1, NSYM 152 ISYMQ = MULD2H(ISYMU,ISYMPT) 153 NASHU = NASH(ISYMU) 154 NORBQ = NORB(ISYMQ) 155 IF (NASHU.GT.0 .AND. NORBQ.GT.0) THEN 156 IOFFU1 = IORB(ISYMU) + NISH(ISYMU) 157 IOFFU2 = IASH(ISYMU) 158 IOFFQ = IORB(ISYMQ) 159 DO 210 U = 1, NASHU 160 U1 = IOFFU1 + U 161 U2 = IOFFU2 + U 162 DO 220 Q = IOFFQ + 1, IOFFQ + NORBQ 163 FOCK(U1,Q) = DDOT(NASHU,DV (IOFFU2 + 1, U2),1, 164 * DMOINT(IOFFU1 + 1, Q ),1) 165 220 CONTINUE 166 210 CONTINUE 167 END IF 168 200 CONTINUE 169 END IF 170C 171 IF (IPRINT .GE. 7) THEN 172 CALL AROUND('Property Fock matrix in OPGFCK') 173 CALL OUTPUT(FOCK,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 174 END IF 175 RETURN 176 END 177C /* Deck opgort */ 178 SUBROUTINE OPGORT(CMO,FOCKMO,ORTOPG,WORK,LWORK,IREP,IPRINT) 179C 180C Dec 89, tuh 181C 182C This subroutine calculates the reorthonormalization contribution 183C to the one-electron property gradient by 184C 185C (1) transforming the property Fock matrices to (contravariant) 186C SO basis 187C (2) contracting these matrices with the differentiated overlap 188C matrices (SO basis) 189C 190#include "implicit.h" 191#include "priunit.h" 192#include "mxcent.h" 193#include "nuclei.h" 194 DIMENSION CMO(NCMOT), FOCKMO(NORBT,NORBT), WORK(LWORK), 195 & ORTOPG(3*NUCDEP) 196#include "inforb.h" 197C 198 CALL QENTER('OPGORT') 199 IF (IPRINT .GT. 5) CALL TITLER('Output from OPGORT','*',103) 200 KSDSP = 1 201 KSDSQ = KSDSP + 3*N2BASX 202 KSDSQT = KSDSQ + N2ORBX 203 KWRK = KSDSQT + N2ORBX 204 LWRK = LWORK - KWRK + 1 205 IF (KWRK .GE. LWORK) CALL STOPIT('OPGORT',' ',KWRK,LWORK) 206 CALL OPGOR1(CMO,FOCKMO,ORTOPG,WORK(KSDSP), 207 & WORK(KSDSQ),WORK(KWRK),LWRK,WORK(KSDSQT), 208 & IREP,IPRINT) 209 CALL QEXIT('OPGORT') 210 RETURN 211 END 212C /* Deck opgor1 */ 213 SUBROUTINE OPGOR1(CMO,FOCKMO,ORTOPG,SDSP,SDSQ,WORK,LWORK, 214 & SDSQT,IREP,IPRINT) 215#include "implicit.h" 216#include "priunit.h" 217#include "maxaqn.h" 218#include "maxorb.h" 219#include "mxcent.h" 220#include "nuclei.h" 221 PARAMETER (D2 = 2.0D0) 222 LOGICAL DERIV(3), DOTD 223 DIMENSION CMO(NCMOT), FOCKMO(NORBT,NORBT), ORTOPG(3*NUCDEP), 224 * SDSP(N2BASX,3), 225 * SDSQ(NORBT,NORBT), SDSQT(NORBT,NORBT), WORK(LWORK) 226 CHARACTER*4 KEY 227#include "abainf.h" 228#include "symmet.h" 229#include "inforb.h" 230C 231 DATA DERIV /3*.TRUE./ 232C 233 DOTD = .FALSE. 234 IF (IPRINT .GT. 10) THEN 235 CALL HEADER('MO Fock matrix in OPGOR1',-1) 236 CALL OUTPUT(FOCKMO,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 237 END IF 238C 239C ***** Calculate reorthonormalization terms ***** 240C 241 CALL DZERO(ORTOPG,3*NUCDEP) 242 DO 100 IATOM = 1, NUCIND 243 IF (IPRINT .GT. 5) WRITE (LUPRI,'(/A,I5)') ' IATOM ', IATOM 244 KEY = 'DMAT' 245 IF (NODIFC) KEY = 'OMAT' 246 CALL GETSD(DERIV,CMO,SDSP,WORK,LWORK,IATOM,.FALSE., 247 * KEY,DOTD,IPRINT,IPRINT) 248 DO 200 ICOOR = 1, 3 249 IF (IPRINT .GT. 5) WRITE (LUPRI,'(/A,I5)') ' ICOOR ', ICOOR 250 ISCOOR = IPTCNT(3*(IATOM - 1) + ICOOR,IREP,1) 251 IF (ISCOOR .GT. 0) THEN 252 IF (IPRINT.GT.5) WRITE (LUPRI,'(/A,I5)') ' ISCOOR',ISCOOR 253 CALL EXFDTD(SDSP(1,ICOOR),SDSQ,IREP+1,IPRINT) 254 CALL TRPMAT(SDSQ,NORBT,NORBT,SDSQT) 255 ORTOPG(ISCOOR) = - DDOT(N2ORBX,SDSQT,1,FOCKMO,1) 256 IF (.NOT. NODIFC) ORTOPG(ISCOOR) = D2*ORTOPG(ISCOOR) 257 END IF 258 200 CONTINUE 259 100 CONTINUE 260C 261C ***** Print ***** 262C 263 IF (IPRINT .GT. 3) THEN 264 CALL HEADER('ORTOPG in OPGOR1',-1) 265 CALL OUTPUT(ORTOPG,1,1,1,3*NUCDEP,1,3*NUCDEP,1,LUPRI) 266 END IF 267 RETURN 268 END 269C /* Deck opgrhs */ 270 SUBROUTINE OPGRHS(RHSOPG,FOCKMO,DMOINT,PRPACT,PRPICT,WORK,LWORK, 271 & IPRINT) 272C 273C Jan 1990 tuh 274C 275C Calculate right-hand side for response equations for 276C one-electron properties 277C 278! module dependencies 279 use lucita_mcscf_ci_cfg 280#include "implicit.h" 281#include "priunit.h" 282#include "iratdef.h" 283#include "mxcent.h" 284#include "maxorb.h" 285 PARAMETER (D0 = 0.0D0, D2 = 2.0D0) 286C 287 DIMENSION RHSOPG(NVARPT), FOCKMO(NORBT,NORBT), 288 & DMOINT(NORBT,NORBT), WORK(LWORK) 289C 290#include "inftap.h" 291#include "infinp.h" 292#include "inforb.h" 293#include "infdim.h" 294#include "inflin.h" 295C 296C 297 CALL QENTER('OPGRHS') 298 CALL TIMER('START ',TIME1,TIME2) 299 IF (IPRINT .GT. 4) CALL TITLER('Output from OPGRHS','*',103) 300C 301C ******************************************* 302C ***** Construction of right-hand side ***** 303C ******************************************* 304C 305 CALL DZERO(RHSOPG,NVARPT) 306C 307C (A) Construct orbital part of right-hand side 308C --------------------------------------------- 309C 310 CALL OPGORB(FOCKMO,RHSOPG(NCONST + 1)) 311C 312C Average rotation gradients if SUPSYM and sym 1: 313C 314 IF (LSYMPT.EQ.1) CALL AVERAG(RHSOPG(NCONST + 1),NWOPPT,1) 315C 316C (B) Construct configuration part of right-hand side 317C --------------------------------------------------- 318C 319 PRPACT = D0 320 IF (NCONST .GT. 1) THEN 321C 322C Work space allocation: 323C 324 KCREF = 1 325 KACAC = KCREF + NCONRF 326 KACACP = KACAC + NASHT*NASHT 327 KCINDX = KACACP + NASHT*NASHT 328 KWRK = KCINDX + LCINDX 329 LWRK = LWORK - KWRK + 1 330 IF (KWRK .GE. LWORK) CALL STOPIT('OPGRHS','GETCIX',KWRK,LWORK) 331C 332C CI index vector: 333C 334 CALL GETCIX(WORK(KCINDX),LSYMRF,LSYMST,WORK(KWRK),LWRK,0) 335C 336C Active part of property integrals: 337C 338 CALL GETAC1(DMOINT,WORK(KACAC)) 339C 340 IF (NCONRF .GT. 0) THEN 341 REWIND LUSIFC 342 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 343 READ (LUSIFC) 344 READ (LUSIFC) 345 READ (LUSIFC) 346 CALL READI(LUSIFC,IRAT*NCONRF,WORK(KCREF)) 347 END IF 348C 349 if(ci_program .eq. 'SIRIUS-CI')then 350 CALL DCOPY(NASHT**2,WORK(KACAC),1,WORK(KACACP),1) 351 else if(ci_program .eq. 'LUCITA ')then 352 cref_is_active_bvec_for_sigma = .true. 353 !... pack matrix in lower triangular form for lucita 354 CALL DSITSP(NASHT,WORK(KACAC),WORK(KACACP)) 355 end if 356 357 CALL CIPRP(1,WORK(KCREF),RHSOPG,NVARPT,WORK(KACACP), 358 * WORK(KCINDX),WORK(KWRK),LWRK) 359C CALL CIPRP(NSIM,CREF,SCVECS,LSCVEC,PRPAC,CINDEX,WORK,LFREE) 360C 361 IF (LSYMPT .EQ. 1) THEN 362 PRPACT = DDOT(NCONST,WORK(KCREF),1,RHSOPG,1) 363 CALL DAXPY(NCONST,(-PRPACT),WORK(KCREF),1,RHSOPG,1) 364 ELSE 365 PRPACT = D0 366 END IF 367 CALL DSCAL(NCONST,D2,RHSOPG,1) 368 END IF 369C 370C ***** Print right-hand side ***** 371C 372 IF (IPRINT .GE. 5) THEN 373 CALL HEADER('Orbital part of right-hand side',-1) 374 IF (IPRINT .GT. 20) THEN 375 PRFAC = 0.0D0 376 ELSE 377 PRFAC = 0.1D0 378 END IF 379 CALL PRKAP(NWOPPT,RHSOPG(NCONST + 1),PRFAC,LUPRI) 380 IF (IPRINT .GT. 20) THEN 381 CALL HEADER('Configuration part of right-hand side',-1) 382 CALL OUTPUT(RHSOPG,1,1,1,NCONST,1,NCONST,1,LUPRI) 383 END IF 384 END IF 385C 386C One-electron property calculated from CI part of RHS 387C 388 IF (LSYMPT .EQ. 1) THEN 389 PRPICT = D0 390 DO 100 ISYM = 1, NSYM 391 IORB1 = IORB(ISYM) + 1 392 NISHI = NISH(ISYM) 393 IF (NISHI .GT. 0) THEN 394 PRPICT = PRPICT + DSUM(NISHI,DMOINT(IORB1,IORB1),NORBT+1) 395 END IF 396 100 CONTINUE 397 PRPICT = D2*PRPICT 398 PRPCI = PRPACT + PRPICT 399 IF (IPRINT .GT. 5) THEN 400 WRITE (LUPRI,'(/,A,D24.12)') ' PRPICT ', PRPICT 401 WRITE (LUPRI,'( A,D24.12)') ' PRPACT ', PRPACT 402 WRITE (LUPRI,'( A,D24.12)') ' PRPCI ', PRPCI 403 END IF 404 ELSE 405 PRPICT = D0 406 PRPACT = D0 407 END IF 408C 409 IF (IPRINT .GT. 1) CALL TIMER('OPGRHS',TIME1,TIME2) 410 CALL QEXIT('OPGRHS') 411 RETURN 412C 413C end of OPGRHS 414C 415 END 416C /* Deck opgorb */ 417 SUBROUTINE OPGORB(FOCKMO,ORBGRD) 418C 419C Dec-1989 hjaaj & tuh 420C 421C Purpose: 422C To add the orbital property gradients in ORBGRD 423C 424C Input: 425C The total property Fock matrix FOCKMO 426C 427C Output: 428C Property orbital gradients in ORBGRD 429C 430#include "implicit.h" 431#include "priunit.h" 432#include "maxorb.h" 433#include "infvar.h" 434 PARAMETER (D2 = 2.0D0) 435 DIMENSION FOCKMO(NORBT,NORBT), ORBGRD(NWOPT) 436C 437C Used from common blocks: 438C INFVAR : NWOPT,JWOP(2,*) 439C INFORB : NORBT 440C 441#include "inforb.h" 442C 443 DO 100 IG = 1, NWOPT 444 K = JWOP(1,IG) 445 L = JWOP(2,IG) 446 ORBGRD(IG) = ORBGRD(IG) + D2*(FOCKMO(K,L) - FOCKMO(L,K)) 447 100 CONTINUE 448 RETURN 449 END 450C /* Deck mo1tra */ 451 SUBROUTINE MO1TRA(CMO,SQMO,SQSO,WORK,LWORK,IREPO,IPRINT) 452C 453C This routine transforms derivative SO matrices to MO basis. 454C 455C tuh 010190 456C 457#include "implicit.h" 458#include "priunit.h" 459#include "maxorb.h" 460#include "inforb.h" 461 DIMENSION CMO(NCMOT), WORK(LWORK), SQMO(NORBT,NORBT), SQSO(N2BASX) 462 PARAMETER (D1 =1.0D0) 463C 464C ***** Print Section ***** 465C 466 IF (IPRINT .GT. 5) THEN 467 CALL HEADER('Output from MO1TRA',-1) 468 WRITE (LUPRI,'(A,I5)') ' IREPO ', IREPO 469 WRITE (LUPRI,'(A,I5)') ' NORBT ', NORBT 470 WRITE (LUPRI,'(A,I5)') ' NBAST ', NBAST 471 WRITE (LUPRI,'(A,I5)') ' NCMOT ', NCMOT 472 WRITE (LUPRI,'(A,I10)') ' LWORK ', LWORK 473 IF (IPRINT .GT. 15) THEN 474 CALL HEADER('Square SO matrix',-1) 475 CALL OUTPUT(SQSO,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 476 END IF 477 END IF 478C 479C ***** Transform matrix from SO to MO basis ***** 480C 481 CALL DZERO(SQMO,N2ORBX) 482C 483C Loop over irreps for orbitals 484C 485 ISYMO = IREPO + 1 486 DO 100 ISYMA = 1, NSYM 487 IF (NORB(ISYMA) .GT. 0) THEN 488 ISYMB = MULD2H(ISYMA,ISYMO) 489 IF ((ISYMA.GE.ISYMB) .AND. (NORB(ISYMB).GT.0)) THEN 490C 491C Print symmetries and orbitals 492C 493 IF (IPRINT.GT.15) THEN 494 WRITE (LUPRI,'(A,3I3)') ' ISYMA/B/O',ISYMA,ISYMB,ISYMO 495 WRITE (LUPRI,'(/A,I5,A)') 496 * ' MO coefficients for symmetry', ISYMA 497 CALL OUTPUT(CMO(ICMO(ISYMA)+1),1,NBAS(ISYMA),1, 498 * NORB(ISYMA),NBAS(ISYMA),NORB(ISYMA),1,LUPRI) 499 WRITE (LUPRI,'(/A,I5,A)') 500 * ' MO coefficients for symmetry', ISYMB 501 CALL OUTPUT(CMO(ICMO(ISYMB)+1),1,NBAS(ISYMB),1, 502 * NORB(ISYMB),NBAS(ISYMB),NORB(ISYMB),1,LUPRI) 503 END IF 504C 505C Transform matrix block(s) 506C 507 CALL UTHV(CMO(ICMO(ISYMA)+1),SQSO,CMO(ICMO(ISYMB)+1), 508 * ISYMA,ISYMB,NBAS(ISYMA),NBAS(ISYMB),SQMO,WORK) 509C 510C Print transformed matrix thus far 511C 512 IF (IPRINT.GT.25) THEN 513 WRITE(LUPRI,'(/4A)')' Unfinished matrix in MO basis' 514 CALL OUTPUT(SQMO,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 515 END IF 516 END IF 517 END IF 518 100 CONTINUE 519 IF (ISYMO .GT. 1) CALL TRANSX(SQMO,SQMO,NORBT,NORBT,D1,IPRINT) 520 IF (IPRINT.GT.15) THEN 521 CALL HEADER('Matrix in MO basis',-1) 522 CALL OUTPUT(SQMO,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 523 END IF 524 RETURN 525 END 526