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