← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/nwpw/band/kbpp/kbpp_e_band.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwpw/band/kbpp/kbpp_e_band.F
 
1
*
 
2
* $Id: kbpp_e_band.F 22548 2012-06-02 21:10:57Z bylaska $
 
3
*
 
4
 
 
5
*     **************************************
 
6
*     *                                    *
 
7
*     *           kbpp_e_band                *
 
8
*     *                                    *
 
9
*     **************************************
 
10
 
 
11
      logical function kbpp_e_band(oprint_in,version,
 
12
     >                  psp_filename,formatted_filename,
 
13
     >                  ngrid,unita,locp,lmax,
 
14
     >                  nbrillioun,kvectors)
 
15
      implicit none
 
16
#include "errquit.fh"
 
17
#include "mafdecls.fh"
 
18
#include "tcgmsg.fh"
 
19
#include "msgtypesf.h"
 
20
#include "util.fh"
 
21
 
 
22
      logical          oprint_in
 
23
      integer          version
 
24
      character*50     psp_filename,formatted_filename
 
25
      integer          ngrid(3)
 
26
      double precision unita(3,3)
 
27
      integer locp,lmax
 
28
      integer nbrillioun
 
29
      real*8  kvectors(3,*)
 
30
 
 
31
 
 
32
*     **** local variables ****
 
33
      character*255 full_filename
 
34
 
 
35
      integer taskid,MASTER,msglen
 
36
      parameter (MASTER=0)
 
37
 
 
38
*     **** 1d pseudopotential data ****
 
39
      character*2 atom
 
40
      character*80 comment
 
41
      double precision zv,amass
 
42
      integer lmax0,lmax1,locp1
 
43
      integer n_extra,n_expansion(0:9),ihasae
 
44
      double precision rc(0:9),rlocal1
 
45
      integer nrho
 
46
      double precision drho
 
47
      integer rho_indx,vp_indx,wp_indx,sc_r_indx,sc_k_indx
 
48
      integer rho_hndl,vp_hndl,wp_hndl,sc_r_hndl,sc_k_hndl
 
49
 
 
50
      integer          isemicore
 
51
      logical          semicore
 
52
      double precision rcore,core_charge
 
53
 
 
54
      integer f_indx,cs_indx,sn_indx
 
55
      integer n_prj_indx,l_prj_indx,m_prj_indx
 
56
      integer f_hndl,cs_hndl,sn_hndl
 
57
      integer n_prj_hndl,l_prj_hndl,m_prj_hndl
 
58
 
 
59
*     ***** ngrid data *****
 
60
      integer vl_indx,vnl_indx,vnlnrm_indx,G_indx
 
61
      integer vl_hndl,vnl_hndl,vnlnrm_hndl,G_hndl
 
62
 
 
63
*     **** ray data ****
 
64
      integer nray,G_ray_hndl,tmp_ray_hndl
 
65
      integer rho_sc_k_ray_hndl,vnl_ray_hndl,vl_ray_hndl
 
66
      integer G_ray_indx,tmp_ray_indx
 
67
      integer rho_sc_k_ray_indx,vnl_ray_indx,vl_ray_indx
 
68
 
 
69
*     **** other variables ****
 
70
      logical value,mprint,hprint,oprint,filter,kray
 
71
      double precision unitg(3,3)
 
72
      integer nsize,i,l,ierr,nb,psp_type
 
73
      integer nfft1,nfft2,nfft3
 
74
      integer nmax,nprj
 
75
 
 
76
*     **** external functions ****
 
77
      logical  control_print
 
78
      external control_print
 
79
      logical  control_kbpp_ray,control_kbpp_filter
 
80
      external control_kbpp_ray,control_kbpp_filter
 
81
      double precision SIMP
 
82
      external         SIMP
 
83
      integer  kbpp_band_calc_nray
 
84
      external kbpp_band_calc_nray
 
85
 
 
86
c      call Parallel_init()
 
