← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/nwdft/zora/dft_zora_utils.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwdft/zora/dft_zora_utils.F
 
1
c
 
2
czora...Scale zora eigenvalues and energy
 
3
          subroutine dft_zora_scale(
 
4
     &                   rtdb,g_dens_at,nexc, ! FA
 
5
     &                   geom, 
 
6
     &                   ao_bas_han,
 
7
     &                   nbf,
 
8
     &                   nbf_ao,
 
9
     &                   g_dens,
 
10
     &                   g_s,
 
11
     &                   g_movecs,
 
12
     &                   g_zora_scale_sf,
 
13
     &                   evals,
 
14
     &                   focc,
 
15
     &                   noc,
 
16
     &                   ipol,
 
17
     &                   ener_scal)
 
18
 
 
19
       implicit none
 
20
#include "errquit.fh"
 
21
#include "mafdecls.fh"
 
22
#include "stdio.fh"
 
23
#include "global.fh"
 
24
#include "msgids.fh"
 
25
#include "zora.fh" 
 
26
#include "rtdb.fh" 
 
27
#include "geom.fh" 
 
28
      integer rtdb,geom,ao_bas_han ! handles
 
29
      integer  ga_create_atom_blocked
 
30
      external ga_create_atom_blocked
 
31
      external get_rhoS,print_EFGZ4_version,
 
32
     &         get_EPRg,print_EPRg_version,
 
33
     &         get_NMRCS_SRZORA,print_NMRCS_SRZORA_version,
 
34
     &         get_NMRHFine_ZORA,print_NMRHypFine_version,
 
35
     &         get_densZ4,
 
36
     &         dft_zorashield_read,dft_zorashield_write      
 
37
      integer dft_zorashield_read,dft_zorashield_write
 
38
      integer g_dens(2),g_movecs(2),g_zora_scale_sf(2)
 
39
      integer g_s,ga_Fji,g_orb,g_dens_sf,g_zora_scr(2),
 
40
     &        g_dens_at(2) 
 
41
      integer l_vecs,k_vecs
 
42
      integer noc(2),iorb,ispin,ipol,i,nexc,noc1
 
43
      integer nbf,nbf_ao,nat_slc  
 
44
      double precision eval_scal,ener_scal,scf_dbl
 
45
      double precision focc(nbf*ipol)  ! occupation no.
 
46
      double precision evals(ipol*nbf) ! eigenvalues                     
 
47
      double precision zora_eint,zora_denom,
 
48
     &                 scale,scale_MO,scale_MO1  
 
49
      logical status,done_Fji,do_zgc_old    
 
50
      integer nospin,nogshift,nogiao,noelfgz4,efgfile,switch_focc
 
51
      data nospin,nogshift,nogiao,noelfgZ4/1,1,1,1/                       
 
52
      data switch_focc/0/ ! =0 not using occ keyword =1 using occ keyword
 
53
      data efgfile/0/     ! =0 no NLMO analysis      =1 doing NLMO analysis
 
54
      character*255 zorafilename
 
55
      integer dims(3),chunk(3) ! for g_Cifull(i) i=1,ipol
 
56
      integer g_densZ4(3),ind_occ
 
57
      logical skip_csAOev,skip_gshiftAOev,
 
58
     &        skip_hypAOev,skip_efgz4AOev
 
59
c ---------- nlmo definitions ----------------- START
 
60
      integer g_munuV6,g_munu_rhoS,g_densx,       
 
61
     &        g_zora_scale_munu(2),ipolmunu 
 
62
      integer ndir
 
63
      logical dft_zoraEFGZ4_NLMOAnalysis_write
 
64
      external dft_zoraEFGZ4_NLMOAnalysis_write
 
65
     &         util_file_name
 
66
c ---------- nlmo definitions ----------------- END
 
67
      status=rtdb_get(rtdb,'prop:efgfile',mt_int,1,efgfile) ! for NLMO analysis
 
68
      status=rtdb_get(rtdb,'focc:occ-switch',mt_int,1,switch_focc)  ! check-occupations-keyword
 
69
      status=rtdb_get(rtdb,'prop:efieldgradZ4',MT_INT,1,noelfgZ4)   
 
70
      status=rtdb_get(rtdb,'prop:giao'        ,MT_INT,1,nogiao)     
 
71
      status=rtdb_get(rtdb,'prop:gshift'      ,MT_INT,1,nogshift)   
 
72
      status=rtdb_get(rtdb,'prop:hyperfine'   ,MT_INT,1,nospin)     
 
73
      do_zgc_old=do_zora_get_correction                             
 
74
      do_zora_get_correction=.true. ! -1  ! FORCE-ZORA-CALC-INTS
 
75
c
 
76
      if (ipol.eq.1) then
 
77
        scf_dbl=2.0d00
 
78
        noc1=noc(1)
 
79
      else if (ipol.eq.2) then 
 
80
        scf_dbl=1.0d00   
 
