File: eveqnit.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (209 lines) | stat: -rw-r--r-- 6,117 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine eveqnit(nmatp,ngp,igpig,vpl,vgpl,vgpc,apwalm,evalfv,evecfv)
use modmain
use modomp
implicit none
! arguments
integer, intent(in) :: nmatp,ngp
integer, intent(in) :: igpig(ngkmax)
real(8), intent(in) :: vpl(3)
real(8), intent(in) :: vgpl(3,ngkmax),vgpc(3,ngkmax)
complex(8), intent(in) :: apwalm(ngkmax,apwordmax,lmmaxapw,natmtot)
real(8), intent(out) :: evalfv(nstfv)
complex(8), intent(out) :: evecfv(nmatmax,nstfv)
! local variables
integer n2,ns,ist,ias
integer it,i,nthd
integer lwork,info
real(8) rmax,t1
real(8) ts1,ts0
! allocatable arrays
real(8), allocatable :: w(:),rwork(:)
complex(8), allocatable :: h(:,:),o(:,:),hv(:,:),ov(:,:)
complex(8), allocatable :: u(:,:),hu(:,:),ou(:,:)
complex(8), allocatable :: hs(:,:),os(:,:),work(:)
! external functions
complex(8) zdotc
external zdotc
n2=2*nmatp
ns=2*nstfv
if (iscl.ge.2) then
! read in the eigenvalues/vectors from file
  call getevalfv(filext,0,vpl,evalfv)
  call getevecfv(filext,0,vpl,vgpl,evecfv)
!**** error with spin-spirals

else
! initialise the eigenvectors to canonical basis vectors
  evecfv(1:nmatp,:)=0.d0
  do ist=1,nstfv
    evecfv(ist,ist)=1.d0
  end do
end if
! compute Hamiltonian and overlap matrices
call timesec(ts0)
allocate(h(nmatp,nmatp),o(nmatp,nmatp))
call omp_hold(2,nthd)
!$OMP PARALLEL SECTIONS DEFAULT(SHARED) &
!$OMP PRIVATE(i,ias) &
!$OMP NUM_THREADS(nthd)
!$OMP SECTION
! Hamiltonian
do i=1,nmatp
  h(1:i,i)=0.d0
end do
do ias=1,natmtot
  call hmlaa(ias,ngp,apwalm(:,:,:,ias),nmatp,h)
  call hmlalo(ias,ngp,apwalm(:,:,:,ias),nmatp,h)
  call hmllolo(ias,ngp,nmatp,h)
end do
call hmlistl(ngp,igpig,vgpc,nmatp,h)
!$OMP SECTION
! overlap
do i=1,nmatp
  o(1:i,i)=0.d0
end do
do ias=1,natmtot
  call olpaa(ias,ngp,apwalm(:,:,:,ias),nmatp,o)
  call olpalo(ias,ngp,apwalm(:,:,:,ias),nmatp,o)
  call olplolo(ias,ngp,nmatp,o)
end do
call olpistl(ngp,igpig,nmatp,o)
!$OMP END PARALLEL SECTIONS
call omp_free(nthd)
call timesec(ts1)
!$OMP CRITICAL(eveqnit_1)
timemat=timemat+ts1-ts0
!$OMP END CRITICAL(eveqnit_1)
call timesec(ts0)
allocate(w(ns),rwork(3*ns))
allocate(hv(nmatp,nstfv),ov(nmatp,nstfv))
allocate(u(nmatp,nstfv),hu(nmatp,nstfv),ou(nmatp,nstfv))
allocate(hs(ns,ns),os(ns,ns))
lwork=2*ns
allocate(work(lwork))
! start iteration loop
do it=1,maxitefv
  rmax=0.d0
  call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(t1) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ist=1,nstfv
! operate with O on the current eigenvector
    call zhemv('U',nmatp,zone,o,nmatp,evecfv(:,ist),1,zzero,ov(:,ist),1)
