File: gndstulr.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 (256 lines) | stat: -rw-r--r-- 6,980 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
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256

! Copyright (C) 2017 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine gndstulr
use modmain
use modulr
use moddftu
use modmpi
use modomp
implicit none
! local variables
logical exist
integer ik0,lp,n,ir
integer nwork,nthd
real(8) dv
! allocatable arrays
integer(8), allocatable :: lock(:)
real(8), allocatable :: v(:),work(:)
complex(8), allocatable :: evecu(:,:)
if (xctype(1).le.1) then
  write(*,*)
  write(*,'("Error(gndstulr): ultra long-range does not work with OEP")')
  write(*,*)
  stop
end if
! initialise global variables
call init0
call init1
! write the q-points to QPOINTS.OUT
call writeqpts
! read the regular Kohn-Sham potential from file
call readstate
! generate the first- and second-variational eigenvectors and eigenvalues for
! the k+kappa-point set
call genvsig
call gencore
call readfermi
call linengy
call genapwlofr
call gensocfr
call genevfsv
call occupy
! initialise the ultra long-range variables
call initulr
! initialise the long-range Kohn-Sham potential and magnetic field
call pksuinit
! size of mixing vector
n=2*nqpt*(npcmtmax*natmtot+ngtot)
if (spinpol) n=n*(1+ndmag)
! allocate mixing array
allocate(v(n))
! determine the size of the mixer work array
nwork=-1
call mixerifc(mixtype,n,v,dv,nwork,v)
allocate(work(nwork))
! initialise the mixer
iscl=0
call mixpacku(.true.,n,v)
call mixerifc(mixtype,n,v,dv,nwork,work)
! initialise the OpenMP locks
allocate(lock(nqpt))
do ir=1,nqpt
  call omp_init_lock(lock(ir))
end do
! set last self-consistent loop flag
tlast=.false.
! begin the self-consistent loop
if (mp_mpi) then
! open ULR_INFO.OUT file
  open(60,file='ULR_INFO.OUT',form='FORMATTED')
! open RMSDVS.OUT
  open(65,file='RMSDVS.OUT',form='FORMATTED')
  call writeinfou(60)
  write(60,*)
  write(60,'("+------------------------------+")')
  write(60,'("| Self-consistent loop started |")')
  write(60,'("+------------------------------+")')
end if
do iscl=1,maxscl
  if (mp_mpi) then
    write(60,*)
    write(60,'("+--------------------+")')
    write(60,'("| Loop number : ",I4," |")') iscl
    write(60,'("+--------------------+")')
  end if
  if (iscl.ge.maxscl) then
    if (mp_mpi) then
      write(60,*)
      write(60,'("Reached self-consistent loops maximum")')
    end if
    write(*,*)
    write(*,'("Warning(gndstulr): failed to reach self-consistency after ",I4,&
     &" loops")') iscl
    tlast=.true.
  end if
! reset the OpenMP thread variables
  call omp_reset
! zero the density
  rhormt(:,:,:)=0.d0
  rhorir(:,:)=0.d0
! zero the magnetisation
  if (spinpol) then
    magrmt(:,:,:,:)=0.d0
    magrir(:,:,:)=0.d0
  end if