81
        noc1=noc(1)+noc(2)    
 
82
      endif
 
83
c
 
84
      if (noelfgZ4.eq.0 .or. nospin.eq.0 .or. 
 
85
     &    nogshift.eq.0 .or. nogiao.eq.0) then
 
86
        if (ga_nodeid().eq.0) then
 
87
          write(LuOut,*)
 
88
          call util_print_centered(LuOut,
 
89
     $   'Commencing ZORA Property Calculations', 23, .true.)
 
90
          write(LuOut,*)
 
91
        end if
 
92
      end if
 
93
c
 
94
c      if (nospin.eq.0 .and. ga_nodeid().eq.0)
 
95
c     &   call print_NMRHypFine_version()
 
96
c      if (nogiao.eq.0 .and. ga_nodeid().eq.0)
 
97
c     &   call print_NMRCS_SRZORA_version()
 
98
c      if (nogshift.eq.0 .and. ga_nodeid().eq.0)
 
99
c     &   call print_EPRg_version()
 
100
c      if (noelfgZ4.eq.0 .and. ga_nodeid().eq.0) 
 
101
c     &   call print_EFGZ4_version()
 
102
c
 
103
      if (noelfgZ4.eq.0 .or. nospin.eq.0 .or. 
 
104
     &    nogshift.eq.0 .or. nogiao.eq.0) then
 
105
       if (.not.ga_create(MT_DBL,ipol,noc1,'dft_zora_scale: g_Ci',
 
106
     $       0,0,g_Ci)) call errquit('dft_zora_scale: g_Ci',0,
 
107
     &       GA_ERR)
 
108
       call ga_zero(g_Ci)
 
109
      endif
 
110
c
 
111
      if (nospin  .eq.0 .or. nogiao  .eq.0 .or.
 
112
     &    nogshift.eq.0 .or. noelfgZ4.eq.0) then
 
113
c ------- create g_Cifull for scaling MO vectors
 
114
       dims(1) =nbf
 
115
       chunk(1)=nbf 
 
116
       do ispin=1,ipol
 
117
        if (.not. nga_create(mt_dbl,1,dims,"g_Cifull",
 
118
     &       chunk,g_Cifull(ispin)))
 
119
     $       call errquit('dft_zora_scale: g_Cifull', 0,
 
120
     &                    GA_ERR)
 
121
        call ga_fill(g_Cifull(ispin),1.0d0)
 
122
       enddo ! end-loop-ispin
 
123
      endif   
 
124
c -- execute this code for EFGSRZ4+NLMO analysis is requested -- START
 
125
        if (noelfgZ4.eq.0 .and. efgfile.eq.1) then ! -- if-noelfgZ4.eq.0---- START
 
126
         ipolmunu=ipol
 
127
         if (.not. ga_create(mt_dbl,nbf,nbf,
 
128
     &        'hnd_elfcon_symm: g_munu',0,0,g_zora_scale_munu(1)))
 
129
     $    call errquit('hnd_elfcon_symm:',0,GA_ERR)
 
130
         call ga_copy(g_zora_scale_sf(1),g_zora_scale_munu(1))
 
131
         if(ipol.gt.1) then
 
132
          if (.not. ga_create(mt_dbl,nbf,nbf,
 
133
     &        'hnd_elfcon_symm: g_munu',0,0,g_zora_scale_munu(2)))
 
134
     $    call errquit('hnd_elfcon_symm:',0,GA_ERR)
 
135
          call ga_copy(g_zora_scale_sf(2),g_zora_scale_munu(2))
 
136
         endif
 
137
c ----- Store efgz4 data in a file ------- START
 
138
c       Note.- lbl_nlmo defined in zora.fh
 
139
        ndir=6
 
140
        call util_file_name(lbl_nlmo,.false.,.false.,zorafilename)
 
141
c ---------> Write NMLO analysis data: 1 of 3 ----- START
 
142
        if (.not. dft_zoraEFGZ4_NLMOAnalysis_write(
 
143
     &       zorafilename, ! in: filename
 
144
     &                nbf, ! in: nr basis functions
 
145
     &               ndir, ! in: nr of directions: 6 = xx yy zz xy xz yz
 
146
     &              nlist, ! in: list of selected atoms 
 
147
     &                  1, ! in: writing order =1,2,3
 
148
     &           ipolmunu, ! in: write for ndada=1
 
149
     &  g_zora_scale_munu, ! in: write for ndada=1
 
150
     &        g_munu_rhoS, ! in: write for ndata=2
 
151
     &            g_densx, ! in: write for ndata=2
 
152
     &           g_munuV6))
 
153
     &     call errquit('get_rhoS: dft_zoraNLMO_write failed',
 
154
     &                  0,DISK_ERR)
 
155
c ---------> Write NMLO analysis data: 1 of 3 ----- END
 
