1      subroutine dft_inpana(rtdb)
2c
3c     $Id$
4c
5c*********************************************************************
6c
7c     inpana (input analysis)
8c     Analyze input to deduce nature of system, and set key flags.
9c     Write pertinent information to user output.
10c
11c*********************************************************************
12      implicit none
13#include "errquit.fh"
14c
15#include "stdio.fh"
16#include "rtdb.fh"
17#include "mafdecls.fh"
18#include "global.fh"
19#include "tcgmsg.fh"
20#include "geom.fh"
21#include "bas.fh"
22#include "cdft.fh"
23#include "util.fh"
24#include "util_params.fh"
25c
26      logical even, cksetd, no_prune
27      Integer rtdb              !  runtime database handle
28      integer me,ichg,nel
29      integer noc1
30      integer n, na, nb, test_sic
31      integer nshells
32      integer ictr
33      double precision bsr, radang
34      double precision anucl_charg, ckfac, t1
35      double precision ictr_coord(3), ictr_chg,smear_sigma
36      character*3 on_off_1, on_off_2
37      logical oprint_general, oprint_grid, oprint_xc,
38     &     oprint_convergence, oprint_tolerances, oprint_sic,
39     ,     lsigma
40c
41      logical odft,rodft
42      character*4 scftype
43c
44      logical cam_exch,cam_srhf
45      double precision cam_alpha,cam_beta,cam_omega
46c
47      character*9 local_c, nonlocal_c, nineb_c
48      character*10 start_10c, NA_10c, asap_10c
49      character*10 strng1, strng2, strng3, strng4, strng5, strng6
50      character*16 tag, theory
51      character*20 rgridnames(10)
52      logical xc_gotxc
53      external xc_gotxc
54      data rgridnames /'Euler-MacLaurin','Mura-Knowles',
55     .     'Treutler-Ahlrichs','Gauss-Legendre','G-C-interv',
56     ,     'Lindh','Chebyshev','Legendre','DE2','DE2D'/
57c
58      logical status  ! FA
59      logical use_nwxc
60      logical dftmp2,out1
61      double precision mp2fac
62c
63      me=ga_nodeid()
64c
65      oprint_general = util_print('general information',print_default)
66      oprint_grid = util_print('grid information',print_default)
67      oprint_xc = util_print('xc information',print_default)
68      oprint_convergence = util_print('convergence information',
69     &                                print_default)
70      oprint_tolerances = util_print('screening tolerance information',
71     &                               print_default)
72      oprint_sic = util_print('sic information',print_default)
73c
74c     Figure out the number of electrons from the required total
75c     charge and the sum of nuclear charges
76c
77      if (.not. rtdb_cget(rtdb, 'dft:theory', 1, theory))
78     $        call errquit('dft_inpana: theory not specified',0,
79     &       RTDB_ERR)
80      if (.not. geom_nuc_charge(geom, anucl_charg))
81     &     call errquit('dft_inpana: geom_nuc_charge failed', 0,
82     &       GEOM_ERR)
83      nel = nint(anucl_charg - rcharge)
84      if (nel .le. 0) call errquit
85     $     ('dft_inpana: negative no. of electrons ?', nel, INPUT_ERR)
86      if (abs(anucl_charg - rcharge - dble(nel)) .gt. 1d-8)
87     $     call errquit('dft_inpana: non-integral # of electrons ?', 0,
88     &       INPUT_ERR)
89c
90c     == Coulomb Attenuation Method (CAM/LC) parameters ==
91      if (.not.rtdb_get(rtdb, 'dft:cam_exch', mt_log, 1, cam_exch))
92     &   cam_exch=.false.
93      if (.not.rtdb_get(rtdb, 'dft:cam_srhf', mt_log, 1, cam_srhf))
94     &   cam_srhf=.false.
95      if (.not.rtdb_get(rtdb, 'dft:cam_omega', mt_dbl, 1,cam_omega))
96     &   cam_omega=0.d0
97      if (.not.rtdb_get(rtdb, 'dft:cam_alpha', mt_dbl, 1,cam_alpha))
98     &   cam_alpha=0.d0
99      if (.not.rtdb_get(rtdb, 'dft:cam_beta', mt_dbl, 1, cam_beta))
100     &   cam_beta=0.d0
101c
102c     Check to see if calculation type is allowed.
103c
104c     Even number of electrons required for RHF.
105c
106      even=mod(nel,2).eq.0
107c
108c     == check for restricted open-shell dft ==
109      if (.not. rtdb_get(rtdb, 'dft:rodft', mt_log, 1,
110     &         rodft))rodft = .false.
111c
112      odft = .false.
113      if (.not. rtdb_cget(rtdb, 'dft:scftype', 1, scftype)) then
114c
115c       If the user did not specify whether the calculation should use
116c       an open-shell or closed-shell formalism then work it out from
117c       the spin multiplicity and potentially from a request for rodft.
118c
119        if (mult.eq.1) then
120          scftype = 'RHF'
121        else
122          odft = .true.
123          if (rodft) then
124            scftype = 'ROHF'
125          else
126            scftype = 'UHF'
127          endif
128        endif
129c
130      else
131c
132c       If the user specified what formalism to use then use what the
133c       user ordered. Even if that means running a closed shell system
134c       with the UHF or ROHF formalism.
135c
136        if (scftype.eq.'UHF') then
137          odft = .true.
138        else if (scftype.eq.'ROHF')  then
139c     fix for closed shell and rodft
140           if(mult.eq.1) then
141              odft=.false.
142              scftype='RHF'
143              ipol=1
144           else
145              odft = .true.
146           endif
147
148       elseif (scftype.eq.'RHF'.and.mult.ne.1)  then
149          odft= .true.
150          scftype='UHF'
151        endif
152c
153      endif
154      if (.not. rtdb_cput(rtdb, 'dft:scftype', 1, scftype))
155     &     call errquit('dft_input: rtdb_put failed', 1500, RTDB_ERR)
156      call dft_cscf_scftype(scftype)
157c
158      if ((.not.even).and.(.not.odft).and.(.not.rodft)) then
159       if (ga_nodeid().eq.0)
160     &    write(luout,"(10x,'Number of electrons: ',i10)") nel
161       call errquit(
162     &  'Odd number of electrons. Please specify a
163     &restricted open-shell or open-shell calculation',nel,INPUT_ERR)
164      end if
165c
166c     odd # of electrons or not a singlet state --> LSD
167c
168      if ((.not.even).or.(mult.ne.1).or.(theory.eq.'sodft').or.
169     &    (scftype.eq.'UHF').or.(scftype.eq.'ROHF')) then
170        ipol=2
171      endif
172      noc(2)=0
173c
174c     Calculate number of occupied orbitals.
175c
176      if (ipol.eq.1)then
177         noc1 = nel/2
178         noc(2)= 0
179         noc(1)= noc1
180      else
181c
182c        check consistency of no. elec and multiplicity
183c
184         even=mod((nel+mult-1),2).eq.0
185         if (.not.even) then
186           write(LuOut,*)' number of electrons :',nel
187           write(LuOut,*)' multiplicity        :',mult
188           call errquit(
189     &         ' no. of electrons and multiplicity not compatible',nel,
190     &       INPUT_ERR)
191         endif
192         if(mult.gt.0) then
193            noc(2) = (nel - mult + 1)/2
194            noc(1) = nel - noc(2)
195            noc1  = noc(1) + noc(2)
196         else
197            noc(1) = (nel + mult + 1)/2
198            noc(2) = nel - noc(1)
199            noc1  = noc(1) + noc(2)
200         endif
201      endif
202      if (.not. rtdb_put(rtdb, 'dft:ipol', mt_int, 1, ipol))
203     $     call errquit('inpana: dft:ipol put failed', 0, RTDB_ERR)
204c
205c     Check to see if there are enough electrons for this
206c     value of the multiplicity.
207c
208      if (noc(2).lt.0)then
209         call errquit('dft: #electrons not valid for multiplicity',mult,
210     &       INPUT_ERR)
211      endif
212c
213c     write noc (consistent with definition in ddscf) to rtdb
214c
215      if (.not. rtdb_put(rtdb, 'dft:noc', mt_int, 2, noc))
216     &   call errquit('inpana: rtdb_put of noc failed', 0, RTDB_ERR)
217cHvD
218      use_nwxc = util_module_avail("nwxc")
219      if (use_nwxc) then
220        call nwxc_rtdb_load(rtdb,"dft",use_nwxc)
221      endif
222      if (use_nwxc) then
223         call nwxc_getvals("nwxc_has_hfx",out1)
224        if (out1) then
225           call nwxc_getwght("nwxc_wght_hfx",xfac(1))
226        endif
227         call nwxc_getvals("nwxc_has_mp2c",dftmp2)
228        if (dftmp2) then
229           call nwxc_getwght("nwxc_wght_mp2c",mp2fac)
230          if (.not.rtdb_put(rtdb,'dft:mp2fac', mt_dbl, 1, mp2fac))
231     &      call errquit('dft_inpana: rtdb_put failed', 2902, RTDB_ERR)
232        endif
233        call nwxc_getvals("nwxc_has_cam",cam_exch)
234        if (cam_exch) then
235          call nwxc_get_cam(cam_alpha,cam_beta,cam_omega,cam_srhf)
236       else
237          cam_alpha = 0.0d0
238          cam_beta  = 0.0d0
239          cam_omega = 0.0d0
240        endif
241        call nwxc_print()
242        if (cam_exch.and.(ga_nodeid().eq.0)) then
243            write(LuOut,*)
244            write(LuOut,8202)
245     &         'Range-Separation Parameters        '
246            write(LuOut,8203)
247            write(LuOut,8201)cam_alpha,cam_beta,cam_omega,cam_srhf
248        end if ! cam_exch
249        call xc_setcamparam(rtdb,cam_exch,cam_srhf,cam_omega,cam_alpha,
250     &       cam_beta)
251      endif
252cHvD
253c
254c     Write new data to checkpoint file.
255c
256c     Analyze any user specified XC functionals ... set if none
257c
258      cksetd = .true.
259      cksetd = cksetd .and. .not. libxcon
260c
261c     Check if user has specified some type of functional,
262c     if so, do not set defaults.
263c
264      do n = 1, numfunc
265         if (lcfac(n).or.nlcfac(n).or.lxfac(n).or.nlxfac(n))
266     &        cksetd = .false.
267      enddo
268      if (cksetd)then
269c
270c        Set functional defaults.
271c
272         cfac(1) = 1.0d0
273         lcfac(1) = .true.
274         xfac(2) = 1.0d0
275         lxfac(2) = .true.
276c
277c        Update rtdb.
278c
279         if (.not. rtdb_put(rtdb, 'dft:cfac', mt_dbl, numfunc, cfac))
280     $        call errquit('dft_input: rtdb_put failed', 210, RTDB_ERR)
281         if (.not. rtdb_put(rtdb, 'dft:xfac', mt_dbl, numfunc, xfac))
282     $        call errquit('dft_input: rtdb_put failed', 211, RTDB_ERR)
283         if (.not. rtdb_put(rtdb, 'dft:lcfac', mt_log, numfunc, lcfac ))
284     $        call errquit('dft_input: rtdb_put failed', 9, RTDB_ERR)
285         if (.not. rtdb_put(rtdb, 'dft:lxfac', mt_log, numfunc, lxfac ))
286     $        call errquit('dft_input: rtdb_put failed', 11, RTDB_ERR)
287      endif
288c
289cc AJL/Begin/FDE
290c Repeat for FDE XC. Is this a safe way to check we are using FDE?
291      if (frozemb_fde) then
292c
293c     Analyze any user specified XC functionals ... set if none
294c
295        cksetd = .true.
296c
297c     Check if user has specified some type of functional,
298c     if so, do not set defaults.
299c
300        do n = 1, numfunc
301           if (lcfac_fde(n).or.nlcfac_fde(n).or.
302     &         lxfac_fde(n).or.nlxfac_fde(n))
303     &          cksetd = .false.
304        enddo
305        if (cksetd)then
306c
307c        Set functional defaults.
308c
309           cfac_fde(1) = 1.0d0
310           lcfac_fde(1) = .true.
311           xfac_fde(2) = 1.0d0
312           lxfac_fde(2) = .true.
313c
314c        Update rtdb.
315c
316           if (.not. rtdb_put(rtdb,
317     &             'dft:frozemb:cfac', mt_dbl, numfunc, cfac_fde))
318     $        call errquit('dft_input: rtdb_put failed', 210, RTDB_ERR)
319           if (.not. rtdb_put(rtdb,
320     &             'dft:frozemb:xfac', mt_dbl, numfunc, xfac_fde))
321     $        call errquit('dft_input: rtdb_put failed', 211, RTDB_ERR)
322           if (.not. rtdb_put(rtdb,
323     &             'dft:frozemb:lcfac', mt_log, numfunc, lcfac_fde))
324     $        call errquit('dft_input: rtdb_put failed', 9, RTDB_ERR)
325           if (.not. rtdb_put(rtdb,
326     &             'dft:frozemb:lxfac', mt_log, numfunc, lxfac_fde))
327     $        call errquit('dft_input: rtdb_put failed', 11, RTDB_ERR)
328        endif
329c
330c Repeat something similar for KE functionals.
331c I'm going to set this up so it can be expanded.
332c
333        cksetd = .true.
334c
335c     Check if user has specified some type of functional,
336c     if so, do not set defaults.
337c
338        do n = 1, numfunc_fde
339           if (tsfac(n).gt.1.d-8) cksetd = .false.
340        enddo
341        if (cksetd)then
342c
343c        Set functional defaults.
344c
345          tsfac(1) = 1.d0
346c
347c        Update rtdb.
348c
349          if (.not. rtdb_put(rtdb,
350     &              'dft:frozemb:ts', mt_dbl, numfunc_fde, tsfac))
351     $        call errquit('dft_input: rtdb_put failed', 210, RTDB_ERR)
352        end if ! tsfac default
353c
354      endif ! frozemb_fde
355c AJL/End
356c
357c     Check/set defaults for convergence schemes
358c
359c     Three types of convergence speedup schemes:
360c     1) based on number of cycles performed,
361c     2) based on differences of total energy less than some threshold,
362c     3) DIIS plus levelshifting if homo-lumo gap is small.
363c
364c     Examples of these might be:
365c
366c     1) based on number of cycles performed,
367c     ncydp = 3
368c     ndamp = 40
369c     ncysh = iterations
370c     rlshift = 0.5
371c     ncyds = iterations
372c
373c     2) based on differences of total energy being less than some threshold,
374c     dampon = 1.d8
375c     dampoff = 1.0d-1
376c     levlon = 1.d-1
377c     levloff = 1.0d-3
378c     rlshift = 0.5
379c     diison = 1.d-1
380c     diisoff = 0.0d0
381c
382c     3) DIIS plus levelshifting if homo-lumo gap is small (default)
383c     nodamping = .true.
384c     ncydp = 0
385c     ndamp = 0
386c     ncysh = iterations
387c     rlshift = 0.5
388c     ncyds = iterations
389c
390      if(nodamping)then
391         damp = .false.
392         ncydp = 0
393         ndamp = 0
394      endif
395      if(nolevelshifting)then
396         levelshift = .false.
397         ncysh = 0
398         rlshift = 0.0
399      endif
400      if(nodiis)then
401         diis = .false.
402         ncyds = 0
403      endif
404      if (damp)then
405c
406c        check to make sure either number of damping iterations or
407c        energy criterion has been specified.
408c
409         if (ncydp.eq.0.and.dampon.eq.0.0d0)then
410            ncydp = iterations
411         endif
412      endif
413      if (levelshift)then
414c
415c        check to make sure either number of levelshifting iterations or
416c        energy criterion has been specified.
417c
418         if (ncysh.eq.0.and.levlon.eq.0.0d0)then
419            ncysh = iterations
420         endif
421      endif
422      if (diis)then
423c
424c        check to make sure either number of diis iterations or
425c        energy criterion has been specified.
426c
427         if (ncyds.eq.0.and.diison.eq.0.0d0)then
428            ncyds = iterations
429         endif
430      endif
431c
432c     If convergence input based upon #cycles then turn off energy constraints.
433c
434      if (ncydp.ne.0)then
435         dampon  = -999.9
436         dampoff = -999.9
437      endif
438      if (ncysh.ne.0)then
439         levlon  = -999.9
440         levloff = -999.9
441      endif
442      if (ncyds.ne.0)then
443         diison  = -999.9
444         diisoff = -999.9
445      endif
446c
447c     check special case with damping - change
448c     default of 2 to "iterations" if no other
449c     convergence control specified
450c
451      if (ncysh.eq.0 .and. ncyds.eq.0 .and. ncydp.eq.2)
452     &   ncydp = iterations
453c
454c     check on 2-e integral and XC grid tolerances
455c
456      call dft_inpanae(rtdb)
457c
458c     Check for no pruning; no_prune
459c
460      if (.not. rtdb_get(rtdb, 'dft:no_prune', mt_log, 1,
461     &         no_prune))no_prune = .false.
462      if (.not. rtdb_get(rtdb, 'dft:test_sic', mt_int, 1,
463     &     test_sic))test_sic = 0
464      lsigma=rtdb_get(rtdb, 'dft:smear_sigma',mt_dbl,1,smear_sigma)
465c
466      if (me.eq.0) then
467c
468c        Write to output.
469c
470         if (oprint_general)then
471            write(LuOut,*)
472            call util_print_centered
473     &         (LuOut,'General Information',20,.true.)
474            write(LuOut,9020)
475            if (scftype.eq."RHF")then
476               write(LuOut,9050)
477            elseif (scftype.eq."ROHF")then
478               write(LuOut,9052)
479            elseif (scftype.eq."UHF")then
480               write(LuOut,9055)
481            endif
482            write(LuOut,8150)ncenters
483            na = noc(1)
484            nb = noc(2)
485            if (ipol.eq.1)nb = noc(1)
486            ichg = nint(rcharge)
487            write(LuOut,8200)nel,na,nb,ichg,mult
488c
489cFA         Writing Number of electrons on rtdb to be used
490c           in create_munu4nbo() defined in hnd_efgmap_Z4.F
491            status = rtdb_parallel(.false.)
492            if (.not. rtdb_put(rtdb, 'prop:Nel',mt_int,
493     $                      1,nel))
494     $      call errquit('prop_input-EFGZ4-nel: rtdb_put failed',
495     $                   555, RTDB_ERR)
496            if (.not. rtdb_put(rtdb, 'prop:Nocc_a',mt_int,
497     $                      1,na))
498     $      call errquit('prop_input-EFGZ4-Nocc_a: rtdb_put failed',
499     $                   555, RTDB_ERR)
500            if (.not. rtdb_put(rtdb, 'prop:Nocc_b',mt_int,
501     $                      1,nb))
502     $      call errquit('prop_input-EFGZ4-Nocc_b: rtdb_put failed',
503     $                   555, RTDB_ERR)
504            status = rtdb_parallel(.true.)
505cFA
506            on_off_1 = 'off'
507            if (oskel)on_off_1 = 'on'
508            on_off_2 = 'off'
509            if (oadapt)on_off_2 = 'on'
510            write(LuOut,9056)on_off_1, on_off_2
511            write(LuOut,9030)iterations
512            if (direct)write(LuOut,9110)
513            if (.not. bas_numcont(AO_bas_han, nshells))
514     &         call errquit('rdinput:rdinput:',86, BASIS_ERR)
515            write(LuOut,4001)nbf, nshells
516            if (nbf_cd.gt.0)then
517               write(LuOut,9130)
518              if (.not. bas_numcont(CD_bas_han, nshells_cd))
519     &           call errquit('rdinput:rdinput:',87, BASIS_ERR)
520               write(LuOut,4002)nbf_cd, nshells_cd
521            endif
522            if (nbf_xc.gt.0)then
523               write(LuOut,9120)
524               if (.not. bas_numcont(XC_bas_han, nshells_xc))
525     &            call errquit('rdinput:rdinput:',88, BASIS_ERR)
526               write(LuOut,4003)nbf_xc, nshells_xc
527            endif
528            write(LuOut,9035)e_conv
529            if (d_conv.gt.0)then
530               write(LuOut,9040)d_conv
531            endif
532            if (g_conv.gt.0)then
533               write(LuOut,9045)g_conv
534            endif
535            call util_flush(LuOut)
536         endif
537         if (oprint_xc)then
538            write(LuOut,*)
539            call util_print_centered
540     &         (LuOut,'XC Information',20,.true.)
541            if (libxcon) call nwchem_libxc_print_header()
542c
543c           Write out XC info. Combo info first, than X components,
544c           than C components.
545c
546            local_c = 'local    '
547            nonlocal_c = 'non-local'
548            nineb_c = '         '
549            do n = 1, numfunc
550               if (xccomb(n))write(LuOut,9223) xcname(n)
551            enddo
552c
553c           Do exact exchange differently.
554c
555            if (lxfac(1).or.nlxfac(1))
556     &          write(LuOut,9224) xname(1), xfac(1), nineb_c
557c
558            do n = 2, numfunc
559               if (lxfac(n).and.nlxfac(n))then
560                  write(LuOut,9224) xname(n), xfac(n), nineb_c
561               elseif (lxfac(n).and.(.not.nlxfac(n)))then
562                  write(LuOut,9224) xname(n), xfac(n), local_c
563               elseif ((.not.lxfac(n)).and.nlxfac(n))then
564                  write(LuOut,9224) xname(n), xfac(n), nonlocal_c
565               endif
566            enddo
567            do n = 1, numfunc
568               if (lcfac(n).and.nlcfac(n))then
569                  write(LuOut,9224) cname(n), cfac(n), nineb_c
570               elseif (lcfac(n).and.(.not.nlcfac(n)))then
571                  write(LuOut,9224) cname(n), cfac(n), local_c
572               elseif ((.not.lcfac(n)).and.nlcfac(n))then
573                  write(LuOut,9224) cname(n), cfac(n), nonlocal_c
574               endif
575            enddo
576
577            if (libxcon) call nwchem_libxc_print()
578c
579c           Check XC coefficients to make sure appropriate components
580c           sum to 1.0
581c
582c            ckfac = 0.0d0
583c            do n = 1, numfunc
584c               if (lcfac(n))ckfac = ckfac + cfac(n)
585c            enddo
586c            if (abs(ckfac-1.0d0).gt.1.d-8)then
587c               write(LuOut,*)
588c     &            ' WARNING: Sum of local correlation is ',ckfac
589c               write(LuOut,*)' Sum of components do not equal unity. '
590c            endif
591c            ckfac = 0.0d0
592c            do n = 1, numfunc
593c               if (nlcfac(n))ckfac = ckfac + cfac(n)
594c            enddo
595c            if (abs(ckfac-1.0d0).gt.1.d-8.and.abs(ckfac).gt.1.d-8)then
596c               write(LuOut,*)
597c     &            ' WARNING: Sum of nonlocal correlation is ',ckfac
598c               write(LuOut,*)
599c     &            ' Sum of components do not equal unity or 0. '
600c            endif
601c            ckfac = 0.0d0
602c            do n = 1, numfunc
603c               if (lxfac(n))ckfac = ckfac + xfac(n)
604c            enddo
605c            if (abs(ckfac-1.0d0).gt.1.d-8)then
606c               write(LuOut,*)
607c     &            ' WARNING: Sum of local exchange is ',ckfac
608c               write(LuOut,*)' Sum of components do not equal unity. '
609c            endif
610c            ckfac = 0.0d0
611c            do n = 1, numfunc
612c               if (nlxfac(n))ckfac = ckfac + xfac(n)
613c            enddo
614c            if (abs(ckfac-1.0d0).gt.1.d-8.and.abs(ckfac).gt.1.d-8)then
615c               write(LuOut,*)
616c     &            ' WARNING: Sum of nonlocal exchange is ',ckfac
617c               write(LuOut,*)
618c     &            ' Sum of components do not equal unity or 0. '
619c            endif
620c
621c           Check if asymptotic correction will be added to potential
622c           If both LB94 and CS00 are .true., it is assumed that the
623c           user meant to use CS00 (since CS00 uses LB94)
624c
625            if (cs00) then
626               if (delta_ac.gt.1.0d90) then
627               write(LuOut,*)
628               write(LuOut,9226)
629     &         '  CS with a Zhan-Nichols-Dixon shift        '
630               else
631               write(LuOut,*)
632               write(LuOut,9227)
633     &         '  Casida-Salahub correction with a shift    ',
634     &         delta_ac,'au'
635               endif
636            else if (lb94) then
637               write(LuOut,*)
638               write(LuOut,9226)
639     &         '         van Leeuwen-Baerends correction    '
640            else if (ncap) then
641               write(LuOut,*)
642               write(LuOut,9227)
643     &         '  NCAP with derivative discontinuity shift    ',
644     &         delta_ac,'au'
645            endif
646c
647c           Print range-separation parameters
648            if (cam_exch) then
649              write(LuOut,*)
650              write(LuOut,8202)
651     &         'Range-Separation Parameters        '
652              write(LuOut,8203)
653              write(LuOut,8201)cam_alpha,cam_beta,cam_omega,cam_srhf
654            end if ! cam_exch
655c
656cc AJL/Begin/FDE
657            if (frozemb_fde) then
658              write(LuOut,*)
659              call util_print_centered
660     &           (LuOut,'XC Information (FDE)',23,.true.)
661c
662c           Write out XC info. Combo info first, than X components,
663c           than C components.
664c
665              do n = 1, numfunc
666                 if (xccomb_fde(n))write(LuOut,9223) xcname(n)
667              enddo
668c
669c           Do exact exchange differently.
670c AJL: Hybrids will give you linear dependencies in the fde/qm combined
671c wavefunction. Therefore, we should never have non-local contributions
672c from exact exchange. I have left the functionality here, as it is
673c lifted from L510 above, but it should never be used.
674c
675c An error message is implemented in dft_rdinput.F
676c If someone thinks it is sensible to re-enable this they can do so themselves.
677c
678              if (lxfac_fde(1).or.nlxfac_fde(1))
679     &            write(LuOut,9224) xname(1), xfac_fde(1), nineb_c
680c
681              do n = 2, numfunc
682                 if (lxfac_fde(n).and.nlxfac_fde(n))then
683                    write(LuOut,9224) xname(n), xfac_fde(n), nineb_c
684                 elseif (lxfac_fde(n).and.(.not.nlxfac_fde(n)))then
685                    write(LuOut,9224) xname(n), xfac_fde(n), local_c
686                 elseif ((.not.lxfac_fde(n)).and.nlxfac_fde(n))then
687                    write(LuOut,9224) xname(n), xfac_fde(n), nonlocal_c
688                 endif
689              enddo
690              do n = 1, numfunc
691                 if (lcfac_fde(n).and.nlcfac_fde(n))then
692                    write(LuOut,9224) cname(n), cfac_fde(n), nineb_c
693                 elseif (lcfac_fde(n).and.(.not.nlcfac_fde(n)))then
694                    write(LuOut,9224) cname(n), cfac_fde(n), local_c
695                 elseif ((.not.lcfac_fde(n)).and.nlcfac_fde(n))then
696                    write(LuOut,9224) cname(n), cfac_fde(n), nonlocal_c
697                 endif
698              enddo
699c
700c The most sensible thing is again to do the error checks for
701c range-separated and self-interaction corrected setups. Due to the
702c complexity of implementation I'm just going to disable this for FDE, but
703c they could in theory be checked and implemented.
704c The error check is in frozemb_fde_init, which is called from dft_scf.
705c
706ccc DISABLED ccc
707c
708cc Lifted from above:
709cc           Check if asymptotic correction will be added to potential
710cc           If both LB94 and CS00 are .true., it is assumed that the
711cc           user meant to use CS00 (since CS00 uses LB94)
712cc
713c              if (cs00) then
714c                 if (delta_ac.gt.1.0d90) then
715c                 write(LuOut,*)
716c                 write(LuOut,9226)
717c     &           '  CS with a Zhan-Nichols-Dixon shift        '
718c                 else
719c                 write(LuOut,*)
720c                 write(LuOut,9227)
721c     &           '  Casida-Salahub correction with a shift    ',
722c     &           delta_ac,'au'
723c                 endif
724c              else if (lb94) then
725c                 write(LuOut,*)
726c                 write(LuOut,9226)
727c     &           '         van Leeuwen-Baerends correction    '
728c              endif
729c
730c           Print range-separation parameters
731c              if (cam_exch) then
732c                write(LuOut,*)
733c                 write(LuOut,8202)
734c     &           'Range-Separation Parameters        '
735c                 write(LuOut,8203)
736c                write(LuOut,8201)cam_alpha,cam_beta,cam_omega,cam_srhf
737c              end if ! cam_exch
738ccc END DISABLED ccc
739c
740c Kinetic Energy Functional
741c
742              write(LuOut,*)
743              call util_print_centered
744     &           (LuOut,'KE Information (FDE)',23,.true.)
745              do n = 1, numfunc_fde
746                 if (tsfac(n).gt.0)
747     &             write(LuOut,9224) tsname(n), tsfac(n), local_c
748              enddo
749            endif ! frozemb_fde
750cc AJL/End
751c
752            call util_flush(LuOut)
753c
754         endif ! oprint_xc
755c
756         if (oprint_sic) then
757           if (test_sic.eq.1) then
758             write(LuOut,'(/,14x,"SIC perturbative approximation")')
759           else
760             if (test_sic.eq.2) then
761               write(LuOut,
762     .               '(/14x,"SIC/OEP without localized orbitals")')
763             else
764               if (test_sic.eq.4) then
765                 write(LuOut,
766     .                 '(/14x,"SIC/OEP with localized orbitals")')
767               end if
768             end if
769           end if
770         end if
771         if (oprint_grid.and.xc_gotxc())then
772            write(LuOut,*)
773            call util_print_centered
774     &           (LuOut,'Grid Information',20,.true.)
775            write(LuOut,9135)gridtype
776            write(LuOut,9136) rgridnames(wradgrid)
777            if (.not.leb)then
778               write(LuOut,9142)
779            else
780               write(LuOut,9143)
781            endif
782            write(LuOut,9144)
783            do n = 1, ntypes
784c
785c              Find an atom of this kind in the complete list.
786c
787               do ictr = 1, ncenters
788                  if (iatype(ictr).eq.n) then
789                     if (.not. geom_cent_get(geom, ictr, tag,
790     &                  ictr_coord, ictr_chg))call errquit
791     &                  ('dft_inpana: geom_cent_get failed', 0,
792     &       GEOM_ERR)
793                     goto 40
794                  endif
795               enddo
796   40          continue
797               bsr = bsrad_atom_type(n)*CAU2ANG
798               radang = dble(nint(cau2ang*dble(rad_cutoff(1,n))))
799               if (leb)then
800                  write(LuOut,9138)tag,bsr,nrad(n),radang,
801     &                nang(n)
802               else
803                  write(LuOut,9137)tag,bsr,nrad(n),radang,
804     &                nang(n),2*nang(n)
805               endif
806            enddo
807            if (no_prune)then
808               on_off_1 = 'off'
809            else
810               on_off_1 = 'on'
811            endif
812            write(LuOut,4005)on_off_1
813            write(LuOut,4004)nqshells
814c
815            if (ldelley)then
816               write(LuOut,9140) 'Delley'
817            elseif(lssw) then
818               if(whichssw.eq.'ssf ') then
819                  write(LuOut,9140) ' Stratmann-Scuseria-Frisch'
820               elseif(whichssw.eq.'erf1') then
821                  write(LuOut,9140) ' Erf1'
822               elseif(whichssw.eq.'erf2') then
823                  write(LuOut,9140) ' Erf2'
824               else
825                  write(LuOut,9140) whichssw
826               endif
827            else
828               write(LuOut,9140) 'Becke'
829            endif
830            if (nquad_task.ne.1)then
831               write(LuOut,9141)nquad_task
832            endif
833            if (nq_chunk.ne.0)then
834               write(LuOut,9145)nq_chunk
835            endif
836            call util_flush(LuOut)
837         endif
838         if (oprint_convergence)then
839            write(LuOut,*)
840            call util_print_centered
841     &         (LuOut,'Convergence Information',20,.true.)
842            write(LuOut,3231)hl_tol, nfock
843            write(LuOut,3232)ndamp, rlshift
844            asap_10c  = '  ASAP    '
845            start_10c = '  start   '
846            NA_10c = '   N/A    '
847            if(ncydp.ne.0)then
848               strng1 = start_10c
849               write(strng4,'(i3,7h iters )')ncydp
850            elseif(nodamping)then
851               strng1 = NA_10c
852               strng4 = NA_10c
853            else
854               write(strng1,'(d10.2)')dampon
855               write(strng4,'(d10.2)')dampoff
856            endif
857c
858            if(ncysh.ne.0)then
859               strng2 = asap_10c
860               write(strng5,'(i3,7h iters )')ncysh
861            elseif(nolevelshifting)then
862               strng2 = NA_10c
863               strng5 = NA_10c
864            else
865               write(strng2,'(d10.2)')levlon
866               write(strng5,'(d10.2)')levloff
867            endif
868c
869            if(ncyds.ne.0)then
870               strng3 = start_10c
871               write(strng6,'(i3,7h iters )')ncyds
872            elseif(nodiis)then
873               strng3 = NA_10c
874               strng6 = NA_10c
875            else
876               write(strng3,'(d10.2)')diison
877               write(strng6,'(d10.2)')diisoff
878            endif
879            write(LuOut,3233)strng1,strng2,strng3,strng4,strng5,strng6
880            call util_flush(LuOut)
881         endif
882         if(lsigma .and. oprint_general) then
883            write(luout,
884     .           "(10x,'Smearing applied: ',d9.2,' (hartree)')"
885     .           ) smear_sigma
886         endif
887
888
889         if (oprint_tolerances)then
890            write(LuOut,*)
891            call util_print_centered
892     &         (LuOut,'Screening Tolerance Information',20,.true.)
893            t1 = 10.d0**(-itol2e)
894c            write(LuOut,9372)tol_rho, iaoacc, icdacc, ixcacc, t1,
895c     &                    r1
896            write(LuOut,9372)tol_rho, iaoacc, icdacc, ixcacc, t1
897            call util_flush(LuOut)
898         endif
899      endif
900c
901      return
902 3231 format(10x,'Convergence aids based upon iterative change in ',/,
903     &       10x,'total energy or number of iterations. ',/,
904     &       10x,'Levelshifting, if invoked, occurs when the ',/,
905     &       10x,'HOMO/LUMO gap drops below (HL_TOL): ',1Pd9.2,/,
906     &       10x,'DIIS, if invoked, will attempt to extrapolate ',/,
907     &       10x,'using up to (NFOCK): ',i2,' stored Fock matrices.',/)
908 3232 format(10x,
909     &       10x,'Damping(',i2,'%)  Levelshifting(',f3.1,
910     &           ')       DIIS',/,
911     &       10x,8x,15('-'),1x,19('-'),1x,15('-'))
912 3233 format(10x,'dE  on:',2x,a10,7x,a10,10x,a10,/,
913     &       10x,'dE off:',2x,a10,7x,a10,10x,a10,/)
914 4001 format(10x,'AO basis - number of functions: ',i5,/,
915     &       10x,'           number of shells: ',i5)
916 4002 format(10x,'CD basis - number of functions: ',i5,/,
917     &       10x,'           number of shells: ',i5)
918 4003 format(10x,'XC basis - number of functions: ',i5,/,
919     &       10x,'           number of shells: ',i5)
920 4004 format(10x,'Number of quadrature shells: ',i5)
921 4005 format(10x,'Grid pruning is: ',a3)
922 9020 format(10x,'SCF calculation type: DFT')
923 9030 format(10x,'Maximum number of iterations: ',I3)
924 9035 format(10x,'Convergence on energy requested: ',1Pd9.2)
925 9040 format(10x,'Convergence on density requested: ',1Pd9.2)
926 9045 format(10x,'Convergence on gradient requested: ',1Pd9.2)
927 9050 format(10x,'Wavefunction type:  closed shell.')
928 9052 format(10x,'Wavefunction type:  restricted open shell.')
929 9055 format(10x,'Wavefunction type:  spin polarized.')
930 9056 format(10x,'Use of symmetry is: ',a3,
931     &           '; symmetry adaption is: ',a3)
932 9110 format(10x,'This is a Direct SCF calculation.')
933 9120 format(10x,'An Exch-Corr fitting basis will be used.')
934 9130 format(10x,'A Charge density fitting basis will be used.')
935 9135 format(10x,'Grid used for XC integration:  ',a)
936cedo 9136 format(10x,'Radial quadrature: Euler-MacLaurin. ')
937 9136 format(10x,'Radial quadrature: ',A)
938 9137 format(10x,a16,2x,f6.2,6x,i3,6x,f6.1,3x,i2,1x,'*',1x,i2)
939 9138 format(10x,a16,2x,f6.2,6x,i3,8x,f6.1,3x,2x,i5)
940c 9139 format(10x,'Spatial weights used: Delley. ')
941 9140 format(10x,'Spatial weights used: ',A)
942 9141 format(10x,'Parallel task size associated with evaluation of ',/,
943     &       10x,'grid based components has been modified to: ',i2)
944 9145 format(10x,'Chunking of the angular grid is being used; ',
945     &           'nq_chunk = ',i4)
946 9142 format(10x,'Angular quadrature: Gauss-Legendre. ')
947 9143 format(10x,'Angular quadrature: Lebedev. ')
948 9144 format(10x,'Tag',14x,'B.-S. Rad.',1x,'Rad. Pts.',1x,'Rad. Cut.',
949     &        1x,'Ang. Pts.',/,
950     &       10x,'---',14x,'----------',1x,'---------',1x,'---------',
951     &        1x,'---------')
952 9223 format(10x,a40)
953 9224 format(10x,a40,1x,f6.3,1x,a9)
954 9225 format(10x,a40,1x,f8.4)
955 9226 format(10x,a44)
956 9227 format(10x,a44,1x,f10.6,1x,a2)
957 8150 format(10x,'No. of atoms     :',2x,i4)
958 8200 format(10x,'No. of electrons :',2x,i4,/,
959     &       10x,' Alpha electrons :',2x,i4,/,
960     &       10x,'  Beta electrons :',2x,i4,/,
961     &       10x,'Charge           :',2x,i4,/,
962     &       10x,'Spin multiplicity:',2x,i4)
963 8201 format(10x,' Alpha           :',f6.2,/,
964     &       10x,' Beta            :',f6.2,/,
965     &       10x,' Gamma           :',f6.2,/,
966     &       10x,' Short-Range HF  :',l6)
967 8202 format(2x,a44)
968 8203 format(11x,'---------------------------')
969 9372 format(10x,'Density screening/tol_rho: ',1Pd9.2,/,
970     &       10x,'AO Gaussian exp screening on grid/accAOfunc: ',i3,/,
971     &       10x,'CD Gaussian exp screening on grid/accCDfunc: ',i3,/,
972     &       10x,'XC Gaussian exp screening on grid/accXCfunc: ',i3,/,
973     &       10x,'Schwarz screening/accCoul: ',1Pd9.2,/)
974      end
975      subroutine dft_inpanae(rtdb)
976      implicit none
977#include "errquit.fh"
978#include "cdft.fh"
979#include "mafdecls.fh"
980#include "global.fh"
981#include "stdio.fh"
982#include "rtdb.fh"
983#include "util.fh"
984      integer rtdb
985c
986      logical oprint_general,oprint_tolgr
987      double precision tol2e
988c
989      oprint_general = util_print('general information',print_default)
990      oprint_tolgr = util_print('grid_tol_info', print_high)
991c
992c     get e_conv from rtdb because hess_init
993c
994      if (.not.rtdb_get(rtdb,'dft:e_conv',mt_dbl,1,e_conv))
995     .     call errquit('dftinpanae: rtdbget econv failed',0, RTDB_ERR)
996c
997c     check integral tolerances to make sure they match
998c     requested convergence tolerances.
999c
1000c     make sure itol2e is less than 0.1*e_conv or (0.01*g_conv**2)
1001c
1002      if (.not. rtdb_get(rtdb, 'dft:itol2e',
1003     &      mt_int, 1, itol2e)) itol2e=8
1004      if (10.d0**(-itol2e).gt.(1.0d-1*e_conv))then
1005         itol2e = -nint(log10(1.0d-1*e_conv))
1006         if (.not. rtdb_put(rtdb, 'dft:itol2e',
1007     &      mt_int, 1, itol2e))
1008     &      call errquit('dft_inpanae: rtdb_put failed', 127, RTDB_ERR)
1009         if (ga_nodeid().eq.0.and.oprint_general)then
1010            write(LuOut,*)' itol2e modified to match energy'
1011            write(LuOut,*)' convergence criterion.'
1012         endif
1013      endif
1014c
1015c     check density tolerance to make sure it matches
1016c     requested convergence tolerances.
1017c
1018c     make sure tol_rho is less than 0.01*e_conv or (0.01*g_conv**2)
1019c
1020      if (tol_rho.gt.(1.0d-3*e_conv))then
1021         tol_rho = 1.0d-3*e_conv
1022         if (.not. rtdb_put(rtdb, 'dft:tol_rho',
1023     &      mt_dbl, 1, tol_rho))
1024     &      call errquit('dft_inpanae: rtdb_put failed', 127, RTDB_ERR)
1025         if (ga_nodeid().eq.0.and.oprint_general)then
1026            write(LuOut,*)' tol_rho modified to match energy'
1027            write(LuOut,*)' convergence criterion.'
1028         endif
1029      endif
1030      if (.not. rtdb_get(rtdb, 'dft:iAOacc', mt_int, 1,
1031     &   iAOacc))then
1032         iAOacc=-nint(log(e_conv))
1033      else
1034        iAOacc=max(iAOacc,-nint(log(e_conv)))
1035      endif
1036       if (.not. rtdb_put(rtdb, 'dft:iAOacc',
1037     &      mt_int, 1, iAOacc))
1038     &      call errquit('dft_inpanae: rtdb_put failed', 124, RTDB_ERR)
1039      tol2e = 10.d0**(-itol2e)
1040      if (.not. rtdb_put(rtdb, 'dft:tol2e', MT_DBL, 1, tol2e))
1041     .     call errquit('dftinpanae:rtdbput failed',0, RTDB_ERR)
1042      if (.not. rtdb_put(rtdb, 'scf:tol2e', MT_DBL, 1, tol2e))
1043     .     call errquit('dftinpanae:rtdbput failed',0, RTDB_ERR)
1044      if(oprint_tolgr.and.ga_nodeid().eq.0) then
1045        write(luout,*) ' dftinpanae: itol2e ',itol2e,
1046     ,     ' iaoacc ',iaoacc,
1047     ,    ' tol_rho ',tol_rho
1048       call util_flush(luout)
1049      endif
1050      return
1051      end
1052