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 CC_EXGR 20 SUBROUTINE CC_EXGR(WORK,LWORK) 21C 22C---------------------------------------------------------------------- 23C 24C Purpose: Direct calculation of Coupled Cluster analytical 25C excited state first order properties and gradients. 26C 27C CIS, CCS, CC2, CCSD 28C 29C Solves for excited state t-bar amplitudes 30C = Lagrangian multipliers. 31C Calculates first order properties: dipole moment, 32C quadrupole moment, electric field gradients, 33C relativistic corrections, electronic moments. 34C 35C Written by Ove Christiansen April 1997. 36C 37C--------------------------------------------------------------------- 38C 39#include "implicit.h" 40#include "priunit.h" 41#include "maxorb.h" 42#include "mxcent.h" 43 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, IZERO = 0 , TWO = 2.0D0) 44#include "codata.h" 45#include "dummy.h" 46#include "iratdef.h" 47#include "ccfop.h" 48#include "ccexgr.h" 49#include "ccexci.h" 50#include "cclr.h" 51#include "ccorb.h" 52#include "ccsdsym.h" 53#include "ccsdio.h" 54#include "ccsdinp.h" 55#include "ccsections.h" 56#include "ccfield.h" 57#include "ccroper.h" 58#include "ccfro.h" 59#include "exeinf.h" 60#include "infvar.h" 61#include "dipole.h" 62#include "quadru.h" 63#include "nqcc.h" 64C 65 LOGICAL LINQCC, LPROJECT, TRIPLET, LREDS,B0SKIP 66 DIMENSION WORK(LWORK), ELSEMO(3,3), SKODE(3,3), SKODN(3,3) 67 CHARACTER*8 FC1AM,FC2AM,FRHO1,FRHO2,FRHO12,FC12AM,FS12AM,FS2AM 68 CHARACTER*3 APROXR12 69 PARAMETER (FC1AM ='CCR_C1AM',FC2AM ='CCR_C2AM') 70 PARAMETER (FRHO1 ='CCR_RHO1',FRHO2 ='CCR_RHO2') 71 PARAMETER (FRHO12='CCR_RH12',FC12AM='CCR_C12M',FS12AM='CCR_S12M') 72 PARAMETER (FS2AM ='CCR_S2DM') 73 CHARACTER MODEL*10 74 CHARACTER MODELPRI*4, MODELPRI2*30 75 CHARACTER LABEL*8, LABR12*4 76 CHARACTER*2 LIST 77 CHARACTER*5 FNDPTIA, FNDPTAB, FNDPTIJ 78 CHARACTER*6 FNDPTIA2 79 PARAMETER(FNDPTIA='DPTIA', FNDPTIA2 = 'DPTIA2', 80 * FNDPTAB='DPTAB' ,FNDPTIJ = 'DPTIJ' ) 81C 82#include "leinf.h" 83C 84 CALL QENTER('CC_EXGR') 85C 86C------------------------------------ 87C Header of Property calculation. 88C------------------------------------ 89C 90C CALL TIMER('START ',TIMEIN,TIMOUT) 91C 92 LUFC1 = -1 93 LUFC2 = -1 94 LUFC12 = -1 95 LUFR1 = -1 96 LUFR2 = -1 97 LUFR12 = -1 98 LUFS12 = -1 99 LUFS2 = -1 100 101 IF (CCR12) THEN 102 LABR12 = '-R12' 103 ELSE 104 LABR12 = ' ' 105 END IF 106C 107 IF (CCR12) CALL QUIT('R12 not yet implemented in CC_EXGR.') 108C 109 WRITE (LUPRI,'(1X,A,/)') ' ' 110 WRITE (LUPRI,'(1X,A)') 111 *'*********************************************************'// 112 *'**********' 113 WRITE (LUPRI,'(1X,A)') 114 *'* '// 115 *' *' 116 WRITE (LUPRI,'(1X,A)') 117 *'*---- OUTPUT FROM COUPLED CLUSTER RESPONSE ----'// 118 *'---------*' 119 IF ( CCFOP ) THEN 120 WRITE (LUPRI,'(1X,A)') 121 * '* '// 122 * ' *' 123 WRITE (LUPRI,'(1X,A)') 124 * '*----- CALCULATION OF EXCITED STATE FIRST ORDER PROPERTI'// 125 * 'ES ------*' 126 ENDIF 127 WRITE (LUPRI,'(1X,A)') 128 *'* '// 129 *' *' 130 WRITE (LUPRI,'(1X,A,/)') 131 *'*********************************************************'// 132 *'**********' 133C 134 MODEL = 'CCSD ' 135 IF (CC2) THEN 136 CALL AROUND('Coupled Cluster model is: CC2'//LABR12) 137 MODEL = 'CC2 ' 138 MODELPRI = ' CC2' 139 ENDIF 140 IF (CCS.AND.(.NOT.CIS)) THEN 141 CALL AROUND('Coupled Cluster model is: CCS') 142 MODEL = 'CCS ' 143 MODELPRI = ' CCS' 144 ENDIF 145 IF (CCS.AND.CIS) THEN 146 CALL AROUND('Model is not CC but CIS crap.') 147 MODEL = 'CCS ' 148 MODELPRI = ' CIS' 149 ENDIF 150 IF (CCD) THEN 151 CALL AROUND('Coupled Cluster model is: CCD'//LABR12) 152 MODEL = 'CCD ' 153 MODELPRI = ' CCD' 154 ENDIF 155 IF (CC3 ) THEN 156 CALL AROUND('Coupled Cluster model is: CC3'//LABR12) 157 MODEL = 'CC3 ' 158 MODELPRI = ' CC3' 159 WRITE(LUPRI,*) 160 * 'CC3 X-state first order properties not implemented' 161 RETURN 162 ENDIF 163 IF (CC1A) THEN 164 CALL AROUND('Coupled Cluster model is: CCSDT-1a'//LABR12) 165 MODEL = 'CCSDT-1a ' 166 WRITE(LUPRI,*) 167 * 'CCSDT-1a X-state first order properties not implemented' 168 RETURN 169 ENDIF 170 IF (CC1B) THEN 171 CALL AROUND('Coupled Cluster model is: CCSDT-1b'//LABR12) 172 MODEL = 'CCSDT-1b ' 173 WRITE(LUPRI,*) 174 * 'CCSDT-1b X-state first order properties not implemented' 175 RETURN 176 ENDIF 177 IF (CCSD) THEN 178 CALL AROUND('Coupled Cluster model is: CCSD'//LABR12) 179 MODEL = 'CCSD ' 180 MODELPRI = 'CCSD' 181 ENDIF 182C 183 NSIDIN = NSIDE 184C 185 LREDS = CCR12 186C 187 IF (IPRINT.GT.10) WRITE(LUPRI,*) 'CC_EXGR-1: Workspace:',LWORK 188C 189C----------------------------- 190C Initialize Variables. 191C----------------------------- 192C 193 ISYM = ISYMOP 194 ISYMTR = ISYMOP 195 IF (CCS) THEN 196 NCCVAR = NT1AM(ISYM) 197 ELSE 198 NCCVAR = NT1AM(ISYM) + NT2AM(ISYM) 199 IF (CCR12) NCCVAR = NCCVAR + NTR12AM(ISYM) 200 END IF 201 TRIPLET = .FALSE. 202 CALL CCLR_LEINFI(TRIPLET) 203 THRLE = THRLEQ 204 LINQCC = .TRUE. 205 NLOAD = 1 206 ISIDE = -1 207 LIST = 'E0' 208 IF (NCCVAR .LE. 0) THEN 209 WRITE(LUPRI,*) 'There are no amplitudes of this symmetry' 210 CALL QUIT('There are no amplitudes of this symmetry') 211 ENDIF 212C 213C------------------------------------------- 214C For CIS jump solution for multipliers. 215C------------------------------------------- 216C 217 IF (CCS.AND.CIS) GOTO 47 218C 219C----------------------------------------------- 220C Calculate appropriate F/B-transformations. 221C----------------------------------------------- 222C 223c IF ((.NOT.E0SKIP).AND.(.NOT.B0SKIP)) THEN 224c hjaaj jan-07: B0SKIP is not defined 225 IF (.NOT.E0SKIP) THEN 226 CALL CC_FRE(WORK,LWORK) 227 ENDIF 228C 229C----------------------------------------------------------------- 230C Start loop over lists over right hand lists to be solved. 231C Find out how many from common block info. 232C----------------------------------------------------------------- 233C 234 IBOFF = 0 235 IBEND = NXGRST 236 NEQSYM = IBEND - IBOFF 237 IF (E0SKIP) THEN 238 NEQSYM = 0 239 WRITE(LUPRI,*) 'Skipping solving for excited state t(0)' 240 ENDIF 241C 242C------------------------------------------------------------------ 243C Make bathching over nr. of simulatneous vectors to be solved. 244C Find out how many simultaneously; 245C all or take the number from input. 246C------------------------------------------------------------------ 247C 248 IF (NEQSYM .EQ. 0 ) THEN 249 NSIM = 0 250 NBAT = 0 251 ELSE 252 IF ( NSIMLE .EQ. 0 ) THEN 253 NSIM = NEQSYM 254 ELSE 255 NSIM = NSIMLE 256 ENDIF 257 NBAT = (NEQSYM-1)/NSIM + 1 258 ENDIF 259C 260 IF (DEBUG ) THEN 261 WRITE(LUPRI,*) ' ' 262 WRITE(LUPRI,*) 'CC_EXGR NEQSYM = ',NEQSYM 263 WRITE(LUPRI,*) 'CC_EXGR IBOFF = ',IBOFF 264 WRITE(LUPRI,*) 'CC_EXGR IBEND = ',IBEND 265 WRITE(LUPRI,*) 'CC_EXGR NSIM = ',NSIM 266 WRITE(LUPRI,*) 'CC_EXGR NBAT = ',NBAT 267 ENDIF 268C 269 DO 2000 IRB = 1, NBAT 270C 271C------------------------------------------------- 272C Find start number for this batch on list. 273C and the number of equations. 274C------------------------------------------------- 275C 276 IRST = IBOFF + (IRB-1)*NSIM + 1 277 NSIMR = MIN(NSIM,NEQSYM - (IRB-1)*NSIM) 278 IREND = IRST + NSIMR - 1 279C 280 IF (DEBUG ) THEN 281 WRITE(LUPRI,*) 'CC_EXGR IRST = ',IRST 282 WRITE(LUPRI,*) 'CC_EXGR IREND = ',IREND 283 WRITE(LUPRI,*) 'CC_EXGR NSIMR = ',NSIMR 284 ENDIF 285C 286 IF (IPRINT .GT. 2) THEN 287 WRITE(LUPRI,'(/,1x,A,I2,A)') 288 * 'Solving ',NSIMR,' response equations for states;' 289 DO 2005 IR = IRST,IREND 290 IEX = IXGRST(IR) 291 WRITE(LUPRI,'(1X,A,F10.6,A,I3,A,I3)') 292 * 'Energy:',EIGVAL(IEX), 293 * ' and nr.: ',IEX-ISYOFE(ISYEXC(IEX)),' of symmetry: ', 294 * ISYEXC(IEX) 295 2005 CONTINUE 296 ENDIF 297C 298 IF (IPRINT.GT.10) 299 * WRITE(LUPRI,*) 'CC_EXGR Workspace:',LWORK 300C 301 CALL FLSHFO(LUPRI) 302C 303C--------------------------------- 304C Allocation of work space. 305C--------------------------------- 306C 307 KIPLAC = 1 308 KREDH = KIPLAC + MAXRED 309 KREDGD = KREDH + MAXRED*MAXRED 310 KEIVAL = KREDGD + MAXRED 311 KSOLEQ = KEIVAL + MAXRED 312 KWRK1 = KSOLEQ + MAXRED*MAXRED 313 IF (LREDS) THEN 314 KREDS = KWRK1 315 KWRK1 = KREDS + MAXRED*MAXRED 316 END IF 317 LWRK1 = LWORK - KWRK1 318C 319 IF (LWRK1.LT. 0 ) 320 * CALL QUIT(' TOO LITTLE WORKSPACE IN CC_RSPSOL') 321C 322C------------------------- 323C Zero frequencies. 324C------------------------- 325C 326 CALL DZERO(WORK(KEIVAL),MAXRED) 327C 328C------------------- 329C Open files. 330C------------------- 331C 332 CALL CC_FILOP(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 333 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 334 * FS12AM,LUFS12,FS2AM,LUFS2) 335C 336C----------------------------- 337C Create start vectors. 338C----------------------------- 339C 340 TRIPLET = .FALSE. 341 LPROJECT = .FALSE. 342 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 343 * FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ, 344 * TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 345 * NREDH,WORK(KEIVAL),WORK(KIPLAC), 346 * WORK(KWRK1),LWRK1,LIST) 347C 348C------------------------------------------ 349C Solve equations by call to solver. 350C------------------------------------------ 351C 352 write(LUPRI,*) 'exgr, irst:',irst 353 ECURR = 0.0D0 354 CALL CCEQ_SOL(LIST,LPROJECT,ISTATPRJ,ECURR, 355 * FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 356 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 357 * LREDS,WORK(KREDS),FS12AM,LUFS12,FS2AM,LUFS2, 358 * LINQCC,TRIPLET,ISIDE,IRST,NSIMR,NUPVEC, 359 * NREDH,WORK(KREDH),WORK(KEIVAL), 360 * WORK(KSOLEQ),WORK(KWRK1),LWRK1,APROXR12) 361C 362C---------------------------------- 363C Analysis of solution vectors. 364C---------------------------------- 365C 366 NVARPT= LETOT + 2*NALLAI(ISYM) 367 KWRK2 = KWRK1 + NVARPT 368 LWRK2 = LWORK - KWRK2 369C 370 THRESH = 0.05 371 MAXLIN = 100 372 NSIMUL = MIN(NSIM,LWRK2/LETOT ) 373 NBATCH = (NSIM-1)/NSIMUL + 1 374 IOFF1 = 1 375 ICOUNT = 0 376 ILSTNR = IRST 377C 378 DO 500 I = 1,NBATCH 379 IOFF2 = MIN(NSIMUL,NSIM - (I-1)*NSIMUL) 380 CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 381 * TRIPLET,CCR12,NREDH,IOFF1,IOFF2,WORK(KSOLEQ), 382 * WORK(KWRK2),WORK(KWRK1)) 383C 384 IF ( IPRINT .GT. 10 ) THEN 385 RHO1N = DDOT(NT1AM(ISYM),WORK(KWRK2),1,WORK(KWRK2),1) 386 RHO2N = DDOT(NT2AM(ISYM),WORK(KWRK2+NT1AMX),1, 387 * WORK(KWRK2+NT1AMX),1) 388 WRITE(LUPRI,*) 'Norm of Lambda vector :',RHO1N 389 WRITE(LUPRI,*) 'Norm of Lambda vector :',RHO2N 390 ENDIF 391C 392 IF ( IPRINT .GT. 30 ) THEN 393 CALL AROUND('CC_EXGR: E0 vector in mo basis' ) 394 CALL OUTPUT(WORK(KWRK2),1,LETOT,1,NSIM, 395 * LETOT,NSIM,1,LUPRI) 396 ENDIF 397C 398 DO 510 J = 1,IOFF2 399 ICOUNT = ICOUNT + 1 400 IF (IPRINT .GT. 1) THEN 401 WRITE(LUPRI,'(//1X,A)') 402 *'Analysis of the X-st 0-th order Lagrangian multipliers: ' 403 WRITE(LUPRI,'(1X,A)') 404 *'--------------------------------------------------------' 405 CALL CC_PRAM(WORK(KWRK2 + (ICOUNT-1)*LETOT), 406 * PT1LOCAL,ISYM,.FALSE.) 407 ENDIF 408 IF (CCSTST) THEN 409 CALL DZERO(WORK(KWRK2+(ICOUNT-1)*LETOT+ 410 * NT1AM(ISYM)),NT2AMX) 411 ENDIF 412C 413C----------------------------------------- 414C Save response vectors on file. 415C----------------------------------------- 416C 417 KT1 = KWRK2 + (ICOUNT-1)*LETOT 418 KT2 = KWRK2 + (ICOUNT-1)*LETOT + NT1AM(ISYM) 419 IOPT = 3 420 CALL CC_WRRSP('E0',ILSTNR,ISYM,IOPT,MODEL,DUMMY, 421 * WORK(KT1),WORK(KT2),WORK(KWRK1),NVARPT) 422 ILSTNR = ILSTNR + 1 423C 424 510 CONTINUE 425 IOFF1 = IOFF1 + NSIMUL 426 500 CONTINUE 427C 428C------------------------------- 429C Close and delete files. 430C------------------------------- 431C 432 CALL CC_FILCL(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 433 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 434 * FS12AM,LUFS12,FS2AM,LUFS2) 435 CALL FLSHFO(LUPRI) 436C 437 2000 CONTINUE 438C 439C-------------------------------- 440C Spaghetti goto end for CIS. 441C-------------------------------- 442C 443 47 CONTINUE 444C 445C----------------------------------------------------------------------- 446C Calculate one electron AO-density and thereby CC nat.occ.num. 447C----------------------------------------------------------------------- 448C 449 DO 3000 IX = 1, NXGRST 450C 451 IEX = IXGRST(IX) 452 ISYMX = ISYEXC(IEX) 453 IXNR = IEX - ISYOFE(ISYMX) 454 WRITE(LUPRI,*) 'Total, actual nr.,loop index',NXGRST,IEX,IX 455 WRITE(LUPRI,*) 'Sym, reduced nr.',ISYMX,IXNR 456 IF (.NOT.CCS) THEN 457 KDENS = 1 458 KWRK2 = KDENS + N2BST(ISYMOP) 459 LWRK2 = LWORK - KWRK2 460C 461 IF (LWRK2 .LT. 0) THEN 462 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK2 463 CALL QUIT('Insufficient memory for one '// 464 & 'e-density in CC_EXGR') 465 ENDIF 466C 467 ILLNR = IX 468C 469Casm ifort does not like LIST. Use 'E0' instead. 470C 471 CALL CC_D1AO(IDUMMY,.FALSE.,WORK(KDENS),DUM,WORK(KWRK2), 472 & LWRK2,MODEL,'E0',ILLNR,.FALSE., 473 & FNDPTIA, FNDPTIA2, FNDPTAB, FNDPTIJ) 474casm CALL CC_D1AO(IDUMMY,.FALSE.,WORK(KDENS),DUM,WORK(KWRK2), 475casm & LWRK2,MODEL,LIST,ILLNR,.FALSE., 476casm & FNDPTIA, FNDPTIA2, FNDPTAB, FNDPTIJ) 477 IF (FROIMP .OR. FROEXP) THEN 478 CALL CC_FCD1AO(WORK(KDENS),WORK(KWRK2),LWRK2,MODELPRI) 479 ENDIF 480C 481 IF (DEBUG) THEN 482 XD = DDOT(N2BST(ISYMOP),WORK(KDENS),1,WORK(KDENS),1) 483 WRITE(LUPRI,*) 'Norm of dens-1: ',XD 484 ENDIF 485C 486 KDENS2 = KWRK2 487 KWRK3 = KDENS2 + N2BST(ISYMOP) 488 LWRK3 = LWORK - KWRK3 489 IF (LWRK2 .LT. 0) THEN 490 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK2 491 CALL QUIT('Insufficient memory for one '// 492 & 'e-density in CC_EXGR') 493 ENDIF 494 CALL DZERO(WORK(KDENS2),N2BST(ISYMOP)) 495 ILLNR = IEX 496 ILRNR = ILLNR 497 CALL CCX_D1AO(WORK(KDENS2),WORK(KWRK3),LWRK3,MODELPRI, 498 * 'LE',ILLNR,'RE',ILRNR) 499 CALL DAXPY(N2BST(ISYMOP),ONE,WORK(KDENS2),1,WORK(KDENS),1) 500C 501 IF (DEBUG) THEN 502 XD = DDOT(N2BST(ISYMOP),WORK(KDENS),1,WORK(KDENS),1) 503 WRITE(LUPRI,*) 'Norm of dens-2: ',XD 504 ENDIF 505C 506 IF (IPRINT .GT. 50) THEN 507 CALL AROUND('One electron unrelaxed density in cc_exgr') 508 CALL CC_PRFCKAO(WORK(KDENS),1) 509 ENDIF 510C 511 CALL FLSHFO(LUPRI) 512 ENDIF 513C 514C------------------------------------------------------------------------ 515C Calculate the simple one electron AO-density in CCS calculation. 516C------------------------------------------------------------------------ 517C 518 IF (CCS) THEN 519C 520 KDENS = 1 521 KWRK3 = KDENS + N2BST(ISYMOP) 522 LWRK3 = LWORK - KWRK3 523C 524 KDENS2= KWRK3 525 KWRK4 = KDENS2+ N2BST(ISYMOP) 526 LWRK4 = LWORK - KWRK3 527C 528 IF (LWRK4 .LT. 0) THEN 529 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK4 530 CALL QUIT('Insufficient memory for CCS '// 531 & 'AO-density in CC_EXGR') 532 ENDIF 533C 534C------------------------------------------ 535C Usual <HF|Emn|HF> contribution. 536C------------------------------------------ 537C 538 CALL CCS_D1AO(WORK(KDENS),WORK(KWRK3),LWRK3) 539 IF (FROIMP .OR. FROEXP) THEN 540 CALL CC_FCD1AO(WORK(KDENS),WORK(KWRK3),LWRK3,MODELPRI) 541 ENDIF 542 XN = DDOT(N2BST(ISYMOP),WORK(KDENS),1,WORK(KDENS),1) 543 WRITE(LUPRI,*) 'Norm of AO density-1' ,XN 544C 545C------------------------------------------------ 546C <LE1|[Emn,RE1]|HF> contribution. 547C For CCS but not CIS also; <L1|Emn|HF> 548C------------------------------------------------ 549C 550 ILLNR = IEX 551 ILRNR = ILLNR 552 CALL CCSX_D1AO(WORK(KDENS2),WORK(KWRK4),LWRK4, 553 * 'LE',ILLNR,'RE',ILRNR,'E0',IX) 554 XN = DDOT(N2BST(ISYMOP),WORK(KDENS2),1,WORK(KDENS2),1) 555 WRITE(LUPRI,*) 'Norm of AO density2' ,XN 556 CALL DAXPY(N2BST(ISYMOP),ONE,WORK(KDENS2),1,WORK(KDENS),1) 557 ELSE 558 KWRK3 = KWRK2 559 LWRK3 = LWRK2 560 ENDIF 561C 562 MODELPRI2 = 'Unrelaxed excited state '//MODELPRI 563 CALL AROUND(MODELPRI2//' First-order properties: ') 564C 565 ISYMX = ISYEXC(IEX) 566 IXNR = IEX - ISYOFE(ISYMX) 567 EXC = EIGVAL(IEX) 568 WRITE(LUPRI,'(/,1X,A,/,1X,A,I5,/,1X,A,I5,/,1X,A,F10.6,F12.6)') 569 * 'Excited state properties for ', 570 * 'Excited state nr.: ',IXNR, 571 * 'Excited state sym: ',ISYMX, 572 * 'Excitation energy (au.,eV): ',EXC,EXC*XTEV 573C 574C======================================= 575C Calculate molecular dipole moment. 576C======================================= 577C 578 IF (XDIPMO) THEN 579C 580 CALL AROUND(' Electric Dipole Moment ') 581C 582C------------------------------------------- 583C Calculate the nuclear contribution. 584C------------------------------------------- 585C 586 IASGER = IPRINT - 9 587 CALL DIPNUC(WORK(KWRK3),WORK(KWRK3),IASGER,.FALSE.) 588C 589 DO 100 IDIP = 1,3 590C 591 IF (IDIP .EQ. 1) LABEL = 'XDIPLEN ' 592 IF (IDIP .EQ. 2) LABEL = 'YDIPLEN ' 593 IF (IDIP .EQ. 3) LABEL = 'ZDIPLEN ' 594C 595C---------------------------------- 596C get property integrals. 597C---------------------------------- 598C 599 KONEP = KWRK3 600 KWRK4 = KONEP + N2BST(ISYMOP) 601 LWRK4 = LWORK - KWRK4 602C 603 IF (LWRK4 .LT. 0) THEN 604 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK4 605 CALL QUIT('Insufficient memory for '// 606 & 'DIPLEN-int. in CC_EXGR') 607 ENDIF 608C 609 CALL DZERO(WORK(KONEP),N2BST(ISYMOP)) 610 FF = 1.0D0 611 ISY = -1 612 CALL CC_ONEP(WORK(KONEP),WORK(KWRK4),LWRK4,FF,ISY,LABEL) 613C 614 IF (IPRINT .GT. 50) THEN 615 CALL AROUND('One electron property integrals in cc_fop') 616 CALL CC_PRFCKAO(WORK(KONEP),ISYMOP) 617 ENDIF 618C 619C---------------------------------------------- 620C Calculate the electronic contribution. 621C---------------------------------------------- 622C 623 IF (ISY .EQ. 1 ) THEN 624 DIPME(IDIP) = -DDOT(N2BST(ISYMOP),WORK(KONEP),1, 625 * WORK(KDENS),1) 626 ELSE 627 DIPME(IDIP) = 0 628 ENDIF 629 DIPMN(IDIP) = DIPMN(IDIP) + DIPME(IDIP) 630C 631C-------------------------------- 632C Store on prpc list. 633C-------------------------------- 634C 635 CALL WRIPRO(DIPMN(IDIP),MODEL,-11,LABEL,LABEL,LABEL,LABEL, 636 * EXC,DUMMY,DUMMY,ISY,ISYMX,1,IXNR) 637C 638 100 CONTINUE 639C 640C--------------------- 641C Print result. 642C--------------------- 643C 644 IF (IASGER .GT. 0) THEN 645 CALL HEADER('Electronic contribution to dipole moment',-1) 646 CALL DP0PRI(DIPME) 647 ENDIF 648 CALL HEADER('Total Molecular Dipole Moment',-1) 649 CALL DP0PRI(DIPMN) 650C 651 CALL FLSHFO(LUPRI) 652C 653 ENDIF 654C 655C=========================================== 656C Calculate molecular quadrupole moment. 657C=========================================== 658C 659 IF (XQUADR) THEN 660C 661 CALL AROUND(' Electric Quadrupole Moment ') 662C 663C------------------------------------------- 664C Calculate the nuclear contribution. 665C------------------------------------------- 666C 667 IOPT = 1 668 IASGER = -1 669 CALL CCNUCQUA(WORK(KWRK3),LWRK3,IOPT,IASGER) 670 CALL DZERO(QDREL,3*3) 671C 672 IJ = 0 673 DO 110 I = 1,3 674 DO 120 J = I,3 675 IJ = IJ + 1 676C 677 IF (IJ .EQ. 1) LABEL = 'XXTHETA ' 678 IF (IJ .EQ. 2) LABEL = 'XYTHETA ' 679 IF (IJ .EQ. 3) LABEL = 'XZTHETA ' 680 IF (IJ .EQ. 4) LABEL = 'YYTHETA ' 681 IF (IJ .EQ. 5) LABEL = 'YZTHETA ' 682 IF (IJ .EQ. 6) LABEL = 'ZZTHETA ' 683C 684C------------------------------------- 685C get property integrals. 686C------------------------------------- 687C 688 KONEP = KWRK3 689 KWRK4 = KONEP + N2BST(ISYMOP) 690 LWRK4 = LWORK - KWRK4 691C 692 IF (LWRK4 .LT. 0) THEN 693 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK4 694 CALL QUIT('Insufficient memory for THETA-int. '// 695 & 'in CC_EXGR') 696 ENDIF 697C 698 CALL DZERO(WORK(KONEP),N2BST(ISYMOP)) 699 FF = 1.0D0 700 ISY = -1 701 CALL CC_ONEP(WORK(KONEP),WORK(KWRK4),LWRK4,FF,ISY,LABEL) 702C 703 IF (IPRINT .GT. 50) THEN 704 CALL AROUND('One electron property int. in cc_fop') 705 CALL CC_PRFCKAO(WORK(KONEP),ISYMOP) 706 ENDIF 707C 708C------------------------------------------------- 709C Calculate the electronic contribution. 710C------------------------------------------------- 711C 712 LENGTH = N2BST(ISYMOP) 713C 714 IF (ISY .EQ. 1) THEN 715 CALL CCELQUA(WORK(KONEP),WORK(KDENS),LENGTH,I,J,QDREL) 716 ENDIF 717C 718 120 CONTINUE 719 110 CONTINUE 720C 721C------------------------ 722C Reorder storing. 723C------------------------ 724C 725 CALL CC_QUAREO(QDREL,SKODE) 726 CALL CC_QUAREO(QDRNUC,SKODN) 727C 728C--------------------- 729C Print result. 730C--------------------- 731C 732 IF (IPRINT .GT. 9) THEN 733 CALL HEADER('Nuclear contr. to quadrupole moment',-1) 734 WRITE(LUPRI,474) 'X','Y','Z' 735 CALL OUTPUT(SKODN,1,3,1,3,3,3,1,LUPRI) 736 CALL HEADER('Electronic contr. to quadrupole moment',-1) 737 WRITE(LUPRI,474) 'X','Y','Z' 738 CALL OUTPUT(SKODE,1,3,1,3,3,3,1,LUPRI) 739 ENDIF 740 CALL DAXPY(9,-ONE,SKODE,1,SKODN,1) 741 CALL HEADER('Total Molecular quadrupole moment',-1) 742 WRITE(LUPRI,474) 'X','Y','Z' 743 CALL OUTPUT(SKODN,1,3,1,3,3,3,1,LUPRI) 744C 745 CALL FLSHFO(LUPRI) 746C 747 ENDIF 748C 749C================================================== 750C Calculate electronic second moment of charge. 751C================================================== 752C 753 IF (XSECMO) THEN 754C 755 CALL AROUND(' Electronic second moment of charge ') 756C 757 CALL DZERO(ELSEMO,9) 758C 759 IJ = 0 760 DO 115 I = 1,3 761 DO 125 J = I,3 762 IJ = IJ + 1 763C 764 IF (IJ .EQ. 1) LABEL = 'XXSECMOM' 765 IF (IJ .EQ. 2) LABEL = 'XYSECMOM' 766 IF (IJ .EQ. 3) LABEL = 'XZSECMOM' 767 IF (IJ .EQ. 4) LABEL = 'YYSECMOM' 768 IF (IJ .EQ. 5) LABEL = 'YZSECMOM' 769 IF (IJ .EQ. 6) LABEL = 'ZZSECMOM' 770C 771C------------------------------------- 772C get property integrals. 773C------------------------------------- 774C 775 KONEP = KWRK3 776 KWRK4 = KONEP + N2BST(ISYMOP) 777 LWRK4 = LWORK - KWRK4 778C 779 IF (LWRK4 .LT. 0) THEN 780 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK4 781 CALL QUIT('Insufficient memory for '// 782 & 'SECMOM-int. in CC_EXGR') 783 ENDIF 784C 785 CALL DZERO(WORK(KONEP),N2BST(ISYMOP)) 786 FF = 1.0D0 787 ISY = -1 788 CALL CC_ONEP(WORK(KONEP),WORK(KWRK4),LWRK4,FF,ISY,LABEL) 789C 790 IF (IPRINT .GT. 50) THEN 791 CALL AROUND('One electron property int. in cc_fop') 792 CALL CC_PRFCKAO(WORK(KONEP),ISYMOP) 793 ENDIF 794C 795C------------------------------------------------- 796C Calculate the electronic contribution. 797C------------------------------------------------- 798C 799 LENGTH = N2BST(ISYMOP) 800C 801 IF (ISY.EQ.1) THEN 802 CALL CCELQUA(WORK(KONEP),WORK(KDENS),LENGTH,I,J,ELSEMO) 803 ENDIF 804C 805C-------------------------------- 806C Store on prpc list. 807C-------------------------------- 808C 809 CALL WRIPRO(ELSEMO(I,J),MODEL,-11,LABEL,LABEL,LABEL,LABEL, 810 * EXC,DUMMY,DUMMY,ISY,ISYMX,1,IXNR) 811C 812 125 CONTINUE 813 115 CONTINUE 814C 815C------------------------ 816C Reorder storing. 817C------------------------ 818C 819 CALL CC_QUAREO(ELSEMO,SKODE) 820C 821C--------------------- 822C Print result. 823C--------------------- 824C 825 WRITE(LUPRI,474) 'X','Y','Z' 826 CALL OUTPUT(SKODE,1,3,1,3,3,3,1,LUPRI) 827 CALL CC_TNSRAN(SKODE,WORK(KWRK3),LWRK3) 828C 829 CALL FLSHFO(LUPRI) 830C 831 ENDIF 832C 833 474 FORMAT(20X,A1,14X,A1,14X,A1) 834C 835C======================================= 836C Calculate electric field gradient. 837C======================================= 838C 839 IF (XNQCC) THEN 840C 841 CALL AROUND(' Electric Field Gradients ') 842C 843C------------------------------------------- 844C Calculate the nuclear contribution. 845C------------------------------------------- 846C 847 IOPT = 2 848 IASGER = IPRINT - 5 849 CALL CCNUCQUA(WORK(KWRK3),LWRK3,IOPT,IASGER) 850C 851C---------------------------------------------- 852C Calculate the electronic contribution. 853C---------------------------------------------- 854C 855 LENGTH = N2BST(ISYMOP) 856 CALL CCELEFG(WORK(KDENS),LENGTH,WORK(KWRK3),LWRK3,IASGER) 857C 858C--------------------- 859C Print result. 860C--------------------- 861C 862 KDIAG = KWRK3 863 KAXIS = KDIAG + 3*MXCENT 864 KWRK4 = KAXIS + 9*MXCENT 865 LWRK4 = LWORK - KWRK4 866C 867 IF (LWRK4 .LT. 0) THEN 868 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK4 869 CALL QUIT('Insufficient memory for EFG-results in CC_EXGR') 870 ENDIF 871C 872 IASGER = 2 873 ICCPRI = 2 874 CALL NQCRES(IASGER,WORK(KDIAG),WORK(KAXIS),ICCPRI) 875C 876 CALL FLSHFO(LUPRI) 877C 878 ENDIF 879C 880C========================================================= 881C Relativistic corrections to the ground-state energy. 882C========================================================= 883C 884 IF (XRELCO) THEN 885C 886 CALL AROUND(' Relativistic corrections to the ground-state' 887 * //' energy ') 888C 889 DO 130 IRC = 1,2 890C 891 IF (IRC .EQ. 1) LABEL = 'DARWIN ' 892 IF (IRC .EQ. 2) LABEL = 'MASSVELO' 893C 894C----------------------------- 895C get the integrals. 896C----------------------------- 897C 898 KONEP = KWRK3 899 KWRK4 = KONEP + N2BST(ISYMOP) 900 LWRK4 = LWORK - KWRK4 901C 902 IF (LWRK4 .LT. 0) THEN 903 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK4 904 CALL QUIT('Insufficient memory for '// 905 & 'Darwin-int. in CC_EXGR') 906 ENDIF 907C 908 CALL DZERO(WORK(KONEP),N2BST(ISYMOP)) 909 FF = 1.0D0 910 ISY = 1 911 CALL CC_ONEP(WORK(KONEP),WORK(KWRK4),LWRK4,FF,ISY,LABEL) 912C 913 IF (IPRINT .GT. 50) THEN 914 CALL AROUND('Relativistic integrals in cc_fop') 915 CALL CC_PRFCKAO(WORK(KONEP),ISYMOP) 916 ENDIF 917C 918C------------------------------------- 919C Calculate the corrections. 920C------------------------------------- 921C 922 IF (IRC .EQ. 1) THEN 923 DARW = DDOT(N2BST(ISYMOP),WORK(KONEP),1,WORK(KDENS),1) 924 ELSE IF (IRC .EQ. 2) THEN 925 VELO = DDOT(N2BST(ISYMOP),WORK(KONEP),1,WORK(KDENS),1) 926 ENDIF 927C 928 130 CONTINUE 929C 930C---------------------- 931C Write out result. 932C---------------------- 933C 934 WRITE(LUPRI,*) ' ' 935 WRITE(LUPRI,131) 'Darwin term:', DARW, 'Mass-Velocity term:', VELO 936 WRITE(LUPRI,132) '----------- ', '------------------ ' 937 WRITE(LUPRI,*) ' ' 938 WRITE(LUPRI,*) ' ' 939 WRITE(LUPRI,133) 'Total relativistic correction:', DARW + VELO 940 WRITE(LUPRI,134) '----------------------------- ' 941 WRITE(LUPRI,*) ' ' 942 WRITE(LUPRI,*) ' ' 943C 944 131 FORMAT(9X,A12,1X,F9.6,7X,A19,1X,F9.6) 945 132 FORMAT(9X,A12,17X,A19) 946 133 FORMAT(19X,A30,1X,F9.6) 947 134 FORMAT(19X,A30) 948C 949 ENDIF 950C 951C-------------------------------------------------------------- 952C Section for general operator APROP represented by LABEL. 953C Note that only the electronic contribution is calculated. 954C-------------------------------------------------------------- 955C 956 DO 140 IOP = 1, NAXGRO 957C 958 LABEL = LBLOPR(IAXGRO(IOP)) 959C 960 IF (IOP .EQ. 1) CALL AROUND( 961 * ' Electronic contribution to operator ') 962C 963C-------------------------- 964C get the integrals. 965C-------------------------- 966C 967 KONEP = KWRK3 968 KWRK4 = KONEP + N2BST(ISYMOP) 969 LWRK4 = LWORK - KWRK4 970C 971 IF (LWRK4 .LT. 0) THEN 972 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KWRK4 973 CALL QUIT('Insufficient memory for property '// 974 & 'integrals in CC_EXGR') 975 ENDIF 976C 977 CALL DZERO(WORK(KONEP),N2BST(ISYMOP)) 978 FF = 1.0D0 979 ISY = -1 980 CALL CC_ONEP(WORK(KONEP),WORK(KWRK4),LWRK4,FF,ISY,LABEL) 981C 982 IF (IPRINT .GT. 50) THEN 983 CALL AROUND('Property integrals in cc_fop') 984 CALL CC_PRFCKAO(WORK(KONEP),ISYMOP) 985 ENDIF 986C 987C-------------------------------------------------------------------- 988C Calculate the electronic contribution to the given property. 989C-------------------------------------------------------------------- 990C 991 PROP = DDOT(N2BST(ISYMOP),WORK(KONEP),1,WORK(KDENS),1) 992 CALL WRIPRO(PROP,MODEL,-11,LABEL,LABEL,LABEL,LABEL, 993 * EXC,DUMMY,DUMMY,ISY,ISYMX,1,IXNR) 994 995C 996C---------------------- 997C Write out result. 998C---------------------- 999C 1000 WRITE(LUPRI,*) ' ' 1001 IF (ISY.EQ.1) WRITE(LUPRI,141) LABEL//':', PROP 1002 IF (ISY.NE.1) WRITE(LUPRI,142) LABEL//':', 'zero by symmetry' 1003 WRITE(LUPRI,*) ' ' 1004 WRITE(LUPRI,*) ' ' 1005C 1006 141 FORMAT(20X,A9,1X,F10.6) 1007 142 FORMAT(20X,A9,1X,A) 1008C 1009 140 CONTINUE 1010C 1011C------------------------------------------------------------------------ 1012C End of loop for write out of excited state one-electron properties. 1013C------------------------------------------------------------------------ 1014C 1015 3000 CONTINUE 1016C 1017C------------------------------------ 1018C Restore NSIDE to input/default. 1019C------------------------------------ 1020C 1021 IF (.NOT. CCS) NSIDE = NSIDIN 1022C 1023 CALL QEXIT('CC_EXGR') 1024C 1025 RETURN 1026 END 1027C /* Deck ccx_d1ao */ 1028 SUBROUTINE CCX_D1AO(AODEN,WORK,LWORK,MODEL, 1029 * LLIST,ILLNR,RLIST,ILRNR) 1030C 1031C Ove Christiansen April 1997 inspired by CC_D1AO 1032C KS Aug 2010: minor modifications to share routines 1033C with CCMM_D2AO 1034C 1035C Purpose: To calculate contributions to the excited state 1036C one electron density matrix and return it backtransformed 1037C to AO-basis in AODEN. 1038C 1039C Current models: CCS, CC2 and CCSD 1040C 1041C 1042#include "implicit.h" 1043#include "priunit.h" 1044#include "dummy.h" 1045#include "maxash.h" 1046#include "maxorb.h" 1047#include "mxcent.h" 1048#include "aovec.h" 1049#include "iratdef.h" 1050 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 1051 DIMENSION AODEN(*), WORK(LWORK) 1052#include "ccorb.h" 1053#include "infind.h" 1054#include "blocks.h" 1055#include "ccsdinp.h" 1056#include "ccsdsym.h" 1057#include "ccexci.h" 1058#include "ccsdio.h" 1059#include "distcl.h" 1060#include "cbieri.h" 1061#include "eribuf.h" 1062#include "cclr.h" 1063C 1064 CHARACTER MODEL*10,MODDUM*10 1065 CHARACTER LLIST*(*),RLIST*(*) 1066C 1067 CALL QENTER('CCX_D1AO') 1068C 1069C--------------------------------------------- 1070C Find symmetry of left and right vectors. 1071C--------------------------------------------- 1072C 1073 ISYMR = ISYEXC(ILRNR) 1074 ISYML = ISYEXC(ILLNR) 1075 IF (ISYMR .NE. ISYML) 1076 & CALL QUIT('CCX_D1AO: Density not total sym.') 1077C 1078C----------------------------------- 1079C Initial work space allocation. 1080C----------------------------------- 1081C 1082 KONEAI = 1 1083 KONEAB = KONEAI + NT1AMX 1084 KONEIJ = KONEAB + NMATAB(1) 1085 KONEIA = KONEIJ + NMATIJ(1) 1086 KL1AM = KONEIA + NT1AMX 1087 KL2AM = KL1AM + NT1AM(ISYML) 1088 KT1AM = KL2AM + NT2SQ(ISYML) 1089 KR1AM = KT1AM + NT1AM(ISYMOP) 1090 KEND0 = KR1AM + NT1AM(ISYMR) 1091C 1092 KR2AM = KEND0 1093 KR2AMT = KR2AM + NT2AM(ISYMR) 1094 KEND1 = KR2AMT + NT2AM(ISYMR) 1095 LWRK1 = LWORK - KEND1 1096C 1097 IF (LWRK1 .LT. 0) THEN 1098 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1099 CALL QUIT('Insufficient memory for initial '// 1100 & 'allocation in CCX_D1AO') 1101 ENDIF 1102C 1103C-------------------------------------------- 1104C Initialize one electron density arrays. 1105C-------------------------------------------- 1106C 1107 CALL DZERO(WORK(KONEAB),NMATAB(1)) 1108 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 1109 CALL DZERO(WORK(KONEAI),NT1AMX) 1110 CALL DZERO(WORK(KONEIA),NT1AMX) 1111C 1112C---------------------- 1113C Read left vector. 1114C---------------------- 1115C 1116 IOPT = 3 1117 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODDUM, 1118 * WORK(KL1AM),WORK(KR2AM)) 1119C 1120C-------------------------------- 1121C Square up zeta2 amplitudes. 1122C-------------------------------- 1123C 1124 CALL CC_T2SQ(WORK(KR2AM),WORK(KL2AM),ISYML) 1125C 1126C---------------------- 1127C Read rigth vector. 1128C---------------------- 1129C 1130 IOPT = 3 1131 CALL CC_RDRSP(RLIST,ILRNR,ISYMR,IOPT,MODDUM, 1132 * WORK(KR1AM),WORK(KR2AM)) 1133 IF (ISYMR.EQ.1) CALL CCLR_DIASCL(WORK(KR2AM),TWO,ISYMR) 1134C 1135C--------------------------------------------------- 1136C Read zero'th order cluster singles amplitudes. 1137C--------------------------------------------------- 1138C 1139 IOPT = 1 1140 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 1141C 1142C--------------------------------------- 1143C Set up 2C-E of cluster amplitudes. 1144C--------------------------------------- 1145C 1146 IF (.NOT. MP2) THEN 1147 CALL DCOPY(NT2AM(ISYMR),WORK(KR2AM),1,WORK(KR2AMT),1) 1148 IOPTTCME = 1 1149 CALL CCSD_TCMEPK(WORK(KR2AMT),1.0D0,ISYMR,IOPTTCME) 1150 ENDIF 1151C 1152C--------------------- 1153C Test amplitudes. 1154C--------------------- 1155C 1156 IF ( DEBUG ) THEN 1157 XLV1 = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 1158 XLV2 = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 1159 WRITE(LUPRI,1) 'Norm of Response vector: L1AM ',XLV1 1160 WRITE(LUPRI,1) 'Norm of Response vector: L2AM ',XLV2 1161 XLV1 = DDOT(NT1AM(ISYML),WORK(KR1AM),1,WORK(KR1AM),1) 1162 XLV2 = DDOT(NT2AM(ISYML),WORK(KR2AM),1,WORK(KR2AM),1) 1163 WRITE(LUPRI,1) 'Norm of Response vector: R1AM ',XLV1 1164 WRITE(LUPRI,1) 'Norm of Response vector: R2AM ',XLV2 1165 XLV2 = DDOT(NT2AM(ISYML),WORK(KR2AMT),1,WORK(KR2AMT),1) 1166 WRITE(LUPRI,1) 'Norm of Response vector: R2AM-Tr.',XLV2 1167 XLV1 = DDOT(NT1AM(ISYMOP),WORK(KT1AM),1,WORK(KT1AM),1) 1168 WRITE(LUPRI,1) 'Norm of referencevector: T1AM ',XLV1 1169 ENDIF 1170C 1171C------------------------------- 1172C Work space allocation one. 1173C Note that D(ai)=ZETA(ai) & 1174C D(ia) is stored transposed 1175C------------------------------- 1176C 1177 KXMAT = KEND1 1178 KYMAT = KXMAT + NMATIJ(1) 1179 KEND2 = KYMAT + NMATAB(1) 1180 LWRK2 = LWORK - KEND1 1181C 1182 IF (LWRK2 .LT. 0) THEN 1183 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1184 CALL QUIT('Insufficient memory for 2. alloc. in CCX_D1AO') 1185 ENDIF 1186C 1187C----------------------------------------------------------- 1188C Calculate X-intermediate of L2AM- and R2AM-amplitudes. 1189C----------------------------------------------------------- 1190C 1191 CALL CC_XI(WORK(KXMAT),WORK(KL2AM),ISYML,WORK(KR2AM),ISYMR, 1192 * WORK(KEND2),LWRK2) 1193C 1194C----------------------------------------------------------- 1195C Calculate Y-intermediate of L2AM- and R2AM-amplitudes. 1196C----------------------------------------------------------- 1197C 1198 CALL CC_YI(WORK(KYMAT),WORK(KL2AM),ISYML,WORK(KR2AM),ISYMR, 1199 * WORK(KEND2),LWRK2) 1200C 1201C-------------------------------------------------------------- 1202C DXai is zero. 1203C-------------------------------------------------------------- 1204C Construct <L2|[Emn,R2]|HF> contribution to DXaa and DXii. 1205C-------------------------------------------------------------- 1206C 1207 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 1208 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND2),LWRK2,1) 1209 CALL DAXPY(NMATIJ(1),-ONE,WORK(KXMAT),1,WORK(KONEIJ),1) 1210C 1211C-------------------------------------------------------------- 1212C Construct <L1|[Emn,R1]|HF> contribution to DXaa and DXii. 1213C-------------------------------------------------------------- 1214C 1215c test 1216 CALL CC_DXIJ(WORK(KONEIJ),WORK(KL1AM),ISYML,WORK(KR1AM),ISYMR) 1217 CALL CC_DXAB(WORK(KONEAB),WORK(KL1AM),ISYML,WORK(KR1AM),ISYMR) 1218C 1219C-------------------------------------------------------------- 1220C Construct <L1|[Eia,R2]|HF> contribution to DXaa and DXii. 1221C-------------------------------------------------------------- 1222C 1223 CALL CC_DXIA12(WORK(KONEIA),WORK(KL1AM),ISYML,WORK(KR2AMT),ISYMR) 1224C 1225C---------------------------------------------------------- 1226C Construct <L2|[[Eia,R1],T2]|HF> contribution to DXia. 1227C---------------------------------------------------------- 1228C 1229 KT2AM = KEND0 1230 KEND3 = KT2AM + NT2AM(ISYMOP) 1231 LWRK3 = LWORK - KEND3 1232 1233C KS: Read amplitudes outside DXIA21 routine 1234 IOPT = 2 1235 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KT2AM)) 1236C 1237 IF (LWRK3 .LT. 0) THEN 1238 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND3 1239 CALL QUIT('Insufficient memory for 3. alloc. in CCX_D1AO') 1240 ENDIF 1241 CALL CC_DXIA21(WORK(KONEIA),WORK(KL2AM),ISYML, 1242 * WORK(KR1AM),ISYMR,WORK(KT2AM),ISYMOP, 1243 * WORK(KEND3),LWRK3) 1244C 1245C-------------------------- 1246C Write out test norms. 1247C-------------------------- 1248C 1249 IF (DEBUG) THEN 1250 XIJ = DDOT(NMATIJ(1),WORK(KONEIJ),1,WORK(KONEIJ),1) 1251 XAB = DDOT(NMATAB(1),WORK(KONEAB),1,WORK(KONEAB),1) 1252 XAI = DDOT(NT1AMX,WORK(KONEAI),1,WORK(KONEAI),1) 1253 XIA = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 1254 WRITE(LUPRI,*) 'Norms: DXIJ',XIJ 1255 WRITE(LUPRI,*) 'Norms: DXAB',XAB 1256 WRITE(LUPRI,*) 'Norms: DXAI',XAI 1257 WRITE(LUPRI,*) 'Norms: DXIA',XIA 1258 ENDIF 1259C 1260C---------------------------------- 1261C Calculate the lamda matrices. 1262C---------------------------------- 1263C 1264 KLAMDP = KEND0 1265 KLAMDH = KLAMDP + NLAMDT 1266 KEND4 = KLAMDH + NLAMDT 1267 LWRK4 = LWORK - KEND4 1268 IF (LWRK4 .LT. 0) THEN 1269 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1270 CALL QUIT('Insufficient memory for 4. alloc. in CCX_D1AO') 1271 ENDIF 1272 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND4), 1273 * LWRK4) 1274C 1275C-------------------------------------------------------- 1276C Add the one electron density in the AO-basis. 1277C-------------------------------------------------------- 1278C 1279 ISDEN = 1 1280 CALL CC_DENAO(AODEN,ISDEN,WORK(KONEAI),WORK(KONEAB), 1281 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 1282 * WORK(KLAMDH),1,WORK(KEND4),LWRK4) 1283C 1284 CALL QEXIT('CCX_D1AO') 1285C 1286 1 FORMAT(1x,A35,1X,E20.10) 1287 RETURN 1288 END 1289C /* Deck cc_dxij */ 1290 SUBROUTINE CC_DXIJ(DIJ,XL1AM,ISYML,XR1AM,ISYMR) 1291C 1292C Written by Ove Christiansen April 1997 1293C 1294C Purpose: To add contributions from the L1 and R1 -amplitudes to 1295C the excited state CC one electron density in MO-basis. 1296C <L1|[Eij,R1]|HF> 1297C 1298#include "implicit.h" 1299 PARAMETER(ZERO = 0.0D0, XMONE = -1.0D0, ONE = 1.0D00) 1300 DIMENSION DIJ(*), XL1AM(*), XR1AM(*) 1301#include "priunit.h" 1302#include "ccorb.h" 1303#include "ccsdsym.h" 1304C 1305 CALL QENTER('CC_DXIJ') 1306C 1307C--------------------------------------------------------- 1308C Add contribution. 1309C Assume ISYMR = ISYML and density is total symmetric. 1310C--------------------------------------------------------- 1311C 1312 ISYRES = MULD2H(ISYMR,ISYML) 1313 IF (ISYRES .NE. ISYMOP) CALL QUIT('CC_DXIJ require tot.sym Dens.') 1314C 1315 DO 100 ISYMI = 1,NSYM 1316C 1317 ISYMJ = MULD2H(ISYRES,ISYMI) 1318 ISYMC = MULD2H(ISYMR,ISYMI) 1319C 1320 KOFF1 = IT1AM(ISYMC,ISYMI) + 1 1321 KOFF2 = IT1AM(ISYMC,ISYMJ) + 1 1322 KOFF3 = IMATIJ(ISYMI,ISYMJ) + 1 1323C 1324 NTOTC = MAX(NVIR(ISYMC),1) 1325 NTOTI = MAX(NRHF(ISYMI),1) 1326C 1327 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC), 1328 * XMONE,XR1AM(KOFF1),NTOTC,XL1AM(KOFF2),NTOTC,ONE, 1329 * DIJ(KOFF3),NTOTI) 1330C 1331 100 CONTINUE 1332C 1333 CALL QEXIT('CC_DXIJ') 1334C 1335 END 1336C /* Deck cc_dxab */ 1337 SUBROUTINE CC_DXAB(DAB,XL1AM,ISYML,XR1AM,ISYMR) 1338C 1339C Written by Ove Christiansen April 1997 1340C 1341C Purpose: To add contributions from the L1 and R1 -amplitudes to 1342C the excited state CC one electron density in MO-basis. 1343C <L1|[Eab,R1]|HF> 1344C 1345#include "implicit.h" 1346 PARAMETER(ZERO = 0.0D0, XMONE = -1.0D0, ONE = 1.0D00) 1347 DIMENSION DAB(*), XL1AM(*), XR1AM(*) 1348#include "priunit.h" 1349#include "ccorb.h" 1350#include "ccsdsym.h" 1351C 1352 CALL QENTER('CC_DXAB') 1353C 1354C--------------------------------------------------------- 1355C Add contribution. 1356C Assume ISYMR = ISYML and density is total symmetric. 1357C--------------------------------------------------------- 1358C 1359 ISYRES = MULD2H(ISYMR,ISYML) 1360 IF (ISYRES .NE. ISYMOP) CALL QUIT('CC_DXAB require tot.sym Dens.') 1361C 1362 DO 100 ISYMA = 1,NSYM 1363C 1364 ISYMB = MULD2H(ISYMA,ISYRES) 1365 ISYMK = MULD2H(ISYMA,ISYMR) 1366C 1367 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1368 KOFF2 = IT1AM(ISYMB,ISYMK) + 1 1369 KOFF3 = IMATAB(ISYMA,ISYMB) + 1 1370C 1371 NTOTA = MAX(NVIR(ISYMA),1) 1372 NTOTB = MAX(NVIR(ISYMB),1) 1373C 1374 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK), 1375 * ONE,XL1AM(KOFF1),NTOTA,XR1AM(KOFF2),NTOTB,ONE, 1376 * DAB(KOFF3),NTOTA) 1377C 1378 100 CONTINUE 1379C 1380 CALL QEXIT('CC_DXAB') 1381C 1382 END 1383C /* Deck cc_dxai12 */ 1384 SUBROUTINE CC_DXIA12(DIA,XL1AM,ISYML,XR2AM,ISYMR) 1385C 1386C Written by Ove Christiansen April 1997 1387C 1388C Purpose: To add contributions from the L1 and R2 -amplitudes to 1389C the excited state CC one electron density in MO-basis. 1390C <L1|[Eia,R2]|HF> 1391C NB - 2c-E XR2AM assumed 1392C and the Dia block is stored transposed, i.e. like a t1-amplitude! 1393C 1394C 1395#include "implicit.h" 1396#include "priunit.h" 1397#include "ccorb.h" 1398#include "ccsdsym.h" 1399 PARAMETER(ZERO = 0.0D0, XMONE = -1.0D0, ONE = 1.0D00) 1400 DIMENSION DIA(*), XL1AM(*), XR2AM(*) 1401C 1402 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1403C 1404 CALL QENTER('CC_DXIA12') 1405C 1406C--------------------------------------------------------- 1407C Assume ISYMR = ISYML and density is total symmetric. 1408C--------------------------------------------------------- 1409C 1410 ISYRES = MULD2H(ISYMR,ISYML) 1411 IF (ISYRES .NE. ISYMOP) CALL QUIT('CC_DXIA12 require '// 1412 & 'tot.sym Dens.') 1413C 1414 ISYMAI = ISYRES 1415 ISYMDK = MULD2H(ISYMR,ISYMAI) 1416C 1417C-------------------------- 1418C Set up density block. 1419C-------------------------- 1420C 1421 DO 100 NAI = 1,NT1AM(ISYMAI) 1422 DO 110 NDK = 1,NT1AM(ISYMDK) 1423C 1424 IF (ISYMAI .EQ. ISYMDK) THEN 1425 NDKAI = IT2AM(ISYMDK,ISYMAI) + INDEX(NDK,NAI) 1426 ELSE IF (ISYMDK .LT. ISYMAI) THEN 1427 NDKAI = IT2AM(ISYMDK,ISYMAI) 1428 * + NT1AM(ISYMDK)*(NAI-1) + NDK 1429 ELSE IF (ISYMDK .GT. ISYMAI) THEN 1430 NDKAI = IT2AM(ISYMDK,ISYMAI) 1431 * + NT1AM(ISYMAI)*(NDK-1) + NAI 1432 ENDIF 1433 DIA(NAI) = DIA(NAI) + XR2AM(NDKAI)*XL1AM(NDK) 1434C 1435 110 CONTINUE 1436 100 CONTINUE 1437C 1438 CALL QEXIT('CC_DXIA12') 1439C 1440 RETURN 1441 END 1442C /* Deck cc_dxia21 */ 1443 SUBROUTINE CC_DXIA21(DIA,XL2AM,ISYML,XR1AM,ISYMR, 1444 * T2AM,ISYMT2,WORK,LWORK) 1445C 1446C Written by Ove Christiansen April 1997 1447C 1448C Purpose: Construct <L2|[[Eia,R1],T2]|HF> contribution to DXia, 1449C the excited state CC one electron density in MO-basis. 1450C The Dia block is stored transposed, i.e. like a t1-amplitude. 1451C KS, Aug 10: 1452C Modified to calculate general <L2|[[Eia,R1],R2]|HF> 1453C 1454C 1455#include "implicit.h" 1456#include "priunit.h" 1457#include "dummy.h" 1458#include "ccorb.h" 1459#include "ccsdsym.h" 1460#include "ccsdinp.h" 1461 PARAMETER(ZERO = 0.0D0, XMONE = -1.0D0, ONE = 1.0D00) 1462 DIMENSION DIA(*), XL2AM(*), XR1AM(*), T2AM(*),WORK(LWORK) 1463 CHARACTER*10 MODEL 1464 LOGICAL LOCDEB 1465 PARAMETER (LOCDEB=.FALSE.) 1466C 1467 CALL QENTER('CC_DXIA21') 1468C 1469C--------------------------------------------------------- 1470C Assume ISYMR = ISYML and density is total symmetric. 1471C--------------------------------------------------------- 1472C 1473 ISYRES = MULD2H(ISYMR,ISYML) 1474 IF (ISYRES .NE. ISYMOP) 1475 & CALL QUIT('CC_DXIA21 require tot.sym Dens.') 1476 IF (ISYRES .NE. ISYMT2) 1477 & CALL QUIT('symmetry mismatch in CC_DXIA21.') 1478C 1479 IF (LOCDEB) THEN 1480 TAL1 = DDOT(NT2AMX,XL2AM,1,XL2AM,1) 1481 WRITE(LUPRI,*) 'CC_EXGR: Norm of L0_2: ', TAL1 1482 TAL1 = DDOT(NT1AMX,XR1AM,1,XR1AM,1) 1483 WRITE(LUPRI,*) 'CC_EXGR: Norm of B_1: ', TAL1 1484 TAL1 = DDOT(NT2AMX,T2AM,1,T2AM,1) 1485 WRITE(LUPRI,*) 'CC_EXGR: Norm of C_2: ', TAL1 1486 END IF 1487C KS schedule delete next 6 lines 1488C since they are now read outside routine and carried 1489C in memory at this time 1490C--------------------------------------------------- 1491C Read zero'th order cluster doubles amplitudes. 1492C--------------------------------------------------- 1493C 1494C IOPT = 2 1495C CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,T2AM) 1496C 1497C---------------------------------------------------------- 1498C Construct <L2|[[Eia,R1],T2]|HF> contribution to DXia. 1499C---------------------------------------------------------- 1500C 1501 KONEAB = 1 1502 KONEIJ = KONEAB + NMATAB(ISYML) 1503 KXMAT = KONEIJ + NMATIJ(ISYML) 1504 KYMAT = KXMAT + NMATIJ(ISYML) 1505 KEND1 = KYMAT + NMATAB(ISYML) 1506 LWRK1 = LWORK - KEND1 1507C 1508 IF (LWRK1 .LT. 0) THEN 1509 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1510 CALL QUIT('Insufficient memory for 1. alloc. in CC_DXIA21.') 1511 ENDIF 1512C 1513C----------------------------------------------------------- 1514C Calculate X-intermediate of L2AM- and T2AM-amplitudes. 1515C----------------------------------------------------------- 1516C 1517 CALL CC_XI(WORK(KXMAT),XL2AM,ISYML,T2AM,ISYMT2, 1518 * WORK(KEND1),LWRK1) 1519C 1520C----------------------------------------------------------- 1521C Calculate Y-intermediate of L2AM- and T2AM-amplitudes. 1522C----------------------------------------------------------- 1523C 1524 CALL CC_YI(WORK(KYMAT),XL2AM,ISYML,T2AM,ISYMT2, 1525 * WORK(KEND1),LWRK1) 1526C TAL1 = DDOT(NMATIJ(ISYML),WORK(KXMAT),1,WORK(KXMAT),1) 1527C TAL2 = DDOT(NMATAB(ISYML),WORK(KYMAT),1,WORK(KYMAT),1) 1528C WRITE(LUPRI,*) 'CC_EXGR: X intermediate: ',TAL1 1529C WRITE(LUPRI,*) 'CC_EXGR: Y intermediate: ',TAL2 1530C 1531 CALL DZERO(WORK(KONEIJ),NMATIJ(ISYML)) 1532 CALL DCOPY(NMATAB(ISYML),WORK(KYMAT),1,WORK(KONEAB),1) 1533 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,ISYML) 1534 CALL DAXPY(NMATIJ(ISYML),XMONE,WORK(KXMAT),1,WORK(KONEIJ),1) 1535 IF (LOCDEB) THEN 1536 XIJ = DDOT(NMATIJ(ISYML),WORK(KONEIJ),1,WORK(KONEIJ),1) 1537 XAB = DDOT(NMATAB(ISYML),WORK(KONEAB),1,WORK(KONEAB),1) 1538 WRITE(LUPRI,*) 'Norms: DXIJ',XIJ 1539 WRITE(LUPRI,*) 'Norms: DXAB',XAB 1540 ENDIF 1541C 1542 CALL CC_IA21(DIA,WORK(KONEAB),WORK(KONEIJ),ISYML,XR1AM,ISYMR, 1543 * WORK(KEND1),LWRK1) 1544C 1545 CALL QEXIT('CC_DXIA21') 1546C 1547 END 1548 SUBROUTINE CC_IA21(DIA,DAB,DIJ,ISYMD,XR1AM,ISYMR,WORK,LWORK) 1549C 1550C Written by Ove Christiansen April 1997 1551C 1552C Version: 1.0 1553C 1554C Purpose: Contract DAB and DIJ with R1AM into DIA 1555C The Dia block is stored transposed, i.e. like a t1-amplitude! 1556C 1557#include "implicit.h" 1558 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, XMONE = -1.0D0) 1559 DIMENSION DIA(*),DAB(*),DIJ(*),XR1AM(*),WORK(LWORK) 1560#include "priunit.h" 1561#include "ccorb.h" 1562#include "ccsdsym.h" 1563C 1564 CALL QENTER('CC_IA21') 1565C 1566C----------------------------------------------------- 1567C Assume ISYMR = ISYMD and Dia is total symmetric. 1568C----------------------------------------------------- 1569C 1570 ISYRES = MULD2H(ISYMR,ISYMD) 1571 IF (ISYRES .NE. ISYMOP) CALL QUIT('CC_IA21 require tot.sym Dens.') 1572 ISYMAI = ISYRES 1573C 1574 DO 100 ISYMA = 1,NSYM 1575C 1576 ISYMI = MULD2H(ISYMA,ISYMAI) 1577 ISYMK = MULD2H(ISYMA,ISYMR) 1578C 1579 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1580 KOFF2 = IMATIJ(ISYMI,ISYMK) + 1 1581 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1582C 1583 NTOTA = MAX(NVIR(ISYMA),1) 1584 NTOTI = MAX(NRHF(ISYMI),1) 1585C 1586 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 1587 * ONE,XR1AM(KOFF1),NTOTA,DIJ(KOFF2),NTOTI,ONE, 1588 * DIA(KOFF3),NTOTA) 1589C 1590 ISYMC = MULD2H(ISYMI,ISYMR) 1591C 1592 KOFF1 = IMATAB(ISYMC,ISYMA) + 1 1593 KOFF2 = IT1AM(ISYMC,ISYMI) + 1 1594 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1595C 1596 NTOTA = MAX(NVIR(ISYMA),1) 1597 NTOTC = MAX(NVIR(ISYMC),1) 1598 NTOTI = MAX(NRHF(ISYMA),1) 1599C 1600 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMC), 1601 * XMONE,DAB(KOFF1),NTOTC,XR1AM(KOFF2),NTOTC,ONE, 1602 * DIA(KOFF3),NTOTA) 1603C 1604 100 CONTINUE 1605C 1606 CALL QEXIT('CC_IA21') 1607C 1608 END 1609 SUBROUTINE CCSX_D1AO(AODEN,WORK,LWORK, 1610 * LLIST,ILLNR,RLIST,ILRNR,L0LIST,ILNRL0) 1611C 1612C Ove Christiansen April 1997 inspired by CC_D1AO 1613C 1614C Purpose: To calculate contributions to the excited state 1615C one electron density matrix and return it backtransformed 1616C to AO-basis in AODEN. 1617C <LE1|[Emn,RE1]|HF> contribution. 1618C For CCS but not CIS also; 1619C <L1|Emn|HF> 1620C 1621C Current models: CCS, CC2 and CCSD 1622C 1623#include "implicit.h" 1624#include "priunit.h" 1625#include "maxash.h" 1626#include "maxorb.h" 1627#include "mxcent.h" 1628#include "aovec.h" 1629#include "iratdef.h" 1630 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 1631 DIMENSION AODEN(*), WORK(LWORK) 1632#include "ccorb.h" 1633#include "infind.h" 1634#include "blocks.h" 1635#include "ccsdinp.h" 1636#include "ccsdsym.h" 1637#include "ccexci.h" 1638#include "ccsdio.h" 1639#include "distcl.h" 1640#include "cbieri.h" 1641#include "eribuf.h" 1642#include "cclr.h" 1643C 1644 CHARACTER MODEL*5,MODDUM*10 1645 CHARACTER LLIST*(*),RLIST*(*), L0LIST*(*) 1646C 1647 CALL QENTER('CCSX_D1AO') 1648C 1649C--------------------------------------------- 1650C Find symmetry of left and right vectors. 1651C--------------------------------------------- 1652C 1653 ISYMR = ISYEXC(ILRNR) 1654 ISYML = ISYEXC(ILLNR) 1655 IF (ISYMR .NE. ISYML) 1656 & CALL QUIT('CCSX_D1AO: Density not total sym.') 1657C 1658C 1659C----------------------------------- 1660C Initial work space allocation. 1661C----------------------------------- 1662C 1663 KONEAI = 1 1664 KONEAB = KONEAI + NT1AMX 1665 KONEIJ = KONEAB + NMATAB(1) 1666 KONEIA = KONEIJ + NMATIJ(1) 1667 KL1AM = KONEIA + NT1AMX 1668 KR1AM = KL1AM + NT1AM(ISYML) 1669 KT1AM = KR1AM + NT1AM(ISYMR) 1670 KLAMDP = KT1AM + NT1AM(ISYMOP) 1671 KLAMDH = KLAMDP + NLAMDT 1672 KEND1 = KLAMDH + NLAMDT 1673 LWRK1 = LWORK - KEND1 1674C 1675 IF (LWRK1 .LT. 0) THEN 1676 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1677 CALL QUIT('Insufficient memory for initial '// 1678 & 'allocation in CCSX_D1AO') 1679 ENDIF 1680C 1681C-------------------------------------------- 1682C Initialize one electron density arrays. 1683C-------------------------------------------- 1684C 1685 CALL DZERO(WORK(KONEAB),NMATAB(1)) 1686 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 1687 CALL DZERO(WORK(KONEAI),NT1AMX) 1688 CALL DZERO(WORK(KONEIA),NT1AMX) 1689C 1690C----------------------------- 1691C Read Left eigen-vector. 1692C----------------------------- 1693C 1694 IOPT = 1 1695 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODDUM, 1696 * WORK(KL1AM),WORK(KEND1)) 1697C 1698C----------------------------- 1699C Read rigth eigen-vector. 1700C----------------------------- 1701C 1702 IOPT = 1 1703 CALL CC_RDRSP(RLIST,ILRNR,ISYMR,IOPT,MODDUM, 1704 * WORK(KR1AM),WORK(KEND1)) 1705C 1706C-------------------------------------------------------------- 1707C Construct <L1|[Emn,R1]|HF> contribution to DXaa and DXii. 1708C-------------------------------------------------------------- 1709C 1710 CALL CC_DXIJ(WORK(KONEIJ),WORK(KL1AM),ISYML,WORK(KR1AM),ISYMR) 1711 CALL CC_DXAB(WORK(KONEAB),WORK(KL1AM),ISYML,WORK(KR1AM),ISYMR) 1712C 1713C------------------------------- 1714C <L1|Emn|HF> Contribution 1715C------------------------------- 1716C 1717 IF (CCS.AND.(.NOT.CIS)) THEN 1718 IOPT = 1 1719 CALL CC_RDRSP(L0LIST,ILNRL0,ISYMOP,IOPT,MODDUM, 1720 * WORK(KONEAI),WORK(KEND1)) 1721 ENDIF 1722C 1723C------------------------------- 1724C Get MO coefficient matrix. 1725C------------------------------- 1726C 1727 CALL DZERO(WORK(KT1AM),NT1AMX) 1728 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 1729 * LWRK1) 1730C 1731C-------------------------------------------------------- 1732C Add the one electron density in the AO-basis. 1733C-------------------------------------------------------- 1734C 1735 CALL DZERO(AODEN,N2BST(1)) 1736 ISDEN = 1 1737 CALL CC_DENAO(AODEN,ISDEN,WORK(KONEAI),WORK(KONEAB), 1738 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 1739 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 1740C 1741 CALL QEXIT('CCSX_D1AO') 1742C 1743 END 1744c*DECK TNSRAN 1745 SUBROUTINE CC_TNSRAN(TNSR,WORK,LWORK) 1746C 1747C------------------------------------------------------------------------ 1748C 1749C Call TNSRAN and write out selected info. 1750C 1751C------------------------------------------------------------------------ 1752C 1753#include "implicit.h" 1754#include "priunit.h" 1755#include "maxorb.h" 1756#include "ccorb.h" 1757#include "iratdef.h" 1758#include "ccsdinp.h" 1759C 1760 PARAMETER (THR = 1.0D-08) 1761 DIMENSION TNSR(3,3),PVAL(3),PAXIS(3,3) 1762C 1763 CALL QENTER('CC_TNSRAN') 1764C 1765 CALL TNSRAN(TNSR,PVAL,PAXIS,ALFSQ,BETSQ,ITST,ITST2, 1766 * APAR,APEN,XKAPPA,IPAR) 1767C 1768 WRITE(LUPRI,'(/,1X,A38,F14.6)') 1769 * 'Alfa**2 Invariant: ' 1770 * //' ',ALFSQ 1771 WRITE(LUPRI,'(1X,A38,F14.6)') 1772 * 'Beta**2 Invariant: ' 1773 * //' ',BETSQ 1774 SHPAL = SQRT(ALFSQ) 1775 ANINV = SQRT(BETSQ) 1776 WRITE(LUPRI,'(/,1X,A42,F10.6,A)') 'Isotropic Property: ' 1777 * //' ',SHPAL,' a.u.' 1778 WRITE(LUPRI,'(1X,A42,F10.6,A)') 'Property anisotropy invariant:' 1779 * //' ',ANINV,' a.u.' 1780 1781 CALL QEXIT('CC_TNSRAN') 1782C 1783 END 1784C /* Deck ccmm_dia */ 1785 SUBROUTINE CCMM_DIA(DIA,DIJ,ISYML,XR1AM,ISYMR) 1786C 1787C Purpose: Calculates DIA block of quadratic density to be use 1788C for a quadratic response CC/MM calculation. It takes as input 1789C a pseudo-density of product coefficients \sum_a t^a_i b_p^a = Dij 1790C calculated in the CC_DXIJ routine. 1791C 1792C KS, Aug 2010 1793C 1794#include "implicit.h" 1795 PARAMETER(ZERO = 0.0D0, XMONE = -1.0D0, ONE = 1.0D00) 1796 DIMENSION DIA(*), DIJ(*), XR1AM(*) 1797#include "priunit.h" 1798#include "ccorb.h" 1799#include "ccsdsym.h" 1800C 1801 CALL QENTER('CCMM_DIA') 1802C 1803C--------------------------------------------------------- 1804C Add contribution. 1805C Assume ISYMR = ISYML and density is total symmetric. 1806C--------------------------------------------------------- 1807C 1808 ISYRES = MULD2H(ISYMR,ISYML) 1809 IF (ISYRES.NE.ISYMOP) CALL QUIT('CCMM_DIA require tot.sym Dens.') 1810C 1811 1812 DO 100 ISYMA = 1,NSYM 1813C 1814 ISYMI = MULD2H(ISYMA,ISYRES) 1815 ISYMK = MULD2H(ISYMA,ISYMR) 1816C 1817 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1818 KOFF2 = IMATIJ(ISYMI,ISYMK) + 1 1819 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1820C 1821 NTOTA = MAX(NVIR(ISYMA),1) 1822 NTOTI = MAX(NRHF(ISYMI),1) 1823 1824 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 1825 * ONE,XR1AM(KOFF1),NTOTA,DIJ(KOFF2),NTOTI,ONE, 1826 * DIA(KOFF3),NTOTA) 1827 1828 100 CONTINUE 1829C 1830 CALL QEXIT('CCMM_DIA') 1831C 1832 END 1833 1834