156
        if (.not. ga_destroy(g_zora_scale_munu(1))) call errquit(
 
157
     &    'dft_zora_utils: ga_destroy failed ',0, GA_ERR)
 
158
         if (ipolmunu.gt.1) then
 
159
         if (.not. ga_destroy(g_zora_scale_munu(2))) call errquit(
 
160
     &    'dft_zora_utils: ga_destroy failed ',0, GA_ERR)
 
161
         endif
 
162
        endif ! -------- if-noelfgZ4.eq.0---- END
 
163
c -- execute this code for EFGSRZ4+NLMO analysis is requested -- END
 
164
      if (.not. ma_alloc_get(mt_dbl,nbf_ao,'vec aux',l_vecs, k_vecs))
 
165
     &    call errquit(
 
166
     &          'dft_zora_utils: cannot allocate vec aux',911,MA_ERR)
 
167
      g_orb     = ga_create_atom_blocked(geom, ao_bas_han, 'orbs')
 
168
      g_dens_sf = ga_create_atom_blocked(geom, ao_bas_han,'orbs dens')
 
169
      do i=1,3 ! for scaled density matrix used in NMR-shieldigns,hyperfine,gshift,EFGZ4
 
170
       g_densZ4(i)=ga_create_atom_blocked(geom,ao_bas_han,'orbs dens')
 
171
       call ga_zero(g_densZ4(i))
 
172
      end do
 
173
c
 
174
c     scale the eigenvalues and energy
 
175
      ener_scal = 0.d0
 
176
      do ispin = 1,ipol
 
177
       do iorb = 1,nbf
 
178
        call ga_get(g_movecs(ispin), 1, nbf, iorb, iorb,
 
179
     &              dbl_mb(k_vecs), 1)
 
180
        call ga_zero(g_orb)
 
181
        call ga_put(g_orb, 1, nbf, iorb, iorb,
 
182
     &              dbl_mb(k_vecs), 1)
 
183
        call ga_zero(g_dens_sf)
 
184
        call ga_dgemm('n','t',nbf,nbf,noc(ispin),
 
185
     &                1.0d00,g_orb,g_orb,
 
186
     &                0.0d00,g_dens_sf)
 
187
        eval_scal = evals(iorb)
 
188
        if (ispin.gt.1) eval_scal = evals(iorb + nbf)
 
189
        zora_eint = ga_ddot(g_dens_sf,g_zora_scale_sf(ispin)) ! FA
 
190
        if ((do_NonRel) .or. (not_zora_scale)) then
 
191
c +++++++++ Create non-scaled density matrix +++++++ START
 
192
         zora_denom=1.0d0
 
193
         ener_scal = ener_scal-eval_scal ! for skipping scaling
 
194
         if (iorb .le. noc(ispin)) then        
 
195
          scale_MO=scf_dbl
 
196
          if (switch_focc.eq.1) then ! adding focc()
 
197
            ind_occ=nbf*(ispin-1)+iorb
 
198
            scale_MO=scf_dbl*focc(ind_occ)
 
199
          endif
 
200
          call ga_add(1.0d0   ,g_densZ4(ispin),
 
201
     &                scale_MO,g_dens_sf,
 
202
     &                g_densZ4(ispin))
 
203
         endif
 
204
c +++++++++ Create non-scaled density matrix +++++++ END
 
205
        else
 
206
c +++++++++ Create scaled density matrix +++++++ START
 
207
         zora_denom= 1.d0+zora_eint                            
 
208
         eval_scal = eval_scal/zora_denom                     
 
209
         ener_scal = ener_scal-eval_scal*zora_eint
 
210
         if (iorb .le. noc(ispin)) then
 
211
          scale=1.0d0/zora_denom
 
212
          call ga_scale(g_dens_sf,scale)
 
213
          scale_MO=scf_dbl
 
214
          if (switch_focc.eq.1) then ! adding focc()
 
215
            ind_occ=nbf*(ispin-1)+iorb
 
216
            scale_MO=scf_dbl*focc(ind_occ)
 
217
          endif
 
218
          call ga_add(1.0d0   ,g_densZ4(ispin),
 
219
     &                scale_MO,g_dens_sf,
 
220
     &                g_densZ4(ispin))
 
221
         endif
 
222
c +++++++++ Create scaled density matrix +++++++ END
 
223
        endif
 
224
        if (ispin.le.1) then
 
225
            evals(iorb)     = eval_scal
 
226
        else
 
227
            evals(iorb+nbf) = eval_scal
 
228
        end if
 
229
c ----- execute this code ONLY when EFGZ4/NMRCSZ4/EPRgshift is requested -- START
 
230
         if(noelfgZ4.eq.0 .or. nospin.eq.0 .or. 
 
231
     &      nogshift.eq.0 .or. nogiao.eq.0) then
 
232
          if ((ipol.eq.1 .and. iorb.le.noc(1)).or.
 
233
     &        (ipol.eq.2 .and.
 
234
     &      (ispin.eq.1 .and. iorb.le.noc(1)).or.
 
235
     &      (ispin.eq.2 .and. iorb.le.noc(2)))) then
 
