1! 2!... Copyright (c) 2015 by the authors of Dalton (see below). 3!... All Rights Reserved. 4!... 5!... The source code in this file is part of 6!... "Dalton, a molecular electronic structure program, 7!... Release DALTON2015 (2015), see http://daltonprogram.org" 8!... 9!... This source code is provided under a written licence and may be 10!... used, copied, transmitted, or stored only in accord with that 11!... written licence. 12!... 13!... In particular, no part of the source code or compiled modules may 14!... be distributed outside the research group of the licence holder. 15!... This means also that persons (e.g. post-docs) leaving the research 16!... group of the licence holder may not take any part of Dalton, 17!... including modified files, with him/her, unless that person has 18!... obtained his/her own licence. 19!... 20!... For further information, including how to get a licence, see: 21!... http://daltonprogram.org 22! 23! 24C 25c /* deck CC_EOM_XOPA */ 26*=====================================================================* 27 SUBROUTINE CC_EOM_XOPA(WORK,LWORK) 28*---------------------------------------------------------------------* 29* 30* Purpose: direct calculation of first-order transition properties 31* (transition moments and oscillator strengths) 32* for transitions between two excited states with the 33* Coupled Cluster models 34* 35* CCS, CC2, CCSD, CC3 36* 37* and partially with SCF and CIS 38* 39* Written by Christof Haettig winter 2002/2003. 40* Modified version by Sonia, to include EOM models, 2015 41*=====================================================================* 42 IMPLICIT NONE 43#include "priunit.h" 44#include "cclists.h" 45#include "ccxopainf.h" 46#include "ccsdinp.h" 47#include "dummy.h" 48#include "second.h" 49#include "ccexcinf.h" 50#include "ccorb.h" 51 52* local parameters: 53 CHARACTER*(16) MSGDBG 54 PARAMETER (MSGDBG = '[debug] CC_EOM_XOPA> ') 55 56#if defined (SYS_CRAY) 57 REAL ZERO 58#else 59 DOUBLE PRECISION ZERO 60#endif 61 PARAMETER (ZERO = 0.0d0) 62 63 CHARACTER*10 MODEL 64 INTEGER LWORK 65 66#if defined (SYS_CRAY) 67 REAL WORK(LWORK) 68 REAL TIM0, TIM1, TIMF, TIMXE1, TIMXE2 69#else 70 DOUBLE PRECISION WORK(LWORK) 71 DOUBLE PRECISION TIM0, TIM1, TIMF, TIMXE1, TIMXE2 72#endif 73 74 LOGICAL LADD 75 INTEGER NBOPA, MXFTRAN, MXATRAN, MXXTRAN, MXFVEC, MXAVEC, MXXVEC, 76 & NFTRAN, NXE1TRAN, NXE2TRAN, NSTATES, 77 & KRESULT, KFTRAN, KFDOTS, KFCONS, KEND0, LEND0, 78 & KE1TRAN, KE1DOTS, KE1CONS, 79 & KX2TRAN, KX2DOTS, KX2CONS, 80 & IOPT, IORDER, ISYM 81 INTEGER KEOMX2TRAN, KEOMX2DOTS, KEOMX2CONS 82 INTEGER KEOML0TRAN, KEOML0DOTS, KEOML0CONS 83 INTEGER NEOMXE2TRAN, MXEOMXTRAN, MXEOMXVEC 84 INTEGER NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC 85 86* external functions: none 87 88*---------------------------------------------------------------------* 89* print header for second-order property section: 90*---------------------------------------------------------------------* 91 WRITE (LUPRI,'(7(/1X,2A),/)') 92 & '************************************', 93 & '******************************', 94 & '* ', 95 & ' *', 96 & '*<<<<<< OUTPUT FROM COUPLED CLUST', 97 & 'ER LINEAR RESPONSE >>>>>>>*', 98 & '*<<<<<< CALCULATION OF ONE-PHOTON A', 99 & 'BSORPTION STRENGTHS >>>>>>>*', 100 & '*<<<<<< FOR EXCITED TO EXCITED S', 101 & 'TATE TRANSITIONS >>>>>>>*', 102 & '* ', 103 & ' *', 104 & '************************************', 105 & '******************************' 106 107*---------------------------------------------------------------------* 108 IF (.NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 109 CALL QUIT('CC_EOM_XOPA called for unknown Coupled Cluster.') 110 END IF 111 112* print some debug/info output 113 IF(IPRINT .GT. 10) WRITE(LUPRI,*) 'CC_EOM_XOPA Workspace:',LWORK 114 115 TIM0 = SECOND() 116 117*---------------------------------------------------------------------* 118* allocate & initialize work space for property contributions: 119*---------------------------------------------------------------------* 120 ! maximum number of transition moments to compute 121 NBOPA = 2 * NQR2OP * NXQR2ST 122 123 ! number of excited states 124 NSTATES = 0 125 DO ISYM = 1, NSYM 126 NSTATES = NSTATES + NCCEXCI(ISYM,1) + NCCEXCI(ISYM,3) 127 END DO 128 129 ! maximum number of transformations or vector calculations 130 ! NSTATES * NQR2OP LE x Eta{X} transformations 131 ! NQR2OP Xi{X} vectors 132 ! 2*NXQR2ST LE x B x RE transformations 133 MXATRAN = NSTATES * NQR2OP 134 MXXTRAN = NQR2OP 135 MXEOMXTRAN = NQR2OP 136 MXEOML0TRAN = 1 137 MXFTRAN = 2*NXQR2ST 138 139 ! maximum number of vectors to dot on 140 ! NSTATES RE vectors dotted on a LE x Eta{X} transformation 141 ! 2*NXQR2ST N2 vectors dotted on a Xi{X} vector 142 ! NQR2OP R1 vectors dotted on a LE x B x RE transformation 143 MXAVEC = NSTATES 144 MXXVEC = 2*NXQR2ST 145 MXEOMXVEC = 2*NXQR2ST !I am not sure I understand why this number... 146 MXFVEC = NQR2OP 147 MXEOML0VEC = NSTATES 148 149 KRESULT = 1 150 KEND0 = KRESULT + NBOPA 151 152 KFTRAN = KEND0 153 KFDOTS = KFTRAN + MXFTRAN * MXDIM_FTRAN 154 KFCONS = KFDOTS + MXFVEC * MXFTRAN 155 KEND0 = KFCONS + MXFVEC * MXFTRAN 156 157 KE1TRAN = KEND0 158 KE1DOTS = KE1TRAN + MXATRAN * MXDIM_XEVEC 159 KE1CONS = KE1DOTS + MXAVEC * MXATRAN 160 KEND0 = KE1CONS + MXAVEC * MXATRAN 161 162 KX2TRAN = KEND0 163 KX2DOTS = KX2TRAN + MXXTRAN * MXDIM_XEVEC 164 KX2CONS = KX2DOTS + MXXVEC * MXXTRAN 165 KEND0 = KX2CONS + MXXVEC * MXXTRAN 166 167 KEOMX2TRAN= KEND0 168 KEOMX2DOTS= KEOMX2TRAN + MXEOMXTRAN * MXDIM_XEVEC 169 KEOMX2CONS= KEOMX2DOTS + MXEOMXVEC * MXEOMXTRAN 170 KEND0 = KEOMX2CONS + MXEOMXVEC * MXEOMXTRAN 171 172 KEOML0TRAN= KEND0 173 KEOML0DOTS= KEOML0TRAN + MXEOML0TRAN 174 KEOML0CONS= KEOML0DOTS + MXEOML0VEC * MXEOML0TRAN 175 KEND0 = KEOML0CONS + MXEOML0VEC * MXEOML0TRAN 176 177 LEND0 = LWORK - KEND0 178 IF (LEND0 .LT. 0) THEN 179 CALL QUIT('Insufficient memory in CC_EOM_XOPA. (1)') 180 END IF 181 182 CALL DZERO(WORK(KRESULT),NBOPA) 183 184 !sonia 185 CALL DZERO(WORK(KEOMX2CONS),MXEOMXVEC * MXEOMXTRAN) 186 187*---------------------------------------------------------------------* 188* set up lists for F transformations, ETA{O} and Xi{O} vectors: 189*---------------------------------------------------------------------* 190 LADD = .FALSE. 191 192 CALL CCXOPA_EOMSETUP(WORK(KFTRAN),WORK(KFDOTS),WORK(KFCONS), 193 & NFTRAN,MXFTRAN,MXFVEC, 194 & WORK(KE1TRAN),WORK(KE1DOTS),WORK(KE1CONS), 195 & NXE1TRAN,MXATRAN,MXAVEC, 196 & WORK(KX2TRAN),WORK(KX2DOTS),WORK(KX2CONS), 197 & NXE2TRAN,MXXTRAN,MXXVEC, 198 & WORK(KEOMX2TRAN),WORK(KEOMX2DOTS), 199 & WORK(KEOMX2CONS), 200 & NEOMXE2TRAN,MXEOMXTRAN,MXEOMXVEC, 201 & WORK(KEOML0TRAN),WORK(KEOML0DOTS), 202 & WORK(KEOML0CONS), 203 & NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC, 204 & WORK(KRESULT),NBOPA,LADD,WORK(KEND0),LEND0) 205 206*---------------------------------------------------------------------* 207* calculate F matrix contributions: 208*---------------------------------------------------------------------* 209 TIM1 = SECOND() 210 211 CALL DZERO(WORK(KFCONS),MXFVEC*NFTRAN) 212 213 IOPT = 5 214 CALL CC_FMATRIX(WORK(KFTRAN),NFTRAN,'LE ','RE ',IOPT,'R1 ', 215 & WORK(KFDOTS),WORK(KFCONS),MXFVEC, 216 & WORK(KEND0), LEND0) 217 218 TIMF = SECOND() - TIM1 219 220 IF (NFTRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 221 & '>>> Time used for',NFTRAN,' F matrix transformations:',TIMF 222 CALL FLSHFO(LUPRI) 223 224*---------------------------------------------------------------------* 225* calculate LE x A{O} x RE contributions: 226*---------------------------------------------------------------------* 227 TIM1 = SECOND() 228 229 CALL DZERO(WORK(KE1CONS),MXAVEC*NXE1TRAN) 230 231 232 233 IOPT = 5 234 IORDER = 1 235 if (leomxopa) then 236 CALL CCEOM_XIETA( WORK(KE1TRAN),NXE1TRAN,IOPT,IORDER,'LE ', 237 & '---',IDUMMY, DUMMY, 238 & 'RE ',WORK(KE1DOTS),WORK(KE1CONS), 239 & MXAVEC, WORK(KEND0), LEND0 ) 240 else 241 CALL CC_XIETA( WORK(KE1TRAN), NXE1TRAN, IOPT, IORDER, 'LE ', 242 & '---',DUMMY,DUMMY, 243 & 'RE ',WORK(KE1DOTS),WORK(KE1CONS), 244 & .FALSE.,MXAVEC, WORK(KEND0), LEND0 ) 245 end if 246 247 TIMXE1 = SECOND() - TIM1 248 IF (NXE1TRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 249 & '>>> Time used for',NXE1TRAN,' A{X} matrix transformations:', 250 & TIMXE1 251 CALL FLSHFO(LUPRI) 252 253*---------------------------------------------------------------------* 254* calculate N2 x Xksi{O} vector contributions: 255*---------------------------------------------------------------------* 256 TIM1 = SECOND() 257 258 CALL DZERO(WORK(KX2CONS),MXXVEC*NXE2TRAN) 259 260 IOPT = 5 261 IORDER = 1 262 CALL CC_XIETA( WORK(KX2TRAN), NXE2TRAN, IOPT, IORDER, '---', 263 & 'N2 ',WORK(KX2DOTS),WORK(KX2CONS), 264 & '---',IDUMMY,DUMMY, 265 & .FALSE.,MXXVEC, WORK(KEND0), LEND0 ) 266 267 TIMXE2 = SECOND() - TIM1 268 IF (NXE2TRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 269 & '>>> Time used for',NXE2TRAN,' O1/X1 vector calculation:',TIMXE2 270 CALL FLSHFO(LUPRI) 271 272*---------------------------------------------------------------------* 273* calculate LE x Xksi{O} vector contributions: 274* calculate RE x Tbar0 vector contributions: 275* multiply them together to get the EOM contrib 276*---------------------------------------------------------------------* 277 TIM1 = SECOND() 278 279 if (leomxopa) then 280 CALL DZERO(WORK(KEOMX2CONS),MXEOMXVEC*NEOMXE2TRAN) 281 CALL DZERO(WORK(KEOML0CONS),MXEOML0VEC*NEOML0TRAN) 282 283 IOPT = 5 284 IORDER = 1 285 CALL CC_XIETA( WORK(KEOMX2TRAN), NEOMXE2TRAN, 286 & IOPT, IORDER, 'L0 ', 287 & 'LE ',WORK(KEOMX2DOTS),WORK(KEOMX2CONS), 288 & '---',IDUMMY,DUMMY, 289 & .FALSE.,MXEOMXVEC, WORK(KEND0), LEND0 ) 290 291 CALL CC_DOTDRV('L0 ','RE ',NEOML0TRAN,MXEOML0VEC, 292 & WORK(KEOML0TRAN), WORK(KEOML0DOTS), 293 & WORK(KEOML0CONS), 294 & WORK(KEND0), LEND0 ) 295 296 TIMXE2 = SECOND() - TIM1 297 IF (NEOMXE2TRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")') 298 & '>>> Time used for',NEOMXE2TRAN, 299 & ' O1/X1 vector calculation:',TIMXE2 300 CALL FLSHFO(LUPRI) 301 end if 302 303*---------------------------------------------------------------------* 304* collect contributions and sum them up to the final results: 305*---------------------------------------------------------------------* 306 LADD = .TRUE. 307 308 CALL CCXOPA_EOMSETUP(WORK(KFTRAN),WORK(KFDOTS),WORK(KFCONS), 309 & NFTRAN,MXFTRAN,MXFVEC, 310 & WORK(KE1TRAN),WORK(KE1DOTS),WORK(KE1CONS), 311 & NXE1TRAN,MXATRAN,MXAVEC, 312 & WORK(KX2TRAN),WORK(KX2DOTS),WORK(KX2CONS), 313 & NXE2TRAN,MXXTRAN,MXXVEC, 314 & WORK(KEOMX2TRAN),WORK(KEOMX2DOTS), 315 & WORK(KEOMX2CONS), 316 & NEOMXE2TRAN,MXEOMXTRAN,MXEOMXVEC, 317 & WORK(KEOML0TRAN),WORK(KEOML0DOTS), 318 & WORK(KEOML0CONS), 319 & NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC, 320 & WORK(KRESULT),NBOPA,LADD,WORK(KEND0),LEND0) 321 322*---------------------------------------------------------------------* 323* print timing: 324*---------------------------------------------------------------------* 325 WRITE (LUPRI,'(/A,I4,A,F12.2," seconds.")') '>>> Total time for', 326 & NBOPA,' quadratic response func.:', SECOND() - TIM0 327 328*---------------------------------------------------------------------* 329* print one-photon absorption properties and return: 330*---------------------------------------------------------------------* 331 CALL CCOPAPRT(WORK(KRESULT),.TRUE.,NQR2OP,NXQR2ST) 332 333 CALL FLSHFO(LUPRI) 334 335 RETURN 336 END 337 338*=====================================================================* 339* END OF SUBROUTINE CC_EOM_XOPA * 340*=====================================================================* 341