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 20C FILE : sirius/sirdiis.F 21C Principal author : Hans Jørgen Aa. Jensen 22C First version : 1993 23C 24!=========================================================================== 25!se efter HJTODO 26!HJTODO 940816: if .FREEZE implemented in FCKEIG (requires new ORDRSS) 27! then check all setting of MXDIIS = 0 and MAXFCK = 0 (I think MAXFCK 28! also will work if .FREEZE implemented in FCKEIG). 29! 30!Revision 1.9 2001/02/01 17:30:12 hjj 31!restart DIIS if DIIS is stalled ... 32! 33!Revision 1.6 2000/05/24 12:28:27 hjj 34!1) fixed solvent calculations (DIFDEN must be false for solvent) 35!2) fixed WR_SIRIFC call for open shell and solvent 36! (required change in SOLFCK also) 37!============================================================================= 38!970326-hjaaj -- adjusted abort DIIS algorithm (allow more NEWOCC) 39!951130-hjaaj 40!DIIS_CTL: moved 'reset to not converged' for solvent and writing SIRIFC to 41! DIIS_CTL from SIRCTL; it is a deficiency in the DIIS routines that 42! it cannot be done in them. Look for MAERKE. 43!950824-hjaaj 44!DIIS_CTL: output for AUTOCC 45!FCKEIG: also AUTOCC for open shell plus smaller changes. 46!9505-k.ruud implemented AUTOCC option 47!940701-hjaaj 48!SRDIIS: let MXERRV to default to MAX(2,MIN(NWOPT,10)) (was MXDIIS). 49! This fixes linear dependency problem for very small debug runs: 50! MXERRV=10 gave lin.dep. for 2 el. in 2 orbitals, 51! MXERRV=2 gave quadratic convergence ! 52!940511-hjaaj 53!SRDIIS,DIIS_CTL: call flshfo so progress can be followed in output file; 54! moved heading printing from DIIS_CTL to SRDIIS 55!940505-hjaaj 56!DIIS_CTL: exit if energy increasing or convergence is too slow. 57!940501-hjaaj 58!FCKEIG: exit if ".FREEZE" not compatible with ORDRSS. 59! ORDRSS now always called (this also fixes error in prev. version). 60!940408-hjaaj 61!DIIS_CTL: s/EMCACT/EACTIV/ to define active energy in /INFOPT/ 62!940311-hjaaj 63!DIIS_CTL: s/EVCNRM/GRDNRM/ (so gradient norm is transferred to SIROUT) 64! nice output to LUW4 if LUW4 .ne. LUPRI 65!931126-hjaaj 66!DIIS_CTL: call NEWORB after each iteration (for restart) 67!930720-hjaaj 68!SRDIIS: new parameter PROCC .FALSE. in CALL PRORB 69!930715-hjaaj 70!------: improved some print 71!SRDIIS: exit MXDGSS .GT. 0 72!930623-hjaaj 73!FCKEIG: inserted PARAMETER(D0=0.0D0) 74!=========================================================================== 75C /* Deck srdiis */ 76 SUBROUTINE SRDIIS (ICONV,CMO,WRK,KFRSAV,LFRSAV) 77C 78C Copyright (c) 1993, 1997 Hans Joergen Aa. Jensen 79C 80C Driver for DIIS iterations in SIRIUS. 81C 82 83#include "implicit.h" 84 DIMENSION CMO(*), WRK(*) 85C 86C 87C Used from common blocks: 88C INFINP : DIRFCK, FLAG(*) 89C SCBRHF : MAXEVC 90C INFORB : N2BAST,NNBAST,N2BASX,NCMOT,NSYM, ...,MXDGSS 91C INFVAR : NWOPT 92C INFPRI : IPRRHF 93C 94#include "maxorb.h" 95#include "priunit.h" 96#include "infinp.h" 97#include "scbrhf.h" 98#include "inforb.h" 99#include "infvar.h" 100#include "infpri.h" 101#include "gnrinf.h" 102C dftcom.h : DFT_SPINDNS 103#include "dftcom.h" 104C 105Cholesky 106#include "ccdeco.h" 107#include "choscf.h" 108C 109 PARAMETER (NBLOC = 1) 110 INTEGER LUBLOC(NBLOC) 111 DATA LUBLOC / 21 / 112Cholesky 113C 114 PARAMETER (D4 = 4.0D0) 115 LOGICAL C2DIIS 116C 117 CALL QENTER('SRDIIS') 118 ICONV = 0 119 IF (MXDGSS .GT. 1) THEN 120 NINFO = NINFO + 1 121 WRITE (LUPRI,'(//A/)') 122 & ' INFO: no DIIS iterations because degenerate'// 123 & 'supersymmetries is not implemented for DIIS' 124 GO TO 9999 125 END IF 126C 127 KFREE = KFRSAV 128 LFREE = LFRSAV 129 130C 131CHJMAERKE: define IPR_DIIS separately? 132C 133 MXERRV = MAXEVC 134 IF (MXERRV .LE. 0) THEN 135 MXERRV = MIN(NWOPT-1,8) 136C ... fixes linear dependency problem for very small debug runs 137C MXERRV=10 gave lin.dep. for 2 el. in 2 orbitals, 138C MXERRV=2 gave quadratic convergence ! (940701-hjaaj) 139C (180705-hjaaj: in a small test case with NWOPT=3, the full space 140C MXERRV=NWOPT gave lin.dep. problem for some solvers, therefore now "NWOPT-1") 141C (130919-hjaaj: decreased default from 10 to 5, because 142C with 10 the span of BMAT may be >1.D15, giving problems 143C in the DIIS solver - an example has been seen !!!) 144C (140307-jmho: increased default from 5 to 8, because 145C with 5 it often resulted in 1-2 extra SCF iterations. At 8 146C we should still avoid the problem mentioned above.) 147 MXERRV = MAX(2,MXERRV) 148 END IF 149 IPR_DIIS = IPRRHF 150 C2DIIS = FLAG(22) 151C 152 LUWOUT = LUPRI 153 10 WRITE (LUWOUT,'(//A/A/A)') 154 & ' ***********************************************', 155 & ' ***** DIIS acceleration of SCF iterations *****', 156 & ' ***********************************************' 157Cholesky 158 IF (CHOINT) THEN 159 WRITE(LUWOUT,'(A)') 160 & ' -- using Cholesky decomposed integrals' 161 IF (.NOT. CCMODSK) THEN 162 WRITE(LUWOUT,'(A,1P,D9.2,A)') 163 & ' -- using Cholesky decomposed density, thr = ',THRDC 164 END IF 165 END IF 166Cholesky 167 IF (C2DIIS) THEN 168 WRITE (LUWOUT,'(/A,I5)') 169 & ' C2-DIIS algorithm; max error vectors =',MXERRV 170 ELSE 171 WRITE (LUWOUT,'(/A,I5)') 172 & ' C1-DIIS algorithm; max error vectors =',MXERRV 173 END IF 174 CALL FLSHFO(LUWOUT) 175 IF (LUWOUT .NE. LUW4) THEN 176 LUWOUT = LUW4 177 GO TO 10 178 END IF 179C 180C Get AO overlap matrix; transform left index to CMO1 basis. 181C 182 CALL MEMGET2('REAL','SMOAO',KSMOAO,N2BAST,WRK,KFREE,LFREE) 183 KREL = KFREE 184 CALL MEMGET2('REAL','SAO' ,KSAO ,NNBAST,WRK,KFREE,LFREE) 185 CALL MEMGET2('REAL','STMP',KSTMP,N2BAST,WRK,KFREE,LFREE) 186 CALL RDONEL('OVERLAP',.TRUE.,WRK(KSAO),NNBAST) 187 IF (IPR_DIIS .GT. 25) THEN 188 WRITE (LUPRI,'(/A)') ' SRDIIS test output of SAO:' 189 CALL OUTPKB(WRK(KSAO),NBAS,NSYM,-1,LUPRI) 190 WRITE (LUPRI,'(/A)') ' SRDIIS test output of CMO1:' 191 CALL PRORB(CMO,.FALSE.,LUPRI) 192C CALL PRORB(CMO,PROCC,IOUT) 193 END IF 194 DO 100 I = 1,NSYM 195 NORBI = NORB(I) 196 IF (NORBI .GT. 0) THEN 197 NBASI = NBAS(I) 198 CALL DSPTSI(NBASI,WRK(KSAO+IIBAS(I)),WRK(KSTMP)) 199 CALL DGEMM('T','N',NORBI,NBASI,NBASI,1.D0, 200 & CMO(1+ICMO(I)),NBASI, 201 & WRK(KSTMP),NBASI,0.D0, 202 & WRK(KSMOAO+I2BAS(I)),NORBI) 203 END IF 204 100 CONTINUE 205C scale by 4 to get correct orbital gradient from DIIS_EVC 206 CALL DSCAL(N2BAST,D4,WRK(KSMOAO),1) 207 CALL MEMREL('SRDIIS.SMOAO',WRK,KFRSAV,KREL,KFREE,LFREE) 208C 209 IF (DFT_SPINDNS) THEN 210 NDMAT = 3 211 ELSE IF (NASHT .EQ. 0) THEN 212 NDMAT = 1 213 ELSE 214 NDMAT = 2 215 END IF 216 CALL MEMGET2('REAL','CMO1',KCMO1 ,NCMOT, WRK,KFREE,LFREE) 217 CALL MEMGET2('REAL','DCAO',KDCAO ,NDMAT*N2BASX,WRK,KFREE,LFREE) 218 CALL MEMGET2('REAL','FCAO',KFCAO ,NDMAT*N2BASX,WRK,KFREE,LFREE) 219 CALL MEMGET2('REAL','H1AO',KH1AO ,NNBAST,WRK,KFREE,LFREE) 220 CALL MEMGET2('REAL','DSAV',KDSAV ,NDMAT*NNBAST,WRK,KFREE,LFREE) 221 IF (FLAG(16)) THEN 222 CALL MEMGET2('REAL','TNLM',KTNLM,2*NLMSOL,WRK,KFREE,LFREE) 223 ELSE 224 CALL MEMGET2('REAL','TNLM',KTNLM,0,WRK,KFREE,LFREE) 225 END IF 226C estimate max value for MXERRV (951130-hjaaj) : 227 IF (DIRFCK) THEN 228C -- direct SCF: 229C NDMAT*N2BASX used in HERFCK(hermit) 230C N2BASX used in SKLFC1(hermit) 231C 1000000 (1 mw estimate for TWOINT(hermit) and the rest) 232 NEVC = (LFREE - (NDMAT+1)*N2BASX - 1000000) / NNBAST 233 ELSE 234C -- conventional SCF: 235C N2BASX used in FCKDEN 236C somewhat arbitrary estimate of memory needed in RDSUPM etc. 237C estimate quite big to be on the safe side /960627-hjaaj 238 NEVC = (LFREE - N2BASX - 100000) / NNBAST 239 END IF 240 NEVC = NEVC / (NDMAT + 1) 241 IF (NEVC .LT. MXERRV) THEN 242 IF (NEVC .LT. 3) THEN 243 NWARN = NWARN + 1 244 WRITE (LUPRI,'(/A/A,I5)') 245 & ' SRDIIS WARNING: insufficient memory for DIIS', 246 & ' max # of error vectors only',NEVC 247 GO TO 9990 248 END IF 249 NINFO = NINFO + 1 250 WRITE (LUPRI,'(/A,I5,A,I5/A)') 251 & ' SRDIIS INFO: MXERRV reduced from',MXERRV,' to',NEVC, 252 & ' SRDIIS INFO: because of memory limitations.' 253 MXERRV = NEVC 254 END IF 255C 256 LBMAT = MXERRV*MXERRV 257 ! was LBMAT = (MXERRV+1)*MXERRV / 2 258 ! but in opt-solvers the square matrix is used 259 ! and memory requirement is small so no case testing here ... 260 CALL MEMGET2('REAL','BMAT',KBMAT ,LBMAT, WRK,KFREE,LFREE) 261 CALL MEMGET2('REAL','CVEC',KCVEC ,MXERRV+1,WRK,KFREE,LFREE) 262 CALL MEMGET2('REAL','GVEC',KGVEC ,MXERRV+1,WRK,KFREE,LFREE) 263 CALL MEMGET2('REAL','FSAV',KFSAV ,NDMAT*MXERRV*NNBAST, 264 & WRK,KFREE,LFREE) 265 CALL MEMGET2('REAL','EVSAV',KEVSAV,MXERRV*NNORBT,WRK,KFREE,LFREE) 266 CALL MEMGET2('REAL','EGSAV',KEGSAV,MXERRV*2 ,WRK,KFREE,LFREE) 267 CALL MEMGET2('REAL','DV',KDV,NNASHX,WRK,KFREE,LFREE) 268 CALL DCOPY(NCMOT,CMO,1,WRK(KCMO1),1) 269 CALL DIIS_CTL(CMO,WRK(KCMO1),WRK(KGVEC),WRK(KBMAT),WRK(KH1AO), 270 & WRK(KFSAV),WRK(KEVSAV),WRK(KEGSAV),MXERRV,C2DIIS, 271 & WRK(KSMOAO),WRK(KDCAO),WRK(KFCAO),WRK(KDSAV),WRK(KDV), 272 & WRK(KCVEC),WRK(KTNLM), 273 & WRK,KFREE,LFREE, IPR_DIIS,ICONV) 274C 275 9990 CALL MEMREL('SRDIIS',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 276 9999 CALL FLSHFO(LUW4) 277 IF (LUPRI .NE. LUW4) CALL FLSHFO(LUPRI) 278 CALL QEXIT('SRDIIS') 279 RETURN 280 END 281C /* Deck DIIS_CTL */ 282 SUBROUTINE DIIS_CTL(CMO,CMO1,GVEC,BMAT,H1AO,FAOSAV,EVCSAV,EGSAV, 283 * MXERRV,C2DIIS,SMOAO,DCAO,FCAO,DAOSAV,DV, 284 * CVEC,TNLM,WRK,KFRSAV,LFRSAV, 285 * IPR_DIIS,ICONV) 286 287C 288C L.r. 5-May-1994 hjaaj / 1-Apr-1997 289C 290C DFT modifications T. Helgaker 291C 292C 293 use qmcmm, only: getdim_relmat, comp_relmat, write_relmat, 294 & getdim_mmmat, comp_mmrelmat 295 use qmcmm_fock, only: qmnpmm_fock 296#ifdef HAS_PCMSOLVER 297 use pcm_scf 298 use pcm_config, only: pcm_configuration, pcm_cfg 299#endif 300 use pelib_interface, only: use_pelib, pelib_ifc_fock, 301 & pelib_ifc_energy, pelib_ifc_mep, 302 & pelib_ifc_do_mep, pelib_ifc_do_savden, 303 & pelib_ifc_save_density 304#include "implicit.h" 305#include "dummy.h" 306 DIMENSION CMO(*), CMO1(*), GVEC(*), BMAT(*), H1AO(*) 307 DIMENSION FAOSAV(NNBAST,MXERRV,2), EVCSAV(NNORBT,MXERRV) 308 DIMENSION SMOAO(*), DCAO(*), FCAO(*), DAOSAV(NNBAST,*), CVEC(*) 309 DIMENSION DV(*) 310 DIMENSION EGSAV(2,MXERRV), TNLM(*), WRK(*) 311 312#ifdef HAS_PCMSOLVER 313 real(8), allocatable :: density_ao(:), density_ao_symm(:) 314 real(8), allocatable :: fock_pcm_ao(:), fock_pcm_ao_symm(:) 315#endif 316 317C 318C Used from common blocks: 319C INFINP : FLAG(*), HSROHF, DODFT; ? 320C INFORB : NCMOT,NNBAST,NSYM 321C INFOPT : EPOT,EMCSCF,EMY,EACTIV,GRDNRM 322C SCBRHF : THRRHF,MXDIIS, ? 323C DFTCOM : HFXFAC, DODFTD 324C INFVAR : NWOPH 325C SCBRHF : IOPRHF,AUTOCC 326C CBIREA : LMULBS 327C R12INT : R12CAL 328C 329#include "maxorb.h" 330#include "priunit.h" 331#include "mxcent.h" 332 333#include "molde.h" 334#include "infinp.h" 335#include "pcmlog.h" 336#include "inforb.h" 337#include "infvar.h" 338#include "infopt.h" 339#include "scbrhf.h" 340#include "dftcom.h" 341#include "cbihr2.h" 342#include "infpri.h" 343#include "gnrinf.h" 344#include "dfterg.h" 345#include "cbirea.h" 346#include "r12int.h" 347#include "qm3.h" 348#include "incore.h" 349#include "qmmm.h" 350#include "infloc.h" 351#include "mmtimes.h" 352#include "infpar.h" 353C 354Cholesky 355#include "ccdeco.h" 356#include "sirftim.h" 357#include "choscf.h" 358C 359#include "center.h" 360#include "nuclei.h" 361C 362Cholesky 363C 364C------------------------ 365#include "pcmdef.h" 366#include "pcm.h" 367 368! QM/NP/MM code 369#include "qmnpmm.h" 370 371 PARAMETER (THDDEF = 1.0D-5) 372 PARAMETER (D0 = 0.0D0, DP5 = 0.5D0, D1 = 1.0D0, DM1 = -1.0D0) 373 PARAMETER (THREVC = 0.01D0, THRINC = 1.0D-10, CNVFAC = 0.8D0) 374C 375 LOGICAL ONLYFC, C2DIIS, DIFDEN, DODIFDEN, NEWTHR, RHFWOP, OLDWOP 376 LOGICAL LCONV, LNEWOCC, HSROHF_save 377 DIMENSION NISH_save(8), NCAN(8) 378 CHARACTER*8 CMOLBL 379 380 REAL(8) :: PE_ENERGY 381 REAL(8), DIMENSION(:), ALLOCATABLE :: PE_DENMAT, PE_FCKMAT 382 383 real(8), allocatable :: relay_matrix(:, :, :) 384 real(8), allocatable :: mm_relay_matrix(:) 385 real(8), allocatable :: qmcmm_sol(:) 386 387!#define EXPERIMENTAL_DMAT_RESTART 388#ifdef EXPERIMENTAL_DMAT_RESTART 389 ! TODO: do not allocate extra matrix, not needed 390 ! but read directly into DCAO 391 ! then DSCAL by 2.0 to match Dalton 392 real(8), allocatable :: dmat_restart(:) 393 logical :: is_gcbasis 394#endif 395C 396 CALL QENTER('DIIS_CTL') 397 KFREE = KFRSAV 398 LFREE = LFRSAV 399 IFTHSV = IFTHRS 400C 401C save current IOPRHF for check below 402 IOPRHF_old = IOPRHF 403 IF (NSYM .EQ. 1) AUTOCC = .FALSE. 404 405! Sep 2013 hjaaj: now HSROHF code here also for NASHT.eq.1 because we 406! have some evidence that the open shell Fock matrix from HSROHF 407! gives better convergence than the one from the older "single open shell". 408! At return we restore HSROHF to .false. for this case, because some 409! ABACUS and other modules require the old "single open shell" 410! to work correctly. 411 412 HSROHF_save = HSROHF 413 HSROHF = HSROHF .OR. (NASHT .EQ. 1 .AND. .NOT. AUTOCC) 414 415C 416C Get one-electron Hamiltonian 417C 418 CALL SIRH1(H1AO,WRK(KFREE),LFREE) 419 420CAMT If requested calculate the DFT-D dispersion energy correction 421 IF (DODFT .AND. DODFTD) THEN 422 NDERIV=0 423 CALL DFT_D_DAL_IFC(EDISP,NDERIV,WRK(KFREE),LFREE) 424 ELSE 425 EDISP = 0.0D0 426 END IF 427 428C 429 IF (LUW4 .NE. LUPRI .OR. IPR_DIIS .LE. 1) THEN 430#ifdef HAS_PCMSOLVER 431 IF (FLAG(16) .OR. PCM .OR. QM3 .OR. QMMM .OR. QMNPMM .OR. 432 & pcm_cfg%do_pcm) THEN 433#else 434 IF (FLAG(16) .OR. PCM .OR. QM3 .OR. QMMM .OR. QMNPMM) THEN 435#endif 436C this also includes the case of MMPCM 437 IF (.NOT.AUTOCC) THEN 438 WRITE (LUW4,'(/A)') 439 & ' Iter Total energy Solvation energy '// 440 & 'Error norm Delta(E)' 441 ELSE IF (IOPRHF .EQ. 0) THEN 442 WRITE (LUW4,'(/A,I4,A//A)') 443 & ' Automatic occupation of symmetries with',NRHFEL, 444 & ' electrons.', 445 & ' Iter Total energy Solvation energy Error norm'// 446 & ' Delta(E) SCF occupation' 447 ELSE 448 WRITE (LUW4,'(/A,I4,A/A//A)') 449 & ' Automatic occupation of symmetries with',NRHFEL, 450 & ' electrons.', 451 & ' "sym." is the symmetry of the one open shell orbital '// 452 & 'in table below.', 453 & ' Iter Total energy Solvation energy Error norm'// 454 & ' Delta(E) Sym. Closed shell occupation' 455 END IF 456 ELSE IF (USE_PELIB()) THEN 457 WRITE(LUW4,'(/A)') 458 & ' Iter Total energy Embedding energy '// 459 & 'Error norm Delta(E)' 460 ELSE ! no solvation energy 461 IF (.NOT.AUTOCC) THEN 462 WRITE (LUW4,'(/A)') 463 & ' Iter Total energy Error norm Delta(E)'// 464 & ' DIIS dim.' 465 ELSE IF (IOPRHF .EQ. 0) THEN 466 WRITE (LUW4,'(/A,I4,A//A)') 467 & ' Automatic occupation of symmetries with',NRHFEL, 468 & ' electrons.', 469 & ' Iter Total energy Error norm Delta(E) '// 470 & ' SCF occupation' 471 ELSE 472 WRITE (LUW4,'(/A,I4,A/A//A)') 473 & ' Automatic occupation of symmetries with',NRHFEL, 474 & ' electrons.', 475 & ' "sym." is the symmetry of the one open shell orbital '// 476 & 'in table below.', 477 & ' Iter Total energy Error norm Delta(E) '// 478 & 'Sym. Closed shell occupation' 479 END IF 480 END IF 481 END IF 482C 483 IF (NISHT.EQ.0) CALL DZERO(DAOSAV,NNBAST) 484 ONLYFC = (NASHT .EQ. 0) 485 IF (ONLYFC) THEN 486C DV(1) = D0 487 NDMAT = 1 488 ELSE IF (NASHT .GT. 0) THEN 489 CALL DZERO(DV,NNASHX) 490 DO I = 1, NASHT 491 II = I*(I+1)/2 492 DV(II) = D1 493 END DO 494 NDMAT = 2 495 ELSE 496 CALL QUIT('DIIS_CTL called with NASHT .gt. 1') 497 END IF 498 THDIIS = THRRHF 499 IF (THDIIS .LE. D0) THDIIS = THDDEF 500C 501C Level shift SHFTLVL ("restricted step" type step reduction) 502C - initial value: see next code block 503C - used to modify FD matrix: FDAO(shifted) = FDAO + SHFTLVL*DCAO 504C - this corresponds to increasing occ-unocc orb.en. gap by SHFTLVL 505C and approximately to add abs(SHFTLVL) to the approximate 506C Hessian diagonal implicit in the Roothaan diagonalization. 507C -- Feb. 2001, Hans Joergen Aa. Jensen. 508C 509c DEFLVL is read from SCF INPUT/.SHIFT 510 IF (DEFLVL .gt. 0.9D0) THEN ! .SHIFT not set in input 511 SHFTLVL = -0.20D0 512 DEFLVL = 0.0D0 ! no default level shift in e.g. next geometry 513 ELSE ! .SHIFT has ben set in input 514 SHFTLVL = DEFLVL 515 END IF 516 SHFTFAC = 0.5D0 517 518C --- hjaaj: find safe screening factor depending on 519C convergence threshold: 520C SCRFAC = 100*NBAST 521 SCRFAC = 100*N2BAST 522 IFTMAX = INT(-LOG10(THDIIS/SCRFAC)) + 1 523 IFTMIN = INT(-LOG10( 0.1D0/SCRFAC)) + 1 524Chjaaj: screening factor for safe Fock matrix construction 525C is chosen to be 100 (we want accuracy to 0.1 times THDIIS) 526C F_ij ~ sum(kl) (ij|kl) * D_kl, where sum(kl) is over 527C max N2BAST elements (for totally symmetric "kl"); 528C statistically we expect error to be ca. sqrt(n2bast) < nbast 529C The 100*NBAST is thus with an extra factor 10 factor for 530C safety, and because IFTHRS only can change screening 531C with a factor 10. /01-Feb-2001 hjaaj 532C 533C Yes, but for the gradient norm we in addition will have a sum 534C of ca. NOCCT*NVIRT/NSYM elements, thus to be safe wrt 535C gradient norm let us use 10*N2BAST in SCRFAC instead 536C of 100*NBAST. 537C Also minimum accuracy corresponding to error in 538C GRDNORM ~ 0.1 (IFTMIN). /31-Jul-2006 hjaaj 539C 540C In iterations we want 0.001 times GRDNRM because we do 541C keep vectors until GRDNRM lowered with a factor 100 542 DAMP = D0 543 EMCLOW = D0 544 GRDNRM = -1.0D0 545 ITNOCC = 0 546 NEWOCC = 0 547 LNEWOCC = AUTOCC 548C if (LNEWOCC) then write new BASINFO to SIRIUS.RST in NEWORB 549 I_STALL = 0 550C------------------ 551 DODIFDEN = .NOT.FLAG(49) .AND. .NOT.HSROHF 552 & .AND. .NOT.DODFT .AND. .NOT.DOHFSRDFT 553 & .AND. .NOT.AOSAVE .AND. .NOT.CHOINT 554 & .AND. .NOT.EMBEDDING 555C 556#ifdef HAS_PCMSOLVER 557 if (pcm_cfg%do_pcm) then 558 dodifden = .false. 559 end if 560#endif 561C------------------ 562Chjaaj: DIFDEN does not work with solvent/PCM/QM3 563C because FC_sol is added to the old FC matrix, but not the new! 564Ctuh: DIFDEN does not work with DFT either 565Cov: DIFDEN disabled for high spin HF 566C-tbp: turn off DIFDEN for Cholesky decomposition (CHOINT true) 567 DIFDEN = DODIFDEN 568C------------------ 569 IF (.NOT. USRSCR) THEN 570C AOSAVE or Cholesky (CHOINT): No screening. 571 IF (AOSAVE .OR. CHOINT .OR. .NOT. DIRFCK) THEN 572 IFTHRS = 20 ! code for no screening 573 ELSE IF (FLAG(13)) THEN 574Chj-aug99: MO's read from input, assume good guess 575 IFTHRS = IFTMAX 576 ELSE 577Chj-aug99: MO's from huckel or H1DIAG, assume poor guess 578 IFTHRS = IFTMIN 579 END IF 580 END IF 581 IF (ICI0 .EQ. 6) THEN 582Chj ... if (GEOWLK) then 583 EMCGEO = EMCOLD 584 DEPGEO = DEPRED 585 END IF 586 587C 588 ITDIIS = 0 589 JTDIIS = 0 590C 591C Projection operators for removing non-rotated orbitals (.FREEZE option) 592C A -> P_f*A*P_f(T), P_f = 1-S*D_f 593C 594 IF (LNOROT) THEN 595 CALL MEMGET('REAL',KPFROZ,N2BASX,WRK,KFREE,LFREE) 596 CALL MEMGET('REAL',KDFROZ,N2BASX,WRK,KFREE,LFREE) 597 CALL MEMGET('REAL',KS,N2BASX,WRK,KFREE,LFREE) 598 CALL GET_H1(WRK(KS),'OVERLAP',WRK(KFREE),LFREE) 599 CALL DUNIT(WRK(KPFROZ),NBAST) 600 CALL DZERO(WRK(KDFROZ),N2BASX) 601 KDFI=KDFROZ 602 DO ISYM=1,NSYM 603 IORBI=IORB(ISYM) 604 NORBI=NORB(ISYM) 605 DO I=1,NORBI 606 IF (NOROT(IORBI+I).NE.0) THEN 607 NBASI=NBAS(ISYM) 608 ICMOI=ICMO(ISYM)+NBASI*(I-1)+1 609 ! C*C.T() 610 CALL DGEMM( 611 & 'N','T',NBASI,NBASI,1, 612 & D1,CMO(ICMOI),NBASI, 613 & CMO(ICMOI),NBASI, 614 & D1,WRK(KDFI),NBAST 615 & ) 616 END IF 617 END DO 618 KDFI=KDFI+NBAST*NBAS(ISYM)+NBAS(ISYM) 619 END DO 620 CALL DGEMM('N','N',NBAST,NBAST,NBAST, 621 & DM1,WRK(KS),NBAST, 622 & WRK(KDFROZ),NBAST, 623 & D1, WRK(KPFROZ),NBAST 624 & ) 625 IF (IPR_DIIS.GT.20) THEN 626 CALL HEADER('Start orbitals',-1) 627 CALL PRORB(CMO,.FALSE.,LUPRI) 628 CALL HEADER('Frozen orbital density',-1) 629 CALL OUTPUT( 630 & WRK(KDFROZ),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI 631 & ) 632 CALL HEADER('Frozen orbital projector',-1) 633 CALL OUTPUT( 634 & WRK(KPFROZ),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI 635 & ) 636 END IF 637 END IF 638 639 ITDIIS = 0 640 JTDIIS = 0 641 642 100 CONTINUE 643 WRITE(LUW4,'(A)') 644 & ' -----------------------------------------------------'// 645 & '------------------------' 646 ITDIIS = ITDIIS + 1 647 JTDIIS = JTDIIS + 1 648 ITNOCC = ITNOCC + 1 649 INDEVC = MOD(JTDIIS-1,MXERRV) + 1 650 NEVC = MIN(MXERRV,JTDIIS) 651C 652C 653 IF (IPR_DIIS .GT. 1) WRITE (LUPRI,'(/A,I5)') 654 & ' *** DIIS iteration number',ITDIIS 655 CALL FLSHFO(LUPRI) 656C 657C ***** Construct folded inactive and active parts of the one-electron 658C density matrix in AO-basis (DCAO and DVAO). DV is the active 659C part of the one-electron density matrix in MO-basis. 660C 661 CALL FCKDEN((NISHT.GT.0.OR.HSROHF),(.NOT.ONLYFC), 662 & DCAO,DCAO(1+N2BASX),CMO,DV,WRK(KFREE),LFREE) 663C CALL FCKDEN(GETDC,GETDV,DCAO,DVAO,CMO,DV,WRK,LWRK) 664 665#ifdef EXPERIMENTAL_DMAT_RESTART 666! radovan: experimental and untested feature, use at own risk 667! reads AO dmat from LSDalton's dens.restart 668 669 if (itdiis == 1) then 670 iunit = -1 671 CALL GPOPEN(iunit,'dens.restart','UNKNOWN',' ', 672 & 'UNFORMATTED',IDUMMY, 673 & .FALSE.) 674 nrow = 0 675 ncol = 0 676 read(iunit) nrow, ncol 677 678 ! check that dimensions match 679 if (nrow*ncol /= n2basx) then 680 call quit('dimensions do not match in restart') 681 end if 682 683 ! check that there is no symmetry 684 if (nsym > 1) then 685 call quit('nsym > 1 in restart') 686 end if 687 688 allocate(dmat_restart(nrow*ncol)) 689 read(iunit) (dmat_restart(I),I=1,nrow*ncol) 690 read(iunit) is_gcbasis 691 692 ! check that gcbasis is .false. 693 if (is_gcbasis == .true.) then 694 call quit('is_gcbasis is true in restart') 695 end if 696 697 CALL GPCLOSE(iunit, 'KEEP') 698 dmat_restart = 2.0*dmat_restart 699 CALL DCOPY(nrow*ncol,dmat_restart,1,DCAO,1) 700 deallocate(dmat_restart) 701 end if 702#endif /* EXPERIMENTAL_DMAT_RESTART */ 703 704C 705C For a restricted open shell DFT calculation, or high spin 706C Let the matrices be Fc=Fi+Fa and Fc-Fo=Fa-Q=-Fa(exch) 707C generated by total and spin (=active) densities 708C This choice avoids the calculation of a third Q matrix 709C 710 IF (HSROHF) THEN 711 CALL DAXPY(N2BASX,D1,DCAO(1+N2BASX),1,DCAO,1) 712 CALL DSCAL(N2BASX,-D1,DCAO(1+N2BASX),1) 713 END IF 714C 715C-tbp: For Cholesky, save a copy of the occupied part of CMO 716C on disk (for use in SIRFCK) 717C 718 IF (CHOINT .AND. CCMODSK) THEN 719 CALL WOPEN2(LUCCMO,FCCMO,64,0) 720 IADR2 = 1 721 DO ISYM = 1,NSYM 722 NDVCS(ISYM) = NISH(ISYM) 723 IF (NISH(ISYM) .GT. 0) THEN 724 JCMO = ICMO(ISYM) + 1 725 LEN = NBAS(ISYM)*NISH(ISYM) 726 CALL PUTWA2(LUCCMO,FCCMO,CMO(JCMO),IADR2,LEN) 727 IADR2 = IADR2 + LEN 728 END IF 729 ENDDO 730 CALL WCLOSE2(LUCCMO,FCCMO,'KEEP') 731 END IF 732C 733 CALL MEMGET2('REAL','TMP1',KTMP1,NNBASX,WRK,KFREE,LFREE) 734 CALL MEMGET2('REAL','TMP2',KTMP2,N2BASX,WRK,KFREE,LFREE) 735 IF (NISHT.GT.0 .OR. HSROHF) THEN 736 CALL DGEFSP(NBAST,DCAO,WRK(KTMP1)) 737 IF (.NOT.DIFDEN .OR. NEVC .EQ. 1) THEN 738 CALL PKSYM1(WRK(KTMP1),DAOSAV,NBAS,NSYM,1) 739C ... save DCAO in DAOSAV 740 DCOVLP = D0 741 ELSE 742 CALL PKSYM1(WRK(KTMP1),WRK(KTMP2),NBAS,NSYM,1) 743 DCOVLP = DDOT(NNBAST,DAOSAV,1,DAOSAV,1) 744 DCOVLP = DDOT(NNBAST,WRK(KTMP2),1,DAOSAV,1) / DCOVLP 745 CALL PKSYM1(WRK(KTMP1),DAOSAV,NBAS,NSYM,-1) 746 CALL DCOPY(NNBAST,WRK(KTMP2),1,DAOSAV,1) 747C ... save new DCAO in DAOSAV 748 CALL DUNFLD(NBAST,WRK(KTMP1),WRK(KTMP2)) 749 CALL DAXPY(N2BASX,(-DCOVLP),WRK(KTMP2),1,DCAO,1) 750C ... difden: DCAO = DCAO - DCOVLP*DAOSAV 751 END IF 752 END IF 753 IF (.NOT.ONLYFC) THEN 754 CALL DGEFSP(NBAST,DCAO(1+N2BASX),WRK(KTMP1)) 755 IF (.NOT.DIFDEN .OR. NEVC .EQ. 1) THEN 756 CALL PKSYM1(WRK(KTMP1),DAOSAV(1,2),NBAS,NSYM,1) 757C ... save "DVAO" in DAOSAV(,2) 758 DVOVLP = D0 759 ELSE 760 CALL PKSYM1(WRK(KTMP1),WRK(KTMP2),NBAS,NSYM,1) 761 DVOVLP = DDOT(NNBAST,DAOSAV(1,2),1,DAOSAV(1,2),1) 762 DVOVLP = DDOT(NNBAST,WRK(KTMP2),1,DAOSAV(1,2),1) 763 & / DVOVLP 764 CALL PKSYM1(WRK(KTMP1),DAOSAV(1,2),NBAS,NSYM,-1) 765 CALL DCOPY(NNBAST,WRK(KTMP2),1,DAOSAV(1,2),1) 766C ... save new "DVAO" in DAOSAV(,2) 767 CALL DUNFLD(NBAST,WRK(KTMP1),WRK(KTMP2)) 768 CALL DAXPY(N2BASX,(-DVOVLP),WRK(KTMP2),1, 769 & DCAO(1+N2BASX),1) 770C ... difden: "DVAO" = "DVAO" - DVOVLP*DVOSAV(,2) 771 END IF 772 END IF 773 CALL MEMREL('DIIS_CTL.DIFDEN',WRK,1,KTMP1,KFREE,LFREE) 774C 775C ***** Calculate 2-electron part of Fock matrix with differential density 776C 777 778 CALL GETTIM(t1_fck,w1_fck) 779 780 CALL FCK2AO(ONLYFC,FCAO,DCAO,WRK,KFREE,LFREE) 781 782 IF (DIRFCK .AND. IFTHRS.LT.20) THEN 783 IF (IPR_DIIS .GE. 0) THEN 784 CALL GETTIM(t2_fck,w2_fck) 785 WRITE (LUPRI,'(I4,A,2I5,L5,1P,2D10.2)') 786 & ITDIIS,' Screening settings (-IFTHRS, JTDIIS, DIFDEN, times)', 787 & -IFTHRS, JTDIIS,DIFDEN,T2_FCK-T1_FCK,W2_FCK-W1_FCK 788 END IF 789 END IF 790 791C 792C FAOSAV(1,INDEVC,1) = FD2AO = FC2AO + FV2AO 793C FAOSAV(1,INDEVC,2) = FV2AO 794C 795 IF (NSYM .GT. 1) THEN 796 CALL MEMGET2('REAL','TMP1',KTMP1,NNBASX,WRK,KFREE,LFREE) 797 CALL DGETSP(NBAST,FCAO,WRK(KTMP1)) 798 CALL PKSYM1(WRK(KTMP1),FAOSAV(1,INDEVC,1),NBAS,NSYM,1) 799 IF (.NOT.ONLYFC) THEN 800 CALL DGETSP(NBAST,FCAO(1+N2BASX),WRK(KTMP1)) 801 CALL PKSYM1(WRK(KTMP1),FAOSAV(1,INDEVC,2),NBAS,NSYM,1) 802 END IF 803 CALL MEMREL('DIIS_CTL.FAOSAV',WRK,1,KTMP1,KFREE,LFREE) 804 ELSE 805 CALL DGETSP(NBAST,FCAO,FAOSAV(1,INDEVC,1)) 806 IF (.NOT.ONLYFC) THEN 807 CALL DGETSP(NBAST,FCAO(1+N2BASX),FAOSAV(1,INDEVC,2)) 808 END IF 809 END IF 810 811 IF (DIFDEN .AND. NEVC .GT. 1) THEN 812 INDEVP = MOD(JTDIIS-2,MXERRV) + 1 813 CALL DAXPY(NNBAST,DCOVLP,FAOSAV(1,INDEVP,1),1, 814 & FAOSAV(1,INDEVC,1),1) 815C .. add DCOVLP*(FCAO_previous + FVAO_previous) 816 IF (.NOT.ONLYFC) THEN 817 IF (.NOT.HSROHF) 818 & CALL DAXPY(NNBAST,-DCOVLP,FAOSAV(1,INDEVP,2),1, 819 & FAOSAV(1,INDEVC,1),1) 820C .. subtract DCOVLP*FVAO_previous to get DCOVLP*FCAO_previous 821 CALL DAXPY(NNBAST, DVOVLP,FAOSAV(1,INDEVP,2),1, 822 & FAOSAV(1,INDEVC,2),1) 823 END IF 824 END IF 825 EMY = DDOT(NNBAST,DAOSAV,1,H1AO,1) 826 & + DDOT(NNBAST,DAOSAV,1,FAOSAV(1,INDEVC,1),1) * DP5 827 IF (DOHFSRDFT) EMY = EMY + EDFTY + ESRDFTY 828C 829 IF (ONLYFC) THEN 830 EACTIV = D0 831 ELSE 832 IF (HSROHF) THEN 833C FAOSAV(:,INDEVC,1) = FDAO = FCAO + FVAO 834C FAOSAV(:,INDEVC,2) = -FVAO_exchange_only 835 EACTIV = DP5 * 836 & DDOT(NNBAST,DAOSAV(1,2),1,FAOSAV(1,INDEVC,2),1) 837 ELSE ! NASHT.eq.1, i.e. PV(1) = 0, i.e. no H2AC contribution. 838C FAOSAV(:,INDEVC,1) = FCAO 839C FAOSAV(:,INDEVC,2) = FVAO 840C EACTIV = DV(u,u)FC(u,u) = FC(u,u) 841C = DVAO(p,q)FC2AO(p,q) + DVAO(p,q)H1AO(p,q) 842 EACTIV = DDOT(NNBAS(IOPRHF),DAOSAV(1+IIBAS(IOPRHF),2),1, 843 & FAOSAV(1+IIBAS(IOPRHF),INDEVC,1),1) 844 & + DDOT(NNBAS(IOPRHF),DAOSAV(1+IIBAS(IOPRHF),2),1, 845 & H1AO(1+IIBAS(IOPRHF)),1) 846 847C make FAOSAV(:,INDEVC,1) = FDAO = FCAO + FVAO 848 CALL DAXPY(NNBAST,D1,FAOSAV(1,INDEVC,2),1, 849 & FAOSAV(1,INDEVC,1),1) 850 END IF 851 END IF 852 853 IF (DODFT) then 854 IF (HSROHF) THEN 855C 856C Reset densities from Dtot, -DV to DC, DV for dft 857C 858 CALL DSCAL(N2BASX,-D1,DCAO(1+N2BASX),1) 859 CALL DAXPY(N2BASX,-D1,DCAO(1+N2BASX),1,DCAO,1) 860 END IF 861 862 CALL DFT_ADD_KS(ONLYFC,EMY,DCAO,NDMAT, 863 & FAOSAV(1,INDEVC,1),FAOSAV(1,INDEVC,2), 864 & WRK(KFREE),LFREE,IPR_DIIS) 865 END IF 866 867CAMT If calculated add the empirical dispersion energy correction 868 EMY = EMY + EDISP 869 870C 871c we need to divide the solvent calls into mmpcm and not mmpcm 872c since we in the latter iterate between the mm and pcm for a given QM 873c density 874 875 ESOLT = D0 876 IF (.NOT. MMPCM) THEN ! do as we usually do ... 877 IF (PCM .AND. (FLAG(16) .OR. QMMM .OR. QMNPMM .OR. QM3)) THEN 878 CALL QUIT("PCM not compatible with other solvent models") 879 END IF 880 IF (QM3 .AND. .NOT. DOCCSD) THEN 881 CALL MEMGET2('REAL','FCQM3',KFCSOL,NNBAST,WRK,KFREE,LFREE) 882 CALL DZERO(WRK(KFCSOL),NNBAST) 883 CALL QM3FCK(DAOSAV(1,1),DAOSAV(1,2), 884 & WRK(KFCSOL),ESOLT,ENSQM3,EPOQM3,EELEL, 885 & WRK(KFREE),LFREE,IPR_DIIS) 886 CALL DAXPY(NNBAST,D1,WRK(KFCSOL),1,FAOSAV(1,INDEVC,1),1) 887 CALL MEMREL('DISCTL.QM3FCK',WRK,1,KFCSOL,KFREE,LFREE) 888 ELSE 889 EELEL = 0.0D0 890 END IF 891 892 IF (QMMM) THEN 893 IF (MMTIME) DTIME = SECOND() 894 CALL MEMGET2('REAL','FQMMM',KFCSOL,NNBAST,WRK,KFREE,LFREE) 895 CALL DZERO(WRK(KFCSOL),NNBAST) 896 CALL QMMMFCK(DAOSAV(1,1),DAOSAV(1,2), 897 & WRK(KFCSOL),ESOLT, 898 & WRK(KFREE),LFREE,IPQMMM) 899 CALL DAXPY(NNBAST,D1,WRK(KFCSOL),1,FAOSAV(1,INDEVC,1),1) 900 CALL MEMREL('DIIS_CTL.QMMMFCK',WRK,1, 901 & KFCSOL,KFREE,LFREE) 902 IF (MMTIME) THEN 903 DTIME = SECOND() - DTIME 904 TMMFCK = TMMFCK + DTIME 905 END IF 906 END IF 907 908 IF (QMNPMM) THEN 909 910 idim = GETDIM_RELMAT(.TRUE.) 911 idim_side = idim**0.5 912 allocate(relay_matrix(idim_side, idim_side, 1)) 913 CALL COMP_RELMAT(relay_matrix,WRK(KFREE),LFREE) 914 CALL WRITE_RELMAT(relay_matrix) 915 916 IF (DOMMSUB.AND.DOMMPOL) THEN 917 CALL GETDIM_MMMAT(IDIM, .TRUE.) 918 allocate(mm_relay_matrix(idim)) 919 CALL COMP_MMRELMAT(mm_relay_matrix,WRK(KFREE),LFREE) 920 END IF 921 922 allocate(qmcmm_sol(nnbast)) 923 qmcmm_sol = 0.0d0 924 CALL QMNPMM_FOCK(DAOSAV(1,1),DAOSAV(1,2),relay_matrix, 925 & mm_relay_matrix, 926 & qmcmm_sol,ESOLT,WRK(KFREE),LFREE) 927 CALL DAXPY(NNBAST,1.0d0,qmcmm_sol,1,FAOSAV(1,INDEVC,1),1) 928 929 deallocate(relay_matrix) 930 if (allocated(mm_relay_matrix)) deallocate(mm_relay_matrix) 931 deallocate(qmcmm_sol) 932 933 END IF 934 935 IF (FLAG(16)) THEN 936 CALL MEMGET2('REAL','FCSOL',KFCSOL,NNBAST,WRK,KFREE,LFREE) 937 CALL SOLFCK(DAOSAV(1,1),DAOSAV(1,2), 938 & WRK(KFCSOL),DUMMY,TNLM,.FALSE.,ESOLT, 939 & WRK(KFREE),LFREE,IPR_DIIS) 940 CALL DAXPY(NNBAST,D1,WRK(KFCSOL),1,FAOSAV(1,INDEVC,1),1) 941 CALL MEMREL('DIIS_CTL.SOLFCK',WRK,1,KFCSOL,KFREE,LFREE) 942 943 ELSE IF (PCM) THEN 944 CALL MEMGET2('REAL','FCPCM',KFCSOL,NNBASX,WRK,KFREE,LFREE) 945 CALL PCMFCK(DAOSAV(1,1),DAOSAV(1,2),WRK(KFCSOL), 946 & DUMMY,.FALSE.,ESOLT,WRK(KFREE), 947 & LFREE,IPR_DIIS) 948! write(lupri,*) 'pcm fock matrix' 949! call outpak(wrk(kfcsol), nbast, 1, lupri) 950 951 CALL DAXPY(NNBAST,DM1,WRK(KFCSOL),1,FAOSAV(1,INDEVC,1),1) 952 CALL MEMREL('DIIS_CTL.PCMFCK',WRK,KFRSAV, 953 & KFCSOL,KFREE,LFREE) 954 END IF 955 ELSE ! the mmpcm interface is only implemented for hf/dft so 956C no reference to doccsd 957 IF (.NOT. QM3 .OR. .NOT. PCM) THEN 958 CALL QUIT("MMPCM have to be used with both PCM and QM3") 959 END IF 960C parameters on input: MXITMP and THRSMP 961 LCONV = .FALSE. 962 XTEST1 = D0 963 CALL MEMGET2('REAL','FCSOL',KFCSOL,NNBAST,WRK,KFREE,LFREE) 964 CALL MEMGET2('REAL','MMPCM',KMMPCM,NNBAST,WRK,KFREE,LFREE) 965 DO ITER = 1, MXITMP 966C First the MM part ... 967 CALL DZERO(WRK(KFCSOL),NNBAST) 968 CALL DZERO(WRK(KMMPCM),NNBAST) 969 CALL QM3FCK(DAOSAV(1,1),DAOSAV(1,2), 970 & WRK(KFCSOL),ESOLT,ENSQM3,EPOQM3,EELEL, 971 & WRK(KFREE),LFREE,IPR_DIIS) 972 CALL DCOPY(NNBAST,WRK(KFCSOL),1,WRK(KMMPCM),1) 973 974 TEMP = ESOLT 975 976C Now the PCM part ... 977 978 CALL DZERO(WRK(KFCSOL),NNBAST) 979 CALL PCMFCK(DAOSAV(1,1),DAOSAV(1,2),WRK(KFCSOL), 980 & DUMMY,.FALSE.,ESOLT,WRK(KFREE), 981 & LFREE,IPR_DIIS) 982 CALL DAXPY(NNBAST,DM1,WRK(KFCSOL),1,WRK(KMMPCM),1) 983 984 ESOLT = ESOLT + TEMP 985 XTEST2 = ESOLT 986 IF ( ((ABS(XTEST1-XTEST2)) .LT. THRSMP) .OR. LOSPC) THEN 987 LCONV = .TRUE. 988 CALL DAXPY(NNBAST,D1,WRK(KMMPCM),1,FAOSAV(1,INDEVC,1),1) 989 write(lupri,*) 'MMPCM converged for given density' 990 write(lupri,*) 'No. of iterations: ',ITER 991 GOTO 999 992 ELSE 993 XTEST1 = XTEST2 994 END IF 995 END DO ! ITER = 1, MXITMP 996 999 CONTINUE 997 IF (.NOT. LCONV) write(lupri,*) 'MMPCM not converged !' 998 CALL MEMREL('DISCTL.PCMQM3FCK',WRK,1,KFCSOL,KFREE,LFREE) 999 END IF 1000C 1001 IF (USE_PELIB()) THEN 1002 IF (FLAG(16) .OR. PCM .OR. QMMM .OR. QMNPMM .OR. QM3) THEN 1003 CALL QUIT('.PELIB cannot be used with other embeddings') 1004 END IF 1005 ALLOCATE(PE_DENMAT(NNBASX), PE_FCKMAT(NNBASX)) 1006 IF (NASHT == 0 .OR. HSROHF) THEN 1007 PE_DENMAT = DAOSAV(:,1) 1008 ELSE 1009 PE_DENMAT = DAOSAV(:,1) + DAOSAV(:,2) 1010 END IF 1011 CALL PELIB_IFC_FOCK(PE_DENMAT, PE_FCKMAT, PE_ENERGY) 1012 ESOLT = PE_ENERGY 1013 FAOSAV(:,INDEVC,1) = FAOSAV(:,INDEVC,1) + PE_FCKMAT 1014 DEALLOCATE(PE_DENMAT, PE_FCKMAT) 1015 END IF 1016C 1017#ifdef HAS_PCMSOLVER 1018 if (pcm_cfg%do_pcm) then 1019! Allocate scratch space for the PCM Fock matrix contribution in AO basis. 1020 allocate(fock_pcm_ao(nnbasx)) 1021 fock_pcm_ao = 0.0d0 1022! Allocate scratch space for the density matrix to be passed to the PCM 1023! subroutines 1024 allocate(density_ao_symm(nnbasx)) 1025 density_ao_symm = 0.0d0 1026! Allocate the PCM Fock matrix contribution in AO basis. 1027! It is of dimension NNBASX and we take the memory from the WRK array. 1028! We zero it out afterwards. 1029! call memget('REAL', kfcsol, nnbasx, wrk, kfree, lfree) 1030! call dzero(wrk(kfcsol), nnbasx) 1031! DAOSAV(:, 1) = DCAO (density matrix for closed shells) 1032! DAOSAV(:, 2) = DVAO (density matrix for open shells) 1033! PCM doesn't care about closed/open shells, so we decide at this 1034! higher level which one of them to pass, based on the wavefunction asked. 1035 if (nasht == 0 .or. (nasht > 0 .and. dodft) .or. hsrohf) then 1036 density_ao_symm(:) = daosav(:, 1) 1037! We pass just DCAO, i.e. density_matrix = DCAO 1038 else 1039! We sum DCAO and DVAO, i.e. density_matrix = DCAO + DVAO 1040 density_ao_symm(:) = daosav(:, 1) + daosav(:, 2) 1041 end if 1042! Call the energy driver with the symmetry-unpacked density matrix 1043 if (nsym > 1) then 1044 allocate(density_ao(nnbasx)) 1045 density_ao = 0.0d0 1046! Symmetry un-pack density_ao_symm to density_ao 1047 call pksym1(density_ao, density_ao_symm, nbas, nsym, -1) 1048 esolt = pcm_scf_driver(density_ao, fock_pcm_ao, 1049 & wrk(kfree), lfree) 1050 else 1051 esolt = pcm_scf_driver(density_ao_symm, fock_pcm_ao, 1052 & wrk(kfree), lfree) 1053 end if 1054! Copy the PCM Fock matrix contribution on top of the AO basis Fock matrix 1055 if (nsym > 1) then 1056 allocate(fock_pcm_ao_symm(nnbasx)) 1057 fock_pcm_ao_symm = 0.0d0 1058! Symmetry-pack fock_pcm_ao to fock_pcm_ao_symm 1059 call pksym1(fock_pcm_ao, fock_pcm_ao_symm, nbas, nsym, +1) 1060 call daxpy(nnbast, -1.0d0, fock_pcm_ao_symm, 1, 1061 & faosav(1, indevc, 1), 1) 1062 else 1063 call daxpy(nnbast, -1.0d0, fock_pcm_ao, 1, 1064 & faosav(1, indevc, 1), 1) 1065 end if 1066! We can print it here, if needed. 1067! write(lupri,*) 'extpcm fock matrix' 1068! if (nsym > 1) then 1069! call outpak(fock_pcm_ao_symm, nbast, 1, lupri) 1070! else 1071! call outpak(fock_pcm_ao, nbast, 1, lupri) 1072! endif 1073 1074! Clean-up 1075 deallocate(fock_pcm_ao) 1076 deallocate(density_ao_symm) 1077 if (nsym > 1) then 1078 deallocate(fock_pcm_ao_symm) 1079 deallocate(density_ao) 1080 end if 1081 end if 1082#endif 1083 1084C 1085 EMCDIF = EMCSCF 1086 EMCSCF = EMY + EACTIV + EPOT + ESOLT 1087 1088 IF (QM3 .AND. (.NOT. OLDTG)) EMCSCF = EMCSCF - EELEL 1089C this also includes the case of MMPCM 1090 1091 EGSAV(1,INDEVC) = EMCSCF 1092 1093 EMCLOW = MIN(EMCLOW,EMCSCF) 1094 EMCDIF = EMCSCF - EMCDIF 1095 EMCDIFERR = THRINC*ABS(EMCLOW) 1096 GRDNRM_old = GRDNRM 1097C 1098 CALL DIIS_EVC(SMOAO,DAOSAV(1,1),DAOSAV(1,2),H1AO, 1099 & FAOSAV(1,INDEVC,1),FAOSAV(1,INDEVC,2),CMO1, 1100 & EVCSAV(1,INDEVC),WRK,KFREE,LFREE,IPR_DIIS) 1101C CALL DIIS_EVC(SMOAO,DCAO,DVAO,H1AO,FDAO,FVAO,CMO1, 1102C & ERRVEC,WRK,KFRSAV,LFRSAV,IPR_DIIS) 1103 IF (IPR_DIIS .GE. 20) THEN 1104 WRITE (LUPRI,*) ' Test output of all of EVCSAV, newest is', 1105 & INDEVC 1106 CALL OUTPUT(EVCSAV,1,NNORBT,1,NEVC,NNORBT,NEVC,-1,LUPRI) 1107 ELSE IF (IPR_DIIS .GE. 10) THEN 1108 WRITE (LUPRI,*) ' Error vector, DIIS iteration',ITDIIS 1109 CALL OUTPKB(EVCSAV(1,INDEVC),NORB,NSYM,-1,LUPRI) 1110 END IF 1111C 1112C Project out .FREEZE'ed orbitals from the gradient matrix 1113C 1114 IF (LNOROT) CALL PHPT(WRK(KPFROZ),WRK(KTMP1),NBAST) 1115 GRDNRM = DNRM2(NNORBT,EVCSAV(1,INDEVC),1) 1116 GRD_CHANGE = GRDNRM_old / GRDNRM 1117 EGSAV(2,INDEVC) = GRDNRM 1118 IF (LUW4 .NE. LUPRI .OR. IPR_DIIS .GT. 1) THEN 1119 WRITE (LUPRI,'(/A/A,I4,1P,G22.12,D15.5,D12.2,I5,L5)') 1120 &' ITDIIS, energy, error norm, change in energy, scr.thr., difden:' 1121 &,'@ DIIS iter.',ITDIIS,EMCSCF,GRDNRM,EMCDIF,-IFTHRS,DIFDEN 1122 IF (AUTOCC) THEN 1123 WRITE (LUPRI,'(/A,8I4)') 1124 & '@ AUTO OCC: current SCF occupation',(NISH(I),I=1,NSYM) 1125 IF (IOPRHF .NE. 0) WRITE (LUPRI,'(A,I4)') 1126 & '@ AUTO OCC: current open shell symmetry',IOPRHF 1127 END IF 1128 CALL FLSHFO(LUPRI) 1129 END IF 1130 IF (LUW4 .NE. LUPRI .OR. IPR_DIIS .LE. 1) THEN 1131#ifdef HAS_PCMSOLVER 1132 IF (EMBEDDING .or. pcm_cfg%do_pcm) then 1133#else 1134 IF (EMBEDDING) THEN 1135#endif 1136C this also includes the case of MMPCM 1137 IF (.NOT.AUTOCC) THEN 1138 WRITE (LUW4,'(A,I3,1P,2G20.12,D15.5,D12.2)') 1139 & '@',ITDIIS,EMCSCF,ESOLT,GRDNRM,EMCDIF 1140 ELSE IF (IOPRHF .EQ. 0) THEN 1141 WRITE (LUW4,'(A,I3,1P,2G20.12,2D11.2,2X,8I4)') 1142 & '@',ITDIIS,EMCSCF,ESOLT,GRDNRM,EMCDIF, 1143 & (NISH(I),I=1,NSYM) 1144 ELSE 1145 WRITE (LUW4,'(A,I3,1P,2G20.12,2D11.2,I4,2X,8I4)') 1146 & '@',ITDIIS,EMCSCF,ESOLT,GRDNRM,EMCDIF, 1147 & IOPRHF,(NISH(I),I=1,NSYM) 1148 END IF 1149 ELSE 1150 IF (.NOT.AUTOCC) THEN 1151 WRITE (LUW4,'(A,I3,1P,G22.12,D15.5,D12.2,I5)') 1152 & '@',ITDIIS,EMCSCF,GRDNRM,EMCDIF,NEVC 1153 ELSE IF (IOPRHF .EQ. 0) THEN 1154 WRITE (LUW4,'(A,I3,1P,G20.12,2D11.2,2X,8I4)') 1155 & '@',ITDIIS,EMCSCF,GRDNRM,EMCDIF, 1156 & (NISH(I),I=1,NSYM) 1157 ELSE 1158 WRITE (LUW4,'(A,I3,1P,G20.12,2D11.2,I4,2X,8I4)') 1159 & '@',ITDIIS,EMCSCF,GRDNRM,EMCDIF, 1160 & IOPRHF,(NISH(I),I=1,NSYM) 1161 END IF 1162 END IF 1163 CALL FLSHFO(LUW4) 1164 END IF 1165 IF (MOLDEN) CALL MOLDEN_SCFCON(ITDIIS,EMCSCF,.FALSE.) 1166 1167 IF (GRDNRM .LE. THDIIS) THEN ! converged 1168 1169 ICONV = 1 1170 1171C No level shift at convergence, otherwise orbital energies 1172C for e.g. CC will not be correct 1173 SHFTLVL = 0.0D0 1174 1175 IF (MOLDEN) CALL MOLDEN_SCFCON(ITDIIS,EMCSCF,.TRUE.) 1176 1177 WRITE (LUPRI,4110) ITDIIS,EMCSCF,GRDNRM 1178 IF (LUW4 .NE. LUPRI) WRITE (LUW4,4110) ITDIIS,EMCSCF,GRDNRM 1179 IF (IPR_DIIS .GE. -1) THEN 1180 IF (CHOINT .AND. (IPR_DIIS .GE. 3)) THEN 1181 ! Cholesky 1182 WRITE(LUPRI,4116) 'Coulomb part ',CSIRF 1183 WRITE(LUPRI,4116) 'Exchange part',XSIRF 1184 WRITE(LUPRI,4116) 'Cholesky I/O ',RSIRF 1185 WRITE(LUPRI,4116) 'MO transform.',OSIRF 1186 WRITE(LUPRI,4116) 'MO sorting ',SSIRF 1187 END IF 1188 WRITE(LUPRI,4115) TSIRF 1189 IF (QMMM) CALL QMMMTIMES('SIRIUS') 1190 IF (LUW4 .NE. LUPRI) THEN 1191 IF (CHOINT .AND. (IPR_DIIS .GE. 3)) THEN 1192 WRITE(LUPRI,4116) 'Coulomb part ',CSIRF 1193 WRITE(LUPRI,4116) 'Exchange part',XSIRF 1194 WRITE(LUPRI,4116) 'Cholesky I/O ',RSIRF 1195 WRITE(LUPRI,4116) 'MO transform.',OSIRF 1196 WRITE(LUPRI,4116) 'MO sorting ',SSIRF 1197 END IF 1198 WRITE(LUW4,4115) TSIRF 1199 END IF 1200 END IF 1201C Save converged CMO 1202 WRITE (CMOLBL,'(A4,I4)') 'DIIS',ITDIIS 1203 CALL NEWORB(CMOLBL,CMO,LNEWOCC) 1204C ... REWIT1 false: do not destroy any GEOWALK information 1205 GO TO 1000 1206 ELSE IF (ITDIIS .GE. MXDIIS) THEN 1207 ICONV = 0 1208 WRITE (LUPRI,4120) 1209 IF (LUW4 .NE. LUPRI) WRITE (LUW4,4120) 1210C 1211C Too many DIIS iterations, save final CMO before exit 1212 WRITE (CMOLBL,'(A4,I4)') 'DIIS',ITDIIS 1213 CALL NEWORB(CMOLBL,CMO,LNEWOCC) 1214C ... REWIT1 false: do not destroy any GEOWALK information 1215 GO TO 1100 1216 END IF 1217 4110 FORMAT(/'@ *** DIIS converged in',I4,' iterations !', 1218 & /'@ Converged SCF energy, gradient:',F20.12,1P,D12.2) 1219 4115 FORMAT( ' - total time used in SIRFCK :' ,F18.2,' seconds') 1220 4116 FORMAT( ' - total time used for ',A13,':',F12.2,' seconds') 1221 4120 FORMAT(/ 1222 & 'WARNING !!! DIIS aborted because max DIIS iterations reached !') 1223C 1224Chj 990819: use DIFDEN more often 1225Chj was: IF (IFTHRS .GE. 12 .AND. .NOT. FLAG(49)) DIFDEN = .TRUE. 1226 NEWTHR = .FALSE. 1227 IF (DIRFCK .AND. .NOT. USRSCR) THEN 1228 ITHRS = INT(-LOG10(GRDNRM/SCRFAC))+3 1229C ... "+3" instead of "+1" because we must 1230C have sufficient accuracy for lowering the 1231C gradient norm with a factor 100 and for 1232C a significant contribution in CVEC later ... 1233C (comment inserted Feb 2001/Jul 2006/hjaaj) 1234 IF (DODIFDEN) THEN 1235 ILIM = 2 1236Chj-aug99: only NEWTHR with skip of 2 for DIFDEN 1237 IF (ITHRS .GE. IFTMAX-2) ITHRS = IFTMAX 1238Chj-aug99: for effective DIFDEN; note this test is OK w. ILIM.eq.2 1239 ELSE 1240 ILIM = 1 1241 END IF 1242 IF (.NOT. AOSAVE .AND. ABS(ITHRS-IFTHRS) .GE. ILIM) THEN 1243 NEWTHR = .TRUE. 1244 IFTHRS = MIN(IFTMAX,MAX(IFTMIN,ITHRS)) 1245 END IF 1246 END IF 1247Chj 990819/000103 new: 1248 DIFDEN = .NOT. NEWTHR .AND. DODIFDEN 1249C DIFDEN true if not new threshold and not "no difden" 1250C -- if we have tightened screening threshold we must 1251C calculate w/o DIFDEN or the error in the Fock matrix 1252C would correspond to the previous screening threshold. 1253C 1254#ifdef USE_OPT_SOLVERS 1255 CALL DIISWEIGHTS(METHOD,ENERG,H1AO,SAOUNP, 1256 & DAOSAV,FAOSAV,EVCSAV, 1257 & INDEVC,NEVC,GVEC,BMAT,CVEC,MXERRV, 1258 & WRK,KFREE,LFREE, WRK(KPFROZ)) 1259 1260#else /* i.e. do not USE_OPT-SOLVERS */ 1261C 1262C Test if convergence is stopping because of numerical problems 1263C (EMCDIF test should make sure we are in local region) 1264C 1265ckr IF (ITNOCC .GT. 3 .AND. ABS(EMCDIF) .LT. THDDEF .AND. 1266ckr & GRDNRM .GT. CNVFAC*GRDNRM_OLD .AND. .NOT. NEWTHR) THEN 1267ckr ICONV = 0 1268ckr WRITE (LUPRI,4140) 1269ckr IF (LUW4 .NE. LUPRI) WRITE (LUW4,4140) 1270ckr GO TO 1100 1271ckr END IF 1272ckr 4140 FORMAT(' DIIS aborted, convergence too slow !') 1273C 1274C hjaaj apr 2002 1275C - check if a previous iteration had a lower gradient and a higher energy, 1276C if yes: eliminate this entry and earlier entries, because they 1277C are not relevant any more and may cause slower convergence 1278C or - even worse - that the convergence is stalled. 1279C Reason: DIIS tries to minimize gradient norm, not energy. 1280C - hjaaj mar 2004: only remove entries with lower gradient and 1281C higher energy, not "earlier entries" 1282C 1283 J_removed = 0 1284 IF (INDEVC.EQ.NEVC .AND. ITDIIS .GT. 3 .AND. 1285 & GRDNRM .LE. 10.0D0) THEN 1286C ... no wrap around yet - INDEVC .lt. NEVC requires more programming 1287C ITDIIS/GRDNRM check: give initial oscillations a chance to settle down 1288C before checking ... experience shows this will often be advantageous 1289C (see below) /hjaaj 1290 IEND = NEVC - 1 1291 DO I = 1,IEND 1292C IF (I .NE. INDEVC) THEN 1293 IF (EGSAV(1,I) .GT. EMCSCF + EMCDIFERR .AND. 1294 & EGSAV(2,I) .LT. GRDNRM) THEN 1295C ... Oops! a previous iter with higher E and lower grd 1296 J_removed = J_removed + 1 1297 EGSAV(1,I) = D0 1298 END IF 1299C END IF 1300 END DO 1301C 1302C remove entries with higher E and lower grd, if any 1303C 1304 IF (J_removed .GT. 0) THEN 1305 WRITE(LUPRI,'(/A/A,I3,A)') 1306 & 'Info: SCF gradient has been lower than now,', 1307 & ' therefore',J_removed, 1308 & ' old iterations are removed from DIIS.' 1309 K = 0 1310 DO J = 1, NEVC 1311 IF (EGSAV(1,J) .NE. D0) THEN 1312 K = K + 1 1313 IF (K .NE. J) THEN 1314 EGSAV(1,K) = EGSAV(1,J) 1315 EGSAV(2,K) = EGSAV(2,J) 1316 CALL DCOPY(NNORBT,EVCSAV(1,J),1, 1317 & EVCSAV(1,K),1) 1318 CALL DCOPY(NNBAST,FAOSAV(1,J,1),1, 1319 & FAOSAV(1,K,1),1) 1320 IF (.NOT. ONLYFC) 1321 & CALL DCOPY(NNBAST,FAOSAV(1,J,2),1, 1322 & FAOSAV(1,K,2),1) 1323 IF (J .LT. NEVC) THEN 1324C ... the new entry (i.e. J.EQ.NEVC) 1325C is added in DIIS_RED below 1326 KROW = K*(K-1) / 2 1327 DO L = 1,K 1328 BMAT(KROW+L) = 1329 & DDOT(NNORBT,EVCSAV(1,L),1,EVCSAV(1,K),1) 1330 END DO 1331 END IF 1332 END IF 1333 END IF 1334 END DO 1335 JTDIIS = K 1336 NEVC = K 1337 INDEVC = K 1338 END IF 1339 END IF 1340C 1341 CALL DIIS_RED(C2DIIS,BMAT,INDEVC,NEVC,DAMP,THREVC, 1342 & CVEC,XLMBDA,EVCSAV,NNORBT, 1343 & WRK,KFREE,LFREE,IPR_DIIS) 1344C CALL DIIS_RED(C2DIIS,BMAT,INDEVC,NEVC,DAMP,THREVC, 1345C & CVEC,XLMBDA,ERRVEC,L_ERRVEC, 1346C & WRK,KFRSAV,LFRSAV,IPR_DIIS) 1347C 1348C 1349C Check to see if DIIS iterations are stalled, 1350C if "yes" then restart DIIS /Feb-2001 hjaaj 1351C JTDIIS check: we need at least two Fock matrices to do this ... 1352C ITDIIS/GRDNRM check: give initial oscillations a chance to settle down 1353C before checking ... experience shows this will often be advantageous 1354C (no. 3 will often be lower in energy than no. 1, even if no. 2 1355C is higher in energy with lower gradient, especially if MOSTART H1DIAG) 1356C Note that we only "stall" if energy goes up while gradient goes 1357C down or vice versa because DIIS minimizes gradient and it is 1358C OK wrt. DIIS if both goes up or both goes down /Feb-2004 hjaaj 1359! Mar. 2011 hjaaj: if abs(EMCDIF) .lt. 1.0D-5 the error can be grid related in DFT 1360 IF (DODFT .AND. ABS(EMCDIF) .LT. 1.0D-5) THEN 1361 EMCDIF_test = 0.0D0 1362 ELSE IF (ABS(EMCDIF) .LT. 1.0D-9) THEN 1363 ! error can be related to integral accuracy or real*8 accuracy 1364 EMCDIF_test = 0.0D0 1365 ELSE 1366 EMCDIF_test = EMCDIF 1367 END IF 1368C Do not try that with open-shell - it does not work./pawsa 1369C Nov. 2018 jmho: disable again for open-shell because it really 1370C does not work. 1371 IF (ONLYFC .AND. JTDIIS.GT.1 .AND. ITDIIS.GT.3 .AND. 1372 & J_removed.EQ.0 .AND. 1373 & GRDNRM .LE. 10.0D0 .AND. 1374 & ( ABS(CVEC(INDEVC)) .LT. 0.1D0 .OR. 1375 & ABS(CVEC(INDEVC)) .GT. 5.0D0 .OR. 1376 & EMCDIF_test*(GRDNRM-GRDNRM_OLD) .LT. -1.0D-10 ) 1377 & .AND. BCKSTP) THEN 1378C .. test EMCDIF*grddif against -1.0D-10 instead of 0.0D0 1379C to avoid round-off problems here /hjaaj Sep 2005 1380 1381 IF (I_STALL .LT. 1) THEN ! do not stop first time energy goes up 1382 I_STALL = I_STALL + 1 1383 GO TO 300 1384 END IF 1385 I_STALL = 0 1386 1387 WRITE(LUPRI,'(1P/A,2(/A,G11.3),/A,I3,A)') 1388 & '!!! Info: DIIS restarted because it was stalled ...', 1389 & ' - energy changed by ',EMCDIF, 1390 & ' - gradient changed by',GRDNRM - GRDNRM_OLD, 1391 & ' - or strange C vector coeff. for current index (=', 1392 & INDEVC,') :' 1393 WRITE (LUPRI,'(10F12.6)') (CVEC(I),I=1,NEVC) 1394 IF (EMCDIF .LE. ABS(EMCSCF)*1.0D-12) THEN 1395 JNDEVC = INDEVC 1396 ELSE 1397C ... use previous Fock matrix if energy increased 1398 WRITE(LUPRI,'(/A)') '! Backstep to previous Fock'// 1399 & ' matrix because energy increased.' 1400 JNDEVC = MOD(JTDIIS-2,MXERRV) + 1 1401C 1402C restore old CMO from LUIT1 (for backstep); 1403C no DIFDEN because DAOSAV is for current, not previous iteration 1404C 1405 CALL READMO(CMO,9) 1406 DIFDEN = .FALSE. 1407 END IF 1408 IF (JNDEVC .NE. 1) THEN 1409 CALL DCOPY(NNORBT,EVCSAV(1,JNDEVC),1, 1410 & EVCSAV(1,1 ),1) 1411 CALL DCOPY(NNBAST,FAOSAV(1,JNDEVC,1),1, 1412 & FAOSAV(1,1 ,1),1) 1413 IF (.NOT. ONLYFC) 1414 & CALL DCOPY(NNBAST,FAOSAV(1,JNDEVC,2),1, 1415 & FAOSAV(1,1 ,2),1) 1416 END IF 1417 EGSAV(1,1) = EGSAV(1,JNDEVC) 1418 EGSAV(2,1) = EGSAV(2,JNDEVC) 1419 CVEC(1) = D1 1420 JTDIIS = 1 1421 NEVC = 1 1422 INDEVC = 1 1423 BMAT(1) = BMAT( (JNDEVC*(JNDEVC+1))/2 ) 1424C also restrict orbital rotation with a level shift 1425c the formula is empirical EMCDIF is positive here. 1426 SHFTLVL = MAX(-12D0,SHFTLVL - MIN(4.0D0,ABS(EMCDIF)*10.D0)) 1427c and slow down the level shifting annealing... 1428 SHFTFAC = MIN(SHFTFAC + (1.D0-SHFTFAC)*0.25,0.85D0) 1429 ELSE IF (ITDIIS .GT. 1) THEN 1430 IF (GRD_CHANGE .GE. 2.0D0) SHFTLVL = SHFTLVL*SHFTFAC 1431c reduce level shift if convergence is good 1432 IF (GRD_CHANGE .GE. 4.0D0) SHFTLVL = SHFTLVL*SHFTFAC 1433c double reduce level shift if convergence is excellent 1434 IF (SHFTLVL .GE. -1.D-2) SHFTLVL = 0.0D0 1435 END IF 1436 300 CONTINUE 1437#endif /* USE_OPT-SOLVER */ 1438 1439C Level shift is only working for closed shell calculations ... 1440C (open shell problems reported by users ...) / hjaaj 22-Jun-2005 1441! IF ( NASHT .GT. 0 ) SHFTLVL = 0.0D0 1442! enabled again 6-Nov-2014 after testing a problematic open-shell DFT 1443! calculation which converged with level shift, but not without. /hjaaj 1444 1445C Iteration accepted; save current CMO for restart 1446 WRITE (CMOLBL,'(A4,I4)') 'DIIS',ITDIIS 1447 CALL NEWORB(CMOLBL,CMO,LNEWOCC) 1448C ... REWIT1 false: do not destroy any GEOWALK information 1449 LNEWOCC = .FALSE. 1450 1451C print atomic populations in each iteration ? /hjaaj 1452 if (lim_poppri .gt. 0) 1453 & call sirpop('DIIS ',DV,wrk(kfree),lfree) 1454 1455C 1456 IF (HSROHF) THEN 1457C 1458C Update FCAO with the FV correction in AO basis 1459C 1460 CALL DAXPY(NNBAST,D1,H1AO,1,FAOSAV(1,INDEVC,1),1) 1461 CALL DCOPY(NNBAST,FAOSAV(1,INDEVC,1),1,FCAO,1) 1462 CALL FVCORR(1,D1,DAOSAV,FAOSAV(1,INDEVC,1),MXERRV,FCAO, 1463 & WRK,KFREE,LFREE ) 1464C 1465C Update the latest FAOSAV with the effective fock matrix 1466C 1467 CALL DCOPY(NNBAST,FCAO,1,FAOSAV(1,INDEVC,1),1) 1468 CALL DZERO(FCAO,NNBAST) 1469 ELSE 1470 CALL DCOPY(NNBAST,H1AO,1,FCAO,1) 1471 END IF 1472 CALL DGEMM('N','N',NNBAST,1,NEVC,1.D0, 1473 & FAOSAV(1,1,1),NNBAST, 1474 & CVEC,NEVC,1.D0, 1475 & FCAO,NNBAST) 1476 IF (.NOT. ONLYFC .AND..NOT.HSROHF) THEN 1477 CALL DGEMM('N','N',NNBAST,1,NEVC,1.D0, 1478 & FAOSAV(1,1,2),NNBAST, 1479 & CVEC,NEVC,0.D0, 1480 & FCAO(1+NNBAST),NNBAST) 1481 END IF 1482C 1483C Transform FDAO to FD (FD saved in DCAO). 1484C Find FV and make open-shell Fock matrix 1485C 1486 CALL UTHUB(FCAO,DCAO,CMO,WRK(KFREE),NSYM,NBAS,NORB) 1487 IF (HSROHF .OR. ONLYFC) THEN 1488 IF (IPR_DIIS .GE. 15) THEN 1489 IF (HSROHF) THEN 1490 WRITE (LUPRI,'(/A)') 1491 & ' Fock matrix after FV correction' 1492 ELSE 1493 WRITE (LUPRI,'(/A)') 1494 & ' Fock matrix before diagonalization' 1495 END IF 1496 CALL OUTPKB(DCAO,NORB,NSYM,-1,LUPRI) 1497 END IF 1498 ELSE 1499 DO 330 ISYM = 1,NSYM 1500 IF (NASH(ISYM) .EQ. 0) GO TO 330 1501 IF (IPR_DIIS .GT. 15) THEN 1502 WRITE (LUPRI,'(/A,I3)') 1503 & ' Fock matrix before FV correction, symmetry',ISYM 1504 CALL OUTPAK(DCAO(1+IIORB(ISYM)),NORB(ISYM),-1,LUPRI) 1505 END IF 1506 CALL UTHU(FCAO(1+NNBAST+IIBAS(ISYM)),DCAO(1+NNBAST), 1507 & CMO(1+ICMO(ISYM)),WRK(KFREE),NBAS(ISYM),NORB(ISYM)) 1508 DO JACT = 1, NASH(ISYM) 1509 KACT = NISH(ISYM) + JACT 1510 KROW = KACT*(KACT-1)/2 1511 DO 310 J = 1, NISH(ISYM) 1512 JFKJ = IIORB(ISYM) + KROW + J 1513 DCAO(JFKJ) = DCAO(JFKJ) + DCAO(NNBAST+KROW+J) 1514 310 CONTINUE 1515 DO 320 J = NISH(ISYM)+NASH(ISYM) + 1,NORB(ISYM) 1516 JROW = J*(J-1)/2 1517 JFJK = IIORB(ISYM) + JROW + KACT 1518 DCAO(JFJK) = DCAO(JFJK) - DCAO(NNBAST+JROW+KACT) 1519 320 CONTINUE 1520 END DO 1521 IF (IPR_DIIS .GE. 15) THEN 1522 WRITE (LUPRI,'(/A,I3)') 1523 & ' Fock matrix after FV correction, symmetry',ISYM 1524 CALL OUTPAK(DCAO(1+IIORB(ISYM)),NORB(ISYM),-1,LUPRI) 1525 END IF 1526 330 CONTINUE 1527 END IF 1528C 1529C If requested (by SHFTLVL .ne. 0) : 1530C Perform "restricted step" step reduction (i.e. level shift) 1531C by adding 0.5*SHFTLVL*DENS(i,i) to all diagonal elements; 1532C /hjaaj-Feb.2001 1533 1534 IF (SHFTLVL .NE. D0) THEN 1535 WRITE(LUPRI,'(I4,A,1P,D10.2)') ITDIIS,' Level shift: '// 1536 & 'doubly occupied orbital energies shifted by',SHFTLVL 1537 IF (NASHT .GT. 0) WRITE(LUPRI,'(T16,A,1P,D10.2)')'and '// 1538 & 'singly occupied orbital energies shifted by',0.5D0*SHFTLVL 1539 END IF 1540chjaaj: improve comment with something about "restricted step" 1541chjaaj: and "0.5*occ(i)* " /23may09 hjaaj 1542C 1543C Diagonalize (level-shifted) FC (saved in DCAO) 1544C 1545 CALL ICOPY(8,NISH,1,NISH_save,1) 1546 JOPRHF = IOPRHF 1547 IF (GRDNRM_old > 0.0D0) THEN 1548 GRDNRM_factor = MAX(12.0D0, GRDNRM_old/GRDNRM) 1549 ELSE 1550 GRDNRM_factor = 12.0D0 1551 END IF 1552 IF ( GRDNRM .LE. MAX(1.D-5, THDIIS*GRDNRM_factor) ) THEN 1553 THR_FCKEIG = 0.8D0*THDIIS 1554 CALL ICOPY(8,NORB,1,NCAN,1) ! we need all orbitals canonical for CC 1555 ELSE 1556 THR_FCKEIG = 0.01D0*GRDNRM 1557 CALL ICOPY(8,NOCC,1,NCAN,1) ! we only need oocupied orbitals canonical in global region 1558 END IF 1559 1560 call gettim(t1,w1) 1561 CALL FCKEIG(CMO,DCAO,SHFTLVL,THR_FCKEIG,NCAN,IPR_DIIS, 1562 & WRK(KFREE),LFREE) 1563 IF (IPR_DIIS .GT. 2) THEN 1564 call gettim(t2,w2) 1565 write(lupri,'(A,1P,3D10.2)') 1566 & ' CPU and wall times used in FCKEIG, THR_FCKEIG: ', 1567 & t2-t1,w2-w1,thr_fckeig 1568 END IF 1569 1570C Test if AUTOCC has changed occupation : 1571 1572 IF (AUTOCC) THEN 1573 NTEST = 0 1574 DO ISYM = 1,8 1575 NTEST = NTEST + ABS(NISH(ISYM)-NISH_save(ISYM)) 1576 END DO 1577 NTEST = NTEST + ABS(IOPRHF-JOPRHF) 1578 IF (NTEST .GT. 0) THEN 1579 NEWOCC = NEWOCC + 1 1580 LNEWOCC = .TRUE. 1581 ITNOCC = 0 1582 IF (NEWOCC .GT. 6) THEN 1583 WRITE (LUPRI,4210) NEWOCC 1584 IF (LUW4 .NE. LUPRI) WRITE (LUW4,4210) NEWOCC 1585 CALL QUIT( 1586 & 'FATAL ERROR: problems with automatic SCF occupation') 1587 END IF 1588Chj-sep99: bugfix for open shell when open shell symmetry changed 1589C (NASH(ISYM) used for test in DIIS_CTL !!!) 1590 IF (NASHT .EQ. 1) THEN 1591 CALL IZERO(NASH,8) 1592 NASH(IOPRHF) = 1 1593 END IF 1594Chj-sep99: and change index arrays to new occupation (for RHFENR and ?) 1595 CALL SETORB 1596 END IF 1597 END IF 1598 4210 FORMAT(/' DIIS has changed occupation numbers',I2,' times now.' 1599 & /' Program aborts because this indicates problems with'// 1600 & ' automatic occupation.') 1601 1602C --> Next DIIS iteration 1603 GO TO 100 1604C 1605C 1000: converged 1606C 1100: not converged 1607 1000 CONTINUE 1608C 1609C QCHF currently needed for solvent and for writing SIRIFC: 1610C reset to not converged in these cases. 951130-hjaaj. 1611C We write to SIRIFC here instead, thus no need to reset ICONV,K.Ruud-May97 1612C hjaaj aug99: revised code for SIRIFC to fix some potential problems, 1613C and call tractl if requested with ITRFIN 1614C 1615C If converged (ICONV .eq. 1) and if this is final level of wave function: 1616C 1617C Print orbital energy analysis for RHF if QCHF not called 1618C because ICONV = 1 and if either IPR_DIIS .gt. 0 or 1619C this is final level of wave function. 1620C 1621C 1622C 1623C Check if we need new RHFWOP because of AUTOCC 1624C (inserted by Kenneth Ruud May 95, open shell hjaaj Aug 95) 1625C Moved from SIRCTL to here since we do not necessarily use 1626C quadratic HF after the DIIS anymore.... 1627C 1628 RHFWOP = .TRUE. 1629 IF (IOPRHF .NE. IOPRHF_old) RHFWOP = .FALSE. 1630 DO ISYM = 1, NSYM 1631 IF (NRHF(ISYM) .NE. NISH(ISYM)) THEN 1632 NRHF(ISYM) = NISH(ISYM) 1633 RHFWOP = .FALSE. 1634 END IF 1635 END DO 1636 IF (.NOT.RHFWOP) THEN 1637 OLDWOP = .FALSE. 1638C 1639 IF (.NOT.DOMC) WRINDX = .TRUE. 1640C ... We need to correct orbital rotation information on file 1641C for response and ABACUS 1642C if occupation has changed and RHF determines orbitals 1643C 1644 JWOPSY = 1 1645 CALL SIRSET(WRK(KFREE),LFREE,OLDWOP) 1646 IAVERR = 0 1647 CALL AVECHK(IAVERR) 1648 IF (IAVERR .NE. 0) CALL QUIT( 1649 & 'SIRCTL DIIS error: inconsistency in sup.sym. after AUTOCC') 1650 RHFWOP = .TRUE. 1651 END IF 1652C 1653 IF ( ICONV .NE. 0. .AND. 1654 & ( IPR_DIIS .GT. 0 .OR. 1655 & CMOPRI .OR. 1656 & .NOT.(DOMP2 .OR. DOCI .OR. DOMC)) ) THEN 1657 IPRENR = 1 1658 ELSE 1659 IPRENR = 0 1660 END IF 1661 1662! IF ( ( FLAG(25) .OR. INERSI ) .AND. 1663! & (ICONV .NE. 0. .AND. .NOT. DOMC) ) THEN 1664 IF (ICONV .NE. 0) THEN 1665 ! hjaaj nov 2010: now always write SIRIFC if converged - so we 1666 ! have it available afterwards for post-SCF programs. This is more 1667 ! modular than making dependent on explicit checking for a 1668 ! post-SCF program, and we can use it in a restart of Dalton. 1669 IWRIFC = 1 1670 ELSE 1671 IWRIFC = 0 1672 END IF 1673C 1674 IF (IPRENR .GT. 0 .OR. IWRIFC .GT. 0) THEN 1675C 1676C Construct FDAO (FDAO saved in DCAO) 1677C Transform FDAO to FD (FD saved in FCAO). 1678C 1679 CALL DCOPY(NNBAST,FAOSAV(1,INDEVC,1),1,DCAO,1) 1680 CALL DAXPY(NNBAST,D1,H1AO,1,DCAO,1) 1681 CALL UTHUB(DCAO,FCAO,CMO,WRK(KFREE), 1682 & NSYM,NBAS,NORB) 1683 END IF 1684 1685C 1686C if qm3 .or. pcm then also mmpcm automatically 1687 IF (IPRENR .GT. 0) THEN 1688C 1689 WRITE (LUW4,'(//A)') ' *** SCF orbital energy analysis ***' 1690 IF (FLAG(16) .OR. PCM .OR. QM3 .OR. QMMM .OR. QMNPMM) 1691 & WRITE (LUW4,'(A)') ' (incl. solvent contribution)' 1692 IF (USE_PELIB()) WRITE(LUW4,'(A)') 1693 & ' (incl. contribution from polarizable embedding potential)' 1694 CALL RHFENR(IPRI4,LUW4,FCAO, WRK(KFREE),LFREE) 1695C------------------- 1696 IF (LUPRI .NE. LUW4) THEN 1697 WRITE (LUPRI,'(//A)') ' *** SCF orbital energy analysis ***' 1698 IF (FLAG(16) .OR. PCM .OR. QM3 .OR. QMMM .OR. QMNPMM) 1699 & WRITE (LUPRI,'(A)') ' (incl. solvent contribution)' 1700 IF (USE_PELIB()) WRITE(LUPRI,'(A)') 1701 & ' (incl. contribution from polarizable embedding potential)' 1702 CALL RHFENR(IPRI6,LUPRI,FCAO, WRK(KFREE),LFREE) 1703 END IF 1704C------------------- 1705C CALL RHFENR(IPRINT,LUPRI,FC,FV,SCRA,LSCRA) 1706C 1707 IF (IPR_DIIS .GT. 10) THEN 1708 WRITE(LUPRI,'(/A)') 1709 & 'Final SCF orbitals after DIIS iterations:' 1710 CALL PRORB(CMO,(IPR_DIIS .LT. 20),LUPRI) 1711 END IF 1712C 1713 END IF 1714 IF (USE_PELIB()) THEN 1715 ALLOCATE(PE_DENMAT(NNBASX)) 1716 IF (NASHT == 0 .OR. HSROHF) THEN 1717 PE_DENMAT = DAOSAV(:,1) 1718 ELSE 1719 PE_DENMAT = DAOSAV(:,1) + DAOSAV(:,2) 1720 END IF 1721 CALL PELIB_IFC_ENERGY(PE_DENMAT) 1722 DEALLOCATE(PE_DENMAT) 1723 END IF 1724 IF (PELIB_IFC_DO_MEP()) THEN 1725 ALLOCATE(PE_DENMAT(NNBASX)) 1726 IF (NASHT == 0 .OR. HSROHF) THEN 1727 PE_DENMAT = DAOSAV(:,1) 1728 ELSE 1729 PE_DENMAT = DAOSAV(:,1) + DAOSAV(:,2) 1730 END IF 1731 CALL PELIB_IFC_MEP(PE_DENMAT) 1732 DEALLOCATE(PE_DENMAT) 1733 END IF 1734 IF (PELIB_IFC_DO_SAVDEN()) THEN 1735 ALLOCATE(PE_DENMAT(NNBASX)) 1736 IF (NASHT == 0 .OR. HSROHF) THEN 1737 PE_DENMAT = DAOSAV(:,1) 1738 ELSE 1739 PE_DENMAT = DAOSAV(:,1) + DAOSAV(:,2) 1740 END IF 1741 CALL PELIB_IFC_SAVE_DENSITY(PE_DENMAT, 1742 & FCAO(1:NORBT*(NORBT+1)/2), 1743 & CMO(1:NBAST*NORBT)) 1744 DEALLOCATE(PE_DENMAT) 1745 END IF 1746C 1747C If (FLAG(25) .OR. INERSI) write to SIRIFC 1748C FLAG(25) is true for .ABACUS (geometry optimization and many properties) 1749C or .RESPONS (response calculation) 1750C or .INTERFACE (request for write of SIRIFC) 1751C or .WESTA 1752C or .CC 1753C INERSI is true if this is initial state calculation 1754C in a solvent calculation with inertial polarization 1755C 1756 IF ( IWRIFC .EQ. 1 ) THEN 1757 IF (IPRI6 .GE. 0) WRITE (LUPRI,'(/A)') 1758 & ' --- Writing SIRIFC interface file' 1759 IF (ICI0 .EQ. 6) THEN 1760Chj ... if (GEOWLK) then 1761 IF (.NOT. OPTNEW) WRITE (LUPRI,'(//A/)') 1762 & ' Check ratio for geometry walk with converged SCF :' 1763 JWSTEP = 1 1764 EMCOLD = EMCGEO 1765 DEPRED = DEPGEO 1766 CALL SIRSTP(JWSTEP,DUMMY,DUMMY,DUMMY,DUMMY,1) 1767Chj ... WLKREJ in gnrinf.h will be true if geom. step rejected. 1768 END IF 1769 1770 LPV = NNASHX*NNASHX 1771 CALL MEMGET2('REAL','PV',KPV,LPV,WRK,KFREE,LFREE) 1772 IF (NASHT .EQ. 1) THEN 1773Chj PV(1) is 0.0D0 for a single open shell RHF 1774 WRK(KPV) = D0 1775 ELSE IF (NASHT .GT. 1) THEN ! HSROHF 1776 ! we fill PV with a very large number, to make sure it gives 1777 ! problems if used later. As it is now, HSROHF is done 1778 ! completely in AO basis, thus PV (and FQ below) are not 1779 ! used. This filling makes it very clear that there are 1780 ! problems if someone programs use of PV for HSROHF without 1781 ! checking code. 1782 WRK(KPV:KPV+LPV-1) = -9.87654321D123 1783 END IF 1784 1785 CALL MEMGET2('REAL','GORB',KGORB,NWOPH,WRK,KFREE,LFREE) 1786 CALL DZERO(WRK(KGORB),NWOPH) 1787Chj ... set orbital gradient to zero (i.e. completely converged!) 1788 1789 LFQ = NASHT*NORBT 1790 CALL MEMGET2('REAL','FQ',KFQ,LFQ,WRK,KFREE,LFREE) 1791 IF (NASHT.GT.0) THEN 1792C 1. put FV in DCAO: 1793 IF (HSROHF .AND. .NOT.HSROHF_save) THEN 1794 ! calculate correct FV for interface SIRIFC 1795 ! (full FV has not been calculated before when HSROHF true, as 1796 ! FAOSAV(:,INDEVC,2) contains minus the exchange part of FVAO) 1797 ! Standard FV is necessary for NASHT.eq.1 because many places 1798 ! in ABACUS and RESPONS the test is on NASHT.EQ.1 instead of IOPRHF.gt.0 /141123-hjaaj 1799 1800 ! retrieve final DVAO in DCAO 1801 IF (NSYM.GT.1) THEN 1802 CALL PKSYM1(DCAO(1+N2BASX),DAOSAV(1,2), 1803 & NBAS,NSYM,-1) 1804 CALL DUNFLD(NBAST,DCAO(1+N2BASX),DCAO) 1805 ELSE 1806 CALL DUNFLD(NBAST,DAOSAV(1,2),DCAO) 1807 ENDIF 1808 DCAO(1:N2BASX) = -DCAO(1:N2BASX) 1809 HSROHF = .FALSE. 1810 CALL FCK2AO(.TRUE.,FCAO(1+N2BASX),DCAO,WRK,KFREE,LFREE) ! calculate FVAO in FCAO(1+N2BASX) 1811 ! ( first parameter .TRUE. tells FCK2AO only one density matrix ) 1812 CALL DGETSP(NBAST,FCAO(1+N2BASX),WRK(KFREE)) 1813 CALL PKSYM1(WRK(KFREE),FCAO(1+N2BASX),NBAS,NSYM,1) 1814 CALL UTHUB(FCAO(1+N2BASX),DCAO,CMO,WRK(KFREE), ! and now FV is in DCAO 1815 & NSYM,NBAS,NORB) 1816 1817 ELSE 1818 CALL UTHUB(FAOSAV(1,INDEVC,2),DCAO,CMO,WRK(KFREE), 1819 & NSYM,NBAS,NORB) 1820 END IF 1821 1822C 2. form FC' = FD - FV' in FCAO, both for standard FV'=FV and for HSROHF FV'=-FV_exch 1823 CALL DAXPY(NNORBT,DM1,DCAO,1,FCAO,1) 1824 1825C 3. FQ matrix is zero matrix for open shell RHF 1826 CALL DZERO(WRK(KFQ),NASHT*NORBT) 1827 END IF ! (NASHT.GT.0) 1828C 1829C *** LOCALIZATION 1830C 1831 IF (BOYORB .OR. PIPORB) CALL SIRLOC(CMO,WRK(KFREE),LFREE) 1832 IF (BOYSEL) CALL LOCCTL(CMO,WRK(KFREE),LFREE) 1833C 1834 IF (LMULBS .OR. R12CAL) THEN 1835C Construct orthonormal auxiliary basis set (WK/UniKA/04-11-2002). 1836 CALL TIMER('START ',TIMSTR,TIMEND) 1837 CALL R12AUX(WRK(KFREE),LFREE) 1838 CALL TIMER('R12AUX',TIMSTR,TIMEND) 1839 END IF 1840C 1841C *** and now we are ready to call WR_SIRIFC to write SIRIFC 1842 CALL MEMGET2('REAL','DUM',KDUM,0,WRK,KFREE,LFREE) 1843 CALL WR_SIRIFC(CMO,DV,WRK(KPV),FCAO,DCAO, 1844 & WRK(KFQ),WRK(KDUM),WRK(KDUM),WRK(KDUM), 1845 & WRK(KFREE),LFREE,WRK(KGORB),TNLM,.FALSE.,WRK(KDUM)) 1846C CALL WR_SIRIFC(CMO,DV,PV,FC,FV,FQ,CREF,FCAC,H2AC,WRK,LFREE, 1847C & GORB,ERLM,ORBHES,XINDX) 1848C 1849 CALL MEMREL('DIIS_CTL.WR_SIRIFC',WRK,1,KPV,KFREE,LFREE) 1850 END IF 1851C 1852C Transform integrals for abacus or response or something else, 1853C if converged and if not followed by higher level, 1854C and if requested in input (.FINAL TRANSFORMATION) 1855C 1856 IF (ICONV .NE. 0 .AND. ABS(ITRFIN) .LE. 10 .AND. 1857 & .NOT.(DOMP2 .OR. DOCI .OR. DOMC)) THEN 1858 IF (IPRI6 .GE. 0) WRITE (LUPRI,'(/A,I3)') 1859 & ' --- Transforming 2-el. integrals acc. to'// 1860 & ' .FINAL TRANSFORMATION =',ITRFIN 1861 JTRLVL = ITRFIN 1862 CALL TRACTL(JTRLVL,CMO,WRK(KFREE),LFREE) 1863 END IF 1864C 1865C *********************************************************** 1866C 1867 1100 CONTINUE 1868C restore original values 1869 IFTHRS = IFTHSV 1870 HSROHF = HSROHF_save 1871 CALL QEXIT('DIIS_CTL') 1872 RETURN 1873 END 1874C /* Deck DIIS_EVC */ 1875 SUBROUTINE DIIS_EVC(SMOAO,DCAO,DVAO,H1AO,FDAO,FVAO,CMO1, 1876 & ERRVEC,WRK,KFRSAV,LFRSAV,IPR_DIIS) 1877C 1878C 1879#include "implicit.h" 1880 DIMENSION DCAO(NNBAST), DVAO(NNBAST), FDAO(NNBAST), FVAO(NNBAST) 1881 DIMENSION H1AO(NNBAST) 1882 DIMENSION SMOAO(N2BAST), CMO1(NCMOT), ERRVEC(NNORBT), WRK(*) 1883C 1884C Used from common blocks: 1885C INFINP : LNOROT,NOROT(), HSROHF 1886C INFORB : N2BAST, NNBAST, NCMOT, NSYM, ..., NFROT, ... 1887C INFDIM : N2BASM 1888C 1889#include "gnrinf.h" 1890#include "maxorb.h" 1891#include "priunit.h" 1892#include "infinp.h" 1893#include "inforb.h" 1894#include "infdim.h" 1895#include "infpri.h" 1896#include "scbrhf.h" 1897C 1898 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, DM1 = -1.0D0) 1899C 1900 CALL QENTER('DIIS_EVC') 1901 KFREE = KFRSAV 1902 LFREE = LFRSAV 1903 CALL MEMGET2('REAL','TMP1',KTMP1,N2BASM,WRK,KFREE,LFREE) 1904 CALL MEMGET2('REAL','TMP2',KTMP2,N2BASM,WRK,KFREE,LFREE) 1905 CALL MEMGET2('REAL','TMP3',KTMP3,N2BASM,WRK,KFREE,LFREE) 1906 IF (HSROHF) THEN 1907 CALL DAXPY(NNBAST,D1,DVAO,1,DCAO,1) 1908 CALL DSCAL(NNBAST,DM1,DVAO,1) 1909 END IF 1910 DO 800 ISYM = 1,NSYM 1911 NOCCI = NOCC(ISYM) 1912 NORBI = NORB(ISYM) 1913 IF (NOCCI .EQ. 0) THEN 1914 IF (NORBI .GT. 0) THEN 1915 CALL DZERO(ERRVEC(1+IIORB(ISYM)),NNORB(ISYM)) 1916 END IF 1917 GO TO 800 1918 END IF 1919 NBASI = NBAS(ISYM) 1920 IF (IPR_DIIS .GE. 35) THEN 1921 WRITE (LUPRI,*) 'Debug output from DIIS_EVC, ISYM=',ISYM 1922 WRITE (LUPRI,*) 'NORBI,NBASI:',NORBI,NBASI 1923 WRITE (LUPRI,*) 'DIIS_EVC CMO1:' 1924 CALL OUTPUT(CMO1(1+ICMO(ISYM)),1,NBASI,1,NORBI, 1925 & NBASI,NORBI,-1,LUPRI) 1926 WRITE (LUPRI,*) 'DIIS_EVC SMOAO = 4 Ct SAO:' 1927 CALL OUTPUT(SMOAO(1+I2BAS(ISYM)),1,NORBI,1,NBASI, 1928 & NORBI,NBASI,-1,LUPRI) 1929 IF (IPR_DIIS .GT. 50) THEN 1930 CALL DGEMM('N','N',NORBI,NORBI,NBASI,1.D0, 1931 & SMOAO(1+I2BAS(ISYM)),NORBI, 1932 & CMO1(1+ICMO(ISYM)),NBASI,0.D0, 1933 & WRK(KTMP3),NORBI) 1934 WRITE (LUPRI,*) 'DIIS_EVC SMO = 4 Ct SAO C:' 1935 CALL OUTPUT(WRK(KTMP3),1,NORBI,1,NORBI, 1936 & NORBI,NORBI,-1,LUPRI) 1937 END IF 1938 END IF 1939 CALL DSPTSI(NBASI,FDAO(1+IIBAS(ISYM)),WRK(KTMP3)) 1940 CALL DSPTSI(NBASI,H1AO(1+IIBAS(ISYM)),WRK(KTMP2)) 1941 CALL DAXPY(N2BAS(ISYM),D1,WRK(KTMP2),1,WRK(KTMP3),1) 1942 IF (NISH(ISYM) .GT. 0) THEN 1943 CALL DUNFLD(NBASI,DCAO(1+IIBAS(ISYM)),WRK(KTMP2)) 1944 CALL DGEMM('N','N',NBASI,NBASI,NBASI,1.D0, 1945 & WRK(KTMP2),NBASI, 1946 & WRK(KTMP3),NBASI,0.D0, 1947 & WRK(KTMP1),NBASI) 1948 IF (IPR_DIIS .GE. 35) THEN 1949 WRITE (LUPRI,*) 'DIIS_EVC DCAO unfolded:' 1950 CALL OUTPUT(WRK(KTMP2),1,NBASI,1,NBASI, 1951 & NBASI,NBASI,-1,LUPRI) 1952 WRITE (LUPRI,*) 'DIIS_EVC FDAO:' 1953 CALL OUTPUT(WRK(KTMP3),1,NBASI,1,NBASI, 1954 & NBASI,NBASI,-1,LUPRI) 1955 WRITE (LUPRI,*) 'DIIS_EVC TMP1 = DCAO FDAO:' 1956 CALL OUTPUT(WRK(KTMP1),1,NBASI,1,NBASI, 1957 & NBASI,NBASI,-1,LUPRI) 1958 END IF 1959 ELSE 1960 CALL DZERO(WRK(KTMP1),NBASI*NBASI) 1961 END IF 1962 IF (NASH(ISYM) .GT. 0) THEN 1963 CALL DSPTSI(NBASI,FVAO(1+IIBAS(ISYM)),WRK(KTMP2)) 1964 CALL DAXPY(NBASI*NBASI,DM1,WRK(KTMP2),1,WRK(KTMP3),1) 1965 CALL DUNFLD(NBASI,DVAO(1+IIBAS(ISYM)),WRK(KTMP2)) 1966 CALL DGEMM('N','N',NBASI,NBASI,NBASI,1.D0, 1967 & WRK(KTMP2),NBASI, 1968 & WRK(KTMP3),NBASI,1.D0, 1969 & WRK(KTMP1),NBASI) 1970 IF (IPR_DIIS .GE. 35) THEN 1971 WRITE (LUPRI,*) 'DIIS_EVC DVAO unfolded:' 1972 CALL OUTPUT(WRK(KTMP2),1,NBASI,1,NBASI, 1973 & NBASI,NBASI,-1,LUPRI) 1974 WRITE (LUPRI,*) 'DIIS_EVC FCAO:' 1975 CALL OUTPUT(WRK(KTMP3),1,NBASI,1,NBASI, 1976 & NBASI,NBASI,-1,LUPRI) 1977 WRITE (LUPRI,*) 'DIIS_EVC TMP1 = DCAO FDAO + DVAO FCAO:' 1978 CALL OUTPUT(WRK(KTMP1),1,NBASI,1,NBASI, 1979 & NBASI,NBASI,-1,LUPRI) 1980 END IF 1981 END IF 1982 CALL DGEMM('N','N',NORBI,NBASI,NBASI,1.D0, 1983 & SMOAO(1+I2BAS(ISYM)),NORBI, 1984 & WRK(KTMP1),NBASI,0.D0, 1985 & WRK(KTMP2),NORBI) 1986 CALL DGEMM('N','N',NORBI,NORBI,NBASI,1.D0, 1987 & WRK(KTMP2),NORBI, 1988 & CMO1(1+ICMO(ISYM)),NBASI,0.D0, 1989 & WRK(KTMP1),NORBI) 1990C 1991C Zero rows and columns corresponding to frozen orbitals 1992C 1993 DO 220 I = 1,NFRO(ISYM) 1994 JOFF = KTMP1 - 1 + (I-1)*NORBI 1995 DO 210 J = 1,NORBI 1996 WRK(KTMP1-1+(I-1)*NORBI+J) = D0 1997 WRK(KTMP1-1+(J-1)*NORBI+I) = D0 1998 210 CONTINUE 1999 220 CONTINUE 2000 IF (LNOROT) THEN 2001 IORBI = IORB(ISYM) 2002 DO 320 I = 1,NORBI 2003 IF (NOROT(IORBI+I) .NE. 0) THEN 2004 DO 310 J = 1,NORBI 2005 WRK(KTMP1-1+(I-1)*NORBI+J) = D0 2006 WRK(KTMP1-1+(J-1)*NORBI+I) = D0 2007 310 CONTINUE 2008 END IF 2009 320 CONTINUE 2010 END IF 2011C 2012C note: now TMP1 = 4 Fmo = 4 Ct SAO (DCAO FDAO + DVAO FCAO) C 2013C (SMOAO contains "4 Ct SAO") 2014C note: DGETAP gives 0.5 (TMP1 - TMP1^t) = 2 ( Fmo - Fmo^t) = gradient 2015C 2016 CALL DGETAP(NORBI,WRK(KTMP1),ERRVEC(1+IIORB(ISYM))) 2017C 2018 IF (IPR_DIIS .GE. 15) THEN 2019 WRITE (LUPRI,'(/A/A,I3/A)') 2020 & ' (DIIS_EVC) 4 Fmo = 4 (DCmo FDmo + DVmo FCmo)', 2021 & ' = 4 Ct SAO (DCAO FDAO + DVAO FCAO) C for symmetry', 2022 & ISYM, 2023 & ' (error vector = gradient = 2 (Fmo - Fmo(transposed))' 2024 CALL OUTPUT(WRK(KTMP1),1,NORBI,1,NORBI,NORBI,NORBI,-1,LUPRI) 2025 END IF 2026 800 CONTINUE 2027 IF (HSROHF) THEN 2028 CALL DAXPY(NNBAST,D1,DVAO,1,DCAO,1) 2029 CALL DSCAL(NNBAST,DM1,DVAO,1) 2030 END IF 2031 CALL MEMREL('DIIS_EVC',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 2032 CALL QEXIT('DIIS_EVC') 2033 RETURN 2034 END 2035C /* Deck DIIS_RED */ 2036 SUBROUTINE DIIS_RED (C2DIIS,BMAT,INDEVC,NEVC,DAMP,THREVC, 2037 & CVEC,XLMBDA,ERRVEC,L_ERRVEC, 2038 & WRK,KFRSAV,LFRSAV,IPR_DIIS) 2039C 2040C 19-May-1993 HJAaJ+HA (based on KAPRED) 2041C 2042C Input: 2043C C2DIIS, to select C1-DIIS or C2-DIIS: 2044C FALSE SOLVE LEVEL SHIFTED LINEAR SET OF EQUATIONS 2045C (LEVEL SHIFT DAMP) TO FIND IMPROVED ORBITAL PARAMETERS 2046C TRUE SOLVE LEVEL SHIFTED LINEAR SET OF EQUATIONS 2047C AS AN EIGENVALUE PROBLEM. ADJUST LEVEL SHIFT 2048C TO OBTAIN STEP LENGTH RTRUST 2049C BMAT, the B matrix 2050C ERRVEC, the NEVC error vectors 2051C NEVC, number of error vectors in ERRVEC 2052C THREVC, threshold for acceptable solution in C2DIIS 2053C DAMP, damping factor IN LINEAR SET OF EQUATIONS 2054C 2055C Output: 2056C BMAT, the new, extended reduced orbital Hessian-matrix 2057C CVEC :.not.C2DIIS - SOLUTION TO LINEAR SET OF EQUATIONS 2058C : C2DIIS - LOWEST acceptable EIGENVECTOR 2059C 2060#include "implicit.h" 2061#include "priunit.h" 2062 DIMENSION BMAT(*), CVEC(*), ERRVEC(L_ERRVEC,*), WRK(*) 2063 LOGICAL C2DIIS 2064C 2065 PARAMETER (D0 = 0.0D0, D1 = 1.0D0) 2066C 2067 IROW(I) = (I*(I-1))/2 2068C 2069 CALL QENTER('DIIS_RED') 2070 KFREE = KFRSAV 2071 LFREE = LFRSAV 2072C 2073 IF (IPR_DIIS .GE. 5) WRITE (LUPRI,7555) C2DIIS,NEVC 2074 7555 FORMAT(/' (DIIS_RED) Control Parameters C2DIIS =',L10, 2075 & ' NEVC =',I5) 2076C 2077C 2078C Section 1: extend BMAT with a new row corresponding to 2079C the new error vector 2080C 2081 K = INDEVC 2082 KROW = IROW(K) 2083 DO, L = 1,K 2084 BMAT(KROW+L) = DDOT(L_ERRVEC,ERRVEC(1,L),1,ERRVEC(1,K),1) 2085 END DO 2086 DO, L = K+1,NEVC 2087 BMAT(IROW(L)+K) = DDOT(L_ERRVEC,ERRVEC(1,L),1,ERRVEC(1,K),1) 2088 END DO 2089 IF ( IPR_DIIS.GE.5 ) THEN 2090 WRITE (LUPRI,'(//A)') ' (DIIS_RED) B matrix:' 2091 CALL OUTPAK(BMAT,NEVC,-1,LUPRI) 2092 END IF 2093C 2094 LBMAT = IROW(NEVC+1) 2095 IF (C2DIIS) THEN 2096 LBAUG = LBMAT 2097 ELSE 2098 LBAUG = IROW(NEVC+2) 2099 END IF 2100 CALL MEMGET2('REAL','BTMP',KBTMP,LBAUG,WRK,KFREE,LFREE) 2101 CALL DCOPY(LBMAT,BMAT,1,WRK(KBTMP),1) 2102 IF (.NOT. C2DIIS) GO TO 2000 2103C 2104C ******************************************************* 2105C 2106C C2DIIS: SOLVE DIIS AS AN EIGENVALUE EQUATION 2107C 2108C 2109 CALL MEMGET2('REAL','EVEC',KEVEC,NEVC*NEVC,WRK,KFREE,LFREE) 2110 CALL MEMGET2('REAL','WJ' ,KWJ ,NEVC ,WRK,KFREE,LFREE) 2111 CALL MEMGET2('INTE','IWJ' ,KIWJ ,NEVC ,WRK,KFREE,LFREE) 2112C 2113 CALL DUNIT(WRK(KEVEC),NEVC) 2114 CALL JACO_THR(WRK(KBTMP),WRK(KEVEC),NEVC,NEVC,NEVC,1.0D-12) 2115 DO 150 I = 1,NEVC 2116 II = KBTMP-1+IROW(I+1) 2117 WRK(KBTMP-1+I) = WRK(II) 2118 150 CONTINUE 2119 CALL ORDER (WRK(KEVEC),WRK(KBTMP),NEVC,NEVC) 2120 IOK = 0 2121 DO 170 I = 1,NEVC 2122 EVCSUM = DSUM(NEVC,WRK(KEVEC+(I-1)*NEVC),1) 2123 IF (ABS(EVCSUM) .GT. THREVC) THEN 2124 IOK = I 2125 XLMBDA = WRK(KBTMP-1+I) 2126 CALL DCOPY(NEVC,WRK(KEVEC+(I-1)*NEVC),1,CVEC,1) 2127 CALL DSCAL(NEVC,(D1/EVCSUM),CVEC,1) 2128 GO TO 171 2129 END IF 2130 170 CONTINUE 2131 171 CONTINUE 2132 IF ( IPR_DIIS.GE.3 .OR. IOK.NE.1) THEN 2133 WRITE(LUPRI,'(//A,I5)') 2134 & ' (DIIS_RED) C2-DIIS B matrix eigenvalues. IOK =',IOK 2135 CALL OUTPUT(WRK(KBTMP),1,1,1,NEVC,1,NEVC,-1,LUPRI) 2136 END IF 2137 IF ( IPR_DIIS.GE.4) THEN 2138 WRITE(LUPRI,'(/A)') ' - and B matrix eigenvectors:' 2139 CALL OUTPUT(WRK(KEVEC),1,NEVC,1,NEVC,NEVC,NEVC,-1,LUPRI) 2140 END IF 2141C 2142 GO TO 9999 2143C 2144C ******************************************** 2145C 2146C DIIS: SOLVE LEVEL SHIFTED linear EQUATIONS 2147C 2148 2000 CONTINUE 2149C 2150 NLEQ = NEVC + 1 2151 CALL MEMGET2('INTE','IPVT',KIPVT,NLEQ,WRK,KFREE,LFREE) 2152C 2153 IF (DAMP .NE. D0) THEN 2154 CALL QUIT('DIIS_RED: DAMP not implemented yet') 2155#if defined (VAR_KAPREDCODE) 2156 NOTE: Hamilton and Pulay uses '* (D1 + DAMP)' 2157 and not ' + DAMP' 2158 DO 510 I = 1,NEVC 2159 II = KBTMP - 1 + IROW(I+1) 2160 WRK(II) = WRK(II) + DAMP 2161 510 CONTINUE 2162 IF ( IPR_DIIS.GE.3 ) THEN 2163 WRITE(LUPRI,'(A,1P,D15.8)') 2164 & ' (DIIS_RED) B matrix damping factor =',DAMP 2165 END IF 2166#endif 2167 END IF 2168 DO 520 I = 1,NEVC 2169 WRK(KBTMP-1+LBMAT+I) = D1 2170 520 CONTINUE 2171 WRK(KBTMP-1+LBAUG) = D0 2172 CALL DZERO(CVEC,NEVC) 2173 CVEC(NLEQ) = D1 2174 CALL DSPSOL(NLEQ,1,WRK(KBTMP),CVEC,WRK(KIPVT),INFO) 2175 XLMBDA = CVEC(NLEQ) 2176 IF ( IPR_DIIS.GE.7 ) THEN 2177 WRITE(LUPRI,'(/A)') 2178 & ' (DIIS_RED) Solutions to C1-DIIS set of linear equations' 2179 WRITE (LUPRI,'(10F12.6)') (CVEC(I),I=1,NLEQ) 2180 END IF 2181 IF (INFO.NE.0) THEN 2182 WRITE (LUPRI,8500) INFO 2183 WRITE (LUPRI,'(//A)') ' (DIIS_RED) B matrix:' 2184 CALL OUTPAK(BMAT,NEVC,-1,LUPRI) 2185 CALL QTRACE(LUPRI) 2186 CALL QUIT('DIIS_RED: no solution to linear equations') 2187 END IF 2188 8500 FORMAT(/' (DIIS_RED) Solution not obtained to linear equations' 2189 & /T11,'Check if matrix is singular, LINPACK DSP code =',I3) 2190C 2191C *** End of subroutine DIIS_RED 2192C 2193 9999 CONTINUE 2194 IF (IPR_DIIS .GT. 1) THEN 2195 WRITE (LUPRI,'(/A,L2,A,1P,D10.2)') 2196 & ' DIIS C vector; C2DIIS =',C2DIIS,', LAMBDA =',XLMBDA 2197 WRITE (LUPRI,'(10F12.6)') (CVEC(I),I=1,NEVC) 2198 END IF 2199 CALL MEMREL('DIIS_RED',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 2200 CALL QEXIT('DIIS_RED') 2201 RETURN 2202 END 2203c /* DFT_ADD_KS */ 2204 SUBROUTINE DFT_ADD_KS(ONLYFC,EMY,DCAO,NDMAT,FCAO,FVAO, 2205 & WRK,LWRK,IPRFCK) 2206c TMP1 is of NNBASX size, TMP2 is of N2BASX size. They are passed 2207c because they are available but they could as well be allocated 2208c here. 2209#include "implicit.h" 2210c 2211#include "inforb.h" 2212#include "dfterg.h" 2213#include "priunit.h" 2214c 2215 PARAMETER (DP5 = 0.5D0, D1 = 1.0D0) 2216c 2217 DIMENSION DCAO(N2BASX*NDMAT),FCAO(N2BASX),FVAO(N2BASX) 2218 DIMENSION WRK(LWRK) 2219 LOGICAL ONLYFC 2220c 2221 CALL QENTER('DFT_ADD_KS') 2222c 2223 KTMP1 = 1 2224 KTMP2 = KTMP1 + NNBASX 2225 KLST = KTMP2 + N2BASX 2226 LFREE = LWRK - KLST +1 2227 IF(LFREE.LT.0) CALL STOPIT('ADDKS1','SIRDIIS',KLST,LFREE) 2228C EXCTRO = DP5*DDOT(NNBASX,DCAO,1,FCAO,1) 2229 IF(ONLYFC) THEN 2230 CALL DZERO(WRK(KTMP2),N2BASX) 2231 CALL DFTKSMb(DCAO,WRK(KTMP2),EDFTY,WRK(KLST),LFREE,IPRFCK) 2232 CALL DGETSP(NBAST,WRK(KTMP2),WRK(KTMP1)) 2233 IF (NSYM .GT. 1) THEN 2234 KKOHN = KLST 2235 KLST = KKOHN + NNBAST 2236 IF (KLST.GT.LWRK) CALL STOPIT('ADDKS','SIRDIIS',KLST,LWRK) 2237 CALL PKSYM1(WRK(KTMP1),WRK(KKOHN),NBAS,NSYM,1) 2238 CALL DAXPY(NNBAST,1D0,WRK(KKOHN),1,FCAO,1) 2239 ELSE 2240 CALL DAXPY(NNBAST,1D0,WRK(KTMP1),1,FCAO,1) 2241 END IF 2242 ELSE 2243C memory allocation 2244 KDENA = KLST 2245 KKSMA = KDENA +N2BASX*2 2246 KKSPA = KKSMA +N2BASX*2 2247 KKSPB = KKSPA +NNBAST 2248 KFREE = KKSPB +NNBAST 2249 IF (KFREE.GT.LFREE) CALL STOPIT('ADDKS1','SIRDIIS',KFREE,LFREE) 2250 CALL DZERO(WRK(KDENA),N2BASX*2) 2251 CALL DZERO(WRK(KKSMA),N2BASX*2) 2252 CALL DZERO(WRK(KKSPA),NNBAST) 2253 CALL DZERO(WRK(KKSPB),NNBAST) 2254C Kohn-Sham contribution evaluation 2255 IF (NISHT .GT. 0) THEN 2256 CALL DAXPY(N2BASX,DP5,DCAO,1,WRK(KDENA),1) 2257 CALL DAXPY(N2BASX,DP5,DCAO,1,WRK(KDENA+N2BASX),1) 2258 END IF 2259 CALL DAXPY(N2BASX,1.D0,DCAO(1+N2BASX),1,WRK(KDENA),1) 2260C 2261C For specific case of molecules bearing only alpha electron(s), 2262C like hydrogen atom, we call "old" open-shell code, which 2263C can handle zero beta density. 2264C 2265 IF (NISHT .EQ. 0) THEN 2266 CALL dft_kohn_shamab(WRK(KDENA),WRK(KKSMA),EDFTY, 2267 & WRK(KFREE),LFREE,IPRFCK) 2268 ELSE 2269 CALL dft_kohn_shamab_b(WRK(KDENA),WRK(KKSMA),EDFTY, 2270 & WRK(KFREE),LFREE,IPRFCK) 2271 END IF 2272C packing 2273 CALL DGETSP(NBAST,WRK(KKSMA),WRK(KTMP1)) 2274 IF (NSYM .GT. 1) THEN 2275 CALL PKSYM1(WRK(KTMP1),WRK(KKSPA),NBAS,NSYM,1) 2276 ELSE 2277 CALL DAXPY(NNBAST,D1,WRK(KTMP1),1,WRK(KKSPA),1) 2278 END IF 2279 CALL DGETSP(NBAST,WRK(KKSMA+N2BASX),WRK(KTMP1)) 2280 IF (NSYM .GT. 1) THEN 2281 CALL PKSYM1(WRK(KTMP1),WRK(KKSPB),NBAS,NSYM,1) 2282 ELSE 2283 CALL DAXPY(NNBAST,D1,WRK(KTMP1),1,WRK(KKSPB),1) 2284 END IF 2285c WRITE(LUPRI,'(/A)') 'ksma: ' 2286c CALL OUTPKB(WRK(KKSPA),NORB,NSYM,-1,LUPRI) 2287c WRITE(LUPRI,'(/A)') 'ksmB: ' 2288c CALL OUTPKB(WRK(KKSPB),NORB,NSYM,-1,LUPRI) 2289c 2290 CALL DAXPY(NNBAST, DP5,WRK(KKSPA),1,FCAO,1) 2291 CALL DAXPY(NNBAST, DP5,WRK(KKSPB),1,FCAO,1) 2292c 2293 CALL DAXPY(NNBAST,-DP5,WRK(KKSPA),1,FVAO,1) 2294 CALL DAXPY(NNBAST, DP5,WRK(KKSPB),1,FVAO,1) 2295 END IF 2296C 2297 EMY = EMY + EDFTY 2298C EDFTY = EDFTY +EXCTRO-DP5*DDOT(N2BASX,DCAO,1,FCAO,1) 2299 CALL QEXIT('DFT_ADD_KS') 2300 RETURN 2301 END 2302C /* Deck fckeig */ 2303 SUBROUTINE FCKEIG(CMO,FC,SHFTLVL,THR_FCKEIG,NCAN,IPRINT_in, 2304 & SCRA,LSCRA) 2305C 2306C Written 20-May-1993 by Hans Jorgen Aa. Jensen 2307C 2308C Purpose: 2309C If requested (by SHFTLVL .ne. 0) : 2310C Perform "restricted step" step reduction (i.e. level shift) 2311C by adding SHFTLVL*DENS(i,i) to all diagonal elements; 2312C Diagonalize Fock matrix 2313C 2314C Input: 2315C CMO; initial molecular orbitals used to build Fock matrix, 2316C assumed to be orthonormal. 2317C FC; the inactive Fock matrix 2318C NCAN(8); number of orbitals in each symmetry to make canonical 2319C 2320C Output: 2321C CMO; molecular orbitals diagonalizing Fock matrix 2322C FC; the orbital energies 2323C 2324C Scratch: 2325C SCRA; general scratch area 2326C 2327#include "implicit.h" 2328 DIMENSION CMO(*),FC(*),NCAN(8),SCRA(*) 2329C 2330C Used from common blocks: 2331C INFINP : SUPSYM, LNOROT, NOROT() 2332C SCBRHF : NFRRHF(*), AUTOCC 2333C INFIND : IROW(*),...,ISSMO(*),? 2334C 2335#include "maxash.h" 2336#include "maxorb.h" 2337#include "priunit.h" 2338#include "infinp.h" 2339#include "inforb.h" 2340#include "scbrhf.h" 2341#include "infind.h" 2342#include "infpri.h" 2343#include "mxcent.h" 2344#include "dftcom.h" 2345#include "dftacb.h" 2346#include "dummy.h" 2347C 2348 LOGICAL LSAVE4,LSAVE6 2349 PARAMETER (D0 = 0.0D0, DP5 = 0.5D0) 2350C 2351 CALL QENTER('FCKEIG') 2352 call gettim(t1,w1) 2353C 2354 IPR_FCKEIG = IPRINT_in 2355 LSAVE4 = P4FLAG(9) 2356 LSAVE6 = P6FLAG(6) 2357 P4FLAG(9) = .FALSE. 2358 P6FLAG(6) = .FALSE. 2359C 2360C Some memory allocation for SCF occupation determination 2361C 2362 KEIG = 1 2363 KSYMS = KEIG + NORBT 2364 KLAST = KSYMS + NORBT 2365C 2366C Step 1: Diagonalize (level-shifted) Fock-matrix: 2367C 2368 K = 0 2369 DO, ISYM = 1,NSYM 2370 K = K + NOCC(ISYM)*(NORB(ISYM)-NOCC(ISYM)) 2371 END DO 2372 FAC_THR_JACO = 16*K 2373 ! factor 16 because gradient ia = 2 F_ia = 4 FC_ia 2374 FAC_THR_JACO = 1.0D0 / SQRT(FAC_THR_JACO) 2375 ! factor to ensure gradient is less than THR_FCKEIG if all 2376 ! occ-vir elements are less than THR_FCKEIG/FAC_THR_JACO /hjaaj 2377 2378 DO 200 ISYM = 1,NSYM 2379 NCANI = NCAN(ISYM) 2380 IF ( .NOT.AUTOCC .AND. NCANI.EQ.0 ) GO TO 200 2381 NORBI = NORB(ISYM) 2382 IF (NORBI.EQ.0) GO TO 200 2383 IORBI = IORB(ISYM) 2384 NBASI = NBAS(ISYM) 2385 NFRZI = NFRRHF(ISYM) 2386 IFSYM = IIORB(ISYM) 2387 ICSYM = ICMO(ISYM) + 1 2388C 2389C If requested (by SHFTLVL .ne. 0) : 2390C Perform "restricted step" step reduction (i.e. level shift) 2391C by adding 0.5*SHFTLVL*DENS(i,i) to all diagonal elements: 2392C /hjaaj March 2001 2393C 2394 IF (SHFTLVL .NE. D0) THEN 2395 DO K = 1,NISH(ISYM) 2396 FC(IFSYM+IROW(K+1)) = FC(IFSYM+IROW(K+1)) + SHFTLVL 2397 END DO 2398 DO K = NISH(ISYM)+1, NISH(ISYM)+NASH(ISYM) 2399 FC(IFSYM+IROW(K+1)) = FC(IFSYM+IROW(K+1)) + DP5*SHFTLVL 2400 END DO 2401 END IF 2402C zero rows and columns corresponding to frozen orbitals 2403 JFRZI = 0 2404 IF (NFRZI .GT. 0 .OR. LNOROT) THEN 2405 JNFRZ = 0 2406 DO 130 K = 1,NORBI 2407 IK=IORBI+K 2408 IF (K .LE. NFRZI .OR. NOROT(IK) .NE. 0) THEN 2409 IF (IOBTYP(IK).NE.JTACT) THEN 2410C ... enable freezing an active orbital 2411 IF (JNFRZ .GT. 0) GO TO 9000 2412C ... exit if free orbital below this frozen orbital 2413C because then ORDRSS cannot be called 2414 JFRZI = JFRZI + 1 2415 END IF 2416 KROW = IROW(K) 2417 DO 110 L = 1,K-1 2418 FC(IFSYM+KROW+L) = D0 2419 110 CONTINUE 2420 DO 120 L = K+1,NORBI 2421 KL = (IFSYM+K)+IROW(L) 2422 FC(KL) = D0 2423 120 CONTINUE 2424 IF (NOROT(IK).NE.0 .AND. IOBTYP(IK).EQ.JTACT) THEN 2425 KK = (IFSYM+K)+IROW(K) 2426 FC(KK) = DUMMY 2427C 2428C This assignment puts it last in the sorting - put back in place later 2429C 2430 END IF 2431 ELSE 2432 JNFRZ = JNFRZ + 1 2433 END IF 2434 130 CONTINUE 2435 END IF 2436C 2437 IF (NCANI .gt. 0) THEN 2438 THR_JACO = THR_FCKEIG * FAC_THR_JACO 2439 THR_JACO = MIN( THR_JACO, 1.D-3) 2440 ELSE IF (AUTOCC) THEN 2441 THR_JACO = 1.D-3 2442 ELSE 2443 GO TO 200 2444 END IF 2445 2446 NMAX = MIN(NORBI, NCANI + 5) 2447 CALL JACO_THR(FC(IFSYM+1),CMO(ICSYM),NORBI,NMAX,NBASI,THR_JACO) 2448C CALL JACO_THR(F,VEC,NB,NMAX,NROWV,THR_CONV) 2449C 2450 IF (IPR_FCKEIG .GT. 20) THEN 2451 write (lupri,*) 'FCKEIG debug: FC after JACO_THR',NMAX,NORBI 2452 call outpak(FC(IFSYM+1),NORBI,-1,LUPRI) 2453 END IF 2454 DO, I=1,NORBI 2455 FC(IORBI+I) = FC(IFSYM+IROW(I+1)) 2456 END DO 2457C 2458 IF (AUTOCC) THEN 2459 CALL DCOPY(NORBI,FC(IORBI + 1),1,SCRA(KEIG + IORBI),1) 2460 DO, IK = 1, NORBI 2461 SCRA(KSYMS + IORBI + IK - 1) = ISYM 2462 END DO 2463 END IF 2464C 2465C 2466 JCSYM = ICSYM + JFRZI*NBASI 2467 JORBI = IORBI + JFRZI + 1 2468 NNOTFR= NORBI - JFRZI 2469 CALL ORDRSS(CMO(JCSYM),FC(JORBI),ISSMO(JORBI),NNOTFR,NBASI) 2470C 2471C CALL HEADER('Sorted orbitals',-1) 2472C CALL PRORB(CMO,.FALSE.,LUPRI) 2473C 2474C If there are frozen open shell orbitals they are last - put back in place 2475C 2476 IF (LNOROT) THEN 2477 NISHI=NISH(ISYM) 2478 NASHI=NASH(ISYM) 2479 JA=NISHI+1 2480 KA=IORBI+NISHI+1 2481 DO IA=1,NASHI 2482 JA=IA+NISHI 2483 KA=JA+IORBI 2484 IF (NOROT(KA).NE.0) THEN 2485 CALL ROTCOL(CMO(ICSYM),NBASI,NORBI,JA,NORBI) 2486 CALL ROTCOL(FC(IORBI+1),1,NORBI,JA,NORBI) 2487 END IF 2488 END DO 2489 END IF 2490C 2491C CALL HEADER('Reordered orbitals',-1) 2492C CALL PRORB(CMO,.FALSE.,LUPRI) 2493 IF (IPR_FCKEIG .GE. 4) WRITE (LUPRI,'(/A,I2/,(5F15.5))') 2494 & ' DIIS Fock eigenvalues in symmetry',ISYM, 2495 & (FC(IORBI+I), I = 1,NORBI) 2496 2497 2498 200 CONTINUE 2499 IF (SUPSYM) CALL AVEORD() 2500C ... remake ISSORD() as ISSMO() may have changed in ORDRSS 2501C 2502C Step 2: Reorthogonalize new mo's 2503C 2504 KSAO = KLAST 2505 KSCR1 = KSAO + NNBAST 2506 LSCR1 = LSCRA - KSCR1 2507 2508 IF (DFTASC) CALL GET_HOMO_ASC(FC,SCRA(KSCR1),LSCR1) 2509 2510 CALL ORTHO(CMO,SCRA(KSAO),SCRA(KSCR1),LSCR1) 2511 IF (SUPSYM) THEN 2512 KFREE = 1 2513 LFREE = LSCRA 2514 CALL AVECPH(IPHCHA,CMO,SCRA(KLAST),KFREE,LFREE) 2515 END IF 2516C 2517C Step 3: Reorder Hartree-Fock occupation just in case last 2518C occupation suggestion was wrong 2519C 2520 IF (AUTOCC) THEN 2521 CALL ORDER(SCRA(KSYMS),SCRA(KEIG),NORBT,1) 2522 CALL IZERO(NISH,8) 2523 MOCC = NRHFEL/2 2524 DO 98 IK = 1, MOCC 2525 ISYM = NINT(SCRA(KSYMS + IK -1)) 2526 NISH(ISYM) = NISH(ISYM) + 1 2527 98 CONTINUE 2528 CALL ICOPY(8,NISH,1,NOCC,1) 2529 IF (2*MOCC .NE. NRHFEL) THEN 2530 IOPRHF = NINT(SCRA(KSYMS + MOCC)) 2531 LSYM = IOPRHF 2532 CALL IZERO(NASH,8) 2533 CALL IZERO(IASH,8) 2534 NASH(IOPRHF) = 1 2535 NOCC(IOPRHF) = NOCC(IOPRHF) + 1 2536 DO 96 ISYM = IOPRHF + 1, 8 2537 IASH(ISYM) = 1 2538 96 CONTINUE 2539 END IF 2540 END IF 2541C 2542C *** end of subroutine FCKEIG 2543C 2544 IF (IPR_FCKEIG > 0) THEN 2545 call gettim(t2,w2) 2546 write(lupri,'(A,1P,3D10.2)') 2547 & ' Time used in FCKEIG, thr_jaco: ',t2-t1,w2-w1,thr_jaco 2548 END IF 2549 CALL QEXIT('FCKEIG') 2550 RETURN 2551C 2552 9000 CONTINUE 2553 WRITE (LUPRI,'(//A/A)') 2554 & ' ERROR in FCKEIG, FCKEIG cannot handle ".FREEZE" for orbitals', 2555 & ' unless they also could have been frozen with ".FROZEN"' 2556 CALL QTRACE(LUPRI) 2557 CALL QUIT('FCKEIG error, cannot handle orbitals '// 2558 & 'frozen with ".FREEZE"') 2559 2560 CONTAINS 2561 2562 SUBROUTINE ROTCOL(A,NR,NC,I1,I2) 2563C 2564C Cyclic permutation of columns in A 2565C I2>I1 is moved to I1 and rest shifted right 2566C 2567 IMPLICIT NONE 2568 DOUBLE PRECISION A,B 2569 INTEGER NR,NC,I1,I2 2570 DIMENSION A(NR,NC),B(NR) 2571 2572 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: TMP 2573 INTEGER I,MC 2574 IF (I2.LE.I1) CALL QEXIT('ROTCOL ERROR: I1<I2 not fulfilled') 2575 IF (I1.LE.0 .OR. I2.LE.0 .OR. I1.GT.NC.OR.I2.GT.NC) 2576 & CALL QEXIT('ROTCOL ERROR: I1 or I2 out of bounds') 2577 MC=(I2-I1+1) 2578 ALLOCATE(TMP(NR,MC)) 2579 TMP(:,1)=A(:,I2) 2580 TMP(:,2:)=A(:,I1:I2-1) 2581 A(:,I1:I2)=TMP(:,:) 2582 DEALLOCATE(TMP) 2583 RETURN 2584 END SUBROUTINE ROTCOL 2585 2586 END ! subroutine FCKEIG 2587 2588 SUBROUTINE FVCORR(NVEC,W,DAO,FAO,MXERRV,F,WRK,KFREE,LFREE) 2589 IMPLICIT NONE 2590C 2591C Form projected (effective) fock matrix in AO basis for each accumulated 2592C iteration and sum up with diis weights. 2593C 2594C Input: NVEC accumulated Fock (FAO) and density (DAO) 2595C matrices from previous DIIS iterations 2596C 2597C Output: sum of weighted effective fock matrices to be diagonalized 2598C 2599C 2600 INTEGER NVEC, MXERRV, KFREE, LFREE 2601 REAL*8 W(NVEC) 2602#include "inforb.h" 2603#include "priunit.h" 2604C REAL*8 DAO(NNBAST,MXERRV,2), FAO(NNBAST,MXERRV,2), F(NNBAST) 2605 REAL*8 DAO(NNBAST,1,2), FAO(NNBAST,MXERRV,2), F(NNBAST) 2606 REAL*8 WRK(*) 2607C 2608C Local 2609C 2610 REAL*8 D0, D1, D2 2611 PARAMETER (D0=0.0D0, D1=1.0D0, D2=2.0D0) 2612 INTEGER ISYM, KS, KTMP1, KTMP2, KDI, KFI, KSI 2613 INTEGER NBASI, NNBASI, IIBASI, N2BASI, IBLK, IVEC, I, II 2614 INTEGER IMXVEC 2615 EXTERNAL IMXVEC 2616 LOGICAL FOUND 2617C 2618 CALL QENTER('FVCORR') 2619 2620 CALL MEMGET2('REAL','S',KS,NNBAST,WRK,KFREE,LFREE) 2621 N2BASI=IMXVEC(N2BAS,NSYM) 2622 CALL MEMGET2('REAL','DI',KDI,N2BASI,WRK,KFREE,LFREE) 2623 CALL MEMGET2('REAL','SI',KSI,N2BASI,WRK,KFREE,LFREE) 2624 CALL MEMGET2('REAL','FI',KFI,N2BASI,WRK,KFREE,LFREE) 2625 CALL MEMGET2('REAL','TMP1',KTMP1,N2BASI,WRK,KFREE,LFREE) 2626 CALL MEMGET2('REAL','TMP2',KTMP2,N2BASI,WRK,KFREE,LFREE) 2627 FOUND = .FALSE. 2628 CALL RDONEL('OVERLAP',FOUND,WRK(KS),NNBAST) 2629 IF (.NOT.FOUND) THEN 2630 CALL QUIT('FVCORR:Error reading overlap') 2631 END IF 2632 DO ISYM=1,NSYM 2633 NBASI=NBAS(ISYM) 2634 NNBASI=NNBAS(ISYM) 2635 N2BASI=N2BAS(ISYM) 2636 IIBASI=IIBAS(ISYM) 2637 IBLK=IIBAS(ISYM)+1 2638 IF (NBASI.EQ.0) GO TO 100 2639 CALL DSPTSI(NBASI,WRK(KS+IIBASI),WRK(KSI)) 2640 DO IVEC=1,NVEC 2641 CALL DUNFLD(NBASI,DAO(IBLK,IVEC,2),WRK(KDI)) 2642 ! -s*(-do) 2643 CALL DGEMM( 2644 & 'N','N',NBASI,NBASI,NBASI, 2645 & -D1,WRK(KSI),NBASI, 2646 & WRK(KDI),NBASI, 2647 & D0,WRK(KTMP1),NBASI 2648 & ) 2649 ! (s*do)*(fc-fo) 2650 CALL DSPTSI(NBASI,FAO(IBLK,IVEC,2),WRK(KFI)) 2651 CALL DGEMM( 2652 & 'N','N',NBASI,NBASI,NBASI, 2653 & D1,WRK(KTMP1),NBASI, 2654 & WRK(KFI),NBASI, 2655 & D0,WRK(KTMP2),NBASI 2656 & ) 2657 ! (do+dc)*s 2658 !CALL DUNFLD(NBASI,DAO(IBLK,IVEC,1),WRK(KTMP1)) 2659 !CALL DAXPY(N2BASI,D1,WRK(KTMP1),1,WRK(KDI),1) 2660 CALL DUNFLD(NBASI,DAO(IBLK,IVEC,1),WRK(KDI)) 2661 CALL DGEMM( 2662 & 'N','N',NBASI,NBASI,NBASI, 2663 & D1,WRK(KDI),NBASI, 2664 & WRK(KSI),NBASI, 2665 & D0,WRK(KTMP1),NBASI 2666 & ) 2667 ! (do+dc)*s-1 2668 DO I=1,NBASI 2669 II=KTMP1+(NBASI+1)*(I-1) 2670 WRK(II) = WRK(II) - D1 2671 END DO 2672 !final correction 2673 ! (s*do)*(fc-fo)*((do+dc)*s-1) 2674 CALL DGEMM( 2675 & 'N','N',NBASI,NBASI,NBASI, 2676 & D1,WRK(KTMP2),NBASI, 2677 & WRK(KTMP1),NBASI, 2678 & D0,WRK(KDI),NBASI 2679 & ) 2680 CALL DGETSP(NBASI,WRK(KDI),WRK(KFI)) 2681 CALL DAXPY(NNBASI,D2*W(IVEC),WRK(KFI),1,F(IBLK),1) 2682 END DO 2683 100 CONTINUE 2684 END DO 2685 CALL MEMREL('FVCORR',WRK,KS,KS,KFREE,LFREE) 2686 CALL QEXIT('FVCORR') 2687 END 2688 2689 SUBROUTINE GET_HOMO_ASC(FD,SCRA,LSCRA) 2690 2691#include "implicit.h" 2692#include "priunit.h" 2693 DIMENSION FD(*),SCRA(*) 2694 PARAMETER (DBIG = 1.D+12) 2695#include "maxash.h" 2696#include "maxorb.h" 2697#include "mxcent.h" 2698#include "inforb.h" 2699#include "infinp.h" 2700#include "infind.h" 2701#include "dftacb.h" 2702#include "dftcom.h" 2703 2704 logical exchanged 2705C 2706 KORBEN = 1 2707 LNEED = KORBEN + NORBT 2708 IF (LNEED .GT. LSCRA) CALL ERRWRK('GET_HOMO_ASC',LNEED,LSCRA) 2709 CALL DZERO(SCRA,NORBT) 2710 2711 call dcopy(NORBT,FD,1,SCRA,1) 2712100 continue 2713 exchanged=.false. 2714 do i=1,norbt-1 2715 if (scra(i+1) .lt. scra(i) ) then 2716 tmp=scra(i+1) 2717 scra(i+1)=scra(i) 2718 scra(i)=tmp 2719 exchanged=.true. 2720 end if 2721 end do 2722 if (exchanged) goto 100 2723 2724 EhomoA = scra(NOCCT) 2725 ehomo = ehomoa 2726 ehomob = scra(NishT) 2727 2728 ElumoA = scra(NOCCT+1) 2729 ElumoB = scra(NishT+1) 2730 2731 SHIFA = DFTIPTA + EHOMOA - SHF 2732 2733 write(LUPRI,10)' E_HOMO Alpha ',EhomoA,' E_HOMO Beta ', 2734 & EhomoB, 'v_xc(inf)', SHIFA 2735 write(LUPRI,10)' E_LUMO Alpha ',ElumoA,' E_LUMO Beta ', 2736 & ElumoB 273710 format(5x,A14,F20.12,A14,F20.12,A14,F20.12) 2738 2739 RETURN 2740 END 2741 2742 SUBROUTINE SDFCOMM(F,D,SAOUNP,COMM,KFRSAV,LFRSAV,WRK) 2743C In: F apk,D apk and S unp 2744C Out: SDF-FDS 2745 2746#include "implicit.h" 2747#include "inforb.h" 2748 DOUBLE PRECISION F(NNBAST),D(NNBAST),SAOUNP(NBAST,NBAST) 2749 DOUBLE PRECISION COMM(NBAST,NBAST),WRK(*) 2750 PARAMETER (D1=1.0D0, D0=0.0D0, DM1=-1.0D0) 2751 2752 KFREE = KFRSAV 2753 LFREE = LFRSAV 2754 CALL MEMGET('REAL',KFTMP2,N2BASX,WRK,KFREE,LFREE) 2755 CALL MEMGET('REAL',KDTMP2,N2BASX,WRK,KFREE,LFREE) 2756 CALL MEMGET('REAL',KDForFD,N2BASX,WRK,KFREE,LFREE) 2757 2758 IF (NSYM.GT.1) THEN 2759 CALL MEMGET('REAL',KFTMP,NNBASX,WRK,KFREE,LFREE) 2760 CALL MEMGET('REAL',KDTMP,NNBASX,WRK,KFREE,LFREE) 2761 ENDIF 2762 2763 IF (NSYM.GT.1) THEN 2764 CALL PKSYM1(WRK(KFTMP),F,NBAS,NSYM,-1) 2765 CALL DSPTSI(NBAST,WRK(KFTMP),WRK(KFTMP2)) 2766 CALL PKSYM1(WRK(KDTMP),D,NBAS,NSYM,-1) 2767 CALL DUNFLD(NBAST,WRK(KDTMP),WRK(KDTMP2)) 2768 ELSE 2769 CALL DSPTSI(NBAST,F,WRK(KFTMP2)) 2770 CALL DUNFLD(NBAST,D,WRK(KDTMP2)) 2771 ENDIF 2772 2773 CALL DGEMM('N','N',NBAST,NBAST,NBAST,D1,WRK(KFTMP2),NBAST, 2774 & WRK(KDTMP2),NBAST,D0,WRK(KDForFD),NBAST) 2775 CALL DGEMM('N','N',NBAST,NBAST,NBAST,D1,WRK(KDForFD),NBAST, 2776 & SAOUNP,NBAST,D0,COMM,NBAST) 2777 CALL DGEMM('N','N',NBAST,NBAST,NBAST,D1,WRK(KDTMP2),NBAST, 2778 & WRK(KFTMP2),NBAST,D0,WRK(KDForFD),NBAST) 2779 CALL DGEMM('N','N',NBAST,NBAST,NBAST,D1,SAOUNP,NBAST, 2780 & WRK(KDForFD),NBAST,DM1,COMM,NBAST) 2781 2782 CALL MEMREL('SDFCOMM',WRK,KFRSAV,KFTMP2,KFREE,LFREE) 2783 END 2784! -- end of sirdiis.F -- 2785