87
      call Parallel_taskid(taskid)
 
88
      hprint = (taskid.eq.MASTER).and.control_print(print_high)
 
89
      mprint = (taskid.eq.MASTER).and.control_print(print_medium)
 
90
      oprint = (oprint_in.or.hprint)
 
91
 
 
92
 
 
93
      value = .false.
 
94
  
 
95
*     ***** read in pseudopotential data ****
 
96
      if (taskid.eq.MASTER) then
 
97
      call util_file_name_noprefix(psp_filename,.false.,.false.,
 
98
     >                    full_filename)
 
99
      l = index(full_filename,' ') - 1
 
100
      open(unit=11,file=full_filename(1:l),
 
101
     >             status='old',form='formatted')
 
102
      read(11,*) ihasae
 
103
      read(11,*) atom
 
104
      read(11,*) zv,amass,lmax0,lmax1,locp1,rlocal1
 
105
      read(11,*) (rc(i),i=0,lmax0)
 
106
 
 
107
      read(11,*) n_extra
 
108
      read(11,*) (n_expansion(i),i=0,lmax0)
 
109
 
 
110
      read(11,*) nrho,drho
 
111
      read(11,'(A)') comment
 
112
      end if
 
113
 
 
114
      msglen = 1
 
115
      psp_type=0
 
116
      call Parallel_Brdcst_values(MASTER,msglen,zv)
 
117
      call Parallel_Brdcst_values(MASTER,msglen,amass)
 
118
      call Parallel_Brdcst_ivalues(MASTER,msglen,lmax0)
 
119
      call Parallel_Brdcst_ivalues(MASTER,msglen,lmax1)
 
120
      call Parallel_Brdcst_ivalues(MASTER,msglen,locp1)
 
121
      call Parallel_Brdcst_ivalues(MASTER,msglen,n_extra)
 
122
      msglen = lmax0+1
 
123
      call Parallel_Brdcst_values(MASTER,msglen,rc)
 
124
      call Parallel_Brdcst_ivalues(MASTER,msglen,n_expansion)
 
125
      msglen = 1
 
126
      call Parallel_Brdcst_ivalues(MASTER,msglen,nrho)
 
127
      call Parallel_Brdcst_values(MASTER,msglen,drho)
 
128
 
 
129
 
 
130
*     **** set the maximum angular momentum ****
 
131
      if (lmax.eq.-1)    lmax = lmax1
 
132
      if (lmax.gt.lmax0) lmax = lmax0
 
133
      if (lmax.lt.0)     lmax = lmax0
 
134
 
 
135
*     **** set the local potential ****
 
136
      if (locp.eq.-1)   locp = locp1
 
137
      if (locp.gt.lmax) locp = lmax
 
138
      if (locp.lt.0)    locp = lmax
 
139
 
 
140
*     **** allocate rho, vp, and wp ****
 
141
      value = MA_alloc_get(mt_dbl,nrho,
 
142
     >                        'rho',rho_hndl,rho_indx)
 
143
      value = value.and.MA_alloc_get(mt_dbl,nrho*(lmax+1),
 
144
     >                        'vp',vp_hndl, vp_indx)
 
145
      value = value.and.MA_alloc_get(mt_dbl,nrho*(lmax+1+n_extra),
 
146
     >                        'wp', wp_hndl, wp_indx)
 
147
      value = value.and.MA_alloc_get(mt_dbl,2*nrho,
 
148
     >                        'sc', sc_r_hndl, sc_r_indx)
 
149
      if (.not.value)
 
150
     > call errquit("kbpp_band: out of stack memory",0,MA_ERR)
 
151
 
 
152
      if (taskid.eq.MASTER) then
 
153
         call read_vpwp_band_extra(11,nrho,lmax,lmax0,n_extra,
 
154
     >                         dbl_mb(rho_indx),
 
155
     >                         dbl_mb(vp_indx),
 
156
     >                         dbl_mb(wp_indx))
 
