1#define TCGMSG
2
3
4*     ***************************
5*     *                         *
6*     *    c_semicore_xc_F      *
7*     *                         *
8*     ***************************
9
10      subroutine c_semicore_xc_F(ispin,xcp,fion)
11      implicit none
12#include "errquit.fh"
13      integer ispin
14      real*8  xcp(*)
15
16      real*8 fion(3,*)
17
18#include "bafdecls.fh"
19#include "c_semicore_common.fh"
20
21*     **** local variables ****
22      logical value
23      integer npack0,nfft3d
24      integer ii,ia,nx,ny,nz
25      real*8  sumx,sumy,sumz
26      real*8  scal1,scal2
27      integer exi(2),vxcG(2)
28      integer Gx(2),Gy(2),Gz(2)
29      integer dng(2)
30      integer dngx(2),dngy(2),dngz(2)
31      integer cngx(2),cngy(2),cngz(2)
32
33*     **** external functions ****
34      integer  ion_nion,ion_katm,c_G_indx
35      real*8   lattice_omega
36      external ion_nion,ion_katm,c_G_indx
37      external lattice_omega
38
39
40      call C3dB_nx(1,nx)
41      call C3dB_ny(1,ny)
42      call C3dB_nz(1,nz)
43      scal1 = 1.0d0/dble(nx*ny*nz)
44      scal2 = 1.0d0/lattice_omega()
45
46      call C3dB_nfft3d(1,nfft3d)
47      call Cram_npack(0,npack0)
48
49      value = BA_push_get(mt_dcpl,nfft3d,'exi', exi(2), exi(1))
50      value = value.and.
51     >        BA_push_get(mt_dbl,nfft3d,'vxcG',vxcG(2),vxcG(1))
52
53      value = value.and.
54     >        BA_push_get(mt_dbl, nfft3d,'Gx',Gx(2),Gx(1))
55      value = value.and.
56     >        BA_push_get(mt_dbl, nfft3d,'Gy',Gy(2),Gy(1))
57      value = value.and.
58     >        BA_push_get(mt_dbl, nfft3d,'Gz',Gz(2),Gz(1))
59
60      value = value.and.
61     >        BA_push_get(mt_dcpl, nfft3d,'dng',dng(2),dng(1))
62      value = value.and.
63     >        BA_push_get(mt_dcpl, nfft3d,'cngx',cngx(2),cngx(1))
64      value = value.and.
65     >        BA_push_get(mt_dcpl, nfft3d,'cngy',cngy(2),cngy(1))
66      value = value.and.
67     >        BA_push_get(mt_dcpl, nfft3d,'cngz',cngz(2),cngz(1))
68
69      value = value.and.
70     >        BA_push_get(mt_dbl, nfft3d,'dngx',dngx(2),dngx(1))
71      value = value.and.
72     >        BA_push_get(mt_dbl, nfft3d,'dngy',dngy(2),dngy(1))
73      value = value.and.
74     >        BA_push_get(mt_dbl, nfft3d,'dngz',dngz(2),dngz(1))
75
76      if (.not. value)
77     > call errquit(' c_semicore_xc_F:out of stack memory',0,MA_ERR)
78
79
80      call C3dB_rr_Sum(1,xcp(1),
81     >                   xcp(1+(ispin-1)*nfft3d),
82     >                   dbl_mb(vxcG(1)))
83c     call C3dB_r_SMul(1,0.5d0,dcpl_mb(vxcG(1)),dcpl_mb(vxcG(1)))
84
85
86*     **** define Gx, Gy, and Gz in packed space ****
87      call C3dB_t_Copy(1,dbl_mb(c_G_indx(1)),dbl_mb(Gx(1)))
88      call C3dB_t_Copy(1,dbl_mb(c_G_indx(2)),dbl_mb(Gy(1)))
89      call C3dB_t_Copy(1,dbl_mb(c_G_indx(3)),dbl_mb(Gz(1)))
90      call Cram_r_pack(0,dbl_mb(Gx(1)))
91      call Cram_r_pack(0,dbl_mb(Gy(1)))
92      call Cram_r_pack(0,dbl_mb(Gz(1)))
93
94
95      do ii=1,ion_nion()
96         ia = ion_katm(ii)
97
98         if (log_mb(semicore(1)+ia)) then
99
100*          **** structure factor and local pseudopotential ****
101           call cstrfac(ii,dcpl_mb(exi(1)))
102           call Cram_c_pack(0,dcpl_mb(exi(1)))
103
104*          **** put sqrt(core-density) at atom position ****
105           call Cram_rc_Mul(0,
106     >               dbl_mb(ncore(1)+(ia-1)*5*npack0),
107     >              dcpl_mb(exi(1)),
108     >              dcpl_mb(dng(1)))
109
110           call Cram_rc_iMul(0,dbl_mb(Gx(1)),dcpl_mb(dng(1)),
111     >                                       dcpl_mb(cngx(1)))
112           call Cram_rc_iMul(0,dbl_mb(Gy(1)),dcpl_mb(dng(1)),
113     >                                       dcpl_mb(cngy(1)))
114           call Cram_rc_iMul(0,dbl_mb(Gz(1)),dcpl_mb(dng(1)),
115     >                                       dcpl_mb(cngz(1)))
116
117*          **** put dng,dngx,dngy,dngz in real space ****
118           call Cram_c_unpack(0,dcpl_mb(dng(1)))
119           call Cram_c_unpack(0,dcpl_mb(cngx(1)))
120           call Cram_c_unpack(0,dcpl_mb(cngy(1)))
121           call Cram_c_unpack(0,dcpl_mb(cngz(1)))
122
123           !call C3dB_cr_fft3b(1,dcpl_mb(dng(1)))
124           !call C3dB_cr_fft3b(1,dcpl_mb(cngx(1)))
125           !call C3dB_cr_fft3b(1,dcpl_mb(cngy(1)))
126           !call C3dB_cr_fft3b(1,dcpl_mb(cngz(1)))
127           call C3dB_cr_pfft3b(1,0,dcpl_mb(dng(1)))
128           call C3dB_cr_pfft3b(1,0,dcpl_mb(cngx(1)))
129           call C3dB_cr_pfft3b(1,0,dcpl_mb(cngy(1)))
130           call C3dB_cr_pfft3b(1,0,dcpl_mb(cngz(1)))
131
132           call C3dB_ccr_Mul(1,dcpl_mb(dng(1)),
133     >                         dcpl_mb(cngx(1)),
134     >                         dbl_mb(dngx(1)))
135           call C3dB_ccr_Mul(1,dcpl_mb(dng(1)),
136     >                         dcpl_mb(cngy(1)),
137     >                         dbl_mb(dngy(1)))
138           call C3dB_ccr_Mul(1,dcpl_mb(dng(1)),
139     >                         dcpl_mb(cngz(1)),
140     >                         dbl_mb(dngz(1)))
141
142           call C3dB_rr_dot(1,dbl_mb(dngx(1)),dbl_mb(vxcG(1)),sumx)
143           call C3dB_rr_dot(1,dbl_mb(dngy(1)),dbl_mb(vxcG(1)),sumy)
144           call C3dB_rr_dot(1,dbl_mb(dngz(1)),dbl_mb(vxcG(1)),sumz)
145
146           fion(1,ii) = fion(1,ii) + sumx*scal2*scal1
147           fion(2,ii) = fion(2,ii) + sumy*scal2*scal1
148           fion(3,ii) = fion(3,ii) + sumz*scal2*scal1
149
150         end if
151
152      end do
153
154      value = BA_pop_stack(dngz(2))
155      value = value.and.BA_pop_stack(dngy(2))
156      value = value.and.BA_pop_stack(dngx(2))
157      value = value.and.BA_pop_stack(cngz(2))
158      value = value.and.BA_pop_stack(cngy(2))
159      value = value.and.BA_pop_stack(cngx(2))
160      value = value.and.BA_pop_stack(dng(2))
161      value = value.and.BA_pop_stack(Gz(2))
162      value = value.and.BA_pop_stack(Gy(2))
163      value = value.and.BA_pop_stack(Gx(2))
164      value = value.and.BA_pop_stack(vxcG(2))
165      value = value.and.BA_pop_stack(exi(2))
166      if (.not. value)
167     > call errquit('c_semicore_force:error popping stack',0, MA_ERR)
168
169
170      return
171      end
172
173
174c $Id$
175