← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/nwpw/nwpwlib/utilities/paw_utilities/nwpw_xc.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwpw/nwpwlib/utilities/paw_utilities/nwpw_xc.F
 
1
c $Id: nwpw_xc.F 21510 2011-11-11 19:59:08Z bylaska $
 
2
 
 
3
*     *********************************
 
4
*     *                               *
 
5
*     *         nwpw_xc_init          *
 
6
*     *                               *
 
7
*     *********************************
 
8
      subroutine nwpw_xc_init(nion0,nkatm0,
 
9
     >                        nprj,nbasis,n1dgrid,psp_type,lmax0,
 
10
     >                        nprj_max0,l_prj,m_prj,b_prj)
 
11
      implicit none
 
12
      integer nion0,nkatm0
 
13
      integer nprj(*),nbasis(*),n1dgrid(*),psp_type(*),lmax0(*)
 
14
      integer nprj_max0
 
15
      integer l_prj(nprj_max0,*),m_prj(nprj_max0,*),b_prj(nprj_max0,*)
 
16
 
 
17
 
 
18
#include "mafdecls.fh"
 
19
#include "errquit.fh"
 
20
#include "nwpw_xc.fh"
 
21
 
 
22
c     **** local variables ****
 
23
      logical ok
 
24
      integer ii,ia,nsize
 
25
      integer ic
 
26
      integer l,m,li,mi,lj,mj,bi,bj,iprj,jprj,lm
 
27
      integer i_p,i_t
 
28
      integer i_tlm,i_plm,i_rlm
 
29
      integer mtr_size
 
30
      integer nb
 
31
      real*8  tmp_theta,cs_theta,tmp_phi,tmp_gaunt
 
32
      real*8  angle_phi
 
33
      integer paw_vxc_size
 
34
 
 
35
c     **** external functions ****
 
36
      integer  control_gga,control_lmax_multipole
 
37
      external control_gga,control_lmax_multipole
 
38
      real*8   rtheta_lm,drtheta_lm,nwpw_gaunt,nwpw_gaunt2,nwpw_gaunt3
 
39
      real*8   rtheta_lm_div
 
40
      external rtheta_lm,drtheta_lm,nwpw_gaunt,nwpw_gaunt2,nwpw_gaunt3
 
41
      external rtheta_lm_div
 
42
 
 
43
 
 
44
      call nwpw_timing_start(4)
 
45
      ok = .true.
 
46
 
 
47
      nion  = nion0
 
48
      nkatm = nkatm0
 
49
      gga   = control_gga()
 
50
 
 
51
      nprj_max    = 0
 
52
      nbasis_max  = 0
 
53
      n1dgrid_max = 0
 
54
      do ia=1,nkatm
 
55
         if (nprj(ia)   .gt.nprj_max)    nprj_max = nprj(ia)
 
56
         if (nbasis(ia) .gt.nbasis_max)  nbasis_max = nbasis(ia)
 
57
         if (n1dgrid(ia).gt.n1dgrid_max) n1dgrid_max = n1dgrid(ia)
 
58
      end do
 
59
 
 
60
      mult_l_max = control_lmax_multipole()
 
61
      if (mult_l_max.lt.0) then
 
62
         mult_l_max = 0
 
63
         do ia=1,nkatm
 
64
            if (psp_type(ia).eq.4) then
 
65
               if (mult_l_max.lt.(2*lmax0(ia))) mult_l_max = 2*lmax0(ia)
 
66
            end if
 
67
         end do
 
68
      end if
 
69
      
 
70
*     ***paw_xc energies ***
 
71
      ok = MA_alloc_get(mt_dbl,nion,"paw_xc_e",paw_xc_e(2),paw_xc_e(1))
 
72
 
 
73
      
 
74
*     *** xc matrix arrays - nbasis_max**2 * (mult_l_max+1)**2 ***
 
75
      mtr_size = 2*(nbasis_max**2)*((mult_l_max+1)**2)
 
76
      ok = ok.and.
 
77
     >     MA_alloc_get(mt_dbl,mtr_size,"paw_xc_matr",
 
78
     >                  paw_xc_matr(2),paw_xc_matr(1))
 
79
 
 
80
      if (gga.ge.10) then
 
81
      ok = ok.and.
 
82
     >     MA_alloc_get(mt_dbl,3*mtr_size,"paw_xc_dmatr_u",
 
83
     >                  paw_xc_dmatr(2),paw_xc_dmatr(1))
 
84
      end if
 
85
      if (.not.ok)
 
86
     > call errquit("init_paw_vxc: error allocating paw_xc_matr",
 
87
     >   mtr_size,0)
 
88
 
 
89
 
 
90
      mtr_size        = 2*nprj_max*nprj_max
 
91
      ok = MA_alloc_get(mt_dbl,mtr_size,"paw_xc_pot",
 
92
     >                  paw_xc_pot(2),paw_xc_pot(1)) 
 
93
      if (.not.ok)
 
94
     > call errquit("init_paw_vxc: error allocating paw_xc_pot",
 
95
     >   mtr_size,0)
 
96
 
 
97
 
 
98
 
 
99
*     *** spherical grid arrays ***
 
100
      if(mult_l_max .eq. 0 ) then
 
101
        paw_xc_nphi   = 1
 
102
        paw_xc_ntheta = 1
 
103
      else
 
104
        paw_xc_nphi   = 3*mult_l_max
 
105
        paw_xc_ntheta = 3*mult_l_max                      
 
106
      end if
 
107
      
 
108
 
 
109
      ok = MA_alloc_get(mt_dbl,paw_xc_nphi,"paw_xc_angle_phi",
 
110
     >                  paw_xc_angle_phi(2),paw_xc_angle_phi(1))
 
111
 
 
112
      ok = ok.and. 
 
113
     >     MA_alloc_get(mt_dbl,paw_xc_ntheta,"paw_xc_cos_theta",
 
114
     >                  paw_xc_cos_theta(2),paw_xc_cos_theta(1))
 
115
 
 
116
      ok = ok.and. 
 
117
     >     MA_alloc_get(mt_dbl,paw_xc_nphi,"paw_xc_w_phi",
 
118
     >                  paw_xc_w_phi(2),paw_xc_w_phi(1))
 
119
 
 
120
      ok = ok.and.
 
121
     >     MA_alloc_get(mt_dbl,paw_xc_ntheta,"paw_xc_w_theta",
 
122
     >                  paw_xc_w_theta(2),paw_xc_w_theta(1))
 
123
 
 
124
      ok = ok.and. 
 
125
     >     MA_alloc_get(mt_dbl,
 
126
     >              paw_xc_ntheta*paw_xc_nphi*(mult_l_max+1)**2,
 
127
     >              "paw_xc_tlm",
 
128
     >              paw_xc_tlm(2),paw_xc_tlm(1))
 
129
 
 
130
 
 
131
*     **** used for generating derivatives of tlm's ****
 
132
      if (gga.ge.10) then
 
133
 
 
134
c     **** derivatives wrt to theta ****
 
135
      ok = ok.and. 
 
136
     >     MA_alloc_get(mt_dbl,
 
137
     >              paw_xc_ntheta*paw_xc_nphi*(mult_l_max+1)**2,
 
138
     >              "paw_xc_dtlm_theta",
 
139
     >              paw_xc_dtlm_theta(2),paw_xc_dtlm_theta(1))
 
140
 
 
141
c     **** derivatives wrt to phi ****
 
142
      ok = ok.and. 
 
143
     >     MA_alloc_get(mt_dbl,
 
144
     >              paw_xc_ntheta*paw_xc_nphi*(mult_l_max+1)**2,
 
145
     >              "paw_xc_dtlm_phi",
 
146
     >              paw_xc_dtlm_phi(2),paw_xc_dtlm_phi(1))
 
147
 
 
148
      end if
 
149
 
 
150
      call nwpw_get_spher_grid(paw_xc_ntheta,paw_xc_nphi,
 
151
     >                         dbl_mb(paw_xc_angle_phi(1)),
 
152
     >                         dbl_mb(paw_xc_cos_theta(1)),
 
153
     >                         dbl_mb(paw_xc_w_theta(1)),
 
154
     >                         dbl_mb(paw_xc_w_phi(1)))
 
155
 
 
156
c     **** define tlm's ****
 
157
      i_tlm = 0
 
158
      do i_t=1,paw_xc_ntheta
 
159
      do i_p=1,paw_xc_nphi
 
160
         do l=0,mult_l_max
 
161
         do m=-l,l
 
162
           tmp_theta = rtheta_lm(l,m,
 
163
     >              dbl_mb(paw_xc_cos_theta(1)+i_t-1))
 
164
           angle_phi=dbl_mb(paw_xc_angle_phi(1)+i_p-1)
 
165
           if (m.lt.0) then
 
166
              tmp_phi = dsin(abs(m)*angle_phi)
 
167
           else if (m.gt.0) then
 
168
              tmp_phi = dcos(abs(m)*angle_phi)
 
169
           else
 
170
              tmp_phi = 1.0d0
 
171
           end if
 
172
           dbl_mb(paw_xc_tlm(1)+i_tlm) = tmp_theta*tmp_phi
 
173
           i_tlm = i_tlm + 1                          
 
174
         end do
 
175
         end do
 
176
      end do
 
177
      end do
 
178
 
 
179
      if (gga.ge.10) then
 
180
 
 
181
c       **** define derivative wrt to theta ****
 
182
        i_tlm = 0
 
183
        do i_t=1,paw_xc_ntheta
 
184
        do i_p=1,paw_xc_nphi
 
185
           do l=0,mult_l_max
 
186
           do m=-l,l
 
187
              tmp_theta = drtheta_lm(l,m,
 
188
     >                               dbl_mb(paw_xc_cos_theta(1)+i_t-1))
 
189
              angle_phi=dbl_mb(paw_xc_angle_phi(1)+i_p-1)
 
190
              if (m.lt.0) then
 
191
                 tmp_phi = dsin(abs(m)*angle_phi)
 
192
              else if (m.gt.0) then
 
193
                 tmp_phi = dcos(abs(m)*angle_phi)
 
194
              else
 
195
                 tmp_phi = 1.0d0
 
196
              end if
 
197
              dbl_mb(paw_xc_dtlm_theta(1)+i_tlm) = tmp_theta*tmp_phi
 
198
              i_tlm = i_tlm + 1 
 
199
           end do
 
200
           end do
 
201
        end do
 
202
        end do
 
203
 
 
204
c       **** define derivative wrt to phi ****
 
205
        i_tlm = 0
 
206
        do i_t=1,paw_xc_ntheta
 
207
        do i_p=1,paw_xc_nphi
 
208
           do l=0,mult_l_max
 
209
           do m=-l,l
 
210
              if (m.eq.0) then
 
211
                 dbl_mb(paw_xc_dtlm_phi(1)+ i_tlm) = 0.0d0
 
212
              else
 
213
                 tmp_theta = rtheta_lm_div(l,m,
 
214
     >                              dbl_mb(paw_xc_cos_theta(1)+i_t-1))
 
215
 
 
216
                 angle_phi=dbl_mb(paw_xc_angle_phi(1)+i_p-1)
 
217
                 if (m.lt.0) then
 
218
                    tmp_phi = abs(m)*dcos(abs(m)*angle_phi)
 
219
                 else if (m.gt.0) then
 
220
                    tmp_phi = -abs(m)*dsin(abs(m)*angle_phi)
 
221
                 else
 
222
                    tmp_phi = 0.0d0
 
223
                 end if
 
224
                 dbl_mb(paw_xc_dtlm_phi(1)+ i_tlm) = tmp_theta*tmp_phi
 
225
              end if
 
226
              i_tlm = i_tlm + 1
 
227
           end do
 
228
           end do
 
229
        end do
 
230
        end do
 
231
 
 
232
      end if
 
233
 
 
234
*     *** temp arrays ***
 
235
      nsize = 2*n1dgrid_max
 
236
      ok = ok.and. 
 
237
     >     MA_alloc_get(mt_dbl,nsize,"rho_ae",rho_ae(2),rho_ae(1))
 
238
      ok = ok.and. 
 
239
     >     MA_alloc_get(mt_dbl,nsize,"rho_ps",rho_ps(2),rho_ps(1))
 
240
      
 
241
*     **** allocate gradient's and agr's ****
 
242
      if (gga.ge.10) then
 
243
       nsize = 2*3*n1dgrid_max
 
244
       ok = ok.and.
 
