1      block data cscfdata
2      implicit none
3#include "cscf.fh"
4c
5      data oinitialized /.false./
6      data g_movecs /0,0/
7      data k_occ /10000000/
8      data k_eval /10000000/
9c
10      end
11      subroutine scf_init(rtdb)
12      implicit none
13#include "errquit.fh"
14#include "mafdecls.fh"
15#include "global.fh"
16#include "bas.fh"
17#include "geom.fh"
18#include "rtdb.fh"
19#include "inp.fh"
20#include "sym.fh"
21#include "util.fh"
22#include "cscf.fh"
23#include "cosmo.fh"
24c
25      integer rtdb              ! database handle
26c
27      double precision nuclear_charge
28      character*255 name
29      integer len_occ
30      external cscfdata ! For T3D linker
31c
32      logical osome
33c
34      logical hf_job
35      character*30 tag
36      character*255 theory
37      integer mult
38      logical xc_gotxc,xc_gothfx
39      external xc_gotxc,xc_gothfx
40c
41      hf_job = (.not. xc_gotxc()).and.(.not. xc_gothfx())
42c
43      if (.not. rtdb_cget(rtdb, 'title', 1, title))
44     $     title = ' '
45c
46c     load geometry and symmetry info
47c
48      if (.not. geom_create(geom, 'geometry'))
49     $     call errquit('scf_init: geom_create?', 0, GEOM_ERR)
50      if (.not. geom_rtdb_load(rtdb, geom, 'geometry'))
51     $     call errquit('scf_init: no geometry ', 0, RTDB_ERR)
52      if (.not.rtdb_get(rtdb, 'scf:skeleton',MT_LOG, 1, oskel)) then
53         oskel = sym_number_ops(geom) .gt. 0
54      endif
55      if (.not.rtdb_get(rtdb, 'scf:adapt',MT_LOG, 1, oadapt)) then
56         oadapt = sym_number_ops(geom) .gt. 0
57      endif
58      if (.not.rtdb_get(rtdb, 'scf:lock',MT_LOG, 1, olock)) then
59         olock = .false.
60      endif
61c
62c     load the basis set and get info about it
63c
64      if (.not. bas_create(basis, 'ao basis'))
65     $     call errquit('scf_init: bas_create?', 0, BASIS_ERR)
66      if (.not. bas_rtdb_load(rtdb, geom, basis, 'ao basis'))
67     $        call errquit('scf_init: no ao basis set', 0, RTDB_ERR)
68c
69c     For debug ... call int_init and do the 2-e
70c
71      if (util_print('texas init debug',print_never)) then
72        call int_init(rtdb, 1, basis)
73        write(6,*) ' DONE INIT'
74*        call schwarz_init(geom, basis)
75*        call schwarz_tidy()
76        call int_terminate
77      endif
78c
79      if (.not. bas_name(basis, name, trans))
80     $     call errquit('scf_init: bas_name?', 0, BASIS_ERR)
81c
82      if (.not. bas_numbf(basis, nbf)) call errquit
83     $     ('scf_init: basis info',0, BASIS_ERR)
84c
85c     Is RI approximation to be used? If so get fitting basis set.
86c
87      if (rtdb_get(rtdb, 'scf:ri', MT_INT, 1, nriscf)) then
88        if (.not. bas_create(riscf_basis, 'riscf basis'))
89     $       call errquit('scf_init: bas_create?', 0, BASIS_ERR)
90        if (.not. bas_rtdb_load(rtdb, geom, riscf_basis, 'riscf basis'))
91     $       call errquit('scf_init: no riscf basis set', 0, RTDB_ERR)
92      else
93        nriscf = 0
94        riscf_basis = 0
95      endif
96c
97c     Figure input/output MO vectors
98c
99      if (hf_job) then
100         tag = 'scf:input vectors'
101      else
102         tag = 'dft:input vectors'
103      endif
104      if (.not. rtdb_cget(rtdb, tag, 1, movecs_in))
105     $     movecs_in = 'atomic'
106      if (hf_job) then
107         tag = 'scf:output vectors'
108      else
109         tag = 'dft:output vectors'
110      endif
111      if (.not. rtdb_cget(rtdb, tag, 1, movecs_out))
112     $     movecs_out = ' '
113      if (movecs_out.eq.' ') then
114         if (movecs_in.eq.'atomic' .or. movecs_in.eq.'hcore' .or.
115     $        movecs_in.eq.'project' .or. movecs_in.eq.'fragment'.or.
116     $        movecs_in.eq.'molden'
117     $          .or.movecs_in.eq.'rotate') then
118            call util_file_name('movecs', .false.,.false.,movecs_out)
119         else
120            movecs_out = movecs_in
121         endif
122      endif
123c
124c     Resolve names of MO files to full paths defaulting to the
125c     permanent directory
126c
127      if (movecs_in.eq.'atomic' .or. movecs_in.eq.'hcore' .or.
128     $        movecs_in.eq.'project' .or. movecs_in.eq.'fragment'.or.
129     $        movecs_in.eq.'molden'
130     $          .or.movecs_in.eq.'rotate') then
131         continue
132      else
133         call util_file_name_resolve(movecs_in, .false.)
134      endif
135      call util_file_name_resolve(movecs_out, .false.)
136c
137c     Figure out the number of electrons from the required total
138c     charge and the sum of nuclear charges
139c
140      if (.not. rtdb_get(rtdb, 'charge', MT_DBL, 1, charge))
141     $     charge = 0.0d0
142      if (.not. geom_nuc_charge(geom, nuclear_charge))
143     $     call errquit('scf: geom_nuc_charge failed', 0, GEOM_ERR)
144      nelec = nint(nuclear_charge - charge)
145      if (nelec .le. 0) call errquit
146     $     ('scf: negative no. of electrons ?', nelec, INPUT_ERR)
147      if (abs(nuclear_charge - charge - dble(nelec)) .gt. 1d-8)
148     $     call errquit('scf: non-integral no. of electrons ?', 0,
149     &       INPUT_ERR)
150c
151c     Determine no. of open and closed shells ... default is to run closed
152c     shell unless told otherwise
153c
154      if(.not.rtdb_cget(rtdb,'task:theory',1,theory))
155     +     call errquit('task: no task input for theory?',0, INPUT_ERR)
156c     if (.not. rtdb_get(rtdb, 'scf:nopen', MT_INT, 1, nopen))
157c     $     nopen = 0
158      if(theory .eq. 'dft')then
159         if (.not. rtdb_get(rtdb, 'dft:mult', MT_INT, 1,mult))
160     *        mult = 1
161         nopen = mult - 1
162      endif
163      if (.not. rtdb_get(rtdb, 'scf:nopen', MT_INT, 1, nopen))
164     $     nopen = 0
165      if (nopen .gt. nelec) call errquit
166     $     ('scf: nopen > nelec ', nopen, INPUT_ERR)
167      if (mod(nelec-nopen,2) .ne. 0) call errquit
168     $     ('scf: no. of closed-shell electrons is not even!',nopen,
169     &       INPUT_ERR)
170      nclosed = (nelec-nopen) / 2
171c
172      if (.not. rtdb_cget(rtdb, 'scf:scftype', 1, scftype)) then
173         if (nopen .eq. 0) then
174            scftype = 'RHF'
175         else
176            scftype = 'ROHF'
177         endif
178      endif
179c
180      call inp_ucase(scftype)
181c
182c     Take care of holes in the input routines
183c
184      if (scftype.eq.'RHF' .and. nopen.gt.0) then
185         scftype = 'ROHF'
186      else if (scftype.eq.'ROHF' .and. nopen.eq.0) then
187         scftype = 'RHF'
188      endif
189c
190      if ( scftype.ne.'ROHF' .and. scftype.ne.'RHF' .and.
191     $     scftype.ne.'UHF' ) call errquit
192     $     ('scf: only ROHF, RHF, and UHF currently supported', 0,
193     &       INPUT_ERR)
194c
195c     Dump lagrangian?  Yes by default now since if the SCF
196c     has converged for an energy it will not be rerun for the gradient.
197c
198      if (.not.rtdb_get(rtdb, 'scf:lagrangian',MT_LOG, 1, olagr))
199     $   olagr = scftype.eq.'ROHF'
200c
201      nalpha = nclosed + nopen
202      nbeta  = nclosed
203c
204c   DIIS toggle
205c
206      if (.not.rtdb_get(rtdb, 'scf:diis',MT_LOG, 1, odiis))
207     $   odiis = .false.
208c
209      call ga_sync()
210c
211c     For now set NMO = NBF, however this may change later when the
212c     linear dependency analysis is done just before the starting guess
213c
214      nmo = nbf
215c
216c     Store the derived info in the database for other wavefunction
217c     modules and/or restart to access
218c
219      if (.not. rtdb_cput(rtdb, 'scf:scftype', 1, scftype))
220     $     call errquit('scf_init: put of scftyp failed',0, RTDB_ERR)
221      if (.not. rtdb_put(rtdb, 'scf:nopen', MT_INT, 1, nopen))
222     $     call errquit('scf_init: put of nopen failed',0, RTDB_ERR)
223      if (.not. rtdb_put(rtdb, 'scf:nclosed', MT_INT, 1, nclosed))
224     $     call errquit('scf_init: put of nclosed failed',0, RTDB_ERR)
225      if (.not. rtdb_put(rtdb, 'scf:nelec', MT_INT, 1, nelec))
226     $     call errquit('scf_init: put of nelec failed',0, RTDB_ERR)
227      if (.not. rtdb_put(rtdb, 'scf:nmo', MT_INT, 1, nmo))
228     $     call errquit('scf_init: put of nmo failed',0, RTDB_ERR)
229      if (scftype .eq. 'UHF') then
230         if (.not. rtdb_put(rtdb, 'scf:nalpha', MT_INT, 1, nalpha))
231     $        call errquit('scf_init: put of nalpha failed',0, RTDB_ERR)
232         if (.not. rtdb_put(rtdb, 'scf:nbeta', MT_INT, 1, nbeta))
233     $        call errquit('scf_init: put of nbeta failed',0, RTDB_ERR)
234      endif
235c
236c     Allocate persistent local and global arrays ... these may
237c     be reallocated later when the dependency analysis is done
238c
239      if (scftype .eq. 'UHF') then
240         if (.not. ga_create(MT_DBL, nbf, nmo, 'scf_init: alpha MOs',
241     $        32, 32, g_movecs)) call errquit('scf_init: alpha MOs', 0,
242     &       GA_ERR)
243         if (.not. ga_create(MT_DBL, nbf, nmo, 'scf_init: beta MOs',
244     $        32, 32, g_movecs(2))) call errquit('scf_init: beta MOs',0,
245     &       GA_ERR)
246      else
247         if (.not. ga_create(MT_DBL, nbf, nmo, 'scf_init: MOs',
248     $        32, 32, g_movecs)) call errquit('scf_init: MOs', 0,
249     &       GA_ERR)
250      endif
251c
252      len_occ = nmo
253      if (scftype .eq. 'UHF') len_occ = nbf * 2
254      if (.not. ma_push_get(mt_dbl, len_occ, 'scf_init: mo evals',
255     $     l_eval, k_eval)) call errquit
256     $     ('scf_init: insufficient memory?', len_occ, MA_ERR)
257c
258      if (.not. ma_push_get(mt_dbl, len_occ, 'scf_init: mo occ',
259     $     l_occ, k_occ)) call errquit
260     $     ('scf_init: insufficient memory?', len_occ, MA_ERR)
261c
262      if (.not. ma_push_get(mt_int, len_occ, 'scf_init: mo irs',
263     $     l_irs, k_irs)) call errquit
264     $     ('scf_init: insufficient memory?', len_occ, MA_ERR)
265c
266      call ifill(len_occ, 1, int_mb(k_irs), 1) ! In case not adapting
267c
268c     Fill in the SCF convergence info
269c
270      call scf_get_conv_info(rtdb)
271      call scf_get_fock_param(rtdb, tol2e)
272c
273c     ----- cosmo initialization ----
274c
275      cosmo_last = .false.
276      if ( rtdb_get(rtdb,'slv:cosmo',mt_log,1,cosmo_on)) then
277         if(cosmo_on) then
278            osome = util_print('information', print_low)
279c
280            call cosmo_initialize(rtdb,geom,basis,osome)
281c
282c           Turn cosmo on, we want to run the calculation
283c           Start with gas_phase run
284c
285            cosmo_last = .true.
286            cosmo_on = .true.
287            if(.not.rtdb_get(rtdb,'cosmo_phase',mt_int,1,cosmo_phase))
288     >         cosmo_phase = 1
289         endif
290      endif
291      oinitialized = .true.
292c
293      end
294      subroutine scf_get_conv_info(rtdb)
295C$Id$
296      implicit none
297#include "mafdecls.fh"
298#include "global.fh"
299#include "rtdb.fh"
300#include "cscf.fh"
301c
302      integer rtdb
303c
304      double precision mone
305      parameter (mone = -1.0d0)
306c
307      ouser_changed_conv = .false.
308c
309      if (.not. rtdb_get(rtdb, 'scf:maxiter', MT_INT, 1, maxiter)) then
310         if (scftype .ne. 'UHF') then
311            maxiter = 30
312         else
313            maxiter = 30
314         endif
315      endif
316      if (.not.rtdb_get(rtdb, 'scf:thresh', MT_DBL, 1, gnorm_tol))
317     $     gnorm_tol = 1.d-4
318c
319c     Ensure that the default integral selection is sufficient
320c     for the request accuracy of the SCF.  However, allow user override.
321c
322      if (.not. rtdb_get(rtdb, 'scf:tol2e', MT_DBL, 1, tol2e))
323     $     tol2e = min(1.0d-7,gnorm_tol*1d-2)
324c
325      if (rtdb_get(rtdb, 'scf:level shift info', MT_DBL, 6,shifts))then
326         ouser_changed_conv = .true.
327      else
328         call dfill(6, mone, shifts, 1)
329      endif
330      if (rtdb_get(rtdb, 'scf:full hessian switch', MT_DBL, 1,
331     $     nr_gswitch)) then
332         ouser_changed_conv = .true.
333      else
334         nr_gswitch = 0.1d0
335      endif
336c
337c     Apply defaults
338c
339      if (shifts(1) .eq. -1.0d0) shifts(1) = 5.0d0
340      if (shifts(2) .eq. -1.0d0) shifts(2) = 0.5d0
341      if (shifts(3) .eq. -1.0d0) shifts(3) = 0.0d0
342      if (shifts(4) .eq. -1.0d0) shifts(4) = 0.0d0
343      if (shifts(5) .eq. -1.0d0) shifts(5) = 0.0d0
344      if (shifts(6) .eq. -1.0d0) shifts(6) = 0.0d0
345c
346      end
347      subroutine scf_get_fock_param(rtdb, tol2e)
348      implicit none
349#include "mafdecls.fh"
350#include "rtdb.fh"
351#include "cfock.fh"
352#include "global.fh"
353#include "util.fh"
354      integer rtdb
355      double precision tol2e
356c
357      logical orestart
358      logical int2e_file_open
359      external int2e_file_open
360c
361c     Load balancing information
362c
363      if (.not. rtdb_get(rtdb, 'fock:task_bf', MT_INT, 1, task_bf))
364     $     task_bf = -1
365c
366c     Dentolmax constrains dentol=tol2e/denmax which is used
367c     to screen integrals before screening integrals*density
368c     against tol2e.  Making dentolmax smaller will, up to a
369c     point, make the hessian-vector products more stable,
370c     though slower.  In strongly convergent cases dentolmax
371c     could be increased to perhaps 1d-3 with some speed gain.
372c     Dentolmax is still used even if density
373c     screening is turned off since otherwise the hessian
374c     vector products would be disastrously expensive.
375c
376c     C4H10 shows worse convergence with dentolmax > 1d-5
377c     C60 shows worsse convergence with dentolmax > 1d-6
378c
379c     Some poor convergence in the CPHF and line search can
380c     be attributed to dentolmax being too high when preconditioning
381c     with Hessian-vector products.  I think there must be a problem
382c     in the screening algorithm somewhere ... but where?
383c
384      if (.not. rtdb_get(rtdb, 'fock:dentolmax', mt_dbl, 1, dentolmax))
385     $     dentolmax = 1d-6
386c
387c     odensityscreen enables/disables density screening inside
388c     the fock build.  Turning it off can increase stability but will
389c     slow things down.
390c
391      if (.not. rtdb_get(rtdb, 'fock:densityscreen', mt_log, 1,
392     $     odensityscreen)) odensityscreen = .true.
393c
394c     Integral caching/file information
395c
396#if defined(CRAYXT) || defined(BGP) || defined(BGQ) || defined(NOIO)
397c     direct by default to avoid io
398      if (.not. rtdb_get(rtdb, 'int2e:filesize',
399     $     MT_INT, 1, filesize)) filesize = -1
400      if (.not. rtdb_get(rtdb, 'int2e:memsize',
401     $     MT_INT, 1, memsize)) memsize = -1
402#else
403      if (.not. rtdb_get(rtdb, 'int2e:filesize',
404     $     MT_INT, 1, filesize)) filesize = 0
405      if (.not. rtdb_get(rtdb, 'int2e:memsize',
406     $     MT_INT, 1, memsize)) memsize = 0
407#endif
408      if (.not. rtdb_get(rtdb, 'int2e:restart',
409     $     MT_LOG, 1, orestart)) orestart = .false.
410c
411c     The opening routine will put the .pid on the integral filename
412c     (hence even tho' parallel file, open as sequential)
413c
414      if (.not. rtdb_cget(rtdb, 'int2e:filename', 1, int2efilename))
415     $     call util_file_name('aoints',.true.,.false.,int2efilename)
416c
417      oreadfile = .false.
418      owritefile= .false.
419c
420      if (filesize.gt.0 .or. memsize.gt.0) then
421         if (.not. int2e_file_open(int2efilename, memsize, filesize,
422     $        tol2e*0.1d0, orestart)) then
423c     Probably memsize/filesize too small to be useful ... just force direct
424            filesize = -1
425            memsize = -1
426         else
427c     Managed to open a file or memory cache
428            if (orestart) then
429               oreadfile = .true.
430            else
431               owritefile = .true.
432            endif
433         endif
434      endif
435c
436c     Blocking integral interface
437c
438      if (.not. rtdb_get(rtdb, 'fock:maxquartet',
439     $     MT_INT, 1, maxquartet)) maxquartet = 10000
440      if (.not. rtdb_get(rtdb, 'fock:maxpair',
441     $     MT_INT, 1, maxquartet)) maxpair = sqrt(dble(maxquartet))
442      if (.not. rtdb_get(rtdb, 'fock:maxeri',
443     $     MT_INT, 1, maxeri)) maxeri = 1296*100 ! 100 d(6) quartets
444      if (.not. rtdb_get(rtdb, 'fock:maxscr',
445     $     MT_INT, 1, maxscr)) maxscr = 0 ! <=0 implies use texas estimate
446      if (.not. rtdb_get(rtdb, 'fock:intacc',
447     $     MT_DBL, 1, intacc)) intacc = 0.0d0 ! Default is variable
448c
449      if (.not. rtdb_get(rtdb, 'fock:replicated',
450     $     mt_log, 1, oreplicated)) oreplicated = .true.
451      if (.not. rtdb_get(rtdb, 'fock:noblock',
452     $     mt_log, 1, onoblock)) onoblock = .false.
453c
454c  Use labels by default
455c
456      if (.not. rtdb_get( rtdb, 'fock:uselabels',  MT_LOG, 1,
457     $                    oerilabel)) oerilabel = .true.
458c
459      if (ga_nodeid().eq.0 .and.
460     $     util_print('fock param',print_high)) then
461         write(6,2) task_bf, maxquartet, maxeri, maxscr, intacc,
462     $        odensityscreen, dentolmax
463 2       format(/
464     $        ' Setting fock-build task_bf   :', i8/
465     $        '                    maxquartet:', i8/
466     $        '                    maxeri    :', i8/
467     $        '                    maxscr    :', i8/
468     $        '                    intacc    :', 1p,d8.1/
469     $        '                    denscreen : ',l1/
470     $        '                    dentol    :', d8.1/)
471         call util_flush(6)
472      endif
473c
474      end
475      subroutine fock_force_direct(rtdb)
476      implicit none
477#include "cfock.fh"
478      integer rtdb
479c
480c     Call after scf_get_fock_param to force a fock build
481c     to be direct
482c
483      call fock_2e_tidy(rtdb)
484c
485      filesize = -1
486      memsize  = -1
487c
488      end
489      subroutine scf_tidy(rtdb)
490      implicit none
491#include "errquit.fh"
492#include "cscf.fh"
493#include "mafdecls.fh"
494#include "global.fh"
495#include "geom.fh"
496#include "bas.fh"
497#include "cfock.fh"
498#include "rtdb.fh"
499#include "cosmo.fh"
500      integer rtdb
501c
502      logical status
503c
504      if (oinitialized) then
505         if (.not. geom_destroy(geom)) call errquit
506     $        ('scf_tidy: geom destroy failed', 0, GEOM_ERR)
507         if (.not. bas_destroy(basis)) call errquit
508     $        ('scf_tidy: basis destroy failed',0, BASIS_ERR)
509         status = ma_pop_stack(l_irs)
510         status = ma_pop_stack(l_occ) .and. status
511         status = ma_pop_stack(l_eval) .and. status
512         if (.not. status) call errquit
513     $        ('scf_tidy: failed to free irs/occupation/evals',0,
514     &       MA_ERR)
515         if (.not. ga_destroy(g_movecs)) call errquit
516     $        ('scf_tidy: failed to free movecs',0, GA_ERR)
517         if (scftype .eq. 'UHF') then
518            if (.not. ga_destroy(g_movecs(2))) call errquit
519     $           ('scf_tidy: failed to free beta movecs',0, GA_ERR)
520         endif
521         oinitialized = .false.
522      endif
523c
524      call fock_2e_tidy(rtdb)
525c
526c     ----- cosmo cleanup and reset -----
527c
528      if ( rtdb_get(rtdb,'slv:cosmo',mt_log,1,cosmo_on)) then
529         if(cosmo_on) then
530           call cosmo_tidy(rtdb)
531           cosmo_on = .false.
532           cosmo_phase = 1
533         endif
534      endif
535c
536      end
537      subroutine fock_2e_tidy(rtdb)
538      implicit none
539#include "errquit.fh"
540#include "cfock.fh"
541#include "rtdb.fh"
542#include "mafdecls.fh"
543      integer rtdb
544c
545      logical int2e_file_close
546      external int2e_file_close
547      logical okeep
548c
549      if (.not. rtdb_get(rtdb, 'int2e:keep',
550     $     MT_LOG, 1, okeep)) okeep = .false.
551      if (.not. int2e_file_close(okeep))
552     $     call errquit('scf_tidy: closing aoints?', 0, INT_ERR)
553c
554      end
555