eig66_da.F90 14.2 KB
Newer Older
1
2
3
4
5
6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7
8
MODULE m_eig66_da
#include "juDFT_env.h"
9
10
11
12
13
14
15
16
17
18
19
  ! Do the IO of the eig-file in fortran direct-access
  ! The eig-file is split into two parts:
  ! eig.bas contains the basis-set information
  ! eig.vec contains the eigenvalues and the eigenvectors
  ! The record number is given by nrec=nk+(jspin-1)*nkpts
  ! each record contains:
  ! eig.bas: el,evac,ello,bkpt,wtkpt,nv,nmat,k1,k2,k3,kveclo
  ! eig.vec: ne,eig,z**
  !**: real or complex depending on calculation type
  USE m_eig66_data
  IMPLICIT NONE
20
21

CONTAINS
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
  SUBROUTINE priv_find_data(id,d)
    INTEGER,INTENT(IN)            :: id
    TYPE(t_data_DA),POINTER,INTENT(out)   :: d

    CLASS(t_data),POINTER   ::dp
    CALL eig66_find_data(dp,id)
    SELECT TYPE(dp)
    TYPE is (t_data_da)
       d=>dp
       CLASS default
       CALL judft_error("BUG: wrong datatype in eig66_da")
    END SELECT
  END SUBROUTINE priv_find_data

  SUBROUTINE open_eig(id,nmat,neig,nkpts,jspins,lmax,nlo,ntype,nlotot,create,l_dos,l_mcd,l_orb,filename,layers,nstars,ncored,nsld,nat)
    INTEGER, INTENT(IN) :: id,nmat,neig,nkpts,jspins,nlo,ntype,lmax,nlotot
    LOGICAL, INTENT(IN) :: create
39
    LOGICAL,INTENT(IN),OPTIONAL ::  l_dos,l_mcd,l_orb
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
    CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: filename
    INTEGER,INTENT(IN),OPTIONAL :: layers,nstars,ncored,nsld,nat
    !locals
    LOGICAL :: l_file
    INTEGER :: i1,recl_z,recl_eig,recl_dos
    REAL    :: r1,r3(3)
    COMPLEX :: c1
    TYPE(t_data_DA),POINTER:: d

    CALL priv_find_data(id,d)

    IF (PRESENT(filename)) d%fname=filename
    CALL eig66_data_storedefault(d,jspins,nkpts,nmat,neig,lmax,nlotot,nlo,ntype,l_dos,l_mcd,l_orb)

    !Allocate the storage for the DATA always read/write
    ALLOCATE(d%el_s(0:lmax,ntype),d%ello_s(nlo,ntype),d%evac_s(2))
    ALLOCATE(d%kvec_s(nmat,3),d%kveclo_s(nlotot))
    !Calculate the record length
    INQUIRE(IOLENGTH=recl_eig) d%el_s,d%evac_s,d%ello_s,r3,r1,i1,i1,d%kvec_s,d%kveclo_s
    d%recl_bas=recl_eig
    INQUIRE(IOLENGTH=recl_eig) r1
    recl_eig=recl_eig*(neig+2) ! add a 2 for integer 'neig'
62
#if ( defined(CPP_INVERSION) && !defined(CPP_SOC) )
63
     INQUIRE(IOLENGTH=recl_z) r1
64
#else
65
    INQUIRE(IOLENGTH=recl_z) c1
66
#endif
67
68
69
70
    recl_z=recl_z*nmat*neig

    d%recl_vec=recl_eig+recl_z

71
    IF (d%l_dos) THEN
72
73
74
75
76
77
       IF (.NOT.(PRESENT(layers).AND.PRESENT(nstars).AND.PRESENT(ncored).AND.PRESENT(nsld).AND.PRESENT(nat))) &
            CALL judft_error("BUG:Could not open file for DOS-data",calledby="eig66_da")
       INQUIRE(IOLENGTH=i1) i1
       recl_dos=i1*2*neig !ksym&jsym
       INQUIRE(IOLENGTH=i1) r1
       recl_dos=recl_dos+i1*3*neig !qvac&qis
78
79
       recl_dos=recl_dos+i1*4*ntype*neig !qal
       recl_dos=recl_dos+i1*neig*2*max(1,layers) !qvlay
80
81
82
83
84
       IF (l_orb) THEN
          recl_dos=recl_dos+i1*2*nsld*neig !qintsl,qmtsl
          recl_dos=recl_dos+i1*24*neig*nat !qmtp,orbcomp
       ENDIF
       INQUIRE(IOLENGTH=i1) c1
