Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2017 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
8 : !!!
9 : !!! This module is a wrapper for the charge density I/O
10 : !!!
11 : !!! GM'17
12 : !!!
13 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
14 :
15 : MODULE m_cdn_io
16 : #ifdef CPP_MPI
17 : use mpi
18 : #endif
19 : USE m_types
20 : USE m_juDFT
21 : USE m_loddop
22 : USE m_wrtdop
23 : USE m_cdnpot_io_hdf
24 : USE m_cdnpot_io_common
25 : USE m_constants
26 : #ifdef CPP_HDF
27 : USE hdf5
28 : #endif
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 : PUBLIC printDensityFileInfo
33 : PUBLIC readDensity, writeDensity
34 : PUBLIC isDensityFilePresent, isCoreDensityPresent
35 : PUBLIC readCoreDensity, writeCoreDensity
36 : PUBLIC readStars, writeStars
37 : PUBLIC readStepfunction, writeStepfunction
38 : PUBLIC setStartingDensity, readPrevEFermi, deleteDensities
39 : PUBLIC storeStructureIfNew,transform_by_moving_atoms, readPrevmmpDistances
40 : PUBLIC getIOMode
41 : PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
42 : PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
43 : PUBLIC CDN_ARCHIVE_TYPE_CDN_const, CDN_ARCHIVE_TYPE_FFN_const
44 :
45 : INTEGER, PARAMETER :: CDN_INPUT_DEN_const = 1
46 : INTEGER, PARAMETER :: CDN_OUTPUT_DEN_const = 2
47 :
48 : INTEGER, PARAMETER :: CDN_ARCHIVE_TYPE_CDN1_const = 1
49 : INTEGER, PARAMETER :: CDN_ARCHIVE_TYPE_NOCO_const = 2
50 : INTEGER, PARAMETER :: CDN_ARCHIVE_TYPE_CDN_const = 3
51 : INTEGER, PARAMETER :: CDN_ARCHIVE_TYPE_FFN_const = 4
52 :
53 : INTEGER, PARAMETER :: CDN_DIRECT_MODE = 1
54 : INTEGER, PARAMETER :: CDN_STREAM_MODE = 2
55 : INTEGER, PARAMETER :: CDN_HDF5_MODE = 3
56 :
57 : CONTAINS
58 :
59 0 : SUBROUTINE printDensityFileInfo()
60 :
61 : INTEGER :: mode, i
62 : LOGICAL :: l_exist
63 :
64 : #ifdef CPP_HDF
65 : INTEGER(HID_T) :: fileID
66 : #endif
67 : INTEGER :: currentStarsIndex,currentLatharmsIndex
68 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
69 : INTEGER :: readDensityIndex, lastDensityIndex
70 : CHARACTER(LEN=30) :: archiveName
71 :
72 : INTEGER :: dateTemp, timeTemp
73 : INTEGER :: iterTemp, starsIndexTemp, latharmsIndexTemp
74 : INTEGER :: structureIndexTemp,stepfunctionIndexTemp
75 : INTEGER :: previousDensityIndex, jspinsTemp
76 : REAL :: fermiEnergyTemp, distanceTemp
77 : LOGICAL :: l_qfixTemp
78 : CHARACTER(LEN=10) :: dateString
79 : CHARACTER(LEN=10) :: timeString
80 : CHARACTER(LEN=19) :: timeStampString
81 : CHARACTER(LEN=15) :: distanceString
82 :
83 0 : CALL getIOMode(mode)
84 :
85 0 : WRITE(*,*) 'Available densities info:'
86 0 : WRITE(*,*)
87 :
88 0 : IF(mode.EQ.CDN_HDF5_MODE) THEN
89 : #ifdef CPP_HDF
90 0 : INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
91 0 : IF (l_exist) THEN
92 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
93 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
94 0 : WRITE(*,*) 'densityIndex iteration prevDensity prevDistance timeStamp'
95 0 : DO i = 1, lastDensityIndex
96 0 : archiveName = ''
97 0 : WRITE(archiveName,'(a,i0)') '/cdn-', i
98 :
99 0 : l_exist = isDensityEntryPresentHDF(fileID,archiveName,DENSITY_TYPE_UNDEFINED_const)
100 0 : IF(.NOT.l_exist) THEN
101 : CYCLE
102 : END IF
103 :
104 : CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_UNDEFINED_const,&
105 : iterTemp, starsIndexTemp, latharmsIndexTemp, structureIndexTemp,&
106 : stepfunctionIndexTemp,previousDensityIndex, jspinsTemp,&
107 0 : dateTemp, timeTemp, distanceTemp, fermiEnergyTemp, l_qfixTemp)
108 :
109 0 : WRITE(dateString,'(i8)') dateTemp
110 0 : WRITE(timeString,'(i6)') timeTemp
111 :
112 0 : distanceString = ''
113 0 : IF (distanceTemp.GE.-1e-10) THEN
114 0 : WRITE(distanceString,'(f15.8)') distanceTemp
115 : END IF
116 :
117 : WRITE(timeStampString,'(a4,a1,a2,a1,a2,1x,a2,a1,a2,a1,a2)') &
118 0 : dateString(1:4),'/',dateString(5:6),'/',dateString(7:8),&
119 0 : timeString(1:2),':',timeString(3:4),':',timeString(5:6)
120 :
121 0 : WRITE(*,'(1x,i7,6x,i7,7x,i7,4x,a15,3x,a)') i, iterTemp, previousDensityIndex, distanceString,&
122 0 : TRIM(ADJUSTL(timeStampString))
123 : END DO
124 0 : CALL closeCDNPOT_HDF(fileID)
125 : ELSE
126 0 : WRITE(*,'(a)') "No cdn.hdf file found. Density file info is not available."
127 : END IF
128 : #else
129 : WRITE(*,'(a)') "Fleur is not compiled with HDF5 support. Density file info is not available."
130 : #endif
131 : ELSE
132 0 : WRITE(*,'(a)') "Density file info is only available if '-hdf_cdn' switch is used."
133 : END IF
134 :
135 0 : END SUBROUTINE printDensityFileInfo
136 :
137 :
138 80 : SUBROUTINE readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,inOrOutCDN,&
139 : relCdnIndex,fermiEnergy,lastDistance,l_qfix,den,inFilename,denIm)
140 :
141 : TYPE(t_stars),INTENT(IN) :: stars
142 : TYPE(t_vacuum),INTENT(IN) :: vacuum
143 : TYPE(t_atoms),INTENT(IN) :: atoms
144 : TYPE(t_cell), INTENT(IN) :: cell
145 : TYPE(t_sphhar),INTENT(IN) :: sphhar
146 : TYPE(t_input),INTENT(IN) :: input
147 : TYPE(t_sym),INTENT(IN) :: sym
148 : TYPE(t_noco),INTENT(IN) :: noco
149 :
150 :
151 : TYPE(t_potden),INTENT(INOUT) :: den
152 :
153 : INTEGER, INTENT (IN) :: inOrOutCDN
154 : INTEGER, INTENT (IN) :: relCdnIndex
155 : INTEGER, INTENT (IN) :: archiveType
156 : REAL, INTENT (OUT) :: fermiEnergy, lastDistance
157 : LOGICAL, INTENT (OUT) :: l_qfix
158 :
159 : CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: inFilename
160 :
161 : TYPE(t_potden), OPTIONAL, INTENT(INOUT) :: denIm
162 :
163 : ! local variables
164 : INTEGER :: mode, datend, k, i, iVac, j, iUnit, l, numLines, ioStatus, iofl
165 : LOGICAL :: l_exist, l_rhomatFile, l_DimChange
166 : CHARACTER(LEN=30) :: filename
167 :
168 : #ifdef CPP_HDF
169 : INTEGER(HID_T) :: fileID
170 : #endif
171 : INTEGER :: currentStarsIndex,currentLatharmsIndex
172 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
173 : INTEGER :: readDensityIndex, lastDensityIndex
174 : INTEGER :: previousDensityIndex, densityType
175 : CHARACTER(LEN=30) :: archiveName
176 : TYPE(t_cell) :: cellTemp
177 :
178 80 : COMPLEX, ALLOCATABLE :: cdomvz(:,:)
179 :
180 80 : fermiEnergy = 0.0
181 80 : lastDistance = -1.0
182 80 : l_qfix = .FALSE.
183 :
184 80 : CALL getIOMode(mode)
185 :
186 : #ifndef CPP_HDF
187 : filename = 'cdn.hdf'
188 : IF(PRESENT(inFilename)) filename = TRIM(ADJUSTL(inFilename))//'.hdf'
189 : INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
190 : IF (l_exist) THEN
191 : CALL juDFT_warn('Fleur not compiled for HDF5, but '//TRIM(ADJUSTL(filename))//' present',calledby='readDensity')
192 : END IF
193 : #endif
194 :
195 : IF(mode.EQ.CDN_HDF5_MODE) THEN
196 : #ifdef CPP_HDF
197 :
198 80 : densityType = 0
199 80 : archiveName = ''
200 :
201 80 : filename = 'cdn.hdf'
202 80 : IF(PRESENT(inFilename)) filename = TRIM(ADJUSTL(inFilename))//'.hdf'
203 :
204 80 : INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
205 80 : IF (l_exist) THEN
206 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
207 160 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
208 :
209 80 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
210 0 : archiveName = 'cdn'
211 : ELSE
212 80 : WRITE(archiveName,'(a,i0)') '/cdn-', readDensityIndex
213 : END IF
214 :
215 160 : SELECT CASE (inOrOutCDN)
216 : CASE (CDN_INPUT_DEN_const)
217 80 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_FFN_const) THEN
218 9 : densityType = DENSITY_TYPE_FFN_IN_const
219 71 : ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
220 17 : densityType = DENSITY_TYPE_NOCO_IN_const
221 : ELSE
222 54 : densityType = DENSITY_TYPE_IN_const
223 : END IF
224 : CASE (CDN_OUTPUT_DEN_const)
225 0 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_FFN_const) THEN
226 0 : densityType = DENSITY_TYPE_FFN_OUT_const
227 0 : ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
228 0 : densityType = DENSITY_TYPE_NOCO_OUT_const
229 : ELSE
230 0 : densityType = DENSITY_TYPE_OUT_const
231 : END IF
232 : CASE DEFAULT
233 0 : WRITE(*,*) 'inOrOutCDN = ', inOrOutCDN
234 80 : CALL juDFT_error("Invalid inOrOutCDN selected.",calledby ="readDensity")
235 : END SELECT
236 80 : l_exist = isDensityEntryPresentHDF(fileID,archiveName,densityType)
237 80 : CALL closeCDNPOT_HDF(fileID)
238 : END IF
239 :
240 80 : IF (l_exist) THEN
241 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
242 160 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
243 :
244 80 : IF (PRESENT(denIm)) THEN
245 : CALL readDensityHDF(fileID, input, stars, sphhar, atoms, vacuum, archiveName, densityType,&
246 0 : fermiEnergy,lastDistance,l_qfix,l_DimChange,den,denIm)
247 : ELSE
248 : CALL readDensityHDF(fileID, input, stars, sphhar, atoms, vacuum, archiveName, densityType,&
249 80 : fermiEnergy,lastDistance,l_qfix,l_DimChange,den)
250 : END IF
251 :
252 80 : CALL closeCDNPOT_HDF(fileID)
253 :
254 80 : IF(l_DimChange) THEN
255 6 : IF (PRESENT(denIm)) THEN
256 : CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,inOrOutCDN,&
257 0 : 1,-1.0,fermiEnergy,-1.0,-1.0,l_qfix,den,denIm=denIm)
258 : ELSE
259 : CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,inOrOutCDN,&
260 6 : 1,-1.0,fermiEnergy,-1.0,-1.0,l_qfix,den)
261 : END IF
262 : END IF
263 : ELSE
264 0 : INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
265 0 : IF(l_exist) THEN
266 0 : WRITE(*,*) 'densityType is ', densityType
267 0 : WRITE(*,*) 'Relevant density entry is '//TRIM(ADJUSTL(archiveName))//'.'
268 0 : WRITE(*,*) 'Entry not found in '//TRIM(ADJUSTL(filename))//'.'
269 : ELSE
270 0 : WRITE(*,*) TRIM(ADJUSTL(filename))//' file not found.'
271 : END IF
272 0 : WRITE(*,*) 'Falling back to stream access.'
273 : mode = CDN_STREAM_MODE
274 : END IF
275 : #endif
276 : END IF
277 :
278 : IF(mode.EQ.CDN_STREAM_MODE) THEN
279 0 : INQUIRE(FILE='cdn.str',EXIST=l_exist)
280 0 : IF (l_exist) THEN
281 : !load density from cdn.str and exit subroutine
282 :
283 : ELSE
284 0 : WRITE(*,*) 'cdn.str file not found.'
285 0 : WRITE(*,*) 'Falling back to direct access file cdn1.'
286 : mode = CDN_DIRECT_MODE
287 : END IF
288 : END IF
289 :
290 80 : IF (mode.EQ.CDN_DIRECT_MODE) THEN
291 0 : if (any(noco%l_unrestrictMT)) call judft_error("the mtNocoPot switch requires HDF5 for the charge density IO")
292 0 : filename = 'cdn1'
293 0 : l_rhomatFile = .FALSE.
294 0 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
295 0 : INQUIRE(file="rhomat_inp",EXIST=l_exist)
296 0 : IF (l_exist) filename = 'rhomat_inp'
297 0 : IF (inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) THEN
298 0 : INQUIRE(file="rhomat_out",EXIST=l_exist)
299 0 : IF (l_exist) filename = 'rhomat_out'
300 : END IF
301 0 : IF(l_exist) l_rhomatFile = .TRUE.
302 : END IF
303 0 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
304 0 : filename = 'cdn'
305 : END IF
306 :
307 0 : IF(PRESENT(inFilename)) filename = inFilename
308 :
309 0 : INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
310 0 : IF(.NOT.l_exist) THEN
311 0 : CALL juDFT_error("charge density file "//TRIM(ADJUSTL(filename))//" missing",calledby ="readDensity")
312 : END IF
313 :
314 0 : iUnit = 93
315 0 : OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='old')
316 :
317 0 : IF ((inOrOutCDN.EQ.CDN_OUTPUT_DEN_const).AND.(archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const)) THEN
318 : ! call loddop to move the file position to the output density
319 : CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
320 0 : iUnit,den%iter,den%mt,den%pw,den%vac)
321 : END IF
322 :
323 : ! read in the density
324 : CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
325 0 : iUnit,den%iter,den%mt,den%pw,den%vac)
326 :
327 : ! read in additional data if l_noco and data is present
328 0 : IF ((archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const).AND.l_rhomatFile) THEN
329 0 : READ (iUnit,iostat=datend) (den%pw(k,3),k=1,stars%ng3)
330 0 : IF (datend == 0) THEN
331 0 : IF (input%film) THEN
332 0 : ALLOCATE(cdomvz(vacuum%nmz,vacuum%nvac))
333 0 : READ (iUnit) ((cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
334 0 : DO iVac = 1, vacuum%nvac
335 0 : DO i = 1, vacuum%nmz
336 0 : den%vac(i,1,iVac,3) = REAL(cdomvz(i,iVac))+ImagUnit*AIMAG(cdomvz(i,iVac))
337 : END DO
338 : END DO
339 0 : DEALLOCATE(cdomvz)
340 0 : READ (iUnit) (((den%vac(i,j,iVac,3),i=1,vacuum%nmzxy),j=2,stars%ng2), iVac=1,vacuum%nvac)
341 : END IF
342 : ELSE
343 : ! (datend < 0) => no off-diagonal magnetisation stored
344 : ! in "rhomat_inp"
345 0 : IF (datend > 0) THEN
346 0 : WRITE(*,*) 'datend = ', datend
347 0 : CALL juDFT_error("density file has illegal state.",calledby ="readDensity")
348 : END IF
349 0 : den%pw(:,3) = CMPLX(0.0,0.0)
350 0 : IF (input%film) THEN
351 0 : den%vac(:,:,:,3) = CMPLX(0.0,0.0)
352 : END IF
353 : END IF
354 0 : ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
355 0 : den%pw(:,3) = CMPLX(0.0,0.0)
356 0 : IF (input%film) THEN
357 0 : den%vac(:,:,:,3) = CMPLX(0.0,0.0)
358 : END IF
359 : END IF
360 0 : CLOSE(iUnit)
361 : END IF
362 :
363 80 : INQUIRE(FILE='n_mmp_mat',EXIST=l_exist)
364 80 : IF(l_exist.AND.atoms%n_denmat.GT.0) THEN
365 2 : OPEN (69,file='n_mmp_mat',status='unknown',form='formatted')
366 2 : READ (69,'(7f20.13)',IOSTAT=ioStatus) den%mmpMat
367 2 : REWIND(69)
368 2 : numLines = 0
369 112 : DO
370 114 : READ (69,*,iostat=iofl)
371 114 : IF (iofl < 0) EXIT
372 112 : numLines = numLines + 1
373 : END DO
374 2 : IF (MOD(numLines,14*SIZE(den%mmpMat,4)).NE.0) THEN
375 0 : WRITE(*,*) 'The n_mmp_mat file could not be read.'
376 0 : WRITE(*,*) 'Was it an old style file with linear mixing parameter specification'
377 0 : WRITE(*,*) 'in the last line? Linear mixing for the density matrix can now be'
378 0 : WRITE(*,*) 'activated and specified in the inp.xml file.'
379 0 : CALL juDFT_error("strange n_mmp_mat-file...", calledby = "readDensity")
380 : END IF
381 2 : CLOSE(69)
382 : END IF
383 :
384 160 : END SUBROUTINE readDensity
385 :
386 482 : SUBROUTINE writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,inOrOutCDN,&
387 : relCdnIndex,distance,fermiEnergy,mmpmatDistance,occDistance,l_qfix,den,inFilename,denIm)
388 :
389 : TYPE(t_noco),INTENT(IN) :: noco
390 : TYPE(t_stars),INTENT(IN) :: stars
391 : TYPE(t_vacuum),INTENT(IN) :: vacuum
392 : TYPE(t_atoms),INTENT(IN) :: atoms
393 : TYPE(t_cell), INTENT(IN) :: cell
394 : TYPE(t_sphhar),INTENT(IN) :: sphhar
395 : TYPE(t_input),INTENT(IN) :: input
396 : TYPE(t_sym),INTENT(IN) :: sym
397 :
398 :
399 : TYPE(t_potden),INTENT(INOUT) :: den
400 :
401 : INTEGER, INTENT (IN) :: inOrOutCDN
402 : INTEGER, INTENT (IN) :: relCdnIndex
403 : INTEGER, INTENT (IN) :: archiveType
404 : REAL, INTENT (IN) :: fermiEnergy, distance
405 : REAL, INTENT (IN) :: mmpmatDistance,occDistance
406 : LOGICAL, INTENT (IN) :: l_qfix
407 :
408 : CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: inFilename
409 :
410 : TYPE(t_potden), OPTIONAL, INTENT(INOUT) :: denIm
411 :
412 1928 : TYPE(t_stars) :: starsTemp
413 482 : TYPE(t_vacuum) :: vacuumTemp
414 482 : TYPE(t_atoms) :: atomsTemp
415 482 : TYPE(t_sphhar) :: sphharTemp
416 : TYPE(t_input) :: inputTemp
417 482 : TYPE(t_sym) :: symTemp
418 : TYPE(t_cell) :: cellTemp
419 :
420 :
421 482 : COMPLEX, ALLOCATABLE :: fpwTemp(:,:), fvacTemp(:,:,:,:)
422 482 : REAL, ALLOCATABLE :: frTemp(:,:,:,:)
423 :
424 : INTEGER :: mode, iterTemp, k, i, iVac, j, iUnit
425 : INTEGER :: d1, d10, asciioffset, iUnitTemp
426 : LOGICAL :: l_exist, l_storeIndices, l_writeNew, l_same
427 : LOGICAL :: l_writeAll
428 : CHARACTER(len=30) :: filename
429 : CHARACTER(len=5) :: cdnfile
430 :
431 : #ifdef CPP_HDF
432 : INTEGER(HID_T) :: fileID
433 : #endif
434 : INTEGER :: currentStarsIndex,currentLatharmsIndex
435 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
436 : INTEGER :: readDensityIndex, writeDensityIndex, lastDensityIndex
437 : INTEGER :: previousDensityIndex, densityType
438 : INTEGER :: starsIndexTemp, latharmsIndexTemp, structureIndexTemp
439 : INTEGER :: stepfunctionIndexTemp
440 : INTEGER :: jspinsTemp
441 : INTEGER :: date, time, dateTemp, timeTemp
442 : INTEGER :: iSpin, iV, iPair, iAtom1, iAtom2
443 : REAL :: fermiEnergyTemp, distanceTemp
444 : LOGICAL :: l_qfixTemp, l_CheckBroyd
445 : CHARACTER(LEN=30) :: archiveName
446 : CHARACTER(LEN=8) :: dateString
447 : CHARACTER(LEN=10) :: timeString
448 : CHARACTER(LEN=10) :: zone
449 :
450 482 : COMPLEX, ALLOCATABLE :: cdomvz(:,:)
451 :
452 482 : CALL getIOMode(mode)
453 482 : CALL DATE_AND_TIME(dateString,timeString,zone)
454 482 : READ(dateString,'(i8)') date
455 482 : READ(timeString,'(i6)') time
456 :
457 482 : l_CheckBroyd = .TRUE.
458 482 : IF(PRESENT(inFilename)) l_CheckBroyd = .FALSE.
459 :
460 482 : IF(mode.EQ.CDN_HDF5_MODE) THEN
461 : #ifdef CPP_HDF
462 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
463 963 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
464 :
465 : CALL checkAndWriteMetadataHDF(fileID, input, atoms, cell, vacuum, stars, sphhar, sym,&
466 : currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
467 482 : currentStepfunctionIndex,l_storeIndices,l_CheckBroyd,.FALSE.)
468 :
469 482 : previousDensityIndex = readDensityIndex
470 482 : writeDensityIndex = readDensityIndex
471 482 : IF(relCdnIndex.NE.0) THEN
472 403 : writeDensityIndex = lastDensityIndex+relCdnIndex
473 403 : lastDensityIndex = writeDensityIndex
474 403 : readDensityIndex = writeDensityIndex
475 403 : l_storeIndices = .TRUE.
476 : END IF
477 :
478 482 : archiveName = ''
479 482 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
480 1 : archiveName = 'cdn'
481 : ELSE
482 481 : WRITE(archiveName,'(a,i0)') '/cdn-', writeDensityIndex
483 : END IF
484 :
485 482 : densityType = 0
486 964 : SELECT CASE (inOrOutCDN)
487 : CASE (CDN_INPUT_DEN_const)
488 482 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_FFN_const) THEN
489 27 : densityType = DENSITY_TYPE_FFN_IN_const
490 455 : ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
491 101 : densityType = DENSITY_TYPE_NOCO_IN_const
492 : ELSE
493 354 : densityType = DENSITY_TYPE_IN_const
494 : END IF
495 : CASE (CDN_OUTPUT_DEN_const)
496 0 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_FFN_const) THEN
497 0 : densityType = DENSITY_TYPE_FFN_OUT_const
498 0 : ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
499 0 : densityType = DENSITY_TYPE_NOCO_OUT_const
500 : ELSE
501 0 : densityType = DENSITY_TYPE_OUT_const
502 : END IF
503 : CASE DEFAULT
504 0 : WRITE(*,*) 'inOrOutCDN = ', inOrOutCDN
505 482 : CALL juDFT_error("Invalid inOrOutCDN selected.",calledby ="writeDensity")
506 : END SELECT
507 :
508 482 : IF(relCdnIndex.EQ.0) THEN
509 79 : l_exist = isDensityEntryPresentHDF(fileID,archiveName,DENSITY_TYPE_UNDEFINED_const)
510 79 : IF(l_exist) THEN
511 : CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_UNDEFINED_const,&
512 : iterTemp, starsIndexTemp, latharmsIndexTemp, structureIndexTemp,&
513 : stepfunctionIndexTemp,previousDensityIndex, jspinsTemp,&
514 78 : dateTemp, timeTemp, distanceTemp, fermiEnergyTemp, l_qfixTemp)
515 : END IF
516 : END IF
517 :
518 482 : IF(vacuum%nvac.EQ.1) THEN
519 45 : IF (sym%invs) THEN
520 1428248 : den%vac(:,:,2,:) = CONJG(den%vac(:,:,1,:))
521 : ELSE
522 3097886 : den%vac(:,:,2,:) = den%vac(:,:,1,:)
523 : END IF
524 : END IF
525 :
526 482 : IF (PRESENT(denIm)) THEN
527 : CALL writeDensityHDF(input, fileID, archiveName, densityType, previousDensityIndex,&
528 : currentStarsIndex, currentLatharmsIndex, currentStructureIndex,&
529 : currentStepfunctionIndex,date,time,distance,fermiEnergy,mmpmatDistance,&
530 0 : occDistance,l_qfix,den%iter+relCdnIndex,den,denIm)
531 : ELSE
532 : CALL writeDensityHDF(input, fileID, archiveName, densityType, previousDensityIndex,&
533 : currentStarsIndex, currentLatharmsIndex, currentStructureIndex,&
534 : currentStepfunctionIndex,date,time,distance,fermiEnergy,mmpmatDistance,&
535 482 : occDistance,l_qfix,den%iter+relCdnIndex,den)
536 : END IF
537 :
538 482 : IF(l_storeIndices) THEN
539 : CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
540 404 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
541 : END IF
542 :
543 482 : CALL closeCDNPOT_HDF(fileID)
544 : #endif
545 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
546 : ! Write density to cdn.str file
547 : STOP 'CDN_STREAM_MODE not yet implemented!'
548 : ELSE
549 0 : filename = 'cdn1'
550 0 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
551 0 : filename = 'rhomat_inp'
552 0 : IF(inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) filename = 'rhomat_out'
553 : END IF
554 0 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
555 0 : filename = 'cdn'
556 : END IF
557 :
558 0 : IF(PRESENT(inFilename)) filename = inFilename
559 :
560 0 : IF ((relCdnIndex.EQ.1).AND.(archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).AND.(den%iter.EQ.0)) THEN
561 0 : INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
562 0 : IF(l_exist) THEN
563 0 : CALL juDFT_error("Trying to generate starting density while a density exists.",calledby ="writeDensity")
564 : END IF
565 : END IF
566 :
567 0 : iUnit = 93
568 0 : OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='unknown')
569 :
570 0 : IF ((relCdnIndex.EQ.1).AND.(archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).AND.(den%iter.GE.1)) THEN
571 0 : inputTemp%jspins = input%jspins
572 0 : vacuumTemp%nmzxyd = vacuum%nmzxyd
573 0 : atomsTemp%jmtd = atoms%jmtd
574 0 : sphharTemp%nlhd = sphhar%nlhd
575 0 : vacuumTemp%nmzd = vacuum%nmzd
576 0 : atomsTemp%ntype = atoms%ntype
577 0 : atomsTemp%firstAtom=atoms%firstAtom
578 0 : ALLOCATE (sphharTemp%nlh(SIZE(sphhar%nlh)))
579 0 : sphharTemp%nlh(:) = sphhar%nlh(:)
580 0 : ALLOCATE (symTemp%ntypsy(SIZE(sym%ntypsy)))
581 0 : symTemp%ntypsy(:) = sym%ntypsy(:)
582 0 : ALLOCATE (atomsTemp%jri(SIZE(atoms%jri)))
583 0 : atomsTemp%jri(:) = atoms%jri(:)
584 0 : ALLOCATE (atomsTemp%neq(SIZE(atoms%neq)))
585 0 : atomsTemp%neq(:) = atoms%neq(:)
586 0 : ALLOCATE (atomsTemp%zatom(SIZE(atoms%zatom)))
587 0 : atomsTemp%zatom(:) = atoms%zatom(:)
588 0 : ALLOCATE (atomsTemp%rmt(SIZE(atoms%rmt)))
589 0 : atomsTemp%rmt(:) = atoms%rmt(:)
590 0 : ALLOCATE (atomsTemp%dx(SIZE(atoms%dx)))
591 0 : atomsTemp%dx(:) = atoms%dx(:)
592 0 : starsTemp%ng3 = stars%ng3
593 0 : symTemp%invs = sym%invs
594 0 : inputTemp%film = input%film
595 0 : vacuumTemp%nvac = vacuum%nvac
596 0 : vacuumTemp%nmzxy = vacuum%nmzxy
597 0 : vacuumTemp%nmz = vacuum%nmz
598 0 : vacuumTemp%dvac = vacuum%dvac
599 0 : vacuumTemp%delz = vacuum%delz
600 0 : starsTemp%ng2 = stars%ng2
601 0 : symTemp%invs2 = sym%invs2
602 0 : ALLOCATE (fpwTemp(stars%ng3,input%jspins))
603 0 : ALLOCATE (fvacTemp(vacuum%nmzd,stars%ng2,2,input%jspins))
604 0 : ALLOCATE (frTemp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins))
605 :
606 : !---> generate name of file to hold the results of this iteration
607 0 : d1 = MOD(den%iter,10)
608 0 : d10 = MOD(INT((den%iter+0.5)/10),10)
609 : asciioffset = IACHAR('1')-1
610 : IF ( d10.GE.10 ) asciioffset = IACHAR('7')
611 0 : cdnfile = 'cdn'//ACHAR(d10+asciioffset)//ACHAR(d1+IACHAR('1')-1)
612 :
613 0 : iUnitTemp = 72
614 0 : OPEN (iUnitTemp,file=cdnfile,form='unformatted',status='unknown')
615 0 : REWIND iUnitTemp
616 :
617 : CALL loddop(starsTemp,vacuumTemp,atomsTemp,sphharTemp,inputTemp,symTemp,&
618 0 : iUnit,iterTemp,frTemp,fpwTemp,fvacTemp)
619 : CALL wrtdop(starsTemp,vacuumTemp,atomsTemp,sphharTemp,inputTemp,symTemp,&
620 0 : iUnitTemp,iterTemp,frTemp,fpwTemp,fvacTemp)
621 :
622 0 : CLOSE(iUnitTemp)
623 0 : REWIND iUnit
624 :
625 0 : DEALLOCATE (frTemp, fpwTemp, fvacTemp)
626 0 : DEALLOCATE (atomsTemp%neq, atomsTemp%jri, atomsTemp%zatom, symTemp%ntypsy, sphharTemp%nlh)
627 0 : DEALLOCATE (atomsTemp%rmt, atomsTemp%dx)
628 : END IF
629 :
630 0 : IF ((inOrOutCDN.EQ.CDN_OUTPUT_DEN_const).AND.(archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const)) THEN
631 :
632 : ! Generate data in temp arrays and variables to be able to perform loddop call.
633 : ! loddop is called to move the file position to the output density position.
634 0 : inputTemp%jspins = input%jspins
635 0 : vacuumTemp%nmzxyd = vacuum%nmzxyd
636 0 : atomsTemp%jmtd = atoms%jmtd
637 0 : sphharTemp%nlhd = sphhar%nlhd
638 0 : vacuumTemp%nmzd = vacuum%nmzd
639 0 : atomsTemp%ntype = atoms%ntype
640 0 : ALLOCATE (sphharTemp%nlh(SIZE(sphhar%nlh)))
641 0 : sphharTemp%nlh(:) = sphhar%nlh(:)
642 0 : ALLOCATE (symTemp%ntypsy(SIZE(sym%ntypsy)))
643 0 : symTemp%ntypsy(:) = sym%ntypsy(:)
644 0 : ALLOCATE (atomsTemp%jri(SIZE(atoms%jri)))
645 0 : atomsTemp%jri(:) = atoms%jri(:)
646 0 : ALLOCATE (atomsTemp%neq(SIZE(atoms%neq)))
647 0 : atomsTemp%neq(:) = atoms%neq(:)
648 0 : starsTemp%ng3 = stars%ng3
649 0 : symTemp%invs = sym%invs
650 0 : inputTemp%film = input%film
651 0 : vacuumTemp%nvac = vacuum%nvac
652 0 : vacuumTemp%nmzxy = vacuum%nmzxy
653 0 : vacuumTemp%nmz = vacuum%nmz
654 0 : vacuumTemp%dvac = vacuum%dvac
655 0 : vacuumTemp%delz = vacuum%delz
656 0 : starsTemp%ng2 = stars%ng2
657 0 : symTemp%invs2 = sym%invs2
658 0 : ALLOCATE (fpwTemp(stars%ng3,input%jspins))
659 0 : ALLOCATE (fvacTemp(vacuum%nmzd,stars%ng2,2,input%jspins))
660 0 : ALLOCATE (frTemp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins))
661 :
662 : CALL loddop(starsTemp,vacuumTemp,atomsTemp,sphharTemp,inputTemp,symTemp,&
663 0 : iUnit,iterTemp,frTemp,fpwTemp,fvacTemp)
664 :
665 0 : DEALLOCATE (frTemp, fpwTemp, fvacTemp)
666 0 : DEALLOCATE (atomsTemp%neq, atomsTemp%jri, symTemp%ntypsy, sphharTemp%nlh)
667 : END IF
668 :
669 : ! Write the density
670 : CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
671 0 : iUnit,den%iter+relCdnIndex,den%mt,den%pw,den%vac)
672 :
673 : ! Write additional data if l_noco
674 0 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
675 0 : WRITE (iUnit) (den%pw(k,3),k=1,stars%ng3)
676 0 : IF (input%film) THEN
677 0 : ALLOCATE(cdomvz(vacuum%nmz,vacuum%nvac))
678 0 : DO iVac = 1, vacuum%nvac
679 0 : DO i = 1, vacuum%nmz
680 0 : cdomvz(i,iVac) = den%vac(i,1,iVac,3)
681 : END DO
682 : END DO
683 0 : WRITE (iUnit) ((cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
684 0 : WRITE (iUnit) (((den%vac(i,j,iVac,3),i=1,vacuum%nmzxy),j=2,stars%ng2), iVac=1,vacuum%nvac)
685 0 : DEALLOCATE(cdomvz)
686 : END IF
687 : END IF
688 :
689 0 : CLOSE(iUnit)
690 : END IF
691 :
692 : !write density matrix to n_mmp_mat_out file
693 482 : IF((inOrOutCDN.EQ.CDN_INPUT_DEN_const).AND.(relCdnIndex.EQ.1).AND.&
694 : ((archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).OR.(archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const).OR.&
695 : (archiveType.EQ.CDN_ARCHIVE_TYPE_FFN_const))) THEN
696 403 : IF(atoms%n_denmat.GT.0) THEN
697 49 : filename = 'n_mmp_mat'
698 49 : IF (mode.EQ.CDN_HDF5_MODE) THEN
699 49 : filename = 'n_mmp_mat_out'
700 : END IF
701 1188 : IF(ANY(den%mmpMat(:,:,:,:).NE.0.0)) THEN
702 46 : IF ((mode.EQ.CDN_HDF5_MODE).AND..NOT.(input%ldauLinMix.AND.(input%ldauMixParam.EQ.0.0))) THEN
703 46 : INQUIRE(file='n_mmp_mat',exist=l_exist)
704 46 : IF(l_exist) THEN
705 1 : CALL system('mv n_mmp_mat n_mmp_mat_old')
706 1 : PRINT *,"n_mmp_mat moved to n_mmp_mat_old"
707 : END IF
708 : END IF
709 46 : OPEN (69,file=TRIM(ADJUSTL(filename)),status='replace',form='formatted')
710 46 : WRITE (69,'(7f20.13)') den%mmpMat(:,:,:,:)
711 46 : CLOSE (69)
712 : END IF
713 : END IF
714 : END IF
715 :
716 : !write density matrix for LDA+V to nIJ_mmp_mat_out file
717 482 : IF((inOrOutCDN.EQ.CDN_INPUT_DEN_const).AND.(relCdnIndex.EQ.1).AND.&
718 : ((archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).OR.(archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const).OR.&
719 : (archiveType.EQ.CDN_ARCHIVE_TYPE_FFN_const))) THEN
720 403 : IF(atoms%n_v.GT.0) THEN
721 0 : filename = 'nIJ_mmp_mat'
722 0 : IF (mode.EQ.CDN_HDF5_MODE) THEN
723 0 : filename = 'nIJ_mmp_mat_out'
724 : END IF
725 :
726 0 : IF ((mode.EQ.CDN_HDF5_MODE).AND..NOT.(input%ldauLinMix.AND.(input%ldauMixParam.EQ.0.0))) THEN
727 0 : INQUIRE(file='nIJ_mmp_mat',exist=l_exist)
728 0 : IF(l_exist) THEN
729 0 : CALL system('mv nIJ_mmp_mat nIJ_mmp_mat_old')
730 0 : PRINT *,"nIJ_mmp_mat moved to nIJ_mmp_mat_old"
731 : END IF
732 : END IF
733 0 : OPEN (79,file=TRIM(ADJUSTL(filename)),status='replace',form='formatted')
734 0 : DO iSpin = 1, input%jspins
735 0 : iPair = 0
736 0 : DO iV = 1, atoms%n_v
737 0 : iAtom1 = atoms%lda_v(iV)%atomIndex
738 0 : DO iAtom2 = 1, atoms%lda_v(iV)%numOtherAtoms
739 0 : iPair = iPair + 1
740 0 : WRITE(79,'(a,4i7)') 'spin,pair_index,atom1,atom2',iSpin, iPair, iAtom1, iAtom2
741 0 : WRITE (79,'(7f20.13)') den%nIJ_llp_mmp(:,:,iPair,iSpin)
742 : END DO
743 : END DO
744 : END DO
745 0 : CLOSE (79)
746 : END IF
747 : END IF
748 :
749 1928 : END SUBROUTINE writeDensity
750 :
751 52 : SUBROUTINE readPrevEFermi(eFermiPrev,l_error)
752 :
753 : REAL, INTENT(OUT) :: eFermiPrev
754 : LOGICAL, INTENT(OUT) :: l_error
755 :
756 : INTEGER :: mode
757 : #ifdef CPP_HDF
758 : INTEGER(HID_T) :: fileID
759 : #endif
760 : INTEGER :: currentStarsIndex,currentLatharmsIndex
761 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
762 : INTEGER :: readDensityIndex, lastDensityIndex
763 :
764 : INTEGER :: starsIndex, latharmsIndex, structureIndex
765 : INTEGER :: stepfunctionIndex
766 : INTEGER :: date, time, iter, jspins, previousDensityIndex
767 : REAL :: fermiEnergy, distance
768 : LOGICAL :: l_qfix, l_exist
769 : CHARACTER(LEN=30) :: archiveName
770 :
771 26 : CALL getIOMode(mode)
772 :
773 26 : eFermiPrev = 0.0
774 26 : l_error = .FALSE.
775 26 : IF(mode.EQ.CDN_HDF5_MODE) THEN
776 : #ifdef CPP_HDF
777 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
778 26 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
779 26 : WRITE(archiveName,'(a,i0)') '/cdn-', readDensityIndex
780 : CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_FFN_IN_const,&
781 : iter, starsIndex, latharmsIndex, structureIndex, stepfunctionIndex,&
782 26 : previousDensityIndex, jspins, date, time, distance, fermiEnergy, l_qfix)
783 26 : eFermiPrev = fermiEnergy
784 26 : CALL closeCDNPOT_HDF(fileID)
785 : #endif
786 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
787 : STOP 'cdn.str not yet implemented!'
788 : ELSE
789 0 : l_error = .TRUE.
790 : END IF
791 :
792 26 : END SUBROUTINE readPrevEFermi
793 :
794 0 : SUBROUTINE readPrevmmpDistances(mmpmatDistancePrev,occDistancePrev,l_error)
795 :
796 : REAL, INTENT(OUT) :: mmpmatDistancePrev,occDistancePrev
797 : LOGICAL, INTENT(OUT) :: l_error
798 :
799 : INTEGER :: mode
800 : #ifdef CPP_HDF
801 : INTEGER(HID_T) :: fileID
802 : #endif
803 : INTEGER :: currentStarsIndex,currentLatharmsIndex
804 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
805 : INTEGER :: readDensityIndex, lastDensityIndex
806 :
807 : INTEGER :: starsIndex, latharmsIndex, structureIndex
808 : INTEGER :: stepfunctionIndex
809 : INTEGER :: date, time, iter, jspins, previousDensityIndex
810 : REAL :: fermiEnergy, distance
811 : REAL :: mmpmatDistance,occDistance
812 : LOGICAL :: l_qfix, l_exist
813 : CHARACTER(LEN=30) :: archiveName
814 :
815 0 : CALL getIOMode(mode)
816 :
817 0 : mmpmatDistancePrev = 0.0
818 0 : occDistancePrev = 0.0
819 0 : l_error = .FALSE.
820 0 : IF(mode.EQ.CDN_HDF5_MODE) THEN
821 : #ifdef CPP_HDF
822 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
823 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
824 0 : WRITE(archiveName,'(a,i0)') '/cdn-', readDensityIndex
825 : CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_FFN_IN_const,&
826 : iter, starsIndex, latharmsIndex, structureIndex, stepfunctionIndex,&
827 : previousDensityIndex, jspins, date, time, distance, fermiEnergy, l_qfix,&
828 0 : mmpmatDistance,occDistance)
829 0 : IF(mmpmatDistance.GE.-1e-10.AND.occDistance.GE.-1e-10) THEN
830 0 : mmpmatDistancePrev = mmpmatDistance
831 0 : occDistancePrev = occDistance
832 : ELSE
833 0 : l_error = .TRUE.
834 : ENDIF
835 0 : CALL closeCDNPOT_HDF(fileID)
836 : #endif
837 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
838 : STOP 'cdn.str not yet implemented!'
839 : ELSE
840 0 : l_error = .TRUE.
841 : END IF
842 :
843 0 : END SUBROUTINE readPrevmmpDistances
844 :
845 74 : SUBROUTINE readCoreDensity(input,atoms,rhcs,tecs,qints)
846 :
847 : TYPE(t_atoms),INTENT(IN) :: atoms
848 : TYPE(t_input),INTENT(IN) :: input
849 :
850 :
851 : REAL, INTENT(OUT) :: rhcs(:,:,:)!(atoms%jmtd,atoms%ntype,input%jspins)
852 : REAL, INTENT(OUT) :: tecs(:,:)!(atoms%ntype,input%jspins)
853 : REAL, INTENT(OUT) :: qints(:,:)!(atoms%ntype,input%jspins)
854 :
855 : INTEGER :: mode, iUnit, iSpin, iAtom, i
856 : LOGICAL :: l_exist
857 : CHARACTER(LEN=30) :: filename
858 :
859 : #ifdef CPP_HDF
860 : INTEGER(HID_T) :: fileID
861 : #endif
862 : INTEGER :: currentStarsIndex,currentLatharmsIndex
863 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
864 : INTEGER :: readDensityIndex, lastDensityIndex
865 :
866 74 : CALL getIOMode(mode)
867 :
868 : IF(mode.EQ.CDN_HDF5_MODE) THEN
869 : #ifdef CPP_HDF
870 74 : l_exist = isCoreDensityPresentHDF()
871 74 : IF (l_exist) THEN
872 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
873 74 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
874 74 : CALL readCoreDensityHDF(fileID,input,atoms,rhcs,tecs,qints)
875 74 : CALL closeCDNPOT_HDF(fileID)
876 74 : RETURN
877 : ELSE
878 0 : WRITE(*,*) 'No core density is available in HDF5 format.'
879 0 : WRITE(*,*) 'Falling back to stream access file cdn.str.'
880 : mode = CDN_STREAM_MODE
881 : END IF
882 : #endif
883 : END IF
884 :
885 : IF(mode.EQ.CDN_STREAM_MODE) THEN
886 0 : INQUIRE(FILE='cdn.str',EXIST=l_exist)
887 0 : IF (l_exist) THEN
888 : !load density from cdn.str and exit subroutine
889 :
890 : RETURN
891 : ELSE
892 0 : WRITE(*,*) 'cdn.str file not found.'
893 0 : WRITE(*,*) 'Falling back to direct access file cdnc.'
894 : mode = CDN_DIRECT_MODE
895 : END IF
896 : END IF
897 :
898 : IF (mode.EQ.CDN_DIRECT_MODE) THEN
899 0 : filename = 'cdnc'
900 0 : INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
901 0 : IF(.NOT.l_exist) THEN
902 0 : CALL juDFT_error("core charge density file "//TRIM(ADJUSTL(filename))//" missing",calledby ="readCoreDensity")
903 : END IF
904 0 : iUnit = 17
905 0 : OPEN (iUnit,file=TRIM(ADJUSTL(filename)),form='unformatted',status='unknown')
906 0 : DO iSpin = 1,input%jspins
907 0 : DO iAtom = 1,atoms%ntype
908 0 : READ (iUnit) (rhcs(i,iAtom,iSpin),i=1,atoms%jri(iAtom))
909 0 : READ (iUnit) tecs(iAtom,iSpin)
910 : END DO
911 0 : READ (iUnit) (qints(iAtom,iSpin),iAtom=1,atoms%ntype)
912 : END DO
913 0 : CLOSE (iUnit)
914 : END IF
915 :
916 : END SUBROUTINE readCoreDensity
917 :
918 331 : SUBROUTINE writeCoreDensity(input,atoms,rhcs,tecs,qints,filename)
919 :
920 : TYPE(t_atoms),INTENT(IN) :: atoms
921 : TYPE(t_input),INTENT(IN) :: input
922 :
923 :
924 : REAL, INTENT(IN) :: rhcs(:,:,:)!(atoms%jmtd,atoms%ntype,input%jspins)
925 : REAL, INTENT(IN) :: tecs(:,:)!(atoms%ntype,input%jspins)
926 : REAL, INTENT(IN) :: qints(:,:)!(atoms%ntype,input%jspins)
927 : CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: filename
928 :
929 : INTEGER :: mode, iUnit, iSpin, iAtom, i
930 :
931 : #ifdef CPP_HDF
932 : INTEGER(HID_T) :: fileID
933 : #endif
934 : INTEGER :: currentStarsIndex,currentLatharmsIndex
935 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
936 : INTEGER :: readDensityIndex, lastDensityIndex
937 :
938 331 : CALL getIOMode(mode)
939 :
940 : IF(mode.EQ.CDN_HDF5_MODE) THEN
941 : #ifdef CPP_HDF
942 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
943 662 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex,filename)
944 331 : CALL writeCoreDensityHDF(fileID,input,atoms,rhcs,tecs,qints)
945 331 : CALL closeCDNPOT_HDF(fileID)
946 : #endif
947 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
948 : ! Write core density to cdn.str file
949 : STOP 'CDN_STREAM_MODE not yet implemented!'
950 : ELSE
951 0 : iUnit = 17
952 0 : OPEN (iUnit,file='cdnc',form='unformatted',status='unknown')
953 0 : DO iSpin = 1,input%jspins
954 0 : DO iAtom = 1,atoms%ntype
955 0 : WRITE (iUnit) (rhcs(i,iAtom,iSpin),i=1,atoms%jri(iAtom))
956 0 : WRITE (iUnit) tecs(iAtom,iSpin)
957 : END DO
958 0 : WRITE (iUnit) (qints(iAtom,iSpin),iAtom=1,atoms%ntype)
959 : END DO
960 0 : CLOSE (iUnit)
961 : END IF
962 :
963 331 : END SUBROUTINE writeCoreDensity
964 :
965 160 : SUBROUTINE storeStructureIfNew(input,stars, atoms, cell, vacuum, sym,fmpi,sphhar,noco)
966 :
967 : TYPE(t_input),INTENT(IN) :: input
968 : TYPE(t_atoms), INTENT(IN) :: atoms
969 : TYPE(t_cell), INTENT(IN) :: cell
970 : TYPE(t_vacuum), INTENT(IN) :: vacuum
971 :
972 : TYPE(t_sym),INTENT(IN) :: sym
973 : TYPE(t_mpi),INTENT(IN) :: fmpi
974 : TYPE(t_sphhar),INTENT(IN) :: sphhar
975 : TYPE(t_noco),INTENT(IN) :: noco
976 : TYPE(t_stars),INTENT(IN) :: stars
977 :
978 :
979 : TYPE(t_input) :: inputTemp
980 160 : TYPE(t_atoms) :: atomsTemp
981 : TYPE(t_cell) :: cellTemp
982 160 : TYPE(t_vacuum) :: vacuumTemp
983 :
984 160 : TYPE(t_sym) :: symTemp
985 :
986 :
987 : INTEGER :: mode
988 : INTEGER :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
989 : INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
990 : LOGICAL :: l_writeStructure, l_same
991 : #ifdef CPP_HDF
992 : INTEGER(HID_T) :: fileID
993 : #endif
994 : #ifdef CPP_MPI
995 : INTEGER :: ierr
996 : #endif
997 :
998 160 : CALL getIOMode(mode)
999 :
1000 160 : IF (fmpi%irank==0) THEN
1001 :
1002 80 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1003 : #ifdef CPP_HDF
1004 80 : l_writeStructure = .FALSE.
1005 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1006 80 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1007 :
1008 80 : IF (currentStructureIndex.EQ.0) THEN
1009 68 : currentStructureIndex = currentStructureIndex + 1
1010 : l_writeStructure = .TRUE.
1011 : ELSE
1012 12 : CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, symTemp,currentStructureIndex)
1013 12 : CALL compareStructure(input, atoms, vacuum, cell, sym, inputTemp, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same)
1014 12 : IF(.NOT.l_same) THEN
1015 0 : currentStructureIndex = currentStructureIndex + 1
1016 : l_writeStructure = .TRUE.
1017 : END IF
1018 : END IF
1019 :
1020 :
1021 : IF (l_writeStructure) THEN
1022 68 : CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, sym, currentStructureIndex,.TRUE.)
1023 : CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1024 68 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1025 : END IF
1026 :
1027 80 : CALL closeCDNPOT_HDF(fileID)
1028 : #endif
1029 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
1030 : ! Write stars to stars file
1031 : STOP 'CDN_STREAM_MODE not yet implemented!'
1032 : ELSE
1033 : ! In direct access mode no structure information is written to any file.
1034 : END IF
1035 : ENDIF
1036 : #ifdef CPP_MPI
1037 160 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr)
1038 : #endif
1039 8800 : END SUBROUTINE storeStructureIfNew
1040 :
1041 108 : SUBROUTINE transform_by_moving_atoms(fmpi,stars,atoms,vacuum, cell, sym, sphhar,input ,noco,nococonv)
1042 : USE m_types
1043 : USE m_constants
1044 : USE m_qfix
1045 : USE m_fix_by_gaussian
1046 : IMPLICIT NONE
1047 : TYPE(t_mpi),INTENT(IN) :: fmpi
1048 : TYPE(t_atoms),INTENT(IN) :: atoms
1049 : TYPE(t_sym),INTENT(IN) :: sym
1050 : TYPE(t_vacuum),INTENT(IN) :: vacuum
1051 : TYPE(t_sphhar),INTENT(IN) :: sphhar
1052 : TYPE(t_input),INTENT(IN) :: input
1053 :
1054 : TYPE(t_cell),INTENT(IN) :: cell
1055 : TYPE(t_noco),INTENT(IN) :: noco
1056 : TYPE(t_nococonv),INTENT(IN) :: nococonv
1057 :
1058 : TYPE(t_stars),INTENT(IN) :: stars
1059 :
1060 : !Locals
1061 : INTEGER :: archiveType
1062 : LOGICAL :: l_qfix
1063 : REAL :: fermiEnergy,fix, tempDistance
1064 108 : REAL :: shifts(3,atoms%nat)
1065 :
1066 108 : TYPE(t_potden):: den
1067 :
1068 : TYPE(t_input) :: inputTemp
1069 108 : TYPE(t_atoms) :: atomsTemp
1070 : TYPE(t_cell) :: cellTemp
1071 108 : TYPE(t_vacuum) :: vacuumTemp
1072 :
1073 108 : TYPE(t_sym) :: symTemp
1074 :
1075 :
1076 :
1077 : INTEGER :: mode,currentStarsIndex,currentLatharmsIndex,currentStructureIndex
1078 : INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex,structureindex
1079 : LOGICAL :: l_same,l_structure_by_shift
1080 :
1081 : #ifdef CPP_HDF
1082 : INTEGER(HID_T) :: fileID
1083 : CHARACTER(len=50) :: archivename
1084 : #endif
1085 : #ifdef CPP_MPI
1086 : INTEGER :: ierr
1087 : #endif
1088 :
1089 108 : l_same=.TRUE.;l_structure_by_shift=.TRUE.
1090 :
1091 108 : CALL getIOMode(mode)
1092 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1093 : #ifdef CPP_HDF
1094 108 : IF (fmpi%irank==0) THEN
1095 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1096 54 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1097 54 : IF (currentStructureIndex>0.AND.lastdensityindex>0) THEN
1098 : !Determine structure of last density
1099 9 : WRITE(archivename,"(a,i0)") "cdn-",lastdensityindex
1100 9 : CALL peekDensityEntryHDF(fileID, archivename, DENSITY_TYPE_IN_const, structureIndex=structureIndex)
1101 : !Read that structure
1102 9 : CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, symTemp,StructureIndex)
1103 9 : CALL compareStructure(input, atoms, vacuum, cell, sym, inputTemp, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same,l_structure_by_shift)
1104 : ENDIF
1105 54 : CALL closeCDNPOT_HDF(fileID)
1106 : ENDIF
1107 : #ifdef CPP_MPI
1108 108 : CALL mpi_bcast(l_same,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
1109 108 : CALL mpi_bcast(l_structure_by_shift,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
1110 : #endif
1111 108 : IF (l_same.OR..NOT.l_structure_by_shift) RETURN ! nothing to do
1112 :
1113 0 : IF (fmpi%irank==0) THEN
1114 0 : WRITE(oUnit,*) "Atomic movement detected, trying to adjust charge density"
1115 :
1116 : !Calculate shifts
1117 0 : shifts=atomsTemp%taual-atoms%taual
1118 :
1119 : !Determine type of charge
1120 0 : IF(any(noco%l_unrestrictMT)) THEN
1121 0 : archiveType=CDN_ARCHIVE_TYPE_FFN_const
1122 0 : ELSE IF (noco%l_noco) THEN
1123 0 : archiveType=CDN_ARCHIVE_TYPE_NOCO_const
1124 : ELSE
1125 0 : archiveType=CDN_ARCHIVE_TYPE_CDN1_const
1126 : END IF
1127 : !read the current density
1128 0 : CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
1129 : CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,CDN_INPUT_DEN_const,&
1130 0 : 0,fermiEnergy,tempDistance,l_qfix,den)
1131 : ENDIF
1132 : !Now fix the density
1133 0 : SELECT CASE(input%qfix)
1134 : CASE (0,1) !just qfix the density
1135 0 : IF (fmpi%irank==0) WRITE(oUnit,*) "Using qfix to adjust density"
1136 0 : IF (fmpi%irank==0) CALL qfix(fmpi,stars,nococonv,atoms,sym,vacuum,sphhar,input,cell ,&
1137 0 : den,noco%l_noco,.FALSE.,.FALSE.,force_fix=.TRUE.,fix=fix)
1138 : CASE(2,3)
1139 0 : IF (fmpi%irank==0) CALL qfix(fmpi,stars,nococonv,atoms,sym,vacuum,sphhar,input,cell ,&
1140 0 : den,noco%l_noco,.FALSE.,.FALSE.,force_fix=.TRUE.,fix=fix,fix_pw_only=.TRUE.)
1141 : CASE(4,5)
1142 0 : IF (fmpi%irank==0) CALL fix_by_gaussian(shifts,atoms,nococonv,stars,fmpi,sym,vacuum,sphhar,input ,cell,noco,den)
1143 : CASE default
1144 0 : CALL judft_error("Wrong choice of qfix in input")
1145 : END SELECT
1146 : !Now write the density to file
1147 0 : IF (fmpi%irank==0) CALL writedensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,CDN_INPUT_DEN_const,&
1148 0 : 0,-1.0,fermiEnergy,-1.0,-1.0,l_qfix,den)
1149 :
1150 : #endif
1151 : END IF
1152 5616 : END SUBROUTINE transform_by_moving_atoms
1153 :
1154 0 : SUBROUTINE writeStars(stars ,l_xcExtended,l_ExtData)
1155 :
1156 : TYPE(t_stars),INTENT(IN) :: stars
1157 :
1158 : LOGICAL, INTENT(IN) :: l_xcExtended, l_ExtData
1159 :
1160 : INTEGER :: mode, ngz, izmin, izmax
1161 :
1162 : INTEGER :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
1163 : INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
1164 : #ifdef CPP_HDF
1165 : INTEGER(HID_T) :: fileID
1166 : #endif
1167 :
1168 0 : INTEGER :: igz(stars%ng3)
1169 :
1170 0 : ngz = 0
1171 0 : izmin = 0
1172 0 : izmax = 0
1173 0 : igz = 0
1174 :
1175 0 : CALL getIOMode(mode)
1176 :
1177 0 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1178 : #ifdef CPP_HDF
1179 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1180 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1181 :
1182 0 : currentStarsIndex = currentStarsIndex + 1
1183 0 : CALL writeStarsHDF(fileID, currentStarsIndex, currentStructureIndex, stars, .TRUE., .FALSE.)
1184 :
1185 : CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1186 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1187 :
1188 0 : CALL closeCDNPOT_HDF(fileID)
1189 : #endif
1190 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
1191 : ! Write stars to stars file
1192 : STOP 'CDN_STREAM_MODE not yet implemented!'
1193 : ELSE
1194 : ! OPEN (51,file='stars',form='unformatted',status='unknown')
1195 : ! WRITE (51) stars%gmax,stars%ng3,stars%ng2,ngz,izmin,izmax,stars%mx1,stars%mx2,stars%mx3
1196 : ! IF(l_ExtData) THEN
1197 : ! IF (.NOT.l_xcExtended) THEN
1198 : ! WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
1199 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1200 : ! stars%igfft2,stars%pgfft2
1201 : ! ELSE
1202 : ! WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
1203 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1204 : ! stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
1205 : ! END IF
1206 : ! ELSE
1207 : ! IF (.NOT.l_xcExtended) THEN
1208 : ! WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
1209 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1210 : ! stars%igfft2,stars%pgfft2
1211 : ! ELSE
1212 : ! WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
1213 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1214 : ! stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
1215 : ! END IF
1216 : ! END IF
1217 : ! CLOSE (51)
1218 : END IF
1219 0 : END SUBROUTINE writeStars
1220 :
1221 0 : SUBROUTINE readStars(stars ,l_xcExtended,l_ExtData,l_error)
1222 :
1223 : TYPE(t_stars),INTENT(INOUT) :: stars
1224 :
1225 : LOGICAL, INTENT(IN) :: l_xcExtended,l_ExtData
1226 : LOGICAL, INTENT(OUT) :: l_error
1227 :
1228 :
1229 0 : TYPE(t_stars) :: starsTemp
1230 :
1231 : INTEGER :: mode, ioStatus, ngz,izmin,izmax
1232 : LOGICAL :: l_exist, l_same
1233 :
1234 : INTEGER :: structureIndexTemp
1235 : INTEGER :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
1236 : INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
1237 : #ifdef CPP_HDF
1238 : INTEGER(HID_T) :: fileID
1239 : #endif
1240 :
1241 0 : INTEGER :: igz(stars%ng3)
1242 :
1243 0 : l_error = .FALSE.
1244 0 : ngz = 0
1245 0 : izmin = 0
1246 0 : izmax = 0
1247 0 : igz = 0
1248 :
1249 0 : CALL getIOMode(mode)
1250 :
1251 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1252 0 : INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
1253 0 : IF (l_exist) THEN
1254 : #ifdef CPP_HDF
1255 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1256 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1257 0 : IF (currentStarsIndex.LT.1) THEN
1258 : mode = CDN_DIRECT_MODE ! (no stars entry found in cdn.hdf file)
1259 : ELSE
1260 0 : CALL peekStarsHDF(fileID, currentStarsIndex, structureIndexTemp)
1261 0 : l_same = structureIndexTemp.EQ.currentStructureIndex
1262 0 : IF(l_same) THEN
1263 0 : CALL readStarsHDF(fileID, currentStarsIndex, starsTemp)
1264 0 : CALL compareStars(stars, starsTemp, l_same)
1265 : END IF
1266 0 : IF(l_same) THEN
1267 0 : CALL readStarsHDF(fileID, currentStarsIndex, stars)
1268 : ELSE
1269 : mode = CDN_DIRECT_MODE ! (no adequate stars entry found in cdn.hdf file)
1270 : END IF
1271 : END IF
1272 0 : CALL closeCDNPOT_HDF(fileID)
1273 : #endif
1274 : END IF
1275 0 : IF(.NOT.l_exist) THEN
1276 : mode = CDN_STREAM_MODE
1277 : END IF
1278 : END IF
1279 :
1280 : IF(mode.EQ.CDN_STREAM_MODE) THEN
1281 0 : INQUIRE(FILE='cdn.str',EXIST=l_exist)
1282 0 : IF (l_exist) THEN
1283 0 : STOP 'cdn.str code path not yet implemented!'
1284 : END IF
1285 : IF (.NOT.l_exist) THEN
1286 : mode = CDN_DIRECT_MODE
1287 : END IF
1288 : END IF
1289 :
1290 0 : IF (mode.EQ.CDN_DIRECT_MODE) THEN
1291 : ! INQUIRE(FILE='stars',EXIST=l_exist)
1292 : ! IF(.NOT.l_exist) THEN
1293 0 : l_error = .TRUE.
1294 0 : RETURN
1295 : ! END IF
1296 : ! OPEN (51,file='stars',form='unformatted',status='unknown')
1297 :
1298 : ! READ (51,IOSTAT=ioStatus) stars%gmax,stars%ng3,stars%ng2,ngz,izmin,izmax,stars%mx1,stars%mx2,stars%mx3
1299 : ! IF (ioStatus.NE.0) THEN
1300 : ! l_error = .TRUE.
1301 : ! RETURN
1302 : ! END IF
1303 :
1304 : ! IF (l_ExtData) THEN
1305 : ! IF (.NOT.l_xcExtended) THEN
1306 : ! READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
1307 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1308 : ! stars%igfft2,stars%pgfft2
1309 : ! stars%ft2_gfx = 0.0
1310 : ! stars%ft2_gfy = 0.0
1311 : ! ELSE
1312 : ! READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
1313 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1314 : ! stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
1315 : ! END IF
1316 : ! ELSE
1317 : ! IF (.NOT.l_xcExtended) THEN
1318 : ! READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
1319 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1320 : ! stars%igfft2,stars%pgfft2
1321 : ! stars%ft2_gfx = 0.0
1322 : ! stars%ft2_gfy = 0.0
1323 : ! ELSE
1324 : ! READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
1325 : ! stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
1326 : ! stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
1327 : ! END IF
1328 : ! END IF
1329 :
1330 : ! IF (ioStatus.NE.0) THEN
1331 : ! l_error = .TRUE.
1332 : ! RETURN
1333 : ! END IF
1334 :
1335 : ! CLOSE (51)
1336 : END IF
1337 :
1338 0 : END SUBROUTINE readStars
1339 :
1340 240 : SUBROUTINE writeStepfunction(stars)
1341 :
1342 : TYPE(t_stars),INTENT(IN) :: stars
1343 :
1344 : INTEGER :: mode, ifftd, i
1345 :
1346 : INTEGER :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
1347 : INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
1348 : #ifdef CPP_HDF
1349 : INTEGER(HID_T) :: fileID
1350 : #endif
1351 :
1352 80 : ifftd=SIZE(stars%ufft)
1353 :
1354 80 : CALL getIOMode(mode)
1355 :
1356 80 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1357 : #ifdef CPP_HDF
1358 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1359 80 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1360 80 : IF((judft_was_argument("-storeSF")).OR.(currentStepfunctionIndex.NE.0)) THEN
1361 0 : currentStepfunctionIndex = currentStepfunctionIndex + 1
1362 0 : CALL writeStepfunctionHDF(fileID, currentStepfunctionIndex, currentStarsIndex, currentStructureIndex, stars,.TRUE.)
1363 : CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1364 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1365 : END IF
1366 :
1367 80 : CALL closeCDNPOT_HDF(fileID)
1368 : #endif
1369 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
1370 : ! Write stars to stars file
1371 : STOP 'CDN_STREAM_MODE not yet implemented!'
1372 : ELSE
1373 : ! OPEN (14,file='wkf2',form='unformatted',status='unknown')
1374 :
1375 : ! WRITE (14) stars%ng3,ifftd
1376 : ! WRITE (14) (stars%ustep(i),i=1,stars%ng3)
1377 : ! WRITE (14) (stars%ufft(i),i=0,ifftd-1)
1378 :
1379 : ! CLOSE (14)
1380 : END IF
1381 :
1382 80 : END SUBROUTINE writeStepfunction
1383 :
1384 240 : SUBROUTINE readStepfunction(stars, atoms, cell, vacuum, l_error)
1385 :
1386 : TYPE(t_stars),INTENT(INOUT) :: stars
1387 : TYPE(t_atoms), INTENT(IN) :: atoms
1388 : TYPE(t_cell), INTENT(IN) :: cell
1389 : TYPE(t_vacuum), INTENT(IN) :: vacuum
1390 : LOGICAL, INTENT(OUT) :: l_error
1391 :
1392 : TYPE(t_stars) :: starsTemp
1393 : TYPE(t_input) :: inputTemp
1394 : TYPE(t_atoms) :: atomsTemp
1395 : TYPE(t_cell) :: cellTemp
1396 : TYPE(t_vacuum) :: vacuumTemp
1397 :
1398 :
1399 : INTEGER :: mode
1400 : INTEGER :: ifftd, ng3Temp, ifftdTemp, ioStatus, i, starsIndexTemp
1401 : INTEGER :: structureIndexTemp
1402 : LOGICAL :: l_exist, l_same, l_sameTemp
1403 :
1404 : INTEGER :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
1405 : INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
1406 : #ifdef CPP_HDF
1407 : INTEGER(HID_T) :: fileID
1408 : #endif
1409 :
1410 80 : l_error = .FALSE.
1411 80 : ioStatus = 0
1412 80 : ifftd = 27*stars%mx1*stars%mx2*stars%mx3
1413 :
1414 80 : CALL getIOMode(mode)
1415 :
1416 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1417 80 : INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
1418 80 : IF (l_exist) THEN
1419 : #ifdef CPP_HDF
1420 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1421 80 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1422 80 : IF(currentStepfunctionIndex.GT.0) THEN
1423 0 : CALL peekStepfunctionHDF(fileID, currentStepfunctionIndex, starsIndexTemp, structureIndexTemp)
1424 0 : IF((starsIndexTemp.EQ.currentStarsIndex).AND.(structureIndexTemp.EQ.currentStructureIndex)) THEN
1425 0 : CALL readStepfunctionHDF(fileID, currentStepfunctionIndex, stars)
1426 : ELSE
1427 : mode = CDN_STREAM_MODE ! No adequate stepfunction entry found. Fall back to other IO modes.
1428 : END IF
1429 : ELSE
1430 : mode = CDN_STREAM_MODE ! No adequate stepfunction entry found. Fall back to other IO modes.
1431 : END IF
1432 80 : CALL closeCDNPOT_HDF(fileID)
1433 : #endif
1434 : END IF
1435 80 : IF(.NOT.l_exist) THEN
1436 : mode = CDN_STREAM_MODE
1437 : END IF
1438 : END IF
1439 :
1440 80 : IF(mode.EQ.CDN_STREAM_MODE) THEN
1441 80 : INQUIRE(FILE='cdn.str',EXIST=l_exist)
1442 80 : IF (l_exist) THEN
1443 0 : STOP 'cdn.str code path not yet implemented!'
1444 : END IF
1445 : IF (.NOT.l_exist) THEN
1446 : mode = CDN_DIRECT_MODE
1447 : END IF
1448 : END IF
1449 :
1450 0 : IF (mode.EQ.CDN_DIRECT_MODE) THEN
1451 : ! INQUIRE(FILE='wkf2',EXIST=l_exist)
1452 : ! IF(.NOT.l_exist) THEN
1453 80 : l_error = .TRUE.
1454 80 : RETURN
1455 : ! END IF
1456 : ! OPEN (14,file='wkf2',form='unformatted',status='unknown')
1457 : ! ng3temp=0;ifftdTemp=0
1458 : ! READ (14,IOSTAT=ioStatus) ng3Temp, ifftdTemp
1459 : ! IF (ng3Temp.NE.stars%ng3) ioStatus = 1
1460 : ! IF (ifftdTemp.NE.ifftd) ioStatus = 1
1461 : ! IF (ioStatus.NE.0) THEN
1462 : ! l_error = .TRUE.
1463 : ! CLOSE (14)
1464 : ! RETURN
1465 : ! END IF
1466 : ! READ (14) (stars%ustep(i),i=1,stars%ng3)
1467 : ! READ (14) (stars%ufft(i),i=0,ifftd-1)
1468 :
1469 : ! CLOSE (14)
1470 : END IF
1471 :
1472 : END SUBROUTINE readStepfunction
1473 :
1474 80 : SUBROUTINE setStartingDensity(l_noco)
1475 :
1476 : LOGICAL,INTENT(IN) :: l_noco
1477 :
1478 : #ifdef CPP_HDF
1479 : INTEGER(HID_T) :: fileID
1480 : #endif
1481 : INTEGER :: currentStarsIndex,currentLatharmsIndex
1482 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
1483 : INTEGER :: readDensityIndex, lastDensityIndex
1484 : INTEGER :: sdIndex, ioStatus, mode
1485 : INTEGER :: densityType
1486 : CHARACTER(LEN=1000) :: numberString
1487 : CHARACTER(LEN=30) :: archiveName
1488 : LOGICAL :: l_exist
1489 :
1490 80 : IF (.NOT.juDFT_was_argument("-sd")) THEN
1491 80 : RETURN
1492 : END IF
1493 0 : numberString = juDFT_string_for_argument("-sd")
1494 0 : IF (TRIM(ADJUSTL(numberString)).EQ.'') THEN
1495 0 : CALL juDFT_error("No number for starting density set in command line.",calledby ="setStartingDensity")
1496 : END IF
1497 0 : ioStatus = 0
1498 0 : READ(numberString,'(i8)',iostat=ioStatus) sdIndex
1499 0 : IF(ioStatus.NE.0) THEN
1500 0 : CALL juDFT_error("Could not convert starting density index string to number.",calledby ="setStartingDensity")
1501 : END IF
1502 :
1503 0 : WRITE(*,'(a,i0,a)') 'Using density ', sdIndex, ' as starting density.'
1504 :
1505 0 : CALL getIOMode(mode)
1506 :
1507 0 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1508 : #ifdef CPP_HDF
1509 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1510 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1511 0 : densityType = DENSITY_TYPE_IN_const
1512 0 : IF(l_noco) THEN
1513 0 : densityType = DENSITY_TYPE_NOCO_IN_const
1514 : END IF
1515 0 : archiveName = ''
1516 0 : WRITE(archiveName,'(a,i0)') '/cdn-', sdIndex
1517 0 : l_exist = isDensityEntryPresentHDF(fileID,archiveName,densityType)
1518 0 : IF(.NOT.l_exist) THEN
1519 0 : WRITE(*,*) 'archiveName: ', TRIM(ADJUSTL(archiveName))
1520 0 : CALL juDFT_error("For selected starting density index no in-density is present.",calledby ="setStartingDensity")
1521 : END IF
1522 : CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1523 0 : currentStepfunctionIndex,sdIndex,lastDensityIndex)
1524 0 : CALL closeCDNPOT_HDF(fileID)
1525 : #endif
1526 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
1527 : STOP 'CDN_STREAM_MODE not yet implemented!'
1528 : ELSE
1529 0 : WRITE(*,*) 'Explicit setting of starting density in direct access mode'
1530 0 : WRITE(*,*) 'not implemented.'
1531 0 : WRITE(*,*) ''
1532 0 : WRITE(*,*) 'Ignoring -sd command line argument.'
1533 : END IF
1534 :
1535 : END SUBROUTINE setStartingDensity
1536 :
1537 80 : SUBROUTINE deleteDensities()
1538 :
1539 : #ifdef CPP_HDF
1540 : INTEGER(HID_T) :: fileID
1541 : #endif
1542 : INTEGER :: currentStarsIndex,currentLatharmsIndex
1543 : INTEGER :: currentStructureIndex,currentStepfunctionIndex
1544 : INTEGER :: readDensityIndex, lastDensityIndex
1545 : INTEGER :: ioStatus, mode, i
1546 : INTEGER :: startNumber, endNumber, separatorIndex
1547 : CHARACTER(LEN=1000) :: ddString
1548 : CHARACTER(LEN=30) :: archiveName
1549 : LOGICAL :: l_exist, l_deleted
1550 :
1551 80 : IF (.NOT.juDFT_was_argument("-delden")) THEN
1552 80 : RETURN
1553 : END IF
1554 0 : ddString = juDFT_string_for_argument("-delden")
1555 0 : IF (TRIM(ADJUSTL(ddString)).EQ.'') THEN
1556 0 : CALL juDFT_error("Densities to be deleted not specified.",calledby ="deleteDensities")
1557 : END IF
1558 :
1559 0 : separatorIndex = -1
1560 0 : startNumber = -1
1561 0 : endNumber = -1
1562 0 : IF (TRIM(ADJUSTL(ddString)).EQ.'allbutlast') THEN
1563 0 : CALL getIOMode(mode)
1564 :
1565 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1566 0 : INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
1567 0 : IF (l_exist) THEN
1568 : #ifdef CPP_HDF
1569 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1570 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1571 0 : CALL closeCDNPOT_HDF(fileID)
1572 0 : startNumber = 1
1573 0 : endNumber = lastDensityIndex - 1
1574 : #endif
1575 : END IF
1576 : END IF
1577 : ELSE
1578 0 : DO i = 1, LEN(TRIM(ADJUSTL(ddString)))
1579 0 : IF(VERIFY(ddString(i:i),'1234567890').NE.0) THEN
1580 0 : IF ((ddString(i:i).EQ.'-').AND.(separatorIndex.EQ.-1)) THEN
1581 : separatorIndex = i
1582 : ELSE
1583 0 : CALL juDFT_error("density deletion string format error",calledby ="deleteDensities")
1584 : END IF
1585 : END IF
1586 : END DO
1587 :
1588 0 : IF(separatorIndex.NE.-1) THEN
1589 0 : READ(ddString(1:separatorIndex-1),'(i7)') startNumber
1590 0 : READ(ddString(separatorIndex+1:LEN(TRIM(ADJUSTL(ddString)))),'(i7)') endNumber
1591 : ELSE
1592 0 : READ(ddString(1:LEN(TRIM(ADJUSTL(ddString)))),'(i7)') startNumber
1593 0 : READ(ddString(1:LEN(TRIM(ADJUSTL(ddString)))),'(i7)') endNumber
1594 : END IF
1595 : END IF
1596 :
1597 0 : CALL getIOMode(mode)
1598 :
1599 : IF(mode.EQ.CDN_HDF5_MODE) THEN
1600 0 : INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
1601 0 : IF (l_exist) THEN
1602 : #ifdef CPP_HDF
1603 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1604 0 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1605 :
1606 0 : DO i = startNumber, endNumber
1607 0 : archiveName = ''
1608 0 : WRITE(archiveName,'(a,i0)') '/cdn-', i
1609 :
1610 0 : l_exist = isDensityEntryPresentHDF(fileID,archiveName,DENSITY_TYPE_UNDEFINED_const)
1611 0 : IF(.NOT.l_exist) THEN
1612 : CYCLE
1613 : END IF
1614 :
1615 0 : l_deleted = deleteDensityEntryHDF(fileID,archiveName)
1616 0 : IF (l_deleted) THEN
1617 0 : WRITE(*,*) 'deleted density entry ', TRIM(ADJUSTL(archiveName))
1618 : END IF
1619 : END DO
1620 :
1621 : CALL deleteObsoleteDensityMetadataHDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1622 0 : currentStepfunctionIndex,lastDensityIndex)
1623 :
1624 0 : CALL closeCDNPOT_HDF(fileID)
1625 : #endif
1626 0 : WRITE(*,*) 'Please note:'
1627 0 : WRITE(*,*) 'The deletion of the densities does not free the associated disk space.'
1628 0 : WRITE(*,*) 'To do this you have to repack the cdn.hdf file.'
1629 0 : WRITE(*,*) 'It can be done by using the tool h5repack, e.g., by invoking'
1630 0 : WRITE(*,*) 'h5repack -i cdn.hdf -o cdn-packed.hdf'
1631 0 : WRITE(*,*) 'mv cdn-packed.hdf cdn.hdf'
1632 : ELSE
1633 0 : WRITE(*,*) "No cdn.hdf file found. No density entry deleted."
1634 : END IF
1635 : ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
1636 : STOP 'CDN_STREAM_MODE not yet implemented!'
1637 : ELSE
1638 0 : WRITE(*,*) 'Explicit deletion of densities in direct access mode'
1639 0 : WRITE(*,*) 'not implemented.'
1640 0 : WRITE(*,*) ''
1641 0 : WRITE(*,*) 'Ignoring -delden command line argument.'
1642 : END IF
1643 :
1644 0 : CALL juDFT_end("Selected densities deleted.")
1645 :
1646 : END SUBROUTINE deleteDensities
1647 :
1648 1501 : SUBROUTINE getIOMode(mode)
1649 : INTEGER, INTENT(OUT) :: mode
1650 :
1651 1501 : mode = CDN_DIRECT_MODE
1652 : !IF (juDFT_was_argument("-stream_cdn")) THEN
1653 : ! mode=CDN_STREAM_MODE
1654 : !END IF
1655 1501 : IF (.NOT.juDFT_was_argument("-no_cdn_hdf")) THEN !juDFT_was_argument("-hdf_cdn")) THEN
1656 : #ifdef CPP_HDF
1657 668 : mode=CDN_HDF5_MODE
1658 : #else
1659 : ! WRITE(*,*) 'Code not compiled with HDF5 support.'
1660 : ! WRITE(*,*) 'Falling back to direct access.'
1661 : #endif
1662 : END IF
1663 0 : END SUBROUTINE getIOMode
1664 :
1665 228 : LOGICAL FUNCTION isDensityFilePresent(archiveType)
1666 :
1667 : INTEGER, INTENT(IN) :: archiveType
1668 :
1669 : LOGICAL :: l_exist
1670 : INTEGER :: mode
1671 :
1672 : INTEGER :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
1673 : INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
1674 : #ifdef CPP_HDF
1675 : INTEGER(HID_T) :: fileID
1676 : #endif
1677 :
1678 80 : CALL getIOMode(mode)
1679 :
1680 : IF (mode.EQ.CDN_HDF5_MODE) THEN
1681 80 : INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
1682 80 : IF(l_exist) THEN
1683 : #ifdef CPP_HDF
1684 : CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
1685 80 : currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
1686 80 : CALL closeCDNPOT_HDF(fileID)
1687 80 : IF(readDensityIndex.GT.0) THEN
1688 12 : isDensityFilePresent = .TRUE.
1689 : RETURN
1690 : END IF
1691 : #endif
1692 : END IF
1693 : END IF
1694 :
1695 68 : IF ((mode.EQ.CDN_STREAM_MODE).OR.(mode.EQ.CDN_HDF5_MODE)) THEN
1696 68 : INQUIRE(FILE='cdn.str',EXIST=l_exist)
1697 68 : IF(l_exist) THEN
1698 12 : isDensityFilePresent = l_exist
1699 : RETURN
1700 : END IF
1701 : END IF
1702 :
1703 : !cdn1 or rhomat_inp should be enough for any mode...
1704 68 : INQUIRE(FILE='cdn1',EXIST=l_exist)
1705 68 : IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const) THEN
1706 45 : isDensityFilePresent = l_exist
1707 45 : RETURN
1708 : END IF
1709 23 : IF ((archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const).AND.(archiveType.NE.CDN_ARCHIVE_TYPE_FFN_const)) THEN
1710 0 : CALL juDFT_error("Illegal archive type selected.",calledby ="isDensityFilePresent")
1711 : END IF
1712 23 : IF (l_exist) THEN
1713 12 : isDensityFilePresent = l_exist
1714 : RETURN
1715 : END IF
1716 23 : INQUIRE(FILE='rhomat_inp',EXIST=l_exist)
1717 23 : isDensityFilePresent = l_exist
1718 23 : END FUNCTION isDensityFilePresent
1719 :
1720 0 : LOGICAL FUNCTION isCoreDensityPresent()
1721 :
1722 : LOGICAL :: l_exist
1723 : INTEGER :: mode
1724 :
1725 0 : CALL getIOMode(mode)
1726 :
1727 : IF (mode.EQ.CDN_HDF5_MODE) THEN
1728 : #ifdef CPP_HDF
1729 0 : l_exist = isCoreDensityPresentHDF()
1730 0 : IF(l_exist) THEN
1731 : isCoreDensityPresent = l_exist
1732 : RETURN
1733 : END IF
1734 : #endif
1735 : END IF
1736 :
1737 0 : IF ((mode.EQ.CDN_STREAM_MODE).OR.(mode.EQ.CDN_HDF5_MODE)) THEN
1738 0 : INQUIRE(FILE='cdn.str',EXIST=l_exist)
1739 0 : IF(l_exist) THEN
1740 0 : STOP 'Not yet implemented!'
1741 : RETURN
1742 : END IF
1743 : END IF
1744 :
1745 : !cdnc should be enough for any mode...
1746 0 : INQUIRE(FILE='cdnc',EXIST=l_exist)
1747 0 : IF (l_exist) THEN
1748 : isCoreDensityPresent = l_exist
1749 : RETURN
1750 : END IF
1751 : isCoreDensityPresent = .FALSE.
1752 : END FUNCTION isCoreDensityPresent
1753 :
1754 : END MODULE m_cdn_io
|
