From 3218c4ff70eeabcd9bf37947f25a98cc965e514c Mon Sep 17 00:00:00 2001 From: ncolonna Date: Thu, 17 Aug 2023 22:03:55 +0200 Subject: [PATCH 01/11] Add n2npm1 PP program. First commit --- PP/Makefile | 7 ++- PP/n2npm1.f90 | 149 ++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 155 insertions(+), 1 deletion(-) create mode 100644 PP/n2npm1.f90 diff --git a/PP/Makefile b/PP/Makefile index 6a4c3ab2..95e656db 100644 --- a/PP/Makefile +++ b/PP/Makefile @@ -77,7 +77,12 @@ MODULES = \ ../Modules/wrappers.o \ ../Modules/compute_dipole.o -all: merge_evc.x bin2xml.x +all: merge_evc.x bin2xml.x n2npm1.x + + +n2npm1.x: n2npm1.o + $(LD) $(LDFLAGS) -o $@ n2npm1.o ../iotk/src/libiotk.a + - ( cd ../bin ; ln -fs ../PP/$@ . ) bin2xml.x: bin2xml.o $(LD) $(LDFLAGS) -o $@ bin2xml.o $(MODULES) $(LIBOBJS) $(LIBS) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 new file mode 100644 index 00000000..2a431cf7 --- /dev/null +++ b/PP/n2npm1.f90 @@ -0,0 +1,149 @@ +! TO DO: +! 1) this is the desired comand with arguments from command line +! --fill tmpdir_in tmpdir_out +! --empty tmpdir_in tmpdir_out +! 2) Add reading routine for empty state (only if --empty) +! 3) Check Makefile and the need of mpif90 wrapper (it loooks like we need it +! as iotk as some dependendence on mpi library. Also we might need a different +! rule to compile this program (no need for any precompilation flags like e.g. -D__MPI) +PROGRAM n2npm1 + !------------------------------------------------------------------------ + ! + IMPLICIT NONE + CHARACTER(LEN=256) :: filename + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + COMPLEX(DP), ALLOCATABLE :: evc(:,:,:) + INTEGER :: ngw, nbnd(2), ispin, nspin + LOGICAL :: is_cmplx + ! + ngw = 0; nbnd=0; nspin=0 + filename='evc01.dat' ! FIXME this will be dependent on the input path + ! + CALL read_sizes (filename, ngw, nbnd(1), nspin, is_cmplx) + ! + IF (nspin==2) THEN + filename='evc02.dat' + CALL read_sizes (filename, ngw, nbnd(2), nspin, is_cmplx) + ENDIF + ! + WRITE(*,*) ngw, nbnd(:), nspin + ALLOCATE (evc(ngw, max(nbnd(1), nbnd(2)), nspin)) + ! + ! reset filename + filename='evc01.dat' + DO ispin = 1, nspin + IF (ispin==2 ) filename="evc02.dat" + CALL read_wf_occ (filename, evc(:,:,ispin), ngw, nbnd(ispin)) + WRITE(*,*) ispin, evc(1:3,1,ispin) ! check debug + ENDDO + ! + CONTAINS + ! + ! + SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx) + !------------------------------------------------------------------------ + ! + USE iotk_module + ! + IMPLICIT NONE + CHARACTER(LEN=256), INTENT(IN) :: filename + INTEGER, INTENT(OUT) :: ngw, nbnd, nspin + LOGICAL, INTENT(OUT) :: is_cmplx + + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + CHARACTER(iotk_attlenx) :: attr + INTEGER :: j + INTEGER :: ierr + INTEGER :: igwx, ik, nk + INTEGER :: iuni + ! + ierr = 0 + iuni=789 + ! + CALL iotk_open_read( iuni, FILE = filename, & + BINARY = .TRUE., IERR = ierr ) + ! + IF (ierr /= 0) THEN + WRITE(*,*) "Something wrong while reading file = ", filename + STOP + ENDIF + ! + CALL iotk_scan_empty( iuni, "INFO", attr ) + ! + WRITE(*,*) "NICOLA", ngw, nbnd, ik, nk, nspin, igwx, is_cmplx + CALL iotk_scan_attr( attr, "ngw", ngw ) + CALL iotk_scan_attr( attr, "nbnd", nbnd ) + CALL iotk_scan_attr( attr, "ik", ik ) + CALL iotk_scan_attr( attr, "nk", nk ) + CALL iotk_scan_attr( attr, "nspin", nspin ) + CALL iotk_scan_attr( attr, "igwx", igwx ) + CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) + ! + WRITE(*,*) "NICOLA", ngw, nbnd, ik, nk, nspin, igwx, is_cmplx + CALL iotk_close_read( iuni ) + RETURN + ! + END subroutine + ! + ! + SUBROUTINE read_wf_occ(filename, wf, ngw, nbnd) + !------------------------------------------------------------------------ + ! + USE iotk_module + ! + IMPLICIT NONE + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + CHARACTER(LEN=256), INTENT(IN) :: filename + COMPLEX(DP) :: wf(ngw, nbnd) + INTEGER, INTENT(IN) :: ngw, nbnd + CHARACTER(iotk_attlenx) :: attr + INTEGER :: j + COMPLEX(DP), ALLOCATABLE :: wtmp(:) + INTEGER :: ierr + INTEGER :: igwx, ngw_, nbnd_ + INTEGER :: iuni + LOGICAL :: is_cmplx + ! + ierr = 0 + iuni=789 + ! + CALL iotk_open_read( iuni, FILE = filename, & + BINARY = .TRUE., IERR = ierr ) + ! + IF (ierr /= 0) THEN + WRITE(*,*) "Something wrong while reading file = ", filename + STOP + ENDIF + ! + CALL iotk_scan_empty( iuni, "INFO", attr ) + ! + ngw_=0; nbnd_=0; + WRITE(*,*) "NICOLA", ngw_, nbnd_ + CALL iotk_scan_attr( attr, "ngw", ngw_ ) + CALL iotk_scan_attr( attr, "nbnd", nbnd_ ) + CALL iotk_scan_attr( attr, "igwx", igwx ) + CALL iotk_scan_attr( attr, "ispin", ispin ) + CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) + ! + WRITE(*,*) "NICOLA", ngw_, nbnd_, ispin, igwx, is_cmplx + ! + IF (nbnd_ /= nbnd .OR. ngw_ /= ngw) THEN + WRITE(*,*) "Size Mismatch. STOP" + STOP + ENDIF + + ALLOCATE( wtmp( igwx ) ) + ! + DO j = 1, nbnd + CALL iotk_scan_dat( iuni, "evc" // iotk_index( j ), wtmp(1:igwx) ) + wf(:,j) = wtmp (:) + ENDDO + ! + DEALLOCATE (wtmp) + ! + CALL iotk_close_read( iuni ) + RETURN + ! + END subroutine + ! +END PROGRAM From eaad7d7e3b2a122d9a5a102d753c087bdcea2b59 Mon Sep 17 00:00:00 2001 From: ncolonna Date: Thu, 17 Aug 2023 22:34:33 +0200 Subject: [PATCH 02/11] proper allocation of the wfc --- PP/n2npm1.f90 | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index 2a431cf7..2637fa7e 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -13,48 +13,48 @@ PROGRAM n2npm1 CHARACTER(LEN=256) :: filename INTEGER, PARAMETER :: DP = selected_real_kind(14,200) COMPLEX(DP), ALLOCATABLE :: evc(:,:,:) - INTEGER :: ngw, nbnd(2), ispin, nspin + INTEGER :: ngw, nbnd(2), ispin, nspin, igwx LOGICAL :: is_cmplx ! ngw = 0; nbnd=0; nspin=0 filename='evc01.dat' ! FIXME this will be dependent on the input path ! - CALL read_sizes (filename, ngw, nbnd(1), nspin, is_cmplx) + CALL read_sizes (filename, ngw, nbnd(1), nspin, is_cmplx, igwx) ! IF (nspin==2) THEN filename='evc02.dat' - CALL read_sizes (filename, ngw, nbnd(2), nspin, is_cmplx) + CALL read_sizes (filename, ngw, nbnd(2), nspin, is_cmplx, igwx) ENDIF ! - WRITE(*,*) ngw, nbnd(:), nspin - ALLOCATE (evc(ngw, max(nbnd(1), nbnd(2)), nspin)) + WRITE(*,*) igwx, ngw, nbnd(:), nspin + ALLOCATE (evc(igwx, max(nbnd(1), nbnd(2)), nspin)) ! ! reset filename filename='evc01.dat' DO ispin = 1, nspin IF (ispin==2 ) filename="evc02.dat" CALL read_wf_occ (filename, evc(:,:,ispin), ngw, nbnd(ispin)) - WRITE(*,*) ispin, evc(1:3,1,ispin) ! check debug + WRITE(*,*) ispin, evc(1,1,ispin), evc(igwx,nbnd(ispin),ispin) ! check debug ENDDO ! CONTAINS ! ! - SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx) + SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) !------------------------------------------------------------------------ ! USE iotk_module ! IMPLICIT NONE CHARACTER(LEN=256), INTENT(IN) :: filename - INTEGER, INTENT(OUT) :: ngw, nbnd, nspin + INTEGER, INTENT(OUT) :: ngw, nbnd, nspin, igwx LOGICAL, INTENT(OUT) :: is_cmplx INTEGER, PARAMETER :: DP = selected_real_kind(14,200) CHARACTER(iotk_attlenx) :: attr - INTEGER :: j + INTEGER :: j, ispin INTEGER :: ierr - INTEGER :: igwx, ik, nk + INTEGER :: ik, nk INTEGER :: iuni ! ierr = 0 @@ -70,16 +70,17 @@ SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx) ! CALL iotk_scan_empty( iuni, "INFO", attr ) ! - WRITE(*,*) "NICOLA", ngw, nbnd, ik, nk, nspin, igwx, is_cmplx + !WRITE(*,*) "NICOLA", ngw, nbnd, ik, nk, nspin, igwx, is_cmplx CALL iotk_scan_attr( attr, "ngw", ngw ) CALL iotk_scan_attr( attr, "nbnd", nbnd ) CALL iotk_scan_attr( attr, "ik", ik ) CALL iotk_scan_attr( attr, "nk", nk ) CALL iotk_scan_attr( attr, "nspin", nspin ) + CALL iotk_scan_attr( attr, "ispin", ispin ) CALL iotk_scan_attr( attr, "igwx", igwx ) CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) ! - WRITE(*,*) "NICOLA", ngw, nbnd, ik, nk, nspin, igwx, is_cmplx + WRITE(*,*) "NICOLA", ispin, ngw, nbnd, ik, nk, nspin, igwx, is_cmplx CALL iotk_close_read( iuni ) RETURN ! @@ -125,7 +126,7 @@ SUBROUTINE read_wf_occ(filename, wf, ngw, nbnd) CALL iotk_scan_attr( attr, "ispin", ispin ) CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) ! - WRITE(*,*) "NICOLA", ngw_, nbnd_, ispin, igwx, is_cmplx + WRITE(*,*) "NICOLA", TRIM(filename), ngw_, nbnd_, ispin, igwx, is_cmplx ! IF (nbnd_ /= nbnd .OR. ngw_ /= ngw) THEN WRITE(*,*) "Size Mismatch. STOP" From dfcf2edc38ba434ec4d7068f056fe9c63921125f Mon Sep 17 00:00:00 2001 From: ncolonna Date: Fri, 18 Aug 2023 09:44:17 +0200 Subject: [PATCH 03/11] Add counting of electron numbers --- PP/n2npm1.f90 | 304 ++++++++++++++++++++++++++++++-------------------- 1 file changed, 186 insertions(+), 118 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index 2637fa7e..01f1fade 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -1,150 +1,218 @@ ! TO DO: ! 1) this is the desired comand with arguments from command line -! --fill tmpdir_in tmpdir_out -! --empty tmpdir_in tmpdir_out -! 2) Add reading routine for empty state (only if --empty) +! --fill tmpdir_in tmpdir_out +! --empty tmpdir_in tmpdir_out +! where is the index of the orbital of the N-electron calculation +! we want to add/remove (NB: need a convention here on how to count empty states +! and spin) +! 2) Add reading routine for empty state (only if --fill) ! 3) Check Makefile and the need of mpif90 wrapper (it loooks like we need it ! as iotk as some dependendence on mpi library. Also we might need a different ! rule to compile this program (no need for any precompilation flags like e.g. -D__MPI) +! RELEVANT observations: +! 1) for the N-1 it seems to me that no matter what the strcuture of the evc files for N-1 +! is identical to that of N. This is because the number of bands is determined by the +! majority spin channel (that does not change in structure going from N to N-1). +! It should be possible to restart a N-1 calculation from the N wfcs. Still it would +! be better to correctly initialize the last band in the minority wfc to ZERO (I guess +! this is what the code will do -it will read only N-1 out of N and set to zero the last +! one- but better to be safe) +! +!-------------------------------------------------------------------------- PROGRAM n2npm1 !------------------------------------------------------------------------ ! IMPLICIT NONE - CHARACTER(LEN=256) :: filename + CHARACTER(LEN=256) :: task, dir_in, filename_in INTEGER, PARAMETER :: DP = selected_real_kind(14,200) COMPLEX(DP), ALLOCATABLE :: evc(:,:,:) - INTEGER :: ngw, nbnd(2), ispin, nspin, igwx - LOGICAL :: is_cmplx - ! - ngw = 0; nbnd=0; nspin=0 - filename='evc01.dat' ! FIXME this will be dependent on the input path + INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd + REAL(DP) :: nel(2) + LOGICAL :: is_cmplx, l_fill + INTEGER :: index + ! + task="empty" ! FIXME this will be dependent on the input + index=1 ! FIXME this will be dependent on the input + dir_in='./' ! FIXME this will be dependent on the input path + ! + ! l_fill determine what to do: remove from an occupied one (l_fill=F) + ! or add to empty (l_fill=T) + l_fill=.false. + IF (task=="fill") l_fill=.true. + ! + IF (l_fill) THEN + WRITE(*, '(/, 3X, A, I5)') "TASK: Going to add one electron from orb ", index + ELSE + WRITE(*, '(/, 3X, A, I5)') "TASK: Going to remove one electron from orb ", index + ENDIF ! - CALL read_sizes (filename, ngw, nbnd(1), nspin, is_cmplx, igwx) + ngw = 0; nbnd=0; nspin=0; nel=0.D0 + ! READ the shape of the wfc object (# of PW, # spin, complx/real, # bands) + filename_in=TRIM(dir_in)//'evc01.dat' + CALL read_sizes (filename_in, ngw, nbnd(1), nspin, is_cmplx, igwx) ! IF (nspin==2) THEN - filename='evc02.dat' - CALL read_sizes (filename, ngw, nbnd(2), nspin, is_cmplx, igwx) + ! If needed READ the shape of the wfc object for spin down + filename_in=TRIM(dir_in)//'evc02.dat' + CALL read_sizes (filename_in, ngw, nbnd(2), nspin, is_cmplx, igwx) ENDIF - ! WRITE(*,*) igwx, ngw, nbnd(:), nspin + ! ALLOCATE (evc(igwx, max(nbnd(1), nbnd(2)), nspin)) ! - ! reset filename - filename='evc01.dat' + ! reset filename + filename_in=TRIM(dir_in)//'evc01.dat' DO ispin = 1, nspin - IF (ispin==2 ) filename="evc02.dat" - CALL read_wf_occ (filename, evc(:,:,ispin), ngw, nbnd(ispin)) - WRITE(*,*) ispin, evc(1,1,ispin), evc(igwx,nbnd(ispin),ispin) ! check debug + IF (ispin==2 ) filename_in=TRIM(dir_in)//'evc02.dat' + CALL read_wf_occ (filename_in, evc(:,:,ispin), ngw, nbnd(ispin)) + CALL check_nele (ngw, nbnd(ispin), evc(:,:,ispin), nel(ispin)) + WRITE(*,*) ispin, nel(ispin), NINT(nel(ispin)), evc(1,1,ispin), evc(igwx,nbnd(ispin),ispin)! check debug ENDDO ! - CONTAINS - ! - ! - SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) - !------------------------------------------------------------------------ - ! - USE iotk_module - ! - IMPLICIT NONE - CHARACTER(LEN=256), INTENT(IN) :: filename - INTEGER, INTENT(OUT) :: ngw, nbnd, nspin, igwx - LOGICAL, INTENT(OUT) :: is_cmplx - - INTEGER, PARAMETER :: DP = selected_real_kind(14,200) - CHARACTER(iotk_attlenx) :: attr - INTEGER :: j, ispin - INTEGER :: ierr - INTEGER :: ik, nk - INTEGER :: iuni - ! - ierr = 0 - iuni=789 - ! - CALL iotk_open_read( iuni, FILE = filename, & - BINARY = .TRUE., IERR = ierr ) - ! - IF (ierr /= 0) THEN - WRITE(*,*) "Something wrong while reading file = ", filename + ! At this stage I read the empty manifold if needed + IF (l_fill) THEN + ! Here we need to read also empty states if + WRITE(*,*) "Case l_fill = T NOT IMPLEMENTED YET: STOP" STOP ENDIF ! - CALL iotk_scan_empty( iuni, "INFO", attr ) - ! - !WRITE(*,*) "NICOLA", ngw, nbnd, ik, nk, nspin, igwx, is_cmplx - CALL iotk_scan_attr( attr, "ngw", ngw ) - CALL iotk_scan_attr( attr, "nbnd", nbnd ) - CALL iotk_scan_attr( attr, "ik", ik ) - CALL iotk_scan_attr( attr, "nk", nk ) - CALL iotk_scan_attr( attr, "nspin", nspin ) - CALL iotk_scan_attr( attr, "ispin", ispin ) - CALL iotk_scan_attr( attr, "igwx", igwx ) - CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) - ! - WRITE(*,*) "NICOLA", ispin, ngw, nbnd, ik, nk, nspin, igwx, is_cmplx - CALL iotk_close_read( iuni ) - RETURN - ! - END subroutine - ! - ! - SUBROUTINE read_wf_occ(filename, wf, ngw, nbnd) - !------------------------------------------------------------------------ - ! - USE iotk_module - ! - IMPLICIT NONE - INTEGER, PARAMETER :: DP = selected_real_kind(14,200) - CHARACTER(LEN=256), INTENT(IN) :: filename - COMPLEX(DP) :: wf(ngw, nbnd) - INTEGER, INTENT(IN) :: ngw, nbnd - CHARACTER(iotk_attlenx) :: attr - INTEGER :: j - COMPLEX(DP), ALLOCATABLE :: wtmp(:) - INTEGER :: ierr - INTEGER :: igwx, ngw_, nbnd_ - INTEGER :: iuni - LOGICAL :: is_cmplx - ! - ierr = 0 - iuni=789 - ! - CALL iotk_open_read( iuni, FILE = filename, & - BINARY = .TRUE., IERR = ierr ) - ! - IF (ierr /= 0) THEN - WRITE(*,*) "Something wrong while reading file = ", filename - STOP + ! check if the index of the orbital is within a valid range + IF (l_fill) THEN + ! + ELSE + IF (index .gt. nbnd(1)+nbnd(2)) THEN + WRITE(*,*) "index out of range ", index, nbnd(1)+nbnd(2) + STOP + ENDIF ENDIF ! - CALL iotk_scan_empty( iuni, "INFO", attr ) - ! - ngw_=0; nbnd_=0; - WRITE(*,*) "NICOLA", ngw_, nbnd_ - CALL iotk_scan_attr( attr, "ngw", ngw_ ) - CALL iotk_scan_attr( attr, "nbnd", nbnd_ ) - CALL iotk_scan_attr( attr, "igwx", igwx ) - CALL iotk_scan_attr( attr, "ispin", ispin ) - CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) - ! - WRITE(*,*) "NICOLA", TRIM(filename), ngw_, nbnd_, ispin, igwx, is_cmplx + CONTAINS ! - IF (nbnd_ /= nbnd .OR. ngw_ /= ngw) THEN - WRITE(*,*) "Size Mismatch. STOP" - STOP - ENDIF + !---------------------------------------------------------------------- + SUBROUTINE check_nele (ngw, nbnd, evc, nele) + !-------------------------------------------------------------------- + IMPLICIT NONE + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + INTEGER, INTENT(IN) :: ngw, nbnd + COMPLEX(DP), INTENT(IN) :: evc(ngw,nbnd) + REAL(DP), INTENT(OUT) :: nele + INTEGER :: ibnd + ! + nele=0.D0 + DO ibnd = 1, nbnd + nele = nele + SUM(evc(:,ibnd)*CONJG(evc(:,ibnd))) + !WRITE(*,'(I5,3x , 2F20.15, 3x, F10.6)') ibnd, SUM(evc(:,ibnd)*CONJG(evc(:,ibnd))), nele + ENDDO + RETURN + ! + END SUBROUTINE + ! + !----------------------------------------------------------------------- + SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) + !------------------------------------------------------------------------ + ! + USE iotk_module + ! + IMPLICIT NONE + CHARACTER(LEN=256), INTENT(IN) :: filename + INTEGER, INTENT(OUT) :: ngw, nbnd, nspin, igwx + LOGICAL, INTENT(OUT) :: is_cmplx - ALLOCATE( wtmp( igwx ) ) - ! - DO j = 1, nbnd - CALL iotk_scan_dat( iuni, "evc" // iotk_index( j ), wtmp(1:igwx) ) - wf(:,j) = wtmp (:) - ENDDO - ! - DEALLOCATE (wtmp) - ! - CALL iotk_close_read( iuni ) - RETURN + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + CHARACTER(iotk_attlenx) :: attr + INTEGER :: j, ispin + INTEGER :: ierr + INTEGER :: ik, nk + INTEGER :: iuni + ! + ierr = 0 + iuni=789 + ! + CALL iotk_open_read( iuni, FILE = filename, & + BINARY = .TRUE., IERR = ierr ) + ! + IF (ierr /= 0) THEN + WRITE(*,*) "Something wrong while reading file = ", filename + STOP + ENDIF + ! + CALL iotk_scan_empty( iuni, "INFO", attr ) + ! + !WRITE(*,*) "NICOLA", ngw, nbnd, ik, nk, nspin, igwx, is_cmplx + CALL iotk_scan_attr( attr, "ngw", ngw ) + CALL iotk_scan_attr( attr, "nbnd", nbnd ) + CALL iotk_scan_attr( attr, "ik", ik ) + CALL iotk_scan_attr( attr, "nk", nk ) + CALL iotk_scan_attr( attr, "nspin", nspin ) + CALL iotk_scan_attr( attr, "ispin", ispin ) + CALL iotk_scan_attr( attr, "igwx", igwx ) + CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) + ! + WRITE(*,*) "NICOLA", ispin, ngw, nbnd, ik, nk, nspin, igwx, is_cmplx + CALL iotk_close_read( iuni ) + RETURN + ! + END subroutine ! + !------------------------------------------------------------------------ + SUBROUTINE read_wf_occ(filename, wf, ngw, nbnd) + !------------------------------------------------------------------------ + ! + USE iotk_module + ! + IMPLICIT NONE + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + CHARACTER(LEN=256), INTENT(IN) :: filename + COMPLEX(DP) :: wf(ngw, nbnd) + INTEGER, INTENT(IN) :: ngw, nbnd + CHARACTER(iotk_attlenx) :: attr + INTEGER :: j + COMPLEX(DP), ALLOCATABLE :: wtmp(:) + INTEGER :: ierr + INTEGER :: igwx, ngw_, nbnd_ + INTEGER :: iuni + LOGICAL :: is_cmplx + ! + ierr = 0 + iuni=789 + ! + CALL iotk_open_read( iuni, FILE = filename, & + BINARY = .TRUE., IERR = ierr ) + ! + IF (ierr /= 0) THEN + WRITE(*,*) "Something wrong while reading file = ", filename + STOP + ENDIF + ! + CALL iotk_scan_empty( iuni, "INFO", attr ) + ! + ngw_=0; nbnd_=0; + WRITE(*,*) "NICOLA", ngw_, nbnd_ + CALL iotk_scan_attr( attr, "ngw", ngw_ ) + CALL iotk_scan_attr( attr, "nbnd", nbnd_ ) + CALL iotk_scan_attr( attr, "igwx", igwx ) + CALL iotk_scan_attr( attr, "ispin", ispin ) + CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) + ! + WRITE(*,*) "NICOLA", TRIM(filename), ngw_, nbnd_, ispin, igwx, is_cmplx + ! + IF (nbnd_ /= nbnd .OR. ngw_ /= ngw) THEN + WRITE(*,*) "Size Mismatch. STOP" + STOP + ENDIF + + ALLOCATE( wtmp( igwx ) ) + ! + DO j = 1, nbnd + CALL iotk_scan_dat( iuni, "evc" // iotk_index( j ), wtmp(1:igwx) ) + wf(:,j) = wtmp (:) + ENDDO + ! + DEALLOCATE (wtmp) + ! + CALL iotk_close_read( iuni ) + RETURN + ! END subroutine ! END PROGRAM From f5de6e0e7e870a3b48d35efac731528144b0da77 Mon Sep 17 00:00:00 2001 From: ncolonna Date: Fri, 18 Aug 2023 11:00:31 +0200 Subject: [PATCH 04/11] prepare the new wcf objects --- PP/n2npm1.f90 | 52 +++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 44 insertions(+), 8 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index 01f1fade..e1ddd342 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -1,10 +1,11 @@ ! TO DO: ! 1) this is the desired comand with arguments from command line -! --fill tmpdir_in tmpdir_out -! --empty tmpdir_in tmpdir_out +! --fill tmpdir_in tmpdir_out +! --empty tmpdir_in tmpdir_out ! where is the index of the orbital of the N-electron calculation ! we want to add/remove (NB: need a convention here on how to count empty states -! and spin) +! and spin. I think it would be convenient to give the spin channel from input as +! we have different evc files for different spin ) ! 2) Add reading routine for empty state (only if --fill) ! 3) Check Makefile and the need of mpif90 wrapper (it loooks like we need it ! as iotk as some dependendence on mpi library. Also we might need a different @@ -26,14 +27,23 @@ PROGRAM n2npm1 CHARACTER(LEN=256) :: task, dir_in, filename_in INTEGER, PARAMETER :: DP = selected_real_kind(14,200) COMPLEX(DP), ALLOCATABLE :: evc(:,:,:) - INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd + COMPLEX(DP), ALLOCATABLE :: evc_empty(:,:,:) + COMPLEX(DP), ALLOCATABLE :: evc_out(:,:,:) + INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd, ibnd_, spin_channel + INTEGER :: nbnd_out(2) REAL(DP) :: nel(2) LOGICAL :: is_cmplx, l_fill INTEGER :: index ! task="empty" ! FIXME this will be dependent on the input index=1 ! FIXME this will be dependent on the input - dir_in='./' ! FIXME this will be dependent on the input path + dir_in='./' ! FIXME this will be dependent on the input path + spin_channel = 1 ! FIXME: this will track on which spin channel we need to add/remove the electron + ! ! Requires a way determine which spin channel from input (either explicite input or + ! ! inferred from index which should go from 1 to nbnd(1)+nbnd(2) + ! + CALL get_command_argument( 1, arg ) + CALL parse_args( io_files, nrtot, output_exst ) ! ! l_fill determine what to do: remove from an occupied one (l_fill=F) ! or add to empty (l_fill=T) @@ -41,7 +51,7 @@ PROGRAM n2npm1 IF (task=="fill") l_fill=.true. ! IF (l_fill) THEN - WRITE(*, '(/, 3X, A, I5)') "TASK: Going to add one electron from orb ", index + WRITE(*, '(/, 3X, A, I5)') "TASK: Going to add one electron taken from orb ", index ELSE WRITE(*, '(/, 3X, A, I5)') "TASK: Going to remove one electron from orb ", index ENDIF @@ -86,6 +96,33 @@ PROGRAM n2npm1 ENDIF ENDIF ! + ! From now on we reshape the wfc object + nbnd_out = nbnd + IF (l_fill) THEN + nbnd_out(spin_channel) = nbnd(spin_channel) + 1 + ALLOCATE (evc_out(igwx, max(nbnd_out(1), nbnd_out(2)), nspin)) + evc_out = CMPLX(0.D0, 0.D0, kind =DP) + DO ispin = 1, nspin + DO ibnd = 1,nbnd(ispin) + evc_out(:,ibnd,ispin) = evc(:,ibnd,ispin) + ENDDO + ENDDO + evc_out(:,nbnd_out(spin_channel),spin_channel) = evc_empty(:,index,spin_channel) ! + ! + ELSE + nbnd_out(spin_channel) = nbnd(spin_channel) - 1 + ALLOCATE (evc_out(igwx, max(nbnd_out(1), nbnd_out(2)), nspin)) + evc_out = CMPLX(0.D0, 0.D0, kind =DP) + DO ispin = 1, nspin + ibnd_=0 + DO ibnd = 1,nbnd(ispin) + IF (ispin == spin_channel .AND. ibnd == index) CYCLE + ibnd_=ibnd_+1 + evc_out(:,ibnd_,ispin) = evc(:,ibnd,ispin) + ENDDO + WRITE(*,*) ispin, evc_out(1,1,ispin), evc_out(igwx,MAX(nbnd_out(1),nbnd_out(2)),ispin)! check debug + ENDDO + ENDIF CONTAINS ! !---------------------------------------------------------------------- @@ -128,8 +165,7 @@ SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) ierr = 0 iuni=789 ! - CALL iotk_open_read( iuni, FILE = filename, & - BINARY = .TRUE., IERR = ierr ) + CALL iotk_open_read( iuni, FILE = filename, BINARY = .TRUE., IERR = ierr ) ! IF (ierr /= 0) THEN WRITE(*,*) "Something wrong while reading file = ", filename From eb5d62e8bd55cfb15fc709ec15c9b371fa4160c2 Mon Sep 17 00:00:00 2001 From: ncolonna Date: Fri, 18 Aug 2023 11:01:26 +0200 Subject: [PATCH 05/11] OOPS --- PP/n2npm1.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index e1ddd342..f568befc 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -42,8 +42,8 @@ PROGRAM n2npm1 ! ! Requires a way determine which spin channel from input (either explicite input or ! ! inferred from index which should go from 1 to nbnd(1)+nbnd(2) ! - CALL get_command_argument( 1, arg ) - CALL parse_args( io_files, nrtot, output_exst ) + !CALL get_command_argument( 1, arg ) + !CALL parse_args( io_files, nrtot, output_exst ) ! ! l_fill determine what to do: remove from an occupied one (l_fill=F) ! or add to empty (l_fill=T) From c616650fe0952545df5a622f0afecf5aa30b127c Mon Sep 17 00:00:00 2001 From: ncolonna Date: Fri, 18 Aug 2023 13:05:08 +0200 Subject: [PATCH 06/11] Print WFC for N-1 case --- PP/n2npm1.f90 | 79 +++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 70 insertions(+), 9 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index f568befc..3e70bd28 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -24,7 +24,7 @@ PROGRAM n2npm1 !------------------------------------------------------------------------ ! IMPLICIT NONE - CHARACTER(LEN=256) :: task, dir_in, filename_in + CHARACTER(LEN=256) :: task, dir_in, filename_in, dir_out, filename_out INTEGER, PARAMETER :: DP = selected_real_kind(14,200) COMPLEX(DP), ALLOCATABLE :: evc(:,:,:) COMPLEX(DP), ALLOCATABLE :: evc_empty(:,:,:) @@ -32,12 +32,14 @@ PROGRAM n2npm1 INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd, ibnd_, spin_channel INTEGER :: nbnd_out(2) REAL(DP) :: nel(2) - LOGICAL :: is_cmplx, l_fill + LOGICAL :: is_cmplx, l_fill, gamma_only INTEGER :: index + REAL(DP) :: scalef ! task="empty" ! FIXME this will be dependent on the input index=1 ! FIXME this will be dependent on the input dir_in='./' ! FIXME this will be dependent on the input path + dir_out='./' ! FIXME this will be dependent on the input path spin_channel = 1 ! FIXME: this will track on which spin channel we need to add/remove the electron ! ! Requires a way determine which spin channel from input (either explicite input or ! ! inferred from index which should go from 1 to nbnd(1)+nbnd(2) @@ -51,20 +53,22 @@ PROGRAM n2npm1 IF (task=="fill") l_fill=.true. ! IF (l_fill) THEN - WRITE(*, '(/, 3X, A, I5)') "TASK: Going to add one electron taken from orb ", index + WRITE(*, '(/, 3X, 2(A, I5))') "TASK: Going to add & + one electron taken from empty orbital= ", index, " spin= ", spin_channel ELSE - WRITE(*, '(/, 3X, A, I5)') "TASK: Going to remove one electron from orb ", index + WRITE(*, '(/, 3X, 2(A, I5))') "TASK: Going to remove & + one electron from occupied orbital= ", index, " spin= ", spin_channel ENDIF ! ngw = 0; nbnd=0; nspin=0; nel=0.D0 ! READ the shape of the wfc object (# of PW, # spin, complx/real, # bands) filename_in=TRIM(dir_in)//'evc01.dat' - CALL read_sizes (filename_in, ngw, nbnd(1), nspin, is_cmplx, igwx) + CALL read_sizes (filename_in, ngw, nbnd(1), nspin, is_cmplx, igwx, scalef, gamma_only) ! IF (nspin==2) THEN ! If needed READ the shape of the wfc object for spin down filename_in=TRIM(dir_in)//'evc02.dat' - CALL read_sizes (filename_in, ngw, nbnd(2), nspin, is_cmplx, igwx) + CALL read_sizes (filename_in, ngw, nbnd(2), nspin, is_cmplx, igwx, scalef, gamma_only) ENDIF WRITE(*,*) igwx, ngw, nbnd(:), nspin ! @@ -123,8 +127,62 @@ PROGRAM n2npm1 WRITE(*,*) ispin, evc_out(1,1,ispin), evc_out(igwx,MAX(nbnd_out(1),nbnd_out(2)),ispin)! check debug ENDDO ENDIF + ! + filename_out=TRIM(dir_out)//'evc01_.dat' + DO ispin = 1, nspin + IF (ispin==2 ) filename_out=TRIM(dir_out)//'evc02_.dat' + CALL write_wf (filename_out, evc_out(:,:,ispin), nspin, ispin, ngw, igwx, & + MAX(nbnd(1), nbnd(2)), is_cmplx, scalef, gamma_only) + ENDDO CONTAINS ! + ! + !---------------------------------------------------------------------- + SUBROUTINE write_wf(filename, wf, nspin, ispin, ngw, igwx, nbnd, is_cmplx, scalef, gamma_only) + !-------------------------------------------------------------------- + ! + USE iotk_module + ! + IMPLICIT NONE + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + CHARACTER(LEN=256), INTENT(IN) :: filename + COMPLEX(DP) :: wf(ngw, nbnd) + INTEGER, INTENT(IN) :: ngw, nbnd, nspin, igwx, ispin + LOGICAL, INTENT(IN) :: is_cmplx, gamma_only + REAL(DP), INTENT(IN) :: scalef + CHARACTER(iotk_attlenx) :: attr + INTEGER :: j + COMPLEX(DP), ALLOCATABLE :: wtmp(:) + INTEGER :: iuni + ! + iuni=789 + ! + CALL iotk_open_write( iuni, FILE = TRIM( filename ), ROOT="WFC", BINARY = .TRUE. ) + ! + CALL iotk_write_attr( attr, "ngw", ngw, FIRST = .TRUE. ) + CALL iotk_write_attr( attr, "igwx", igwx ) + CALL iotk_write_attr( attr, "do_wf_cmplx", is_cmplx ) !added:giovanni + CALL iotk_write_attr( attr, "gamma_only", gamma_only.and..not.is_cmplx ) + CALL iotk_write_attr( attr, "nbnd", nbnd ) + CALL iotk_write_attr( attr, "ik", ispin ) + CALL iotk_write_attr( attr, "nk", nspin ) + CALL iotk_write_attr( attr, "ispin", ispin ) + CALL iotk_write_attr( attr, "nspin", nspin ) + CALL iotk_write_attr( attr, "scale_factor", scalef ) + ! + CALL iotk_write_empty( iuni, "INFO", attr ) + ! + ALLOCATE( wtmp( MAX( igwx, 1 ) ) ) + wtmp = 0.0_DP + DO j = 1, nbnd + wtmp(:)=wf(:,j) + CALL iotk_write_dat( iuni, "evc" // iotk_index( j ), wtmp(1:igwx) ) + ENDDO + ! + CALL iotk_close_write( iuni ) + ! + END SUBROUTINE + ! !---------------------------------------------------------------------- SUBROUTINE check_nele (ngw, nbnd, evc, nele) !-------------------------------------------------------------------- @@ -145,7 +203,7 @@ SUBROUTINE check_nele (ngw, nbnd, evc, nele) END SUBROUTINE ! !----------------------------------------------------------------------- - SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) + SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx, scalef, gamma_only) !------------------------------------------------------------------------ ! USE iotk_module @@ -153,7 +211,7 @@ SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) IMPLICIT NONE CHARACTER(LEN=256), INTENT(IN) :: filename INTEGER, INTENT(OUT) :: ngw, nbnd, nspin, igwx - LOGICAL, INTENT(OUT) :: is_cmplx + LOGICAL, INTENT(OUT) :: is_cmplx, gamma_only INTEGER, PARAMETER :: DP = selected_real_kind(14,200) CHARACTER(iotk_attlenx) :: attr @@ -161,6 +219,7 @@ SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) INTEGER :: ierr INTEGER :: ik, nk INTEGER :: iuni + REAL(DP) :: scalef ! ierr = 0 iuni=789 @@ -183,8 +242,10 @@ SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx) CALL iotk_scan_attr( attr, "ispin", ispin ) CALL iotk_scan_attr( attr, "igwx", igwx ) CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) + CALL iotk_scan_attr( attr, "gamma_only", gamma_only ) + CALL iotk_scan_attr( attr, "scale_factor", scalef ) ! - WRITE(*,*) "NICOLA", ispin, ngw, nbnd, ik, nk, nspin, igwx, is_cmplx + WRITE(*,*) "NICOLA", ispin, ngw, nbnd, ik, nk, nspin, igwx, is_cmplx, scalef CALL iotk_close_read( iuni ) RETURN ! From 9fd6e32441c8500d610c43e1dc03307dcdc30bf3 Mon Sep 17 00:00:00 2001 From: ncolonna Date: Fri, 18 Aug 2023 15:54:44 +0200 Subject: [PATCH 07/11] Basic implementation done. --- PP/n2npm1.f90 | 175 +++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 153 insertions(+), 22 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index 3e70bd28..ab321484 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -1,4 +1,4 @@ -! TO DO: +! DONE: ! 1) this is the desired comand with arguments from command line ! --fill tmpdir_in tmpdir_out ! --empty tmpdir_in tmpdir_out @@ -7,6 +7,7 @@ ! and spin. I think it would be convenient to give the spin channel from input as ! we have different evc files for different spin ) ! 2) Add reading routine for empty state (only if --fill) +! TODO: ! 3) Check Makefile and the need of mpif90 wrapper (it loooks like we need it ! as iotk as some dependendence on mpi library. Also we might need a different ! rule to compile this program (no need for any precompilation flags like e.g. -D__MPI) @@ -29,14 +30,17 @@ PROGRAM n2npm1 COMPLEX(DP), ALLOCATABLE :: evc(:,:,:) COMPLEX(DP), ALLOCATABLE :: evc_empty(:,:,:) COMPLEX(DP), ALLOCATABLE :: evc_out(:,:,:) - INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd, ibnd_, spin_channel + INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd, ibnd_, spin_channel, ig + INTEGER :: nbnd_emp(2) INTEGER :: nbnd_out(2) REAL(DP) :: nel(2) - LOGICAL :: is_cmplx, l_fill, gamma_only + LOGICAL :: is_cmplx, l_fill, gamma_only, exst INTEGER :: index REAL(DP) :: scalef + INTEGER :: iuni ! - task="empty" ! FIXME this will be dependent on the input + !Defauls + task="fill" ! FIXME this will be dependent on the input index=1 ! FIXME this will be dependent on the input dir_in='./' ! FIXME this will be dependent on the input path dir_out='./' ! FIXME this will be dependent on the input path @@ -44,8 +48,7 @@ PROGRAM n2npm1 ! ! Requires a way determine which spin channel from input (either explicite input or ! ! inferred from index which should go from 1 to nbnd(1)+nbnd(2) ! - !CALL get_command_argument( 1, arg ) - !CALL parse_args( io_files, nrtot, output_exst ) + CALL read_command_line (task, index, spin_channel, dir_in, dir_out) ! ! l_fill determine what to do: remove from an occupied one (l_fill=F) ! or add to empty (l_fill=T) @@ -53,13 +56,23 @@ PROGRAM n2npm1 IF (task=="fill") l_fill=.true. ! IF (l_fill) THEN - WRITE(*, '(/, 3X, 2(A, I5))') "TASK: Going to add & + WRITE(*, '(/, 3X, 2(A, I5))') "TASK: N --> N+1; Adding & one electron taken from empty orbital= ", index, " spin= ", spin_channel ELSE - WRITE(*, '(/, 3X, 2(A, I5))') "TASK: Going to remove & + WRITE(*, '(/, 3X, 2(A, I5))') "TASK: N --> N-1; Removing & one electron from occupied orbital= ", index, " spin= ", spin_channel ENDIF + WRITE(*,'(/, 5X, "SUMMARY: task =", A)') TRIM(task) + WRITE(*,'( 5X, " index =", I5)') index + WRITE(*,'( 5X, " spin =", I5)') spin_channel + WRITE(*,'( 5X, " dir_in =", A)') TRIM(dir_in) + WRITE(*,'( 5X, " dir_out =", A)') TRIM(dir_out) ! + IF (dir_in == dir_out) THEN + WRITE(*,'(3X, 3A)') "ERROR: dir_out = dir_in" + STOP + ENDIF + WRITE(*,'(/ 5X, A)') "Reading OCC manifold ..." ngw = 0; nbnd=0; nspin=0; nel=0.D0 ! READ the shape of the wfc object (# of PW, # spin, complx/real, # bands) filename_in=TRIM(dir_in)//'evc01.dat' @@ -70,7 +83,14 @@ PROGRAM n2npm1 filename_in=TRIM(dir_in)//'evc02.dat' CALL read_sizes (filename_in, ngw, nbnd(2), nspin, is_cmplx, igwx, scalef, gamma_only) ENDIF - WRITE(*,*) igwx, ngw, nbnd(:), nspin + ! + WRITE(*,'(/, 7X, "SUMMARY: igwx =", I5)') igwx + WRITE(*,'( 7X, " ngw =", I5)') ngw + WRITE(*,'( 7X, " nspin =", I5)') nspin + WRITE(*,'( 7X, " nbnd(1) =", I5)') nbnd(1) + WRITE(*,'( 7X, " nbnd(2) =", I5)') nbnd(2) + WRITE(*,'( 7X, " is_cmplx =", L5)') is_cmplx + WRITE(*,'( 7X, " G_trick =", L5, /)') gamma_only ! ALLOCATE (evc(igwx, max(nbnd(1), nbnd(2)), nspin)) ! @@ -80,27 +100,78 @@ PROGRAM n2npm1 IF (ispin==2 ) filename_in=TRIM(dir_in)//'evc02.dat' CALL read_wf_occ (filename_in, evc(:,:,ispin), ngw, nbnd(ispin)) CALL check_nele (ngw, nbnd(ispin), evc(:,:,ispin), nel(ispin)) - WRITE(*,*) ispin, nel(ispin), NINT(nel(ispin)), evc(1,1,ispin), evc(igwx,nbnd(ispin),ispin)! check debug + WRITE(*,'(7X, "CHECK: ispin =", I5, " nel =", F18.12, " INT(nel)=", I5)') & + ispin, nel(ispin), NINT(nel(ispin)) + !WRITE(*,) ispin, nel(ispin), NINT(nel(ispin)) ENDDO ! ! At this stage I read the empty manifold if needed IF (l_fill) THEN - ! Here we need to read also empty states if - WRITE(*,*) "Case l_fill = T NOT IMPLEMENTED YET: STOP" - STOP + WRITE(*,'(/ 5X, A)') "Reading EMP manifold ..." + ! Here we need to read also empty states if needed. + ! First the shape of the arrey ngw, nbnd_emp for spin up + iuni=789 + filename_in=TRIM(dir_in)//'evc0_empty1.dat' + INQUIRE (FILE=TRIM(filename_in), EXIST=EXST) + IF (.NOT. exst) THEN + WRITE(*,'(3A)') "File ", TRIM(filename_in), " do not exist" + STOP + ENDIF + OPEN (UNIT=iuni, FILE=TRIM(filename_in), status='unknown', FORM='UNFORMATTED') + READ (iuni) ngw, nbnd_emp(1) + ! then for spin down channel + IF (nspin == 2) THEN + iuni=iuni+1 + filename_in=TRIM(dir_in)//'evc0_empty2.dat' + INQUIRE (FILE=TRIM(filename_in), EXIST=EXST) + IF (.NOT. exst) THEN + WRITE(*,'(3A)') "File ", TRIM(filename_in), " do not exist" + STOP + ENDIF + OPEN (UNIT=iuni, FILE=TRIM(filename_in), status='unknown', FORM='UNFORMATTED') + READ (iuni) ngw, nbnd_emp(2) + ENDIF + WRITE(*,'(/, 7X, "SUMMARY: ngw =", I5)') ngw + WRITE(*,'( 7X, " nbnd(1) =", I5)') nbnd_emp(1) + WRITE(*,'( 7X, " nbnd(2) =", I5)') nbnd_emp(2) + ! allocation + ALLOCATE ( evc_empty (ngw, MAX(nbnd_emp(1), nbnd_emp(2)),nspin) ) + ! + ! read the empty-state wfcs + iuni=789 + filename_in=TRIM(dir_in)//'evc0_empty1.dat' + DO ispin = 1, nspin + IF (ispin == 2) filename_in=TRIM(dir_in)//'evc0_empty2.dat' + IF (ispin == 2) iuni=iuni+1 + DO ibnd = 1, nbnd_emp(ispin) + READ (iuni) (evc_empty(ig,ibnd,ispin), ig=1, ngw) + ENDDO + ENDDO ENDIF ! ! check if the index of the orbital is within a valid range IF (l_fill) THEN ! + IF (spin_channel == 1 .AND. index .gt. nbnd_emp(1)) THEN + WRITE(*,*) "index out of range ", index, nbnd_emp(1); STOP + ENDIF + IF (spin_channel == 2 .AND. index .gt. nbnd_emp(2)) THEN + WRITE(*,*) "index out of range ", index, nbnd_emp(2); STOP + ENDIF + ! ELSE - IF (index .gt. nbnd(1)+nbnd(2)) THEN - WRITE(*,*) "index out of range ", index, nbnd(1)+nbnd(2) - STOP + ! + IF (spin_channel == 1 .AND. index .gt. nbnd(1)) THEN + WRITE(*,*) "index out of range ", index, nbnd(1); STOP + ENDIF + IF (spin_channel == 2 .AND. index .gt. nbnd(2)) THEN + WRITE(*,*) "index out of range ", index, nbnd(2); STOP ENDIF + ! ENDIF ! ! From now on we reshape the wfc object + WRITE(*,'(/ 5X, A)') "Re-shaping OCC manifold ..." nbnd_out = nbnd IF (l_fill) THEN nbnd_out(spin_channel) = nbnd(spin_channel) + 1 @@ -128,14 +199,74 @@ PROGRAM n2npm1 ENDDO ENDIF ! - filename_out=TRIM(dir_out)//'evc01_.dat' + WRITE(*,'(/, 7X, "SUMMARY: igwx =", I5)') igwx + WRITE(*,'( 7X, " ngw =", I5)') ngw + WRITE(*,'( 7X, " nspin =", I5)') nspin + WRITE(*,'( 7X, " nbnd(1) =", I5)') MAX(nbnd_out(1),nbnd_out(2)) + WRITE(*,'( 7X, " nbnd(2) =", I5)') MAX(nbnd_out(1),nbnd_out(2)) + WRITE(*,'( 7X, " is_cmplx =", L5)') is_cmplx + WRITE(*,'( 7X, " G_trick =", L5, /)') gamma_only + ! + WRITE(*,'(/ 5X, A)') "Writing new WFCs to disk ..." + ! + filename_out=TRIM(dir_out)//'evc01.dat' DO ispin = 1, nspin - IF (ispin==2 ) filename_out=TRIM(dir_out)//'evc02_.dat' + IF (ispin==2 ) filename_out=TRIM(dir_out)//'evc02.dat' CALL write_wf (filename_out, evc_out(:,:,ispin), nspin, ispin, ngw, igwx, & - MAX(nbnd(1), nbnd(2)), is_cmplx, scalef, gamma_only) + MAX(nbnd_out(1), nbnd_out(2)), is_cmplx, scalef, gamma_only) ENDDO + ! + WRITE(*,'(/ 3X, A, /)') "GAME OVER" CONTAINS ! + !---------------------------------------------------------------------- + SUBROUTINE read_command_line (task, index, spin_channel, dir_in, dir_out) + !------------------------------------------------------------------ + ! + IMPLICIT NONE + INTEGER :: index, spin_channel + CHARACTER(LEN=256) :: task, dir_in, dir_out + INTEGER nargs, iarg, counter + CHARACTER(LEN=256) :: arg + ! + nargs = command_argument_count() + !WRITE(*,*) nargs + ! + IF ( nargs /= 5 ) THEN + WRITE(*,'(/, 5X, A, I5)') "Wrong number of arguments: STOP", nargs + WRITE(*,'( 5X, A, /)') "Usage: n2npm1.x -- " + STOP + ENDIF + ! + iarg = 1 + counter = 0 + ! + CALL get_command_argument( iarg, arg ) + iarg = iarg + 1 + ! + SELECT CASE ( trim(arg) ) + CASE ( '--fill' ) + task="fill" + CASE ( '--empty') + task="empty" + CASE DEFAULT + WRITE(*,'(A)') 'unrecognised argument option' + STOP + END SELECT + ! + CALL get_command_argument( iarg, arg ) + READ(arg,*) index + iarg = iarg + 1 + CALL get_command_argument( iarg, arg ) + READ(arg,*) spin_channel + iarg = iarg + 1 + CALL get_command_argument( iarg, dir_in ) + iarg = iarg + 1 + CALL get_command_argument( iarg, dir_out ) + ! + RETURN + ! + END SUBROUTINE ! !---------------------------------------------------------------------- SUBROUTINE write_wf(filename, wf, nspin, ispin, ngw, igwx, nbnd, is_cmplx, scalef, gamma_only) @@ -245,7 +376,7 @@ SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx, scalef, gamma_ CALL iotk_scan_attr( attr, "gamma_only", gamma_only ) CALL iotk_scan_attr( attr, "scale_factor", scalef ) ! - WRITE(*,*) "NICOLA", ispin, ngw, nbnd, ik, nk, nspin, igwx, is_cmplx, scalef + !WRITE(*,*) "NICOLA", ispin, ngw, nbnd, ik, nk, nspin, igwx, is_cmplx, scalef CALL iotk_close_read( iuni ) RETURN ! @@ -284,14 +415,14 @@ SUBROUTINE read_wf_occ(filename, wf, ngw, nbnd) CALL iotk_scan_empty( iuni, "INFO", attr ) ! ngw_=0; nbnd_=0; - WRITE(*,*) "NICOLA", ngw_, nbnd_ + !WRITE(*,*) "NICOLA", ngw_, nbnd_ CALL iotk_scan_attr( attr, "ngw", ngw_ ) CALL iotk_scan_attr( attr, "nbnd", nbnd_ ) CALL iotk_scan_attr( attr, "igwx", igwx ) CALL iotk_scan_attr( attr, "ispin", ispin ) CALL iotk_scan_attr( attr, "do_wf_cmplx", is_cmplx ) ! - WRITE(*,*) "NICOLA", TRIM(filename), ngw_, nbnd_, ispin, igwx, is_cmplx + !WRITE(*,*) "NICOLA", TRIM(filename), ngw_, nbnd_, ispin, igwx, is_cmplx ! IF (nbnd_ /= nbnd .OR. ngw_ /= ngw) THEN WRITE(*,*) "Size Mismatch. STOP" From a10d82403982253e9e533267009b25da6c7fafb3 Mon Sep 17 00:00:00 2001 From: ncolonna Date: Fri, 18 Aug 2023 22:26:45 +0200 Subject: [PATCH 08/11] Minor changes in formatting --- PP/n2npm1.f90 | 36 +++++++++++++++++++----------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index ab321484..608b03df 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -84,13 +84,13 @@ PROGRAM n2npm1 CALL read_sizes (filename_in, ngw, nbnd(2), nspin, is_cmplx, igwx, scalef, gamma_only) ENDIF ! - WRITE(*,'(/, 7X, "SUMMARY: igwx =", I5)') igwx - WRITE(*,'( 7X, " ngw =", I5)') ngw - WRITE(*,'( 7X, " nspin =", I5)') nspin - WRITE(*,'( 7X, " nbnd(1) =", I5)') nbnd(1) - WRITE(*,'( 7X, " nbnd(2) =", I5)') nbnd(2) - WRITE(*,'( 7X, " is_cmplx =", L5)') is_cmplx - WRITE(*,'( 7X, " G_trick =", L5, /)') gamma_only + WRITE(*,'(/, 7X, "SUMMARY: igwx =", I8)') igwx + WRITE(*,'( 7X, " ngw =", I8)') ngw + WRITE(*,'( 7X, " nspin =", I8)') nspin + WRITE(*,'( 7X, " nbnd(1) =", I8)') nbnd(1) + WRITE(*,'( 7X, " nbnd(2) =", I8)') nbnd(2) + WRITE(*,'( 7X, " is_cmplx =", L8)') is_cmplx + WRITE(*,'( 7X, " G_trick =", L8, /)') gamma_only ! ALLOCATE (evc(igwx, max(nbnd(1), nbnd(2)), nspin)) ! @@ -131,9 +131,9 @@ PROGRAM n2npm1 OPEN (UNIT=iuni, FILE=TRIM(filename_in), status='unknown', FORM='UNFORMATTED') READ (iuni) ngw, nbnd_emp(2) ENDIF - WRITE(*,'(/, 7X, "SUMMARY: ngw =", I5)') ngw - WRITE(*,'( 7X, " nbnd(1) =", I5)') nbnd_emp(1) - WRITE(*,'( 7X, " nbnd(2) =", I5)') nbnd_emp(2) + WRITE(*,'(/, 7X, "SUMMARY: ngw =", I8)') ngw + WRITE(*,'( 7X, " nbnd(1) =", I8)') nbnd_emp(1) + WRITE(*,'( 7X, " nbnd(2) =", I8)') nbnd_emp(2) ! allocation ALLOCATE ( evc_empty (ngw, MAX(nbnd_emp(1), nbnd_emp(2)),nspin) ) ! @@ -199,13 +199,13 @@ PROGRAM n2npm1 ENDDO ENDIF ! - WRITE(*,'(/, 7X, "SUMMARY: igwx =", I5)') igwx - WRITE(*,'( 7X, " ngw =", I5)') ngw - WRITE(*,'( 7X, " nspin =", I5)') nspin - WRITE(*,'( 7X, " nbnd(1) =", I5)') MAX(nbnd_out(1),nbnd_out(2)) - WRITE(*,'( 7X, " nbnd(2) =", I5)') MAX(nbnd_out(1),nbnd_out(2)) - WRITE(*,'( 7X, " is_cmplx =", L5)') is_cmplx - WRITE(*,'( 7X, " G_trick =", L5, /)') gamma_only + WRITE(*,'(/, 7X, "SUMMARY: igwx =", I8)') igwx + WRITE(*,'( 7X, " ngw =", I8)') ngw + WRITE(*,'( 7X, " nspin =", I8)') nspin + WRITE(*,'( 7X, " nbnd(1) =", I8)') MAX(nbnd_out(1),nbnd_out(2)) + WRITE(*,'( 7X, " nbnd(2) =", I8)') MAX(nbnd_out(1),nbnd_out(2)) + WRITE(*,'( 7X, " is_cmplx =", L8)') is_cmplx + WRITE(*,'( 7X, " G_trick =", L8, /)') gamma_only ! WRITE(*,'(/ 5X, A)') "Writing new WFCs to disk ..." ! @@ -222,6 +222,8 @@ PROGRAM n2npm1 !---------------------------------------------------------------------- SUBROUTINE read_command_line (task, index, spin_channel, dir_in, dir_out) !------------------------------------------------------------------ + ! arg# 1 2 3 4 5 + ! n2npm1.x -- ! IMPLICIT NONE INTEGER :: index, spin_channel From b3ecf6988b561363f8a2facc372060bcef23a0a6 Mon Sep 17 00:00:00 2001 From: ncolonna Date: Fri, 18 Aug 2023 22:59:05 +0200 Subject: [PATCH 09/11] more --- PP/n2npm1.f90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index 608b03df..64fd4d0e 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -99,6 +99,8 @@ PROGRAM n2npm1 DO ispin = 1, nspin IF (ispin==2 ) filename_in=TRIM(dir_in)//'evc02.dat' CALL read_wf_occ (filename_in, evc(:,:,ispin), ngw, nbnd(ispin)) + ! TODO: if KS state are used at initialization, also empty states are written to evc + ! then this countig does not make sense CALL check_nele (ngw, nbnd(ispin), evc(:,:,ispin), nel(ispin)) WRITE(*,'(7X, "CHECK: ispin =", I5, " nel =", F18.12, " INT(nel)=", I5)') & ispin, nel(ispin), NINT(nel(ispin)) @@ -175,6 +177,10 @@ PROGRAM n2npm1 nbnd_out = nbnd IF (l_fill) THEN nbnd_out(spin_channel) = nbnd(spin_channel) + 1 + nel(spin_channel) = nel(spin_channel)+1 + IF (nel(1) .lt. nel(2)) THEN + WRITE(*,*) "nel up must be greater than or equal to nel dw ", NINT(nel(:)); STOP + ENDIF ALLOCATE (evc_out(igwx, max(nbnd_out(1), nbnd_out(2)), nspin)) evc_out = CMPLX(0.D0, 0.D0, kind =DP) DO ispin = 1, nspin @@ -186,6 +192,10 @@ PROGRAM n2npm1 ! ELSE nbnd_out(spin_channel) = nbnd(spin_channel) - 1 + nel(spin_channel) = nel(spin_channel)-1 + IF (nel(1) .lt. nel(2)) THEN + WRITE(*,*) "nel up must be greater than or equal to nel dw ", NINT(nel(:)); STOP + ENDIF ALLOCATE (evc_out(igwx, max(nbnd_out(1), nbnd_out(2)), nspin)) evc_out = CMPLX(0.D0, 0.D0, kind =DP) DO ispin = 1, nspin From 82181a354ead48ab3667dedcbab1b504374ea622 Mon Sep 17 00:00:00 2001 From: ncolonna Date: Mon, 21 Aug 2023 09:06:36 +0200 Subject: [PATCH 10/11] minor changes description --- PP/n2npm1.f90 | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index 64fd4d0e..8bdd6c0d 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -1,7 +1,7 @@ ! DONE: ! 1) this is the desired comand with arguments from command line -! --fill tmpdir_in tmpdir_out -! --empty tmpdir_in tmpdir_out +! --fill tmpdir_in tmpdir_out +! --empty tmpdir_in tmpdir_out ! where is the index of the orbital of the N-electron calculation ! we want to add/remove (NB: need a convention here on how to count empty states ! and spin. I think it would be convenient to give the spin channel from input as @@ -11,14 +11,10 @@ ! 3) Check Makefile and the need of mpif90 wrapper (it loooks like we need it ! as iotk as some dependendence on mpi library. Also we might need a different ! rule to compile this program (no need for any precompilation flags like e.g. -D__MPI) -! RELEVANT observations: -! 1) for the N-1 it seems to me that no matter what the strcuture of the evc files for N-1 -! is identical to that of N. This is because the number of bands is determined by the -! majority spin channel (that does not change in structure going from N to N-1). -! It should be possible to restart a N-1 calculation from the N wfcs. Still it would -! be better to correctly initialize the last band in the minority wfc to ZERO (I guess -! this is what the code will do -it will read only N-1 out of N and set to zero the last -! one- but better to be safe) +! Alternatively, keep mpi and use ionode only to perform operations +! 4) Add the possibility to initialize the extra electron/hole wfc with a completely delocalzed +! state (vector full of 1/npw in G-space). This is to prepare the starting point for the +! calculation of the reference state in a polaron simulation ! !-------------------------------------------------------------------------- PROGRAM n2npm1 @@ -48,6 +44,7 @@ PROGRAM n2npm1 ! ! Requires a way determine which spin channel from input (either explicite input or ! ! inferred from index which should go from 1 to nbnd(1)+nbnd(2) ! + ! Here read command line CALL read_command_line (task, index, spin_channel, dir_in, dir_out) ! ! l_fill determine what to do: remove from an occupied one (l_fill=F) From 6239cce0aeb9fa2770242a1f10e54210b581542f Mon Sep 17 00:00:00 2001 From: ncolonna Date: Tue, 22 Aug 2023 09:17:24 +0200 Subject: [PATCH 11/11] Modularization of the main program --- PP/n2npm1.f90 | 170 +++++++++++++++++++++++++++++++++----------------- 1 file changed, 113 insertions(+), 57 deletions(-) diff --git a/PP/n2npm1.f90 b/PP/n2npm1.f90 index 8bdd6c0d..991de003 100644 --- a/PP/n2npm1.f90 +++ b/PP/n2npm1.f90 @@ -26,13 +26,14 @@ PROGRAM n2npm1 COMPLEX(DP), ALLOCATABLE :: evc(:,:,:) COMPLEX(DP), ALLOCATABLE :: evc_empty(:,:,:) COMPLEX(DP), ALLOCATABLE :: evc_out(:,:,:) - INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd, ibnd_, spin_channel, ig + INTEGER :: ngw, nbnd(2), ispin, nspin, igwx, ibnd, spin_channel, ig + INTEGER :: ngw_, nspin_, igwx_, ibnd_ INTEGER :: nbnd_emp(2) INTEGER :: nbnd_out(2) REAL(DP) :: nel(2) - LOGICAL :: is_cmplx, l_fill, gamma_only, exst + LOGICAL :: is_cmplx, is_cmplx_, l_fill, gamma_only, gamma_only_, exst INTEGER :: index - REAL(DP) :: scalef + REAL(DP) :: scalef, scalef_ INTEGER :: iuni ! !Defauls @@ -54,31 +55,32 @@ PROGRAM n2npm1 ! IF (l_fill) THEN WRITE(*, '(/, 3X, 2(A, I5))') "TASK: N --> N+1; Adding & - one electron taken from empty orbital= ", index, " spin= ", spin_channel + & one electron taken from empty orbital= ", index, " spin= ", spin_channel ELSE WRITE(*, '(/, 3X, 2(A, I5))') "TASK: N --> N-1; Removing & - one electron from occupied orbital= ", index, " spin= ", spin_channel + & one electron from occupied orbital= ", index, " spin= ", spin_channel ENDIF - WRITE(*,'(/, 5X, "SUMMARY: task =", A)') TRIM(task) - WRITE(*,'( 5X, " index =", I5)') index - WRITE(*,'( 5X, " spin =", I5)') spin_channel - WRITE(*,'( 5X, " dir_in =", A)') TRIM(dir_in) - WRITE(*,'( 5X, " dir_out =", A)') TRIM(dir_out) ! - IF (dir_in == dir_out) THEN - WRITE(*,'(3X, 3A)') "ERROR: dir_out = dir_in" - STOP - ENDIF WRITE(*,'(/ 5X, A)') "Reading OCC manifold ..." + ! ngw = 0; nbnd=0; nspin=0; nel=0.D0 ! READ the shape of the wfc object (# of PW, # spin, complx/real, # bands) filename_in=TRIM(dir_in)//'evc01.dat' CALL read_sizes (filename_in, ngw, nbnd(1), nspin, is_cmplx, igwx, scalef, gamma_only) ! IF (nspin==2) THEN + ! ! If needed READ the shape of the wfc object for spin down filename_in=TRIM(dir_in)//'evc02.dat' - CALL read_sizes (filename_in, ngw, nbnd(2), nspin, is_cmplx, igwx, scalef, gamma_only) + CALL read_sizes (filename_in, ngw_, nbnd(2), nspin_, is_cmplx_, igwx_, scalef_, gamma_only_) + ! + ! The following should never happen + IF (ngw_ /= ngw) CALL error("ngw from evc01 and evc02 does not match", ngw, ngw_) + IF (igwx_ /= igwx) CALL error("ngw from evc01 and evc02 does not match", igwx, igwx_) + IF (nspin_ /= nspin) CALL error("nspin from evc01 and evc02 does not match", nspin, nspin_) + IF (is_cmplx_ .NEQV. is_cmplx_) CALL error("is_cmplx from evc01 and evc02 does not match", -1) + IF (gamma_only_ .NEQV. gamma_only) CALL error("gamma_only from evc01 and evc02 does not match", -1) + ! ENDIF ! WRITE(*,'(/, 7X, "SUMMARY: igwx =", I8)') igwx @@ -91,48 +93,47 @@ PROGRAM n2npm1 ! ALLOCATE (evc(igwx, max(nbnd(1), nbnd(2)), nspin)) ! - ! reset filename + ! Read the wfcs + ! Reset filename first filename_in=TRIM(dir_in)//'evc01.dat' DO ispin = 1, nspin + ! IF (ispin==2 ) filename_in=TRIM(dir_in)//'evc02.dat' CALL read_wf_occ (filename_in, evc(:,:,ispin), ngw, nbnd(ispin)) - ! TODO: if KS state are used at initialization, also empty states are written to evc + ! + ! FIXME: if KS state are used at initialization, also empty states are written to evc ! then this countig does not make sense CALL check_nele (ngw, nbnd(ispin), evc(:,:,ispin), nel(ispin)) + ! WRITE(*,'(7X, "CHECK: ispin =", I5, " nel =", F18.12, " INT(nel)=", I5)') & ispin, nel(ispin), NINT(nel(ispin)) - !WRITE(*,) ispin, nel(ispin), NINT(nel(ispin)) + ! ENDDO ! ! At this stage I read the empty manifold if needed IF (l_fill) THEN + ! + ! If needed we read here empty_state WFCs WRITE(*,'(/ 5X, A)') "Reading EMP manifold ..." - ! Here we need to read also empty states if needed. - ! First the shape of the arrey ngw, nbnd_emp for spin up + ! + ! First the shape of the arrray (ngw, nbnd_emp) for spin up iuni=789 filename_in=TRIM(dir_in)//'evc0_empty1.dat' - INQUIRE (FILE=TRIM(filename_in), EXIST=EXST) - IF (.NOT. exst) THEN - WRITE(*,'(3A)') "File ", TRIM(filename_in), " do not exist" - STOP - ENDIF - OPEN (UNIT=iuni, FILE=TRIM(filename_in), status='unknown', FORM='UNFORMATTED') - READ (iuni) ngw, nbnd_emp(1) + CALL read_size_emp(filename_in, iuni, ngw_, nbnd_emp(1)) + ! ! then for spin down channel IF (nspin == 2) THEN + ! iuni=iuni+1 filename_in=TRIM(dir_in)//'evc0_empty2.dat' - INQUIRE (FILE=TRIM(filename_in), EXIST=EXST) - IF (.NOT. exst) THEN - WRITE(*,'(3A)') "File ", TRIM(filename_in), " do not exist" - STOP - ENDIF - OPEN (UNIT=iuni, FILE=TRIM(filename_in), status='unknown', FORM='UNFORMATTED') - READ (iuni) ngw, nbnd_emp(2) + CALL read_size_emp(filename_in, iuni, ngw_, nbnd_emp(2)) + ! ENDIF + ! WRITE(*,'(/, 7X, "SUMMARY: ngw =", I8)') ngw WRITE(*,'( 7X, " nbnd(1) =", I8)') nbnd_emp(1) WRITE(*,'( 7X, " nbnd(2) =", I8)') nbnd_emp(2) + ! ! allocation ALLOCATE ( evc_empty (ngw, MAX(nbnd_emp(1), nbnd_emp(2)),nspin) ) ! @@ -147,25 +148,24 @@ PROGRAM n2npm1 ENDDO ENDDO ENDIF + CLOSE (iuni) + IF (nspin==2) CLOSE (iuni+1) ! ! check if the index of the orbital is within a valid range IF (l_fill) THEN ! - IF (spin_channel == 1 .AND. index .gt. nbnd_emp(1)) THEN - WRITE(*,*) "index out of range ", index, nbnd_emp(1); STOP - ENDIF - IF (spin_channel == 2 .AND. index .gt. nbnd_emp(2)) THEN - WRITE(*,*) "index out of range ", index, nbnd_emp(2); STOP - ENDIF + IF (spin_channel == 1 .AND. index .gt. nbnd_emp(1)) & + & CALL error ("index out of range", index, nbnd_emp(1)) + ! + IF (spin_channel == 2 .AND. index .gt. nbnd_emp(2)) & + & CALL error ("index out of range", index, nbnd_emp(2)) ! ELSE ! - IF (spin_channel == 1 .AND. index .gt. nbnd(1)) THEN - WRITE(*,*) "index out of range ", index, nbnd(1); STOP - ENDIF - IF (spin_channel == 2 .AND. index .gt. nbnd(2)) THEN - WRITE(*,*) "index out of range ", index, nbnd(2); STOP - ENDIF + IF (spin_channel == 1 .AND. index .gt. nbnd(1)) & + & CALL error ("index out of range", index, nbnd(1)) + IF (spin_channel == 2 .AND. index .gt. nbnd(2)) & + & CALL error ("index out of range", index, nbnd(2)) ! ENDIF ! @@ -175,9 +175,8 @@ PROGRAM n2npm1 IF (l_fill) THEN nbnd_out(spin_channel) = nbnd(spin_channel) + 1 nel(spin_channel) = nel(spin_channel)+1 - IF (nel(1) .lt. nel(2)) THEN - WRITE(*,*) "nel up must be greater than or equal to nel dw ", NINT(nel(:)); STOP - ENDIF + IF (nel(1) .lt. nel(2)) CALL error("nel up must be greater than or equal to nel dw", & + & NINT(nel(1)),NINT(nel(2))) ALLOCATE (evc_out(igwx, max(nbnd_out(1), nbnd_out(2)), nspin)) evc_out = CMPLX(0.D0, 0.D0, kind =DP) DO ispin = 1, nspin @@ -190,9 +189,8 @@ PROGRAM n2npm1 ELSE nbnd_out(spin_channel) = nbnd(spin_channel) - 1 nel(spin_channel) = nel(spin_channel)-1 - IF (nel(1) .lt. nel(2)) THEN - WRITE(*,*) "nel up must be greater than or equal to nel dw ", NINT(nel(:)); STOP - ENDIF + IF (nel(1) .lt. nel(2)) CALL error("nel up must be greater than or equal to nel dw", & + & NINT(nel(1)),NINT(nel(2))) ALLOCATE (evc_out(igwx, max(nbnd_out(1), nbnd_out(2)), nspin)) evc_out = CMPLX(0.D0, 0.D0, kind =DP) DO ispin = 1, nspin @@ -202,7 +200,7 @@ PROGRAM n2npm1 ibnd_=ibnd_+1 evc_out(:,ibnd_,ispin) = evc(:,ibnd,ispin) ENDDO - WRITE(*,*) ispin, evc_out(1,1,ispin), evc_out(igwx,MAX(nbnd_out(1),nbnd_out(2)),ispin)! check debug + !WRITE(*,*) ispin, evc_out(1,1,ispin), evc_out(igwx,MAX(nbnd_out(1),nbnd_out(2)),ispin)! check debug ENDDO ENDIF ! @@ -227,6 +225,47 @@ PROGRAM n2npm1 CONTAINS ! !---------------------------------------------------------------------- + SUBROUTINE read_size_emp(filename, iuni, ngw, nbnd) + !------------------------------------------------------------------ + IMPLICIT NONE + INTEGER, INTENT(IN) :: iuni + CHARACTER (256), INTENT(IN):: filename + INTEGER, INTENT(OUT) :: ngw, nbnd + LOGICAL :: exst + ! + ! + INQUIRE (FILE=TRIM(filename), EXIST=EXST) + IF (.NOT. exst) CALL error ("File not found", -1) + ! + OPEN (UNIT=iuni, FILE=TRIM(filename), status='unknown', FORM='UNFORMATTED') + READ (iuni) ngw_, nbnd + IF (ngw_ /= ngw) CALL error("ngw from evc0 and evc0_empty does not match", ngw, ngw_) + RETURN + ! + END SUBROUTINE + ! + !---------------------------------------------------------------------- + SUBROUTINE error (mess, ierr, ierr2) + !------------------------------------------------------------------ + ! + IMPLICIT NONE + CHARACTER (LEN=*), INTENT(IN) :: mess + INTEGER, INTENT(IN) :: ierr + INTEGER, INTENT(IN), OPTIONAL :: ierr2 + ! + IF (PRESENT (ierr2)) THEN + WRITE(*,'(/,3X, 60("%"), /, 5X, 2A, 2I5, /, 3X, 60("%"),/)') "ERROR: ", TRIM(mess), ierr, ierr2 + ELSE + WRITE(*,'(/,3X, 60("%"), /, 5X, 2A, 1I5, /, 3X, 60("%"),/)') "ERROR: ", TRIM(mess), ierr + ENDIF + ! + STOP + RETURN + ! + END SUBROUTINE + ! + ! + !---------------------------------------------------------------------- SUBROUTINE read_command_line (task, index, spin_channel, dir_in, dir_out) !------------------------------------------------------------------ ! arg# 1 2 3 4 5 @@ -244,6 +283,12 @@ SUBROUTINE read_command_line (task, index, spin_channel, dir_in, dir_out) IF ( nargs /= 5 ) THEN WRITE(*,'(/, 5X, A, I5)') "Wrong number of arguments: STOP", nargs WRITE(*,'( 5X, A, /)') "Usage: n2npm1.x -- " + WRITE(*,'( 5X, A, /)') " = fill/empty" + WRITE(*,'( 5X, A )') " index of the occ orbital to remove form the occ manifold (==empty)" + WRITE(*,'( 5X, A, /)') " index of the empty orbital to add to the occ manifold (==fill)" + WRITE(*,'( 5X, A, /)') " = 1/2" + WRITE(*,'( 5X, A, /)') " Directory where to read N-electron WFC files" + WRITE(*,'( 5X, A, /)') " Directory where to write N+1(N-1)-electron WFC files" STOP ENDIF ! @@ -259,8 +304,7 @@ SUBROUTINE read_command_line (task, index, spin_channel, dir_in, dir_out) CASE ( '--empty') task="empty" CASE DEFAULT - WRITE(*,'(A)') 'unrecognised argument option' - STOP + CALL error('unrecognised argument option', -1) END SELECT ! CALL get_command_argument( iarg, arg ) @@ -273,6 +317,18 @@ SUBROUTINE read_command_line (task, index, spin_channel, dir_in, dir_out) iarg = iarg + 1 CALL get_command_argument( iarg, dir_out ) ! + ! + WRITE(*,'(/, 3X, A )') "INPUT SUMMARY:" + WRITE(*,'( 3X, " task = ", A )') TRIM(task) + WRITE(*,'( 3X, " index = ", I5)') index + WRITE(*,'( 3X, " spin = ", I5)') spin_channel + WRITE(*,'( 3X, " dir_in = ", A )') TRIM(dir_in) + WRITE(*,'( 3X, " dir_out = ", A )') TRIM(dir_out) + ! + IF (dir_in == dir_out) CALL error ("dir_out = dir_in", -1) + IF (spin_channel >2 .OR. spin_channel < 0) CALL error("spin channel out of range", spin_channel) + IF (index < 0 ) CALL error("index out of range", index) + ! RETURN ! END SUBROUTINE @@ -335,7 +391,7 @@ SUBROUTINE check_nele (ngw, nbnd, evc, nele) ! nele=0.D0 DO ibnd = 1, nbnd - nele = nele + SUM(evc(:,ibnd)*CONJG(evc(:,ibnd))) + nele = nele + REAL(SUM(evc(:,ibnd)*CONJG(evc(:,ibnd)))) !WRITE(*,'(I5,3x , 2F20.15, 3x, F10.6)') ibnd, SUM(evc(:,ibnd)*CONJG(evc(:,ibnd))), nele ENDDO RETURN @@ -355,7 +411,7 @@ SUBROUTINE read_sizes(filename, ngw, nbnd, nspin, is_cmplx, igwx, scalef, gamma_ INTEGER, PARAMETER :: DP = selected_real_kind(14,200) CHARACTER(iotk_attlenx) :: attr - INTEGER :: j, ispin + INTEGER :: ispin INTEGER :: ierr INTEGER :: ik, nk INTEGER :: iuni