236
           call ga_fill_patch(g_Ci,ispin,ispin,iorb,iorb,
 
237
     &                        zora_denom)
 
238
          endif
 
239
         endif
 
240
         if(nogiao  .eq.0 .or.
 
241
     &      nogshift.eq.0 .or. nospin.eq.0) then
 
242
          if ((ipol.eq.1 .and. iorb.le.noc(1)).or.
 
243
     &        (ipol.eq.2 .and.
 
244
     &      (ispin.eq.1 .and. iorb.le.noc(1)).or.
 
245
     &      (ispin.eq.2 .and. iorb.le.noc(2)))) then
 
246
            if (nogiao.eq.0 .or. nogshift.eq.0 .or.
 
247
     &          nospin .eq. 0) then
 
248
             scale_MO1=1.0d0/zora_denom
 
249
             if (switch_focc.eq.1) then ! adding focc()
 
250
              ind_occ=nbf*(ispin-1)+iorb
 
251
              scale_MO1=1.0d0/zora_denom*focc(ind_occ)
 
252
c              if (ga_nodeid().eq.0) then
 
253
c               write(*,5) 1.0d0/zora_denom,ind_occ,focc(ind_occ)
 
254
c5              format('In dft_zora_utils:: (scl_MO,ind_occ,focc)=(',
 
255
c    &                  f15.8,',',i4,',',f15.8,')') 
 
256
c              endif
 
257
             endif     
 
258
             call ga_fill_patch(g_Cifull(ispin),iorb,iorb,1,1,
 
259
     &                          scale_MO1)
 
260
            endif
 
261
          end if
 
262
         endif
 
263
c ----- execute this code ONLY when EFGZ4/NMRCSZ4/EPRgshift is requested -- END
 
264
       end do ! orbital loop
 
265
      end do ! polarization loop
 
266
c--- Create total scaled-density matrix ---- START
 
267
      call ga_zero(g_densZ4(3))
 
268
      call ga_copy(g_densZ4(1),g_densZ4(3))
 
269
      if (ipol.gt.1) then 
 
270
       call ga_add(1.0d00,g_densZ4(1),
 
271
     &             1.0d00,g_densZ4(2),
 
272
     &                    g_densZ4(3))      
 
273
      end if
 
274
c--- Create total scaled-density matrix ---- END
 
275
c     double the energy for closed-shell calculations
 
276
      if (ipol.eq.1) ener_scal = 2.d0*ener_scal     
 
277
c      write(luout,*) "scaled_energy:",ener_scal
 
278
c ----- execute this code ONLY when NMR-CS-Z4 is requested -- START
 
279
      done_Fji=.false.
 
280
      if(.not.rtdb_get(rtdb,'zora:skip_csAOev',
 
281
     &                 mt_log,1,skip_csAOev))        
 
282
     &  skip_csAOev = .false.       
 
283
      if (nogiao.eq.0 .and. .not.(skip_csAOev)) then
 
284
       call get_NMRCS_SRZORA(rtdb,g_dens_at,nexc,
 
285
     &                       geom,ao_bas_han,
 
286
     &                       nbf,noc,ipol,
 
287
     &                       ga_Fji,
 
288
     &                       g_densZ4)
 
289
       done_Fji=.true. ! I do not need to calculate again 
 
290
                       ! if calling get_EPRg
 
291
      endif
 
292
c ----- execute this code ONLY when NMR-CS-Z4 is requested -- END
 
293
c ----- execute this code ONLY when NMR-spinspin is requested -- START
 
294
      if(.not.rtdb_get(rtdb,'zora:skip_hypAOev',
 
295
     &                 mt_log,1,skip_hypAOev))        
 
296
     &  skip_hypAOev = .false.    
 
297
      if (nospin.eq.0 .and. .not.(skip_hypAOev)) then
 
298
        call get_NMRHFine_ZORA(rtdb,g_dens_at,nexc,
 
299
     &                         geom,ao_bas_han,
 
300
     &                         nbf,focc,noc,ipol,
 
301
     &                         g_densZ4) 
 
302
      endif
 
303
c ----- execute this code ONLY when NMR-spinspin is requested -- END
 
304
c ----- execute this code ONLY when EPR-gshift is requested -- START
 
305
      skip_gshiftAOev = .false. 
 
306
      if(.not.rtdb_get(rtdb,'zora:skip_gshiftAOev',
 
307
     &                 mt_log,1,skip_gshiftAOev))        
 
308
     &  skip_gshiftAOev = .false.  
 
309
      
 
310
      if (ga_nodeid().eq.0)
 
311
     &  write(*,22) nogshift,skip_gshiftAOev,done_Fji
 
312
 22   format('(nogshift,skip_gshiftAOev,done_Fji)=(',
 
313
     &        i2,',',L,',',L,')')
 
