1 subroutine hnd_sdfc(rtdb,geom,basis,vectors,nclosed,nopen,nvirt, 2 & nbf,nmo,pairlist,translate,ipairs,unique, 3 & i_pert,i_resp,tensor,coords,nfc,nsd,nsdfc) 4c $Id$ 5 implicit none 6#include "errquit.fh" 7#include "global.fh" 8#include "mafdecls.fh" 9#include "msgids.fh" 10#include "geom.fh" 11#include "rtdb.fh" 12#include "bas.fh" 13#include "nwc_const.fh" 14#include "stdio.fh" 15#include "apiP.fh" 16c 17 integer rtdb ! [input] rtdb handle 18 integer basis ! [input] basis handle 19 integer geom ! [input] geometry handle 20 integer vectors(2) ! [input] vectors 21 integer nclosed(2), nopen(2), nvirt(2) ! [input] occupation info 22 integer nbf, nmo ! [input] orbital info 23 integer ipairs ! [input] number of spin-spin pairs 24 integer pairlist(2*ipairs) ! [input] list of the pairs 25 integer translate(2*ipairs) ! [input] translation to unique list 26 integer i_pert, i_resp ! [input] # of unique responding and perturbing atoms 27 integer unique(i_pert+i_resp) ! [input] list of unique atoms 28 double precision tensor(3,3,5,ipairs) ! [output] spin-spin tensor, one for each term 29 double precision coords(3,i_pert+i_resp) ! [input] coordinates of unique atoms 30 double precision nfc, nsd, nsdfc ! [input] prefactors for the three terms 31c 32 integer ixy, ix, iy, iz, ifld 33 integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3) 34 integer dlo(3), dhi(3) 35 integer g_fca, g_fcb, g_sda, g_sdb, ii 36 integer g_d1(3),g_rhs,g_u(2) 37 integer i, j, i_total 38 double precision tol2e 39 character*256 cphf_rhs, cphf_sol 40c 41 logical cphf2, file_write_ga, file_read_ga, cphf 42 external cphf2, file_write_ga, file_read_ga, cphf 43c 44 logical oskel 45 double precision valuea, valueb 46 data tol2e /1.0d-16/ 47c 48 double precision pifac, froth 49c 50 integer ilist(3,3) 51 data ilist /1,4,5, 4,2,6, 5,6,3/ 52c 53 parameter(froth=4.0d0/3.0d0) 54 parameter(pifac=froth*3.14159265358979323846264338327950288419d0) 55c 56 oskel = .false. 57c 58c Integral initialization 59c 60 call int_init(rtdb,1,basis) 61 call schwarz_init(geom,basis) 62 call hnd_giao_init(basis,1) 63 call scf_get_fock_param(rtdb,tol2e) 64c 65c Total number of responses to be calculated 66c i_pert FC + 6 * i_pert SD + 6 * i_resp SD-FC 67c 68 i_total = 7*i_pert + 6*i_resp 69c 70c Create U matrix of dimension (nbf,nmo,3) and zero 71c Use ahi for dimension and ahi array for chunking/blocking 72c 73 alo(1) = nbf 74 alo(2) = -1 75 alo(3) = -1 76 ahi(1) = nbf 77 ahi(2) = nopen(1) 78 ahi(3) = i_total 79 if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u(1))) call 80 & errquit('hnd_sdfc: nga_create failed g_u',0,GA_ERR) 81 call ga_zero(g_u(1)) 82 if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u(2))) call 83 & errquit('hnd_sdfc: nga_create failed g_u',0,GA_ERR) 84 call ga_zero(g_u(2)) 85c 86c Construction of right-hand side CPHF 87c Create CPHF array of proper dimension : (nocc*nvirt,atom) 88c 89 if(.not.ga_create(MT_DBL,2*nopen(1)*nvirt(1),i_total,'RHS', 90 & -1,-1,g_rhs)) 91 & call errquit('hnd_sdfc: ga_create failed g_rhs',0,GA_ERR) 92 call ga_zero(g_rhs) 93c 94c Get FC and SD+FC in GA 95c 96 alo(1) = nbf 97 alo(2) = -1 98 alo(3) = -1 99 ahi(1) = nbf 100 ahi(2) = nbf 101 ahi(3) = i_pert+i_resp 102c 103c Allocate memory 104c 105 if (.not.nga_create(MT_DBL,3,ahi,'fca matrix',alo,g_fca)) call 106 & errquit('hnd_giaox: nga_create failed g_fca',0,GA_ERR) 107 if (.not.nga_create(MT_DBL,3,ahi,'fcb matrix',alo,g_fcb)) call 108 & errquit('hnd_giaox: nga_create failed g_fca',0,GA_ERR) 109 ahi(3) = 6*(i_pert+i_resp) 110 if (.not.nga_create(MT_DBL,3,ahi,'sda matrix',alo,g_sda)) call 111 & errquit('hnd_giaox: nga_create failed g_sda',0,GA_ERR) 112 if (.not.nga_create(MT_DBL,3,ahi,'sda matrix',alo,g_sdb)) call 113 & errquit('hnd_giaox: nga_create failed g_sda',0,GA_ERR) 114c 115c Calculate integrals 116c 117 call ga_zero(g_fca) 118 call ga_zero(g_sda) 119 call int_giao_1ega(basis,basis,g_fca,'fc',coords, 120 & i_pert+i_resp,oskel) 121 call int_giao_1ega(basis,basis,g_sda,'sd+fc',coords, 122 & i_pert+i_resp,oskel) 123c 124c Take out FC part from SD+FC (i.e. from xx, yy and zz components) 125c 126 alo(1) = 1 127 ahi(1) = nbf 128 alo(2) = 1 129 ahi(2) = nbf 130 blo(1) = 1 131 bhi(1) = nbf 132 blo(2) = 1 133 bhi(2) = nbf 134 do ixy = 1, i_pert+i_resp 135 alo(3) = 1+(ixy-1)*6 136 ahi(3) = 1+(ixy-1)*6 137 blo(3) = ixy 138 bhi(3) = ixy 139 do i = 1, 3 140 call nga_add_patch(1.0d0,g_sda,alo,ahi,pifac,g_fca,blo,bhi, 141 & g_sda,alo,ahi) 142 alo(3) = alo(3)+1 143 ahi(3) = ahi(3)+1 144 enddo 145 enddo 146c 147c NGA dimension arrays for copying will be the same every time 148c Also third NGA dimension for any of the three dimensional 149c arrays will be the same everytime (running from 1 to 3) 150c So, lets define them once and for all in blo and bhi 151c 152 blo(1) = 1 153 bhi(1) = nopen(1)*nvirt(1) 154 blo(2) = 1 155 bhi(2) = i_pert 156c 157c Now that all the work is done on the integrals in Alpha, copy to Beta 158c 159 call ga_copy(g_fca,g_fcb) 160 call ga_copy(g_sda,g_sdb) 161c 162c ga_rhs(a,i) = ga_rhs(a,i) + FC(a,i) 163c Transform FC and SD to MO and add to g_rhs 164c 165 call giao_aotomo(g_fca,vectors(1),nopen(1),nvirt(1),1,i_pert,nbf) 166 call giao_aotomo(g_fcb,vectors(2),nopen(2),nvirt(2),1,i_pert,nbf) 167 call giao_aotomo(g_sda,vectors(1),nopen(1),nvirt(1),1, 168 & 6*(i_pert+i_resp),nbf) 169 call giao_aotomo(g_sdb,vectors(2),nopen(2),nvirt(2),1, 170 & 6*(i_pert+i_resp),nbf) 171c 172c Add to g_rhs 173c 174 alo(1) = nopen(1)+1 175 ahi(1) = nmo 176 alo(2) = 1 177 ahi(2) = nopen(1) 178 alo(3) = 1 179 ahi(3) = i_pert 180 call nga_add_patch(1.0d0,g_rhs,blo,bhi,0.5d0,g_fca,alo,ahi, 181 & g_rhs,blo,bhi) 182 blo(1) = blo(1)+nopen(1)*nvirt(1) 183 bhi(1) = bhi(1)+nopen(2)*nvirt(2) 184 call nga_add_patch(1.0d0,g_rhs,blo,bhi,-0.5d0,g_fcb,alo,ahi, 185 & g_rhs,blo,bhi) 186c 187 blo(1) = 1 188 bhi(1) = nopen(1)*nvirt(1) 189 blo(2) = i_pert+1 190 bhi(2) = i_pert+6*(i_pert+i_resp) 191 ahi(3) = 6*(i_pert+i_resp) 192 call nga_add_patch(1.0d0,g_rhs,blo,bhi,0.5d0,g_sda,alo,ahi, 193 & g_rhs,blo,bhi) 194 blo(1) = blo(1)+nopen(1)*nvirt(1) 195 bhi(1) = bhi(1)+nopen(2)*nvirt(2) 196 call nga_add_patch(1.0d0,g_rhs,blo,bhi,-0.5d0,g_sdb,alo,ahi, 197 & g_rhs,blo,bhi) 198 call ga_scale(g_rhs,-4.0d0) 199c 200c Some memory cleanup here 201c 202 if (.not.ga_destroy(g_fcb)) call 203 & errquit('hnd_sdfc: ga_destroy failed g_fcb',0,GA_ERR) 204 if (.not.ga_destroy(g_sdb)) call 205 & errquit('hnd_sdfc: ga_destroy failed g_sdb',0,GA_ERR) 206c 207c Write ga_rhs to disk 208c 209 call cphf_fname('cphf_rhs',cphf_rhs) 210 call cphf_fname('cphf_sol',cphf_sol) 211 if(.not.file_write_ga(cphf_rhs,g_rhs)) call errquit 212 $ ('hnd_sdfc: could not write cphf_rhs',0, DISK_ERR) 213c 214 call schwarz_tidy() 215 call int_terminate() 216c 217c Call the CPHF routine 218c 219 if (.not.cphf2(rtdb)) call errquit 220 $ ('hnd_sdfc: failure in cphf ',0, RTDB_ERR) 221c 222c Occ-virt blocks are the solution pieces of the CPHF 223c Read solution vector from disk and put solutions in U matrices 224c 225 call int_init(rtdb,1,basis) 226 call schwarz_init(geom,basis) 227 call ga_zero(g_rhs) 228 if(.not.file_read_ga(cphf_sol,g_rhs)) call errquit 229 $ ('hnd_sdfc: could not read cphf_rhs',0, DISK_ERR) 230 alo(1) = nopen(1)+1 231 ahi(1) = nmo 232 alo(2) = 1 233 ahi(2) = nopen(1) 234 alo(3) = 1 235 ahi(3) = i_total 236 blo(1) = 1 237 bhi(1) = nopen(1)*nvirt(1) 238 blo(2) = 1 239 bhi(2) = i_total 240 call nga_copy_patch('n',g_rhs,blo,bhi,g_u(1),alo,ahi) 241 blo(1) = blo(1)+nopen(1)*nvirt(1) 242 bhi(1) = bhi(1)+nopen(2)*nvirt(2) 243 call nga_copy_patch('n',g_rhs,blo,bhi,g_u(2),alo,ahi) 244c 245 if (.not.ga_destroy(g_rhs)) call 246 & errquit('hnd_sdfc: ga_destroy failed g_rhs',0,GA_ERR) 247c 248c From U matrices, generate the perturbed density matrices D1 249c C1 = C0 * U10 250c D1 = [(C1*C0+) + (C0*C1+)] 251c 252 alo(1) = nbf 253 alo(2) = -1 254 alo(3) = -1 255 ahi(1) = nbf 256 ahi(2) = nbf 257 ahi(3) = i_total 258 if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1(1))) call 259 & errquit('hnd_sdfc: nga_create failed g_d1(1)',0,GA_ERR) 260 if (.not.nga_create(MT_DBL,3,ahi,'D10 matrix',alo,g_d1(2))) call 261 & errquit('hnd_sdfc: nga_create failed g_d1(2)',0,GA_ERR) 262c 263 alo(1) = 1 264 alo(2) = 1 265 blo(1) = 1 266 blo(2) = 1 267 clo(1) = 1 268 chi(1) = nbf 269 clo(2) = 1 270 chi(2) = nbf 271 dlo(1) = 1 272 dhi(1) = nbf 273 dlo(2) = 1 274 do ii = 1, 2 275 dhi(2) = nopen(ii) 276 do ifld = 1, i_total 277 alo(3) = ifld 278 ahi(3) = ifld 279 blo(3) = ifld 280 bhi(3) = ifld 281 clo(3) = ifld 282 chi(3) = ifld 283 dlo(3) = ifld 284 dhi(3) = ifld 285 bhi(1) = nbf 286 bhi(2) = nmo 287 ahi(1) = nmo 288 ahi(2) = nopen(ii) 289c 290c Make C1 291c 292 call nga_matmul_patch('n','n',1.0d0,0.0d0,vectors(ii),blo,bhi, 293 & g_u(ii),alo,ahi,g_d1(ii),dlo,dhi) 294 call nga_copy_patch('n',g_d1(ii),dlo,dhi,g_u(ii),dlo,dhi) 295 bhi(1) = nbf 296 bhi(2) = nopen(ii) 297 ahi(1) = nopen(ii) 298 ahi(2) = nbf 299c 300c Make D1 301c 302 call nga_matmul_patch('n','t',1.0d0,0.0d0,vectors(ii),blo,bhi, 303 & g_u(ii),alo,ahi,g_d1(ii),clo,chi) 304 call nga_matmul_patch('n','t',1.0d0,0.0d0,g_u(ii),blo,bhi, 305 & vectors(ii),alo,ahi,g_d1(ii),clo,chi) 306 enddo 307 enddo 308c 309 call ga_sync() 310c 311 if (.not.ga_destroy(g_u(1))) call 312 & errquit('hnd_sdfc: ga_destroy failed g_u(1)',0,GA_ERR) 313 if (.not.ga_destroy(g_u(2))) call 314 & errquit('hnd_sdfc: ga_destroy failed g_u(2)',0,GA_ERR) 315c 316c Now we have in g_d1(nmo,nmo,3) the derivative densities and 317c hence we can calculate the contributions to the FC term 318c 319 call ga_zero(g_fca) 320 call int_giao_1ega(basis,basis,g_fca,'fc',coords(1,i_pert+1), 321 & i_resp,oskel) 322 alo(1) = 1 323 ahi(1) = nbf 324 alo(2) = 1 325 ahi(2) = nbf 326 blo(1) = 1 327 bhi(1) = nbf 328 blo(2) = 1 329 bhi(2) = nbf 330 do i = 1, ipairs 331c alo(3) = unique(translate(i)) 332c ahi(3) = unique(translate(i)) 333c blo(3) = unique(translate(i+ipairs)) 334c bhi(3) = unique(translate(i+ipairs)) 335 alo(3) = translate(i) 336 ahi(3) = translate(i) 337 blo(3) = translate(i+ipairs) 338 bhi(3) = translate(i+ipairs) 339 valuea=nga_ddot_patch(g_d1(1),'n',alo,ahi,g_fca,'n',blo,bhi) 340 valueb=nga_ddot_patch(g_d1(2),'n',alo,ahi,g_fca,'n',blo,bhi) 341c if (ga_nodeid().eq.0) 342c & write(6,'(A,2I4,4F12.6)') 'FC tens i j com', 343c & pairlist(ixy),pairlist(ixy+ipairs), 344c & tensor(1,1,1,ixy),(valuea-valueb)*0.5d0*nfc, 345c & valuea,valueb 346 do ixy = 1, 3 347 tensor(ixy,ixy,1,i)=(valuea-valueb)*nfc*0.5d0 348 enddo 349 enddo 350c 351c Calculate the SD-SD tensor. First calculate the sd+fc integrals 352c and subtract the fc contribution 353c 354 call ga_zero(g_sda) 355 call int_giao_1ega(basis,basis,g_sda,'sd+fc',coords(1,i_pert+1), 356 & i_resp,oskel) 357 alo(1) = 1 358 alo(2) = 1 359 blo(1) = 1 360 bhi(1) = nbf 361 blo(2) = 1 362 bhi(2) = nbf 363 clo(1) = 1 364 chi(1) = nbf 365 clo(2) = 1 366 chi(2) = nbf 367 do i = 1, i_resp 368 alo(3) = 1+(i-1)*6 369 ahi(3) = 1+(i-1)*6 370 blo(3) = i 371 bhi(3) = i 372 do ixy = 1, 3 373 call nga_add_patch(1.0d0,g_sda,alo,ahi,pifac,g_fca,blo,bhi, 374 & g_sda,alo,ahi) 375 alo(3) = alo(3)+1 376 ahi(3) = ahi(3)+1 377 enddo 378 enddo 379c 380c Calculate the tensor contributions to the SD-SD term 381c 382 do ixy = 1, ipairs 383c i = unique(translate(ixy)) 384c j = unique(translate(ixy+ipairs)) 385 i = translate(ixy) 386 j = translate(ixy+ipairs) 387c RAW contributions for reference 388c do ix = 1, 3 389c do iy = 1, 3 390c alo(3) = i_pert + (i-1)*6 + ilist(ix,iy) 391c ahi(3) = i_pert + (i-1)*6 + ilist(ix,iy) 392c blo(3) = (j-1)*6 + ilist(ix,iy) 393c bhi(3) = (j-1)*6 + ilist(ix,iy) 394c valuea = nga_ddot_patch(g_sda,'n',blo,bhi, 395c & g_d1(1),'n',alo,ahi) 396c valueb = nga_ddot_patch(g_sda,'n',blo,bhi, 397c & g_d1(2),'n',alo,ahi) 398c if (ga_nodeid().eq.0) 399c & write(6,'(A,4I4,5F12.6)') 'SD raw i j com', 400c & pairlist(ixy),pairlist(ixy+ipairs),ix,iy, 401c & (valuea-valueb)*nsd*0.5d0, 402c & valuea,valueb,nsd 403c enddo 404c enddo 405 do ix = 1, 3 406 do iy = 1, 3 407 valuea = 0.0d0 408 valueb = 0.0d0 409 do iz = 1, 3 410 alo(3) = i_pert + (i-1)*6 + ilist(iz,iy) 411 ahi(3) = i_pert + (i-1)*6 + ilist(iz,iy) 412 blo(3) = (j-1)*6 + ilist(iz,ix) 413 bhi(3) = (j-1)*6 + ilist(iz,ix) 414 valuea = valuea + nga_ddot_patch(g_sda,'n',blo,bhi, 415 & g_d1(1),'n',alo,ahi) 416 valueb = valueb + nga_ddot_patch(g_sda,'n',blo,bhi, 417 & g_d1(2),'n',alo,ahi) 418 enddo 419 tensor(iy,ix,2,ixy)=(valuea-valueb)*nsd*0.5d0 420 enddo 421 enddo 422 enddo 423c 424c The final term to be calculated is the SD-FC tensorterm , which consists 425c of two contributions: 426c D1(SD)_N * FC_M + D1(SD)_M * FC_N N=perturbing, M=responding 427c We have all the contributions for the first term in the matrices, 428c hence lets do the ddots 429c 430 do ixy = 1, ipairs 431c i = unique(translate(ixy)) 432c j = unique(translate(ixy+ipairs)) 433 i = translate(ixy) 434 j = translate(ixy+ipairs) 435 do ix = 1, 3 436 do iy = 1, 3 437 alo(3) = i_pert + (i-1)*6 + ilist(ix,iy) 438 ahi(3) = i_pert + (i-1)*6 + ilist(ix,iy) 439 blo(3) = j 440 bhi(3) = j 441 valuea = nga_ddot_patch(g_fca,'n',blo,bhi, 442 & g_d1(1),'n',alo,ahi) 443 valueb = nga_ddot_patch(g_fca,'n',blo,bhi, 444 & g_d1(2),'n',alo,ahi) 445 tensor(ix,iy,3,ixy)=(valuea-valueb)*0.5d0 446c if (ga_nodeid().eq.0) 447c & write(6,'(A,4I4,3F12.6)') 'SDFC NM tens i j com', 448c & pairlist(ixy),pairlist(ixy+ipairs), 449c & ix,iy,tensor(ix,iy,3,ixy)*nsdfc,valuea,valueb 450 enddo 451 enddo 452c if (ga_nodeid().eq.0) print*,'' 453 enddo 454c 455c For the D1(SD)_M * FC_N N=perturbing, M=responding we need to calculate the 456c fc integrals for the perturbing atoms. We already calculated the D1(SD) for 457c the responding atoms as extra terms. 458c 459 call ga_zero(g_fca) 460 call int_giao_1ega(basis,basis,g_fca,'fc',coords, 461 & i_pert,oskel) 462 do ixy = 1, ipairs 463c i = unique(translate(ixy)) 464c j = unique(translate(ixy+ipairs)) 465 i = translate(ixy) 466 j = translate(ixy+ipairs) 467 do ix = 1, 3 468 do iy = 1, 3 469 valuea = 0.0d0 470 valueb = 0.0d0 471 alo(3) = 7*i_pert + (j-1)*6 + ilist(ix,iy) 472 ahi(3) = 7*i_pert + (j-1)*6 + ilist(ix,iy) 473 blo(3) = i 474 bhi(3) = i 475 valuea = valuea + nga_ddot_patch(g_fca,'n',blo,bhi, 476 & g_d1(1),'n',alo,ahi) 477 valueb = valueb + nga_ddot_patch(g_fca,'n',blo,bhi, 478 & g_d1(2),'n',alo,ahi) 479 tensor(ix,iy,3,ixy)=tensor(ix,iy,3,ixy)+ 480 & (valuea-valueb)*0.5d0 481 tensor(ix,iy,3,ixy)=tensor(ix,iy,3,ixy)*nsdfc 482c if (ga_nodeid().eq.0) 483c & write(6,'(A,4I4,4F12.6)') 'SDFC MN tens i j com', 484c & pairlist(ixy),pairlist(ixy+ipairs), 485c & ix,iy,tensor(ix,iy,3,ixy),(valuea-valueb)*0.5d0*nsdfc, 486c & valuea,valueb 487 enddo 488 enddo 489c if (ga_nodeid().eq.0) print*,'' 490 enddo 491c 492 call ga_sync() 493c 494c Clean up all remaining memory 495c 496 if (.not.ga_destroy(g_d1(1))) call 497 & errquit('hnd_sdfc: ga_destroy failed g_d1(1)',0,GA_ERR) 498 if (.not.ga_destroy(g_d1(2))) call 499 & errquit('hnd_sdfc: ga_destroy failed g_d1(2)',0,GA_ERR) 500 if (.not.ga_destroy(g_sda)) call 501 & errquit('hnd_sdfc: ga_destroy failed g_sda',0,GA_ERR) 502 if (.not.ga_destroy(g_fca)) call 503 & errquit('hnd_sdfc: ga_destroy failed g_fca',0,GA_ERR) 504c 505 call schwarz_tidy() 506 call int_terminate() 507c 508 return 509 end 510