1*
2* $Id$
3*
4
5*     ***********************************
6*     *             			*
7*     *           v_psi_read		*
8*     *             			*
9*     ***********************************
10
11      subroutine v_psi_read(ispin0,ne0,psi2)
12      implicit none
13      integer ispin0,ne0(2)
14      double complex psi2(*)
15      integer occupation
16
17#include "bafdecls.fh"
18#include "util.fh"
19#include "errquit.fh"
20
21*    *** local variables ***
22      integer ispin,ne(2)
23      integer version,l
24      integer nfft3d,npack1
25      integer nfft(3)
26      real*8  unita(3,3)
27      character*50 filename
28      character*255 full_filename
29
30      integer MASTER,taskid,taskid_i,taskid_j,taskid_p,com_p
31      parameter(MASTER=0)
32      integer n,q,pj
33      integer msglen
34
35c     complex*16 tmp(*)
36      integer tmp(2),tmp2(2)
37      logical value,readblock,lprint,pio
38
39*     ***** local functions ****
40      character*50 control_input_v_psi
41      external     control_input_v_psi
42      double precision control_unita
43      external         control_unita
44      integer  control_ngrid
45      external control_ngrid
46      logical  control_print,control_parallel_io
47      external control_print,control_parallel_io
48
49      call nwpw_timing_start(50)
50      call Parallel_taskid(taskid)
51      call Parallel2d_taskid_i(taskid_i)
52      call Parallel2d_taskid_j(taskid_j)
53      call D3dB_nfft3d(1,nfft3d)
54      call Pack_npack(1,npack1)
55
56      pio = control_parallel_io()
57      if (pio) then
58         taskid_p = taskid_i
59         com_p = 1
60      else
61         taskid_p = taskid
62         com_p = 0
63      end if
64
65      lprint= ((taskid.eq.MASTER).and.control_print(print_low))
66
67      value = BA_push_get(mt_dcpl,nfft3d,'tmp',tmp(2),tmp(1))
68      value = value.and.
69     >        BA_push_get(mt_dcpl,nfft3d,'tmp2',tmp2(2),tmp2(1))
70      if (.not. value)
71     > call errquit('v_psi_read:out of stack memory',0,MA_ERR)
72
73*     **** open ELCIN binary file ****
74      if (taskid_p.eq.MASTER) then
75         filename = control_input_v_psi()
76         full_filename = filename
77         call util_file_name_resolve(full_filename, .false.)
78c        call util_file_name_noprefix(filename,.false.,
79c    >                                .false.,
80c    >                        full_filename)
81         l = index(full_filename,' ') -1
82         if (lprint) write(*,1210) full_filename(1:l)
83 1210    FORMAT(/' input vpsi filename:',A)
84
85
86         l = index(full_filename,' ') -1
87         call openfile(5,full_filename,l,'r',l)
88         call iread(5,version,1)
89         call iread(5,nfft,3)
90         call dread(5,unita,9)
91         call iread(5,ispin,1)
92         call iread(5,ne,2)
93         call iread(5,occupation,1)
94      end if
95
96c     **** send header to all nodes ****
97      msglen = 1
98      call Parallela_Brdcst_ivalues(com_p,MASTER,msglen,version)
99      msglen = 3
100      call Parallela_Brdcst_ivalues(com_p,MASTER,msglen,nfft)
101      msglen = 9
102      call Parallela_Brdcst_values(com_p,MASTER,msglen,unita)
103      msglen = 1
104      call Parallela_Brdcst_ivalues(com_p,MASTER,msglen,ispin)
105      call Parallela_Brdcst_ivalues(com_p,MASTER,msglen,occupation)
106      msglen = 2
107      call Parallela_Brdcst_ivalues(com_p,MASTER,msglen,ne)
108
109
110*     ***** Error checking ****
111      readblock=.true.
112      if ( (nfft(1).ne.control_ngrid(1)) .or.
113     >     (nfft(2).ne.control_ngrid(2)) .or.
114     >     (nfft(3).ne.control_ngrid(3)) ) then
115         readblock = .false.
116         if (taskid_p.eq.MASTER) then
117            write(*,*) "Error reading v_psi, v_psi set to zero"
118            write(*,*) " - bad ngrid "
119            write(*,*) " - vpsi_ngrid(1)=",nfft(1),
120     >                 " ngrid(1)=",control_ngrid(1)
121            write(*,*) " - vpsi_ngrid(2)=",nfft(2),
122     >                 " ngrid(2)=",control_ngrid(2)
123            write(*,*) " - vpsi_ngrid(3)=",nfft(3),
124     >                 " ngrid(3)=",control_ngrid(3)
125         end if
126      end if
127
128      if ( (unita(1,1).ne.control_unita(1,1)) .or.
129     >     (unita(2,1).ne.control_unita(2,1)) .or.
130     >     (unita(3,1).ne.control_unita(3,1)) .or.
131     >     (unita(1,2).ne.control_unita(1,2)) .or.
132     >     (unita(2,2).ne.control_unita(2,2)) .or.
133     >     (unita(3,2).ne.control_unita(3,2)) .or.
134     >     (unita(1,3).ne.control_unita(1,3)) .or.
135     >     (unita(2,3).ne.control_unita(2,3)) .or.
136     >     (unita(3,3).ne.control_unita(3,3)) .or.
137     >     (ispin     .ne.ispin0)             .or.
138     >     (ne(1)     .ne.ne0(1))             .or.
139     >     (ne(2)     .ne.ne0(2))  ) then
140         readblock = .false.
141         if (taskid.eq.MASTER) then
142            write(6,*) "Error reading v_psi, v_psi set to zero"
143            write(6,*) " - incorrect unitcell,ispin, or ne"
144         end if
145      end if
146
147
148*     *** read in 3d blocks ***
149      if (readblock) then
150      do n=1,(ne(1)+ne(2))
151         call Dneall_ntoqp(n,q,pj)
152         if (pio) then
153            call D3dB_c_read_pio(1,5,dcpl_mb(tmp2(1)),
154     >                               dcpl_mb(tmp(1)),pj)
155         else
156            call D3dB_c_read(1,5,dcpl_mb(tmp2(1)),
157     >                           dcpl_mb(tmp(1)),pj)
158         end if
159         if (pj.eq.taskid_j) then
160           call Pack_c_pack(1,dcpl_mb(tmp2(1)))
161           call Pack_c_Copy(1,dcpl_mb(tmp2(1)),psi2(1+(q-1)*npack1))
162         end if
163      end do
164      end if
165
166*     *** close ELCIN binary file ***
167      if (taskid_p.eq.MASTER) then
168        call closefile(5)
169      end if
170
171      value =           BA_pop_stack(tmp2(2))
172      value = value.and.BA_pop_stack(tmp(2))
173      if (.not. value)
174     >  call errquit('v_psi_read:error popping stack',0,MA_ERR)
175
176      call nwpw_timing_end(50)
177      return
178      end
179