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_AAMAT */ 20*=====================================================================* 21 SUBROUTINE CC_AAMAT(IATRAN,NATRAN,LISTO,LISTR,IOPTRES,FILAAM, 22 & IADOTS,ACONS,MXVEC,WORK,LWORK) 23*---------------------------------------------------------------------* 24* 25* Purpose: driver for a list of A{O} matrix transformations 26* 27* IOPTRES = 3 : each result vector is written to its own 28* file by a call to CC_WRRSP, FILAAM is used 29* as list type and IATRAN(3,*) as list index 30* 31* IOPTRES = 4 : each result vector is added to a vector on 32* file by a call to CC_WARSP, FILAAM is used 33* as list type and IATRAN(3,*) as list index 34* 35* IOPTRES = 5 : the result vectors are dotted on an array 36* of vectors, the type of the arrays is given 37* by FILAAM and the indices are taken from 38* IADOTS. the results are returned in ACONS 39* 40* N.B.: this is a very first quick hack solution... 41* should finally be incorporated in the CCCR_AA routine... 42* 43* Written by Christof Haettig, maj 1996. 44* 45*=====================================================================* 46 IMPLICIT NONE 47#include "priunit.h" 48#include "ccorb.h" 49#include "ccsdinp.h" 50#include "ccsdsym.h" 51#include "ccroper.h" 52#include "cclists.h" 53#include "ccnoddy.h" 54#include "dummy.h" 55 56 CHARACTER*1 CDUMMY 57 PARAMETER (CDUMMY = ' ') 58 59* local parameters: 60 CHARACTER*(19) MSGDBG 61 PARAMETER (MSGDBG = '[debug] CC_AAMAT> ') 62 LOGICAL LOCDBG 63 PARAMETER (LOCDBG = .FALSE.) 64 65 INTEGER LWORK, NATRAN, MXVEC 66 67#if defined (SYS_CRAY) 68 REAL WORK(LWORK), ACONS(MXVEC,NATRAN) 69 REAL ZERO, TWO, FREQ 70#else 71 DOUBLE PRECISION WORK(LWORK), ACONS(MXVEC,NATRAN) 72 DOUBLE PRECISION ZERO, TWO, FREQ 73#endif 74 PARAMETER (ZERO = 0.0d0, TWO = 2.0D0) 75 76 CHARACTER LISTO*(*), LISTR*(*), FILAAM*(*) 77 CHARACTER*(10) MODEL 78 CHARACTER*(8) LABEL 79 CHARACTER APROXR12*3 80 LOGICAL LTWOEL 81 INTEGER IATRAN(MXDIM_AATRAN,NATRAN), IADOTS(MXVEC,NATRAN) 82 INTEGER IOPER, IDLSTR, IFILE, IVEC, ITRAN, IOPT, IOPTRES 83 INTEGER ISYOPE, ISYTAM, ISYRES, IRELAX, IOPTE 84 INTEGER KEND1, LWRK1, KRES1, KRES2, LEN, KRES1EFF, KRES2EFF 85 INTEGER LENMOD,IOPTWR12,KRESR12 86 87 88* external functions: 89 INTEGER ILSTSYM 90 REAL*8 DDOT, FREQLST 91 92*---------------------------------------------------------------------* 93* initializations: 94*---------------------------------------------------------------------* 95 CALL QENTER('CC_AAMAT') 96 IF (LOCDBG) THEN 97 WRITE(LUPRI,*) 'Entered cc_aamat...' 98 CALL FLSHFO(LUPRI) 99 END IF 100 101 IF ( .NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 102 WRITE (LUPRI,'(/1x,a)') 'CC_AAMAT Called for a Coupled Cluster ' 103 & //'method not implemented in CC_AAMAT...' 104 CALL QUIT('Unknown CC method in CC_AAMAT.') 105 END IF 106 107* check list combinations: 108 IF ( (LISTO(1:2).NE.'o1' .AND. LISTO(1:2).NE.'o2') 109 & .OR. ( LISTR(1:1).NE.'R'.AND. LISTR(1:2).NE.'ER')) THEN 110 WRITE (LUPRI,*) 111 & 'A{O} matrix driver called with a strange list set:' 112 WRITE (LUPRI,*) 'LISTO:',LISTO 113 WRITE (LUPRI,*) 'LISTR:',LISTR 114 CALL QUIT('A{O} mat. driver called with a strange '// 115 & 'set of list types.') 116 END IF 117 118 IF ( LISTO(1:2).EQ.'O2' .AND. NATRAN.NE.0) THEN 119 WRITE (LUPRI,*) 120 & 'LISTO="o2" case not yet implemented in CC_AAMAT.' 121 CALL QUIT('LISTO="o2" case not yet implemented in CC_AAMAT.') 122 END IF 123 124*---------------------------------------------------------------------* 125* start loop over all requested A{O} matrix transformations 126*---------------------------------------------------------------------* 127 DO ITRAN = 1, NATRAN 128 129 IOPER = IATRAN(1,ITRAN) 130 IDLSTR = IATRAN(2,ITRAN) 131 IFILE = IATRAN(3,ITRAN) 132 IRELAX = IATRAN(4,ITRAN) 133 134 LABEL = LBLOPR(IOPER) 135 LTWOEL = LPDBSOP(IOPER) 136 137 ISYOPE = ILSTSYM(LISTO,IOPER) 138 ISYTAM = ILSTSYM(LISTR,IDLSTR) 139 ISYRES = MULD2H(ISYOPE,ISYTAM) 140 141 IF ((IOPTRES.EQ.3.OR.IOPTRES.EQ.4) .AND. 142 & ILSTSYM(FILAAM,IFILE).NE.ISYRES) THEN 143 CALL QUIT('Symmetry mismatch in CC_AAMAT.') 144 END IF 145 146*---------------------------------------------------------------------* 147* calculate A{O} matrix times a response amplitude vector: 148*---------------------------------------------------------------------* 149 150 ! allocate memory for results vector 151 KRES1 = 1 152 KRES2 = KRES1 + NT1AM(ISYRES) 153 KEND1 = KRES2 + NT2AM(ISYRES) 154 IF (CCR12) THEN 155 KRESR12 = KRES2 + NT2AM(ISYRES) 156 KEND1 = KRESR12 + NTR12AM(ISYRES) 157 END IF 158 LWRK1 = LWORK - KEND1 159 IF (LWRK1.LT.0) CALL QUIT('Insufficient memory in CC_AAMAT.') 160 161 IF ( (IRELAX.GT.0) .OR. LTWOEL ) THEN 162 163 IF (CCSDT) CALL QUIT('relaxed with triples not available!') 164 165 CALL CC_FDAAMAT(LISTR,IDLSTR,WORK(KRES1),WORK(KRES2), 166 & LABEL,IRELAX,WORK(KEND1),LWRK1) 167 168 ELSE 169 170c -------------------------------------------------- 171c do the CCS/CC2/CCSD parts: 172c -------------------------------------------------- 173 CALL CCCR_AA(LABEL,ISYOPE,LISTR,IDLSTR,DUMMY,WORK,LWORK) 174 CALL CCLR_DIASCL(WORK(KRES2),2.0d0,ISYRES) 175 176c -------------------------------------------------- 177c do the triples parts: 178c -------------------------------------------------- 179 IF (CCSDT) THEN 180 ! find frequency associates with result vector 181 IF (IOPTRES.EQ.5) THEN 182 FREQ = 0.0D0 183 ELSE 184 FREQ = FREQLST(FILAAM,IFILE) 185 END IF 186 187 ! allocate memory for effective results vector 188 KRES1EFF = KEND1 189 KRES2EFF = KRES1EFF + NT1AM(ISYRES) 190 KEND1 = KRES2EFF + NT2AM(ISYRES) 191 LWRK1 = LWORK - KEND1 192 IF (LWRK1.LT.0) 193 & CALL QUIT('Insufficient memory in CC_AAMAT. (2)') 194 195 CALL DZERO(WORK(KRES1EFF),NT1AM(ISYRES)) 196 CALL DZERO(WORK(KRES2EFF),NT2AM(ISYRES)) 197 198 IF (NODDY_AAMAT) THEN 199 CALL CCSDT_AAMAT_NODDY(IOPTRES,FREQ,LABEL,ISYOPE, 200 & LISTR,IDLSTR,.TRUE., 201 & WORK(KRES1),WORK(KRES2), 202 & WORK(KRES1EFF),WORK(KRES2EFF), 203 & IADOTS,ACONS,FILAAM,ITRAN, 204 & NATRAN,MXVEC,WORK(KEND1),LWRK1) 205 206 ELSE 207 208 WRITE(LUPRI,*)'AAMATSD should be called from ' 209 WRITE(LUPRI,*)'CC_BMAT module. Check NEW_RHS flag ' 210 WRITE(LUPRI,*)'in CCRHSVEC.' 211 CALL QUIT('Real code for AAMATSD not called here...') 212 213 END IF 214 215 216 END IF 217 218 END IF 219 220*---------------------------------------------------------------------* 221* write/add to vector on file: 222*---------------------------------------------------------------------* 223 IF (CCS) THEN 224 MODEL = 'CCS ' 225 IOPT = 1 226 ELSE IF (CC2) THEN 227 MODEL = 'CC2 ' 228 IOPT = 3 229 ELSE IF (CCSD) THEN 230 MODEL = 'CCSD ' 231 IOPT = 3 232 ELSE IF (CC3) THEN 233 MODEL = 'CC3 ' 234 IOPT = 3 235 IOPTE = 24 236 ELSE 237 CALL QUIT('Unknown coupled cluster model in CC_AAMAT.') 238 END IF 239 IF (CCR12) THEN 240 APROXR12 = ' ' 241 CALL CCSD_MODEL(MODEL,LENMOD,10,MODEL,10,APROXR12) 242 IOPTWR12 = 32 243 END IF 244 245 IF (IOPTRES.EQ.3) THEN 246 CALL CC_WRRSP(FILAAM, IFILE, ISYRES, IOPT, MODEL, 247 & DUMMY,WORK(KRES1),WORK(KRES2), 248 & WORK(KEND1),LWRK1) 249 IF (CCR12) THEN 250 CALL CC_WRRSP(FILAAM,IFILE,ISYRES,IOPTWR12, 251 & MODEL,DUMMY,DUMMY,WORK(KRESR12), 252 & WORK(KEND1),LWRK1) 253 END IF 254 IF (CCSDT) THEN 255 CALL CC_WRRSP(FILAAM,IFILE,ISYRES,IOPTE,MODEL,DUMMY, 256 & WORK(KRES1EFF),WORK(KRES2EFF),WORK(KEND1),LWRK1) 257 END IF 258 ELSE IF (IOPTRES.EQ.4) THEN 259 CALL CC_WARSP(FILAAM, IFILE, ISYRES, IOPT, MODEL, 260 & DUMMY,WORK(KRES1),WORK(KRES2), 261 & WORK(KEND1),LWRK1) 262 IF (CCR12) THEN 263 CALL CC_WARSP(FILAAM,IFILE,ISYRES,IOPTWR12, 264 & MODEL,DUMMY,DUMMY,WORK(KRESR12), 265 & WORK(KEND1),LWRK1) 266 END IF 267 IF (CCSDT) THEN 268 CALL CC_WARSP(FILAAM,IFILE,ISYRES,IOPTE,MODEL,DUMMY, 269 & WORK(KRES1EFF),WORK(KRES2EFF),WORK(KEND1),LWRK1) 270 END IF 271 ELSE IF (IOPTRES.EQ.5) THEN 272 IF (LOCDBG) THEN 273 IVEC = 1 274 WRITE(LUPRI,*) 'ACONS TRIPLES CONTRIBUTION:' 275 DO WHILE (IADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 276 WRITE (LUPRI,*) 277 & 'ACONS:',IVEC,ITRAN,ACONS(IVEC,ITRAN),IOPT 278 IVEC = IVEC + 1 279 END DO 280 END IF 281 282 CALL CCDOTRSP(IADOTS,ACONS,IOPT,FILAAM,ITRAN,NATRAN,MXVEC, 283 & WORK(KRES1),WORK(KRES2),ISYRES, 284 & WORK(KEND1),LWRK1) 285 IF (CCR12) THEN 286 CALL CCDOTRSP(IADOTS,ACONS,IOPTWR12,FILAAM,ITRAN,NATRAN, 287 & MXVEC,DUMMY,WORK(KRESR12),ISYRES, 288 & WORK(KEND1),LWRK1) 289 END IF 290 291 IF (LOCDBG) THEN 292 IVEC = 1 293 DO WHILE (IADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 294 WRITE (LUPRI,*) 295 & 'ACONS:',IVEC,ITRAN,ACONS(IVEC,ITRAN),IOPT 296 IVEC = IVEC + 1 297 END DO 298 END IF 299 ELSE 300 WRITE (LUPRI,*) 'Illegal option IOPTRES in CC_AAMAT.' 301 CALL QUIT('Illegal option IOPTRES in CC_AAMAT.') 302 END IF 303 304 IF (LOCDBG .AND. (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4)) THEN 305 WRITE (LUPRI,*) MSGDBG, 'wrote ',FILAAM,':',IFILE,' to disk.' 306 LEN = NT1AM(ISYRES) 307 IF (.NOT.CCS) LEN = LEN + NT2AM(ISYRES) 308 IF (CCR12) LEN = LEN + NTR12AM(ISYRES) 309 WRITE (LUPRI,*) MSGDBG, 'NORM^2 for result vector = ', 310 & DDOT(LEN,WORK(KRES1),1,WORK(KRES1),1) 311 WRITE (LUPRI,*) 'Listing of the result vector:' 312 CALL CC_PRP(WORK(KRES1),WORK(KRES2),ISYRES,1,1) 313 IF (CCR12) CALL CC_PRPR12(WORK(KRESR12),ISYRES,1,.TRUE.) 314 IF (CCSDT) THEN 315 WRITE (LUPRI,*) 'norm^2 of eff. result vector:', 316 & DDOT(LEN,WORK(KRES1EFF),1,WORK(KRES1EFF),1) 317 WRITE (LUPRI,*) 'Listing of the eff. result vector:' 318 CALL CC_PRP(WORK(KRES1EFF),WORK(KRES2EFF),ISYRES,1,1) 319 END IF 320 END IF 321 322*---------------------------------------------------------------------* 323* end of loop over A{O} transformations: 324*---------------------------------------------------------------------* 325 END DO 326 327 CALL QEXIT('CC_AAMAT') 328 RETURN 329 END 330*---------------------------------------------------------------------* 331* END OF SUBROUTINE CC_AAMAT * 332*---------------------------------------------------------------------* 333