1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19c /* deck cc_cphf */ 20*=====================================================================* 21 SUBROUTINE CC_CPHF(TYPE,LABEL,ISYMS,ISTAT,EIGV, 22 & ISYMO,FREQS,ICAU,NVEC,MAXVEC, 23 & WORK,LWORK) 24*---------------------------------------------------------------------* 25* 26* Purpose: solve CPHF equations needed in CC program 27* 28* implemented types: R1 29* 30* Written by Christof Haettig, november 1998 31* 32*=====================================================================* 33#if defined (IMPLICIT_NONE) 34 IMPLICIT NONE 35#else 36#include "implicit.h" 37#endif 38#include "priunit.h" 39#include "dummy.h" 40#include "exeinf.h" 41#include "maxorb.h" 42#include "infvar.h" 43#include "inftap.h" 44#include "iratdef.h" 45#include "ccsdinp.h" 46#include "ccsdsym.h" 47#include "ccorb.h" 48#include "cclr.h" 49#include "ccfro.h" 50#include "inflin.h" 51Cholesky 52#include "ccdeco.h" 53Cholesky 54 55* local parameters: 56 CHARACTER*(18) MSGDBG 57 PARAMETER (MSGDBG = '[debug] CC_CPHF> ') 58 LOGICAL LOCDBG 59 PARAMETER (LOCDBG = .true. ) 60 CHARACTER*8 FNGDVE, FNSOVE, FNREVE, FILCPHF 61 PARAMETER(FNGDVE='CCPHFRHS',FNSOVE='CCPHFSOL') 62 PARAMETER(FNREVE='CCPHFRES') 63 INTEGER MXFRVAL 64 PARAMETER (MXFRVAL = 100) 65C PARAMETER (MXFRVAL = 1) ! switch off simultaneous equations 66 67 CHARACTER TYPE*(3) 68 69 INTEGER NVEC, MAXVEC, LWORK 70 INTEGER ISYMS(MAXVEC,*), ISYMO(MAXVEC,*) 71 INTEGER ISTAT(MAXVEC,*), ICAU(MAXVEC,*) 72 73 CHARACTER*8 LABEL(MAXVEC,*) 74 75 REAL*8 FREQS(MAXVEC,*), EIGV(MAXVEC,*) 76 REAL*8 WORK(LWORK) 77 REAL*8 ZERO, RDUM 78 REAL*8 XNORM, DDOT, DSQRT, GDNORM 79 80 PARAMETER (ZERO = 0.0d0) 81 82 CHARACTER MODEL*(10), MODWR*(10) 83 84 LOGICAL MAYBE_MORE, RELAX, CICLC, HFCLC, TRPCLC,OOTV 85 LOGICAL EXCLC, NEWCMO_SAVE, OPTST, EX 86 INTEGER IOPT, ISYM, IVEC, NSTAT, ORDER, IDUM, IPERT,LWRK0,LWRK1 87 INTEGER NCMOT,NASHT,N2ASHX,LCINDX,KCMO,KUDV,KXINDX,KEND0,KEND1 88 INTEGER LUGDVE, LUSOVE, MXRM, MXPHP, NCOSAV, IOPTWR, KSLVEC 89 INTEGER IREAL, NFRVAL, KFRVAL, NABAOP, NABATY, IFRVAL, IDX 90 INTEGER IPRABA, LUREVE, LUCPHF, MALLAI, INUM 91 92* external functions: 93 INTEGER IR1TAMP 94 95*---------------------------------------------------------------------* 96* do some checks: 97*---------------------------------------------------------------------* 98 CALL QENTER('CC_CPHF') 99 100 IF (NVEC.LT.1) THEN 101 CALL QEXIT('CC_CPHF') 102 RETURN 103 END IF 104 105 IF (LOCDBG) THEN 106 WRITE(LUPRI,*) 'Entered CC_CPHF. NVEC =',NVEC 107 WRITE(LUPRI,*) 'TYPE ',TYPE 108 CALL FLSHFO(LUPRI) 109 END IF 110 111* check vector type: 112 IF (TYPE(1:3).EQ.'R1 ') THEN 113 ORDER = 1 114 NSTAT = 0 115 MODWR = 'SCF? ' 116 ELSE 117 WRITE (LUPRI,*) 'CPHF vectors ',TYPE(1:3),' not implemented.' 118 CALL QUIT('required CPHF vectors not implemented.') 119 END IF 120 121* Refuse to do anything for Cholesky. 122 IF (CHOINT) THEN 123 NWARN = NWARN + 1 124 WRITE(LUPRI,'(/A/A/A/A/)') '*** WARNING from CC_CPHF:', 125 & ' WARNING: Refusing to do CPHF for Cholesky, because', 126 & ' WARNING: ABACUS has not been modified yet!', 127 & ' WARNING: Program continues nevertheless...' 128 RETURN 129 ENDIF 130 131* get some variables from SIRIUS common blocks 132 CALL CC_SIRINF(NCMOT,NASHT,N2ASHX,LCINDX) 133 134* check number of active shells from SIRIUS: 135 IF (NASHT.NE.0) THEN 136 WRITE (LUPRI,*) 'non-zero number of active shells:',NASHT 137 CALL QUIT('Non-zero number of active shells in CC_CPHF.') 138 END IF 139 140*---------------------------------------------------------------------* 141* print header for rhs vector section 142*---------------------------------------------------------------------* 143 WRITE (LUPRI,'(7(/1X,2A),/)') 144 & '------------------------------------', 145 & '-------------------------------', 146 & '| OUTPUT FROM CPHF', 147 & ' SECTION |', 148 & '------------------------------------', 149 & '-------------------------------' 150 151* print some debug/info output 152 IF (IPRINT .GT. 10 .OR. LOCDBG) THEN 153 WRITE(LUPRI,*) 'CC_CPHF Workspace:',LWORK 154 WRITE(LUPRI,*) '1 MODWR ',MODWR 155 END IF 156 157 CALL FLSHFO(LUPRI) 158 159*---------------------------------------------------------------------* 160* some initilizations: 161*---------------------------------------------------------------------* 162 163* maximum of nallai over isym (used as fixed record lengths for LUCPHF: 164 MALLAI = NALLAI(1) 165 DO ISYM = 2, NSYM 166 MALLAI = MAX(MALLAI,NALLAI(ISYM)) 167 END DO 168 169* arrays for GETGPV and ABARSP: 170 KCMO = 1 171 KUDV = KCMO + NCMOT 172 KXINDX = KUDV + N2ASHX 173 KFRVAL = KXINDX + LCINDX 174 KEND0 = KFRVAL + MXFRVAL 175 LWRK0 = LWORK - KEND0 176 177 IF (LWRK0 .LT. 0) THEN 178 CALL QUIT('Insufficient memory in CC_CPHF.') 179 END IF 180 181* read MO coefficients from file: 182 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 183 & .FALSE.) 184 REWIND LUSIFC 185 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 186 READ (LUSIFC) 187 READ (LUSIFC) 188 CALL READI(LUSIFC,IRAT*NCMOT,WORK(KCMO)) 189 CALL GPCLOSE(LUSIFC,'KEEP') 190 191* flags for ABARSP: 192 CICLC = .FALSE. 193 HFCLC = .TRUE. 194 TRPCLC = .FALSE. 195 OOTV = .FALSE. 196 EXCLC = .FALSE. 197 MXRM = 40 198 MXPHP = 1 199 200 NEWCMO_SAVE = NEWCMO 201 NCOSAV = NCONF 202 203 NEWCMO = .TRUE. 204 NCONF = 1 205 206 IF (DIRECT) CALL CCDFFOP 207 208* flags for CC_WRRSP routine: 209 IOPTWR = 4 210 211* open file for right hand side and solution vectors: 212 213 LUGDVE = -1 214 LUSOVE = -1 215 LUREVE = -1 216 CALL GPOPEN(LUGDVE,FNGDVE,'UNKNOWN',' ','UNFORMATTED',IDUMMY, 217 & .FALSE.) 218 CALL GPOPEN(LUSOVE,FNSOVE,'UNKNOWN',' ','UNFORMATTED',IDUMMY, 219 & .FALSE.) 220 CALL GPOPEN(LUREVE,FNREVE,'UNKNOWN',' ','UNFORMATTED',IDUMMY, 221 & .FALSE.) 222* open property integral file: 223 IF (LUPROP .LE. 0) CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN',' ', 224 & 'UNFORMATTED',IDUMMY,.FALSE.) 225 226* close twoelectron integral file 'AOTWOINT': 227 IF (LUINTA.GT.0) CALL GPCLOSE(LUINTA,'KEEP') 228 229* file for CPHF vectors: 230 WRITE(FILCPHF,'(A5,A3)') 'CPHF_',TYPE(1:3) 231 DO I = 6,8 232 IF (FILCPHF(I:I).EQ.' ') FILCPHF(I:I) = '_' 233 END DO 234 235 LUCPHF = -1 236 CALL WOPEN2(LUCPHF,FILCPHF,64,0) 237 238*---------------------------------------------------------------------* 239* loop over vectors, set up rhs vectors and call ABARSP to solve eqs.: 240*---------------------------------------------------------------------* 241 242 IVEC = 0 243 DO WHILE(IVEC.LT.NVEC) 244 IVEC = IVEC + 1 245 246 ISYM = 1 247 DO IDX = 1, NSTAT 248 ISYM = MULD2H(ISYM,ISYMS(IVEC,IDX)) 249 END DO 250 DO IDX = 1, ORDER 251 ISYM = MULD2H(ISYM,ISYMO(IVEC,IDX)) 252 END DO 253 254 IF (LWRK0 .LT. NALLAI(ISYM)) THEN 255 CALL QUIT('Insufficient memory in CC_CPHF.') 256 END IF 257 258 IF (LOCDBG) THEN 259 WRITE (LUPRI,*) 'CC_CPHF> GP vector, label, symmetry:', 260 & LABEL(IVEC,1), ISYM 261 WRITE(LUPRI,*) '2 MODWR ',MODWR 262 CALL FLSHFO(LUPRI) 263 END IF 264 265C --------------------------- 266C get right hand side vector: 267C --------------------------- 268 CALL CC_GETHFGD(IVEC,TYPE,LABEL,ISYMS,ISTAT,EIGV,ISYMO, 269 & FREQS,ICAU,NVEC,MAXVEC,IREAL, 270 & WORK(KCMO),WORK(KUDV),WORK(KXINDX), 271 & WORK(KFRVAL),WORK(KEND0),LWRK0) 272 273 IF (LOCDBG) THEN 274 WRITE (LUPRI,'(5X,I5,F12.8)') (I,WORK(KEND0-1+I), 275 & I=1,NALLAI(ISYM)) 276 WRITE(LUPRI,*) '3 MODWR ',MODWR 277 END IF 278 279 280C ------------------------------ 281C write gradient vector to file: 282C ------------------------------ 283 REWIND LUGDVE 284 CALL WRITT(LUGDVE,NALLAI(ISYM),WORK(KEND0)) 285 286 287C ---------------------------------- 288C check norm of the gradient vector: 289C ---------------------------------- 290 GDNORM=DSQRT(DDOT(NALLAI(ISYM),WORK(KEND0),1,WORK(KEND0),1)) 291 IF (LOCDBG) WRITE (LUPRI,*) 'GDNORM:',GDNORM 292 IF (LOCDBG) WRITE (LUPRI,*) '4 MODWR ',MODWR 293 294 295C -------------------------------------------------------- 296C for 'R1' check if several frequencies for same operator: 297C -------------------------------------------------------- 298 NFRVAL = 1 299 MAYBE_MORE = .TRUE. 300 DO WHILE (MAYBE_MORE) 301 IF (TYPE(1:3).EQ.'R1 ' 302 & .AND. IVEC.LT.NVEC .AND. NFRVAL.LT.MXFRVAL 303 & .AND. LABEL(IVEC,1).EQ.LABEL(IVEC+1,1) ) THEN 304 WORK(KFRVAL+NFRVAL) = FREQS(IVEC+1,1) 305 IVEC = IVEC + 1 306 NFRVAL = NFRVAL + 1 307 ELSE 308 MAYBE_MORE = .FALSE. 309 END IF 310 END DO 311 312 313 IF (GDNORM.GT.THRLDPHF) THEN 314C 315 CALL SETSIR(.FALSE.,WORK(KEND0),LWRK0) 316 317 IF (LOCDBG) WRITE (LUPRI,*) 'going to call ABARSP' 318 IF (LOCDBG) WRITE (LUPRI,*) '5 MODWR ',MODWR 319C ----------------------- 320C solve CPHF equation(s): 321C ----------------------- 322 IPRABA = IPRINT 323 NABAOP = 1 324 NABATY = IREAL ! flag for real/imaginary operators 325 CALL ABARSP(CICLC,HFCLC,TRPCLC,OOTV,ISYM,EXCLC, 326 & WORK(KFRVAL),NFRVAL,NABATY,NABAOP, 327 & LABEL(IVEC,1),LUGDVE,LUSOVE,LUREVE,THRLEQ, 328 & MAXITE,IPRABA,MXRM,MXPHP,WORK(KEND0),LWRK0) 329 330 IF (LOCDBG) WRITE (LUPRI,*) 'returned from ABARSP' 331 IF (LOCDBG) WRITE (LUPRI,*) '6 MODWR ',MODWR 332 333C -------------------------------------------- 334C clean up `left overs' from ABARSP: 335C -------------------------------------------- 336 CALL GPCLOSE(LUSIFC,'KEEP') 337 IF (LUONEL .GT. 0) CALL GPCLOSE(LUONEL,'KEEP') 338 IF (LUPROP .GT. 0) CALL GPCLOSE(LUPROP,'KEEP') 339C 340C -------------------------------------------- 341C read solution vector(s) and save on CC file: 342C -------------------------------------------- 343 REWIND LUSOVE 344 345 END IF 346 347 DO IFRVAL = 1, NFRVAL 348 KSLVEC = KEND0 349 KEND1 = KSLVEC + 2*MALLAI 350 LWRK1 = LWORK - KEND1 351 IF (LWRK1.LT.0) THEN 352 CALL QUIT('Insufficient memory in CC_CPHF.') 353 END IF 354 355 CALL DZERO(WORK(KSLVEC),2*MALLAI) 356 IF (GDNORM.GT.THRLDPHF) THEN 357 CALL READT(LUSOVE,2*NALLAI(ISYM),WORK(KSLVEC)) 358 END IF 359 360 361 ! save on CPHF vector file 362 IDX = IVEC - NFRVAL + IFRVAL 363 CALL PUTWA2(LUCPHF,FILCPHF,WORK(KSLVEC), 364 & 1+2*MALLAI*(IDX-1),2*MALLAI) 365 366 367 ! check if a corresponding CC vector exists 368 INUM = -1 369 IF (TYPE(1:3).EQ.'R1 ') THEN 370 INUM=IR1TAMP(LABEL(IDX,1),.TRUE.,WORK(KFRVAL-1+IFRVAL),ISYM) 371 END IF 372 IF (LOCDBG) WRITE (LUPRI,*) 'IFRVAL, TYPE, MODWR ', 373 & IFRVAL,TYPE,MODWR 374 375 ! if yes put the CPHF also on the CC vector file 376 IF (INUM.GT.0) THEN 377 CALL CC_WRRSP(TYPE,INUM,ISYM,IOPTWR,MODWR, 378 & WORK(KSLVEC),DUMMY,DUMMY,WORK(KEND1),LWRK1) 379 END IF 380 381 IF (LOCDBG) THEN 382 WRITE (LUPRI,*) 383 & 'CC_CPHF> solution vector, label, freq:', 384 & LABEL(IDX,1),WORK(KFRVAL-1+IFRVAL) 385 WRITE (LUPRI,'(5X,I5,F12.8)') 386 & (I,WORK(KSLVEC-1+I),I=1,2*NALLAI(ISYM)) 387 WRITE (LUPRI,*) 388 & 'CC_CPHF> saved CPHF solution for ',TYPE, 389 & ' equation nb. ',IDX,INUM 390 END IF 391 IF (LOCDBG) WRITE (LUPRI,*) '7 MODWR ',MODWR 392 393 END DO 394 395 END DO 396*---------------------------------------------------------------------* 397* that's it: close files, restore variables and return 398*---------------------------------------------------------------------* 399 CALL WCLOSE2(LUCPHF,FILCPHF,'KEEP') 400 401 IF (LUINTM .GT. 0) CALL GPCLOSE(LUINTM,'DELETE') 402 CALL GPCLOSE(LUGDVE,'DELETE') 403 CALL GPCLOSE(LUSOVE,'DELETE') 404 CALL GPCLOSE(LUREVE,'DELETE') 405 406 IF (LUINTA .LE. 0) THEN 407 CALL MAKE_AOTWOINT(WORK,LWORK) 408 CALL GPOPEN(LUINTA,'AOTWOINT','UNKNOWN',' ','UNFORMATTED', 409 * IDUMMY,.FALSE.) 410 END IF 411 412 NEWCMO = NEWCMO_SAVE 413 NCONF = NCOSAV 414 415 IF (LOCDBG) THEN 416 WRITE(LUPRI,*) 'leave CC_CPHF' 417 CALL FLSHFO(LUPRI) 418 END IF 419 IF (LOCDBG) WRITE (LUPRI,*) '8 MODWR ',MODWR 420 421 CALL QEXIT('CC_CPHF') 422 RETURN 423 END 424*=====================================================================* 425* END OF SUBROUTINE CC_CPHF * 426*=====================================================================* 427c /* deck cc_sirinf */ 428*=====================================================================* 429 SUBROUTINE CC_SIRINF(NCMOT1,NASHT1,N2ASHX1,LCINDX1) 430*---------------------------------------------------------------------* 431* Purpose: read some variables from SIRIUS common blocks 432*---------------------------------------------------------------------* 433#include "implicit.h" 434#include "inforb.h" 435#include "inftap.h" 436#include "infdim.h" 437 438 NCMOT1 = NCMOT 439 NASHT1 = NASHT 440 N2ASHX1 = N2ASHX 441 LCINDX1 = LCINDX 442 RETURN 443 END 444*=====================================================================* 445c /* deck cc_rdhfrsp */ 446*=====================================================================* 447 SUBROUTINE CC_RDHFRSP(LIST,IDXLST,ISYM,XKAPPA) 448*---------------------------------------------------------------------* 449C 450C Purpose: Read a CPHF response vector from file 451C for explanation of LIST, IDXLIST & MODFIL see CC_RDRSP 452C 453C Christof Haettig, summer 2003 454*---------------------------------------------------------------------* 455 IMPLICIT NONE 456#include "priunit.h" 457#include "dummy.h" 458#include "ccorb.h" 459#include "ccsdsym.h" 460#include "ccfro.h" 461 462#if defined (SYS_CRAY) 463 REAL XKAPPA(*) 464#else 465 DOUBLE PRECISION XKAPPA(*) 466#endif 467 468 CHARACTER FILCPHF*8, LIST*3 469 INTEGER MALLAI, ISYM, JSYM, LUCPHF, IDXLST 470 471 CALL QENTER('CCRDHFRSP') 472 473* set file name: 474 WRITE(FILCPHF,'(A5,A3)') 'CPHF_',LIST(1:3) 475 DO I = 6,8 476 IF (FILCPHF(I:I).EQ.' ') FILCPHF(I:I) = '_' 477 END DO 478 479* calculate record lengths: 480 MALLAI = NALLAI(1) 481 DO JSYM = 2, NSYM 482 MALLAI = MAX(MALLAI,NALLAI(JSYM)) 483 END DO 484 485* open direct access file: 486 LUCPHF = -1 487 CALL WOPEN2(LUCPHF,FILCPHF,64,0) 488 489* read vector number IDXLST from file: 490 CALL GETWA2(LUCPHF,FILCPHF,XKAPPA, 491 & 1+2*MALLAI*(IDXLST-1),2*NALLAI(ISYM)) 492 493* close file and return: 494 CALL WCLOSE2(LUCPHF,FILCPHF,'KEEP') 495 496 CALL QEXIT('CCRDHFRSP') 497 RETURN 498 END 499*=====================================================================* 500* END OF SUBROUTINE CC_RDHFRSP * 501*=====================================================================* 502