1 subroutine hnd_efgmap(rtdb,basis,geom) 2c 3c $Id$ 4c 5c This routine calculates the electric field gradient and 6c the orientation of the EFG for a given density at the 7c atomic positions. 8c 9 implicit none 10#include "nwc_const.fh" 11#include "errquit.fh" 12#include "global.fh" 13#include "bas.fh" 14#include "mafdecls.fh" 15#include "geom.fh" 16#include "stdio.fh" 17#include "rtdb.fh" 18#include "cosmo.fh" 19c 20 integer rtdb ! [Input] rtdb 21 integer basis ! [Input] Basis set 22 integer geom ! [Input] Geometry 23c 24 character*2 symbol 25 character*16 element, at_tag 26 integer iat, atn, nat, i, j, ij 27 integer l_xyzpt, k_xyzpt, l_zanpt, k_zanpt, l_efgs, k_efgs 28 integer g_dens(3),ndens,nclosed(2),nopen(2),nvirt(2) 29 integer nefc, l_efcc, k_efcc, l_efcz, k_efcz 30 character*3 scftyp 31 double precision xp, yp, zp, xn, yn, zn, zan 32 double precision vec(3,3), eig(3), a(6) 33 double precision pi, deg, efgxx, efgyy, efgzz, efgxy, efgxz, efgyz 34 double precision rr, rr5, eta 35 logical status 36c bq variables (MV) 37 logical dobq 38 integer bq_ncent 39 integer i_cbq 40 integer i_qbq 41 double precision elpotbq 42c 43 logical docosmo 44 integer ncosbq 45c 46c Initialize integrals 47c 48 call int_init(rtdb,1, basis) 49 call schwarz_init(geom, basis) 50c 51c Get density matrix 52c 53 call hnd_prp_get_dens(rtdb,geom,basis,g_dens,ndens,scftyp, 54 & nclosed,nopen,nvirt) 55c 56c ----- calculate electric field gradient ----- 57c 58 if (ga_nodeid().eq.0) write(luout,9999) 59 if (ga_nodeid().eq.0) write(luout,9994) 60c 61 pi = acos(-1.0d0) 62 deg = 180.0d0/pi 63c 64 call ecce_print_module_entry('efg') 65c 66c 67c ----- define points for calculation ----- 68c 1. nuclei 69c 70 status=geom_ncent(geom,nat) 71c 72 if (.not. ma_push_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 73 & call errquit('hnd_efgmap: ma failed',911,MA_ERR) 74 if (.not. ma_push_get(mt_dbl,6*nat,'efg pnt',l_efgs,k_efgs)) 75 & call errquit('hnd_efgmap: ma failed',911,MA_ERR) 76 if (.not. ma_push_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 77 & call errquit('hnd_efgmap: ma failed',911,MA_ERR) 78c 79 do 30 iat=1,nat 80 status=geom_cent_get(geom,iat,at_tag,dbl_mb(k_xyzpt+3*(iat-1)), 81 & dbl_mb(k_zanpt+iat-1)) 82 30 continue 83c 84 call hnd_elfcon(basis,geom,g_dens(ndens),dbl_mb(k_xyzpt),nat, 85 & dbl_mb(k_efgs),2) 86c 87c get bq structures if any (MV) 88c ----------------------------- 89 dobq = .false. 90 if(geom_extbq_on()) then 91 dobq = .true. 92 bq_ncent = geom_extbq_ncenter() 93 i_cbq = geom_extbq_coord() 94 i_qbq = geom_extbq_charge() 95 end if 96c 97 docosmo = .false. 98 ncosbq = 0 99 if (rtdb_get(rtdb,'cosmo:nefc',mt_int,1,ncosbq).and.(ncosbq.gt.0)) 100 & docosmo = .true. 101c 102c ----- collect and output results of all points ----- 103c 104cc if (ga_nodeid().gt.0) goto 300 105c 106 if (docosmo) then 107 if (.not.rtdb_get(rtdb,'cosmo:nefc',mt_int,1,nefc)) 108 & call errquit('hnd_efgmap: rtdb get failed for nefc ',911, 109 & RTDB_ERR) 110 if (.not.ma_push_get(mt_dbl,nefc*3,'efcc',l_efcc,k_efcc)) 111 & call errquit('hnd_efgmap: malloc k_efcc fail',911,ma_err) 112 if (.not.ma_push_get(mt_dbl,nefc,'efcz',l_efcz,k_efcz)) 113 & call errquit('hnd_efgmap: malloc k_efcz fail',911,ma_err) 114 if (.not.rtdb_get(rtdb,'cosmo:efcc',mt_dbl,3*nefc, 115 & dbl_mb(k_efcc))) call 116 & errquit('hnd_efgmap: rtdb get failed efcc',912,rtdb_err) 117 if (.not.rtdb_get(rtdb,'cosmo:efcz',mt_dbl,nefc, 118 & dbl_mb(k_efcz))) call 119 & errquit('hnd_efgmap: rtdb get failed efcz',913,rtdb_err) 120 end if ! docosmo 121c 122 do 230 iat=1,nat 123 xp = dbl_mb(k_xyzpt +3*(iat-1)) 124 yp = dbl_mb(k_xyzpt+1+3*(iat-1)) 125 zp = dbl_mb(k_xyzpt+2+3*(iat-1)) 126c 127c ----- add nuclear contribution ----- 128c 129 efgxx = dbl_mb(k_efgs +6*(iat-1))/3.0d0 130 efgyy = dbl_mb(k_efgs+1+6*(iat-1))/3.0d0 131 efgzz = dbl_mb(k_efgs+2+6*(iat-1))/3.0d0 132 efgxy = dbl_mb(k_efgs+3+6*(iat-1))/3.0d0 133 efgxz = dbl_mb(k_efgs+4+6*(iat-1))/3.0d0 134 efgyz = dbl_mb(k_efgs+5+6*(iat-1))/3.0d0 135 do 210 i = 1,nat 136 xn = dbl_mb(k_xyzpt +3*(i-1)) - xp 137 yn = dbl_mb(k_xyzpt+1+3*(i-1)) - yp 138 zn = dbl_mb(k_xyzpt+2+3*(i-1)) - zp 139 zan = dbl_mb(k_zanpt+i-1) 140 rr = sqrt(xn*xn + yn*yn + zn*zn) 141 if (rr.lt.1.0d-3) go to 210 142 rr5=rr*rr*rr*rr*rr 143 efgxx = efgxx - zan*xn*xn/rr5 144 efgyy = efgyy - zan*yn*yn/rr5 145 efgzz = efgzz - zan*zn*zn/rr5 146 efgxy = efgxy - zan*xn*yn/rr5 147 efgxz = efgxz - zan*xn*zn/rr5 148 efgyz = efgyz - zan*yn*zn/rr5 149 210 continue 150c 151c ----- form -efc- contribution ----- 152c from cosmo point charges !!!! 153c 154 if (docosmo) then 155 do i = 1,nefc 156 xn = dbl_mb(k_efcc+3*(i-1) ) - xp 157 yn = dbl_mb(k_efcc+3*(i-1)+1) - yp 158 zn = dbl_mb(k_efcc+3*(i-1)+2) - zp 159 rr = sqrt(xn*xn + yn*yn + zn*zn) 160 if (rr.lt.1.0d-3) then 161 if (ga_nodeid().eq.0) write(luout,9993) xp,yp,zp,i 162 else 163 rr5=rr*rr*rr*rr*rr 164 efgxx = efgxx - dbl_mb(k_efcz+i-1)*xn*xn/rr5 165 efgyy = efgyy - dbl_mb(k_efcz+i-1)*yn*yn/rr5 166 efgzz = efgzz - dbl_mb(k_efcz+i-1)*zn*zn/rr5 167 efgxy = efgxy - dbl_mb(k_efcz+i-1)*xn*yn/rr5 168 efgxz = efgxz - dbl_mb(k_efcz+i-1)*xn*zn/rr5 169 efgyz = efgyz - dbl_mb(k_efcz+i-1)*yn*zn/rr5 170 endif 171 enddo 172 end if ! docosmo 173c 174c adding external bq contributions(MV) 175c ---------------------------------- 176 if (dobq) then 177 do i = 1,bq_ncent 178 xn = dbl_mb(i_cbq+3*(i-1) ) - xp 179 yn = dbl_mb(i_cbq+3*(i-1)+1) - yp 180 zn = dbl_mb(i_cbq+3*(i-1)+2) - zp 181 rr = sqrt(xn*xn + yn*yn + zn*zn) 182 if (rr.lt.1.0d-3) then 183 write(luout,9993) xp,yp,zp,i 184 else 185 rr5=rr*rr*rr*rr*rr 186 efgxx = efgxx - dbl_mb(i_qbq+i-1)*xn*xn/rr5 187 efgyy = efgyy - dbl_mb(i_qbq+i-1)*yn*yn/rr5 188 efgzz = efgzz - dbl_mb(i_qbq+i-1)*zn*zn/rr5 189 efgxy = efgxy - dbl_mb(i_qbq+i-1)*xn*yn/rr5 190 efgxz = efgxz - dbl_mb(i_qbq+i-1)*xn*zn/rr5 191 efgyz = efgyz - dbl_mb(i_qbq+i-1)*yn*zn/rr5 192 endif 193 end do 194 end if 195c 196 dbl_mb(k_efgs +6*(iat-1)) = 2.0d0*efgxx - efgyy - efgzz 197 dbl_mb(k_efgs+1+6*(iat-1)) = 2.0d0*efgyy - efgxx - efgzz 198 dbl_mb(k_efgs+2+6*(iat-1)) = 2.0d0*efgzz - efgxx - efgyy 199 dbl_mb(k_efgs+3+6*(iat-1)) = 3.0d0*efgxy 200 dbl_mb(k_efgs+4+6*(iat-1)) = 3.0d0*efgxz 201 dbl_mb(k_efgs+5+6*(iat-1)) = 3.0d0*efgyz 202c 203c ----- reorder into a as xx xy yy xz yz zz to form matrix ----- 204c 205 a(1) = dbl_mb(k_efgs +6*(iat-1)) 206 a(2) = dbl_mb(k_efgs+3+6*(iat-1)) 207 a(3) = dbl_mb(k_efgs+1+6*(iat-1)) 208 a(4) = dbl_mb(k_efgs+4+6*(iat-1)) 209 a(5) = dbl_mb(k_efgs+5+6*(iat-1)) 210 a(6) = dbl_mb(k_efgs+2+6*(iat-1)) 211 ij=0 212 do 241 i = 1, 3 213 do 241 j = 1, i 214 ij = ij + 1 215 vec(i,j) = a(ij) 216 vec(j,i) = a(ij) 217 241 continue 218c 219c ----- store ecce data ----- 220c 221 if (.not. geom_cent_tag(geom,iat,at_tag)) call 222 & errquit('hnd_efgmap: geom_cent_tag failed',0,GEOM_ERR) 223c geom_tag_to_element returns false for Bq elements (MV) 224c ----------------------------------------------------- 225 if (.not. geom_tag_to_element(at_tag,symbol,element,atn)) then 226 if(symbol.ne."bq") call 227 & errquit('hnd_efgmap: geom_tag_to_element failed',0,GEOM_ERR) 228 end if 229c 230c if (.not. geom_tag_to_element(at_tag,symbol,element,atn)) call 231c & errquit('hnd_efgmap: geom_tag_to_element failed',0,GEOM_ERR) 232 call ecce_print1_char('atom name',symbol,1) 233 call ecce_print2('EFG tensor',MT_DBL,vec,3,3,3) 234c 235c ----- print tensor components ----- 236c 237 if (ga_nodeid().eq.0) then 238 write(luout,9998) iat,symbol,xp,yp,zp 239 write(luout,9997) 240 write(luout,9995) (dbl_mb(k_efgs+6*(iat-1)+i),i=0,5) 241 end if 242c 243c ----- diagonalize to get principal components and vectors ----- 244c 245 call hnd_diag(vec,eig,3,.true.,.false.) 246 eta = abs( (eig(3)-eig(2)) / eig(1) ) 247c 248 call ecce_print1('EFG eigenvalues',MT_DBL,eig,3) 249 call ecce_print2('EFG eigenvectors',MT_DBL,vec,3,3,3) 250 call ecce_print1('EFG asymmetry',MT_DBL,eta,1) 251c 252 if (ga_nodeid().eq.0) then 253 write(luout,9992) 254 write(luout,9991) eig(1),eig(2),eig(3),eta 255 write(luout,9988) ((vec(i,j),j=1,3),i=1,3) 256 write(luout,*) ' ' 257 end if 258c 259 230 continue ! Assemblin and printing next atom 260c 261 if (docosmo) then 262 if (.not.ma_chop_stack(l_efcc)) 263 & call errquit('hnd_efgmap: chop stack l_efcc',913,ma_err) 264 end if ! docosmo 265c 266 call ecce_print_module_exit('EFG','ok') 267 call util_flush(luout) 268c 269c ----- release memory block ----- 270c 271 300 call ga_sync() 272c 273c ------- Deallocate MA memory ------ 274c 275 if (.not.ma_pop_stack(l_zanpt)) call errquit 276 & ('hnd_efgmap, ma_pop_stack of l_zanpt failed',911,MA_ERR) 277 if (.not.ma_pop_stack(l_efgs)) call errquit 278 & ('hnd_efgmap, ma_pop_stack of l_efgs failed',911,MA_ERR) 279 if (.not.ma_pop_stack(l_xyzpt)) call errquit 280 & ('hnd_efgmap, ma_pop_stack of l_xyzpt failed',911,MA_ERR) 281c 282 do i = 1, ndens 283 if (.not.ga_destroy(g_dens(i))) call 284 & errquit('efgmap: ga_destroy failed g_dens',0,GA_ERR) 285 enddo 286c 287c Terminate integrals 288c 289 call schwarz_tidy() 290 call int_terminate() 291c 292 return 293 9999 format(/,10x,23(1h-),/,10x,'Electric field gradient', 294 1 /,10x,23(1h-),/) 295 9998 format(/,1x,60(1h-),/,3x,'Atom',6x,'X',9x,'Y',9x,'Z',/,1x,60(1h-), 296 1 /,i5,1x,a2,3f10.5,/,1x,60(1h-),/) 297 9997 format(1x,'Electric field gradient in molecular frame (a.u.)',/, 298 2 9x,'XX',13x,'YY',13x,'ZZ',13x,'XY',13x,'XZ',13x,'YZ',/, 299 3 1x,90(1h-)) 300 9996 format(' --- Warning - electric field gradient at ', 301 1 3F10.5,' . contribution from nucleus ',i3,' ignored') 302 9995 format(1x,6f15.6,/) 303 9994 format(' 1 a.u. = 0.324123 10**(16) esu/cm**3 ', 304 1 ' ( or statvolts/cm**2 )',' = 0.97174 10**(22) v/m**2 ',/) 305 9993 format(' --- Warning - electric field gradient at ', 306 1 3f10.5,' . contribution from -efc- ',i3,' ignored') 307 9992 format(1x,'Principal components (a.u.) and orientation ', 308 1 /,' of principal axis w.r.t. absolute frame', 309 2 22x,'Asymmetry parameter eta',/,1x,86(1h-)) 310 9991 format(1x,3f15.6,16x,f15.6,/) 311 9988 format(1X,3F15.6) 312 end 313