1czora...
2czora...Scale zora eigenvalues and energy
3czora...
4          subroutine dft_zora_scale_so(
5     &                   rtdb,
6     &                   g_dens_ini,nexc,
7     &                   geom,
8     &                   ao_bas_han,
9     &                   nbf,
10     &                   nbf_ao,
11     &                   nbf_mo,
12     &                   g_dens,
13     &                   g_s,
14     &                   g_moso,
15     &                   g_zora_scale_sf,
16     &                   g_zora_scale_so,
17     &                   evals,
18     &                   focc,
19     &                   noc,
20     &                   nocc, ! FA
21     &                   ipol,
22     &                   ener_scal)
23
24       implicit none
25#include "errquit.fh"
26#include "mafdecls.fh"
27#include "stdio.fh"
28#include "global.fh"
29#include "msgids.fh"
30#include "zora.fh"
31#include "rtdb.fh"
32
33      integer  ga_create_atom_blocked
34      external ga_create_atom_blocked
35      external get_rhoS_so         ! FA
36      external hnd_efgmap_Z4_so    ! FA
37      external print_EFGZ4_version ! FA
38      integer g_dens(2)
39      integer g_moso(2)
40      integer g_orb(2)
41      integer g_dens_so(2)
42      integer g_scr(2)
43      integer g_s
44      integer g_zora_scale_sf(2)
45      integer g_zora_scale_so(3)
46      double precision numelecs, ener_kin
47      integer l_vecsre, k_vecsre
48      integer l_vecsim, k_vecsim
49      integer iorb
50      double precision eval_scal
51      double precision ener_scal
52      double precision ener_sf
53      double precision ener_so_x
54      double precision ener_so_y
55      double precision ener_so_z
56      double precision ener_tot
57      integer noc
58      integer ispin
59      integer ipol
60
61      integer geom
62      integer ao_bas_han
63      integer nbf
64      integer nbf_ao
65      integer nbf_mo
66
67      double precision focc(*)  ! occupation no.
68      double precision evals(*) ! eigenvalues
69      integer do_replnxn,mem_replnxn,val_in
70      integer l_evals_in, k_evals_in,
71     &          l_zora_scale_sf, k_zora_scale_sf,
72     &          l_zora_scale_so, k_zora_scale_so,
73     &          l_dens_so, k_dens_so,
74     &          l_vtmp, k_vtmp
75      double precision zora_denom          ! Added by FA
76      integer rtdb                         ! Added by FA
77      integer nexc,g_dens_ini(2)           ! Added by FA
78      integer j,j1,iorb1,nocc(2)           ! Added by FA
79      logical do_zgc_old                   ! Added by FA
80      integer ac_AB(2),iter                ! Added by FA
81      character*1  soxyz(3)                ! Added by FA
82      double precision ener_sf_1,ener_sf_2 ! Added by FA
83      logical status                       ! Added by FA
84      integer noelfgz4                     ! Added by FA
85      data noelfgZ4/1/                     ! Added by FA
86      soxyz(1)='z'
87      soxyz(2)='y'
88      soxyz(3)='x'
89c +++++++++++++++++++++++++++++++++++++++++++++++++++++
90      status=rtdb_get(rtdb,'prop:efieldgradZ4',MT_INT,1,noelfgZ4) ! FA
91      do_zgc_old=do_zora_get_correction
92      do_zora_get_correction= .true. !-1  ! FORCE-ZORA-CALC-INTS
93c +++++++++++++++++++++++++++++++++++++++++++++++++++++
94c ----- execute this code ONLY when EFGZ4 is requested -- START
95      if(noelfgZ4.eq.0) then
96       if (ga_nodeid().eq.0) call print_EFGZ4_version()
97       if (.not.ga_create(MT_DBL,2,nocc(1)+nocc(2),      ! FA
98     &      'dft_zora_scale: g_Ci',0,0,g_Ci))            ! FA
99     &  call errquit('dft_zora_scale_so: g_Ci',0,GA_ERR) ! FA
100       call ga_zero(g_Ci)                                ! FA
101      endif
102c ----- execute this code ONLY when EFGZ4 is requested -- END
103c     allocate memory
104
105      if (.not.MA_Push_Get(MT_Dbl, 2*nbf_mo, 'real vec aux',
106     &          l_vecsre, k_vecsre))
107     & call errquit('dft_zora_scale_so: cannot allocate vec',0, MA_ERR)
108c
109#if 0
110      if (.not.MA_Push_Get(MT_Dbl, nbf_mo, 'imag vec aux',
111     &          l_vecsim, k_vecsim))
112     & call errquit('dft_zora_scale_so: cannot allocate vec',0, MA_ERR)
113#else
114      k_vecsim=k_vecsre+nbf_mo
115#endif
116c
117c     spin-orbit vector - real
118      if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'orbs re',0,0,g_orb(1)))
119     & call errquit('dft_zora_scale_so: orb real error',0, GA_ERR)
120      call ga_zero(g_orb(1))
121c
122c     spin-orbit vector - imag
123      if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'orbs im',0,0,g_orb(2)))
124     & call errquit('dft_zora_scale_so: orb imag 2 error',0, GA_ERR)
125      call ga_zero(g_orb(2))
126c
127c     spin-orbit density matrix - real
128      if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'denmxre',0,0,
129     & g_dens_so(1)))
130     & call errquit('dft_zora_scale_so: dens real error',0, GA_ERR)
131      call ga_zero(g_dens_so(1))
132c
133c     spin-orbit density matrix - imag
134      if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'denmxim',0,0,
135     & g_dens_so(2)))
136     & call errquit('dft_zora_scale_so: dens imag error',0, GA_ERR)
137      call ga_zero(g_dens_so(2))
138c
139c     scratch array
140       if(.not.ga_duplicate(g_dens(1),g_scr(1),'scratch 1'))
141     &  call errquit('dft_zora_scale_so: ga_duplicate failed',1, GA_ERR)
142        call ga_zero(g_scr(1))
143       if(.not.ga_duplicate(g_dens(2),g_scr(2),'scratch 2'))
144     &  call errquit('dft_zora_scale_so: ga_duplicate failed',1, GA_ERR)
145        call ga_zero(g_scr(2))
146      j=0 ! FA
147      ac_AB(1)=1
148      ac_AB(2)=1
149      ener_scal = 0.d0
150cwhy??      ener_tot  = 0.d0
151c
152c     check if enough memory is available to allocate a few nxn objects
153c
154      do_replnxn=0
155      mem_replnxn=
156     C     nbf_mo*nbf_mo*2  +   ! dens_so
157     C     nbf_ao*nbf_ao    +   ! vtmp
158     C     nbf_ao*nbf_ao*2  +   ! zora_scale_sf
159     C     nbf_ao*nbf_ao*3  +   ! zora_scale_so
160     C     nbf_mo               ! evals_in
161      if (MA_inquire_avail(mt_dbl)*8/10.gt.mem_replnxn) then
162         do_replnxn=1
163      else
164         if(ga_nodeid().eq.0) then
165            write(6,*) ' mem needed    in MB',mem_replnxn*
166     / 8/1024/1024
167            write(6,*) ' mem available in MB ',MA_inquire_avail(mt_dbl)*
168     / 8/1024/1024
169         endif
170      endif
171      if(noelfgZ4.eq.0) then
172         write(6,*) ' need to compute elfgZ4, do_replnxn=0 '
173         do_replnxn=0
174      endif
175      if(rtdb_get(rtdb,'zora:do_replnxn',MT_INT,1,val_in))
176     c    do_replnxn=val_in
177      call ga_igop(19481984, do_replnxn, 1, 'min')
178      if(do_replnxn.eq.0) then
179         if(ga_nodeid().eq.0)
180     c        write(6,*) ' zoraso: NOT enough mem for repl'
181      do iorb = 1, noc
182       call ga_get(g_moso(1),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsre),1)
183       call ga_zero(g_orb(1))
184       call  ga_put(g_orb(1),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsre),1)
185       call ga_get(g_moso(2),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsim),1)
186       call ga_zero(g_orb(2))
187       call  ga_put(g_orb(2),1,nbf_mo,iorb,iorb,dbl_mb(k_vecsim),1)
188       call dft_densm_so(g_dens_so, g_orb, nbf_ao, noc)
189       call ga_zero(g_scr(1))
190       call ga_zero(g_scr(2))
191       call ga_dens_sf(g_scr, g_dens_so(1), nbf_ao)
192        ener_sf_1=ga_ddot(g_scr(1),g_zora_scale_sf(1))
193        ener_sf_2=ga_ddot(g_scr(2),g_zora_scale_sf(2))
194        ener_sf=ener_sf_1+ener_sf_2
195        ener_tot=ener_sf
196        do iter=1,3
197         call ga_zero(g_scr(1))
198         call ga_dens_so(g_scr(1),g_dens_so,nbf_ao,soxyz(iter))
199         ener_tot=ener_tot+
200     &            ga_ddot(g_scr(1),g_zora_scale_so(iter))
201        end do
202       zora_denom=1.d0+ener_tot
203       eval_scal = evals(iorb)
204       eval_scal = eval_scal/zora_denom
205       ener_scal = ener_scal - eval_scal*ener_tot
206       evals(iorb) = eval_scal
207c ----- execute this code ONLY when EFGZ4 is requested -- START
208        if(noelfgZ4.eq.0) then
209         j1=j+1                          ! FA
210         if (ac_AB(j1).le.nocc(j1)) then ! FA
211          call ga_fill_patch(g_Ci,j1,j1,ac_AB(j1),ac_AB(j1),
212     &                       zora_denom) ! FA
213          ac_AB(j1)=ac_AB(j1)+1          ! FA
214         end if                          ! FA
215        j=1-j                            ! FA
216        endif
217c ----- execute this code ONLY when EFGZ4 is requested -- END
218      enddo ! end-loop-iorb
219      else
220         if(ga_nodeid().eq.0)
221     c        write(6,*) ' zoraso: enough mem for repl'
222      if (.not.ma_push_get(mt_dbl, nbf_mo, 'evals_in',
223     &          l_evals_in, k_evals_in))
224     & call errquit('dft_zora_scale_so: cannot allocate evals',0,MA_ERR)
225      if (.not.ma_push_get(mt_dbl, nbf_ao*nbf_ao*2, 'zora_scale_sf',
226     &          l_zora_scale_sf, k_zora_scale_sf))
227     & call errquit('dft_zora_scale_so:cant all zora_scale_sf',0,MA_ERR)
228      if (.not.ma_push_get(mt_dbl, nbf_ao*nbf_ao*3, 'zora_scale_so',
229     &          l_zora_scale_so, k_zora_scale_so))
230     & call errquit('dft_zora_scale_so:cant all zora_scale_so',0,MA_ERR)
231      if (.not.ma_push_get(mt_dbl, nbf_mo*nbf_mo*2, 'dens_so',
232     &          l_dens_so, k_dens_so))
233     & call errquit('dft_zora_scale_so:cant all dens_so',0,MA_ERR)
234      if (.not.ma_push_get(mt_dbl, nbf_ao*nbf_ao, 'vtmp',
235     &          l_vtmp, k_vtmp))
236     & call errquit('dft_zora_scale_so:cant all vtmp',0,MA_ERR)
237c
238c     replicate zora_scale_sx
239c
240      call util_mygabcast(g_zora_scale_sf(1),
241     M  nbf_ao,nbf_ao, dbl_mb(k_zora_scale_sf),nbf_ao)
242      call util_mygabcast(g_zora_scale_sf(2),
243     M  nbf_ao,nbf_ao, dbl_mb(k_zora_scale_sf+nbf_ao*nbf_ao),nbf_ao)
244
245      call util_mygabcast(g_zora_scale_so(1),
246     M  nbf_ao,nbf_ao, dbl_mb(k_zora_scale_so),nbf_ao)
247      call util_mygabcast(g_zora_scale_so(2),
248     M  nbf_ao,nbf_ao, dbl_mb(k_zora_scale_so+nbf_ao*nbf_ao),nbf_ao)
249      call util_mygabcast(g_zora_scale_so(3),
250     M  nbf_ao,nbf_ao, dbl_mb(k_zora_scale_so+2*nbf_ao*nbf_ao),nbf_ao)
251
252c     copy input evals
253      call dcopy(nbf_mo,evals,1,dbl_mb(k_evals_in),1)
254c     zero evals [output] to get dgop working
255      call dcopy(noc,0d0,0,evals,1)
256      do iorb = ga_nodeid()+1, noc, ga_nnodes()
257      call dft_zora_nogaga(iorb,nbf_ao,nbf_mo,noc,g_moso,
258     C     soxyz,
259     D     dbl_mb(k_zora_scale_sf),dbl_mb(k_zora_scale_so),
260     E     ener_tot,
261     D     dbl_mb(k_dens_so),
262     O     dbl_mb(k_vtmp),dbl_mb(k_vecsre))
263       zora_denom=1.d0+ener_tot
264       eval_scal = dbl_mb(k_evals_in+iorb-1)/zora_denom
265       ener_scal = ener_scal - eval_scal*ener_tot
266       evals(iorb) = eval_scal
267#if 0
268NOT WORKING YET
269c ----- execute this code ONLY when EFGZ4 is requested -- START
270        if(noelfgZ4.eq.0) then
271         j1=j+1                          ! FA
272         if (ac_AB(j1).le.nocc(j1)) then ! FA
273          call ga_fill_patch(g_Ci,j1,j1,ac_AB(j1),ac_AB(j1),
274     &                       zora_denom) ! FA
275          ac_AB(j1)=ac_AB(j1)+1          ! FA
276         end if                          ! FA
277        j=1-j                            ! FA
278        endif
279#endif
280c ----- execute this code ONLY when EFGZ4 is requested -- END
281      enddo ! end-loop-iorb
282c dgop on evals
283      call ga_dgop(20166102, evals, noc, '+')
284      call ga_dgop(20166103, ener_scal, 1, '+')
285      endif
286c ----- execute this code ONLY when EFGZ4 is requested -- START
287       if(noelfgZ4.eq.0) then
288c       if (ga_nodeid().eq.0) write(*,*) "---g_Ci_so -- START"
289c       call ga_print(g_Ci)
290c       if (ga_nodeid().eq.0) write(*,*) "---g_Ci_so -- END"
291c       write(luout,*) "ener_scal:",ener_scal
292        call get_rhoS_so(rtdb,g_dens_ini,nexc,
293     &                   geom,ao_bas_han,
294     &                   nbf,nbf_ao,nbf_mo,
295     &                   g_dens,g_moso,
296     &                   noc,nocc,ipol)
297c -------- computing EFGs --------- START
298        call hnd_efgmap_Z4_so(rtdb,ao_bas_han,geom,
299     &                        nbf,nbf_ao,nbf_mo,
300     &                        g_dens,g_moso,
301     &                        noc,nocc)
302c -------- computing EFGs --------- END
303       endif
304c ----- execute this code ONLY when EFGZ4 is requested -- START
305
306c     deallocate memory
307#if 0
308      if (.not. ma_chop_stack(l_vecsim))
309     & call errquit('dft_zora_scale_so:l_vecsim', 0, MA_ERR)
310#endif
311      if (.not. ma_chop_stack(l_vecsre))
312     & call errquit('dft_zora_scale_so:l_vecsre', 0, MA_ERR)
313
314      if (.not. ga_destroy(g_orb(1)))
315     & call errquit('dft_zora_scale_so: ga_destroy failed',0, GA_ERR)
316      if (.not. ga_destroy(g_orb(2)))
317     & call errquit('dft_zora_scale_so: ga_destroy failed',0, GA_ERR)
318
319      if (.not. ga_destroy(g_dens_so(1)))
320     & call errquit('dft_zora_scale_so: ga_destroy failed',0, GA_ERR)
321      if (.not. ga_destroy(g_dens_so(2)))
322     & call errquit('dft_zora_scale_so: ga_destroy failed',0, GA_ERR)
323
324      if (.not. ga_destroy(g_scr(1)))
325     & call errquit('dft_zora_scale_so: ga_destroy failed',0, GA_ERR)
326      if (.not. ga_destroy(g_scr(2)))
327     & call errquit('dft_zora_scale_so: ga_destroy failed',0, GA_ERR)
328c +++++++++++++++ RESTORE VALUE +++++++++++++START
329      do_zora_get_correction=do_zgc_old
330c +++++++++++++++ RESTORE VALUE +++++++++++++END
331      return
332      end
333c
334czora...Inquire if the zora correction file is present
335c
336      logical function dft_zora_inquire_file_so(filename)
337c
338      implicit none
339c
340#include "errquit.fh"
341#include "global.fh"
342#include "tcgmsg.fh"
343#include "msgtypesf.h"
344#include "mafdecls.fh"
345#include "msgids.fh"
346#include "cscfps.fh"
347#include "inp.fh"
348#include "util.fh"
349#include "stdio.fh"
350c
351      character*(*) filename
352      logical found
353c
354      call ga_sync()
355c
356c     Inquire if file is present
357      inquire(file=filename,exist=found)
358      dft_zora_inquire_file_so = found
359c
360      call ga_sync()
361c
362      return
363      end
364c
365czora...Read in the zora atomic corrections from disk
366c
367      logical function dft_zora_read_so(filename, nbf, nsets, nmo,
368     &               mult, g_zora_sf, g_zora_scale_sf,
369     &               g_zora_so, g_zora_scale_so)
370c
371      implicit none
372c
373#include "errquit.fh"
374#include "global.fh"
375#include "tcgmsg.fh"
376#include "msgtypesf.h"
377#include "mafdecls.fh"
378#include "msgids.fh"
379#include "cscfps.fh"
380#include "inp.fh"
381#include "util.fh"
382#include "stdio.fh"
383c
384      character*(*) filename
385      integer iset              ! restricted or unrestricted
386c
387      integer g_zora_sf(2)
388      integer g_zora_scale_sf(2)
389      integer g_zora_so(3)
390      integer g_zora_scale_so(3)
391c
392      integer nsets             ! restricted or unrestricted
393      integer nbf               ! No. of functions in basis
394      integer nmo(nsets)
395      integer mult
396      integer ok, jset, i, j
397      integer l_zora, k_zora
398c
399      integer unitno
400      parameter (unitno = 77)
401      integer inntsize,ddblsize
402c
403      integer nsets_read
404      integer nbf_read
405      integer mult_read
406c
407      l_zora = -1               ! An invalid MA handle
408c
409      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
410      ddblsize=MA_sizeof(MT_DBL,1,MT_BYTE)
411c
412      call ga_sync()
413      ok = 0
414      if (ga_nodeid() .eq. 0) then
415c
416c     Print a message indicating the file being read
417      write(6,22) filename(1:inp_strlen(filename))
418 22   format(/' Read atomic ZORA corrections from ',a/)
419      call util_flush(luout)
420c
421c     Open the file
422      open(unitno, status='old', form='unformatted', file=filename,
423     $        err=1000)
424c
425c      Read in some basics to check if they are consistent with the calculation
426       read(unitno, err=1001, end=1001) nsets_read
427       read(unitno, err=1001, end=1001) nbf_read
428       read(unitno, err=1001, end=1001) mult_read
429c
430c      Error checks
431       if ((nsets_read .ne. nsets)
432     &  .or. (nbf_read .ne. nbf)
433     &  .or. (mult_read .ne. mult) ) goto 1003
434c
435c     Allocate the temporary buffer
436      if (.not.ma_push_get(mt_dbl,nbf,'dft_zora_read_so',l_zora,k_zora))
437     $        call errquit('dft_zora_read_so: ma failed', nbf, MA_ERR)
438c
439c     Read in g_zora_sf
440      do iset = 1,2  ! 2 components
441       do i = 1,nbf
442        call sread(unitno, dbl_mb(k_zora), nbf)
443        call ga_put(g_zora_sf(iset), 1, nbf, i, i, dbl_mb(k_zora), 1)
444       end do
445      end do
446c
447c     Read in g_zora_scale_sf
448      do iset = 1,2  ! 2 components
449       do i = 1,nbf
450        call sread(unitno, dbl_mb(k_zora), nbf)
451        call ga_put(g_zora_scale_sf(iset), 1, nbf, i, i,
452     &   dbl_mb(k_zora), 1)
453       end do
454      end do
455c
456c     Read in g_zora_so
457      do iset = 1,3  ! 3 components: x,y,z
458       do i = 1,nbf
459        call sread(unitno, dbl_mb(k_zora), nbf)
460        call ga_put(g_zora_so(iset), 1, nbf, i, i, dbl_mb(k_zora), 1)
461       end do
462      end do
463c
464c     Read in g_zora_scale_so
465      do iset = 1,3  ! 3 components: x,y,z
466       do i = 1,nbf
467        call sread(unitno, dbl_mb(k_zora), nbf)
468        call ga_put(g_zora_scale_so(iset), 1, nbf, i, i,
469     &   dbl_mb(k_zora), 1)
470       end do
471      end do
472c
473c     Close the file
474      close(unitno,err=1002)
475      ok = 1
476c
477c     Deallocate the temporary buffer
478      if (.not. ma_pop_stack(l_zora)) call errquit
479     $      ('dft_zora_read_so: pop failed', l_zora, MA_ERR)
480c
481      end if
482c
483c     Broadcast status to other nodes
484 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
485      call ga_sync()
486c
487      dft_zora_read_so = ok .eq. 1
488      return
489c
490 1000 write(6,*) 'dft_zora_read_so: failed to open ',
491     $     filename(1:inp_strlen(filename))
492      call util_flush(luout)
493      ok = 0
494      goto 10
495c
496 1001 write(6,*) 'dft_zora_read_so: failed to read ',
497     $     filename(1:inp_strlen(filename))
498      call util_flush(luout)
499      ok = 0
500      close(unitno,err=1002)
501      goto 10
502c
503 1003 write(6,*) 'dft_zora_read_so: file inconsistent with calculation',
504     $     filename(1:inp_strlen(filename))
505      call util_flush(luout)
506      ok = 0
507      close(unitno,err=1002)
508      goto 10
509c
510 1002 write(6,*) 'dft_zora_read_so: failed to close',
511     $     filename(1:inp_strlen(filename))
512      call util_flush(luout)
513      ok = 0
514      goto 10
515c
516      end
517c
518czora...Write out the zora atomic corrections to disk
519c
520      logical function dft_zora_write_so(rtdb, basis, filename,
521     &     nbf, nsets, nmo, mult,
522     &     g_zora_sf, g_zora_scale_sf,
523     &     g_zora_so, g_zora_scale_so)
524c
525      implicit none
526c
527#include "errquit.fh"
528#include "mafdecls.fh"
529#include "global.fh"
530#include "tcgmsg.fh"
531#include "msgtypesf.h"
532#include "inp.fh"
533#include "msgids.fh"
534#include "cscfps.fh"
535#include "util.fh"
536#include "bas.fh"
537#include "geom.fh"
538#include "rtdb.fh"
539#include "stdio.fh"
540c
541c     Temporary routine
542c
543      integer rtdb              ! [input] RTDB handle (-1 if not accessible)
544      integer basis             ! [input] Basis handle(-1 if not accessible)
545      character*(*) filename    ! [input] File to write to
546      integer nbf               ! [input] No. of functions in basis
547      integer nsets             ! [input] No. of sets of matrices
548      integer nmo(nsets)        ! [input] No. of mos in each set
549      integer mult
550c
551      integer g_zora_sf(2)
552      integer g_zora_scale_sf(2)
553      integer g_zora_so(3)
554      integer g_zora_scale_so(3)
555c
556      integer unitno
557      parameter (unitno = 77)
558c
559      integer lentit
560      integer lenbas
561      integer l_zora, k_zora
562      integer ok, iset, i, j
563      integer geom, ma_type, nelem
564      character*26 date
565      character*32 geomsum, basissum, key
566      character*20 scftype20
567      character*128 basis_name, trans_name
568      double precision energy, enrep
569      integer inntsize
570c
571      l_zora = -1               ! An invalid MA handle
572c
573      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
574      call ga_sync()
575c
576      ok = 0
577c
578c     Read routines should be consistent with this
579c
580c     Write out the atomic zora corrections
581c
582      if (ga_nodeid() .eq. 0) then
583c
584c     Open the file
585      open(unitno, status='unknown', form='unformatted',
586     $        file=filename, err=1000)
587c
588c     Write out the number of sets and basis functions
589      write(unitno, err=1001) nsets
590      write(unitno, err=1001) nbf
591      write(unitno, err=1001) mult
592c
593c     Allocate the temporary buffer
594      if (.not.ma_push_get(mt_dbl,nbf,'dft_zora_write_so',
595     & l_zora,k_zora))
596     &   call errquit('dft_zora_write_so: ma failed', nbf, MA_ERR)
597c
598c     Write out g_zora_sf
599      do iset = 1,2  ! 2 components
600       do i = 1, nbf
601        call ga_get(g_zora_sf(iset), 1, nbf, i, i, dbl_mb(k_zora),1)
602        call swrite(unitno, dbl_mb(k_zora), nbf)
603       end do
604      end do
605c
606c     Write out g_zora_scale_sf
607      do iset = 1,2  ! 2 components
608       do i = 1, nbf
609        call ga_get(g_zora_scale_sf(iset), 1, nbf, i, i,
610     &              dbl_mb(k_zora),1)
611        call swrite(unitno, dbl_mb(k_zora), nbf)
612       end do
613      end do
614c
615c     Write out g_zora_so
616      do iset = 1,3  ! 3 components: x,y,z
617       do i = 1, nbf
618        call ga_get(g_zora_so(iset), 1, nbf, i, i, dbl_mb(k_zora),1)
619        call swrite(unitno, dbl_mb(k_zora), nbf)
620       end do
621      end do
622c
623c     Write out g_zora_scale_so
624      do iset = 1,3  ! 3 components: x,y,z
625       do i = 1, nbf
626        call ga_get(g_zora_scale_so(iset), 1, nbf, i, i,
627     &              dbl_mb(k_zora),1)
628        call swrite(unitno, dbl_mb(k_zora), nbf)
629       end do
630      end do
631c
632c     Deallocate the temporary buffer
633      if (.not. ma_pop_stack(l_zora))
634     $  call errquit('dft_zora_write_so: ma pop failed', l_zora, MA_ERR)
635c
636c     Close the file
637      close(unitno,err=1002)
638c
639      ok = 1
640      end if
641c
642c     Broadcast status to other nodes
643 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
644      call ga_sync()
645c
646      dft_zora_write_so = (ok .eq. 1)
647      if (ga_nodeid() .eq. 0) then
648         write(6,22) filename(1:inp_strlen(filename))
649 22      format(/' Wrote atomic ZORA corrections to ',a/)
650         call util_flush(luout)
651      endif
652c
653      return
654c
655 1000 write(6,*) 'dft_zora_write_so: failed to open ',
656     $     filename(1:inp_strlen(filename))
657      call util_flush(luout)
658      ok = 0
659      goto 10
660c
661 1001 write(6,*) 'dft_zora_write_so: failed to write ',
662     $     filename(1:inp_strlen(filename))
663      call util_flush(luout)
664      ok = 0
665      close(unitno,err=1002)
666      goto 10
667c
668 1002 write(6,*) 'dft_zora_write_so: failed to close',
669     $     filename(1:inp_strlen(filename))
670      call util_flush(luout)
671      ok = 0
672      goto 10
673c
674      end
675c $Id$
676      subroutine dft_zora_nogaga(iorb,nbf_ao,nbf_mo,noc,g_moso,
677     C     soxyz,
678     D     zora_scale_sf,zora_scale_so,
679     E     ener_tot,
680     D     dens_so,
681     O     vtmp,vecs)
682      implicit none
683      integer iorb
684      integer nbf_mo,nbf_ao ! [in] nbf_mo=2*nbf_ao
685      integer noc
686      double precision zora_scale_sf(nbf_ao,nbf_ao,2)! replicated[in]
687      double precision zora_scale_so(nbf_ao,nbf_ao,3)! replicated[in]
688      integer g_moso(2) ! [in]
689      double precision vecs(nbf_mo,2) ! [scratch]
690      double precision dens_so(nbf_mo,nbf_mo,2) ! repl [scratch]
691      double precision vtmp(nbf_ao,nbf_ao) ! repl [scratch]
692      double precision ener_tot ! [out]
693      character*1  soxyz(3)
694
695c
696      integer xyz
697      integer j,i
698c
699      external ddot
700      double precision ddot
701c
702
703       call ga_get(g_moso(1),1,nbf_mo,iorb,iorb,vecs(1,1),1)
704       call ga_get(g_moso(2),1,nbf_mo,iorb,iorb,vecs(1,2),1)
705       call dftv_densm_so(dens_so,vecs,nbf_mo,noc)
706c        ener_sf_1
707c strided copy
708        do j=1,nbf_ao
709           call dcopy(nbf_ao,dens_so(1,j,1),1,vtmp(1,j),1)
710        enddo
711        ener_tot=ddot(nbf_ao*nbf_ao,
712     D      vtmp,1,zora_scale_sf(1,1,1),1)
713cstrided copy
714        do j=1+nbf_ao,2*nbf_ao
715           call dcopy(nbf_ao,dens_so(1+nbf_ao,j,1),1,vtmp(1,j-nbf_ao),1)
716        enddo
717c        ener_sf_2
718        ener_tot=ener_tot+ddot(nbf_ao*nbf_ao,
719     D      vtmp,1,zora_scale_sf(1,1,2),1)
720        do xyz=1,3
721         call vec_dens_so(vtmp,dens_so,nbf_ao,soxyz(xyz))
722         ener_tot=ener_tot+ddot(nbf_ao*nbf_ao,
723     &            vtmp,1,zora_scale_so(1,1,xyz),1)
724        end do
725        return
726        end
727