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