245
     >      MA_alloc_get(mt_dbl,nsize,
 
246
     >              "paw_tmp_drho_ae",
 
247
     >              rho_ae_prime(2),rho_ae_prime(1))
 
248
       ok = ok.and.
 
249
     >      MA_alloc_get(mt_dbl,nsize,
 
250
     >                   "paw_tmp_drho_ps",
 
251
     >                   rho_ps_prime(2),rho_ps_prime(1))
 
252
       nsize = (2*2-1)*n1dgrid_max
 
253
       ok = ok.and.
 
254
     >      MA_alloc_get(mt_dbl,nsize,"agr_ae",agr_ae(2),agr_ae(1))
 
255
       ok = ok.and.
 
256
     >      MA_alloc_get(mt_dbl,nsize,"agr_ps",agr_ps(2),agr_ps(1))
 
257
       ok = ok.and.
 
258
     >      MA_alloc_get(mt_dbl,nsize,"fdn_ae",fdn_ae(2),fdn_ae(1))
 
259
       ok = ok.and.
 
260
     >      MA_alloc_get(mt_dbl,nsize,"fdn_ps",fdn_ps(2),fdn_ps(1))
 
261
      end if
 
262
 
 
263
      nsize = 2*n1dgrid_max
 
264
      ok = ok.and. 
 
265
     >     MA_alloc_get(mt_dbl,nsize,"vxc_ae",vxc_ae(2),vxc_ae(1))
 
266
      ok = ok.and.
 
267
     >     MA_alloc_get(mt_dbl,nsize,"vxc_ps",vxc_ps(2),vxc_ps(1))
 
268
      ok = ok.and.
 
269
     >     MA_alloc_get(mt_dbl,nsize,"exc_ae",exc_ae(2),exc_ae(1))
 
270
      ok = ok.and.
 
271
     >     MA_alloc_get(mt_dbl,nsize,"exc_ps",exc_ps(2),exc_ps(1))
 
272
      ok = ok.and.
 
273
     >     MA_alloc_get(mt_dbl,nsize,"xc_temp",
 
274
     >                   xc_temp(2),xc_temp(1))
 
275
      ok = ok.and.
 
276
     >     MA_alloc_get(mt_dbl,nsize,"xc_temp2",
 
277
     >                  xc_temp2(2),xc_temp2(1))
 
278
 
 
279
 
 
280
*     *** allocate vxclm multipole expansion  arrays ****
 
281
      nsize = 2*n1dgrid_max*(mult_l_max+1)**2
 
282
 
 
283
      ok = ok.and.
 
284
     >     MA_alloc_get(mt_dbl,nsize,"paw_vxc_ae",
 
285
     >                  paw_vxc_ae(2),paw_vxc_ae(1))
 
286
      ok = ok.and.
 
287
     >     MA_alloc_get(mt_dbl,nsize,"paw_vxc_ps",
 
288
     >                  paw_vxc_ps(2),paw_vxc_ps(1))
 
289
 
 
290
*     *** allocate dvxclm multipole expansion  arrays ****
 
291
      if (gga.ge.10) then
 
292
      ok = ok.and.
 
293
     >     MA_alloc_get(mt_dbl,3*nsize,
 
294
     >                  "paw_dvxc_u_ae",
 
295
     >                  paw_dvxc_ae(2),paw_dvxc_ae(1))
 
296
      ok = ok.and.
 
297
     >     MA_alloc_get(mt_dbl,3*nsize,
 
298
     >                  "paw_dvxc_u_ps",
 
299
     >                  paw_dvxc_ps(2),paw_dvxc_ps(1))
 
300
      end if
 
301
 
 
302
*     *** allocate rholm multipole expansion  arrays ****
 
303
      ok = ok.and.
 
304
     >     MA_alloc_get(mt_dbl,nsize,"paw_rho2_ae",
 
305
     >                  paw_rho2_ae(2),paw_rho2_ae(1))
 
306
      ok = ok.and.
 
307
     >     MA_alloc_get(mt_dbl,nsize,"paw_rho2_ps",
 
308
     >                  paw_rho2_ps(2),paw_rho2_ps(1))
 
309
 
 
310
*     *** allocate rholm_prime multipole expansion arrays - need for GGA's****
 
311
      if (gga.ge.10) then
 
312
      ok = ok.and.
 
313
     >     MA_alloc_get(mt_dbl,nsize,
 
314
     >                  "paw_rho2_ae_prime",
 
315
     >                   paw_rho2_ae_prime(2),paw_rho2_ae_prime(1))
 
316
      ok = ok.and.
 
317
     >     MA_alloc_get(mt_dbl,nsize,
 
318
     >                  "paw_rho2_ps_prime",
 
319
     >                  paw_rho2_ps_prime(2),paw_rho2_ps_prime(1))
 
320
      end if
 
321
      if (.not.ok)
 
322
     > call errquit("nwpw_xc_init: error allocating work arrays",0,0)
 
323
      
 
324
c      call nwpw_init_gntxc(mult_l_max)
 
325
c      if ((gga.ge.10).and.(mult_l_max.ge.1)) then
 
326
c        call nwpw_init_gntxc2(mult_l_max)
 
327
c        call nwpw_init_gntxc3(mult_l_max)
 
328
c      end if
 
329
 
 
330
 
 
331
*      ***** set up index arrays for nwpw_xc_density_solve2 *****
 
332
      ok = MA_alloc_get(mt_int,nkatm,"nindx_rholm",
 
333
     >                  nindx_rholm(2),nindx_rholm(1))
 
334
      ok = ok.and.
 
335
     >     MA_alloc_get(mt_int,nkatm,"shift_rholm",
 
336
     >                  shift_rholm(2),shift_rholm(1))
 
337
      if (.not.ok)
 
338
     > call errquit("nwpw_xc_init: error allocating work arrays",0,0)
 
339
 
 
340
      nsize = 0
 
341
      do ia=1,nkatm
 
342
         int_mb(shift_rholm(1)+ia-1) = nsize
 
343
         do jprj = 1,nprj(ia)
 
344
            lj = l_prj(jprj,ia)
 
345
            mj = m_prj(jprj,ia)
 
346
            do iprj = 1,nprj(ia)
 
347
               li = l_prj(iprj,ia)
 
348
               mi = m_prj(iprj,ia)
 
349
               do l=0,mult_l_max
 
350
                  do m=-l,l
 
351
                      if ((l.le.(li+lj)) .and. (l.ge.abs(li-lj))) then
 
352
                        tmp_gaunt = nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
 
353
                        if (dabs(tmp_gaunt).gt.1.0d-11) then
 
354
                           nsize = nsize + 1
 
355
                        end if
 
356
                      end if
 
357
                  end do
 
358
               end do
 
359
 
 
360
            end do
 
361
         end do
 
362
        int_mb(nindx_rholm(1)+ia-1) = nsize-int_mb(shift_rholm(1)+ia-1)
 
363
      end do
 
364
 
 
365
      ok = MA_alloc_get(mt_int,nsize,"bi_rholm",bi_rholm(2),bi_rholm(1))
 
366
      ok = ok.and.
 
367
     >     MA_alloc_get(mt_int,nsize,"bj_rholm",bj_rholm(2),bj_rholm(1))
 
368
      ok = ok.and.
 
369
     >     MA_alloc_get(mt_int,nsize,"lm_rholm",lm_rholm(2),lm_rholm(1))
 
370
      ok = ok.and.
 
371
     >     MA_alloc_get(mt_int,nsize,"iprj_rholm",
 
372
     >                 iprj_rholm(2),iprj_rholm(1))
 
373
      ok = ok.and.
 
374
     >     MA_alloc_get(mt_int,nsize,"jprj_rholm",
 
375
     >                 jprj_rholm(2),jprj_rholm(1))
 
376
      ok = ok.and.
 
377
     >     MA_alloc_get(mt_dbl,nsize,"coeff_rholm",
 
378
     >                 coeff_rholm(2),coeff_rholm(1))
 
379
      if (.not.ok)
 
380
     > call errquit("nwpw_xc_init: error allocating work arrays",0,0)
 
381
 
 
382
      nsize = 0
 
383
      do ia=1,nkatm
 
384
         do jprj = 1,nprj(ia)
 
385
            lj = l_prj(jprj,ia)
 
386
            mj = m_prj(jprj,ia)
 
387
            bj = b_prj(jprj,ia)
 
388
            do iprj = 1,nprj(ia)
 
389
               li = l_prj(iprj,ia)
 
390
               mi = m_prj(iprj,ia)
 
391
               bi = b_prj(iprj,ia)
 
392
               lm = 0
 
393
               do l=0,mult_l_max
 
394
                  do m=-l,l
 
395
                     if ((l.le.(li+lj)) .and. (l.ge.abs(li-lj))) then
 
396
                        tmp_gaunt = nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
 
397
                        if (dabs(tmp_gaunt).gt.1.0d-11) then
 
398
                           dbl_mb(coeff_rholm(1)+nsize) = tmp_gaunt
 
399
                           int_mb(lm_rholm(1)+nsize)   = lm
 
400
                           int_mb(iprj_rholm(1)+nsize) = iprj
 
401
                           int_mb(jprj_rholm(1)+nsize) = jprj
 
402
                           int_mb(bi_rholm(1)+nsize)   = bi
 
403
                           int_mb(bj_rholm(1)+nsize)   = bj
 
404
                           nsize = nsize + 1
 
405
                        end if
 
406
                     end if
 
407
                     lm = lm + 1
 
408
                  end do
 
409
               end do
 
410
            end do
 
411
         end do
 
412
      end do
 
413
 
 
414
      if ((gga.ge.10).and.(mult_l_max.gt.0)) then
 
415
         ok = MA_alloc_get(mt_int,nkatm,"nindx_2rholm",
 
416
     >                  nindx_2rholm(2),nindx_2rholm(1))
 
417
         ok = ok.and.
 
418
     >        MA_alloc_get(mt_int,nkatm,"shift_2rholm",
 
419
     >                  shift_2rholm(2),shift_2rholm(1))
 
420
         if (.not.ok)
 
421
     >    call errquit("nwpw_xc_init: error allocating work arrays",0,0)
 
422
 
 
423
         nsize = 0
 
424
         do ia=1,nkatm
 
425
            int_mb(shift_2rholm(1)+ia-1) = nsize
 
426
            do jprj = 1,nprj(ia)
 
427
               lj = l_prj(jprj,ia)
 
428
               mj = m_prj(jprj,ia)
 
429
               do iprj = 1,nprj(ia)
 
430
                  li = l_prj(iprj,ia)
 
431
                  mi = m_prj(iprj,ia)
 
432
                  do l=0,mult_l_max
 
433
                     do m=-l,l
 
434
                        tmp_gaunt=nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
 
435
                        if (dabs(tmp_gaunt).gt.1.0d-11) then
 
436
                           nsize = nsize + 1
 
437
                        end if
 
438
                     end do
 
439
                  end do
 
440
               end do
 
441
            end do
 
442
            int_mb(nindx_2rholm(1)+ia-1)=nsize
 
443
     >                                  -int_mb(shift_2rholm(1)+ia-1)
 
444
         end do
 
445
         ok = MA_alloc_get(mt_int,nsize,"bi_2rholm",
 
446
     >                     bi_2rholm(2),bi_2rholm(1))
 
447
         ok = ok.and.
 
448
     >        MA_alloc_get(mt_int,nsize,"bj_2rholm",
 
449
     >                     bj_2rholm(2),bj_2rholm(1))
 
450
         ok = ok.and.
 
451
     >        MA_alloc_get(mt_int,nsize,"lm_2rholm",
 
452
     >                     lm_2rholm(2),lm_2rholm(1))
 
453
         ok = ok.and.
 
454
     >        MA_alloc_get(mt_int,nsize,"iprj_2rholm",
 
455
     >                    iprj_2rholm(2),iprj_2rholm(1))
 
456
         ok = ok.and.
 
457
     >        MA_alloc_get(mt_int,nsize,"jprj_2rholm",
 
458
     >                     jprj_2rholm(2),jprj_2rholm(1))
 
459
         ok = ok.and.
 
460
     >        MA_alloc_get(mt_dbl,nsize,"coeff_2rholm",
 
461
     >                     coeff_2rholm(2),coeff_2rholm(1))
 
462
         if (.not.ok)
 
463
     >    call errquit("nwpw_xc_init: error allocating work arrays",0,0)
 
464
 
 
465
 
 
466
         nsize = 0
 
467
         do ia=1,nkatm
 
