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 FILE : siropt.F 20C 21C SIROPT controls second order MCSCF or SCF optimization in SIRIUS. 22C RHFENR print RHF orbital energies 23C OPTST starting the optimization 24C SIRSTP step control 25C TIMOPT timing statistics of optimization 26C WR_SIRIFC write SIRIFC interface file from SIRIUS to RESPONS and ABACUS 27C 28C=========================================================================== 29C /* Deck siropt */ 30 SUBROUTINE SIROPT(WRK,LFREE,ICONV) 31C 32C 31-Oct-1983-V1/6-May-1984-V2 -- Hans Jorgen Aa. Jensen 33C Revisions 34C 6-Aug-1984/16-May-1984 - hjaaj 35C 14-May-1985 - hjaaj (update -- sirupd) 36C 10-Nov-1985 - hjaaj (absorption -- sirorb) 37C 14-May-1986 - hjaaj (nr eq.s -- sirnr) 38C 39C Purpose: Control MCSCF optimization. 40C 41C MOTECC-90: The purpose of this module, SIROPT, and the algorithms used 42C are described in Chapter 8 Section C.0 of MOTECC-90 43C "Implementation" 44C 45C Output: 46C ICONV .eq.0 MCSCF not converged. 47C .gt.0 MCSCF converged. 48C ----- ICONV can be used for conditional call of 49C subsequent steps requiring a converged w.f. 50C 51#ifdef VAR_MPI 52 use par_mcci_io 53#endif 54 use pelib_interface, only: use_pelib, pelib_ifc_fock, 55 & pelib_ifc_grad, pelib_ifc_energy 56 57#include "implicit.h" 58 59C 60 DIMENSION WRK(LFREE) 61C 62#include "thrldp.h" 63 PARAMETER (MAXBCK = 5) 64 PARAMETER (FAKABS = 0.40D0, STPABS = 0.20D0) 65 PARAMETER (STPLCL = 0.15D0, SHFLCL = 1.D-3) 66 PARAMETER (D0 = 0.D0, D1 = 1.D0, D2 = 2.0D0) 67#include "iratdef.h" 68#include "dummy.h" 69C 70C used from common blocks: 71C INFINP: LSYM,FLAG(*),...,ITRLVL,ITRFIN,THRPWF,THRCGR,MAXUIT, 72C INERSI,... 73C INFVAR: NCONF,NWOPT,NVAR,NVARH 74C INFORB: NCMOT,NNASHX,... 75C INFOPT: EMCSCF,... 76C INFDIM: LCINDX,MAXRL,? 77C INFTAP: LUIT1,? 78C CBIHR2: IFTHRS,USRSCR 79C CBIREA: LMULBS 80C R12INT: R12CAL 81C 82#include "maxorb.h" 83#include "priunit.h" 84#include "infinp.h" 85#include "pcmlog.h" 86#include "infvar.h" 87#include "inforb.h" 88#include "infopt.h" 89#include "infdim.h" 90#include "inftap.h" 91#include "infpri.h" 92#include "itinfo.h" 93#include "gnrinf.h" 94#include "cbihr2.h" 95#include "cbirea.h" 96#include "r12int.h" 97C------------------------ 98#include "dfterg.h" 99#include "dftcom.h" 100C --- 101C 102C *** local: 103 CHARACTER*8 STAR8, RTNLBL(2) 104 CHARACTER*8 WF_TYPE 105 LOGICAL ABSORP, UPDATE, SOLVNT, GEOWLK, RHF_CALC 106 LOGICAL CHKSTP, DONR, NRALW, LOCAL 107C *** data: 108 DATA STAR8/'********'/ 109C 110 CALL QENTER('SIROPT') 111C 112C *** Start timing: 113C 114 IF (IPRSTAT .GT. 0) CALL TIMOPT('START',LUSTAT,WRK,LFREE) 115 CALL GETTIM(T_start,W_start) 116C 117C *** WF_TYPE 118C 119 RHF_CALC = (NASHT .LE. 1 .OR. HSROHF) 120 IF ( (MCTYPE.LE.0) .NEQV. RHF_CALC ) THEN 121 WRITE (LUPRI,'(//A,I5)') 122 & 'WARNING siropt: RHF_CALC neqv MCTYPE =',MCTYPE 123 END IF 124 IF (RHF_CALC) THEN 125 IF (DOHFSRDFT) THEN 126 WF_TYPE = 'HF-srDFT' 127 len_WF_TYPE = 8 128 ELSE IF (DODFT) THEN 129 ! radovan: added block 2014-12-09 130 ! my information is that this is not fully implemented 131 call quit('QC-SCF not implemented/tested for DFT') 132 WF_TYPE = 'QC-DFT' 133 len_WF_TYPE = 6 134 ELSE IF (HSROHF) THEN 135 WF_TYPE = 'HSROHF' 136 len_WF_TYPE = 6 137 ELSE 138 WF_TYPE = 'QC-HF' 139 len_WF_TYPE = 5 140 END IF 141 ELSE IF (DOMCSRDFT) THEN 142 WF_TYPE = 'MC-srDFT' 143 len_WF_TYPE = 8 144 ELSE 145 WF_TYPE = 'MCSCF' 146 len_WF_TYPE = 5 147 END IF 148 149 WRITE (LUW4,960) WF_TYPE(1:len_WF_TYPE) 150 flush(LUW4) 151 IF (LUPRI .NE. LUW4) THEN 152 WRITE (LUPRI,960) WF_TYPE(1:len_WF_TYPE) 153 flush(LUPRI) 154 END IF 155 960 FORMAT(//T9, 'SIRIUS ',A,' optimization (SIROPT)' 156 * /' ================================================'/) 157C 158C *** Initializations 159C convergence flag ICONV to "not converged" 160C UPDATE , if use UPDATE for MC optimization 161C 162 ICONV = 0 163 SOLVNT = FLAG(16) 164 165 ABSORP = ( FLAG(51) .OR. FLAG(52) .OR. FLAG(53) ) 166#if defined (VAR_NOEXCABS) 167 IF ( ABSORP .AND. ISTATE .GT. 1 ) THEN 168C WRITE (LUW4,'(//A/A/)') ' *** Orbital absorption requested,', 169C * ' orbital absorption is not currently defined for'// 170C * ' excited states and request is ignored.' 171C ABSORP = .FALSE. 172C FLAG(51) = .FALSE. 173C FLAG(52) = .FALSE. 174C FLAG(53) = .FALSE. 175 NWARN = NWARN + 1 176 WRITE (LUW4,'(//A/)') ' *** WARNING *** ABSORPTION FOR'// 177 * ' AN EXCITED STATE IS EXPERIMENTAL' 178 END IF 179#endif 180 IF (ABSORP .AND. MAXAPM .LE. 0) THEN 181 WRITE (LUW4,'(//A/A/)') ' *** Orbital absorption requested,', 182 * ' request is ignored because zero iterations specified.' 183 ABSORP = .FALSE. 184 FLAG(51) = .FALSE. 185 FLAG(52) = .FALSE. 186 FLAG(53) = .FALSE. 187 END IF 188CHJAaJ Nov 09: KTRLVL=0 does not work in TR1H2M in current Dalton 189CHJ IF (FLAG(51) .OR. FLAG(52)) THEN 190CHJ KTRLVL = 0 191CHJ ELSE 192 KTRLVL = ITRLVL 193CHJ END IF 194C If FLAG(39) always NR (e.g. for core hole calculations /860515-hj) 195C 921028-hjaaj: FLAG(39) copied to local variable, because 196C it may be reset to true in local region 197 NRALW = FLAG(39) 198 DONR = NRALW 199 UPDATE = FLAG(19) 200 ITRLSV = ITRLVL 201 LTRLVL = -999 202 NWOPS = NWOPH 203 NVARS = NVARH 204C 205C dec 99 hjaaj: 206C 1) save initial screening threshold for Direct Fock matrix construction 207C 2) set screening threshold to accuracy for converged gradient for first 208C iteration (only used if DIRFCK true) 209C 210 IFTHSV = IFTHRS 211 GNORM(3) = THRGRD 212C 213C *** Core allocation: 214C 215C.. 1.0 .. (OPTST, READMO, GETCIX, GRAD and SIRNEO) 216 KCINDX = 1 217 KCMO = KCINDX + LCINDX 218C.. 1.1.1 .. (GETCIX, TRACTL, READMO) 219 KWRK10 = KCMO + NCMOT 220 LWRK10 = LFREE + 1 - KWRK10 221C.. 1.2 .. (GRAD, SIRORB and SIRNEO) 222 IF (DFT_SPINDNS) THEN 223 NDV = 2 224 ELSE 225 NDV = 1 226 END IF 227 KGTOT = KWRK10 228 KDV = KGTOT + NVARH 229 KPV = KDV + NDV*NNASHX 230 KFC = KPV + NNASHX*NNASHX 231 KFV = KFC + NNORBT 232!HJJ - FRAN KFS added for the SR_DFT Spin Dens. contribution 233 KFS = KFV + NNORBT 234 IF (DFT_SPINDNS) THEN 235 KSRHYB = KFS + NNORBT 236 ELSE 237 KSRHYB = KFS 238 END IF 239!HJJ - FRAN END 240 IF (SRHYBR .AND. MCTYPE.GT.0) THEN 241 KFQ = KSRHYB + NNORBT 242 ELSE 243 KFQ = KSRHYB 244 END IF 245 KFCAC = KFQ + NASHT*NORBT 246 KCREF = KFCAC + NDV*NNASHX ! FCAC, FSAC for DFT_SPINDNS 247 KH2AC = KCREF + NCONF 248 KWRK12 = KH2AC + NNASHX*NNASHX 249C.. 1.2.1 .. (GRAD,SOLGRD,WR_SIRIFC) 250 KERLM = KWRK12 251 IF (SOLVNT) THEN 252 KWRKG = KERLM + 2*NLMSOL 253 ELSE 254 KWRKG = KWRK12 255 END IF 256 LWRKG = LFREE - (KWRKG-1) 257 IF (LWRKG.LT.0) CALL ERRWRK('SIROPT.GRAD',-KWRKG,LFREE) 258C.. 1.2.2 .. (SIRORB) 259 KWOCTL = KFCAC 260 LWOCTL = LFREE - KWOCTL 261 IF (ABSORP .AND. LWOCTL .LE. 0) 262 * CALL ERRWRK('SIROPT.SIRORB',-KWOCTL,LFREE) 263 IF (.NOT.FLAG(20)) THEN 264C.. 1.2.3 .. (SIRSTP, SIRNEO) 265 KIBNDX = KWRK12 266 KWNEO = KIBNDX + (MAXRL + MAXRTS*3)/IRAT + 1 267 LWRK2 = KWNEO + 3*NVAR 268C (KWNEO space for vectors as above; 269C ?????? estimate of space for DTV,PTV,FXC,...; MAERKE 270C 3*NVAR space for one Bvec, one Svec, and diagonal) 271 IF (LWRK2 .GT. LFREE) CALL ERRWRK('SIROPT.SIRNEO',-LWRK2,LFREE) 272 LWNEO = LFREE + 1 - KWNEO 273 END IF 274C.. 1.3 .. (SIRNEO) 275C.. 1.4 .. (SIRUPD) 276 KREFU = KWRK10 277 KGRDU = KREFU + NVAR 278 KWU = KGRDU + NVAR 279 LWU = LFREE - KWU 280C.. 281 NCONF4 = MAX(4,NCONF) 282C 283C Get CI information if active orbitals (and ci_program .eq. SIRIUS-CI) 284C 285 IF (.NOT.RHF_CALC) 286 & CALL GETCIX(WRK(KCINDX),LSYM,LSYM,WRK(KWRK10),LWRK10,0) 287C 288C **************************** 289C *** GET STARTING VECTORS *** 290C **************************** 291C 292C OPTST constructs LROOTS starting vectors for 293C the very first macroiteration. 294C The C-vectors are constructed in OPTST or from the result of 295C another run (another geometry). Later the resulting vectors 296C from the previous macroiteration is read in. 297C ITMAC and EMCOLD are initialized before CALL OPTST 298C so they can be defined by OPTST when restart. 299C 300C ITRLVL is temporarily set to KTRLVL, so we don't perform larger 301C integral transformation in OPTST than needed. 302C 303 ITMAC = 0 304 EMCOLD = D0 305 DEPRED = D0 306 STPLEN = D0 307 JTRLVL = ITRLVL 308 ITRLVL = KTRLVL 309C 310C jan00 hjaaj: set screening threshold for SIRCNO in OPTST for direct 311C SCF (canonical orbitals are not needed with high precision). 312C 313 IF (DIRFCK .AND. .NOT. USRSCR) THEN 314 IFTHRS = MIN( IFTHSV, 9 ) 315 IF (IPRI6 .GE. 0) THEN 316 WRITE (LUPRI,'(/A,I5)') 317 & ' Fock matrix screening setting for SIRCNO: IFTHRS =', 318 & IFTHRS 319 flush(LUPRI) 320 END IF 321 END IF 322 CALL OPTST(WRK(KCINDX),WRK(KWRK10),LWRK10) 323C CALL OPTST(INDXCI,WRK,LFREE) 324 ITRLVL = JTRLVL 325C 326C set parameters 327C ABSORP may be reset by OPTST, if restart 328C DONR may be reset by OPTST, if restart 329C CHKSTP false for first macro iteration, unless restart 330C 331 ABSORP = ( FLAG(51) .OR. FLAG(52) .OR. FLAG(53) ) 332 IF (.NOT. FLAG(54)) ABSORP = ABSORP .AND. (MAXABS .GT. 0) 333 DONR = FLAG(39) 334 GEOWLK = (ICI0 .EQ. 6) 335 IF ( GEOWLK ) THEN 336 JWSTEP = 1 337 CHKSTP = .FALSE. 338 ISTEP = 0 339 EMCGEO = EMCOLD 340 DEPGEO = DEPRED 341 THRGEO = MIN(1.D-2,0.1D0*SQRT( ABS(DEPGEO) )) 342 ELSE 343 JWSTEP = 0 344 IF (DEPRED .NE. D0) THEN 345 CHKSTP = .TRUE. 346 ISTEP = -1 347 ELSE 348 CHKSTP = .FALSE. 349 ISTEP = 0 350 END IF 351 END IF 352C 353C-841209-hj: force beta positive 354 BETA = ABS(BETA) 355 IF (BETA .EQ. D0) BETA = D1 356 IF (BETMIN .LE. D0 .OR. 357 * BETMIN .GT. BETMAX) BETMIN = D1 358 IF (BETMAX .LT. D1) BETMAX = D1 359C 360C 361 IF (IPRI6 .GT. 0) THEN 362 CALL GETTIM(T_1,W_1) 363 T_1 = T_1 - T_start 364 W_1 = W_1 - W_start 365 WRITE (LUPRI,'(/A,2F15.2)') 366 & ' --- CPU and WALL time used in SIROPT startup (seconds):', 367 & T_1,W_1 368 flush(LUPRI) 369 END IF 370C 371C *********************************************** 372C *** START LOOP OVER MACRO ITERATIONS ********** 373C *********************************************** 374C 375 ITMACN = 0 376 ITMICT = 0 377 ITBCK = 0 378 ITABS = 0 379C*****NON-STANDARD LOOP (MACRO ITERATIONS) 380 100 CONTINUE 381 IF (.NOT.ABSORP) ITABS = -1 382C 383 TIMMAC = SECOND() 384 IF (ITBCK .GT. 0) THEN 385 WRITE (LUW4,975) ITMAC,ITBCK 386 IF (LUPRI.NE.LUW4) WRITE (LUPRI,975) ITMAC,ITBCK 387 ELSE IF (ITABS .GT. 0) THEN 388 WRITE (LUW4,977) ITMAC 389 IF (LUPRI.NE.LUW4) WRITE (LUPRI,977) ITMAC 390 ELSE 391 ITMAC = ITMAC + 1 392 ITMACN = ITMACN + 1 393 WRITE (LUW4,970) ITMAC 394 IF (LUPRI.NE.LUW4) WRITE (LUPRI,970) ITMAC 395 END IF 396 970 FORMAT(//' --- MACRO ITERATION',I3,' ---' 397 * /' --------------------------') 398 975 FORMAT(//' --- MACRO ITERATION',I3,' BACKSTEP NO.',I3,' ---') 399 977 FORMAT(//' --- MACRO ITERATION',I3,' AFTER ABSORPTION ---') 400C 401 CALL FLSHFO(LUW4) 402 IF (LUPRI.NE.LUW4) CALL FLSHFO(LUPRI) 403 CALL DZERO(DINFO,LDINFO) 404 CALL IZERO(IINFO,LIINFO) 405 SHFLVL = D0 406C 407C First retrieve orbitals : 408C 409 CALL READMO (WRK(KCMO),9) 410C 411C *** Transform two-electron integrals 412C 413 TIMITR = SECOND() 414 IF (.NOT.FLAG(14)) THEN 415 JTRLVL = ITRLVL 416 IF (ITABS .EQ. 0) JTRLVL = KTRLVL 417 IF (FLAG(20) .OR. ITMACN .GE. MAXMAC) JTRLVL = 0 418C ... stop after gradient: only 1. order transform. needed 419 CALL TRACTL(JTRLVL,WRK(KCMO),WRK(KWRK10),LWRK10) 420C CALL TRACTL(ITRLVL,CMO,WRK,LFREE) 421 LTRLVL = JTRLVL 422 ELSE 423 LTRLVL = -999 424 END IF 425 TIMITR = SECOND() - TIMITR 426 IF (NWOPT.EQ.0 .OR. FLAG(34)) THEN 427C ... We do not need the integrals for optimization. 428C ... flag(34) = (int.transf. not needed in optimization) 429 FLAG(14) = .TRUE. 430 ELSE 431 FLAG(14) = .FALSE. 432 END IF 433C 434C *** Retrieve CI coefficients: 435C 436 CALL RDCREF(WRK(KCREF),.true.) 437C 438C --- Normalize CI vector (i.e. CREF) 439C 6-Dec-1984-hjaaj: SIRSAV changed so CREF normally normalized now 440C 441 IF (NCONF .GT. 1) THEN 442 C0NORM = DNRM2(NCONF,WRK(KCREF),1) 443 THREQL = SQRT(NCONF*THRLDP) 444 IF (ABS(C0NORM-D1) .GT. THREQL) THEN 445 IF (C0NORM.EQ.D0) THEN 446 WRITE (LUW4,1110) 447 IF (LUPRI.NE.LUW4) WRITE (LUPRI,1110) 448 WRITE (LUERR,1110) 449 CALL QUIT('ERROR (SIROPT) CI vector has zero norm') 450 END IF 451 C0NORM = D1/C0NORM 452 NWARN = NWARN + 1 453 WRITE (LUPRI,1100) C0NORM 454 CALL DSCAL(NCONF,C0NORM,WRK(KCREF),1) 455 END IF 456 END IF 457 1100 FORMAT(/' WARNING: Reference CI vector normalized, factor:', 458 & F15.10) 459 1110 FORMAT(/' Optimization control: FATAL ERROR', 460 * ' - zero norm of reference CI vector') 461C 462 IF (P6FLAG(32) .AND. NCONF.GT.1) THEN 463C use PRWF here ?? 890925-hjaaj 464 WRITE (LUPRI,2010) ITMAC,THRPWF 465 2010 FORMAT(/' Reference CI vector in iteration',I3, 466 * /' -----------------------------------', 467 * //' Cutoff for print:',F7.4,/) 468 CALL PRVEC(NCONF,WRK(KCREF),1,THRPWF,200,LUPRI) 469 CALL PRMGN(NCONF,WRK(KCREF),1,12,LUPRI) 470#if defined (VAR_OLDCODE) || defined (VAR_SECSEC) 471 DO 202 I = 1,NCONF 472 CREFI = WRK((KCREF-1)+I) 473 IF (ABS(CREFI) .GE. THRPWF) WRITE (LUPRI,2020) I,CREFI 474 202 CONTINUE 475 2010 FORMAT(/' Reference CI vector in iteration',I3, 476 * /' -----------------------------------', 477 * //' Cutoff for print:',F7.4, 478 * //' Configuration no. value', 479 * /' ----------------- -----') 480 2020 FORMAT(I16,F20.10) 481#endif 482 END IF 483C 484C *** Calculate gradient for CREF and some matrices/arrays needed later. 485C 486C Dec99 hjaaj: implemented dynamic Fock matrix screening for QCHF and DIRFCK 487C (based on code in sirdiis.F) : 488C set screening threshold for direct SCF gradient 489C GNORM(3) = grdnrm from prev.iter. (THRGRD in iter 1) 490C Nov07 hjaaj: we expect quadratic convergence and GNORM(3) is from 491C prev.iter., thus GNORM(3)**2 492C NORBT factor: statistical error is ca. sqrt(# elements) ~ NORBT 493C 494 IF (DIRFCK .AND. .NOT. USRSCR) THEN 495 IF (UPDATE) THEN 496 ! Hessian update gives max a factor 0.1, so 0.01 should be safe 497 GRDACC = MAX(THRGRD,0.01D0*GNORM(3)) 498 ELSE 499 ! quadratic convergence; IFTHRS .ge. 9 takes care of linear region 500 GRDACC = MAX(THRGRD,GNORM(3)**2) 501 END IF 502 GRDACC = GRDACC / NORBT 503 IFTHRS = MIN( IFTHSV, NINT(-LOG10(GRDACC)) + 1) 504 IFTHRS = MAX( 9, IFTHRS ) 505 IF (IPRI6 .GE. 0) WRITE (LUPRI,'(/A,I3,A)') 506 & ' Fock matrix screening setting for this macro iteration: 10^(' 507 & ,-IFTHRS,')' 508 END IF 509 IF (.NOT.UPDATE) THEN 510 CALL GRAD (WRK(KCREF),WRK(KCMO),WRK(KCINDX),WRK(KDV), 511 & WRK(KPV),WRK(KFC),WRK(KFV),WRK(KFQ),WRK(KFCAC), 512 & WRK(KH2AC),WRK(KGTOT),EMY,EACTIV, 513 & WRK(KWRKG),LWRKG) 514C CALL GRAD (CREF,CMO,INDXCI,DV,PV,FC,FV,FQ, 515C * FCAC,H2AC,G,EMY,EACTIV,WRK,LFREE) 516 ELSE 517 CALL UPDGRD(0,WRK(KCREF),WRK(KCMO),WRK(KDV),WRK(KPV), 518 & WRK(KGTOT),EMY,EACTIV,WRK(KFC),WRK(KFV),WRK(KFQ), 519 & DUMMY,WRK(KCINDX),WRK(KWRKG),LWRKG) 520C CALL UPDGRD(IGRTYP,REF1,CMO1,DV,PV,GRD,EMY, 521C * EACTIV,FC,FV,DIAOR,INDXCI,WRK,LFREE) 522 END IF 523C 524C hj-861130 solvent contribution, ESOLT is saved in /INFOPT/ 525C 526 IF (SOLVNT) THEN 527 CALL SOLGRD(WRK(KCREF),WRK(KCMO),WRK(KCINDX),WRK(KDV), 528 & WRK(KGTOT),ESOLT,WRK(KERLM),WRK(KWRKG),LWRKG) 529C CALL SOLGRD(CREF,CMO,INDXCI,DV,G,ESOLT,ERLM,WRK,LFREE) 530 ELSEIF (PCM) THEN 531 CALL PCMGRD(WRK(KCREF),WRK(KCMO),WRK(KCINDX),WRK(KDV), 532 % WRK(KGTOT),ESOLT,WRK(KWRKG),LWRKG) 533 ELSEIF (QMMM) THEN 534 CALL PEGRD(WRK(KCREF),WRK(KCMO),WRK(KCINDX),WRK(KDV), 535 & WRK(KGTOT),ESOLT,WRK(KWRKG),LWRKG) 536 ELSE IF (USE_PELIB()) THEN 537 CALL PELIB_IFC_GRAD(WRK(KCREF), WRK(KCMO), WRK(KCINDX), 538 & WRK(KDV), WRK(KGTOT), ESOLT, 539 & WRK(KWRKG), LWRKG) 540 ELSE 541 ESOLT = D0 542 END IF 543C---------------- 544C CBN+JK 03.01.06 545C---------------- 546 IF (QM3) THEN 547 CALL QM3GRAD(WRK(KCREF),WRK(KCMO), 548 & WRK(KCINDX),WRK(KGTOT), 549 & EQM3,WRK(KWRKG),LWRKG,IPRSOL) 550 ELSE 551 EQM3=D0 552 END IF 553C---------------- 554C CBN+JK 03.01.06 555C---------------- 556 GNRMSV = GNORM(3) 557 CALL GRDINF(GNORM,WRK(KGTOT)) 558C 559C HJ-860808 no absorption if orb.grd. .lt. fakabs*ci.grd. 560C 561 IF ( ITABS .EQ. 0 .AND. GNORM(2) .LT. FAKABS*GNORM(1) ) THEN 562 IF (IPRI6.GE.5) WRITE (LUPRI,'(/A/A/)') 563 * ' *** no orbital absorption this macro', 564 * ' orbital gradient significantly smaller than CI gradient.' 565 ITABS = -1 566 END IF 567C 568C---------------- 569 EMCSCF = POTNUC + EMY + EACTIV + ESOLT + EQM3 570C---------------- 571C 572 IF (ITABS .GT. 0) THEN 573 WRITE (LUW4,1005) WF_TYPE(1:len_WF_TYPE),EMCSCF,ITABS 574 IF (LUPRI.NE.LUW4) 575 & WRITE (LUPRI,1005) WF_TYPE(1:len_WF_TYPE),EMCSCF,ITABS 576 ELSE 577 WRITE (LUW4,1010) WF_TYPE(1:len_WF_TYPE),EMCSCF,ITMAC 578 IF (LUPRI.NE.LUW4) 579 & WRITE (LUPRI,1010) WF_TYPE(1:len_WF_TYPE),EMCSCF,ITMAC 580 END IF 581 IF (DOHFSRDFT.OR.DOMCSRDFT) THEN 582 EMYTMP = EMY - ESRDFTY 583 ELSE 584 EMYTMP = EMY 585 ENDIF 586 IF (P4FLAG(2)) THEN 587 WRITE (LUW4,1015) POTNUC,EMYTMP,EACTIV 588 IF (DOHFSRDFT.OR.DOMCSRDFT) WRITE(LUW4,1018) ESRDFTY 589 END IF 590chj IF (P6FLAG(1)) THEN 591 IF (LUPRI.NE.LUW4 .OR. .NOT.P4FLAG(2)) THEN 592 WRITE (LUPRI,1015) POTNUC,EMYTMP,EACTIV 593 IF (DOHFSRDFT.OR.DOMCSRDFT) WRITE(LUPRI,1018) ESRDFTY 594 END IF 595chj END IF 596 IF (SOLVNT) WRITE (LUW4,1016) ESOLT 597 IF (PCM) WRITE (LUW4,1019) ESOLT 598 IF (QM3) WRITE (LUW4,1017) EQM3 599 IF (QMMM) WRITE (LUW4,1017) ESOLT 600 IF (USE_PELIB()) WRITE (LUW4,1021) ESOLT 601 WRITE (LUW4,1020) GNORM(3),GNORM(1),GNORM(2) 602 IF (LUPRI.NE.LUW4) THEN 603 IF (SOLVNT) WRITE (LUPRI,1016) ESOLT 604 IF (PCM) WRITE (LUPRI,1019) ESOLT 605 IF (QM3) WRITE (LUPRI,1017) EQM3 606 IF (QMMM) WRITE (LUPRI,1017) ESOLT 607 IF (USE_PELIB()) WRITE(LUPRI,1021) ESOLT 608 WRITE (LUPRI,1020) GNORM(3),GNORM(1),GNORM(2) 609 END IF 610 1005 FORMAT(/' Total ',A,' energy',T27,':',F30.15,T54, 611 * '(after absorp',I4,')' ) 612 1010 FORMAT(/' Total ',A,' energy',T27,':',F30.15,T60, 613 * '(MACRO ',I4,')' ) 614 1015 FORMAT(/' - Nuclear repulsion :',F30.15, 615 * /' - Inactive energy :',F30.15, 616 * /' - Active energy :',F30.15) 617 1016 FORMAT( ' - Dielec. solvation ener.:',F30.15) 618 1017 FORMAT( ' - QM/MM solvation energy :',F30.15) 619 1021 FORMAT( ' - Embedding energy :',F30.15) 620 1018 FORMAT( ' - srDFT effective energy :',F30.15) 621 1019 FORMAT( ' - PCM solvation energy :',F30.15) 622 1020 FORMAT(/' Norm of total gradient :',F27.12, 623 * /' - of CI gradient :',F27.12, 624 * /' - of orbital gradient :',F27.12) 625 626 IF (lim_poppri .gt. 0) THEN 627 call sirpop('MCITER',WRK(KDV),WRK(KWRKG),LWRKG) 628 END IF 629 630 ! next line is debug print for impatient developers ;-) 631 ! print *,'MCITER',ITMAC,EMCSCF,GNORM(1:3) 632 633C 634C 635 IF (SOLVNT .AND. IPRI6 .GT. 1) THEN 636 WRITE (LUPRI,'(//A/)') 637 * ' Multipole analysis of dielectric solvation energy :' 638 JERLM = KERLM 639 DO 244 LVAL = 0,LSOLMX 640 NMVAL = 2*LVAL + 1 641 ERLTOT = DSUM(NMVAL,WRK(JERLM),1) 642 WRITE (LUPRI,'(A5,I3,T12,F12.8)') 643 * ' L =',LVAL,ERLTOT 644 IF (IPRI6 .GT. 5) THEN 645 WRITE(LUPRI,'(A11,5F12.8/,(T12,5F12.8))') 646 * ' Components', (WRK(JERLM+MVAL),MVAL=0,(NMVAL-1)) 647 WRITE(LUPRI,'()') 648 END IF 649 JERLM = JERLM + NMVAL 650 244 CONTINUE 651 END IF 652C 653 IF (P4FLAG(15) .AND. NCONF .GT. 1) THEN 654CHJ: in future version maybe CALL PRWF here (and above for CREF) 655 PTEST = MAX(1.D-10,GNORM(1)*THRCGR) 656 WRITE (LUW4,3000) ITMAC 657 CALL PRVEC(NCONF,WRK(KGTOT),1,-PTEST,200,LUW4) 658 CALL PRMGN(NCONF,WRK(KGTOT),1,12,LUW4) 659#if defined (VAR_OLDCODE) || defined (VAR_SECSEC) 660 WRITE (LUW4,3010) ITMAC, PTEST 661 DO 301 I = 1,NCONF 662 GCI = WRK((KGTOT-1)+I) 663 IF (ABS(GCI) .GE. PTEST) WRITE (LUW4,3011) I,GCI 664 301 CONTINUE 665 3010 FORMAT(//' Configuration gradient iteration no.',I3, 666 * /' ---------------------------------------' 667 * //' Cutoff for print:',1P,D10.2, 668 * //' Configuration no. value' 669 * /' ----------------- -----') 670 3011 FORMAT(I16,F20.10) 671#endif 672 END IF 673 IF (P4FLAG(14) .AND. NWOPT.GT.0) THEN 674 WRITE (LUW4,3020) ITMAC 675 IF (MPRI4.GT.100) THEN 676 PRFAC = 0.0D0 677 ELSE IF (MPRI4.GT.10) THEN 678 PRFAC = 0.1D0 679 ELSE 680 PRFAC = 0.2D0 681 END IF 682 CALL PRKAP (NWOPT,WRK(KGTOT+NCONF),PRFAC,LUW4) 683 END IF 684 3000 FORMAT(//' Configuration gradient iteration no.',I3, 685 * /' ---------------------------------------') 686 3020 FORMAT(//' Orbital gradient iteration no.',I3, 687 * /' ---------------------------------') 688C 689 CALL FLSHFO(LUW4) 690 IF (LUPRI.NE.LUW4) CALL FLSHFO(LUPRI) 691C 692C *** Test for convergence and trust region size 693C 694 IF (GNORM(3) .LE. THRGRD) THEN 695 ICONV = 1 696 ISTEP = -2 697 END IF 698C 699 IF (FLAG(20)) GO TO 700 700C 701 IF (ITMACN .GE. MAXMAC) ISTEP = -2 702C 703 IF ( CHKSTP ) THEN 704C istep = -2 : converged or stop after gradient, 705C get step info for analysis 706C istep = -1 : restart prediction available 707C istep = 0 : normal second-order step check 708 CALL SIRSTP(ISTEP,WRK(KCMO),WRK(KIBNDX),WRK(KCINDX), 709 * WRK(KWNEO),LWNEO) 710C CALL SIRSTP(ISTEP,CMO,IBNDX,INDXCI,WRK,LFREE) 711C 712 IF ( ICONV .NE. 1) THEN 713C Check if RATIO is OK for local UPDATE/DONR 714C NB! Only disable DONR if not flag(39) ".NR ALWAYS" 715 IF ( UPDATE .OR. (DONR .AND. .NOT. NRALW) ) THEN 716C if ( update ) then this is local update, 717C otherwise CHKSTP would be false 718 RATCHK = ABS(DINFO(6) - D1) 719 IF (RATCHK .GT. 0.02) THEN 720 IF (UPDATE) THEN 721 WRITE (LUPRI,'(/A,F10.5/A)') 722 & ' Local UPDATE disabled because 0.02 .lt. abs(ratio-1) =', 723 & RATCHK, ' We must therefore recalculate gradient.' 724 UPDATE = .FALSE. 725 DONR = NRALW 726Chj-sep99: next 3 lines work, the rest from go to 100 to 'ELSE' does not, why??? 727 CHKSTP = .FALSE. 728 itrlvl = itrlsv 729 go to 100 730Chj-sep99: previous 3 lines work, the rest from go to 100 to 'ELSE' does not, why??? 731#ifdef HJ_CHECK_LATER 732 IF (ITRLSV .GT. ITRLVL) THEN 733C need higher int. transf. level for NEO/NR 734 ITRLVL = ITRLSV 735 JTRLVL = ITRLVL 736 CALL TRACTL(JTRLVL,WRK(KCMO),WRK(KWRK10),LWRK10) 737 LTRLVL = JTRLVL 738 END IF 739 CALL GRAD (WRK(KCREF),WRK(KCMO),WRK(KCINDX),WRK(KDV), 740 * WRK(KPV),WRK(KFC),WRK(KFV),WRK(KFQ),WRK(KFCAC), 741 * WRK(KH2AC),WRK(KGTOT),EMY,EACTIV, 742 & WRK(KWRKG),LWRKG) 743 IF (SOLVNT) THEN 744 CALL SOLGRD(WRK(KCREF),WRK(KCMO),WRK(KCINDX), 745 * WRK(KDV),WRK(KGTOT),ESOLT, 746 * WRK(KERLM),WRK(KWRKG),LWRKG) 747 ELSEIF (PCM) THEN 748 CALL PCMGRD(WRK(KCREF),WRK(KCMO),WRK(KCINDX), 749 * WRK(KDV),WRK(KGTOT),ESOLT, 750 * WRK(KWRKG),LWRKG) 751 ELSEIF (QMMM) THEN 752 CALL PEGRD(WRK(KCREF),WRK(KCMO),WRK(KCINDX), 753 * WRK(KDV),WRK(KGTOT),ESOLT, 754 * WRK(KWRKG),LWRKG) 755 ELSE IF (USE_PELIB()) THEN 756 CALL PELIB_IFC_GRAD(WRK(KCREF), WRK(KCMO), 757 & WRK(KCINDX), WRK(KDV), 758 & WRK(KGTOT), ESOLT, 759 & WRK(KWRKG), LWRKG) 760 END IF 761#endif 762 ELSE 763 WRITE (LUPRI,'(/A,F10.5)') 764 & ' Local NR disabled because 0.02 .lt. abs(ratio-1) =', 765 & RATCHK 766 DONR = .FALSE. 767 END IF 768 END IF 769 END IF 770 END IF 771 ELSE 772 EMCOLD = EMCSCF 773 ISTEP = 0 774 END IF 775 DINFO(9) = RTRUST 776 CALL FLSHFO(LUW4) 777 IF (LUPRI.NE.LUW4) CALL FLSHFO(LUPRI) 778C --- Exit if converged --- 779 IF (ICONV .EQ. 1) GO TO 700 780C --- Exit if max macro iter reached --- 781 IF (ITMACN .GE. MAXMAC) GO TO 700 782C 783 IF (ISTEP.GT.0) THEN 784 ISTEP = 0 785 UPDATE = .FALSE. 786 DONR = NRALW 787 ITRLVL = ITRLSV 788 ITBCK = ITBCK + 1 789 IF (ITBCK.LE.MAXBCK) THEN 790 WRITE(LUW4,9000) 791 IF (LUPRI.NE.LUW4) WRITE(LUPRI,9000) 792 ELSE 793 WRITE (LUERR,9100) ITBCK-1 794 WRITE (LUW4 ,9100) ITBCK-1 795 IF (LUPRI.NE.LUW4) WRITE (LUPRI,9100) ITBCK-1 796 END IF 797 GO TO 700 798 ELSE IF (ISTEP.LT.0) THEN 799C close to convergence, no step control performed 800 IF (ITMACN .GT. 1 .AND. GNORM(3) .GT. GNRMSV) THEN 801C ... hjaaj: no test if restart (ITMACN .eq. 1) 802 WRITE (LUW4,'(/A/A)') 803 & ' ERROR EXIT: gradient norm increasing in local region.', 804 & ' Probable cause: numerical round-off errors'// 805 & ' or too many integrals neglected.' 806 IF (LUPRI .NE. LUW4) WRITE (LUPRI,'(/A/A)') 807 & ' ERROR EXIT: gradient norm increasing in local region.', 808 & ' Probable cause: numerical round-off errors'// 809 & ' or too many integrals neglected.' 810c GO TO 750 811 write(lupri,*) 'Well, let us try anyway ...' 812 END IF 813 ITBCK = 0 814 ELSE 815 ITBCK = 0 816 END IF 817 9000 FORMAT(//' Optimization control: Step was too large, we backstep') 818 9100 FORMAT(//' Optimization control: ERROR EXIT NO CONVERGENCE' 819 * //' Step was too large and maximum number of backsteps,', 820 * I2,', reached'/) 821C 822C Check ratio for geometry optimization (ABACUS interface) 823C 824 IF (GEOWLK .AND. JWSTEP .EQ. 1) THEN 825 IF (GNORM(3) .LE. THRGEO) THEN 826 DEPSAV = DEPRED 827 EMCOLD = EMCGEO 828 DEPRED = DEPGEO 829 CALL SIRSTP(JWSTEP,VDUMMY,DUMMY,DUMMY,DUMMY,1) 830 EMCOLD = EMCSCF 831 DEPRED = DEPSAV 832 IF (WLKREJ) GO TO 700 833 END IF 834 END IF 835C 836C ********** 837C We have not converged (yet), 838C 839C 840C if (ABSORP) then call SIRORB for orbital optimization. 841C 842 CHKSTP = .TRUE. 843 IF (ITABS .EQ. 0) THEN 844 CHKSTP = .FALSE. 845 CALL SIRORB(MAXAPM,ITABS,WRK(KCMO),WRK(KGTOT+NCONF),WRK(KDV), 846 & WRK(KPV),WRK(KFC),WRK(KFV),WRK(KFQ), WRK(KWOCTL),LWOCTL) 847C CALL SIRORB(MAXITO,NITORB,CMO,GORB,DV,PV,FC,FV,FQ,WRK,LFREE) 848 GO TO 700 849C ... change 870607/hjaaj, was GO TO 100 850C use GO TO 700 to get timing statistics. 851 ELSE IF (ABSORP) THEN 852 ITABS = 0 853 END IF 854C 855C If (UPDATE) then call SIRUPD to do UPDATE optimization, 856C else call SIRNEO to find next step vector (using the NEO 857C optimization technique). 858C 859 IF (UPDATE) THEN 860 CHKSTP = .FALSE. 861 ITRLVL = 0 862 NWOPH = NWOPT 863 NVARH = NVAR 864 CALL DCOPY(NCONF,WRK(KCREF),1,WRK(KREFU),1) 865 CALL DCOPY(NVAR,WRK(KGTOT),1,WRK(KGRDU),1) 866 CALL SIRUPD(ICONV,MAXUIT,WRK(KREFU),WRK(KGRDU),WRK(KCMO), 867 & WRK(KCINDX),WRK(KWU),LWU) 868C CALL SIRUPD(ICONV,MAXIT,REF1,GRD,CMO1,INDXCI,WRK,LFREE) 869 LTRLVL = ITRLVL 870 ITRLVL = ITRLSV 871 NWOPH = NWOPS 872 NVARH = NVARS 873 NREDL = 0 874 UPDATE = .FALSE. 875 IF (ICONV .GT. 0) THEN 876 IF (FLAG(25)) THEN 877 WRITE (LUPRI,'(//A/)') 878 & ' --- after UPDATE: recalculating gradient for SIRIFC' 879 ICONV = 0 880 UPDATE = .FALSE. 881 IF (.NOT. RHF_CALC) FLAG(14) = .TRUE. 882 GO TO 700 883 END IF 884 GO TO 800 885 ELSE 886 IF (ISTATE .EQ. 1 .AND. NROOTS .EQ. 1) THEN 887 GO TO 700 888C ... change 870607/hjaaj, was GO TO 100 889C use GO TO 700 to get timing statistics. 890 ELSE 891 WRITE (LUW4,'(/A)') 892 & ' ERROR EXIT: NEO cannot proceed after UPDATE for NROOTS .gt. 1' 893 IF (LUPRI .NE. LUW4) WRITE (LUPRI,'(/A)') 894 & ' ERROR EXIT: NEO cannot proceed after UPDATE for NROOTS .gt. 1' 895 GO TO 750 896C ... too few start vectors on LUIT1 897 END IF 898 END IF 899 END IF 900C 901C 902C ************************************************************ 903C second order step 904C ************************************************************ 905C 906C 907 IF (DONR) THEN 908 CALL SIRNR (WRK(KCREF),WRK(KGTOT),WRK(KCMO),WRK(KCINDX), 909 * WRK(KDV),WRK(KPV), 910 * WRK(KFC),WRK(KFV),WRK(KFCAC), 911 * WRK(KH2AC),EACTIV,WRK(KIBNDX), 912 & WRK(KWNEO),LWNEO) 913C CALL SIRNR (CREF,G,CMO,INDXCI,DV,PV, 914C * FC,FV,FCAC,H2AC,EACTIV,IBNDX,WRK,LFREE) 915C 916C Absorption should not be used after SIRNR 917C (because SIRNR can not check if root flip has occurred 918C for excited states, and anyway SIRNR is normally only 919C used in the local region). 920C 921 ABSORP = .FALSE. 922 FLAG(51) = .FALSE. 923 FLAG(52) = .FALSE. 924 FLAG(53) = .FALSE. 925 ELSE 926 CALL SIRNEO(WRK(KCREF),WRK(KGTOT),WRK(KCMO),WRK(KCINDX), 927 * WRK(KDV),WRK(KPV),WRK(KFC),WRK(KFV),WRK(KFCAC), 928 * WRK(KH2AC),EACTIV,WRK(KIBNDX), 929 & WRK(KWNEO),LWNEO) 930 931C CALL SIRNEO(CREF,G,CMO,INDXCI,DV,PV, 932C * FC,FV,FCAC,H2AC,EACTIV,IBNDX,WRK,LFREE) 933C 934C Test if SIRNEO has disabled absorption 935C (because of root flipping for excited state optimization) 936C 937 IF (ABSORP) ABSORP = ( FLAG(51) .OR. FLAG(52) .OR. FLAG(53) ) 938 END IF 939C 940C Check if local region where absorption normally is 941C switched off (unless FLAG(54) has been set in input) 942C 943 IF (ABSORP .AND. .NOT.FLAG(54)) THEN 944 LOCAL = (STPLEN .LT. STPABS*RTRUST) .AND. ITBCK .EQ. 0 945 IF (LOCAL) ITABS = -1 946C ... no absorption next iteration 947 LOCAL = (LOCAL .AND. STPC .LT. D2*STPO) 948C ... only switch absorption completely off if also 949C STPC .lt. 2*STPO (else we may get a large orbital step 950C again). 951 IF (LOCAL .OR. ITMACN .GE. MAXABS) THEN 952 WRITE (LUPRI,'(//A)') ' *** Orbital absorption switched off' 953 IF (LOCAL) THEN 954 WRITE (LUPRI,'(A/)') 955 * ' because local region has been reached.' 956 ELSE 957 WRITE (LUPRI,'(A/)') 958 * ' because max macro iterations w/ absorption reached.' 959 END IF 960 ABSORP = .FALSE. 961 FLAG(51) = .FALSE. 962 FLAG(52) = .FALSE. 963 FLAG(53) = .FALSE. 964 END IF 965 END IF 966C 967C Test if switch to update method or switch to NR linear equations 968C (condition: local region, defined as step .le. STPLCL*RTRUST) 969C 970 IF (.NOT. DOMCSRDFT) THEN 971 ! Oct 2019 hjaaj: DONR is not working for excited state MCsrDFT ... 972 IF (.NOT.FLAG(37) .OR. .NOT.FLAG(38) .OR. .NOT.NRALW) THEN 973C If (not noupdate or not neo_always or not nr_always) then 974 LOCAL = (STPLEN .LT. STPLCL*RTRUST .AND. STPC .LT. STPO) 975 I = ISTATE + 1 - IINFO(13) 976 SHFLVL = DINFO(20+I) 977 IF (IPRSTAT .GT. 5) WRITE (LUSTAT,'(/A,F20.10)') 978 * ' Level shift used for UPDATE/SIRNR test:',SHFLVL 979 LOCAL = LOCAL .AND. (ABS(SHFLVL) .LT. SHFLCL) 980 LOCAL = LOCAL .AND. ITBCK .EQ. 0 981 IF ( LOCAL ) THEN 982 IF (.NOT. FLAG(37)) THEN 983 WRITE (LUPRI,'(/A)') ' Local region, UPDATE activated' 984 UPDATE = .TRUE. 985 DONR = .FALSE. 986 ITRLVL = 0 987 ELSE IF (.NOT. FLAG(38) .AND. ISTATE .GT. 1) THEN 988C ... no reason to switch to NR for ground state 989 WRITE (LUPRI,'(/A)') ' Local region, DONR activated' 990 DONR = .TRUE. 991 FLAG(39) = .TRUE. 992 END IF 993 END IF 994 END IF 995 END IF ! IF (.NOT. DOMCSRDFT) THEN 996C 997C save information for final summary print-out 998C 999 700 CONTINUE 1000 TIMMAC = SECOND() - TIMMAC 1001 IF (IPRI6 .GT. 0) WRITE (LUPRI,'(/A,F15.2,A)') 1002 * ' --- Time used in this macro iteration :',TIMMAC,' seconds.' 1003C 1004 IINFO(1) = ITMAC 1005 DINFO(1) = EMY 1006 DINFO(2) = EACTIV 1007 DINFO(3) = EMCSCF 1008 DINFO(16) = TIMMAC 1009 DINFO(18) = TIMITR 1010 WRITE (LUINF) DINFO,IINFO 1011C 1012 IF (IPRSTAT .GT. 1) THEN 1013 WRITE (LUSTAT,'(//A/A/)') 1014 * ' Iteration no., total time, integral transformation time,', 1015 * ' micro it. time, and lin.trn. time' 1016 WRITE (LUSTAT,'(I5,4F10.2)') 1017 * ITMAC,TIMMAC,TIMITR,DINFO(17),DINFO(19) 1018 END IF 1019C ( DINFO(17) = TIMMIC and DINFO(19) = TIMLIN ) 1020C 1021C CI case special output 1022C 1023 IF (NWOPT.EQ.0 .AND. NCONF.GT.1) THEN 1024 WRITE (LUW4,7011) DEPRED,(EMCSCF+DEPRED-POTNUC),EMCSCF+DEPRED 1025 IF (LUPRI.NE.LUW4) 1026 * WRITE (LUPRI,7011) DEPRED,(EMCSCF+DEPRED-POTNUC),EMCSCF+DEPRED 1027 END IF 1028 7011 FORMAT(/' (SIROPT) CI energy lowering: ',F30.15, 1029 * /T11,'CI electronic energy:',F30.15, 1030 * /T11,'CI total energy: ',F30.15) 1031C 1032 IF (ICONV .GT. 0) GO TO 800 1033 IF (FLAG(20)) GO TO 890 1034 IF (ITBCK .GT. 0 .AND. ITBCK.LE.MAXBCK) GO TO 100 1035 IF (ITBCK .EQ. 0 .AND. ITMACN.LT.MAXMAC .AND. 1036 * ITMICT.LT.MAXMIC) GO TO 100 1037C 1038C *********************************** 1039C *** END OF MACRO ITERATION LOOP *** 1040C *********************************** 1041C 1042C We did not converge... 1043C 1044 NWARN = NWARN + 1 1045 CALL GETTIM(T_opt, W_opt) 1046 T_opt = T_opt - T_start 1047 W_opt = W_opt - W_start 1048 WRITE (LUSTAT,7050) WF_TYPE(1:len_WF_TYPE), 1049 & ITMACN,MAXMAC,ITMICT,MAXMIC,ITBCK,MAXBCK,T_opt,W_opt 1050 WRITE (LUW4,7050) WF_TYPE(1:len_WF_TYPE), 1051 & ITMACN,MAXMAC,ITMICT,MAXMIC,ITBCK,MAXBCK,T_opt,W_opt 1052 IF (LUPRI.NE.LUW4) WRITE (LUPRI,7050) WF_TYPE(1:len_WF_TYPE), 1053 & ITMACN,MAXMAC,ITMICT,MAXMIC,ITBCK,MAXBCK,T_opt,W_opt 1054 7050 FORMAT( 1055 * /' *** Optimization control WARNING: ',A, 1056 & ' not converged ***', 1057 * /5X, 'Maximum number of iterations or backsteps reached:', 1058 * /5X, 'Number of macro iterations used',T45,I5, 1059 * /5X, ' Maximum',T45,I5, 1060 * /5X, 'Number of micro iterations used',T45,I5, 1061 * /5X, ' Maximum',T45,I5, 1062 * /5X, 'Number of back steps this macro',T45,I5, 1063 * /5X, ' Maximum',T45,I5, 1064 * /5X, 'Total number of CPU seconds used',F13.2, 1065 * /5X, 'Total WALL time used (seconds) ',F13.2) 1066 GO TO 890 1067 750 CONTINUE 1068 NWARN = NWARN + 1 1069 CALL GETTIM(T_opt, W_opt) 1070 T_opt = T_opt - T_start 1071 W_opt = W_opt - W_start 1072 WRITE (LUSTAT,7051) 1073 & WF_TYPE(1:len_WF_TYPE),ITMACN,ITMICT,T_opt,W_opt 1074 WRITE (LUW4 ,7051) 1075 & WF_TYPE(1:len_WF_TYPE),ITMACN,ITMICT,T_opt,W_opt 1076 IF (LUPRI.NE.LUW4) WRITE (LUPRI ,7051) 1077 & WF_TYPE(1:len_WF_TYPE),ITMACN,ITMICT,T_opt,W_opt 1078 7051 FORMAT( 1079 * /' *** Optimization control WARNING: ',A, 1080 & ' not converged ***', 1081 * /5X, 'Number of macro iterations used',T45,I5, 1082 * /5X, 'Number of micro iterations used',T45,I5, 1083 * /5X, 'Total number of CPU seconds used',F13.2, 1084 * /5X, 'Total WALL time used (seconds) ',F13.2) 1085 GO TO 890 1086C 1087C *********************************************************** 1088C 1089C We have converged... 1090C 1091 800 CONTINUE 1092 CALL GETTIM(T_opt, W_opt) 1093 T_opt = T_opt - T_start 1094 W_opt = W_opt - W_start 1095 WRITE (LUW4 ,8000) WF_TYPE(1:len_WF_TYPE),ITMAC,ITMICT,T_opt,W_opt 1096 IF (LUPRI.NE.LUW4) 1097 &WRITE (LUPRI,8000) WF_TYPE(1:len_WF_TYPE),ITMAC,ITMICT,T_opt,W_opt 1098 8000 FORMAT(/' *** Optimization control: ',A,' converged ***', 1099 * /5X, 'Number of macro iterations used',T45,I5, 1100 * /5X, 'Number of micro iterations used',T45,I5, 1101 * /5X, 'Total number of CPU seconds used',F13.2, 1102 * /5X, 'Total WALL time used (seconds) ',F13.2) 1103C 1104 IF (SOLVNT) THEN 1105 KFCSOL = KWRKG 1106 KWRKG = KFCSOL + NNORBT 1107 LWRKG = LFREE - (KWRKG-1) 1108 CALL SOLFCK(DUMMY,DUMMY,WRK(KFCSOL),WRK(KCMO), 1109 & WRK(KERLM),.TRUE.,DUMSOL, 1110 & WRK(KWRKG),LWRKG,IPRSOL) 1111 CALL DAXPY(NNORBT,D1,WRK(KFC),1,WRK(KFCSOL),1) 1112 ELSE IF (PCM) THEN 1113 KFCSOL = KWRKG 1114 KWRKG = KFCSOL + NNORBT 1115 LWRKG = LFREE - (KWRKG-1) 1116 CALL PCMFCK(DUMMY,DUMMY,WRK(KFCSOL),WRK(KCMO), 1117 & .TRUE.,DUMSOL, WRK(KWRKG),LWRKG,IPRSOL) 1118 CALL DSCAL(NNORBT,-1.0D0,WRK(KFCSOL),1) 1119 CALL DAXPY(NNORBT,D1,WRK(KFC),1,WRK(KFCSOL),1) 1120C CBN+JK 03.01.06 1121 ELSE IF (QM3) THEN 1122 KFCSOL = KWRKG 1123 KWRKG = KFCSOL + NNORBT 1124 LWRKG = LFREE - (KWRKG-1) 1125 CALL QM3FCKMO(WRK(KCMO),WRK(KFCSOL),WRK(KWRKG), 1126 & LWRKG,IPRSOL) 1127 CALL DAXPY(NNORBT,D1,WRK(KFC),1,WRK(KFCSOL),1) 1128 ELSE IF (QMMM) THEN 1129 KFCSOL = KWRKG 1130 KWRKG = KFCSOL + NNORBT 1131 LWRKG = LFREE - (KWRKG-1) 1132 ! NB: only works without symmetry 1133 CALL PEFCMO(WRK(KCMO),WRK(KFCSOL),WRK(KDV), 1134 & WRK(KWRKG),LWRKG,IPRINT) 1135 CALL DAXPY(NNORBT,D1,WRK(KFC),1,WRK(KFCSOL),1) 1136 ELSE IF (USE_PELIB()) THEN 1137 KFCSOL = KWRKG 1138 KENR = KFCSOL + NNORBX 1139 KDCAO = KENR + 1 1140 KDVAO = KDCAO + N2BASX 1141 KFDTAO = KDVAO + N2BASX 1142 KFCKAO = KFDTAO + NNBASX 1143 KWRKG = KFCKAO + NNBASX 1144 LWRKG = LFREE - KWRKG + 1 1145 IF (LWRKG > LFREE) CALL ERRWRK('SIROPT-PELIB', -KWRKG, LFREE) 1146 CALL FCKDEN((NISHT>0), (NASHT>0), WRK(KDCAO), WRK(KDVAO), 1147 & WRK(KCMO), WRK(KDV), WRK(KWRKG), LWRKG) 1148 IF (NISHT == 0) THEN 1149 CALL DZERO(WRK(KDCAO), N2BASX) 1150 END IF 1151 IF (NASHT > 0) THEN 1152 CALL DAXPY(N2BASX, 1.0D0, WRK(KDVAO), 1, WRK(KDCAO), 1) 1153 END IF 1154 CALL DGEFSP(NBAST, WRK(KDCAO), WRK(KFDTAO)) 1155 CALL PELIB_IFC_FOCK(WRK(KFDTAO), WRK(KFCKAO), WRK(KENR)) 1156 CALL UTHU(WRK(KFCKAO), WRK(KFCSOL), WRK(KCMO), WRK(KWRKG), 1157 & NBAST, NORBT) 1158 CALL DAXPY(NNORBX, 1.0D0, WRK(KFC), 1, WRK(KFCSOL), 1) 1159 CALL PELIB_IFC_ENERGY(WRK(KFDTAO)) 1160 ELSE 1161 KFCSOL = KFC 1162 END IF 1163 1164 IF (NASHT .EQ. 0) THEN 1165C hjaaj dec 2002: if extended to cover nasht.eq.1 as in sirdiis.F 1166C then remember to add FCSOL and FV to FD for RHFENR call 1167 WRITE (LUW4,'(//A)') ' *** SCF orbital energy analysis ***' 1168 IF (SOLVNT .OR. QM3) 1169 & WRITE (LUW4,'(A)') ' (incl. solvent contribution)' 1170 IF (USE_PELIB()) THEN 1171 WRITE(LUW4,'(A)') 1172 & ' (incl. contribution from polarizable embedding potential)' 1173 END IF 1174 CALL RHFENR(IPRI4,LUW4,WRK(KFCSOL), WRK(KWRKG),LWRKG) 1175 IF (LUPRI .NE. LUW4) THEN 1176 WRITE (LUPRI,'(//A)') ' *** SCF orbital energy analysis ***' 1177 IF (SOLVNT .OR. QM3) 1178 & WRITE (LUPRI,'(A)') ' (incl. solvent contribution)' 1179 IF (USE_PELIB()) THEN 1180 WRITE(LUW4,'(A)') 1181 & ' (incl. contribution from polarizable embedding potential)' 1182 END IF 1183 CALL RHFENR(IPRI6,LUPRI,WRK(KFCSOL), WRK(KWRKG),LWRKG) 1184 END IF 1185C CALL RHFENR(IPRINT,LUPRI,FD,SCRA,LSCRA) 1186 END IF 1187C 1188C 1189C *** LOCALIZATION 1190C 1191 IF (BOYORB.OR.PIPORB) CALL SIRLOC(WRK(KCMO),WRK(KWRKG),LWRKG) 1192C 1193 IF (LMULBS .OR. R12CAL) THEN 1194C Construct orthonormal auxiliary basis sets (WK/UniKA/04-11-2002). 1195 LWRKG = LFREE - (KWRKG-1) 1196 CALL TIMER('START ',TIMSTR,TIMEND) 1197 CALL R12AUX(WRK(KWRKG),LWRKG) 1198 CALL TIMER('R12AUX',TIMSTR,TIMEND) 1199 END IF 1200C 1201C If (FLAG(25) .OR. INERSI) write LUSIFC 1202C modified hjaaj-aug99: only if this is final MO level 1203C FLAG(25) is true for .ABACUS (geometry optimization) 1204C or .RESPONSE (response calculation) 1205C or .INTERFACE (request for LUSIFC) 1206C INERSI is true if this is initial state calculation 1207C in a solvent calculation with inertial polarization 1208C 1209 LSTLVL = 1 1210 IF (FLAG(21) .AND. DOMC) LSTLVL = 0 1211C ... this is RHF which will be followed by MCSCF 1212 IF (LSTLVL .EQ. 1 .AND. (FLAG(25) .OR. INERSI)) THEN 1213 IF (GEOWLK) THEN 1214 IF (.NOT. OPTNEW) WRITE (LUW4,'(//A/)') 1215 & ' Check ratio for geometry walk with converged MC :' 1216 JWSTEP = 1 1217 EMCOLD = EMCGEO 1218 DEPRED = DEPGEO 1219 CALL SIRSTP(JWSTEP,VDUMMY,DUMMY,DUMMY,DUMMY,1) 1220 ELSE 1221 WLKREJ = .FALSE. 1222 END IF 1223C 1224Chj if walk is rejected: do not update information on SIRIFC 1225 IF (.NOT.WLKREJ) 1226 & CALL WR_SIRIFC(WRK(KCMO),WRK(KDV),WRK(KPV), 1227 & WRK(KFCSOL),WRK(KFV),WRK(KFQ), 1228 & WRK(KCREF),WRK(KFCAC),WRK(KH2AC), 1229 & WRK(KWRKG),LWRKG,WRK(KGTOT+NCONF),WRK(KERLM), 1230 & .TRUE.,WRK(KCINDX)) 1231C CALL WR_SIRIFC(CMO,DV,PV,FC,FV,FQ,CREF,FCAC,H2AC,WRK,LFREE, 1232C & GORB,ERLM,ORBHES,XINDX) 1233C ... note: GORB will in a normal macro iteration be 1234C overwritten in SIRNEO or SIRNR, but none of 1235C these routines have been called after GRAD 1236C when the wave function is converged. 1237C 1238 END IF 1239C 1240C Transform integrals for abacus or response or 1241C for other programs, if requested. 1242C modified hjaaj-aug99: only if this is final wf level 1243C 1244 LSTLVL = 1 1245 IF (FLAG(21) .AND. (DOMP2 .OR. DOCI .OR. DOMC)) LSTLVL = 0 1246C ... if (true) this QCHF will be followed by higher wf level 1247 IF (LSTLVL .EQ. 1 .AND. ABS(ITRFIN) .LE. 10) THEN 1248 CALL FLSHFO(LUW4) 1249 IF (LUPRI.NE.LUW4) CALL FLSHFO(LUPRI) 1250 IF (IPRI6 .GE. 5) WRITE (LUPRI,'(/2A,I3,A)') 1251 * ' --- Transforming integrals for Abacus/response with', 1252 * ' ITRLVL =',ITRFIN,' ---' 1253 IF (ABS(ITRFIN) .EQ. LTRLVL) THEN 1254 IF (IPRI6 .GE. 5) WRITE (LUPRI,'(/A)') 1255 * ' Transformation skipped, was already performed.' 1256 ELSE 1257 JTRLVL = ITRFIN 1258 CALL TRACTL(JTRLVL,WRK(KCMO),WRK(KWRK10),LWRK10) 1259 END IF 1260 END IF 1261C 1262C *********************************************************** 1263C 1264 890 CONTINUE 1265 1266 !> store 1-particle density matrix in AO-basis on file 1267!#define AO_DENSITIES_ON_FILE 1268#ifdef AO_DENSITIES_ON_FILE 1269 call rdcref(wrk(kcref)) 1270 call makdv(wrk(kcref),wrk(kdv),wrk(kcindx),wrk(kwrkg),lwrkg) 1271 call write_aodens(wrk(kdv),nisht,nasht,nbast,nnashx, 1272 & ncmot,nsym,nbas,ibas,lupri,.false.,"MC") 1273#undef AO_DENSITIES_ON_FILE 1274#endif 1275 1276C restore control variables 1277 FLAG(39) = NRALW 1278 IFTHRS = IFTHSV 1279 CALL FLSHFO(LUW4) 1280 IF (LUPRI.NE.LUW4) CALL FLSHFO(LUPRI) 1281C 1282C write timing information for optimization to LUSTAT 1283C 1284 IF (IPRSTAT.GT.0) CALL TIMOPT('END',LUSTAT,WRK,LFREE) 1285C *** end of subroutine SIROPT 1286 CALL QEXIT('SIROPT') 1287 RETURN 1288 END 1289C /* Deck rhfenr */ 1290 SUBROUTINE RHFENR(IPRINT,LUPRI,FD,SCRA,LSCRA) 1291C 1292C 27-Oct-1990 Hans Joergen Aa. Jensen (revised 17-Dec-2002 hjaaj) 1293C 1294C Purpose: 1295C Check if orbitals are canonical Hartree-Fock orbitals 1296C and print orbital energies 1297C 1298C Input: 1299C FD; the total inactive + active Fock matrix 1300C 1301C Scratch: 1302C SCRA; general scratch area 1303C 1304#include "implicit.h" 1305 DIMENSION FD(*),SCRA(*) 1306C 1307C 1308 PARAMETER (D0 = 0.0D0, EMYCNV = 1.D-4, 1309 & DBIG = 1.D+12, GAPMIN = 0.1D0) 1310#include "dummy.h" 1311C 1312C Used from common blocks: 1313C pgroup : REP 1314C INFORB : NISH(*), NISHT, ... 1315Chj-sep99: do not use NRHF(*) and NRHFT, because they are not 1316C correct in DIIS if AUTOCC and occupations have changed. 1317C SCBRHF : IOPRHF,NFRRHF(*) 1318C INFIND : IROW(*),ISW(*) 1319C 1320#include "maxash.h" 1321#include "maxorb.h" 1322#include "pgroup.h" 1323#include "molde.h" 1324#include "infinp.h" 1325#include "inforb.h" 1326#include "scbrhf.h" 1327#include "infind.h" 1328#include "mxcent.h" 1329C 1330 CALL QENTER('RHFENR') 1331C 1332 IF (NASHT .GT. 0) THEN 1333C hjaaj Dec 2002: 1334C allow printing of orbital energies for open-shell systems. 1335 WRITE (LUPRI,'(/A/A/A/A)') 1336 &' Orbital energy analysis for an open-shell system.', 1337 &' Orbital energies are not uniquely defined'// 1338 &' for open-shell systems', 1339 &' here is used block diagonalization of'// 1340 & ' the FD=FC+FV Fock matrix.', 1341 &" NOTE that Koopmans' theorem is not fulfilled for this case." 1342C WRITE (LUPRI,'(/A)') 1343C * ' Orbital energy analysis skipped for open-shell systems' 1344C GO TO 9999 1345C hjaaj TODO: EKT orbital energies may be defined later. 1346 END IF 1347 IF (IPRINT .LT. 10 .AND. NSSHT .GT. 20) 1348 & WRITE (LUPRI,'(/A)') ' Only the 20'// 1349 & ' lowest virtual orbital energies printed in each symmetry.' 1350 IF (IPRINT .GE. 2) WRITE (LUPRI,'(/A,I5/A,8I5)') 1351 & ' Number of electrons :',2*NISHT, 1352 & ' Orbital occupations :',(NISH(I),I=1,NSYM) 1353C 1354 KORBEN = 1 1355 LNEED = KORBEN + NORBT 1356 IF (LNEED .GT. LSCRA) CALL ERRWRK('RHFENR',LNEED,LSCRA) 1357 CALL DZERO(SCRA,NORBT) 1358C 1359C Print orbital energies. 1360C Check if orbitals are canonical orbitals 1361C 1362 IF (SUPSYM) THEN 1363 IF (DODFT) THEN 1364 WRITE (LUPRI,1841) 1365 ELSE 1366 WRITE (LUPRI,1741) 1367 END IF 1368 ELSE 1369 IF (DODFT) THEN 1370 WRITE (LUPRI,1840) 1371 ELSE 1372 WRITE (LUPRI,1740) 1373 END IF 1374 END IF 1375 VALCON = D0 1376 EHOMO = -DBIG 1377 IHOMO = 0 1378 ESOMO = -DBIG 1379 ISOMO = 0 1380 ELUMO = DBIG 1381 ILUMO = 0 1382 DO 200 ISYM = 1,NSYM 1383 NORBI = NORB(ISYM) 1384 IF (NORBI.EQ.0) GO TO 200 1385 NFRHFI= NFRRHF(ISYM) 1386 IORBI = IORB(ISYM) 1387 ISSYM = IIORB(ISYM) 1388 EHOMOI= EHOMO 1389 ESOMOI= ESOMO 1390 ELUMOI= ELUMO 1391 DO 300 I = 1,NORBI 1392 IJOFF = ISSYM + IROW(I) 1393 SCRA(IORBI+I) = FD(IJOFF+I) 1394 IF (I .GT. NFRHFI) THEN 1395C ... not frozen 1396 IW = ISW(IORBI+I) 1397 IF (IW.LE.NISHT) THEN 1398C ... closed shell (inactive) orbital 1399 EHOMOI = MAX(EHOMOI,SCRA(IORBI+I)) 1400 DO J = NFRHFI+1,I-1 1401 VALCON = MAX(VALCON,ABS(FD(IJOFF+J))) 1402 END DO 1403 ELSE IF (IW.LE.NOCCT) THEN 1404C ... open shell (active) orbital 1405 ESOMOI = MAX(ESOMOI,SCRA(IORBI+I)) 1406C No 'canonical' test for the open shell active orbitals 1407C because the Fock matrix is not zero for these 1408 ELSE 1409C ... virtual (secondary) orbital 1410 ELUMOI = MIN(ELUMOI,SCRA(IORBI+I)) 1411 DO J = NFRHFI+1,I-1 1412 VALCON = MAX(VALCON,ABS(FD(IJOFF+J))) 1413 END DO 1414C ... skipping the open shell active orbitals 1415C because the Fock matrix is not zero for these 1416 DO J = NOCCT+1,I-1 1417 VALCON = MAX(VALCON,ABS(FD(IJOFF+J))) 1418 END DO 1419 END IF 1420 END IF 1421 300 CONTINUE 1422 IF (EHOMOI .GT. EHOMO) THEN 1423 EHOMO = EHOMOI 1424 IHOMO = ISYM 1425 END IF 1426 IF (ESOMOI .GT. ESOMO) THEN 1427 ESOMO = ESOMOI 1428 ISOMO = ISYM 1429 END IF 1430 IF (ELUMOI .LT. ELUMO) THEN 1431 ELUMO = ELUMOI 1432 ILUMO = ISYM 1433 END IF 1434 IF (IPRINT .GE. 10) THEN 1435 NLAST = NORBI 1436 ELSE 1437 NLAST = MIN(NORBI,NOCC(ISYM)+20) 1438 END IF 1439 IF (NFRHFI .GT. 0) WRITE (LUPRI,'(A,I3,A)') 1440 &' Note: the first',NFRHFI,' orbitals in next symmetry are frozen' 1441 IF (SUPSYM) THEN 1442 WRITE (LUPRI,1751) ISYM,REP(ISYM-1), 1443 & (SCRA(IORBI+I),ISSMO(IORBI+I),I=1,NLAST) 1444 ELSE 1445 WRITE (LUPRI,1750) ISYM,REP(ISYM-1), 1446 & (SCRA(IORBI+I),I=1,NLAST) 1447 END IF 1448 200 CONTINUE 1449 1450! save orbital energies for MOLDEN 1451 IF (MOLDEN) CALL MOLDEN_MOS(2,SCRA,DUMMY,DUMMY,DUMMY,DUMMY) 1452 1453C take care of 1-electron case and no virtuals case 1454 IF (EHOMO .EQ. -DBIG) EHOMO = D0 1455 IF (ELUMO .EQ. DBIG) ELUMO = D0 1456 WRITE (LUPRI,'(2(/A,F15.8,A,I2,A)/A/A,F15.8,A)') 1457 & ' E(LUMO) :',ELUMO,' au (symmetry',ILUMO,')', 1458 & ' - E(HOMO) :',EHOMO,' au (symmetry',IHOMO,')', 1459 & ' ------------------------------------------', 1460 & ' gap :',ELUMO-EHOMO,' au' 1461 IF (NASHT .EQ. 1) THEN 1462 WRITE (LUPRI,'(/A,F15.8,A,I2,A)') 1463 & 'and E(SOMO) :',ESOMO,' au (symmetry',ISOMO,')' 1464 ELSE IF (HSROHF) THEN 1465 WRITE (LUPRI,'(/A)') 1466 & 'and energy and symmetry of HSROHF open shell orbitals:' 1467 DO IW = NISHT+1,NOCCT 1468 I = ISX(IW) 1469 ISYM = ISMO(I) 1470 WRITE (LUPRI,'(F15.8,I10)') SCRA(I), ISYM 1471 END DO 1472 END IF 1473 IF (ELUMO-EHOMO .LT. GAPMIN) THEN 1474 WRITE(LUPRI,'(/A)') 1475 & ' INFO: E(LUMO) - E(HOMO) small or negative.' 1476 END IF 1477 IF (EHOMO .GT. 0.0d0) THEN 1478 WRITE(LUPRI,'(/A)') ' WARNING: E(HOMO) positive!' 1479 END IF 1480 IF ( VALCON .GT. EMYCNV ) THEN 1481 WRITE(LUPRI,'(/A//A,1P,D10.2)') ' NOTE:'// 1482 * ' MOLECULAR ORBITALS ARE NOT CANONICAL HARTREE-FOCK ORBITALS', 1483 * ' Largest off-diagonal Fock matrix element is',VALCON 1484 IF (IPRINT .GT. 10) THEN 1485 WRITE(LUPRI,'(/A)') ' The total Fock matrix :' 1486 CALL OUTPKB(FD,NORB,NSYM,1,LUPRI) 1487 END IF 1488 ELSE IF (IPRINT .GT. 10) THEN 1489 WRITE(LUPRI,'(/A/A,1P,D10.2)') 1490 * ' Deviation from canonical Hartree-Fock orbitals :', 1491 * ' Largest off-diagonal Fock matrix element is',VALCON 1492 END IF 1493 9999 CALL QEXIT('RHFENR') 1494 RETURN 1495C 1496C1730 FORMAT(//' Hartree-Fock electronic energy:',F30.15) 1497C1732 FORMAT( /' Hartree-Fock total energy:',F30.15) 1498 1740 FORMAT( /' Sym Hartree-Fock orbital energies') 1499 1741 FORMAT( /' Sym Hartree-Fock orbital energies', 1500 & ' (supersymmetry in parenthesis)') 1501 1840 FORMAT( /' Sym Kohn-Sham orbital energies') 1502 1841 FORMAT( /' Sym Kohn-Sham orbital energies', 1503 & ' (supersymmetry in parenthesis)') 1504 1750 FORMAT(/I1,1X,A3,5F15.8/,(5X,5F15.8) ) 1505 1751 FORMAT(/I1,1X,A3,4(F15.8,' (',I2,')')/,(5X,4(F15.8,' (',I2,')')) ) 1506C 1507C 1508C *** end of subroutine RHFENR 1509C 1510 END 1511C /* Deck optst */ 1512 SUBROUTINE OPTST (INDXCI,WRK,LFREE) 1513C 1514C Written december 83 by Hans Joergen Aa. Jensen and Hans Agren 1515C Revisions: 1516C 29-Jul-1984 hjaaj // 19-May-1984 hjaaj/ha 1517C 7-Jan-1985 hjaaj (ICI0=5 option) 1518C 22-Apr-1985 hjaaj (ICI0=6 option) 1519C 6-Jul-1985 hjaaj (CICTL for ICI0=1 option) 1520C 4-Aug-1986 hjaaj (removed obsolete ICI0=2, =3 options) 1521C 1522C MOTECC-90: Some of the algorithms used in OPTST are 1523C described in Chapter 8 Section E.6 of MOTECC-90 1524C "Initial CI iterations for Frozen Orbitals" 1525C 1526C Purpose: 1527C Obtain CI vector for the expansion point in the first macro 1528C iteration and a start guess for the other vectors in a 1529C simultaneous expansion algorithm for solving the eigenvalue/ 1530C eigenvector problem. To write molecular orbitals on LUIT1. 1531C The routine is controlled by the ICIO parameter. 1532C 1533C if (ICI0<0) then 1534C select -ICIO as start configuration 1535C else 1536C 1537C ICI0>100 : transform to natural orbitals 1538C Should not be used, use FLAG instead. 1539C 1540C mod(ICI0,100) -- 1541C =1 : Compute start guess of CI-vector from a new H-diagonal. 1542C The diagonal is computed in and returned from the HDIAG 1543C routine. This call is preceded by a transformation call. 1544C The CI-vector is afterwards written on LUIT1. 1545C If MAXCIT .gt. 0, or the ISTATE configuration is degenerate, 1546C then call CICTL to do MAXCIT iterations (3 if deg.). 1547C If MAXCIT .lt. 0, then decide here how many CI iterations 1548C 1549C =2 : *OBSOLETE* (860804 hjaaj); now as ICI0 = 4 (880616 hjaaj) 1550C 1551C =3 : *OBSOLETE* (860804 hjaaj); now as ICI0 = 4 (880616 hjaaj) 1552C 1553C =4 : Start CI-vector is already on LUIT1 (which is checked). 1554C 1555C =5 : Restart, as ICI0=4 plus read EMCOLD,BETA,RTRUST from LUIT1. 1556C ('NEOSAVE' in SIRSAV) 1557C 1558C =6 : Geometry walk, new geometry ('GEOSAVE' in SIRSAV) 1559C 1560C =7 : Start using MC update information ('UPDSAVE' in SIRSAV) 1561C end if 1562C 1563 use lucita_mcscf_ci_cfg 1564#include "implicit.h" 1565C 1566 DIMENSION INDXCI(*),WRK(LFREE) 1567C 1568 PARAMETER ( D0=0.0D0, DP5=0.5D0, D1=1.0D0 ) 1569#include "dummy.h" 1570C 1571C Used from common blocks: 1572C exeinf.h : NEWCMO 1573C INFINP: ISTATE,LROOTS,FLAG(*),MAXCIT,THRGRD,ITRLVL,? 1574C INFVAR: NCONF (p.t. only NCONF) 1575C INFIND: ISX() 1576C INFORB: NCMOT,? 1577C INFOPT: EMCOLD,BETA,RTRUST (read from LUIT1 when restart, ICI0=5) 1578C INFDIM: IVEX4,INTOT 1579C INFTAP: LUIT1,? 1580C INFPRI: ? 1581C 1582#include "maxorb.h" 1583#include "maxash.h" 1584#include "priunit.h" 1585#include "exeinf.h" 1586#include "infinp.h" 1587#include "infvar.h" 1588#include "infind.h" 1589#include "inforb.h" 1590#include "infopt.h" 1591#include "infdim.h" 1592#include "inftap.h" 1593#include "infpri.h" 1594#include "gnrinf.h" 1595#include "dftcom.h" 1596C 1597C *** Externals: 1598C 1599 LOGICAL FNDLAB, FNDLB2 1600C 1601C *** Local variables: 1602C 1603 LOGICAL HFCALC, RSTFLG(4), TRCNO, LSAV15 1604 LOGICAL DFT_SPINDNS_SAVE, DFT_LOCALSPIN_SAVE 1605 CHARACTER*8 TABLE(8), LBLINF(2), LBLDAT(2), STAR8 1606 CHARACTER*6 CNOLBL 1607C 1608C *** data: 1609C 1610 DATA TABLE/'EODATA ', 'OLDORB ', 'STARTVEC', 'CIDIAG2 ', 1611 * 'RESTART ', 'NEWORB ', 'GEOWALK ', 'NATOCC '/ 1612 DATA STAR8/'********'/ 1613C 1614 CALL QENTER('OPTST ') 1615 1616 IF (ISTATE .GT. NCONF) THEN 1617 WRITE(LUPRI,'(//A,I0,A,I0)') 1618 & 'FATAL ERROR: requested .STATE ',ISTATE, 1619 & ' is greater than number of configurations: ',NCONF 1620 CALL QUIT('ERROR: '// 1621 & 'requested .STATE is greater than number of configurations') 1622 END IF 1623 1624C Transfer ICI0 to local variable ICSTRT 1625 ICSTRT = ICI0 1626C 1627 CALL GETDAT(LBLDAT(1),LBLDAT(2)) 1628 NC4 = MAX(4,NCONF) 1629 NCMOT4 = MAX(4,NCMOT) 1630 TRCNO = .FALSE. 1631 HFCALC = (NCONF .EQ. 1) 1632 10 CONTINUE 1633C 1634C *** Go to appropriate section to obtain start guess for CI vectors 1635C (as specified with the input parameter ICI0 (here ICSTRT)) 1636C However, if RHF or NCONF.eq.1, simply write CREF(1) = D1 to LUIT1 1637C 1638 IF (HFCALC) THEN 1639 ICSTRT = MOD(ICSTRT,100) 1640 IF (ICSTRT .EQ. 1) ICSTRT = -1 1641 TRCNO = FLAG(15) 1642 CNOLBL = 'ONLYFD' 1643 ELSE IF (FLAG(15)) THEN 1644 TRCNO = .TRUE. 1645 IF (FLAG(48)) THEN 1646 CNOLBL = 'ONLYFD' 1647 ELSE IF (FLAG(46)) THEN 1648 CNOLBL = 'ONLYNO' 1649 ELSE 1650 CNOLBL = 'FD+NO ' 1651 END IF 1652 ICSTRT = MOD(ICSTRT,100) 1653 ELSE IF (ICSTRT .GE. 100) THEN 1654 TRCNO = .TRUE. 1655 CNOLBL = 'ONLYNO' 1656 ICSTRT = MOD(ICSTRT,100) 1657 ELSE 1658 TRCNO = .FALSE. 1659 END IF 1660C 1661C *** ICSTRT < 0 *** 1662C === Start guess = CSF no. -ICSTRT 1663C 1664 IF (ICSTRT.LT.0) THEN 1665 JCSTRT = -ICSTRT 1666C ... 930521-hjaaj: JCSTRT introduced to 1667C fix AIX xlf v.2.3 error for 'xlf -O'. 1668 IF (JCSTRT.GT.NCONF) THEN 1669 WRITE (LUPRI,5010) JCSTRT,NCONF 1670 WRITE (LUERR,5010) JCSTRT,NCONF 1671 CALL QTRACE(LUERR) 1672 CALL QUIT( 1673 & 'ERROR (OPTST) non-existent configuration specified') 1674 END IF 1675 IF (LROOTS.GT.1) THEN 1676 WRITE (LUPRI,5020) ICSTRT 1677 WRITE (LUERR,5020) ICSTRT 1678 CALL QTRACE(LUERR) 1679 CALL QUIT('(OPTST) ICSTRT must be pos. when LROOTS .gt. 1') 1680 END IF 1681 KCMO = 1 1682 KCREF = KCMO + NCMOT 1683 KW11 = KCREF + NCONF 1684 IF (KW11 .GT. LFREE) CALL ERRWRK('OPTST.CISTRT 0',KW11,LFREE) 1685C 1686C Read NEWORB and write as OLDORB 1687C 1688 REWIND LUIT1 1689 CALL MOLLB2(TABLE(6),LBLINF,LUIT1,LUPRI) 1690 CALL READT(LUIT1,NCMOT,WRK(KCMO)) 1691 CALL NEWIT1 1692 WRITE (LUIT1) STAR8,LBLINF,TABLE(2) 1693 CALL WRITT(LUIT1,NCMOT4,WRK(KCMO)) 1694C 1695C Write STARTVEC 1696C 1697 CALL DZERO(WRK(KCREF),NCONF) 1698 WRK(KCREF-1+JCSTRT) = D1 1699 LBLINF(1) = ' 1 1' 1700 LBLINF(2) = 'CISTRT 0' 1701 WRITE (LUIT1) STAR8,LBLINF,TABLE(3) 1702 CALL WRITT(LUIT1,NC4,WRK(KCREF)) 1703C 1704C Write NEWORB 1705C 1706 WRITE (LUIT1) STAR8,LBLDAT(1),LBLINF(2),TABLE(6) 1707 CALL WRITT(LUIT1,NCMOT4,WRK(KCMO)) 1708 NEWCMO = .TRUE. 1709C single configuration has always natural orbitals 1710 IF (CNOLBL .EQ. 'ONLYNO') THEN 1711 TRCNO = .FALSE. 1712 END IF 1713 GO TO 9000 1714 END IF 1715 5010 FORMAT(/' OPTST FATAL ERROR,' 1716 * /' specified start configuration no.',I10, 1717 * ' is non-existent (last is no.',I10,')') 1718 5020 FORMAT(/' OPTST FATAL ERROR,' 1719 * /' negative CISTRT option not allowed',I10, 1720 * /' when LROOTS is greater then one.') 1721C 1722C 1723C *** ICSTRT > 0 *** 1724C 1725C 1726 IF (ICSTRT .EQ. 2 .OR. ICSTRT .EQ. 3) ICSTRT = 4 1727C ... 880616 hjaaj (CISTRT = 2 has now same meaning as in CIST). 1728 GO TO (100,50,50,400,400,400,400) ICSTRT 1729 50 CONTINUE 1730 WRITE (LUPRI,5030) ICSTRT 1731 WRITE (LUERR ,5030) ICSTRT 1732 CALL QTRACE(LUERR) 1733 CALL QUIT('ERROR, illegal control option in OPTST') 1734 5030 FORMAT(/' OPTST FATAL ERROR, ICSTRT =',I4,' is not defined') 1735C 1736C *** ICSTRT = 1 *** 1737C *** Compute H-diagonal and then select starting configurations 1738C (If MAXCIT .eq. 0 then CICTL selects smallest diag. elements of H; 1739C if MAXCIT .gt. 0 then CICTL follows up with MAXCIT CI iterations) 1740C 1741 100 CONTINUE 1742C 1743 KCMO = 1 1744 KW21 = KCMO + NCMOT 1745 LW21 = LFREE - KW21 1746 KECI = KW21 1747 KICROO = KECI + LROOTS 1748 KW22 = KICROO + LROOTS 1749 LW22 = LFREE - KW22 1750 IF (LW22 .LT. 0) CALL ERRWRK('OPTST.CICTL',-KW22,LFREE) 1751C 1752C A single configuration will already have natural orbitals 1753 IF (MAXCIT .LE. 0) THEN 1754 IF (CNOLBL .EQ. 'ONLYNO') THEN 1755 TRCNO = .FALSE. 1756 ELSE IF (CNOLBL .NE. 'TRTOFC') THEN 1757 CNOLBL = 'ONLYFD' 1758 END IF 1759 END IF 1760C 1761C Read NEWORB and write as OLDORB 1762C 1763 REWIND LUIT1 1764 CALL MOLLB2(TABLE(6),LBLINF,LUIT1,LUPRI) 1765 CALL READT(LUIT1,NCMOT,WRK(KCMO)) 1766 CALL NEWIT1 1767 WRITE (LUIT1) STAR8,LBLINF,TABLE(2) 1768 CALL WRITT(LUIT1,NCMOT4,WRK(KCMO)) 1769C 1770C --- Call TRACTL first, if needed 1771C (do it out here unless we later will transform to NO, 1772C CICTL only does an ITRLVL = 0 transformation) 1773C 1774 IF (.NOT.FLAG(14) .AND. .NOT.TRCNO) THEN 1775 CALL TRACTL(ITRLVL,WRK(KCMO),WRK(KW21),LW21) 1776C CALL TRACTL(ITRLVL,CMO,WRK,LFREE) 1777 FLAG(14) = .TRUE. 1778 END IF 1779C 1780 IF (MAXCIT .LT. 0) THEN ! user not specified .MAX CI 1781 MAXCIT = 8 1782 THRCIX = 1.D-2 ! Usually initial orbital gradient is larger ... 1783 ELSE 1784 THRCIX = DP5*THRGRD 1785 END IF 1786 IF (ICHECK .GT. 0) THEN 1787 IF (MAXCIT .LE. 2) THEN 1788 WRITE (LUPRI,'(//A/A)') 1789 & ' --- INFO OPTST, symmetry check of CI vectors'// 1790 & ' only works after CI iterations.', 1791 & ' ".MAX CI ITERATIONS" reset to 3' 1792 MAXCIT = 3 1793 END IF 1794 IF (LROOTS .EQ. NROOTS) THEN 1795C ... if (LROOTS .gt. NROOTS) then LROOTS has been 1796C specified in input; use this value 1797 LROOTS = MAX(LROOTS,2*ISTATE+2) 1798 LROOTS = MIN(LROOTS,NCONF) 1799 END IF 1800 END IF 1801C 1802 LSAV15 = FLAG(15) 1803 FLAG(15) = .FALSE. 1804 ISTASV = ISTACI 1805 IF (ICHECK .EQ. 1 .AND. ISTATE .GT. 1) THEN 1806C ... we do not know which root will become ISTATE 1807 ISTACI = 0 1808 ELSE 1809 ISTACI = ISTATE 1810 END IF 1811C 1812 IF (DOMCSRDFT) THEN 1813 ADDSRI = .FALSE. 1814 DOCISRDFT = .TRUE. 1815 IF (DFT_SPINDNS .OR. DFT_LOCALSPIN)THEN 1816 WRITE(LUPRI,'(/A)') 1817 & 'INFO : spin density ignored in initial CI iterations.' 1818 END IF 1819 DFT_SPINDNS_SAVE = DFT_SPINDNS 1820 DFT_SPINDNS = .FALSE. 1821 DFT_LOCALSPIN_SAVE = DFT_LOCALSPIN 1822 DFT_LOCALSPIN = .FALSE. 1823 END IF 1824 CALL CICTL(1,LROOTS,MAXCIT,THRCIX,WRK(KCMO),INDXCI, 1825 * WRK(KECI),ICONV,WRK(KICROO),WRK(KW22),LW22) 1826C CALL CICTL(ICICTL,NCROOT,MAXITC,THRCIX,CMO,INDXCI, 1827C * ECI,ICONV,ICROOT,WRK,LFREE) 1828 IF (DOMCSRDFT) THEN 1829 ADDSRI = .TRUE. 1830 DOCISRDFT = .FALSE. 1831 DFT_SPINDNS = DFT_SPINDNS_SAVE 1832 DFT_LOCALSPIN = DFT_LOCALSPIN_SAVE 1833 END IF 1834 FLAG(15) = LSAV15 1835 ISTACI = ISTASV 1836C 1837C --- Call CICHCK to remove CI vectors that do not have the same 1838C symmetry as ISTATE (ICHECK =1) or as the state of the lowest 1839C symmetry (ICHECK=2). (This call can only be made after some 1840C CI-iterations have been done.) 1841C 1842 IF (ICHECK .GE. 0 .AND. MAXCIT .GT. 0) THEN 1843 CALL CICHCK(WRK,LFREE,ICHECK) 1844C CALL CICHCK(WRK,LWRK,CCHECK) 1845 END IF 1846C 1847 IF (MAXCIT .GT. 0) LROOTS = NROOTS 1848C 1849C 1850 GO TO 9000 1851C (See what codes means at beginning of this subroutine) 1852C 1853C === start guess should already be on LUIT1 (SIRIUS.RST) 1854C (but we test to be sure!) 1855C 1856 400 CONTINUE 1857 REWIND LUIT1 1858 IF ( .NOT.FNDLB2(TABLE(3),LBLINF,LUIT1) ) GO TO 490 ! label 'STARTVEC' 1859 IF ( ICSTRT.EQ.5 ) THEN ! .RESTART 1860 REWIND LUIT1 1861 IF ( .NOT.FNDLAB(TABLE(2),LUIT1) ) GO TO 490 ! label 'OLDORB ' for backstep after restart 1862 IF ( LBLINF(2).EQ.'NRSAVE ') FLAG(39) = .TRUE. ! if user has specified .RESTART, then respect shift to NR 1863 END IF 1864 IF (FLAG(39)) GO TO 495 1865C ... if (ONLYNR) one vector is all that is needed, 1866C so no reason to do more checking. 1867C We now check if sufficient vectors available using label information. 1868 READ (LBLINF(1),'(2I4)',ERR=405) LROIT1,JSTA 1869 IF (LROIT1 .LT. LROOTS) GO TO 493 1870 GO TO 495 1871C error reading LBLINF(1) -- use the old check 1872 405 LROIT1 = 0 1873 DO 410 I = 1,LROOTS 1874 READ (LUIT1,END=493,ERR=493) LBLINF 1875 IF (LBLINF(1) .EQ. STAR8) GO TO 493 1876 LROIT1 = LROIT1 + 1 1877 410 CONTINUE 1878 GO TO 495 1879C 1880 490 CONTINUE 1881 REWIND LUIT1 1882 CALL DMPLAB(LUIT1, LUPRI) 1883 IF (ICSTRT .EQ. 4) THEN 1884 NWARN = NWARN + 1 1885 WRITE (LUPRI,'(//A/A//)') 1886 & '@@*** (OPTST) MC start vectors not found on SIRIUS.RST', 1887 & '@@*** WARNING lowest diagonal elements used.' 1888 ICSTRT = 1 1889 GO TO 10 1890 ELSE IF (ICSTRT .EQ. 5) THEN 1891 WRITE (LUPRI,'(//A/A//)') 1892 & '@@*** (OPTST) MC start vectors or OLDORB'// 1893 & ' not found on SIRIUS.RST,', 1894 & '@@ .RESTART not possible.' 1895 CALL QTRACE(LUPRI) 1896 CALL QUIT( 'ERROR: '// 1897 & ' All necessary info for .RESTART not found on SIRIUS.RST.') 1898 ELSE 1899 WRITE (LUPRI,'(//A/A//)') 1900 & '@@*** (OPTST) Old CI vectors not found on SIRIUS.RST,', 1901 & '@@ no recover when CISTRT 6 or 7,'// 1902 & ' optimization abandoned.' 1903 CALL QTRACE(LUPRI) 1904 CALL QUIT('ERROR: '// 1905 & 'Old CI vectors not found on SIRIUS.RST as expected.') 1906 END IF 1907 493 CONTINUE 1908 REWIND LUIT1 1909 IF ( LROIT1 .GE. NROOTS) THEN 1910 WRITE (LUPRI,'(//A,I3,A/T14,A,I3,A,I3/)') 1911 * ' *** (OPTST) only',LROIT1,' start vectors found ', 1912 * 'on SIRIUS.RST,', 'LROOTS reset from',LROOTS,' to',LROIT1 1913 LROOTS = LROIT1 1914 ELSE 1915 GO TO 490 1916 END IF 1917 495 CONTINUE 1918C 1919C =================== 1920C 1921C READ RESTART INFORMATION 1922C 1923C =================== 1924C 1925 REWIND LUIT1 1926 IF (ICSTRT.EQ.5) THEN 1927C "NEOSAVE" or "NRSAVE" restart; 1928C 1) get trust radius and beta 1929C 2) if (absorption) check if absorption switched off 1930 IF ( .NOT.FNDLB2(TABLE(5),LBLINF,LUIT1) ) THEN 1931 WRITE (LUPRI,'(//A)') 1932 & ' (OPTST) Error, restart information not found on SIRIUS.RST.' 1933 CALL QUIT( 1934 & '(OPTST) Error, restart information not found on SIRIUS.RST.') 1935 END IF 1936 READ (LUIT1) EMCOLD,DEPRED,BETA,RDUMMY,RTRUST,ITMAC 1937 WRITE (LUPRI,5000) EMCOLD,ITMAC,BETA,RTRUST,EMCOLD+DEPRED 1938 IF (LBLINF(2) .NE. STAR8) THEN 1939C this is a new restart file (6-Aug-1986) 1940 READ (LUIT1) DUM,DUM,DUM,DUM,RSTFLG 1941 IF ((FLAG(51) .OR. FLAG(52) .OR. FLAG(53)) .AND. 1942 * .NOT.FLAG(54)) THEN 1943 FLAG(51) = RSTFLG(1) 1944 FLAG(52) = RSTFLG(2) 1945 FLAG(53) = RSTFLG(3) 1946 IF (.NOT. (FLAG(51) .OR. FLAG(52) .OR. FLAG(53)) ) 1947 * WRITE (LUPRI,'(/A/)') 1948 * ' *** absorption switched off in a previous iteration.' 1949 END IF 1950 IF (.NOT.FLAG(38) .AND. .NOT. FLAG(39)) THEN 1951C if (not always neo and not always nr) then 1952 FLAG(39) = RSTFLG(4) 1953 IF (FLAG(39)) WRITE (LUPRI,'(/A/)') 1954 * ' *** Switched from NEO to NR in a previous iteration.' 1955 END IF 1956 END IF 1957 TRCNO = .FALSE. 1958C ... no TRCNO when restart, presumably done in prev. iteration 1959 ELSE IF (ICSTRT.EQ.6) THEN 1960C "GEOSAVE" restart 1961 IF (ISTATE .GT. 1) THEN 1962 IF (FLAG(39)) THEN 1963 WRITE (LUPRI,6008) ISTATE 1964 ELSE 1965 WRITE (LUPRI,6010) 1966 CALL QTRACE(LUPRI) 1967 CALL QUIT('(OPTST) NEO and ISTATE .gt. 1 for GEO WALK') 1968 END IF 1969 END IF 1970 IF ( FNDLAB(TABLE(7),LUIT1) ) THEN 1971 READ (LUIT1) EMCOLD,DEPRED,REJWMI,REJWMA 1972 WRITE (LUPRI,6000)EMCOLD,DEPRED,EMCOLD+DEPRED,REJWMI,REJWMA 1973 ELSE 1974C 1975C We do not complain about lack of GEOWALK if .OPTIMI, K.Ruud, Nov.29-96 1976C OK! But we must anyway reset ICSTRT and ICI0 to 4 hjaaj/961213. 1977C 1978 IF (.NOT. OPTNEW) THEN 1979 WRITE (LUPRI,'(//A/A/)') 1980 * '@@ INFO : Label "GEOWALK " not found on SIRIUS.RST,', 1981 * '@@ INFO : calculation continues without'// 1982 * ' test for rejection of geometry change.' 1983 END IF 1984 ICSTRT = 4 1985 ICI0 = 4 1986C ICI0 modified such that SIROPT does not perform 1987C walk rejection check 1988 END IF 1989 ELSE IF (ICSTRT.EQ.7) THEN 1990C "UPDSAVE" restart 1991 IF (ISTATE .GT. 1) THEN 1992 WRITE (LUPRI,7010) 1993 CALL QTRACE(LUPRI) 1994 CALL QUIT('ERROR (OPTST) ISTATE .gt. 1 for ICSTRT .eq. 7') 1995 END IF 1996 END IF 1997 REWIND LUIT1 1998 GO TO 9000 1999 5000 FORMAT(/' (OPTST) RESTART, Old MCSCF energy:',T50,F30.15, 2000 * /T18,'for macro iteration no.',T49,I10 2001 * /T23,'BETA value:',T49,F10.2, 2002 * /T23,'Trust radius:',T44,F15.5, 2003 * /T18,'Predicted new MCSCF energy:',T50,F30.15) 2004 6000 FORMAT(//' (OPTST) Geometry walk', 2005 * /T16,'Old MCSCF energy: ',F30.15, 2006 * /T16,'Predicted energy change:',F30.15, 2007 * /T16,'Predicted new energy: ',F30.15, 2008 * //T16,'Ratios for rejecting walk,', 2009 * /T16,' Minimum :',F14.5, 2010 * /T16,' Maximum :',F14.5) 2011 6008 FORMAT(//' (OPTST) INFO: Newton-Raphson used in geometry walk', 2012 * ' for state no.',I3/ 2013 * ' INFO: no absolute guarantee for staying', 2014 * ' on same electronic surface'//) 2015 6010 FORMAT(//' (OPTST) ERROR, ICSTRT = 6 is not', 2016 * ' allowed for NEO excited state optimization') 2017 7010 FORMAT(//' (OPTST) ERROR, ICSTRT = 7 is not', 2018 * ' allowed for excited state optimization') 2019C 2020C 2021C 2022 9000 CONTINUE 2023 REWIND LUIT1 2024 IF ( .NOT.FNDLAB(TABLE(6),LUIT1) ) THEN 2025 CALL READMO(WRK,7) 2026 CALL NEWORB(LBLDAT(2),WRK,.FALSE.) 2027 END IF 2028 IF (TRCNO) THEN 2029C -- check if already done, e.g. by CISAVE 900720/hjaaj 2030 REWIND (LUIT1) 2031 IF (FNDLB2(TABLE(6),LBLINF,LUIT1)) THEN 2032 IF (LBLINF(2) .EQ. '(CNOORB)') TRCNO = .FALSE. 2033 END IF 2034 END IF 2035 IF (TRCNO) THEN 2036C 2037C -- Transform to natural orbitals -- 2038C 2039 KCMO = 1 2040 KCVECS = KCMO + NCMOT 2041 KAOCC = KCVECS + LROOTS*NCONF 2042 KW31 = KAOCC + NASHT 2043 LW31 = LFREE - KW31 2044 IF (LW31 .LE. 0) CALL ERRWRK('OPTST.SIRCNO',-KW31,LFREE) 2045 CALL READMO(WRK(KCMO),9) 2046 IF (HFCALC) THEN 2047 WRK(KCVECS) = D1 2048 ICREF = 1 2049 NCVECS = 1 2050 JSTA = 1 2051 ELSE 2052 REWIND (LUIT1) 2053C ... search for "STARTVEC" 2054 CALL MOLLB2(TABLE(3),LBLINF,LUIT1,LUPRI) 2055 READ (LBLINF(1),'(2I4)',ERR=405) LROIT1,JSTA 2056 IF (JSTA .LT. 0) THEN 2057 ICREF = 1 2058 NCVECS = 1 2059 JSTA = -ISTATE 2060 ELSE IF (LROOTS .LT. ISTATE) THEN 2061C ... to handle .NR ALWAYS calculation following 2062C NEO calculation, e.g. for double core hole. 2063 ICREF = 1 2064 NCVECS = 1 2065 JSTA = -ISTATE 2066 DO 9005 I = 1,ISTATE-1 2067 READ (LUIT1) 2068 9005 CONTINUE 2069 ELSE 2070 ICREF = ISTATE 2071 NCVECS = LROOTS 2072 JSTA = ISTATE 2073 END IF 2074 JCVECS = KCVECS 2075 DO 9010 I = 1,NCVECS 2076 CALL READT(LUIT1,NCONF,WRK(JCVECS)) 2077 JCVECS = JCVECS + NCONF 2078 9010 CONTINUE 2079 END IF 2080! set logical flag for new reference CI vector in WRK(KCREF) - 2081! in the interface to lucita we use this info to save/broadcast the new 2082! reference vector on lucita internal sequential/parallel MPI-I/O files 2083 vector_exchange_type1 = 1 2084 vector_update_mc2lu_lu2mc((2-1)*vector_exchange_types+ 2085 & vector_exchange_type1) = .true. 2086 2087 CALL SIRCNO(CNOLBL,NCVECS,ICREF,WRK(KCVECS), 2088 * WRK(KCMO),WRK(KAOCC),INDXCI,WRK(KW31),1,LW31) 2089C CALL SIRCNO(KEYCNO,NCVECS,ICREF,CREF,CMO,AOCC, 2090C * INDXCI,WRK,KFRSAV,LFRSAV) 2091C 2092 REWIND (LUIT1) 2093 CALL MOLLAB(TABLE(3),LUIT1,LUPRI) 2094C search for "STARTVEC" 2095 WRITE (LBLINF(1),'(2I4)') NCVECS,JSTA 2096 LBLINF(2) = 'CNOSAVE ' 2097 BACKSPACE LUIT1 2098C write startvec with new vectors 2099 WRITE (LUIT1) STAR8,LBLINF,TABLE(3) 2100 JCVECS = KCVECS 2101 DO 9020 I = 1,NCVECS 2102 CALL WRITT(LUIT1,NC4,WRK(JCVECS)) 2103 JCVECS = JCVECS + NCONF 2104 9020 CONTINUE 2105C write the natural orbitals to LUIT1 label NEWORB 2106 LBLINF(2) = '(CNOORB)' 2107 WRITE (LUIT1) STAR8,LBLDAT(1),LBLINF(2),TABLE(6) 2108 CALL WRITT(LUIT1,NCMOT4,WRK(KCMO)) 2109 NEWCMO = .TRUE. 2110 IF (NASHT .EQ. 1 .OR. 2111 & (CNOLBL .NE. 'ONLYFD' .AND. NASHT .GT. 1)) THEN 2112 IF (NASHT.EQ.1) WRK(KAOCC) = NACTEL 2113 WRITE (LUIT1) STAR8,LBLDAT(1),LBLINF(2),TABLE(8) 2114 KOCCUP = KW31 2115 DO, I = 1, NISHT 2116 IX = ISX(I) 2117 WRK(KOCCUP-1+IX) = 2.0D0 2118 END DO 2119 DO, I = 1, NASHT 2120 IX = ISX(NISHT+I) 2121 WRK(KOCCUP-1+IX) = WRK(KAOCC-1+I) 2122 END DO 2123 NORBT4 = MAX(NORBT,4) 2124 CALL WRITT(LUIT1,NORBT4,WRK(KOCCUP)) 2125 END IF 2126 IF (.NOT.FLAG(34)) FLAG(14) = .FALSE. 2127C if (int.transf. needed in optimization) 2128 END IF 2129C 2130C *** Append label EODATA, if not already there 2131C 2132 REWIND LUIT1 2133 IF ( .NOT.FNDLAB(TABLE(1),LUIT1) ) THEN 2134 WRITE (LUIT1) STAR8,LBLDAT,TABLE(1) 2135 END IF 2136 REWIND LUIT1 2137C 2138C 2139C *** End of subroutine OPTST. 2140C 2141C Set LROOTS negative to tell SIRNEO to use LROOTS and not NROOTS. 2142C 2143 LROOTS = -LROOTS 2144 CALL QEXIT('OPTST ') 2145 RETURN 2146 END 2147C /* Deck sirstp */ 2148 SUBROUTINE SIRSTP(ISTEP,CMO,IBNDX,INDXCI,WRK,LFREE) 2149C 2150C Written 8-Apr-1984 by Hans Jorgen Aa. Jensen and Hans Agren 2151C Last revision 14-May-1984 hjaaj 2152C 21-May-1986 hjaaj (geowlk step control) 2153C 2154C Purpose: 2155C STEP CONTROL -- check if step was too large, if actual 2156C energy deviates too much from the one predicted by last 2157C macro iteration. 2158C If step is too large, read back reduced L-matrix from 2159C last macro iteration and find the level-shift (beta value) 2160C giving the step length corresponding to the revised trust 2161C radius. Then read the b-vectors from last macro iteration 2162C to find the revised eigenvectors and write those on LUIT1 2163C (in SIRSAV). 2164C Signal to calling program with ISTEP = 0 if step is OK, 2165C with ISTEP = 1 if step is to large. 2166C 2167C MOTECC-90: The purpose of this module, SIRSTP, and the algorithms 2168C used are described in Chapter 8 Section C.5 of MOTECC-90 2169C "Step-Control Algorithm" 2170C 2171C 2172C Input: 2173C ISTEP, =1 check geometry walk step 2174C =-1 restart step (step cannot be rejected, 2175C but trust radius will be adjusted) 2176C =-2 wave function is converged, only generate info 2177C else normal MC step check 2178C Output: 2179C ISTEP, = 0 if step is OK 2180C = 1 if step is too large 2181C =-1 if no check because close to convergence 2182C 2183C Scratch: 2184C WRK 2185C 2186#include "implicit.h" 2187 DIMENSION CMO(*), IBNDX(*), INDXCI(*), WRK(LFREE) 2188 PARAMETER ( DP1 = 0.1D0, DP5 = 0.5D0, 2189 * D0 = 0.0D0, D1 = 1.0D0, D2 = 2.0D0 ) 2190 PARAMETER ( THDE = 1.D-10 ) 2191#include "dummy.h" 2192C 2193C Used from common blocks: 2194C INFINP: ISTATE,FLAG(*),? 2195C INFORB: NCMOT 2196C INFVAR: NWOPT 2197C INFOPT: NREDL,EMCSCF,EMCOLD,DEPRED,BETA,GAMMA,SHFLVL, 2198C STPLEN,STPMAX,STPINC,STPRED, 2199C RTRUST,RTTOL,RATMIN,RATGOD,RATREJ,REJWMI,REJWMA, GRDNRM 2200C ITINFO: DINFO(*) 2201C INFTAP: LUIT1 2202C INFPRI: P4FLAG(*),P6FLAG(*) 2203C 2204#include "maxorb.h" 2205#include "priunit.h" 2206#include "infinp.h" 2207#include "inforb.h" 2208#include "infvar.h" 2209#include "infopt.h" 2210#include "itinfo.h" 2211#include "inftap.h" 2212#include "infpri.h" 2213#include "gnrinf.h" 2214C 2215C -- local: 2216C 2217 CHARACTER*8 TABLE(3), RTNLBL(2) 2218 LOGICAL GEOWLK, NOREJT 2219C 2220C 2221C -- data: 2222C 2223 DATA TABLE /'OLDORB ','RESTART ','LREDUCED'/ 2224C 2225C 2226 CALL QENTER('SIRSTP') 2227C 2228C ***************************************************************** 2229C *** Trust region control 2230C 2231 GEOWLK = .FALSE. 2232 NOREJT = .FALSE. 2233 IF (ISTEP .EQ. 1 .AND. .NOT. OPTNEW) THEN 2234 GEOWLK = .TRUE. 2235 WRITE (LUW4,'(//A/)') 2236 * ' --- step control for geometry walk (ABACUS interface)' 2237 ELSE IF (ISTEP .LT. 0) THEN 2238 NOREJT = .TRUE. 2239 END IF 2240 ISTEP = 0 2241 DEACT = EMCSCF - EMCOLD 2242 2243 IF ( ABS(DEPRED/EMCSCF).GT.THDE ) THEN 2244 RATIO = DEACT / DEPRED 2245 IF (GEOWLK .OR. P4FLAG(7)) WRITE (LUW4,1100) DEACT,DEPRED,RATIO 2246 ELSE 2247 IF (DEPRED .NE. 0.0D0) THEN 2248 RATIO = DEACT / DEPRED 2249 WRITE(LUW4,1100) DEACT, DEPRED, RATIO 2250 WRITE(LUW4,'(A)') ' Close to convergence, ratio set to one.' 2251 ELSE 2252 IF (GEOWLK .OR. P4FLAG(7)) WRITE (LUW4,1110) DEACT,DEPRED 2253 END IF 2254 RATIO = D1 2255 ISTEP = -1 2256 END IF 2257 1100 FORMAT (/' (SIRSTP) Energy difference;', 2258 * /' actual, predicted and ratio:',2F20.15,F12.6) 2259 1110 FORMAT (/' (SIRSTP) Close to convergence, ratio set to one.', 2260 * /' Energy difference; actual and predicted:',1P,2D15.5) 2261C 2262C special test if geometry walk 2263C 2264 IF (GEOWLK) THEN 2265 IF (RATIO .LT. REJWMI .OR. RATIO. GT. REJWMA) THEN 2266 WLKREJ = .TRUE. 2267 WRITE (LUW4,1210) 2268 IF (LUPRI .NE. LUW4) WRITE(LUPRI,1210) 2269 ICI0 = 0 2270 GOTO 9999 2271C 2272C We don't want to abort the program just because of a long step 2273C 2274C WRITE (LUERR,1210) 2275C CALL QTRACE(LUERR) 2276C CALL QUIT('***step control:geometry walk step too long***') 2277 ELSE 2278 WRITE (LUW4,'(//A/)') ' *** GEOMETRY WALK STEP ACCEPTED.' 2279 END IF 2280 EMCOLD = EMCSCF 2281 ICI0 = 0 2282 GO TO 9999 2283 END IF 2284 1210 FORMAT(//' *** GEOMETRY WALK STEP TOO LONG,', 2285 * /' DECREASE WALK TRUST RADIUS AND TRY AGAIN') 2286C 2287C save information for final summary print-out 2288C 2289 DINFO(4) = DEPRED 2290 DINFO(5) = DEACT 2291 DINFO(6) = RATIO 2292C 2293C 2294C Normal MCSCF optimization (not first iteration in geometry walk) 2295C 2296 RATVGD = D1 - DP1*(D1 - RATGOD) 2297 RTSAVE = RTRUST 2298 IF (ISTATE .EQ. 1 .AND. .NOT. FLAG(35)) THEN 2299C (case 1: ground state optimization) 2300C flag(35): tight step control also for ground state optimization 2301 IF (RATIO .LT. RATMIN) THEN 2302 RTRUST = STPRED*RTRUST 2303 IF (RATIO .LT. RATREJ .AND. .NOT. NOREJT) THEN 2304 ISTEP = 1 2305 RTRUST = MIN(RTRUST, STPRED*STPLEN) 2306 END IF 2307 ELSE IF (RATIO.GT.RATGOD) THEN 2308 RTRUST = MIN(STPMAX,STPINC*RTRUST) 2309 IF (RATIO.GT.RATVGD) RTRUST = MIN(STPMAX,STPINC*RTRUST) 2310C extra increase if ratio is very good 2311 END IF 2312 ELSE 2313C (case 2: excited state optimization) 2314 IF (RATIO .LT. RATMIN .OR. RATIO .GT. (D2-RATMIN)) THEN 2315 RTRUST = STPRED*RTRUST 2316 IF (.NOT. NOREJT) THEN 2317 IF (RATIO .LT. RATREJ .OR. RATIO .GT. (D2-RATREJ)) THEN 2318 ISTEP = 1 2319 RTRUST = MIN(RTRUST, STPRED*STPLEN) 2320 END IF 2321 END IF 2322 ELSE IF (RATIO .GT. RATGOD .AND. RATIO .LT. (D2-RATGOD)) THEN 2323 RTRUST = MIN(STPMAX,STPINC*RTRUST) 2324 IF (RATIO .GT. RATVGD .AND. RATIO .LT. (D2-RATVGD)) 2325 * RTRUST = MIN(STPMAX,STPINC*RTRUST) 2326C extra increase if ratio is very good 2327 END IF 2328 END IF 2329C 2330C ***************************************************************** 2331C *** Step is acceptable ( or restart ) : 2332C 2333 IF (ISTEP.LE.0) THEN 2334 EMCOLD = EMCSCF 2335 GO TO 9000 2336 END IF 2337C 2338C ***************************************************************** 2339C *** Step is too large : 2340C 2341 IF (.NOT.P4FLAG(7)) WRITE (LUW4,1100) DEACT,DEPRED,RATIO 2342 IF (LUPRI.NE.LUW4) WRITE (LUPRI,1100) DEACT,DEPRED,RATIO 2343 WRITE (LUW4,7000) 2344 IF (LUPRI.NE.LUW4 .AND. .NOT.P6FLAG(10)) WRITE (LUPRI,7000) 2345 7000 FORMAT(/' (SIRSTP) step is too large -- backstep') 2346C 2347C Core allocation: 2348C 2349 NNREDL = NREDL*(NREDL+1)/2 2350 KREDL = 1 2351 KEVEC = KREDL + NNREDL 2352 KEVAL = KEVEC + NREDL*NREDL 2353 KWRK11 = KEVAL + NREDL 2354 KWRK12 = KWRK11 + NNREDL 2355 KWRK13 = KWRK12 + NREDL 2356 LWRK11 = LFREE + 1 - KWRK11 2357 KXKAP = KWRK11 2358 KWRK21 = KXKAP + NWOPT 2359 LWRK21 = LFREE + 1 - KWRK21 2360C 2361C Recover old orbitals, IBNDX and REDL from previous macro iteration 2362C 2363 REWIND LUIT1 2364 CALL MOLLAB(TABLE(1),LUIT1,LUPRI) 2365 CALL READT(LUIT1,NCMOT,CMO) 2366 CALL MOLLB2(TABLE(2),RTNLBL,LUIT1,LUPRI) 2367 IF (RTNLBL(2) .NE. 'NEOSAVE ') THEN 2368 WRITE (LUW4,'(/A/2A/A)') 2369 * ' Sorry, backstep only implemented for NEO', 2370 * ' Label from last iterations : ',RTNLBL(2), 2371 * ' Program cannot continue (hint: try ".NEO ALWAYS" option).' 2372 CALL QTRACE(LUW4) 2373 CALL QUIT('Sorry, backstep only implemented for NEO') 2374 END IF 2375 READ (LUIT1) DUM,DUM,DUM,DUM,DUM,IDUM, 2376 * MREDL,(IBNDX(I),I=1,MREDL) 2377 IF (MREDL.NE.NREDL) THEN 2378C this should never occur ... 2379 WRITE (LUERR,'(/1X,A,2I5)') 2380 * 'SIRSTP: NREDL on LUIT1 inconsistent with NREDL in common', 2381 * MREDL,NREDL 2382 CALL QTRACE(LUERR) 2383 CALL QUIT('ERROR (SIRSTP) unexpected value of NREDL on LUIT1') 2384 END IF 2385 CALL MOLLAB(TABLE(3),LUIT1,LUPRI) 2386 CALL READT(LUIT1,NNREDL,WRK(KREDL)) 2387C 2388C 2389 CALL NEORED (2,0,0,WRK(KWRK11),DUMMY,BETA,DUMMY,IBNDX, 2390 * NREDL,WRK(KREDL),WRK(KEVAL),WRK(KEVEC),LWRK11) 2391C CALL NEORED (ICTL,NCSIM,NOSIM,BVECS,SVECS,BETVAL,G,IBNDX, 2392C * NREDL,REDL,EVAL,EVEC, LBVECS) 2393C 2394C 2395 CALL UPDBET(WRK(KREDL),WRK(KEVEC),WRK(KEVAL),WRK(KWRK11),LWRK11) 2396C CALL UPDBET(REDL,EVEC,EVAL,WRK,LFREE) 2397 IF (P6FLAG(10)) THEN 2398 WRITE (LUPRI,8010) BETA,RTRUST 2399 CALL OUTPAK(WRK(KREDL),NREDL,1,LUPRI) 2400 WRITE (LUPRI,8012) 2401 WRITE (LUPRI,8015) (WRK(KEVAL-1+I),I=1,NREDL) 2402 LIMP = MIN(((ISTATE+5)/4) * 4, NREDL) 2403 CALL OUTPUT(WRK(KEVEC),1,NREDL,1,LIMP,NREDL,NREDL,1,LUPRI) 2404 END IF 2405 8010 FORMAT(/' (SIRSTP) step was too long;', 2406 * /T11,'new beta and trust radius: ',2F15.8, 2407 * /' - The updated reduced L matrix:') 2408 8012 FORMAT(/' - Eigenvalues and eigenvectors of ', 2409 * 'the updated reduced L matrix:') 2410 8015 FORMAT(/,(10X,1P,4D15.6)) 2411C 2412C *** Save the revised (restricted step) eigenvectors (CI part only) 2413C and rotate the orbitals corresponding to the revised eigenvector. 2414C 2415C 2416 SHFLVL = WRK(KEVAL-1+ISTATE) 2417 CALL SIRSAV ('NEOSAVE',CMO,IBNDX,WRK(KREDL),WRK(KEVEC), 2418 * WRK(KXKAP),INDXCI,WRK(KWRK21),LWRK21) 2419C 2420C CALL SIRSAV (KEYWRD,CMO,IBNDX,REDL,EVEC,XKAP,INDXCI,WRK,LFREE) 2421C 2422C === Revise predicted energy difference 2423C between next macro iteration and this one. 2424C 2425 DEPRED = ( GAMMA*(D1-GAMMA) + DP5*GAMMA*GAMMA* 2426 * (WRK(KEVEC+(ISTATE-1)*NREDL) ** (-2)) ) * 2427 * WRK(KEVAL-1+ISTATE) / (BETA*BETA) 2428C 2429C 2430C ***************************************************************** 2431C *** End of subroutine SIRSTP 2432C 2433 9000 CONTINUE 2434 IF (P4FLAG(7)) WRITE (LUW4,1200) RTSAVE,RTRUST 2435 1200 FORMAT(/' (SIRSTP) Old and new MC trust radius:',2F15.10) 2436 9999 CALL QEXIT('SIRSTP') 2437 RETURN 2438 END 2439C /* Deck timopt */ 2440 SUBROUTINE TIMOPT(KEY,LUPRI,WRK,LFREE) 2441C 2442C Written 18-Dec-1984 Hans Joergen Aa. Jensen. 2443C Revised 891201-hjaaj: LUPRI in parameter list. 2444C 2445C Purpose: 2446C Write timing statistics to LUPRI for central (time-consuming) 2447C routines used in MC optimization, that is from SIROPT. 2448C 2449#include "implicit.h" 2450 CHARACTER*(*) KEY 2451C 2452 PARAMETER ( NTYP = 7 ) 2453 DIMENSION WRK(7,2,NTYP) 2454C 2455 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 ) 2456C 2457C 2458C Used from common blocks: 2459C INFTIM : NCALLS,TIMCPU,TIMWAL 2460C 2461#include "inftim.h" 2462C 2463C Local variables: 2464C 2465 INTEGER INAME(NTYP),IEND(NTYP) 2466 CHARACTER*8 NAME(7,NTYP) 2467 DATA INAME /1,2,6,7,10,8,3/, IEND /7,4,5,5,5,5,7/ 2468 DATA NAME /' GRAD',' MAKDM',' FCKMAT', 2469 * ' GETH2',' ORDIAG',' ORBGRD',' CIGRAD', 2470 * ' CISIGC',' diag',' ONEELC', 2471 * ' TWOELC',' ',' ',' ', 2472 * ' LINTRN',' SIRTR1',' ORBSIG', 2473 * ' CISIG',' SOLLIN',' ',' ', 2474 * ' MAKTDM',' "diag"',' "oneel"', 2475 * ' "twoel"','"sym PV"',' ',' ', 2476 * ' CISIGO',' DIAG1X',' ONEELX', 2477 * ' TWOELX',' other',' ',' ', 2478 * ' ORBLIN',' SIRTR1',' ORBSIG', 2479 * ' CISIG',' SOLLIN',' ',' ', 2480 * ' UPDGRAD',' MAKDM',' FCKMAT', 2481 * ' GETH2',' ORDIAG',' ORBGRD',' CIGRAD'/ 2482C 2483C 2484 IF (KEY(1:3) .EQ. 'STA') THEN 2485 DO 140 ITYP = 1,20 2486 NCALLS(ITYP) = 0 2487 NVECS (ITYP) = 0 2488 DO 120 IBLCK = 1,7 2489 TIMCPU(IBLCK,ITYP) = D0 2490 TIMWAL(IBLCK,ITYP) = D0 2491 120 CONTINUE 2492 140 CONTINUE 2493 RETURN 2494 END IF 2495C 2496 WRITE (LUPRI,10) 2497 WRITE (LUPRI,15) 2498C 2499 DO 400 I = 1,NTYP 2500 II = INAME(I) 2501 JEND = IEND(I) 2502 MCALLS = NCALLS(II) 2503 IF (MCALLS .EQ. 0) THEN 2504 WRITE (LUPRI,25) NAME(1,I) 2505 WRITE (LUPRI,15) 2506 GO TO 400 2507 END IF 2508 IF (NVECS(II) .GT. 0) THEN 2509 FAC = NVECS(II) 2510 ELSE 2511 FAC = MCALLS 2512 END IF 2513 FAC = D1 / FAC 2514 DO 300 J = 1,JEND 2515 WRK(J,1,I) = FAC * TIMCPU(J,II) 2516 WRK(J,2,I) = FAC * TIMWAL(J,II) 2517 300 CONTINUE 2518C 2519 IF (NVECS(II) .LE. 0) THEN 2520 WRITE (LUPRI,20) NAME(1,I),MCALLS,TIMCPU(1,II),WRK(1,1,I), 2521 * TIMWAL(1,II),WRK(1,2,I) 2522 ELSE 2523 WRITE (LUPRI,21) NAME(1,I),MCALLS,NVECS(II), 2524 * TIMCPU(1,II),WRK(1,1,I), 2525 * TIMWAL(1,II),WRK(1,2,I) 2526 END IF 2527 IF (JEND.GT.1) 2528 * WRITE (LUPRI,30) (NAME(J,I),TIMCPU(J,II),WRK(J,1,I), 2529 * TIMWAL(J,II),WRK(J,2,I),J=2,JEND) 2530 WRITE (LUPRI,15) 2531 400 CONTINUE 2532C 2533 RETURN 2534C 2535 10 FORMAT(//' (TIMOPT) Timing of MC optimization:', 2536 * //' Name . # calls . total CPU . CPU / call .', 2537 * ' total wall . wall / call .') 2538 15 FORMAT(1X,78('-')) 2539 20 FORMAT(/1X,A8,I10,4F15.2) 2540 21 FORMAT(/1X,A8,I10/' # vecs',I10,4F15.2) 2541 25 FORMAT(/1X,A8,' has not been called in this calculation.') 2542 30 FORMAT(/1X,A8,10X,4F15.2) 2543 END 2544C /* Deck wr_sirifc */ 2545 SUBROUTINE WR_SIRIFC(CMO,DV,PV,FC,FV,FQ,CREF,FCAC,H2AC,WRK,LFREE, 2546 & GORB,ERLM,ORBHES,XINDX) 2547C 2548C MOTECC: this routine writes the interface file SIRIFC for post-programs. 2549C The structure of the SIRIFC file is described here. 2550C 2551C Written 14-Feb-1985 Hans Joergen Aa. Jensen 2552C Revisions: 2553C 21-Nov-1985 hjaaj (for RHF) 2554C 8-Feb-1987 hjaaj (full Fock matrix) 2555C Oct-1988 hjaaj (for non-guga mcscf) 2556C 6-Sep-1989 hjaaj (corrected FCDIA calculation) 2557C 2-Sep-1999 hjaaj (do not use FV for NASHT.eq.0) 2558C 2559C Purpose: 2560C Write information to LUSIFC needed for (1) Trygve Helgaker's 2561C first and second order geometry derivatives program and/or 2562C (2) response calculation. 2563C 2564C Suggestions: 2565C 890906-hjaaj: save parameter telling if orbitals are canonical HF 2566C orbitals or natural orbitals or ? 2567C 2568C The following records are written: 2569C 2570C 0) label LBSIFC ("SIR IPH ") 2571C 1) POTNUC,EMY,EACTIV,EMCSCF,ISTATE,ISPIN,NACTEL,LSYM,MS2 2572C 2) NISHT,NASHT,... 2573C 3) CMO 2574C 4) CREF 2575C 5) DV 2576C 6) FOCK 2577C 7) PV 2578C 8) FC 2579C 9) FV 2580C 10) FCAC 2581C 11) H2AC 2582C 12) GORB 2583C If (GUGA) then 2584C 13) label "CIDIAG1 " 2585C 14) CIDIAG1 2586C end if 2587C 15) label "CIDIAG2 " 2588C 16) CIDIAG2 2589C 17) label "ORBDIAG " 2590C 18) ORBDIAG 2591C 2592C *) label "SIRFLAGS" 2593C *) FLAG(1:NFLAG) 2594C *) GRDNRM,POTNUC,EMY,EACTIV,ESOLT,EMCSCF 2595C If (fields) then 2596C *) label "EXTFIELD" 2597C *) NFIELD 2598C *) (EFIELD(I),I=1,NFIEL4) 2599C *) (LFIELD(I),I=1,NFIEL4) 2600C end if 2601C If (solvent) then 2602C *) label "SOLVINFO" 2603C *) EPSOL,EPSTAT,EPPN,RSOL(1:3),LSOLMX,INERSI,INERSF 2604C *) GRDNRM,POTNUC,EMY,EACTIV,ESOLT,EMCSCF 2605C *) ERLM(LM,1), LM = 1,NLMSOL) 2606C *) (TRLM(LM), LM = 1,NLMSOL) where TRLM(i) = ERLM(i,2) 2607C *) NSYM, NBAS 2608C *) label "SOLVTMAT" 2609C *) TMAT(1:NNORBX) 2610C end if 2611C If (CSF's) then 2612C *) label "CREFDETS" 2613C *) CREF in dets 2614C end if 2615C *) label 'TRCCINT" 2616C *) NSYM, NORBT, ... 2617C *) FCDIA, ISMO 2618C *) CMO 2619C 2620 use pelib_interface, only: use_pelib, pelib_ifc_fock 2621#if defined (HAS_PCMSOLVER) 2622 use pcm_config, only: pcm_configuration, pcm_cfg 2623 use pcm_scf 2624#endif 2625#include "implicit.h" 2626 DIMENSION CMO(*),DV(*),PV(*),FC(*),FV(*),FQ(NORBT,*),XINDX(*) 2627 DIMENSION ERLM(NLMSOL,2),GORB(*),CREF(*),FCAC(*),H2AC(*),WRK(*) 2628C 2629 PARAMETER ( D0=0.0D0, D2=2.0D0 ) 2630#include "dummy.h" 2631C 2632C Used from common blocks: 2633C INFINP : ABAIPH,ISTATE,ISPIN,NACTEL,FLAG(*),INERSI 2634C SPINFO : MS2 2635C INFOPT : EMY, EACTIV, EMCSCF, ESOLT 2636C INFVAR : NCONF,NWOPT 2637C INFORB : NSYM,NISHT,NASHT,NNASHX,NNASHY,NOCCT,NNOCCX,NCMOT,... 2638C INFIND : IROW(*) 2639C INFTAP : LUSIFC 2640C INFDIM : NASHDI,NASHMA,NORBMA 2641C INFPRI : ? 2642C CBGETDIS : IADH2 2643C CBIREA : LMULBS 2644C R12INT : R12CAL 2645C dftcom.h : DFT_SPINDNS 2646C 2647#include "maxash.h" 2648#include "maxorb.h" 2649#include "mxcent.h" 2650#include "pcmdef.h" 2651#include "priunit.h" 2652#include "pcm.h" 2653#include "pcmlog.h" 2654#include "infinp.h" 2655#include "spinfo.h" 2656#include "infopt.h" 2657#include "infvar.h" 2658#include "inforb.h" 2659#include "infind.h" 2660#include "inftap.h" 2661#include "infdim.h" 2662#include "infpri.h" 2663#include "cbgetdis.h" 2664#include "cbirea.h" 2665#include "r12int.h" 2666#include "qmmm.h" 2667#include "qm3.h" 2668#include "gnrinf.h" 2669#include "dftcom.h" 2670C 2671 LOGICAL FNDLAB, ORBHES, LGLOCAL, EXP1VL, TOFILE 2672C 2673 CHARACTER*8 LAB123(3), TABLE(12), LABEOD, LABCC, LABABA 2674 CHARACTER*8 LABAUX, LABUNI 2675 DATA LAB123/'********','********','********'/ 2676 DATA TABLE /'CIDIAG1 ','CIDIAG2 ','ORBDIAG ','SIRFLAGS', 2677 * 'EXTFIELD','SOLVINFO','SOLVTMAT','CREFDETS', 2678 * 'PCMINFO ','PCMJXMAT','QMMMTMAT','PEFMAT '/ 2679 DATA LABEOD, LABCC /'EODATA ', 'TRCCINT '/ 2680 DATA LABUNI, LABAUX /'FULLBAS ', 'AUXILBAS'/ 2681#if defined (HAS_PCMSOLVER) 2682 real(8), allocatable :: fock_pcm_mo(:) 2683#endif 2684C 2685 CALL QENTER('WR_SIRIFC') 2686 LGLOCAL = BOYORB.OR.PIPORB 2687 CALL GETDAT(LAB123(2),LAB123(3)) 2688C ... place date in LAB123(2) and time in LAB123(3) 2689 2690C 2691C *** Now create LUSIFC and write ... 2692C 2693 IF (LUSIFC .LE. 0) CALL GPOPEN(LUSIFC,FNSIFC, 2694 & 'UNKNOWN',' ','UNFORMATTED',IDUMMY,.FALSE.) 2695 REWIND (LUSIFC) 2696 2697 IF (IPRI6 .GE. 5) WRITE (LUPRI,'(/A)') 2698 * ' --- Writing Sirius interface file for'// 2699 * ' post-programs in WR_SIRIFC. ---' 2700C 2701C Calculate Fock matrix 2702C 2703 KFCDIA = 1 2704 KFOCK = KFCDIA + NORBT 2705 KWRK1 = KFOCK + N2ORBT 2706 KFCI = KWRK1 2707 KUDV = KFCI + NORBMA*NORBMA 2708 KFCFV = KUDV + N2ASHX 2709 KEND = KFCFV + NNORBT 2710 KEND = MAX(KEND,KWRK1+NCONF,KWRK1+NWOPT) 2711 IF (KEND .GT. LFREE) CALL ERRWRK('WR_SIRIFC',KEND,LFREE) 2712C 2713 IF (NASHT .EQ. 0) THEN 2714 DO I = 1,NNORBT 2715 WRK((KFCFV-1)+I) = D2*FC(I) 2716 END DO 2717 ELSE 2718 CALL DSPTSI(NASHT,DV,WRK(KUDV)) 2719 DO I = 1,NNORBT 2720 WRK((KFCFV-1)+I) = D2*(FC(I) + FV(I)) 2721 END DO 2722 END IF 2723C 2724 CALL DZERO(WRK(KFOCK),N2ORBT) 2725 DO 500 ISYM = 1,NSYM 2726 NORBI = NORB(ISYM) 2727 IORBI = IORB(ISYM) 2728C 2729 JFCDIA = KFCDIA - 1 + IORBI 2730 DO I = 1,NORBI 2731 WRK(JFCDIA+I) = FC( IIORB(ISYM) + IROW(I+1) ) 2732 END DO 2733C 2734 JFOCK = KFOCK + I2ORB(ISYM) 2735 NOCCI = NOCC(ISYM) 2736 IF (NOCCI .EQ. 0) THEN 2737 IF (NORBI .GT. 0) CALL DZERO(WRK(JFOCK),N2ORB(ISYM)) 2738 ELSE 2739 NISHI = NISH(ISYM) 2740 NASHI = NASH(ISYM) 2741 NSSHI = NSSH(ISYM) 2742C inactive-general + junk in active-general and secondary-general 2743 IF (NISHI .GT. 0) THEN 2744 JFCFV = KFCFV + IIORB(ISYM) 2745 CALL DSPTSI(NORBI,WRK(JFCFV),WRK(JFOCK)) 2746 END IF 2747C active-general 2748 IF (NASHI .GT. 0) THEN 2749 CALL DSPTSI(NORBI,FC(1+IIORB(ISYM)),WRK(KFCI)) 2750 IASHI = IASH(ISYM) 2751 JUDV = KUDV + IASHI + IASHI*NASHT 2752 JFCI = KFCI + NISHI 2753 IFOCKJ = JFOCK + NISHI 2754 CALL DGEMM('N','N',NASHI,NORBI,NASHI,1.D0, 2755 & WRK(JUDV),NASHT, 2756 & WRK(JFCI),NORBI,0.D0, 2757 & WRK(IFOCKJ),NORBI) 2758 DO 360 J = 1,NORBI 2759 IFOCKJ = JFOCK - 1 + NISHI + (J-1)*NORBI 2760 DO 340 NI = 1,NASHI 2761 WRK(IFOCKJ + NI) = WRK(IFOCKJ + NI) 2762 & + FQ(IORBI + J, IASHI + NI) 2763 340 CONTINUE 2764 360 CONTINUE 2765 END IF 2766C secondary-general is zero 2767 IF (NSSHI .GT. 0) THEN 2768 DO 460 J = 1,NORBI 2769 IFOCKJ = JFOCK - 1 + NOCCI + (J-1)*NORBI 2770 DO 440 NI = 1,NSSHI 2771 WRK(IFOCKJ + NI) = D0 2772 440 CONTINUE 2773 460 CONTINUE 2774 END IF 2775 END IF 2776C 2777 500 CONTINUE 2778C 2779 IF (P6FLAG(14)) THEN 2780 WRITE (LUPRI,'(/A)') ' Full Fock matrix from WR_SIRIFC :' 2781 CALL OUTPTB(WRK(KFOCK),NORB,NSYM,1,LUPRI) 2782 WRITE (LUPRI,'(/A)') ' FC diagonal elememts from WR_SIRIFC :' 2783 WRITE (LUPRI,'(/,(1X,5F15.8))') (WRK(KFCDIA-1+I),I=1,NORBT) 2784 END IF 2785C 2786C 2787 WRITE (LUSIFC) LAB123,LBSIFC 2788 WRITE (LUSIFC) POTNUC,EMY,EACTIV,EMCSCF, 2789 & ISTATE,ISPIN,NACTEL,LSYM,MS2 2790 IF (MCTYPE .EQ. 0 .AND. HSROHF) THEN 2791 MCTYPE_IFC = -1 2792 ELSE 2793 MCTYPE_IFC = MCTYPE 2794 END IF 2795 IF (DFT_SPINDNS) MCTYPE_IFC = MCTYPE_IFC + 1000 2796 WRITE (LUSIFC) NISHT,NASHT,NOCCT,NORBT,NBAST,NCONF,NWOPT,NWOPH, 2797 & NCDETS, NCMOT,NNASHX,NNASHY,NNORBT,N2ORBT, 2798 & NSYM,MULD2H, NRHF,NFRO, 2799 & NISH,NASH,NORB,NBAS, 2800 & NELMN1, NELMX1, NELMN3, NELMX3, MCTYPE_IFC, 2801 & NAS1, NAS2, NAS3 2802C 2803C 880920-hjaaj: later write label here for orbitals 2804C 2805 IF (DFT_SPINDNS) THEN 2806 NDV = 2 ! DV, DS; singlet and triplet active density matrices 2807 ELSE 2808 NDV = 1 2809 END IF 2810 NC4 = MAX(4,NCONF) 2811 NW4 = MAX(4,NWOPT) 2812 NWH4 = MAX(4,NWOPH) 2813 NCMOT4 = MAX(4,NCMOT) 2814 MMORBT = MAX(4,NNORBT) 2815 M2ORBT = MAX(4,N2ORBT) 2816 MMASHX = MAX(4,NNASHX) 2817 MMASHY = MAX(4,NNASHY) 2818 M2ASHY = MAX(4,NNASHX*NNASHX) 2819C 2820 CALL WRITT (LUSIFC,NCMOT4,CMO) ! rec no. 3 2821 CALL WRITT (LUSIFC,NC4,CREF) ! rec no. 4 2822 CALL WRITT (LUSIFC,NDV*MMASHX,DV) ! rec no. 5; DV, (DS) 2823C Write Fock matrix. 2824 CALL WRITT (LUSIFC,M2ORBT,WRK(KFOCK)) ! rec no. 6 2825 CALL WRITT (LUSIFC,M2ASHY,PV) ! rec no. 7 2826 CALL WRITT (LUSIFC,MMORBT,FC) ! rec no. 8 2827 IF (NASHT .EQ. 0) THEN 2828 CALL DZERO(WRK(KFOCK),NNORBT) 2829 CALL WRITT (LUSIFC,MMORBT,WRK(KFOCK)) 2830 ELSE 2831 CALL WRITT (LUSIFC,NDV*MMORBT,FV) ! rec. no. 9; FC, (FS) 2832 END IF 2833C 2834 CALL WRITT (LUSIFC,NDV*MMASHX,FCAC) ! rec. no. 10; FCAC, (FSAC) 2835 2836 IF (IADH2 .GE. 0) THEN 2837#if defined (VAR_NEWCODE) 2838... 27-Jun-1987/hjaaj -- noget lignende : MAERKE 2839 DO 100 IJ = 1,NNASHX 2840 CALL READDI(LUDA2,IADH2+IJ,IRAT*NNASHX,H2AC) 2841 CALL WRITT (LUSIFC,MMASHX,H2AC) 2842 100 CONTINUE 2843... men der mangler en label som fortaeller response at H2AC on disk. 2844#else 2845 WRITE (LUPRI,*) 2846 & 'SIRCI.WR_SIRIFC: ".DISKH2" not implemented here yet.' 2847 CALL QTRACE(LUPRI) 2848 CALL QUIT('SIRCI.WR_SIRIFC: ".DISKH2" not implemented here.') 2849#endif 2850 ELSE 2851 CALL WRITT (LUSIFC,(MMASHX*MMASHX),H2AC) ! rec no. 11 2852 END IF 2853 CALL WRITT (LUSIFC,NWH4,GORB) ! rec no. 12 2854C 2855C *** now write diagonals of L-matrix needed for LINTRN 2856C and for conjugate gradient "next" algorithm 2857C 2858 IF (NCONF .GT. 1) THEN 2859 REWIND LUIT2 2860 IF ( FNDLAB(TABLE(1),LUIT2) ) THEN 2861C ... CIDIAG1 is only used if CASGUGA 2862 CALL READT (LUIT2,NCONF,WRK(KWRK1)) 2863 WRITE (LUSIFC) LAB123,TABLE(1) 2864 CALL WRITT (LUSIFC,NC4,WRK(KWRK1)) 2865 ELSE 2866 REWIND LUIT2 2867 END IF 2868C 2869 CALL MOLLAB(TABLE(2),LUIT2,LUPRI) 2870 CALL READT (LUIT2,NCONF,WRK(KWRK1)) 2871 WRITE (LUSIFC) LAB123,TABLE(2) 2872 CALL WRITT (LUSIFC,NC4,WRK(KWRK1)) 2873 END IF 2874C 2875 IF (NWOPT .GT. 0 .AND. ORBHES) THEN 2876 REWIND LUIT2 2877 CALL MOLLAB(TABLE(3),LUIT2,LUPRI) 2878 CALL READT (LUIT2,NWOPT,WRK(KWRK1)) 2879 WRITE (LUSIFC) LAB123,TABLE(3) 2880 CALL WRITT (LUSIFC,NW4,WRK(KWRK1)) 2881 END IF 2882 WRITE (LUSIFC) LAB123,TABLE(4) 2883 WRITE (LUSIFC) (FLAG(I),I=1,NFLAG) 2884 WRITE (LUSIFC) GRDNRM,POTNUC,EMY,EACTIV,ESOLT,EMCSCF 2885C Write External fields, if any 2886 IF (NFIELD .GT. 0) THEN 2887 WRITE (LUSIFC) LAB123,TABLE(5) 2888 WRITE (LUSIFC) NFIELD,DUMMY,DUMMY,DUMMY,DUMMY 2889 NFIEL4 = MAX(4,NFIELD) 2890 WRITE (LUSIFC) (EFIELD(I),I=1,NFIEL4) 2891 WRITE (LUSIFC) (LFIELD(I),I=1,NFIEL4) 2892 END IF 2893 2894C Write Solvent information, if Solvent calculation. 2895 IF (FLAG(16)) THEN 2896 KTMAT = KWRK1 2897 KWRK2 = KTMAT + NNORBT 2898 LWRK2 = LFREE - KWRK2 + 1 2899 IF (KWRK2 .GT. LFREE) 2900 & CALL ERRWRK('WR_SIRIFC-SOLVENT',KWRK2,LFREE) 2901 WRITE (LUSIFC) LAB123,TABLE(6) 2902 WRITE (LUSIFC) EPSOL,EPSTAT,EPPN,RSOL,LSOLMX,INERSI,INERSF 2903 WRITE (LUSIFC) GRDNRM,POTNUC,EMY,EACTIV,ESOLT,EMCSCF 2904 NLMSO4 = MAX(4,NLMSOL) 2905 CALL WRITT (LUSIFC,NLMSO4,ERLM(1,1)) 2906 CALL WRITT (LUSIFC,NLMSO4,ERLM(1,2)) 2907 WRITE (LUSIFC) NSYM, NBAS 2908 CALL SOLFCK(DUMMY,DUMMY,WRK(KTMAT),CMO,ERLM(1,1),.TRUE., 2909 & DUMSOL,WRK(KWRK2),LWRK2,IPRSOL) 2910 WRITE (LUSIFC) LAB123,TABLE(7) 2911 CALL WRITT (LUSIFC,MMORBT,WRK(KTMAT)) 2912 END IF 2913C-------------- 2914C CBN+JK, 03.01.06 2915C------------------- 2916C Write qm3 solvent part of the Fock matrix to interface file 2917 2918 IF (QM3) THEN 2919 KTMAT = KWRK1 2920 KFSOLAO = KTMAT + NNORBT 2921 KUCMO = KFSOLAO + NNBASX 2922 KWRK2 = KUCMO + NORBT*NBAST 2923 LWRK2 = LFREE - KWRK2 + 1 2924 IF (KWRK2 .GT. LFREE) CALL ERRWRK('WRSIFC-2',KWRK2,LFREE) 2925 CALL QM3_FMO(WRK(KFSOLAO),WRK(KWRK2),LWRK2,IQM3PR) 2926 CALL UPKCMO(CMO,WRK(KUCMO)) 2927 CALL UTHU(WRK(KFSOLAO),WRK(KTMAT),WRK(KUCMO),WRK(KWRK2), 2928 & NBAST,NORBT) 2929 WRITE (LUSIFC) LAB123,TABLE(7) 2930 CALL WRITT (LUSIFC,MMORBT,WRK(KTMAT)) 2931 END IF 2932C-------------- 2933C CBN+JK, 03.01.06 2934C------------------- 2935C Write PCM solvent info if PCM calculation 2936 IF (PCM) THEN 2937 KTMAT = 1 2938 KDCAO = KTMAT + NNORBT 2939 KDVAO = KDCAO + N2BASX 2940 KPOT = KDVAO + N2BASX 2941 KWRK2 = KPOT + NTS 2942 LWRK2 = LFREE - KWRK2 + 1 2943 IF (KWRK2 .GT. LFREE) CALL ERRWRK('WR_SIRIFC-PCM',KWRK2,LFREE) 2944C 2945 CALL PCMFCK(DUMMY,DUMMY,WRK(KTMAT),CMO,.TRUE., 2946 * DUMSOL,WRK(KWRK2),LWRK2,IPRSOL) 2947 WRITE (LUSIFC) LAB123,TABLE(10) 2948 CALL WRITT (LUSIFC,MMORBT,WRK(KTMAT)) 2949 WRITE (LUSIFC) LAB123,TABLE(9) 2950 WRITE (LUSIFC) NTS 2951 WRITE (LUSIFC) (QSE(ITS),ITS = 1, NTS) 2952 WRITE (LUSIFC) (QSN(ITS),ITS = 1, NTS) 2953 WRITE (LUSIFC) (-QSE(ITS)-QSN(ITS),ITS = 1, NTS) 2954C 2955Clf WRITE (LUSIFC) (QSE(ITS)-QSN(ITS),ITS = 1, NTS) 2956Clf 2957C 2958C Calculate expectation value of the electronic potential 2959C at the tessera 2960C 2961 EXP1VL = .TRUE. 2962 TOFILE = .FALSE. 2963 CALL DZERO(WRK(KDCAO),2*N2BASX) 2964 CALL FCKDEN((NISHT.GT.0),(NASHT.GT.0),WRK(KDCAO), 2965 & WRK(KDVAO),CMO,DV,WRK(KWRK2),LWRK2) 2966 CALL DGEFSP(NBAST,WRK(KDCAO),WRK(KWRK2)) 2967 CALL DCOPY(NNORBX,WRK(KWRK2),1,WRK(KDCAO),1) 2968ckr CALL PKSYM1(WRK(KWRK2),WRK(KDCAO),NBAS,NSYM,1) 2969 IF (NASHT .GT. 0) THEN 2970 CALL DGEFSP(NBAST,WRK(KDVAO),WRK(KWRK2)) 2971 CALL DCOPY(NNORBX,WRK(KWRK2),1,WRK(KDVAO),1) 2972ckr CALL PKSYM1(WRK(KWRK2),WRK(KDVAO),NBAS,NSYM,1) 2973 END IF 2974Ckr 2975Ckr Check whether DV+DC is good enough for open-shell systems (should 2976Ckr be the case). 2977Ckr 2978 CALL DAXPY(NNBASX,1.0D0,WRK(KDVAO),1,WRK(KDCAO),1) 2979 CALL J1INT(WRK(KPOT),EXP1VL,WRK(KDCAO),1,TOFILE,'NPETES ', 2980 & 1,WRK(KWRK2),LWRK2) 2981 WRITE (LUSIFC) (WRK(KPOT + I - 1), I = 1, NTS) 2982 IF (NEQRSP) CALL V2Q(WRK(KWRK2),WRK(KPOT),QSENEQ,QETNEQ,NEQRSP) 2983 END IF 2984 2985#if defined (HAS_PCMSOLVER) 2986! Write PCM contribution to Fock matrix in MO basis to SIRIFC 2987 if(pcm_cfg%do_pcm) then 2988 kwrk2 = 1 2989 lwrk2 = lfree - kwrk2 +1 2990! Allocate scratch space for the PCM Fock matrix contribution in AO basis. 2991 allocate(fock_pcm_mo(nnorbt)) 2992 fock_pcm_mo = 0.0d0 2993! Calculate PCM contribution to Fock matrix in MO basis 2994 call pcm_mo_fock(fock_pcm_mo, cmo, wrk(kwrk2), lwrk2) 2995 WRITE (LUSIFC) LAB123,'EXTPCMMAT' 2996 CALL WRITT (LUSIFC,MMORBT,fock_pcm_mo) 2997 end if 2998#endif 2999 3000C ****** edh: 24/11/2011 ****** 3001 3002C Write qm/mm solvent part of the Fock matrix to interface file 3003 3004 IF (QMMM) THEN 3005 3006 IF (IPQMMM .GE. 15) THEN 3007 WRITE(LUPRI,*)'WRITING QMMM FOCK MAT TO SIRIUS INTERFACE FILE' 3008 ENDIF 3009 3010 KDCAO = KWRK1 3011 KDVAO = KDCAO + N2BASX 3012 KDENSF = KDVAO + N2BASX 3013 KMAT = KDENSF + NNBASX 3014 KUCMO = KMAT + NNBASX 3015 KMATMO = KUCMO + NORBT*NBAST 3016 KWRK2 = KMATMO + NNORBT 3017 LWRK2 = LFREE - KWRK2 + 1 3018 3019 IF (LWRK2 .LT. 0) CALL ERRWRK('WR_SIRIFC-QMMM',-KWRK2,LFREE) 3020 3021 CALL DZERO(WRK(KDCAO),N2BASX) 3022 CALL DZERO(WRK(KDVAO),N2BASX) 3023 CALL DZERO(WRK(KDENSF),NNBASX) 3024 CALL DZERO(WRK(KMAT),NNBASX) 3025 CALL DZERO(WRK(KUCMO),NORBT*NBAST) 3026 CALL DZERO(WRK(KMATMO),NNORBT) 3027 3028 CALL FCKDEN((NISHT.GT.0),(NASHT.GT.0),WRK(KDCAO),WRK(KDVAO), 3029 & CMO,DV,WRK(KWRK2),LWRK2) 3030 IF (NASHT .GT. 0) THEN 3031 CALL DAXPY(N2BASX,1.0D0,WRK(KDVAO),1,WRK(KDCAO),1) 3032 ENDIF 3033 3034 CALL DGEFSP(NBAST,WRK(KDCAO),WRK(KDENSF)) 3035 3036 CALL QMMM_FCK_AO(WRK(KMAT),WRK(KDENSF),EQMMM,WRK(KWRK2),LWRK2, 3037 & IPQMMM) 3038 CALL UPKCMO(CMO,WRK(KUCMO)) 3039 CALL UTHU(WRK(KMAT),WRK(KMATMO),WRK(KUCMO),WRK(KWRK2), 3040 & NBAST,NORBT) 3041 3042 IF (IPQMMM .GE. 15) THEN 3043 WRITE (LUPRI,'(/A)') ' QMMM_ao matrix from WR_SIRIFC:' 3044 CALL OUTPAK(WRK(KMAT),NBAST,1,LUPRI) 3045 WRITE (LUPRI,'(/A)') ' QMMM_mo matrix from WR_SIRIFC' 3046 CALL OUTPAK(WRK(KMATMO),NORBT,1,LUPRI) 3047 ENDIF 3048 3049 WRITE (LUSIFC) LAB123,TABLE(11) 3050 CALL WRITT (LUSIFC,MMORBT,WRK(KMATMO)) 3051 END IF 3052 3053 IF (USE_PELIB()) THEN 3054 ! Calculate ground-state PE contribution which is subtracted from 3055 ! the full Fock matrix in response for MCSCF and CI wavefunction 3056 KFCKMO = KWRK1 3057 KENR = KFCKMO + NNORBX 3058 KDCAO = KENR + 1 3059 KDVAO = KDCAO + N2BASX 3060 KFDTAO = KDVAO + N2BASX 3061 KFCKAO = KFDTAO + NNBASX 3062 KWRK2 = KFCKAO + NNBASX 3063 LWRK2 = LFREE - KWRK2 + 1 3064 IF (LWRK2 > LFREE) CALL ERRWRK('WR_SIRIFC-PELIB',-KWRK2,LFREE) 3065 CALL FCKDEN((NISHT>0), (NASHT>0), WRK(KDCAO), WRK(KDVAO), 3066 & CMO, DV, WRK(KWRK2), LWRK2) 3067 IF (NISHT == 0) THEN 3068 CALL DZERO(WRK(KDCAO), N2BASX) 3069 END IF 3070 IF (NASHT > 0) THEN 3071 CALL DAXPY(N2BASX, 1.0D0, WRK(KDVAO), 1, WRK(KDCAO), 1) 3072 END IF 3073 CALL DGEFSP(NBAST, WRK(KDCAO), WRK(KFDTAO)) 3074 CALL PELIB_IFC_FOCK(WRK(KFDTAO), WRK(KFCKAO), WRK(KENR)) 3075 CALL AROUND('Writing FOCKMAT in siropt') 3076 IF (.NOT.HFFLD) CALL PUT_TO_FILE('FOCKMAT',NNBASX, 3077 & WRK(KFCKAO)) 3078 CALL PUT_TO_FILE('FOCKMHF',NNBASX,WRK(KFCKAO)) 3079 CALL UTHU(WRK(KFCKAO), WRK(KFCKMO), CMO, WRK(KWRK2), 3080 & NBAST, NORBT) 3081 WRITE (LUSIFC) LAB123, TABLE(12) 3082 CALL WRITT(LUSIFC, NNORBX, WRK(KFCKMO)) 3083 END IF 3084C 3085C Revision 2000/05/24 12:43:25 hjj 3086C Write CREF in determinants here, if NCDETS .gt. NCONF 3087C 3088 IF (NCDETS .GT. NCONF) THEN 3089 IF (FLAG(27)) THEN 3090 WRITE (LUPRI,*) 3091 & 'WR_SIRIFC ERROR, .DETERMINANTS but NCDETS.gt.NCONF:', 3092 & NCDETS, NCONF 3093 CALL QUIT('WR_SIRIFC: NCDETS.gt. NCONF for .DETERMINANTS') 3094 END IF 3095 WRITE (LUSIFC) LAB123,TABLE(8) ! label "CREFDETS" 3096 KWRK1 = 1 + NCDETS 3097 LWRK1 = LFREE - KWRK1 3098 CALL GETDETS(LSYM,NCONF,NCDETS,XINDX,CREF,WRK,WRK(KWRK1),LWRK1) 3099 NC4 = MAX(4,NCDETS) 3100 CALL WRITT(LUSIFC,NC4,WRK) 3101 END IF 3102C 3103C Interface for CC modules 3104C 3105 IF (LMULBS) THEN 3106 IOFF1 = 1 3107 IOFF2 = 1 + NORBT 3108 NBAST1 = 0 3109 NCMOT1 = 0 3110 DO ISYM = 1, NSYM 3111 NBAST1 = NBAST1 + MBAS1(ISYM) 3112 NCMOT1 = NCMOT1 + MBAS1(ISYM) * NORB(ISYM) 3113 DO I = 1, NORB(ISYM) 3114 CALL DCOPY(MBAS1(ISYM),CMO(IOFF1),1,WRK(IOFF2),1) 3115 IOFF1 = IOFF1 + NBAS(ISYM) 3116 IOFF2 = IOFF2 + MBAS1(ISYM) 3117 ENDDO 3118 ENDDO 3119 NCMOTX = MAX(4,NCMOT1) 3120 WRITE (LUSIFC) LAB123, LABCC ! label "TRCCINT " 3121 WRITE (LUSIFC) NSYM,NORBT,NBAST1,NCMOT1,(NOCC(I),I=1,NSYM), 3122 * (NORB(I),I=1,NSYM),(MBAS1(I),I=1,NSYM),POTNUC,EMCSCF 3123 WRITE (LUSIFC) (WRK(KFCDIA-1+I),I=1,NORBT),(ISMO(I),I=1,NORBT), 3124 * DUMMY,DUMMY,DUMMY 3125 WRITE (LUSIFC) (WRK(I),I=1+NORBT,NCMOTX+NORBT) 3126 ELSE 3127 WRITE (LUSIFC) LAB123, LABCC 3128 WRITE (LUSIFC) NSYM,NORBT,NBAST,NCMOT,(NOCC(I),I=1,NSYM), 3129 * (NORB(I),I=1,NSYM),(NBAS(I),I=1,NSYM),POTNUC,EMCSCF 3130 WRITE (LUSIFC) (WRK(KFCDIA-1+I),I=1,NORBT),(ISMO(I),I=1,NORBT), 3131 * DUMMY,DUMMY,DUMMY 3132 WRITE (LUSIFC) (CMO(I),I=1,NCMOT4) 3133 END IF 3134C 3135 IF (LMULBS .OR. R12CAL) THEN 3136C Interface for auxiliary basis set (WK/UniKA/04-11-2002). 3137 DO ISYM = 1, NSYM 3138 IF (NRXR12(ISYM).GT.(NORB1(ISYM)-NOCC(ISYM))) THEN 3139 WRITE(LUPRI,*) 'You specified more virtual R12 orbitals '// 3140 & 'than there are virtual orbitals in this '// 3141 & 'symmery:' 3142 WRITE(LUPRI,*) 'NVIR(',ISYM,') = ',NORB1(ISYM)-NOCC(ISYM) 3143 WRITE(LUPRI,*) 'NRXR12(',ISYM,') = ',NRXR12(ISYM) 3144 CALL QUIT('Too many virtual R12 orbitals') 3145 ELSE 3146 !redefine NOCC 3147 NOCC(ISYM) = NOCC(ISYM) + NRXR12(ISYM) 3148 END IF 3149 END DO 3150 LUMULB = 34 3151 CALL GPOPEN(LUMULB,'AUXBAS','UNKNOWN','SEQUENTIAL', 3152 & 'FORMATTED',IDUMMY,LDUMMY) 3153 READ (LUMULB,*) NDUMM 3154 NAUXT = 0 3155 MOAUX = 0 3156 NFULLT = 0 3157 NMOFLT = 0 3158 NOCCT = 0 3159 NORBFT = 0 3160 DO ISYM = 1, NSYM 3161 NOCCT = NOCCT + NOCC(ISYM) 3162 NORBFT = NORBFT + NORB1(ISYM) + NORB2(ISYM) 3163 END DO 3164 NBASC = NBAST + NOCCT 3165 NORBFT = NORBFT + NOCCT 3166 IOFFEV = KFCDIA + NBASC 3167 IOFFMO = IOFFEV + NBASC * 2 3168 DO ISYM = 1, NSYM 3169 READ (LUMULB,*) IDUMM,NORBI,NAUXI,NBASI 3170 NAUXT = NAUXT + NAUXI 3171 MOAUX = MOAUX + NAUXI * NBASI 3172 NFULLT = NFULLT + NBASI 3173 NMOFLT = NMOFLT + NBASI * (NORBI + NAUXI + NOCC(ISYM)) 3174 DO I = 1, NAUXI 3175 READ (LUMULB,*) NDUMM 3176 READ (LUMULB,'(4E30.20)') WRK(IOFFEV) 3177 IOFFEV = IOFFEV + 1 3178 READ (LUMULB,'(4E30.20)') (WRK(IOFFMO+K), K = 0, NBASI-1) 3179 IOFFMO = IOFFMO + NBASI 3180 END DO 3181 END DO 3182 CALL GPCLOSE(LUMULB,'KEEP') 3183 WRITE (LUSIFC) LAB123, LABAUX ! label "AUXILBAS" 3184 WRITE (LUSIFC) NSYM,NAUXT,NBAST,MOAUX, 3185 * (0,I=1,NSYM), 3186 * (NORB2(I),I=1,NSYM), 3187 * (NBAS(I),I=1,NSYM),POTNUC,EMCSCF 3188 WRITE (LUSIFC) (WRK(NBASC+I),I=1,MAX(NAUXT,4)) 3189 WRITE (LUSIFC) (WRK(3*NBASC+I),I=1,MAX(MOAUX,4)) 3190C Interface for primary and secondary basis sets (WK/UniKA/04-11-2002). 3191 IF (LGLOCAL) THEN 3192 NOCL = NOCC(1) 3193 NBASL = NBAS(1) 3194 NORB1L = NORB1(1) 3195 NORB2L = NORB2(1) 3196 NORBFTL = NORB1L + NORB2L + NOCL 3197 NMOFLTL = NORBFTL * NBASL 3198 IEIGL = IOFFMO 3199 IOCCL = IEIGL + NORBFTL 3200 ICMOL = IOCCL + NBASL * NOCL 3201 CALL DZERO(WRK(IEIGL),NORBFTL) 3202 OPEN(99,FILE='LOCMO',FORM='FORMATTED') 3203 IJ = 0 3204 DO IJ = 1, NBASL * NORB1L 3205 READ(99,'(D30.20)') WRK(ICMOL+IJ-1) 3206 END DO 3207 CLOSE(99) 3208 CALL DCOPY(NBASL*NOCL,WRK(ICMOL),1,WRK(IOCCL),1) 3209 CALL DCOPY(NORB2L*NBASL,WRK(3*NBASC+1),1, 3210 * WRK(IOCCL+NBASL*(NORB1L+NOCL)),1) 3211 WRITE (LUSIFC) LAB123, LABUNI 3212 WRITE (LUSIFC) 1,NORBFTL,NBAST,NMOFLTL,NOCL,NORB1L+ 3213 * NORB2L+NOCL,NBASL,POTNUC,EMCSCF 3214 WRITE (LUSIFC) (WRK(IEIGL+I-1),I=1,MAX(NORBFTL,4)) 3215 WRITE (LUSIFC) (WRK(IOCCL+I-1),I=1,MAX(NMOFLTL,4)) 3216 ELSE 3217 IOFF1 = KFCDIA 3218 IOFF2 = KFCDIA + NBASC 3219 IOFFEV = IOFF2 + NBASC 3220 IOFM1 = 1 3221 IOFM2 = IOFFEV + NBASC 3222 IOFMO = IOFM2 + NBASC*NBASC 3223 NFBAST = 0 3224 DO ISYM = 1, NSYM 3225 CALL DCOPY(NOCC(ISYM),WRK(IOFF1),1,WRK(IOFFEV),1) 3226 IOFFEV = IOFFEV + NOCC(ISYM) 3227 CALL DCOPY(NORB1(ISYM),WRK(IOFF1),1,WRK(IOFFEV),1) 3228 IOFFEV = IOFFEV + NORB1(ISYM) 3229 IOFF1 = IOFF1 + NORB1(ISYM) 3230 CALL DCOPY(NORB2(ISYM),WRK(IOFF2),1,WRK(IOFFEV),1) 3231 IOFF2 = IOFF2 + NORB2(ISYM) 3232 IOFFEV = IOFFEV + NORB2(ISYM) 3233 NBASI = NBAS(ISYM) 3234 CALL DCOPY(NOCC(ISYM)*NBASI,CMO(IOFM1),1,WRK(IOFMO),1) 3235 IOFMO = IOFMO + NOCC(ISYM)*NBASI 3236 CALL DCOPY(NORB1(ISYM)*NBASI,CMO(IOFM1),1,WRK(IOFMO),1) 3237 IOFMO = IOFMO + NORB1(ISYM)*NBASI 3238 IOFM1 = IOFM1 + NORB1(ISYM)*NBASI 3239 CALL DCOPY(NORB2(ISYM)*NBASI,WRK(IOFM2),1,WRK(IOFMO),1) 3240 IOFM2 = IOFM2 + NORB2(ISYM)*NBASI 3241 IOFMO = IOFMO + NORB2(ISYM)*NBASI 3242 NFBAST = NFBAST + 3243 * (NORB1(ISYM) + NORB2(ISYM) + NOCC(ISYM)) * NBASI 3244 END DO 3245 WRITE (LUSIFC) LAB123, LABUNI 3246 WRITE (LUSIFC) NSYM,NORBFT,NBAST,NMOFLT,(NOCC(I),I=1,NSYM), 3247 * (NORB1(I)+NORB2(I)+NOCC(I),I=1,NSYM), 3248 * (NBAS(I),I=1,NSYM),POTNUC,EMCSCF 3249 WRITE (LUSIFC) (WRK(2*NBASC+I),I=1,MAX(NORBFT,4)) 3250 WRITE (LUSIFC) (WRK((3+NBASC)*NBASC+I),I=1,MAX(NMOFLT,4)) 3251 END IF 3252 END IF 3253C 3254C *** end of subroutine WR_SIRIFC 3255C 3256 WRITE (LUSIFC) LAB123,LABEOD 3257 CALL GPCLOSE(LUSIFC,'KEEP') 3258 CALL QEXIT('WR_SIRIFC') 3259 RETURN 3260 END 3261C --- end of siropt.F --- 3262