1*
2* $Id$
3*
4
5
6*     ***********************************
7*     *                                 *
8*     *          paw_psi_MakeOrtho      *
9*     *                                 *
10*     ***********************************
11
12*   This routine orthonormalizes the orbitals using a modified
13* Gram-Schmidt algorithm.
14*
15      subroutine paw_psi_MakeOrtho(npack,ne,psi)
16      implicit none
17      integer npack,ne
18      double complex psi(npack,ne)
19
20*     **** local variables ****
21      integer j,k
22      real*8  w
23
24      do k=1,ne
25         call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,k),w)
26         w = 1.0d0/dsqrt(w)
27c         call Pack_c_SMul(1,w,psi(1,k),psi(1,k))
28         call Pack_c_SMul1(1,w,psi(1,k))
29
30         do j=k+1,ne
31            call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,j),w)
32            w = -w
33            call Pack_cc_daxpy(1,w,psi(1,k),psi(1,j))
34         end do
35      end do
36
37      return
38      end
39
40*     ***********************************
41*     *                                 *
42*     *          paw_psi_CheckOrtho     *
43*     *                                 *
44*     ***********************************
45
46*   This routine return true if the orbitals are
47* orthonormal.
48
49      real*8 function paw_psi_CheckOrtho(npack,ne,psi)
50      implicit none
51      integer npack,ne
52      double complex psi(npack,ne)
53
54
55*     **** local variables ****
56      integer j,k
57      real*8  w,error
58
59      error = 0.0d0
60      do k=1,ne
61        call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,k),w)
62        error = error + dabs(1.0d0-w)
63
64        do j=k+1,ne
65           call paw_overlap_matrix_gen(1,1,psi(1,j),psi(1,k),w)
66           error = error + dabs(w)
67        end do
68      end do
69
70      paw_psi_CheckOrtho = error
71      return
72      end
73
74
75*     ***********************************
76*     *                                 *
77*     *          paw_psi_CheckOrtho2    *
78*     *                                 *
79*     ***********************************
80
81*   This routine return true if the orbitals are
82* orthonormal.
83
84      subroutine paw_psi_CheckOrtho2(npack,ne,psi)
85      implicit none
86      integer npack,ne
87      double complex psi(npack,ne)
88
89
90*     **** local variables ****
91      integer j,k
92      real*8  w,error
93
94
95      write(*,*)
96      do k=1,ne
97        call paw_overlap_matrix_gen(1,1,psi(1,k),psi(1,k),w)
98        write(*,*) "CheckOrtho2:",k,k,w
99
100        do j=k+1,ne
101           call paw_overlap_matrix_gen(1,1,psi(1,j),psi(1,k),w)
102        write(*,*) "CheckOrtho2:",j,k,w
103        end do
104      end do
105      write(*,*)
106
107      return
108      end
109
110
111
112*     ***********************************
113*     *                                 *
114*     *          paw_psi_lmbda          *
115*     *                                 *
116*     ***********************************
117
118      subroutine paw_psi_lmbda(ispin,ne,nemax,npack1,
119     >                     psi1,psi2,
120     >                     dte,
121     >                     lmbda,tmp,ierr)
122
123      implicit none
124      integer ispin,ne(2),nemax,npack1
125      complex*16 psi1(npack1,nemax)
126      complex*16 psi2(npack1,nemax)
127      real*8     dte
128      real*8     lmbda(*)
129      real*8     tmp(*)
130      integer    ierr
131
132      integer MASTER
133      parameter (MASTER=0)
134
135*     **** local variables ****
136      logical failed
137      integer taskid
138      integer n1(2),n2(2)
139      integer i,j,ii,jj,ms
140      integer n,nn,index
141      integer st1,st2
142      integer A,B,C,U,D,Ba,Bs,fnm
143      integer sl(2)
144      real*8  alpha
145
146
147
148      call nwpw_timing_start(3)
149
150      call Parallel_taskid(taskid)
151
152      n    = ne(1)
153      nn   = n**2
154
155      A    = 0*nn + 1
156      B    = 1*nn + 1
157      C    = 2*nn + 1
158      Ba   = 3*nn + 1
159      Bs   = 4*nn + 1
160      fnm  = 5*nn + 1
161      st1  = 6*nn + 1
162      D    = 7*nn + 1
163
164      U    = Bs
165      st2  = B
166
167      call dcopy(8*nn,0.0d0,0,tmp,1)
168
169      sl(1)  = 0*nn + 1
170      sl(2)  = 1*nn + 1
171      call dcopy(2*nn,0.0d0,0,lmbda,1)
172
173      n1(1)=1
174      n2(1)=ne(1)
175      n1(2)=ne(1)+1
176      n2(2)=ne(1)+ne(2)
177
178
179      do ms=1,ispin
180        IF(ne(ms).le.0) go to 640
181
182
183*       ***** compute the overlap matrices ****
184        call paw_overlap_sym_matrix_gen(n,ne(ms),
185     >                          psi2(1,n1(ms)),
186     >                          psi2(1,n1(ms)),
187     >                          tmp(A))
188        call paw_overlap_matrix_gen(n,ne(ms),
189     >                          psi1(1,n1(ms)),
190     >                          psi2(1,n1(ms)),
191     >                          tmp(B))
192        call paw_overlap_sym_matrix_gen(n,ne(ms),
193     >                          psi1(1,n1(ms)),
194     >                          psi1(1,n1(ms)),
195     >                          tmp(C))
196
197
198        call paw_psi_gen_Ba_Bs(n,ne(ms),tmp(B),tmp(Bs),tmp(Ba))
199        call paw_psi_gen_UD(n,ne(ms),tmp(Bs),tmp(D),lmbda)
200
201
202        call paw_psi_gen_X(n,ne(ms),tmp(st1),tmp(st2),
203     >                     tmp(A),tmp(Ba),tmp(C),
204     >                     tmp(U),tmp(D),tmp(fnm),
205     >                     failed)
206
207        if (failed) then
208          if (taskid.eq.MASTER) then
209            write(6,*)
210     >     'Warning: Lagrange Multiplier generation failed.'
211            write(6,*) '        +Try using a smaller time step'
212            write(6,*) '        +Gram-Schmidt being performed, spin:',ms
213          end if
214          call paw_psi_MakeOrtho(npack1,ne(ms),psi2(1,n1(ms)))
215        else
216          call dcopy(n*ne(ms),tmp(st1),1,lmbda(sl(ms)),1)
217          call dscal(n*ne(ms),(1.0d0/dte),lmbda(sl(ms)),1)
218
219
220*         ****  correction due to the constraint ****
221          call dgemm('N','N',2*npack1,ne(ms),ne(ms),
222     >              (1.0d0),
223     >              psi1(1,n1(ms)),2*npack1,
224     >              tmp(st1),n,
225     >              (1.0d0),
226     >              psi2(1,n1(ms)),2*npack1)
227
228        end if
229  640   continue
230      end do !*ms*
231      call nwpw_timing_end(3)
232
233      return
234      end
235
236
237
238*     ***********************************
239*     *                                 *
240*     *        paw_psi_gen_Ba_Bs        *
241*     *                                 *
242*     ***********************************
243      subroutine paw_psi_gen_Ba_Bs(n_max,n,B,Bs,Ba)
244      implicit none
245      integer n_max,n
246      real*8 B(n_max,n)
247      real*8 Bs(n_max,n)
248      real*8 Ba(n_max,n)
249
250      !*** local variables ***
251      integer i,j
252
253      do i=1,n
254      do j=1,n
255         Bs(i,j) = 0.5d0*(B(i,j)+B(j,i))
256         Ba(i,j) = 0.5d0*(B(i,j)-B(j,i))
257      end do
258      end do
259      return
260      end
261
262*     ***********************************
263*     *                                 *
264*     *        paw_psi_gen_UD           *
265*     *                                 *
266*     ***********************************
267      subroutine paw_psi_gen_UD(n_max,n,Bs,D,work)
268      implicit none
269      integer n_max,n
270      real*8 Bs(n_max,n)
271      real*8 D(n_max,n)
272      real*8 Work(n_max,n)
273
274      !*** local variables ***
275      integer ierr
276
277      !call eigen(n_max,n,Bs,D,D(1,2))
278      call dsyev('V','U',n,Bs,n_max, D,Work,2*n_max*n_max,ierr)
279      return
280      end
281
282
283
284
285*     ***********************************
286*     *                                 *
287*     *        paw_psi_gen_X            *
288*     *                                 *
289*     ***********************************
290      subroutine paw_psi_gen_X(n_max,n,
291     >                     X1,tmp,
292     >                     A,Ba,C,
293     >                     U,D,fnm,
294     >                     failed)
295     >
296      implicit none
297      integer n_max,n
298      real*8 X1(n_max,n)
299      real*8 tmp(*)
300      real*8 A(n_max,n)
301      real*8 Ba(n_max,n)
302      real*8 C(n_max,n)
303      real*8 U(n_max,n)
304      real*8 D(n_max,n)
305      real*8 fnm(n_max,n)
306      logical failed
307
308      !**** local variables ****
309      integer itrlmd
310      real*8  convg
311      parameter (itrlmd=40, convg=1.0d-15)
312
313      integer i,it
314      real*8  adiff
315
316      !**** external functions ****
317      integer  idamax
318      external idamax
319
320
321      !**** A = I-A ***
322      call dscal(n_max*n,(-1.0d0),A,1)
323      do i=1,n
324         A(i,i) = A(i,i) + 1.0d0
325      end do
326
327      !*** fnm = I-A ****
328      call dcopy(n_max*n,A,1,fnm,1)
329
330      !*** solve U*D*Ut*X + X*U*D*Ut = fnm for X ***
331      call paw_psi_fnm_to_X(n_max,n,fnm,U,D,tmp)
332      call dcopy(n_max*n,fnm,1,X1,1)
333
334
335      it     = 0
336      failed = .true.
337      do while (failed .and. (it.lt.itrlmd))
338        it = it + 1
339
340        !*** fnm = X*C*X ***
341        call DMMUL(n_max,n,C,X1,tmp)
342        call DMMUL(n_max,n,X1,tmp,fnm)
343
344
345        !*** fnm = Ba*X - X*C*X ***
346        call DMMUL(n_max,n,Ba,X1,tmp)
347        call DMSUB(n_max,n,tmp,fnm,fnm)
348
349
350        !*** fnm = Ba*X - X*Ba - X*C*X ***
351        call DMMUL(n_max,n,X1,Ba,tmp)
352        call DMSUB(n_max,n,fnm,tmp,fnm)
353
354        !*** fnm = I-A + Ba*X - X*Ba - X*C*X ***
355        call DMADD(n_max,n,fnm,A,fnm)
356
357
358        !*** solve U*D*Ut*X + X*U*D*Ut = fnm for X ***
359        call paw_psi_fnm_to_X(n_max,n,fnm,U,D,tmp)
360
361        call DMSUB(n_max,n,fnm,X1,tmp)
362        adiff = tmp(idamax(n_max*n,tmp,1))
363        call dcopy(n_max*n,fnm,1,X1,1)
364
365        adiff = dabs(adiff)
366
367        if (adiff.lt.convg) failed = .false.
368      end do
369
370      return
371      end
372
373
374*     ***********************************
375*     *                                 *
376*     *        paw_psi_fnm_to_X         *
377*     *                                 *
378*     ***********************************
379      subroutine paw_psi_fnm_to_X(n_max,n,
380     >                            fnm,U,D,tmp)
381      implicit none
382      integer n_max,n
383      real*8 fnm(n_max,n)
384      real*8 U(n_max,n)
385      real*8 D(n_max,n)
386      real*8 tmp(n_max,n)
387
388      !**** local variables ****
389      integer i,j
390      real*8  d2
391
392
393      !**** fnm = Ut*fnm*U ***
394      call dgemm('N','N',n,n,n,1.0d0,
395     >           fnm,n_max,
396     >           U,n_max,
397     >           0.0d0,
398     >           tmp,n_max)
399      call dgemm('T','N',n,n,n,1.0d0,
400     >           U,n_max,
401     >           tmp,n_max,
402     >           0.0d0,
403     >           fnm,n_max)
404
405
406      !**** fnm = (Ut*fnm*U)_nm/(d_n+d_m) ***
407      do j=1,n
408      do i=1,n
409        d2 = D(i,1)+D(j,1)
410        fnm(i,j) = fnm(i,j)/d2
411      end do
412      end do
413
414      !**** fnm = X = U*{(Ut*fnm*U)_nm/(d_n+d_m)}*Ut ***
415      call dgemm('N','N',n,n,n,1.0d0,
416     >           U,n_max,
417     >           fnm,n_max,
418     >           0.0d0,
419     >           tmp,n_max)
420      call dgemm('N','T',n,n,n,1.0d0,
421     >           tmp,n_max,
422     >           U,n_max,
423     >           0.0d0,
424     >           fnm,n_max)
425
426      return
427      end
428