468
            do jprj = 1,nprj(ia)
 
469
               lj = l_prj(jprj,ia)
 
470
               mj = m_prj(jprj,ia)
 
471
               bj = b_prj(jprj,ia)
 
472
               do iprj = 1,nprj(ia)
 
473
                  li = l_prj(iprj,ia)
 
474
                  mi = m_prj(iprj,ia)
 
475
                  bi = b_prj(iprj,ia)
 
476
                  lm = 0
 
477
                  do l=0,mult_l_max
 
478
                     do m=-l,l
 
479
                        tmp_gaunt=nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
 
480
                        if (dabs(tmp_gaunt).gt.1.0d-11) then
 
481
                           dbl_mb(coeff_2rholm(1)+nsize) = tmp_gaunt
 
482
                           int_mb(lm_2rholm(1)+nsize)   = lm
 
483
                           int_mb(iprj_2rholm(1)+nsize) = iprj
 
484
                           int_mb(jprj_2rholm(1)+nsize) = jprj
 
485
                           int_mb(bi_2rholm(1)+nsize)   = bi
 
486
                           int_mb(bj_2rholm(1)+nsize)   = bj
 
487
                           nsize = nsize + 1
 
488
                        end if
 
489
                        lm = lm + 1
 
490
                     end do
 
491
                  end do
 
492
               end do
 
493
            end do
 
494
         end do
 
495
 
 
496
         ok = MA_alloc_get(mt_int,nkatm,"nindx_3rholm",
 
497
     >                  nindx_3rholm(2),nindx_3rholm(1))
 
498
         ok = ok.and.
 
499
     >        MA_alloc_get(mt_int,nkatm,"shift_3rholm",
 
500
     >                  shift_3rholm(2),shift_3rholm(1))
 
501
         if (.not.ok)
 
502
     >    call errquit("nwpw_xc_init: error allocating work arrays",0,0)
 
503
 
 
504
         nsize = 0
 
505
         do ia=1,nkatm
 
506
            int_mb(shift_3rholm(1)+ia-1) = nsize
 
507
            do jprj = 1,nprj(ia)
 
508
               lj = l_prj(jprj,ia)
 
509
               mj = m_prj(jprj,ia)
 
510
               do iprj = 1,nprj(ia)
 
511
                  li = l_prj(iprj,ia)
 
512
                  mi = m_prj(iprj,ia)
 
513
                  do l=0,mult_l_max
 
514
                     do m=-l,l
 
515
                        tmp_gaunt=nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
 
516
                        if (dabs(tmp_gaunt).gt.1.0d-11) then
 
517
                           nsize = nsize + 1
 
518
                        end if
 
519
                     end do
 
520
                  end do
 
521
               end do
 
522
            end do
 
523
            int_mb(nindx_3rholm(1)+ia-1)=nsize
 
524
     >                                  -int_mb(shift_3rholm(1)+ia-1)
 
525
         end do
 
526
         ok = MA_alloc_get(mt_int,nsize,"bi_3rholm",
 
527
     >                     bi_3rholm(2),bi_3rholm(1))
 
528
         ok = ok.and.
 
529
     >        MA_alloc_get(mt_int,nsize,"bj_3rholm",
 
530
     >                     bj_3rholm(2),bj_3rholm(1))
 
531
         ok = ok.and.
 
532
     >        MA_alloc_get(mt_int,nsize,"lm_3rholm",
 
533
     >                     lm_3rholm(2),lm_3rholm(1))
 
534
         ok = ok.and.
 
535
     >        MA_alloc_get(mt_int,nsize,"iprj_3rholm",
 
536
     >                    iprj_3rholm(2),iprj_3rholm(1))
 
537
         ok = ok.and.
 
538
     >        MA_alloc_get(mt_int,nsize,"jprj_3rholm",
 
539
     >                     jprj_3rholm(2),jprj_3rholm(1))
 
540
         ok = ok.and.
 
541
     >        MA_alloc_get(mt_dbl,nsize,"coeff_3rholm",
 
542
     >                     coeff_3rholm(2),coeff_3rholm(1))
 
543
         if (.not.ok)
 
544
     >    call errquit("nwpw_xc_init: error allocating work arrays",0,0)
 
545
 
 
546
         nsize = 0
 
547
         do ia=1,nkatm
 
548
            do jprj = 1,nprj(ia)
 
549
               lj = l_prj(jprj,ia)
 
550
               mj = m_prj(jprj,ia)
 
551
               bj = b_prj(jprj,ia)
 
552
               do iprj = 1,nprj(ia)
 
553
                  li = l_prj(iprj,ia)
 
554
                  mi = m_prj(iprj,ia)
 
555
                  bi = b_prj(iprj,ia)
 
556
                  lm = 0
 
557
                  do l=0,mult_l_max
 
558
                     do m=-l,l
 
559
                        tmp_gaunt=nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
 
560
                        if (dabs(tmp_gaunt).gt.1.0d-11) then
 
561
                           dbl_mb(coeff_3rholm(1)+nsize) = tmp_gaunt
 
562
                           int_mb(lm_3rholm(1)+nsize)   = lm
 
563
                           int_mb(iprj_3rholm(1)+nsize) = iprj
 
564
                           int_mb(jprj_3rholm(1)+nsize) = jprj
 
565
                           int_mb(bi_3rholm(1)+nsize)   = bi
 
566
                           int_mb(bj_3rholm(1)+nsize)   = bj
 
567
                           nsize = nsize + 1
 
568
                        end if
 
569
                        lm = lm + 1
 
570
                     end do
 
571
                  end do
 
572
               end do
 
573
            end do
 
574
         end do
 
575
 
 
576
      end if
 
577
 
 
578
 
 
579
      call nwpw_timing_end(4)
 
580
      return
 
581
      end 
 
582
 
 
583
*     ********************************************
 
584
*     *                                          *
 
585
*     *             nwpw_xc_solve                *
 
586
*     *                                          *
 
587
*     ********************************************
 
588
      subroutine nwpw_xc_solve(ii,ia,n1dgrid,nbasis,
 
589
     >               phi_ae,phi_ps,phi_ae_prime,phi_ps_prime,
 
590
     >               core_ae,core_ps,core_ae_prime,core_ps_prime,
 
591
     >               rgrid,log_amesh,
 
592
     >               ispin,ne,nprj,sw1,sw2)
 
593
      implicit none
 
594
      integer ii,ia
 
595
      integer n1dgrid,nbasis
 
596
      real*8  phi_ae(n1dgrid,nbasis)
 
597
      real*8  phi_ps(n1dgrid,nbasis)
 
598
      real*8  phi_ae_prime(n1dgrid,nbasis)
 
599
      real*8  phi_ps_prime(n1dgrid,nbasis)
 
600
      real*8  core_ae(n1dgrid)
 
601
      real*8  core_ps(n1dgrid)
 
602
      real*8  core_ae_prime(n1dgrid)
 
603
      real*8  core_ps_prime(n1dgrid)
 
604
      real*8  rgrid(n1dgrid)
 
605
      real*8  log_amesh
 
606
 
 
607
      integer ispin,ne(2),nprj
 
608
      real*8  sw1(ne(1)+ne(2),nprj)
 
609
      real*8  sw2(ne(1)+ne(2),nprj)
 
610
 
 
611
#include "mafdecls.fh"
 
612
#include "errquit.fh"
 
613
#include "nwpw_xc.fh"
 
614
 
 
615
      logical value
 
616
      integer i,j,np,i_tlm,i_tlm1,lmax2,nsize,shift,ig
 
617
      integer i_t,i_p,lm
 
618
      real*8 exc_tmp
 
619
 
 
620
*     **** external functions ****
 
621
      real*8   log_integrate_def
 
622
      external log_integrate_def
 
623
 
 
624
      call nwpw_timing_start(4)
 
625
      call nwpw_timing_start(22)
 
626
 
 
627
*     *** index for spin down in temp arrays ***      
 
628
      lmax2  = (mult_l_max+1)**2
 
629
      dbl_mb(paw_xc_e(1)+ii-1)=0.0d0
 
630
     
 
631
*     *** zero out vxclm multipole arrays ***
 
632
      nsize = ispin*n1dgrid_max*lmax2
 
633
      call dcopy(nsize,0.0d0,0,dbl_mb(paw_vxc_ae(1)),1)
 
634
      call dcopy(nsize,0.0d0,0,dbl_mb(paw_vxc_ps(1)),1)
 
635
 
 
636
*     *** zero out dvxclm multipole arrays ***
 
637
      if (gga.ge.10) then
 
638
         call dcopy(3*nsize,0.0d0,0,dbl_mb(paw_dvxc_ae(1)),1)
 
639
         call dcopy(3*nsize,0.0d0,0,dbl_mb(paw_dvxc_ps(1)),1)
 
640
      end if
 
641
 
 
642
      i_tlm = 0
 
643
      i_tlm1 = 0
 
644
      do i_t = 1,paw_xc_ntheta
 
645
      do i_p = 1,paw_xc_nphi
 
646
      
 
647
*       *** zero out temp arrays ***
 
648
        nsize = ispin*n1dgrid_max
 
649
        call dcopy(nsize,0.0d0,0,dbl_mb(rho_ae(1)),1)
 
650
        call dcopy(nsize,0.0d0,0,dbl_mb(rho_ps(1)),1)
 
651
        call dcopy(nsize,0.0d0,0,dbl_mb(vxc_ae(1)),1)
 
652
        call dcopy(nsize,0.0d0,0,dbl_mb(vxc_ps(1)),1)
 
653
        call dcopy(nsize,0.0d0,0,dbl_mb(exc_ae(1)),1)
 
654
        call dcopy(nsize,0.0d0,0,dbl_mb(exc_ps(1)),1)
 
655
 
 
656
*       *** find rholm multipole expansion on spherical grid ***          
 
657
        shift = int_mb(shift_rholm(1)+ia-1)
 
658
        call nwpw_xc_density_solve2(ispin,ne,nprj,sw1,
 
659
     >                               n1dgrid,nbasis,lmax2,
 
660
     >                               int_mb(nindx_rholm(1)+ia-1),
 
661
     >                               int_mb(lm_rholm(1)+shift),
 
662
     >                               int_mb(iprj_rholm(1)+shift),
 
663
     >                               int_mb(jprj_rholm(1)+shift),
 
664
     >                               int_mb(bi_rholm(1)+shift),
 
665
     >                               int_mb(bj_rholm(1)+shift),
 
666
     >                               dbl_mb(coeff_rholm(1)+shift),
 
667
     >                               phi_ae,phi_ps,rgrid,
 
668
     >                               dbl_mb(paw_rho2_ae(1)),
 
669
     >                               dbl_mb(paw_rho2_ps(1)))
 
670
 
 
671
 
 
672
*       *** generate atomic densities on spherical grid ***          
 
673
        call nwpw_xc_gen_atomic_densities(n1dgrid,lmax2,ispin,
 
674
     >                                    dbl_mb(paw_rho2_ae(1)),
 
675
     >                                    dbl_mb(paw_xc_tlm(1)+i_tlm),
 
676
     >                                    core_ae,
 
677
     >                                    dbl_mb(rho_ae(1)))
 
678
        call nwpw_xc_gen_atomic_densities(n1dgrid,lmax2,ispin,
 
679
     >                                    dbl_mb(paw_rho2_ps(1)),
 
680
     >                                    dbl_mb(paw_xc_tlm(1)+i_tlm),
 
681
     >                                    core_ps,
 
682
     >                                    dbl_mb(rho_ps(1)))
 
683
 
 
684
*       ***************************************************
 
685
*       *** find exchange-correlation on spherical grid ***
 
686
*       ***************************************************
 
687
 
 
688
 
 
689
*       **** LDA functionals on spherical grid ****
 
690
        if (gga.lt.10) then
 
691
           
 
692
*          **** LDA functionals of ae and ps atomic densities ****
 
693
           call nwpw_vosko(n1dgrid,n1dgrid,ispin,
 
694
     >              dbl_mb(rho_ae(1)),        
 
695
     >              dbl_mb(exc_ae(1)),        
 
696
     >              dbl_mb(vxc_ae(1)),        
 
697
     >              dbl_mb(xc_temp(1)))       
 
698
           call nwpw_vosko(n1dgrid,n1dgrid,ispin,
 
699
     >              dbl_mb(rho_ps(1)),        
 
700
     >              dbl_mb(exc_ps(1)),        
 
701
     >              dbl_mb(vxc_ps(1)),        
 
702
     >              dbl_mb(xc_temp(1)))       
 
