1C
2C     Initialization for charge density Coulomb fitting.
3C
4      subroutine rt_tddft_init_coulcdfit (params)
5      implicit none
6
7#include "errquit.fh"
8#include "mafdecls.fh"
9#include "stdio.fh"
10#include "global.fh"
11#include "msgids.fh"
12#include "util.fh"
13#include "cdft.fh"
14#include "rt_tddft.fh"
15#ifdef SCALAPACK
16            integer ga_cholesky,ga_llt_i
17            external ga_cholesky,ga_llt_i
18#endif
19
20
21C     == In/out ==
22      type(rt_params_t) params  !cdfit params stored in here
23
24
25C     == Parameters ==
26      integer, parameter   :: npol = 1
27      character(*), parameter :: pname = "rt_tddft_init_coulcdfit: "
28
29
30C     == External ==
31      logical dft_mem3c
32      external dft_mem3c
33
34
35C     == Variables ==
36      logical new_incore
37      integer n_batch
38      integer n3c_int, n3c_dbl, n_semi_bufs
39      integer iwhat_max
40      integer l_3ceri
41      integer k_3ceri
42      integer l_3cwhat
43      integer k_3cwhat
44      integer fd
45
46C     (for part from dft_main0d.F)
47      integer g_tmpb
48      integer info
49      logical IOLGC
50      integer lmiss
51      integer me
52
53
54      me = ga_nodeid ()
55
56
57C
58C     (ripped from dft_main0d.F)
59C
60      IOLGC = .TRUE.
61
62C      if (noio.eq.1) IOLGC = .FALSE.
63
64      if (params%nodisk) then
65         call errquit (pname//"not working with nodisk yet",0,0)
66         IOLGC = .false.
67      endif
68
69
70C      call errquit (pname//"XXX check npol, always 1?",0,0)
71C      call rt_tddft_print_warning ("CHECK THAT CD FITTING WORKING")
72
73
74      if (CDFIT)then
75c
76c        Determine the characteristics of the CD Gaussian basis set.
77c
78c
79c        Compute the matrix inverse of the CD 2-ctr ERIs.
80c
81         if (.not. ga_create(mt_dbl, nbf_cd, nbf_cd, 'CD 2cERI',
82     &                       0, nbf_cd, g_2ceri))
83     &      call errquit(pname//'Error creating g_2ceri',0,
84     &       GA_ERR)
85         call ga_zero(g_2ceri)
86         call dft_get2eri(CD_bas_han, g_2ceri,oskel)
87         if (oskel)call
88     .        sym_symmetrize(geom,cd_bas_han,.false.,g_2ceri)
89         call ga_sync()
90         if (.not. ga_duplicate(g_2ceri, g_cdinv, 'CD 2cERInv'))
91     &    call errquit(pname//'Error creating g_cdinv',0, GA_ERR)
92         call ga_zero(g_cdinv)
93         lmiss = 1
94c         if (lmiss.eq.1) then
95            call ga_zero(g_cdinv)
96            info = 0
97#if defined(PARALLEL_DIAG)
98#ifdef SCALAPACK
99            call ga_copy(g_2ceri, g_cdinv)
100            call ga_sync()
101            info= ga_cholesky('U',g_cdinv)
102#else
103            call ga_chol(g_2ceri, g_cdinv, info)
104#endif
105#else
106            call ga_chol_seq(g_2ceri, g_cdinv, info)
107#endif
108            if (info.ne.0)then
109               if (me.eq.0)then
110                  write(LuOut,*)' Problem in performing a Cholesky '
111                  write(LuOut,*)' decomposition of the 2-ctr ERI '
112                  write(LuOut,*)' matrix using CD fitting basis. '
113                  write(LuOut,*)' Attempting a diag/inverse. '
114               endif
115            endif
116            if (info.eq.0) then
117               g_tmpb = g_2ceri
118#if defined(PARALLEL_DIAG)
119#ifdef SCALAPACK
120               info = ga_llt_i('U',g_cdinv,-1)
121            if (info.ne.0)then
122               if (me.eq.0)then
123                  write(LuOut,*)' Problem in performing a Invers. '
124                  write(LuOut,*)' of the 2-ctr ERI '
125               endif
126               call ga_sync
127               call errquit(pname//'Inverse failed ',0,0)
128            endif
129
130#else
131               call ga_inverse(g_cdinv, g_tmpb)
132#endif
133#ifndef SCALAPACK
134               call ga_dgemm('T', 'N', nbf_cd, nbf_cd, nbf_cd, 1.d0,
135     &              g_tmpb, g_tmpb, 0.d0, g_cdinv)
136#endif
137#else
138               call ga_copy(g_cdinv, g_tmpb)
139               call ga_inv_seq(g_tmpb, g_cdinv)
140#endif
141            else
142               call dft_invdiag(g_2ceri, g_cdinv,
143     &                          nbf_cd)
144            endif
145#ifndef SCALAPACK
146c
147c     second build of g_2ceri needed becuase previous calls destroyed it
148c
149            call ga_zero(g_2ceri)
150            call dft_get2eri(CD_bas_han, g_2ceri,oskel)
151            if (oskel)call
152     .           sym_symmetrize(geom,cd_bas_han,.false.,g_2ceri)
153#endif
154            if (IOLGC.and.(me.eq.0)) then
155               lmiss = 0
156               call dft_invio('CDI', g_cdinv, nbf_cd, 'WRITE', lmiss)
157               if (lmiss.ne.0)call errquit
158     &         (pname//'dft_invio - abnormal write of CDI ', 0,
159     &       DISK_ERR)
160               lmiss = 0
161               call dft_invio('CD', g_2ceri, nbf_cd, 'WRITE', lmiss)
162               if (lmiss.ne.0)call errquit
163     &         (pname//'dft_invio - abnormal write of CD ', 0,
164     &       DISK_ERR)
165            endif
166c         endif
167         if (IOLGC) then
168            if (.not. ga_destroy(g_cdinv)) call errquit
169     &         (pname//'Could not destroy g_xcinv', 0, GA_ERR)
170            if (.not. ga_destroy(g_2ceri)) call errquit
171     &         (pname//'Could not destroy g_xcinv', 0, GA_ERR)
172         endif
173      endif
174C
175C     (end from dft_main0d.F)
176C
177
178
179C
180C     Set up three center integrals.  This was already done in SCF, but
181C     we have to redo it here.
182C
183      if (dft_mem3c(params%rtdb,
184     $     params%natoms, npol, .false., .false.,
185     $     n3c_int, n3c_dbl, !! n_semi_bufs,
186     $     l_3ceri, k_3ceri, l_3cwhat, k_3cwhat)) then
187
188         incore=.false.
189         call dft_3cincor (n_batch, n3c_int, int_mb(k_3cwhat),
190     $        dbl_mb(k_3cERI), n3c_dbl)
191
192         new_incore = .true.
193      else
194         new_incore = .false.
195      endif
196
197
198C
199C     Replace the "incore" flag in cdft.fh
200C
201      incore = new_incore
202
203C      if (new_incore .neqv. incore)
204C     $     call errquit (pname//"incore different from SCF", 0, 0)
205
206
207C
208C     Load params necessary for CD fitting into "params" struct.
209C
210      params%n_batch = n_batch
211      params%k_3ceri = k_3ceri
212      params%l_3ceri = l_3ceri
213      params%n3c_int = n3c_int
214      params%k_3cwhat = k_3cwhat
215      params%l_3cwhat = l_3cwhat
216      params%n3c_dbl = n3c_dbl
217      params%iwhat_max = iwhat_max
218      params%n_semi_bufs = n_semi_bufs
219      params%fd = fd
220
221
222      end subroutine rt_tddft_init_coulcdfit
223
224
225
226c $Id$
227