1c     *****************************************
2c     *                                       *
3c     *         psp_cmp_dE_Qlm_pw             *
4c     *                                       *
5c     *****************************************
6c    0.5*ncmp*ncmp + ncmp*n
7      subroutine psp_cmp_dE_Qlm_pw(ispin,move,fion)
8      implicit none
9      integer ispin
10      logical move
11      real*8  fion(3,*)
12
13#include "bafdecls.fh"
14#include "psp.fh"
15#include "errquit.fh"
16
17*     **** local variables ****
18      logical periodic,value
19      integer npack0,nfft3d,nx,ny,nz
20      integer dng_cmp(2),vcmp(2)
21      real*8  dv,scal1,omega
22
23*     **** external functions ****
24      real*8   lattice_omega
25      external lattice_omega
26      integer  control_version
27      external control_version
28
29      periodic = (control_version().eq.3)
30
31      !****************
32      !*** periodic ***
33      !****************
34      if (periodic) then
35         omega = lattice_omega()
36         call Pack_npack(0,npack0)
37         value = BA_push_get(mt_dcpl,npack0,'dng_cmp',
38     >                    dng_cmp(2),dng_cmp(1))
39         value = value.and.
40     >           BA_push_get(mt_dcpl,npack0,'vcmp',vcmp(2),vcmp(1))
41         if (.not.value)
42     >     call errquit('psp_cmp_dE_Qlm_pw:out stack',0,MA_ERR)
43
44         !*** generate vcmp ****
45         call nwpw_compcharge_gen_dn_cmp(ispin,dcpl_mb(dng_cmp(1)))
46         call coulomb_v(dcpl_mb(dng_cmp(1)),dcpl_mb(vcmp(1)))
47         call Pack_c_Copy(0,dcpl_mb(vcmp(1)),dcpl_mb(vcmp_tmp(1)))
48
49         !*** generate vcmp = vcmp+vc ****
50         call Pack_cc_Sum2(0,
51     >                     dcpl_mb(vc_tmp(1)),
52     >                     dcpl_mb(vcmp(1)))
53
54         !*** generate dEcmp/dQlm ****
55         call Pack_c_SMul1(0,omega,dcpl_mb(vcmp(1)))
56         call nwpw_compcharge_gen_dE_Qlm(1,ispin,
57     >                                   dcpl_mb(vcmp(1)),
58     >                                   dcpl_mb(vcmp(1)),
59     >                                   move,fion)
60
61         value =           BA_pop_stack(vcmp(2))
62         value = value.and.BA_pop_stack(dng_cmp(2))
63         if (.not.value)
64     >      call errquit('psp_cmp_dE_Qlm_pw:popping stack',1,MA_ERR)
65
66
67      !*****************
68      !*** aperiodic ***
69      !*****************
70      else
71         call D3dB_nx(1,nx)
72         call D3dB_ny(1,ny)
73         call D3dB_nz(1,nz)
74         scal1 = 1.0d0/dble(nx*ny*nz)
75         dv = lattice_omega()*scal1
76
77         call D3dB_nfft3d(1,nfft3d)
78         value = BA_push_get(mt_dcpl,nfft3d,'dng_cmp',
79     >                    dng_cmp(2),dng_cmp(1))
80         value = value.and.
81     >           BA_push_get(mt_dcpl,nfft3d,'vcmp',vcmp(2),vcmp(1))
82         if (.not.value)
83     >     call errquit('psp_cmp_dE_Qlm_pw:out stack',2,MA_ERR)
84         call nwpw_compcharge_gen_dn_cmp(ispin,dcpl_mb(dng_cmp(1)))
85         call Pack_c_unpack(0,dcpl_mb(dng_cmp(1)))
86         call D3dB_cr_pfft3b(1,0,dcpl_mb(dng_cmp(1)))
87
88         !*** generate vcmp ****
89         call coulomb2_v(dcpl_mb(dng_cmp(1)),dcpl_mb(vcmp(1)))
90         call D3dB_r_Copy(0,dcpl_mb(vcmp(1)),dcpl_mb(vcmp_tmp(1)))
91
92         !*** generate vcmp = vcmp+vc ****
93         call D3dB_rr_Sum2(0,
94     >                     dcpl_mb(vc_tmp(1)),
95     >                     dcpl_mb(vcmp(1)))
96         call D3dB_r_SMul1(1,dv,dcpl_mb(vcmp(1)))
97         call D3dB_rc_pfft3f(1,0,dcpl_mb(vcmp(1)))
98         call Pack_c_pack(0,dcpl_mb(vcmp(1)))
99
100         !*** generate dEcmp/dQlm ***
101         call nwpw_compcharge_gen_dE_Qlm(1,ispin,
102     >                                   dcpl_mb(vcmp(1)),
103     >                                   dcpl_mb(vcmp(1)),
104     >                                   move,fion)
105
106         value =           BA_pop_stack(vcmp(2))
107         value = value.and.BA_pop_stack(dng_cmp(2))
108         if (.not.value)
109     >      call errquit('psp_cmp_dE_Qlm_pw:popping stack',3,MA_ERR)
110      end if
111
112      return
113      end
114
115c     *****************************************
116c     *                                       *
117c     *         psp_cmp_dE_Qlm                *
118c     *                                       *
119c     *****************************************
120
121      subroutine psp_cmp_dE_Qlm(ispin,move,fion)
122      implicit none
123      integer ispin
124      logical move
125      real*8  fion(3,*)
126
127#include "bafdecls.fh"
128#include "psp.fh"
129#include "errquit.fh"
130
131*     **** local variables ****
132      logical periodic,value
133      integer npack0,nfft3d,nx,ny,nz
134      integer dng_cmp(2),dng_cmp_smooth(2),vcmp(2),vcmp_smooth(2)
135      real*8  dv,scal1,omega
136
137*     **** external functions ****
138      real*8   lattice_omega
139      external lattice_omega
140      integer  control_version
141      external control_version
142
143      periodic = (control_version().eq.3)
144
145      !****************
146      !*** periodic ***
147      !****************
148      if (periodic) then
149
150         omega = lattice_omega()
151         call Pack_npack(0,npack0)
152         value = BA_push_get(mt_dcpl,npack0,'dng_cmp',
153     >                       dng_cmp(2),dng_cmp(1))
154         value = value.and.
155     >           BA_push_get(mt_dcpl,npack0,'dng_cmp_smooth',
156     >                       dng_cmp_smooth(2),dng_cmp_smooth(1))
157         value = value.and.
158     >           BA_push_get(mt_dcpl,npack0,'vcmp',
159     >                       vcmp(2),vcmp(1))
160         value = value.and.
161     >           BA_push_get(mt_dcpl,npack0,'vcmp_smooth',
162     >                       vcmp_smooth(2),vcmp_smooth(1))
163         if (.not.value)
164     >      call errquit('psp_cmp_dE_Qlm:out stack',0,MA_ERR)
165
166         call nwpw_compcharge_gen_dn_cmp2(ispin,
167     >                                    dcpl_mb(dng_cmp(1)),
168     >                                    dcpl_mb(dng_cmp_smooth(1)))
169
170         !*** compute hartree potential of ntilde + ncmp_tilde ***
171         !*** compute hartree potential ncmp   - ncmp_tilde ***
172         call coulomb_v(dcpl_mb(dng_cmp(1)),dcpl_mb(vcmp(1)))
173         call coulomb_v(dcpl_mb(dng_cmp_smooth(1)),
174     >                  dcpl_mb(vcmp_smooth(1)))
175         call Pack_c_Copy(0,dcpl_mb(vcmp(1)),dcpl_mb(vcmp_tmp(1)))
176
177         call Pack_cc_Sub2(0,
178     >                     dcpl_mb(vcmp_smooth(1)),
179     >                     dcpl_mb(vcmp(1)))
180
181         call Pack_cc_Sum2(0,dcpl_mb(vc_tmp(1)),
182     >                       dcpl_mb(vcmp_smooth(1)))
183
184         call Pack_c_SMul1(0,omega,dcpl_mb(vcmp(1)))
185         call Pack_c_SMul1(0,omega,dcpl_mb(vcmp_smooth(1)))
186
187         call nwpw_compcharge_gen_dE_Qlm(3,ispin,
188     >                                   dcpl_mb(vcmp_smooth(1)),
189     >                                   dcpl_mb(vcmp(1)),
190     >                                   move,fion)
191
192         value =           BA_pop_stack(vcmp_smooth(2))
193         value = value.and.BA_pop_stack(vcmp(2))
194         value = value.and.BA_pop_stack(dng_cmp_smooth(2))
195         value = value.and.BA_pop_stack(dng_cmp(2))
196         if (.not.value)
197     >     call errquit('psp_cmp_dE_Qlm:pop stack',1,MA_ERR)
198
199
200      !******************
201      !*** aperiodic ****
202      !******************
203      else
204         call D3dB_nx(1,nx)
205         call D3dB_ny(1,ny)
206         call D3dB_nz(1,nz)
207         scal1 = 1.0d0/dble(nx*ny*nz)
208         dv = lattice_omega()*scal1
209
210         call D3dB_nfft3d(1,nfft3d)
211         value = BA_push_get(mt_dcpl,nfft3d,'dng_cmp',
212     >                       dng_cmp(2),dng_cmp(1))
213         value = value.and.
214     >           BA_push_get(mt_dcpl,nfft3d,'dng_cmp_smooth',
215     >                       dng_cmp_smooth(2),dng_cmp_smooth(1))
216         value = value.and.
217     >           BA_push_get(mt_dcpl,nfft3d,'vcmp',
218     >                       vcmp(2),vcmp(1))
219         value = value.and.
220     >           BA_push_get(mt_dcpl,nfft3d,'vcmp_smooth',
221     >                       vcmp_smooth(2),vcmp_smooth(1))
222         if (.not.value)
223     >      call errquit('psp_cmp_dE_Qlm2:out stack',2,MA_ERR)
224
225         call nwpw_compcharge_gen_dn_cmp2(ispin,
226     >                                    dcpl_mb(dng_cmp(1)),
227     >                                    dcpl_mb(dng_cmp_smooth(1)))
228
229         !*** dng_cmp(G),dng_cmp_smooth(G) --> dng_cmp(r),dng_cmp_smooth(r) ***
230         call Pack_c_unpack(0,dcpl_mb(dng_cmp(1)))
231         call Pack_c_unpack(0,dcpl_mb(dng_cmp_smooth(1)))
232         call D3dB_cr_pfft3b(1,0,dcpl_mb(dng_cmp(1)))
233         call D3dB_cr_pfft3b(1,0,dcpl_mb(dng_cmp_smooth(1)))
234
235         !*** generate vcmp and vcmp_smooth ***
236         call coulomb2_v(dcpl_mb(dng_cmp(1)),
237     >                   dcpl_mb(vcmp(1)))
238         call coulomb2_v(dcpl_mb(dng_cmp_smooth(1)),
239     >                   dcpl_mb(vcmp_smooth(1)))
240         call D3dB_r_Copy(1,dcpl_mb(vcmp(1)),dcpl_mb(vcmp_tmp(1)))
241
242
243         !*** vcmp        = vcmp-vcmp_smooth ***
244         !*** vcmp_smooth = vcmp_smooth + vc ***
245         call D3dB_rr_Sub2(1,
246     >                     dcpl_mb(vcmp_smooth(1)),
247     >                     dcpl_mb(vcmp(1)))
248         call D3dB_rr_Sum2(1,
249     >                     dcpl_mb(vc_tmp(1)),
250     >                     dcpl_mb(vcmp_smooth(1)))
251
252         call D3dB_r_SMul1(1,dv,dcpl_mb(vcmp(1)))
253         call D3dB_r_SMul1(1,dv,dcpl_mb(vcmp_smooth(1)))
254         call D3dB_rc_pfft3f(1,0,dcpl_mb(vcmp(1)))
255         call D3dB_rc_pfft3f(1,0,dcpl_mb(vcmp_smooth(1)))
256         call Pack_c_pack(0,dcpl_mb(vcmp(1)))
257         call Pack_c_pack(0,dcpl_mb(vcmp_smooth(1)))
258
259         call nwpw_compcharge_gen_dE_Qlm(3,ispin,
260     >                                   dcpl_mb(vcmp_smooth(1)),
261     >                                   dcpl_mb(vcmp(1)),move,fion)
262
263         value =           BA_pop_stack(vcmp_smooth(2))
264         value = value.and.BA_pop_stack(vcmp(2))
265         value = value.and.BA_pop_stack(dng_cmp_smooth(2))
266         value = value.and.BA_pop_stack(dng_cmp(2))
267         if (.not.value)
268     >      call errquit('psp_cmp_dE_Qlm:pop stack',3,MA_ERR)
269
270      end if
271
272      call nwpw_compcharge_gen_dEmult_Qlm(ispin)
273      call nwpw_compcharge_add_dEmult_Qlm(ispin)
274
275      return
276      end
277
278
279