File: readstate.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 (354 lines) | stat: -rw-r--r-- 9,565 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
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354

! 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.

!BOP
! !ROUTINE: readstate
! !INTERFACE:
subroutine readstate
! !USES:
use modmain
use moddftu
! !DESCRIPTION:
!   Reads in the charge density and other relevant variables from the file
!   {\tt STATE.OUT}. Checks for version and parameter compatibility.
!
! !REVISION HISTORY:
!   Created May 2003 (JKD)
!EOP
!BOC
implicit none
! local variables
logical spinpol_
integer iostat
integer is,ia,ias,lmmax,lm,ir,jr
integer idm,jdm,mapidm(3)
integer i1,i2,i3,j1,j2,j3,n
integer version_(3)
integer nspecies_,natoms_,lmmaxo_
integer nrmt_(maxspecies),nrmtmax_
integer nrcmt_(maxspecies),nrcmtmax_
integer ngridg_(3),ngtot_,ngvec_
integer ndmag_,nspinor_,fsmtype_,ftmtype_
integer dftu_,lmmaxdm_,xcgrad_
real(8) t1
! allocatable arrays
integer, allocatable :: mapir(:)
real(8), allocatable :: rsp_(:,:),rcmt_(:,:)
real(8), allocatable :: rfmt_(:,:,:),rfir_(:)
real(8), allocatable :: rvfmt_(:,:,:,:),rvfir_(:,:)
real(8), allocatable :: rvfcmt_(:,:,:,:),rfmt(:,:)
real(8), allocatable :: bfsmcmt_(:,:),fi(:),fo(:)
complex(8), allocatable :: vsig_(:)
complex(8), allocatable :: vmatmt_(:,:,:,:,:),vmftm_(:,:,:,:,:)
open(40,file='STATE'//trim(filext),form='UNFORMATTED',status='OLD', &
 iostat=iostat)
if (iostat.ne.0) then
  write(*,*)
  write(*,'("Error(readstate): error opening ",A)') 'STATE'//trim(filext)
  write(*,*)
  stop
