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