1c======================================================= 2c================ Fredy Aquino's routines ======== START 3c======================================================= 4 subroutine get_rhoS_so( 5 & rtdb, 6 & g_dens_ini,nexc, 7 & geom, 8 & ao_bas_han, 9 & nbf,nbf_ao,nbf_mo, 10 & g_dens, 11 & g_moso, 12 & noc, 13 & nocc, 14 & ipol) 15c -- Purpose: Calculation of small component density, \rho_S 16c for Electric Field Gradient calculation 17c Source: van Lenthe, et.al.,JCP, V112,N19,Y2000 18c -- Author : Fredy Aquino 10-30-10 19 20 implicit none 21 22#include "errquit.fh" 23#include "mafdecls.fh" 24#include "stdio.fh" 25#include "global.fh" 26#include "msgids.fh" 27#include "rtdb.fh" 28#include "geom.fh" 29#include "zora.fh" 30 31 integer g_dens(2) 32 integer g_moso(2) 33 integer g_dens_ini(2) 34 integer noc,nocc(2) 35 integer geom 36 integer ao_bas_han 37 integer i,pos_dat, 38 & do_zgc_old,nexc 39 logical dft_zoraEFGZ4_write 40 integer ipol,iat,nat,noc1,ac 41 integer l_xyzpt,k_xyzpt, 42 & l_zanpt,k_zanpt 43 logical status 44 character*16 element, at_tag 45 character*255 zorafilename 46 double precision zora_eint,toac(5) 47 character*100 fname_RSLC 48 integer stat_read 49 integer iat1,l_AtNr,k_AtNr 50 integer rtdb 51 integer g_densZ4so(5),g_rhoS 52 integer read_SLCTD_EFG_Atoms 53 external get_densZ4_so,read_SLCTD_EFG_Atoms, 54 & zora_getv_EFGZ4_SO,dft_zoraEFGZ4_write, 55 & util_file_name 56 integer nbf,nbf_ao,nbf_mo 57 integer zora_Qpq 58 double precision xyz_EFGcoords(3) 59 integer g_efgz4_sf(2) 60 integer g_efgz4_so(3) 61c 62c---- Input global variables (defined in zora.fh): 63c 1. zora_calc_type ! =3 -> ZORA4-EFG 64c 2. so_term ! =0 -> ZORA4-spin-free 65c 3. xyz_EFGcoords(i) i=1,2,3 66c 4. zora_Qpq=1,6 : xx,yy,zz,xy,xz,yz 67 if (ipol.eq.1) then 68 noc1=nocc(1) 69 else if (ipol.eq.2) then 70 noc1=nocc(1)+nocc(2) 71 endif 72 status=geom_ncent(geom,nat) 73c----- Allocate memory - FA 74 if (.not. ma_alloc_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 75 & call errquit('get_rhoS: ma failed',911,MA_ERR) 76 if (.not. ma_alloc_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 77 & call errquit('get_rhoS: ma failed',911,MA_ERR) 78c----- Allocate global arrays - FA 79 do i=1,2 ! up,down 80 if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_scale_sf 1',0,0, 81 & g_efgz4_sf(i))) 82 & call errquit('dft_scf_so: error creating g_zora_scale_sf 1',0, 83 & GA_ERR) 84 enddo 85 do i=1,3 ! x,y,z 86 if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_scale_so 1',0,0, 87 & g_efgz4_so(i))) 88 & call errquit('dft_scf_so: error creating g_zora_scale_so 1',0, 89 & GA_ERR) 90 enddo 91c ++++++++++++++++++++++++++++++++++ 92c +++++ Read Atom Nr for EFG calc ++ 93 if (.not. ga_create(mt_dbl,1,nat, 94 & 'get_rhoS: g_AtNr',0,0,g_AtNr)) 95 $ call errquit('get_rhoS: g_AtNr', 0,GA_ERR) 96 call ga_zero(g_AtNr) 97 stat_read=read_SLCTD_EFG_Atoms 98 & (rtdb,nat,nlist,g_AtNr) 99 if (ga_nodeid().eq.0) 100 & write(*,*) "==> nlist=",nlist 101c Allocate memory for l_AtNr,k_AtNr 102 if (.not.ma_alloc_get(mt_dbl,nat,'AtNr', 103 & l_AtNr,k_AtNr)) 104 & call errquit('get_rhoS: ma failed',0,MA_ERR) 105 call ga_get(g_AtNr,1,1,1,nat,dbl_mb(k_AtNr),1) 106c ++++++++++++++++++++++++++++++++++ 107c--- About content of g_rhoS: 108c--- 1st set of nat*6 elements corresponds to ZORA4-EFG 109c--- 2nd set of nat*6 elements corresponds to NUM-EFG 110 if (.not. ga_create(mt_dbl,1,nlist*6*2, 111 & 'get_rhoS: g_rhoS',0,0,g_rhoS)) 112 & call errquit('get_rhoS: g_rhoS', 0,GA_ERR) 113 call ga_zero(g_rhoS) 114 pos_dat=1 115c do zora_calc_type =3,4 ! ZORA4-EFG,NUM-EFG 116 do zora_calc_type =4,3,-1 ! NUM-EFG,ZORA4-EFG 117 call get_densZ4_so(rtdb,geom,ao_bas_han, 118 & nbf,nbf_ao,nbf_mo, 119 & g_dens,g_moso, 120 & noc,nocc,g_densZ4so) 121 do iat1=1,nlist ! nlist<=nat 122 iat=dbl_mb(k_AtNr+iat1-1) 123 status=geom_cent_get(geom,iat,at_tag, 124 & dbl_mb(k_xyzpt+3*(iat-1)), 125 & dbl_mb(k_zanpt+iat-1)) 126 xyz_EFGcoords(1)= dbl_mb(k_xyzpt +3*(iat-1)) 127 xyz_EFGcoords(2)= dbl_mb(k_xyzpt+1+3*(iat-1)) 128 xyz_EFGcoords(3)= dbl_mb(k_xyzpt+2+3*(iat-1)) 129 if (ga_nodeid().eq.0) then 130 write(*,19) iat,xyz_EFGcoords(1),xyz_EFGcoords(2), 131 & xyz_EFGcoords(3) 132 19 format('xyz_EFG(',i2,')=(',f15.8,',',f15.8,',', 133 & f15.8,')') 134 endif 135 do zora_Qpq=1,6 ! xx,yy,zz,xy,xz,yz - FA 136c------Generate munu A^{pq}_r ----- START 137 so_term=0 ! ZORA-spin-free 138 do i=1,ipol 139 call ga_zero(g_efgz4_sf(i)) 140 enddo 141 do i=1,3 142 call ga_zero(g_efgz4_so(i)) 143 enddo 144 call zora_getv_EFGZ4_SO(rtdb,g_dens_ini, 145 & zora_calc_type, 146 & zora_Qpq,xyz_EFGcoords, 147 & g_efgz4_sf, ! out: munu matrix 148 & g_efgz4_so, ! out: munu matrix 149 & nexc) 150 do i=1,ipol 151 toac(i)=ga_ddot(g_densZ4so(i),g_efgz4_sf(i)) 152 enddo ! end-loop-ipol 153 if (zora_calc_type.eq.3) then 154 ac=3 155 do i=1,3 156 toac(ac)=ga_ddot(g_densZ4so(ac),g_efgz4_so(i)) 157 ac=ac+1 158 enddo 159 zora_eint=toac(1)+toac(2)+ 160 & toac(3)+toac(4)+toac(5) 161 if (ga_nodeid().eq.0) then 162 write(*,155) zora_calc_type,iat,zora_Qpq,pos_dat, 163 & toac(1),toac(2),toac(1)+toac(2), 164 & toac(3),toac(4),toac(5),zora_eint 165 155 format('zora-efg-so(',i3,',',i3,',',i3,',',i3,')=(', 166 & f15.8,',',f15.8,',',f15.8,',',f15.8,',', 167 & f15.8,',',f15.8,',',f15.8,')') 168 endif 169 else 170 zora_eint=toac(1)+toac(2) 171 if (ga_nodeid().eq.0) then 172 write(*,15) zora_calc_type,iat,zora_Qpq,pos_dat, 173 & toac(1),toac(2),toac(1)+toac(2), 174 & zora_eint 175 15 format('zora-efg-so(',i3,',',i3,',',i3,',',i3,')=(', 176 & f15.8,',',f15.8,',', f15.8,',',f15.8,')') 177 endif 178 endif 179 call ga_fill_patch(g_rhoS,1,1,pos_dat,pos_dat,zora_eint) 180c------Generate munu A^{pq}_r ----- END 181 pos_dat=pos_dat+1 182 end do ! zora_Qpq loop 183 end do ! iat-loop 184 do i=1,5 185 if (.not. ga_destroy(g_densZ4so(i))) call errquit( 186 & 'dft_zora_rhos_so: ga_destroy failed ',0, GA_ERR) 187 enddo 188 end do ! zora_calc_type loop 189c ----- Store efgz4 data in a file ------- START 190 call ga_sync() 191c Note.- lbl_efgz4 defined in zora.fh 192 call util_file_name(lbl_efgz4,.false.,.false.,zorafilename) 193 if (.not.dft_zoraEFGZ4_write( 194 & zorafilename, 195 & nlist, 196 & nat, 197 & g_AtNr, 198 & g_rhoS)) 199 & call errquit('get_rhoS_so: dft_zoraEFGZ4_write failed', 200 & 0,DISK_ERR) 201c ----- Store efgz4 data in a file ------- END 202c----deallocate memory - FA 203 if (.not.ma_free_heap(l_zanpt)) call errquit 204 & ('get_rhoS_so, ma_free_heap of l_zanpt failed',911,MA_ERR) 205 if (.not.ma_free_heap(l_xyzpt)) call errquit 206 & ('get_rhoS_so, ma_free_heap of l_xyzpt failed',911,MA_ERR) 207 if (.not.ma_free_heap(l_AtNr)) call 208 & errquit('dft_zora_utils: ma_free_heap l_AtNr',0, MA_ERR) 209 do i=1,ipol 210 if (.not. ga_destroy(g_efgz4_sf(i))) call errquit( 211 & 'dft_zora_rhos_so: ga_destroy failed ',0, GA_ERR) 212 enddo 213 do i=1,3 214 if (.not. ga_destroy(g_efgz4_so(i))) call errquit( 215 & 'dft_zora_rhos_so: ga_destroy failed ',0, GA_ERR) 216 enddo 217 if (.not. ga_destroy(g_rhoS)) call errquit( 218 & 'dft_zora_rhos_so: ga_destroy failed ',0, GA_ERR) 219 return 220 end 221 222 subroutine get_densZ4_so( 223 & rtdb, 224 & geom, 225 & ao_bas_han, 226 & nbf,nbf_ao,nbf_mo, 227 & g_dens,g_moso, 228 & noc,nocc, 229 & g_densZ4so) 230c Purpose : Calculation of density matrix-like 231c for EFG-ZORA-4 232c Output : g_densZ4so 233c Author : Fredy Aquino 234c Note : Modified from dft_zora_scale_so() 235c Date : 12-09-09 236 implicit none 237#include "nwc_const.fh" 238#include "errquit.fh" 239#include "mafdecls.fh" 240#include "stdio.fh" 241#include "global.fh" 242#include "msgids.fh" 243#include "rtdb.fh" 244#include "zora.fh" 245 246 integer ga_create_atom_blocked 247 external ga_create_atom_blocked 248 249 integer g_dens(2) 250 integer g_moso(2) 251 integer g_orb(2) 252 integer g_dens_so(2) 253 integer g_scr(2) 254 integer l_vecsre, k_vecsre 255 integer l_vecsim, k_vecsim 256 integer iorb 257 integer noc 258 integer ispin 259 integer geom 260 integer ao_bas_han 261 integer nbf 262 integer nbf_ao 263 integer nbf_mo 264 integer g_densZ4so(5) ! FA - Main output 265 integer rtdb ! FA 266 integer i,j,ac,noc1,count ! FA 267 integer nocc(2) ! FA 268 integer l_Ci,k_Ci ! FA 269 character*1 soxyz(3) ! FA 270 character*7 lbls(6) ! FA 271 character*34 lbls1(6) ! FA 272 double precision scale ! FA 273c Using global input : g_Ci 274c defined in zora.fh and calculated 275c in dft_zora_scale_so() 276 soxyz(1)='z' 277 soxyz(2)='y' 278 soxyz(3)='x' 279 lbls(1)='orbs re' 280 lbls(2)='denmxre' 281 lbls(3)='scrch 1' 282 lbls(4)='orbs im' 283 lbls(5)='denmxim' 284 lbls(6)='scrch 2' 285 lbls1(1)='get_densZ4_so: orb real error' 286 lbls1(2)='get_densZ4_so: dens real error' 287 lbls1(3)='get_densZ4_so: ga_duplicate failed' 288 lbls1(4)='get_densZ4_so: orb imag error' 289 lbls1(5)='get_densZ4_so: dens imag error' 290 lbls1(6)='get_densZ4_so: ga_duplicate failed' 291 292 noc1=nocc(1)+nocc(2) 293 294cc allocate memory 295 if (.not.MA_Push_Get(MT_Dbl, nbf_mo, 'real vec aux', 296 & l_vecsre, k_vecsre)) 297 & call errquit('get_densZ4_so: cannot allocate vec',0, MA_ERR) 298c 299 if (.not.MA_Push_Get(MT_Dbl, nbf_mo, 'imag vec aux', 300 & l_vecsim, k_vecsim)) 301 & call errquit('get_densZ4_so: cannot allocate vec',0, MA_ERR) 302 303 if (.not.ma_alloc_get(MT_Dbl, 2*noc1, 'Ci', 304 & l_Ci, k_Ci)) 305 & call errquit('get_densZ4_so: ma failed',0,MA_ERR) 306 307 call ga_get(g_Ci,1,2,1,noc1,dbl_mb(k_Ci),2) 308 309 ac=1 310 do i=1,2 311c spin-orbit vector - i=1(real),2(imaginary) 312 if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,lbls(ac),0,0, 313 & g_orb(i))) 314 & call errquit(lbls1(ac),0, GA_ERR) 315 call ga_zero(g_orb(i)) 316 ac=ac+1 317c spin-orbit density matrix - i=1(real),2(imaginary) 318 if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,lbls(ac),0,0, 319 & g_dens_so(i))) 320 & call errquit(lbls1(ac),0, GA_ERR) 321 call ga_zero(g_dens_so(i)) 322 ac=ac+1 323c scratch array i=1(real),2(imaginary) 324 if(.not.ga_duplicate(g_dens(i),g_scr(i),lbls(ac))) 325 & call errquit(lbls1(ac),1, GA_ERR) 326 call ga_zero(g_scr(i)) 327 ac=ac+1 328 end do ! end i-loop 329 330 do i=1,5 ! Added by FA 331 if(.not.ga_duplicate(g_dens(1),g_densZ4so(i),'densZ4so')) 332 & call errquit('get_densZ4_so: ga_duplicate failed',1, GA_ERR) 333 call ga_zero(g_densZ4so(i)) 334 end do 335 336 do ispin=1,2 337 iorb=ispin 338 do count=1,nocc(ispin) 339 call ga_get(g_moso(1),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsre),1) 340 call ga_zero(g_orb(1)) 341 call ga_put(g_orb(1),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsre),1) 342 call ga_get(g_moso(2),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsim),1) 343 call ga_zero(g_orb(2)) 344 call ga_put(g_orb(2),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsim),1) 345 call dft_densm_so(g_dens_so,g_orb,nbf_ao,noc) 346 call ga_zero(g_scr(1)) 347 call ga_zero(g_scr(2)) 348 call ga_dens_sf(g_scr, g_dens_so(1), nbf_ao) 349 if (zora_calc_type.eq.3) scale=1.0d0/dbl_mb(k_Ci+iorb-1) 350 if (zora_calc_type.eq.4 .or. zora_calc_type.eq.5) scale=1.0d0 ! NUM-EFG,NUM-NMR-K=1 351 do i=1,2 352 call ga_scale(g_scr(i),scale) 353 call ga_add(1.0d0,g_densZ4so(i),1.0d0,g_scr(i),g_densZ4so(i)) 354 end do ! end-loop-i 355 j=3 356 do i=1,3 357 call ga_zero(g_scr(1)) 358 call ga_dens_so(g_scr(1),g_dens_so,nbf_ao,soxyz(i)) 359 call ga_scale(g_scr(1),scale) 360 call ga_add(1.0d0,g_densZ4so(j),1.0d0,g_scr(1),g_densZ4so(j)) 361 j=j+1 362 end do ! end-loop-i 363 iorb=iorb+2 364 end do ! count-loop 365 end do ! ispin-loop 366 367c deallocate memory 368 if (.not. ma_chop_stack(l_vecsim)) 369 & call errquit('get_densZ4_so:l_vecsim', 0, MA_ERR) 370 if (.not. ma_chop_stack(l_vecsre)) 371 & call errquit('get_densZ4_so:l_vecsre', 0, MA_ERR) 372 if (.not. MA_free_heap(l_Ci)) 373 & call errquit('get_densZ4_so:cannot free heap',111, MA_ERR) 374 do i=1,2 375 if (.not. ga_destroy(g_orb(i))) 376 & call errquit('get_densZ4_so: ga_destroy failed',0, GA_ERR) 377 if (.not. ga_destroy(g_dens_so(i))) 378 & call errquit('get_densZ4_so: ga_destroy failed',0, GA_ERR) 379 if (.not. ga_destroy(g_scr(i))) 380 & call errquit('dft_zora_scale_so: ga_destroy failed',0, GA_ERR) 381 end do 382 383 return 384 end 385 386 subroutine hnd_efgmap_Z4_so(rtdb,basis,geom, 387 & nbf,nbf_ao,nbf_mo, 388 & g_dens_4SO,g_moso, 389 & noc,nocc) 390c 391c $Id$ 392c 393c This routine calculates the electric field gradient and 394c the orientation of the EFG for a given density at the 395c atomic positions. 396 397 implicit none 398#include "nwc_const.fh" 399#include "errquit.fh" 400#include "global.fh" 401#include "bas.fh" 402#include "mafdecls.fh" 403#include "geom.fh" 404#include "stdio.fh" 405#include "rtdb.fh" 406#include "cosmo.fh" 407#include "zora.fh" ! Added by FA 408 409 integer rtdb ! [Input] rtdb 410 integer basis ! [Input] Basis set 411 integer geom ! [Input] Geometry 412 character*255 zorafilename 413 character*2 symbol 414 character*16 element, at_tag 415 integer iat, atn, nat, i, j, ij 416 integer l_xyzpt , k_xyzpt, 417 & l_xyzpt1, k_xyzpt1, 418 & l_zanpt , k_zanpt, 419 & l_efgs , k_efgs 420 integer l_tmp,k_tmp ! to store indices of selected atoms 421 integer g_dens(3),ndens,nclosed(2),nopen(2),nvirt(2) 422 integer nefc, l_efcc, k_efcc, l_efcz, k_efcz 423 double precision xp, yp, zp, xn, yn, zn, zan 424 double precision vec(3,3), eig(3), a(6) 425 double precision pi, deg, efgxx, efgyy, efgzz, efgxy, efgxz, efgyz 426 double precision rr, rr5, eta 427 logical status 428c bq variables (MV) 429 logical dobq 430 integer bq_ncent 431 integer i_cbq 432 integer i_qbq 433 double precision elpotbq 434 logical dft_zoraEFGZ4_read 435 integer g_rhoS,g_Atnr1,icalczora 436 integer nder 437 integer l_rhoS,k_rhoS 438 integer ii,jj 439 integer g_densZ4(3) 440 double precision sum_efgs 441 integer ipol,pos,indx,indy 442 integer count_efgtyp 443 integer g_densZ4so(5) 444 integer g_dens_4SO(2) 445 integer g_moso(2) 446 integer nbf,nbf_ao,nbf_mo, 447 & noc,nocc(*),efgfile,typeprop, 448 & indx1,indx2,nat_slc 449 external hnd_elfcon_symm,hnd_elfcon 450 external get_densZ4_so, 451 & dft_zoraEFGZ4_read,util_file_name 452 integer iat1 453c 454c Initialize integrals 455 456 call int_init(rtdb,1, basis) 457 call schwarz_init(geom, basis) 458c 459c ----- calculate electric field gradient ----- 460 461 if (ga_nodeid().eq.0) write(luout,9999) 462 if (ga_nodeid().eq.0) write(luout,9994) 463 464 pi = acos(-1.0d0) 465 deg = 180.0d0/pi 466c 467 call ecce_print_module_entry('efg') 468c 469c ----- define points for calculation ----- 470c 1. nuclei 471c ------- Read (nat,atmnr) --------- START 472 status=geom_ncent(geom,nat) 473 if (.not.ma_alloc_get( 474 & mt_int,nat,'nmt tmp',l_tmp,k_tmp)) 475 & call errquit('hnd_efgmap_Z4: ma failed',0,MA_ERR) 476 typeprop=1 ! =1 EFG =2 Shieldings =3 Hyperfine 477 call get_slctd_atoms(nat_slc, ! out: selected atoms 478 & int_mb(k_tmp), ! out: list of selected atom nr. 479 & nat, ! in : total nr atoms in molecule 480 & rtdb, ! in : rdt handle 481 & typeprop) ! in : =1,2,3=EFG,Shieldings,Hyperfine 482 if (ga_nodeid().eq.0) then 483 write(*,*) 'nat_slc=',nat_slc 484 do i=1,nat_slc 485 write(*,7) i,int_mb(k_tmp+i-1) 486 7 format('In hnd_efgmap_Z4:: atomnr(',i3,')=',i5) 487 enddo 488 endif 489c ------- Read (nat,atmnr) --------- END 490 if (.not.ma_alloc_get( 491 & mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 492 & call errquit('hnd_efgmap_Z4: ma failed',0,MA_ERR) 493 if (.not.ma_alloc_get( 494 & mt_dbl,6*nat,'efg pnt',l_efgs,k_efgs)) 495 & call errquit('hnd_efgmap_Z4: ma failed',0,MA_ERR) 496 if (.not.ma_alloc_get( 497 & mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 498 & call errquit('hnd_efgmap_Z4: ma failed',0,MA_ERR) 499 500 do 30 iat=1,nat 501 status=geom_cent_get(geom,iat,at_tag, 502 & dbl_mb(k_xyzpt+3*(iat-1)), 503 & dbl_mb(k_zanpt+iat-1)) 504 30 continue 505c ++++++ Reading EFGZ4 data from file ++++++ START 506c Note.- lbl_efgz4 defined in zora.fh 507 call util_file_name(lbl_efgz4,.false.,.false.,zorafilename) 508 icalczora = 0 ! initialize the flag 509 if (.not.dft_zoraEFGZ4_read( 510 & zorafilename, 511 & nat_slc, 512 & nat, 513 & g_AtNr1, 514 & g_rhoS)) icalczora=1 515c Note.- If I print the GAs here it gets freezed 516c ++++++ Reading EFGZ4 data from file ++++++ END 517 if (ga_nodeid().eq.0) 518 & write(*,*) '-------hnd_efgmat_Z4: g_rhoS ---------- START' 519 call ga_print(g_AtNr1) 520 call ga_print(g_rhoS) 521 if (ga_nodeid().eq.0) 522 & write(*,*) '-------hnd_efgmat_Z4: g_rhoS ---------- END' 523c 524c Allocate memory for l_rhoS,k_rhoS 525 if (.not.ma_alloc_get( 526 & mt_dbl,nat_slc*6*2,'rhoS',l_rhoS,k_rhoS)) 527 & call errquit('hnd_efgmap_Z4_so: ma failed',0,MA_ERR) 528 call ga_get(g_rhoS,1,1,1,nat_slc*6*2,dbl_mb(k_rhoS),1) 529 530c count_efgtyp=1 531 count_efgtyp=0 ! NUM-EFG, Z4-EFG 532 do zora_calc_type=4,3,-1 533 do_NonRel=.false. 534 if (zora_calc_type.eq.4) do_NonRel=.true. 535 call get_densZ4_so(rtdb,geom,basis, 536 & nbf,nbf_ao,nbf_mo, 537 & g_dens_4SO,g_moso, 538 & noc,nocc,g_densZ4so) 539 540c if (zora_calc_type.eq.4) then 541c if (ga_nodeid().eq.0) 542c & write(*,*) '----g_densZ4so(1)----- START' 543c call ga_print(g_densZ4so(1)) 544c if (ga_nodeid().eq.0) 545c & write(*,*) '----g_densZ4so(1)----- END' 546c if (ga_nodeid().eq.0) 547c & write(*,*) '----g_densZ4so(2)----- START' 548c call ga_print(g_densZ4so(2)) 549c if (ga_nodeid().eq.0) 550c & write(*,*) '----g_densZ4so(2)----- END' 551c endif 552 553c ----- compute g_densZ4_all=g_densZ4so(1)+g_densZ4so(2) 554 call ga_add(1.0d00,g_densZ4so(1), 555 & 1.0d00,g_densZ4so(2), 556 & g_densZ4so(1)) 557 nder=2 558 if (.not.ma_alloc_get( 559 & mt_dbl,3*nat_slc,'xyz pnt1',l_xyzpt1,k_xyzpt1)) 560 & call errquit('hnd_efgmap_Z4: ma failed',0,MA_ERR) 561 do iat1=1,nat_slc 562 iat=int_mb(k_tmp+iat1-1) 563 indx1=k_xyzpt1+3*(iat1-1) 564 indx2=k_xyzpt +3*(iat -1) 565 dbl_mb(indx1 )= dbl_mb(indx2 ) 566 dbl_mb(indx1+1)= dbl_mb(indx2+1) 567 dbl_mb(indx1+2)= dbl_mb(indx2+2) 568 if (ga_nodeid().eq.0) then 569 write(*,12) iat1,iat, 570 & dbl_mb(indx1), 571 & dbl_mb(indx1+1), 572 & dbl_mb(indx1+2) 573 12 format('Atom(',i3,',',i3,')=(', 574 & f15.8,',',f15.8,',',f15.8,')') 575 endif 576 enddo 577c ---- extract selected atoms coordinates ----- END 578c goto 13 579 efgfile=0 ! NO NLMO/NBO analysis 580 call hnd_elfcon_symm(basis, ! in: basis handle 581 & geom, ! in: geom handle 582 & g_densZ4so(1), ! in: electron density 583 & dbl_mb(k_xyzpt1), ! in: (x,y,z) centers 584 & nat_slc, ! in: number of centers 585 & dbl_mb(k_efgs), !out: EFG values at centers 586 & nder, ! in: =2 for second derivative 587 & efgfile) ! in: efgfile=0,1= NO,YES NLMONBO analysis 588c 13 continue 589 590c call hnd_elfcon(basis, ! in: basis handle 591c & geom, ! in: geom handle 592c & g_densZ4so(1), ! in: electron density 593c & dbl_mb(k_xyzpt1), ! in: (x,y,z) centers 594c & nat_slc, ! in: number of centers 595c & dbl_mb(k_efgs), !out: EFG values at centers 596c & nder) ! in: =2 for second derivative 597 598c Note.- Getting NaN in dbl_mb(k_efgs) for CuI + HiraoBS_uncontracted + BLYP 599c filename: CuI_HiraoBS_uCNT_Z4_BLYP.pbs 600c if (ga_nodeid().eq.0) then 601c write(*,*) '--------check k_efgs ------ START' 602c do i=1,6*nat 603c write(*,171) i,dbl_mb(k_efgs+i-1) 604c 171 format('efgs(',i3,')=',f15.8) 605c enddo 606c write(*,*) '--------check k_efgs ------ END' 607c endif 608 609 if (.not.ma_free_heap(l_xyzpt1)) call 610 & errquit('hnd_efgmap_Z4_so: ma_free_heap l_xyzpt1',0, MA_ERR) 611 612 if (ga_nodeid().eq.0) then ! START-if-do-it-once 613 write(*,112) 614 do iat1=1,nat_slc 615 iat=int_mb(k_tmp+iat1-1) 616c ------- get Atom name: symbol ----------- START 617 if (.not. geom_cent_tag(geom,iat,at_tag)) call 618 & errquit('hnd_efgmap_Z4_so: geom_cent_tag failed',0,GEOM_ERR) 619 if (.not. geom_tag_to_element(at_tag,symbol,element,atn)) then 620 if(symbol.ne."bq") call 621 & errquit('hnd_efgmap_Z4_so: geom_tag_to_element failed', 622 & 0,GEOM_ERR) 623 end if 624c ------- get Atom name: symbol ----------- END 625 indx=k_efgs+6*(iat1-1) 626 efgxx = dbl_mb(indx) 627 efgyy = dbl_mb(indx+1) 628 efgzz = dbl_mb(indx+2) 629 efgxy = dbl_mb(indx+3) 630 efgxz = dbl_mb(indx+4) 631 efgyz = dbl_mb(indx+5) 632 sum_efgs=(efgxx+efgyy+efgzz)/3.0d0 633 efgxx=efgxx-sum_efgs 634 efgyy=efgyy-sum_efgs 635 efgzz=efgzz-sum_efgs 636 indx=k_efgs+6*(iat1-1) 637 indy=k_rhoS+6*(iat1-1)+6*nlist*count_efgtyp 638 if (zora_calc_type.eq.3) then 639 dbl_mb(indx )=efgxx+dbl_mb(indy ) 640 dbl_mb(indx+1)=efgyy+dbl_mb(indy+1) 641 dbl_mb(indx+2)=efgzz+dbl_mb(indy+2) 642 dbl_mb(indx+3)=efgxy+dbl_mb(indy+3) 643 dbl_mb(indx+4)=efgxz+dbl_mb(indy+4) 644 dbl_mb(indx+5)=efgyz+dbl_mb(indy+5) 645 write(*,113) symbol,efgxx,efgyy,efgzz,efgxy, 646 & efgxz,efgyz 647 write(*,114) symbol,dbl_mb(indy),dbl_mb(indy+1), 648 & dbl_mb(indy+2),dbl_mb(indy+3), 649 & dbl_mb(indy+4),dbl_mb(indy+5) 650 write(*,115) symbol,dbl_mb(indx),dbl_mb(indx+1), 651 & dbl_mb(indx+2),dbl_mb(indx+3), 652 & dbl_mb(indx+4),dbl_mb(indx+5) 653 endif 654 if (zora_calc_type.eq.4) then 655 dbl_mb(indx )=efgxx 656 dbl_mb(indx+1)=efgyy 657 dbl_mb(indx+2)=efgzz 658 dbl_mb(indx+3)=efgxy 659 dbl_mb(indx+4)=efgxz 660 dbl_mb(indx+5)=efgyz 661 write(*,116) symbol,efgxx,efgyy,efgzz,efgxy, 662 & efgxz,efgyz 663 write(*,117) symbol,dbl_mb(indy),dbl_mb(indy+1), 664 & dbl_mb(indy+2),dbl_mb(indy+3), 665 & dbl_mb(indy+4),dbl_mb(indy+5) 666 endif 667 end do ! iat-loop 668 end if ! END-if-do-it-once 669c count_efgtyp=count_efgtyp-1 670 count_efgtyp=count_efgtyp+1 ! NUM-EFG, Z4-EFG 671 do i=1,5 672 if (.not. ga_destroy(g_densZ4so(i))) call errquit( 673 & 'dft_zora_rhos_so: ga_destroy failed ',0, GA_ERR) 674 enddo 675 end do ! zora_calc_type loop 676c ------- All-FA-formats------------------------------------ START 677 112 format('====> Electronic contribution to EFG', 678 & ' in molecular frame (a.u.)',/, 679 & 21x,'XX',12x,'YY',12x,'ZZ',12x,'XY',12x,'XZ',12x,'YZ',/, 680 & 16x,82(1h-)) 681 113 format('EFG-elec(',a2,')=(',f13.8,',',f13.8,',', 682 & f13.8,',',f13.8,',',f13.8,',',f13.8,')') 683 114 format('EFG-rhoS(',a2,')=(',f13.8,',',f13.8,',', 684 & f13.8,',',f13.8,',',f13.8,',',f13.8,')') 685 115 format('EFG-tot (',a2,')=(',f13.8,',',f13.8,',', 686 & f13.8,',',f13.8,',',f13.8,',',f13.8,')') 687 116 format(' ANALYT(',a2,')=(',f13.8,',',f13.8,',', 688 & f13.8,',',f13.8,',',f13.8,',',f13.8,')') 689 117 format(' NUMERI(',a2,')=(',f13.8,',',f13.8,',', 690 & f13.8,',',f13.8,',',f13.8,',',f13.8,')') 691c ------- All-FA-formats------------------------------------ END 692c 693c get bq structures if any (MV) 694c ----------------------------- 695 dobq = .false. 696 if(geom_extbq_on()) then 697 dobq = .true. 698 bq_ncent = geom_extbq_ncenter() 699 i_cbq = geom_extbq_coord() 700 i_qbq = geom_extbq_charge() 701 end if 702c 703c ----- collect and output results of all points ----- 704 705 status = rtdb_parallel(.false.) ! FA-04-23-10 706 if (ga_nodeid().gt.0) goto 300 707 do 230 iat=1,nat_slc 708 iat1=int_mb(k_tmp+iat-1) 709 xp = dbl_mb(k_xyzpt +3*(iat1-1)) 710 yp = dbl_mb(k_xyzpt+1+3*(iat1-1)) 711 zp = dbl_mb(k_xyzpt+2+3*(iat1-1)) 712c ----- add nuclear contribution ----- 713 efgxx = 0.0d0 ! FA 714 efgyy = 0.0d0 ! FA 715 efgzz = 0.0d0 ! FA 716 efgxy = 0.0d0 ! FA 717 efgxz = 0.0d0 ! FA 718 efgyz = 0.0d0 ! FA 719 do 210 i = 1,nat 720 xn = dbl_mb(k_xyzpt +3*(i-1)) - xp 721 yn = dbl_mb(k_xyzpt+1+3*(i-1)) - yp 722 zn = dbl_mb(k_xyzpt+2+3*(i-1)) - zp 723 zan = dbl_mb(k_zanpt+i-1) 724 rr = sqrt(xn*xn + yn*yn + zn*zn) 725 if (rr.lt.1.0d-3) go to 210 726 rr5=rr*rr*rr*rr*rr 727 efgxx = efgxx - zan*xn*xn/rr5 728 efgyy = efgyy - zan*yn*yn/rr5 729 efgzz = efgzz - zan*zn*zn/rr5 730 efgxy = efgxy - zan*xn*yn/rr5 731 efgxz = efgxz - zan*xn*zn/rr5 732 efgyz = efgyz - zan*yn*zn/rr5 733 210 continue 734c 735c ----- form -efc- contribution ----- 736c from cosmo point charges !!!! 737 if (cosmo_last) then 738 if (.not.rtdb_get(rtdb,'cosmo:nefc',mt_int,1,nefc)) 739 & call errquit('hnd_efgmap: rtdb get failed for nefc ',911, 740 & RTDB_ERR) 741 if (.not.ma_push_get(mt_dbl,nefc*3,'efcc',l_efcc,k_efcc)) 742 & call errquit('hnd_efgmap: malloc k_efcc fail',911,ma_err) 743 if (.not.ma_push_get(mt_dbl,nefc,'efcz',l_efcz,k_efcz)) 744 & call errquit('hnd_efgmap: malloc k_efcz fail',911,ma_err) 745 if (.not.rtdb_get(rtdb,'cosmo:efcc',mt_dbl,3*nefc, 746 & dbl_mb(k_efcc))) call 747 & errquit('hnd_efgmap: rtdb get failed efcc',912,rtdb_err) 748 if (.not.rtdb_get(rtdb,'cosmo:efcz',mt_dbl,nefc, 749 & dbl_mb(k_efcz))) call 750 & errquit('hnd_efgmap: rtdb get failed efcz',913,rtdb_err) 751 do i = 1,nefc 752 xn = dbl_mb(k_efcc+3*(i-1) ) - xp 753 yn = dbl_mb(k_efcc+3*(i-1)+1) - yp 754 zn = dbl_mb(k_efcc+3*(i-1)+2) - zp 755 rr = sqrt(xn*xn + yn*yn + zn*zn) 756 if (rr.lt.1.0d-3) then 757 write(luout,9993) xp,yp,zp,i 758 else 759 rr5=rr*rr*rr*rr*rr 760 efgxx = efgxx - dbl_mb(k_efcz+i-1)*xn*xn/rr5 761 efgyy = efgyy - dbl_mb(k_efcz+i-1)*yn*yn/rr5 762 efgzz = efgzz - dbl_mb(k_efcz+i-1)*zn*zn/rr5 763 efgxy = efgxy - dbl_mb(k_efcz+i-1)*xn*yn/rr5 764 efgxz = efgxz - dbl_mb(k_efcz+i-1)*xn*zn/rr5 765 efgyz = efgyz - dbl_mb(k_efcz+i-1)*yn*zn/rr5 766 endif 767 enddo 768 220 continue 769 if (.not.ma_chop_stack(l_efcc)) call 770 & errquit('hnd_efgmap: chop stack l_efcc',913,ma_err) 771 endif 772c 773c adding external bq contributions(MV) 774c ---------------------------------- 775 if (dobq) then 776 do i = 1,bq_ncent 777 xn = dbl_mb(i_cbq+3*(i-1) ) - xp 778 yn = dbl_mb(i_cbq+3*(i-1)+1) - yp 779 zn = dbl_mb(i_cbq+3*(i-1)+2) - zp 780 rr = sqrt(xn*xn + yn*yn + zn*zn) 781 if (rr.lt.1.0d-3) then 782 write(luout,9993) xp,yp,zp,i 783 else 784 rr5=rr*rr*rr*rr*rr 785 efgxx = efgxx - dbl_mb(i_qbq+i-1)*xn*xn/rr5 786 efgyy = efgyy - dbl_mb(i_qbq+i-1)*yn*yn/rr5 787 efgzz = efgzz - dbl_mb(i_qbq+i-1)*zn*zn/rr5 788 efgxy = efgxy - dbl_mb(i_qbq+i-1)*xn*yn/rr5 789 efgxz = efgxz - dbl_mb(i_qbq+i-1)*xn*zn/rr5 790 efgyz = efgyz - dbl_mb(i_qbq+i-1)*yn*zn/rr5 791 endif 792 end do 793 end if 794c ------- Adding modified electronic part + nuclear contribution 795 indx=k_efgs+6*(iat-1) 796 dbl_mb(indx )=dbl_mb(indx )+2.0d0*efgxx - efgyy - efgzz 797 dbl_mb(indx+1)=dbl_mb(indx+1)+2.0d0*efgyy - efgxx - efgzz 798 dbl_mb(indx+2)=dbl_mb(indx+2)+2.0d0*efgzz - efgxx - efgyy 799 dbl_mb(indx+3)=dbl_mb(indx+3)+3.0d0*efgxy 800 dbl_mb(indx+4)=dbl_mb(indx+4)+3.0d0*efgxz 801 dbl_mb(indx+5)=dbl_mb(indx+5)+3.0d0*efgyz 802c 803c ----- reorder into a as xx xy yy xz yz zz to form matrix ----- 804 a(1) = dbl_mb(k_efgs +6*(iat-1)) 805 a(2) = dbl_mb(k_efgs+3+6*(iat-1)) 806 a(3) = dbl_mb(k_efgs+1+6*(iat-1)) 807 a(4) = dbl_mb(k_efgs+4+6*(iat-1)) 808 a(5) = dbl_mb(k_efgs+5+6*(iat-1)) 809 a(6) = dbl_mb(k_efgs+2+6*(iat-1)) 810 ij=0 811 do 241 i = 1, 3 812 do 241 j = 1, i 813 ij = ij + 1 814 vec(i,j) = a(ij) 815 vec(j,i) = a(ij) 816 241 continue 817c 818c ----- store ecce data ----- 819 if (.not. geom_cent_tag(geom,iat1,at_tag)) call 820 & errquit('hnd_efgmap: geom_cent_tag failed',0,GEOM_ERR) 821c geom_tag_to_element returns false for Bq elements (MV) 822c ----------------------------------------------------- 823 if (.not. geom_tag_to_element(at_tag,symbol,element,atn)) then 824 if(symbol.ne."bq") call 825 & errquit('hnd_efgmap: geom_tag_to_element failed',0,GEOM_ERR) 826 end if 827c 828c if (.not. geom_tag_to_element(at_tag,symbol,element,atn)) call 829c & errquit('hnd_efgmap: geom_tag_to_element failed',0,GEOM_ERR) 830 831 call ecce_print1_char('atom name',symbol,1) 832 call ecce_print2('EFG tensor',MT_DBL,vec,3,3,3) 833c 834c ----- print tensor components ----- 835 836 write(luout,9998) iat,symbol,xp,yp,zp 837 write(luout,9997) 838 write(luout,9995) (dbl_mb(k_efgs+6*(iat-1)+i),i=0,5) 839c 840c ----- diagonalize to get principal components and vectors ----- 841c 842c FA: I found that for some few cases it halts when calling hnd_diag() 843c I don't know why. 844c Example: AuI_SARC1-ZORA_BS_CNT_Z4_BLYP_4.out (04-16-10) 845 846 call hnd_diag(vec,eig,3,.true.,.false.) 847 eta = abs( (eig(3)-eig(2)) / eig(1) ) 848 849 call ecce_print1('EFG eigenvalues',MT_DBL,eig,3) 850 call ecce_print2('EFG eigenvectors',MT_DBL,vec,3,3,3) 851 call ecce_print1('EFG asymmetry',MT_DBL,eta,1) 852 853 write(luout,9992) 854 write(luout,9991) eig(1),eig(2),eig(3),eta 855 write(luout,9988) ((vec(i,j),j=1,3),i=1,3) 856 write(luout,*) ' ' 857c 858 230 continue ! Assemblin and printing next atom 859 call ecce_print_module_exit('EFG','ok') 860 call util_flush(luout) 861 862c ----- release memory block ----- 863 300 call ga_sync() 864 status = rtdb_parallel(.true.) ! FA-04-23-10 865c 866c ------- Deallocate MA memory ------ 867 if (.not.ma_free_heap(l_rhoS)) call 868 & errquit('hnd_efgmap_Z4: ma_free_heap l_rhoS',0, MA_ERR) 869 if (.not.ma_free_heap(l_zanpt)) call 870 & errquit('hnd_efgmap_Z4: ma_free_heap l_zanpt',0, MA_ERR) 871 if (.not.ma_free_heap(l_efgs)) call 872 & errquit('hnd_efgmap_Z4: ma_free_heap l_efgs',0, MA_ERR) 873 if (.not.ma_free_heap(l_xyzpt)) call 874 & errquit('hnd_efgmap_Z4: ma_free_heap l_xyzpt',0, MA_ERR) 875 if (.not.ma_free_heap(l_tmp)) call 876 & errquit('hnd_efgmap_Z4: ma_free_heap l_tmp',0, MA_ERR) 877c 878 return 879 9999 format(/,10x,23(1h-),/,10x,'Electric field gradient', 880 1 /,10x,23(1h-),/) 881 9998 format(/,1x,60(1h-),/,3x,'Atom',6x,'X',9x,'Y',9x,'Z',/,1x,60(1h-), 882 1 /,i5,1x,a2,3f10.5,/,1x,60(1h-),/) 883 9997 format(1x,'Electric field gradient in molecular frame (a.u.)',/, 884 2 9x,'XX',13x,'YY',13x,'ZZ',13x,'XY',13x,'XZ',13x,'YZ',/, 885 3 1x,90(1h-)) 886 9996 format(' --- Warning - electric field gradient at ', 887 1 3F10.5,' . contribution from nucleus ',i3,' ignored') 888 9995 format(1x,6f15.6,/) 889 9994 format(' 1 a.u. = 0.324123 10**(16) esu/cm**3 ', 890 1 ' ( or statvolts/cm**2 )',' = 0.97174 10**(22) v/m**2 ',/) 891 9993 format(' --- Warning - electric field gradient at ', 892 1 3f10.5,' . contribution from -efc- ',i3,' ignored') 893 9992 format(1x,'Principal components (a.u.) and orientation ', 894 1 /,' of principal axis w.r.t. absolute frame', 895 2 22x,'Asymmetry parameter eta',/,1x,86(1h-)) 896 9991 format(1x,3f15.6,16x,f15.6,/) 897 9988 format(1X,3F15.6) 898 end 899c======================================================= 900c================ Fredy Aquino's routines ======== END 901c======================================================= 902