157
         call read_semicore_band(11,isemicore,rcore,nrho,
 
158
     >                           dbl_mb(sc_r_indx))
 
159
         close(11)
 
160
      end if
 
161
 
 
162
      msglen = nrho
 
163
      call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(rho_indx))
 
164
      msglen = nrho*(lmax+1)
 
165
      call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(vp_indx))
 
166
      msglen = nrho*(lmax+1+n_extra)
 
167
      call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(wp_indx))
 
168
 
 
169
      msglen = 1
 
170
      call Parallel_Brdcst_ivalues(MASTER,msglen,isemicore)
 
171
      semicore = (isemicore.eq.1)
 
172
      if (semicore) then
 
173
      msglen = 2*nrho
 
174
      call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(sc_r_indx))
 
175
      else
 
176
         rcore = 0.0d0
 
177
      end if
 
178
 
 
179
 
 
180
*    **** more temporary space ****
 
181
      value = MA_alloc_get(mt_dbl,nrho,
 
182
     >                        'f',f_hndl,f_indx)
 
183
      value = MA_alloc_get(mt_dbl,nrho,
 
184
     >                        'cs',cs_hndl,cs_indx)
 
185
      value = MA_alloc_get(mt_dbl,nrho,
 
186
     >                        'sn',sn_hndl,sn_indx)
 
187
      if (.not.value)
 
188
     > call errquit("kbpp_band: out of stack memory",0,MA_ERR)
 
189
 
 
190
*     **** allocate vl,vnl,vnlnrm G ****
 
191
      nprj = 0
 
192
      nmax = 1
 
193
      do l=0,lmax
 
194
         if (l.ne.locp) then
 
195
            nprj = nprj + n_expansion(l)*(2*l+1)
 
196
            if (n_expansion(l).gt.nmax) nmax = n_expansion(l)
 
197
         end if
 
198
      end do
 
199
      !lmmax = (lmax+1)**2 - (2*locp+1)
 
200
 
 
201
      nsize = (ngrid(1))*ngrid(2)*ngrid(3)
 
202
      value = MA_alloc_get(mt_dbl,nsize,
 
203
     >                        'vl',vl_hndl,vl_indx)
 
204
      value = value.and.MA_alloc_get(mt_dbl,nsize*(nprj),
 
205
     >                        'vnl',vnl_hndl, vnl_indx)
 
206
      value = value.and.MA_alloc_get(mt_dbl,nmax*nmax*(lmax+1),
 
207
     >                        'vnlnrm', vnlnrm_hndl, vnlnrm_indx)
 
208
      value = value.and.MA_alloc_get(mt_dbl,nsize*(3),
 
209
     >                        'G',G_hndl, G_indx)
 
210
      value = value.and.MA_alloc_get(mt_dbl,4*nsize,
 
211
     >                        'sc_k',sc_k_hndl,sc_k_indx)
 
212
      value = value.and.MA_alloc_get(mt_int,nprj,
 
213
     >                        'n_prj', n_prj_hndl, n_prj_indx)
 
214
      value = value.and.MA_alloc_get(mt_int,nprj,
 
215
     >                        'l_prj', l_prj_hndl, l_prj_indx)
 
216
      value = value.and.MA_alloc_get(mt_int,nprj,
 
217
     >                        'm_prj', m_prj_hndl, m_prj_indx)
 
218
      if (.not.value)
 
219
     > call errquit("kbpp_band: out of stack memory",0,MA_ERR)
 
220
 
 
221
 
 
222
 
 
223
*     **** preparation of constants ****
 
224
      nfft1=ngrid(1)
 
225
      nfft2=ngrid(2)
 
226
      nfft3=ngrid(3)
 
227
      call setup_kbpp_band(nfft1,nfft2,nfft3,unita,unitg,dbl_mb(G_indx))
 
228
      filter = control_kbpp_filter()
 
229
      kray   = control_kbpp_ray()
 
