First steps towards wrapping charge density I/O

The interface is by far not yet complete. Up to now the charge
density I/O operations in main/vgen.F90 are wrapped.

io/CMakeLists.txt

@@ -14,6 +14,7 @@ io/eig66_mpi.F90
 io/nocoInputCheck.F90
 io/inpnoco.F90
 io/loddop.f90
+io/cdn_io.f90
 io/rw_inp.f90
 io/rw_noco.f90
 io/r_inpXML.F90

io/cdn_io.f90

+!--------------------------------------------------------------------------------
+! 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.
+!--------------------------------------------------------------------------------
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! This module is a wrapper for the charge density I/O
+!!!
+!!!                             GM'17
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+MODULE m_cdn_io
+
+   USE m_types
+   USE m_juDFT
+   USE m_loddop
+   USE m_wrtdop
+   IMPLICIT NONE
+
+   PRIVATE
+   PUBLIC readDensity, writeDensity
+   PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
+   PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
+   PUBLIC CDN_ARCHIVE_TYPE_CDN_const
+
+   INTEGER,          PARAMETER :: CDN_INPUT_DEN_const = 1
+   INTEGER,          PARAMETER :: CDN_OUTPUT_DEN_const = 2
+
+   INTEGER,          PARAMETER :: CDN_ARCHIVE_TYPE_CDN1_const = 1
+   INTEGER,          PARAMETER :: CDN_ARCHIVE_TYPE_NOCO_const = 2
+   INTEGER,          PARAMETER :: CDN_ARCHIVE_TYPE_CDN_const  = 3
+
+   INTEGER,          PARAMETER :: CDN_DIRECT_MODE = 1
+   INTEGER,          PARAMETER :: CDN_STREAM_MODE = 2
+   INTEGER,          PARAMETER :: CDN_HDF5_MODE   = 3
+
+   CONTAINS
+
+   SUBROUTINE readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
+                          relCdnIndex,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
+
+      TYPE(t_stars),INTENT(IN)  :: stars
+      TYPE(t_vacuum),INTENT(IN) :: vacuum
+      TYPE(t_atoms),INTENT(IN)  :: atoms
+      TYPE(t_sphhar),INTENT(IN) :: sphhar
+      TYPE(t_input),INTENT(IN)  :: input
+      TYPE(t_sym),INTENT(IN)    :: sym
+      TYPE(t_oneD),INTENT(IN)   :: oneD
+
+      INTEGER, INTENT (IN)      :: inOrOutCDN
+      INTEGER, INTENT (IN)      :: relCdnIndex
+      INTEGER, INTENT (OUT)     :: iter
+      INTEGER, INTENT (IN)      :: archiveType
+
+      !     ..
+      !     .. Array Arguments ..
+      COMPLEX, INTENT (OUT) :: fpw(stars%n3d,input%jspins), fzxy(vacuum%nmzxyd,stars%n2d-1,2,input%jspins)
+      COMPLEX, INTENT (OUT) :: cdom(:), cdomvz(:,:), cdomvxy(:,:,:)
+      REAL,    INTENT (OUT) :: fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins), fz(vacuum%nmzd,2,input%jspins)
+
+      ! local variables
+      INTEGER           :: mode, datend, k, i, iVac, j, iUnit
+      LOGICAL           :: l_exist
+      CHARACTER(len=30) :: filename
+
+      CALL getMode(mode)
+
+      IF(mode.EQ.CDN_HDF5_MODE) THEN
+         INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
+         IF (l_exist) THEN
+            !load density from cdn.hdf and exit subroutine
+
+            RETURN
+         ELSE
+            WRITE(*,*) 'cdn.hdf file not found.'
+            WRITE(*,*) 'Falling back to stream access file cdn.str.'
+            mode = CDN_STREAM_MODE
+         END IF
+      END IF
+
+      IF(mode.EQ.CDN_STREAM_MODE) THEN
+         INQUIRE(FILE='cdn.str',EXIST=l_exist)
+         IF (l_exist) THEN
+            !load density from cdn.str and exit subroutine
+
+            RETURN
+         ELSE
+            WRITE(*,*) 'cdn.str file not found.'
+            WRITE(*,*) 'Falling back to direct access file cdn1.'
+            mode = CDN_DIRECT_MODE
+         END IF
+      END IF
+
+      IF (mode.EQ.CDN_DIRECT_MODE) THEN
+         filename = 'cdn1'
+         IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
+            INQUIRE(file="rhomat_inp",EXIST=l_exist)
+            IF (l_exist) filename = 'rhomat_inp'
+         END IF
+         IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
+            filename = 'cdn'
+         END IF
+
+         INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
+         IF(.NOT.l_exist) THEN
+            CALL juDFT_error("charge density file missing",calledby ="readDensity")
+         END IF
+
+         iUnit = 93
+         OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='old')
+
+         IF (inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) THEN
+            IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
+               CALL juDFT_error("inOrOutCDN.EQ.CDN_OUTPUT_DEN_const incompatible to l_noco",calledby ="readDensity")
+            END IF
+            ! call loddop to move the file position to the output density
+            CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
+                        iUnit,iter,fr,fpw,fz,fzxy)
+         END IF
+
+         ! read in the density
+         CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
+                     iUnit,iter,fr,fpw,fz,fzxy)
+
+         ! read in additional data if l_noco and data is present
+         IF ((archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const).AND.(TRIM(ADJUSTL(filename)).EQ.'rhomat_inp')) THEN
+            READ (iUnit,iostat=datend) (cdom(k),k=1,stars%ng3)
+            IF (datend == 0) THEN
+               IF (input%film) THEN
+                  READ (iUnit) ((cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
+                  READ (iUnit) (((cdomvxy(i,j-1,iVac),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), iVac=1,vacuum%nvac)
+               END IF
+            ELSE
+               ! (datend < 0)  =>  no off-diagonal magnetisation stored
+               !                   in "rhomat_inp"
+               IF (datend > 0) THEN
+                  WRITE(*,*) 'datend = ', datend
+                  CALL juDFT_error("density file has illegal state.",calledby ="readDensity")
+               END IF
+               cdom = CMPLX(0.0,0.0)
+               IF (input%film) THEN
+                  cdomvz = CMPLX(0.0,0.0)
+                  cdomvxy = CMPLX(0.0,0.0)
+               END IF
+            END IF
+         ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
+            cdom = CMPLX(0.0,0.0)
+            IF (input%film) THEN
+               cdomvz = CMPLX(0.0,0.0)
+               cdomvxy = CMPLX(0.0,0.0)
+            END IF
+         END IF
+      END IF
+
+      CLOSE(iUnit)
+
+   END SUBROUTINE readDensity
+
+   SUBROUTINE writeDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
+                           iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
+
+      TYPE(t_stars),INTENT(IN)  :: stars
+      TYPE(t_vacuum),INTENT(IN) :: vacuum
+      TYPE(t_atoms),INTENT(IN)  :: atoms
+      TYPE(t_sphhar),INTENT(IN) :: sphhar
+      TYPE(t_input),INTENT(IN)  :: input
+      TYPE(t_sym),INTENT(IN)    :: sym
+      TYPE(t_oneD),INTENT(IN)   :: oneD
+
+      INTEGER, INTENT (IN)      :: inOrOutCDN
+      INTEGER, INTENT (IN)      :: iter
+      INTEGER, INTENT (IN)      :: archiveType
+      !     ..
+      !     .. Array Arguments ..
+      COMPLEX, INTENT (IN) :: fpw(stars%n3d,input%jspins), fzxy(vacuum%nmzxyd,stars%n2d-1,2,input%jspins)
+      COMPLEX, INTENT (IN) :: cdom(:), cdomvz(:,:), cdomvxy(:,:,:)
+      REAL,    INTENT (IN) :: fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins), fz(vacuum%nmzd,2,input%jspins)
+
+      TYPE(t_stars)        :: starsTemp
+      TYPE(t_vacuum)       :: vacuumTemp
+      TYPE(t_atoms)        :: atomsTemp
+      TYPE(t_sphhar)       :: sphharTemp
+      TYPE(t_input)        :: inputTemp
+      TYPE(t_sym)          :: symTemp
+
+      COMPLEX, ALLOCATABLE :: fpwTemp(:,:), fzxyTemp(:,:,:,:)
+      REAL, ALLOCATABLE    :: frTemp(:,:,:,:), fzTemp(:,:,:)
+
+      INTEGER           :: mode, iterTemp, k, i, iVac, j, iUnit
+      CHARACTER(len=30) :: filename
+
+      CALL getMode(mode)
+
+      IF(mode.EQ.CDN_HDF5_MODE) THEN
+         ! Write density to cdn.hdf file
+         STOP 'CDN_HDF5_MODE not yet implemented!'
+      ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
+         ! Write density to cdn.str file
+         STOP 'CDN_STREAM_MODE not yet implemented!'
+      ELSE
+         filename = 'cdn1'
+         IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
+            filename = 'rhomat_inp'
+         END IF
+         IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
+            filename = 'cdn'
+         END IF
+
+         iUnit = 93
+         OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='unknown')
+
+         IF (inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) THEN
+            IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
+               CALL juDFT_error("inOrOutCDN.EQ.CDN_OUTPUT_DEN_const incompatible to l_noco",calledby ="writeDensity")
+            END IF
+
+            ! Generate data in temp arrays and variables to be able to perform loddop call.
+            ! loddop is called to move the file position to the output density position.
+            starsTemp%n3d = stars%n3d
+            inputTemp%jspins = input%jspins
+            vacuumTemp%nmzxyd = vacuum%nmzxyd
+            starsTemp%n2d = stars%n2d
+            atomsTemp%jmtd = atoms%jmtd
+            sphharTemp%nlhd = sphhar%nlhd
+            vacuumTemp%nmzd = vacuum%nmzd
+            atomsTemp%ntype = atoms%ntype
+            ALLOCATE (sphharTemp%nlh(SIZE(sphhar%nlh)))
+            sphharTemp%nlh(:) = sphhar%nlh(:)
+            ALLOCATE (atomsTemp%ntypsy(SIZE(atoms%ntypsy)))
+            atomsTemp%ntypsy(:) = atoms%ntypsy(:)
+            ALLOCATE (atomsTemp%jri(SIZE(atoms%jri)))
+            atomsTemp%jri(:) = atoms%jri(:)
+            ALLOCATE (atomsTemp%neq(SIZE(atoms%neq)))
+            atomsTemp%neq(:) = atoms%neq(:)
+            starsTemp%ng3 = stars%ng3
+            symTemp%invs = sym%invs
+            inputTemp%film = input%film
+            vacuumTemp%nvac = vacuum%nvac
+            starsTemp%ng2 = stars%ng2
+            symTemp%invs2 = sym%invs2
+            ALLOCATE (fpwTemp(stars%n3d,input%jspins))
+            ALLOCATE (fzxyTemp(vacuum%nmzxyd,stars%n2d-1,2,input%jspins))
+            ALLOCATE (frTemp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins))
+            ALLOCATE (fzTemp(vacuum%nmzd,2,input%jspins))
+
+            CALL loddop(starsTemp,vacuumTemp,atomsTemp,sphharTemp,inputTemp,symTemp,&
+                        iUnit,iterTemp,frTemp,fpwTemp,fzTemp,fzxyTemp)
+
+            DEALLOCATE (fzTemp, frTemp, fzxyTemp, fpwTemp)
+            DEALLOCATE (atomsTemp%neq, atomsTemp%jri, atomsTemp%ntypsy, sphharTemp%nlh)
+         END IF
+
+         ! Write the density
+         CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
+                     iUnit,iter,fr,fpw,fz,fzxy)
+
+         ! Write additional data if l_noco
+         IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
+            WRITE (iUnit) (cdom(k),k=1,stars%ng3)
+            IF (input%film) THEN
+               WRITE (iUnit) ((cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
+               WRITE (iUnit) (((cdomvxy(i,j-1,iVac),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), iVac=1,vacuum%nvac)
+            END IF
+         END IF
+
+         CLOSE(iUnit)
+      END IF
+
+   END SUBROUTINE writeDensity
+
+   SUBROUTINE getMode(mode)
+      INTEGER, INTENT(OUT) :: mode
+
+      mode = CDN_DIRECT_MODE
+      IF (juDFT_was_argument("-stream_cdn")) mode=CDN_STREAM_MODE
+      IF (juDFT_was_argument("-hdf_cdn")) mode=CDN_HDF5_MODE
+   END SUBROUTINE getMode
+
+END MODULE m_cdn_io

main/optional.f90

@@ -137,9 +137,10 @@ CONTAINS
     IF (.NOT.(strho.OR.obsolete%l_f2u.OR.obsolete%l_u2f.OR.sliceplot%iplot)) THEN
        IF (noco%l_noco) THEN
           INQUIRE (file='rhomat_inp',exist=strho) ! if no density (rhoma
-       ELSE
+       END IF
+       IF(.NOT.strho) THEN
           INQUIRE (file='cdn1',exist=strho)       ! if no density (cdn1)
-       ENDIF
+       END IF
        strho = .NOT.strho                ! create a starting density
     ENDIF
 

main/vgen.F90

@@ -32,8 +32,8 @@ CONTAINS
     USE m_vvacxy
     USE m_vintcz
     USE m_checkdop
-    USE m_loddop
     USE m_wrtdop
+    USE m_cdn_io
     USE m_qfix
     USE m_types
     USE m_od_vvac
@@ -73,9 +73,9 @@ CONTAINS
     !     .. Local Scalars ..
     COMPLEX vintcza,xint,rhobar
     INTEGER i,i3,irec2,irec3,ivac,j,js,k,k3,lh,n,nzst1
-    INTEGER imz,imzxy,ichsmrg,ivfft,nt,npd 
+    INTEGER imz,imzxy,ichsmrg,ivfft,npd 
     INTEGER ifftd,ifftd2, ifftxc3d,iter,datend
-    INTEGER itypsym,itype,jsp,l,nat
+    INTEGER itypsym,itype,jsp,l,nat,archiveType
     !      INTEGER i_sm,n_sm,i_sta,i_end
     REAL ani,g3,signum,z,rhmn,fix,mfie
     REAL sig1dh,vz1dh,zat_l(atoms%ntype),rdum,dpdot ! ,delta,deltb,corr
@@ -133,26 +133,14 @@ CONTAINS
 
     IF (noco%l_noco) THEN
        ALLOCATE ( cdom(stars%n3d), cdomvz(vacuum%nmzd,2),cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2) )
+       archiveType = CDN_ARCHIVE_TYPE_NOCO_const
     ELSE
        ALLOCATE ( cdom(1),cdomvz(1,1),cdomvxy(1,1,1) )
+       archiveType = CDN_ARCHIVE_TYPE_CDN1_const
     END IF
     !
 
     IF (mpi%irank == 0) THEN
-       ! ff
-       IF (noco%l_noco) THEN
-          ! nt = 70
-          ! OPEN (nt,file='cdn',form='unformatted',status='old')
-          ! In the non-collinear case |m| is calculated from mx,my,mz
-          ! in "visxc(g)","vvacxc(g)" instead of using the "cdn"-file.
-          ! It is done this way as accuracy is lost when transforming
-          ! |m| back to g-space in "rhodirgen" .
-          nt = 26
-          OPEN (nt,file='rhomat_inp',form='unformatted',status='old')
-       ELSE
-          nt = 71
-          OPEN (nt,file='cdn1',form='unformatted',status='old')
-       ENDIF
        !
        ! --  total = .false. and reap = .false. means, that we now calculate
        !     the potential from the output density
@@ -161,33 +149,11 @@ CONTAINS
           IF (noco%l_noco) THEN
              CALL juDFT_error("vgen:1",calledby ="vgen")
           ENDIF
-          CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
-               nt, iter,rho,qpw,rht,rhtxy)
-       ENDIF
-       !
-       ! --> load density and fix
-       !
-       CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
-            nt, iter,rho,qpw,rht,rhtxy)
-       IF (noco%l_noco) THEN
-          READ (nt,iostat=datend) (cdom(k),k=1,stars%ng3)
-          IF (datend == 0) THEN
-             IF (input%film) THEN
-                READ (nt) ((cdomvz(i,ivac),i=1,vacuum%nmz),ivac=1,vacuum%nvac)
-                READ (nt) (((cdomvxy(i,j-1,ivac),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), ivac=1,vacuum%nvac)
-             END IF
-          ELSE
-             ! (datend < 0)  =>  no off-diagonal magnetisation stored
-             !                   in "rhomat_inp"
-             IF (datend > 0) THEN
-                CALL juDFT_error("vgen: error reading ",calledby ="vgen")
-             END IF
-             cdom= CMPLX(0.,0.)
-             IF (input%film) THEN
-                cdomvz= CMPLX(0.,0.)
-                cdomvxy= CMPLX(0.,0.)
-             END IF
-          END IF
+          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,CDN_OUTPUT_DEN_const,&
+                           0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
+       ELSE
+          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
+                           0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
        END IF
 
        IF (.NOT.l_xyav) THEN
@@ -195,17 +161,12 @@ CONTAINS
           CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD, qpw,rhtxy,rho,rht,.FALSE., fix)
           CALL timestop("Qfix")
        ENDIF
-       IF (input%total.OR.reap) REWIND nt
-       CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
-            nt, iter,rho,qpw,rht,rhtxy)
-       IF (noco%l_noco) THEN
-          WRITE (nt) (cdom(k),k=1,stars%ng3)
-          IF (input%film) THEN
-             WRITE (nt) ((cdomvz(i,ivac),i=1,vacuum%nmz),ivac=1,vacuum%nvac)
-             WRITE (nt) (((cdomvxy(i,j-1,ivac),i=1,vacuum%nmzxy),j=2,oneD%odi%nq2), ivac=1,vacuum%nvac)
-          ENDIF
-       ENDIF
-       !
+
+       IF (input%total.OR.reap) THEN
+          CALL writeDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
+                            iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
+       END IF
+
        WRITE (6,FMT=8000)
 8000   FORMAT (/,/,t10,' p o t e n t i a l   g e n e r a t o r',/)
        vpw  = CMPLX(0.,0.)
