1
2*
3* $Id$
4*
5
6*     ***********************************
7*     *					*
8*     *	 	     borbs_init		*
9*     *					*
10*     ***********************************
11      logical function borbs_init()
12      implicit none
13
14#include "bafdecls.fh"
15#include "borbs.fh"
16
17*     **** local variables ****
18      integer npack1,nion
19      logical value
20
21*     **** external functions *****
22      integer  ion_nkatm_qm,ion_nion_qm
23      external ion_nkatm_qm,ion_nion_qm
24
25
26      call Pack_npack(1,npack1)
27      nion   = ion_nion_qm()
28      nkatmx = ion_nkatm_qm()
29
30*     **** allocate borb datastructure  and borb tag lists****
31      call borb_projector_init(2*nkatmx)
32      value = BA_alloc_get(mt_int,(nkatmx),'borbs',borbs(2),borbs(1))
33
34      value = value.and.
35     >        BA_alloc_get(mt_int,(nkatmx),'lmmax',lmmax(2),lmmax(1))
36      value = value.and.
37     >        BA_alloc_get(mt_int,(nkatmx),'lmax',lmax(2),lmax(1))
38      value = value.and.
39     >        BA_alloc_get(mt_int,(nkatmx),'locp',locp(2),locp(1))
40      value = value.and.
41     >        BA_alloc_get(mt_dbl,(nkatmx),'rcut',rcut(2),rcut(1))
42      value = value.and.
43     >        BA_alloc_get(mt_dbl,(nkatmx),'lmbda',lmbda(2),lmbda(1))
44
45      value = value.and.
46     >        BA_alloc_get(mt_int,(nion*norbs_max),
47     >                    'lmborb',lmborb(2),lmborb(1))
48      value = value.and.
49     >        BA_alloc_get(mt_int,(nion*norbs_max),
50     >                    'iaborb',iaborb(2),iaborb(1))
51      value = value.and.
52     >        BA_alloc_get(mt_int,(nion*norbs_max),
53     >                    'iiborb',iiborb(2),iiborb(1))
54      value = value.and.
55     >        BA_alloc_get(mt_int,(nion*norbs_max),
56     >                    'basisborb',basisborb(2),basisborb(1))
57
58      borbs_init = value
59      return
60      end
61
62
63*     ***********************************
64*     *					*
65*     *	 	     borbs_end		*
66*     *					*
67*     ***********************************
68      subroutine borbs_end()
69      implicit none
70
71#include "bafdecls.fh"
72#include "errquit.fh"
73#include "borbs.fh"
74
75*     **** local variables ****
76      logical value
77
78*     **** deallocate projector data ****
79      call borb_projector_end()
80
81      value = BA_free_heap(borbs(2))
82      value = value.and.BA_free_heap(lmmax(2))
83      value = value.and.BA_free_heap(lmax(2))
84      value = value.and.BA_free_heap(locp(2))
85      value = value.and.BA_free_heap(rcut(2))
86      value = value.and.BA_free_heap(lmbda(2))
87      value = value.and.BA_free_heap(lmborb(2))
88      value = value.and.BA_free_heap(iaborb(2))
89      value = value.and.BA_free_heap(iiborb(2))
90      value = value.and.BA_free_heap(basisborb(2))
91      if (.not. value) call errquit('borbs_end:freeing heap memory',0,
92     &       MA_ERR)
93
94      return
95      end
96
97
98*     ***********************************
99*     *					*
100*     *	 	     borbs_norbs	*
101*     *					*
102*     ***********************************
103
104      integer function borbs_norbs(ia)
105      implicit none
106      integer  ia
107
108#include "bafdecls.fh"
109#include "borbs.fh"
110
111      borbs_norbs = int_mb(lmmax(1)+ia-1)
112      return
113      end
114
115
116*     ***********************************
117*     *					*
118*     *	 	   borbs_nbasis		*
119*     *					*
120*     ***********************************
121      integer function borbs_nbasis()
122      implicit none
123
124#include "borbs.fh"
125
126      borbs_nbasis = ibasis
127      return
128      end
129
130*     ***********************************
131*     *                                 *
132*     *              borbs_lmax         *
133*     *                                 *
134*     ***********************************
135
136      integer function borbs_lmax(ia)
137      implicit none
138      integer  ia
139
140#include "bafdecls.fh"
141#include "borbs.fh"
142
143      borbs_lmax = int_mb(lmax(1)+ia-1)-1
144      return
145      end
146
147
148*     ***********************************
149*     *                                 *
150*     *              borbs_rcut         *
151*     *                                 *
152*     ***********************************
153      real*8 function borbs_rcut(ia)
154      implicit none
155      integer  ia
156
157#include "bafdecls.fh"
158#include "borbs.fh"
159
160      borbs_rcut = dbl_mb(rcut(1)+ia-1)
161      return
162      end
163
164*     ***********************************
165*     *                                 *
166*     *           borbs_lmbda           *
167*     *                                 *
168*     ***********************************
169      real*8 function borbs_lmbda(ia)
170      implicit none
171      integer  ia
172
173#include "bafdecls.fh"
174#include "borbs.fh"
175
176      borbs_lmbda = dbl_mb(lmbda(1)+ia-1)
177      return
178      end
179
180
181*     ***********************************
182*     *                                 *
183*     *              borbs_l            *
184*     *                                 *
185*     ***********************************
186
187      integer function borbs_l(ia,n)
188      implicit none
189      integer  ia
190      integer n          ! basis number
191
192#include "bafdecls.fh"
193#include "borbs.fh"
194
195*     *** local variables ***
196      integer l,m,lm
197
198      lm = int_mb(lmborb(1)+n-1)
199      l = 0
200      if (lm.eq.1) l = 0
201
202      if (lm.eq.2) l = 1
203      if (lm.eq.3) l = 1
204      if (lm.eq.4) l = 1
205
206      if (lm.eq.5) l = 2
207      if (lm.eq.6) l = 2
208      if (lm.eq.7) l = 2
209      if (lm.eq.8) l = 2
210      if (lm.eq.9) l = 2
211
212      if (lm.eq.10) l = 3
213      if (lm.eq.11) l = 3
214      if (lm.eq.12) l = 3
215      if (lm.eq.13) l = 3
216      if (lm.eq.14) l = 3
217      if (lm.eq.15) l = 3
218      if (lm.eq.16) l = 3
219
220      borbs_l = l
221      return
222      end
223
224
225*     ***********************************
226*     *                                 *
227*     *     borbs_get_basis_number      *
228*     *                                 *
229*     ***********************************
230
231      integer function borbs_get_basis_number(ii,lm)
232      implicit none
233      integer ii,lm
234
235#include "bafdecls.fh"
236#include "borbs.fh"
237
238      borbs_get_basis_number=int_mb(basisborb(1)+lm-1+(ii-1)*norbs_max)
239
240      return
241      end
242
243
244
245*     ***********************************
246*     *					*
247*     *	 	   borbs_normalize  	*
248*     *					*
249*     ***********************************
250      subroutine borbs_normalize()
251      implicit none
252
253#include "bafdecls.fh"
254#include "borbs.fh"
255
256*     **** local variables ****
257      integer lm,ia,npack1,nbrillq,nbq,shift
258      real*8  sum
259
260*     **** external functions ****
261      integer  Pneb_nbrillq,borb_projector_get_ptr
262      external Pneb_nbrillq,borb_projector_get_ptr
263
264
265      nbrillq = Pneb_nbrillq()
266
267*     **** Normalize atomic orbitals in k space *****
268      do nbq=1,nbrillq
269         do ia=1,nkatmx
270         do lm=1,int_mb(lmmax(1)+ia-1)
271
272            shift = borb_projector_get_ptr(int_mb(borbs(1)+ia-1),nbq,lm)
273            call Cram_rr_dot(nbq,dbl_mb(shift),dbl_mb(shift),sum)
274            sum = 1.0d0/dsqrt(sum)
275            call Cram_r_SMul1(nbq,sum,dbl_mb(shift))
276         end do
277         end do
278      end do
279
280      return
281      end
282
283*     ***********************************
284*     *					*
285*     *	 	   borbs_weight		*
286*     *					*
287*     ***********************************
288      real*8 function borbs_weight(n)
289      implicit none
290      integer n          ! basis number
291
292#include "bafdecls.fh"
293#include "borbs.fh"
294
295*     **** local variables ****
296      integer ia
297      real*8 zv,zcount
298
299*     **** external functions ****
300      real*8   psp_zv
301      external psp_zv
302
303      ia     = int_mb(iaborb(1)+n-1)
304      zcount = int_mb(lmmax(1)+ia-1)
305      zv     = psp_zv(ia)
306
307      borbs_weight = (zv/zcount)
308      return
309      end
310
311*     ***********************************
312*     *					*
313*     *	 	   borbs_borb  		*
314*     *					*
315*     ***********************************
316
317      subroutine borbs_borb(nbq,n,borb)
318      implicit none
319      integer nbq
320      integer n          ! basis number
321      complex*16 borb(*) ! return orbital
322
323#include "bafdecls.fh"
324#include "errquit.fh"
325#include "borbs.fh"
326
327*     **** local variables ****
328      logical value
329      integer lm,ia,ii
330      integer nfft3d,npack1,shift
331      integer exi(2)
332      complex*16 cxr
333
334*     **** external functions ****
335      logical  is_sORd
336      external is_sORd
337      integer  borb_projector_get_ptr
338      external borb_projector_get_ptr
339
340      call C3dB_nfft3d(1,nfft3d)
341      call Cram_npack(nbq,npack1)
342
343      value = BA_push_get(mt_dcpl,npack1,'exi', exi(2), exi(1))
344      if (.not. value) call errquit('borbs_borb:out of heap memory',0,
345     &       MA_ERR)
346
347*     **** structure factor ****
348      lm = int_mb(lmborb(1)+n-1)
349      ia = int_mb(iaborb(1)+n-1)
350      ii = int_mb(iiborb(1)+n-1)
351      call cstrfac_pack(nbq,ii,dcpl_mb(exi(1)))
352      call cstrfac_k(ii,nbq,cxr)
353      call zscal(npack1,cxr,dcpl_mb(exi(1)),1)
354
355      shift = borb_projector_get_ptr(int_mb(borbs(1)+ia-1),nbq,lm)
356
357
358*     **** phase factor does not matter therefore ****
359*     **** (-i)^l is the same as (i)^l in the     ****
360*     **** Rayleigh scattering formula            ****
361
362*     *** current function is s or d ****
363      if (is_sORd(lm,int_mb(lmax(1)+ia-1),
364     >                   int_mb(locp(1)+ia-1))
365     >        ) then
366         call Cram_rc_Mul(nbq,dbl_mb(shift),dcpl_mb(exi(1)),borb)
367
368*     *** current function is p or f ****
369      else
370         call Cram_irc_Mul(nbq,dbl_mb(shift),dcpl_mb(exi(1)),borb)
371      end if
372
373
374      value = BA_pop_stack(exi(2))
375      if (.not. value) call errquit('borbs_borb:freeing heap memory',0,
376     &       MA_ERR)
377
378      return
379      end
380
381
382*     ***********************************
383*     *					*
384*     *	 	   borbs_read 		*
385*     *					*
386*     ***********************************
387      subroutine borbs_read(fname,
388     >                      version,
389     >                       nfft,unita,
390     >                       atom,
391     >                       lmmax,lmax,locp,rcut,lmbda,
392     >                       npack1,borbs_tag,
393     >                       tmp,tmp2,
394     >                       ierr)
395      implicit none
396      character*50 fname
397      integer version
398      integer nfft(3)
399      real*8  unita(3,3)
400      character*2 atom
401      integer lmmax,lmax,locp
402      real*8 rcut,lmbda
403      integer npack1
404      integer borbs_tag
405      complex*16 tmp(*)
406      real*8     tmp2(*)
407      integer ierr
408
409#include "bafdecls.fh"
410#include "btdb.fh"
411#include "util.fh"
412
413*    *** local variables ***
414      logical mprint,value
415      integer MASTER,taskid,taskid_k
416      parameter(MASTER=0)
417      integer n,l,nbrillioun,nb,nbq,pk
418      integer msglen
419      integer iatom(2)
420      character*255 full_filename
421      real*8 kv(3)
422
423
424*     ***** external functions ****
425      integer  brillioun_nbrillioun,borb_projector_alloc
426      external brillioun_nbrillioun,borb_projector_alloc
427      integer  brillioun_nbrillq
428      external brillioun_nbrillq
429      real*8   brillioun_all_k
430      external brillioun_all_k
431      logical  control_print
432      external control_print
433
434      call Parallel_taskid(taskid)
435      call Parallel3d_taskid_k(taskid_k)
436      mprint = (taskid.eq.MASTER).and.control_print(print_medium)
437
438
439*     **** open fname binary file ****
440      if (taskid.eq.MASTER) then
441         call util_file_name_noprefix(fname,.false.,
442     >                             .false.,
443     >                       full_filename)
444         l = index(full_filename,' ') - 1
445         call openfile(5,full_filename,l,'r',l)
446         call iread(5,version,1)
447         call iread(5,nfft,3)
448         call dread(5,unita,9)
449         call cread(5,atom,2)
450         call iread(5,lmax,1)
451         call iread(5,locp,1)
452         call dread(5,rcut,1)
453         call dread(5,lmbda,1)
454
455         call iread(5,nbrillioun,1)
456         ierr = 0
457         if (nbrillioun.eq.brillioun_nbrillioun()) then
458           do nb=1,nbrillioun
459               call dread(5,kv,3)
460               if ((brillioun_all_k(1,nb).ne.kv(1)).or.
461     >             (brillioun_all_k(2,nb).ne.kv(2)).or.
462     >             (brillioun_all_k(3,nb).ne.kv(3)))
463     >          ierr = 1
464           end do
465           if (ierr.eq.1) then
466              if (mprint) then
467                 write(*,*)"Brillioun Zone Vectors do not match!"
468                 call flush(6)
469              end if
470           end if
471         else
472           if (mprint) then
473           write(*,*)"Brillioun Zone Points do not match!"
474           write(*,*)"NB = ",nbrillioun," not equal ",
475     >       brillioun_nbrillioun()
476           call flush(6)
477           end if
478           ierr = 1
479         end if
480      end if
481
482      msglen = 1
483      call Parallel_Brdcst_ivalues(MASTER,msglen,ierr)
484      if (ierr.ne.0) then
485         if (taskid.eq.MASTER) call closefile(5)
486         return
487      end if
488
489
490*     **** send header data to all processors ****
491      msglen = 1
492      call Parallel_Brdcst_ivalues(MASTER,msglen,version)
493      msglen = 3
494      call Parallel_Brdcst_ivalues(MASTER,msglen,nfft)
495      msglen = 9
496      call Parallel_Brdcst_values(MASTER,msglen,unita)
497
498      iatom(1) = ichar(atom(1:1))
499      iatom(2) = ichar(atom(2:2))
500      msglen = 2
501      call Parallel_Brdcst_ivalues(MASTER,msglen,iatom)
502      atom(1:1) = char(iatom(1))
503      atom(2:2) = char(iatom(2))
504
505      msglen = 1
506      call Parallel_Brdcst_ivalues(MASTER,msglen,lmax)
507      call Parallel_Brdcst_ivalues(MASTER,msglen,locp)
508      call Parallel_Brdcst_ivalues(MASTER,msglen,nbrillioun)
509      call Parallel_Brdcst_values(MASTER,msglen,rcut)
510      call Parallel_Brdcst_values(MASTER,msglen,lmbda)
511      lmmax=(lmax+1)**2 - (2*locp+1)
512
513
514
515*     **** read in borb 3d blocks ****
516      borbs_tag = borb_projector_alloc(brillioun_nbrillq(),lmmax,npack1)
517      do nb=1,brillioun_nbrillioun()
518         call K1dB_ktoqp(nb,nbq,pk)
519         do n=1,lmmax
520            call C3dB_r_read(1,5,tmp2,tmp,-1,pk,.true.)
521            if (pk.eq.taskid_k) then
522               call Cram_r_pack(nbq,tmp2)
523               call borb_projector_add(borbs_tag,nbq,n,tmp2)
524            end if
525         end do
526      end do
527
528*     *** close fname binary file ***
529      if (taskid.eq.MASTER) then
530         call closefile(5)
531      end if
532
533      ierr = 0
534      return
535      end
536
537*     ***********************************
538*     *					*
539*     *	 	  borbs_readall		*
540*     *					*
541*     ***********************************
542      logical function borbs_readall()
543      implicit none
544
545#include "bafdecls.fh"
546#include "borbs.fh"
547
548*     **** local variables ****
549      integer ngp(3),version,nfft3d,npack1
550      integer ia,l,lm,ii,icount
551      real*8 unita(3,3)
552      integer tmp(2),tmp2(2),ierr
553      logical value,found,correct_box,value2
554      character*2  atom
555      character*4  element
556      character*50 fname
557
558*     **** parallel i/o variable ****
559      integer MASTER,taskid
560      parameter(MASTER=0)
561
562*     **** external functions ****
563      logical      nwpw_filefind
564      integer      control_ngrid
565      real*8       control_unita
566      character*12 control_boundry
567      character*4  ion_atom_qm
568      external     nwpw_filefind
569      external     control_ngrid
570      external     control_unita
571      external     control_boundry
572      external     ion_atom_qm
573      integer      ion_nion_qm,ion_katm_qm
574      external     ion_nion_qm,ion_katm_qm
575
576      call C3dB_nfft3d(1,nfft3d)
577      call Cram_max_npack(npack1)
578      call Parallel_taskid(taskid)
579
580      value = BA_push_get(mt_dbl,(4*nfft3d),'tmp',tmp(2),tmp(1))
581      if (.not. value) go to 9000
582
583      value = BA_push_get(mt_dbl,(2*nfft3d),'tmp2',tmp2(2),tmp2(1))
584      if (.not. value) go to 9000
585
586*     **** read pseudopotentials ****
587      do ia=1,nkatmx
588
589*      **** define formatted borb name ****
590       element = '    '
591       element = ion_atom_qm(ia)
592       l = index(element,' ') - 1
593       fname = element(1:l)//'.borb'
594
595
596       found = .false.
597       do while (.not.found)
598
599         if (nwpw_filefind(fname)) then
600            call borbs_read(fname,
601     >                  version,
602     >                  ngp,unita,
603     >                  atom,
604     >                  int_mb(lmmax(1)+ia-1),
605     >                  int_mb(lmax(1)+ia-1),
606     >                  int_mb(locp(1)+ia-1),
607     >                  dbl_mb(rcut(1)+ia-1),
608     >                  dbl_mb(lmbda(1)+ia-1),
609     >                  npack1,
610     >                  int_mb(borbs(1)+ ia-1),
611     >                  dbl_mb(tmp(1)),dbl_mb(tmp2(1)),
612     >                  ierr)
613
614           if (ierr.gt.0) then
615              value = .false.
616              go to 9000
617           end if
618
619*          **************************************************************
620*          ***** logic for finding out if psp is correctly formatted ****
621*          **************************************************************
622           correct_box = .true.
623           if ( (ngp(1).ne.control_ngrid(1)) .or.
624     >       (ngp(2).ne.control_ngrid(2)) .or.
625     >       (ngp(3).ne.control_ngrid(3)) .or.
626     >       (unita(1,1).ne.control_unita(1,1)) .or.
627     >       (unita(2,1).ne.control_unita(2,1)) .or.
628     >       (unita(3,1).ne.control_unita(3,1)) .or.
629     >       (unita(1,2).ne.control_unita(1,2)) .or.
630     >       (unita(2,2).ne.control_unita(2,2)) .or.
631     >       (unita(3,2).ne.control_unita(3,2)) .or.
632     >       (unita(1,3).ne.control_unita(1,3)) .or.
633     >       (unita(2,3).ne.control_unita(2,3)) .or.
634     >       (unita(3,3).ne.control_unita(3,3))) then
635              correct_box = .false.
636              if (taskid.eq.MASTER) then
637              write(6,*) "atomic orbitals are not correctly formatted:",
638     >                    fname
639              end if
640           end if
641           if (correct_box) found = .true.
642
643         end if
644
645*        **** generate formatted pseudopotential atom.borb *****
646         if (.not.found) then
647             call borbs_formatter_auto(ion_atom_qm(ia),0.0d0,0.0d0)
648         end if
649
650       end do !***do while ****
651
652
653      end do
654
655*     ***********************************************
656*     **** set up the index for the atomic basis ****
657*     ***********************************************
658      icount = 0
659      do ii=1,ion_nion_qm()
660         ia = ion_katm_qm(ii)
661
662         do lm=1,int_mb(lmmax(1)+ia-1)
663            icount = icount + 1
664
665            int_mb(lmborb(1)+icount-1)  = lm
666            int_mb(iaborb(1)+icount-1)  = ia
667            int_mb(iiborb(1)+icount-1)  = ii
668            int_mb(basisborb(1)+lm-1+(ii-1)*norbs_max) = icount
669         end do
670      end do
671      ibasis = icount
672      call borbs_normalize()
673
674 9000 value2 = BA_pop_stack(tmp2(2))
675      value2 = BA_pop_stack(tmp(2))
676
677
678      borbs_readall = value
679      return
680      end
681
682