703
 
 
704
 
 
705
*       **** GGA functionals on spherical grid ****
 
706
        else
 
707
 
 
708
*         *** find rholm_prime multipole expansion on spherical grid ***          
 
709
          shift = int_mb(shift_rholm(1)+ia-1)
 
710
          call nwpw_xc_density_prime_solve2(ispin,ne,nprj,sw1,
 
711
     >                                  n1dgrid,nbasis,lmax2,
 
712
     >                                  int_mb(nindx_rholm(1)+ia-1),
 
713
     >                                  int_mb(lm_rholm(1)+shift),
 
714
     >                                  int_mb(iprj_rholm(1)+shift),
 
715
     >                                  int_mb(jprj_rholm(1)+shift),
 
716
     >                                  int_mb(bi_rholm(1)+shift),
 
717
     >                                  int_mb(bj_rholm(1)+shift),
 
718
     >                                  dbl_mb(coeff_rholm(1)+shift),
 
719
     >                                  phi_ae,phi_ps,
 
720
     >                                  phi_ae_prime,phi_ps_prime,rgrid,
 
721
     >                                  dbl_mb(paw_rho2_ae_prime(1)),
 
722
     >                                  dbl_mb(paw_rho2_ps_prime(1)))
 
723
 
 
724
*       *** generate the gradient of atomic densities in spherical ****
 
725
*       **** coordinates  on spherical grid                        ****          
 
726
          call nwpw_xc_gen_atomic_gradients(n1dgrid,lmax2,ispin,
 
727
     >                 dbl_mb(paw_rho2_ae(1)),
 
728
     >                 dbl_mb(paw_rho2_ae_prime(1)),
 
729
     >                 dbl_mb(paw_xc_tlm(1)+i_tlm),
 
730
     >                 dbl_mb(paw_xc_dtlm_theta(1)+i_tlm),
 
731
     >                 dbl_mb(paw_xc_dtlm_phi(1)+i_tlm),
 
732
     >                 rgrid,core_ae,core_ae_prime,
 
733
     >                 dbl_mb(rho_ae_prime(1)))
 
734
 
 
735
          call nwpw_xc_gen_atomic_gradients(n1dgrid,lmax2,ispin,
 
736
     >                 dbl_mb(paw_rho2_ps(1)),
 
737
     >                 dbl_mb(paw_rho2_ps_prime(1)),
 
738
     >                 dbl_mb(paw_xc_tlm(1)+i_tlm),
 
739
     >                 dbl_mb(paw_xc_dtlm_theta(1)+i_tlm),
 
740
     >                 dbl_mb(paw_xc_dtlm_phi(1)+i_tlm),
 
741
     >                 rgrid,core_ps,core_ps_prime,
 
742
     >                 dbl_mb(rho_ps_prime(1)))
 
743
 
 
744
*         **** generate agr ****
 
745
          call nwpw_xc_gen_atomic_agr(n1dgrid,lmax2,ispin,
 
746
     >                               dbl_mb(rho_ae_prime(1)),
 
747
     >                               dbl_mb(agr_ae(1)))
 
748
          call nwpw_xc_gen_atomic_agr(n1dgrid,lmax2,ispin,
 
749
     >                               dbl_mb(rho_ps_prime(1)),
 
750
     >                               dbl_mb(agr_ps(1)))
 
751
 
 
752
*          **** GGA functionals of ae and ps atomic densities ****
 
753
           call nwpw_gga(gga,n1dgrid,ispin,
 
754
     >                    dbl_mb(rho_ae(1)),        
 
755
     >                    dbl_mb(agr_ae(1)),        
 
756
     >                    dbl_mb(exc_ae(1)),
 
757
     >                    dbl_mb(vxc_ae(1)),    !* df/dnup, df/dndn
 
758
     >                    dbl_mb(fdn_ae(1)),    !* df/d|grad nup|,df/d|grad ndn|, df/d|grad n|
 
759
     >                    dbl_mb(xc_temp(1)))       
 
760
           call nwpw_gga(gga,n1dgrid,ispin,
 
761
     >                    dbl_mb(rho_ps(1)),        
 
762
     >                    dbl_mb(agr_ps(1)),        
 
763
     >                    dbl_mb(exc_ps(1)),
 
764
     >                    dbl_mb(vxc_ps(1)),   !* df/dnup, df/dndn
 
765
     >                    dbl_mb(fdn_ps(1)),   !* df/d|grad nup|,df/d|grad ndn|, df/d|grad n|
 
766
     >                    dbl_mb(xc_temp(1)))       
 
767
 
 
768
 
 
769
c          **** if restricted calculate df/d|grad n| *(grad n)/|grad n|  ****
 
770
*          *****   and put it in rho_prime                               ****
 
771
 
 
772
*          **** if unrestricted calculate (df/d|grad nup|* (grad nup)/|grad nup|)  ****
 
773
*          ****                         + (df/d|grad n|  * (grad n)/|grad n|)      ****
 
774
*          ****             and calculate (df/d|grad ndn|* (grad ndn)/|grad ndn|)  ****
 
775
*          ****                         + (df/d|grad n|  * (grad n)/|grad n|)      ****
 
776
*          *****   and put it in rho_prime                                         ****
 
777
           call nwpw_xc_gen_dvxc(n1dgrid,lmax2,ispin,
 
778
     >                          dbl_mb(fdn_ae(1)),
 
779
     >                          dbl_mb(agr_ae(1)),
 
780
     >                          dbl_mb(rho_ae_prime(1)))
 
781
           call nwpw_xc_gen_dvxc(n1dgrid,lmax2,ispin,
 
782
     >                          dbl_mb(fdn_ps(1)),
 
783
     >                          dbl_mb(agr_ps(1)),
 
784
     >                          dbl_mb(rho_ps_prime(1)))
 
785
 
 
786
c          **** find dvxclm multipole expansion ****
 
787
           call nwpw_xc_addto_dvxclm(n1dgrid,lmax2,ispin,
 
788
     >                          (dbl_mb(paw_xc_w_theta(1)+i_t-1)
 
789
     >                          *dbl_mb(paw_xc_w_phi(1)  +i_p-1)),
 
790
     >                          dbl_mb(paw_xc_tlm(1)+i_tlm),
 
791
     >                          dbl_mb(rho_ae_prime(1)),
 
792
     >                          dbl_mb(rho_ps_prime(1)),
 
793
     >                          dbl_mb(paw_dvxc_ae(1)),
 
794
     >                          dbl_mb(paw_dvxc_ps(1)))
 
795
 
 
796
        end if
 
797
        i_tlm = i_tlm + lmax2
 
798
 
 
799
 
 
800
*       ************************************************************************************
 
801
*       *** find vxclm the multipole expansion of atomic vxc=df/dn  or (df/dnup,df/dndn) ***          
 
802
*       ************************************************************************************
 
803
         call nwpw_xc_addto_vxclm(n1dgrid,lmax2,ispin,
 
804
     >                          (dbl_mb(paw_xc_w_theta(1)+i_t-1)
 
805
     >                          *dbl_mb(paw_xc_w_phi(1)  +i_p-1)),
 
806
     >                          dbl_mb(paw_xc_tlm(1)+i_tlm1),
 
807
     >                          dbl_mb(vxc_ae(1)),
 
808
     >                          dbl_mb(vxc_ps(1)),
 
809
     >                          dbl_mb(paw_vxc_ae(1)),
 
810
     >                          dbl_mb(paw_vxc_ps(1)))
 
811
         i_tlm1 = i_tlm1 + lmax2
 
812
 
 
813
*        ******************************************************
 
814
*        *** compute the atomic exchange correlation energy ***
 
815
*        ******************************************************
 
816
         do ig=1,n1dgrid
 
817
            dbl_mb(xc_temp(1)+ig-1)=
 
818
     >             (
 
819
     >             dbl_mb(rho_ae(1)+ig-1)  + 
 
820
     >             dbl_mb(rho_ae(1)+(ispin-1)*n1dgrid+ig-1)
 
821
     >             )*dbl_mb(exc_ae(1)+ig-1)
 
822
     >             -
 
823
     >             (
 
824
     >             dbl_mb(rho_ps(1)+ig-1)+ 
 
825
     >             dbl_mb(rho_ps(1)+(ispin-1)*n1dgrid+ig-1)
 
826
     >             )*
 
827
     >             dbl_mb(exc_ps(1)+ig-1)
 
828
         end do
 
829
         exc_tmp = log_integrate_def(0,dbl_mb(xc_temp(1)),
 
830
     >                               2,rgrid,log_amesh,n1dgrid)
 
831
         dbl_mb(paw_xc_e(1)+ii-1)=dbl_mb(paw_xc_e(1)+ii-1)+
 
832
     >                           exc_tmp*
 
833
     >                           dbl_mb(paw_xc_w_theta(1)+i_t-1)*
 
834
     >                           dbl_mb(paw_xc_w_phi(1)+i_p-1)
 
835
 
 
836
      end do !i_p
 
837
      end do !i_t
 
838
 
 
839
 
 
840
*       *********************************************************
 
841
*       ****    non-local operator computation - LDA part    ****
 
842
*       *********************************************************
 
843
 
 
844
*       **** compute (vxc^a)_nln'l'^lm radial integrals *****
 
845
        call nwpw_xc_gen_matr(n1dgrid,nbasis,lmax2,ispin,
 
846
     >                        log_amesh,
 
847
     >                        phi_ae,phi_ps,rgrid,
 
848
     >                        dbl_mb(paw_vxc_ae(1)),
 
849
     >                        dbl_mb(paw_vxc_ps(1)),
 
850
     >                        dbl_mb(xc_temp(1)),
 
851
     >                        dbl_mb(paw_xc_matr(1)))
 
852
 
 
853
 
 
854
*       **** xc potential  non-local matrix elements ****
 
855
        shift = int_mb(shift_rholm(1)+ia-1)
 
856
        call nwpw_xc_sw1tosw2(ispin,ne,nprj,sw1,sw2,
 
857
     >                        n1dgrid,nbasis,lmax2,
 
858
     >                        int_mb(nindx_rholm(1)+ia-1),
 
859
     >                        int_mb(lm_rholm(1)+shift),
 
860
     >                        int_mb(iprj_rholm(1)+shift),
 
861
     >                        int_mb(jprj_rholm(1)+shift),
 
862
     >                        int_mb(bi_rholm(1)+shift),
 
863
     >                        int_mb(bj_rholm(1)+shift),
 
864
     >                         dbl_mb(coeff_rholm(1)+shift),
 
865
     >                         dbl_mb(paw_xc_matr(1)))
 
866
 
 
867
 
 
868
 
 
869
c*       *********************************************************
 
870
c*       ****    non-local operator computation - GGA part    ****
 
871
c*       *********************************************************
 
872
c        if (gga.ge.10) then
 
873
c
 
874
c          i_phi_ae0_prime = paw_basis_i_phi_ae_prime(ia)
 
875
c          i_phi_ps0_prime = paw_basis_i_phi_ps_prime(ia)
 
876
c
 
877
c*         **** compute (dvxc^a)_nln'l'^lm radial integrals *****
 
878
c          call paw_xc_gen_dmatr(ng,nb,ic,istart,lmax2,ispin,
 
879
c     >                       log_amesh,
 
880
c     >                       dbl_mb(i_phi_ae0),
 
881
c     >                       dbl_mb(i_phi_ps0),
 
882
c     >                       dbl_mb(i_phi_ae0_prime),
 
883
c     >                       dbl_mb(i_phi_ps0_prime),
 
884
c     >                       dbl_mb(i_r+istart-1),
 
885
c     >                       dbl_mb(paw_dvxc_ae(1)),
 
886
c     >                       dbl_mb(paw_dvxc_ps(1)),
 
887
c     >                       dbl_mb(xc_c_temp(1)),
 
888
c     >                       dbl_mb(paw_xc_dmatr(1)))
 
889
c
 
890
c
 
891
c*         **** xc potential  non-local matrix elements ****
 
892
c          !i_pot0 = int_mb(i_paw_xc_pot(1) + in - 1)
 
893
c          do ms=1,ispin
 
894
c           call paw_addto_potential_gntxc(ia,
 
895
c     >           dbl_mb(paw_xc_dmatr(1)+(ms-1)*nb2*lmax2*3),
 
896
c     >           dbl_mb(paw_xc_pot(1)+i_pot0+(ms-1)*paw_xc_pot_size))
 
897
c          end do
 