314
 
 
315
      if (nogshift.eq.0 .and. .not.(skip_gshiftAOev)) then
 
316
       call get_EPRg(rtdb,g_dens_at,nexc,
 
317
     &               geom,ao_bas_han,
 
318
     &               nbf,focc,noc,ipol,
 
319
     &               done_Fji,ga_Fji,
 
320
     &               g_densZ4)
 
321
      endif   
 
322
c ----- execute this code ONLY when EPR-gshift is requested -- END
 
323
      if ((nogiao  .eq.0 .and. .not.(skip_csAOev)).or. 
 
324
     &    (nogshift.eq.0 .and. .not.(skip_gshiftAOev))) then
 
325
        if (.not.ga_destroy(ga_Fji)) call 
 
326
     &    errquit('dft_zora_utils: ga_destroy failed ga_Fji',
 
327
     &            0,GA_ERR)
 
328
      endif
 
329
c ----- execute this code ONLY when EFGZ4 is requested -- START
 
330
      if(.not.rtdb_get(rtdb,'zora:skip_efgz4AOev',
 
331
     &                 mt_log,1,skip_efgz4AOev))        
 
332
     &  skip_efgz4AOev = .false.  
 
333
      if (noelfgZ4.eq.0 .and. .not.(skip_efgz4AOev)) then
 
334
       call get_rhoS(rtdb,g_dens_at,nexc,
 
335
     &               geom,ao_bas_han,nbf,nbf_ao,
 
336
     &               noc,ipol)     
 
337
      endif 
 
338
c ----- execute this code ONLY when EFGZ4 is requested -- END
 
339
c     deallocate memory
 
340
      if (.not.ma_free_heap(l_vecs)) call errquit
 
341
     &   ('dft_zora_utils, ma_free_heap of l_vecs failed',911,MA_ERR)
 
342
      if (.not. ga_destroy(g_orb)) call errquit(
 
343
     &  'zora_scale_evals: ga_destroy failed ',0, GA_ERR)
 
344
      if (.not. ga_destroy(g_dens_sf)) call errquit(
 
345
     &  'zora_scale_evals: ga_destroy failed ',0, GA_ERR)
 
346
      do i=1,ipol
 
347
        if (.not. ga_destroy(g_densZ4(i))) call errquit(
 
348
     &    'dft_zora_utils: ga_destroy failed ',0, GA_ERR)
 
349
      enddo
 
350
      if (ipol.eq.2) then
 
351
        if (.not. ga_destroy(g_densZ4(3))) call errquit(
 
352
     &    'dft_zora_utils: ga_destroy failed ',0, GA_ERR)
 
353
      endif
 
354
c
 
355
      do_zora_get_correction=do_zgc_old   ! restore value
 
356
c
 
357
      return
 
358
      end
 
359
c
 
360
czora...Inquire if the zora correction file is present
 
361
c
 
362
      logical function dft_zora_inquire_file(filename)
 
363
c
 
364
      implicit none
 
365
c
 
366
#include "errquit.fh"
 
367
#include "global.fh"
 
368
#include "tcgmsg.fh"
 
369
#include "msgtypesf.h"
 
370
#include "mafdecls.fh"
 
371
#include "msgids.fh"
 
372
#include "cscfps.fh"
 
373
#include "inp.fh"
 
374
#include "util.fh"
 
375
#include "stdio.fh"
 
376
c
 
377
      character*(*) filename
 
378
      logical found
 
379
c
 
380
      call ga_sync()
 
381
c
 
382
c     Inquire if file is present
 
383
      inquire(file=filename,exist=found)
 
384
      dft_zora_inquire_file = found 
 
385
c
 
386
      call ga_sync()
 
387
c
 
388
      return
 
389
      end
 
390
c
 
391
czora...Read in the zora atomic corrections from disk
 
392
c
 
393
      logical function dft_zora_read(filename, nbf, nsets, nmo, 
 
394
     &               mult, g_zora_sf, g_zora_scale_sf)
 
395
c
 
396
      implicit none
 
397
c
 
398
#include "errquit.fh"
 
399
#include "global.fh"
 
400
#include "tcgmsg.fh"
 
401
#include "msgtypesf.h"
 
402
#include "mafdecls.fh"
 
403
#include "msgids.fh"
 
404
#include "cscfps.fh"
 
405
#include "inp.fh"
 
406
#include "util.fh"
 
407
#include "stdio.fh"
 
408
c
 
409
      integer nsets             ! restricted or unrestricted
 
410
      character*(*) filename
 
411
      integer iset              ! restricted or unrestricted
 
412
      integer g_zora_sf(nsets)
 
413
      integer g_zora_scale_sf(nsets)
 
414
c
 
415
      integer nbf               ! No. of functions in basis
 
416
      integer nmo(nsets)       
 
417
      integer mult              ! multiplicity
 
418
      integer ok, jset, i, j
 
