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_drv */ 20 SUBROUTINE CC_DRV(WORK,LWORK) 21C 22C 23C Ove Christiansen 14-11-1995 24C 25C 26C Purpose: Driver routine for coupled cluster modules in 27C calculation of 28C ground state energies, excitation energies, 29C oscilator strenghts, 30C frequency dependent polarizabilities and first order 31C properties without orbital relaxation. 32C 33C Pass control to; 34C 35C CCSD_ENERGY for (integral direct) calculation of 36C groundstate energies and amplitudes. 37C 38C CC_EXCI for (integral direct) calculation of 39C excitation energies and amplitudes. 40C 41C CC_FOP for integral direct calculation of first 42C order properties. 43C Also calculates the left amplitudes. 44C 45C CC_EXGR for (integral direct) calculation of excited 46C state gradient properties. 47C Currently only first order properties. 48C 49C CC_LR & Linear Response or Second-Order Properties: 50C CC_SOP polarizabilities, property gradients etc. 51C 52C CC_LRESID & Control calculation of residues and 53C CC_SOPR oscillator strenths. 54C 55C CC_HYPPOL calculation of quadratic response properties. 56C 57C CC_2HYP calculation of cubic response properties. 58C 59C CC_3HYP calculation of quartic response properties. 60C 61C CC_TPA Two-Photon Absorption strengths and moments. 62C 63C CC_TMCAL calculation of third order moments 64C 65C CC_MCDCAL calculation of B terms in MCD 66C 67C CC_LANCZOS_DRV tridiagonal matrix build and LR Lanczos 68C 69C Loop over all models specified. The order is supposed 70C to be the optimal one from a computational point of view. 71C Energy calculations are restarted from the preceding 72C calculation. Excitation vectors is also saved 73C and the excitation energy calculation is restarted 74C from preceding model solution as well. Similar for 75C response vectors in property calculations. 76C 77C (ground state) geometry optimizations using the 78C analytic gradient are done for the first (i.e. lowest 79C order) model specified in the input. 80C 81C Order: CIS,CCS,MP2,CC(2)(=CIS(D)),CC2,CCD,CCSD, 82C CCSD(T),CCSDR(T),CCSDR(1a),CCSDR(1b), 83C CC(3),CCSDR(3), CC3,CC1A,CC1B 84C (see list below for name of logical flags.) 85C 86C 87 USE PELIB_INTERFACE, ONLY: USE_PELIB 88#include "implicit.h" 89#include "priunit.h" 90#include "mxcent.h" 91#include "iratdef.h" 92#include "dummy.h" 93#include "ccorb.h" 94#include "ccsdinp.h" 95#include "ccsections.h" 96#include "ccinftap.h" 97#include "inftap.h" 98#include "ccgr.h" 99#include "cclrinf.h" 100#include "cclr.h" 101#include "ccsdsym.h" 102#include "ccfdgeo.h" 103#include "ccslvinf.h" 104#include "ccexci.h" 105#include "prpc.h" 106#include "ccfop.h" 107#include "ccfro.h" 108#include "eritap.h" 109#include "cbirea.h" 110#include "r12int.h" 111#include "grdccpt.h" 112#include "ccfield.h" 113#include "qm3.h" 114C 115Cholesky 116#include "maxorb.h" 117#include "ccdeco.h" 118Cholesky 119C 120#include "ccexcicvs.h" 121#include "ccxscvs.h" 122C 123 PARAMETER(NMODEL = 22) 124 !sonia: prevent projection in L0 125 logical lcvsbck,lionbck,lrmcorebck 126 ! 127 LOGICAL LCCSAV(0:NMODEL), LTCEXC, LPTEXC 128 LOGICAL LTBAR0, LRSPSYM, LLEFTEIG 129 LOGICAL MOLGRD, LHTF, NOAUXBSAVE 130 LOGICAL CC3RSP, CCR12RSP, CCR12LIM 131 PARAMETER (CC3RSP = .TRUE.) 132 DIMENSION WORK(LWORK) 133 CHARACTER FILPQIM*(6), FILXYM*(6) 134 PARAMETER (FILPQIM = 'CCPQ00', FILXYM = 'CC_XYM') 135 CHARACTER*3 APROXR12 136 CHARACTER*8 NAME(8) 137 DATA NAME /'CCAOIN_1','CCAOIN_2','CCAOIN_3','CCAOIN_4', 138 * 'CCAOIN_5','CCAOIN_6','CCAOIN_7','CCAOIN_8'/ 139 COMMON/SORTIO/LUAOIN(8) 140 141 INTEGER NCCEXCI2(8,3),ISYOFE2(8),ITROFE2(8) 142C 143 INTEGER ISYEXC2(MAXEXCI),IMULTE2(MAXEXCI) 144C 145 CALL QENTER('CC_DRV') 146 CALL GETTIM(TSTART,WSTART) 147C 148 IF (.NOT. (ONLYMO)) THEN 149 CALL CCTAPINI 150 WRITE (LUPRI,'(A,/)') ' ' 151 WRITE (LUPRI,'(1x,A)') 152 *'*********************************************************'// 153 *'**********************' 154 WRITE (LUPRI,'(1x,A)') 155 *'*********************************************************'// 156 *'**********************' 157 WRITE (LUPRI,'(1x,A)') 158 *'* '// 159 *' *' 160 WRITE (LUPRI,'(1x,A)') 161 *'* '// 162 *' *' 163 WRITE (LUPRI,'(1x,A)') 164 *'* START OF COUPLED CLUSTER CALCULATION '// 165 *' *' 166 WRITE (LUPRI,'(1x,A)') 167 *'* '// 168 *' *' 169 WRITE (LUPRI,'(1x,A)') 170 *'* '// 171 *' *' 172 WRITE (LUPRI,'(1x,A)') 173 *'*********************************************************'// 174 *'**********************' 175 WRITE (LUPRI,'(1x,A,/)') 176 *'*********************************************************'// 177 *'**********************' 178C 179 CALL FLSHFO(LUPRI) 180C 181 IDUM = 0 182C 183 IF (IPRINT.GE.5) WRITE(LUPRI,'(/1X,A,I15,/)') 'LWORK = ',LWORK 184 END IF 185C 186 IF (ONLYMO) THEN 187 LUSIFC = -1 188 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 189 & .FALSE.) 190 REWIND LUSIFC 191C 192 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 193 READ (LUSIFC) NSYM, NORBTS, NBAST, NLAMDS, (NRHFS(I),I=1,NSYM), 194 & (NORBS(I),I=1,NSYM), (NBAS(I),I=1,NSYM), DUM, DUM 195C 196 CALL GPCLOSE(LUSIFC,'KEEP') 197C 198 CALL CCMM_REP1(WORK,LWORK) 199C 200C Unit no. for info file set to zero 201C 202 LUMMMO = -1 203C 204 CALL QUIT('SCF orbitals and energies calculated') 205C 206 END IF 207C---------------------------------------------------- 208C read some flags from /ABAINF/ common block 209C and close SIRIFC in case of gradient requested. 210C---------------------------------------------------- 211C 212 CALL CCRDABAINF(MOLGRD) 213 IF (LUSIFC .GT. 0) CALL GPCLOSE(LUSIFC,'KEEP') 214C 215C---------------------------------------------------------------------- 216C initialize common block for Cray I/O: 217C---------------------------------------------------------------------- 218C 219 CALL INITWIO() 220C 221C---------------------------------------------------------------------- 222C Open files needed open throughout entire CC calculation. 223C---------------------------------------------------------------------- 224C 225 LURES = -1 226 CALL GPOPEN(LURES,'CC_RES','UNKNOWN',' ','FORMATTED',IDUMMY, 227 & .FALSE.) 228 REWIND(LURES) 229 230 LUPRPC = -1 231 CALL GPOPEN(LUPRPC,'CC_PRPC','UNKNOWN',' ','FORMATTED',IDUMMY, 232 & .FALSE.) 233 REWIND(LUPRPC) 234 235 LUIAJB = -1 236 CALL GPOPEN(LUIAJB,'CCSD_IAJB','UNKNOWN',' ','UNFORMATTED',IDUMMY, 237 & .FALSE.) 238 REWIND(LUIAJB) 239C 240C--------------------------------- 241C Initialize energy variables. 242C--------------------------------- 243C 244 ECCGRS = 0.0D0 245 OMECCX = 0.0D0 246 ECCXST = 0.0D0 247C 248C----------------------------------------------------------------------* 249C Save information on the number of states etc. for subsequent 250C calls in numerical differentiations. 251C----------------------------------------------------------------------* 252C 253 NEXCI2 = NEXCI 254 DO IMULT = 1, 3, 2 255 DO ISYM =1, 8 256 NCCEXCI2(ISYM,IMULT) = NCCEXCI(ISYM,IMULT) 257 IF (IMULT.EQ.1) ISYOFE2(ISYM) = ISYOFE(ISYM) 258 IF (IMULT.EQ.3) ITROFE2(ISYM) = ITROFE(ISYM) 259 ENDDO 260 ENDDO 261 DO IEX=1,MAXEXCI 262 ISYEXC2(IEX) = ISYEXC(IEX) 263 IMULTE2(IEX) = IMULTE(IEX) 264 ENDDO 265 NPRPC = 0 266 NPRMI = 0 267C 268C--------------------------------------------------- 269C Set some logicals needed for R12 calculations: 270C--------------------------------------------------- 271C 272 R12CAL = R12CAL .OR. LMULBS 273 274 ! request that some intermediates needed for CC2-R12 275 ! are calculated during MP2-R12 calculation 276 CC2R12INT = R12CAL .AND. CC2 277 CCSDR12INT = R12CAL .AND. (CCSD.OR.CCPT.OR.CCP3.OR.CC3) 278 279 ! switch on R12 flag used in CC routines 280 CCR12 = R12CAL 281 282 ! CCR12 not possible as CCS: 283 IF (CCR12.AND.CCS) CALL QUIT('CCS not reasonable with R12') 284C 285C------------------------------------- 286C Set logicals: 4,6,8 is obsolete. 287C------------------------------------- 288C 289 LCCSAV(0) = CIS 290 LCCSAV(1) = CCS 291 LCCSAV(2) = MP2 292 LCCSAV(3) = CCP2 293 LCCSAV(4) = .FALSE. 294 LCCSAV(5) = CC2 295 LCCSAV(6) = .FALSE. 296 LCCSAV(7) = CCD 297 LCCSAV(8) = .FALSE. 298 LCCSAV(9) = CCSD 299 LCCSAV(10) = CCPT 300 LCCSAV(11) = CCRT 301 LCCSAV(12) = CCR1A 302 LCCSAV(13) = CCR1B 303 LCCSAV(14) = CCP3 304 LCCSAV(15) = CCR3 305 LCCSAV(16) = CC3 306 LCCSAV(17) = CC1A 307 LCCSAV(18) = CC1B 308! FRAN/SONIA 309 LCCSAV(19) = RCCD 310 LCCSAV(20) = DRCCD 311 LCCSAV(21) = RTCCD 312! AMT 313 LCCSAV(22) = DCPT2 314 315 APROXR12 = ' ' 316C 317 LTCEXC = .FALSE. 318C 319 MMOD = 0 320 DO I = 0, NMODEL 321 IF (LCCSAV(I)) MMOD = MMOD + 1 322 END DO 323C 324C------------------------------------------ 325C Call the CCSD initialization routine. 326C------------------------------------------ 327C 328 ISYMOP = 1 329 LABEL = 'TRCCINT ' 330 IF (LMULBS) THEN 331 NOAUXBSAVE = NOAUXB 332 NOAUXB = .TRUE. 333 END IF 334 CALL CCSD_INIT1(WORK,LWORK) 335C 336C----------------------------------------------------------------- 337C flag for ground state zeroth-order Lagrange multiplier: 338C Lanczos added. Sonia 339C----------------------------------------------------------------- 340 LTBAR0 = (CCLRSD.OR.CCLR.OR.CCQR.OR.CCCR.OR.CC4R.OR. 341 & CC5R.OR.CCQR2R.OR.CCEXGR.OR.CCTM.OR. 342 & CCMCD.OR.CCEXLR.OR.MOLGRD.OR.CCSLV.OR.CCDERI.OR. 343 & CCLRLCZ.OR. 344 & CCOPA.OR.CCTPA.OR.CCXOPA.OR.USE_PELIB()) 345 346C----------------------------------------------------------------- 347C flag for kappa-bar-0 Lagrange multiplier: 348C----------------------------------------------------------------- 349 RELORB = MOLGRD .OR. CCDERI .OR. RELORB 350 351C----------------------------------------------------------------- 352C consistency check: no gradient for solvent... 353C----------------------------------------------------------------- 354 IF ((MOLGRD.OR.CCDERI) .AND. (CCSLV.OR.USE_PELIB())) THEN 355 WRITE(LUPRI,*) 'No analytic derivatives in solvent available!' 356 IF (MOLGRD) THEN 357 CALL QUIT("No optimiz. with anal. gradient in solvent.") 358 END IF 359 IF (CCDERI) THEN 360 WRITE(LUPRI,*) 'CCDERI flag will be ignored.' 361 CCDERI = .FALSE. 362 END IF 363 END IF 364 365C----------------------------------------------------------------- 366C flag for deleting AOTWOINT integral file: 367C KEEPAOTWO = 0 --> delete AOTWOINT in CCSD_SORTAO 368C KEEPAOTWO < 2 --> delete AOTWOINT in CC_GRAD 369C KEEPAOTWO >= 2 --> keep AOTWOINT file until the end 370C----------------------------------------------------------------- 371 IF (RELORB .AND. MMOD.LE.1) KEEPAOTWO = MAX(KEEPAOTWO,1) 372 IF (RELORB .AND. MMOD.GT.1) KEEPAOTWO = MAX(KEEPAOTWO,2) 373 374C----------------------------------------------------------------- 375C flag for response calculation (-> execute CCRSPSYM routine): 376C Lanczos added. Sonia 377C----------------------------------------------------------------- 378 LRSPSYM = (CCLRSD.OR.CCLR.OR.CCQR.OR.CCCR.OR.CC4R.OR. 379 & CC5R.OR.CCQR2R.OR.CCEXGR.OR.CCTM.OR. 380 & CCLRLCZ.OR. 381 & CCMCD.OR.CCEXLR.OR.CCOPA.OR.CCTPA.OR.CCXOPA) 382 383 RSPIM = RSPIM .OR. LRSPSYM .OR. LTBAR0 384 385C----------------------------------------------------------------- 386C set flag for CCR12 response beyond excitation energies 387C (controls the calculation of certain integrals) 388C----------------------------------------------------------------- 389 CCR12RSP = CCR12 .AND. (LRSPSYM .OR. NONHF .OR. CCFOP) .AND. 390 & .NOT. R12PRP 391 IF (CCR12RSP.AND.(IANCC2.NE.1)) THEN 392 WRITE(LUPRI,*) 'CCR12RSP = ',CCR12RSP 393 WRITE(LUPRI,*) 'IANR12 = ',IANCC2 394 CALL QUIT('Sorry, only Ansatz 1 implemented '// 395 & 'for CC response') 396 END IF 397 398C----------------------------------------------------------------- 399C flag for left eigenvector calculation: 400C----------------------------------------------------------------- 401 LLEFTEIG = CCOPA .OR. CCTPA .OR. CCXOPA .OR. 402 & CCMCD .OR. CCTM .OR. CCEXLR 403 404C----------------------------------------------------------------- 405C LPTEXC - to control if any perturbative triples corrections 406C to excitation energies and thus calculated left vectors. 407C----------------------------------------------------------------- 408C 409 LPTEXC = .FALSE. 410 IF (CCRT.OR.CCR3.OR.CCR1A.OR.CCR1B) LPTEXC = .TRUE. 411C 412C---------------------------------------------------------------- 413C CCR12LIM - controls if intermediates for CCR12 left transf. 414C are needed 415C---------------------------------------------------------------- 416C 417 CCR12LIM = CCR12.AND.(CCR12RSP.OR.LLEFTEIG.OR.LHTR) 418C 419C------------------------------------------ 420C Sort AO-integrals into distributions. 421C------------------------------------------ 422C 423 IF ((.NOT. CHOINT) .AND. (.NOT. DIRECT)) THEN 424 CALL CCSD_SORTAO(WORK,LWORK) 425 ENDIF 426C 427C----------------------------------------------------------- 428C Calculate the ( ia | jb ) integrals and write to disk. 429C----------------------------------------------------------- 430C 431 IF ((.NOT. CHOINT) .AND. (.NOT.LISKIP)) THEN 432 IF (LMULBS.AND.MMOD.LE.1.AND.MP2.AND.(.NOT.NATVIR).AND. 433 * (.NOT.R12PRP)) THEN 434 ! if only MP2-R12 the (IA|JB) integrals are not needed 435 ! on file LUIAJB 436 CONTINUE 437 ELSE IF (LMULBS.AND.HERDIR) THEN 438 CALL QUIT('CC-R12 with HERDIR not implemented') 439 ELSE 440 KT1AM = 1 441 KT2AM = KT1AM + NT1AMX 442 KEND = KT2AM + NT2AMX 443 LWRK = LWORK - KEND 444 IF (LWRK .LT. 0) CALL QUIT ('Insufficient memory in CC_DRV') 445 CALL DZERO(WORK(KT1AM),NT1AMX) 446 LHTF = .FALSE. 447 CALL CCSD_IAJB(WORK(KT2AM),WORK(KT1AM),LHTF, 448 * .FALSE.,.FALSE.,WORK(KEND),LWRK) 449 REWIND(LUIAJB) 450 CALL WRITI(LUIAJB,IRAT*NT2AM(ISYMOP),WORK(KT2AM)) 451 ENDIF 452 ENDIF 453C 454C-------------------------------------------------------------------- 455C Calculate semi-natural virtual orbitals from MP2-Guess for R12: 456C-------------------------------------------------------------------- 457C 458 IF ((.NOT.LISKIP).AND.NATVIR) THEN 459 CALL CC_R12NO(WORK(KT1AM),WORK(KT2AM),WORK(KEND),LWRK) 460 END IF 461C 462 IF (LMULBS) NOAUXB = NOAUXBSAVE 463C 464C---------------------- 465C Loop over models. 466C---------------------- 467C 468 INTTR = .FALSE. 469 IMOD = 0 470C 471 DO 100 I = 0, NMODEL 472C 473 CALL CC_FALSE() 474C 475 IF (.NOT.LCCSAV(I)) THEN 476C 477 IF ((I .EQ. NMODEL).AND.(IMOD.EQ.0)) THEN 478 WRITE(LUPRI,'(1x,A)') 'No CC model specified - '// 479 * 'default is CCSD' 480 CCSD = .TRUE. 481 ELSE 482 GOTO 100 483 ENDIF 484C 485 ENDIF 486C 487 IF (I.EQ.0) CIS = LCCSAV(0) 488 IF (I.EQ.0) CCS = LCCSAV(0) 489 IF (I.EQ.1) CCS = LCCSAV(1) 490 IF (I.EQ.2) MP2 = LCCSAV(2) 491 IF (I.EQ.3) CCP2 = LCCSAV(3) 492 IF (I.EQ.4) GOTO 100 493 IF (I.EQ.4) CC2 = LCCSAV(4) 494 IF (I.EQ.5) CC2 = LCCSAV(5) 495 IF (I.EQ.6) GOTO 100 496 IF (I.EQ.7) CCD = LCCSAV(7) 497 IF (I.EQ.8) GOTO 100 498 IF (I.EQ.9) CCSD = LCCSAV(9) 499 IF (I.EQ.10) CCPT = LCCSAV(10) 500 IF (I.EQ.11) CCRT = LCCSAV(11) 501 IF (I.EQ.12) CCR1A = LCCSAV(12) 502 IF (I.EQ.13) CCR1B = LCCSAV(13) 503 IF (I.EQ.14) CCP3 = LCCSAV(14) 504 IF (I.EQ.15) CCR3 = LCCSAV(15) 505 IF (I.EQ.16) CC3 = LCCSAV(16) 506 IF (I.EQ.17) CC1A = LCCSAV(17) 507 IF (I.EQ.18) CC1B = LCCSAV(18) 508!FRAN/SONIA 509 IF (I.EQ.19) RCCD = LCCSAV(19) 510 IF (I.EQ.20) DRCCD = LCCSAV(20) 511 IF (I.EQ.21) RTCCD = LCCSAV(21) 512!AMT 513 IF (I.EQ.22) DCPT2 = LCCSAV(22) 514C 515 IF (CCS .AND. LCCSAV(3)) GOTO 100 516 IF (MP2 .AND. LCCSAV(3)) GOTO 100 517 IF (CCP2 ) MP2 = .TRUE. 518C IF (CCSD .AND. LCCSAV(10)) GOTO 100 519 IF (CCSD .AND. LCCSAV(11)) GOTO 100 520 IF (CCSD .AND. LCCSAV(12)) GOTO 100 521 IF (CCSD .AND. LCCSAV(13)) GOTO 100 522 IF (CCSD .AND. LCCSAV(14) .AND.(.NOT.LCCSAV(15))) GOTO 100 523 IF (CCP3 .AND. LCCSAV(15)) GOTO 100 524C 525 IMOD = IMOD + 1 526C 527C ----------------------------------------------------------- 528C Save integral transformation (if not triples) and 529C restart from previous cc amplitudes and excitation vectors. 530C ----------------------------------------------------------- 531 IF (IMOD.GT.1) THEN 532 STOLD = .TRUE. 533 IF ((I.GT.3).AND..NOT.(LCCSAV(1).AND.(IMOD.EQ.2))) 534 * CCRSTR = .TRUE. 535 ENDIF 536 537 ! in geometry optimizations we can restart from the 538 ! amplitudes of the previous points 539 IF (MOLGRD) CCRSTR = .TRUE. 540 541 IF (I .GT. 9) INTTR = .TRUE. 542C 543 !IF ((I.GT.15).AND.(.NOT.CCSD)) THEN 544 !FRAN/SONIA, 19,20 and 21 are not CCSDT methods! 545 IF (((I.GT.15).AND.(I.LT.19)).AND.(.NOT.CCSD)) THEN 546 547 CCSDT = .TRUE. 548 ENDIF 549C 550 IF (CCR12 .AND. .NOT.(CCS .OR. MP2 .OR. CC2.OR. 551 & CCSD.OR.CCPT.OR.CCP3.OR.CC3)) THEN 552 CALL QUIT('R12 not yet implemented for this CC model') 553 END IF 554 555C ----------------------------------------------------------- 556C in geometry optimization runs compute the gradient only 557C for the first model 558C ----------------------------------------------------------- 559 IF (IMOD.GT.1 .AND. MOLGRD) GOTO 100 560C 561C---------------------------------------------- 562C ************************************** 563C ***** Zeroth-order CC amplitudes ***** 564C ************************************** 565C---------------------------------------------- 566C 567 IF (I .GT. 13) LTCEXC = .FALSE. 568 IF (LTCEXC) GO TO 100 569 IF (CCP2) CCS = .FALSE. 570 IF (CCR3) CCP3 = .TRUE. 571C 572C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 573C Start loop over R12 ansaetze and approximations 574C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 575C 576 IANR12 = IANCC2 577 IAPR12 = IAPCC2 578C 579 WRITE(LUPRI,'(//A,I3)') ' CCR12 ANSATZ = ',IANR12 580 WRITE(LUPRI,'(/A,I3/)') ' CCR12 APPROX = ',IAPR12 581 CALL FLSHFO(LUPRI) 582 583C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 584C 1. Start loop over different solvents. 585C 2. Start loop for self-consistent CC solvent 586C solution. 587C Do at least one(which could be vacuum). 588C SLV98,OC 589C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 590C 591 NSLV = MAX(1,NCCSLV) 592C 593 DO 9999 ISLV = 1, NSLV 594C 595 LMAXCU = LMAXCC(ISLV) 596 NLMCU = (LMAXCU+1)*(LMAXCU+1) 597 EPSTCU = EPSTCC(ISLV) 598 EPOPCU = EPOPCC(ISLV) 599 RCAVCU = RCAVCC(ISLV) 600 LSTBTR = .FALSE. 601 LSLECVG= .FALSE. 602 LSLTCVG= .FALSE. 603 LSLLCVG= .FALSE. 604 ICCSLIT= 0 605C 606C Spaghetti GOTO for finding CCSCRF 607C 608 31 CONTINUE 609C 610C ! no energy in gradient step for geometry optimization 611 IF (.NOT.MOLGRD) THEN 612 613 CALL CCSD_ENERGY(WORK,LWORK,APROXR12,CCR12RSP,CCR12LIM) 614c 615c write(lupri,*) 'after CCSD_ENERGY: APROXR12 = ',APROXR12 616c 617 618 END IF 619C 620C switch off flags for sort and transformation 621C done in CCSD_ENERGY 622C 623 INTTR = .FALSE. 624C 625 IF (CCP2) CCS = .TRUE. 626C KS: For CCSDR(3)/MM model we use vacuum triples amplitudes, i.e. we skip the 627C call to CC_FOP by a dirty goto 628 IF (CCR3 .AND. (CCSLV.OR.USE_PELIB())) THEN 629 LSTBTR=.TRUE. 630 GOTO 32 631 ENDIF 632C 633C-------------------------------------- 634C **************************** 635C ***** CC Lambda vector ***** 636C **************************** 637C-------------------------------------- 638C 639 640 IF ( (CCFOP.OR.LTBAR0) .AND. 641 & (.NOT. (CCR3.OR.CCR1A.OR.CCR1B.OR.CCRT.OR.CCR3))) THEN 642C 643 IF (CCR12.AND.MP2.AND..NOT.R12PRP) THEN 644 NWARN = NWARN + 1 645 WRITE(LUPRI,'(//A/A)') 646 & 'WARNING: Please specify R12PRP in your input...', 647 & 'WARNING: CCFOP will be ignored...' 648 ELSEIF (.NOT. (CCSDT .AND. (.NOT. CC3))) THEN 649 !................................... 650 !SONIA: I am not sure but let us try 651 !................................... 652 !IF ( I .EQ. 10 ) CALL CC3_FILOP() 653 !................................... 654 IF (R12PRP) THEN 655 DO IPDD = 1,5 656 IF (IPDD .EQ.2.OR.IPDD.EQ.3.OR.IPDD.EQ.5)THEN 657 CALL CC_FOP(IPDD,WORK,LWORK,APROXR12) 658 ENDIF 659 ENDDO 660 ELSE 661 lcvsbck = LCVSEXCI 662 lionbck = LIONIZEXCI 663 lrmcorebck = LRMCORE 664 LCVSEXCI = .false. 665 LIONIZEXCI = .false. 666 LRMCORE = .false. 667 CALL CC_FOP(IDUMMY,WORK,LWORK,APROXR12) 668 LCVSEXCI = lcvsbck 669 LIONIZEXCI = lionbck 670 LRMCORE = lrmcorebck 671 ENDIF 672 ENDIF 673C 674 ENDIF 675C 676C------------------------------------------------------------ 677C SLV98,OC 678C End loop for self-consistent CC solvent solution. 679C Force restart after first run. 680C After CC SCRF convergence then left 681C transformations are no longer static. 682C------------------------------------------------------------ 683C 684 IF (CCSLV.OR.USE_PELIB()) THEN 685 CCRSTR = .TRUE. 686 IF (LSLECVG.AND.LSLTCVG.AND.LSLLCVG) THEN 687 LSTBTR = .TRUE. 688 ELSE 689 GOTO 31 690 ENDIF 691 ENDIF 692C 693 32 CONTINUE 694 IF ((CCSLV.OR.USE_PELIB()) .AND. DISCEX) THEN 695C IF (CCEXCI .OR. CCEXGR .OR. LRSPSYM .OR. CCLRSD .OR. 696C * CCQR2R .OR. CCSM .OR. CCLR .OR. CCQR ) THEN 697 IF (.NOT. (CCMM.OR.USE_PELIB())) THEN 698 CALL CCSL_XIETA(WORK,LWORK) 699 ELSE IF (CCMM.OR.USE_PELIB()) THEN 700 CALL AROUND('ENTERING CCMM_XIETA') 701 CALL CCMM_XIETA(WORK,LWORK) 702 END IF 703C END IF 704 END IF 705C 706C----------------------------------------------------------------------* 707C 708C CCDERI -- gradient calculation forced in CC input 709C 710C MOLGRD -- gradient required from DALTON (via /ABAINF/) 711C some short cuts are made in the latter case: 712C 1) the energy calculation and all property/excitation 713C energy calculations are saved because they were 714C done in a previous call to the CC program 715C 2) we do the optimization (and thus the gradient) 716C for the first (i.e. lowest order) model 717c 718C----------------------------------------------------------------------* 719C 720 IF (CCDERI .OR. MOLGRD) THEN 721 CALL CC_GRADIENT(WORK,LWORK) 722 IGRDCCPT = IGRDCCPT + 1 723 724 IF (MOLGRD) GOTO 100 725 END IF 726C 727C--------------------------------------------------------------------- 728C If triples corrections have been calculated goto end of loop. 729C or do all triple corrections in one call. 730C But if both ccsd(t) and cc(3) ground state energies are wanted 731C wait to next time. 732C--------------------------------------------------------------------- 733C 734 IF (LTCEXC) GO TO 90 735C 736 IF ((I.GE.10).AND.(I.LE.13)) THEN 737C 738 LTCEXC = .TRUE. 739 CCRT = LCCSAV(11) 740 CCR1A = LCCSAV(12) 741 CCR1B = LCCSAV(13) 742C 743 ENDIF 744C 745 IF ( CCR3 ) CCT = .TRUE. 746C 747C--------------------------------------------- 748C *********************************** 749C ***** CC Excitation energies ***** 750C *********************************** 751C--------------------------------------------- 752C 753 IF (CCEXCI) THEN 754C 755 IF (((.NOT. MP2).OR.CCP2).AND.(.NOT.(CCD.OR.CCPT))) THEN 756C 757 IF ((.NOT. LHTR).AND.(.NOT.RESKIP)) THEN 758 CALL CC_EXCI(WORK,LWORK,'RE ',APROXR12) 759 ENDIF 760C 761 IF (CCSPIC.AND.CCS) CALL CC_REDEIG(WORK,LWORK,OMPCCS) 762 IF (CC2PIC.AND.CC2) CALL CC_REDEIG(WORK,LWORK,OMPCC2) 763 IF (CCSDPI.AND.CCSD) CALL CC_REDEIG(WORK,LWORK,OMPCCSD) 764C 765 IF (LHTR.OR.CCLRSD.OR.CCQR2R.OR.LPTEXC.OR.CCP2.OR. 766 * LLEFTEIG.OR.CCEXGR) THEN 767 IF ((I .LT. 17).AND.(.NOT.LESKIP)) THEN 768 CALL CC_EXCI(WORK,LWORK,'LE ',APROXR12) 769 ELSE 770 WRITE(LUPRI,*) 'No triples left excitation energies ' 771 ENDIF 772 ENDIF 773 774 IF ((.NOT. LHTR).AND.(.NOT.RESKIP)) THEN 775 CALL CC_MOLDEN_NTO(WORK,LWORK) 776 ENDIF 777 ELSE 778 WRITE(LUPRI,'(/,1X,A)') 779 & ' No MP2 and CCD excitation energies' 780 ENDIF 781 782C 783 ENDIF 784C 785C------------------------------------------------------------ 786C Settings for numerical gradient (driven from dalton) 787C------------------------------------------------------------ 788C 789 CALL CC_FDGD 790C 791C------------------------------------------------ 792C The rest is not implemented for triples. 793C------------------------------------------------ 794C 795 IF (I.GT.9 .AND. I.NE.16) GOTO 90 796C 797C------------------------------------------------------------ 798C **************************************************** 799C ***** general set up for response calculations ***** 800C **************************************************** 801C------------------------------------------------------------ 802C 803 IF (LRSPSYM) THEN 804 Call CCRSPSYM(MOLGRD,WORK,LWORK) 805 END IF 806C 807C--------------------------------------------------------------- 808C precalculate expectatation values & eff. Fock matrices: 809C--------------------------------------------------------------- 810C 811 IF (CCS .or. CC2 .or. CCSD) THEN 812 I1DXORD = 0 813 IOPREL = 0 814 CALL CC_GRAD2(I1DXORD,WORK,LWORK) 815C CALL CCEFFFOCK(I1DXORD,IOPREL,WORK,LWORK) 816 END IF 817 818C 819C----------------------------------------------- 820C *************************************** 821C ***** Solve CC-response equations.***** 822C *************************************** 823C----------------------------------------------- 824C 825 !IF (LRSPSYM) THEN 826 !Do not enter if Lanczos calculation. Sonia 827 !IF ((LRSPSYM).and.(.NOT.CCLRLCZ)) THEN 828 !reactivate for Excited state calculation 829 IF (LRSPSYM) THEN 830 if (LXSCVS.or.LXRMCORE) then 831 !transfer info for projection 832 call cc_cvs_interface(MSYM) 833 end if 834 CALL CC_SOLRSP(WORK,LWORK,APROXR12) 835 END IF 836C 837C----------------------------------------------------- 838C ********************************************* 839C ***** Excited state gradient properties ***** 840C ********************************************* 841C----------------------------------------------------- 842C 843 IF (CCEXGR) THEN 844 IF (.NOT. (CCSLV.OR.USE_PELIB())) THEN 845 IF ((CCS.or.CC2.or.CCSD.or.CC3).and.(.not.CCR12)) THEN 846 CALL CC_EXGR(WORK,LWORK) 847 ELSE 848 WRITE(LUPRI,*) 'Excited state gradient properties'// 849 & 'are not implemented for the requested model.' 850 END IF 851 ELSE 852 WRITE(LUPRI,*) 'Solvent model not implemented for'// 853 & 'excited state gradients.' 854 ENDIF 855 ENDIF 856C 857C----------------------------------------------- 858C *************************************** 859C ***** CC linear response residues ***** 860C *************************************** 861C----------------------------------------------- 862C 863 IF (CCLRSD) THEN 864 IF ((CCS .or. CC2 .or. CCSD) .and. (.not. CCR12)) THEN 865 IF (OLDLRS) THEN 866 CALL CC_LRESID(WORK,LWORK) 867 ELSE 868 CALL CC_SOPR(WORK,LWORK) 869 END IF 870 ELSE 871 WRITE(LUPRI,*) 'Linear response residues'// 872 & ' are not implemented for the requested model.' 873 END IF 874 END IF 875 876 IF (CCOPA) THEN 877 IF ((CCS .or. CC2 .or. CCSD .or. (CC3.AND.CC3RSP)) 878 & .and. (.not. CCR12)) THEN 879 CALL CC_OPA(WORK,LWORK) 880 ELSE 881 WRITE(LUPRI,*) 'One-photon absorption strengths'// 882 & ' are not implemented for the requested model.' 883 END IF 884 END IF 885C 886C-------------------------------------------------------- 887C ************************************************ 888C **** CC quadratic response second residues ***** 889C ************************************************ 890C-------------------------------------------------------- 891C 892 IF (CCQR2R) THEN 893 IF ((CCS .or. CC2 .or. CCSD) .and. (.not. CCR12)) THEN 894 CALL CC_QR2RSD(WORK,LWORK) 895 ELSE 896 WRITE(LUPRI,*) 'Quadratic response residues'// 897 & ' are not implemented for the requested model.' 898 END IF 899 END IF 900 901 IF (CCXOPA) THEN 902 IF ((CCS .or. CC2 .or. CCSD .or. (CC3.AND.CC3RSP)) 903 & .and. (.not. CCR12)) THEN 904 write(lupri,*) "SONIA cc_drv: CALLING CC_eom_XOPA" 905 if (.false.) then 906 CALL CC_XOPA(WORK,LWORK) 907 else 908 CALL CC_eom_XOPA(WORK,LWORK) 909 end if 910 ELSE 911 WRITE(LUPRI,*) 'One-photon absorption strengths'// 912 & ' are not implemented for the requested model.' 913 END IF 914 END IF 915C 916C------------------------------------------------------- 917C *********************************************** 918C ***** CC two-photon transition moments ***** 919C *********************************************** 920C------------------------------------------------------- 921C 922 IF (CCTPA) THEN 923 IF ((CCS .or. CC2 .or. CCSD .or. (CC3.AND.CC3RSP)) 924 & .and. (.not. CCR12)) THEN 925 CALL CC_TPA(WORK,LWORK) 926 ELSE 927 WRITE(LUPRI,*) 'two-photon absorption strengths'// 928 & ' are not implemented for the requested model.' 929 END IF 930 ENDIF 931C 932C------------------------------------------------------ 933C ********************************************** 934C ***** CC third-order transition moments ***** 935C ********************************************** 936C------------------------------------------------------ 937C 938 IF (CCTM) THEN 939 IF ((CCS .or. CC2 .or. CCSD) .and. (.not. CCR12) 940 & .AND. (.NOT. (CCSLV.OR.USE_PELIB()))) THEN 941 CALL CC_TMCAL(WORK,LWORK) 942 ELSE 943 WRITE(LUPRI,*) 'third-order transition moments'// 944 & ' are not implemented for the requested model.' 945 END IF 946 ENDIF 947C 948C------------------------------------------------------- 949C *********************************************** 950C ***** CC magnetic circular dichroism ***** 951C *********************************************** 952C------------------------------------------------------- 953C 954 IF (CCMCD) THEN 955 IF ((CCS .or. CC2 .or. CCSD) .and. (.not. CCR12) 956 & .AND. (.NOT. (CCSLV.OR.USE_PELIB()))) THEN 957 CALL CC_MCDCAL(WORK,LWORK) 958 ELSE 959 WRITE(LUPRI,*) 'magnetic circular dichroism'// 960 & ' is not implemented for the requested model.' 961 END IF 962 ENDIF 963C 964C----------------------------------------- 965C ******************************* 966C ***** CC Polarizabilities ***** 967C ******************************* 968C----------------------------------------- 969C 970 IF ( CCLR .AND. NBLRFR.GT.0) THEN 971C 972 IF (CIS .OR. CCS .OR. CC2 .OR. CCSD .OR. (CC3.AND.CC3RSP)) 973 & THEN 974 IF ((OLDLR .AND. .NOT.CC3) .OR. CIS) THEN 975 CALL CC_LR(WORK,LWORK) 976 ELSE 977 CALL CC_SOP(WORK,LWORK) 978 END IF 979 ELSE 980 WRITE(LUPRI,'(/,1X,A)') 981 & ' No polarizabilities for present model' 982 ENDIF 983C 984 ENDIF 985C 986C----------------------------------------- 987C ***************************** 988C ***** CC Cauchy Moments ***** 989C ***************************** 990C----------------------------------------- 991C 992 IF ( CCLR .AND. CAUCHY) THEN 993 IF ((CCS .or. CC2 .or. CCSD .OR. (CC3.AND.CC3RSP)) 994 & .AND. (.NOT.CCR12) .AND. (.NOT. CCSLV) .AND. 995 & (.NOT.USE_PELIB()) ) THEN 996 CALL CC_CAUCHY(WORK,LWORK) 997 ELSE 998 WRITE(LUPRI,'(/,1X,2A)') 'Cauchy moments', 999 & ' are not yet implemented for the requested model' 1000 ENDIF 1001 ENDIF 1002C 1003C---------------------------------------------------- 1004C ******************************************** 1005C ***** CC excited state linear response ***** 1006C ******************************************** 1007C---------------------------------------------------- 1008C 1009 IF (CCEXLR) THEN 1010 IF ((CCS .or. CC2 .or. CCSD .or. (CC3.AND.CC3RSP)) 1011 & .and. (.not. CCR12) .AND. (.NOT. CCSLV) 1012 & .AND. (.NOT. USE_PELIB())) THEN 1013 CALL CC_EXLR(WORK,LWORK) 1014 ELSE 1015 WRITE(LUPRI,*) 'Excited state linear response'// 1016 & ' is not implemented for the requested model.' 1017 END IF 1018 END IF 1019C 1020C-------------------------------------------------- 1021C ****************************************** 1022C ***** CC first Hyperpolarizabilities ***** 1023C ****************************************** 1024C-------------------------------------------------- 1025C 1026 IF ( CCQR ) THEN 1027C 1028 IF (CCS .or. CC2 .or. CCSD .or. (CC3.AND.CC3RSP)) THEN 1029 CALL CC_HYPPOL(WORK,LWORK) 1030 ELSE 1031 WRITE(LUPRI,'(/,1X,2A)') 'first hyperpolarizabilities', 1032 & ' are not yet implemented for the requested model' 1033 ENDIF 1034C 1035 ENDIF 1036C 1037C--------------------------------------------------- 1038C ******************************************* 1039C ***** CC second Hyperpolarizabilities ***** 1040C ******************************************* 1041C--------------------------------------------------- 1042C 1043 IF ( CCCR ) THEN 1044C 1045 IF ( (CCS .or. CC2 .or. CCSD .or. (CC3.AND.CC3RSP)) 1046 & .AND. (.NOT. CCSLV) .AND. (.NOT. USE_PELIB()) ) THEN 1047 CALL CC_2HYP(WORK,LWORK) 1048 ELSE 1049 WRITE(LUPRI,'(/,1X,A)') 1050 & 'requested model not yet implemented' 1051 ENDIF 1052C 1053 ENDIF 1054C 1055C-------------------------------------------------- 1056C ****************************************** 1057C ***** CC third Hyperpolarizabilities ***** 1058C ****************************************** 1059C-------------------------------------------------- 1060C 1061 IF ( CC4R ) THEN 1062C 1063 IF ( (CCS .or. CC2 .or. CCSD) .AND. (.NOT. CCSLV) .AND. 1064 & (.NOT. USE_PELIB()) ) THEN 1065 CALL CC_3HYP(WORK,LWORK) 1066 ELSE 1067 WRITE(LUPRI,'(/,1X,A)') 1068 & 'requested model not yet implemented' 1069 ENDIF 1070C 1071 ENDIF 1072C 1073C--------------------------------------------------- 1074C ******************************************* 1075C ***** CC fourth Hyperpolarizabilities ***** 1076C ******************************************* 1077C--------------------------------------------------- 1078C 1079 IF ( CC5R ) THEN 1080C 1081 IF ( (CCS .or. CC2 .or. CCSD) .AND. (.NOT. CCSLV) .AND. 1082 & (.NOT. USE_PELIB()) ) THEN 1083 CALL CC_4HYP(WORK,LWORK) 1084 ELSE 1085 WRITE(LUPRI,'(/,1X,A)') 1086 & 'requested model not yet implemented' 1087 ENDIF 1088C 1089 ENDIF 1090C--------------------------------------------------------- 1091C ************************************************* 1092C ***** CC LANCZOS: damped response via Lanczos *** 1093C ************************************************* 1094C--------------------------------------------------------- 1095 IF (CCLRLCZ) THEN 1096 IF ((.not. CCR12).AND. (.NOT. CCSLV).AND. 1097 & (.NOT. USE_PELIB())) THEN 1098 write(lupri,*)'CCLRLANCZOS: building the '// 1099 & 'Lanczos Chain' 1100 !if (.false.) then 1101 !IF (REDMEML) then 1102 ! if (Debug_Sym) then 1103 CALL CC_LANCZOS_DRV(WORK,LWORK,APROXR12) 1104 ! else 1105 ! CALL CC_LANCZOS_drv2(WORK,LWORK,APROXR12) 1106 ! end if 1107 !ELSE 1108 ! CALL CC_LANCZOS_drv(WORK,LWORK) 1109 !END IF 1110 !else 1111 ! WRITE(LUPRI,*) 'Lanczos Using ATRR' 1112 ! call cc_lanczos_drv3(work,lwork) 1113 !end if 1114 ELSE 1115 WRITE(LUPRI,*) 'Lanczos Damped CC response theory'// 1116 & ' is not implemented for the requested model.' 1117 END IF 1118 ENDIF 1119C 1120C-------------------------------------------------- 1121C ****************************************** 1122C ***** B,F,C,D... matrix test calls ***** 1123C ****************************************** 1124C-------------------------------------------------- 1125C 1126 IF ( .FALSE. ) THEN 1127C 1128 IF (CCS .or. CC2 .or. CCSD) THEN 1129C !NOTE: Only CC_FTSTNEW, CC_BTST are compatible with CC-R12! 1130C CALL CC_BTST(WORK,LWORK,APROXR12) 1131C CALL CC_CTST(WORK,LWORK) 1132C CALL CC_DTST(WORK,LWORK) 1133C CALL CC_FTSTNEW(WORK,LWORK,APROXR12) 1134C CALL CC_GTST(WORK,LWORK) 1135C CALL CC_HTST(WORK,LWORK) 1136C CALL CC_XETST(WORK,LWORK) 1137C CALL CC_FBTST(WORK,LWORK) 1138C CALL CC_AATST(WORK,LWORK) 1139C CALL CC_BATST(WORK,LWORK) 1140 ELSE 1141 WRITE(LUPRI,'(/,1X,A)') 1142 & 'requested model not yet implemented' 1143 ENDIF 1144C 1145 ENDIF 1146C 1147 90 CONTINUE 1148C 1149C SLV98,OC End of solvent loop 1150C 1151 9999 CONTINUE 1152C 1153C------------------------------------------------- 1154C Close files used in triples calculations. 1155C------------------------------------------------- 1156C 1157!SONIA SONIA SONIA 1158!SONIA SONIA SONIA: take care. Closed files we need for geoopt???? 1159!SONIA SONIA SONIA 1160 1161 IF ( I .GT. 9 ) CALL CC3_FILCL() 1162C 1163 100 CONTINUE 1164C 1165C------------------------------------- 1166C Make summary of CC calculations. 1167C------------------------------------- 1168C 1169 CALL CC_RESUME(WORK,LWORK) 1170C 1171 CCS = LCCSAV(1) 1172 MP2 = LCCSAV(2) 1173 CCP2 = LCCSAV(3) 1174 CC2 = LCCSAV(5) 1175 CCD = LCCSAV(7) 1176 CCSD = LCCSAV(9) 1177 CCPT = LCCSAV(10) 1178 CCRT = LCCSAV(11) 1179 CCR1A = LCCSAV(12) 1180 CCR1B = LCCSAV(13) 1181 CCP3 = LCCSAV(14) 1182 CCR3 = LCCSAV(15) 1183 CC3 = LCCSAV(16) 1184 CC1A = LCCSAV(17) 1185 CC1B = LCCSAV(18) 1186! FRAN/SONIA 1187 RCCD = LCCSAV(19) 1188 DRCCD = LCCSAV(20) 1189 RTCCD = LCCSAV(21) 1190! AMT 1191 DCPT2 = LCCSAV(22) 1192C 1193C Close files. 1194C 1195C 1196 CALL GPCLOSE(LURES,'KEEP') 1197 CALL GPCLOSE(LUPRPC,'KEEP') 1198 CALL GPCLOSE(LUIAJB,'KEEP') 1199 IF (LUAORC(0).GT.0) CALL GPCLOSE(LUAORC(0),'DELETE') 1200 IF (LUINTR.GT.0) CALL GPCLOSE(LUINTR,'DELETE') 1201 IF (LUINTA.GT.0) THEN 1202 IF (MOLGRD.AND.(.NOT.DIRECT)) THEN 1203 CALL GPCLOSE(LUINTA,'DELETE') 1204 ELSE 1205 CALL GPCLOSE(LUINTA,'KEEP') 1206 END IF 1207 END IF 1208C 1209 DO I = 1, 8 1210 IF (LUAOIN(I) .GT. 0) THEN 1211 IF ((.NOT.KEEPAOIN).OR.MOLGRD) THEN 1212 CALL WCLOSE2(LUAOIN(I),NAME(I),'DELETE') 1213 ELSE 1214 ! needed for gradient calculation... 1215 CALL WCLOSE2(LUAOIN(I),NAME(I),'KEEP') 1216 END IF 1217 END IF 1218 END DO 1219 1220 CALL WCLOSEALL() 1221C 1222C------------------------------------------------ 1223C Restore information on number of roots etc. 1224C if there is an additional call to cc_drv 1225C------------------------------------------------ 1226C 1227 NEXCI = NEXCI2 1228 DO IMULT = 1, 3, 2 1229 DO ISYM =1, 8 1230 NCCEXCI(ISYM,IMULT) = NCCEXCI2(ISYM,IMULT) 1231 IF (IMULT.EQ.1) ISYOFE(ISYM) = ISYOFE2(ISYM) 1232 IF (IMULT.EQ.3) ITROFE(ISYM) = ITROFE2(ISYM) 1233 ENDDO 1234 ENDDO 1235 DO IEX=1,MAXEXCI 1236 ISYEXC(IEX) = ISYEXC2(IEX) 1237 IMULTE(IEX) = IMULTE2(IEX) 1238 ENDDO 1239C 1240 CALL GETTIM(TEND,WEND) 1241 TUSED = TEND - TSTART 1242 WUSED = WEND - WSTART 1243 IF (IPRSTAT .GT. 0) THEN 1244 WRITE (LUSTAT,'(/A,2F12.3)') 1245 & ' CPU and wall time for CC :',TUSED,WUSED 1246 CALL FLSHFO(LUSTAT) 1247 END IF 1248 WRITE (LUPRI,'(/A,2F12.3)') 1249 & ' CPU and wall time for CC :',TUSED,WUSED 1250 CALL TSTAMP(' ',LUPRI) 1251 CALL FLSHFO(LUPRI) 1252C 1253 CALL QEXIT('CC_DRV') 1254C 1255 END 1256C /* Deck cc_resume */ 1257 SUBROUTINE CC_RESUME(WORK,LWORK) 1258C 1259C Ove Christiansen 14-11-1995 1260C 1261C Purpose: Resume all results obtained in this call to coupled 1262C cluster routines: ground state energies, excitation energies etc. 1263C 1264C 1265C 1266#include "implicit.h" 1267#include "priunit.h" 1268#include "mxcent.h" 1269#include "iratdef.h" 1270#include "dummy.h" 1271 DIMENSION WORK(LWORK) 1272 CHARACTER*80 STR 1273 CHARACTER STHELP*10 1274 CHARACTER LABEL*10, LABEL1*8 1275 1276#include "codata.h" 1277#include "ccorb.h" 1278#include "maxorb.h" 1279#include "ccsdinp.h" 1280#include "ccsections.h" 1281#include "cclr.h" 1282#include "prpc.h" 1283#include "ccpack.h" 1284#include "inftap.h" 1285#include "ccinftap.h" 1286#include "numder.h" 1287C 1288 LOGICAL FIRSTCALL 1289 SAVE FIRSTCALL 1290 DATA FIRSTCALL /.TRUE./ 1291C 1292 CALL QENTER('CC_RESUME') 1293C 1294 REWIND(LURES) 1295C 1296 WRITE (LUPRI,'(A,/)') ' ' 1297 WRITE (LUPRI,'(1x,A)') 1298 *'*********************************************************'// 1299 *'**********************' 1300 WRITE (LUPRI,'(1x,A)') 1301 *'*********************************************************'// 1302 *'**********************' 1303 WRITE (LUPRI,'(1x,A)') 1304 *'* '// 1305 *' *' 1306 WRITE (LUPRI,'(1x,A)') 1307 *'* '// 1308 *' *' 1309 WRITE (LUPRI,'(1x,A)') 1310 *'* SUMMARY OF COUPLED CLUSTER CALCULATION '// 1311 *' *' 1312 WRITE (LUPRI,'(1x,A)') 1313 *'* '// 1314 *' *' 1315 WRITE (LUPRI,'(1x,A)') 1316 *'* '// 1317 *' *' 1318 WRITE (LUPRI,'(1x,A)') 1319 *'*********************************************************'// 1320 *'**********************' 1321 WRITE (LUPRI,'(1x,A,/)') 1322 *'*********************************************************'// 1323 *'**********************' 1324C 1325 100 READ(LURES,'(A80)',END=200) STR 1326C 1327 WRITE(LUPRI,'(A80)') STR 1328 GOTO 100 1329C 1330 200 CONTINUE 1331C 1332 IF ( CCEXCI ) THEN 1333 WRITE(LUPRI,'(/A11,F15.5,A6)') ' 1 a.u. = ',XTEV,' eV. ' 1334 WRITE(LUPRI,'(A11,F15.5,A6)') ' 1 a.u. = ',XTKAYS,' cm-1.' 1335 ENDIF 1336C 1337 IF (FIRSTCALL) THEN 1338 FIRSTCALL = .FALSE. 1339 CALL CC_WPRPC 1340C 1341C----------------------------------------------- 1342C Save the list ("old prpc") on file. 1343C Count the number of properties in 1344C each symmetry class 1345C----------------------------------------------- 1346C 1347 DO ISYM=1,8 1348 NPRPINSYM(ISYM) = 0 1349 ENDDO 1350 LUPRPCO = -1 1351 CALL GPOPEN(LUPRPCO,'CC_PRPC_O','UNKNOWN',' ','FORMATTED', 1352 * IDUMMY,.FALSE.) 1353 REWIND(LUPRPCO) 1354 REWIND(LUPRPC) 1355 1356 NPRPCO = NPRPC 1357 DO I = 1, NPRPC 1358 READ(LUPRPC, 1359 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1360 * IPRPC,ISYMIN,NORD,LABEL,PROP, 1361 * LABX,LABY,LABZ,LABU,FRQY,FRQZ,FRQU,ISYEX,ISPEX,IEX 1362 WRITE(LUPRPCO, 1363 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1364 * IPRPC,ISYMIN,NORD,LABEL,PROP, 1365 * LABX,LABY,LABZ,LABU,FRQY,FRQZ,FRQU,ISYEX,ISPEX,IEX 1366 NPRPINSYM(ISYMIN) = NPRPINSYM(ISYMIN) + 1 1367 END DO 1368 CALL GPCLOSE(LUPRPCO,'KEEP') 1369C 1370 ELSE 1371C 1372C Reorder the list of properties to match the firstcall order 1373C 1374 CALL CC_PRPCREORDER 1375C 1376C *** For numerical differentiation of cc-properties, write the *** 1377C *** property to file. *** 1378C 1379 IF (NMRDRP.GT.0) THEN 1380 KPRPC = 1 1381 KEND = KPRPC + NPRPC 1382 LEND = LWORK - KEND 1383 IF (LEND.LE.0)CALL QUIT('Too little workspace in cc_resume') 1384 REWIND(LUPRPC) 1385 1386 DO I = 1, NPRPC 1387 READ(LUPRPC, 1388 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1389 * IPRPC,ISYMIN,NORD,LABEL,WORK(KPRPC+I-1), 1390 * LABX,LABY,LABZ,LABU,FRQY,FRQZ,FRQU,ISYEX,ISPEX,IEX 1391 ENDDO 1392 CALL WRAVFL(WORK(KPRPC),NPRPC,1,1,'CC PROPER',IPRINT) 1393 END IF 1394 ENDIF 1395C 1396C Writeout the list of properties to summary section of output file. 1397C 1398 1399 CALL FLSHFO(LUPRI) 1400 IF ((NPRPC.GT.0).AND.(IPRINT.GT.100)) THEN 1401 CALL PRPRPC(LUPRPC,1,DUMMY,NPRPC) 1402 ENDIF 1403C 1404 IF ( LPACKINT ) THEN 1405 WRITE (LUPRI,'(//10X,A,G9.2)') 1406 & 'Threshold used for packing of integrals:',THRPCKINT 1407 WRITE (LUPRI,'(10X,A,F8.2,A)') 1408 & 'Time used for packing and unpacking: ', 1409 & PCKTIME, ' seconds' 1410 END IF 1411 WRITE (LUPRI,'(A,/)') ' ' 1412 WRITE (LUPRI,'(1x,A)') 1413 *'*********************************************************'// 1414 *'**********************' 1415 WRITE (LUPRI,'(1x,A)') 1416 *'*********************************************************'// 1417 *'**********************' 1418 WRITE (LUPRI,'(1x,A)') 1419 *'* '// 1420 *' *' 1421 WRITE (LUPRI,'(1x,A)') 1422 *'* '// 1423 *' *' 1424 WRITE (LUPRI,'(1x,A)') 1425 *'* END OF COUPLED CLUSTER CALCULATION '// 1426 *' *' 1427 WRITE (LUPRI,'(1x,A)') 1428 *'* '// 1429 *' *' 1430 WRITE (LUPRI,'(1x,A)') 1431 *'* '// 1432 *' *' 1433 WRITE (LUPRI,'(1x,A)') 1434 *'*********************************************************'// 1435 *'**********************' 1436 WRITE (LUPRI,'(1x,A,/)') 1437 *'*********************************************************'// 1438 *'**********************' 1439C 1440 LABEL1 = 'ALL_DONE' 1441 STHELP = 'THE_END ' 1442 ZERO=0.0D0 1443 CALL CC_PRPC(ZERO,STHELP,666, 1444 * LABEL1,LABEL1,LABEL1,LABEL1, 1445 * ZERO,ZERO,ZERO,1,0,0,0) 1446 1447 CALL QEXIT('CC_RESUME') 1448C 1449 END 1450C /* Deck cc_false */ 1451 SUBROUTINE CC_FALSE() 1452C 1453C Ove Christiansen 14-11-1995 1454C 1455C Purpose: set all model flags to false. 1456C 1457#include "implicit.h" 1458#include "ccsdinp.h" 1459#include "cclr.h" 1460C 1461 CALL QENTER('CC_FALSE') 1462C 1463 CIS = .FALSE. 1464 CCS = .FALSE. 1465 MP2 = .FALSE. 1466 CCP2 = .FALSE. 1467 CC2 = .FALSE. 1468 CCD = .FALSE. 1469 CCSD = .FALSE. 1470 CCPT = .FALSE. 1471 CCP3 = .FALSE. 1472 CCRT = .FALSE. 1473 CCR3 = .FALSE. 1474 CCR1A = .FALSE. 1475 CCR1B = .FALSE. 1476 CC3 = .FALSE. 1477 CC1A = .FALSE. 1478 CC1B = .FALSE. 1479 CCSDT = .FALSE. 1480C 1481 CCT = .FALSE. 1482 STCCS = .FALSE. 1483CSonia 1484 RCCD = .false. 1485 DRCCD = .false. 1486 RTCCD = .false. 1487 SOSEX = .false. 1488C 1489 CALL QEXIT('CC_FALSE') 1490C 1491 END 1492C /* Deck cc_false */ 1493 SUBROUTINE CC_FDGD 1494C 1495C Ove Christiansen 13-2-1997 1496C 1497C Purpose: Set energy for finite difference gradient calculation. 1498C 1499#include "implicit.h" 1500#include "priunit.h" 1501#include "ccsdinp.h" 1502#include "ccgr.h" 1503#include "ccfdgeo.h" 1504#include "infopt.h" 1505C 1506 CALL QENTER('CC_FDGD') 1507C 1508 IF (IPRINT.GT.2) THEN 1509 WRITE (LUPRI,'(/,1X,A)') 1510 * '*********************************************************'// 1511 * '**********' 1512 END IF 1513 IF ((IXSTSY.GT.0).AND.(IXSTAT.GT.0)) THEN 1514 ECORR = ECCXST 1515 IF (IPRINT.GT.2) THEN 1516 WRITE(LUPRI,'(A,I3,A,I3)') ' Excited state nr.',IXSTAT, 1517 * ' of symmetry ',IXSTSY 1518 WRITE(LUPRI,'(A,F20.10)') ' Ground state energy: ', ECCGRS 1519 WRITE(LUPRI,'(A,F20.10)') ' Excitation energy: ', OMECCX 1520 WRITE(LUPRI,'(A,F20.10)') ' Excited state energy:',ECCXST 1521 ENDIF 1522 ELSE 1523 ECORR = ECCGRS 1524 IF (IPRINT.GT.2) THEN 1525 WRITE(LUPRI,'(/,1X,A)') 'Ground state energy is target ' 1526 ENDIF 1527 ENDIF 1528C 1529 IF (IPRINT.GT.2) THEN 1530 WRITE(LUPRI,'(/,A,F30.20)') ' Ecorr = ',Ecorr 1531 WRITE(LUPRI,'(/,1X,A)') 1532 * '*********************************************************'// 1533 * '**********' 1534 END IF 1535C 1536 CALL QEXIT('CC_FDGD') 1537C 1538 END 1539*=====================================================================* 1540 SUBROUTINE CCRDABAINF(MOLGRDCC) 1541*---------------------------------------------------------------------* 1542* 1543* Purpose: read flags from /ABAINF/ common block used in CC 1544* 1545*=====================================================================* 1546#include "implicit.h" 1547#include "mxcent.h" 1548#include "abainf.h" 1549 1550 LOGICAL MOLGRDCC 1551 1552 CALL QENTER('CCRDABAINF') 1553 1554 MOLGRDCC = MOLGRD 1555 1556 CALL QEXIT('CCRDABAINF') 1557 1558 RETURN 1559 END 1560 1561*=====================================================================* 1562* END OF SUBROUTINE CCRDABAINF * 1563*=====================================================================* 1564*=====================================================================* 1565 SUBROUTINE CCTAPINI 1566*---------------------------------------------------------------------* 1567* 1568* Purpose: Initialize coupled cluster unit numbers and file names 1569* 1570*=====================================================================* 1571#include "implicit.h" 1572#include "ccinftap.h" 1573C 1574 CALL QENTER('CCTAPINI') 1575C 1576 LUBFDN = -1 1577 LURES = -1 1578 LUCSOL = -1 1579C 1580 CALL QEXIT('CCTAPINI') 1581C 1582 RETURN 1583 END 1584 1585*=====================================================================* 1586 SUBROUTINE CC_WPRPC 1587C 1588C Ove Christiansen Nov. 1999 1589C 1590C Purpose: Write out operators. 1591C 1592#include "implicit.h" 1593#include "ccsdinp.h" 1594#include "ccroper.h" 1595#include "dummy.h" 1596C 1597 CALL QENTER('CC_WPRPC ') 1598C 1599 LUPRP = -1 1600 CALL GPOPEN(LUPRP,'PRPCOP','UNKNOWN',' ','FORMATTED', 1601 * IDUMMY,.FALSE.) 1602C 1603 DO IROPER = 1,NRSOLBL 1604 WRITE(LUPRP,'(I3,A8,I3)') IROPER,LBLOPR(IROPER),ISYOPR(IROPER) 1605 END DO 1606C 1607 CALL GPCLOSE(LUPRP,'KEEP') 1608C 1609 CALL QEXIT('CC_WPRPC ') 1610C 1611 END 1612 1613*=====================================================================* 1614 SUBROUTINE CC_PRPCREORDER 1615C 1616C Ove Christiansen Nov. 1999 1617C 1618C Purpose: Reorder PRPC list to fit old PRPC2 lists. 1619C Property list for numerical differentiation. 1620C 1621#include "implicit.h" 1622#include "prpc.h" 1623#include "priunit.h" 1624#include "ccinftap.h" 1625 CHARACTER LABELI*10, LABXI*8, LABYI*8, LABZI*8, LABUI*8 1626 CHARACTER LABELJ*10, LABXJ*8, LABYJ*8, LABZJ*8, LABUJ*8 1627C 1628c DIMENSION FRQ12(MXPRPC),FRQ22(MXPRPC),FRQ32(MXPRPC) 1629c DIMENSION PRPC2(MXPRPC) 1630c CHARACTER*8 LAB02(MXPRPC), LAB12(MXPRPC), LAB22(MXPRPC), 1631c * LAB32(MXPRPC) 1632c CHARACTER*10 CPRPC2(MXPRPC) 1633c DIMENSION IOPRPC2(MXPRPC), ISAVSY2(MXPRPC) 1634 PARAMETER (TOLPRPFD=1.0D-10,TOLE=1.0D-02) 1635C 1636 CALL QENTER('CC_WPRPCREORDER') 1637C 1638C TEST TEST 1639C IF (NPRPC.NE.NPRPCO) CALL QUIT('Mistmatch in WPRPC') 1640C TEST TEST 1641C 1642C----------------------------------------------- 1643C Save the current list on ("prpc_tmp") on file. 1644C----------------------------------------------- 1645C 1646 LUPRPCTMP = -1 1647 CALL GPOPEN(LUPRPCTMP,'CC_PRPC_TMP','UNKNOWN',' ','FORMATTED', 1648 * IDUMMY,.FALSE.) 1649 REWIND(LUPRPCTMP) 1650 REWIND(LUPRPC) 1651 DO I = 1, NPRPC 1652 READ(LUPRPC, 1653 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1654 * IPRPCI,ISYMINI,NORDI,LABELI,PROPI, 1655 * LABXI,LABYI,LABZI,LABUI,FRQYI,FRQZI,FRQUI, 1656 * ISYEXI,ISPEXI,IEXI 1657 WRITE(LUPRPCTMP, 1658 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1659 * IPRPCI,ISYMINI,NORDI,LABELI,PROPI, 1660 * LABXI,LABYI,LABZI,LABUI,FRQYI,FRQZI,FRQUI, 1661 * ISYEXI,ISPEXI,IEXI 1662 ENDDO 1663C 1664C----------------------------------------------- 1665C Open the saved list ("prpc_o") on file. 1666C----------------------------------------------- 1667C 1668 LUPRPCO = -1 1669 CALL GPOPEN(LUPRPCO,'CC_PRPC_O','UNKNOWN',' ','FORMATTED', 1670 * IDUMMY,.FALSE.) 1671C 1672C Rewind luprpc to make ready for the ordered set of data. 1673C 1674 REWIND(LUPRPC) 1675C 1676C Rewind luprpco to prepare the reference set of data. 1677C 1678 REWIND(LUPRPCO) 1679 DO 300 I = 1, NPRPC 1680 READ(LUPRPCO, 1681 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1682 * IPRPCI,ISYMINI,NORDI,LABELI,PROPI, 1683 * LABXI,LABYI,LABZI,LABUI,FRQYI,FRQZI,FRQUI, 1684 * ISYEXI,ISPEXI,IEXI 1685 1686C 1687C Rewind luprpc_tmp to prepare the data that is reordered 1688C 1689 REWIND(LUPRPCTMP) 1690 DO 400 J = 1, NPRPC 1691 READ(LUPRPCTMP, 1692 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1693 * IPRPCJ,ISYMINJ,NORDJ,LABELJ,PROPJ, 1694 * LABXJ,LABYJ,LABZJ,LABUJ,FRQYJ,FRQZJ,FRQUJ, 1695 * ISYEXJ,ISPEXJ,IEXJ 1696 IF (LABELJ.EQ.LABELI) THEN 1697 IF (NORDI.EQ.NORDJ) THEN 1698 IF ((LABXJ.EQ.LABXI).AND. 1699 * (LABYJ.EQ.LABYI).AND. 1700 * (LABZJ.EQ.LABZI).AND. 1701 * (LABUJ.EQ.LABUI)) THEN 1702C Take care with excitation energies !!!! 1703 IF ( 1704 * ((NORDI.GT.0).AND. 1705 * ((ABS(FRQYJ-FRQYI).LT.TOLPRPFD).AND. 1706 * (ABS(FRQZJ-FRQZI).LT.TOLPRPFD).AND. 1707 * (ABS(FRQUJ-FRQUI).LT.TOLPRPFD))).OR. 1708 * ((NORDI.LE.0).AND. 1709 * ((ABS(FRQYJ-FRQYI).LT.TOLE).AND. 1710 * (ABS(FRQZJ-FRQZI).LT.TOLE).AND. 1711 * (ABS(FRQUJ-FRQUI).LT.TOLE)))) THEN 1712 WRITE(LUPRPC, 1713 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1714 * IPRPCI,ISYMINJ,NORDJ,LABELJ,PROPJ, 1715 * LABXJ,LABYJ,LABZJ,LABUJ,FRQYJ,FRQZJ,FRQUJ, 1716 * ISYEXJ,ISPEXJ,IEXJ 1717 GOTO 300 1718 ENDIF 1719 ENDIF 1720 ENDIF 1721 ENDIF 1722 400 CONTINUE 1723 300 CONTINUE 1724 1725 CALL GPCLOSE(LUPRPCTMP,'KEEP') 1726 CALL GPCLOSE(LUPRPCO,'KEEP') 1727 CALL QEXIT('CC_WPRPCREORDER') 1728C 1729 END 1730C 1731C 1732 SUBROUTINE PRPRPC(LU,IOPT,PROPIN,NUOFPR) 1733C 1734C Outputs the prpc to lupri the list of properties 1735C that is stored on the file opened with unit number LU 1736C NB: File should be open on entry 1737C IF IOPT = 1 the properties on file is output. 1738C IF IOPT = 2 the properties in vector propin is output with 1739C labels etc. comming from the prpc list. 1740C 1741#include "implicit.h" 1742#include "priunit.h" 1743#include "prpc.h" 1744 CHARACTER STHELP*10 1745 CHARACTER LABEL*10, LABX*8, LABY*8, LABZ*8, LABU*8 1746 DIMENSION PROPIN(*) 1747 1748 WRITE (LUPRI,'(A,/)') ' ' 1749 WRITE (LUPRI,'(A)') 1750 *' +---------------------------------------------------------+ ' 1751 WRITE (LUPRI,'(A)') 1752 *' | List of selected calculated molecular properties (a.u.) | ' 1753 WRITE (LUPRI,'(A,/)') 1754 *' +---------------------------------------------------------+ ' 1755 STHELP = 'Dummy ' 1756 WRITE (LUPRI,*) ' Properties are read from unit number ', LU 1757 REWIND(LU) 1758 DO I = 1, NUOFPR 1759 READ(LU, 1760 * '(I5,I3,I4,1X,A10,E23.16,4(1X,A8),3E23.16,3I4)') 1761 * IPRPC,ISYMIN,NORD,LABEL,PROP, 1762 * LABX,LABY,LABZ,LABU,FRQY,FRQZ,FRQU,ISYEX,ISPEX,IEX 1763 IF (LABEL .EQ. 'THE_END ') THEN 1764 WRITE(LUPRI,'(A)') 1765 WRITE(LUPRI,'(1X,A)') LABEL 1766 return 1767 END IF 1768 IF (I.NE.IPRPC) WRITE(LUPRI,*) ' PROBLEM OM PROPERTY LIST?' 1769 IF (IOPT.EQ.2) PROP = PROPIN(I) 1770 1771 IF (STHELP.NE.LABEL) WRITE(LUPRI,'(/,1X,A10,A)') LABEL, 1772 * ' properties :' 1773 STHELP = LABEL 1774 IF (NORD.EQ.0) THEN 1775 IF (LABX.EQ."Excita ") THEN 1776 IF (IOPT.NE.2) THEN 1777 WRITE(LUPRI,'(1X,A,I3,I2,2(A,I1),A,I3,3(A,F18.10))') 1778 * '#',I,ISYMIN, 1779 * ' Sy: ',ISYEX,' Sp: ',ISPEX,' Exc:',IEX, 1780 * ' E_tot:',FRQU,' w:',PROP,' E_grs:',FRQY 1781 ELSE 1782 WRITE(LUPRI,'(1X,A,I3,I2,2(A,I1),A,I3,1(A,F18.10))') 1783 * '#',I,ISYMIN, 1784 * ' Sy: ',ISYEX,' Sp: ',ISPEX,' Exc:',IEX, 1785 * ' E_tot:',PROP 1786 ENDIF 1787 ELSE IF (LABX.EQ."Exctot ") THEN 1788 WRITE(LUPRI,'(1X,A,I3,I2,A,A8,A,F18.10)') 1789 * '#',I,ISYMIN,' ',LABX,' Excited state Energy = ',PROP 1790 ELSE 1791 WRITE(LUPRI,'(1X,A,I3,I2,A,A8,A,F18.10)') 1792 * '#',I,ISYMIN,' ',LABX,' Ground state Energy = ',PROP 1793 ENDIF 1794 ENDIF 1795 IF (NORD.EQ.1) WRITE(LUPRI,'(1X,A,I3,I2,A,A8,A3,F20.12)') 1796 * '#',I,ISYMIN,' <',LABX,'> = ',PROP 1797 IF (NORD.EQ.-11) 1798 * WRITE(LUPRI,'(1X,A,I3,I2,A,A8,A,I1,I2,I3,F9.6,A,F20.12)') 1799 * '#',I,ISYMIN,' <',LABX,'>(',ISYEX,ISPEX,IEX,FRQY,') = ',PROP 1800 IF (NORD.EQ.2) WRITE(LUPRI, 1801 * '(1X,A,I3,I2,A,A8,A1,A8,A3,F9.6,A,F20.12)') 1802 * '#',I,ISYMIN,' <<',LABX,',',LABY, 1803 * '>>(',FRQY,') =',PROP 1804 IF (NORD.EQ.3) WRITE(LUPRI, 1805 * '(1X,A,I3,I2,A,A8,A1,A8,A1,A8,A3,F9.6,A1,F9.6,A,F20.12)') 1806 * '#',I,ISYMIN,' <<',LABX,',',LABY,',',LABZ, 1807 * '>>(',FRQY,',',FRQZ,') =',PROP 1808 IF (NORD.EQ.4) WRITE(LUPRI, 1809 * '(1X,A,I3,I2,A,A8,A1,A8,A1,A8,A1,A8,A3,F9.6,A1,F9.6,A1,F9.6, 1810 * A,F20.12)') 1811 * '#',I,ISYMIN,' <<',LABX,',',LABY,',',LABZ, 1812 * ',',LABU, 1813 * '>>(',FRQY,',',FRQZ,',',FRQU,') =',PROP 1814 IF (NORD.EQ.-1) WRITE(LUPRI, 1815 * '(1X,A,I3,I2,A,A8,A3,F9.6,A,F15.8,A,F12.8,A)') 1816 * '#',I,ISYMIN,' |<O|',LABX,'|i(',FRQY,')>| =',PROP 1817 IF (NORD.EQ.-2) WRITE(LUPRI, 1818 * '(1X,A,I3,I2,A,F9.6,A,A8,A3,F9.6,A,F15.8,A,F12.8,A)') 1819 * '#',I,ISYMIN,' |<i(',FRQY,')|',LABX,'|f(',FRQZ, 1820 * ')>| =',PROP 1821 IF (NORD.EQ.-20) WRITE(LUPRI, 1822 * '(1X,A,I3,I2,A,A8,A1,A8,A3,F9.6,A,F15.8,A,F12.8,A)') 1823 * '#',I,ISYMIN,' RES{<<',LABX,',',LABY, 1824 * '>>(',FRQY,')} =', 1825 * PROP,' ( ',SQRT(ABS(PROP)),')' 1826 ENDDO 1827 END 1828