898
c
 
899
cc         **** theta and phi parts ****
 
900
c          if (lmax2.gt.1) then
 
901
c
 
902
cc           **** theta derivatives ****
 
903
c            do ms=1,ispin
 
904
c             call paw_addto_potential_gntxc2(ia,
 
905
c     >         dbl_mb(paw_xc_dmatr(1)+nb2*lmax2+(ms-1)*nb2*lmax2*3),
 
906
c     >         dbl_mb(paw_xc_pot(1)+i_pot0+(ms-1)*paw_xc_pot_size))
 
907
c            end do
 
908
c
 
909
cc           **** phi derivatives ****
 
910
c            do ms=1,ispin
 
911
c             call paw_addto_potential_gntxc3(ia,
 
912
c     >         dbl_mb(paw_xc_dmatr(1)+nb2*lmax2*2+(ms-1)*nb2*lmax2*3),
 
913
c     >         dbl_mb(paw_xc_pot(1)+i_pot0+(ms-1)*paw_xc_pot_size))
 
914
c            end do
 
915
c
 
916
c
 
917
c         end if !**lmax2>1**
 
918
c
 
919
c        
 
920
c        end if !**gga**
 
921
c 
 
922
c
 
923
c      if (np.gt.1)  then
 
924
c
 
925
cc       **** xc non-local matrix elements ****
 
926
c        call D3dB_Vector_Sumall(2*paw_xc_pot(3),
 
927
c     >                          dbl_mb(paw_xc_pot(1)))
 
928
c
 
929
cc       **** atomic exchange-correlation energies ****
 
930
c        call D3dB_Vector_Sumall(paw_xc_e(3),dbl_mb(paw_xc_e(1)))
 
931
c      end if
 
932
c     
 
933
c
 
934
      call nwpw_timing_end(4)
 
935
      call nwpw_timing_end(22)
 
936
 
 
937
      return
 
938
      end 
 
939
 
 
940
 
 
941
*     *****************************************
 
942
*     *                                       *
 
943
*     *             nwpw_xc_end               *
 
944
*     *                                       *
 
945
*     *****************************************
 
946
      subroutine nwpw_xc_end()
 
947
      implicit none
 
948
 
 
949
#include "mafdecls.fh"
 
950
#include "errquit.fh"
 
951
#include "nwpw_xc.fh"
 
952
 
 
953
      logical ok
 
954
      integer ms,i
 
955
      
 
956
      call nwpw_timing_start(4)
 
957
      ok = .true.
 
958
      ok = ok.and.MA_free_heap(coeff_rholm(2))
 
959
      ok = ok.and.MA_free_heap(jprj_rholm(2))
 
960
      ok = ok.and.MA_free_heap(iprj_rholm(2))
 
961
      ok = ok.and.MA_free_heap(lm_rholm(2))
 
962
      ok = ok.and.MA_free_heap(bj_rholm(2))
 
963
      ok = ok.and.MA_free_heap(bi_rholm(2))
 
964
      ok = ok.and.MA_free_heap(shift_rholm(2))
 
965
      ok = ok.and.MA_free_heap(nindx_rholm(2))
 
966
 
 
967
 
 
968
      ok = ok.and.MA_free_heap(exc_ps(2))
 
969
      ok = ok.and.MA_free_heap(exc_ae(2))
 
970
      ok = ok.and.MA_free_heap(vxc_ps(2))
 
971
      ok = ok.and.MA_free_heap(vxc_ae(2))
 
972
      ok = ok.and.MA_free_heap(rho_ps(2))
 
973
      ok = ok.and.MA_free_heap(rho_ae(2))
 
974
      ok = ok.and.MA_free_heap(xc_temp(2))
 
975
      ok = ok.and.MA_free_heap(xc_temp2(2))
 
976
 
 
977
      ok = ok.and.MA_free_heap(paw_xc_e(2))
 
978
      ok = ok.and.MA_free_heap(paw_xc_tlm(2)) 
 
979
      ok = ok.and.MA_free_heap(paw_xc_w_theta(2) )
 
980
      ok = ok.and.MA_free_heap(paw_xc_w_phi(2))
 
981
      ok = ok.and.MA_free_heap(paw_xc_cos_theta(2)) 
 
982
      ok = ok.and.MA_free_heap(paw_xc_angle_phi(2))
 
983
 
 
984
      ok = ok.and.MA_free_heap(paw_xc_pot(2))
 
985
      ok = ok.and.MA_free_heap(paw_xc_matr(2))
 
986
      ok = ok.and.MA_free_heap(paw_vxc_ae(2))
 
987
      ok = ok.and.MA_free_heap(paw_vxc_ps(2))
 
988
      ok = ok.and.MA_free_heap(paw_rho2_ae(2))
 
989
      ok = ok.and.MA_free_heap(paw_rho2_ps(2))
 
990
 
 
991
      if (gga.ge.10) then
 
992
         ok = ok.and.MA_free_heap(paw_rho2_ae_prime(2))
 
993
         ok = ok.and.MA_free_heap(paw_rho2_ps_prime(2))
 
994
         ok = ok.and.MA_free_heap(rho_ps_prime(2))
 
995
         ok = ok.and.MA_free_heap(rho_ae_prime(2))
 
996
         ok = ok.and.MA_free_heap(agr_ae(2))
 
997
         ok = ok.and.MA_free_heap(agr_ps(2))
 
998
         ok = ok.and.MA_free_heap(fdn_ae(2))
 
999
         ok = ok.and.MA_free_heap(fdn_ps(2))
 
1000
         ok = ok.and.MA_free_heap(paw_xc_dtlm_phi(2))
 
1001
         ok = ok.and.MA_free_heap(paw_xc_dtlm_theta(2))
 
1002
         ok = ok.and.MA_free_heap(paw_dvxc_ae(2))
 
1003
         ok = ok.and.MA_free_heap(paw_dvxc_ps(2))
 
1004
         ok = ok.and.MA_free_heap(paw_xc_dmatr(2))
 
1005
      end if
 
1006
 
 
1007
      if ((gga.ge.10).and.(mult_l_max.gt.0)) then
 
1008
         ok = ok.and.MA_free_heap(coeff_2rholm(2))
 
1009
         ok = ok.and.MA_free_heap(jprj_2rholm(2))
 
1010
         ok = ok.and.MA_free_heap(iprj_2rholm(2))
 
1011
         ok = ok.and.MA_free_heap(lm_2rholm(2))
 
1012
         ok = ok.and.MA_free_heap(bj_2rholm(2))
 
1013
         ok = ok.and.MA_free_heap(bi_2rholm(2))
 
1014
         ok = ok.and.MA_free_heap(shift_2rholm(2))
 
1015
         ok = ok.and.MA_free_heap(nindx_2rholm(2))
 
1016
 
 
1017
         ok = ok.and.MA_free_heap(coeff_3rholm(2))
 
1018
         ok = ok.and.MA_free_heap(jprj_3rholm(2))
 
1019
         ok = ok.and.MA_free_heap(iprj_3rholm(2))
 
1020
         ok = ok.and.MA_free_heap(lm_3rholm(2))
 
1021
         ok = ok.and.MA_free_heap(bj_3rholm(2))
 
1022
         ok = ok.and.MA_free_heap(bi_3rholm(2))
 
1023
         ok = ok.and.MA_free_heap(shift_3rholm(2))
 
1024
         ok = ok.and.MA_free_heap(nindx_3rholm(2))
 
1025
      end if
 
1026
 
 
1027
      if (.not.ok)
 
1028
     > call errquit("nwpw_xc_end: error freeing heap",0,MA_ERR)
 
1029
 
 
1030
 
 
1031
      call nwpw_timing_end(4)
 
1032
      return
 
1033
      end 
 
1034
 
 
1035
 
 
1036
*      ******************************************
 
1037
*      *                                        *
 
1038
*      *          nwpw_xc_energy_atom           *
 
1039
*      *                                        *
 
1040
*      ******************************************
 
1041
      real*8 function nwpw_xc_energy_atom()
 
1042
      implicit none
 
1043
 
 
1044
#include "mafdecls.fh"
 
1045
#include "errquit.fh"
 
1046
#include "nwpw_xc.fh"
 
1047
 
 
1048
      integer ii
 
1049
      real*8 e
 
1050
 
 
1051
      e = 0.0d0
 
1052
      do ii=1,nion
 
1053
         e = e + dbl_mb(paw_xc_e(1)+ii-1)
 
1054
      end do
 
1055
      nwpw_xc_energy_atom = e
 
1056
      return
 
1057
      end
 
1058
 
 
1059
 
 
1060
 
 
1061
*     ******************************************
 
1062
*     *                                        *
 
1063
*     *          nwpw_xc_gen_sphere_rho        *
 
1064
*     *                                        *
 
1065
*     ******************************************
 
1066
      subroutine nwpw_xc_gen_sphere_rho(ng,lmax2,rholm,Tlm,rho)
 
1067
      implicit none
 
1068
      integer ng,lmax2
 
1069
      real*8  rholm(ng,lmax2)
 
1070
      real*8  Tlm(lmax2)
 
1071
      real*8  rho(ng)
 
1072
 
 
1073
      integer    lm,ig
 
1074
 
 
1075
      do lm=1,lmax2
 
1076
         do ig=1,ng
 
1077
           rho(ig) = rho(ig) + rholm(ig,lm)*Tlm(lm)
 
1078
         end do
 
1079
      end do
 
1080
      
 
1081
      return
 
1082
      end
 
1083
 
 
1084
 
 
1085
*     ******************************************
 
1086
*     *                                        *
 
1087
*     *          nwpw_xc_rho_div_r             *
 
1088
*     *                                        *
 
1089
*     ******************************************
 
1090
      subroutine nwpw_xc_rho_div_r(ic,r,rho)
 
1091
      implicit none
 
1092
      integer ic
 
1093
      real*8     r(*)
 
1094
      real*8     rho(*)
 
1095
 
 
1096
      integer    ig
 
1097
      do ig=1,ic
 
1098
        rho(ig) = rho(ig)/r(ig)
 
1099
      end do
 
1100
      return
 
1101
      end
 
1102
 
 
1103
 
 
1104
 
 
1105
*     ********************************************************
 
1106
*     *                                                      *
 
1107
*     *             nwpw_xc_addto_vxclm                      *
 
1108
*     *                                                      *
 
1109
*     ********************************************************
 
1110
      subroutine nwpw_xc_addto_vxclm(ic,lmax2,ispin,
 
1111
     >                              alpha,
 
1112
     >                              tlm,
 
1113
     >                              vxc_ae,vxc_ps,
 
1114
     >                              vxclm_ae,vxclm_ps)
 
1115
      implicit none
 
1116
      integer ic,lmax2,ispin
 
1117
      real*8 alpha
 
1118
      real*8 tlm(*)
 
1119
      real*8 vxc_ae(ic,ispin)
 
1120
      real*8 vxc_ps(ic,ispin)
 
1121
      real*8 vxclm_ae(ic,lmax2,ispin)
 
1122
      real*8 vxclm_ps(ic,lmax2,ispin)
 
1123
 
 
1124
*     **** local variables ****
 
1125
      integer    i,lm,ms
 
1126
 
 
1127
      do ms = 1,ispin
 
1128
         do lm = 1,lmax2
 
1129
            do i  = 1,ic
 
1130
               vxclm_ae(i,lm,ms) = vxclm_ae(i,lm,ms) 
 
1131
     >                           + vxc_ae(i,ms)*(tlm(lm)*alpha)
 
1132
               vxclm_ps(i,lm,ms) = vxclm_ps(i,lm,ms) 
 
1133
     >                          + vxc_ps(i,ms)*(tlm(lm)*alpha)
 
1134
            end do
 
1135
         end do
 
1136
      end do
 
1137
      return
 
1138
      end
 
1139
 
 
1140
*     ********************************************************
 
1141
*     *                                                      *
 
1142
*     *             nwpw_xc_gen_matr                         *
 
1143
*     *                                                      *
 
1144
*     ********************************************************
 
1145
      subroutine nwpw_xc_gen_matr(ng,nb,lmax2,ispin,log_amesh,
 
1146
     >                           phi_ae,phi_ps,r,
 
1147
     >                           vxclm_ae,vxclm_ps,
 
1148
     >                           tmpC,
 
1149
     >                           matr) 
 
1150
      implicit none
 
1151
      integer ng,nb,lmax2,ispin
 
1152
      real*8  log_amesh
 
1153
      real*8  phi_ae(ng,nb)
 
