1*
2* $Id$
3*
4
5*     ***********************************************
6*     *                                             *
7*     *               wvfnc_adjust                  *
8*     *                                             *
9*     ***********************************************
10
11      subroutine wvfnc_adjust(wavefunction_filename,ispin,nein)
12      implicit none
13      character*50 wavefunction_filename
14      integer      ispin,nein(2)
15
16#include "bafdecls.fh"
17#include "errquit.fh"
18
19*     **** local variables ****
20      logical value,fractional
21      integer MASTER,taskid
22      parameter (MASTER=0)
23
24      integer NMAX
25      integer filling(2)
26      integer fractional_orbitals(2),ne(2)
27      character*255 new_filename,old_filename,emo_filename
28
29*     **** external functions ****
30      logical  control_fractional
31      integer  control_fractional_orbitals
32      external control_fractional
33      external control_fractional_orbitals
34      character*50 control_input_epsi
35      external     control_input_epsi
36
37      ne(1) = nein(1)
38      ne(2) = nein(2)
39      fractional = control_fractional()
40      if (fractional) then
41         fractional_orbitals(1) = control_fractional_orbitals(1)
42         ne(1) = nein(1) + fractional_orbitals(1)
43         if (ispin.eq.2) then
44           fractional_orbitals(2) = control_fractional_orbitals(2)
45           ne(2) = nein(2) + fractional_orbitals(2)
46         end if
47      end if
48
49      NMAX = ne(1)+ne(2)
50      call Parallel_taskid(taskid)
51      if (taskid.eq.MASTER) then
52        value = BA_push_get(mt_int,8*NMAX,
53     >                    'filling',filling(2),filling(1))
54        if (.not. value)
55     >    call errquit('wvfnc_adjust:out of stack memory',0,MA_ERR)
56
57        call util_file_name_noprefix('wvfnc_adjust',
58     >                    .false.,
59     >                    .false.,
60     >                     old_filename)
61        call util_file_name_noprefix(wavefunction_filename,
62     >                    .false.,
63     >                    .false.,
64     >                    new_filename)
65        call util_file_copy(new_filename,old_filename)
66
67        !*** fetch the emovecs filename ***
68        call util_file_name_noprefix(control_input_epsi(),
69     >                    .false.,
70     >                    .false.,
71     >                    emo_filename)
72
73        !*** adjust wavefunctions ***
74        call sub_wvfnc_adjust(NMAX,int_mb(filling(1)),
75     >                     new_filename,
76     >                     old_filename,
77     >                     emo_filename,
78     >                     ispin,
79     >                     ne,
80     >                     fractional,
81     >                     fractional_orbitals)
82
83        !*** remove temporary wvfnc_adjust file ***
84        call util_file_unlink(old_filename)
85
86
87        write(*,*) "wavefunction adjust, new psi:",
88     >             wavefunction_filename
89        write(*,*) "-   spin, nalpha, nbeta:",ispin,ne
90        value = BA_pop_stack(filling(2))
91        if (.not. value) call errquit('popping stack memory',0, MA_ERR)
92      end if
93      call ga_sync()
94
95      return
96      end
97
98
99      subroutine sub_wvfnc_adjust(NMAX,filling,
100     >                         new_filename,
101     >                         old_filename,
102     >                         emo_filename,
103     >                         ispin,
104     >                         ne,
105     >                         fractional,
106     >                         frac_orb)
107      implicit none
108      integer NMAX
109      integer filling(4,NMAX,2)
110      character*255 new_filename
111      character*255 old_filename
112      character*255 emo_filename
113      integer      ispin,ne(2)
114      logical      fractional
115      integer      frac_orb(2)
116
117#include "bafdecls.fh"
118#include "errquit.fh"
119
120      logical value,emo_found,emo_used
121      character*255 full_filename
122
123      integer      version
124      integer      ngrid(3)
125      real*8       unita(3,3)
126
127      integer nfft1,nfft2,nfft3,nfft3d,n2ft3d
128      integer inc2c,inc3c
129      integer cfull_indx,cfull_hndl,l,l1,l2
130      integer i,j,k,ms,n,occupation
131      integer n0,ms0,n0max,ispin0,ne0(2)
132      integer nx,msx,nxmax,ispinx,nex(2)
133
134      double precision p,scale
135      double complex cc,cx,sx,zx,zc,rx,ry
136
137*     **** external functions ****
138      double precision GCDOTC,util_random
139      external         GCDOTC,util_random
140
141
142      p = util_random(5291999) !*** initialize the random sequence ****
143
144      call getfilling(.true.,ne(1),filling)
145      if (ispin.eq.2) call getfilling(.true.,ne(2),filling(1,1,2))
146
147
148      scale=1.0d0/dsqrt(2.0d0)
149      zx=(1.0d0,0.0d0)
150      sx=(0.0d0,1.0d0)*scale
151      cx=(1.0d0,0.0d0)*scale
152
153
154*     **** write wavefunction in CPMDV3 format ****
155      l = index(old_filename,' ') - 1
156      call openfile(5,old_filename,l,'r',l)
157      call iread(5,version,1)
158      call iread(5,ngrid,3)
159      call dread(5,unita,9)
160      call iread(5,ispin0,1)
161      call iread(5,ne0,2)
162      call iread(5,occupation,1)
163
164      l = index(new_filename,' ') - 1
165      call openfile(6,new_filename,l,'w',l)
166      call iwrite(6,version,1)
167      call iwrite(6,ngrid,3)
168      call dwrite(6,unita,9)
169      call iwrite(6,ispin,1)
170      call iwrite(6,ne,2)
171      if (fractional) then
172         occupation = ispin
173      else
174         occupation = -1
175      end if
176      call iwrite(6,occupation,1)
177
178
179*     ***** constants *****
180      nfft1=ngrid(1)
181      nfft2=ngrid(2)
182      nfft3=ngrid(3)
183      nfft3d=(nfft1/2+1)*nfft2*nfft3
184      n2ft3d=2*nfft3d
185      inc2c = nfft1/2+1
186      inc3c =inc2c*nfft2
187
188*     ***** allocate wavefunction memory ****
189      value = BA_push_get(mt_dcpl,2*nfft3d,
190     >                     'cfull',cfull_hndl,cfull_indx)
191      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
192
193*     **** modularize the filling ****
194      do ms=1,ispin
195        do n=1,ne(ms)
196           i = filling(1,n,ms)
197           j = filling(2,n,ms)
198           k = filling(3,n,ms)
199           filling(1,n,ms) = mod(i+inc2c,inc2c)
200           filling(2,n,ms) = mod(j+nfft2,nfft2)
201           filling(3,n,ms) = mod(k+nfft3,nfft3)
202        end do
203      end do
204
205*     **** try to read emo_filename ***
206      nex(1) = 0
207      nex(2) = 0
208      emo_found = .false.
209      emo_used  = .false.
210      inquire(file=emo_filename,exist=emo_found)
211      if (emo_found) then
212         l = index(emo_filename,' ') - 1
213         call openfile(3,emo_filename,l,'r',l)
214         call iread(3,version,1)
215         call iread(3,ngrid,3)
216         call dread(3,unita,9)
217         call iread(3,ispinx,1)
218         call iread(3,nex,2)
219         call iread(3,occupation,1)
220      end if
221
222
223
224      ms0 = 1
225      msx = 1
226      do ms=1,ispin
227         n0 = 1
228         nx = 1
229         n0max = ne0(ms0)
230         nxmax = nex(ms0)
231         do n=1,ne(ms)
232
233            !*** read from old filename ***
234            if (n.le.n0max) then
235               call dread(5,dcpl_mb(cfull_indx),n2ft3d)
236               n0 = n0 + 1
237
238            !*** read from emo_filename ***
239            else if (n.le.(n0max+nxmax)) then
240               call dread(3,dcpl_mb(cfull_indx),n2ft3d)
241               nx = nx + 1
242               emo_used  = .true.
243               call sub_wvfnc_project_out(n2ft3d,
244     >                                  dcpl_mb(cfull_indx),
245     >                                  dcpl_mb(cfull_indx+nfft3d),
246     >                                  old_filename)
247               P=GCDOTC(nfft1,nfft2,nfft3,
248     >                     dcpl_mb(cfull_indx),
249     >                     dcpl_mb(cfull_indx))
250               P=1.0d0/dsqrt(P)
251               call dscal(n2ft3d,P,dcpl_mb(cfull_indx),1)
252
253            !*** generate new random wavefunction ***
254            else
255               call dcopy(n2ft3d,0.0d0,0,dcpl_mb(cfull_indx),1)
256               l1= inc3c*filling(3,n,ms)
257     >           + inc2c*filling(2,n,ms)
258     >           +       filling(1,n,ms)
259               if (filling(4,n,ms).lt.0) cc=sx
260               if (filling(4,n,ms).eq.0) cc=zx
261               if (filling(4,n,ms).gt.0) cc=cx
262                l2=l1
263                dcpl_mb(cfull_indx+l1)=cc
264                if (filling(1,n,ms).eq.0) then
265                  l2 = inc3c*mod(nfft3-filling(3,n,ms),nfft3)
266     >               + inc2c*mod(nfft2-filling(2,n,ms),nfft2)
267     >               +       filling(1,n,ms)
268                  dcpl_mb(cfull_indx+l2)=dconjg(cc)
269                end if
270c                if((ABS(filling(4,n,ms)).gt.1)) then
271                  do 125 k=0,nfft3d-1
272                    dcpl_mb(cfull_indx+k) = dcpl_mb(cfull_indx+k)
273     >                 + dcmplx((0.5d0-util_random(0)),
274     >                          (0.5d0-util_random(0)))
275     >                    /dsqrt(dble(nfft3d))
276  125             continue
277                  zc = dcpl_mb(cfull_indx)
278                  dcpl_mb(cfull_indx) = dcmplx(dble(zc),0.0d0)
279                  call gctimereverse(nfft1,nfft2,nfft3,
280     >                               dcpl_mb(cfull_indx))
281
282                  call sub_wvfnc_project_out(n2ft3d,
283     >                                  dcpl_mb(cfull_indx),
284     >                                  dcpl_mb(cfull_indx+nfft3d),
285     >                                  old_filename)
286
287                  P=gcdotc(nfft1,nfft2,nfft3,
288     >                     dcpl_mb(cfull_indx),
289     >                     dcpl_mb(cfull_indx))
290                  P=1.0d0/dsqrt(P)
291                  call dscal(n2ft3d,P,dcpl_mb(cfull_indx),1)
292c                end if
293            end if
294
295            call dwrite(6,dcpl_mb(cfull_indx),n2ft3d)
296            !n0 = n0 + 1
297         end do
298
299         ms0 = ms0 + 1
300         msx = msx + 1
301
302*        **** rewind the wavefunction read ****
303         if (ms0.gt.ispin0) then
304            ms0 = 1
305            call closefile(5)
306            l = index(old_filename,' ') - 1
307            call openfile(5,old_filename,l,'r',l)
308            call iread(5,version,1)
309            call iread(5,ngrid,3)
310            call dread(5,unita,9)
311            call iread(5,ispin0,1)
312            call iread(5,ne0,2)
313            call iread(5,occupation,1)
314         end if
315
316*        **** read remaining wvfunctions in spin ****
317         if ((ispin0.eq.2).and.(ispin.eq.2).and.(n0.le.n0max)) then
318           do n=n0,n0max
319             call dread(5,dcpl_mb(cfull_indx),n2ft3d)
320           end do
321         end if
322
323*        **** read remaining emo_filename wavefunctions in spin ***
324         if (emo_found) then
325            if (msx.gt.ispinx) then
326               msx = 1
327               call closefile(3)
328               l = index(emo_filename,' ') - 1
329               call openfile(3,emo_filename,l,'r',l)
330               call iread(3,version,1)
331               call iread(3,ngrid,3)
332               call dread(3,unita,9)
333               call iread(3,ispinx,1)
334               call iread(3,nex,2)
335               call iread(3,occupation,1)
336            end if
337
338*           **** read remaining wvfunctions in spin ****
339            if ((ispinx.eq.2).and.(ispin.eq.2).and.(nx.le.nxmax)) then
340              do n=nx,nxmax
341                call dread(3,dcpl_mb(cfull_indx),n2ft3d)
342              end do
343            end if
344         end if
345
346      end do
347
348c     **** add occupation - don't use previous start from scratch ****
349      if (fractional) then
350         rx = 1.0d0
351         ry = 0.0d0
352         do ms=1,ispin
353           do n=1,ne(ms)
354             if (n.le.ne0(ms)) then
355                call dwrite(6,rx,1)
356             else
357                call dwrite(6,ry,1)
358             end if
359           end do
360         end do
361      end if
362
363      if (emo_found) then
364         call closefile(3)
365         if (emo_used) call util_file_unlink(emo_filename) !*** remove emo_filename ***
366      end if
367      call closefile(5)
368      call closefile(6)
369
370      value = BA_pop_stack(cfull_hndl)
371      if (.not. value) call errquit('popping stack memory',0, MA_ERR)
372
373      return
374      end
375
376      subroutine sub_wvfnc_project_out(n2ft3d,psi,tpsi,old_filename)
377      implicit none
378      integer n2ft3d
379      complex*16 psi(*)
380      complex*16 tpsi(*)
381      character*255 old_filename
382
383      integer l,version,ngrid(3)
384      integer i,occupation
385      integer ispin0,ne0(2)
386
387      double precision p,unita(9)
388
389*     **** external functions ****
390      double precision GCDOTC
391      external         GCDOTC
392
393      l = index(old_filename,' ') - 1
394      call openfile(4,old_filename,l,'r',l)
395      call iread(4,version,1)
396      call iread(4,ngrid,3)
397      call dread(4,unita,9)
398      call iread(4,ispin0,1)
399      call iread(4,ne0,2)
400      call iread(4,occupation,1)
401      do i=1,ne0(1)
402         call dread(4,tpsi,n2ft3d)
403         p = GCDOTC(ngrid(1),ngrid(2),ngrid(3),tpsi,psi)
404         call daxpy(n2ft3d,-p,tpsi,1,psi,1)
405      end do
406      call closefile(4)
407
408      return
409      end
410
411