419
      integer l_zora, k_zora
 
420
c
 
421
      integer unitno
 
422
      parameter (unitno = 77)
 
423
      integer inntsize,ddblsize
 
424
c
 
425
      integer nsets_read       ! No. of sets  from file
 
426
      integer nbf_read         ! No. of functions from file 
 
427
      integer mult_read
 
428
c     
 
429
c     Initialise to invalid MA handle
 
430
      l_zora = -1
 
431
c
 
432
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
433
      ddblsize=MA_sizeof(MT_DBL,1,MT_BYTE)
 
434
c
 
435
      call ga_sync()
 
436
      ok = 0
 
437
      if (ga_nodeid() .eq. 0) then
 
438
c
 
439
c      Print a message indicating the file being read
 
440
       write(6,22) filename(1:inp_strlen(filename))
 
441
 22    format(/' Read atomic ZORA corrections from ',a/)
 
442
       call util_flush(luout)
 
443
c
 
444
c      Open the file
 
445
       open(unitno, status='old', form='unformatted', file=filename,
 
446
     $        err=1000)
 
447
c
 
448
c      Read in some basics to check if they are consistent with the calculation
 
449
       read(unitno, err=1001, end=1001) nsets_read
 
450
       read(unitno, err=1001, end=1001) nbf_read
 
451
       read(unitno, err=1001, end=1001) mult_read
 
452
c
 
453
c      Error checks
 
454
       if ((nsets_read .ne. nsets)
 
455
     &  .or. (nbf_read .ne. nbf)
 
456
     &  .or. (mult_read .ne. mult) ) goto 1003
 
457
c
 
458
c      Allocate the temporary buffer
 
459
       if (.not. ma_push_get(mt_dbl,nbf,'dft_zora_read',l_zora,k_zora))
 
460
     $        call errquit('dft_zora_read: ma failed', nbf, MA_ERR)
 
461
c
 
462
c      Read in g_zora_sf
 
463
       do iset = 1, nsets
 
464
        do i = 1, nmo(iset)
 
465
         call sread(unitno, dbl_mb(k_zora), nbf)
 
466
         call ga_put(g_zora_sf(iset), 1, nbf, i, i, dbl_mb(k_zora), 1)
 
467
        end do
 
468
       end do
 
469
c
 
470
c      Read in g_zora_scale_sf
 
471
       do iset = 1, nsets
 
472
        do i = 1, nmo(iset)
 
473
         call sread(unitno, dbl_mb(k_zora), nbf)
 
474
         call ga_put(g_zora_scale_sf(iset), 1, nbf, i, i, 
 
475
     & dbl_mb(k_zora), 1)
 
476
        end do
 
477
       end do
 
478
c
 
479
c      Close the file
 
480
       close(unitno,err=1002)
 
481
       ok = 1
 
482
c
 
483
c      Deallocate the temporary buffer
 
484
       if (.not. ma_pop_stack(l_zora)) call errquit
 
485
     $      ('dft_zora_read: pop failed', l_zora, MA_ERR)
 
486
c
 
487
      end if
 
488
c
 
489
c     Broadcast status to other nodes
 
490
 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
 
491
      call ga_sync()
 
492
c
 
493
c       write(6,*)' g_zora_sf(1) from dft_scf'
 
494
c       call ga_print(g_zora_sf(1))
 
495
c       write(6,*)' g_zora_scale_sf(1) from dft_scf'
 
496
c       call ga_print(g_zora_scale_sf(1))
 
497
c
 
498
      dft_zora_read = ok .eq. 1
 
499
      return
 
500
c
 
501
 1000 write(6,*) 'dft_zora_read: failed to open',
 
502
     $     filename(1:inp_strlen(filename))
 
503
      call util_flush(luout)
 
504
      ok = 0
 
505
      goto 10
 
506
c
 
507
 1001 write(6,*) 'dft_zora_read: failed to read',
 
508
     $     filename(1:inp_strlen(filename))
 
509
      call util_flush(luout)
 
510
      ok = 0
 
511
      close(unitno,err=1002)
 
512
      goto 10
 
513
c
 
514
 1003 write(6,*) 'dft_zora_read: file inconsistent with calculation',
 
515
     $     filename(1:inp_strlen(filename))
 
516
      call util_flush(luout)
 
517
      ok = 0
 
518
      close(unitno,err=1002)
 
519
      goto 10
 
520
c
 
521
 1002 write(6,*) 'dft_zora_read: failed to close',
 
522
     $     filename(1:inp_strlen(filename))
 
523
      call util_flush(luout)
 
524
      ok = 0
 
525
      goto 10
 
526
c
 
527
      end
 
528
c
 
529
czora...Write out the zora atomic corrections to disk
 
530
c
 
531
      logical function dft_zora_write(rtdb, basis, filename,
 
532
     &     nbf, nsets, nmo, mult, g_zora_sf, g_zora_scale_sf)
 