1154
      real*8  phi_ps(ng,nb)
 
1155
      real*8  r(ng)
 
1156
      real*8  vxclm_ae(ng,lmax2,ispin)
 
1157
      real*8  vxclm_ps(ng,lmax2,ispin)
 
1158
      real*8  tmpC(ng)
 
1159
      real*8 matr(nb,nb,lmax2,ispin)
 
1160
 
 
1161
*     **** local varialbles ****
 
1162
      integer ig,i,j,lm,ms
 
1163
      real*8  tmp_ae,tmp_ps
 
1164
 
 
1165
*     **** external functions ****
 
1166
      real*8   log_integrate_def
 
1167
      external log_integrate_def
 
1168
 
 
1169
      do ms=1,ispin
 
1170
         do lm=1,lmax2
 
1171
           do i=1,nb
 
1172
           do j=i,nb
 
1173
              do ig=1,ng
 
1174
                 tmp_ae = phi_ae(ig,i)
 
1175
     >                   *phi_ae(ig,j)/(r(ig)**2)
 
1176
                 tmp_ps = phi_ps(ig,i)
 
1177
     >                   *phi_ps(ig,j)/(r(ig)**2)
 
1178
                 tmpC(ig) = vxclm_ae(ig,lm,ms)*tmp_ae
 
1179
     >                    - vxclm_ps(ig,lm,ms)*tmp_ps
 
1180
              end do
 
1181
              matr(i,j,lm,ms)=log_integrate_def(0,tmpC,2,r,log_amesh,ng)
 
1182
              if (i.ne.j) matr(j,i,lm,ms) = matr(i,j,lm,ms)
 
1183
           end do
 
1184
           end do
 
1185
         end do
 
1186
      end do
 
1187
      return
 
1188
      end
 
1189
 
 
1190
 
 
1191
*     ********************************************************
 
1192
*     *                                                      *
 
1193
*     *             nwpw_xc_gen_atomic_densities             *
 
1194
*     *                                                      *
 
1195
*     ********************************************************
 
1196
      subroutine nwpw_xc_gen_atomic_densities(ng,lmax2,ispin,
 
1197
     >                                 rholm,
 
1198
     >                                 tlm,
 
1199
     >                                 rhocore,
 
1200
     >                                 rho)
 
1201
      implicit none
 
1202
      integer ng,lmax2,ispin
 
1203
      real*8 rholm(ng,lmax2,ispin)
 
1204
      real*8 tlm(lmax2)
 
1205
      real*8 rhocore(ng)
 
1206
      real*8 rho(ng,ispin)
 
1207
 
 
1208
      integer ms,i
 
1209
 
 
1210
      do ms=1,ispin
 
1211
         call nwpw_xc_gen_sphere_rho(ng,lmax2,
 
1212
     >                              rholm(1,1,ms),
 
1213
     >                              tlm,
 
1214
     >                              rho(1,ms))
 
1215
         call daxpy(ng,0.5d0,rhocore,1,rho(1,ms),1)
 
1216
         do i=1,ng
 
1217
            rho(i,ms) = dabs(rho(i,ms))
 
1218
         end do
 
1219
      end do
 
1220
 
 
1221
      return
 
1222
      end
 
1223
 
 
1224
 
 
1225
 
 
1226
*     ********************************************************
 
1227
*     *                                                      *
 
1228
*     *             nwpw_xc_gen_atomic_gradients             *
 
1229
*     *                                                      *
 
1230
*     ********************************************************
 
1231
      subroutine nwpw_xc_gen_atomic_gradients(ic,lmax2,ispin,
 
1232
     >                                       rholm,rholm_prime,
 
1233
     >                                       tlm,dtlm_theta,dtlm_phi,
 
1234
     >                                       r,rhocore,rhocore_prime,
 
1235
     >                                       rho_prime)
 
1236
      implicit none
 
1237
      integer ic,lmax2,ispin
 
1238
      real*8 rholm(ic,lmax2,ispin)
 
1239
      real*8 rholm_prime(ic,lmax2,ispin)
 
1240
      real*8 tlm(lmax2)
 
1241
      real*8 dtlm_theta(lmax2)
 
1242
      real*8 dtlm_phi(lmax2)
 
1243
      real*8 r(ic)
 
1244
      real*8 rhocore(ic)
 
1245
      real*8 rhocore_prime(ic)
 
1246
 
 
1247
      real*8 rho_prime(ic,3,ispin)
 
1248
 
 
1249
 
 
1250
      integer ms,i
 
1251
 
 
1252
      call dcopy(3*ic*ispin,0.0d0,0,rho_prime,1)
 
1253
 
 
1254
      do ms=1,ispin
 
1255
 
 
1256
*        *** find [d/dr paw_rho_ae] on spherical grid ***
 
1257
         call nwpw_xc_gen_sphere_rho(ic,lmax2,
 
1258
     >                              rholm_prime(1,1,ms),
 
1259
     >                              tlm,
 
1260
     >                              rho_prime(1,1,ms))
 
1261
 
 
1262
*        *** add core d/dr densities***
 
1263
         call daxpy(ic,0.5d0,rhocore_prime,1,rho_prime(1,1,ms),1)
 
1264
 
 
1265
 
 
1266
         !*** only computate radial derivatives if lmax2==1 ****
 
1267
         if (lmax2.gt.1) then
 
1268
 
 
1269
*          *** find  (1/r) * d/dtheta paw_rho_ae on spherical grid ***
 
1270
           call nwpw_xc_gen_sphere_rho(ic,lmax2,
 
1271
     >               rholm(1,1,ms),
 
1272
     >               dtlm_theta,
 
1273
     >               rho_prime(1,2,ms))
 
1274
           call nwpw_xc_rho_div_r(ic,r,rho_prime(1,2,ms))
 
1275
 
 
1276
*          *** find  [(r*sin(theta)) * d/dphi paw_rho_ae] on spherical grid ***
 
1277
           call nwpw_xc_gen_sphere_rho(ic,lmax2,
 
1278
     >               rholm(1,1,ms),
 
1279
     >               dtlm_phi,
 
1280
     >               rho_prime(1,3,ms))
 
1281
           call nwpw_xc_rho_div_r(ic,r,rho_prime(1,3,ms))
 
1282
 
 
1283
         end if
 
1284
 
 
1285
 
 
1286
      end do
 
1287
 
 
1288
      return
 
1289
      end
 
1290
 
 
1291
 
 
1292
*    ************************************
 
1293
*    *                                  *
 
1294
*    *      nwpw_xc_gen_atomic_agr      *
 
1295
*    *                                  *
 
1296
*    ************************************
 
1297
*
 
1298
*   This function returns  the absolute values of the gradient.
 
1299
*
 
1300
*   Entry - ic     : number of grid points
 
1301
*           ispin  : restricted/unrestricted
 
1302
*           rho_prime(ic,3,ispin) : gradient in spherical coordinates
 
1303
*                                   of atomic spin densites nup and ndn
 
1304
*
 
1305
*   Exit - agr_in(*,1): |grad n| if restricted
 
1306
*   Exit - agr_in(*,3): |grad nup|, |grad ndn|, and |grad n| if unrestricted
 
1307
 
 
1308
      subroutine nwpw_xc_gen_atomic_agr(ic,lmax2,ispin,
 
1309
     >                                 rho_prime,
 
1310
     >                                 agr)
 
1311
      implicit none
 
1312
      integer ic,lmax2,ispin
 
1313
      real*8  rho_prime(ic,3,ispin)
 
1314
      real*8  agr(ic,*)                !*(ic,2*ispin-1)
 
1315
 
 
1316
      integer ig,ms
 
1317
 
 
1318
 
 
1319
      !*** only computate radial derivatives if lmax2==1 ****
 
1320
      if (lmax2.eq.1) then
 
1321
 
 
1322
c        **** compute |grad n| ****
 
1323
         do ig=1,ic
 
1324
            agr(ig,2*ispin-1)
 
1325
     >       = dsqrt( (rho_prime(ig,1,1)+rho_prime(ig,1,ispin))**2)
 
1326
         end do
 
1327
 
 
1328
c        **** compute |grad nup| and |grad ndn| ****
 
1329
         if (ispin.eq.2) then
 
1330
           do ms=1,ispin
 
1331
           do ig=1,ic
 
1332
              agr(ig,ms) = dsqrt(rho_prime(ig,1,ms)**2)
 
1333
           end do
 
1334
           end do
 
1335
         end if
 
1336
 
 
1337
 
 
1338
 
 
1339
      !***  computate theta and phi derivatives if lmax2>1 ****
 
1340
      else
 
1341
c        **** compute |grad n| ****
 
1342
         do ig=1,ic   
 
1343
            agr(ig,2*ispin-1) 
 
1344
     >       = dsqrt( (rho_prime(ig,1,1)+rho_prime(ig,1,ispin))**2
 
1345
     >              + (rho_prime(ig,2,1)+rho_prime(ig,2,ispin))**2
 
1346
     >              + (rho_prime(ig,3,1)+rho_prime(ig,3,ispin))**2)
 
1347
         end do
 
1348
 
 
1349
c        **** compute |grad nup| and |grad ndn| ****
 
1350
         if (ispin.eq.2) then
 
1351
           do ms=1,ispin
 
1352
           do ig=1,ic
 
1353
              agr(ig,ms) = dsqrt( rho_prime(ig,1,ms)**2
 
1354
     >                          + rho_prime(ig,2,ms)**2
 
1355
     >                          + rho_prime(ig,3,ms)**2)
 
1356
           end do
 
1357
           end do
 
1358
         end if
 
1359
 
 
1360
      end if
 
1361
      return
 
1362
      end
 
1363
 
 
1364
 
 
1365
*    ************************************
 
1366
*    *                                  *
 
1367
*    *      nwpw_xc_gen_dvxc            *
 
1368
*    *                                  *
 
1369
*    ************************************
 
1370
      subroutine nwpw_xc_gen_dvxc(ic,lmax2,ispin,
 
1371
     >                           fdn,agr,rho_prime)
 
1372
      implicit none
 
1373
      integer ic,lmax2,ispin
 
1374
      real*8  fdn(ic,*)
 
1375
      real*8  agr(ic,*)
 
1376
      real*8  rho_prime(ic,3,ispin)
 
1377
 
 
1378
      integer i,j,lm,ms,jmax
 
1379
      real*8  drho1,drho2,drhoa
 
1380
 
 
1381
      !*** only computate radial derivatives if lmax2=1 ****
 
1382
      if (lmax2.eq.1) then
 
1383
         jmax = 1
 
1384
      else
 
1385
         jmax = 3
 
1386
      end if
 
1387
 
 
1388
*     *** restricted ****
 
1389
      if (ispin.eq.1) then
 
1390
         do j=1,jmax
 
1391
          do i=1,ic
 
1392
           rho_prime(i,j,1) = (rho_prime(i,j,1)+rho_prime(i,j,1))
 
1393
     >                       *(fdn(i,1)/agr(i,1))
 
1394
          end do
 
1395
         end do
 
1396
 
 
1397
*     *** unrestricted ****
 
1398
      else
 
1399
         do j=1,jmax
 
1400
          do i=1,ic
 
1401
            drho1 = rho_prime(i,j,1)
 
1402
            drho2 = rho_prime(i,j,2)
 
1403
            drhoa = drho1+drho2
 
1404
            rho_prime(i,j,1) = (drho1/agr(i,1))*fdn(i,1)
 
1405
     >                       + (drhoa/agr(i,3))*fdn(i,3)
 
1406
            rho_prime(i,j,2) = (drho2/agr(i,2))*fdn(i,2)
 
1407
     >                       + (drhoa/agr(i,3))*fdn(i,3)
 
1408
          end do
 
1409
         end do
 
1410
      end if
 
1411
      return
 
1412
      end
 
1413
 
 
1414
*    ************************************
 
1415
*    *                                  *
 
1416
*    *      nwpw_xc_addto_dvxclm        *
 
1417
*    *                                  *
 
1418
*    ************************************
 
1419
      subroutine nwpw_xc_addto_dvxclm(ic,lmax2,ispin,
 
1420
     >                              alpha,
 
1421
     >                              tlm,
 
1422
     >                              dvxc_ae,dvxc_ps,
 
1423
     >                              dvxclm_ae,dvxclm_ps)
 
1424
      implicit none
 
1425
      integer ic,lmax2,ispin
 
1426
      real*8  alpha
 
1427
      real*8  tlm(*)
 