85
       recl_dos=recl_dos+i1*nstars*neig*max(1,layers)*2 !qstars
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
       IF (l_mcd) recl_dos=recl_dos+i1*3*ntype*ncored*neig !mcd
    ELSE
       recl_dos=-1
    ENDIF
    d%recl_dos=recl_dos


    IF (create) THEN
       d%file_io_id_bas=priv_free_uid()
       INQUIRE(file=TRIM(d%fname)//".bas",opened=l_file)
       DO WHILE(l_file)
          PRINT *,"eig66_open_da:",d%fname," in use"
          d%fname=TRIM(d%fname)//"6"
          INQUIRE(file=TRIM(d%fname)//".bas",opened=l_file)
       ENDDO
       OPEN(d%file_io_id_bas,FILE=TRIM(d%fname)//".bas",ACCESS='direct',FORM='unformatted',RECL=d%recl_bas,STATUS='unknown')
       d%file_io_id_vec=priv_free_uid()
       OPEN(d%file_io_id_vec,FILE=TRIM(d%fname)//".vec",ACCESS='direct',FORM='unformatted',RECL=d%recl_vec,STATUS='unknown')
       IF(d%recl_dos>0) THEN
          d%file_io_id_dos=priv_free_uid()
          OPEN(d%file_io_id_dos,FILE=TRIM(d%fname)//".dos",ACCESS='direct',FORM='unformatted',RECL=d%recl_dos,STATUS='unknown')
       ENDIF

    ELSE
       d%file_io_id_bas=priv_free_uid()
       OPEN(d%file_io_id_bas,FILE=TRIM(d%fname)//".bas",ACCESS='direct',FORM='unformatted',RECL=d%recl_bas,STATUS='old')
       d%file_io_id_vec=priv_free_uid()
       OPEN(d%file_io_id_vec,FILE=TRIM(d%fname)//".vec",ACCESS='direct',FORM='unformatted',RECL=d%recl_vec,STATUS='old')
       IF(d%recl_dos>0) THEN
          d%file_io_id_dos=priv_free_uid()
          OPEN(d%file_io_id_dos,FILE=TRIM(d%fname)//".dos",ACCESS='direct',FORM='unformatted',RECL=d%recl_dos,STATUS='old')
       ENDIF
    ENDIF
  CONTAINS
    INTEGER FUNCTION priv_free_uid() RESULT(uid)
      IMPLICIT NONE
      LOGICAL::used
      used=.TRUE.
      uid=665
      DO WHILE(used)
         uid=uid+1
         INQUIRE(UNIT=uid,OPENED=used)
      END DO
    END FUNCTION priv_free_uid
  END SUBROUTINE open_eig
  SUBROUTINE close_eig(id,filename)
    INTEGER,INTENT(IN)::id
    CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: filename
    TYPE(t_data_DA),POINTER:: d

    CALL priv_find_data(id,d)

    DEALLOCATE(d%el_s,d%ello_s,d%evac_s,d%kvec_s,d%kveclo_s)
    CLOSE(d%file_io_id_bas)
    CLOSE(d%file_io_id_vec)
    d%recl_vec=0
    d%recl_bas=0

    !If a filename was given and the name is not the current filename then rename
    IF (PRESENT(filename)) THEN
       IF (filename.NE.d%fname) THEN
          CALL system("mv "//TRIM(d%fname)//".bas "//TRIM(filename)//".bas")
          CALL system("mv "//TRIM(d%fname)//".vec "//TRIM(filename)//".vec")
       ENDIF
    ENDIF
    d%fname="eig"
    CALL eig66_remove_data(id)
  END SUBROUTINE close_eig
  SUBROUTINE read_eig(id,nk,jspin,nv,nmat,k1,k2,k3,bk,wk,neig,eig,el,ello,evac,kveclo,n_start,n_end,z)
    IMPLICIT NONE
    INTEGER, INTENT(IN)            :: id,nk,jspin
    INTEGER, INTENT(OUT),OPTIONAL  :: nv,nmat
    INTEGER, INTENT(OUT),OPTIONAL  :: neig
    REAL,    INTENT(OUT),OPTIONAL  :: eig(:)
    INTEGER, INTENT(OUT),OPTIONAL  :: k1(:),k2(:),k3(:),kveclo(:)
    REAL,    INTENT(OUT),OPTIONAL  :: evac(:),ello(:,:),el(:,:)
    REAL,    INTENT(OUT),OPTIONAL  :: bk(:),wk
    INTEGER, INTENT(IN),OPTIONAL   :: n_start,n_end
    CLASS(*),TARGET,INTENT(OUT),OPTIONAL  :: z(:,:)

    !Local variables
    INTEGER:: nv_s,nmat_s,n,nrec,neig_s
    REAL   :: bkpt(3),wtkpt
    REAL,ALLOCATABLE::eig_s(:),zr_s(:,:)
    COMPLEX,ALLOCATABLE::zc_s(:,:)
    REAL,POINTER :: zr(:,:)
    COMPLEX,POINTER:: zc(:,:)
    TYPE(t_data_DA),POINTER:: d



    CALL priv_find_data(id,d)
    ! check if io is performed correctly
    IF (PRESENT(n_start)) THEN
       IF (n_start/=1) &
            CALL juDFT_error("In direct access mode only all eigenstates can be read")
    ENDIF

    nrec=nk+(jspin-1)*d%nkpts
    IF (PRESENT(el).OR.PRESENT(ello).OR.PRESENT(evac).OR.PRESENT(bk).OR.PRESENT(wk).OR.&
         PRESENT(nv).OR.PRESENT(nmat).OR.PRESENT(k1).OR.PRESENT(k2).OR.PRESENT(k3).OR.&
         PRESENT(kveclo)) THEN
       !IO of basis-set information
       READ(d%file_io_id_bas,REC=nrec) nmat_s,d%el_s,d%evac_s,d%ello_s,bkpt,wtkpt,nv_s,d%kvec_s,d%kveclo_s
       IF (PRESENT(el)) el=d%el_s
       IF (PRESENT(evac)) evac=d%evac_s
       IF (PRESENT(ello)) ello=d%ello_s
       IF (PRESENT(bk)) bk=bkpt
       IF (PRESENT(wk)) wk=wtkpt
       IF (PRESENT(nv)) nv=nv_s
       IF (PRESENT(nmat)) nmat=nmat_s
       IF (PRESENT(k1)) k1=d%kvec_s(:,1)
       IF (PRESENT(k2)) k2=d%kvec_s(:,2)
       IF (PRESENT(k3)) k3=d%kvec_s(:,3)
       IF (PRESENT(kveclo)) kveclo=d%kveclo_s
    ENDIF

    IF (.NOT.(PRESENT(eig).OR.PRESENT(neig).OR.PRESENT(z))) RETURN
204
    READ(d%file_io_id_vec,REC=nrec) neig_s
205
    IF (PRESENT(neig)) neig=neig_s
206
207
    IF (.NOT.(PRESENT(eig).OR.PRESENT(z))) RETURN
    ALLOCATE(eig_s(neig_s))
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
    IF (PRESENT(z)) THEN
       SELECT TYPE(z)
       TYPE IS (REAL)
          INQUIRE(IOLENGTH=n) neig_s,eig_s,REAL(z)
          IF (n>d%recl_vec) CALL juDFT_error("BUG: Too long record")
          zr=>z
          READ(d%file_io_id_vec,REC=nrec) neig_s,eig_s,zr
       TYPE IS (COMPLEX)
          INQUIRE(IOLENGTH=n) neig_s,eig_s,CMPLX(z)
          IF (n>d%recl_vec) CALL juDFT_error("BUG: Too long record")
          zc=>z
          READ(d%file_io_id_vec,REC=nrec) neig_s,eig_s,zc
       END SELECT
    ELSE
       INQUIRE(IOLENGTH=n) neig_s,eig_s
       IF (n>d%recl_vec) CALL juDFT_error("BUG: Too long record")
       READ(d%file_io_id_vec,REC=nrec) neig_s,eig_s
    ENDIF
226
    IF (PRESENT(eig)) eig(:min(size(eig),neig_s))=eig_s(:min(size(eig),neig_s))
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

  END SUBROUTINE read_eig

  SUBROUTINE write_eig(id,nk,jspin,neig,neig_total,nv,nmat,k1,k2,k3,bk,wk, &
       eig,el,ello,evac,nlotot,kveclo,n_size,n_rank,z)
    INTEGER, INTENT(IN)          :: id,nk,jspin
    INTEGER, INTENT(IN),OPTIONAL :: n_size,n_rank
    REAL,    INTENT(IN),OPTIONAL :: wk
    INTEGER, INTENT(IN),OPTIONAL :: neig,nv,nmat,nlotot,neig_total
    INTEGER, INTENT(IN),OPTIONAL :: k1(:),k2(:),k3(:),kveclo(:)
    REAL,    INTENT(IN),OPTIONAL :: bk(3),eig(:),el(:,:)
    REAL,    INTENT(IN),OPTIONAL :: evac(:),ello(:,:)
    CLASS(*),INTENT(IN),OPTIONAL :: z(:,:)

    INTEGER:: nrec,r_len
    INTEGER:: nv_s,nmat_s
    REAL   :: bkpt(3),wtkpt
    TYPE(t_data_DA),POINTER:: d

    CALL priv_find_data(id,d)
    !This mode requires all data to be written at once!!

    IF (PRESENT(n_size).AND.PRESENT(n_rank)) THEN
       IF (n_size/=1.OR.n_rank/=0) &
            CALL juDFT_error("Direct Access IO not possible in eigenvalue parallel code")
    ENDIF
    !check record length
    !INQUIRE(iolength=r_len) nmat,el,evac,ello,bk,wk,nv,d%kvec_s,kveclo
    !if (r_len>recl_bas) call juDFT_error("BUG: too long record")

    !Now it is time for the IO :-)
    nrec=nk+(jspin-1)*d%nkpts
    IF (PRESENT(nmat).AND..NOT.PRESENT(el)) THEN
       !IO of basis-set information
       READ(d%file_io_id_bas,REC=nrec,ERR=88) nmat_s,d%el_s,d%evac_s,d%ello_s,bkpt,wtkpt,nv_s,d%kvec_s,d%kveclo_s
88     WRITE(d%file_io_id_bas,REC=nrec)       nmat  ,d%el_s,d%evac_s,d%ello_s,bkpt,wtkpt,nv_s,d%kvec_s,d%kveclo_s
       IF (PRESENT(wk).OR.PRESENT(nv).OR.PRESENT(nlotot) &
            .OR.PRESENT(k1).OR.PRESENT(k2).OR.PRESENT(k3).OR.PRESENT(kveclo).OR.&
            PRESENT(bk).OR.PRESENT(ello).OR.PRESENT(evac)) THEN
          CALL juDFT_error("BUG:Direct access IO of eig-file only with all scalar data")
       ENDIF
    ELSE
       IF (.NOT.(PRESENT(wk).AND.PRESENT(nv).AND.PRESENT(nmat).AND.PRESENT(nlotot) &
            .AND.PRESENT(k1).AND.PRESENT(k2).AND.PRESENT(k3).AND.PRESENT(kveclo).AND.&
            PRESENT(bk).AND.PRESENT(el).AND.PRESENT(ello).AND.PRESENT(evac))) THEN
          CALL juDFT_error("BUG:Direct access IO of eig-file only with all data")
       ENDIF
       d%kvec_s(:,1)=k1
       d%kvec_s(:,2)=k2
       d%kvec_s(:,3)=k3
       IF ((SIZE(el).NE.SIZE(d%el_s)).OR.(SIZE(ello).NE.SIZE(d%ello_s).OR.(SIZE(evac).NE.SIZE(d%evac_s)))) THEN
          WRITE(*,*) SHAPE(el),SHAPE(d%el_s)
          WRITE(*,*) SHAPE(ello),SHAPE(d%ello_s)
          WRITE(*,*) SHAPE(evac),SHAPE(d%evac_s)
          CALL juDFT_error("Mismatch of sizes")
       ENDIF
       WRITE(d%file_io_id_bas,REC=nrec) nmat,el,evac,ello,bk,wk,nv,d%kvec_s,kveclo

    ENDIF
    IF (PRESENT(neig).AND.PRESENT(neig_total)) THEN
       IF (neig.NE.neig_total) THEN
          CALL juDFT_error("Neig and neig_total have to be equal in DA mode",calledby="eig66_da")
       ENDIF
    ENDIF
    IF (.NOT.PRESENT(eig).OR..NOT.PRESENT(neig)) RETURN
    !Now the IO of the eigenvalues/vectors
    IF (PRESENT(z)) THEN
       SELECT TYPE(z)
       TYPE IS (REAL)
          INQUIRE(IOLENGTH=r_len) neig,eig,REAL(z)
          IF (r_len>d%recl_vec) CALL juDFT_error("BUG: too long record")
          WRITE(d%file_io_id_vec,REC=nrec) neig,eig,REAL(z)
       TYPE IS (COMPLEX)
300
          INQUIRE(IOLENGTH=r_len) neig,eig(:neig),CMPLX(z)
301
          IF (r_len>d%recl_vec) CALL juDFT_error("BUG: too long record")
302
          WRITE(d%file_io_id_vec,REC=nrec) neig,eig(:neig),CMPLX(z)
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
355
356
357
358
359
360
361
       END SELECT
    ELSE
       INQUIRE(IOLENGTH=r_len) neig,eig
       IF (r_len>d%recl_vec) CALL juDFT_error("BUG: too long record")
       WRITE(d%file_io_id_vec,REC=nrec) neig,eig
    ENDIF

  END SUBROUTINE write_eig


  SUBROUTINE write_dos(id,nk,jspin,qal,qvac,qis,qvlay,qstars,ksym,jsym,mcd,qintsl,qmtsl,qmtp,orbcomp)
    IMPLICIT NONE
    INTEGER, INTENT(IN)          :: id,nk,jspin
    REAL,INTENT(IN)              :: qal(:,:,:),qvac(:,:),qis(:),qvlay(:,:,:)
    COMPLEX,INTENT(IN)           :: qstars(:,:,:,:)
    INTEGER,INTENT(IN)           :: ksym(:),jsym(:)
    REAL,INTENT(IN),OPTIONAL  :: mcd(:,:,:)
    REAL,INTENT(IN),OPTIONAL     :: qintsl(:,:),qmtsl(:,:),qmtp(:,:),orbcomp(:,:,:)
    TYPE(t_data_DA),POINTER:: d
    INTEGER:: nrec
    CALL priv_find_data(id,d)

    nrec=nk+(jspin-1)*d%nkpts

    IF (d%l_orb.AND.PRESENT(qmtsl)) THEN
       IF (d%l_mcd) CPP_error("mcd & orbital decomposition not implemented in IO")
       WRITE(d%file_io_id_dos,REC=nrec) qal,qvac,qis,qvlay,qstars,ksym,jsym,qintsl,qmtsl,qmtp,orbcomp
    ELSEIF(d%l_mcd.AND.PRESENT(mcd)) THEN
       WRITE(d%file_io_id_dos,REC=nrec) qal,qvac,qis,qvlay,qstars,ksym,jsym,mcd
    ELSE
       WRITE(d%file_io_id_dos,REC=nrec) qal,qvac,qis,qvlay,qstars,ksym,jsym
    END IF
  END SUBROUTINE write_dos

  SUBROUTINE read_dos(id,nk,jspin,qal,qvac,qis,qvlay,qstars,ksym,jsym,mcd,qintsl,qmtsl,qmtp,orbcomp)
    IMPLICIT NONE
    INTEGER, INTENT(IN)          :: id,nk,jspin
    REAL,INTENT(OUT)              :: qal(:,:,:),qvac(:,:),qis(:),qvlay(:,:,:)
    COMPLEX,INTENT(OUT)           :: qstars(:,:,:,:)
    INTEGER,INTENT(OUT)           :: ksym(:),jsym(:)
    REAL,INTENT(OUT),OPTIONAL  :: mcd(:,:,:)
    REAL,INTENT(OUT),OPTIONAL     :: qintsl(:,:),qmtsl(:,:),qmtp(:,:),orbcomp(:,:,:)
    TYPE(t_data_DA),POINTER:: d
    INTEGER:: nrec
    CALL priv_find_data(id,d)

    nrec=nk+(jspin-1)*d%nkpts

    IF (d%l_orb.AND.PRESENT(qmtsl)) THEN
       IF (d%l_mcd) CPP_error("mcd & orbital decomposition not implemented in IO")
       READ(d%file_io_id_dos,REC=nrec) qal,qvac,qis,qvlay,qstars,ksym,jsym,qintsl,qmtsl,qmtp,orbcomp
    ELSEIF(d%l_mcd.AND.PRESENT(mcd)) THEN
       READ(d%file_io_id_dos,REC=nrec) qal,qvac,qis,qvlay,qstars,ksym,jsym,mcd
    ELSE
       READ(d%file_io_id_dos,REC=nrec) qal,qvac,qis,qvlay,qstars,ksym,jsym
    END IF
  END SUBROUTINE read_dos


362
END MODULE m_eig66_da