1c---------------- get_chi_centers_ga() ------------- START 2 subroutine get_chi_centers(chi_cntr, ! out 3 & basis, ! in : basis handle 4 & nbf, ! in : nr basis functions 5 & geom, ! in : geometry handle 6 & mcenters) ! in : nr. atoms 7c 8c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 9 10 implicit none 11 12#include "rtdb.fh" 13#include "nwc_const.fh" 14#include "errquit.fh" 15#include "global.fh" 16#include "bas.fh" 17#include "geom.fh" 18#include "mafdecls.fh" 19#include "stdio.fh" 20#include "util.fh" 21#include "msgids.fh" 22 integer basis,geom,nbf 23 double precision chi_cntr(3,nbf) ! OUTPUT 24 double precision cnt(3),valZ 25 integer ictr,ic1,ic2,icset 26 integer l,nprim,ncontr,isphere,nshbf 27 integer mcenters,i,n1 28 integer iniz,ifin 29 character*16 at_tag 30 integer lo1(3),hi1(3),ld(2) 31 logical status 32 n1=0 33 do ictr=1,mcenters 34 iniz=0 35 ifin=0 36 if (.not.bas_ce2cnr(basis,ictr,ic1,ic2)) 37 & call errquit('Exiting in get_chi_centers_ga.', 38 & 11, BASIS_ERR) 39 do icset = ic1,ic2 40c ----- get info about current contraction set 41 if (.not. bas_continfo(basis,icset,l,nprim, 42 & ncontr,isphere)) 43 & call errquit('Exiting in get_chi_centers_ga.', 44 & 5, BASIS_ERR) 45 nshbf=ncontr*(((l+1)*(l+2))/2) 46 if(isphere.eq.1) then 47 nshbf=ncontr*(2*l+1) 48 endif 49 if (iniz.eq.0) iniz=n1+1 50 n1=n1+nshbf 51 enddo ! end loop icset 52 ifin= n1 53 status=geom_cent_get(geom,ictr,at_tag, 54 & cnt,valZ) 55c if iniz=0, must be a bq or bare ecp 56 if(iniz.ne.0) then 57 do i=iniz,ifin 58 chi_cntr(1,i)=cnt(1) 59 chi_cntr(2,i)=cnt(2) 60 chi_cntr(3,i)=cnt(3) 61 enddo ! end loop i 62 endif 63 enddo ! end loop ictr 64 return 65 end 66 subroutine get_chi_centers_ga(g_chi_cntr, ! out 67 & basis, ! in : basis handle 68 & nbf,geom,mcenters) 69c 70c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 71 72 implicit none 73 74#include "rtdb.fh" 75#include "nwc_const.fh" 76#include "errquit.fh" 77#include "global.fh" 78#include "bas.fh" 79#include "geom.fh" 80#include "mafdecls.fh" 81#include "stdio.fh" 82#include "util.fh" 83#include "msgids.fh" 84 integer basis,geom,nbf 85 integer g_chi_cntr(3) ! OUTPUT 86 double precision cnt(3),valZ 87 integer ictr,ic1,ic2,icset 88 integer l,nprim,ncontr,isphere,nshbf 89 integer mcenters,i,n1 90 integer iniz,ifin 91 character*16 at_tag 92 integer l_buf,k_buf 93 integer lo1(3),hi1(3),ld(2) 94 logical status 95c ----- allocate array to store centers ---- START 96 if(.not.MA_push_get(MT_DBL,3*nbf,'get_chi_centers_ga:buf', 97 & l_buf,k_buf)) 98 $ call errquit('get_chi_centers_ga: ma failed', 99 & 3*nbf, MA_ERR) 100c ----- allocate array to store centers ---- END 101 n1=0 102 do ictr=1,mcenters 103 iniz=0 104 ifin=0 105 if (.not.bas_ce2cnr(basis,ictr,ic1,ic2)) 106 & call errquit('Exiting in get_chi_centers_ga.', 107 & 11, BASIS_ERR) 108 do icset = ic1,ic2 109c ----- get info about current contraction set 110 if (.not. bas_continfo(basis,icset,l,nprim, 111 & ncontr,isphere)) 112 & call errquit('Exiting in get_chi_centers_ga.', 113 & 5, BASIS_ERR) 114 nshbf=ncontr*(((l+1)*(l+2))/2) 115 if(isphere.eq.1) then 116 nshbf=ncontr*(2*l+1) 117 endif 118 if (iniz.eq.0) iniz=n1+1 119 n1=n1+nshbf 120 enddo ! end loop icset 121 ifin= n1 122 status=geom_cent_get(geom,ictr,at_tag, 123 & cnt,valZ) 124c if iniz=0, must be a bq or bare ecp 125 if(iniz.ne.0) then 126 do i=iniz,ifin 127 dbl_mb(k_buf +i-1)=cnt(1) 128 dbl_mb(k_buf+nbf +i-1)=cnt(2) 129 dbl_mb(k_buf+2*nbf+i-1)=cnt(3) 130 enddo ! end loop i 131 endif 132 enddo ! end loop ictr 133c ----- store in g_chi_cntr() --- START 134c dbl_mb() ---> g_chi_cntr() 135 do i=1,3 136 ld(1)=nbf 137 lo1(1)=1 138 hi1(1)=nbf 139 lo1(2)=i 140 hi1(3)=i 141 call nga_put(g_chi_cntr(i), 142 & lo1,hi1,dbl_mb(k_buf+(i-1)*nbf),ld) 143 enddo ! end-loop-i 144c ----- store in g_chi_cntr() --- END 145c --- Free memory 146 if (.not. MA_pop_stack(l_buf)) call errquit 147 $ ('get_chi_centers_ga: pop failed', 0, GA_ERR) 148 return 149 end 150c---------------- get_chi_centers_ga() ------------- END 151 subroutine get_3rdterm_R(g_N, ! to be scaled 152 & g_R, ! scaling 153 & ind_a, ! from kab=123,231,312 154 & ind_b, ! from kab=123,231,312 155 & g_tmp2, ! scratch 156 & g_N_scld)! output 157c 158c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 159 160 implicit none 161#include "errquit.fh" 162#include "mafdecls.fh" 163#include "stdio.fh" 164#include "global.fh" 165#include "msgids.fh" 166 integer g_N,g_R(3),g_M 167 integer g_tmp2,g_N_scld 168 integer ind_a,ind_b 169 call ga_copy(g_N,g_tmp2) 170 call ga_copy(g_N,g_N_scld) 171 call ga_scale_cols(g_tmp2,g_R(ind_b)) ! R_{nu,b} g_N 172 call ga_scale_rows(g_tmp2,g_R(ind_a)) ! R_{mu,a} [R_{nu,b} g_N] -> g_tmp2 173 call ga_scale_cols(g_N_scld,g_R(ind_a)) ! R_{nu,a} g_N 174 call ga_scale_rows(g_N_scld,g_R(ind_b)) ! R_{mu,b} [R_{nu,a} g_N] -> g_N_scld 175 call ga_add(1.0d0,g_tmp2,-1.0d0,g_N_scld,g_N_scld) 176 return 177 end 178 179 subroutine get_scld_A(g_A, ! ga-arr to scale - OUT 180 & g_R, ! scaling arr 181 & g_tmp)! scratch arr 182c 183c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 184 185 implicit none 186#include "errquit.fh" 187#include "mafdecls.fh" 188#include "stdio.fh" 189#include "global.fh" 190#include "msgids.fh" 191 integer g_A,g_R 192 integer g_tmp 193 integer nbf 194c Purpose: Compute R_{nu} U_{munu} - R_{mu} U_{munu} 195c g_R -> R_{mu} 196c g_A -> U_{munu} 197 call ga_copy(g_A,g_tmp) 198 call ga_scale_cols(g_A ,g_R) 199 call ga_scale_rows(g_tmp,g_R) 200 call ga_add(1.0d0,g_A,-1.0d0,g_tmp,g_A) 201 return 202 end 203c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 204c +++++++++++ READ/WRITE NMR-ZORA data +++++++++++++ START 205c Note.- Using modified versions of 206c dft_zora_read() and dft_zoraNMR_write() 207c --> located in dft_zora_utils.F 208czora...Write out the zora NMR shieldings to disk 209 210 logical function dft_zoraNMR_write( 211 & filename, ! in: filename 212 & type_nmrdata, ! in: =1,2,3=shieldings,hyperfine,gshift 213 & nbf, ! in: number of basis functions 214 & nlist, ! in: number of selected atoms 215 & g_AtNr1, ! in: list of atoms to calc. shieldings 216 & g_dia, ! in: dia tensor 217 & g_para1, ! in: para1 tensor 218 & g_h01, ! in: h01 AO matrix 219 & g_Fji) ! in: Perturbed Fock matrix 220c 221c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 222 223 implicit none 224#include "errquit.fh" 225#include "mafdecls.fh" 226#include "global.fh" 227#include "tcgmsg.fh" 228#include "msgtypesf.h" 229#include "inp.fh" 230#include "msgids.fh" 231#include "cscfps.fh" 232#include "util.fh" 233#include "stdio.fh" 234 character*(*) filename ! [input] File to write to 235 integer nbf ! [input] No. of functions in basis 236 integer nlist ! nr. slc atoms 237 integer g_AtNr1, ! in: list of atoms to calc. shieldings 238 & g_dia, ! in: dia tensor 239 & g_para1, ! in: para1 tensor 240 & g_h01, ! in: h01 AO matrix 241 & g_Fji 242 integer unitno 243 parameter (unitno = 77) 244 integer l_AtNr,k_AtNr, 245 & l_tens,k_tens, 246 & l_h01 ,k_h01, 247 & l_Fji ,k_Fji 248 integer ok, iset, i, j 249 integer inntsize 250 integer n9,nxyz,nh01,nFji,type_nmrdata 251 integer alo(3),ahi(3),ld(2) 252 nxyz=3 ! =x,y,z 253 l_tens = -1 ! An invalid MA handle 254 l_h01 = -1 ! An invalid MA handle 255 l_Fji = -1 ! An invalid MA handle 256 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 257 call ga_sync() 258 ok = 0 259c Read routines should be consistent with this 260c Write out the atomic zora corrections 261 if (ga_nodeid() .eq. 0) then 262c Open the file 263 open(unitno, status='unknown', form='unformatted', 264 $ file=filename, err=1000) 265c Write out the number of sets and basis functions 266 write(unitno, err=1001) nbf 267 write(unitno, err=1001) nxyz 268 write(unitno, err=1001) nlist 269c Allocate the temporary buffer 270c ++++++++ using ma_alloc_get +++++++++++++++++++ START 271c ---> ma_alloc_get: to allocate memory 272c ---> ma_free_heap: to release allocated memory 273 if (.not. ma_alloc_get(mt_dbl,nlist,'dft_zoraNMR_write', 274 & l_AtNr,k_AtNr)) 275 $ call errquit('dft_zoraNMR_write: ma failed', 276 & nlist, MA_ERR) 277 n9=nxyz*nxyz*nlist 278 if (.not. ma_alloc_get(mt_dbl,n9,'dft_zoraNMR_write', 279 & l_tens,k_tens)) 280 $ call errquit('dft_zoraNMR_write: ma failed', 281 & n9, MA_ERR) 282 nh01=nbf*nbf*nxyz*nlist 283 if (.not. ma_alloc_get(mt_dbl,nh01,'dft_zoraNMR_write', 284 & l_h01,k_h01)) 285 $ call errquit('dft_zoraNMR_write: ma failed', 286 & nh01, MA_ERR) 287 nFji=nbf*nbf*nxyz 288 if (.not. ma_alloc_get(mt_dbl,nFji,'dft_zoraNMR_write', 289 & l_Fji,k_Fji)) 290 $ call errquit('dft_zoraNMR_write: ma failed', 291 & nFji, MA_ERR) 292c ++++++++ using ma_alloc_get +++++++++++++++++++ END 293 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.2) then 294 call ga_get(g_AtNr1,1,1,1,nlist,dbl_mb(k_AtNr),1) 295 call swrite(unitno,dbl_mb(k_AtNr),nlist) 296 endif 297 alo(1)=1 298 ahi(1)=3 299 alo(2)=1 300 ahi(2)=3 301 alo(3)=1 302 ahi(3)=nlist 303 ld(1)=3 304 ld(2)=3 305 call nga_get(g_dia,alo,ahi,dbl_mb(k_tens),ld) 306 call swrite(unitno,dbl_mb(k_tens),n9) 307 if (type_nmrdata .ne. 1 .and. type_nmrdata .ne. 2 .and. 308 & type_nmrdata .ne. 3) then 309 write(*,*) 'Error in dft_zoraNMR_read:: ', 310 & 'type_nmrdata not correct.', 311 & 'It should be equal 1 or 2 or 3' 312 endif 313 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.3) then 314 call nga_get(g_para1,alo,ahi,dbl_mb(k_tens),ld) 315 call swrite(unitno,dbl_mb(k_tens),n9) 316 endif 317 alo(1)=1 318 ahi(1)=nbf 319 alo(2)=1 320 ahi(2)=nbf 321 alo(3)=1 322 ahi(3)=nxyz*nlist 323 ld(1)=nbf 324 ld(2)=nbf 325 call nga_get(g_h01,alo,ahi,dbl_mb(k_h01),ld) 326 call swrite(unitno,dbl_mb(k_h01),nh01) 327 alo(1)=1 328 ahi(1)=nbf 329 alo(2)=1 330 ahi(2)=nbf 331 alo(3)=1 332 ahi(3)=nxyz 333 ld(1)=nbf 334 ld(2)=nbf 335 call nga_get(g_Fji,alo,ahi,dbl_mb(k_Fji),ld) 336 call swrite(unitno,dbl_mb(k_Fji),nFji) 337c 338c Deallocate the temporary buffer 339c ----- Using ma_free_heap ------------ START 340 if (.not. ma_free_heap(l_AtNr)) 341 $ call errquit('dft_zoraNMR_write: ma free_heap failed', 342 & 911, MA_ERR) 343 if (.not. ma_free_heap(l_tens)) 344 $ call errquit('dft_zoraNMR_write: ma free_heap failed', 345 & 911, MA_ERR) 346 if (.not. ma_free_heap(l_h01)) 347 $ call errquit('dft_zoraNMR_write: ma free_heap failed', 348 & 911, MA_ERR) 349 if (.not. ma_free_heap(l_Fji)) 350 $ call errquit('dft_zoraNMR_write: ma free_heap failed', 351 & 911, MA_ERR) 352c ----- Using ma_free_heap ------------ END 353c Close the file 354 close(unitno,err=1002) 355 ok = 1 356 end if 357c Broadcast status to other nodes 358 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 359 call ga_sync() 360 dft_zoraNMR_write = (ok .eq. 1) 361 if (ga_nodeid() .eq. 0) then 362 write(6,22) filename(1:inp_strlen(filename)) 363 22 format(/' Wrote ZORA NMR data to ',a/) 364 call util_flush(luout) 365 endif 366 return 367 1000 write(6,*) 'dft_zoraNMR_write: failed to open ', 368 $ filename(1:inp_strlen(filename)) 369 call util_flush(luout) 370 ok = 0 371 goto 10 372 1001 write(6,*) 'dft_zoraNMR_write: failed to write ', 373 $ filename(1:inp_strlen(filename)) 374 call util_flush(luout) 375 ok = 0 376 close(unitno,err=1002) 377 goto 10 378 1002 write(6,*) 'dft_zoraNMR_write: failed to close', 379 $ filename(1:inp_strlen(filename)) 380 call util_flush(luout) 381 ok = 0 382 goto 10 383 end 384 385 logical function dft_zoraNMR_read( 386 & filename, ! in: filename 387 & type_nmrdata, ! in: =1,2,3=shieldings,hyperfine,gshift 388 & nbf, ! in: number of basis functions 389 & nlist, ! in: number of selected atoms 390 & g_AtNr1, ! out: list of atoms to calc. shieldings 391 & g_dia, ! out: dia tensor 392 & g_para1, ! out: para1 tensor 393 & g_h01, ! out: h01 AO matrix 394 & g_Fji) ! out: Perturbed Fock matrix 395c 396c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 397 398 implicit none 399#include "errquit.fh" 400#include "global.fh" 401#include "tcgmsg.fh" 402#include "msgtypesf.h" 403#include "mafdecls.fh" 404#include "msgids.fh" 405#include "cscfps.fh" 406#include "inp.fh" 407#include "util.fh" 408#include "stdio.fh" 409 character*(*) filename ! [input] File to write to 410 integer nbf ! [input] No. of functions in basis 411 integer nlist ! nr. slc atoms 412 integer g_AtNr1, ! in: list of atoms to calc. shieldings 413 & g_dia, ! in: dia tensor 414 & g_para1, ! in: para1 tensor 415 & g_h01, ! in: h01 AO matrix 416 & g_Fji 417 integer unitno 418 parameter (unitno = 77) 419 integer l_AtNr,k_AtNr, 420 & l_tens,k_tens, 421 & l_h01 ,k_h01, 422 & l_Fji ,k_Fji 423 integer n9,nxyz,nh01,nFji,type_nmrdata 424 integer alo(3),ahi(3),ld(2) 425 integer ok,inntsize 426 integer nxyz_read,nlist_read, 427 & nbf_read 428 if (type_nmrdata .ne. 1 .and. type_nmrdata .ne. 2 .and. 429 & type_nmrdata .ne. 3) then 430 write(*,*) 'Error in dft_zoraNMR_read::', 431 & ' type_nmrdata not correct.', 432 & 'It should be equal 1 or 2 or 3' 433 stop 434 endif 435 nxyz=3 ! =x,y,z 436c Initialise to invalid MA handle 437 l_tens = -1 ! An invalid MA handle 438 l_h01 = -1 ! An invalid MA handle 439 l_Fji = -1 ! An invalid MA handle 440 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 441 call ga_sync() 442 ok = 0 443c---- Create GA arrays: g_AtNr1,g_dia,g_para1,g_h01,g_Fji --- START 444 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.2) then 445 if (.not. ga_create(mt_dbl,1,nlist, 446 & 'dft_zoraNMR_read: g_AtNr1',0,0,g_AtNr1)) 447 $ call errquit('gCSSR: g_AtNr1',0,GA_ERR) 448 call ga_zero(g_AtNr1) 449 endif 450 alo(1) = 3 451 alo(2) = -1 452 alo(3) = -1 453 ahi(1) = 3 454 ahi(2) = 3 455 ahi(3) = nlist 456 if (.not.nga_create(MT_DBL,3,ahi,'g_DIA matrix',alo,g_dia)) call 457 & errquit('dft_zoraNMR_read: nga_create failed g_dia', 458 & 0,GA_ERR) 459 call ga_zero(g_dia) 460 if (type_nmrdata .eq. 1 .or. type_nmrdata .eq. 3) then 461 if (.not.nga_create(MT_DBL,3,ahi,'gPAR1 matrix', 462 & alo,g_para1)) 463 & call errquit('dft_zoraNMR_read: nga_create failed gpar1', 464 & 0,GA_ERR) 465 call ga_zero(g_para1) 466 endif 467 alo(1) = nbf 468 alo(2) = -1 469 alo(3) = -1 470 ahi(1) = nbf 471 ahi(2) = nbf 472 ahi(3) = nxyz 473 if (.not.nga_create(MT_DBL,3,ahi,'Fji matrix',alo,g_Fji)) 474 & call 475 & errquit('dft_zoraNMR_read: nga_create failed g_Fji', 476 & 0,GA_ERR) 477 call ga_zero(g_Fji) 478 alo(1) = nbf 479 alo(2) = -1 480 alo(3) = -1 481 ahi(1) = nbf 482 ahi(2) = nbf 483 ahi(3) = 3*nlist 484 if (.not.nga_create(MT_DBL,3,ahi,'h01 matrix',alo,g_h01)) call 485 & errquit('dft_zoraNMR_read: nga_create failed g_h01_num', 486 & 0,GA_ERR) 487 call ga_zero(g_h01) 488c---- Create GA arrays: g_AtNr1,g_dia,g_para1,g_h01,g_Fji --- END 489 if (ga_nodeid() .eq. 0) then 490c Print a message indicating the file being read 491 write(6,22) filename(1:inp_strlen(filename)) 492 22 format(/' Read ZORA NMR data from ',a/) 493 call util_flush(luout) 494c Open the file 495 open(unitno, status='old', form='unformatted', file=filename, 496 $ err=1000) 497c Read in some basics to check if they are consistent with the calculation 498 read(unitno, err=1001, end=1001) nbf_read 499 read(unitno, err=1001, end=1001) nxyz_read 500 read(unitno, err=1001, end=1001) nlist_read 501c Error checks 502 if ((nxyz_read .ne. nxyz) .or. 503 & (nbf_read .ne. nbf) .or. 504 & (nlist_read .ne. nlist) ) goto 1003 505c ++++++++ using ma_alloc_get +++++++++++++++++++ START 506c ---> ma_alloc_get: to allocate memory 507c ---> ma_free_heap: to release allocated memory 508 if (.not. ma_alloc_get(mt_dbl,nlist,'dft_zoraNMR_read', 509 & l_AtNr,k_AtNr)) 510 $ call errquit('dft_zoraNMR_read: ma failed', 511 & nlist, MA_ERR) 512 n9=nxyz*nxyz*nlist 513 if (.not. ma_alloc_get(mt_dbl,n9,'dft_zoraNMR_read', 514 & l_tens,k_tens)) 515 $ call errquit('dft_zoraNMR_read: ma failed', 516 & n9, MA_ERR) 517 nh01=nbf*nbf*nxyz*nlist 518 if (.not. ma_alloc_get(mt_dbl,nh01,'dft_zoraNMR_read', 519 & l_h01,k_h01)) 520 $ call errquit('dft_zoraNMR_read: ma failed', 521 & nh01, MA_ERR) 522 nFji=nbf*nbf*nxyz 523 if (.not. ma_alloc_get(mt_dbl,nFji,'dft_zoraNMR_read', 524 & l_Fji,k_Fji)) 525 $ call errquit('dft_zoraNMR_read: ma failed', 526 & nFji, MA_ERR) 527c ++++++++ using ma_alloc_get +++++++++++++++++++ END 528 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.2) then 529 call sread(unitno,dbl_mb(k_AtNr),nlist) 530 call ga_put(g_AtNr1,1,1,1,nlist,dbl_mb(k_AtNr),1) 531 endif 532 alo(1)=1 533 ahi(1)=3 534 alo(2)=1 535 ahi(2)=3 536 alo(3)=1 537 ahi(3)=nlist 538 ld(1)=3 539 ld(2)=3 540 call sread(unitno,dbl_mb(k_tens),n9) 541 call nga_put(g_dia,alo,ahi,dbl_mb(k_tens),ld) 542 if (type_nmrdata .eq. 1 .or. type_nmrdata .eq. 3) then 543 call sread(unitno,dbl_mb(k_tens),n9) 544 call nga_put(g_para1,alo,ahi,dbl_mb(k_tens),ld) 545 endif 546 alo(1)=1 547 ahi(1)=nbf 548 alo(2)=1 549 ahi(2)=nbf 550 alo(3)=1 551 ahi(3)=nxyz*nlist 552 ld(1)=nbf 553 ld(2)=nbf 554 call sread(unitno,dbl_mb(k_h01),nh01) 555 call nga_put(g_h01,alo,ahi,dbl_mb(k_h01),ld) 556 alo(1)=1 557 ahi(1)=nbf 558 alo(2)=1 559 ahi(2)=nbf 560 alo(3)=1 561 ahi(3)=nxyz 562 ld(1)=nbf 563 ld(2)=nbf 564 call sread(unitno,dbl_mb(k_Fji),nFji) 565 call nga_put(g_Fji,alo,ahi,dbl_mb(k_Fji),ld) 566c 567c Deallocate the temporary buffer 568c ----- Using ma_free_heap ------------ START 569 if (.not. ma_free_heap(l_AtNr)) 570 $ call errquit('dft_zoraNMR_read: ma free_heap failed', 571 & 911, MA_ERR) 572 if (.not. ma_free_heap(l_tens)) 573 $ call errquit('dft_zoraNMR_read: ma free_heap failed', 574 & 911, MA_ERR) 575 if (.not. ma_free_heap(l_h01)) 576 $ call errquit('dft_zoraNMR_read: ma free_heap failed', 577 & 911, MA_ERR) 578 if (.not. ma_free_heap(l_Fji)) 579 $ call errquit('dft_zoraNMR_read: ma free_heap failed', 580 & 911, MA_ERR) 581c ----- Using ma_free_heap ------------ END 582c Close the file 583 close(unitno,err=1002) 584 ok = 1 585 end if 586c 587c Broadcast status to other nodes 588 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 589 call ga_sync() 590 dft_zoraNMR_read = ok .eq. 1 591 return 592 1000 write(6,*) 'dft_zoraNMR_read: failed to open', 593 $ filename(1:inp_strlen(filename)) 594 call util_flush(luout) 595 ok = 0 596 goto 10 597 1001 write(6,*) 'dft_zoraNMR_read: failed to read', 598 $ filename(1:inp_strlen(filename)) 599 call util_flush(luout) 600 ok = 0 601 close(unitno,err=1002) 602 goto 10 603 1003 write(6,*) 'dft_zshield_read: file inconsistent with calculation', 604 $ filename(1:inp_strlen(filename)) 605 call util_flush(luout) 606 ok = 0 607 close(unitno,err=1002) 608 goto 10 609 1002 write(6,*) 'dft_zoraNMR_read: failed to close', 610 $ filename(1:inp_strlen(filename)) 611 call util_flush(luout) 612 ok = 0 613 goto 10 614 end 615c --- 05-02-11 ------- writing/reading A,B contributions ----- START 616 logical function dft_zoraNMR_write_AB( 617 & filename, ! in: filename 618 & type_nmrdata, ! in: =1,2,3=shieldings,hyperfine,gshift 619 & nbf, ! in: number of basis functions 620 & nlist, ! in: number of selected atoms 621 & g_AtNr1, ! in: list of atoms to calc. shieldings 622 & g_dia, ! in: dia A,B tensor 623 & g_para1) ! in: par A,B tensor 624c 625c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 626 627 implicit none 628#include "errquit.fh" 629#include "mafdecls.fh" 630#include "global.fh" 631#include "tcgmsg.fh" 632#include "msgtypesf.h" 633#include "inp.fh" 634#include "msgids.fh" 635#include "cscfps.fh" 636#include "util.fh" 637#include "stdio.fh" 638 character*(*) filename ! [input] File to write to 639 integer nbf ! [input] No. of functions in basis 640 integer nlist ! nr. slc atoms 641 integer g_AtNr1, ! in: list of atoms to calc. shieldings 642 & g_dia, ! in: dia tensor 643 & g_para1 ! in: para1 tensor 644 integer unitno 645 parameter (unitno = 77) 646 integer l_AtNr,k_AtNr, 647 & l_tens,k_tens 648 integer ok, iset, i, j 649 integer inntsize 650 integer n9,nxyz,nh01,nFji,type_nmrdata 651 integer alo(3),ahi(3),ld(2) 652 nxyz=3 ! =x,y,z 653 l_tens = -1 ! An invalid MA handle 654 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 655 call ga_sync() 656 ok = 0 657c Read routines should be consistent with this 658c Write out the atomic zora corrections 659 if (ga_nodeid() .eq. 0) then 660c Open the file 661 open(unitno, status='unknown', form='unformatted', 662 $ file=filename, err=1000) 663c Write out the number of sets and basis functions 664 write(unitno, err=1001) nbf 665 write(unitno, err=1001) nxyz 666 write(unitno, err=1001) nlist 667c Allocate the temporary buffer 668c ++++++++ using ma_alloc_get +++++++++++++++++++ START 669c ---> ma_alloc_get: to allocate memory 670c ---> ma_free_heap: to release allocated memory 671 if (.not. ma_alloc_get(mt_dbl,nlist,'dft_zoraNMR_write', 672 & l_AtNr,k_AtNr)) 673 $ call errquit('dft_zoraNMR_write: ma failed', 674 & nlist, MA_ERR) 675 n9=2*nxyz*nxyz*nlist 676 if (.not. ma_alloc_get(mt_dbl,n9,'dft_zoraNMR_write', 677 & l_tens,k_tens)) 678 $ call errquit('dft_zoraNMR_write: ma failed', 679 & n9, MA_ERR) 680c ++++++++ using ma_alloc_get +++++++++++++++++++ END 681 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.2) then 682 call ga_get(g_AtNr1,1,1,1,nlist,dbl_mb(k_AtNr),1) 683 call swrite(unitno,dbl_mb(k_AtNr),nlist) 684 endif 685 alo(1)=1 686 ahi(1)=3 687 alo(2)=1 688 ahi(2)=3 689 alo(3)=1 690 ahi(3)=2*nlist 691 ld(1)=3 692 ld(2)=3 693 call nga_get(g_dia,alo,ahi,dbl_mb(k_tens),ld) 694 call swrite(unitno,dbl_mb(k_tens),n9) 695 if (type_nmrdata .ne. 1 .and. type_nmrdata .ne. 2 .and. 696 & type_nmrdata .ne. 3) then 697 write(*,*) 'Error in dft_zoraNMR_read:: ', 698 & 'type_nmrdata not correct.', 699 & 'It should be equal 1 or 2 or 3' 700 endif 701 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.3) then 702 call nga_get(g_para1,alo,ahi,dbl_mb(k_tens),ld) 703 call swrite(unitno,dbl_mb(k_tens),n9) 704 endif 705c 706c Deallocate the temporary buffer 707c ----- Using ma_free_heap ------------ START 708 if (.not. ma_free_heap(l_AtNr)) 709 $ call errquit('dft_zoraNMR_write: ma free_heap failed', 710 & 911, MA_ERR) 711 if (.not. ma_free_heap(l_tens)) 712 $ call errquit('dft_zoraNMR_write: ma free_heap failed', 713 & 911, MA_ERR) 714c ----- Using ma_free_heap ------------ END 715c Close the file 716 close(unitno,err=1002) 717 ok = 1 718 end if 719c Broadcast status to other nodes 720 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 721 call ga_sync() 722 dft_zoraNMR_write_AB = (ok .eq. 1) 723 if (ga_nodeid() .eq. 0) then 724 write(6,22) filename(1:inp_strlen(filename)) 725 22 format(/' Wrote ZORA NMR data to ',a/) 726 call util_flush(luout) 727 endif 728 return 729 1000 write(6,*) 'dft_zoraNMR_write: failed to open ', 730 $ filename(1:inp_strlen(filename)) 731 call util_flush(luout) 732 ok = 0 733 goto 10 734 1001 write(6,*) 'dft_zoraNMR_write: failed to write ', 735 $ filename(1:inp_strlen(filename)) 736 call util_flush(luout) 737 ok = 0 738 close(unitno,err=1002) 739 goto 10 740 1002 write(6,*) 'dft_zoraNMR_write: failed to close', 741 $ filename(1:inp_strlen(filename)) 742 call util_flush(luout) 743 ok = 0 744 goto 10 745 end 746 747 logical function dft_zoraNMR_read_AB( 748 & filename, ! in: filename 749 & type_nmrdata, ! in: =1,2,3=shieldings,hyperfine,gshift 750 & nbf, ! in: number of basis functions 751 & nlist, ! in: number of selected atoms 752 & g_AtNr1, ! out: list of atoms to calc. shieldings 753 & g_dia, ! out: dia-A,B tensor 754 & g_para1) ! out: par-A,B tensor 755c 756c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 757 758 implicit none 759#include "errquit.fh" 760#include "global.fh" 761#include "tcgmsg.fh" 762#include "msgtypesf.h" 763#include "mafdecls.fh" 764#include "msgids.fh" 765#include "cscfps.fh" 766#include "inp.fh" 767#include "util.fh" 768#include "stdio.fh" 769 character*(*) filename ! [input] File to write to 770 integer nbf ! [input] No. of functions in basis 771 integer nlist ! nr. slc atoms 772 integer g_AtNr1, ! in: list of atoms to calc. shieldings 773 & g_dia, ! in: dia tensor 774 & g_para1 ! in: para1 tensor 775 integer unitno 776 parameter (unitno = 77) 777 integer l_AtNr,k_AtNr, 778 & l_tens,k_tens 779 integer n9,nxyz,nh01,nFji,type_nmrdata 780 integer alo(3),ahi(3),ld(2) 781 integer ok,inntsize 782 integer nxyz_read,nlist_read, 783 & nbf_read 784 if (type_nmrdata .ne. 1 .and. type_nmrdata .ne. 2 .and. 785 & type_nmrdata .ne. 3) then 786 write(*,*) 'Error in dft_zoraNMR_read::', 787 & ' type_nmrdata not correct.', 788 & 'It should be equal 1 or 2 or 3' 789 stop 790 endif 791 nxyz=3 ! =x,y,z 792c Initialise to invalid MA handle 793 l_tens = -1 ! An invalid MA handle 794 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 795 call ga_sync() 796 ok = 0 797c---- Create GA arrays: g_AtNr1,g_dia,g_para1,g_h01,g_Fji --- START 798 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.2) then 799 if (.not. ga_create(mt_dbl,1,nlist, 800 & 'dft_zoraNMR_read: g_AtNr1',0,0,g_AtNr1)) 801 $ call errquit('gCSSR: g_AtNr1',0,GA_ERR) 802 call ga_zero(g_AtNr1) 803 endif 804 alo(1) = 3 805 alo(2) = -1 806 alo(3) = -1 807 ahi(1) = 3 808 ahi(2) = 3 809 ahi(3) = 2*nlist 810 if (.not.nga_create(MT_DBL,3,ahi,'g_DIA matrix',alo,g_dia)) call 811 & errquit('dft_zoraNMR_read: nga_create failed g_dia', 812 & 0,GA_ERR) 813 call ga_zero(g_dia) 814 if (type_nmrdata .eq. 1 .or. type_nmrdata .eq. 3) then 815 if (.not.nga_create(MT_DBL,3,ahi,'gPAR1 matrix', 816 & alo,g_para1)) 817 & call errquit('dft_zoraNMR_read: nga_create failed gpar1', 818 & 0,GA_ERR) 819 call ga_zero(g_para1) 820 endif 821c---- Create GA arrays: g_AtNr1,g_dia,g_para1,g_h01,g_Fji --- END 822 if (ga_nodeid() .eq. 0) then 823c Print a message indicating the file being read 824 write(6,22) filename(1:inp_strlen(filename)) 825 22 format(/' Read ZORA NMR data from ',a/) 826 call util_flush(luout) 827c Open the file 828 open(unitno, status='old', form='unformatted', file=filename, 829 $ err=1000) 830c Read in some basics to check if they are consistent with the calculation 831 read(unitno, err=1001, end=1001) nbf_read 832 read(unitno, err=1001, end=1001) nxyz_read 833 read(unitno, err=1001, end=1001) nlist_read 834c Error checks 835 if ((nxyz_read .ne. nxyz) .or. 836 & (nbf_read .ne. nbf) .or. 837 & (nlist_read .ne. nlist) ) goto 1003 838c ++++++++ using ma_alloc_get +++++++++++++++++++ START 839c ---> ma_alloc_get: to allocate memory 840c ---> ma_free_heap: to release allocated memory 841 if (.not. ma_alloc_get(mt_dbl,nlist,'dft_zoraNMR_read', 842 & l_AtNr,k_AtNr)) 843 $ call errquit('dft_zoraNMR_read: ma failed', 844 & nlist, MA_ERR) 845 n9=2*nxyz*nxyz*nlist 846 if (.not. ma_alloc_get(mt_dbl,n9,'dft_zoraNMR_read', 847 & l_tens,k_tens)) 848 $ call errquit('dft_zoraNMR_read: ma failed', 849 & n9, MA_ERR) 850c ++++++++ using ma_alloc_get +++++++++++++++++++ END 851 if (type_nmrdata.eq.1 .or. type_nmrdata.eq.2) then 852 call sread(unitno,dbl_mb(k_AtNr),nlist) 853 call ga_put(g_AtNr1,1,1,1,nlist,dbl_mb(k_AtNr),1) 854 endif 855 alo(1)=1 856 ahi(1)=3 857 alo(2)=1 858 ahi(2)=3 859 alo(3)=1 860 ahi(3)=2*nlist 861 ld(1)=3 862 ld(2)=3 863 call sread(unitno,dbl_mb(k_tens),n9) 864 call nga_put(g_dia,alo,ahi,dbl_mb(k_tens),ld) 865 if (type_nmrdata .eq. 1 .or. type_nmrdata .eq. 3) then 866 call sread(unitno,dbl_mb(k_tens),n9) 867 call nga_put(g_para1,alo,ahi,dbl_mb(k_tens),ld) 868 endif 869c 870c Deallocate the temporary buffer 871c ----- Using ma_free_heap ------------ START 872 if (.not. ma_free_heap(l_AtNr)) 873 $ call errquit('dft_zoraNMR_read: ma free_heap failed', 874 & 911, MA_ERR) 875 if (.not. ma_free_heap(l_tens)) 876 $ call errquit('dft_zoraNMR_read: ma free_heap failed', 877 & 911, MA_ERR) 878c ----- Using ma_free_heap ------------ END 879c Close the file 880 close(unitno,err=1002) 881 ok = 1 882 end if 883c 884c Broadcast status to other nodes 885 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 886 call ga_sync() 887 dft_zoraNMR_read_AB = ok .eq. 1 888 return 889 1000 write(6,*) 'dft_zoraNMR_read: failed to open', 890 $ filename(1:inp_strlen(filename)) 891 call util_flush(luout) 892 ok = 0 893 goto 10 894 1001 write(6,*) 'dft_zoraNMR_read: failed to read', 895 $ filename(1:inp_strlen(filename)) 896 call util_flush(luout) 897 ok = 0 898 close(unitno,err=1002) 899 goto 10 900 1003 write(6,*) 'dft_zshield_read: file inconsistent with calculation', 901 $ filename(1:inp_strlen(filename)) 902 call util_flush(luout) 903 ok = 0 904 close(unitno,err=1002) 905 goto 10 906 1002 write(6,*) 'dft_zoraNMR_read: failed to close', 907 $ filename(1:inp_strlen(filename)) 908 call util_flush(luout) 909 ok = 0 910 goto 10 911 end 912c --- 05-02-11 ------- writing/reading A,B contributions ----- END 913c +++++++++++ READ/WRITE NMR-ZORA data +++++++++++++ END 914c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 915c ======================================================== 916c =========== READ/WRITE CPHF (g_rhs) data ==========START 917 logical function dft_zoraCPHF_write( 918 & filename, ! in: filename 919 & npol, ! in: nr polarization 920 & nocc, ! in: nr occupied MOs 921 & nvirt, ! in: nr virtual MOs 922 & nbf, ! in: nr basis functions 923 & vectors, ! in: MOs 924 & g_rhs0, ! in: (ntot,3) GA matrix 925 & g_rhs) ! in: (nocc*nvirt,3) GA matrix 926c 927c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 928 929 implicit none 930#include "errquit.fh" 931#include "mafdecls.fh" 932#include "global.fh" 933#include "tcgmsg.fh" 934#include "msgtypesf.h" 935#include "inp.fh" 936#include "msgids.fh" 937#include "cscfps.fh" 938#include "util.fh" 939#include "stdio.fh" 940 character*(*) filename ! [input] File to write to 941 integer npol,nbf, 942 & nocc(npol),nvirt(npol), 943 & vectors(npol), 944 & ispin,ntot, 945 & g_rhs0,g_rhs 946 integer unitno 947 parameter (unitno = 77) 948 integer l_rhs0,k_rhs0, 949 & l_rhs,k_rhs, 950 & l_mo,k_mo 951 integer ok,i,j 952 integer inntsize 953 integer nrhs,nxyz 954 955 nxyz=3 ! =x,y,z 956 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 957 call ga_sync() 958 ok = 0 959c Read routines should be consistent with this 960c Write out the atomic zora corrections 961 if (ga_nodeid() .eq. 0) then 962c Open the file 963 open(unitno, status='unknown', form='unformatted', 964 $ file=filename, err=1000) 965c Write out the number of sets and basis functions 966 write(unitno, err=1001) npol 967 do i=1,npol 968 write(unitno, err=1001) nocc(i) 969 enddo 970 do i=1,npol 971 write(unitno, err=1001) nvirt(i) 972 enddo 973 write(unitno, err=1001) nxyz 974 write(unitno, err=1001) nbf 975c Allocate the temporary buffer 976c ++++++++ using ma_alloc_get +++++++++++++++++++ START 977c ---> ma_alloc_get: to allocate memory 978c ---> ma_free_heap: to release allocated memory 979c ----- Add MOs in file ----- START 980 if (.not. ma_alloc_get( 981 & mt_dbl,nbf,'dft_zoraNLMO_writehyp', 982 & l_mo,k_mo)) 983 $ call errquit('dft_zoraCPHF_write: k_mo failed', 984 & nbf,MA_ERR) 985 do i=1,npol 986 do j=1,nbf 987 call dcopy(nbf,0.0d0,0,dbl_mb(k_mo),1) ! reset 988 call ga_get(vectors(i),1,nbf,j,j,dbl_mb(k_mo),nbf) 989 call swrite(unitno,dbl_mb(k_mo),nbf) 990 enddo ! end-loop-j 991 enddo ! end-loop-i 992c ----- Add MOs in file ----- END 993c ++++++++ using ma_alloc_get +++++++++++++++++++ END 994 ntot=0 995 do ispin=1,npol 996 ntot=ntot+nocc(ispin)*nocc(ispin) 997 enddo 998 write(unitno, err=1001) ntot 999 if (.not. ma_alloc_get(mt_dbl,ntot, 1000 & 'dft_zoraCPHF_write',l_rhs0,k_rhs0)) 1001 $ call errquit('dft_zoraCPHF_write: ma failed', 1002 & ntot, MA_ERR) 1003 do i=1,nxyz 1004 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs0),1) 1005 call ga_get(g_rhs0,1,ntot,i,i,dbl_mb(k_rhs0),ntot) 1006 call swrite(unitno,dbl_mb(k_rhs0),ntot) 1007 enddo 1008 ntot=0 1009 do ispin=1,npol 1010 ntot=ntot+nocc(ispin)*nvirt(ispin) 1011 enddo 1012 write(unitno, err=1001) ntot 1013 if (.not. ma_alloc_get(mt_dbl,ntot, 1014 & 'dft_zoraCPHF_write',l_rhs,k_rhs)) 1015 $ call errquit('dft_zoraCPHF_write: ma failed', 1016 & ntot, MA_ERR) 1017 do i=1,nxyz 1018 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1019 call ga_get(g_rhs,1,ntot,i,i,dbl_mb(k_rhs),ntot) 1020 call swrite(unitno,dbl_mb(k_rhs),ntot) 1021 enddo 1022c 1023c Deallocate the temporary buffer 1024c ----- Using ma_free_heap ------------ START 1025 if (.not. ma_free_heap(l_rhs0)) 1026 $ call errquit('dft_zoraCPHF_write: ma free_heap failed', 1027 & 911, MA_ERR) 1028 if (.not. ma_free_heap(l_rhs)) 1029 $ call errquit('dft_zoraCPHF_write: ma free_heap failed', 1030 & 911, MA_ERR) 1031 if (.not. ma_free_heap(l_mo)) 1032 $ call errquit('dft_zoraCPHF_write: ma free_heap failed', 1033 & 911, MA_ERR) 1034c ----- Using ma_free_heap ------------ END 1035c Close the file 1036 close(unitno,err=1002) 1037 ok = 1 1038 end if 1039c Broadcast status to other nodes 1040 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1041 call ga_sync() 1042 dft_zoraCPHF_write = (ok .eq. 1) 1043 if (ga_nodeid() .eq. 0) then 1044 write(6,22) filename(1:inp_strlen(filename)) 1045 22 format(/' Wrote CPHF data to ',a/) 1046 call util_flush(luout) 1047 endif 1048 return 1049 1000 write(6,*) 'dft_zoraCPHF_write: failed to open ', 1050 $ filename(1:inp_strlen(filename)) 1051 call util_flush(luout) 1052 ok = 0 1053 goto 10 1054 1001 write(6,*) 'dft_zoraCPHF_write: failed to write ', 1055 $ filename(1:inp_strlen(filename)) 1056 call util_flush(luout) 1057 ok = 0 1058 close(unitno,err=1002) 1059 goto 10 1060 1002 write(6,*) 'dft_zoraCPHF_write: failed to close', 1061 $ filename(1:inp_strlen(filename)) 1062 call util_flush(luout) 1063 ok = 0 1064 goto 10 1065 end 1066 1067 logical function dft_zoraCPHF_read( 1068 & filename, ! in: filename 1069 & npol, ! in: nr polarization 1070 & nocc, ! in: nr occupied MOs 1071 & nvirt, ! in: nr virtual MOs 1072 & nbf, ! in: nr basis functions 1073 & vectors, ! out: MOs 1074 & g_rhs0, ! out: (ntot,3) GA matrix 1075 & g_rhs) ! out: (nocc*nvirt,3) GA matrix 1076c 1077c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1078 1079 implicit none 1080#include "errquit.fh" 1081#include "global.fh" 1082#include "tcgmsg.fh" 1083#include "msgtypesf.h" 1084#include "mafdecls.fh" 1085#include "msgids.fh" 1086#include "cscfps.fh" 1087#include "inp.fh" 1088#include "util.fh" 1089#include "stdio.fh" 1090 character*(*) filename ! [input] File to write to 1091 integer npol,nbf, 1092 & nocc(npol),nvirt(npol), 1093 & vectors(npol), 1094 & ispin,ntot, 1095 & g_rhs0,g_rhs 1096 integer unitno 1097 parameter (unitno = 77) 1098 integer l_rhs0,k_rhs0, 1099 & l_rhs,k_rhs, 1100 & l_mo,k_mo 1101 integer ok,i,j 1102 integer inntsize 1103 integer nrhs,nxyz 1104 integer npol_read,nxyz_read,ntot_read, 1105 & nbf_read, 1106 & nocc_read(2),nvirt_read(2) 1107 1108 nxyz=3 ! =x,y,z 1109c Initialise to invalid MA handle 1110 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 1111 call ga_sync() 1112 ok = 0 1113 if (ga_nodeid() .eq. 0) then 1114c Print a message indicating the file being read 1115 write(6,22) filename(1:inp_strlen(filename)) 1116 22 format(/' Read ZORA NMR data from ',a/) 1117 call util_flush(luout) 1118c Open the file 1119 open(unitno, status='old', form='unformatted', file=filename, 1120 $ err=1000) 1121c Read in some basics to check if they are consistent with the calculation 1122 read(unitno, err=1001, end=1001) npol_read 1123 do i=1,npol_read 1124 read(unitno, err=1001, end=1001) nocc_read(i) 1125 enddo 1126 do i=1,npol_read 1127 read(unitno, err=1001, end=1001) nvirt_read(i) 1128 enddo 1129 read(unitno, err=1001, end=1001) nxyz_read 1130 read(unitno, err=1001, end=1001) nbf_read 1131c Error checks 1132 if ((nxyz_read .ne. nxyz) .or. 1133 & (npol_read .ne. npol) .or. 1134 & (nbf_read .ne. nbf) ) goto 1003 1135c ----- Read MOs ----- START 1136 if (.not. ma_alloc_get( 1137 & mt_dbl,nbf,'dft_zoraNLMO_readhyp', 1138 & l_mo,k_mo)) 1139 $ call errquit('dft_zoraNLMO_readhyp: ma failed', 1140 & nbf,MA_ERR) 1141 do i=1,npol 1142 do j=1,nbf 1143 call dcopy(nbf,0.0d0,0,dbl_mb(k_mo),1) ! reset 1144 call sread(unitno,dbl_mb(k_mo),nbf) 1145 call ga_put(vectors(i),1,nbf,j,j,dbl_mb(k_mo),nbf) 1146 enddo ! end-loop-j 1147 enddo ! end-loop-i 1148c ----- Read MOs ----- END 1149 ntot=0 1150 do ispin=1,npol 1151 ntot=ntot+nocc(ispin)*nocc(ispin) 1152 enddo 1153 read(unitno, err=1001, end=1001) ntot_read 1154 if (.not. ma_alloc_get(mt_dbl,ntot, 1155 & 'dft_zoraCPHF_write',l_rhs0,k_rhs0)) 1156 $ call errquit('dft_zoraCPHF_write: ma failed', 1157 & ntot, MA_ERR) 1158 do i=1,nxyz 1159 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs0),1) 1160 call sread(unitno,dbl_mb(k_rhs0),ntot) 1161 call ga_put(g_rhs0,1,ntot,i,i,dbl_mb(k_rhs0),ntot) 1162 enddo 1163 ntot=0 1164 do ispin=1,npol 1165 ntot=ntot+nocc(ispin)*nvirt(ispin) 1166 enddo 1167 read(unitno, err=1001, end=1001) ntot_read 1168 if (.not. ma_alloc_get(mt_dbl,ntot, 1169 & 'dft_zoraCPHF_write',l_rhs,k_rhs)) 1170 $ call errquit('dft_zoraCPHF_write: ma failed', 1171 & ntot, MA_ERR) 1172 do i=1,nxyz 1173 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1174 call sread(unitno,dbl_mb(k_rhs),ntot) 1175 call ga_put(g_rhs,1,ntot,i,i,dbl_mb(k_rhs),ntot) 1176 enddo 1177c 1178c Deallocate the temporary buffer 1179c ----- Using ma_free_heap ------------ START 1180 if (.not. ma_free_heap(l_mo)) ! deallocate memory 1181 $ call errquit('dft_zoraCPHF_read: ma free_heap failed', 1182 & 911, MA_ERR) 1183 if (.not. ma_free_heap(l_rhs0)) 1184 $ call errquit('dft_zoraCPHF_read: ma free_heap failed', 1185 & 911, MA_ERR) 1186 if (.not. ma_free_heap(l_rhs)) 1187 $ call errquit('dft_zoraCPHF_read: ma free_heap failed', 1188 & 911, MA_ERR) 1189c ----- Using ma_free_heap ------------ END 1190c Close the file 1191 close(unitno,err=1002) 1192 ok = 1 1193 end if 1194c 1195c Broadcast status to other nodes 1196 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1197 call ga_sync() 1198 dft_zoraCPHF_read = ok .eq. 1 1199 return 1200 1000 write(6,*) 'dft_zoraCPHF_read: failed to open', 1201 $ filename(1:inp_strlen(filename)) 1202 call util_flush(luout) 1203 ok = 0 1204 goto 10 1205 1001 write(6,*) 'dft_zoraCPHF_read: failed to read', 1206 $ filename(1:inp_strlen(filename)) 1207 call util_flush(luout) 1208 ok = 0 1209 close(unitno,err=1002) 1210 goto 10 1211 1003 write(6,*) 1212 & 'dft_zoraCPHF_read: file inconsistent with calculation', 1213 $ filename(1:inp_strlen(filename)) 1214 call util_flush(luout) 1215 ok = 0 1216 close(unitno,err=1002) 1217 goto 10 1218 1002 write(6,*) 'dft_zoraCPHF_read: failed to close', 1219 $ filename(1:inp_strlen(filename)) 1220 call util_flush(luout) 1221 ok = 0 1222 goto 10 1223 end 1224c =========== READ/WRITE CPHF (g_rhs) data ==========END 1225c $Id$ 1226c =========== READ/WRITE CPHF-1 (g_rhs) data ==========START 1227c To be used in aoresponse module: fiao_f1_movecs 1228 logical function dft_CPHF1_write( 1229 & filename, ! in: filename 1230 & npol, ! in: nr polarization 1231 & nocc, ! in: nr occupied MOs 1232 & nvirt, ! in: nr virtual MOs 1233 & ncomp, ! in: nr. components 1234 & g_rhs_re, ! in: (nocc*nvirt,3)(ipm) GA matrix 1235 & g_rhs_im, ! in: (nocc*nvirt,3)(ipm) GA matrix 1236 & lifetime) ! in: =T if (RE,IM) =F if RE 1237c 1238c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1239c 1240c Note.- nmo ne nbf if linear dependencies appear. 1241 1242 implicit none 1243#include "errquit.fh" 1244#include "mafdecls.fh" 1245#include "global.fh" 1246#include "tcgmsg.fh" 1247#include "msgtypesf.h" 1248#include "inp.fh" 1249#include "msgids.fh" 1250#include "cscfps.fh" 1251#include "util.fh" 1252#include "stdio.fh" 1253 character*(*) filename ! [input] File to write to 1254 integer npol, 1255 & nocc(npol),nvirt(npol), 1256 & ispin,ntot, 1257 & ipm,ncomp, 1258 & g_rhs_re(ncomp), 1259 & g_rhs_im(ncomp) 1260 integer unitno 1261 parameter (unitno = 77) 1262 integer l_rhs0,k_rhs0, 1263 & l_rhs,k_rhs 1264 integer ok,i,j 1265 integer inntsize 1266 integer nrhs,nxyz 1267 logical lifetime 1268 1269 nxyz=3 ! =x,y,z 1270 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 1271 call ga_sync() 1272 ok = 0 1273c Read routines should be consistent with this 1274c Write out the atomic zora corrections 1275 if (ga_nodeid() .eq. 0) then 1276c Open the file 1277 open(unitno, status='unknown', form='unformatted', 1278 $ file=filename, err=1000) 1279c Write out the number of sets and basis functions 1280 write(unitno, err=1001) npol 1281 do i=1,npol 1282 write(unitno, err=1001) nocc(i) 1283 enddo 1284 do i=1,npol 1285 write(unitno, err=1001) nvirt(i) 1286 enddo 1287 write(unitno, err=1001) nxyz 1288c Allocate the temporary buffer 1289 ntot=0 1290 do ispin=1,npol 1291 ntot=ntot+nocc(ispin)*nvirt(ispin) 1292 enddo 1293 write(unitno, err=1001) ntot 1294 if (.not. ma_alloc_get(mt_dbl,ntot, 1295 & 'dft_CPHF_write',l_rhs,k_rhs)) 1296 $ call errquit('dft_CPHF_write: ma failed', 1297 & ntot, MA_ERR) 1298 do ipm=1,ncomp 1299 do i=1,2*nxyz ! write (g_b,g_rhs_sol) 1300 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1301 call ga_get(g_rhs_re(ipm),1,ntot,i,i,dbl_mb(k_rhs),ntot) 1302 call swrite(unitno,dbl_mb(k_rhs),ntot) 1303 enddo ! end-loop-i 1304 if (lifetime) then 1305 do i=1,2*nxyz ! write (g_b,g_rhs_sol) 1306 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1307 call ga_get(g_rhs_im(ipm),1,ntot,i,i,dbl_mb(k_rhs),ntot) 1308 call swrite(unitno,dbl_mb(k_rhs),ntot) 1309 enddo ! end-loop-i 1310 endif ! end-if-lifetime 1311 enddo ! end-loop-ipm 1312c 1313c Deallocate the temporary buffer 1314c ----- Using ma_free_heap ------------ START 1315 if (.not. ma_free_heap(l_rhs)) 1316 $ call errquit('dft_CPHF_write: ma free_heap failed', 1317 & 911, MA_ERR) 1318c ----- Using ma_free_heap ------------ END 1319c Close the file 1320 close(unitno,err=1002) 1321 ok = 1 1322 end if 1323c Broadcast status to other nodes 1324 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1325 call ga_sync() 1326 dft_CPHF1_write = (ok .eq. 1) 1327 if (ga_nodeid() .eq. 0) then 1328 write(6,22) filename(1:inp_strlen(filename)) 1329 22 format(/' dft_CPHF_write: Wrote aoresponse g_rhs data to ',a/) 1330 call util_flush(luout) 1331 endif 1332 return 1333 1000 write(6,*) 'dft_CPHF_write: failed to open ', 1334 $ filename(1:inp_strlen(filename)) 1335 call util_flush(luout) 1336 ok = 0 1337 goto 10 1338 1001 write(6,*) 'dft_CPHF_write: failed to write ', 1339 $ filename(1:inp_strlen(filename)) 1340 call util_flush(luout) 1341 ok = 0 1342 close(unitno,err=1002) 1343 goto 10 1344 1002 write(6,*) 'dft_CPHF_write: failed to close', 1345 $ filename(1:inp_strlen(filename)) 1346 call util_flush(luout) 1347 ok = 0 1348 goto 10 1349 end 1350 1351 logical function dft_CPHF1_read( 1352 & filename, ! in: filename 1353 & npol, ! in: nr polarization 1354 & nocc, ! in: nr occupied MOs 1355 & nvirt, ! in: nr virtual MOs 1356 & ncomp, ! out: nr. components 1357 & g_rhs_re, ! out: (nocc*nvirt,3)(ipm) GA matrix 1358 & g_rhs_im, ! out: (nocc*nvirt,3)(ipm) GA matrix 1359 & lifetime) ! out: =T if (RE,IM) =F if RE 1360c 1361c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1362 1363 implicit none 1364#include "errquit.fh" 1365#include "global.fh" 1366#include "tcgmsg.fh" 1367#include "msgtypesf.h" 1368#include "mafdecls.fh" 1369#include "msgids.fh" 1370#include "cscfps.fh" 1371#include "inp.fh" 1372#include "util.fh" 1373#include "stdio.fh" 1374 character*(*) filename ! [input] File to write to 1375 integer npol, 1376 & nocc(npol),nvirt(npol), 1377 & ispin,ntot, 1378 & ipm,ncomp, 1379 & g_rhs_re(ncomp), 1380 & g_rhs_im(ncomp) 1381 integer unitno 1382 parameter (unitno = 77) 1383 integer l_rhs,k_rhs 1384 integer ok,i,j 1385 integer inntsize 1386 integer nrhs,nxyz 1387 integer npol_read,nxyz_read,ntot_read, 1388 & nocc_read(2),nvirt_read(2) 1389 logical lifetime 1390 nxyz=3 ! =x,y,z 1391c Initialise to invalid MA handle 1392 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 1393 call ga_sync() 1394 ok = 0 1395 if (ga_nodeid() .eq. 0) then 1396c Print a message indicating the file being read 1397 write(6,22) filename(1:inp_strlen(filename)) 1398 22 format(/' dft_CPHF_read:Read aoresponse g_rhs data from ',a/) 1399 call util_flush(luout) 1400c Open the file 1401 open(unitno, status='old', form='unformatted', file=filename, 1402 $ err=1000) 1403c Read in some basics to check if they are consistent with the calculation 1404 read(unitno, err=1001, end=1001) npol_read 1405 do i=1,npol_read 1406 read(unitno, err=1001, end=1001) nocc_read(i) 1407 enddo 1408 do i=1,npol_read 1409 read(unitno, err=1001, end=1001) nvirt_read(i) 1410 enddo 1411 read(unitno, err=1001, end=1001) nxyz_read 1412c Error checks 1413 if ((nxyz_read .ne. nxyz) .or. 1414 & (npol_read .ne. npol)) goto 1003 1415 ntot=0 1416 do ispin=1,npol 1417 ntot=ntot+nocc(ispin)*nvirt(ispin) 1418 enddo 1419 read(unitno, err=1001, end=1001) ntot_read 1420 if (.not. ma_alloc_get(mt_dbl,ntot, 1421 & 'dft_CPHF_read',l_rhs,k_rhs)) 1422 $ call errquit('dft_CPHF_read: ma failed', 1423 & ntot, MA_ERR) 1424 1425 do ipm=1,ncomp 1426 do i=1,nxyz ! skip 1st subspace 1427 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1428 call sread(unitno,dbl_mb(k_rhs),ntot) 1429 enddo 1430 do i=nxyz+1,2*nxyz ! read 2nd subspace and copy to g_rhs_re 1431 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1432 call sread(unitno,dbl_mb(k_rhs),ntot) 1433 call ga_put(g_rhs_re(ipm),1,ntot,i,i,dbl_mb(k_rhs),ntot) 1434 enddo ! end-loop-i 1435 if (lifetime) then 1436 do i=1,nxyz ! skip 1st subspace 1437 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1438 call sread(unitno,dbl_mb(k_rhs),ntot) 1439 enddo ! end-loop-i 1440 do i=nxyz+1,2*nxyz ! read 2nd subspace and copy to g_rhs_im 1441 call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1) 1442 call sread(unitno,dbl_mb(k_rhs),ntot) 1443 call ga_put(g_rhs_im(ipm),1,ntot,i,i,dbl_mb(k_rhs),ntot) 1444 enddo ! end-loop-i 1445 endif ! end-if-lifetime 1446 enddo ! end-loop-ipm 1447c 1448c Deallocate the temporary buffer 1449c ----- Using ma_free_heap ------------ START 1450 if (.not. ma_free_heap(l_rhs)) 1451 $ call errquit('dft_CPHF_read: ma free_heap failed', 1452 & 911, MA_ERR) 1453c ----- Using ma_free_heap ------------ END 1454c Close the file 1455 close(unitno,err=1002) 1456 ok = 1 1457 end if 1458c 1459c Broadcast status to other nodes 1460 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1461 call ga_sync() 1462 dft_CPHF1_read = ok .eq. 1 1463 return 1464 1000 write(6,*) 'dft_CPHF_read: failed to open', 1465 $ filename(1:inp_strlen(filename)) 1466 call util_flush(luout) 1467 ok = 0 1468 goto 10 1469 1001 write(6,*) 'dft_CPHF_read: failed to read', 1470 $ filename(1:inp_strlen(filename)) 1471 call util_flush(luout) 1472 ok = 0 1473 close(unitno,err=1002) 1474 goto 10 1475 1003 write(6,*) 1476 & 'dft_CPHF_read: file inconsistent with calculation', 1477 $ filename(1:inp_strlen(filename)) 1478 call util_flush(luout) 1479 ok = 0 1480 close(unitno,err=1002) 1481 goto 10 1482 1002 write(6,*) 'dft_CPHF_read: failed to close', 1483 $ filename(1:inp_strlen(filename)) 1484 call util_flush(luout) 1485 ok = 0 1486 goto 10 1487 end 1488c 000000000000000000000000000000000000000000000000000000 1489c 000000000000000000000000000000000000000000000000000000 1490 logical function dft_CPHF2_write( 1491 & filename, ! in: filename 1492 & n, ! in: sum_{i=1,npol} nocc(i)*nvirt(i) 1493 & ncomp, ! in: nr. components 1494 & nvec, ! in: nr. of directions = 3 1495 & n1, ! in: =n*ncomp 1496 & nsub, ! in: last subspace index (nsub+1)= nr of subspaces stored 1497 & nsub_file,! ub: subspace counter 1498 & g_z1, ! in: history matrix z 1499 & g_Az1) ! in: history matrix Az 1500c 1501c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1502c 1503c Note.- nsub = cc * nvec (nvec=3=x,y,z) multiple of nvec 1504c dim(g_z1) =(n1,maxsub) maxsub=nvec*11 default 1505c dim(g_Az1)=(n1,maxsub) maxsub=nvec*11 default 1506c nsub <= maxsub 1507c -->It will write only subscape nsub for (z1,Az1) 1508 implicit none 1509#include "errquit.fh" 1510#include "mafdecls.fh" 1511#include "global.fh" 1512#include "tcgmsg.fh" 1513#include "msgtypesf.h" 1514#include "inp.fh" 1515#include "msgids.fh" 1516#include "cscfps.fh" 1517#include "util.fh" 1518#include "stdio.fh" 1519 character*(*) filename ! [input] File to write to 1520 character*255 filename_mini ! only to store nblocks 1521 integer g_z1,g_Az1 1522 integer unitno 1523 parameter (unitno = 77) 1524 integer l_dat,k_dat, 1525 & l_z,k_z 1526 integer ok,i,j,nblock,nblock_file,idat 1527 integer inntsize,g_xre,g_xim,m1,iter 1528 integer n,ncomp,nvec,n1,nsub,nsub_file 1529 double precision val_re,val_im 1530 external conv2reim_z1 ! defined in ga_lkain_2cpl3.F 1531 1532 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 1533 call ga_sync() 1534 ok = 0 1535c Read routines should be consistent with this 1536c Write out the atomic zora corrections 1537 if (ga_nodeid() .eq. 0) then 1538c +++++++++ store nsub +++++++++ START 1539 write(filename_mini,30) trim(filename),'_nblock' 1540 30 format(a,4a) 1541c write(*,*) 'Writing mini:',filename_mini 1542 nblock=nsub/3+1 1543 nblock_file=nsub_file/3+1 1544c write(*,2) nblock,nblock_file 1545c 2 format('Writing (nblock,nblock_file)=(',i5,',',i5,')') 1546 open(unitno, status='unknown', form='unformatted', 1547 $ file=filename_mini, err=1000) 1548 write(unitno, err=1001) nblock_file 1549 close(unitno,err=1002) 1550c +++++++++ store nsub +++++++++ END 1551c Open the file 1552 open(unitno, status='unknown', form='unformatted', 1553 $ file=filename, err=1000,position='append') 1554c Write out the number of sets and basis functions 1555 write(unitno, err=1001) n 1556 write(unitno, err=1001) ncomp 1557 write(unitno, err=1001) nvec 1558 write(unitno, err=1001) n1 1559 write(unitno, err=1001) nsub_file 1560c Allocate the temporary buffer 1561 if (.not. ma_alloc_get(mt_dcpl,n1, 1562 & 'dft_CPHF2_write',l_z,k_z)) 1563 $ call errquit('dft_CPHF2_write: ma failed', 1564 & n1, MA_ERR) 1565 if (.not. ma_alloc_get(mt_dbl,n1, 1566 & 'dft_CPHF2_write',l_dat,k_dat)) 1567 $ call errquit('dft_CPHF2_write: ma failed', 1568 & n1, MA_ERR) 1569c ------- write g_z1 ---------------------- START 1570 do i=1,nvec 1571 m1=(nblock-1)*nvec+i 1572 call ga_get(g_z1,1,n1,m1,m1,dcpl_mb(k_z),n1) 1573 do idat=1,n1 1574 val_re=dreal(dcpl_mb(k_z+idat-1)) 1575 dbl_mb(k_dat+idat-1)=val_re 1576 enddo ! end-loop-idat 1577 call swrite(unitno,dbl_mb(k_dat),n1) 1578 do idat=1,n1 1579 val_im=dimag(dcpl_mb(k_z+idat-1)) 1580 dbl_mb(k_dat+idat-1)=val_im 1581 enddo ! end-loop-idat 1582 call swrite(unitno,dbl_mb(k_dat),n1) 1583 enddo ! end-loop-i 1584c ------- write g_z1 ---------------------- END 1585c ------- write g_Az1 ---------------------- START 1586 do i=1,nvec 1587 m1=(nblock-1)*nvec+i 1588 call ga_get(g_Az1,1,n1,m1,m1,dcpl_mb(k_z),n1) 1589 do idat=1,n1 1590 val_re=dreal(dcpl_mb(k_z+idat-1)) 1591 dbl_mb(k_dat+idat-1)=val_re 1592 enddo ! end-loop-idat 1593 call swrite(unitno,dbl_mb(k_dat),n1) 1594 do idat=1,n1 1595 val_im=dimag(dcpl_mb(k_z+idat-1)) 1596 dbl_mb(k_dat+idat-1)=val_im 1597 enddo ! end-loop-idat 1598 call swrite(unitno,dbl_mb(k_dat),n1) 1599 enddo ! end-loop-i 1600c ------- write g_Az1 ---------------------- END 1601c Deallocate the temporary buffer 1602c ----- Using ma_free_heap ------------ START 1603 if (.not. ma_free_heap(l_dat)) 1604 $ call errquit('dft_CPHF2_write: ma free_heap failed', 1605 & 911, MA_ERR) 1606 if (.not. ma_free_heap(l_z)) 1607 $ call errquit('dft_CPHF2_write: ma free_heap failed', 1608 & 911, MA_ERR) 1609c ----- Using ma_free_heap ------------ END 1610c Close the file 1611 close(unitno,err=1002) 1612 ok = 1 1613 end if 1614c Broadcast status to other nodes 1615 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1616 call ga_sync() 1617 dft_CPHF2_write = (ok .eq. 1) 1618 if (ga_nodeid() .eq. 0) then 1619 iter=nsub/3+1 ! estimate iter from nsub 1620 write(6,22) filename(1:inp_strlen(filename)), 1621 & iter,nsub,nsub_file 1622 22 format(' dft_CPHF2_write: Wrote ',a, 1623 & ' (iter,nsub,nsub_file)=(', 1624 & i4,',',i5,',',i5,')') 1625 call util_flush(luout) 1626 endif 1627 return 1628 1000 write(6,*) 'dft_CPHF2_write: failed to open ', 1629 $ filename(1:inp_strlen(filename)) 1630 call util_flush(luout) 1631 ok = 0 1632 goto 10 1633 1001 write(6,*) 'dft_CPHF2_write: failed to write ', 1634 $ filename(1:inp_strlen(filename)) 1635 call util_flush(luout) 1636 ok = 0 1637 close(unitno,err=1002) 1638 goto 10 1639 1002 write(6,*) 'dft_CPHF2_write: failed to close', 1640 $ filename(1:inp_strlen(filename)) 1641 call util_flush(luout) 1642 ok = 0 1643 goto 10 1644 end 1645 1646 logical function dft_CPHF2_read( 1647 & filename, ! in: filename 1648 & n, ! in: sum_{i=1,npol} nocc(i)*nvirt(i) 1649 & ncomp, ! in: nr. components 1650 & nvec, ! in: nr. of directions = 3 1651 & n1, ! in: =n*ncomp 1652 & nsub, ! ou: last subspace index (nsub+1)= nr of subspaces stored 1653 & nsub_read,! ou: last subspace read from file 1654 & maxsub, ! in: max subspace of (g_z1,g_Az1) 1655 & g_z1, ! ou: history matrix z 1656 & g_Az1) ! ou: history matrix Az 1657c 1658c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1659 1660 implicit none 1661#include "errquit.fh" 1662#include "global.fh" 1663#include "tcgmsg.fh" 1664#include "msgtypesf.h" 1665#include "mafdecls.fh" 1666#include "msgids.fh" 1667#include "cscfps.fh" 1668#include "inp.fh" 1669#include "util.fh" 1670#include "stdio.fh" 1671 character*(*) filename ! [input] File to write to 1672 character*255 filename_mini ! only to store nsub 1673 integer ivec,idat, 1674 & g_z1,g_Az1 1675 integer unitno 1676 parameter (unitno = 77) 1677 integer l_zre,k_zre, 1678 & l_zim,k_zim 1679 integer ok,i,j,m1,m2, 1680 & imin,imax,iskip,nskip 1681 integer inntsize 1682 integer n,ncomp,nvec,n1,nsub,maxsub,iblock, 1683 & n_read,n1_read,nvec_red,ncomp_read, 1684 & nsub_read,nvec_read, 1685 & nblocks_ga,nblocks_file 1686 double precision val_re,val_im 1687 double complex val_cmplx 1688 external conv2reim_z1 ! defined in ga_lkain_2cpl3.F 1689 1690 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 1691 call ga_sync() 1692 ok = 0 1693 imin=0 1694 imax=0 1695 if (ga_nodeid() .eq. 0) then ! ----- if-ga_nodeid-eq-0 ----- START 1696c +++++++++ read nsub +++++++++ START 1697 write(filename_mini,30) trim(filename),'_nblock' 1698 30 format(a,4a) 1699c write(*,*) 'Reading mini:',filename_mini 1700 open(unitno, status='old', form='unformatted', 1701 $ file=filename_mini, err=1000) 1702 read(unitno, err=1001, end=1001) nblocks_file 1703 close(unitno,err=1002) 1704 nblocks_ga=maxsub/nvec 1705 write(*,2) nblocks_file,nblocks_ga 1706 2 format('(nblocks_file,nblocks_ga)=(',i4,',',i4,')') 1707c --- Compare (nblocks_ga,nblocks_file) 1708 if (nblocks_file .le. nblocks_ga-1) then 1709 imin=1 1710 imax=nblocks_file 1711 nskip=0 1712 else 1713 imin=1 1714 imax=nblocks_ga-1 1715 nskip=nblocks_file-(nblocks_ga-1) 1716 endif 1717 write(*,3) imin,imax,nskip 1718 3 format('(imin,imax,nskip)=(',i4,',',i4,',',i4,')') 1719c +++++++++ read nsub +++++++++ END 1720c Print a message indicating the file being read 1721 write(6,22) filename(1:inp_strlen(filename)) 1722 22 format(/' dft_CPHF_read:Read aoresponse g_rhs data from ',a/) 1723 call util_flush(luout) 1724c Open the file 1725 open(unitno, status='old', form='unformatted', file=filename, 1726 $ err=1000) 1727c Read in some basics to check if they are consistent with the calculation 1728 if (.not.MA_Push_Get(mt_dbl,n1,'hessv jfacs',l_zre,k_zre)) 1729 & call errquit('conv2complex: cannot allocate zre', 1730 & n1, MA_ERR) 1731 if (.not.MA_Push_Get(mt_dbl,n1,'hessv kfacs',l_zim,k_zim)) 1732 & call errquit('conv2complex: cannot allocate zim', 1733 & n1, MA_ERR) 1734c ======= skip blocks ======================= START 1735 do iskip=1,nskip 1736 read(unitno, err=1001, end=1001) n_read 1737 read(unitno, err=1001, end=1001) ncomp_read 1738 read(unitno, err=1001, end=1001) nvec_read 1739 read(unitno, err=1001, end=1001) n1_read 1740 read(unitno, err=1001, end=1001) nsub_read 1741c ------- skip g_z1 ---------------------- START 1742 do ivec=1,nvec 1743 call sread(unitno,dbl_mb(k_zre),n1) 1744 call sread(unitno,dbl_mb(k_zim),n1) 1745 enddo ! end-loop-i 1746c ------- skip g_z1 ---------------------- END 1747c ------- skip g_Az1 ---------------------- START 1748 do ivec=1,nvec 1749 call sread(unitno,dbl_mb(k_zre),n1) 1750 call sread(unitno,dbl_mb(k_zim),n1) 1751 enddo ! end-loop-i 1752c ------- skip g_Az1 ---------------------- END 1753 enddo ! end-loop-iskip 1754c ======= skip blocks ======================= END 1755c ======= Loop in subspaces ===== START 1756 do iblock=imin,imax 1757 read(unitno, err=1001, end=1001) n_read 1758 read(unitno, err=1001, end=1001) ncomp_read 1759 read(unitno, err=1001, end=1001) nvec_read 1760 read(unitno, err=1001, end=1001) n1_read 1761 read(unitno, err=1001, end=1001) nsub_read 1762 write(*,14) iblock,nsub_read 1763 14 format('(iblock,nsub_read)=(',i4,',',i4,')') 1764c ------- read g_z1 ---------------------- START 1765 m1=(iblock-1)*nvec+1 1766 m2=m1+nvec-1 1767c write(*,4) m1,m2,nvec 1768c 4 format('(m1,m2,nvec)=(',i4,',',i4,',',i4,')') 1769 do ivec=m1,m2 1770 call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1) 1771 call sread(unitno,dbl_mb(k_zre),n1) 1772 call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1) 1773 call sread(unitno,dbl_mb(k_zim),n1) 1774 do idat=1,n1 1775 val_cmplx=dcmplx(dbl_mb(k_zre+idat-1), 1776 & dbl_mb(k_zim+idat-1)) 1777 call ga_put(g_z1,idat,idat,ivec,ivec,val_cmplx,1) 1778 enddo ! end-loop-idat 1779 enddo ! end-loop-ivec 1780c ------- read g_z1 ---------------------- END 1781c ------- read g_Az1 ---------------------- START 1782 do ivec=m1,m2 1783 call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1) 1784 call sread(unitno,dbl_mb(k_zre),n1) 1785 call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1) 1786 call sread(unitno,dbl_mb(k_zim),n1) 1787 do idat=1,n1 1788 val_cmplx=dcmplx(dbl_mb(k_zre+idat-1), 1789 & dbl_mb(k_zim+idat-1)) 1790 call ga_put(g_Az1,idat,idat,ivec,ivec,val_cmplx,1) 1791 enddo ! end-loop-idat 1792 enddo ! end-loop-i 1793c ------- read g_Az1 ---------------------- END 1794 enddo ! end-loop-iblock 1795 nsub=(imax-1)*3 1796 write(*,5) nsub,nblocks_file,nblocks_ga 1797 5 format('dft_CPHF2_read: (nsub,nblocks_file,nblocks_ga)=(', 1798 & i12,',',i12,',',i12,')') 1799c ======= Loop in subspaces ===== END 1800c Deallocate the temporary buffer 1801c ----- Using ma_free_heap ------------ START 1802 if (.not.ma_pop_stack(l_zim)) 1803 $ call errquit('dft_CPHF2_read: pop problem with l_zim', 1804 & 555,MA_ERR) 1805 if (.not.ma_pop_stack(l_zre)) 1806 $ call errquit('dft_CPHF2_read: pop problem with l_zre', 1807 & 555,MA_ERR) 1808c ----- Using ma_free_heap ------------ END 1809c Close the file 1810 close(unitno,err=1002) 1811 ok = 1 1812 end if ! ----- if-ga_nodeid-eq-0 ----- END 1813c Broadcast status to other nodes 1814 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1815 call ga_sync() 1816 dft_CPHF2_read = ok .eq. 1 1817 return 1818 1000 write(6,*) 'dft_CPHF2_read: failed to open', 1819 $ filename(1:inp_strlen(filename)) 1820 call util_flush(luout) 1821 ok = 0 1822 goto 10 1823 1001 write(6,*) 'dft_CPHF2_read: failed to read', 1824 $ filename(1:inp_strlen(filename)) 1825 call util_flush(luout) 1826 ok = 0 1827 close(unitno,err=1002) 1828 goto 10 1829 1003 write(6,*) 1830 & 'dft_CPHF2_read: file inconsistent with calculation', 1831 $ filename(1:inp_strlen(filename)) 1832 call util_flush(luout) 1833 ok = 0 1834 close(unitno,err=1002) 1835 goto 10 1836 1002 write(6,*) 'dft_CPHF2_read: failed to close', 1837 $ filename(1:inp_strlen(filename)) 1838 call util_flush(luout) 1839 ok = 0 1840 goto 10 1841 end 1842 1843 logical function dft_CPHF2_read2fix( 1844 & filename, ! in: filename 1845 & n, ! in: sum_{i=1,npol} nocc(i)*nvirt(i) 1846 & ncomp, ! in: nr. components 1847 & nvec, ! in: nr. of directions = 3 1848 & n1, ! in: =n*ncomp 1849 & nsub, ! ou: last subspace index (nsub+1)= nr of subspaces stored 1850 & maxsub, ! in: max subspace of (g_z1,g_Az1) 1851 & g_z1, ! ou: history matrix z 1852 & g_Az1) ! ou: history matrix Az 1853c 1854c Author : Fredy W. Aquino, Northwestern University (Oct 2012) 1855 1856 implicit none 1857#include "errquit.fh" 1858#include "global.fh" 1859#include "tcgmsg.fh" 1860#include "msgtypesf.h" 1861#include "mafdecls.fh" 1862#include "msgids.fh" 1863#include "cscfps.fh" 1864#include "inp.fh" 1865#include "util.fh" 1866#include "stdio.fh" 1867 character*(*) filename ! [input] File to write to 1868 character*255 filename_mini ! only to store nsub 1869 character*255 filename_mini_1 ! only to store nsub 1870 character*255 filename_1 ! only to store nsub 1871 integer ivec,idat, 1872 & g_z1,g_Az1 1873 integer unitno 1874 parameter (unitno = 77) 1875 1876 integer unitno1 1877 parameter (unitno1 = 78) 1878 1879 integer l_zre,k_zre, 1880 & l_zim,k_zim 1881 integer ok,i,j,m1,m2, 1882 & imin,imax,iskip,nskip,nsub1 1883 integer inntsize 1884 integer n,ncomp,nvec,n1,nsub,maxsub,iblock, 1885 & n_read,n1_read,nvec_red,ncomp_read, 1886 & nsub_read,nvec_read, 1887 & nblocks_ga,nblocks_file, 1888 & ntotblock_true 1889 double complex val_cmplx 1890 external conv2reim_z1 ! defined in ga_lkain_2cpl3.F 1891c 00000000000000000000000000 1892 ntotblock_true=75 ! FePc 1893c 00000000000000000000000000 1894 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 1895 call ga_sync() 1896 ok = 0 1897 if (ga_nodeid() .eq. 0) then 1898 write(filename_1,31) trim(filename),'_1' 1899 31 format(a,2a) 1900 write(filename_mini_1,32) trim(filename),'_nblock_1' 1901 32 format(a,9a) 1902 write(*,33) filename_1(1:inp_strlen(filename_1)), 1903 & filename_mini_1(1:inp_strlen(filename_mini_1)) 1904 33 format('Creating files: filename_1=',a, 1905 & ' filename_mini_1=',a) 1906c +++++++++ read nsub +++++++++ START 1907 write(filename_mini,30) trim(filename),'_nblock' 1908 30 format(a,4a) 1909c write(*,*) 'Reading mini:',filename_mini 1910 open(unitno, status='unknown', form='unformatted', 1911 $ file=filename_mini, err=1000) 1912 read(unitno, err=1001, end=1001) nblocks_file 1913 close(unitno,err=1002) 1914 1915 write(*,*) 'Writing nblock=',ntotblock_true 1916 open(unitno1, status='unknown', form='unformatted', 1917 $ file=filename_mini_1, err=1000) 1918 write(unitno1, err=1001) ntotblock_true 1919 close(unitno1,err=1002) 1920 1921 nblocks_ga=maxsub/nvec 1922 write(*,2) nblocks_file,nblocks_ga 1923 2 format('(nblocks_file,nblocks_ga)=(',i4,',',i4,')') 1924c --- Compare (nblocks_ga,nblocks_file) 1925 if (nblocks_file .le. nblocks_ga) then 1926 imin=1 1927 imax=nblocks_file 1928 nskip=0 1929 else 1930 imin=1 1931 imax=nblocks_ga 1932 nskip=nblocks_file-nblocks_ga 1933 endif 1934 write(*,3) imin,imax,nskip 1935 3 format('(imin,imax,nskip)=(',i4,',',i4,',',i4,')') 1936c +++++++++ read nsub +++++++++ END 1937c Print a message indicating the file being read 1938 write(6,22) filename(1:inp_strlen(filename)) 1939 22 format(/' dft_CPHF_read:Read aoresponse g_rhs data from ',a/) 1940 call util_flush(luout) 1941c Open the file 1942 open(unitno, status='old', form='unformatted', file=filename, 1943 $ err=1000) 1944 1945 open(unitno1, status='unknown', form='unformatted', 1946 $ file=filename_1, err=1000,position='append') 1947 1948c Read in some basics to check if they are consistent with the calculation 1949 if (.not.MA_Push_Get(mt_dbl,n1,'hessv jfacs',l_zre,k_zre)) 1950 & call errquit('conv2complex: cannot allocate zre', 1951 & n1, MA_ERR) 1952 if (.not.MA_Push_Get(mt_dbl,n1,'hessv kfacs',l_zim,k_zim)) 1953 & call errquit('conv2complex: cannot allocate zim', 1954 & n1, MA_ERR) 1955c ======= Loop in subspaces ===== START 1956 nsub1=0 1957 do iblock=1,ntotblock_true 1958 read(unitno, err=1001, end=1001) n_read 1959 read(unitno, err=1001, end=1001) ncomp_read 1960 read(unitno, err=1001, end=1001) nvec_read 1961 read(unitno, err=1001, end=1001) n1_read 1962 read(unitno, err=1001, end=1001) nsub_read 1963 1964 write(unitno1, err=1001) n_read 1965 write(unitno1, err=1001) ncomp_read 1966 write(unitno1, err=1001) nvec_read 1967 write(unitno1, err=1001) n1_read 1968 write(unitno1, err=1001) nsub1 1969 nsub1=nsub1+nvec 1970 write(*,14) iblock,nsub1,nsub_read 1971 14 format('(iblock,nsub1,nsub_read)=(', 1972 & i4,',',i4,',',i4,')') 1973c ------- read g_z1 ---------------------- START 1974 m1=(iblock-1)*nvec+1 1975 m2=m1+nvec-1 1976c write(*,4) m1,m2,nvec 1977c 4 format('(m1,m2,nvec)=(',i4,',',i4,',',i4,')') 1978 do ivec=m1,m2 1979 call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1) 1980 call sread(unitno,dbl_mb(k_zre),n1) 1981 call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1) 1982 call sread(unitno,dbl_mb(k_zim),n1) 1983 call swrite(unitno1,dbl_mb(k_zre),n1) 1984 call swrite(unitno1,dbl_mb(k_zim),n1) 1985 enddo ! end-loop-i 1986c ------- read g_z1 ---------------------- END 1987c ------- read g_Az1 ---------------------- START 1988 do ivec=m1,m2 1989 call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1) 1990 call sread(unitno,dbl_mb(k_zre),n1) 1991 call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1) 1992 call sread(unitno,dbl_mb(k_zim),n1) 1993 call swrite(unitno1,dbl_mb(k_zre),n1) 1994 call swrite(unitno1,dbl_mb(k_zim),n1) 1995 enddo ! end-loop-i 1996c ------- read g_Az1 ---------------------- END 1997 enddo ! end-loop-iblock 1998c ======= Loop in subspaces ===== END 1999c Deallocate the temporary buffer 2000c ----- Using ma_free_heap ------------ START 2001 if (.not.ma_pop_stack(l_zim)) 2002 $ call errquit('dft_CPHF2_read: pop problem with l_zim', 2003 & 555,MA_ERR) 2004 if (.not.ma_pop_stack(l_zre)) 2005 $ call errquit('dft_CPHF2_read: pop problem with l_zre', 2006 & 555,MA_ERR) 2007c ----- Using ma_free_heap ------------ END 2008c Close the file 2009 close(unitno,err=1002) 2010 close(unitno1,err=1002) 2011 ok = 1 2012 end if 2013c 2014c Broadcast status to other nodes 2015 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 2016 call ga_sync() 2017 dft_CPHF2_read2fix = ok .eq. 1 2018 return 2019 1000 write(6,*) 'dft_CPHF2_read: failed to open', 2020 $ filename(1:inp_strlen(filename)) 2021 call util_flush(luout) 2022 ok = 0 2023 goto 10 2024 1001 write(6,*) 'dft_CPHF2_read: failed to read', 2025 $ filename(1:inp_strlen(filename)) 2026 call util_flush(luout) 2027 ok = 0 2028 close(unitno,err=1002) 2029 goto 10 2030 1003 write(6,*) 2031 & 'dft_CPHF2_read: file inconsistent with calculation', 2032 $ filename(1:inp_strlen(filename)) 2033 call util_flush(luout) 2034 ok = 0 2035 close(unitno,err=1002) 2036 goto 10 2037 1002 write(6,*) 'dft_CPHF2_read: failed to close', 2038 $ filename(1:inp_strlen(filename)) 2039 call util_flush(luout) 2040 ok = 0 2041 goto 10 2042 end 2043c =========== READ/WRITE CPHF-1 g_rhs) data ==========END 2044