230
     
 
231
      if (kray) then
 
232
        !**** allocate memory for rays ****
 
233
        nray = kbpp_band_calc_nray(nfft1,nfft2,nfft3,dbl_mb(G_indx))
 
234
 
 
235
        value =           MA_alloc_get(mt_dbl,nray,
 
236
     >                  'G_ray',G_ray_hndl,G_ray_indx)
 
237
        value = value.and.MA_alloc_get(mt_dbl,2*nray,
 
238
     >                  'vl_ray',vl_ray_hndl,vl_ray_indx)
 
239
        value = value.and.MA_alloc_get(mt_dbl,2*nray*(lmax+1+n_extra),
 
240
     >                  'vnl_ray',vnl_ray_hndl,vnl_ray_indx)
 
241
        value = value.and.MA_alloc_get(mt_dbl,2*nray*2,
 
242
     >              'rho_sc_k_ray',rho_sc_k_ray_hndl,rho_sc_k_ray_indx)
 
243
        value = value.and.MA_alloc_get(mt_dbl,nray,
 
244
     >                  'tmp_ray',tmp_ray_hndl,tmp_ray_indx)
 
245
        if (.not.value)
 
246
     >   call errquit('kbpp_band:out of heap memory',0,MA_ERR)
 
247
 
 
248
        call kbpp_band_generate_G_ray(nfft1,nfft2,nfft3,
 
249
     >                         dbl_mb(G_indx),
 
250
     >                         dbl_mb(G_ray_indx))
 
251
        call integrate_kbpp_e_band_local_new(version,
 
252
     >                      nrho,drho,lmax,locp,nmax,
 
253
     >                      n_extra,n_expansion,zv,
 
254
     >                                dbl_mb(vp_indx),
 
255
     >                                dbl_mb(wp_indx),
 
256
     >                                dbl_mb(rho_indx),
 
257
     >                                dbl_mb(f_indx),
 
258
     >                                dbl_mb(cs_indx),
 
259
     >                                dbl_mb(sn_indx),
 
260
     >                      nfft1,nfft2,nfft3,nprj,
 
261
     >                                dbl_mb(G_indx),
 
262
     >                                dbl_mb(vl_indx),
 
263
     >                                int_mb(n_prj_indx),
 
264
     >                                int_mb(l_prj_indx),
 
265
     >                                int_mb(m_prj_indx),
 
266
     >                                dbl_mb(vnlnrm_indx),
 
267
     >                                semicore,
 
268
     >                                dbl_mb(sc_r_indx),
 
269
     >                                dbl_mb(sc_k_indx),
 
270
     >                      nray,
 
271
     >                                dbl_mb(G_ray_indx),
 
272
     >                                dbl_mb(vl_ray_indx),
 
273
     >                                dbl_mb(vnl_ray_indx),
 
274
     >                                dbl_mb(rho_sc_k_ray_indx),
 
275
     >                                dbl_mb(tmp_ray_indx),
 
276
     >                                filter,
 
277
     >                      ierr)
 
278
      else
 
279
      call integrate_kbpp_e_band_local(version,
 
280
     >                      nrho,drho,lmax,locp,nmax,
 
281
     >                      n_extra,n_expansion,zv,
 
282
     >                                dbl_mb(vp_indx),
 
283
     >                                dbl_mb(wp_indx),
 
284
     >                                dbl_mb(rho_indx),
 
285
     >                                dbl_mb(f_indx),
 
286
     >                                dbl_mb(cs_indx),
 
287
     >                                dbl_mb(sn_indx),
 
288
     >                      nfft1,nfft2,nfft3,nprj,
 
289
     >                                dbl_mb(G_indx),
 
290
     >                                dbl_mb(vl_indx),
 
291
     >                                int_mb(n_prj_indx),
 
292
     >                                int_mb(l_prj_indx),
 
293
     >                                int_mb(m_prj_indx),
 
294
     >                                dbl_mb(vnlnrm_indx),
 
295
     >                                semicore,
 
296
     >                                dbl_mb(sc_r_indx),
 
297
     >                                dbl_mb(sc_k_indx),
 
298
     >                      ierr)
 
