1*
2* $Id$
3*
4
5*     **************************************
6*     *                                    *
7*     *           hghppv1                  *
8*     *                                    *
9*     **************************************
10
11      logical function hghppv1(oprint_in,version,
12     >                  psp_filename,formatted_filename,
13     >                  ngrid,unita,rlocal)
14      implicit none
15
16#include "bafdecls.fh"
17#include "tcgmsg.fh"
18#include "msgtypesf.h"
19#include "errquit.fh"
20#include "util.fh"
21#include "stdio.fh"
22
23      logical          oprint_in,filter
24      integer          version
25      character*50     psp_filename,formatted_filename
26      integer          ngrid(3)
27      double precision unita(3,3)
28      real*8  rlocal
29      character*255 full_filename
30
31      logical value,mprint,hprint,oprint
32
33      integer taskid,MASTER,msglen
34      parameter (MASTER=0)
35
36*     **** 1d pseudopotential data ****
37      character*80 comment
38
39      character*2 atom
40      integer Zion
41      double precision zv,amass
42
43      double precision rloc,C1,C2,C3,C4
44      double precision r(0:3),H(3,0:3),K(3,0:3)
45      double precision Gijl(9*4),Gtmp(3,3,0:3)
46      double precision rcore
47      double precision ecut,wcut
48
49      integer lmax,locp
50
51      integer n_prj_indx,l_prj_indx,m_prj_indx,b_prj_indx
52      integer n_prj_hndl,l_prj_hndl,m_prj_hndl,b_prj_hndl
53
54*     ***** ngrid data *****
55      integer Ylm_indx,vl_indx,vnl_indx,G_indx
56      integer Ylm_hndl,vl_hndl,vnl_hndl,G_hndl
57
58*     **** other variables ****
59      double precision unitg(3,3)
60      integer nsize,i,j,n,l,m,psp_type,indx
61      integer nfft1,nfft2,nfft3
62      integer nmax,nprj
63
64*     **** external functions ****
65      logical  control_print,control_kbpp_filter
66      external control_print,control_kbpp_filter
67      real*8   control_ecut,control_wcut
68      external control_ecut,control_wcut
69
70
71
72
73c      call Parallel_init()
74      call Parallel_taskid(taskid)
75      hprint = (taskid.eq.MASTER).and.control_print(print_high)
76      mprint = (taskid.eq.MASTER).and.control_print(print_medium)
77      oprint = (oprint_in.or.hprint)
78
79      value = .false.
80
81*     ***** read in pseudopotential data ****
82      if (taskid.eq.MASTER) then
83         call util_file_name_noprefix(psp_filename,.false.,.false.,
84     >                       full_filename)
85         l = index(full_filename,' ') - 1
86         open(unit=11,file=full_filename(1:l),
87     >                status='old',form='formatted')
88         read(11,*) psp_type
89         read(11,'(A2)') atom
90         read(11,*) Zion
91         read(11,*) lmax
92         read(11,*) rloc,C1,C2,C3,C4
93
94         if (lmax.ge.0) then
95            read(11,*) r(0),H(1,0),H(2,0),H(3,0)
96            do l=1,lmax
97              read(11,*) r(l),H(1,l),H(2,l),H(3,l)
98              read(11,*)      K(1,l),K(2,l),K(3,l)
99            end do
100         end if
101         read(11,'(A)') comment
102         close(11)
103
104
105         !**** determine nmax ****
106         nmax = 0
107         do l=0,lmax
108            do i=1,3
109              if ((H(i,l).ne.0.0d0) .and. (i.gt.nmax)) nmax = i
110            end do
111         end do
112      end if
113
114
115      msglen = 1
116      call BRDCST(9+MSGINT,psp_type,mitob(msglen),MASTER)
117      call BRDCST(9+MSGINT,Zion,mitob(msglen),MASTER)
118      call BRDCST(9+MSGINT,lmax,mitob(msglen),MASTER)
119      call BRDCST(9+MSGINT,nmax,mitob(msglen),MASTER)
120
121      call BRDCST(9+MSGDBL,rloc,mdtob(msglen),MASTER)
122      call BRDCST(9+MSGDBL,C1,mdtob(msglen),MASTER)
123      call BRDCST(9+MSGDBL,C2,mdtob(msglen),MASTER)
124      call BRDCST(9+MSGDBL,C3,mdtob(msglen),MASTER)
125      call BRDCST(9+MSGDBL,C4,mdtob(msglen),MASTER)
126
127      msglen = 4
128      call BRDCST(9+MSGDBL,r,mdtob(msglen),MASTER)
129
130      msglen = 12
131      call BRDCST(9+MSGDBL,H,mdtob(msglen),MASTER)
132      call BRDCST(9+MSGDBL,K,mdtob(msglen),MASTER)
133
134
135
136*     **** set the maximum angular momentum ****
137
138*     **** set the local potential ****
139      locp = lmax+1
140
141*     **** set the local potential ****
142      rlocal = 1.0d0
143
144      !**** determine nprj ****
145      nprj= 0
146      do l=0,lmax
147            !write(luout,*) "???H :",l,(H(i,l),i=1,3)
148      do i=1,3
149        if ((H(i,l).ne.0.0d0)) nprj = nprj + (2*l+1)
150      end do
151      end do
152      !write(luout,*) "???nprj:", nprj
153
154
155
156*     **** set the projector coeficients ****
157      call dcopy(9*4,0.0d0,0,Gtmp,1)
158      call dcopy(9*4,0.0d0,0,Gijl,1)
159      Gtmp(1,1,0) = H(1,0)
160      Gtmp(2,2,0) = H(2,0)
161      Gtmp(3,3,0) = H(3,0)
162      Gtmp(1,2,0) = -0.5d0*dsqrt(3.0d0/5.0d0)   *Gtmp(2,2,0)
163      Gtmp(2,1,0) = Gtmp(1,2,0)
164      Gtmp(1,3,0) =  0.5d0*dsqrt(5.0d0/21.0d0)  *Gtmp(3,3,0)
165      Gtmp(3,1,0) = Gtmp(1,3,0)
166      Gtmp(2,3,0) = -0.5d0*dsqrt(100.0d0/63.0d0)*Gtmp(3,3,0)
167      Gtmp(3,2,0) = Gtmp(2,3,0)
168
169
170      Gtmp(1,1,1) = H(1,1)
171      Gtmp(2,2,1) = H(2,1)
172      Gtmp(3,3,1) = H(3,1)
173      Gtmp(1,2,1) = -0.5d0*dsqrt(5.0d0/7.0d0)   *Gtmp(2,2,1)
174      Gtmp(2,1,1) = Gtmp(1,2,1)
175      Gtmp(1,3,1) =  (1.0d0/6.0d0)*dsqrt(35.0d0/11.0d0)*Gtmp(3,3,1)
176      Gtmp(3,1,1) = Gtmp(1,3,1)
177      Gtmp(2,3,1) = -(14.0d0/6.0d0)*dsqrt(1.0d0/11.0d0)*Gtmp(3,3,1)
178      Gtmp(3,2,1) = Gtmp(2,3,1)
179
180      Gtmp(1,1,2) = H(1,2)
181      Gtmp(2,2,2) = H(2,2)
182      Gtmp(3,3,2) = H(3,2)
183      Gtmp(1,2,2) = -0.5d0*dsqrt(7.0d0/9.0d0)   *Gtmp(2,2,2)
184      Gtmp(2,1,2) = Gtmp(1,2,2)
185      Gtmp(1,3,2) =  0.5d0*dsqrt(63.0d0/143.0d0)*Gtmp(3,3,2)
186      Gtmp(3,1,2) = Gtmp(1,3,2)
187      Gtmp(2,3,2) = -0.5d0*18.0d0*dsqrt(1.0d0/143.0d0)*Gtmp(3,3,2)
188      Gtmp(3,2,2) = Gtmp(2,3,2)
189
190      Gtmp(1,1,3) = H(1,3)
191      Gtmp(2,2,3) = H(2,3)
192      Gtmp(3,3,3) = H(3,3)
193
194
195
196      do l=0,lmax
197        do i=1,nmax
198        do j=1,nmax
199           Gijl(i+(j-1)*nmax+l*nmax*nmax) = Gtmp(i,j,l)
200        end do
201        end do
202      end do
203
204
205*     **** allocate vl, vnl, G ****
206      nsize = (ngrid(1)/2+1)*ngrid(2)*ngrid(3)
207      value = BA_alloc_get(mt_dbl,nsize,
208     >                        'vl',vl_hndl,vl_indx)
209      value = value.and.BA_alloc_get(mt_dbl,nsize*(nprj),
210     >                        'vnl',vnl_hndl, vnl_indx)
211      value = value.and.BA_alloc_get(mt_dbl,nsize,
212     >                        'Ylm',Ylm_hndl, Ylm_indx)
213
214      value = value.and.BA_alloc_get(mt_dbl,nsize*(3),
215     >                        'G',G_hndl, G_indx)
216      value = value.and.BA_alloc_get(mt_int,nprj,
217     >                        'n_prj', n_prj_hndl, n_prj_indx)
218      value = value.and.BA_alloc_get(mt_int,nprj,
219     >                        'l_prj', l_prj_hndl, l_prj_indx)
220      value = value.and.BA_alloc_get(mt_int,nprj,
221     >                        'm_prj', m_prj_hndl, m_prj_indx)
222      value = value.and.BA_alloc_get(mt_int,nprj,
223     >                        'b_prj', b_prj_hndl, b_prj_indx)
224      if(.not.value)
225     >    call errquit('hghppv1: out of heap memory', 0, MA_ERR)
226
227
228      !**** determine n_prj, l_prj, and m_prj arrays ****
229      indx  = 0
230      nfft1 = 1
231      do l=0,lmax
232      do i=1,3
233        if ((H(i,l).ne.0.0d0)) then
234          do m=-l,l
235            int_mb(n_prj_indx+indx) = i
236            int_mb(l_prj_indx+indx) = l
237            int_mb(m_prj_indx+indx) = m
238            int_mb(b_prj_indx+indx) = nfft1
239            indx = indx + 1
240          end do
241          nfft1=nfft1+1
242        end if
243      end do
244      end do
245
246      filter = control_kbpp_filter()
247      ecut   = control_ecut()
248      wcut   = control_wcut()
249
250*     **** preparation of constants ****
251      nfft1=ngrid(1)
252      nfft2=ngrid(2)
253      nfft3=ngrid(3)
254      call setup_kbpp(nfft1,nfft2,nfft3,unita,unitg,dbl_mb(G_indx))
255      zv = dble(Zion)
256
257      call HGH_local(version,rlocal,
258     >               zv,rloc,C1,C2,C3,C4,
259     >               nfft1,nfft2,nfft3,
260     >               dbl_mb(G_indx),
261     >               dbl_mb(vl_indx))
262      if (filter)
263     >   call kbpp_filter(nfft1,nfft2,nfft3,
264     >                    dbl_mb(G_indx),
265     >                    ecut,
266     >                    dbl_mb(vl_indx))
267
268      do i=1,nprj
269         n=int_mb(n_prj_indx+i-1)
270         l=int_mb(l_prj_indx+i-1)
271         m=int_mb(m_prj_indx+i-1)
272         call Tesseral(l,m,
273     >                 nfft1,nfft2,nfft3,
274     >                 dbl_mb(G_indx),
275     >                 dbl_mb(Ylm_indx))
276         call HGH_nonlocal(n,l,
277     >                     r(l),
278     >                     nfft1,nfft2,nfft3,
279     >                     dbl_mb(G_indx),
280     >                     dbl_mb(vnl_indx+(i-1)*nsize))
281         if ((taskid.eq.MASTER).and.(oprint))
282     >      write(luout,*) "creating projector:",n,l,m
283         do j=1,nsize
284          dbl_mb(vnl_indx+(i-1)*nsize+j-1)
285     >    = dbl_mb(vnl_indx+(i-1)*nsize+j-1)*dbl_mb(Ylm_indx+j-1)
286         end do
287
288         if (filter)
289     >      call kbpp_filter(nfft1,nfft2,nfft3,
290     >                    dbl_mb(G_indx),
291     >                    wcut,
292     >                    dbl_mb(vnl_indx+(i-1)*nsize))
293      end do
294
295
296      if ((taskid.eq.MASTER).and.(oprint)) then
297      write(luout,*) "     ********************************************"
298      write(luout,*) "     *                                          *"
299      write(luout,*) "     *    HGHPPV1 - Pseudopotential Formatter   *"
300      write(luout,*) "     *                                          *"
301      write(luout,*) "     *      version last updated 11/13/03       *"
302      write(luout,*) "     *                                          *"
303      write(luout,*) "     *       developed by Eric J. Bylaska       *"
304      write(luout,*) "     *                                          *"
305      write(luout,*) "     ********************************************"
306      call nwpw_message(1)
307      write(luout,*)
308      write(luout,*) "Pseudpotential Data"
309      write(luout,*) "-------------------"
310      write(luout,*) "  atom     :",atom
311      write(luout,*) "  charge   :",Zion
312      write(luout,*) "  highest angular component used :",lmax
313      write(luout,*) "  highest radial  component used :",nmax
314      write(luout,*) "  number of non-local projectors :",nprj
315      write(luout,111) "   projector cutoffs: ",(r(i), i=0,lmax)
316      write(luout,*)
317      write(luout,111) " local psp cutoff       : ",rloc
318      write(luout,111) " local psp coefficients : ",C1,C2,C3,C4
319      if (version.eq.4)
320     >   write(luout,*) "  aperiodic cutoff radius        :",rlocal
321
322      write(luout,*)
323      write(luout,*) "Simulation Cell"
324      write(luout,*) "---------------"
325      if (version.eq.3) write(luout,112) "  boundry: periodic"
326      if (version.eq.4) write(luout,112) "  boundry: aperiodic"
327      write(luout,113) "  ngrid  :",ngrid
328      write(luout,114) "  unita  :",unita(1,1),unita(2,1),unita(3,1)
329      write(luout,114) "          ",unita(1,2),unita(2,2),unita(3,2)
330      write(luout,114) "          ",unita(1,3),unita(2,3),unita(3,3)
331      write(luout,*)
332  111 format(a,10f10.3)
333  112 format(a)
334  113 format(a,3I4)
335  114 format(a,3F10.3)
336  115 format(a,2E14.6)
337      end if
338
339
340
341
342      if (taskid.eq.MASTER) then
343      call util_file_name_noprefix(formatted_filename,
344     >                    .false.,
345     >                    .false.,
346     >                    full_filename)
347      l = index(full_filename,' ') - 1
348      if (mprint) then
349      write(luout,*)
350      write(luout,*) "Generated formatted_filename: ",full_filename(1:l)
351      if (filter) write(luout,*) "- filtering pseudopotential -"
352      !write(luout,*)
353      end if
354      call openfile(2,full_filename,l,'w',l)
355
356         call cwrite(2,comment,80)
357         call iwrite(2,psp_type,1)
358         call iwrite(2,version,1)
359         call iwrite(2,ngrid,3)
360         call dwrite(2,unita,9)
361         call cwrite(2,atom,2)
362         call dwrite(2,amass,1)
363         call dwrite(2,zv,1)
364         call iwrite(2,lmax,1)
365         call iwrite(2,locp,1)
366
367         call iwrite(2,nmax,1)
368         call dwrite(2,r,lmax+1)
369
370
371         call iwrite(2,nprj,1)
372         if (nprj.gt.0) then
373          call iwrite(2,int_mb(n_prj_indx),nprj)
374          call iwrite(2,int_mb(l_prj_indx),nprj)
375          call iwrite(2,int_mb(m_prj_indx),nprj)
376          call iwrite(2,int_mb(b_prj_indx),nprj)
377          call dwrite(2,Gijl,(nmax*nmax*(lmax+1)))
378         end if
379
380         if (version.eq.4) call dwrite(2,rlocal,1)
381         rcore = 0.0d0
382         call dwrite(2,rcore,1)
383
384         call dwrite(2,dbl_mb(vl_indx),nsize)
385         if (nprj.gt.0) then
386           call dwrite(2,dbl_mb(vnl_indx),nsize*nprj)
387         end if
388
389      call closefile(2)
390      end if
391
392*     **** free heap space ****
393      value = BA_free_heap(vl_hndl)
394      value = value.and.BA_free_heap(vnl_hndl)
395      value = value.and.BA_free_heap(Ylm_hndl)
396      value = value.and.BA_free_heap(G_hndl)
397      value = value.and.BA_free_heap(n_prj_hndl)
398      value = value.and.BA_free_heap(l_prj_hndl)
399      value = value.and.BA_free_heap(m_prj_hndl)
400      value = value.and.BA_free_heap(b_prj_hndl)
401      if(.not.value)
402     >  call errquit('hghppv1: deallocatin heap memory', 0, MA_ERR)
403
404
405      if ((taskid.eq.MASTER).and.(oprint)) call nwpw_message(4)
406      hghppv1 = value
407      return
408
409 9999 call errquit('hghppv1:Error reading psp_filename',0, DISK_ERR)
410
411      END
412
413
414
415
416