1C 2C /* Deck dc_calc */ 3 SUBROUTINE DC_CALC(MODEL,ISYMTR,NEXCI,EXVAL,LEXVAL, 4 & DENSIJ,LDENSIJ,DENSAB,LDENSAB, 5 & DENSAI,LDENSAI, 6 & T2MP,LT2MP,FOCKD,LFOCKD, 7 & WORK,LWORK) 8C 9C This routine is part of the atomic integral direct SOPPA program. 10C 11C Keld Bak, June 1997 12C Stephan P. A. Sauer: 10.11.2003: merge with Dalton 2.0 13C 14C PURPOSE: Determine "double corrected RPA" excitation energies. 15C 16#ifdef VAR_MPI 17 use so_parutils, only: parsoppa_do_eres, my_mpi_integer, 18 & soppa_nint, sop_master, one_mpi 19#endif 20 use so_info, only: sop_excita, so_full_name, so_model_number 21C 22#include "implicit.h" 23#ifdef VAR_MPI 24#include "mpif.h" 25C IRAT in order to assign space for load-balancing 26#include "iratdef.h" 27#endif 28#include "priunit.h" 29C 30#include "soppinf.h" 31#include "ccsdsym.h" 32C 33 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 34C 35 LOGICAL NONEWT 36 CHARACTER*8 PDENS_LABEL 37CPi 38 CHARACTER*5 MODEL 39 CHARACTER*11 FULL_NAME 40C 41 DIMENSION EXVAL(LEXVAL) 42 DIMENSION DENSIJ(LDENSIJ), DENSAB(LDENSAB), DENSAI(LDENSAI) 43 DIMENSION T2MP(LT2MP), FOCKD(LFOCKD), WORK(LWORK) 44C 45 INTEGER NIT, NOLDTR, NNEWTR, IDTYPE, IMODEL 46#ifdef VAR_MPI 47 INTEGER CP_ISYMTR 48 INTEGER MAXNUMJOBS 49C This array is only there to ensure that the four above variables 50C are allocated consecutively, so that it can be send together. Only 51C use it for this purpose. 52C The definition must match that in soppa_nodedriver 53 INTEGER INFO_ARRAY(6) 54 EQUIVALENCE (info_array(1), cp_isymtr), (info_array(2),nit), 55 & (info_array(3), nnewtr), (info_array(4),noldtr), 56 & (info_array(5), idtype), (info_array(6),imodel) 57 INTEGER(MPI_INTEGER_KIND) :: ierr_mpi, numprocs_mpi 58C IF (numprocs_mpi .GT. 1) CALL PARQUIT('RPA(D)') 59 IDTYPE = 0 60#endif 61C 62C------------------ 63C Add to trace. 64C------------------ 65C 66 CALL QENTER('DC_CALC') 67C 68 FULL_NAME = SO_FULL_NAME(MODEL) 69 IMODEL = SO_MODEL_NUMBER(MODEL) 70C 71 KEND1 = 1 72 LWORK1 = LWORK 73C 74C-------------------------------------------------------------- 75C Initialize iteration counter, number of old trialvectors, 76C and number of new trialvectors. 77C-------------------------------------------------------------- 78C 79 NIT = 1 80C 81 NOLDTR = 0 82C 83 NNEWTR = NEXCI 84C 85#ifdef VAR_MPI 86C------------------------------------------------------------------ 87C For MPI, we need some space in which to store the indices each 88C process is to work with in so_eres. 89C------------------------------------------------------------------ 90C 91 call mpi_comm_size( mpi_comm_world, numprocs_mpi, ierr_mpi) 92 maxnumjobs = soppa_nint - min(soppa_nint, numprocs_mpi) + 1 93 if ( numprocs_mpi .eq. 1 ) then 94C Not a real parallel job, don't bother 95 lAssignedIndices = 1 96 kAssignedIndices = 0 97 else 98 lAssignedIndices = (maxnumjobs + irat-1) /IRAT 99 kAssignedIndices = KEND1 100 KEND1 = kAssignedIndices + lAssignedIndices 101 LWORK1 = LWORK - KEND1 102 CALL SO_MEMMAX ('DC_CALC.1A',LWORK1) 103 IF (LWORK1 .LT. 0) CALL STOPIT('DC_CALC.1A',' ',KEND1,LWORK) 104 endif 105#endif 106C------------------ 107C Write banner. 108C------------------ 109CPi 110 WRITE(LUPRI,'(/,2X,A)') '*********************************'// 111 & '*********************************' 112 WRITE(LUPRI,'(11X,A,A,I1)') ADJUSTR(FULL_NAME), 113 & ' calculation, Excitation symmetry ',ISYMTR 114 WRITE(LUPRI,'(2X,A)') '*********************************'// 115 & '*********************************' 116C 117C--------------------------------------------------------- 118C Make E[2] linear transformation of RPA eigenvectors. 119C--------------------------------------------------------- 120C 121#ifdef VAR_MPI 122C In parallel, send slaves to so_eres 123C 124 call mpi_bcast( parsoppa_do_eres, one_mpi, 125 & my_mpi_integer, sop_master, 126 & mpi_comm_world, ierr_mpi ) 127C ISYMTR is a non-local parameter, we need to copy it to the info-array 128 CP_ISYMTR = ISYMTR 129 CALL MPI_BCAST( INFO_ARRAY, 6_mpi_integer_kind, MY_MPI_INTEGER, 130 & sop_master, MPI_COMM_WORLD, ierr_mpi) 131#endif 132CPi 18.08.16: Copy from so_lrsolv 133 CALL GETTIM (DUMMY,WTIMES) 134 DTIME = SECOND() 135 CALL SO_ERES(MODEL,0,NEXCI,DENSIJ,LDENSIJ,DENSAB,LDENSAB, 136 & T2MP,LT2MP,FOCKD,LFOCKD,DENSAI,LDENSAI,1,ISYMTR, 137 & 0, 138#ifdef VAR_MPI 139 & WORK(kAssignedIndices), maxnumjobs, 140#endif 141 & WORK(KEND1),LWORK1) 142 DTIME = SECOND() - DTIME 143 SOTIME(35) = SOTIME(35) + DTIME 144 CALL GETTIM (DUMMY,WTIMEE) 145 SOWTIM(1) = SOWTIM(1) + WTIMEE - WTIMES 146C 147C----------------------------------------------------------- 148C Make S[2] linear transformation of trialvectors giving 149C resultvectors. 150C----------------------------------------------------------- 151C 152 DTIME = SECOND() 153 CALL DC_SRES(0,NEXCI,DENSIJ,LDENSIJ,DENSAB,LDENSAB, 154 & ISYMTR,WORK,LWORK) 155 DTIME = SECOND() - DTIME 156 SOTIME(40) = SOTIME(40) + DTIME 157C 158C----------------------------------------------------------------- 159C Calculate diagonal parts of E[2] and S[2] and write to disk. 160C----------------------------------------------------------------- 161C 162 DTIME = SECOND() 163 CALL SO_DIAG(MODEL,FOCKD,LFOCKD,DENSIJ,LDENSIJ,DENSAB,LDENSAB, 164 & ISYMTR,WORK,LWORK) 165 DTIME = SECOND() - DTIME 166 SOTIME(31) = SOTIME(31) + DTIME 167C 168C------------------------------ 169C Allocation of work space. 170C------------------------------ 171C 172 LTR1E = NT1AM(ISYMTR) 173 LTR1D = NT1AM(ISYMTR) 174 LTR2E = N2P2HOP(ISYMTR) 175 LTR2D = N2P2HOP(ISYMTR) 176C 177 KTR1 = 1 178 KRES1 = KTR1 + LTR1E 179 KRESO1 = KRES1 + LTR1E 180 KRES2 = KRESO1 + LTR1E 181 KEND1 = KRES2 + LTR2E 182 LWORK1 = LWORK - KEND1 183C 184 CALL SO_MEMMAX ('DC_CALC',LWORK1) 185 IF (LWORK1 .LT. 0) CALL STOPIT('DC_CALC',' ',KEND1,LWORK) 186C 187C---------------- 188C Open files. 189C---------------- 190C 191 CALL SO_OPEN(LUTR1E,FNTR1E,LTR1E) 192 CALL SO_OPEN(LUTR1D,FNTR1D,LTR1D) 193 CALL SO_OPEN(LUTR2E,FNTR2E,LTR2E) 194 CALL SO_OPEN(LUTR2D,FNTR2D,LTR2D) 195 CALL SO_OPEN(LURS1E,FNRS1E,LTR1E) 196 CALL SO_OPEN(LURS1D,FNRS1D,LTR1D) 197 CALL SO_OPEN(LURO1E,FNRO1E,LTR1E) 198 CALL SO_OPEN(LURO1D,FNRO1D,LTR1D) 199 CALL SO_OPEN(LURS2E,FNRS2E,LTR2E) 200 CALL SO_OPEN(LURS2D,FNRS2D,LTR2D) 201C 202C--------------------------- 203C Loop over excitations. 204C--------------------------- 205 WRITE(LUPRI,9010) 206 WRITE(LUPRI,9011) 207 WRITE(LUPRI,9010) 208C 209 DO 100 IEXCI = 1,NEXCI 210C 211C----------------------------------------------------------------------- 212C Calculate 1p-1h part of double corrected RPA excitation energy. 213C----------------------------------------------------------------------- 214C 215 CALL SO_READ(WORK(KTR1), LTR1E,LUTR1E,FNTR1E,IEXCI) 216 CALL SO_READ(WORK(KRES1), LTR1E,LURS1E,FNRS1E,IEXCI) 217 CALL SO_READ(WORK(KRESO1),LTR1E,LURO1E,FNRO1E,IEXCI) 218C 219 IF (MODEL .EQ. 'DCRPA' .OR. MODEL .EQ. 'SDCHR') THEN 220C 221 CALL DAXPY(LTR1E,-EXVAL(IEXCI),WORK(KRESO1),1,WORK(KRES1),1) 222C 223 END IF 224C 225 EX1P1H = DDOT(LTR1E,WORK(KRES1),1,WORK(KTR1),1) 226C 227C----------------------------------------------------------------------- 228C Calculate 1h-1p part of double corrected RPA excitation energy. 229C----------------------------------------------------------------------- 230C 231 CALL SO_READ(WORK(KTR1), LTR1D,LUTR1D,FNTR1D,IEXCI) 232 CALL SO_READ(WORK(KRES1), LTR1D,LURS1D,FNRS1D,IEXCI) 233 CALL SO_READ(WORK(KRESO1),LTR1D,LURO1D,FNRO1D,IEXCI) 234C 235 IF (MODEL .EQ. 'DCRPA' .OR. MODEL .EQ. 'SDCHR') THEN 236C 237 CALL DAXPY(LTR1D,-EXVAL(IEXCI),WORK(KRESO1),1,WORK(KRES1),1) 238C 239 END IF 240C 241 EX1H1P = DDOT(LTR1D,WORK(KRES1),1,WORK(KTR1),1) 242C 243C----------------------------------------------------------------------- 244C Calculate 2p-2h part of double corrected RPA excitation energy. 245C----------------------------------------------------------------------- 246C 247 CALL SO_READ(WORK(KRES2), LTR2E,LURS2E,FNRS2E,IEXCI) 248C 249 IF (TRIPLET) THEN 250C 251 CALL DC_OMEC(EX2P2H,WORK(KRES2),LTR2E,EXVAL(IEXCI), 252 & ISYMTR,WORK(KEND1),LWORK1) 253C 254 ELSE 255C 256 CALL CC_OMEC(EX2P2H,WORK(KRES2),EXVAL(IEXCI), 257 & WORK(KEND1),LWORK1,ISYMTR) 258C 259C-------------------------------------------------------------------- 260C Calculate the first order 2p-2h eigenvector in RPA(D) theory 261C and write to output. 262C-------------------------------------------------------------------- 263C 264 CALL DC_R1VEC(WORK(KRES2), LTR2E,EXVAL(IEXCI),ISYMTR, 265 & WORK(KEND1),LWORK1) 266C 267 END IF 268C 269 CALL SO_WRITE(WORK(KRES2), LTR2E, LUTR2E,FNTR2E, IEXCI) 270C 271C----------------------------------------------------------------------- 272C Calculate 2h-2p part of double corrected RPA excitation energy. 273C----------------------------------------------------------------------- 274C 275 CALL SO_READ(WORK(KRES2), LTR2D,LURS2D,FNRS2D,IEXCI) 276C 277 IF (TRIPLET) THEN 278C 279 CALL DC_OMEC(EX2H2P,WORK(KRES2),LTR2D,EXVAL(IEXCI), 280 & ISYMTR,WORK(KEND1),LWORK1) 281C 282 ELSE 283C 284 CALL CC_OMEC(EX2H2P,WORK(KRES2),EXVAL(IEXCI), 285 & WORK(KEND1),LWORK1,ISYMTR) 286C 287C-------------------------------------------------------------------- 288C Calculate the first order 2h-2p eigenvector in RPA(D) theory 289C and write to output. 290C-------------------------------------------------------------------- 291C 292 CALL DC_R1VEC(WORK(KRES2), LTR2D,EXVAL(IEXCI),ISYMTR, 293 & WORK(KEND1),LWORK1) 294C 295 END IF 296C 297 CALL SO_WRITE(WORK(KRES2), LTR2D, LUTR2D,FNTR2D, IEXCI) 298C 299C----------------------------------------------------------- 300C Add contributions to give RPA(D) excitation energy. 301C----------------------------------------------------------- 302C 303 IF (TRIPLET) THEN 304C 305 WRITE(LUPRI,9012) IEXCI, 306 & EX1P1H + EX1H1P + EX2P2H + EX2H2P, 307 & EXVAL(IEXCI), 308 & (EX1P1H + EX1H1P) - EXVAL(IEXCI), 309 & EX2P2H + EX2H2P 310C 311 EXVAL(IEXCI) = EX1P1H + EX1H1P + EX2P2H + EX2H2P 312C 313 ELSE 314C 315 WRITE(LUPRI,9012) IEXCI, 316 & EX1P1H + EX1H1P + (EX2P2H + EX2H2P)/TWO, 317 & EXVAL(IEXCI), 318 & (EX1P1H + EX1H1P)- EXVAL(IEXCI), 319 & (EX2P2H + EX2H2P)/TWO 320C 321 EXVAL(IEXCI) = EX1P1H + EX1H1P + EX2P2H/TWO + EX2H2P/TWO 322C 323 END IF 324C 325 100 CONTINUE 326C 327 WRITE(LUPRI,9010) 328C 329C----------------- 330C Close files. 331C----------------- 332C 333 CALL SO_CLOSE(LUTR1E,FNTR1E,'KEEP') 334 CALL SO_CLOSE(LUTR1D,FNTR1D,'KEEP') 335 CALL SO_CLOSE(LUTR2E,FNTR2E,'KEEP') 336 CALL SO_CLOSE(LUTR2D,FNTR2D,'KEEP') 337 CALL SO_CLOSE(LURS1E,FNRS1E,'KEEP') 338 CALL SO_CLOSE(LURS1D,FNRS1D,'KEEP') 339 CALL SO_CLOSE(LURO1E,FNRO1E,'KEEP') 340 CALL SO_CLOSE(LURO1D,FNRO1D,'KEEP') 341 CALL SO_CLOSE(LURS2E,FNRS2E,'KEEP') 342 CALL SO_CLOSE(LURS2D,FNRS2D,'KEEP') 343C 344C-------------------------------------------------------------- 345C Calculate and save the pertubed density matrix 346C------------------------------------------------------------- 347CPi Solution vectors are not normalized as in so_rpslex based on 348CPi excitation weights !!! 349C 350 WRITE(PDENS_LABEL,'(A7,I1)') 'EXCITA ',ISYMTR 351 CALL SO_PERTDENS('DCRPA',SOP_EXCITA,NEXCI, 352 & EXVAL,NEXCI,PDENS_LABEL, 353 & ISYMTR,.FALSE.,1.0D0, 354 & T2MP,LT2MP,DENSIJ,LDENSIJ, 355 & DENSAB,LDENSAB,DENSAI,LDENSAI, 356 & WORK(KEND1),LWORK1) 357C 358C 359C------------------------------------ 360C Flush the standard output file. 361C------------------------------------ 362C 363 CALL FLSHFO(LUPRI) 364C 365C----------------------- 366C Remove from trace. 367C----------------------- 368C 369 CALL QEXIT('DC_CALC') 370C 371 RETURN 372C 373 9010 FORMAT(1X,'----------------------------------------', 374 & '-------------------------') 375 9011 FORMAT(1X,'Excitation Energy (au) ph(0) ', 376 & ' ph(2) 2p2h(2)') 377 9012 FORMAT(2X,I5,6X,F12.8,1X,F12.8,1X,F12.8,1X,F12.8) 378C 379 END 380