299
      end if
 
300
 
 
301
      if ((taskid.eq.MASTER).and.(oprint)) then
 
302
      write(*,*) "     ********************************************"
 
303
      write(*,*) "     *                                          *"
 
304
      write(*,*) "     *   KBPP_BAND - Pseudopotential Formatter  *"
 
305
      write(*,*) "     *                                          *"
 
306
      write(*,*) "     *        version last updated 01/01/02     *"
 
307
      write(*,*) "     *                                          *"
 
308
      write(*,*) "     *        developed by Eric J. Bylaska      *"
 
309
      write(*,*) "     *                                          *"
 
310
      write(*,*) "     ********************************************"
 
311
      call nwpw_message(1)
 
312
      write(*,*)
 
313
      write(*,*) "Pseudpotential Data"
 
314
      write(*,*) "-------------------"
 
315
      write(*,*) "  atom     :",atom
 
316
      write(*,*) "  charge   :",zv
 
317
      write(*,*) "  mass no. :",amass
 
318
      write(*,*) "  highest angular component      :",lmax0
 
319
      write(*,*) "  highest angular component used :",lmax
 
320
      write(*,*) "  local potential used           :",locp
 
321
      write(*,*) "  number of projectors           :",nprj
 
322
      write(*,111) "  cutoffs: ",(rc(i), i=0,lmax)
 
323
      if (semicore) then
 
324
        write(*,*)
 
325
        write(*,115) "  semi-core charge included, rcore:",rcore
 
326
        do i=1,nrho
 
327
           dbl_mb(f_indx+i-1) = dbl_mb(sc_r_indx+i-1)
 
328
     >                        * dbl_mb(rho_indx+i-1)**2
 
329
        end do
 
330
        core_charge=16.0d0*datan(1.0d0)*SIMP(nrho,dbl_mb(f_indx),drho)
 
331
        write(*,115) "  semi-core charge                :",core_charge,
 
332
     >                dbl_mb(sc_k_indx)
 
333
         do i=1,nrho
 
334
           dbl_mb(f_indx+i-1) = dbl_mb(sc_r_indx+i-1+nrho)
 
335
     >                        * dbl_mb(rho_indx+i-1)**2
 
336
         end do
 
337
         core_charge=16.0d0*datan(1.0d0)*SIMP(nrho,dbl_mb(f_indx),drho)
 
338
         write(*,115) "  Semi-core charge gradient       :",
 
339
     >                core_charge
 
340
 
 
341
      end if
 
342
      write(*,*)
 
343
      write(*,*) "Simulation Cell"
 
344
      write(*,*) "---------------"
 
345
      if (version.eq.3) write(*,112) "  boundry: periodic"
 
346
      write(*,113) "  ngrid  :",ngrid
 
347
      write(*,114) "  unita  :",unita(1,1),unita(2,1),unita(3,1)
 
348
      write(*,114) "          ",unita(1,2),unita(2,2),unita(3,2)
 
349
      write(*,114) "          ",unita(1,3),unita(2,3),unita(3,3)
 
350
      write(*,*)
 
351
  111 format(a,10f10.3)
 
352
  112 format(a)
 
353
  113 format(a,3I4)
 
354
  114 format(a,3F10.3)
 
355
  115 format(a,2E14.6)
 
356
      end if
 
357
 
 
358
 
 
359
      if (taskid.eq.MASTER) then
 
360
        call util_file_name_noprefix(formatted_filename,
 
361
     >                    .false.,
 
362
     >                    .false.,
 
363
     >                    full_filename)
 
364
          l = index(full_filename,' ') - 1
 
365
         if (mprint) then
 
366
         write(*,*)
 
