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 cimain */ 20 SUBROUTINE CIMAIN(WORK,LWORK,ICONV) 21C 22C 28-Jun-1985 hjaaj 23C 24#ifdef MOD_SRDFT 25 use lucita_mcscf_srdftci_cfg 26#endif 27#include "implicit.h" 28#include "priunit.h" 29 DIMENSION WORK(LWORK) 30#include "dummy.h" 31C 32 PARAMETER (MXCIDFT = 20) 33C 34C Used from common blocks: 35C INFINP: ICI0,NROOCI,MXCIMA,THRCI,LSYM 36C INFORB: NCMOT,NASHT 37C INFDIM: LCINDX 38C dftcom: DFTADD 39C 40#include "maxorb.h" 41#include "gnrinf.h" 42#include "infinp.h" 43#include "inforb.h" 44#include "infdim.h" 45#include "dftcom.h" 46C 47 LOGICAL ADDSRI_SAVE 48C 49 CALL QENTER('CIMAIN') 50 ICIST = ICI0 51 NCROOT = NROOCI 52 THRCIX = THRCI 53C 54C 55 KFREE = 1 56 LFREE = LWORK 57 CALL MEMGET('REAL',KCINDX,LCINDX,WORK,KFREE,LFREE) 58 CALL MEMGET('REAL',KCMO ,NCMOT ,WORK,KFREE,LFREE) 59 CALL MEMGET('REAL',KECI ,NCROOT,WORK,KFREE,LFREE) 60 CALL MEMGET('INTE',KICROO,NCROOT,WORK,KFREE,LFREE) 61C 62 ADDSRI_SAVE = ADDSRI 63 IF (DOCISRDFT) THEN 64#ifndef MOD_SRDFT 65 call quit('CI-srdft not implemented in this version') 66#else 67C ... do not add SR integrals to Fock matrices for CI-DFT 68C (will be true on entry if CI-srDFT preceded by normal RHF) 69 ADDSRI = .FALSE. 70C ... Dont risk adding DFT contributions in later call of FCKMAT 71C The SR-DFT terms needed in CI-DFT is handled by SRFMAT. 72C Maybe this should be added to the generel code and not 73C just the CIDFT code if one wants to run a CI with KS orbitals ! 74C (if the CISRDFT was specified together with a presceding 75C HFSRDFT call, we assume the HFSRDFT is done by now!) 76 DFTADD = .FALSE. 77 DOHFSRDFT = .FALSE. 78#endif 79 ENDIF 80C 81 IF (NASHT .GT. 1 .AND. .NOT. HSROHF) 82 * CALL GETCIX(WORK(KCINDX),LSYM,LSYM,WORK(KFREE),LFREE,0) 83 CALL READMO(WORK(KCMO),9) 84C 85 if(do_sc_ensemble_dft)then 86#ifndef MOD_SRDFT 87 call quit('srdft not implemented in this version') 88#else 89 CALL MEMGET('REAL',KECID,MXCIDFT,WORK,KFREE,LFREE) 90 IMCITR = 0 91 DO I = 1,MXCIDFT 92 IMCITR = IMCITR + 1 93 CALL CICTL(ICIST,NCROOT,MXCIMA,THRCIX,WORK(KCMO), 94 * WORK(KCINDX),WORK(KECI),ICONV,WORK(KICROO), 95 * WORK(KFREE),LFREE) 96 ICIST = 5 97C ... restart CI in next macro ensemble-DFT iteration ... 98 ensemble_energy = WORK(KECI) 99 WORK(KECID + I - 1) = ensemble_energy 100 IF(IMCITR.NE.1) EDIFF = ensemble_energy-WORK(KECID+I-2) 101 IF (IMCITR .GE. MXCIDFT) THEN 102 WRITE(LUPRI,'(///A//A)') 103 & ' ----- ensemble-DFT optimization -- Summary -----', 104 & ' ** Maximum number of Macro ensemble-DFT iterations reached **' 105 GO TO 3 106 ELSE IF ((ABS(EDIFF).LT.THRCIDFT).AND.(IMCITR.NE.1)) THEN 107 WRITE(LUPRI,'(///A//A,I3,A)') 108 & ' ----- ensemble-DFT optimization -- Summary -----', 109 & ' ** Convergence reached in',IMCITR, 110 & ' ensemble-DFT macro iterations.' 111 GO TO 3 112 ENDIF 113 ENDDO 1143 CONTINUE 115 WRITE(LUPRI,'(3(/2X,A))') 116 & '---------------------------------------------------------', 117 & ' Itr. Energy Delta(E)', 118 & '---------------------------------------------------------' 119 WRITE(LUPRI,'(9X,I2,3X,F20.12,A)') 120 & 1,WORK(KECID),' ------' 121 DO I = 2,IMCITR 122 DE = WORK(KECID+I-2) - WORK(KECID+I-1) 123 WRITE(LUPRI,'(9X,I2,2(3X,F20.12))') 124 & I,WORK(KECID + I - 1),DE 125 ENDDO 126 WRITE(LUPRI,'(2X,A/)') 127 & '---------------------------------------------------------' 128 WRITE(LUPRI,'(/2X,A,F20.12/)') 129 & '@ Final variational ensemble-DFT energy: ', 130 & WORK(KECID + IMCITR -1) 131#endif 132 else 133C ================= 134 IF (.NOT.DOCISRDFT .OR. ISTACI .LE. 0) THEN 135C ================= 136 CALL CICTL(ICIST,NCROOT,MXCIMA,THRCIX,WORK(KCMO),WORK(KCINDX), 137 * WORK(KECI),ICONV,WORK(KICROO),WORK(KFREE),LFREE) 138C ================= 139 ELSE 140C ================= 141C 142C regular CIDFT Macro-iteration control 143C 144#ifndef MOD_SRDFT 145 call quit('srdft not implemented in this version') 146#else 147 CALL MEMGET('REAL',KECID,NCROOT,WORK,KFREE,LFREE) 148 IMCITR = 0 149 DO I = 1,MXCIDFT 150 IMCITR = IMCITR + 1 151 CALL CICTL(ICIST,NCROOT,MXCIMA,THRCIX,WORK(KCMO), 152 * WORK(KCINDX),WORK(KECI),ICONV,WORK(KICROO), 153 * WORK(KFREE),LFREE) 154 ICIST = 5 155C ... restart CI in next macro CIDFT iteration ... 156 ECIDFT = WORK(KECI + (ISTACI-1) ) 157 WORK(KECID + I - 1) = ECIDFT 158 IF(IMCITR.NE.1) EDIFF = ECIDFT - WORK(KECID + I - 2) 159 IF (IMCITR .GE. MXCIDFT) THEN 160 WRITE(LUPRI,'(///A//A)') 161 & ' ----- CI-srDFT optimization -- Summary -----', 162 & ' ** Maximum number of Macro CI-srDFT iterations reached **' 163 GO TO 5 164 ELSE IF ((ABS(EDIFF).LT.THRCIDFT).AND.(IMCITR.NE.1)) THEN 165 WRITE(LUPRI,'(///A//A,I3,A)') 166 & ' ----- CI-DFT optimization -- Summary -----', 167 & ' ** Convergence reached in',IMCITR, 168 & ' CI-DFT macro iterations.' 169 GO TO 5 170 ENDIF 171 ENDDO 1725 CONTINUE 173 WRITE(LUPRI,'(3(/2X,A))') 174 & '---------------------------------------------------------', 175 & ' Itr. Energy Delta(E)', 176 & '---------------------------------------------------------' 177 WRITE(LUPRI,'(9X,I2,3X,F20.12,A)') 178 & 1,WORK(KECID),' ------' 179 DO I = 2,IMCITR 180 DE = WORK(KECID+I-2) - WORK(KECID+I-1) 181 WRITE(LUPRI,'(9X,I2,2(3X,F20.12))') 182 & I,WORK(KECID + I - 1),DE 183 ENDDO 184 WRITE(LUPRI,'(2X,A/)') 185 & '---------------------------------------------------------' 186 WRITE(LUPRI,'(/2X,A,F20.12/)') 187 & '@ Final CI-srDFT energy: ',WORK(KECID + IMCITR -1) 188#endif 189C 190C ================= 191 ENDIF 192C ================= 193 end if ! ensemble-DFT 194 CALL MEMREL('CIMAIN',WORK,1,1,KFREE,LFREE) 195 ADDSRI = ADDSRI_SAVE 196C 197 CALL QEXIT('CIMAIN') 198 RETURN 199 END 200C /* Deck cictl */ 201 SUBROUTINE CICTL(ICICTL,NCROOT,MAXITC,THRCIX,CMO,INDXCI, 202 * ECI,ICONV,ICROOT,WRK,LWRK) 203C 204C 20-Jun-1985 hjaaj; l.r. 18-Sep-1990 hjaaj 205C 206C Purpose: 207C Perform CI optimization. 208C If MAXITC .le. 0 then return after CIST 209C (e.g. when called from OPTST to find lowest diagonal elements). 210C 211C Input: 212C ICICTL control of methods used 213C NCROOT number of CI states to converge (from below) 214C MAXITC maximum number of iterations 215C THRCIX threshold for convergence 216C CMO(*) m.o. coefficients 217C 218C Output: 219C ECI(*) final CI energies. 220C ICONV number of CI states converged 221C (NCROOT-ICONV states are not converged) 222C ICROOT(*) index of CI states not converged 223C 224! module dependencies 225 use lucita_mcscf_ci_interface_procedures 226 use mcscf_or_gasci_2_define_cfg, only : 227 & define_lucita_cfg_dynamic 228#ifdef MOD_SRDFT 229 use lucita_mcscf_srdftci_cfg 230#endif 231 232#include "implicit.h" 233 DIMENSION CMO(*), INDXCI(*), ECI(*), ICROOT(*), WRK(*) 234C 235#include "iratdef.h" 236#include "thrldp.h" 237 PARAMETER ( D0 = 0.0D0 , DP5 = 0.5D0, D1 = 1.0D0 ) 238#include "dummy.h" 239 LOGICAL CINMEM 240C 241C Used from common blocks: 242C pgroup : REP 243C CICB01 : NCIRED,NJCR,JCROOT(*),THRCIT 244C INFINP : POTNUC,FLAG(*),LSYM,ISTACI 245C INFVAR : NCONF 246C INFORB : NNASHX 247C cbespn.h : ispin1,ispin2 248C INFDIM : LCINDX,LACIMX,LBCIMX,LPHPMX,MAXPHP 249C INFTAP : LUIT1,LUIT3,LUIT5 250C 251#include "maxorb.h" 252#include "priunit.h" 253#include "pgroup.h" 254#include "cicb01.h" 255#include "infinp.h" 256#include "infvar.h" 257#include "inforb.h" 258#include "cbespn.h" 259#include "infdim.h" 260#include "inftap.h" 261#include "infpri.h" 262C 263C Local variables 264C 265 logical :: diskh2, natorb, moisno, calc_2p_dens 266 character (len=6) :: keywrd 267 integer :: xctype1, xctype2, spindens_lvl_save 268C 269 real*8 :: EMYDFTAUX, EENSEMB, ESRDV_ens, ESRDFT(3) 270 character(len=12) :: myciruntype 271 EMYDFTAUX = 0.0D0 ! initialization 272 EENSEMB = 0.0D0 273 ESRDV_ens = 0.0D0 274C 275C ***************************************************************** 276C Analyze input control specification 277C ***************************************************************** 278C 279C 280 CALL QENTER('CICTL ') 281 TIMTOT = SECOND() 282 WRITE (LUW4,'(//A///)') 283 * ' ----- Output from SIRIUS CI module (CICTL) -----' 284 IF (LUPRI .NE. LUW4) WRITE (LUPRI,'(//A///)') 285 * ' ----- Output from SIRIUS CI module (CICTL) -----' 286 IF (NASHT .LE. 1) THEN 287 WRITE(LUPRI,'(//A,I8/)') 288 * ' *** ERROR *** CICTL CALLED, BUT NASHT =',NASHT 289 CALL QTRACE(LUPRI) 290 CALL QUIT('*** ERROR (CICTL) NASHT .LE. 1') 291 END IF 292 IF (NCROOT .LE. 0) THEN 293 NWARN = NWARN + 1 294 WRITE (LUPRI,'(//A,I4,A/)') 295 * ' *** WARNING,',NCROOT,' CI roots requested, nothing done.' 296 ICONV = 0 297 CALL QTRACE(LUPRI) 298 GO TO 9999 299 END IF 300 IF (NCROOT .GT. NCONF) THEN 301 WRITE (LUPRI,'(//A,I5/A,I5)') 302 * ' *** number of CI roots to calculate reduced from',NCROOT, 303 * ' to the number of configurations =',NCONF 304 NCROOT = NCONF 305 END IF 306 IF (NCROOT .GT. MXCIRT) THEN 307 WRITE (LUPRI,'(//A,I5,A/A,I5,A)') 308 * ' *** ERROR in sirci,',NCROOT,' CI roots requested', 309 * ' which exceeds maximum of',MXCIRT,', increase MXCIRT.' 310 CALL QUIT('SIRCI, number of CI roots exceeds predefined '// 311 & 'maximum') 312 END IF 313C 314C 315C ***************************************************************** 316C Core allocation 317C ***************************************************************** 318C 319C 320 KFCAC = 1 321 IF (DOCISRDFT) THEN 322C allocate space for Short Range ACtive effective 1-el. integrals 323 KSRAC = KFCAC + NNASHX 324 KH2AC = KSRAC + NNASHX 325 ELSE 326 KSRAC = KFCAC 327 KH2AC = KFCAC + NNASHX 328 END IF 329C ... first attempt to keep H2AC in core : 330 DISKH2 = .FALSE. 331 IF (FLAG(69)) DISKH2 = .TRUE. 332 10 CONTINUE 333C 334 IF (DISKH2) THEN 335 KCDIAG = KH2AC + NNASHX 336C ... allocate space for one distribution (* * / u v) of H2AC 337 ELSE 338 KCDIAG = KH2AC + NNASHX*NNASHX 339C ... allocate space for whole H2AC 340 END IF 341 CALL PHPINI(LPHPMX,NCONF,0,MAXPHP) 342C CALL PHPINI(LPHPT,NCVAR,NOVAR,MAXPHP) 343 KW1 = KCDIAG + LPHPMX 344 LW1 = LWRK - KW1 345C 346C 347 KREDCI = KW1 348 KEVEC = KREDCI + (MAXRC*MAXRC+MAXRC) /2 349 KW2 = KEVEC + MAXRC * MAXRC 350 LW2 = LWRK - KW2 351C 352C 353C work space needed for CILIN is .lt. LA + NCSIM*LB 354C 355 LA = LACIMX 356 LB = 2*NCONF + LBCIMX 357C 358 MAXSIM = MIN(NCROOT, (LW2 - LA) / LB ) 359 IF (MAXSIM .LE. 0) THEN 360 IF (.NOT. DISKH2) THEN 361C ... try put H2AC on disk : 362 DISKH2 = .TRUE. 363 GO TO 10 364 END IF 365 LNEED = KW2 + LA + LB 366 CALL ERRWRK('CICTL.CILIN',-LNEED,LWRK ) 367 END IF 368 IF (IPRSTAT .GT. 1) WRITE (LUSTAT,'(/A,I8/A,4I12/)') 369 & ' CICTL : MAXSIM =',MAXSIM, 370 & ' KW1, KW2, LWRK =',KW1,KW2,LWRK 371C 372 11 KCVECS = KW2 373 KSVECS = KCVECS + MAXSIM*NCONF 374 LSVECS = LWRK - KSVECS 375 KW3 = KSVECS + MAXSIM*NCONF 376 LW3 = LWRK - KW3 377 MW3 = MAX(KW3 + 2 * MAXRC, KSVECS + MAXRC*MAXRC) 378 IF (MW3 .GT. LWRK ) THEN 379 MAXSIM = MAXSIM - 1 - (MW3 - LWRK ) / (2*NCONF) 380 WRITE (LUSTAT,'(/A,I8/)') 381 & ' CICTL: MAXSIM REDUCED TO',MAXSIM 382 IF (MAXSIM .GT. 0) GO TO 11 383 CALL ERRWRK('CICTL.CIRED',MW3,LWRK ) 384 END IF 385C 386 NCONF4 = MAX(4,NCONF) 387C 388C 389C ***************************************************************** 390C Initialization 391C Set up initial guess for CI vectors 392C ***************************************************************** 393C 394C 395!#define LUCI_DEBUG 396#ifdef LUCI_DEBUG 397 thrcix = 1.0d-11 398#undef LUCI_DEBUG 399#endif 400 NCIRED = 0 401 THRCIT = THRCIX 402 NJCR = NCROOT 403 JCONV = NCROOT 404 TIMITR = SECOND() 405 IF (.NOT. FLAG(14)) THEN 406 CALL TRACTL(0,CMO,WRK,LWRK) 407C CALL TRACTL(ITRLVL,CMO,WRK,LWRK) 408 FLAG(14) = .TRUE. 409 END IF 410 TIMITR = SECOND() - TIMITR 411C 412C *** Obtain FCAC and H2AC integral matrices 413C 414 CALL CIINTS(DISKH2,CMO, 415 * EMY_CI,WRK(KFCAC),WRK(KH2AC),WRK(KW1),LW1) 416C CALL CIINTS(DISKH2,CMO,EMY,FCAC,H2AC,WRK,LFREE) 417C 418! set keyword for saving the final vector(s) in SIRSAV (called from CISAVE) 419 keywrd = 'CISAVE' 420! if(DOCISRDFT)then 421! thrcit = 5.0d-4 422! maxitc = 50 423! end if 424 425! thrcit = 1.0d-11 426! write(lupri,*) ' incoming thrcit (updated) ==> ',thrcit 427 428 IF(CI_PROGRAM .eq. 'LUCITA ')THEN 429! 430! note for stefan/hjaaj: kcdiag is not used for lucita; we might want to NOT allocate it then above!!! 431! 432 xctype1 = -1 433 xctype2 = -1 434 435 calc_2p_dens = .false. 436 437! restart CI if in iteration > 1 438 if(icictl == 5) srdft_restart_ci = 1 439 440 if(docisrdft)then 441 srdft_ci_with_lucita = .true. 442! calc_2p_dens = spindens_lvl > 1 443 spindens_lvl_save = spindens_lvl 444 spindens_lvl = -1 445 end if 446 447 448 call define_lucita_cfg_dynamic(lsym, 449 & lsym, 450 & istaci, 451 & ispin1, 452 & ispin2, 453 & ncroot, ! # of CI roots 454 & ispin, 455 & maxitc, ! max davidson CI iterations 456 & spindens_lvl, ! 1/2-particle density calcs according to spindens_lvl 457 & iprcix, ! print level 458 & 1.0d-10, ! screening threshold 459 & thrcit, ! default davidson threshold 460 & thrcit, ! default davidson threshold 461 & .false., ! do not calculate nat. orb. occ. num. 462! & .true., ! do not calculate nat. orb. occ. num. 463 & P6FLAG(32).and.FLAG(67), ! wave function analysis 464! & .true., ! wave function analysis 465 & .false., ! no 2-el operators active 466 & .true., ! integrals provided by the MCSCF environment 467 & calc_2p_dens, ! calculate 2-particle density matrix 468 & .false., ! calculate transition density matrix 469 & xctype1, ! vector exchange type1 470 & xctype2, ! vector exchange type2 471 & .false., ! vector exchange active in parallel runs in mc2lu interface (both mc-->lu and lu-->mc) 472 & .true.) ! vector exchange between lucita/mc using io2io, e.g. vectors are written to luit3 473 474 IF(DOCISRDFT)THEN 475 476 477 if(.not.allocated(srdft_srac_lucita))then 478 allocate(srdft_srac_lucita(nnashx)) 479 end if 480 481 srdft_srac_lucita = 0.0d0 482 483 if(.not.allocated(srdft_cmo_lucita))then 484 allocate(srdft_cmo_lucita(ncmot)) 485 end if 486 487 call dcopy(ncmot,cmo,1,srdft_cmo_lucita,1) 488 489 myciruntype = 'srdft ci ' 490 491 else 492 myciruntype = 'initial ci ' 493 end if 494 call mcscf_lucita_interface(wrk(kcvecs),wrk(ksvecs), 495 & wrk(kfcac), 496 & wrk(kh2ac),myciruntype, 497 & wrk(kw3),lw3,iprcix) 498 499 jconv = jconv_c ! jconv_c is set in the interface procedure 500 kvec_generic = kcvecs ! for a unique interface to CISAVE (--> final CI vectors are stored in wrk(kcvecs) NOT ACTIVE AT PRESENT) 501 502 ELSE IF(CI_PROGRAM .eq. 'SIRIUS-CI')THEN 503 504 kvec_generic = KEVEC ! for a unique interface to CISAVE 505 506 CALL CIOPT(EMY_CI,JCONV,ICICTL,NCROOT,MAXITC,CMO,INDXCI, 507 & EJCSR,EJVSR,EDSR,ESRDFT,EMYDFT,EMYDFTAUX, 508 & NCLIN,ITCI, 509 & TIMLIN,TIMCIT,NCONF4,LSVECS,MAXSIM, 510 & residual_croot,energy_root, 511 & KFCAC,KH2AC,KCDIAG,KW1,KCIDIA,KSRAC, 512 & KCVECS,KSVECS,KW3,KW2,KREDCI,KEVEC, 513 & LW1,LW2,LW3,WRK,LWRK) 514C write(lupri,*) 'just after calling ciopt ...' 515C write(lupri,*) 'EMYDFTAUX = ', EMYDFTAUX 516 517 END IF ! CI_PROGRAM switch 518 519 DO 810 I = 1,NCROOT 520 ECI(I) = EMY_CI + energy_root(i) 521 810 CONTINUE 522 523 IF (MAXITC .LE. 0) THEN 524 GO TO 9999 525 END IF 526 527 IF(DOCISRDFT)THEN 528 529 if(ci_program .eq. 'LUCITA ')then 530 call dcopy(nnashx,srdft_srac_lucita,1,wrk(ksrac),1) 531 EMYDFT = emydft_mc2lu 532 EJCSR = ejcsr_mc2lu 533 EJVSR = ejvsr_mc2lu 534 EDSR = edsr_mc2lu 535 ESRDFT(1:3) = edft_mc2lu(1:3) 536 EMYDFTAUX = emydftaux_mc2lu 537 spindens_lvl = spindens_lvl_save 538 if(spindens_lvl > 0)then 539 !> compute the spin-density and S**2 540 DO IST = 1,NCROOT 541 KDVREF = KW2 542 KPVREF = KDVREF + NNASHX 543 KW2A = KPVREF + NNASHX**2 544 LW2A = LW2 - KW2A 545 CALL GETS2REF(WRK(KPVREF),WRK(KDVREF),IST, 546 & dummy,WRK(KW2A),LW2A) 547 end do 548 end if 549 end if 550 551 WRITE(LUPRI,'(/A,I5)') 'SR 0-el. energy for input .STATE', 552 & ISTACI 553 WRITE(LUPRI,'(A/)') '-------------------------------------' 554 WRITE(LUPRI,805) ' SR core Hartree energy :',EJCSR 555 WRITE(LUPRI,805) '- SR valence Hartree energy :',EJVSR 556 WRITE(LUPRI,805) '+ SR Exchange-correlation :',ESRDFT(1) 557 WRITE(LUPRI,805) '- SR Exchange-correlation pot.:',EDSR 558 WRITE(LUPRI,806) '= Total eff. SR 0-el. energy :',EMYDFT 559C 560 561 562 DO IST = 1,NCROOT 563 KDVREF = KW2 564 KW2A = KDVREF + NNASHX 565 LW2A = LW2 - KW2A 566 CALL GETDVREF(ci_program,WRK(KDVREF),IST, 567 & INDXCI,WRK(KW2A),LW2A) 568C ... DVREF is equal to DV for state abs(ISTI), 569C i.e. 1-el. energy can be calculated 570 ESRDV = DDOT(NNASHX,WRK(KDVREF),1,WRK(KSRAC),1) 571C write(lupri,*) 'before printing final results ...' 572C write(lupri,*) 'ESRDV = ', ESRDV 573 WRITE(LUPRI,'(/A,I5)') 'CI-DFT energy for state no.',IST 574 WRITE(LUPRI,'(A/)') '-------------------------------------' 575 WRITE(LUPRI,805) 'SR eff. 1-el. energy :',ESRDV 576 EJTOT = ESRDV + EDSR + EJVSR + EJCSR 577 WRITE(LUPRI,805) 'SR total Hartree energy :',EJTOT 578 ETDFT = EMYDFT + ESRDV 579 WRITE(LUPRI,805) 'SR eff. total DFT energy :',ETDFT 580 ELRCI = ECI(IST) - ETDFT 581 WRITE(LUPRI,805) 'LR total CI energy :',ELRCI 582 WRITE(LUPRI,806) 'Total CI-DFT energy :',ECI(IST) 583 ECIAUX = ELRCI+EMYDFTAUX+ESRDV-POTNUC 584! 585 if (IST .EQ. 1) then 586 ECIAUX_gs = ECIAUX ! saving the GS auxiliary energy 587 end if 588! write(lupri,*) 'eciaux_gs =', ECIAUX_gs 589! 590 write(lupri,'(/A)') 591 & ' Decomposition of the auxiliary CI-srDFT energy:' 592 WRITE(LUPRI,805) 'ELRCI :',ELRCI 593 WRITE(LUPRI,805) 'EMYDFTAUX :',EMYDFTAUX 594 WRITE(LUPRI,805) 'ESRDV :',ESRDV 595 WRITE(LUPRI,805) 'POTNUC :',POTNUC 596 write(lupri,*) '' 597 WRITE(LUPRI,807) 'Auxiliary CI-srDFT energy for root', 598 & IST, ': ', ECIAUX 599 if (IST .GT. 1) then 600! WRITE(LUPRI,806) 'Auxiliary excitation energy :', 601! & ECIAUX-ECIAUX_gs 602 WRITE(LUPRI,807) 'Auxiliary excitation energy for root', 603 & IST, ': ', ECIAUX-ECIAUX_gs 604 endif 605! compute the ensemble CI-srDFT energy 606 if(do_sc_ensemble_dft) 607 & EENSEMB = EENSEMB + weights(IST) * ECI(IST) 608! compute the ensemble ESRDV energy contribution 609! ESRDV_ens = ESRDV_ens + weights(IST) * ESRDV 610 END DO 611 805 FORMAT( 1X,A,F25.12) 612 806 FORMAT(/1X,A,F25.12) 613 807 FORMAT(/1X,A,I2,A,F20.12) 614 615 if(do_sc_ensemble_dft)then 616 617 WRITE(LUPRI,'(/A/A,F25.12/A)') 618 & ' *****************************************', 619 & ' Non-variational ensemble CI-srDFT energy:', EENSEMB, 620 & ' *****************************************' 621 622 veensemb = eensemb 623 call get_sc_ensemble_energy(veensemb,ESRDFT, 624 & cmo,nnashx,nnorbt, 625 & ncroot,weights,lupri, 626 & wrk(kw2),lw2) 627 end if 628 ENDIF 629 630 J = LUW4 631 815 CONTINUE 632 633 WRITE (J,'(//A,I2,3A/,(A,I5,0P,F25.15,1P,D15.2))') 634 & '@ Final CI energies and residuals in symmetry',LSYM, 635 & ' (irrep ',REP(LSYM-1),')', 636 & ('@',I,ECI(I),residual_croot(I), I = 1,NCROOT) 637 638 IF (J .NE. LUPRI) THEN 639 J = LUPRI 640 GO TO 815 641 END IF 642C 643C ***************************************************************** 644C Save final result 645C ***************************************************************** 646C 647C Save CI vectors and orbitals on LUIT1 648C 649 if(maxitc .gt. 0 .or. ci_program .eq. 'LUCITA ')then 650 CALL CISAVE(KEYWRD,NCROOT,NCIRED,ECI,residual_croot, 651 & WRK(kvec_generic),CMO,INDXCI,WRK(KW2),LW2) 652C CALL CISAVE(NCROOT,NCIRED,ECI,RESID,EVEC,CMO,INDXCI,WRK,LFREE) 653 end if 654 ICONV = 0 655 J = 0 656 DO 820 I = 1,NCROOT 657 ICROOT(I) = 0 658 IF (JCROOT(I) .EQ. 0) THEN 659 ICONV = ICONV + 1 660 ELSE 661 J = J + 1 662 ICROOT(J) = JCROOT(I) 663 END IF 664 820 CONTINUE 665C 666 IF (FLAG(67)) THEN 667 668 WRITE (LUW4,'(//T6,A/T6,A/T6,A)') 669 & '***********************************', 670 & '*** CI natural orbital analysis ***', 671 & '***********************************' 672 673 IPRNO = MAX(IPRI4-2,3) ! default: IPRI4 = IPRUSR + 5 674 KCVEC = KW2 675 KOCCNO= KCVEC + NCONF 676 KUNO = KOCCNO+ NORBT 677 KW3 = KUNO + N2ASHX 678 LW3 = LWRK - KW3 679C 680 REWIND LUIT1 681 CALL MOLLAB('STARTVEC',LUIT1,LUPRI) 682 DO 900 JSTATE = 1,NCROOT 683 CALL READT(LUIT1,NCONF,WRK(KCVEC)) 684C 685 IF (JSTATE .EQ. ISTACI) THEN 686 WRITE (LUW4,2000) JSTATE 687 IF (LUPRI .NE. LUW4) WRITE (LUPRI,2000) JSTATE 688 NATORB = FLAG(68) 689 ELSE 690 WRITE (LUW4,2001) JSTATE 691 IF (LUPRI .NE. LUW4) WRITE (LUPRI,2001) JSTATE 692 NATORB = .FALSE. 693 END IF 694 2000 FORMAT(/' ------------------------------------------------', 695 & /' CI eigenvector no.',I4,' (= the reference state)', 696 & /' ------------------------------------------------') 697 2001 FORMAT(/' ----------------------', 698 & /' CI eigenvector no.',I4, 699 & /' ----------------------') 700C 701 IF (P6FLAG(32) .OR. NATORB) THEN 702C NATORB is true if CINO and reference state 703 WRITE (LUPRI,'(/A,F8.5)') 704 & ' Printout of CI-coefficients abs greater than:',THRPWF 705 CALL PRWF(WRK(KCVEC),INDXCI,LUPRI,THRPWF,WRK(KW3),LW3) 706C CALL PRWF(C,INDXCI,IOUT,THRPWF,WORK,LFREE) 707 END IF 708C 709C TODO MAERKE insert calculation of spin-density matrix here, 710C if ISPIN .ne. 1 711C 712C insert option for calculation of two-electron density matrix, 713C using MAKDM, for deeper analysis of state vectors. 714C 715 IF (NATORB) THEN 716 KCMO = KW3 717 KW4 = KCMO + NCMOT 718 CALL DCOPY(NCMOT,CMO,1,WRK(KCMO),1) 719 ELSE 720 KCMO = KW3 721 KW4 = KCMO 722 END IF 723 KWNO = 1 724 LWNO = LWRK - KW4 + 1 725 CALL GETNO(WRK(KCVEC),LSYM,WRK(KOCCNO),WRK(KUNO),WRK(KCMO), 726 * .TRUE.,NATORB,.FALSE.,MOISNO,INDXCI, 727 * WRK(KW4),KWNO,LWNO,LUW4,IPRNO) 728C CALL GETNO(CVEC,ICSYM,OCCNO,UNO,CMO,ORDER,NATORB,NOAVER, 729C * MOISNO,INDXCI,WRK,KFRSAV,LFRSAV,LUPRI,IPRINT) 730 IF (NATORB) THEN 731 WRITE (LUW4,'(/A/)') 732 * ' The natural orbitals are saved with label "NEWORB ".' 733 CALL NEWORB('CINOSAVE',WRK(KCMO),.TRUE.) 734 END IF 735 900 CONTINUE 736 END IF 737C 738C Generate density and dump to file. 739C 740C CALL MCRHO(CMO,INDXCI,WRK(KW2),LW2,IPRSIR) 741C 742 743 IF (FLAG(4) .AND. FLAG(25)) THEN 744C ... if (cicalc and save for response) then 745 IF (DOCISRDFT) THEN 746 NWARN = NWARN + 1 747 WRITE (LUPRI,'(//1X,A/)') 748 & 'WARNING, CI RESPONSE is not possible after CI-srDFT' 749 ELSE IF (FLAG(67) .AND. FLAG(68)) THEN 750 NWARN = NWARN + 1 751 WRITE (LUPRI,'(//1X,A/)') 752 & 'WARNING, CI RESPONSE is not possible after CI-NO' 753 ELSE 754 IF (.NOT. FLAG(14)) THEN 755C ... orbitals have been changed by natural orbital 756C transformation in CISAVE 757 WRITE (LUPRI,'(//A/A/A)') ' Molecular orbitals'// 758 * ' were modified after CI convergence', 759 * ' integrals are transformed to new orbital basis', 760 * ' before saving information for post-programs' 761 CALL TRACTL(0,CMO,WRK,LWRK) 762C CALL TRACTL(ITRLVL,CMO,WRK,LWRK) 763 FLAG(14) = .TRUE. 764C 765C *** Obtain FCAC and H2AC integral matrices 766C 767 CALL CIINTS(DISKH2,CMO, 768 * EMY_CI,WRK(KFCAC),WRK(KH2AC),WRK(KW1),LW1) 769C CALL CIINTS(DISKH2,CMO,EMY,FCAC,H2AC,WRK,LFREE) 770C 771C *** Calculate the diagonal of the CI matrix. 772C 773 CALL CIDIAG(LSYM,.FALSE.,WRK(KFCAC),WRK(KH2AC),INDXCI, 774 & WRK(KCDIAG),WRK(KW1),LW1) 775C CALL CIDIAG(ICSYM,NOH2,FCAC,H2AC,XNDXCI,DIAGC,WRK,LFREE) 776C 777 EMY_CI = EMY_CI + POTNUC 778 ELSE 779C subtract EMY from diagonal again, as expected in CIRSPS 780 DO 950 I = 0,NCONF-1 781 WRK(KCDIAG+I) = WRK(KCDIAG+I) - EMY_CI 782 950 CONTINUE 783 END IF 784C 785 CALL CIRSPS(ECI(ISTACI),EMY_CI,CMO,WRK(KFCAC),WRK(KH2AC), 786 * WRK(KCDIAG),INDXCI,WRK(KW2),LW2) 787C CALL CIRSPS(ECIREF,EMY,CMO,FCAC,H2AC,HDIAG,INDXCI,WRK,LFREE) 788 END IF 789 END IF 790C 791 IF (IPRI6 .GT. 0) THEN 792 TIMTOT = SECOND() - TIMTOT 793 WRITE (LUPRI,'(2(/A,F10.2),2(/A,I3,A,F10.2))') 794 * ' Total time used in CICTL :',TIMTOT, 795 * ' Integral transformation :',TIMITR, 796 * ' Total for',ITCI,' CI iterations :',TIMCIT, 797 * ' hereof used in',NCLIN,' linear transformations :',TIMLIN 798 END IF 799C 800 9999 CONTINUE 801 802 !> save the ensemble-DFT energy for self-consistent convergence 803 if(do_sc_ensemble_dft)then 804 eci(1) = VEENSEMB 805 end if 806C 807 if(ci_program .eq. 'LUCITA ')then 808 if(allocated(srdft_srac_lucita)) 809 & deallocate(srdft_srac_lucita) 810 if(allocated(srdft_cmo_lucita)) 811 & deallocate(srdft_cmo_lucita) 812 end if 813 814 CALL FLSHFO(LUW4) 815 IF (LUPRI.NE.LUW4) CALL FLSHFO(LUPRI) 816 CALL QEXIT('CICTL ') 817 RETURN 818C 819C 820C 821C ***************************************************************** 822C End of CICTL 823C ***************************************************************** 824C 825 END 826C /* Deck ciints */ 827 SUBROUTINE CIINTS(DISKH2,CMO,EMY,FCAC,H2AC,WRK,LFREE) 828C 829C Written by Hans Joergen Jensen 13-Dec-1984 830C Revision for LUCITA Dec 2010 /HJAaJ (moved CALL CIDIAG outside) 831C 832C Purpose: 833C CIIN*: Calculate integrals needed for CI. EMY, FCAC and H2AC 834C are needed by calling program. 835C 836#include "implicit.h" 837 838 DIMENSION CMO(*), FCAC(*),H2AC(*),WRK(*) 839 LOGICAL DISKH2 840 841#include "priunit.h" 842#include "dummy.h" 843C 844C Used from common blocks: 845C INFORB : NNASHX,NNBAST,NNORBT 846C INFINP : SRHYBR 847C 848#include "maxorb.h" 849#include "inforb.h" 850#include "infinp.h" 851 852 LOGICAL SRHYBR_SAVE 853C 854 CALL QENTER('CIINTS') 855C 856C *** Section 0 *** 857C 858C work space allocation: 859C 860 NF = 1 861 SRHYBR_SAVE = SRHYBR 862 IF (SRHYBR) THEN 863 SRHYBR = .FALSE. 864 WRITE(LUPRI,'(//A)') ' INFO: .SRHYBR ignored in CIINTS' 865! NF = 3 866 END IF 867 KFC = 1 868 KW0 = KFC + NF*NNORBT 869 LW0 = LFREE - KW0 870C 871 KW1 = 1 872 LW1 = LFREE - KW1 873C 874 KEND = MAX(KW0,KW1) 875 IF (KEND.GT.LFREE) CALL ERRWRK('CIINTS',KEND,LFREE) 876C 877C *** Section 1 *** 878C 879C Calculate F-matrix FCAC, 880C and integrals H2AC. 881C 882 CALL FCKMAT(.TRUE.,DUMMY,CMO,EMY,WRK(KFC),DUMMY, 883 & WRK(KW0),LW0) 884C CALL FCKMAT(ONLYFC,DV,CMO,EMCMY,FC,FV,WRK,LFREE) 885 CALL GETAC(WRK(KFC),FCAC) 886 IF (DISKH2) THEN 887 CALL CIIN2B(H2AC,WRK(KW1),LW1) 888 ELSE 889 CALL CIIN2A(H2AC,WRK(KW1),LW1) 890 END IF 891C 892C *** End of subroutine CIINTS 893C 894 SRHYBR = SRHYBR_SAVE 895 CALL QEXIT('CIINTS') 896 RETURN 897 END 898C /* Deck ciin2a */ 899 SUBROUTINE CIIN2A(H2AC,WRK,LWRK) 900C 901C Written by Hans Agren 27-Jun-1987 902C Rewritten 891017 hjaaj using NXTH2M 903C 904C Revisions: 905C Purpose : To obtain the necessary two-electron 906C integrals to perform a CI calculation. 907C 908#include "implicit.h" 909C 910 DIMENSION H2AC(NNASHX,NNASHX),WRK(LWRK) 911C 912C Used from common blocks: 913C INFORB : NNASHX,NNBAST,? 914C INFIND : ISMO,? 915C CBGETDIS : IADH2 916C 917#include "maxash.h" 918#include "maxorb.h" 919#include "priunit.h" 920#include "inforb.h" 921#include "infind.h" 922#include "infpri.h" 923#include "cbgetdis.h" 924C 925 DIMENSION NEEDMU(-4:6) 926 CALL QENTER('CIIN2A') 927C 928C get array of 2-electron integrals with active indices 929C for ci-use. 930C 931C 932 IF (NASHT .LE. 1) THEN 933 CALL QTRACE(LUPRI) 934 CALL QUIT('PROGRAM ERROR, CIIN2A called with NASHT .le. 1') 935 END IF 936C 937C *** Initialize H2AC *** 938C 939 CALL DZERO(H2AC,NNASHX*NNASHX) 940C 941 KFREE = 1 942 LFREE = LWRK 943 CALL MEMGET('REAL',KH2CD,N2ORBX,WRK,KFREE,LFREE) 944C 945C **************************************************************** 946C Loop over active-active distributions 947 NEEDMU(-4:6) = 0 948 NEEDMU(3) = 1 ! only active-active distributions needed (type 3) 949C 950 IDIST = 0 951 100 CALL NXTH2M(IC,ID,WRK(KH2CD),NEEDMU,WRK,KFREE,LFREE,IDIST) 952 IF (IDIST .LT. 0) GO TO 800 953C ... if idist .lt. 0 then no more distributions 954C 955C We know that IC and ID are both active orbitals 956C because we have specified so in NEEDMU. 957C 958#if defined (VAR_DEBUG) 959 IF (IOBTYP(IC).NE.JTACT .OR. IOBTYP(ID).NE.JTACT) THEN 960 CALL QUIT('CIIN2A: NXTH2M did not return '// 961 & 'active-active distribution') 962 END IF 963#endif 964 IF (IC.LT.ID) THEN 965 ISWAP = IC 966 IC = ID 967 ID = ISWAP 968 END IF 969 ICSYM = ISMO(IC) 970 IDSYM = ISMO(ID) 971 ICDSYM = MULD2H(ICSYM,IDSYM) 972 NCW = ICH(IC) 973 NDW = ICH(ID) 974 NCDW = IROW(NCW) + NDW 975C 976C Add active-active elements to H2AC 977C 978 CALL ADH2AC(H2AC(1,NCDW),WRK(KH2CD),ICDSYM) 979C 980C Go get next active-active distributions 981C 982 GO TO 100 983C 984C arrive at 800 when finished with all active-active distributions 985C 986 800 CONTINUE 987C 988C *** IADH2 .lt. 0 means that H2AC is in core. 989C 990 IADH2 = -1 991C 992C 993 IF (P6FLAG(30)) THEN 994 NCDW = 0 995 DO 300 NCW = 1,NASHT 996 DO 300 NDW = 1,NCW 997 NCDW = NCDW + 1 998 WRITE(LUPRI,3100)NCW,NDW 999 CALL OUTPAK(H2AC(1,NCDW),NASHT,-1,LUPRI) 1000 300 CONTINUE 1001 END IF 1002 3100 FORMAT(/' CIIN2A: Two-electron integrals H2AC, distribution', 1003 * 2I4) 1004C end of CIIN2A 1005C 1006 CALL QEXIT('CIIN2A') 1007 RETURN 1008 END 1009C /* Deck ciin2b */ 1010 SUBROUTINE CIIN2B(H2ACCD,WRK,LWRK) 1011C 1012C 27-6 1987 Hans Agren 1013C 1014C Purpose: 1015C 1016C To get array of 2-electron integrals with active indices 1017C and to store them on direct access file (LUDA2=24), one 1018C CD-distribution per record. 1019C 1020C Special SIRIUS ordering of MO integrals assumed (only 1021C one distribution in each record). 1022C 1023C 1024C 1025C Scratch: H2ACCD (active 2-el integrals) 1026C WRK 1027C 1028C **************************************************************** 1029#include "implicit.h" 1030 DIMENSION H2ACCD(NNASHX), WRK(LWRK) 1031#include "iratdef.h" 1032C 1033C 1034C Used from common blocks: 1035C INFORB : NNASHX,... 1036C INFIND : ISMO,ICH 1037C INFTAP : LUH2AC 1038C INFPRI : P6FLAG 1039C CBGETDIS : IADH2 1040C 1041#include "maxash.h" 1042#include "maxorb.h" 1043#include "priunit.h" 1044 LOGICAL OLDDX 1045#include "inforb.h" 1046#include "infind.h" 1047#include "inftap.h" 1048#include "infpri.h" 1049#include "cbgetdis.h" 1050 1051 DIMENSION NEEDMU(-4:6) 1052 1053 CALL QENTER('CIIN2B') 1054 IF (NASHT .LE. 1) THEN 1055 CALL QTRACE(LUPRI) 1056 CALL QUIT('PROGRAM ERROR, CIIN2B called with NASHT .le. 1') 1057 END IF 1058C 1059C Define off-set on LUH2AC for H2AC. 1060C (IADH2 .lt. 0 means that H2AC is in core.) 1061C 1062C MAERKE 900302: lav mekanisme til at samordne allok. af 1063C H2AC, H2XAC, H2DAC. 1064C 1065 IADH2 = 0 1066C 1067C open direct access file 1068C 1069 CALL GPOPEN(LUH2AC,'MOH2AC.DA','NEW','DIRECT',' ',NNASHX*IRAT, 1070 & OLDDX) 1071C 1072C 1073C **************************************************************** 1074C 1075 KFREE = 1 1076 LFREE = LWRK 1077 CALL MEMGET('REAL',KH2CD,N2ORBX,WRK,KFREE,LFREE) 1078C 1079C **************************************************************** 1080C Loop over Mulliken distributions allowed in NEEDMU(*) 1081C 1082 NEEDMU(-4:6) = 0 1083 NEEDMU(3) = 1 ! only active-active distributions needed (type 3) 1084 1085 IDIST = 0 1086 100 CALL NXTH2M(IC,ID,WRK(KH2CD),NEEDMU,WRK,KFREE,LFREE,IDIST) 1087 IF (IDIST .LT. 0) GO TO 800 1088C ... if idist .lt. 0 then no more distributions 1089C 1090 IF (IC.LT.ID) THEN 1091 ISWAP = IC 1092 IC = ID 1093 ID = ISWAP 1094 END IF 1095 ICSYM = ISMO(IC) 1096 IDSYM = ISMO(ID) 1097 ICDSYM = MULD2H(ICSYM,IDSYM) 1098 NCW = ICH(IC) 1099 NDW = ICH(ID) 1100 NCDW = IROW(NCW) + NDW 1101C 1102C Collect active-active elements in H2ACCD 1103C 1104 CALL DZERO(H2ACCD,NNASHX) 1105 CALL ADH2AC(H2ACCD,WRK(KH2CD),ICDSYM) 1106C 1107C Write this CD-distribution on disk: 1108C 1109 CALL WRITDX(LUH2AC,IADH2+NCDW,IRAT*NNASHX,H2ACCD) 1110 IF (P6FLAG(30)) THEN 1111 WRITE (LUPRI,3100) NCW,NDW 1112 CALL OUTPAK(H2ACCD,NASHT,-1,LUPRI) 1113 END IF 1114 3100 FORMAT(/' CIIN2B: Two-electron integrals H2AC, distribution', 1115 * 2I4) 1116C 1117C Go get next active-active Mulliken distribution 1118C 1119 GO TO 100 1120C 1121C arrive at 800 when finished with all needed Mulliken distributions 1122C 1123 800 CONTINUE 1124C 1125C *** end of subroutine CIIN2B 1126C 1127 CALL QEXIT('CIIN2B') 1128 RETURN 1129 END 1130C /* Deck cilin */ 1131 SUBROUTINE CILIN(NCSIM,BCVECS,FCAC,H2AC,INDXCI,WRK,LFREE) 1132C 1133C Written 30-Apr-1985 by Hans Jorgen Aa. Jensen 1134C (based on LINTRN) 1135C 1136C Purpose: 1137C Do direct CI: calculate the simultaneous linear transformation 1138C defined by the CI matrix H of the NCSIM(ultaneous) B CI vectors : 1139C 1140C SCVECS(i,j) = sum(k=1,NCONF) H(i,k) * BCVECS(k,j) 1141C (i = 1,NCONF; j = 1,NCSIM) 1142C 1143C Input: 1144C 1145C Output: 1146C SCVECS; will start at WRK(1) on output; followed by 1147C 1148#include "implicit.h" 1149#include "infvar.h" 1150 DIMENSION BCVECS(NCONF,*), FCAC(*), H2AC(*) 1151 DIMENSION INDXCI(*),WRK(LFREE) 1152C 1153C *** local constants 1154C 1155 PARAMETER (D1 = 1.0D0, D2 = 2.0D0) 1156C 1157C Used from common blocks: 1158C INFVAR : NCONF 1159C CBGETDIS : DISTYP,IADINT,IADH2 1160C 1161#include "maxorb.h" 1162#include "cbgetdis.h" 1163C 1164 CALL QENTER('CILIN ') 1165C 1166C *** Allocate work area for CISIGC 1167C 1168 KSCVEC = 1 1169 KCW = KSCVEC + NCSIM*NCONF 1170 LCW = LFREE - KCW 1171C-- 1172 IF (LCW .LT. 0) CALL ERRWRK('CILIN',-KCW,LFREE) 1173C 1174C *** call CISIGC to calculate CI sigma vector 1175C 1176 DISTYP = 1 1177 IADINT = IADH2 1178 CALL CISIGC(NCSIM,BCVECS,WRK(KSCVEC),NCONF, 1179 * FCAC,H2AC,INDXCI,WRK(KCW),LCW) 1180C CALL CISIGC(NSIM,BCVECS,SCVECS,LSCVEC,FCAC,H2AC,INDXCI,WRK,LFREE) 1181C 1182C *** MAERKE implement solvent contribution. 1183C 1184C *** end of SUBROUTINE CILIN 1185C 1186 CALL QEXIT('CILIN ') 1187 RETURN 1188 END 1189C /* Deck cinex */ 1190 SUBROUTINE CINEX(EVAL,EVEC,EMY,ESHIFT,CDIAG, 1191 * NCSIM,JCONV,RESID,WRK,LFREE) 1192C 1193C Written 21-Jun-1985 by Hans Jorgen Aa. Jensen 1194C (based on NEONEX version 2) 1195C Revisions: 1196C 1-Jul-1987 hjaaj: CTRUNC and SYMTRZ. 1197C 1198C Purpose: 1199C Find the next trial vectors for the CI optimization. 1200C 1201C Input: 1202C EVAL(NCIRED), eigenvalues of reduced CI matrix. 1203C EVEC(NCIRED,NCIRED), corresponding eigenvectors. 1204C ESHIFT, shift of "Davidson" denominator (usually zero) 1205C CDIAG(*), diagonal of CI matrix + any PHP information 1206C NCSIM = NCROOT, number of roots we want to converge to 1207C 1208C 1209C Output: 1210C JCONV number of roots not converged 1211C NCSIM abs value is number of CI trial vectors written to LUIT3 1212C if positive they are returned in WRK(1) 1213C 1214C note: if NCSIM .eq. 0 and JCONV .gt. 0, then not converged 1215C ---- and no new trial vectors could be found with the algorithm 1216C used (e.g. all removed because of linear dependency). 1217C 1218C Scratch: 1219C WRK(LFREE) 1220C 1221#include "implicit.h" 1222 DIMENSION EVAL(*), EVEC(NCIRED,*), CDIAG(*), RESID(*), WRK(*) 1223C 1224C ********** Define parameters, common blocks, and local variables. 1225C 1226C 1227 PARAMETER ( THRAUX = 0.1D0 ) 1228 PARAMETER ( THKC = 0.10D0, THKO = 0.001D0 ) 1229 PARAMETER ( THLIU1 = 0.01D0, THLIU2 = 0.0D0 ) 1230 PARAMETER ( THRSIM = 1.D-6 , DTEST = 0.01D0 ) 1231#include "thrldp.h" 1232 PARAMETER ( D0=0.0D0, D1 = 1.0D0, DHALF=0.5D0 ) 1233C 1234C Common blocks: 1235C NJCR, number of residual vectors to construct 1236C (number of eigenvalues/-vectors of REDL to use) 1237C Used from common blocks: 1238C CICB01 : NJCR,NKCR,JCROOT(*),KCROOT(*),NCIRED,THRCIT 1239C INFINP : FLAG(*),JOLSEN,RESPHP,ISTACI 1240C INFVAR : NCONF 1241C INFTAP : LUIT3,LUIT5 1242C INFPRI : P4FLAG(*) 1243C PRIUNIT: IPRSTAT 1244C 1245#include "maxorb.h" 1246#include "priunit.h" 1247#include "cicb01.h" 1248#include "infinp.h" 1249#include "infvar.h" 1250#include "inftap.h" 1251#include "infpri.h" 1252C 1253C local variables: 1254C 1255 LOGICAL FIRST, ALLRES, WCRESD, CTRUNC, SYMTRZ 1256C 1257C 1258C data statement: 1259C 1260 SAVE FIRST 1261 DATA FIRST /.TRUE./ 1262C 1263C ********** 1264C ********** End of declarations, start execution. 1265C ********** First, define some entities we will need. 1266C ********** 1267C 1268 CALL QENTER('CINEX ') 1269 IF (FIRST) THEN 1270 WRITE (LUW4,'(//A,F15.8/)') 1271 * ' Convergence threshold for CI optimization :',THRCIT 1272 IF (LUPRI .NE. LUW4) WRITE (LUPRI,'(//A,F15.8/)') 1273 * ' Convergence threshold for CI optimization :',THRCIT 1274 IF (.NOT.DOCISRDFT .AND. ISTACI.GT.0 .AND. NROOCI .GT. 1) THEN 1275 WRITE (LUW4,1010) THRAUX 1276 IF (LUPRI.NE.LUW4) WRITE (LUPRI,1010) THRAUX 1277 END IF 1278 TLINDP = 100*NCONF*THRLDP 1279 TLINDP = SQRT(TLINDP) 1280 IF (THRCIT .LT. TLINDP) THEN 1281 WRITE(LUPRI,'(//A/A,2(/A,D10.2))') 1282 * ' CINEX ERROR: Convergence threshold is less than the', 1283 * ' numerical linear dependence limit', 1284 * ' Convergence threshold ',THRGRD, 1285 * ' Linear dependence limit ',TLINDP 1286 CALL QUIT('CINEX:convergence threshold below lin.dep.limit') 1287 END IF 1288 FIRST = .FALSE. 1289 END IF 1290 1010 FORMAT(' Convergence threshold for auxiliary roots:',F15.8/) 1291 NCIT3 = NCIRED 1292C 1293 WCRESD = FLAG(76) 1294C ... use energy weighted residuals for convergence 1295C this gives more uniform accuracy in state vectors 1296 CTRUNC = FLAG(77) 1297C ... "truncate ci trial vectors", i.e. zero small elements. 1298 SYMTRZ = FLAG(78) 1299C ... symmetrize CI vector 1300C 1301C ALLRES = FLAG(59) 1302 ALLRES = .TRUE. 1303C ... if flag(59) check convergence of all roots,, 1304C roots may have switched order if CI matrix has symmetry 1305C 890915-hjaaj: ALLRES true always now. 1306C 1307 IF (.NOT.ALLRES) NCSIM = NJCR 1308C ... else NCSIM = NCROOT (set by calling routine) 1309 DO 1 I = 1,NCSIM 1310 KCROOT(I) = JCROOT(I) 1311 IF (ALLRES) JCROOT(I) = I 1312 1 CONTINUE 1313C 1314C ********** 1315C ********** Find MAXRCV: how many we can treat simultaneously. 1316C ********** 1317C MWRK0 is space needed for arrays independent of MAXRCV. 1318C MWRK1 is space needed per residual vector. 1319C 1320 MWRK0 = NCONF * 2 1321 IF (JOLSEN) THEN 1322 MWRK1 = NCONF * 2 1323C solution vector plus residual vector 1324 ELSE 1325 MWRK1 = NCONF 1326C only residual vector 1327 END IF 1328 MAXRCV = MIN(NCSIM, (LFREE-MWRK0) / MWRK1) 1329 IF (IPRSTAT .GE. 5) THEN 1330 WRITE (LUSTAT,'(/,(A,I10))') 1331 * ' (CINEX) number of residual vectors',NCSIM, 1332 * ' number of residuals/batch ',MAXRCV 1333 WRITE (LUSTAT,'(/A,I4,A/,(10I5))') ' Index of',NJCR, 1334 * ' non-conv. vectors from last iteration:', 1335 * (KCROOT(I),I=1,NJCR) 1336 END IF 1337 IF (MAXRCV.LE.0) CALL ERRWRK('CINEX',(MWRK0+MWRK1),LFREE) 1338C 1339 NJCR = NCSIM 1340C ... set NJCR to actual number of residual vectors 1341C (if ALLRES it is larger than no. non-conv. vectors) 1342C 1343C CHECK IF WE SHOULD INCREASE NCSIM (Liu simultaneous expansion) 1344C BECAUSE EVAL(last root to converge) 1345C IS CLOSE TO EVAL(first root not to converge). 1346C 1347C CRIT1 / THLIU1 : for absolute test 1348C CRIT2 / THLIU2 : maybe for relative test, we haven't decided yet 1349C what the algorithm should look like. 851213/hj 1350C 1351 MJCR = NJCR 1352 LSTJCR = 0 1353 DO 2 I = 1,NJCR 1354 LSTJCR = MAX(LSTJCR,JCROOT(I)) 1355 2 CONTINUE 1356 IJCR = LSTJCR 1357 3 IF (IJCR .LT. NCIRED .AND. IJCR .LT. MXCIRT) THEN 1358 CRIT1 = EVAL(IJCR+1) - EVAL(IJCR) 1359C CRIT2 = ABS(EVAL(IJCR)) / (EVAL(NCIRED) - EVAL(1)) ??? 1360 IF ( CRIT1 .LT. THLIU1 ) THEN 1361C here ought to be a symmetry test of IJCR vs. IJCR + 1 1362 MJCR = MJCR + 1 1363 IJCR = IJCR + 1 1364 JCROOT(MJCR) = IJCR 1365 GO TO 3 1366 END IF 1367 END IF 1368 IF (IPRSTAT .GT. 0 .AND. MJCR .GT. NJCR) THEN 1369 WRITE (LUSTAT,'(//A,2(/A,I4),//A)') 1370 * ' -- CINEX, number of trial vectors increased --', 1371 * ' Number of CI roots :',NJCR, 1372 * ' Number of simultaneous trial vectors :',MJCR, 1373 * ' Lowest eigenvalues of reduced CI matrix:' 1374 I = MIN(NCIRED, MJCR+2) 1375 WRITE (LUSTAT,'(5X,5F12.8)') (EVAL(J), J = 1,I) 1376 END IF 1377C 1378C ********** 1379C ********** Loop over batches 1380C ********** 1381C 1382 ISTCNV = 0 1383 NCSIM = 0 1384 JSIMX = 0 1385 NSIML = MJCR 1386 IBATCH = 0 1387 10 CONTINUE 1388 IBATCH = IBATCH + 1 1389 NSIMX = MIN(MAXRCV,NSIML) 1390C 1391C ********** 1392C ********** Calculate residual vectors. 1393C ********** 1394C 1395C ********** a) core allocation 1396C RVCS for residual vectors 1397C 1398 MSIMX = 0 1399 DO 205 ISIMX = JSIMX+1, JSIMX+NSIMX 1400 IF (JCROOT(ISIMX) .NE. 0) THEN 1401 MSIMX = MSIMX + 1 1402 JCROOT(JSIMX+MSIMX) = JCROOT(ISIMX) 1403 END IF 1404 205 CONTINUE 1405 DO 206 ISIMX = JSIMX+MSIMX+1,JSIMX+NSIMX 1406 JCROOT(ISIMX) = 0 1407 206 CONTINUE 1408C 1409 KRVECS = 1 1410 KSCR = KRVECS + MSIMX*NCONF 1411 IF (JOLSEN) THEN 1412 KXVECS = KSCR 1413 KSCR = KXVECS + MSIMX*NCONF 1414 ELSE 1415 KXVECS = KRVECS 1416 END IF 1417 CALL DZERO (WRK(KRVECS),(KSCR-KRVECS)) 1418C 1419C ********** b) calculate RVCS 1420C 1421C ********** sigma vector part 1422C 1423 REWIND LUIT5 1424 DO 220 ICIRED = 1,NCIRED 1425 JRVECS = KRVECS 1426 CALL READT(LUIT5,NCONF,WRK(KSCR)) 1427 DO 210 ISIMX = JSIMX+1,JSIMX+MSIMX 1428 FAC = EVEC(ICIRED,JCROOT(ISIMX)) 1429 CALL DAXPY(NCONF,FAC,WRK(KSCR),1,WRK(JRVECS),1) 1430 JRVECS = JRVECS + NCONF 1431 210 CONTINUE 1432 220 CONTINUE 1433C 1434C ********** subtract E times eigenvector to finish residual 1435C note that if JOLSEN false, then KXVECS = KRVECS. 1436C 1437 REWIND LUIT3 1438 DO 330 ICIRED = 1,NCIRED 1439 CALL READT(LUIT3,NCONF,WRK(KSCR)) 1440 JXVECS = KXVECS 1441 DO 310 ISIMX = JSIMX+1, JSIMX+MSIMX 1442 JCROTX = JCROOT(ISIMX) 1443 IF (JOLSEN) THEN 1444 FAC = EVEC(ICIRED,JCROTX) 1445 ELSE 1446 FAC = - EVAL(JCROTX) * EVEC(ICIRED,JCROTX) 1447 END IF 1448 CALL DAXPY(NCONF,FAC,WRK(KSCR),1,WRK(JXVECS),1) 1449 JXVECS = JXVECS + NCONF 1450 310 CONTINUE 1451 330 CONTINUE 1452 IF (JOLSEN) THEN 1453 JXVECS = KXVECS 1454 JRVECS = KRVECS 1455 DO 340 ISIMX = JSIMX+1, JSIMX+MSIMX 1456 JCROTX = JCROOT(ISIMX) 1457 FAC = - EVAL(JCROTX) 1458 CALL DAXPY(NCONF,FAC,WRK(JXVECS),1,WRK(JRVECS),1) 1459 JXVECS = JXVECS + NCONF 1460 JRVECS = JRVECS + NCONF 1461 340 CONTINUE 1462 END IF 1463 IF (P6FLAG(43)) THEN 1464 WRITE (LUPRI,9001) (JCROOT(ISIMX),ISIMX=JSIMX+1,JSIMX+MSIMX) 1465 CALL OUTPUT(WRK(KRVECS),1,NCONF,1,MSIMX,NCONF,MSIMX,-1,LUPRI) 1466 END IF 1467 9001 FORMAT(/' (CINEX) residual vectors for roots',(10I4)) 1468C 1469C ********** 1470C ********** Convergence check, contract RVCS 1471C ********** 1472C 1473C 1474 IF (P4FLAG(6)) THEN 1475 IF (WCRESD) THEN 1476 WRITE (LUW4,9012) 1477 ELSE 1478 WRITE (LUW4,9010) 1479 END IF 1480 END IF 1481 JRVECS = KRVECS 1482 DO 1100 ISIMX = JSIMX+1,JSIMX+MSIMX 1483 JCROTX = JCROOT(ISIMX) 1484 IF (JCROTX .EQ. 0) GO TO 1100 1485 ECROOT = EMY + EVAL(JCROTX) 1486 IF (JCROTX .LE. LSTJCR) THEN 1487C this is one of the roots we want to converge 1488 QCNRM = DNRM2(NCONF,WRK(JRVECS),1) 1489C 1490 IF (WCRESD) THEN 1491C 1492C Weighted residual gives more uniform vector 1493C (wave function) accuracy, however, unweighted 1494C is equivalent to CI part of MC gradient in 1495C SIRIUS.SIROPT MCSCF optimization 1496C 1497 QCFAC = MAX( DTEST, ABS(EVAL(JCROTX)) ) 1498 QCNRM = QCNRM / QCFAC 1499 END IF 1500C 1501C save residual for final print 1502C 1503 RESID(JCROTX) = QCNRM 1504C 1505 IF (DOCISRDFT .OR. ISTACI .LE. 0) THEN 1506 THRQC = THRCIT 1507 ELSE 1508 IF (JCROTX .EQ. ISTACI) THEN 1509 THRQC = THRCIT 1510 ELSE 1511 THRQC = THRAUX 1512 END IF 1513 END IF 1514 IF (QCNRM .LE. THRQC) THEN 1515 IF (JCROTX .EQ. LSTJCR) THEN 1516 DO 1095 IJCR = NJCR+1,MJCR 1517 JCROOT(IJCR) = 0 1518 1095 CONTINUE 1519 END IF 1520 JCROOT(ISIMX) = 0 1521 IF (P4FLAG(6)) THEN 1522 WRITE (LUW4,9030) JCROTX,ECROOT,QCNRM 1523 ELSE IF (IPRI4 .GT. 5) THEN 1524 IF (WCRESD) THEN 1525 WRITE (LUW4,9012) 1526 ELSE 1527 WRITE (LUW4,9010) 1528 END IF 1529 WRITE (LUW4,9030) JCROTX,ECROOT,QCNRM 1530 END IF 1531 IF (JCROTX .EQ. ISTACI) THEN 1532 ISTCNV = 1 1533 WRITE (LUW4,9014) 1534 END IF 1535 ELSE 1536 IF (P4FLAG(6)) WRITE(LUW4,9020) JCROTX,ECROOT,QCNRM 1537 END IF 1538 END IF 1539 JRVECS = JRVECS + NCONF 1540 1100 CONTINUE 1541C 1542 9010 FORMAT(/' CI root eigenvalue residual norm'/) 1543 9012 FORMAT(/' CI root eigenvalue weighted residual norm'/) 1544 9014 FORMAT(/' The requested root number is now converged.') 1545 9020 FORMAT(I5,2F20.10) 1546 9030 FORMAT(I5,2F20.10,' converged') 1547C 1548C set all to converged if ISTACI converged 1549C 1550 IF (ISTCNV .NE. 0) THEN 1551 DO 1120 IJCR = 1,NJCR 1552 JCROOT(IJCR) = 0 1553 1120 CONTINUE 1554 END IF 1555C 1556 KCVECS = KRVECS 1557 JRVECS = KRVECS 1558 JCVECS = KCVECS 1559 JXVECS = KXVECS 1560 JYVECS = KXVECS 1561 NCVEC = 0 1562 DO 1200 ISIMX = JSIMX+1,JSIMX+MSIMX 1563 JCROTX = JCROOT(ISIMX) 1564 IF (JCROTX .NE. 0) THEN 1565 NCVEC = NCVEC + 1 1566 KCROOT(NCVEC) = JCROTX 1567 IF (JRVECS.NE.JCVECS) THEN 1568 CALL DCOPY(NCONF,WRK(JRVECS),1,WRK(JCVECS),1) 1569 IF (JOLSEN) THEN 1570 CALL DCOPY(NCONF,WRK(JXVECS),1,WRK(JYVECS),1) 1571 END IF 1572 END IF 1573 JCVECS = JCVECS + NCONF 1574 JYVECS = JYVECS + NCONF 1575 END IF 1576 JRVECS = JRVECS + NCONF 1577 JXVECS = JXVECS + NCONF 1578 1200 CONTINUE 1579 KCPREV = KCVECS + NCVEC*NCONF 1580C ... KCPREV is used in CIORT, it is OK that it overlaps with KXVECS 1581 KWPHP = KXVECS + NCVEC*NCONF 1582 LWPHP = LFREE + 1 - KWPHP 1583C 1584C 1585C ********** 1586C ********** Use Davidson's shift to get CI trial vectors 1587C ********** 1588C 1589C D = ( CIDIAG(I) - EVAL ) + ESHIFT 1590C 1591 IF (ESHIFT.NE.D0 .AND. IPRI4.GT.5) WRITE (LUW4,9120) ESHIFT 1592 9120 FORMAT(/' (CINEX) Davidson level shift shifted by ESHIFT =', 1593 * F20.15) 1594C 1595C 1596 ICRV = KCVECS - 1 1597 DO 2400 ICVEC = 1,NCVEC 1598 KCROTX = KCROOT(ICVEC) 1599#if defined (VAR_NOPHP) 1600 DSHIFT = ESHIFT - EVAL(KCROTX) - EMY 1601 DO 2200 ICONF = 1,NCONF 1602 D = CDIAG(ICONF) + DSHIFT 1603 IF (D.LT.DTEST) THEN 1604 IF (D.GT.(-DTEST)) D = SIGN(DTEST,D) 1605 END IF 1606 WRK(ICRV+ICONF) = WRK(ICRV+ICONF) / D 1607 2200 CONTINUE 1608#else 1609 DSHIFT = EMY + EVAL(KCROTX) - ESHIFT 1610 IF (JOLSEN) THEN 1611 JXVEC = KXVECS + (ICVEC-1)*NCONF 1612 ELSE 1613 JXVEC = 999 999 999 1614 END IF 1615 IPRPHP = MAX(0,IPRI6-5) 1616 CALL NEXCI(JOLSEN,DSHIFT,NCONF,WRK(ICRV+1),WRK(JXVEC), 1617 * WRK(ICRV+1),CDIAG,IPRPHP,WRK(KWPHP),LWPHP) 1618C CALL NEXCI(OLSEN,ENER,NCVAR,D,XVEC, 1619C & RES,DIAG,IPRPHP,WRK,LWRK) 1620#endif 1621 IF (CTRUNC) THEN 1622C ... increase efficiency by zeroing elements in 1623C trial vector abs less than 1.d-3 abs largest element. 1624 FACTRN = 1.D-3 1625 ICMAX = IDAMAX(NCONF,WRK(ICRV+1),1) 1626 CTEST = FACTRN * ABS(WRK(ICRV+ICMAX)) 1627 DO 2300 ICONF = ICRV+1,ICRV+NCONF 1628 IF (ABS(WRK(ICONF)) .LT. CTEST) WRK(ICONF) = D0 1629 2300 CONTINUE 1630 END IF 1631 IF (SYMTRZ) THEN 1632C ... symmetrize CI vector so that numerical 1633C inaccuracy won't break symmetry; 1634C WARNING : if symmetry were broken in residual vectors, 1635C and thus in sigma vectors, then this will 1636C lead to constant residual. 1637C 1638 CNORM = DNRM2(NCONF,WRK(ICRV+1),1) 1639 TLINDP = SQRT(NCONF*THRLDP) 1640 IF (CNORM .GT. TLINDP) THEN 1641 CNORM = D1 / CNORM 1642 CALL DSCAL(NCONF,CNORM,WRK(ICRV+1),1) 1643 CALL VSYMTR(NCONF,WRK(ICRV+1),THRSIM) 1644 END IF 1645 END IF 1646 ICRV = ICRV + NCONF 1647 2400 CONTINUE 1648C 1649C ********** 1650C ********** 1651C ********** 1652C 1653C 1654C ********** 1655C ********** Orthogonalize these NCVEC new b-vectors 1656C ********** against previous b-vectors and each other. 1657C ********** 1658C 1659C CIORT updates NCIT3. 1660C LUIT3 is positioned, ready for the new trial vectors. 1661C 1662C 1663 THRLDV = NCONF * THRLDP 1664 CALL CIORT (NCVEC,NCIT3,LUIT3,KCROOT,WRK(KCVECS),THRLDV, 1665 * WRK(KCPREV)) 1666C CALL CIORT (NCNEW,NCPREV,LUCVEC,KCROOT,CVEC,THRLDP,CPREV) 1667C 1668C ********** 1669C ********** Check if finished all NJCR residuals, 1670C ********** if not go back to 10 for next batch. 1671C ********** 1672C 1673 NCSIM = NCSIM + NCVEC 1674 NSIML = NSIML - NSIMX 1675 JSIMX = JSIMX + NSIMX 1676 IF (NSIML.GT.0) GO TO 10 1677C ^----------------------- 1678C 1679C 1680C ********** 1681C ********** Test if "in memory" 1682C ********** 1683C 1684C NECgh971126 Fopp has problems with this kind of IF 1685C NECgh971126 IF (IBATCH .EQ. 1) THEN 1686C NECgh971126 ELSE 1687 IF (IBATCH .NE. 1) THEN 1688 NCSIM = -NCSIM 1689 END IF 1690C 1691C ********** 1692C ********** Update NJCR and JCROOT by removing converged vectors. 1693C ********** 1694C 1695 NJ = 0 1696 DO 8000 ISIMX = 1,NJCR 1697 IF (JCROOT(ISIMX).GT.0 .AND. JCROOT(ISIMX).LE.LSTJCR) THEN 1698 NJ = NJ + 1 1699 JCROOT(NJ) = JCROOT(ISIMX) 1700 END IF 1701 8000 CONTINUE 1702 NJCR = NJ 1703 JCONV = NJCR 1704C 1705C ********** 1706C ********** End of subroutine CINEX 1707C ********** 1708C 1709 CALL QEXIT('CINEX ') 1710 RETURN 1711 END 1712C /* Deck ciort */ 1713 SUBROUTINE CIORT (NCNEW,NCPREV,LUCVEC,KCROOT, 1714 * CVEC,THRLDP,CPREV) 1715C 1716C Written 23-Jun-1985 by Hans Jorgen Aa. Jensen 1717C 1718C Purpose: 1719C Orthogonalize new CI vectors against all previous CI vectors 1720C and among themselves, and renormalize. 1721C (Orthogonalization is performed twice if round-off is large, 1722C if larger than THRRND). 1723C 1724C Input: 1725C NCNEW, number of new CVEC 1726C NCPREV, number of previous, orthonormal vectors on LUCVEC 1727C LUCVEC, file unit with previous vectors 1728C CVEC(*,*), non-orthogonal configurational vectors 1729C KCROOT(*), index of input CVEC 1730C (if KCROOT(i) = 0 then vector no. i is skipped) 1731C THRLDP, threshold for linear dependence 1732C 1733C Output: 1734C NCNEW, number of linear indep. vectors in CVEC 1735C NCPREV, number of vectors now on LUCVEC 1736C KCROOT(*), index of output CVEC 1737C CVEC, NCNEW orthonormal vectors 1738C 1739C Scratch: 1740C CPREV(*), scratch array for old vectors on LUCVEC 1741C 1742#include "implicit.h" 1743#include "infvar.h" 1744 DIMENSION KCROOT(*), CVEC(NCONF,*), CPREV(*) 1745C 1746 PARAMETER ( THRRND=1.D-2, THRTT=1.D-4, D1=1.0D0 ) 1747C 1748C Used from common blocks: 1749C INFVAR : NCONF 1750C 1751#include "maxorb.h" 1752#include "priunit.h" 1753#include "infpri.h" 1754C 1755C 1756 CALL QENTER('CIORT ') 1757 IF (NCNEW .EQ. 0) GO TO 9999 1758C 1759 NCONF4 = MAX(4,NCONF) 1760 TLINDP = SQRT(THRLDP) 1761C 1762C Normalize input CVEC 1763C 1764 DO 100 I = 1,NCNEW 1765 TT = DNRM2(NCONF,CVEC(1,I),1) 1766 IF (TT.LE.TLINDP) THEN 1767 WRITE (LUPRI,3240) I,KCROOT(I),TT 1768 KCROOT(I) = 0 1769 ELSE 1770 IF (TT.LT.THRTT) THEN 1771 TT = D1 / TT 1772 CALL DSCAL(NCONF,TT,CVEC(1,I),1) 1773 TT = DNRM2(NCONF,CVEC(1,I),1) 1774 END IF 1775 TT = D1 / TT 1776 CALL DSCAL(NCONF,TT,CVEC(1,I),1) 1777 END IF 1778 100 CONTINUE 1779C 1780 ITURN=0 1781 1500 ITURN=ITURN+1 1782 IROUND=0 1783C 1784C Orthogonalize new vectors agains previous vectors 1785C 1786 REWIND LUCVEC 1787 DO 2000 I=1,NCPREV 1788 CALL READT(LUCVEC,NCONF,CPREV) 1789 DO 1910 J = 1,NCNEW 1790 IF (KCROOT(J) .NE. 0) THEN 1791 TT = -DDOT(NCONF,CPREV,1,CVEC(1,J),1) 1792 CALL DAXPY(NCONF,TT,CPREV,1,CVEC(1,J),1) 1793 END IF 1794 1910 CONTINUE 1795 2000 CONTINUE 1796C 1797C Orthogonalize new vectors against each other and normalization 1798C 1799 DO 3200 I = 1,NCNEW 1800 IF (KCROOT(I) .EQ. 0) GO TO 3200 1801C Orthogonalize against previous new vectors 1802C (which have been normalized in the code following 2300) 1803 DO 2300 J = 1,(I-1) 1804 IF (KCROOT(J) .NE. 0) THEN 1805 TT = -DDOT(NCONF,CVEC(1,J),1,CVEC(1,I),1) 1806 CALL DAXPY(NCONF,TT,CVEC(1,J),1,CVEC(1,I),1) 1807 END IF 1808 2300 CONTINUE 1809C Normalize 1810 TT = DNRM2(NCONF,CVEC(1,I),1) 1811 IF (TT .LE. TLINDP) THEN 1812 WRITE (LUPRI,3250) I,KCROOT(I),TT 1813 KCROOT(I) = 0 1814 ELSE 1815 IF (TT .LT. THRRND) IROUND=IROUND+1 1816 IF (TT.LT.THRTT) THEN 1817 TT = D1 / TT 1818 CALL DSCAL(NCONF,TT,CVEC(1,I),1) 1819 TT = DNRM2(NCONF,CVEC(1,I),1) 1820 END IF 1821 TT = D1 / TT 1822 CALL DSCAL(NCONF,TT,CVEC(1,I),1) 1823 END IF 1824 3200 CONTINUE 1825C 1826C 1827 IF (IROUND.GT.0 .AND. ITURN.EQ.1) GO TO 1500 1828 IF (IROUND.GT.0) CALL QUIT('CIORT: round-off errors, see code.') 1829 3240 FORMAT(/' (CIORT) trial vector no.',I4,' for root',I3, 1830 * ' is removed because of linear dependence;', 1831 * /' norm of vector on input',1P,D12.5) 1832 3250 FORMAT(/' (CIORT) trial vector no.',I4,' for root',I3, 1833 * ' is removed because of linear dependence;', 1834 * /' norm of vector after Gram-Schmidt''s orthogonalization', 1835 * 1P,D12.5) 1836C 1837C Save new vector together with old ones on LUCVEC 1838C 1839 J = 0 1840 DO 4300 I = 1,NCNEW 1841 IF (KCROOT(I) .NE. 0) THEN 1842 J = J + 1 1843 IF (J .LT. I) THEN 1844 KCROOT(J) = KCROOT(I) 1845 CALL DCOPY(NCONF,CVEC(1,I),1,CVEC(1,J),1) 1846 END IF 1847 CALL WRITT(LUCVEC,NCONF4,CVEC(1,J)) 1848 END IF 1849 4300 CONTINUE 1850 IF (J .LT. NCNEW) THEN 1851 WRITE (LUPRI,4310) NCNEW,J 1852 NCNEW = J 1853 END IF 1854C 1855C update NCPREV 1856C 1857 NCPREV = NCPREV + NCNEW 1858C 1859 4310 FORMAT(/' (CIORT) NCNEW reduced from',I3,' to',I3) 1860C 1861C *** End of subroutine CIORT 1862C 1863 9999 CALL QEXIT('CIORT ') 1864 RETURN 1865 END 1866C /* Deck circhk */ 1867 SUBROUTINE CIRCHK(WRK,LFREE) 1868C 1869C 1-Aug-1987 Hans Joergen Aa. Jensen, based on REDCHK from SIRNEO. 1870C 1871C Purpose: 1872C Calculate full reduced CI-matrix for symmetry check in 1873C CI calculation; 1874C asymmetry would indicate probable bug in CILIN module. 1875C 1876#include "implicit.h" 1877 DIMENSION WRK(LFREE) 1878 PARAMETER ( D1 = 1.0D0, DIFTST = 1.0D-8 ) 1879C 1880C Used from common blocks: 1881C INFVAR : NCONF 1882C CICB01 : NCIRED 1883C INFTAP : LUIT3, LUIT5 1884C 1885#include "maxorb.h" 1886#include "priunit.h" 1887#include "infvar.h" 1888#include "cicb01.h" 1889#include "inftap.h" 1890C 1891C 1892 REDH(I,J) = WRK((J-1)*NCIRED+I) 1893 REDS(I,J) = WRK((KREDS-1)+(J-1)*NCIRED+I) 1894C 1895 CALL QENTER('CIRCHK') 1896 KREDH = 1 1897 KREDS = KREDH + NCIRED*NCIRED 1898 KSVEC = KREDS + NCIRED*NCIRED 1899 KBVECS = KSVEC + NCONF 1900 LBVECS = LFREE - KBVECS 1901 IF (LBVECS .LT. NCIRED*NCONF) THEN 1902C insufficient space for the algorithm used 1903 WRITE (LUPRI,'(//2A,/A/)') ' %%% Insufficient space', 1904 * ' for symmetry check of reduced CI matrix.', 1905 * ' Nothing done.' 1906 GO TO 9999 1907 END IF 1908C 1909 JBVECI = KBVECS 1910 REWIND LUIT3 1911 DO 100 I = 1,NCIRED 1912 CALL READT(LUIT3,NCONF,WRK(JBVECI)) 1913 JBVECI = JBVECI + NCONF 1914 100 CONTINUE 1915C 1916C REDH(*,I) = B(*,*)t S(*,I) 1917C 1918 I1 = KREDH 1919 REWIND LUIT5 1920 DO 200 I = 1,NCIRED 1921 CALL READT(LUIT5,NCONF,WRK(KSVEC)) 1922 CALL DGEMM('T','N',NCIRED,1,NCONF,1.D0, 1923 & WRK(KBVECS),NCONF, 1924 & WRK(KSVEC),NCONF,0.D0, 1925 & WRK(I1),NCIRED) 1926 I1 = I1 + NCIRED 1927 200 CONTINUE 1928C 1929C Reduced overlap matrix, get b vectors 1930C 1931C 1932C REDS(I,J) = B(*,I)t B(*,J) 1933C 1934 CALL DGEMM('T','N',NCIRED,NCIRED,NCONF,1.D0, 1935 & WRK(KBVECS),NCONF, 1936 & WRK(KBVECS),NCONF,0.D0, 1937 & WRK(KREDS),NCIRED) 1938C 1939C check reduced matrices 1940C 1941 WRITE (LUPRI,'(//A/A/A,1P,D10.2,A)') 1942 * ' %%% Non-symmetric elements in redH, reduced CI matrix,', 1943 * ' and wrong elements in redS, reduced CI overlap', 1944 * ' ( threshold is',DIFTST,' ) :' 1945 NERROR = 0 1946 DO 340 I = 2,NCIRED 1947 DO 320 J = 1,(I-1) 1948 DIFF = ABS( REDH(I,J) - REDH(J,I) ) 1949 IF (DIFF .GT. DIFTST) THEN 1950 NERROR = NERROR + 1 1951 WRITE (LUPRI,'(A,2I5,2F20.14)') 1952 * ' i, j, redh(i,j), redh(j,i) :', 1953 * I,J,REDH(I,J),REDH(J,I) 1954 END IF 1955 IF (ABS(REDS(I,J)) .GT. DIFTST .OR. 1956 * ABS(REDS(J,I)) .GT. DIFTST) THEN 1957 NERROR = NERROR + 1 1958 WRITE (LUPRI,'(A,2I5,2F20.14)') 1959 * ' i, j, reds(i,j), reds(j,i) :', 1960 * I,J,REDS(I,J),REDS(J,I) 1961 END IF 1962 320 CONTINUE 1963 IF (ABS(REDS(I,I) - D1) .GT. DIFTST) THEN 1964 NERROR = NERROR + 1 1965 WRITE (LUPRI,'(A,2I5,2F20.14)') 1966 * ' i, i, reds(i,i) :', 1967 * I,I,REDS(I,I) 1968 END IF 1969 340 CONTINUE 1970 IF (NERROR .EQ. 0) WRITE (LUPRI,'(A)') ' None, redH is symmetric.' 1971C 1972 WRITE (LUPRI,'(//A)') ' Reduced CI matrix (as bT*s) :' 1973 CALL OUTPUT(WRK(KREDH),1,NCIRED,1,NCIRED,NCIRED,NCIRED,-1,LUPRI) 1974 WRITE (LUPRI,'(//A)') ' Reduced overlap (as bT*b) :' 1975 CALL OUTPUT(WRK(KREDS),1,NCIRED,1,NCIRED,NCIRED,NCIRED,-1,LUPRI) 1976C 1977 IF (NERROR .GT. 0) THEN 1978 CALL QTRACE(LUPRI) 1979 CALL QUIT('ERROR detected in CIRCHK.') 1980 END IF 1981C 1982 9999 CALL QEXIT('CIRCHK') 1983 RETURN 1984 END 1985C /* Deck cired */ 1986 SUBROUTINE CIRED (ICTL,NREDCI,REDCI,NCNEW, 1987 * CVECS,SVECS,EVAL,EVEC,WRK,LWRK) 1988C 1989C Written 25-Jun-1985 by Hans Jorgen Aa. Jensen 1990C 1991C Input: 1992C ICTL, flow control 1993C =1, extend the reduced CI-matrix 1994C =2, find the new eigenvalues and -vectors 1995C =3, both 1996C REDCI, the old reduced CI-matrix (dimension: NREDCI-NCNEW) 1997C NREDCI, dimension of new reduced CI-matrix 1998C NCNEW, number of new CI-vectors 1999C CVECS, the NCNEW new CI-vectors (is destroyed) 2000C SVECS, the sigma vectors of the NCNEW new CI-vectors 2001C 2002C (The reduced CI-matrix is the projection of the CI-matrix 2003C on the basis of CI-vectors.) 2004C 2005C Output: 2006C REDCI, the new, extended reduced CI-matrix (dimension: NREDCI) 2007C EVAL, eigenvalues of the new reduced CI matrix 2008C EVEC, eigenvectors of the new reduced CI matrix 2009C 2010C Scratch: 2011C CVECS, is used for CI vectors from unit LUIT3 (section 1) 2012C and for the reduced CI matrix (section 2) 2013C dimension: NCONF and NREDCI*(NREDCI+1)/2, rsp. 2014C WRK, real scratch array 2015C 2016#include "implicit.h" 2017#include "infvar.h" 2018 DIMENSION CVECS(*), SVECS(NCONF,*) 2019 DIMENSION REDCI(*), EVEC(NREDCI,*),EVAL(*), WRK(LWRK) 2020C 2021C Used from common blocks: 2022C INFVAR : NCONF 2023C INFIND : IROW(*) 2024C 2025#include "maxash.h" 2026#include "maxorb.h" 2027#include "infinp.h" 2028#include "infind.h" 2029#include "inftap.h" 2030#include "infpri.h" 2031C 2032C 2033 CALL QENTER('CIRED ') 2034C 2035C Section 1: extend reduced CI matrix with NCNEW new CI-vectors 2036C 2037 IF (MOD(ICTL,2) .NE. 1) GO TO 1000 2038 IF (NCNEW.LT.1) GO TO 1000 2039 IF (NREDCI .GE. LIROW) 2040 & CALL QUIT('CIRED error: NREDCI exceeds LIROW') 2041C IROW(NREDCI+1) used below 2042 IREDCI = NREDCI - NCNEW 2043C 2044C New CI-vectors (are in CVECS) 2045C 2046 K1 = 1 2047 DO 140 K = 1,NCNEW 2048 KL = IROW(IREDCI+K) + IREDCI 2049 DO 120 L=1,K 2050 REDCI(KL+L) = DDOT(NCONF,CVECS(K1),1,SVECS(1,L),1) 2051 120 CONTINUE 2052 K1 = K1 + NCONF 2053 140 CONTINUE 2054C 2055C Old CI-vectors. 2056C Rewind LUIT3 2057C 2058 IF (IREDCI .GT. 0) THEN 2059 REWIND LUIT3 2060 LL = 0 2061 DO 240 J = 1,IREDCI 2062 CALL READT(LUIT3,NCONF,CVECS) 2063 LL = LL + 1 2064 DO 220 K=1,NCNEW 2065 REDCI(IROW(IREDCI+K) + LL ) = 2066 * DDOT(NCONF,CVECS,1,SVECS(1,K),1) 2067 220 CONTINUE 2068 240 CONTINUE 2069 END IF 2070C 2071C ********************************************************* 2072C Section 2: find eigenvalues and -vectors of reduced CI matrix 2073C 2074 1000 CONTINUE 2075 IF (MOD(ICTL,4) .LE. 1) GO TO 2000 2076C 2077C Diagonalize reduced CI matrix. 2078C Copy REDCI to CVECS for diagonalization. 2079C 2080 LMA = IROW(NREDCI+1) 2081 KFREE = 1 2082 LFREE = LWRK 2083 CALL MEMGET('REAL',KREDCI,LMA ,WRK,KFREE,LFREE) 2084 CALL DCOPY(LMA,REDCI,1,WRK(KREDCI),1) 2085 CALL DUNIT(EVEC,NREDCI) 2086 CALL JACO_THR(WRK(KREDCI),EVEC,NREDCI,NREDCI,NREDCI,0.0D0) 2087 II = KREDCI - 1 2088 DO 1100 I = 1,NREDCI 2089 II = II + I 2090 EVAL(I) = WRK(II) 2091 1100 CONTINUE 2092 CALL MEMREL('CIRED.JACO_THR',WRK,1,1,KFREE,LFREE) 2093 CALL ORDER (EVEC,EVAL,NREDCI,NREDCI) 2094C 2095C *** End of subroutine CIRED 2096C 2097 2000 CONTINUE 2098 CALL QEXIT('CIRED ') 2099 RETURN 2100 END 2101C /* Deck cisave */ 2102 SUBROUTINE CISAVE(KEYWRD,NCROOT,NCIRED,ECI,RESID,EVEC,CMO,INDXCI, 2103 * WRK,LFREE) 2104C 2105C 27-Sep-1986 Hans Joergen Aa. Jensen 2106C 2107C revised 870627-hjaaj, write ci energies to LUIT1 for OCE program 2108C 890925-hjaaj, write ci residuals to LUIT1 also 2109C 2110#include "implicit.h" 2111#include "priunit.h" 2112 DIMENSION ECI(*), RESID(*), EVEC(*), CMO(*), INDXCI(*), WRK(*) 2113 character (len=6), intent(in) :: keywrd 2114#include "dummy.h" 2115C 2116C Used from common blocks : 2117C INFINP : NROOTS 2118C INFOPT : NREDL 2119C INFTAP : LUIT1 2120C 2121#include "maxorb.h" 2122#include "infinp.h" 2123#include "infopt.h" 2124#include "inftap.h" 2125#include "infpri.h" 2126C 2127 CHARACTER*8 LABENR(4), LABEOD(4) 2128 DATA LABENR/'********','********','********','ENERGIES'/ 2129 DATA LABEOD/'********','********','********','EODATA '/ 2130C 2131 CALL QENTER('CISAVE') 2132C 2133 NREDL = NCIRED 2134 ISTSAV = ISTATE 2135 NRTSAV = NROOTS 2136 ISTATE = ISTACI 2137 NROOTS = NCROOT 2138 CALL SIRSAV(keywrd,CMO,DUMMY,DUMMY,EVEC,DUMMY,INDXCI,WRK,LFREE) 2139C CALL SIRSAV(KEYWRD,CMO,IBNDX,REDL,EVEC,XKAP,INDXCI,WRK,LFREE) 2140 ISTATE = ISTSAV 2141 NROOTS = NRTSAV 2142C 2143C save ECI on LUIT1 2144C 2145 REWIND LUIT1 2146 CALL MOLLAB('NEWORB ',LUIT1,LUPRI) 2147 READ (LUIT1) 2148C ... skip "new orbitals", and write energies label: 2149C 2150 WRITE (LABENR(3),'(I8)') NCROOT 2151 WRITE (LUIT1) LABENR 2152 NCR4 = MAX(4,NCROOT) 2153 WRITE (LUIT1) (ECI(I),I=1,NCR4) 2154 WRITE (LUIT1) (RESID(I),I=1,NCR4) 2155 CALL GETDAT(LABEOD(2),LABEOD(3)) 2156 WRITE (LUIT1) LABEOD 2157C 2158 CALL QEXIT('CISAVE') 2159 RETURN 2160 END 2161C /* Deck cist */ 2162 SUBROUTINE CIST (ICI1,NCSIM,CMO,HDIAG,EMY,WRK,LFREE) 2163C 2164C Written 20-Jun-1985 by Hans Joergen Aa. Jensen 2165C 2166C Purpose: 2167C Obtain initial CI trial vectors for a CI calculation. 2168C To write molecular orbitals on LUIT1. 2169C The routine is controlled by the ICI1 parameter 2170C (similar to the ICI0 parameter for OPTST). 2171C 2172C ICI1<0 : select -ICI1 as start configuration 2173C 2174C ICI1=1 : Compute startguess of CI-vector from H-diagonal (HDIAG) 2175C 2176C ICI1=2 : Start CI-vector(s) are already on LUIT1 (this is checked). 2177C 2178C ICI1=4 : Same as ICI1=2, for compatibility with MCSCF CISTRT 2179C 2180#include "implicit.h" 2181C 2182 DIMENSION CMO(*),HDIAG(*),WRK(*) 2183C 2184 PARAMETER (D0 = 0.0D0, D1 = 1.0D0) 2185C 2186C Used from common blocks: 2187C INFINP : FLAG(),ISTACI,? 2188C INFVAR : NCONF,? 2189C INFORB : NCMOT,? 2190C INFIND : ? 2191C INFDIM : ? 2192C INFTAP : LUIT1,? 2193C 2194#include "maxash.h" 2195#include "maxorb.h" 2196#include "priunit.h" 2197#include "infinp.h" 2198#include "infvar.h" 2199#include "inforb.h" 2200#include "infind.h" 2201#include "infdim.h" 2202#include "inftap.h" 2203#include "infpri.h" 2204C 2205C *** Externals: 2206C 2207 LOGICAL FNDLAB 2208C 2209C *** Local variables: 2210C 2211 CHARACTER*8 TABLE(4), LAB123(3), CHRDUM 2212C 2213C *** data: 2214C 2215 DATA LAB123(1)/'********'/ 2216 DATA TABLE/'********','OLDORB ','STARTVEC','CIDIAG2 '/ 2217C 2218 CALL QENTER('CIST ') 2219 CALL GETDAT(LAB123(2),LAB123(3)) 2220 LAB123(3) = '( CIST )' 2221 IF (IPRSIR .GT. 1) WRITE(lupri,*) 'CIST called with ICI1 = ',ICI1 2222 KFREE = 1 2223 NC4 = MAX(4,NCONF) 2224 MCSIM = 0 2225 NCSIM0 = NCSIM 2226 JCI1 = ICI1 2227 IF (JCI1 .EQ. 4 .OR. JCI1 .EQ. 5) JCI1 = 2 2228 10 CONTINUE 2229 IF (JCI1 .NE. 2 .AND. MCSIM .EQ. 0) THEN 2230C 2231 CALL NEWIT1 2232 WRITE (LUIT1) LAB123,TABLE(2) 2233 CALL WRITT(LUIT1,NCMOT,CMO) 2234C 2235C If RHF or NCONF.eq.1, simply write C0(1) = D1 to LUIT1 2236C 2237 IF (FLAG(21) .OR. NCONF.EQ.1) THEN 2238 WRK(1) = D1 2239 WRITE (LAB123(2),'(2I4)') 1,1 2240 WRITE (LUIT1) LAB123,TABLE(3) 2241 CALL WRITT(LUIT1,4,WRK) 2242 GO TO 9000 2243 END IF 2244 END IF 2245C 2246C *** Go to appropriate section to obtain start guess for 2247C CI vectors (as specified with the input parameter JCI1) 2248C 2249C *** JCI1 < 0 *** 2250C === start guess = CSF no. -JCI1 2251C 2252 IF (JCI1.LT.0) THEN 2253 IF (-JCI1.GT.NCONF) THEN 2254 WRITE (LUPRI,5010) -JCI1,NCONF 2255 CALL QTRACE(LUPRI) 2256 CALL QUIT('*** ERROR (CIST) illegal JCI1 parameter') 2257 END IF 2258 IF (NCSIM.GT.1) THEN 2259 WRITE (LUPRI,5020) 2260 CALL QTRACE(LUPRI) 2261 CALL QUIT('*** ERROR (CIST) illegal JCI1 parameter') 2262 END IF 2263 CALL DZERO(WRK,NCONF) 2264 WRK(-JCI1) = D1 2265 WRITE (LAB123(2),'(2I4)') 1,-1 2266 WRITE (LUIT1) LAB123,TABLE(3) 2267 CALL WRITT(LUIT1,NC4,WRK) 2268 GO TO 9000 2269 END IF 2270 5010 FORMAT(/' (CIST) ERROR, specified start configuration no.', 2271 * I10,' is non-existent (last CSF is no.',I10,')') 2272 5020 FORMAT(/' (CIST) ERROR, negative JCI1 option not allowed ', 2273 * 'when NCSIM is greater then one.') 2274C 2275 IF (JCI1.LT.1 .OR. JCI1.GT.2) THEN 2276 WRITE (LUPRI,5030) JCI1 2277 CALL QTRACE(LUPRI) 2278 CALL QUIT('*** ERROR (CIST) illegal JCI1 parameter') 2279 END IF 2280 5030 FORMAT(/' CIST FATAL ERROR, JCI1 =',I4,' is not defined') 2281C 2282C 2283C 2284C 2285 NCSIM = NCSIM0 + MCSIM 2286C ... (870804-hjaaj) Until now we used IFRST = MCSIM + 1 2287C and no change in NCSIM, however, it gave bad convergence 2288C in roots MCSIM+1 to NCSIM not to include the lowest 2289C MCSIM diagonal elements explicitly. 2290C 940511-hjaaj MAERKE: Maybe this can be done with .OLSEN ?? 2291 2292 GO TO (100,200) JCI1 2293C 2294C *** JCI1 = 1 *** 2295C *** Use PHPGET and PHPVEC to get lowest eigenvectors of subspace 2296C *** or find smallest diag. elements of H to select start config. 2297C 2298 100 CONTINUE 2299 IF (MAXPHP .GE. 2*NCSIM) THEN 2300 IF (MCSIM .EQ. 0) THEN 2301 WRITE (LAB123(2),'(2I4)') NCSIM,-1 2302 LAB123(3) = 'CIST-PHP' 2303 WRITE (LUIT1) LAB123,TABLE(3) 2304 ELSE 2305C ... some, but not all, requested CI trial vectors available 2306 REWIND LUIT1 2307 CALL MOLLAB(TABLE(3),LUIT1,LUPRI) 2308 DO 161 I = 1,MCSIM 2309 READ (LUIT1) 2310 161 CONTINUE 2311 END IF 2312 DO 110 ICSIM = 1,NCSIM 2313 CALL PHPVEC(1,ICSIM,WRK,NCONF,HDIAG) 2314C CALL PHPVEC(NVEC,I1VEC,CVECS,NCONF,DIAG) 2315 CALL WRITT(LUIT1,NC4,WRK) 2316 110 CONTINUE 2317 NCSIM = MCSIM + NCSIM 2318 ELSE 2319 CALL CIST1(ISTACI,FLAG(58),NCONF,NCSIM,MCSIM,HDIAG,EMY, 2320 & WRK,LFREE) 2321C CALL CIST1(ISTACI,PLUSCB,NCONF,NCSIM,MCSIM,HDIAG,EMY,WRK,LWRK) 2322 END IF 2323 GO TO 9000 2324C 2325C *** JCI1 = 2 *** 2326C === start guess should already be on LUIT1 2327C (but we test to be sure!) 2328C 2329 200 CONTINUE 2330 MCSIM = 0 2331 REWIND LUIT1 2332 IF ( .NOT.FNDLAB(TABLE(2),LUIT1) ) GO TO 290 2333 READ (LUIT1,END=290,ERR=290) 2334 IF ( .NOT.FNDLAB(TABLE(3),LUIT1) ) GO TO 290 2335 DO 210 I = 1,NCSIM 2336 READ (LUIT1,END=290,ERR=290) CHRDUM 2337 IF (CHRDUM .EQ. TABLE(1)) GO TO 290 2338 MCSIM = I 2339 210 CONTINUE 2340 REWIND LUIT1 2341 GO TO 9000 2342C 2343 290 CONTINUE 2344 REWIND LUIT1 2345 IF (MCSIM .EQ. 0) THEN 2346 WRITE (LUW4,'(//A//)') 2347 * ' *** CI control: start vectors not found on LUIT1,'// 2348 * ' lowest diagonal elements used.' 2349 ELSE 2350 WRITE (LUW4,'(//A,I4,A,/A/A//)') 2351 * ' *** CI control: only',MCSIM, 2352 * ' start vectors found on LUIT1,', 2353 * ' In addition to these start vectors, lowest diagonal', 2354 * ' elements used for requested trial vectors.' 2355 END IF 2356 JCI1 = 1 2357 GO TO 10 2358C 2359C *** End of subroutine CIST 2360C 2361C Write NEWORB at end of file (so e.g. SIRORB won't overwrite 2362C STARTVEC). 2363C 2364 9000 CONTINUE 2365 IF (JCI1 .NE. 2) THEN 2366 CALL NEWORB('CIST ',CMO,.FALSE.) 2367 END IF 2368 CALL QEXIT('CIST ') 2369 RETURN 2370 END 2371C /* Deck cist1 */ 2372 SUBROUTINE CIST1(ISTACI,PLUSCB,NCONF,NCSIM,MCSIM, 2373 * HDIAG,EMY,WRK,LWRK) 2374C 2375C 30-Oct-1990 hjaaj 2376C 2377C Block of old CIST, purpose: find start vectors corresponding 2378C to lowest diagonal elements and write these to LUIT1. 2379C 2380#include "implicit.h" 2381 LOGICAL PLUSCB 2382 REAL*8 HDIAG(*), WRK(*) 2383 PARAMETER (THRSIM = 0.05D0, THRDEG = 1.0D-6) 2384C ... THRSIM is threshold for two values are similar 2385C THRDEG is threshold for two values are degenerate 2386C 2387C Used from common blocks: 2388C INFTAP : LUIT1 2389C INFPRI : LUW*, IPRI* 2390C 2391#include "infpri.h" 2392#include "inftap.h" 2393#include "priunit.h" 2394C 2395C *** Local variables: 2396C 2397 INTEGER, ALLOCATABLE :: NLIST(:) 2398 CHARACTER*8 LBLIT1(4) 2399C 2400C *** data: 2401C 2402 DATA LBLIT1/'********',' ',' ','STARTVEC'/ 2403C 2404 CALL QENTER('CIST1 ') 2405C 2406 NC4 = MAX(4,NCONF) 2407C 2408 IF (PLUSCB) THEN 2409C ... take plus combination of degenerate diagonal elements 2410C we thus need more ... 2411 NRS = MIN( NCONF, 2 * NCSIM + 6) 2412 WRITE (LUPRI,'(//A/A/)') 2413 * ' --- SIRCI.CIST1: plus combination of all degenerate', 2414 * ' configurations is used as start vectors.' 2415 ELSE 2416 NRS = MIN(NCSIM + 7,NCONF) 2417C ... we add 5 so we can expand if degeneracy. 2418 END IF 2419 IF (LWRK .LT. NRS+NCONF) CALL ERRWRK('CIST1',(NRS+NCONF),LWRK) 2420 2421 allocate(NLIST(NRS)) 2422 CALL CIFNMN(NRS,NLIST,HDIAG,NCONF,WRK,LWRK) 2423C CALL CIFNMN(NELMNT,IPLACE,VEC,NVEC,WRK,LWRK) 2424 IF (IPRI4 .GT. 2) WRITE (LUW4,1410) 2425 * NRS,(I,NLIST(I),WRK(I)-EMY,WRK(I),I=1,NRS) 2426 IF (LUPRI .NE. LUW4 .AND. IPRSIR .GT. 0) WRITE (LUPRI,1410) 2427 * NRS,(I,NLIST(I),WRK(I)-EMY,WRK(I),I=1,NRS) 2428 1410 FORMAT(/' (CIST1) ',I0,' lowest diagonal elements:', 2429 * //' Element no. Config.no. Active energy', 2430 * ' Total energy', 2431 * //,(I10,' : ',I10,2F18.10)) 2432 2433 IFRST = 1 2434 IF (.NOT.PLUSCB) THEN 2435C ... include the degenerate configurations separately, 2436C check if last element is degenerate with next element(s) 2437 LAST = NCSIM 2438 DO 150 I = NCSIM+1,NRS 2439 IF ((WRK(I)-WRK(NCSIM)).LT.THRSIM) LAST = I 2440 150 CONTINUE 2441 IF (LAST .GT. NCSIM) THEN 2442 WRITE (LUPRI,1500) NCSIM,LAST 2443 NCSIM = LAST 2444 END IF 2445 1500 FORMAT(//' (CIST1) number of start vectors is increased from', 2446 * I3,' to',I3,' because of (near) degeneracy') 2447 END IF 2448C 2449C === Write start guess on LUIT1 2450C 2451 LAST = NCSIM 2452 IF (MCSIM .EQ. 0) THEN 2453 WRITE(LBLIT1(2),'(2I4)') NCSIM,ISTACI 2454 LBLIT1(3) = '(CIST1 )' 2455 WRITE (LUIT1) LBLIT1 2456 ELSE 2457C ... some, but not all, requested CI trial vectors available 2458 REWIND LUIT1 2459 CALL MOLLAB(LBLIT1(4),LUIT1,LUPRI) 2460 DO 161 I = 1,MCSIM 2461 READ (LUIT1) 2462 161 CONTINUE 2463 END IF 2464 2465#ifdef LUCI_DEBUG_test_case 2466 CALL DZERO(WRK(NRS+1),NCONF) 2467 WRK(NRS+2) = 1.0d0 2468 CALL WRITT(LUIT1,NC4,WRK(NRS+1)) 2469 MCSIM = MCSIM + 1 2470 LAST = IFRST 2471#endif 2472 CALL DZERO(WRK(NRS+1),NCONF) 2473 2474 IF (PLUSCB) THEN 2475C ... take plus combination of degenerate diagonal elements 2476 170 CONTINUE 2477C ... 170: repeat until (MCSIM .GE. NCSIM .OR. NRS .GE. LAST) 2478 LAST = IFRST 2479 172 IF ( (WRK(LAST+1)-WRK(IFRST)) .LT. THRDEG ) THEN 2480 LAST = LAST + 1 2481 IF (LAST .LT. NRS) GO TO 172 2482 END IF 2483 FAC = 1 + LAST - IFRST ! yields degeneracy 2484 FAC = 1.0D0 / SQRT(FAC) 2485 IF (LAST .EQ. IFRST) THEN 2486 WRITE(LUPRI,'(/A,6I10)') 2487 & 'Start CI vector is configuration no.',NLIST(IFRST:LAST) 2488 ELSE 2489 WRITE(LUPRI,'(/A,6I10)') 2490 & 'Start CI vector is sum of',NLIST(IFRST:LAST) 2491 END IF 2492 DO 174 I = IFRST,LAST 2493 WRK(NRS + NLIST(I)) = FAC 2494 174 CONTINUE 2495 CALL WRITT(LUIT1,NC4,WRK(NRS+1)) 2496 MCSIM = MCSIM + 1 2497 DO 176 I = IFRST,LAST 2498 WRK(NRS + NLIST(I)) = 0.0D0 2499 176 CONTINUE 2500 IFRST = LAST + 1 2501 IF (MCSIM .LT. NCSIM .AND. LAST .LT. NRS) GOTO 170 2502C ... end repeat 2503 ELSE 2504 WRITE(LUPRI,'(A,I0,A)') 2505 & 'First ',LAST,' elements used as separate start vectors.' 2506 DO 185 I = IFRST,LAST 2507 MCSIM = MCSIM + 1 2508 WRK(NRS + NLIST(I)) = 1.0D0 2509 CALL WRITT(LUIT1,NC4,WRK(NRS+1)) 2510 WRK(NRS + NLIST(I)) = 0.0D0 2511 185 CONTINUE 2512 END IF 2513 IF (NCSIM .GT. MCSIM .AND. LAST .LT. NCONF) THEN 2514C ... NRS too small compared to NCSIM. 2515 WRITE (LUPRI,1850) MCSIM, NCSIM 2516 1850 FORMAT( 2517 * //' (CIST1) INTERNAL ERROR, we only got',I4,' start vectors', 2518 * /' ',I4,' were requested', 2519 * /' --- action: revise code in SIRCI.CIST1 or modify problem.') 2520 CALL QUIT('CIST1 internal error, found too few start vectors.') 2521 END IF 2522 NCSIM = MCSIM 2523 deallocate(NLIST) 2524 CALL QEXIT('CIST1 ') 2525 RETURN 2526 END 2527C /* Deck cifnmn */ 2528 SUBROUTINE CIFNMN(NELMNT,IPLACE,VEC,NVEC,WRK,LWRK) 2529C 2530C 12-Aug-1990 hjaaj, based on CIST 2531C (Find minimum elemnts) 2532C 2533C Purpose: 2534C Return in IPLACE addresses on lowest NELMNT elements in VEC. 2535C 2536#include "implicit.h" 2537 DIMENSION VEC(NVEC), IPLACE(NELMNT), WRK(LWRK) 2538C 2539 IF (LWRK .LT. NELMNT) THEN 2540 CALL QUIT('CIFNMN: Insufficient memory (LWRK .lt. NELMNT)') 2541 END IF 2542 IF (NELMNT .GT. NVEC) THEN 2543 CALL QUIT('CIFNMN ERROR: NELMNT .gt. NVEC') 2544 END IF 2545C 2546C Sort first NELMNT elements of VEC 2547C 2548 WRK(1) = VEC(1) 2549 IPLACE(1) = 1 2550 DO 120 I = 2,NELMNT 2551 VECI = VEC(I) 2552 DO 115 J = 1,(I-1) 2553 IF (WRK(J) .GT. VECI) THEN 2554 DO 112 K = I,(J+1),-1 2555 WRK(K) = WRK(K-1) 2556 IPLACE(K) = IPLACE(K-1) 2557 112 CONTINUE 2558 IPLACE(J) = I 2559 WRK(J) = VECI 2560 GO TO 120 2561 END IF 2562 115 CONTINUE 2563 IPLACE(I) = I 2564 WRK(I) = VECI 2565 120 CONTINUE 2566C 2567C Find lowest elements by insertion sort 2568C 2569 XHGH = WRK(NELMNT) 2570 DO 140 I = NELMNT+1,NVEC 2571 IF (VEC(I).GE.XHGH) GO TO 140 2572 DO 130 J = 1,NELMNT 2573 IF (VEC(I) .LT. WRK(J)) THEN 2574 DO 135 K = NELMNT,(J+1),-1 2575 WRK(K) = WRK(K-1) 2576 IPLACE(K) = IPLACE(K-1) 2577 135 CONTINUE 2578 IPLACE(J) = I 2579 WRK(J) = VEC(I) 2580 XHGH = WRK(NELMNT) 2581 GO TO 140 2582 END IF 2583 130 CONTINUE 2584 140 CONTINUE 2585 RETURN 2586 END 2587C /* Deck cirsps */ 2588 SUBROUTINE CIRSPS(ECIREF,EMY,CMO,FCAC,H2AC,HDIAG,INDXCI,WRK,LFREE) 2589C 2590C 17-Jun-1987 hjaaj 2591C Revised 13-Jun-1991 hjaaj : also calculate and write DV 2592C 2593C Purpose: 2594C Interface to RESPONSE for CI response calculations: 2595C write information to LUSIFC ("SIRIFC") needed for response calculation. 2596C 2597C The following records are written: 2598C 2599C 1) POTNUC,EMY,EACTIV,ECIREF,ISTACI,ISPIN,NACTEL,LSYM,MS2 2600C 2) NISHT,NASHT,... 2601C 3) CMO 2602C 4) CREF 2603C 5) DV 2604C 6) FCAC 2605C 7) H2AC 2606C 8) label "CIDIAG2 " 2607C 9) HDIAG (diagonal of CI matrix) 2608C 10) label "SIRFLAGS" 2609C 11) FLAG(1:NFLAG) 2610C 2611#include "implicit.h" 2612#include "priunit.h" 2613 DIMENSION CMO(*), FCAC(*), H2AC(*), HDIAG(*), INDXCI(*), WRK(*) 2614#include "iratdef.h" 2615#include "dummy.h" 2616C 2617C Used from common blocks: 2618C 2619C INFINP : ISTACI,ISPIN,NACTEL,FLAG(1:NFLAG) 2620C INFORB : NSYM,NISHT,NASHT,NNASHX,NNASHY,NOCCT,NNOCCX,NCMOT,... 2621C INFVAR : NCONF 2622C INFTAP : LUIT1, LUSIFC 2623C CBGETD : IADH2 2624C SPINFO : MS2 2625C 2626#include "maxorb.h" 2627#include "infinp.h" 2628#include "inforb.h" 2629#include "infvar.h" 2630#include "inftap.h" 2631#include "infpri.h" 2632#include "cbgetdis.h" 2633#include "spinfo.h" 2634C 2635 CHARACTER*8 LAB123(3), TABLE(5) 2636 DATA LAB123/'********','********','********'/ 2637 DATA TABLE /'CIDIAG2 ','CIRESPON','SIRFLAGS','EODATA ', 2638 * 'CREFDETS'/ 2639C 2640 CALL QENTER('CIRSPS') 2641 CALL GETDAT(LAB123(2),LAB123(3)) 2642C ... place date in LAB123(2) and time in LAB123(3) 2643C 2644C *** Now create LUSIFC and write ... 2645C 2646 CALL GPOPEN(LUSIFC,FNSIFC,'UNKNOWN',' ','UNFORMATTED',IDUMMY, 2647 & .FALSE.) 2648C 2649C Retrieve CREF 2650C Calculate DV 2651C 2652 KCREF = 1 2653 KDV = KCREF + NCONF 2654 KWRK1 = KDV + NNASHX 2655 LWRK1 = LFREE + 1 - KWRK1 2656C 2657 REWIND LUIT1 2658 CALL MOLLAB('STARTVEC ',LUIT1,LUPRI) 2659 DO 100 I = 1, (ISTACI-1) 2660 READ (LUIT1) 2661 100 CONTINUE 2662 CALL READT(LUIT1,NCONF,WRK(KCREF)) 2663C 2664 CALL MAKDV(WRK(KCREF),WRK(KDV),INDXCI,WRK(KWRK1),LWRK1) 2665C 2666C Write SIRIFC 2667C 2668 LBSIFC = 'CIRESPON' 2669 WRITE (LUSIFC) LAB123,LBSIFC 2670 ECACTV = ECIREF - EMY 2671 WRITE (LUSIFC) POTNUC,EMY,ECACTV,ECIREF,ISTACI,ISPIN,NACTEL, 2672 & LSYM,MS2 2673 WRITE (LUSIFC) NISHT,NASHT,NOCCT,NORBT,NBAST,NCONF,NWOPT,NWOPH, 2674 * NCDETS,NCMOT,NNASHX,NNASHY,NNORBT,N2ORBT, 2675 & NSYM, MULD2H, NRHF, NFRO, 2676 & NISH, NASH, NORB, NBAS, 2677 & NELMN1, NELMX1, NELMN3, NELMX3, MCTYPE, 2678 & NAS1, NAS2, NAS3 2679C 2680 NC4 = MAX(4,NCONF) 2681 NCMOT4 = MAX(4,NCMOT) 2682 MMASHX = MAX(4,NNASHX) 2683C 2684 CALL WRITT (LUSIFC,NCMOT4,CMO) 2685 CALL WRITT (LUSIFC,NC4,WRK(KCREF)) 2686 CALL WRITT (LUSIFC,MMASHX,WRK(KDV)) 2687 CALL WRITT (LUSIFC,MMASHX,FCAC) 2688 IF (IADH2 .GE. 0) THEN 2689#if defined (VAR_NEWCODE) 2690... 27-Jun-1987/hjaaj -- noget lignende : MAERKE 2691 DO 100 IJ = 1,NNASHX 2692 CALL READDX(LUH2AC,IADH2+IJ,IRAT*NNASHX,H2AC) 2693 CALL WRITT (LUSIFC,MMASHX,H2AC) 2694 100 CONTINUE 2695... men der mangler en label som fortaeller response at H2AC on disk. 2696#else 2697 CALL QUIT('SIRCI.CIRSPS: ".DISKH2" not implemented here yet.') 2698#endif 2699 ELSE 2700 CALL WRITT (LUSIFC,(MMASHX*MMASHX),H2AC) 2701 END IF 2702 WRITE (LUSIFC) LAB123,TABLE(1) 2703 CALL WRITT (LUSIFC,NC4,HDIAG) 2704 WRITE (LUSIFC) LAB123,TABLE(3) 2705 WRITE (LUSIFC) (FLAG(I),I=1,NFLAG) 2706C 2707C Write CREF in determinants here, if NCDETS .gt. NCONF 2708C 2709 IF (NCDETS .GT. NCONF) THEN 2710 IF (FLAG(27)) THEN 2711 WRITE (LUPRI,*) 2712 & 'WR_SIRIFC ERROR, .DETERMINANTS but NCDETS.gt.NCONF:', 2713 & NCDETS, NCONF 2714 CALL QUIT('CIRSPS: NCDETS .gt. NCONF for .DETERMINANTS') 2715 END IF 2716 WRITE (LUSIFC) LAB123,TABLE(5) 2717 KCDET = KDV 2718 KWRK1 = KCDET + NCDETS 2719 LWRK1 = LFREE - KWRK1 2720 CALL GETDETS(LSYM,NCONF,NCDETS,INDXCI,WRK(KCREF),WRK(KCDET), 2721 & WRK(KWRK1),LWRK1) 2722 NC4 = MAX(4,NCDETS) 2723 CALL WRITT(LUSIFC,NC4,WRK(KCDET)) 2724 END IF 2725C 2726C Write EODATA 2727C 2728 WRITE (LUSIFC) LAB123,TABLE(4) 2729C 2730C *** end of subroutine CIRSPS 2731C 2732 CALL GPCLOSE(LUSIFC,'KEEP') 2733 CALL QEXIT('CIRSPS') 2734 RETURN 2735 END 2736 SUBROUTINE GETDVREF(ci_program,DVREF,IST,INDXCI,WRK,LWRK) 2737C 2738C Get DV_ref for CI-DFT and MC-DFT models. 2739C 2740C (c) jkp + hjaaj July 2003 2741C 2742C ===================================================== 2743C Construct active one-electron density matrix 2744C DVREF for CREF vector no. IST 2745C ===================================================== 2746C 2747 use lucita_mcscf_ci_cfg 2748 use lucita_mcscf_srdftci_cfg 2749#include "implicit.h" 2750C 2751 DIMENSION DVREF(*), INDXCI(*), WRK(LWRK) 2752C 2753C Used from common blocks: 2754C infvar: NCONF 2755C inftap: LUIT1 2756C 2757#include "infvar.h" 2758#include "inftap.h" 2759#include "priunit.h" 2760 character (len=9) :: ci_program 2761 integer, parameter :: luci_cvec = 99 2762C 2763C 2764 IF (IST .LE. 0) THEN 2765 CALL QUIT('GETDVREF called with IST .le. 0') 2766 END IF 2767 KCREF = 1 2768 KW1 = KCREF + NCONF 2769 LW1 = LWRK - KW1 2770 IF (LW1 .LE. 0) CALL ERRWRK('GETDVREF',-KW1,LWRK) 2771C 2772 if(ci_program .eq. 'LUCITA ')then 2773 open(file='srdft-lucita.cvec',unit=luci_cvec,status='old', 2774 & form='unformatted',action='read',position='rewind') 2775 call skpvcd(luci_cvec,ist-1,wrk(kcref),-1,-1) 2776 call dzero(wrk(kcref),nconf) 2777 call readvcd_exp(luci_cvec,wrk(kcref),0,-1) 2778! call wrtmatmn(wrk(kcref),1,nconf,1,nconf,lupri) 2779 close(luci_cvec,status='keep') 2780! set logical flag for new reference CI vector in WRK(KCREF) - 2781! in the interface to lucita we use this info to save/broadcast the new 2782! reference vector on lucita internal sequential/parallel MPI-I/O files 2783 vector_exchange_type1 = 1 2784 vector_update_mc2lu_lu2mc((2-1)*vector_exchange_types+ 2785 & vector_exchange_type1) = .true. 2786 2787 srdft_ci_1pdens_cref_restore = .true. 2788 2789 else 2790 REWIND LUIT1 2791 CALL MOLLAB('STARTVEC ',LUIT1,LUPRI) 2792 DO I = 1, (IST-1) 2793 READ (LUIT1) 2794 END DO 2795 CALL READT(LUIT1,NCONF,WRK(KCREF)) 2796 end if 2797 2798 CALL MAKDV(WRK(KCREF),DVREF,INDXCI,WRK(KW1),LW1) 2799 RETURN 2800 END 2801!****************************************************************************** 2802 2803 SUBROUTINE GETUDVREF(ci_program,UDVREF,IST,INDXCI,WRK,LWRK) 2804C 2805C Get DV_ref for CI-DFT and MC-DFT models. 2806C 2807C (c) jkp + hjaaj July 2003 2808C 2809C ===================================================== 2810C Construct active one-electron density matrix 2811C DVREF for CREF vector no. IST 2812C ===================================================== 2813C 2814 use lucita_mcscf_ci_cfg 2815 use lucita_mcscf_srdftci_cfg 2816#include "implicit.h" 2817C 2818 DIMENSION UDVREF(*), INDXCI(*), WRK(LWRK) 2819C 2820C Used from common blocks: 2821C infvar: NCONF 2822C inftap: LUIT1 2823C 2824#include "infvar.h" 2825#include "inftap.h" 2826#include "priunit.h" 2827 character (len=9) :: ci_program 2828 integer, parameter :: luci_cvec = 99 2829C 2830C 2831 IF (IST .LE. 0) THEN 2832 CALL QUIT('GETUDVREF called with IST .le. 0') 2833 END IF 2834 KCREF = 1 2835 KW1 = KCREF + NCONF 2836 LW1 = LWRK - KW1 2837 IF (LW1 .LE. 0) CALL ERRWRK('GETUDVREF',-KW1,LWRK) 2838C 2839 if(ci_program .eq. 'LUCITA ')then 2840 open(file='srdft-lucita-final.cvec',unit=luci_cvec,status='old', 2841 & form='unformatted',action='read',position='rewind') 2842 call skpvcd(luci_cvec,ist-1,wrk(kcref),-1,-1) 2843 call dzero(wrk(kcref),nconf) 2844 call readvcd_exp(luci_cvec,wrk(kcref),0,-1) 2845! call wrtmatmn(wrk(kcref),1,nconf,1,nconf,lupri) 2846 close(luci_cvec,status='keep') 2847! set logical flag for new reference CI vector in WRK(KCREF) - 2848! in the interface to lucita we use this info to save/broadcast the new 2849! reference vector on lucita internal sequential/parallel MPI-I/O files 2850 vector_exchange_type1 = 1 2851 vector_update_mc2lu_lu2mc((2-1)*vector_exchange_types+ 2852 & vector_exchange_type1) = .true. 2853 2854 srdft_ci_1pdens_cref_restore = .true. 2855 2856 else 2857 REWIND LUIT1 2858 CALL MOLLAB('STARTVEC ',LUIT1,LUPRI) 2859 DO I = 1, (IST-1) 2860 READ (LUIT1) 2861 END DO 2862 CALL READT(LUIT1,NCONF,WRK(KCREF)) 2863 end if 2864 2865 CALL MAKDV(WRK(KCREF),UDVREF,INDXCI,WRK(KW1),LW1) 2866 RETURN 2867 END 2868!****************************************************************************** 2869 2870 SUBROUTINE GETS2REF(PVREF,DVREF,IST,INDXCI,WRK,LWRK) 2871C 2872C Get DV_ref for CI-DFT and MC-DFT models. 2873C 2874C (c) jkp + hjaaj July 2003 2875C 2876C ===================================================== 2877C Construct active one- and two-electron density matrix 2878C for CREF vector no. IST 2879C ===================================================== 2880C 2881 use lucita_mcscf_ci_cfg 2882 use lucita_mcscf_srdftci_cfg 2883#include "implicit.h" 2884C 2885 DIMENSION DVREF(*), PVREF(*), INDXCI(*), WRK(LWRK) 2886C 2887C Used from common blocks: 2888C infvar: NCONF 2889C inftap: LUIT1 2890C 2891#include "infvar.h" 2892#include "inftap.h" 2893#include "priunit.h" 2894 integer, parameter :: luci_cvec = 99 2895C 2896C 2897 IF (IST .LE. 0) THEN 2898 CALL QUIT('GETS2REF called with IST .le. 0') 2899 END IF 2900 KCREF = 1 2901 KW1 = KCREF + NCONF 2902 LW1 = LWRK - KW1 2903 IF (LW1 .LE. 0) CALL ERRWRK('GETS2REF',-KW1,LWRK) 2904C 2905 open(file='srdft-lucita-final.cvec',unit=luci_cvec,status='old', 2906 & form='unformatted',action='read',position='rewind') 2907 call skpvcd(luci_cvec,ist-1,wrk(kcref),-1,-1) 2908 call dzero(wrk(kcref),nconf) 2909 call readvcd_exp(luci_cvec,wrk(kcref),0,-1) 2910! call wrtmatmn(wrk(kcref),1,nconf,1,nconf,lupri) 2911 close(luci_cvec,status='keep') 2912! set logical flag for new reference CI vector in WRK(KCREF) - 2913! in the interface to lucita we use this info to save/broadcast the new 2914! reference vector on lucita internal sequential/parallel MPI-I/O files 2915 vector_exchange_type1 = 1 2916 vector_update_mc2lu_lu2mc((2-1)*vector_exchange_types+ 2917 & vector_exchange_type1) = .true. 2918 2919 srdft_ci_1pdens_cref_restore = .true. 2920 2921 CALL MAKDM(WRK(KCREF),DVREF,PVREF,INDXCI,WRK(KW1),LW1) 2922 2923 END 2924!****************************************************************************** 2925 2926 SUBROUTINE CIOPT(EMY,JCONV,ICICTL,NCROOT,MAXITC,CMO,INDXCI, 2927 & EJCSR,EJVSR,EDSR,ESRDFT,EMYDFT,EMYDFTAUX, 2928 & NCLIN,ITCI, 2929 & TIMLIN,TIMCIT,NCONF4,LSVECS,MAXSIM, 2930 & residual_c,energy_c, 2931 & KFCAC,KH2AC,KCDIAG,KW1,KCIDIA,KSRAC, 2932 & KCVECS,KSVECS,KW3,KW2,KREDCI,KEVEC, 2933 & LW1,LW2,LW3,WRK,LWRK) 2934 2935C 2936C 2937C 23-08-2012 Manu: added one more argument, EMYDFTAUX, in order to compute 2938C auxiliary long-range interacting energies 2939 2940#include <implicit.h> 2941 DIMENSION CMO(*), INDXCI(*), WRK(*), residual_c(*), energy_c(*) 2942 REAL*8 ESRDFT(3) 2943C 2944#include <iratdef.h> 2945#include <thrldp.h> 2946 PARAMETER ( D0 = 0.0D0 , DP5 = 0.5D0, D1 = 1.0D0 ) 2947#include <dummy.h> 2948C 2949C Used from common blocks: 2950C CICB01 : NCIRED,NJCR,JCROOT(*),THRCIT 2951C INFINP : POTNUC,FLAG(*),LSYM 2952C INFVAR : NCONF 2953C INFORB : NNASHX 2954C INFDIM : LCINDX,LACIMX,LBCIMX,LPHPMX,MAXPHP 2955C INFTAP : LUIT1,LUIT3,LUIT5 2956C 2957#include <maxorb.h> 2958#include <priunit.h> 2959#include <cicb01.h> 2960#include <infinp.h> 2961#include <infvar.h> 2962#include <inforb.h> 2963#include <infdim.h> 2964#include <inftap.h> 2965#include <infpri.h> 2966 LOGICAL CINMEM 2967 2968C *** Calculate the diagonal of the CI matrix. 2969C 2970 IF (IPRSIR .GT. 10) WRITE (LUPRI,'(/a)') ' --- Calling CIDIAG ...' 2971 CALL CIDIAG(LSYM,.FALSE.,WRK(KFCAC),WRK(KH2AC),INDXCI, 2972 & WRK(KCDIAG),WRK(KW1),LW1) 2973C CALL CIDIAG(ICSYM,NOH2,FCAC,H2AC,XNDXCI,DIAGC,WRK,LFREE) 2974C 2975 EMY = EMY + POTNUC 2976 DO 50 I = 0,NCONF-1 2977 WRK(KCDIAG+I) = WRK(KCDIAG+I) + EMY 2978 50 CONTINUE 2979 IF (IPRSIR .GT. 10) WRITE(LUPRI,'(/a)') ' --- CIDIAG finished ...' 2980 2981#if defined (VAR_CIHAMTST) 2982C 2983C --- 890131 -- explicit CI Hamiltonian matrix --- 2984C 2985 IF ( FLAG(43) ) THEN 2986 WRITE (LUPRI,'(//A/)') ' CICTL: test call of CIHAM' 2987 IWAY = 1 2988 MXPRDT = MIN(20,NCONF) 2989 IPRHAM = 10 2990C 2991 KIDET = KW1 2992 KCIDIA = KIDET + MXPRDT/IRAT + 1 2993 KWHAM = KCIDIA + NCONF 2994 LWHAM = LWRK - KWHAM 2995 CALL DCOPY(NCONF,WRK(KCDIAG),1,WRK(KCIDIA),1) 2996 CALL CIHAM(IWAY,MXPRDT,LSYM,WRK(KIDET),WRK(KCIDIA),NCONF, 2997 * EMY,WRK(KFCAC),WRK(KH2AC),INDXCI, 2998 * WRK(KWHAM),LWHAM,IPRHAM) 2999C CALL CIHAM(IWAY,MXPRDT,ICSYM,IDET,CIDIA,NCONF, 3000C * ECORE,FCAC,H2AC,INDXCI,WRK,LWRK,IPRINT) 3001 END IF 3002#endif 3003C 3004C Get PHP CI matrix (with correct CI energies) 3005C 3006 IF (MAXPHP .GE. 1) THEN 3007 IPWAY = 1 3008 IPRPHP = IPRI6 - 3 3009 IF (IPRSIR .GT. 10) WRITE (LUPRI,'(/a,i0)') 3010 & ' --- Calling PHPGET with MAXPHP = ',MAXPHP 3011 CALL PHPGET(LSYM,NCONF,INDXCI,WRK(KFCAC),WRK(KH2AC), 3012 & WRK(KCDIAG),EMY,IPWAY,IPRPHP,WRK(KW1),LW1) 3013 END IF 3014C 3015C 3016C ****************************** 3017C *** GET INITIAL CI VECTORS *** 3018C ****************************** 3019C 3020C CIST constructs NCROOT start vectors for CI optimization. 3021C The C-vectors are constructed in CIST or from the result of 3022C another run (another geometry). 3023C 3024 NCSIM = NCROOT 3025 IF (IPRSIR .GT. 10) WRITE (LUPRI,'(/a)') ' --- Calling CIST ...' 3026 CALL CIST(ICICTL,NCSIM,CMO,WRK(KCDIAG),EMY,WRK(KW1),LW1) 3027 IF (NCSIM .LE. 0) THEN 3028 WRITE (LUPRI,'(//A,I5)') 3029 * ' *** ERROR after CIST, NCSIM =',NCSIM 3030 CALL QTRACE(LUPRI) 3031 CALL QUIT('*** ERROR (CICTL) NCSIM .le. 0 after CIST') 3032 END IF 3033 IF (IPRSIR .GT. 10) WRITE (LUPRI,'(/a)') ' --- CIST finished ...' 3034C 3035C ******************************* 3036 IF (MAXITC .LE. 0) THEN 3037 ITCI = 0 3038 RETURN 3039 END IF 3040C 3041C ***************************************************************** 3042C CI-DFT contributions 3043C ***************************************************************** 3044C 3045 IF (DOCISRDFT) THEN 3046#ifndef MOD_SRDFT 3047 call quit('srdft not implemented in this version') 3048#else 3049 KDVREF = KW2 3050 KFSR = KDVREF + NNASHX 3051 KW2A = KFSR + NNORBT 3052 LW2A = LWRK - KW2A 3053 IF (ISTACI .EQ. 0) THEN 3054 CALL DZERO(WRK(KDVREF),NNASHX) 3055 ELSE 3056 IST = ABS(ISTACI) 3057 CALL GETDVREF('SIRIUS-CI',WRK(KDVREF),IST, 3058 & INDXCI,WRK(KW2A),LW2A) 3059 END IF 3060C 3061 CALL SRFMAT(WRK(KFSR),CMO,WRK(KDVREF),EJCSR,EJVSR,EDSR,ESRDFT, 3062 & EMYDFTAUX,UEJCVSR,WRK(KW2A),LW2A,IPRI6) 3063C ... Correct EMY and CI diagonal 3064 EMYDFT = EJCSR + EJVSR + EDSR + ESRDFT(1) 3065C write(lupri,*) 'just after srfmat' 3066C write(lupri,*) 'EMYDFTAUX = ', EMYDFTAUX 3067 EMY = EMY + EMYDFT 3068C ... extract and save SR over active indices 3069 CALL GETAC(WRK(KFSR),WRK(KSRAC)) 3070! call wrtmatmn(wrk(ksrac),1,nnashx,1,nnashx,lupri) 3071! call wrtmatmn(wrk(kfcac),1,nnashx,1,nnashx,lupri) 3072C ... add the contributions over the active indices 3073 CALL DAXPY(NNASHX,D1,WRK(KSRAC),1,WRK(KFCAC),1) 3074 !call wrtmatmn(wrk(kfcac),1,nnashx,1,nnashx,lupri) 3075C ... fold SRAC for 3076C 1) Adding it to CI-diagonal to get better pre-conditioning 3077C 2) Later printing of CIDFT specific energy contributions 3078 CALL DSPTGE(NASHT,WRK(KSRAC),WRK(KFSR)) 3079 CALL DGEFSP(NASHT,WRK(KFSR),WRK(KSRAC)) 3080 ESRDV = DDOT(NNASHX,WRK(KDVREF),1,WRK(KSRAC),1) 3081C write(lupri,*) 'esrdv... ',esrdv 3082 DO I = 0,NCONF-1 3083 WRK(KCDIAG+I) = WRK(KCDIAG+I) + EMYDFT + ESRDV 3084 END DO 3085#endif 3086 ENDIF 3087C ******************************* 3088 CINMEM = (NCSIM .LE. MAXSIM) 3089C 3090 IF (NCSIM .GT. MXCIRT) THEN 3091 WRITE (LUPRI,'(//A,I5,A,I5)') 3092 * ' *** ERROR after CIST, NCSIM =',NCSIM, 3093 * ', maximum (MXCIRT) =',MXCIRT 3094 CALL QTRACE(LUPRI) 3095 CALL QUIT('*** ERROR (CICTL) MXCIRT parameter .lt. NCSIM') 3096 END IF 3097 DO 100 I = 1,NCSIM 3098 JCROOT(I) = I 3099 residual_c(i) = D0 3100 100 CONTINUE 3101C 3102 REWIND LUIT1 3103 CALL MOLLAB('STARTVEC',LUIT1,LUPRI) 3104 REWIND LUIT3 3105 MXCSIM = 2*MAXSIM - 1 3106C reuse space for orthonormalization, 3107C except reserving last NCONF elements for CPREV 3108 KCPREV = KCVECS + MXCSIM*NCONF 3109 NCPREV = 0 3110 DO 490 I = 1,NCSIM,MXCSIM 3111 NXSIM = MIN(MXCSIM,(NCSIM-I+1)) 3112 JCVECS = KCVECS 3113 DO 420 J = 1,NXSIM 3114 CALL READT(LUIT1,NCONF,WRK(JCVECS)) 3115 JCVECS = JCVECS + NCONF 3116 420 CONTINUE 3117 THRLDV = NCONF * THRLDP 3118 CALL CIORT(NXSIM,NCPREV,LUIT3,JCROOT(NCPREV+1), 3119 * WRK(KCVECS),THRLDV,WRK(KCPREV)) 3120 490 CONTINUE 3121 IF (NCPREV .NE. NCSIM) THEN 3122 WRITE (LUPRI,'(//A,I8,A)') 3123 * ' *** ERROR ***',(NCSIM-NCPREV),' OF THE TRIAL CI VECTORS'// 3124 * ' WERE LINEAR DEPENDENT (CICTL).' 3125 CALL QUIT('***ERROR (CICTL) LIN.DEP. AMONG INITIAL CI '// 3126 & 'TRIAL VEC.S') 3127 END IF 3128C 3129C ***************************************************************** 3130C CI iteration loop 3131C ***************************************************************** 3132C 3133C 3134 TIMLIN = D0 3135 TIMCIT = SECOND() 3136 NCLIN = 0 3137 ITCI = 0 3138C 3139 500 ITCI = ITCI + 1 3140 IF (IPRI4 .GT. 5) 3141 * WRITE (LUW4,'(//A,I4,/)') ' *** CI ITERATION NO.',ITCI 3142 IF (LUPRI.NE.LUW4 .AND. IPRI6 .GE. 5) WRITE (LUPRI,'(//A,I4,/)') 3143 * ' *** CI ITERATION NO.',ITCI 3144 CALL FLSHFO(LUW4) 3145 IF (LUPRI.NE.LUW4) CALL FLSHFO(LUPRI) 3146 REWIND LUIT5 3147 DO 220 I = 1,NCIRED 3148 220 READ (LUIT5) 3149C 3150C REPEAT UNTIL (NCSIM .LE. 0) 3151C 3152 240 CONTINUE 3153 NXSIM = MIN(MAXSIM,NCSIM) 3154 IF (.NOT. CINMEM) THEN 3155 REWIND LUIT3 3156 DO 310 I = 1,NCIRED 3157 310 READ (LUIT3) 3158 DO 320 I = 1,NXSIM 3159 JCVECS = KCVECS + (I-1)*NCONF 3160 CALL READT(LUIT3,NCONF,WRK(JCVECS)) 3161 320 CONTINUE 3162 END IF 3163 TIMLIN = TIMLIN - SECOND() 3164 NCLIN = NCLIN + NXSIM 3165 CALL CILIN(NXSIM,WRK(KCVECS),WRK(KFCAC),WRK(KH2AC), 3166 * INDXCI,WRK(KSVECS),LSVECS) 3167 TIMLIN = TIMLIN + SECOND() 3168 DO 340 I = 1,NXSIM 3169 JSVECS = KSVECS + (I-1)*NCONF 3170 CALL WRITT(LUIT5,NCONF4,WRK(JSVECS)) 3171 340 CONTINUE 3172 NCIRED = NCIRED + NXSIM 3173 CALL CIRED(1,NCIRED,WRK(KREDCI),NXSIM,WRK(KCVECS), 3174 * WRK(KSVECS),DUMMY,DUMMY,WRK(KW3),LW3) 3175 NCSIM = NCSIM - NXSIM 3176 IF (NCSIM .GT. 0) GO TO 240 3177C ^----------------------------- END REPEAT 3178C 3179 CALL CIRED(2,NCIRED,WRK(KREDCI),0,WRK(KCVECS), 3180 * DUMMY,energy_c,WRK(KEVEC),WRK(KSVECS),LSVECS) 3181C CALL CIRED(ICTL,NCIRED,REDCI,NCNEW,CVECS,SVECS,EVAL,EVEC,WRK,LWRK) 3182 IF (P6FLAG(10) .OR. P6FLAG(34)) THEN 3183 WRITE (LUPRI,'(/A,I3)') 3184 * ' The reduced CI matrix, iteration',ITCI 3185 IF (P6FLAG(34)) CALL OUTPAK(WRK(KREDCI),NCIRED,-1,LUPRI) 3186 MCROOT = MIN(NCROOT+2,NCIRED) 3187 WRITE (LUPRI,'(/I4,A//,(10X,1P,6D15.6))') 3188 * MCROOT,' lowest eigenvalues of reduced CI matrix :', 3189 * (energy_c(i),I=1,MCROOT) 3190 IF (P6FLAG(34) .OR. IPRI6 .GT. 10) THEN 3191 WRITE (LUPRI,'(/A,I3,A)') 3192 * ' - and the lowest',NCROOT,' eigenvectors :' 3193 CALL OUTPUT(WRK(KEVEC),1,NCIRED,1,NCROOT, 3194 * NCIRED,NCIRED,-1,LUPRI) 3195 END IF 3196 END IF 3197 IF (P6FLAG(51)) THEN 3198C ... check symmetry of reduced CI matrix and CI overlap 3199 CALL CIRCHK(WRK(KSVECS),LSVECS) 3200 END IF 3201 ESHIFT = D0 3202 NCSIM = NCROOT 3203 CALL CINEX(energy_c,WRK(KEVEC),EMY,ESHIFT,WRK(KCDIAG), 3204 * NCSIM,JCONV,residual_c,WRK(KW2),LW2) 3205 IF (JCONV .GT. 0) THEN 3206 IF (NCSIM .EQ. 0) THEN 3207 GO TO 8100 3208 ELSE IF (NCSIM .LT. 0) THEN 3209 CINMEM = .FALSE. 3210 NCSIM = -NCSIM 3211 ELSE 3212 CINMEM = (NCSIM .LE. MAXSIM) 3213 END IF 3214 IF (ITCI .GE. MAXITC) GO TO 8000 3215 IF (NCIRED+NCSIM .GE. MAXRC ) GO TO 8200 3216 GO TO 500 3217 END IF 3218 WRITE (LUW4,'(//A,I3,A/)') 3219 * ' *** CI converged in',ITCI,' iterations.' 3220 IF (LUPRI .NE. LUW4) WRITE (LUPRI,'(//A,I3,A/)') 3221 * ' *** CI converged in',ITCI,' iterations.' 3222 GO TO 9000 3223 8000 CONTINUE 3224 WRITE (LUPRI,'(//A,I5/I15,A)') 3225 * ' *** Reached maximum number of CI iterations:',MAXITC, 3226 * JCONV,' CI roots are not converged.' 3227 GO TO 9000 3228 8100 CONTINUE 3229 NWARN = NWARN + 1 3230 WRITE (LUPRI,'(//A/I15,A)') 3231 * ' *** WARNING, all CI trial vectors linear dependent', 3232 * JCONV,' CI roots are not converged.' 3233 GO TO 9000 3234 8200 CONTINUE 3235 NWARN = NWARN + 1 3236 WRITE (LUPRI,'(//A,I5/I15,A)') 3237 * ' *** WARNING, reached maximum dimension of reduced matrix:', 3238 * MAXRC,JCONV,' CI roots are not converged.' 3239 GO TO 9000 3240 3241 9000 CONTINUE 3242 3243! finish timing 3244 TIMCIT = SECOND() - TIMCIT 3245 3246C write(lupri,*) 'right at the end of ine ciopt' 3247C write(lupri,*) 'EMYDFTAUX = ', EMYDFTAUX 3248 END SUBROUTINE CIOPT 3249! - end of file sirci.F - 3250