1428
      real*8  dvxc_ae(ic,3,ispin)
 
1429
      real*8  dvxc_ps(ic,3,ispin)
 
1430
 
 
1431
      real*8  dvxclm_ae(ic,lmax2,3,ispin)
 
1432
      real*8  dvxclm_ps(ic,lmax2,3,ispin)
 
1433
 
 
1434
      integer    i,j,lm,ms,jmax
 
1435
 
 
1436
      !*** only computate radial derivatives if lmax2=1 ****
 
1437
      if (lmax2.eq.1) then
 
1438
         jmax = 1
 
1439
      else
 
1440
         jmax = 3
 
1441
      end if
 
1442
 
 
1443
      do ms = 1,ispin
 
1444
         do j  = 1,jmax
 
1445
            do lm = 1,lmax2
 
1446
               do i  = 1,ic
 
1447
                 dvxclm_ae(i,lm,j,ms) = dvxclm_ae(i,lm,j,ms)
 
1448
     >                                + dvxc_ae(i,j,ms)*(tlm(lm)*alpha)
 
1449
                 dvxclm_ps(i,lm,j,ms) = dvxclm_ps(i,lm,j,ms)
 
1450
     >                                + dvxc_ps(i,j,ms)*(tlm(lm)*alpha)
 
1451
               end do
 
1452
            end do
 
1453
         end do
 
1454
      end do
 
1455
      return
 
1456
      end
 
1457
 
 
1458
c*    ************************************
 
1459
c*    *                                  *
 
1460
c*    *      nwpw_xc_gen_dmatr           *
 
1461
c*    *                                  *
 
1462
c*    ************************************
 
1463
c
 
1464
c      subroutine nwpw_xc_gen_dmatr(ng,nb,ic,istart,lmax2,
 
1465
c     >                           ispin,log_amesh,
 
1466
c     >                           phi_ae,phi_ps,
 
1467
c     >                           phi_ae_prime,phi_ps_prime,
 
1468
c     >                           r,
 
1469
c     >                           dvxclm_ae,dvxclm_ps,
 
1470
c     >                           tmpC,
 
1471
c     >                           dmatr)
 
1472
c      implicit none
 
1473
c      integer ng,nb,ic,istart,lmax2,ispin
 
1474
c      real*8  log_amesh
 
1475
c
 
1476
c      real*8  phi_ae(ng,nb)
 
1477
c      real*8  phi_ps(ng,nb)
 
1478
c      real*8  phi_ae_prime(ng,nb)
 
1479
c      real*8  phi_ps_prime(ng,nb)
 
1480
c      real*8  r(ic)
 
1481
c      real*8  dvxclm_ae(ic,lmax2,3,ispin)
 
1482
c      real*8  dvxclm_ps(ic,lmax2,3,ispin)
 
1483
c      real*8  tmpC(ic)
 
1484
c
 
1485
c      real*8   dmatr(nb,nb,lmax2,3,ispin)
 
1486
c
 
1487
c*     **** local varialbles ****
 
1488
c      integer ig,i,j,lm,ms
 
1489
c      double precision tmp_ae,tmp_ps
 
1490
c
 
1491
c
 
1492
c*     **** radial derivative integral ****
 
1493
c      do ms=1,ispin
 
1494
c      do lm=1,lmax2
 
1495
c        do i=1,nb
 
1496
c        do j=i,nb
 
1497
c
 
1498
c          do ig=1,ic
 
1499
c            tmp_ae 
 
1500
c     >        = ( phi_ae_prime(ig+istart-1,i)*phi_ae(ig+istart-1,j)
 
1501
c     >          + phi_ae(ig+istart-1,i)*phi_ae_prime(ig+istart-1,j))
 
1502
c     >          /(r(ig)**2)
 
1503
c     >        - 2.0d0*phi_ae(ig+istart-1,i)*phi_ae(ig+istart-1,j)
 
1504
c     >          /(r(ig)**3)
 
1505
c            tmp_ps 
 
1506
c     >        = ( phi_ps_prime(ig+istart-1,i)*phi_ps(ig+istart-1,j)
 
1507
c     >          + phi_ps(ig+istart-1,i)*phi_ps_prime(ig+istart-1,j))
 
1508
c     >          /(r(ig)**2)
 
1509
c     >        - 2.0d0*phi_ps(ig+istart-1,i)*phi_ps(ig+istart-1,j)
 
1510
c     >          /(r(ig)**3)
 
1511
c
 
1512
c
 
1513
c             tmpC(ig) = dvxclm_ae(ig,lm,1,ms)*tmp_ae
 
1514
c     >                - dvxclm_ps(ig,lm,1,ms)*tmp_ps
 
1515
c
 
1516
c          end do
 
1517
c          dmatr(i,j,lm,1,ms)=log_integrate_def(0,tmpC,2,r,log_amesh,ic)
 
1518
c          if (i.ne.j) dmatr(j,i,lm,1,ms) = dmatr(i,j,lm,1,ms)
 
1519
c        end do
 
1520
c        end do
 
1521
c      end do
 
1522
c      end do
 
1523
c
 
1524
c
 
1525
c*     *** only comnpute theta and phi integrals if lmax2>1 ****
 
1526
c      if (lmax2.gt.1) then
 
1527
c
 
1528
c*     **** theta derivative integral ****
 
1529
c      do ms=1,ispin
 
1530
c      do lm=1,lmax2
 
1531
c        do i=1,nb
 
1532
c        do j=i,nb
 
1533
c
 
1534
c          do ig=1,ic
 
1535
c             tmp_ae = phi_ae(ig+istart-1,i)
 
1536
c     >               *phi_ae(ig+istart-1,j)
 
1537
c     >               /(r(ig)**3)
 
1538
c             tmp_ps = phi_ps(ig+istart-1,i)
 
1539
c     >               *phi_ps(ig+istart-1,j)
 
1540
c     >               /(r(ig)**3)
 
1541
c             tmpC(ig) = dvxclm_ae(ig,lm,2,ms)*tmp_ae
 
1542
c     >                - dvxclm_ps(ig,lm,2,ms)*tmp_ps
 
1543
c          end do
 
1544
c          dmatr(i,j,lm,2,ms) =
 
1545
c     >         log_integrate_def(0,tmpC,2,r,log_amesh,ic,istart)
 
1546
c          if (i.ne.j) dmatr(j,i,lm,2,ms) = dmatr(i,j,lm,2,ms)
 
1547
c        end do
 
1548
c        end do
 
1549
c      end do
 
1550
c      end do
 
1551
c
 
1552
c*     **** phi   derivative integral ****
 
1553
c      do ms=1,ispin
 
1554
c      do lm=1,lmax2
 
1555
c        do i=1,nb
 
1556
c        do j=i,nb
 
1557
c
 
1558
c          do ig=1,ic
 
1559
c             tmp_ae = phi_ae(ig+istart-1,i)
 
1560
c     >               *phi_ae(ig+istart-1,j)
 
1561
c     >               /(r(ig)**3)
 
1562
c             tmp_ps = phi_ps(ig+istart-1,i)
 
1563
c     >               *phi_ps(ig+istart-1,j)
 
1564
c     >               /(r(ig)**3)
 
1565
c             tmpC(ig) = dvxclm_ae(ig,lm,3,ms)*tmp_ae
 
1566
c     >                - dvxclm_ps(ig,lm,3,ms)*tmp_ps
 
1567
c          end do
 
1568
c          dmatr(i,j,lm,3,ms)
 
1569
c     >      = log_integrate_def(0,tmpC,2,r,log_amesh,ic)
 
1570
c          if (i.ne.j) dmatr(j,i,lm,3,ms) = dmatr(i,j,lm,3,ms)
 
1571
c        end do
 
1572
c        end do
 
1573
c      end do
 
1574
c      end do
 
1575
c
 
1576
c      end if
 
1577
c
 
1578
c      return
 
1579
c      end
 
1580
 
 
1581
 
 
1582
*     *****************************************************
 
1583
*     *                                                   *
 
1584
*     *          nwpw_get_spher_grid                      *
 
1585
*     *                                                   *
 
1586
*     *****************************************************
 
1587
      subroutine nwpw_get_spher_grid(ntheta,nphi,angle_phi,
 
1588
     >                          cos_theta,w_theta,w_phi)
 
1589
      implicit none
 
1590
 
 
1591
      integer ntheta
 
1592
      integer nphi 
 
1593
      double precision cos_theta(ntheta)
 
1594
      double precision angle_phi(nphi)
 
1595
      double precision w_theta(ntheta)
 
1596
      double precision w_phi(nphi)
 
1597
 
 
1598
*     *** local variables ***
 
1599
      integer i
 
1600
      real*8 pi
 
1601
 
 
1602
      pi = 4.0d0*datan(1.0d0)
 
1603
 
 
1604
*     *** gaussian quadrature angular grid for cos_theta ***
 
1605
      call nwpw_gauss_weights(-1.0d0,1.0d0,cos_theta,w_theta,ntheta)
 
1606
 
 
1607
      if (nphi.gt.1) then
 
1608
*       *** linear angular grid for angle_phi***
 
1609
        do i=1,nphi
 
1610
         angle_phi(i) = 2.0d0*pi*(i-1)/dble(nphi-1)
 
1611
         w_phi(i) = 2.0d0*pi/dble(nphi-1)
 
1612
        end do
 
1613
        w_phi(1)    = 0.5d0*w_phi(1)
 
1614
        w_phi(nphi) = w_phi(1)
 
1615
      else
 
1616
        angle_phi(1) = 0.0d0
 
1617
        w_phi(1)     = 2.0d0*pi
 
1618
      end if
 
1619
      return
 
1620
      end  
 
1621
 
 
1622
*     *****************************************************
 
1623
*     *                                                   *
 
1624
*     *          nwpw_gauss_weights                       *
 
1625
*     *                                                   *
 
1626
*     *****************************************************
 
1627
      subroutine nwpw_gauss_weights(x1, x2, x, w, n)
 
1628
      implicit none
 
1629
      integer n
 
1630
      double precision x1, x2
 
1631
      double precision x(*), w(*)
 
1632
 
 
1633
!    *** local variables ***
 
1634
      integer i, j, m, niter
 
1635
      double precision eps
 
1636
      parameter (eps = 3.0d-14)
 
1637
      double precision p1, p2, p3, pp, xl, xm, z, z1,pi
 
1638
 
 
1639
      pi = 4.0d0*datan(1.0d0)
 
1640
      m = (n + 1)/2  
 
1641
      xm = 0.5d0*(x2 + x1)
 
1642
      xl = 0.5d0*(x2 - x1)
 
1643
 
 
1644
      do i = 1, m
 
1645
         z = cos(pi*(i - 0.25d0)/(n + 0.5d0))
 
1646
         niter = 0
 
1647
    1    continue
 
1648
 
 
1649
         niter = niter + 1
 
1650
         if (niter .ge. 1000000)
 
1651
     >     call errquit('cannot converge in gauss_weights',0,1)
 
1652
 
 
1653
         p1 = 1.0d0
 
1654
         p2 = 0.0d0
 
1655
         do j = 1, n
 
1656
            p3 = p2
 
1657
            p2 = p1
 
1658
            p1 = ((2.0d0*j - 1.0d0)*z*p2 - (j - 1.0d0)*p3)/j
 
1659
         end do
 
1660
 
 
1661
         pp = n*(z*p1 - p2)/(z*z - 1.0d0)
 
1662
         z1 = z
 
1663
         z = z1 - p1/pp
 
1664
 
 
1665
         if (abs(z - z1) .gt. eps) go to 1
 
1666
         x(i) = xm - xl*z
 
1667
         x(n+1-i) = xm + xl*z
 
1668
         w(i) = 2.0d0*xl/((1.0d0 - z*z)*pp*pp)
 
1669
         w(n+1-i) = w(i)
 
1670
      end do
 
1671
      return
 
1672
      end
 
1673
 
 
1674
 
 
1675
c     *********************************************
 
1676
c     *                                           *
 
1677
c     *           nwpw_xc_density_solve2          *
 
1678
c     *                                           *
 
1679
c     *********************************************
 
1680
c
 
1681
c   Calculates the atomic density lm expansions  from the
 
1682
c   overlap coefficients.
 