367
         write(*,*) "Generated formatted_filename: ",full_filename(1:l)
 
368
         if (kray)write(*,'(A,I5,A)')" - Spline fitted, nray=",nray," -"
 
369
         if (filter) write(*,'(A)') " - filtered -"
 
370
         end if
 
371
         call openfile(2,full_filename,l,'w',l)     
 
372
         call cwrite(2,comment,80)
 
373
         call iwrite(2,psp_type,1)
 
374
         call iwrite(2,version,1)
 
375
         call iwrite(2,ngrid,3)
 
376
         call dwrite(2,unita,9)
 
377
         call cwrite(2,atom,2)
 
378
         call dwrite(2,amass,1)
 
379
         call dwrite(2,zv,1)
 
380
         call iwrite(2,lmax,1)
 
381
         call iwrite(2,locp,1)
 
382
 
 
383
         call iwrite(2,nmax,1)
 
384
         call dwrite(2,rc,lmax+1)
 
385
 
 
386
         call iwrite(2,nprj,1)
 
387
         if (nprj.gt.0) then
 
388
         call iwrite(2,int_mb(n_prj_indx),nprj)
 
389
         call iwrite(2,int_mb(l_prj_indx),nprj)
 
390
         call iwrite(2,int_mb(m_prj_indx),nprj)
 
391
         call dwrite(2,dbl_mb(vnlnrm_indx),nmax*nmax*(lmax+1))
 
392
         end if
 
393
         
 
394
         call dwrite(2,rcore,1)
 
395
         call iwrite(2,nbrillioun,1)
 
396
         call dwrite(2,kvectors,3*nbrillioun)
 
397
         call dwrite(2,dbl_mb(vl_indx),nsize)
 
398
      end if
 
399
 
 
400
      do nb=1,nbrillioun
 
401
 
 
402
         if ((oprint).and.(taskid.eq.MASTER))
 
403
     >      write(*,*) "generating brillioun #",nb
 
404
 
 
405
            if (kray) then 
 
406
            call integrate_kbpp_e_band_nonlocal_new(version,
 
407
     >                      kvectors(1,nb),
 
408
     >                      nrho,drho,lmax,locp,nmax,
 
409
     >                      n_extra,n_expansion,zv,
 
410
     >                                dbl_mb(vp_indx),
 
411
     >                                dbl_mb(wp_indx),
 
412
     >                                dbl_mb(rho_indx),
 
413
     >                                dbl_mb(f_indx),
 
414
     >                                dbl_mb(cs_indx),
 
415
     >                                dbl_mb(sn_indx),
 
416
     >                      nfft1,nfft2,nfft3,nprj,
 
417
     >                                dbl_mb(G_indx),
 
418
     >                                dbl_mb(vnl_indx),
 
419
     >                      nray,
 
420
     >                                dbl_mb(G_ray_indx),
 
421
     >                                dbl_mb(vnl_ray_indx),
 
422
     >                      ierr)
 
423
            else
 
424
            call integrate_kbpp_e_band_nonlocal(version,
 
425
     >                      kvectors(1,nb),
 
426
     >                      nrho,drho,lmax,locp,nmax,
 
427
     >                      n_extra,n_expansion,zv,
 
428
     >                                dbl_mb(vp_indx),
 
429
     >                                dbl_mb(wp_indx),
 
430
     >                                dbl_mb(rho_indx),
 
431
     >                                dbl_mb(f_indx),
 
432
     >                                dbl_mb(cs_indx),
 
433
     >                                dbl_mb(sn_indx),
 
434
     >                      nfft1,nfft2,nfft3,nprj,
 
435
     >                                dbl_mb(G_indx),
 
436
     >                                dbl_mb(vnl_indx),
 
437
     >                      ierr)
 
438
            end if
 
439
 
 
440
         if (taskid.eq.MASTER) 
 
441
     >     call dwrite(2,dbl_mb(vnl_indx),nsize*nprj)
 
442
      end do
 
