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_symmback */ 20 SUBROUTINE SYMMBACK(MODEL, 21 & LUIAJK,FNDIAJK, 22 & LUAIJK,FNDAIJK, 23 & LUABI1,FNDABI1, 24 & LUABI3,FNDABI3, 25 & LUPTIAJB,FNDIAJB, 26 & ISYRES,WORK,LWORK) 27C 28C----------------------------------------------------------------------------- 29C Purpose: construct/resort the symmetrized (T) densities for CCSD(T) 30C resort the other (T) densities 31C and drive backtranformation of last index to AO 32C 33C S. Coriani, December 2001/January 2002 34C----------------------------------------------------------------------------- 35C 36 IMPLICIT NONE 37#include "priunit.h" 38#include "dummy.h" 39#include "iratdef.h" 40#include "ccsdsym.h" 41#include "inftap.h" 42#include "ccinftap.h" 43#include "ccorb.h" 44#include "ccsdinp.h" 45 46 INTEGER LUPTIAJB, LUAIJK, LUIAJK, LUABI1, LUABI3, LWORK 47 INTEGER LUIJKAS, LUIAJKS, LUAIJKS, LUABICS, LUAIBCS, LUIABCS 48 INTEGER LUIJKDS, LUIAJDS, LUAIJDS, LUABIDS 49 INTEGER LUIAJDE, LUAIJDE, LUAIBDE, LUIABDE 50 INTEGER ISYRES, KSTART 51 INTEGER KOCC1,KOCC2,KOCCR1,KOCCR2,KOCCS1,LWRK1,KEND1 52 INTEGER IOFF, ISYMA, ISYMI, ISYMJ, ISYMK 53 INTEGER ISYKJI,ISYMKJ 54 integer ISYMIJ, ISYIJK, KIJKA, KIJAK, KIJKAS 55 integer LENGTH, KAIB, KBIA, ISYM, KABIS 56 integer KDEN4, KDEN3 57 integer ISYBIA, ISYAIB, KAIB_C, KBIA_C, KABIS_C, ISYMB 58 integer ISYBIC, ISYCIB 59 integer KBICS_A 60 integer KOFF2, KOFF3, KOFF1, ISYMAI, ISYMC, ISYMBI 61 integer ISYMCI 62 integer IOPT, IOPTX 63 integer kaijk, kiajk, isyaij, isyden, isymab, isyabi 64 integer kcmo,kend0,lwrk0,koff 65 integer knoddy1,knoddy2,kaibd,kbida 66 integer kdaijb, kdaibj, kscrpk, isymbj, nai, naij, nbj 67 integer kaibj, kaijb, kiajb,kdiajb 68 integer kd1kjia, kd2jkia, kdiajk, kdaijk,isycmo 69 integer kbicr4,ioffbic_a,ioffcib_a 70 71#if defined (SYS_CRAY) 72 REAL WORK(LWORK) 73 REAL ZERO, ONE, TWO, FACT 74 REAL DDOT,XNORM2,XNORM,XNABI1,XNABI3 75#else 76 DOUBLE PRECISION WORK(LWORK) 77 DOUBLE PRECISION ZERO, ONE, TWO, FACT 78 DOUBLE PRECISION DDOT,XNORM2,XNORM,XNABI1,XNABI3 79#endif 80 81 82 CHARACTER MODEL*10 83 CHARACTER MODELPRI*4, MODELPRI2*14 84 CHARACTER*(*) FNDIAJK, FNDAIJK, FNDABI1, FNDABI3, FNDIAJB 85 LOGICAL LOCDBG 86C 87C Locally define density file names for easier debug 88C 89 CHARACTER FNDSIAJK*9, FNDSAIJK*9, FNDSIJKA*9 90 CHARACTER FNDSABIC*9, FNDSAIBC*9, FNDSIABC*9 91 CHARACTER FNDSIJKD*9, FNDSIAJD*9, FNDSAIJD*9, FNDSABID*9 92 CHARACTER FNDSAIBD*9, FNDSIABD*9 93 PARAMETER (FNDSIAJK = 'PT_DSIAJK', FNDSAIJK = 'PT_DSAIJK') 94 PARAMETER (FNDSIJKA = 'PT_DSIJKA', FNDSABIC = 'PT_DSABIC') 95 PARAMETER (FNDSAIBC = 'PT_DSAIBC', FNDSIABC = 'PT_DSIABC') 96c 97 PARAMETER (FNDSIJKD = 'PT_DSIJKD', FNDSIAJD = 'PT_DSIAJD') 98 PARAMETER (FNDSAIJD = 'PT_DSAIJD', FNDSABID = 'PT_DSABID') 99 PARAMETER (FNDSAIBD = 'PT_DSAIBD', FNDSIABD = 'PT_DSIABD') 100C 101 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 102 PARAMETER (LOCDBG = .FALSE.) 103C 104C#include "leinf.h" 105C 106 CALL QENTER('SYMMBACK') 107C 108C------------------------------- 109C First allocation 110C------------------------------- 111C 112 KCMO = 1 113 KEND0 = KCMO + NLAMDS 114 LWRK0 = LWORK - KEND0 115C 116 IF (LWRK0 .LT. 0) THEN 117 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND0 118 CALL QUIT( 119 * 'Insufficient memory for initial allocation') 120 ENDIF 121C 122C---------------------------------------------------------- 123C Read MO-coefficients from interface file and reorder. 124C---------------------------------------------------------- 125C 126 LUSIFC = -1 127 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 128 * IDUMMY,.FALSE.) 129 REWIND LUSIFC 130 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 131 READ (LUSIFC) 132 READ (LUSIFC) 133 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 134 CALL GPCLOSE (LUSIFC,'KEEP') 135C 136 CALL CMO_REORDER(WORK(KCMO),WORK(KEND0),LWRK0) 137C 138C--------------------------------------------------------- 139C Open files for backtransformed/symmetrized densities 140C--------------------------------------------------------- 141C 142 LUIJKAS = -1 143 LUIJKDS = -1 144 LUABICS = -1 145 LUABIDS = -1 146C 147C Temporarily keep the other "S" files with the resorted densities 148C (instead of symmetrized) 149C 150 LUIAJKS = -1 151 LUAIJKS = -1 152 LUAIBCS = -1 153 LUIABCS = -1 154C 155 LUIAJDE = -1 156 LUAIJDE = -1 157 LUAIBDE = -1 158 LUIABDE = -1 159C 160C Symmetrized and backtransformed 161C 162C d^s_{ijka} 163 CALL WOPEN2(LUIJKAS,FNDSIJKA,64,0) 164C d^s_{ijkdelta} 165 CALL WOPEN2(LUIJKDS,FNDSIJKD,64,0) 166C d^s_{abic} 167 CALL WOPEN2(LUABICS,FNDSABIC,64,0) 168C d^s_{abidelta} 169 CALL WOPEN2(LUABIDS,FNDSABID,64,0) 170C 171C Only backtransformed 172C 173C d^s_{iajdelta} 174 CALL WOPEN2(LUIAJDE,FNDSIAJD,64,0) 175C d^s_{aijdelta} 176 CALL WOPEN2(LUAIJDE,FNDSAIJD,64,0) 177C d^s_{aibdelta} 178 CALL WOPEN2(LUAIBDE,FNDSAIBD,64,0) 179C d^s_{iabdelta} 180 CALL WOPEN2(LUIABDE,FNDSIABD,64,0) 181C 182C Only symmetrized 183C 184C d^s_{aibc} 185 CALL WOPEN2(LUAIBCS,FNDSAIBC,64,0) 186C d^s_{iabc} 187 CALL WOPEN2(LUIABCS,FNDSIABC,64,0) 188C 189C--------------------------------------------------------------------- 190C - Read in iajk, aijk, iajb (the latter read-in inside backtransform) 191C - Resort iajk, aijk 192C - Backtransform last index and save on file 193C--------------------------------------------------------------------- 194 195 KSTART = KEND0 196 KOCC1 = KSTART 197 KOCC2 = KOCC1 + NCKIJ(ISYRES) 198 KOCCR1 = KOCC2 + NCKIJ(ISYRES) !resorted density 1 diajk 199 KOCCR2 = KOCCR1 + NCKIJ(ISYRES) !resorted density 2 daijk 200 KEND1 = KOCCR2 + NCKIJ(ISYRES) 201 LWRK1 = LWORK - KEND1 202C 203C----------- 204C Initialize 205C----------- 206C 207 CALL DZERO(WORK(KOCC1),NCKIJ(ISYRES)) 208 CALL DZERO(WORK(KOCC2),NCKIJ(ISYRES)) 209 CALL DZERO(WORK(KOCCR1),NCKIJ(ISYRES)) 210 CALL DZERO(WORK(KOCCR2),NCKIJ(ISYRES)) 211C 212C-------------------------------------------------------------- 213C Read in the entire diajk(kjia) and daijk(jkia) distributions, 214C and keep in WORK(KOCC1) and WORK(KOCC2) for later use 215C-------------------------------------------------------------- 216C 217 IOFF = 1 218 CALL GETWA2(LUIAJK,FNDIAJK,WORK(KOCC1),IOFF,NCKIJ(ISYRES)) 219C 220!OBS: scale -1 221C 222 CALL DSCAL(NCKIJ(ISYRES),-ONE,WORK(KOCC1),1) 223 224 CALL GETWA2(LUAIJK,FNDAIJK,WORK(KOCC2),IOFF,NCKIJ(ISYRES)) 225C 226!OBS: scale -1 227C 228 CALL DSCAL(NCKIJ(ISYRES),-ONE,WORK(KOCC2),1) 229 230C 231C------------------------------------ 232C Resort diajk(kjia) to diajk(ai,j,k) 233C Resort daijk(jkia) to daijk(ai,j,k) 234C------------------------------------ 235C 236 DO ISYMA = 1, NSYM 237 ISYKJI = MULD2H(ISYRES,ISYMA) 238 DO A = 1, NVIR(ISYMA) 239 DO ISYMI = 1, NSYM 240 ISYMKJ = MULD2H(ISYKJI,ISYMI) 241 DO I = 1, NRHF(ISYMI) 242 DO ISYMJ = 1, NSYM 243 ISYMAI = MULD2H(ISYMA,ISYMI) 244 ISYAIJ = MULD2H(ISYMAI,ISYMJ) 245 ISYMK = MULD2H(ISYMKJ,ISYMJ) 246 DO J = 1, NRHF(ISYMJ) 247 DO K = 1, NRHF(ISYMK) 248 KD1KJIA = I3OVIR(ISYKJI,ISYMA) 249 & + NMAIJK(ISYKJI)*(A-1) 250 & + IMAIJK(ISYMKJ,ISYMI) 251 & + NMATIJ(ISYMKJ)*(I-1) 252 & + IMATIJ(ISYMK,ISYMJ) 253 & + NRHF(ISYMK)*(J-1) 254 & + K 255 & + KOCC1 - 1 256 257 KD2JKIA = I3OVIR(ISYKJI,ISYMA) 258 & + NMAIJK(ISYKJI)*(A-1) 259 & + IMAIJK(ISYMKJ,ISYMI) 260 & + NMATIJ(ISYMKJ)*(I-1) 261 & + IMATIJ(ISYMJ,ISYMK) 262 & + NRHF(ISYMJ)*(K-1) 263 & + J 264 & + KOCC2 - 1 265 266 KDAIJK = ICKITR(ISYAIJ,ISYMK) 267 & + NCKI(ISYAIJ)*(K-1) 268 & + ICKI(ISYMAI,ISYMJ) 269 & + NT1AM(ISYMAI)*(J-1) 270 & + IT1AM(ISYMA,ISYMI) 271 & + NVIR(ISYMA)*(I-1) + A 272 273 KDIAJK = KDAIJK !d_iajk is stored aijk 274 275 WORK(KDIAJK + KOCCR1 - 1) = WORK(KD1KJIA) 276 WORK(KDAIJK + KOCCR2 - 1) = WORK(KD2JKIA) 277 278 END DO !K 279 END DO !J 280 END DO !ISYMJ 281 END DO !I 282 END DO !ISYMI 283 END DO !A 284 END DO !ISYMA 285C 286C-------------------------------------------------------- 287C Read-in d_iajb(ai,bj), unpack and reorder to ai,j,b 288C-------------------------------------------------------- 289C 290 KDIAJB = KEND1 291 KDAIBJ = KDIAJB + NT2SQ(ISYRES) 292 KSCRPK = KDAIBJ + NT2SQ(ISYRES) 293 KEND1 = KSCRPK + NT2AM(ISYRES) 294 LWRK1 = LWORK - KEND1 295 296 CALL DZERO(WORK(KDIAJB),NT2SQ(ISYRES)) 297 CALL DZERO(WORK(KDAIBJ),NT2SQ(ISYRES)) 298 299 IF (NT2AM(ISYRES) .GT. 0) THEN 300 IOFF = 1 301 CALL GETWA2(LUPTIAJB,FNDIAJB,WORK(KSCRPK), 302 & IOFF,NT2AM(ISYRES)) 303 !square up 304 CALL CC_T2SQ(WORK(KSCRPK),WORK(KDAIBJ),1) 305 !Husk: scale by 2* 306 CALL DSCAL(NT2SQ(ISYRES),TWO,WORK(KDAIBJ),1) 307 ENDIF 308 309 CALL CC3_T2TP(WORK(KDIAJB),WORK(KDAIBJ),ISYRES) 310C 311C----------------------------------------------- 312C Backtransform 313C----------------------------------------------- 314C 315 IOPTX=2 316 CALL BTRIAJDEL(IOPTX,WORK(KCMO),WORK(KOCCR1),WORK(KDIAJB), 317 & ISYRES,LUIAJDE,FNDSIAJD,WORK(KEND1),LWRK1) 318 319 IOPTX=1 320 CALL BTRIAJDEL(IOPTX,WORK(KCMO),WORK(KOCCR2),WORK(KDIAJB), 321 & ISYRES,LUAIJDE,FNDSAIJD,WORK(KEND1),LWRK1) 322C 323C------------------------------------------------------------------- 324C STEP A.2 325C Construct symmetrized occ.occ.occ.vir and backtransform last index 326C d^s_ijka (ijka) = d2_ijak(ijka) + d1_ijka(jika) 327C------------------------------------------------------------------- 328C 329 KOCCS1 = KOCCR1 330 KEND1 = KOCCS1 + NCKIJ(ISYRES) 331 LWRK1 = LWORK - KEND1 332 333 CALL DZERO(WORK(KOCCS1),NCKIJ(ISYRES)) 334 335 DO ISYMA = 1, NSYM 336 ISYIJK = MULD2H(ISYRES,ISYMA) 337 DO A = 1, NVIR(ISYMA) 338 DO ISYMK = 1, NSYM 339 ISYMIJ = MULD2H(ISYIJK,ISYMK) 340 DO K = 1, NRHF(ISYMK) 341 DO ISYMJ = 1, NSYM 342 ISYMI = MULD2H(ISYMIJ,ISYMJ) 343 DO J = 1, NRHF(ISYMJ) 344 DO I = 1, NRHF(ISYMI) 345 346 KIJAK = I3OVIR(ISYIJK,ISYMA) 347 & + NMAIJK(ISYIJK)*(A-1) 348 & + IMAIJK(ISYMIJ,ISYMK) 349 & + NMATIJ(ISYMIJ)*(K-1) 350 & + IMATIJ(ISYMI,ISYMJ) 351 & + NRHF(ISYMI)*(J-1) 352 & + I 353 & + KOCC2 - 1 354 355 356 KIJKA = I3OVIR(ISYIJK,ISYMA) 357 & + NMAIJK(ISYIJK)*(A-1) 358 & + IMAIJK(ISYMIJ,ISYMK) 359 & + NMATIJ(ISYMIJ)*(K-1) 360 & + IMATIJ(ISYMJ,ISYMI) 361 & + NRHF(ISYMJ)*(I-1) 362 & + J 363 & + KOCC1 - 1 364 365 KIJKAS = I3OVIR(ISYIJK,ISYMA) 366 & + NMAIJK(ISYIJK)*(A-1) 367 & + IMAIJK(ISYMIJ,ISYMK) 368 & + NMATIJ(ISYMIJ)*(K-1) 369 & + IMATIJ(ISYMI,ISYMJ) 370 & + NRHF(ISYMI)*(J-1) 371 & + I 372 & + KOCCS1 - 1 373 374 WORK(KIJKAS) = WORK(KIJAK) + WORK(KIJKA) 375 376 END DO !I 377 END DO !J 378 END DO !ISYMJ 379 END DO !K 380 END DO !ISYMK 381 END DO !A 382 END DO !ISYMA 383C 384C-------------- 385C Backtransform 386C-------------- 387C 388 if (locdbg) then 389 xnorm = ddot(NCKIJ(ISYRES),WORK(KOCC1),1,WORK(KOCC1),1) 390 write(lupri,*) 'Symmback: norm of d_ijka', xnorm 391 xnorm = ddot(NCKIJ(ISYRES),WORK(KOCC2),1,WORK(KOCC2),1) 392 write(lupri,*) 'Symmback: norm of d_ijak', xnorm 393 xnorm = ddot(NCKIJ(ISYRES),WORK(KOCCS1),1,WORK(KOCCS1),1) 394 write(lupri,*) 'Symmback: norm of ds_ijka', xnorm 395 end if 396 397 ISYCMO = 1 398 CALL BTRIJKDEL(WORK(KCMO),ISYCMO,WORK(KOCCS1),ISYRES, 399 & LUIJKDS,FNDSIJKD,WORK(KEND1),LWRK1) 400C 401C----------------------------------------------------------------- 402C STEP A.4 403C Construct symmetrized vir.vir.occ.vir density and backtransform 404C d^s_abic (abic) = d3_abic(aibc) + d4_abci(biac) 405C----------------------------------------------------------------- 406C 407 LENGTH = 0 408 DO ISYM = 1, NSYM 409 LENGTH = MAX(LENGTH,NCKATR(ISYM)) 410 END DO 411! 412 KAIB = KSTART !til den3 413 KBIA = KAIB + LENGTH !til den4 414 KABIS = KBIA + LENGTH 415 KEND1 = KABIS + LENGTH 416 LWRK1 = LWORK - KEND1 417! 418 CALL DZERO(WORK(KAIB),LENGTH) 419 CALL DZERO(WORK(KBIA),LENGTH) 420 CALL DZERO(WORK(KABIS),LENGTH) 421 422 xnorm = zero 423 xnabi1 = zero 424 xnabi3 = zero 425 426 DO ISYMC = 1, NSYM 427 ISYAIB = MULD2H(ISYRES,ISYMC) 428 ISYBIA = ISYAIB 429 ISYABI = ISYAIB 430 431 CALL DZERO(WORK(KAIB),NCKATR(ISYAIB)) 432 CALL DZERO(WORK(KBIA),NCKATR(ISYAIB)) 433 CALL DZERO(WORK(KABIS),NCKASR(ISYABI)) 434 435 DO C = 1, NVIR(ISYMC) 436 437 !determine the offset on file to the given block aib;c 438 !and bia;c 439 440 KAIB_C = ICKBD(ISYAIB,ISYMC) 441 & + NCKATR(ISYAIB)*(C-1) + 1 442 443 KBIA_C = ICKBD(ISYBIA,ISYMC) 444 & + NCKATR(ISYBIA)*(C-1) + 1 445 446C Read in aib and bia distributions for a given c. Put in Work restarting 447C from the beginning for each C (probabilmente poi per quella dopo devo 448C rileggere dal file) 449 450 CALL GETWA2(LUABI1,FNDABI1,WORK(KAIB),KAIB_C, 451 & NCKATR(ISYAIB)) 452!OBS: scale -1 453 CALL DSCAL(NCKATR(ISYAIB),-ONE,WORK(KAIB),1) 454 455 xnabi1 = xnabi1 + ddot(nckatr(isyaib),work(kaib), 456 & 1,work(kaib),1) 457C 458 CALL GETWA2(LUABI3,FNDABI3,WORK(KBIA),KBIA_C, 459 & NCKATR(ISYBIA)) 460 461 xnabi3 = xnabi3 + ddot(nckatr(isyaib),work(kbia), 462 & 1,work(kbia),1) 463 464 DO ISYMB = 1, NSYM 465 ISYMAI = MULD2H(ISYAIB,ISYMB) 466 DO B = 1, NVIR(ISYMB) 467 DO ISYMI = 1, NSYM 468 ISYMA = MULD2H(ISYMAI,ISYMI) 469 ISYMBI = MULD2H(ISYMB,ISYMI) 470 ISYMAB = MULD2H(ISYMA,ISYMB) 471 DO I = 1, NRHF(ISYMI) 472 DO A = 1, NVIR(ISYMA) 473 474 KOFF1 = KAIB - 1 475 & + ICKATR(ISYMAI,ISYMB) 476 & + NT1AM(ISYMAI)*(B-1) 477 & + IT1AM(ISYMA,ISYMI) 478 & + NVIR(ISYMA)*(I-1) + A 479 480 KOFF2 = KBIA - 1 481 & + ICKATR(ISYMBI,ISYMA) 482 & + NT1AM(ISYMBI)*(A-1) 483 & + IT1AM(ISYMB,ISYMI) 484 & + NVIR(ISYMB)*(I-1) + B 485 486 KOFF3 = KABIS - 1 487 & + ICKASR(ISYMAB,ISYMI) 488 & + NMATAB(ISYMAB)*(I-1) 489 & + IMATAB(ISYMA,ISYMB) 490 & + NVIR(ISYMA)*(B-1) + A 491C 492 WORK(KOFF3) = WORK(KOFF1) + WORK(KOFF2) 493 494 END DO !A 495 END DO !I 496 END DO !ISYMI 497 END DO !B 498 END DO 499 500 KABIS_C = ICDKVI(ISYABI,ISYMC) 501 & + NCKASR(ISYABI)*(C-1) + 1 502 xnorm = xnorm + 503 & ddot(NCKASR(ISYABI),WORK(KABIS),1,WORK(KABIS),1) 504 505 CALL PUTWA2(LUABICS,FNDSABIC,WORK(KABIS),KABIS_C, 506 & NCKASR(ISYABI)) 507C 508 ENDDO !C 509 ENDDO !ISYMC 510 511 if (locdbg) then 512 write(lupri,*) '--------------------------' 513 write(lupri,*) 'symmback: norm of ds_abic : ', xnorm 514 write(lupri,*) 'symmback: norm of d_abic : ', xnabi1 515 write(lupri,*) 'symmback: norm of d_abci : ', xnabi3 516 end if 517C 518C------------------------------- 519C Backtransform and dump on file 520C------------------------------- 521C 522 CALL BTRABIDEL(LUABIDS,FNDSABID,LUABICS,FNDSABIC,ISYRES, 523 & WORK(KCMO),WORK(KEND1),LWRK1) 524 525C----------------------------------------------------------------- 526C Resort 527C d3_iabc(bica) = d3_iabc(bica) (nothing to do) ABI1 528C d4_aibc(bica) = d4_aibc(ciba) ABI3 529C----------------------------------------------------------------- 530!1) resort d4 ciba->bica and put on file 531!2) resort d3 and d4 bica->biac and then biac->aibc 532!3) backtransform last index to delta 533 534! iabc 535 536 IOPT = 2 537 FACT = -ONE 538 CALL RESORT_DAIBC(IOPT,ISYRES,FACT,LUABI1,FNDABI1, 539 & LUIABCS,FNDSIABC,WORK(KEND1),LWRK1) 540C 541C------------------------------- 542C Backtransform and dump on file 543C------------------------------- 544C 545 CALL BTRAIBDEL(LUIABDE,FNDSIABD,LUIABCS,FNDSIABC, 546 & ISYRES,WORK(KCMO),WORK(KEND1),LWRK1) 547 548!aibc ------------------ 549 550 !first resort ciba-->bica 551 ! 552 IOPT = 1 553 FACT = ONE 554 CALL RESORT_DAIBC(IOPT,ISYRES,FACT,LUABI3,FNDABI3, 555 & LUABI3,FNDABI3,WORK(KEND1),LWRK1) 556 557 IOPT = 2 558 FACT = ONE 559 CALL RESORT_DAIBC(IOPT,ISYRES,FACT,LUABI3,FNDABI3, 560 & LUAIBCS,FNDSAIBC,WORK(KEND1),LWRK1) 561 562C 563C------------------------------- 564C Backtransform and dump on file 565C------------------------------- 566C 567 CALL BTRAIBDEL(LUAIBDE,FNDSAIBD,LUAIBCS,FNDSAIBC, 568 & ISYRES,WORK(KCMO),WORK(KEND1),LWRK1) 569C 570C----------------------------------------------------------------- 571C Close files and delete densities no longer required 572C----------------------------------------------------------------- 573C 574 CALL WCLOSE2(LUIAJDE,FNDSIAJD,'KEEP') 575 CALL WCLOSE2(LUAIJDE,FNDSAIJD,'KEEP') 576 CALL WCLOSE2(LUIJKDS,FNDSIJKD,'KEEP') 577 CALL WCLOSE2(LUABIDS,FNDSABID,'KEEP') 578 CALL WCLOSE2(LUAIBDE,FNDSAIBD,'KEEP') 579 CALL WCLOSE2(LUIABDE,FNDSIABD,'KEEP') 580 581 CALL WCLOSE2(LUABICS,FNDSABIC,'DELETE') 582 CALL WCLOSE2(LUAIBCS,FNDSAIBC,'DELETE') 583 CALL WCLOSE2(LUIABCS,FNDSIABC,'DELETE') 584 CALL WCLOSE2(LUIJKAS,FNDSIJKA,'DELETE') 585C 586C----------------------------------------------------------------- 587C Finished, return 588C----------------------------------------------------------------- 589C 590 CALL QEXIT('SYMMBACK') 591 RETURN 592 END 593C------------------------------------------------------------------ 594C /* Deck cc_symmback_1 */ 595 SUBROUTINE SYMMBACK_1(MODEL, 596 & LUIAJK,FNDIAJK, 597 & LUAIJK,FNDAIJK, 598 & LUABI1,FNDABI1, 599 & LUABI3,FNDABI3, 600 & LUPTIAJB,FNDIAJB, 601 & ISYRES,WORK,LWORK) 602C 603C----------------------------------------------------------------------------- 604C Purpose: construct/resort the symmetrized (T) densities for CCSD(T) 605C resort the other (T) densities 606C and drive backtranformation of last index to AO 607C 608C S. Coriani, December 2001/January 2002 609C Special Modified version for Orbit-Orbit corrections. TO be merged back!!! 610C----------------------------------------------------------------------------- 611C 612 IMPLICIT NONE 613#include "priunit.h" 614#include "dummy.h" 615#include "iratdef.h" 616#include "ccsdsym.h" 617#include "inftap.h" 618#include "ccinftap.h" 619#include "ccorb.h" 620#include "ccsdinp.h" 621Csonia 622#include "ccfop.h" 623 624 INTEGER LUPTIAJB, LUAIJK, LUIAJK, LUABI1, LUABI3, LWORK 625 INTEGER LUIJKAS, LUIAJKS, LUAIJKS, LUABICS, LUAIBCS, LUIABCS 626 INTEGER LUIJKDS, LUIAJDS, LUAIJDS, LUABIDS 627 INTEGER LUIAJDE, LUAIJDE, LUAIBDE, LUIABDE 628 INTEGER ISYRES, KSTART 629 INTEGER KOCC1,KOCC2,KOCCR1,KOCCR2,KOCCS1,LWRK1,KEND1 630 INTEGER IOFF, ISYMA, ISYMI, ISYMJ, ISYMK 631 INTEGER ISYKJI,ISYMKJ 632 integer iii 633 integer ISYMIJ, ISYIJK, KIJKA, KIJAK, KIJKAS 634 integer LENGTH, KAIB, KBIA, ISYM, KABIS 635 integer KDEN4, KDEN3 636 integer ISYBIA, ISYAIB, KAIB_C, KBIA_C, KABIS_C, ISYMB 637 integer ISYBIC, ISYCIB 638 integer KBICS_A 639 integer KOFF2, KOFF3, KOFF1, ISYMAI, ISYMC, ISYMBI 640 integer ISYMCI 641 integer IOPT,ioptx 642 integer kaijk, kiajk, isyaij, isyden, isymab, isyabi 643 integer kcmo,kend0,lwrk0,koff 644 integer knoddy1,knoddy2,kaibd,kbida 645 integer kdaijb, kdaibj, kscrpk, isymbj, nai, naij, nbj 646 integer kaibj, kaijb, kiajb,kdiajb 647 integer kd1kjia, kd2jkia, kdiajk, kdaijk,isycmo 648 integer kbicr4,ioffbic_a,ioffcib_a 649 650 integer LUABICA, LUABIDA,LUIJKDA, KOCCA1, KABIA 651 652#if defined (SYS_CRAY) 653 REAL WORK(LWORK) 654 REAL ZERO, ONE, TWO, FACT 655 real ddot,xnorm2,xnorm,xnabi1,xnabi3 656#else 657 DOUBLE PRECISION WORK(LWORK) 658 DOUBLE PRECISION ZERO, ONE, TWO, FACT 659 double precision ddot,xnorm2,xnorm,xnabi1,xnabi3 660#endif 661 662 663 CHARACTER MODEL*10 664 CHARACTER MODELPRI*4, MODELPRI2*14 665 CHARACTER*(*) FNDIAJK, FNDAIJK, FNDABI1, FNDABI3, FNDIAJB 666 LOGICAL LOCDBG 667C 668C Locally define density file names for easier debug 669C 670 CHARACTER FNDSIAJK*9, FNDSAIJK*9, FNDSIJKA*9 671 CHARACTER FNDSABIC*9, FNDSAIBC*9, FNDSIABC*9 672 CHARACTER FNDSIJKD*9, FNDSIAJD*9, FNDSAIJD*9, FNDSABID*9 673 CHARACTER FNDSAIBD*9, FNDSIABD*9 674c 675c Special extra antisymmetrized densities for BP2EOO 676c 677 CHARACTER FNDAABIC*9, FNDAIJKD*9, FNDAABID*9 678 PARAMETER (FNDAABIC = 'PT_DAABIC') 679 PARAMETER (FNDAABID = 'PT_DAABID',FNDAIJKD='PT_DAIJKD') 680c 681 PARAMETER (FNDSIAJK = 'PT_DSIAJK', FNDSAIJK = 'PT_DSAIJK') 682 PARAMETER (FNDSIJKA = 'PT_DSIJKA', FNDSABIC = 'PT_DSABIC') 683 PARAMETER (FNDSAIBC = 'PT_DSAIBC', FNDSIABC = 'PT_DSIABC') 684c 685 PARAMETER (FNDSIJKD = 'PT_DSIJKD', FNDSIAJD = 'PT_DSIAJD') 686 PARAMETER (FNDSAIJD = 'PT_DSAIJD', FNDSABID = 'PT_DSABID') 687 PARAMETER (FNDSAIBD = 'PT_DSAIBD', FNDSIABD = 'PT_DSIABD') 688C 689 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 690 PARAMETER (LOCDBG = .FALSE.) 691C 692C#include "leinf.h" 693C 694 CALL QENTER('SYMMBACK_1') 695C 696C------------------------------- 697C First allocation 698C------------------------------- 699C 700 KCMO = 1 701 KEND0 = KCMO + NLAMDS 702 LWRK0 = LWORK - KEND0 703C 704 IF (LWRK0 .LT. 0) THEN 705 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND0 706 CALL QUIT( 707 * 'Insufficient memory for initial allocation') 708 ENDIF 709C 710C---------------------------------------------------------- 711C Read MO-coefficients from interface file and reorder. 712C---------------------------------------------------------- 713C 714 LUSIFC = -1 715 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 716 * IDUMMY,.FALSE.) 717 REWIND LUSIFC 718 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 719 READ (LUSIFC) 720 READ (LUSIFC) 721 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 722 CALL GPCLOSE (LUSIFC,'KEEP') 723C 724 CALL CMO_REORDER(WORK(KCMO),WORK(KEND0),LWRK0) 725C 726C--------------------------------------------------------- 727C Open files for backtransformed/symmetrized densities 728C--------------------------------------------------------- 729C 730! LUIJKAS = -1 731 LUABICS = -1 732 LUIJKDS = -1 733 LUABIDS = -1 734C 735C Temporarily keep the other "S" files with the resorted densities 736C (instead of symmetrized) 737C 738! LUIAJKS = -1 739! LUAIJKS = -1 740 LUAIBCS = -1 741 LUIABCS = -1 742C 743 LUIAJDE = -1 744 LUAIJDE = -1 745 LUAIBDE = -1 746 LUIABDE = -1 747C 748C Symmetrized and backtransformed 749C 750C d^s_{ijka} 751! CALL WOPEN2(LUIJKAS,FNDSIJKA,64,0) 752C d^s_{ijkdelta} 753 CALL WOPEN2(LUIJKDS,FNDSIJKD,64,0) 754C d^s_{abic} 755 CALL WOPEN2(LUABICS,FNDSABIC,64,0) 756C d^s_{abidelta} 757 CALL WOPEN2(LUABIDS,FNDSABID,64,0) 758C 759C Only backtransformed 760C 761C d^s_{iajdelta} 762 CALL WOPEN2(LUIAJDE,FNDSIAJD,64,0) 763C d^s_{aijdelta} 764 CALL WOPEN2(LUAIJDE,FNDSAIJD,64,0) 765C d^s_{aibdelta} 766 CALL WOPEN2(LUAIBDE,FNDSAIBD,64,0) 767C d^s_{iabdelta} 768 CALL WOPEN2(LUIABDE,FNDSIABD,64,0) 769C 770C Only symmetrized 771C 772C d^s_{aibc} 773 CALL WOPEN2(LUAIBCS,FNDSAIBC,64,0) 774C d^s_{iabc} 775 CALL WOPEN2(LUIABCS,FNDSIABC,64,0) 776 777 if (BP2EOO) then 778 779! LUIJKAA = -1 780 LUABICA = -1 781 LUIJKDA = -1 782 LUABIDA = -1 783C 784C Antisymmetrized and backtransformed 785C 786C d^a_{ijka} 787! CALL WOPEN2(LUIJKAA,FNDAIJKA,64,0) 788C d^a_{ijkdelta} 789 CALL WOPEN2(LUIJKDA,FNDAIJKD,64,0) 790C d^a_{abic} 791 CALL WOPEN2(LUABICA,FNDAABIC,64,0) 792C d^a_{abidelta} 793 CALL WOPEN2(LUABIDA,FNDAABID,64,0) 794 795 end if 796C 797C--------------------------------------------------------------------- 798C - Read in iajk, aijk, iajb (the latter read-in inside backtransform) 799C - Resort iajk, aijk 800C - Backtransform last index and save on file 801C--------------------------------------------------------------------- 802 803 KSTART = KEND0 804 KOCC1 = KSTART 805 KOCC2 = KOCC1 + NCKIJ(ISYRES) 806 KOCCR1 = KOCC2 + NCKIJ(ISYRES) !resorted density 1 diajk 807 KOCCR2 = KOCCR1 + NCKIJ(ISYRES) !resorted density 2 daijk 808 KEND1 = KOCCR2 + NCKIJ(ISYRES) 809 LWRK1 = LWORK - KEND1 810C 811C----------- 812C Initialize 813C----------- 814C 815 CALL DZERO(WORK(KOCC1),NCKIJ(ISYRES)) 816 CALL DZERO(WORK(KOCC2),NCKIJ(ISYRES)) 817 CALL DZERO(WORK(KOCCR1),NCKIJ(ISYRES)) 818 CALL DZERO(WORK(KOCCR2),NCKIJ(ISYRES)) 819C 820C-------------------------------------------------------------- 821C Read in the entire diajk(kjia) and daijk(jkia) distributions, 822C and keep in WORK(KOCC1) and WORK(KOCC2) for later use 823C-------------------------------------------------------------- 824C 825 IOFF = 1 826 CALL GETWA2(LUIAJK,FNDIAJK,WORK(KOCC1),IOFF,NCKIJ(ISYRES)) 827C 828!OBS: scale -1 829C 830 CALL DSCAL(NCKIJ(ISYRES),-ONE,WORK(KOCC1),1) 831 832 CALL GETWA2(LUAIJK,FNDAIJK,WORK(KOCC2),IOFF,NCKIJ(ISYRES)) 833C 834!OBS: scale -1 835C 836 CALL DSCAL(NCKIJ(ISYRES),-ONE,WORK(KOCC2),1) 837 838C 839C------------------------------------ 840C Resort diajk(kjia) to diajk(ai,j,k) 841C Resort daijk(jkia) to daijk(ai,j,k) 842C------------------------------------ 843C 844 DO ISYMA = 1, NSYM 845 ISYKJI = MULD2H(ISYRES,ISYMA) 846 DO A = 1, NVIR(ISYMA) 847 DO ISYMI = 1, NSYM 848 ISYMKJ = MULD2H(ISYKJI,ISYMI) 849 DO I = 1, NRHF(ISYMI) 850 DO ISYMJ = 1, NSYM 851 ISYMAI = MULD2H(ISYMA,ISYMI) 852 ISYAIJ = MULD2H(ISYMAI,ISYMJ) 853 ISYMK = MULD2H(ISYMKJ,ISYMJ) 854 DO J = 1, NRHF(ISYMJ) 855 DO K = 1, NRHF(ISYMK) 856 KD1KJIA = I3OVIR(ISYKJI,ISYMA) 857 & + NMAIJK(ISYKJI)*(A-1) 858 & + IMAIJK(ISYMKJ,ISYMI) 859 & + NMATIJ(ISYMKJ)*(I-1) 860 & + IMATIJ(ISYMK,ISYMJ) 861 & + NRHF(ISYMK)*(J-1) 862 & + K 863 & + KOCC1 - 1 864 865 KD2JKIA = I3OVIR(ISYKJI,ISYMA) 866 & + NMAIJK(ISYKJI)*(A-1) 867 & + IMAIJK(ISYMKJ,ISYMI) 868 & + NMATIJ(ISYMKJ)*(I-1) 869 & + IMATIJ(ISYMJ,ISYMK) 870 & + NRHF(ISYMJ)*(K-1) 871 & + J 872 & + KOCC2 - 1 873 874 KDAIJK = ICKITR(ISYAIJ,ISYMK) 875 & + NCKI(ISYAIJ)*(K-1) 876 & + ICKI(ISYMAI,ISYMJ) 877 & + NT1AM(ISYMAI)*(J-1) 878 & + IT1AM(ISYMA,ISYMI) 879 & + NVIR(ISYMA)*(I-1) + A 880 881 KDIAJK = KDAIJK !d_iajk is stored aijk 882 883 WORK(KDIAJK + KOCCR1 - 1) = WORK(KD1KJIA) 884 WORK(KDAIJK + KOCCR2 - 1) = WORK(KD2JKIA) 885 886 END DO !K 887 END DO !J 888 END DO !ISYMJ 889 END DO !I 890 END DO !ISYMI 891 END DO !A 892 END DO !ISYMA 893 894C 895C Save on file 896C 897! IOFF = 1 898! CALL PUTWA2(LUIAJKS,FNDSIAJK,WORK(KOCCR1),IOFF,NCKIJ(ISYRES)) 899! CALL PUTWA2(LUAIJKS,FNDSAIJK,WORK(KOCCR2),IOFF,NCKIJ(ISYRES)) 900C 901C-------------------------------------------------------- 902C Read-in d_iajb(ai,bj), unpack and reorder to ai,j,b 903C-------------------------------------------------------- 904 905 KDIAJB = KEND1 906 KDAIBJ = KDIAJB + NT2SQ(ISYRES) 907 KSCRPK = KDAIBJ + NT2SQ(ISYRES) 908 KEND1 = KSCRPK + NT2AM(ISYRES) 909 LWRK1 = LWORK - KEND1 910 911 CALL DZERO(WORK(KDIAJB),NT2SQ(ISYRES)) 912 CALL DZERO(WORK(KDAIBJ),NT2SQ(ISYRES)) 913 914 IF (NT2AM(ISYRES) .GT. 0) THEN 915 IOFF = 1 916 CALL GETWA2(LUPTIAJB,FNDIAJB,WORK(KSCRPK), 917 & IOFF,NT2AM(ISYRES)) 918 !square up 919 CALL CC_T2SQ(WORK(KSCRPK),WORK(KDAIBJ),1) 920 !Husk: scale by 2* 921 CALL DSCAL(NT2SQ(ISYRES),TWO,WORK(KDAIBJ),1) 922 ENDIF 923 924 CALL CC3_T2TP(WORK(KDIAJB),WORK(KDAIBJ),ISYRES) 925C 926C----------------------------------------------- 927C Backtransform 928C----------------------------------------------- 929C 930 IOPTX=2 931 CALL BTRIAJDEL(IOPTX,WORK(KCMO),WORK(KOCCR1),WORK(KDIAJB), 932 & ISYRES,LUIAJDE,FNDSIAJD,WORK(KEND1),LWRK1) 933! -------------------------------------------------------- 934! write(lupri,*)'Symmback: call BTRAIJDEL (aijdel) ' 935! CALL BTRAIJDEL(WORK(KCMO),WORK(KOCCR2),ISYRES,LUAIJDE, 936! & FNDSAIJD,WORK(KEND1),LWRK1) 937! write(lupri,*) 'Symmback: exit BTRAIJDEL (aijdelta)' 938! call flshfo(lupri) 939! -------------------------------------------------------- 940 IOPTX=1 941 CALL BTRIAJDEL(IOPTX,WORK(KCMO),WORK(KOCCR2),WORK(KDIAJB), 942 & ISYRES,LUAIJDE,FNDSAIJD,WORK(KEND1),LWRK1) 943C 944C------------------------------------------------------------------- 945C STEP A.2 946C 947C Construct symmetrized occ.occ.occ.vir and backtransform last index 948C d^s_ijka (ijka) = d2_ijak(ijka) + d1_ijka(jika) 949C 950C If (BP2EOO) antisymmetrize as well 951C 952C------------------------------------------------------------------- 953C 954 KOCCS1 = KOCCR1 955 KEND1 = KOCCS1 + NCKIJ(ISYRES) 956 LWRK1 = LWORK - KEND1 957 958 CALL DZERO(WORK(KOCCS1),NCKIJ(ISYRES)) 959 960 if (BP2EOO) then 961 962 KOCCA1 = KEND1 963 KEND1 = KOCCA1 + NCKIJ(ISYRES) 964 LWRK1 = LWORK - KEND1 965 966 end if 967 968 DO ISYMA = 1, NSYM 969 ISYIJK = MULD2H(ISYRES,ISYMA) 970 DO A = 1, NVIR(ISYMA) 971 DO ISYMK = 1, NSYM 972 ISYMIJ = MULD2H(ISYIJK,ISYMK) 973 DO K = 1, NRHF(ISYMK) 974 DO ISYMJ = 1, NSYM 975 ISYMI = MULD2H(ISYMIJ,ISYMJ) 976 DO J = 1, NRHF(ISYMJ) 977 DO I = 1, NRHF(ISYMI) 978 979 KIJAK = I3OVIR(ISYIJK,ISYMA) 980 & + NMAIJK(ISYIJK)*(A-1) 981 & + IMAIJK(ISYMIJ,ISYMK) 982 & + NMATIJ(ISYMIJ)*(K-1) 983 & + IMATIJ(ISYMI,ISYMJ) 984 & + NRHF(ISYMI)*(J-1) 985 & + I 986 & + KOCC2 - 1 987 988 989 KIJKA = I3OVIR(ISYIJK,ISYMA) 990 & + NMAIJK(ISYIJK)*(A-1) 991 & + IMAIJK(ISYMIJ,ISYMK) 992 & + NMATIJ(ISYMIJ)*(K-1) 993 & + IMATIJ(ISYMJ,ISYMI) 994 & + NRHF(ISYMJ)*(I-1) 995 & + J 996 & + KOCC1 - 1 997 998 KIJKAS = I3OVIR(ISYIJK,ISYMA) 999 & + NMAIJK(ISYIJK)*(A-1) 1000 & + IMAIJK(ISYMIJ,ISYMK) 1001 & + NMATIJ(ISYMIJ)*(K-1) 1002 & + IMATIJ(ISYMI,ISYMJ) 1003 & + NRHF(ISYMI)*(J-1) 1004 & + I 1005 & + KOCCS1 - 1 1006 1007 1008 WORK(KIJKAS) = WORK(KIJAK) + WORK(KIJKA) 1009 1010 if (BP2EOO) then 1011 WORK(KIJKAS - KOCCS1 + KOCCA1) = 1012 & - WORK(KIJAK) + WORK(KIJKA) 1013 end if 1014 1015 END DO !I 1016 END DO !J 1017 END DO !ISYMJ 1018 END DO !K 1019 END DO !ISYMK 1020 END DO !A 1021 END DO !ISYMA 1022 1023C 1024C 1025! IOFF = 1 1026! CALL PUTWA2(LUIJKAS,FNDSIJKA,WORK(KOCCS1),IOFF,NCKIJ(ISYRES)) 1027! CALL WCLOSE2(LUIJKAS,FNDSIJKA,'KEEP') 1028C 1029C Backtransform 1030 1031 if (locdbg) then 1032 xnorm = ddot(NCKIJ(ISYRES),WORK(KOCC1),1,WORK(KOCC1),1) 1033 write(lupri,*) 'Symmback: norm of d_ijka', xnorm 1034 xnorm = ddot(NCKIJ(ISYRES),WORK(KOCC2),1,WORK(KOCC2),1) 1035 write(lupri,*) 'Symmback: norm of d_ijak', xnorm 1036 xnorm = ddot(NCKIJ(ISYRES),WORK(KOCCS1),1,WORK(KOCCS1),1) 1037 write(lupri,*) 'Symmback: norm of ds_ijka', xnorm 1038 end if 1039 1040 ISYCMO = 1 1041 CALL BTRIJKDEL(WORK(KCMO),ISYCMO,WORK(KOCCS1),ISYRES, 1042 & LUIJKDS,FNDSIJKD,WORK(KEND1),LWRK1) 1043 1044 if (BP2EOO) then 1045 ISYCMO = 1 1046 CALL BTRIJKDEL(WORK(KCMO),ISYCMO,WORK(KOCCA1),ISYRES, 1047 & LUIJKDA,FNDAIJKD,WORK(KEND1),LWRK1) 1048 end if 1049C 1050C----------------------------------------------------------------- 1051C STEP A.4 1052C Construct symmetrized vir.vir.occ.vir density and backtransform 1053C d^s_abic (abic) = d3_abic(aibc) + d4_abci(biac) 1054C----------------------------------------------------------------- 1055C 1056 LENGTH = 0 1057 DO ISYM = 1, NSYM 1058 LENGTH = MAX(LENGTH,NCKATR(ISYM)) 1059 END DO 1060! 1061 KAIB = KSTART !til den3 1062 KBIA = KAIB + LENGTH !til den4 1063 KABIS = KBIA + LENGTH 1064 KEND1 = KABIS + LENGTH 1065 LWRK1 = LWORK - KEND1 1066! 1067 CALL DZERO(WORK(KAIB),LENGTH) 1068 CALL DZERO(WORK(KBIA),LENGTH) 1069 CALL DZERO(WORK(KABIS),LENGTH) 1070 1071 1072 if (BP2EOO) then 1073 1074 KABIA = KABIS + LENGTH 1075 KEND1 = KABIA + LENGTH 1076 LWRK1 = LWORK - KEND1 1077 CALL DZERO(WORK(KABIA),LENGTH) 1078 1079 end if 1080 1081 xnorm = zero 1082 xnabi1 = zero 1083 xnabi3 = zero 1084 1085 DO ISYMC = 1, NSYM 1086 ISYAIB = MULD2H(ISYRES,ISYMC) 1087 ISYBIA = ISYAIB 1088 ISYABI = ISYAIB 1089 1090 CALL DZERO(WORK(KAIB),NCKATR(ISYAIB)) 1091 CALL DZERO(WORK(KBIA),NCKATR(ISYAIB)) 1092 CALL DZERO(WORK(KABIS),NCKASR(ISYABI)) 1093 1094 DO C = 1, NVIR(ISYMC) 1095 1096 !determine the offset on file to the given block aib;c 1097 !and bia;c 1098 1099 KAIB_C = ICKBD(ISYAIB,ISYMC) 1100 & + NCKATR(ISYAIB)*(C-1) + 1 1101 1102 KBIA_C = ICKBD(ISYBIA,ISYMC) 1103 & + NCKATR(ISYBIA)*(C-1) + 1 1104 1105C Read in aib and bia distributions for a given c. Put in Work restarting 1106C from the beginning for each C (probabilmente poi per quella dopo devo 1107C rileggere dal file) 1108 1109 CALL GETWA2(LUABI1,FNDABI1,WORK(KAIB),KAIB_C, 1110 & NCKATR(ISYAIB)) 1111!OBS: scale -1 1112 CALL DSCAL(NCKATR(ISYAIB),-ONE,WORK(KAIB),1) 1113 1114 xnabi1 = xnabi1 + ddot(nckatr(isyaib),work(kaib), 1115 & 1,work(kaib),1) 1116C 1117 CALL GETWA2(LUABI3,FNDABI3,WORK(KBIA),KBIA_C, 1118 & NCKATR(ISYBIA)) 1119 1120 xnabi3 = xnabi3 + ddot(nckatr(isyaib),work(kbia), 1121 & 1,work(kbia),1) 1122 1123 DO ISYMB = 1, NSYM 1124 ISYMAI = MULD2H(ISYAIB,ISYMB) 1125 DO B = 1, NVIR(ISYMB) 1126 DO ISYMI = 1, NSYM 1127 ISYMA = MULD2H(ISYMAI,ISYMI) 1128 ISYMBI = MULD2H(ISYMB,ISYMI) 1129 ISYMAB = MULD2H(ISYMA,ISYMB) 1130 DO I = 1, NRHF(ISYMI) 1131 DO A = 1, NVIR(ISYMA) 1132 1133 KOFF1 = KAIB - 1 1134 & + ICKATR(ISYMAI,ISYMB) 1135 & + NT1AM(ISYMAI)*(B-1) 1136 & + IT1AM(ISYMA,ISYMI) 1137 & + NVIR(ISYMA)*(I-1) + A 1138 1139 KOFF2 = KBIA - 1 1140 & + ICKATR(ISYMBI,ISYMA) 1141 & + NT1AM(ISYMBI)*(A-1) 1142 & + IT1AM(ISYMB,ISYMI) 1143 & + NVIR(ISYMB)*(I-1) + B 1144 1145 KOFF3 = KABIS - 1 1146 & + ICKASR(ISYMAB,ISYMI) 1147 & + NMATAB(ISYMAB)*(I-1) 1148 & + IMATAB(ISYMA,ISYMB) 1149 & + NVIR(ISYMA)*(B-1) + A 1150C 1151 WORK(KOFF3) = WORK(KOFF1) + WORK(KOFF2) 1152 1153 if (BP2EOO) then 1154 1155 WORK(KOFF3 - KABIS + KABIA) = 1156 & WORK(KOFF1) - WORK(KOFF2) 1157 1158 end if 1159 1160 END DO !A 1161 END DO !I 1162 END DO !ISYMI 1163 END DO !B 1164 END DO 1165 1166 KABIS_C = ICDKVI(ISYABI,ISYMC) 1167 & + NCKASR(ISYABI)*(C-1) + 1 1168 xnorm = xnorm + 1169 & ddot(NCKASR(ISYABI),WORK(KABIS),1,WORK(KABIS),1) 1170 1171 CALL PUTWA2(LUABICS,FNDSABIC,WORK(KABIS),KABIS_C, 1172 & NCKASR(ISYABI)) 1173 1174 if (BP2EOO) then 1175 CALL PUTWA2(LUABICA,FNDAABIC,WORK(KABIA),KABIS_C, 1176 & NCKASR(ISYABI)) 1177 end if 1178 1179C 1180 ENDDO !C 1181 !CALL DZERO(WORK(KABIS),NCKASR(ISYABI)) 1182 ENDDO !ISYMC 1183 if (locdbg) then 1184 write(lupri,*) '--------------------------' 1185 write(lupri,*) 'symmback: norm of ds_abic : ', xnorm 1186 write(lupri,*) 'symmback: norm of d_abic : ', xnabi1 1187 write(lupri,*) 'symmback: norm of d_abci : ', xnabi3 1188 end if 1189C 1190C------------------------------- 1191C Backtransform and dump on file 1192C------------------------------- 1193C 1194 CALL BTRABIDEL(LUABIDS,FNDSABID,LUABICS,FNDSABIC,ISYRES, 1195 & WORK(KCMO),WORK(KEND1),LWRK1) 1196 if (BP2EOO) then 1197 CALL BTRABIDEL(LUABIDA,FNDAABID,LUABICA,FNDAABIC,ISYRES, 1198 & WORK(KCMO),WORK(KEND1),LWRK1) 1199 end if 1200 1201C----------------------------------------------------------------- 1202C Resort 1203C d3_iabc(bica) = d3_iabc(bica) (nothing to do) ABI1 1204C d4_aibc(bica) = d4_aibc(ciba) ABI3 1205C----------------------------------------------------------------- 1206!1) resort d4 ciba->bica and put on file 1207!2) resort d3 and d4 bica->biac and then biac->aibc 1208!3) backtransform last index to delta 1209 1210! iabc 1211 1212 IOPT = 2 1213 FACT = -ONE 1214 CALL RESORT_DAIBC(IOPT,ISYRES,FACT,LUABI1,FNDABI1, 1215 & LUIABCS,FNDSIABC,WORK(KEND1),LWRK1) 1216C 1217C------------------------------- 1218C Backtransform and dump on file 1219C------------------------------- 1220C 1221 CALL BTRAIBDEL(LUIABDE,FNDSIABD,LUIABCS,FNDSIABC, 1222 & ISYRES,WORK(KCMO),WORK(KEND1),LWRK1) 1223 1224!aibc ------------------ 1225 1226 !first resort ciba-->bica 1227 ! 1228 IOPT = 1 1229 FACT = ONE 1230 CALL RESORT_DAIBC(IOPT,ISYRES,FACT,LUABI3,FNDABI3, 1231 & LUABI3,FNDABI3,WORK(KEND1),LWRK1) 1232 1233 IOPT = 2 1234 FACT = ONE 1235 CALL RESORT_DAIBC(IOPT,ISYRES,FACT,LUABI3,FNDABI3, 1236 & LUAIBCS,FNDSAIBC,WORK(KEND1),LWRK1) 1237 1238C 1239C------------------------------- 1240C Backtransform and dump on file 1241C------------------------------- 1242C 1243 CALL BTRAIBDEL(LUAIBDE,FNDSAIBD,LUAIBCS,FNDSAIBC, 1244 & ISYRES,WORK(KCMO),WORK(KEND1),LWRK1) 1245C 1246C----------------------------------------------------------------- 1247C Close files and delete densities no longer required 1248C----------------------------------------------------------------- 1249C 1250 CALL WCLOSE2(LUIAJDE,FNDSIAJD,'KEEP') 1251 CALL WCLOSE2(LUAIJDE,FNDSAIJD,'KEEP') 1252 CALL WCLOSE2(LUIJKDS,FNDSIJKD,'KEEP') 1253 CALL WCLOSE2(LUABIDS,FNDSABID,'KEEP') 1254 CALL WCLOSE2(LUAIBDE,FNDSAIBD,'KEEP') 1255 CALL WCLOSE2(LUIABDE,FNDSIABD,'KEEP') 1256 1257 IF (BP2EOO) THEN 1258 CALL WCLOSE2(LUIJKDA,FNDAIJKD,'KEEP') 1259 CALL WCLOSE2(LUABIDA,FNDAABID,'KEEP') 1260 END IF 1261 1262 CALL WCLOSE2(LUABICS,FNDSABIC,'DELETE') 1263 CALL WCLOSE2(LUAIBCS,FNDSAIBC,'DELETE') 1264 CALL WCLOSE2(LUIABCS,FNDSIABC,'DELETE') 1265! CALL WCLOSE2(LUIJKAS,FNDSIJKA,'DELETE') 1266C 1267C----------------------------------------------------------------- 1268C Finished, return 1269C----------------------------------------------------------------- 1270C 1271 CALL QEXIT('SYMMBACK_1') 1272 RETURN 1273 END 1274