1*
2* $Id$
3*
4
5
6*     **************************************
7*     *                                    *
8*     *           kbpp_band_stress         *
9*     *                                    *
10*     **************************************
11
12      logical function kbpp_band_stress(oprint_in,version,
13     >                  psp_filename,formatted_filename,
14     >                  ngrid,unita,locp,lmax,
15     >                  nbrillioun,kvectors)
16      implicit none
17#include "errquit.fh"
18#include "bafdecls.fh"
19#include "tcgmsg.fh"
20#include "msgtypesf.h"
21#include "util.fh"
22
23      logical          oprint_in
24      integer          version
25      character*50     psp_filename,formatted_filename
26      integer          ngrid(3)
27      double precision unita(3,3)
28      integer locp,lmax
29      integer nbrillioun
30      real*8  kvectors(3,*)
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 nb,lmax0,lmmax,lmax1,locp1
43      double precision rc(0:9),rlocal1
44      integer nrho
45      double precision drho
46      integer rho_indx,vp_indx,wp_indx,sc_r_indx,sc_k_indx
47      integer rho_hndl,vp_hndl,wp_hndl,sc_r_hndl,sc_k_hndl
48      integer          isemicore
49      logical          semicore
50      double precision rcore
51
52      integer f_indx,cs_indx,sn_indx
53      integer f_hndl,cs_hndl,sn_hndl
54
55*     ***** ngrid data *****
56      integer vl_indx,vnl_indx,G_indx
57      integer vl_hndl,vnl_hndl,G_hndl
58
59*     **** ray data ****
60      integer nray,G_ray_hndl,tmp_ray_hndl
61      integer rho_sc_k_ray_hndl,dvnl_ray_hndl,dvl_ray_hndl
62      integer G_ray_indx,tmp_ray_indx
63      integer rho_sc_k_ray_indx,dvnl_ray_indx,dvl_ray_indx
64
65*     **** other variables ****
66      logical value,mprint,hprint,oprint,filter,kray
67      double precision unitg(3,3)
68      integer nsize,i,l,ierr,psp_type
69      integer nfft1,nfft2,nfft3
70
71*     **** external functions ****
72      logical  control_print,control_kbpp_ray,control_kbpp_filter
73      external control_print,control_kbpp_ray,control_kbpp_filter
74      integer  kbpp_band_calc_nray
75      external kbpp_band_calc_nray
76
77
78c      call Parallel_init()
79      call Parallel_taskid(taskid)
80      hprint = (taskid.eq.MASTER).and.control_print(print_high)
81      mprint = (taskid.eq.MASTER).and.control_print(print_medium)
82      oprint = (oprint_in.or.hprint)
83
84
85      value = .false.
86
87*     ***** read in pseudopotential data ****
88
89      if (taskid.eq.MASTER) then
90      call util_file_name_noprefix(psp_filename,.false.,.false.,
91     >                    full_filename)
92      l = index(full_filename,' ') - 1
93      open(unit=11,file=full_filename(1:l),
94     >             status='old',form='formatted')
95
96      read(11,'(A2)') atom
97      read(11,*) zv,amass,lmax0,lmax1,locp1,rlocal1
98      read(11,*) (rc(i),i=0,lmax0)
99      read(11,*) nrho,drho
100      read(11,'(A)') comment
101      end if
102
103      msglen = 1
104      call BRDCST(9+MSGDBL,zv,mdtob(msglen),MASTER)
105      call BRDCST(9+MSGDBL,amass,mdtob(msglen),MASTER)
106      call BRDCST(9+MSGINT,lmax0,mitob(msglen),MASTER)
107      call BRDCST(9+MSGINT,lmax1,mitob(msglen),MASTER)
108      call BRDCST(9+MSGINT,locp1,mitob(msglen),MASTER)
109      msglen = lmax0+1
110      call BRDCST(9+MSGDBL,rc,mdtob(msglen),MASTER)
111      msglen = 1
112      call BRDCST(9+MSGINT,nrho,mitob(msglen),MASTER)
113      call BRDCST(9+MSGDBL,drho,mdtob(msglen),MASTER)
114
115
116*     **** set the maximum angular momentum ****
117      if (lmax.eq.-1)    lmax = lmax1
118      if (lmax.gt.lmax0) lmax = lmax0
119      if (lmax.lt.0)     lmax = lmax0
120
121*     **** set the local potential ****
122      if (locp.eq.-1)   locp = locp1
123      if (locp.gt.lmax) locp = lmax
124      if (locp.lt.0)    locp = lmax
125
126*     **** allocate rho, vp, and wp ****
127      value = BA_alloc_get(mt_dbl,nrho,
128     >                        'rho',rho_hndl,rho_indx)
129      value = value.and.BA_alloc_get(mt_dbl,nrho*(lmax+1),
130     >                        'vp',vp_hndl, vp_indx)
131      value = value.and.BA_alloc_get(mt_dbl,nrho*(lmax+1),
132     >                        'wp', wp_hndl, wp_indx)
133      value = value.and.BA_alloc_get(mt_dbl,2*nrho,
134     >                        'sc', sc_r_hndl, sc_r_indx)
135      if (.not.value)
136     > call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR)
137
138      if (taskid.eq.MASTER) then
139      call read_vpwp_band(11,nrho,lmax,dbl_mb(rho_indx),
140     >                         dbl_mb(vp_indx),
141     >                         dbl_mb(wp_indx))
142      call read_semicore_band(11,isemicore,rcore,nrho,dbl_mb(sc_r_indx))
143      close(11)
144      end if
145
146      msglen = nrho
147      call BRDCST(9+MSGDBL,dbl_mb(rho_indx),mdtob(msglen),MASTER)
148      msglen = nrho*(lmax+1)
149      call BRDCST(9+MSGDBL,dbl_mb(vp_indx),mdtob(msglen),MASTER)
150      call BRDCST(9+MSGDBL,dbl_mb(wp_indx),mdtob(msglen),MASTER)
151
152      msglen = 1
153      call BRDCST(9+MSGINT,isemicore,mitob(msglen),MASTER)
154      semicore = (isemicore.eq.1)
155      if (semicore) then
156      msglen = 2*nrho
157      call BRDCST(9+MSGDBL,dbl_mb(sc_r_indx),mdtob(msglen),MASTER)
158      else
159         rcore = 0.0d0
160      end if
161
162
163*    **** more temporary space ****
164      value = BA_alloc_get(mt_dbl,nrho,
165     >                        'f',f_hndl,f_indx)
166      value = value.and.BA_alloc_get(mt_dbl,nrho,
167     >                        'cs',cs_hndl,cs_indx)
168      value = value.and.BA_alloc_get(mt_dbl,nrho,
169     >                        'sn',sn_hndl,sn_indx)
170      if (.not.value)
171     > call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR)
172
173*     **** allocate vl,vnl,vnlnrm G ****
174      lmmax = (lmax+1)**2 - (2*locp+1)
175      nsize = ngrid(1)*ngrid(2)*ngrid(3)
176      value = BA_alloc_get(mt_dbl,nsize,
177     >                        'vl',vl_hndl,vl_indx)
178      value = value.and.BA_alloc_get(mt_dbl,nsize*(lmmax)*3,
179     >                        'vnl',vnl_hndl, vnl_indx)
180      value = value.and.BA_alloc_get(mt_dbl,nsize*(3),
181     >                        'G',G_hndl, G_indx)
182      value = value.and.BA_alloc_get(mt_dbl,4*nsize,
183     >                        'sc_k',sc_k_hndl,sc_k_indx)
184      if (.not.value)
185     > call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR)
186
187
188*     **** preparation of constants ****
189      nfft1=ngrid(1)
190      nfft2=ngrid(2)
191      nfft3=ngrid(3)
192      call setup_kbpp_band(nfft1,nfft2,nfft3,unita,unitg,dbl_mb(G_indx))
193      filter = control_kbpp_filter()
194      kray   = control_kbpp_ray()
195
196      if (kray) then
197        !**** allocate memory for rays ****
198        nray = kbpp_band_calc_nray(nfft1,nfft2,nfft3,dbl_mb(G_indx))
199
200        value =           BA_alloc_get(mt_dbl,nray,
201     >                  'G_ray2',G_ray_hndl,G_ray_indx)
202        value = value.and.BA_alloc_get(mt_dbl,2*nray,
203     >                  'dvl_ray2',dvl_ray_hndl,dvl_ray_indx)
204        value = value.and.BA_alloc_get(mt_dbl,4*nray*(lmax+1),
205     >                  'dvnl_ray2',dvnl_ray_hndl,dvnl_ray_indx)
206        value = value.and.BA_alloc_get(mt_dbl,2*nray,
207     >              'rho_sc_k_ray2',rho_sc_k_ray_hndl,rho_sc_k_ray_indx)
208        value = value.and.BA_alloc_get(mt_dbl,nray,
209     >                  'tmp_ray2',tmp_ray_hndl,tmp_ray_indx)
210        if (.not.value)
211     >   call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR)
212
213        call kbpp_band_generate_G_ray(nfft1,nfft2,nfft3,
214     >                         dbl_mb(G_indx),
215     >                         dbl_mb(G_ray_indx))
216
217        call integrate_kbpp_band_stress_local_new(version,
218     >                      nrho,drho,lmax,locp,zv,
219     >                                dbl_mb(vp_indx),
220     >                                dbl_mb(wp_indx),
221     >                                dbl_mb(rho_indx),
222     >                                dbl_mb(f_indx),
223     >                                dbl_mb(cs_indx),
224     >                                dbl_mb(sn_indx),
225     >                      nfft1,nfft2,nfft3,lmmax,
226     >                                dbl_mb(G_indx),
227     >                                dbl_mb(vl_indx),
228     >                                semicore,
229     >                                dbl_mb(sc_r_indx),
230     >                                dbl_mb(sc_k_indx),
231     >                      nray,
232     >                                dbl_mb(G_ray_indx),
233     >                                dbl_mb(dvl_ray_indx),
234     >                                dbl_mb(dvnl_ray_indx),
235     >                                dbl_mb(rho_sc_k_ray_indx),
236     >                                dbl_mb(tmp_ray_indx),
237     >                                filter,
238     >                      ierr)
239
240
241      else
242      call integrate_kbpp_band_stress_local(version,
243     >                      nrho,drho,lmax,locp,zv,
244     >                                dbl_mb(vp_indx),
245     >                                dbl_mb(wp_indx),
246     >                                dbl_mb(rho_indx),
247     >                                dbl_mb(f_indx),
248     >                                dbl_mb(cs_indx),
249     >                                dbl_mb(sn_indx),
250     >                      nfft1,nfft2,nfft3,lmmax,
251     >                                dbl_mb(G_indx),
252     >                                dbl_mb(vl_indx),
253     >                                semicore,
254     >                                dbl_mb(sc_r_indx),
255     >                                dbl_mb(sc_k_indx),
256     >                      ierr)
257      end if
258
259
260      if ((taskid.eq.MASTER).and.(oprint)) then
261      write(*,*) "     ************************************************"
262      write(*,*) "     *                                              *"
263      write(*,*) "     * KBPP_Band_Stress - Pseudopotential Formatter *"
264      write(*,*) "     *                                              *"
265      write(*,*) "     *      version last updated 4/15/99            *"
266      write(*,*) "     *                                              *"
267      write(*,*) "     ************************************************"
268      call nwpw_message(1)
269      write(*,*)
270      write(*,*) "Pseudpotential Data"
271      write(*,*) "-------------------"
272      write(*,*) "  atom     :",atom
273      write(*,*) "  charge   :",zv
274      write(*,*) "  mass no. :",amass
275      write(*,*) "  highest angular component      :",lmax0
276      write(*,*) "  highest angular component used :",lmax
277      write(*,*) "  local potential used           :",locp
278
279      if (semicore) then
280        write(*,*)
281        write(*,115) "  semi-core stress derivative added"
282      end if
283      write(*,*)
284      write(*,*) "Simulation Cell"
285      write(*,*) "---------------"
286      if (version.eq.3) write(*,112) "  boundry: periodic"
287      write(*,113) "  ngrid  :",ngrid
288      write(*,114) "  unita  :",unita(1,1),unita(2,1),unita(3,1)
289      write(*,114) "          ",unita(1,2),unita(2,2),unita(3,2)
290      write(*,114) "          ",unita(1,3),unita(2,3),unita(3,3)
291      write(*,*)
292  111 format(a,10f10.3)
293  112 format(a)
294  113 format(a,3I4)
295  114 format(a,3F10.3)
296  115 format(a,2E14.6)
297      end if
298
299
300
301      if (taskid.eq.MASTER) then
302        call util_file_name_noprefix(formatted_filename,
303     >                    .false.,
304     >                    .false.,
305     >                    full_filename)
306         l = index(full_filename,' ') - 1
307         if (mprint) then
308         write(*,*)
309         write(*,*) "Generated formatted_stress_filename: ",
310     >               full_filename(1:l)
311         if (kray)write(*,'(A,I5,A)')" - Spline fitted, nray=",nray," -"
312         if (filter) write(*,'(A)') " - filtered -"
313         end if
314         call openfile(2,full_filename,l,'w',l)
315         call iwrite(2,version,1)
316         call iwrite(2,ngrid,3)
317         call dwrite(2,unita,9)
318
319         call iwrite(2,nbrillioun,1)
320         call dwrite(2,kvectors,3*nbrillioun)
321
322         call dwrite(2,dbl_mb(vl_indx),nsize)
323      end if
324
325
326      do nb=1,nbrillioun
327
328         if ((oprint).and.(taskid.eq.MASTER))
329     >      write(*,*) "generating brillioun #",nb,
330     >                  kvectors(1,nb),
331     >                  kvectors(2,nb),
332     >                  kvectors(3,nb)
333
334            if (kray) then
335            call integrate_kbpp_band_stress_nonlocal_new(version,
336     >                      kvectors(1,nb),
337     >                      nrho,drho,lmax,locp,zv,
338     >                                dbl_mb(vp_indx),
339     >                                dbl_mb(wp_indx),
340     >                                dbl_mb(rho_indx),
341     >                                dbl_mb(f_indx),
342     >                                dbl_mb(cs_indx),
343     >                                dbl_mb(sn_indx),
344     >                      nfft1,nfft2,nfft3,lmmax,
345     >                                dbl_mb(G_indx),
346     >                                dbl_mb(vnl_indx),
347     >                      nray,
348     >                                dbl_mb(G_ray_indx),
349     >                                dbl_mb(dvnl_ray_indx),
350     >                      ierr)
351
352            else
353            call integrate_kbpp_band_stress_nonlocal(version,
354     >                      kvectors(1,nb),
355     >                      nrho,drho,lmax,locp,zv,
356     >                                dbl_mb(vp_indx),
357     >                                dbl_mb(wp_indx),
358     >                                dbl_mb(rho_indx),
359     >                                dbl_mb(f_indx),
360     >                                dbl_mb(cs_indx),
361     >                                dbl_mb(sn_indx),
362     >                      nfft1,nfft2,nfft3,lmmax,
363     >                                dbl_mb(G_indx),
364     >                                dbl_mb(vnl_indx),
365     >                      ierr)
366
367         end if
368         if (taskid.eq.MASTER)
369     >     call dwrite(2,dbl_mb(vnl_indx),nsize*lmmax*3)
370      end do
371
372
373      if (taskid.eq.MASTER) then
374         if (semicore) then
375           call dwrite(2,dbl_mb(sc_k_indx),4*nsize)
376         end if
377      call closefile(2)
378      end if
379
380*     **** free ray heap space ****
381      if (kray) then
382        value = BA_free_heap(tmp_ray_hndl)
383        value = value.and.BA_free_heap(rho_sc_k_ray_hndl)
384        value = value.and.BA_free_heap(dvnl_ray_hndl)
385        value = value.and.BA_free_heap(dvl_ray_hndl)
386        value = value.and.BA_free_heap(G_ray_hndl)
387        if (.not.value)
388     >   call errquit('kbpp_band_stress:Error freeing memory',0,MA_ERR)
389      end if
390
391*     **** free heap space ****
392      value = BA_free_heap(rho_hndl)
393      value = value.and.BA_free_heap(vp_hndl)
394      value = value.and.BA_free_heap(wp_hndl)
395      value = value.and.BA_free_heap(sc_r_hndl)
396      value = value.and.BA_free_heap(sc_k_hndl)
397      value = value.and.BA_free_heap(f_hndl)
398      value = value.and.BA_free_heap(cs_hndl)
399      value = value.and.BA_free_heap(sn_hndl)
400
401      value = value.and.BA_free_heap(vl_hndl)
402      value = value.and.BA_free_heap(vnl_hndl)
403      value = value.and.BA_free_heap(G_hndl)
404      if (.not.value)
405     > call errquit('kbpp_band_stress:Error freeing heap',0,MA_ERR)
406
407c      call Parallel_Finalize()
408
409      if ((taskid.eq.MASTER).and.(oprint)) call nwpw_message(4)
410
411      kbpp_band_stress = value
412      return
413
414 9999 kbpp_band_stress = .false.
415      call errquit('kbpp_band_stress:Error reading psp_filename',
416     > 0,DISK_ERR)
417      END
418
419
420
421
422