443
 
 
444
 
 
445
 
 
446
      if (taskid.eq.MASTER) then
 
447
         if (semicore) then
 
448
           call dwrite(2,dbl_mb(sc_k_indx),4*nsize)
 
449
         end if
 
450
      call closefile(2)
 
451
      end if
 
452
     
 
453
*     **** free ray heap space ****
 
454
      if (kray) then
 
455
        value = MA_free_heap(tmp_ray_hndl)
 
456
        value = value.and.MA_free_heap(rho_sc_k_ray_hndl)
 
457
        value = value.and.MA_free_heap(vnl_ray_hndl)
 
458
        value = value.and.MA_free_heap(vl_ray_hndl)
 
459
        value = value.and.MA_free_heap(G_ray_hndl)
 
460
        if (.not.value)
 
461
     >   call errquit('kbpp_band:Error freeing memory',0,MA_ERR)
 
462
      end if
 
463
 
 
464
*     **** free heap space ****
 
465
      value = MA_free_heap(rho_hndl)
 
466
      value = value.and.MA_free_heap(vp_hndl)
 
467
      value = value.and.MA_free_heap(wp_hndl)
 
468
      value = value.and.MA_free_heap(sc_r_hndl)
 
469
      value = value.and.MA_free_heap(sc_k_hndl)
 
470
      value = value.and.MA_free_heap(f_hndl)
 
471
      value = value.and.MA_free_heap(cs_hndl)
 
472
      value = value.and.MA_free_heap(sn_hndl)
 
473
 
 
474
      value = value.and.MA_free_heap(vl_hndl)
 
475
      value = value.and.MA_free_heap(vnl_hndl)
 
476
      value = value.and.MA_free_heap(vnlnrm_hndl)
 
477
      value = value.and.MA_free_heap(G_hndl)
 
478
      value = value.and.MA_free_heap(n_prj_hndl)
 
479
      value = value.and.MA_free_heap(l_prj_hndl)
 
480
      value = value.and.MA_free_heap(m_prj_hndl)
 
481
      if (.not.value)
 
482
     > call errquit("kbpp_band: error freeing memory",0,MA_ERR)
 
483
 
 
484
c      call Parallel_Finalize()
 
485
      
 
486
      if ((taskid.eq.MASTER).and.(oprint)) call nwpw_message(4)
 
487
      kbpp_e_band = value
 
488
      return
 
489
 
 
490
 9999 call errquit('kbpp_e_band:Error reading psp_filename',0,DISK_ERR)
 
491
      kbpp_e_band = value
 
492
      END
 
493
 
 
494
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
495
      subroutine read_vpwp_band_extra(unit,nrho,lmax,lmax0,n_extra,
 
496
     >                                rho,vp,wp)
 
497
      implicit none
 
498
      integer unit
 
499
      integer nrho,lmax,lmax0,n_extra
 
500
      double precision rho(nrho)
 
501
      double precision vp(nrho,0:lmax)
 
502
      double precision wp(nrho,0:(lmax+n_extra))
 
503
 
 
504
#include "errquit.fh"
 
505
 
 
506
      integer i,j
 
507
      real*8 tmp(0:20)
 
508
 
 
509
      do i=1,nrho
 
510
        read(unit,*,ERR=9999,END=9999) rho(i),(vp(i,j),j=0,lmax)
 
511
      end do
 
512
      do i=1,nrho
 
513
        read(unit,*,ERR=9999,END=9999) rho(i),(tmp(j),j=0,lmax0+n_extra)
 
514
        do j=0,lmax
 
515
           wp(i,j) = tmp(j)
 
516
        end do
 
517
        do j=1,n_extra
 
518
           wp(i,j+lmax) = tmp(j+lmax0)
 
519
        end do
 
520
      end do
 
521
      return
 
522
 9999 call errquit('Error reading psp_filename',0, DISK_ERR)
 
523
      end
 
524
 
 
525
 
 
526