1      Subroutine open_xc_exact_pot(g_rho_exact, g_tmp_exact,
2     &                             g_addit_exact, g_movecs_aux,g_dens,
3     &                             dif_lamda, lamda, p_lamda,
4     &                             lamda_old, Ecoul_aux)
5c
6      implicit none
7c
8#include "global.fh"
9#include "cdft.fh"
10      integer g_rho_exact(2), g_tmp_exact(2), g_addit_exact(2),
11     &        g_movecs_aux(2), g_dens(2)
12      integer i
13      double precision dif_lamda, lamda, p_lamda, lamda_old,
14     &                 Ecoul_aux(2)
15      integer  ga_create_atom_blocked
16      external ga_create_atom_blocked
17      g_rho_exact(1) = ga_create_atom_blocked(geom, AO_bas_han,
18     &                                       'ga_rho_exact_1')
19      g_tmp_exact(1) = ga_create_atom_blocked(geom, AO_bas_han,
20     &                                       'ga_tmp_exact_1')
21      g_addit_exact(1) = ga_create_atom_blocked(geom, AO_bas_han,
22     &                                     'ga_addit_exact_1')
23      g_movecs_aux(1) = ga_create_atom_blocked(geom, AO_bas_han,
24     &                                        'ga_movecs_aux_1')
25      if (ipol.eq.2) then
26        g_rho_exact(2) = ga_create_atom_blocked(geom, AO_bas_han,
27     &                                         'ga_rho_exact_2')
28        g_tmp_exact(2) = ga_create_atom_blocked(geom, AO_bas_han,
29     &                                         'ga_tmp_exact_2')
30        g_addit_exact(2) = ga_create_atom_blocked(geom,
31     &                          AO_bas_han, 'ga_addit_exact_2')
32        g_movecs_aux(2) = ga_create_atom_blocked(geom, AO_bas_han,
33     &                                        'ga_movecs_aux_2')
34      end if
35      do i = 1, ipol
36        call ga_zero(g_rho_exact(i))
37        call ga_zero(g_tmp_exact(i))
38        call ga_zero(g_addit_exact(i))
39        call ga_zero(g_movecs_aux(i))
40        call ga_copy(g_dens(i),g_rho_exact(i))
41        call ga_scale(g_rho_exact(i),-1.0d00)
42      end do
43      do i = 1, numfunc
44        xfac(i) = 0.0d00
45        cfac(i) = 0.0d00
46      end do
47      dif_lamda = 0.d00
48      lamda = 0.d00
49      p_lamda = 0.d00
50      lamda_old = 0.d00
51      do i = 1 , ipol
52        Ecoul_aux(i) = 100.0d0
53      end do
54      return
55      end
56c
57cccc
58      Subroutine xc_exact_pot(ecoul_aux, g_dens, g_vxc,
59     &                        g_rho_exact, dif_lamda, g_tmp_exact,
60     &                        g_addit_exact)
61c
62      implicit none
63#include "errquit.fh"
64c
65c     integer noc(2)
66      integer g_dens(2), g_vxc(4),
67     &        g_vxc_aux(4), g_rho_exact(2), g_rho_diff(2),
68     &        g_addit_exact(2), g_tmp_exact(2), g_rho_zero(2)
69      integer i_dif
70c
71#include "global.fh"
72#include "cdft.fh"
73c
74      double precision jfac(4),kfac(4)
75c
76      integer  ga_create_atom_blocked
77      external ga_create_atom_blocked
78c
79      double precision fer_amal, dif_lamda
80      double precision ecoul, ecoul_aux(2), ocup
81c
82      integer me,nproc
83c     double precision zero,one,onem
84c     parameter(zero=0.d0,one=1.d0,onem=-1.d0)
85      double precision tol2e
86c
87      me=ga_nodeid()
88      nproc=ga_nnodes()
89c
90      g_vxc_aux(1) = ga_create_atom_blocked(geom,ao_bas_han,
91     &                                      'g_vxc_aux_1')
92      g_rho_diff(1) = ga_create_atom_blocked(geom,ao_bas_han,
93     &                                      'g_rho_diff_1')
94      tol2e=10.d0**(-itol2e)
95      call ga_sync
96      kfac(1) = 0.0d00
97      jfac(1) = 0.0d0
98      jfac(2) = 1.0d0
99      kfac(2) = 0d0
100      ocup = 2.0d00
101      if (ipol.eq.2) ocup = 1.0d00
102      do i_dif = 1 , ipol
103        fer_amal = -1.0d00/(ocup*noc(i_dif))
104        g_vxc_aux(2) = ga_create_atom_blocked(geom,ao_bas_han,'jk')
105        g_rho_zero(1) = ga_create_atom_blocked(geom, ao_bas_han,
106     &                                        'g_rho_zero')
107        call ga_copy(g_dens(i_dif),g_rho_zero(1))
108        g_rho_zero(2) = g_rho_zero(1)
109        call fock_2e(geom, AO_bas_han, 2, jfac, kfac,
110     &               tol2e, oskel, g_rho_zero(1), g_vxc_aux(1),.false.)
111        ecoul = 0.5d00*ga_ddot(g_rho_zero(1),g_vxc_aux(2))
112        call ga_scale(g_vxc_aux(2),fer_amal)
113        call ga_dadd(1.0d00, g_vxc(i_dif), 1.0d00, g_vxc_aux(2),
114     &                 g_vxc(i_dif))
115        if (dif_lamda.gt.0.0d00) then
116           call ga_add(1.0d00, g_vxc(i_dif), 1.0d00,
117     &                 g_tmp_exact(i_dif), g_vxc(i_dif))
118           call ga_zero(g_vxc_aux(1))
119           call ga_zero(g_vxc_aux(2))
120           call ga_dadd(1.0d00, g_dens(i_dif), 1.0d00,
121     &                  g_rho_exact(i_dif), g_rho_diff(1))
122           g_rho_diff(2) = g_rho_diff(1)
123           call fock_2e(geom, AO_bas_han, 2, jfac, kfac,
124     &          tol2e, oskel, g_rho_diff(1), g_vxc_aux(1),.false.)
125           ecoul_aux(i_dif) = 0.5d0*ga_ddot(g_rho_diff(1),
126     &                                      g_vxc_aux(2))
127           call ga_scale(g_vxc_aux(2),dif_lamda)
128           call ga_copy(g_vxc_aux(2),g_addit_exact(i_dif))
129           call ga_dadd(1.0d00, g_vxc(i_dif),1.0d00,
130     &                  g_addit_exact(i_dif), g_vxc(i_dif))
131        end if
132        if (.not. ga_destroy(g_rho_zero(1))) call errquit
133     $    ('xc_exact_pot: could not detroy g_rho_zero',0, GA_ERR)
134        if (.not. ga_destroy(g_vxc_aux(2))) call errquit
135     $        ('xc_exact_pot: could not detroy jk',0, GA_ERR)
136      end do
137      if (.not. ga_destroy(g_vxc_aux(1))) call errquit
138     $             ('xc_exact_pot: could not detroy g_vxc_aux_1',0,
139     &       GA_ERR)
140      if (.not. ga_destroy(g_rho_diff(1))) call errquit
141     $             ('xc_exact_pot: could not detroy g_rho_diff_1',0,
142     &       GA_ERR)
143      return
144      end
145c
146ccccc
147      Subroutine parlam_xc_exact_pot(lamda, lamda_old, dif_lamda,
148     &                               p_lamda, incre, g_tmp_exact,
149     &                               g_addit_exact, g_movecs_aux,
150     &                               g_movecs, iter, k_eval,
151     &                               Ecoul_aux, split)
152c
153      implicit none
154c
155#include "mafdecls.fh"
156#include "global.fh"
157#include "cdft.fh"
158#include "util.fh"
159#include "stdio.fh"
160      integer g_addit_exact(2), g_movecs_aux(2), g_movecs(2),
161     &        g_tmp_exact(2), k_eval(2)
162c     integer noc(2)
163      integer split, i, incre, iter, me, nproc, isp
164      double precision reference, lamda, p_lamda, dif_lamda,
165     &                 lamda_old, Ecoul_aux(2)
166      me=ga_nodeid()
167      nproc=ga_nnodes()
168      if (iter.ge.iterations) then
169        if (me.eq.0) write(LuOut,*) ' Lamda, NO SCF CONVERGENCE',
170     &                                  lamda
171      else
172        if (me.eq.0)then
173          write(LuOut,*) ' Lamda = ',lamda
174          do isp = 1, ipol
175            write(LuOut,*) ' Restriction = ',Ecoul_aux(isp)
176            write(LuOut,*) (dbl_mb(k_eval(isp)+i-1),i=1,noc(1)+3)
177          end do
178        end if
179      end if
180      reference = 64.0d00
181      if (iter.lt.iterations) then
182        split = 0
183        call ga_copy(g_movecs(1),g_movecs_aux(1))
184        if (ipol.eq.2) call ga_copy(g_movecs(2),g_movecs_aux(2))
185        if (lamda.lt.reference) then
186          lamda_old = lamda
187          p_lamda = p_lamda + 1.0d00
188          lamda = 2.0d00**p_lamda
189          incre = 0
190        else
191          lamda_old = lamda
192          incre = incre + 1
193          lamda = reference + dble(incre)*32.0d00
194          if (lamda.eq.896.d00) lamda = 900.0d00
195        end if
196      else
197        split = split + 1
198        do i = 1, ipol
199          call ga_zero(g_addit_exact(i))
200          call ga_copy(g_movecs_aux(i),g_movecs(i))
201        end do
202        if (lamda.le.reference) then
203          if (split.eq.1) p_lamda = p_lamda - 1.0d00
204        else
205          if (split.eq.1) incre = incre - 1
206        end if
207        lamda = lamda_old + dif_lamda/dble(2*split)
208      end if
209      dif_lamda = lamda - lamda_old
210      do i = 1, ipol
211        call ga_dadd(1.0d0, g_tmp_exact(i), 1.0d0, g_addit_exact(i),
212     &               g_tmp_exact(i))
213      end do
214      return
215      end
216c
217cccccc
218      Subroutine close_xc_exact_pot(g_rho_exact, g_tmp_exact,
219     &                             g_addit_exact, g_movecs_aux)
220c
221      implicit none
222#include "errquit.fh"
223c
224#include "global.fh"
225#include "cdft.fh"
226      integer g_rho_exact(2), g_tmp_exact(2), g_addit_exact(2),
227     &        g_movecs_aux(2)
228
229      if (.not. ga_destroy(g_rho_exact(1))) call errquit
230     &       ('xc_exact_pot: could not destroy g_rho_exact(1)', 0,
231     &       GA_ERR)
232      if (.not. ga_destroy(g_tmp_exact(1))) call errquit
233     &       ('xc_exact_pot: could not destroy g_tmp_exact_1', 0,
234     &       GA_ERR)
235      if (.not. ga_destroy(g_addit_exact(1))) call errquit
236     &       ('xc_exact_pot: could not destroy g_addit_exact_1', 0,
237     &       GA_ERR)
238      if (.not. ga_destroy(g_movecs_aux(1))) call errquit
239     &       ('xc_exact_pot: could not destroy g_movecs_aux_1', 0,
240     &       GA_ERR)
241      if (ipol.eq.2) then
242        if (.not. ga_destroy(g_rho_exact(2))) call errquit
243     &     ('xc_exact_pot: could not destroy g_rho_exact(2)', 0,
244     &       GA_ERR)
245        if (.not. ga_destroy(g_tmp_exact(2))) call errquit
246     &     ('xc_exact_pot: could not destroy g_tmp_exact_2', 0,
247     &       GA_ERR)
248        if (.not. ga_destroy(g_addit_exact(2))) call errquit
249     &     ('xc_exact_pot: could not destroy g_addit_exact_2', 0,
250     &       GA_ERR)
251        if (.not. ga_destroy(g_movecs_aux(2))) call errquit
252     &       ('xc_exact_pot: could not destroy g_movecs_aux_2', 0,
253     &       GA_ERR)
254      end if
255      return
256      end
257
258
259
260
261c $Id$
262