533
c
 
534
      implicit none
 
535
c
 
536
#include "errquit.fh"
 
537
#include "mafdecls.fh"
 
538
#include "global.fh"
 
539
#include "tcgmsg.fh"
 
540
#include "msgtypesf.h"
 
541
#include "inp.fh"
 
542
#include "msgids.fh"
 
543
#include "cscfps.fh"
 
544
#include "util.fh"
 
545
#include "bas.fh"
 
546
#include "geom.fh"
 
547
#include "rtdb.fh"
 
548
#include "stdio.fh"
 
549
c
 
550
c     Temporary routine
 
551
c
 
552
      integer rtdb              ! [input] RTDB handle (-1 if not accessible)
 
553
      integer basis             ! [input] Basis handle(-1 if not accessible)
 
554
      character*(*) filename    ! [input] File to write to
 
555
      integer nbf               ! [input] No. of functions in basis
 
556
      integer nsets             ! [input] No. of sets of matrices
 
557
      integer nmo(nsets)        ! [input] No. of mos in each set
 
558
      integer mult
 
559
c
 
560
      integer g_zora_sf(nsets)    
 
561
      integer g_zora_scale_sf(nsets)  
 
562
c
 
563
      integer unitno
 
564
      parameter (unitno = 77)
 
565
c
 
566
      integer lentit
 
567
      integer lenbas
 
568
      integer l_zora, k_zora
 
569
      integer ok, iset, i, j
 
570
      integer geom, ma_type, nelem
 
571
      character*26 date
 
572
      character*32 geomsum, basissum, key
 
573
      character*20 scftype20 
 
574
      character*128 basis_name, trans_name
 
575
      double precision energy, enrep
 
576
      integer inntsize
 
577
c
 
578
      l_zora = -1               ! An invalid MA handle
 
579
c
 
580
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
581
      call ga_sync()
 
582
c
 
583
      ok = 0
 
584
c
 
585
c     Read routines should be consistent with this
 
586
c
 
587
c     Write out the atomic zora corrections
 
588
c
 
589
      if (ga_nodeid() .eq. 0) then
 
590
c
 
591
c     Open the file
 
592
      open(unitno, status='unknown', form='unformatted',
 
593
     $        file=filename, err=1000)
 
594
c
 
595
c     Write out the number of sets and basis functions
 
596
      write(unitno, err=1001) nsets
 
597
      write(unitno, err=1001) nbf
 
598
      write(unitno, err=1001) mult
 
599
c
 
600
c     Allocate the temporary buffer
 
601
      if (.not. ma_push_get(mt_dbl,nbf,'dft_zora_write',l_zora,k_zora))
 
602
     $        call errquit('dft_zora_write: ma failed', nbf, MA_ERR)
 
603
c
 
604
c     Write out g_zora_sf
 
605
      do iset = 1, nsets
 
606
       do i = 1, nmo(iset)
 
607
        call ga_get(g_zora_sf(iset), 1, nbf, i, i, dbl_mb(k_zora),1)
 
608
        call swrite(unitno, dbl_mb(k_zora), nbf)
 
609
       end do
 
610
      end do
 
611
c
 
612
c     Write out g_zora_scale_sf
 
613
      do iset = 1, nsets
 
614
       do i = 1, nmo(iset)
 
615
        call ga_get(g_zora_scale_sf(iset), 1, nbf, i, i, 
 
616
     &            dbl_mb(k_zora),1)
 
617
        call swrite(unitno, dbl_mb(k_zora), nbf)
 
618
       end do
 
619
      end do
 
620
c
 
621
c     Deallocate the temporary buffer
 
622
      if (.not. ma_pop_stack(l_zora))
 
623
     $  call errquit('dft_zora_write: ma pop failed', l_zora, MA_ERR)
 
624
c
 
625
c     Close the file
 
626
      close(unitno,err=1002)
 
627
c
 
628
      ok = 1
 
629
      end if
 
630
c
 
631
c     Broadcast status to other nodes
 
632
 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
 
633
      call ga_sync()
 
634
c
 
635
c       write(6,*)' g_zora_sf(1) from dft_scf'
 
636
c       call ga_print(g_zora_sf(1))
 
637
c       write(6,*)' g_zora_scale_sf(1) from dft_scf'
 
638
c       call ga_print(g_zora_scale_sf(1))
 
639
c
 
640
      dft_zora_write = (ok .eq. 1)
 
641
      if (ga_nodeid() .eq. 0) then
 
642
         write(6,22) filename(1:inp_strlen(filename))
 
643
 22      format(/' Wrote atomic ZORA corrections to ',a/)
 
644
         call util_flush(luout)
 
645
      endif
 
646
c
 
647
      return
 
648
c
 
649
 1000 write(6,*) 'dft_zora_write: failed to open ',
 
650
     $     filename(1:inp_strlen(filename))
 
651
      call util_flush(luout)
 
652
      ok = 0
 