end if
read(40) version_
if (version_(1).lt.2) then
  write(*,*)
  write(*,'("Error(readstate): unable to read STATE.OUT from versions earlier &
   &than 2.0.0")')
  write(*,*)
  stop
end if
if ((version(1).ne.version_(1)).or.(version(2).ne.version_(2)).or. &
    (version(3).ne.version_(3))) then
  write(*,*)
  write(*,'("Warning(readstate): different versions")')
  write(*,'(" current   : ",I3.3,".",I3.3,".",I3.3)') version
  write(*,'(" STATE.OUT : ",I3.3,".",I3.3,".",I3.3)') version_
end if
read(40) spinpol_
read(40) nspecies_
if (nspecies.ne.nspecies_) then
  write(*,*)
  write(*,'("Error(readstate): differing nspecies")')
  write(*,'(" current   : ",I4)') nspecies
  write(*,'(" STATE.OUT : ",I4)') nspecies_
  write(*,*)
  stop
end if
read(40) lmmaxo_
lmmax=min(lmmaxo,lmmaxo_)
read(40) nrmtmax_
read(40) nrcmtmax_
allocate(rsp_(nrmtmax_,nspecies))
allocate(rcmt_(nrcmtmax_,nspecies))
do is=1,nspecies
  read(40) natoms_
  if (natoms(is).ne.natoms_) then
    write(*,*)
    write(*,'("Error(readstate): differing natoms for species ",I4)') is
    write(*,'(" current   : ",I4)') natoms(is)
    write(*,'(" STATE.OUT : ",I4)') natoms_
    write(*,*)
    stop
  end if
  read(40) nrmt_(is)
  read(40) rsp_(1:nrmt_(is),is)
  read(40) nrcmt_(is)
  read(40) rcmt_(1:nrcmt_(is),is)
end do
read(40) ngridg_
read(40) ngvec_
read(40) ndmag_
if ((spinpol_).and.(ndmag_.ne.1).and.(ndmag_.ne.3)) then
  write(*,*)
  write(*,'("Error(readstate): invalid ndmag in STATE.OUT : ",I8)') ndmag_
  write(*,*)
  stop
end if
read(40) nspinor_
read(40) fsmtype_
if ((version_(1).gt.2).or.(version_(2).ge.3)) then
  read(40) ftmtype_
else
  ftmtype_=0
end if
read(40) dftu_
read(40) lmmaxdm_
if ((version_(1).gt.5).or.((version_(1).eq.5).and.(version_(2).ge.1))) then
  read(40) xcgrad_
else
  xcgrad_=0
end if
ngtot_=ngridg_(1)*ngridg_(2)*ngridg_(3)
! map from old interstitial grid to new
allocate(mapir(ngtot))
ir=0
do i3=0,ngridg(3)-1
  t1=dble(i3*ngridg_(3))/dble(ngridg(3))
  j3=modulo(nint(t1),ngridg_(3))
  do i2=0,ngridg(2)-1
    t1=dble(i2*ngridg_(2))/dble(ngridg(2))
    j2=modulo(nint(t1),ngridg_(2))
    do i1=0,ngridg(1)-1
      t1=dble(i1*ngridg_(1))/dble(ngridg(1))
      j1=modulo(nint(t1),ngridg_(1))
      ir=ir+1
      jr=j3*ngridg_(2)*ngridg_(1)+j2*ngridg_(1)+j1+1
      mapir(ir)=jr
    end do
  end do
end do
allocate(rfmt_(lmmaxo_,nrmtmax_,natmtot),rfir_(ngtot_))
allocate(rfmt(lmmaxo,nrmtmax))
n=max(nrmtmax,nrmtmax_)
allocate(fi(n),fo(n))
! read the muffin-tin density
read(40) rfmt_,rfir_
! regrid and pack the muffin-tin function
call rgfmt(rhomt)
! regrid the interstitial function
rhoir(:)=rfir_(mapir(:))
! read the Coulomb potential, regrid and pack
read(40) rfmt_,rfir_
call rgfmt(vclmt)
vclir(:)=rfir_(mapir(:))
! read the exchange-correlation potential, regrid and pack
read(40) rfmt_,rfir_
call rgfmt(vxcmt)
vxcir(:)=rfir_(mapir(:))
! read the Kohn-Sham effective potential, regrid and pack
if ((version_(1).gt.2).or.(version_(2).ge.2)) then
  read(40) rfmt_,rfir_
else
  allocate(vsig_(ngvec_))
  read(40) rfmt_,rfir_,vsig_
  deallocate(vsig_)
end if
call rgfmt(vsmt)
vsir(:)=rfir_(mapir(:))
! read the magnetisation, exchange-correlation and effective magnetic fields
if (spinpol_) then
! component map for spin-polarised case
  mapidm(:)=0
  if (ndmag.eq.ndmag_) then
    do idm=1,ndmag
      mapidm(idm)=idm
    end do
  else
    mapidm(ndmag)=ndmag_
  end if
  allocate(rvfmt_(lmmaxo_,nrmtmax_,natmtot,ndmag_))
  allocate(rvfir_(ngtot_,ndmag_))
  allocate(rvfcmt_(lmmaxo_,nrcmtmax_,natmtot,ndmag_))
  read(40) rvfmt_,rvfir_
  call rgvfmt(magmt)
  call rgvir(magir)
  read(40) rvfmt_,rvfir_
  call rgvfmt(bxcmt)
  call rgvir(bxcir)
  read(40) rvfcmt_,rvfir_
  call rgvfcmt(bsmt)
  call rgvir(bsir)
  deallocate(rvfmt_,rvfir_,rvfcmt_)
! read fixed spin moment effective fields
  if (fsmtype_.ne.0) then
    allocate(bfsmcmt_(3,natmtot))
    read(40) bfsmc
    read(40) bfsmcmt_
    if (fsmtype.ne.0) bfsmcmt(:,:)=bfsmcmt_(:,:)
! make sure that the constraining fields are perpendicular to the fixed moments
! for fixed direction calculations (Y. Kvashnin and LN)
    if (fsmtype.lt.0) then
      if (ncmag) then
        call r3vo(momfix,bfsmc)
        do is=1,nspecies
          do ia=1,natoms(is)
            ias=idxas(ia,is)
            call r3vo(mommtfix(:,ia,is),bfsmcmt(:,ias))
          end do
        end do
      else
        bfsmc(:)=0.d0
        bfsmcmt(:,:)=0.d0
      end if
    end if
    deallocate(bfsmcmt_)
  end if
end if
if (xcgrad.eq.4) then
  if (xcgrad_.eq.4) then
    read(40) rfmt_,rfir_
    call rgfmt(wxcmt)
    wxcir(:)=rfir_(mapir(:))
  else
    wxcmt(:,:)=0.d0
    wxcir(:)=0.d0
  end if
  call genws
end if
deallocate(rfmt_,rfir_,rfmt,fi,fo)
! read DFT+U potential matrix in each muffin-tin
if (((dftu.ne.0).and.(dftu_.ne.0)).or. &
    ((ftmtype.ne.0).and.(ftmtype_.ne.0))) then
  allocate(vmatmt_(lmmaxdm_,nspinor_,lmmaxdm_,nspinor_,natmtot))
  read(40) vmatmt_
  lmmax=min(lmmaxdm,lmmaxdm_)
  vmatmt(:,:,:,:,:)=0.d0
  if (nspinor.eq.nspinor_) then
    vmatmt(1:lmmax,:,1:lmmax,:,:)=vmatmt_(1:lmmax,:,1:lmmax,:,:)
  else if ((nspinor.eq.1).and.(nspinor_.eq.2)) then
    vmatmt(1:lmmax,1,1:lmmax,1,:)=0.5d0*(vmatmt_(1:lmmax,1,1:lmmax,1,:) &
     +vmatmt_(1:lmmax,2,1:lmmax,2,:))
  else
    vmatmt(1:lmmax,1,1:lmmax,1,:)=vmatmt_(1:lmmax,1,1:lmmax,1,:)
    vmatmt(1:lmmax,2,1:lmmax,2,:)=vmatmt_(1:lmmax,1,1:lmmax,1,:)
  end if
  deallocate(vmatmt_)
end if
! read fixed tensor moment potential matrix elements
if ((ftmtype.ne.0).and.(ftmtype_.ne.0)) then
  allocate(vmftm_(lmmaxdm_,nspinor_,lmmaxdm_,nspinor_,natmtot))
  read(40) vmftm_
  lmmax=min(lmmaxdm,lmmaxdm_)
  vmftm_(:,:,:,:,:)=0.d0
  if (nspinor.eq.nspinor_) then
    vmftm(1:lmmax,:,1:lmmax,:,:)=vmftm_(1:lmmax,:,1:lmmax,:,:)
  else if ((nspinor.eq.1).and.(nspinor_.eq.2)) then
    vmftm(1:lmmax,1,1:lmmax,1,:)=0.5d0*(vmftm_(1:lmmax,1,1:lmmax,1,:) &
     +vmftm_(1:lmmax,2,1:lmmax,2,:))
  else
    vmftm(1:lmmax,1,1:lmmax,1,:)=vmftm_(1:lmmax,1,1:lmmax,1,:)
    vmftm(1:lmmax,2,1:lmmax,2,:)=vmftm_(1:lmmax,1,1:lmmax,1,:)
  end if
  deallocate(vmftm_)
end if
close(40)
return

contains

subroutine rgfmt(rfmtp)
implicit none
! arguments
real(8), intent(out) :: rfmtp(npmtmax,natmtot)
do ias=1,natmtot
  is=idxis(ias)
! regrid the muffin-tin function
  do lm=1,lmmax
    fi(1:nrmt_(is))=rfmt_(lm,1:nrmt_(is),ias)
    call rfinterp(nrmt_(is),rsp_(:,is),fi,nrmt(is),rsp(:,is),fo)
    rfmt(lm,1:nrmt(is))=fo(1:nrmt(is))
  end do
  rfmt(lmmax+1:lmmaxo,1:nrmt(is))=0.d0
! pack the muffin-tin function
  call rfmtpack(.true.,nrmt(is),nrmti(is),rfmt,rfmtp(:,ias))
end do
return
end subroutine

subroutine rgvfmt(rvfmt)
implicit none
! arguments
real(8), intent(out) :: rvfmt(npmtmax,natmtot,ndmag)
do idm=1,ndmag
  jdm=mapidm(idm)
  if (jdm.eq.0) then
    rvfmt(:,:,idm)=0.d0
    cycle
  end if
  do ias=1,natmtot
    is=idxis(ias)
    do lm=1,lmmax
      fi(1:nrmt_(is))=rvfmt_(lm,1:nrmt_(is),ias,jdm)
      call rfinterp(nrmt_(is),rsp_(:,is),fi,nrmt(is),rsp(:,is),fo)
      rfmt(lm,1:nrmt(is))=fo(1:nrmt(is))
    end do
    rfmt(lmmax+1:lmmaxo,1:nrmt(is))=0.d0
    call rfmtpack(.true.,nrmt(is),nrmti(is),rfmt,rvfmt(:,ias,idm))
  end do
end do
return
end subroutine

subroutine rgvfcmt(rvfcmt)
implicit none
! arguments
real(8), intent(out) :: rvfcmt(npcmtmax,natmtot,ndmag)
do idm=1,ndmag
  jdm=mapidm(idm)
  if (jdm.eq.0) then
    rvfcmt(:,:,idm)=0.d0
    cycle
  end if
  do ias=1,natmtot
    is=idxis(ias)
    do lm=1,lmmax
      fi(1:nrcmt_(is))=rvfcmt_(lm,1:nrcmt_(is),ias,jdm)
      call rfinterp(nrcmt_(is),rcmt_(:,is),fi,nrcmt(is),rcmt(:,is),fo)
      rfmt(lm,1:nrcmt(is))=fo(1:nrcmt(is))
    end do
    rfmt(lmmax+1:lmmaxo,1:nrcmt(is))=0.d0
    call rfmtpack(.true.,nrcmt(is),nrcmti(is),rfmt,rvfcmt(:,ias,idm))
  end do
end do
return
end subroutine

subroutine rgvir(rvfir)
implicit none
! arguments
real(8), intent(out) :: rvfir(ngtot,ndmag)
do idm=1,ndmag
  jdm=mapidm(idm)
  if (jdm.eq.0) then
    rvfir(:,idm)=0.d0
    cycle
  end if
  rvfir(:,idm)=rvfir_(mapir(:),jdm)
end do
return
end subroutine

end subroutine
!EOC