! normalise the eigenvector
    t1=dble(zdotc(nmatp,evecfv(:,ist),1,ov(:,ist),1))
    if (t1.gt.0.d0) then
      t1=1.d0/sqrt(t1)
      call dscal(n2,t1,evecfv(:,ist),1)
      call dscal(n2,t1,ov(:,ist),1)
    end if
! operate with H on the current eigenvector
    call zhemv('U',nmatp,zone,h,nmatp,evecfv(:,ist),1,zzero,hv(:,ist),1)
! estimate the eigenvalue
    t1=dble(zdotc(nmatp,evecfv(:,ist),1,hv(:,ist),1))
    if ((iscl.le.1).and.(it.eq.1)) then
      evalfv(ist)=t1
    else
      evalfv(ist)=(1.d0-befvit)*evalfv(ist)+befvit*t1
    end if
! compute the residual |u> = (H - eO)|v>
    call zcopy(nmatp,hv(:,ist),1,u(:,ist),1)
    call daxpy(n2,-evalfv(ist),ov(:,ist),1,u(:,ist),1)
! apply the overlap matrix to the residual
    call zhemv('U',nmatp,zone,o,nmatp,u(:,ist),1,zzero,ou(:,ist),1)
! compute the overlap of the residual with itself
    t1=dble(zdotc(nmatp,u(:,ist),1,ou(:,ist),1))
!$OMP CRITICAL(eveqnit_2)
    if (t1.gt.rmax) rmax=t1
!$OMP END CRITICAL(eveqnit_2)
! normalise the residual
    if (t1.gt.0.d0) then
      t1=1.d0/sqrt(t1)
      call dscal(n2,t1,u(:,ist),1)
      call dscal(n2,t1,ou(:,ist),1)
    end if
! apply the Hamiltonian matrix to the residual
    call zhemv('U',nmatp,zone,h,nmatp,u(:,ist),1,zzero,hu(:,ist),1)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
! compute the Hamiltonian and overlap matrices in the subspace formed by the
! eigenvectors and their residuals
  call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ist=1,nstfv
    call zgemv('C',nmatp,nstfv,zone,evecfv,nmatmax,hv(:,ist),1,zzero, &
     hs(1,ist),1)
    call zgemv('C',nmatp,nstfv,zone,evecfv,nmatmax,hu(:,ist),1,zzero, &
     hs(1,nstfv+ist),1)
    call zgemv('C',nmatp,nstfv,zone,u,nmatp,hu(:,ist),1,zzero, &
     hs(nstfv+1,nstfv+ist),1)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
  call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ist=1,nstfv
    call zgemv('C',nmatp,nstfv,zone,evecfv,nmatmax,ov(:,ist),1,zzero, &
     os(1,ist),1)
    call zgemv('C',nmatp,nstfv,zone,evecfv,nmatmax,ou(:,ist),1,zzero, &
     os(1,nstfv+ist),1)
    call zgemv('C',nmatp,nstfv,zone,u,nmatp,ou(:,ist),1,zzero, &
     os(nstfv+1,nstfv+ist),1)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
! solve the generalised eigenvalue problem in the subspace
  call zhegv(1,'V','U',ns,hs,ns,os,ns,w,work,lwork,rwork,info)
  if (info.ne.0) exit
! construct the new eigenvectors
  call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ist=1,nstfv
    call zgemv('N',nmatp,nstfv,zone,evecfv,nmatmax,hs(1,ist),1,zzero, &
     ov(:,ist),1)
    call zgemv('N',nmatp,nstfv,zone,u,nmatp,hs(nstfv+1,ist),1,zone,ov(:,ist),1)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
  call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ist=1,nstfv
    call zcopy(nmatp,ov(:,ist),1,evecfv(:,ist),1)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
! check for convergence
  rmax=sqrt(abs(rmax)/dble(nmatp))
  if ((it.ge.minitefv).and.(rmax.lt.epsefvit)) exit
! end iteration loop
end do
deallocate(w,rwork,h,o,hv,ov)
deallocate(u,hu,ou,hs,os,work)
call timesec(ts1)
!$OMP CRITICAL(eveqnit_3)
timefv=timefv+ts1-ts0
!$OMP END CRITICAL(eveqnit_3)
return
end subroutine