1c
2c
3c  Form
4c
5c       F    = h                                core Hamiltonian
6c        ij     ij
7c
8c                              rs    rs
9c       F    =  h   +  sum ( 2J   - K   )       inactive Fock element
10c        rs      rs       i    ii    ii
11c
12c
13       subroutine mcscf_fcore( basis, nbf, nclosed, nact,
14     $                         g_movecs, g_coul, g_exch, g_fcore )
15*
16* $Id$
17*
18       implicit none
19#include "errquit.fh"
20#include "mafdecls.fh"
21#include "global.fh"
22#include "util.fh"
23       integer basis
24       integer nbf
25       integer nclosed
26       integer nact
27       integer g_movecs
28       integer g_coul
29       integer g_exch
30       integer g_fcore
31       integer nn, t, u, tu, jlo, jhi, k_j, k_k, ld1
32       logical ga_check_JKblocked
33       external ga_check_JKblocked
34c
35c  1e Hamiltononian
36c
37       if (util_print('fcore', print_never)) then
38          write(6,*) ' MOVECS in fcore '
39          call ga_print(g_movecs)
40       endif
41       call moints_1e( nbf, basis, g_movecs, g_fcore)
42       if (util_print('fcore', print_never))  call ga_print(g_fcore)
43c
44c  Inactive-active Coulomb interation
45c
46       nn = nbf*nbf
47       if (.not.ga_check_JKblocked(g_coul,nact,nbf,jlo,jhi))
48     $      call errquit('mcscf_fcore: wrong distrib. for Coulomb',0,
49     &       GA_ERR)
50       do t=1,nact
51         do u=1,t
52           tu = ((t-1)*t)/2 + u
53           if ((tu.ge.jlo).and.(tu.le.jhi)) then
54             call ga_access(g_coul,1,nn,tu,tu,k_j,ld1)
55             call mcscf_fcore01( nbf, nclosed, nact, t, u, 2.d0,
56     $                           dbl_mb(k_j), g_fcore)
57             call ga_release(g_coul,1,nn,tu,tu)
58           endif
59         enddo
60       enddo
61c
62c  Inactive-active Exchange interaction
63c
64       if (.not.ga_check_JKblocked(g_exch,nact,nbf,jlo,jhi))
65     $      call errquit('mcscf_fcore: wrong distrib. for exchange',0,
66     &       GA_ERR)
67       do t=1,nact
68         do u=1,t
69           tu = ((t-1)*t)/2 + u
70           if ((tu.ge.jlo).and.(tu.le.jhi)) then
71             call ga_access(g_exch,1,nn,tu,tu,k_k,ld1)
72             call mcscf_fcore01( nbf, nclosed, nact, t, u, -1.d0,
73     $                           dbl_mb(k_k), g_fcore)
74             call ga_release(g_exch,1,nn,tu,tu)
75           endif
76         enddo
77       enddo
78c
79c  Properly synchronize
80c
81       call ga_sync()
82c
83c  Done
84c
85       return
86       end
87
88
89
90
91
92
93       subroutine mcscf_fcore01( nbf, nclosed, nact, t, u, fact,
94     $                           xmo, g_fcore )
95       implicit none
96#include "global.fh"
97       integer nbf
98       integer nclosed
99       integer nact
100       integer t, u
101       double precision fact
102       double precision xmo(nbf,nbf)
103       integer g_fcore
104       double precision xx
105       integer i, tt, uu
106
107       xx = 0.d0
108       do i=1,nclosed
109         xx = xx + xmo(i,i)
110       enddo
111       tt = nclosed + t
112       uu = nclosed + u
113       call ga_acc(g_fcore, tt, tt, uu, uu, xx, 1, fact )
114       if (t.ne.u) call ga_acc(g_fcore, uu, uu, tt, tt, xx, 1, fact )
115       return
116       end
117
118
119
120
121
122
123
124
125
126
127
128
129
130c
131c  Evaluate inactive Fock piece separately
132c  (see mcscf_fock)
133c
134c
135       subroutine mcscf_ifock( geom, basis, nbf, nclosed, nact,
136     $                         oskel, tol2e, g_movecs, eone, etwo,
137     $                         ecore, g_ifock )
138       implicit none
139#include "errquit.fh"
140#include "mafdecls.fh"
141#include "global.fh"
142#include "msgids.fh"
143#include "rtdb.fh"
144#include "bas.fh"
145#include "util.fh"
146#include "geom.fh"
147       integer geom, basis                            ! [input] Geometry and basis handles
148       integer nbf                                    ! [input] Number of basis functions
149       integer nclosed                                ! [input] Number of closed shells
150       integer nact                                   ! [input] Number of open shells
151       logical oskel                                  ! [input] Symmetry toggle
152       double precision tol2e                         ! [input] Integral tolerance
153       integer g_movecs                               ! [input] MO coefficients
154       double precision eone, etwo, ecore             ! [output] Energy components
155       integer g_ifock                                ! [output] Inactive Fock matrix
156c
157       integer nset
158       parameter(nset=1)
159c
160       integer g_cdens, g_tmp
161       double precision e2ii
162       integer i
163       double precision xx
164       integer iv_dens(2), iv_fock(2)
165       double precision jfac(2), kfac(2)
166       integer ga_create_atom_blocked
167       external ga_create_atom_blocked
168       data jfac/1.0d0, 1.d0/, kfac/-0.5d0, -0.5d0/
169c
170c
171c
172       g_tmp = ga_create_atom_blocked(geom, basis, 'temp1')
173       g_cdens = ga_create_atom_blocked(geom, basis, 'closed dens')
174       call ga_dgemm( 'n', 't', nbf, nbf, nclosed, 2.d0,
175     $                g_movecs, g_movecs, 0.d0, g_cdens )
176c
177c One-electron component
178c
179       call ga_zero(g_ifock)
180       call int_1e_ga( basis, basis, g_ifock, 'kinetic', oskel)
181       call int_1e_ga( basis, basis, g_ifock, 'potential', oskel)
182       if (oskel) call sym_symmetrize(geom, basis, .false., g_ifock)
183       eone = ga_ddot(g_ifock,g_cdens)
184c
185c Two-electron component of the AO Fock matrices
186c
187       call ga_zero(g_tmp)
188       iv_dens(1) = g_cdens
189       iv_fock(1) = g_tmp
190       call fock_2e( geom, basis, nset, jfac, kfac, tol2e, oskel,
191     $               iv_dens, iv_fock, .false. )
192c
193c  Symmetrize Fock AO components
194c
195       if (oskel)
196     $    call sym_symmetrize(geom, basis, .false., g_tmp )
197       e2ii = ga_ddot(g_tmp,g_cdens)
198c
199c  MO integral contribution
200c
201       etwo = e2ii*0.5d0
202       if ((ga_nodeid().eq.0).and.
203     $     (util_print('energy trace',print_debug))) then
204         write(6,911) e2ii*0.5d0
205 911     format(  'energy components  eii: ',f12.6)
206       endif
207c
208c  Inactive core energy
209c
210       call ga_dadd( 1.d0, g_tmp, 1.d0, g_ifock, g_cdens )
211       call two_index_transf( g_cdens, g_movecs, g_movecs,
212     $                       g_tmp, g_ifock )
213       ecore = 0.d0
214       do i=ga_nodeid()+1,nclosed,ga_nnodes()
215         call ga_get(g_ifock,i,i,i,i,xx,1)
216         ecore = ecore + 2.d0*xx
217       enddo
218       call ga_sync()
219       call ga_dgop(msg_mcscf_ifocktrace,ecore,1,'+')
220       ecore = ecore - 0.5d0*e2ii
221c
222c  Clean up
223c
224       if (.not.ga_destroy(g_tmp))
225     $   call errquit('mcscf_etr: cannot destroy',0, GA_ERR)
226       if (.not.ga_destroy(g_cdens))
227     $   call errquit('mcscf_etr: cannot destroy',0, GA_ERR)
228       return
229       end
230