1c======================================================= 2c================ Fredy Aquino's routines ======== START 3c======================================================= 4 subroutine get_rhoS(rtdb,g_dens_ini,nexc, 5 & geom, 6 & ao_bas_han, 7 & nbf,nbf_ao, 8 & noc, 9 & ipol) 10c -- Purpose: Calculation of small component density, \rho_S 11c for Electric Field Gradient calculation 12c Source: van Lenthe, et.al.,JCP, V112,N19,Y2000 13c -- Author : Fredy Aquino 12-07-09 14 15 implicit none 16#include "errquit.fh" 17#include "mafdecls.fh" 18#include "stdio.fh" 19#include "global.fh" 20#include "msgids.fh" 21#include "rtdb.fh" 22#include "geom.fh" 23#include "zora.fh" 24 25 integer rtdb,geom,ao_bas_han ! handles 26 integer g_efgz4(2),g_dens_ini(2),g_densZ4(3), 27 & g_zora_scale_munu(2) 28 integer noc(2),noc1,ipol,ispin,pos_dat, 29 & nbf,nbf_ao,pos,ndir, 30 & do_zgc_old,nexc 31 integer iat,iat1,nat,ipolmunu,zora_Qpq 32 integer l_xyzpt,k_xyzpt, 33 & l_zanpt,k_zanpt, 34 & l_ilst,k_ilst, 35 & l_jlst,k_jlst, 36 & l_AtNr,k_AtNr, 37 & nlst,i,j,count 38 integer g_rhoS,g_munu_rhoS,g_munuV6, 39 & efgfile 40 logical dft_zoraEFGZ4_write 41c g_munu_rhoS, contains xx,yy,zz,xy,xz,yz 42c munu half-triangle matrices 43c including main diagonal 44 logical status 45 character*16 element, at_tag 46 character*255 zorafilename 47 double precision scf_dbl,zora_eint,toac, 48 & xyz_EFGcoords(3) 49 integer stat_read,read_SLCTD_EFG_Atoms 50 logical dft_zoraEFGZ4_NLMOAnalysis_write 51 external get_densZ4,read_SLCTD_EFG_Atoms, 52 & zora_getv_EFGZ4_SR,dft_zoraEFGZ4_write, 53 & util_file_name, 54 & dft_zoraEFGZ4_NLMOAnalysis_write 55c Note.- g_munu_rhoS is created in dft_zora_rhos.F 56c g_munuV6 is created in hnd_elfcon_symm.F 57c g_munuV6 is not created/defined here yet 58c 59c---- Input global variables (defined in zora.fh): 60c 1. zora_calc_type ! =3 -> ZORA4-EFG 61c 2. so_term ! =0 -> ZORA4-spin-free 62c 3. xyz_EFGcoords(i) i=1,2,3 63c 4. zora_Qpq=1,6 : xx,yy,zz,xy,xz,yz 64 if (ipol.eq.1) then 65 scf_dbl=2.0d00 66 noc1=noc(1) 67 else if (ipol.eq.2) then 68 scf_dbl=1.0d00 69 noc1=noc(1)+noc(2) 70 endif 71 status=geom_ncent(geom,nat) 72c----- Allocate memory - FA 73 if (.not. ma_alloc_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 74 & call errquit('get_rhoS: ma failed',911,MA_ERR) 75 if (.not. ma_alloc_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 76 & call errquit('get_rhoS: ma failed',911,MA_ERR) 77c----- Allocate global arrays - FA 78 do i=1,ipol 79 if (.not. ga_create(mt_dbl,nbf,nbf, 80 & 'get_rhoS: g_efgz4',0,0,g_efgz4(i))) 81 $ call errquit('get_rhoS: g_efgz4', 0,GA_ERR) 82 call ga_zero(g_efgz4(i)) 83 enddo 84c +++++ Read Atom Nr for EFG calc ++ 85 if (.not. ga_create(mt_dbl,1,nat, 86 & 'get_rhoS: g_AtNr',0,0,g_AtNr)) 87 $ call errquit('get_rhoS: g_AtNr', 0,GA_ERR) 88 call ga_zero(g_AtNr) 89 stat_read=read_SLCTD_EFG_Atoms 90 & (rtdb,nat,nlist,g_AtNr) 91 if (ga_nodeid().eq.0) 92 & write(*,*) "==> nlist=",nlist 93c Allocate memory for l_AtNr,k_AtNr 94 if (.not.ma_alloc_get(mt_dbl,nat,'AtNr',l_AtNr,k_AtNr)) 95 & call errquit('get_rhoS: ma failed',0,MA_ERR) 96 call ga_get(g_AtNr,1,1,1,nat,dbl_mb(k_AtNr),1) 97c ------------ for NLMO analysis --------- START 98 efgfile=0 ! not doing NLMO analysis by default 99 status=rtdb_get(rtdb,'prop:efgfile',mt_int,1,efgfile) ! for NLMO analysis 100 if (efgfile.eq.1) then ! ============= efgfile-if-START 101c Allocate memory for (l_ilst,k_ilst) (l_jlst,k_jlst) 102 nlst=nbf*(nbf+1)/2 ! single triangle matrix 103 if (.not.ma_alloc_get(mt_int,nlst,'ijlst', 104 & l_ilst,k_ilst)) 105 & call errquit('get_ijlst: ma failed',0,MA_ERR) 106 if (.not.ma_alloc_get(mt_int,nlst,'ijlst', 107 & l_jlst,k_jlst)) 108 & call errquit('get_ijlst: ma failed',0,MA_ERR) 109c ------- create ga_munu_rhoS --- START 110c WARNING : ONLY for one atom 111 ndir=6 ! Nr. directions: xx,yy,zz,xy,xz,yz 112 if (.not. ga_create(mt_dbl,1,nlst*ndir*nlist, 113 & 'get_munu_rhoS: g_munu_rhoS', 0,0,g_munu_rhoS)) 114 $ call errquit('get_rhoS: g_AtNr', 0,GA_ERR) 115 call ga_zero(g_munu_rhoS) 116c ------- create ga_munu_rhoS --- END 117c Define (ilst,jlst) indices to store munu-rhoS elem. 118 count=0 119 do i=1,nbf 120 int_mb(k_ilst+count)=i 121 int_mb(k_jlst+count)=i 122 count=count+1 123 enddo 124 do i=2,nbf 125 do j=1,i-1 126 int_mb(k_ilst+count)=i 127 int_mb(k_jlst+count)=j 128 count=count+1 129 enddo 130 enddo 131 endif !=============================== efgfile-if-END 132c ------------ for NLMO analysis --------- END 133c ++++++++++++++++++++++++++++++++++ 134c--- About content of g_rhoS: 135c--- 1st set of nat*6 elements corresponds to ZORA4-EFG 136c--- 2nd set of nat*6 elements corresponds to NUM-EFG 137 if (.not. ga_create(mt_dbl,1,nlist*6*2, 138 & 'get_rhoS: g_rhoS', 139 $ 0,0,g_rhoS)) 140 & call errquit('get_rhoS: g_rhoS', 0, 141 & GA_ERR) 142 efgfile=0 ! not doing NLMO analysis by default 143 status=rtdb_get(rtdb,'prop:efgfile',mt_int,1,efgfile) ! for NLMO analysis 144 145 pos_dat=1 ! data-counter in g_rhoS 146c do zora_calc_type =3,4 ! ZORA4-EFG,NUM-EFG 147 do zora_calc_type =4,3,-1 ! NUM-EFG,ZORA4-EFG 148 do_NonRel=.false. 149 if (zora_calc_type.eq.4) do_NonRel=.true. 150 call get_densZ4(rtdb,ao_bas_han,geom,g_densZ4) 151 so_term=0 ! ZORA-spin-free 152 pos=0 ! chunk-index counter for g_munu_rhoS 153 do iat1=1,nlist ! nlist <= nat 154 iat=dbl_mb(k_AtNr+iat1-1) 155 status=geom_cent_get(geom,iat,at_tag, 156 & dbl_mb(k_xyzpt+3*(iat-1)), 157 & dbl_mb(k_zanpt+iat-1)) 158 xyz_EFGcoords(1)= dbl_mb(k_xyzpt +3*(iat-1)) 159 xyz_EFGcoords(2)= dbl_mb(k_xyzpt+1+3*(iat-1)) 160 xyz_EFGcoords(3)= dbl_mb(k_xyzpt+2+3*(iat-1)) 161 if (ga_nodeid().eq.0) then 162 write(*,19) iat,xyz_EFGcoords(1),xyz_EFGcoords(2), 163 & xyz_EFGcoords(3) 164 19 format('xyz_EFG(',i2,')=(',f15.8,',',f15.8,',', 165 & f15.8,')') 166 endif 167 do zora_Qpq=1,6 ! xx,yy,zz,xy,xz,yz 168c------Generate munu A^{pq}_r ----- START 169 do i=1,ipol 170 call ga_zero(g_efgz4(i)) 171 enddo 172 call zora_getv_EFGZ4_SR(rtdb,g_dens_ini, 173 & zora_calc_type, 174 & zora_Qpq,xyz_EFGcoords, 175 & g_efgz4, ! out: munu matrix 176 & nexc) 177 zora_eint=0.0d0 178 do ispin=1,ipol 179 toac=ga_ddot(g_densZ4(ispin),g_efgz4(ispin)) 180 zora_eint=zora_eint+toac 181 if (ga_nodeid().eq.0) then 182 write(*,15) zora_calc_type,iat,zora_Qpq, 183 & ispin,pos_dat,toac,zora_eint 184 15 format('zora-efg(',i3,',',i3,',',i3,',',i3,',',i3,')=(', 185 & f15.8,',',f15.8,')') 186 endif 187 end do ! ispin-loop 188c ++++++ for NLMO analysis ++++++++++++++++++ START 189 if (zora_calc_type.eq.3 .and. efgfile.eq.1) then 190 call get_munu_rhos_symm(g_efgz4, ! input 191 & ipol, 192 & l_ilst,k_ilst, 193 & l_jlst,k_jlst, 194 & nlst, 195 & g_munu_rhoS, ! output 196 & pos) 197 endif 198 pos=pos+1 199c ++++++ for NLMO analysis ++++++++++++++++++ END 200 call ga_fill_patch(g_rhoS,1,1,pos_dat,pos_dat, 201 & zora_eint) ! zora_eint --> g_rhoS(pos_dat) 202c------Generate munu A^{pq}_r ----- END 203 pos_dat=pos_dat+1 204 end do ! zora_Qpq loop 205 end do ! iat loop 206c ------------- destroy g_densZ4() ------------ START 207 if (zora_calc_type.eq.4) then 208 do i=1,3 209 if (.not. ga_destroy(g_densZ4(i))) call errquit( 210 & 'dft_zora_rhos: ga_destroy failed ',0, GA_ERR) 211 enddo 212 endif 213c ------------- destroy g_densZ4() ------------ END 214 end do ! zora_calc_type loop 215c ----- Store efgz4 data in a file ------- START 216c Note.- lbl_efgz4 defined in zora.fh 217 call util_file_name(lbl_efgz4,.false.,.false.,zorafilename) 218 if (.not.dft_zoraEFGZ4_write( 219 & zorafilename, 220 & nlist, 221 & nat, 222 & g_AtNr, 223 & g_rhoS)) 224 & call errquit('get_rhoS: dft_zoraNMR_write failed', 225 & 0,DISK_ERR) 226c ----- Store efgz4 data in a file ------- END 227 efgfile=0 ! not doing NLMO analysis by default 228 status=rtdb_get(rtdb,'prop:efgfile',mt_int,1,efgfile) ! for NLMO analysis 229 if (efgfile.eq.1) then ! ============ efgfile-if---START 230c if (ga_nodeid().eq.0) 231c & write(*,*) '--------- g_munu_rhoS ---------- START' 232c call ga_print(g_munu_rhoS) 233c if (ga_nodeid().eq.0) 234c & write(*,*) '--------- g_munu_rhoS ---------- END' 235c ---------> Write NMLO analysis data: 2 of 3 ----- START 236 call util_file_name(lbl_nlmo,.false.,.false.,zorafilename) 237 if (.not.dft_zoraEFGZ4_NLMOAnalysis_write( 238 & zorafilename, ! in: filename 239 & nbf, ! in: nr basis functions 240 & ndir, ! in: nr of directions: 6 = xx yy zz xy xz yz 241 & nlist, ! in: list of selected atoms 242 & 2, ! in: writing order =1,2,3 243 & ipolmunu, ! in: write for ndata=1 244 & g_zora_scale_munu, ! in: write for ndata=1 245 & g_munu_rhoS, ! in: write for ndata=2 246 & g_densZ4(3), ! in: write for ndata=2 247 & g_munuV6))! in: write for ndata=3 248 & call errquit('get_rhoS: dft_zoraNLMO_write failed', 249 & 0,DISK_ERR) 250 if (.not. ga_destroy(g_munu_rhoS)) call errquit( 251 & 'dft_zora_utils: ga_destroy failed ',0, GA_ERR) 252c ---------> Write NMLO analysis data: 2 of 3 ----- END 253 if (.not.ma_free_heap(l_ilst)) call 254 & errquit('dft_zora_rhos: ma_free_heap l_ilst',0, MA_ERR) 255 if (.not.ma_free_heap(l_jlst)) call 256 & errquit('dft_zora_rhos: ma_free_heap l_jlst',0, MA_ERR) 257 endif! ============================= efgfile-if---END 258c ------ destroy remaining from zora_calc_type=3 ------ START 259 do i=1,3 260 if (.not. ga_destroy(g_densZ4(i))) call errquit( 261 & 'dft_zora_rhos: ga_destroy failed ',0, GA_ERR) 262 enddo 263c ------ destroy remaining from zora_calc_type=3 ------ END 264c----deallocate memory 265 if (.not.ma_free_heap(l_zanpt)) call errquit 266 & ('dft_zora_utils, ma_free_heap of l_zanpt failed',911,MA_ERR) 267 if (.not.ma_free_heap(l_xyzpt)) call errquit 268 & ('dft_zora_utils, ma_free_heap of l_xyzpt failed',911,MA_ERR) 269 if (.not.ma_free_heap(l_AtNr)) call 270 & errquit('dft_zora_utils: ma_free_heap l_AtNr',0, MA_ERR) 271 do i=1,ipol 272 if (.not. ga_destroy(g_efgz4(i))) call errquit( 273 & 'dft_zora_rhos: ga_destroy failed ',0, GA_ERR) 274 enddo 275 if (.not. ga_destroy(g_rhoS)) call errquit( 276 & 'dft_zora_rhos: ga_destroy failed ',0, GA_ERR) 277 return 278 end 279 280 subroutine get_gdensZ4_ith(basis,geom,iorb,nbf,ispin, 281 & nocc,scftyp_dbl) 282c Purpose: Calculate Cmunu_iorb matrix 283c iorb, index for MO 284c mu,nu indices for AO or Basis Set (BS) 285c Author : Fredy Aquino 286#include "nwc_const.fh" 287#include "errquit.fh" 288#include "global.fh" 289#include "bas.fh" 290#include "mafdecls.fh" 291#include "geom.fh" 292#include "stdio.fh" 293#include "rtdb.fh" 294 integer basis ! [Input] Basis set 295 integer geom ! [Input] Geometry 296 integer ispin,iorb,nbf 297 integer nocc 298 integer g_orb,prpvectors(2) 299 integer k_prpocc,l_prpocc 300 double precision scftyp_dbl 301 integer g_dens_sf 302 integer ga_create_atom_blocked 303 external ga_create_atom_blocked 304 common /Cmunu/g_dens_sf,prpvectors 305c --- Cmunu --> g_dens_sf 306 if (.not.ma_push_get(mt_dbl,nbf*2,'MO occ',l_prpocc,k_prpocc)) 307 & call errquit('get_gdensZ4_ith:ma_push_get l_prpocc', 308 & 0,MA_ERR) 309 g_orb = ga_create_atom_blocked(geom,basis,'orbs') 310 call ga_get(prpvectors(ispin),1,nbf,iorb,iorb, 311 & dbl_mb(k_prpocc),1) 312 call ga_zero(g_orb) 313 call ga_put(g_orb,1,nbf,iorb,iorb, 314 & dbl_mb(k_prpocc),1) 315 call ga_zero(g_dens_sf) 316 call ga_dgemm('n','t',nbf,nbf,nocc,scftyp_dbl, 317 & g_orb,g_orb,0.d00,g_dens_sf) 318 if (.not.ma_pop_stack(l_prpocc)) call 319 & errquit('get_gdensZ4_ith: ma_pop_stack l_occ',0, MA_ERR) 320 if (.not. ga_destroy(g_orb)) call errquit( 321 & 'get_gdensZ4_ith: ga_destroy failed ',0, GA_ERR) 322 return 323 end 324c ------ calculate get_gdensZ4_jith -------------------- START 325c --> To be used in evaluation of F_{ji}^{1k} NMRZ4 326 subroutine get_gdensZ4_jith(basis,geom, 327 & jorb,iorb,nbf,ispin, 328 & nocc,scftyp_dbl) 329c Purpose: Calculate Cmunu_jiorb matrix 330c jorb, can be index of virtual MO 331c iorb, can be index of occupied MO 332c mu,nu indices for AO or Basis Set (BS) 333c Output : g_dens_sf 334c Author : Fredy Aquino 335#include "nwc_const.fh" 336#include "errquit.fh" 337#include "global.fh" 338#include "bas.fh" 339#include "mafdecls.fh" 340#include "geom.fh" 341#include "stdio.fh" 342#include "rtdb.fh" 343 integer basis ! [Input] Basis set 344 integer geom ! [Input] Geometry 345 integer ispin,iorb,jorb,nbf 346 integer nocc 347 integer g_orbj,gorbi,g_dens_sf,prpvectors(2) 348 integer k_prpocc,l_prpocc 349 double precision scftyp_dbl 350 integer ga_create_atom_blocked 351 external ga_create_atom_blocked 352 common /Cmunu/g_dens_sf,prpvectors 353c --- Cmunu --> g_dens_sf 354 if (.not.ma_push_get(mt_dbl,nbf*2,'MO occ',l_prpocc,k_prpocc)) 355 & call errquit('get_gdensZ4_jith:ma_push_get l_prpocc', 356 & 0,MA_ERR) 357 g_orbj = ga_create_atom_blocked(geom,basis,'orbs') 358 g_orbi = ga_create_atom_blocked(geom,basis,'orbs') 359c ---- Getting Cj (g_orbj) MO coeffs 360 call ga_get(prpvectors(ispin),1,nbf,jorb,jorb, 361 & dbl_mb(k_prpocc),1) 362 call ga_zero(g_orbj) 363 call ga_put(g_orbj,1,nbf,jorb,jorb, 364 & dbl_mb(k_prpocc),1) 365c ---- Getting Ci (g_orbi) MO coeffs 366 call ga_get(prpvectors(ispin),1,nbf,iorb,iorb, 367 & dbl_mb(k_prpocc),1) 368 call ga_zero(g_orbi) 369 call ga_put(g_orbi,1,nbf,iorb,iorb, 370 & dbl_mb(k_prpocc),1) 371 call ga_zero(g_dens_sf) 372 call ga_dgemm('n','t',nbf,nbf,nocc,scftyp_dbl, 373 & g_orbj,g_orbi,0.d00,g_dens_sf) 374 if (.not.ma_pop_stack(l_prpocc)) call 375 & errquit('get_gdensZ4_jith: ma_pop_stack l_occ',0, MA_ERR) 376 if (.not. ga_destroy(g_orbi)) call errquit( 377 & 'get_gdensZ4_jith: ga_destroy failed ',0, GA_ERR) 378 if (.not. ga_destroy(g_orbj)) call errquit( 379 & 'get_gdensZ4_jith: ga_destroy failed ',0, GA_ERR) 380 return 381 end 382c ------ calculate get_gdensZ4_jith -------------------- END 383 384 subroutine hnd_prp_get_vecs(rtdb,geom,basis, 385 2 prpvectors, 386 2 scftyp, 387 2 nclosed,nopen,nvirt) 388 implicit none 389#include "errquit.fh" 390#include "mafdecls.fh" 391#include "global.fh" 392#include "bas.fh" 393#include "util.fh" 394c 395c Assumes energy has been completed, MO vectors stored 396c and all information is still in the RTDB 397c 398 integer rtdb ! [input] database handle 399 integer geom ! [input] geometry handle 400 integer basis ! [input] handles to basis 401 integer ndens ! [output] number of active density handles (RHF=1, UHF=3) 402 character*3 scftyp ! [output] type of wave function 403 integer nclosed(2),nopen(2),nvirt(2) ! [output] occupation info 404c 405 integer nbf, nmo, k_prpocc, l_prpocc, k_prpeval, l_prpeval 406 integer prpvectors(2) 407c 408c Get vectors and other information 409c Arrays occ(nbf*2) and evals(nbf*2) are needed 410c 411 if (.not. bas_numbf(basis,nbf)) call 412 & errquit('hnd_prp_get_dens: could not get nbf',0, BASIS_ERR) 413 if (.not.ma_push_get(MT_DBL,nbf*2,'MO eval',l_prpeval,k_prpeval)) 414 & call errquit('hnd_prp_get_dens:ma_push_get l_prpeval',0,MA_ERR) 415 if (.not.ma_push_get(MT_DBL,nbf*2,'MO occ',l_prpocc,k_prpocc)) 416 & call errquit('hnd_prp_get_dens:ma_push_get l_prpocc',0,MA_ERR) 417 call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen,nvirt, 418 & scftyp,prpvectors,dbl_mb(k_prpocc), 419 & dbl_mb(k_prpeval),nmo) 420c 421c Make the density matrix 422c 423c Cleanup of MA arrays that are not needed 424c 425 if (.not.ma_pop_stack(l_prpocc)) call 426 & errquit('hnd_prp_get_dens: ma_pop_stack l_occ',0, MA_ERR) 427 if (.not.ma_pop_stack(l_prpeval)) call 428 & errquit('hnd_prp_get_dens: ma_pop_stack l_eval',0, MA_ERR) 429c 430 return 431 end 432 433 subroutine get_densZ4(rtdb,basis,geom,g_densZ4) 434c 435c Purpose: Calculate Cmunu Z4 density matrix 436c mu,nu indices for AO or Basis Set (BS) 437c Author : Fredy Aquino 438 439 implicit none 440#include "nwc_const.fh" 441#include "errquit.fh" 442#include "global.fh" 443#include "bas.fh" 444#include "mafdecls.fh" 445#include "geom.fh" 446#include "stdio.fh" 447#include "rtdb.fh" 448#include "zora.fh" 449 character*3 scftyp ! output from hnd_prp_get_vecs() 450 integer rtdb 451 integer basis ! [Input] Basis set 452 integer geom ! [Input] Geometry 453 integer ispin,ipol,iorb,porb,nbf 454 integer nocc(2) 455 integer prpvectors(2) 456 integer g_dens_sf,g_densZ4(3) 457 integer l_Ci,k_Ci 458 integer noc1 ! Total # of occupied MOs 459 integer nclosed(2),nopen(2),nvirt(2) ! [output] occupation info 460 double precision scftyp_dbl,scale 461 integer ga_create_atom_blocked 462 external ga_create_atom_blocked 463 integer i,ii 464 465 common /Cmunu/g_dens_sf,prpvectors 466 467 if (.not. bas_numbf(basis,nbf)) call 468 & errquit('get_densZ4: could not get nbf',0, BASIS_ERR) 469 g_dens_sf=ga_create_atom_blocked(geom,basis,'sf-1') 470 do ii=1,3 471 g_densZ4(ii)=ga_create_atom_blocked(geom,basis,'densZ4-1') 472 end do 473 call hnd_prp_get_vecs(rtdb,geom,basis,prpvectors,scftyp, 474 & nclosed,nopen,nvirt) 475 if (scftyp.eq.'RHF') then 476 ipol=1 477 scftyp_dbl=2.0d00 478 nocc(1)=nclosed(1)+nopen(1) ! is that right?? 479 noc1=nocc(1) 480 else if (scftyp.eq.'UHF') then 481 ipol=2 482 scftyp_dbl=1.0d00 483 nocc(1)=nopen(1) 484 nocc(2)=nopen(2) 485 noc1=nocc(1)+nocc(2) 486 endif 487 if (.not.ma_alloc_get(mt_dbl,ipol*noc1,'Ci',l_Ci,k_Ci)) 488 & call errquit('hnd_efgmap_Z4: ma failed',0,MA_ERR) 489 call ga_get(g_Ci,1,ipol,1,noc1,dbl_mb(k_Ci),ipol) 490 do ispin=1,ipol 491 porb=ispin 492 call ga_zero(g_densZ4(ispin)) 493 do iorb=1,nocc(ispin) 494 call get_gdensZ4_ith(basis,geom,iorb,nbf,ispin, 495 & nocc(ispin),scftyp_dbl) 496 if ((do_NonRel) .or. (not_zora_scale)) then 497 scale=1.0d0 498 else 499 scale=1.0d0/dbl_mb(k_Ci+porb-1) 500 endif 501 call ga_scale(g_dens_sf,scale) 502 call ga_add(1.0d00,g_densZ4(ispin), 503 & 1.0d00,g_dens_sf,g_densZ4(ispin)) 504 porb=porb+ipol 505 end do ! iorb 506 end do ! ispin 507 call ga_zero(g_densZ4(3)) 508 call ga_copy(g_densZ4(1),g_densZ4(3)) 509 if (ipol.gt.1) then 510 call ga_add(1.0d00,g_densZ4(1), 511 & 1.0d00,g_densZ4(2), 512 & g_densZ4(3)) 513 end if 514 if (.not. ga_destroy(g_dens_sf)) call errquit( 515 & 'get_densZ4: ga_destroy failed ',0, GA_ERR) 516 if (.not.ga_destroy(prpvectors(1))) call 517 & errquit('get_densZ4: ga_destroy vecs 1',0, GA_ERR) 518 if (scftyp.eq.'UHF') then 519 if (.not.ga_destroy(prpvectors(2))) call 520 & errquit('get_densZ4: ga_destroy vecs 2',0, GA_ERR) 521 endif 522 if (.not. MA_free_heap(l_Ci)) 523 & call errquit('get_densZ4:cannot free heap',111, MA_ERR) 524 return 525 end 526 527 integer function read_SLCTD_EFG_Atoms 528 & (rtdb,nat,nlist,g_AtNr) 529c---- GA output: g_AtNr 530 531 implicit none 532#include "errquit.fh" 533#include "global.fh" 534#include "tcgmsg.fh" 535#include "msgtypesf.h" 536#include "mafdecls.fh" 537#include "msgids.fh" 538#include "cscfps.fh" 539#include "inp.fh" 540#include "util.fh" 541#include "stdio.fh" 542#include "rtdb.fh" 543 544#include "context.fh" 545 546 integer rtdb,ii,nlist,nat,g_AtNr 547 integer atomnr(nat) 548 integer efgz4atoms 549 double precision AtNr_dbl 550 551 read_SLCTD_EFG_Atoms = 0 552 if (.not. rtdb_get(rtdb, 'efgz4:natoms',mt_int, 553 $ 1,efgz4atoms)) 554 & efgz4atoms=0 ! reset 555 if (efgz4atoms.eq.0) then 556 efgz4atoms=nat 557 nlist=efgz4atoms 558 do ii=1,efgz4atoms 559 AtNr_dbl=ii 560 call ga_put(g_AtNr,1,1,ii,ii,AtNr_dbl,1) 561 enddo 562 else 563 if (.not. rtdb_get(rtdb, 'efgz4:atom list',mt_int, 564 $ efgz4atoms,atomnr)) 565 $ call errquit('prop_input-EFGZ4: rtdb_get failed', 566 $ 555, RTDB_ERR) 567 nlist=efgz4atoms 568 do ii=1,efgz4atoms 569 AtNr_dbl=atomnr(ii) 570 call ga_put(g_AtNr,1,1,ii,ii,AtNr_dbl,1) 571 enddo 572 endif 573 read_SLCTD_EFG_Atoms = 1 574 return 575 end 576 subroutine get_munu_rhos_symm(g_munu, ! in : munu matrix 577 & ipol, ! in : nr. polarizations 578 & l_ilst,k_ilst, ! in : array of indices ilst 579 & l_jlst,k_jlst, ! in : array of indices jlst 580 & nlst, ! in : = nbf*(nbf+1)/2 581 & g_munu_rhoS, ! out : accumulate main-diag,off-diag munu elements 582 & pos) 583c g_munu_rhoS accumulates xx,yy,zz,xy,xz,yz 584c matrices in the following format: 585c 11 22 ... (nbf nbf) -> main diagonal first 586c 21 31 32 41 42 43 -> off-diagonal later 587c nbf, number of basis functions 588c nlst=nbf*(nbf*+1)/2 589c 1st chunk of nlst numbers corresponds to xx 590c 2nd chunk --> yy, 3rd chunk --> zz, 591c 4th chunk --> xy, 5th chunk --> xz, 6th chunk --> yz 592c Input: pos=0,1,..., chunk-index 593c --> the size of one chunk is nlst 594 implicit none 595#include "errquit.fh" 596#include "mafdecls.fh" 597#include "stdio.fh" 598#include "global.fh" 599#include "msgids.fh" 600#include "rtdb.fh" 601#include "geom.fh" 602#include "zora.fh" 603 604 integer g_munu_rhoS 605 integer g_munu(2) 606 integer i,ipol 607 integer nlst,pos 608 integer l_ilst,k_ilst, 609 & l_jlst,k_jlst 610 integer jlo,jhi 611 double precision v(nlst) 612 call ga_gather(g_munu(1),v, ! g_munu ---> v 613 & int_mb(k_ilst),int_mb(k_jlst), 614 & nlst) 615c Now accumulates ith-chunk on g_munu_rhoS 616 jlo=pos*nlst+1 617 jhi=pos*nlst+nlst 618 call ga_put(g_munu_rhoS,1,1,jlo,jhi,v,1) ! v --> g_munu_rhoS 619 return 620 end 621c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 622c +++++++++++ READ/WRITE EFGZ4 data +++++++++++++ START 623c Note.- Using modified versions of 624c dft_zora_read() and dft_zoraNMR_write() 625c --> located in dft_zora_utils.F 626czora...Write out the zora NMR shieldings to disk 627 628 logical function dft_zoraEFGZ4_write( 629 & filename, ! in: filename 630 & nlist, ! in: number of selected atoms 631 & nat, ! in: total number of atoms 632 & g_AtNr1, ! in: list of atoms to calc. shieldings 633 & g_rhoS) ! in: rhoS values 634 implicit none 635#include "errquit.fh" 636#include "mafdecls.fh" 637#include "global.fh" 638#include "tcgmsg.fh" 639#include "msgtypesf.h" 640#include "inp.fh" 641#include "msgids.fh" 642#include "cscfps.fh" 643#include "util.fh" 644#include "stdio.fh" 645 character*(*) filename ! [input] File to write to 646 integer nlist, ! nr. slc atoms 647 & nat, ! total nr. of atoms 648 & ntot_efg ! = ntat*6*2 6=xx,yy,zz,xy,xz,yz 2=NonRel, EFGZ4 649 integer g_AtNr1, ! in: list of atoms to calc. shieldings 650 & g_rhoS ! rhoS values 651 integer unitno 652 parameter (unitno = 77) 653 integer l_AtNr,k_AtNr, 654 & l_rhoS,k_rhoS 655 integer ok, iset, i, j 656 integer inntsize 657 integer alo(3),ahi(3),ld(2) 658 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 659 call ga_sync() 660 ok = 0 661c Read routines should be consistent with this 662c Write out the atomic zora corrections 663 if (ga_nodeid() .eq. 0) then 664c Open the file 665 open(unitno, status='unknown', form='unformatted', 666 $ file=filename, err=1000) 667c Write out the number of sets and basis functions 668 write(unitno, err=1001) nlist 669 write(unitno, err=1001) nat 670c Allocate the temporary buffer 671c ++++++++ using ma_alloc_get +++++++++++++++++++ START 672c ---> ma_alloc_get: to allocate memory 673c ---> ma_free_heap: to release allocated memory 674 if (.not. ma_alloc_get(mt_dbl,nlist,'dft_zoraEFGZ4_write', 675 & l_AtNr,k_AtNr)) 676 $ call errquit('dft_zoraEFGZ4_write: ma failed', 677 & nlist, MA_ERR) 678 ntot_efg=nlist*6*2 679 if (.not. ma_alloc_get(mt_dbl,ntot_efg, 680 & 'dft_zoraEFGZ4_write', 681 & l_rhoS,k_rhoS)) 682 $ call errquit('dft_zoraEFGZ4_write: ma failed', 683 & ntot_efg, MA_ERR) 684c ++++++++ using ma_alloc_get +++++++++++++++++++ END 685 call ga_get(g_AtNr1,1,1,1,nlist,dbl_mb(k_AtNr),1) 686 call swrite(unitno,dbl_mb(k_AtNr),nlist) 687 call ga_get(g_rhoS,1,1,1,ntot_efg,dbl_mb(k_rhoS),1) 688 call swrite(unitno,dbl_mb(k_rhoS),ntot_efg) 689c 690c Deallocate the temporary buffer 691c ----- Using ma_free_heap ------------ START 692 if (.not. ma_free_heap(l_AtNr)) 693 $ call errquit('dft_zoraEFGZ4_write: ma free_heap failed', 694 & 911, MA_ERR) 695 if (.not. ma_free_heap(l_rhoS)) 696 $ call errquit('dft_zoraEFGZ4_write: ma free_heap failed', 697 & 911, MA_ERR) 698c ----- Using ma_free_heap ------------ END 699c Close the file 700 close(unitno,err=1002) 701 ok = 1 702 end if 703c Broadcast status to other nodes 704 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 705 call ga_sync() 706 dft_zoraEFGZ4_write = (ok .eq. 1) 707 if (ga_nodeid() .eq. 0) then 708 write(6,22) filename(1:inp_strlen(filename)) 709 22 format(/' Wrote ZORA EFGZ4 data to ',a/) 710 call util_flush(luout) 711 endif 712 return 713 1000 write(6,*) 'dft_zoraEFGZ4_write: failed to open ', 714 $ filename(1:inp_strlen(filename)) 715 call util_flush(luout) 716 ok = 0 717 goto 10 718 1001 write(6,*) 'dft_zoraEFGZ4_write: failed to write ', 719 $ filename(1:inp_strlen(filename)) 720 call util_flush(luout) 721 ok = 0 722 close(unitno,err=1002) 723 goto 10 724 1002 write(6,*) 'dft_zoraEFGZ4_write: failed to close', 725 $ filename(1:inp_strlen(filename)) 726 call util_flush(luout) 727 ok = 0 728 goto 10 729 end 730 731 logical function dft_zoraEFGZ4_read( 732 & filename, ! in : filename 733 & nlist, ! in : number of selected atoms 734 & nat, ! in : total number of atoms 735 & g_AtNr1, ! out: list of atoms to calc. shieldings 736 & g_rhoS) ! out: rhoS values 737 implicit none 738#include "errquit.fh" 739#include "global.fh" 740#include "tcgmsg.fh" 741#include "msgtypesf.h" 742#include "mafdecls.fh" 743#include "msgids.fh" 744#include "cscfps.fh" 745#include "inp.fh" 746#include "util.fh" 747#include "stdio.fh" 748 character*(*) filename ! [input] File to write to 749 integer nlist, ! nr. slc atoms 750 & nat, ! total nr. of atoms 751 & ntot_efg ! = ntat*6*2 6=xx,yy,zz,xy,xz,yz 2=NonRel, EFGZ4 752 integer g_AtNr1, ! in: list of atoms to calc. shieldings 753 & g_rhoS 754 integer unitno 755 parameter (unitno = 77) 756 integer l_AtNr,k_AtNr, 757 & l_rhoS,k_rhoS 758 integer ok,inntsize 759 integer nat_read,nlist_read 760c Initialise to invalid MA handle 761 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 762 call ga_sync() 763 ok = 0 764c---- Create GA arrays: g_AtNr1,g_rhoS --- START 765 if (.not. ga_create(mt_dbl,1,nlist, 766 & 'dft_zoraEFGZ4_read: g_AtNr1',0,0,g_AtNr1)) 767 $ call errquit('dft_zoraEFGZ4_read: g_AtNr1',0,GA_ERR) 768 call ga_zero(g_AtNr1) 769 ntot_efg=nlist*6*2 770 if (.not. ga_create(mt_dbl,1,ntot_efg, 771 & 'dft_zoraEFGZ4_read: g_rhoS',0,0,g_rhoS)) 772 $ call errquit('dft_zoraEFGZ4_read: g_rhoS',0,GA_ERR) 773 call ga_zero(g_rhoS) 774c---- Create GA arrays: g_AtNr1,g_rhoS --- END 775 if (ga_nodeid() .eq. 0) then 776c Print a message indicating the file being read 777 write(6,22) filename(1:inp_strlen(filename)) 778 22 format(/' Read ZORA EFGZ4 data from ',a/) 779 call util_flush(luout) 780c Open the file 781 open(unitno, status='old', form='unformatted', file=filename, 782 $ err=1000) 783c Read in some basics to check if they are consistent with the calculation 784 read(unitno, err=1001, end=1001) nlist_read 785 read(unitno, err=1001, end=1001) nat_read 786c Error checks 787 if ((nat_read .ne. nat) .or. 788 & (nlist_read .ne. nlist) ) goto 1003 789c ++++++++ using ma_alloc_get +++++++++++++++++++ START 790c ---> ma_alloc_get: to allocate memory 791c ---> ma_free_heap: to release allocated memory 792 if (.not. ma_alloc_get(mt_dbl,nlist,'dft_zoraEFGZ4_read', 793 & l_AtNr,k_AtNr)) 794 $ call errquit('dft_zoraEFGZ4_read: ma failed', 795 & nlist, MA_ERR) 796 if (.not. ma_alloc_get(mt_dbl,ntot_efg,'dft_zoraEFGZ4_read', 797 & l_rhoS,k_rhoS)) 798 $ call errquit('dft_zoraEFGZ4_read: ma failed', 799 & ntot_efg, MA_ERR) 800c ++++++++ using ma_alloc_get +++++++++++++++++++ END 801 call sread(unitno,dbl_mb(k_AtNr),nlist) 802 call ga_put(g_AtNr1,1,1,1,nlist,dbl_mb(k_AtNr),1) 803 call sread(unitno,dbl_mb(k_rhoS),ntot_efg) 804 call ga_put(g_rhoS,1,1,1,ntot_efg,dbl_mb(k_rhoS),1) 805c 806c Deallocate the temporary buffer 807c ----- Using ma_free_heap ------------ START 808 if (.not. ma_free_heap(l_AtNr)) 809 $ call errquit('dft_zoraEFGZ4_read: ma free_heap failed', 810 & 911, MA_ERR) 811 if (.not. ma_free_heap(l_rhoS)) 812 $ call errquit('dft_zoraEFGZ4_read: ma free_heap failed', 813 & 911, MA_ERR) 814c ----- Using ma_free_heap ------------ END 815c Close the file 816 close(unitno,err=1002) 817 ok = 1 818 end if 819c 820c Broadcast status to other nodes 821 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 822 call ga_sync() 823 dft_zoraEFGZ4_read = ok .eq. 1 824 return 825 1000 write(6,*) 'dft_zoraEFGZ4_read: failed to open', 826 $ filename(1:inp_strlen(filename)) 827 call util_flush(luout) 828 ok = 0 829 goto 10 830 1001 write(6,*) 'dft_zoraEFGZ4_read: failed to read', 831 $ filename(1:inp_strlen(filename)) 832 call util_flush(luout) 833 ok = 0 834 close(unitno,err=1002) 835 goto 10 836 1003 write(6,*) 'dft_zoraEFGZ4_read: file inconsistent', 837 & ' with calculation', 838 $ filename(1:inp_strlen(filename)) 839 call util_flush(luout) 840 ok = 0 841 close(unitno,err=1002) 842 goto 10 843 1002 write(6,*) 'dft_zoraEFGZ4_read: failed to close', 844 $ filename(1:inp_strlen(filename)) 845 call util_flush(luout) 846 ok = 0 847 goto 10 848 end 849 850 logical function dft_zoraEFGZ4_NLMOAnalysis_write( 851 & filename, ! in: filename 852 & nbf, ! in: nr basis functions 853 & ndir, ! in: nr of directions: 6 = xx yy zz xy xz yz 854 & nlist, ! in: list of selected atoms 855 & ndata, ! in: writing order =1,2,3 856 & ipolmunu, ! in: write for ndata=1 857 & g_zora_scale_munu, ! in: write for ndata=1 858 & g_munu_rhoS, ! in: write for ndata=2 859 & g_dens, ! in: write for ndata=2 860 & g_munuV6) ! in: write for ndata=3 861c Description: Collecting three matrices to be used 862c in wefgfile(rtdb) and wnbofile(rtdb) 863c Those routines are called in prop.F 864c after hnd_property(rtdb) 865c The info collected is: 866c 1. ipolmunu, nr. of polarizations 867c 2. g_zora_scale_munu, = g_zora_scale_sf 868c --> Stored in dft_zora_scale() [dft_zora_utils.F] 869c 3. g_munu_rhoS, small component density munu matrix 870c about the format check info in header 871c of get_munu_rhos_symm() 872c --> Stored in get_rhoS() [dft_zora_rhos.F] 873c 4. g_munuV6, electronic EFG munu matrix (without the 874c small component density correction) 875c --> Stored in get_munuV6() [hnd_elfcon_symm.F] 876 implicit none 877#include "errquit.fh" 878#include "mafdecls.fh" 879#include "global.fh" 880#include "tcgmsg.fh" 881#include "msgtypesf.h" 882#include "inp.fh" 883#include "msgids.fh" 884#include "cscfps.fh" 885#include "util.fh" 886#include "stdio.fh" 887 character*(*) filename ! [input] File to write to 888 integer nlist, ! = nr. slc atoms 889 & nlst, ! = nbf*(nbf+1)/2 890 & nmat, ! = nlst*ndir*nlist 891 & nmat1 ! = nbf*nbf 892 integer ndata ! = 1,2,3 893 integer ipolmunu,g_zora_scale_munu(2) ! write for ndata=1 894 integer g_munu_rhoS,g_dens ! write for ndata=2 895 integer g_munuV6 ! write for ndata=3 896 integer ndir,nbf 897 integer unitno 898 parameter (unitno = 77) 899 integer l_mat,k_mat,l_mat1,k_mat1 900 integer ok, iset, i, j 901 integer inntsize 902c WARNING: The sequence in which is called this routine is: 1->2->3. 903c ---- check ndata=1,2,3 ----- START 904 if (ndata.ne.1 .and. 905 & ndata.ne.2 .and. ndata.ne.3) then 906 write(*,*) 'Error in dft_zoraEFGZ4_NLMOAnalysis_write:', 907 & 'ndata=',ndata,' It should be 1, 2 or 3' 908 stop 909 endif 910c ---- check ndata=1,2,3 ----- END 911 nlst=nbf*(nbf+1)/2 912 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 913 call ga_sync() 914 ok = 0 915c Read routines should be consistent with this 916c Write out the atomic zora corrections 917 if (ga_nodeid() .eq. 0) then 918c Allocate the temporary buffer 919c ++++++++ using ma_alloc_get +++++++++++++++++++ START 920c ---> ma_alloc_get: to allocate memory 921c ---> ma_free_heap: to release allocated memory 922 if (ndata.eq.1) then 923 nmat=nbf 924 if (.not. ma_alloc_get( 925 & mt_dbl,nmat,'dft_zoraNLMO_write', 926 & l_mat,k_mat)) 927 $ call errquit('dft_zoraNLMO_write: ma failed', 928 & nmat, MA_ERR) 929 else if (ndata.eq.2) then 930 nmat=nlst*ndir*nlist 931 nmat1=nbf*nbf 932 if (.not. ma_alloc_get( 933 & mt_dbl,nmat,'dft_zoraNLMO_write', 934 & l_mat,k_mat)) 935 $ call errquit('dft_zoraNLMO_write: ma failed', 936 & nmat, MA_ERR) 937 if (.not. ma_alloc_get( 938 & mt_dbl,nmat1,'dft_zoraNLMO_write', 939 & l_mat1,k_mat1)) 940 $ call errquit('dft_zoraNLMO_write: ma failed', 941 & nmat1, MA_ERR) 942 else if (ndata.eq.3) then 943 nmat=nlst*ndir*nlist 944 if (.not. ma_alloc_get( 945 & mt_dbl,nmat,'dft_zoraNLMO_write', 946 & l_mat,k_mat)) 947 $ call errquit('dft_zoraNLMO_write: ma failed', 948 & nmat, MA_ERR) 949 endif 950c ++++++++ using ma_alloc_get +++++++++++++++++++ END 951 if (ndata.eq.1) then 952c Open the file - 1st time 953 open(unitno, status='unknown', form='unformatted', 954 $ file=filename, err=1000) 955c Write out the number of sets and basis functions 956 write(unitno, err=1001) nbf 957 write(unitno, err=1001) nlst 958 write(unitno, err=1001) ndir 959 write(unitno, err=1001) ipolmunu 960c Write out g_zora_scale_sf 961 do iset = 1, ipolmunu 962 do i = 1, nbf 963 call ga_get(g_zora_scale_munu(iset), 1, nbf, i, i, 964 & dbl_mb(k_mat),1) 965 call swrite(unitno, dbl_mb(k_mat), nbf) 966 end do 967 end do 968 else if (ndata.eq.2) then 969c Open the file - 2nd time 970 open(unitno, status='unknown', form='unformatted', 971 $ file=filename, err=1000,position='append') 972 write(unitno, err=1001) nlist 973 call ga_get(g_munu_rhoS,1,1,1,nmat, 974 & dbl_mb(k_mat),1) 975 call swrite(unitno,dbl_mb(k_mat),nmat) 976 call ga_get(g_dens,1,nbf,1,nbf, 977 & dbl_mb(k_mat1),nbf) 978 call swrite(unitno,dbl_mb(k_mat1),nmat1) 979 else if (ndata.eq.3) then 980c Open the file - 3rd time 981 open(unitno, status='unknown', form='unformatted', 982 $ file=filename, err=1000,position='append') 983 call ga_get(g_munuV6,1,1,1,nmat, 984 & dbl_mb(k_mat),1) 985 call swrite(unitno,dbl_mb(k_mat),nmat) 986 endif 987c 988c Deallocate the temporary buffer 989c ----- Using ma_free_heap ------------ START 990 if (ndata.eq.1) then 991 if (.not. ma_free_heap(l_mat)) 992 $ call errquit('dft_zoraNLMO_write: ma free_heap failed', 993 & 911, MA_ERR) 994 else if (ndata.eq.2) then 995 if (.not. ma_free_heap(l_mat)) 996 $ call errquit('dft_zoraNLMO_write: ma free_heap failed', 997 & 911, MA_ERR) 998 if (.not. ma_free_heap(l_mat1)) 999 $ call errquit('dft_zoraNLMO_write: ma free_heap1 failed', 1000 & 911, MA_ERR) 1001 else if (ndata.eq.3) then 1002 if (.not. ma_free_heap(l_mat)) 1003 $ call errquit('dft_zoraNLMO_write: ma free_heap failed', 1004 & 911, MA_ERR) 1005 endif 1006c ----- Using ma_free_heap ------------ END 1007c Close the file 1008 close(unitno,err=1002) 1009 ok = 1 1010 end if 1011c Broadcast status to other nodes 1012 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1013 call ga_sync() 1014 dft_zoraEFGZ4_NLMOAnalysis_write = (ok .eq. 1) 1015 if (ga_nodeid() .eq. 0) then 1016 write(6,22) filename(1:inp_strlen(filename)) 1017 22 format(/' Wrote ZORA NLMO EFGZ4 data to ',a/) 1018 call util_flush(luout) 1019 endif 1020 return 1021 1000 write(6,*) 'dft_zoraNLMO_write: failed to open ', 1022 $ filename(1:inp_strlen(filename)) 1023 call util_flush(luout) 1024 ok = 0 1025 goto 10 1026 1001 write(6,*) 'dft_zoraNLMO_write: failed to write ', 1027 $ filename(1:inp_strlen(filename)) 1028 call util_flush(luout) 1029 ok = 0 1030 close(unitno,err=1002) 1031 goto 10 1032 1002 write(6,*) 'dft_zoraNLMO_write: failed to close', 1033 $ filename(1:inp_strlen(filename)) 1034 call util_flush(luout) 1035 ok = 0 1036 goto 10 1037 end 1038 1039 logical function dft_zoraEFGZ4_NLMOAnalysis_read( 1040 & filename, ! in: filename 1041 & nbf, ! in: nr basis functions 1042 & ndir, ! in: nr of directions: 6 = xx yy zz xy xz yz 1043 & nlist, ! in: list of selected atoms 1044 & ipolmunu, ! in: 1045 & g_zora_scale_munu, ! out: 1046 & g_munu_rhoS, ! out: 1047 & g_dens, ! out: 1048 & g_munuV6) ! out: 1049 implicit none 1050#include "errquit.fh" 1051#include "global.fh" 1052#include "tcgmsg.fh" 1053#include "msgtypesf.h" 1054#include "mafdecls.fh" 1055#include "msgids.fh" 1056#include "cscfps.fh" 1057#include "inp.fh" 1058#include "util.fh" 1059#include "stdio.fh" 1060 1061 character*(*) filename ! [input] File to write to 1062 integer nmat,nmat1 1063 integer g_zora_scale_munu(2), 1064 & g_munu_rhoS, 1065 & g_dens, 1066 & g_munuV6 1067 integer nbf,nbf_read,nlst,nlst_read, 1068 & ndir,ndir_read,nlist,nlist_read, 1069 & ipolmunu,ipolmunu_read 1070 integer unitno 1071 parameter (unitno = 77) 1072 integer l_mat,k_mat, 1073 & l_mat1,k_mat1 1074 integer ok, iset, i, j 1075 integer inntsize 1076c Initialise to invalid MA handle 1077 nlst=nbf*(nbf+1)/2 1078 inntsize=MA_sizeof(MT_INT,1,MT_BYTE) 1079 call ga_sync() 1080 ok = 0 1081c---- Create GA arrays: g_AtNr1,g_rhoS --- START 1082 if (.not. ga_create(mt_dbl,nbf,nbf, 1083 & 'dft_zoraNLMO_read: g_zora_scale_munu', 1084 & 0,0,g_zora_scale_munu(1))) 1085 $ call errquit('dft_zoraNLMO_read: g_rhoS',0,GA_ERR) 1086 call ga_zero(g_zora_scale_munu(1)) 1087 if (ipolmunu.gt.1) then 1088 if (.not. ga_create(mt_dbl,nbf,nbf, 1089 & 'dft_zoraNLMO_read: g_zora_scale_munu', 1090 & 0,0,g_zora_scale_munu(2))) 1091 $ call errquit('dft_zoraNLMO_read: g_rhoS',0,GA_ERR) 1092 call ga_zero(g_zora_scale_munu(2)) 1093 endif 1094 nmat=nlst*ndir*nlist 1095 if (.not. ga_create(mt_dbl,1,nmat, 1096 & 'dft_zoraNLMO_read: g_rhoS',0,0,g_munu_rhoS)) 1097 $ call errquit('dft_zoraNLMO_read: g_rhoS',0,GA_ERR) 1098 call ga_zero(g_munu_rhoS) 1099 if (.not. ga_create(mt_dbl,nbf,nbf, 1100 & 'dft_zoraNLMO_read: g_sdens',0,0,g_dens)) 1101 $ call errquit('dft_zoraNLMO_read: g_dens',0,GA_ERR) 1102 call ga_zero(g_dens) 1103 if (.not. ga_create(mt_dbl,1,nmat, 1104 & 'dft_zoraNLMO_read: g_munuV6',0,0,g_munuV6)) 1105 $ call errquit('dft_zoraNLMO_read: g_munuV6',0,GA_ERR) 1106 call ga_zero(g_munuV6) 1107c---- Create GA arrays: g_AtNr1,g_rhoS --- END 1108 if (ga_nodeid() .eq. 0) then 1109c Print a message indicating the file being read 1110 write(6,22) filename(1:inp_strlen(filename)) 1111 22 format(/' Read ZORA EFGZ4 data from ',a/) 1112 call util_flush(luout) 1113c Open the file 1114 open(unitno, status='old', form='unformatted', file=filename, 1115 $ err=1000) 1116c Read in some basics to check if they are consistent with the calculation 1117 read(unitno, err=1001, end=1001) nbf_read 1118 read(unitno, err=1001, end=1001) nlst_read 1119 read(unitno, err=1001, end=1001) ndir_read 1120 read(unitno, err=1001, end=1001) ipolmunu_read 1121c Error checks 1122 if ((nbf_read .ne. nbf) .or. 1123 & (nlst_read .ne. nlst) .or. 1124 & (ndir_read .ne. ndir) .or. 1125 & (ipolmunu_read .ne. ipolmunu)) goto 1003 1126c ---> ma_alloc_get: to allocate memory 1127c ---> ma_free_heap: to release allocated memory 1128c ------- Reading 1st matrix: g_zora_scale_munu ----- START 1129 nmat=nbf 1130 if (.not. ma_alloc_get(mt_dbl,nmat, ! allocate memory 1131 & 'dft_zoraNLMO_read',l_mat,k_mat)) 1132 $ call errquit('dft_zoraNLMO_read: ma failed', 1133 & nmat, MA_ERR) 1134 do iset = 1, ipolmunu 1135 do i = 1, nbf 1136 call sread(unitno, dbl_mb(k_mat), nbf) 1137 call ga_put(g_zora_scale_munu(iset), 1, nbf, i, i, 1138 & dbl_mb(k_mat), 1) 1139 end do 1140 end do 1141 if (.not. ma_free_heap(l_mat)) ! deallocate memory 1142 $ call errquit('dft_zoraNLMO_read: ma free_heap failed', 1143 & 911, MA_ERR) 1144c ------- Reading 1st matrix: g_zora_scale_munu ----- END 1145c ------- Reading 2nd matrix: g_munu_rhoS ----- START 1146 read(unitno, err=1001, end=1001) nlist_read 1147 if (nlist_read .ne. nlist) goto 1003 1148 nmat=nlst*ndir*nlist 1149 if (.not. ma_alloc_get(mt_dbl,nmat, ! allocate memory 1150 & 'dft_zoraNLMO_read',l_mat,k_mat)) 1151 $ call errquit('dft_zoraNLMO_read: ma failed', 1152 & nmat, MA_ERR) 1153 nmat1=nbf*nbf 1154 if (.not. ma_alloc_get(mt_dbl,nmat1, ! allocate memory 1155 & 'dft_zoraNLMO_read',l_mat1,k_mat1)) 1156 $ call errquit('dft_zoraNLMO_read: ma failed', 1157 & nmat1, MA_ERR) 1158 call sread(unitno,dbl_mb(k_mat),nmat) 1159 call ga_put(g_munu_rhoS,1,1,1,nmat,dbl_mb(k_mat),1) 1160 call sread(unitno,dbl_mb(k_mat1),nmat1) 1161 call ga_put(g_dens,1,nbf,1,nbf,dbl_mb(k_mat1),nbf) 1162 if (.not. ma_free_heap(l_mat)) ! deallocate memory 1163 $ call errquit('dft_zoraNLMO_read: ma free_heap failed', 1164 & 911, MA_ERR) 1165 if (.not. ma_free_heap(l_mat1)) ! deallocate memory 1166 $ call errquit('dft_zoraNLMO_read: ma free_heap failed', 1167 & 911, MA_ERR) 1168c ------- Reading 2nd matrix: g_munu_rhoS ----- END 1169c ------- Reading 3rd matrix: g_munuV6 ----- START 1170 nmat=nlst*ndir*nlist 1171 if (.not. ma_alloc_get(mt_dbl,nmat, ! allocate memory 1172 & 'dft_zoraNLMO_read',l_mat,k_mat)) 1173 $ call errquit('dft_zoraNLMO_read: ma failed', 1174 & nmat, MA_ERR) 1175 call sread(unitno,dbl_mb(k_mat),nmat) 1176 call ga_put(g_munuV6,1,1,1,nmat,dbl_mb(k_mat),1) 1177 if (.not. ma_free_heap(l_mat)) ! deallocate memory 1178 $ call errquit('dft_zoraNLMO_read: ma free_heap failed', 1179 & 911, MA_ERR) 1180c ------- Reading 3rd matrix: g_munuV6 ----- END 1181c Close the file 1182 close(unitno,err=1002) 1183 ok = 1 1184 end if 1185c 1186c Broadcast status to other nodes 1187 10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status 1188 call ga_sync() 1189 dft_zoraEFGZ4_NLMOAnalysis_read = ok .eq. 1 1190 return 1191 1000 write(6,*) 'dft_zoraNLMO_read: failed to open', 1192 $ filename(1:inp_strlen(filename)) 1193 call util_flush(luout) 1194 ok = 0 1195 goto 10 1196 1001 write(6,*) 'dft_zoraNLMO_read: failed to read', 1197 $ filename(1:inp_strlen(filename)) 1198 call util_flush(luout) 1199 ok = 0 1200 close(unitno,err=1002) 1201 goto 10 1202 1003 write(6,*) 'dft_zoraNLMO_read: file inconsistent', 1203 & ' with calculation', 1204 $ filename(1:inp_strlen(filename)) 1205 call util_flush(luout) 1206 ok = 0 1207 close(unitno,err=1002) 1208 goto 10 1209 1002 write(6,*) 'dft_zoraNLMO_read: failed to close', 1210 $ filename(1:inp_strlen(filename)) 1211 call util_flush(luout) 1212 ok = 0 1213 goto 10 1214 end 1215c +++++++++++ READ/WRITE EFGZ4 data +++++++++++++ END 1216c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1217c 1218 subroutine print_EFGZ4_version() 1219c 1220 implicit none 1221c 1222#include "stdio.fh" 1223c 1224 write(LuOut,*) 1225 call util_print_centered(LuOut, 1226 $ 'EPR Z4', 23, .true.) 1227 write(LuOut,*) 1228c 1229 return 1230 end 1231c $Id$ 1232