1683
 
 
1684
      subroutine nwpw_xc_density_solve2(ispin,ne,nprj,sw1,
 
1685
     >                                  n1dgrid,nbasis,lmax2,
 
1686
     >                                  nindx,lm_rholm,
 
1687
     >                                  iprj_rholm,jprj_rholm,
 
1688
     >                                  bi_rholm,bj_rholm,
 
1689
     >                                  coeff_rholm,
 
1690
     >                                  phi_ae,phi_ps,
 
1691
     >                                  rgrid,
 
1692
     >                                  rholm_ae,rholm_ps)
 
1693
      implicit none
 
1694
      integer ispin,ne(2),nprj
 
1695
      real*8  sw1(ne(1)+ne(2),nprj)
 
1696
      integer n1dgrid,nbasis,lmax2
 
1697
      integer nindx,lm_rholm(*),iprj_rholm(*),jprj_rholm(*)
 
1698
      integer bi_rholm(*),bj_rholm(*)
 
1699
      real*8  coeff_rholm(*)
 
1700
      real*8  phi_ae(n1dgrid,nbasis),phi_ps(n1dgrid,nbasis)
 
1701
      real*8  rgrid(n1dgrid)
 
1702
      real*8  rholm_ae(n1dgrid,lmax2,ispin)
 
1703
      real*8  rholm_ps(n1dgrid,lmax2,ispin)
 
1704
 
 
1705
      integer ms,n,i,lm,iprj,jprj,bi,bj,n1(2),n2(2)
 
1706
      real*8  tmp_ms,coeff,w,scal
 
1707
 
 
1708
*     **** external functions ****
 
1709
      real*8   lattice_omega
 
1710
      external lattice_omega
 
1711
      
 
1712
      call nwpw_timing_start(21)
 
1713
      n1(1) = 1 
 
1714
      n1(2) = ne(1)+1 
 
1715
      n2(1) = ne(1)
 
1716
      n2(2) = ne(1)+ne(2)
 
1717
      scal = 1.0d0/lattice_omega()
 
1718
 
 
1719
*     ***init to zero***
 
1720
      call dcopy(ispin*lmax2*n1dgrid,0.0d0,0,rholm_ae,1)
 
1721
      call dcopy(ispin*lmax2*n1dgrid,0.0d0,0,rholm_ps,1)
 
1722
 
 
1723
      do i=1,nindx
 
1724
         lm    = lm_rholm(i)
 
1725
         iprj  = iprj_rholm(i)
 
1726
         jprj  = jprj_rholm(i)
 
1727
         bi    = bi_rholm(i)
 
1728
         bj    = bj_rholm(i)
 
1729
         coeff = coeff_rholm(i)*scal
 
1730
         do ms=1,ispin
 
1731
            if (ne(ms).gt.0) then
 
1732
               w = 0.0d0
 
1733
               do n=n1(ms),n2(ms)
 
1734
                  w = w + sw1(n,iprj)*sw1(n,jprj)
 
1735
               end do
 
1736
               tmp_ms = w*coeff
 
1737
               call nwpw_xc_density_gen_rho(n1dgrid,tmp_ms,
 
1738
     >                  rgrid,phi_ae(1,bi),phi_ae(1,bj),
 
1739
     >                  rholm_ae(1,lm+1,ms))
 
1740
               call nwpw_xc_density_gen_rho(n1dgrid,tmp_ms,
 
1741
     >                  rgrid,phi_ps(1,bi),phi_ps(1,bj),
 
1742
     >                  rholm_ps(1,lm+1,ms))
 
1743
            end if
 
1744
         end do
 
1745
      end do
 
1746
 
 
1747
 
 
1748
      call nwpw_timing_end(21)
 
1749
      return
 
1750
      end
 
1751
 
 
1752
      subroutine nwpw_xc_density_gen_rho(ic,alpha,r,phi1,phi2,rho)
 
1753
      implicit none
 
1754
      integer ic
 
1755
      double precision alpha
 
1756
      double precision r(ic)
 
1757
      double precision phi1(ic)
 
1758
      double precision phi2(ic)
 
1759
      double precision rho(ic)
 
1760
 
 
1761
*     **** local variables ****
 
1762
      integer i
 
1763
      double precision tmp
 
1764
 
 
1765
      do i=1,ic
 
1766
         tmp = (phi1(i)*phi2(i))/(r(i)**2)
 
1767
         rho(i) = rho(i) + alpha*tmp
 
1768
      end do
 
1769
 
 
1770
      return
 
1771
      end
 
1772
 
 
1773
 
 
1774
c     *********************************************
 
1775
c     *                                           *
 
1776
c     *        nwpw_xc_density_prime_solve2       *
 
1777
c     *                                           *
 
1778
c     *********************************************
 
1779
c
 
1780
c   Calculates the atomic density lm expansions  from the
 
1781
c   overlap coefficients.
 
1782
      subroutine nwpw_xc_density_prime_solve2(ispin,ne,nprj,sw1,
 
1783
     >                                  n1dgrid,nbasis,lmax2,
 
1784
     >                                  nindx,lm_rholm,
 
1785
     >                                  iprj_rholm,jprj_rholm,
 
1786
     >                                  bi_rholm,bj_rholm,
 
1787
     >                                  coeff_rholm,
 
1788
     >                                  phi_ae,phi_ps,
 
1789
     >                                  dphi_ae,dphi_ps,rgrid,
 
1790
     >                                  rholm_ae_prime,rholm_ps_prime)
 
1791
      implicit none
 
1792
      integer ispin,ne(2),nprj
 
1793
      real*8  sw1(ne(1)+ne(2),nprj)
 
1794
      integer n1dgrid,nbasis,lmax2
 
1795
      integer nindx,lm_rholm(*),iprj_rholm(*),jprj_rholm(*)
 
1796
      integer bi_rholm(*),bj_rholm(*)
 
1797
      real*8  coeff_rholm(*)
 
1798
      real*8  phi_ae(n1dgrid,nbasis),phi_ps(n1dgrid,nbasis)
 
1799
      real*8  dphi_ae(n1dgrid,nbasis),dphi_ps(n1dgrid,nbasis)
 
1800
      real*8  rgrid(n1dgrid)
 
1801
      real*8  rholm_ae_prime(n1dgrid,lmax2,ispin)
 
1802
      real*8  rholm_ps_prime(n1dgrid,lmax2,ispin)
 
1803
 
 
1804
      integer ms,n,i,lm,iprj,jprj,bi,bj,n1(2),n2(2)
 
1805
      real*8  tmp_ms,coeff,w,scal
 
1806
 
 
1807
*     **** external functions ****
 
1808
      real*8   lattice_omega
 
1809
      external lattice_omega
 
1810
 
 
1811
      call nwpw_timing_start(21)
 
1812
      n1(1) = 1
 
1813
      n1(2) = ne(1)+1
 
1814
      n2(1) = ne(1)
 
1815
      n2(2) = ne(1)+ne(2)
 
1816
      scal = 1.0d0/lattice_omega()
 
1817
 
 
1818
*     ***init to zero***
 
1819
      call dcopy(n1dgrid*lmax2*ispin,0.0d0,0,rholm_ae_prime,1)
 
1820
      call dcopy(n1dgrid*lmax2*ispin,0.0d0,0,rholm_ps_prime,1)
 
1821
 
 
1822
      do i=1,nindx
 
1823
         lm    = lm_rholm(i)
 
1824
         iprj  = iprj_rholm(i)
 
1825
         jprj  = jprj_rholm(i)
 
1826
         bi    = bi_rholm(i)
 
1827
         bj    = bj_rholm(i)
 
1828
         coeff = coeff_rholm(i)*scal
 
1829
         do ms=1,ispin
 
1830
            if (ne(ms).gt.0) then
 
1831
               w = 0.0d0
 
1832
               do n=n1(ms),n2(ms)
 
1833
                  w = w + sw1(n,iprj)*sw1(n,jprj)
 
1834
               end do
 
1835
               tmp_ms = w*coeff
 
1836
               call nwpw_xc_density_gen_drho(n1dgrid,tmp_ms,
 
1837
     >                  rgrid,
 
1838
     >                  phi_ae(1,bi),dphi_ae(1,bi),
 
1839
     >                  phi_ae(1,bj),dphi_ae(1,bj),
 
1840
     >                  rholm_ae_prime(1,lm+1,ms))
 
1841
               call nwpw_xc_density_gen_drho(n1dgrid,tmp_ms,
 
1842
     >                  rgrid,
 
1843
     >                  phi_ps(1,bi),dphi_ps(1,bi),
 
1844
     >                  phi_ps(1,bj),dphi_ps(1,bj),
 
1845
     >                  rholm_ps_prime(1,lm+1,ms))
 
1846
            end if
 
1847
         end do
 
1848
      end do
 
1849
 
 
1850
      call nwpw_timing_end(21)
 
1851
      return
 
1852
      end
 
1853
 
 
1854
      subroutine nwpw_xc_density_gen_drho(ic,alpha,r,
 
1855
     >                                phi1,dphi1,
 
1856
     >                                phi2,dphi2,
 
1857
     >                                drho)
 
1858
      implicit none
 
1859
      integer ic
 
1860
      double precision alpha
 
1861
      double precision r(ic)
 
1862
      double precision phi1(ic),dphi1(ic)
 
1863
      double precision phi2(ic),dphi2(ic)
 
1864
      double precision drho(ic)
 
1865
 
 
1866
*     **** local variables ****
 
1867
      integer i
 
1868
      double precision tmp
 
1869
 
 
1870
      do i=1,ic
 
1871
         tmp = (dphi1(i)*phi2(i)+phi1(i)*dphi2(i))/(r(i)**2)
 
1872
     >       - 2.0d0*(phi1(i)*phi2(i))/(r(i)**3)
 
1873
         drho(i) = drho(i) + alpha*tmp
 
1874
      end do
 
1875
      return
 
1876
      end
 
1877
 
 
1878
 
 
1879
c     *********************************************
 
1880
c     *                                           *
 
1881
c     *           nwpw_xc_sw1tosw2                *
 
1882
c     *                                           *
 
1883
c     *********************************************
 
1884
c
 
1885
      subroutine nwpw_xc_sw1tosw2(ispin,ne,nprj,sw1,sw2,
 
1886
     >                            n1dgrid,nbasis,lmax2,
 
1887
     >                            nindx,lm_rholm,
 
1888
     >                            iprj_rholm,jprj_rholm,
 
1889
     >                            bi_rholm,bj_rholm,
 
1890
     >                            coeff_rholm,
 
1891
     >                            matr)
 
1892
      implicit none
 
1893
      integer ispin,ne(2),nprj
 
1894
      real*8  sw1(ne(1)+ne(2),nprj)
 
1895
      real*8  sw2(ne(1)+ne(2),nprj)
 
1896
      integer n1dgrid,nbasis,lmax2
 
1897
      integer nindx,lm_rholm(*),iprj_rholm(*),jprj_rholm(*)
 
1898
      integer bi_rholm(*),bj_rholm(*)
 
1899
      real*8  coeff_rholm(*)
 
1900
      real*8  matr(nbasis,nbasis,lmax2,ispin)
 
1901
 
 
1902
      integer ms,n,i,lm,iprj,jprj,bi,bj,n1(2),n2(2)
 
1903
      real*8  tmp_ms,coeff,w,scal
 
1904
 
 
1905
*     **** external functions ****
 
1906
      real*8   lattice_omega
 
1907
      external lattice_omega
 
1908
 
 
1909
      call nwpw_timing_start(21)
 
1910
      n1(1) = 1 
 
1911
      n1(2) = ne(1)+1
 
1912
      n2(1) = ne(1)
 
1913
      n2(2) = ne(1)+ne(2)
 
1914
      scal = 1.0d0/dsqrt(lattice_omega())
 
1915
 
 
1916
*     ***init to zero***
 
1917
      do i=1,nindx
 
1918
         lm    = lm_rholm(i)
 
1919
         iprj  = iprj_rholm(i)
 
1920
         jprj  = jprj_rholm(i)
 
1921
         bi    = bi_rholm(i)
 
1922
         bj    = bj_rholm(i)
 
1923
         coeff = coeff_rholm(i)*scal
 
1924
         do ms=1,ispin
 
1925
            if (ne(ms).gt.0) then
 
1926
               do n=n1(ms),n2(ms)
 
1927
                  sw2(n,iprj) = sw2(n,iprj) 
 
1928
     >                        + coeff*matr(bi,bj,lm+1,ms)*sw1(n,jprj)
 
1929
               end do 
 
1930
            end if
 
1931
         end do
 
1932
      end do
 
1933
      call nwpw_timing_end(21)
 
1934
      return
 
1935
      end 
 
1936
 
 
1937