! loop over original k-points
  call omp_hold(nkpt0/np_mpi,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(evecu) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ik0=1,nkpt0
! distribute among MPI processes
    if (mod(ik0,np_mpi).ne.lp_mpi) cycle
    allocate(evecu(nstulr,nstulr))
! solve the ultra long-range eigenvalue equation
    call eveqnulr(ik0,evecu)
! add to the density, magnetisation and current
    call rhomaguk(ik0,lock,evecu)
    deallocate(evecu)
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
  if (np_mpi.gt.1) then
! broadcast eigenvalue array to every process
    do ik0=1,nkpt0
      lp=mod(ik0,np_mpi)
      call mpi_bcast(evalu(:,ik0),nstulr,mpi_double_precision,lp,mpicom,ierror)
    end do
! add densities from each process and redistribute
    n=npcmtmax*natmtot*nqpt
    call mpi_allreduce(mpi_in_place,rhormt,n,mpi_double_precision,mpi_sum, &
     mpicom,ierror)
    n=ngtot*nqpt
    call mpi_allreduce(mpi_in_place,rhorir,n,mpi_double_precision,mpi_sum, &
     mpicom,ierror)
    if (spinpol) then
      n=npcmtmax*natmtot*ndmag*nqpt
      call mpi_allreduce(mpi_in_place,magrmt,n,mpi_double_precision,mpi_sum, &
       mpicom,ierror)
      n=ngtot*ndmag*nqpt
      call mpi_allreduce(mpi_in_place,magrir,n,mpi_double_precision,mpi_sum, &
       mpicom,ierror)
    end if
  end if
! find the occupation numbers and Fermi energy
  call occupyulr
! perform partial Fourier transform to Q-space
  call rhomagq
! compute the ultra long-range Kohn-Sham potential
  call potksulr
! pack interstitial and muffin-tin potential and field into one array
  call mixpacku(.true.,n,v)
! mix in the old potential and field with the new
  call mixerifc(mixtype,n,v,dv,nwork,work)
! make sure every MPI process has a numerically identical potential
  if (np_mpi.gt.1) then
    call mpi_bcast(v,n,mpi_double_precision,0,mpicom,ierror)
  end if
! unpack potential and field
  call mixpacku(.false.,n,v)
  if (mp_mpi) then
! write eigenvalues to file
    call writeevalu
! output energy components
    call writeengyu(60)
! output charges and moments
    call writechgu(60)
! write muffin-tin moments for each R-vector
    call writemomrmt
  end if
! exit self-consistent loop if required
  if (tlast) goto 10
! check for convergence
  if (iscl.ge.2) then
    if (mp_mpi) then
      write(60,*)
      write(60,'("RMS change in Kohn-Sham potential (target) : ",G18.10," (",&
       &G18.10,")")') dv,epspot
      flush(60)
      write(65,'(G18.10)') dv
      flush(65)
    end if
    if (dv.lt.epspot) then
      if (mp_mpi) then
        write(60,*)
        write(60,'("Convergence targets achieved")')
      end if
      tlast=.true.
    end if
  end if
  if (mp_mpi) then
! check for WRITE file
    inquire(file='WRITE',exist=exist)
    if (exist) then
      write(60,*)
      write(60,'("WRITE file exists - writing STATE_ULR.OUT")')
      call writestulr
      open(50,file='WRITE')
      close(50,status='DELETE')
    end if
! write STATE_ULR.OUT file if required
    if (nwrite.ge.1) then
      if (mod(iscl,nwrite).eq.0) then
        call writestulr
        write(60,*)
        write(60,'("Wrote STATE_ULR.OUT")')
      end if
    end if
  end if
! check for STOP file (only master process)
  if (mp_mpi) then
    inquire(file='STOP',exist=exist)
    if (exist) then
      write(60,*)
      write(60,'("STOP file exists - stopping self-consistent loop")')
      open(50,file='STOP')
      close(50,status='DELETE')
      tlast=.true.
    end if
  end if
! broadcast tlast from master process to all other processes
  call mpi_bcast(tlast,1,mpi_logical,0,mpicom,ierror)
! reset the OpenMP thread variables
  call omp_reset
end do
10 continue
if (mp_mpi) then
  write(60,*)
  write(60,'("+------------------------------+")')
  write(60,'("| Self-consistent loop stopped |")')
  write(60,'("+------------------------------+")')
  if (maxscl.gt.1) then
    call writestulr
    write(60,*)
    write(60,'("Wrote STATE_ULR.OUT")')
  end if
! close the ULR_INFO.OUT file
  close(60)
! close the RMSDVS.OUT file
  close(65)
end if
! destroy the OpenMP locks
do ir=1,nqpt
  call omp_destroy_lock(lock(ir))
end do
deallocate(lock)
deallocate(v,work)
! synchronise MPI processes
call mpi_barrier(mpicom,ierror)
return
end subroutine