653
      goto 10
 
654
c
 
655
 1001 write(6,*) 'dft_zora_write: failed to write ',
 
656
     $     filename(1:inp_strlen(filename))
 
657
      call util_flush(luout)
 
658
      ok = 0
 
659
      close(unitno,err=1002)
 
660
      goto 10
 
661
c
 
662
 1002 write(6,*) 'dft_zora_write: failed to close',
 
663
     $     filename(1:inp_strlen(filename))
 
664
      call util_flush(luout)
 
665
      ok = 0
 
666
      goto 10
 
667
c
 
668
      end
 
669
c
 
670
c     Convergence check for zora calculations
 
671
c
 
672
      subroutine dft_zora_scfcvg(rms, derr, etold, etnew, e_conv, 
 
673
     &                      d_conv, g_conv, ipol, iter, iterations, 
 
674
     &                      idone, rtdb, converged, diising)
 
675
c
 
676
c     $Id: dft_zora_utils.F 21851 2012-01-25 00:02:19Z niri $
 
677
c
 
678
      implicit none
 
679
#include "errquit.fh"
 
680
c
 
681
      double precision rms(2)   ! [input]
 
682
      double precision derr(2)  ! [input]
 
683
      double precision etold    ! [input]
 
684
      double precision etnew    ! [input]
 
685
      double precision e_conv   ! [input]
 
686
      double precision d_conv   ! [input]
 
687
      double precision g_conv   ! [input]
 
688
      integer ipol              ! [input]
 
689
      integer iter              ! [input]
 
690
      integer iterations        ! [input]
 
691
      integer idone             ! [output]
 
692
      integer rtdb              ! [input]
 
693
      logical converged         ! [output]
 
694
      logical diising           ! [input]
 
695
c
 
696
#include "mafdecls.fh"
 
697
#include "rtdb.fh"
 
698
c     
 
699
      logical e_conv_logical, d_conv_logical, g_conv_logical
 
700
      logical ENERGY, DENSITY, GRADIENT
 
701
      double precision de, abde
 
702
c
 
703
      converged = .false.
 
704
c
 
705
      e_conv_logical = .false.
 
706
      d_conv_logical = .false.
 
707
      g_conv_logical = .false.
 
708
c
 
709
      ENERGY = e_conv.gt.0
 
710
      DENSITY = d_conv.gt.0
 
711
      GRADIENT = g_conv.gt.0
 
712
c
 
713
      idone = 0
 
714
c
 
715
c     Evaluate change in energy.
 
716
c     
 
717
      de = etnew - etold
 
718
      etold=etnew
 
719
      abde = dabs(de)
 
720
c     
 
721
c     Check to see if energy is converged.
 
722
c
 
723
      if (ENERGY)then     
 
724
         if (abde.lt.e_conv)e_conv_logical = .true.
 
725
      else
 
726
         e_conv_logical = .true.
 
727
      endif
 
728
c     
 
729
c     Check for density matrix convergence.
 
730
c     
 
731
      if (DENSITY)then
 
732
         if (dsqrt(rms(1)).le.d_conv)d_conv_logical = .true.
 
733
         if (ipol.eq.2)then
 
734
           if (dsqrt(rms(2)).le.d_conv) then
 
735
             d_conv_logical = d_conv_logical.and..true.
 
736
          else
 
737
             d_conv_logical = d_conv_logical.and..false.
 
738
          endif
 
739
         endif
 
740
      else
 
741
         d_conv_logical = .true.
 
742
      endif
 
743
c     
 
744
c     Check for gradient convergence.
 
745
c     
 
746
      if (GRADIENT.and.diising)then
 
747
         if (derr(1).le.g_conv)g_conv_logical = .true.
 
748
         if (ipol.eq.2)then
 
749
           if (derr(2).le.g_conv) then
 
750
              g_conv_logical = g_conv_logical.and..true.
 
751
             else
 
752
                g_conv_logical = g_conv_logical.and..false.
 
753
             endif
 
754
         endif
 
755
      else
 
756
         g_conv_logical = .true.
 
757
      endif
 
758
c
 
759
c     Check over-all convergence.
 
760
c
 
761
      converged = e_conv_logical
 
762
      if (converged)idone = 1
 
763
c     
 
764
c     Check iteration value.
 
765
c     
 
766
      if (iter.lt.1)then
 
767
         return
 
768
      elseif (iter.eq.iterations)then
 
769
         idone = 1
 
770
      endif
 
771
c
 
772
c     If all convergence criterion met or number of iterations has been
 
773
c     exceeded, write "converged" to RTDB.
 
774
c
 
775
      if (idone.eq.1)then
 
776
         if (.not.rtdb_put(rtdb, 'dft:converged', MT_LOG, 1, converged))
 
777
     &      call errquit('dft_scfcvg: rtdb_put failed', 1, RTDB_ERR)
 
778
      endif
 
779
      return
 
780
      end