@@ -439,19 +400,14 @@ CONTAINS
        !     ----> reload the density for calculating vxc (for spin-pol. case)
        !
        IF (input%jspins.EQ.2) THEN
-          !
-          REWIND nt
-          CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
-               nt, iter,rho,qpw,rht,rhtxy)
-          CLOSE (nt)
+          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
+                           0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
           vr(:,0:,:,2) = vr(:,0:,:,1)
           vpw(:,2) = vpw(:,1)
           IF (input%film) THEN
              vxy(:,:,:,2) = vxy(:,:,:,1)
              vz(:,:,2)=vz(:,:,1)
           END IF
-       ELSE
-          CLOSE (nt)
        END IF
        IF (input%total) THEN
           OPEN (11,file='potcoul',form='unformatted',status='unknown')
@@ -707,11 +663,8 @@ CONTAINS
        IF (input%total) THEN
 
           IF (noco%l_noco) THEN ! load qpw,rht,rhtxy from 'cdn'-file
-             nt = 70
-             OPEN (nt,file='cdn',form='unformatted',status='old')
-             CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
-                  nt, iter,rho,qpw,rht,rhtxy)
-             CLOSE (nt)
+             CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
+                              0,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
           ENDIF
           !
           !     CALCULATE THE INTEGRAL OF n1*Veff1 + n2*Veff2