From 13e61c21106f75a797ff6f270bd274dc5e72b771 Mon Sep 17 00:00:00 2001 From: giannozz Date: Fri, 11 Jun 2010 09:00:04 +0000 Subject: [PATCH] Source "normalization" using Norbert's script dev-tools/src-normal. Changes should not affect functionalities (please verify!) git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6825 c92efa57-630b-4861-b058-cf58834340f0 --- upftools/any2upf.f90 | 80 ++-- upftools/casino2upf.f90 | 162 ++++---- upftools/cpmd2upf.f90 | 426 +++++++++---------- upftools/fhi2upf.f90 | 392 +++++++++--------- upftools/fpmd2upf.f90 | 526 ++++++++++++------------ upftools/ncpp2upf.f90 | 406 +++++++++---------- upftools/oldcp2upf.f90 | 230 +++++------ upftools/read_upf.f90 | 448 ++++++++++---------- upftools/read_upf_tofile.f90 | 32 +- upftools/rrkj2upf.f90 | 308 +++++++------- upftools/uspp2upf.f90 | 32 +- upftools/vanderbilt.f90 | 414 +++++++++---------- upftools/vdb2upf.f90 | 28 +- upftools/virtual.f90 | 766 +++++++++++++++++------------------ upftools/write_upf.f90 | 506 +++++++++++------------ 15 files changed, 2378 insertions(+), 2378 deletions(-) diff --git a/upftools/any2upf.f90 b/upftools/any2upf.f90 index db67a42c0..0a75dafea 100644 --- a/upftools/any2upf.f90 +++ b/upftools/any2upf.f90 @@ -7,72 +7,72 @@ ! ! !--------------------------------------------------------------------- -program mypp2upf +PROGRAM mypp2upf !--------------------------------------------------------------------- ! - ! Convert a pseudopotential written in a user-supplied format + ! Convert a pseudopotential written in a user-supplied format ! to unified pseudopotential format - sample program ! - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open (unit = 1, file = filein, status = 'old', form = 'formatted') - call read_mypp(1) - close (1) + CALL get_file ( filein ) + OPEN (unit = 1, file = filein, status = 'old', form = 'formatted') + CALL read_mypp(1) + CLOSE (1) ! convert variables read from user-supplied format into those needed ! by the upf format - add missing quantities - call convert_mypp + CALL convert_mypp fileout=trim(filein)//'.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) -stop -20 write (6,'("mypp2upf: error reading pseudopotential file name")') - stop -end program mypp2upf +STOP +20 WRITE (6,'("mypp2upf: error reading pseudopotential file name")') + STOP +END PROGRAM mypp2upf -module mypp +MODULE mypp ! ! All variables read from user-supplied file format ! Must have distinct names from variables in the "upf" module ! -end module mypp -! +END MODULE mypp +! ! ---------------------------------------------------------- -subroutine read_mypp(iunps) - ! ---------------------------------------------------------- - ! - use mypp - implicit none - integer :: iunps - ! - ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully read' +SUBROUTINE read_mypp(iunps) ! ---------------------------------------------------------- ! - return -100 write (6,'("read_mypp: error reading pseudopotential file")') - stop -end subroutine read_mypp + USE mypp + IMPLICIT NONE + INTEGER :: iunps + ! + ! ---------------------------------------------------------- + WRITE (6,'(a)') 'Pseudopotential successfully read' + ! ---------------------------------------------------------- + ! + RETURN +100 WRITE (6,'("read_mypp: error reading pseudopotential file")') + STOP +END SUBROUTINE read_mypp ! ---------------------------------------------------------- -subroutine convert_mypp +SUBROUTINE convert_mypp ! ---------------------------------------------------------- ! - use mypp - use upf - implicit none + USE mypp + USE upf + IMPLICIT NONE ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully converted' + WRITE (6,'(a)') 'Pseudopotential successfully converted' ! ---------------------------------------------------------- - return -end subroutine convert_mypp + RETURN +END SUBROUTINE convert_mypp diff --git a/upftools/casino2upf.f90 b/upftools/casino2upf.f90 index 7023d41ae..b10fa08a2 100644 --- a/upftools/casino2upf.f90 +++ b/upftools/casino2upf.f90 @@ -29,9 +29,9 @@ PROGRAM casino2upf PRINT*, 'Enter wavefunction files, starting with the ground state:' DO i=1,nofiles CALL get_file ( wavefile(i) ) - OPEN(unit=i,file=wavefile(i),status='old',form='formatted') + OPEN(unit=i,file=wavefile(i),status='old',form='formatted') ENDDO - OPEN(unit=99,file=filein,status='old',form='formatted') + OPEN(unit=99,file=filein,status='old',form='formatted') CALL read_casino(99,nofiles) CLOSE (unit=99) @@ -44,7 +44,7 @@ PROGRAM casino2upf CALL convert_casino - fileout=TRIM(filein)//'.UPF' + fileout=trim(filein)//'.UPF' PRINT '(''Output PP file in US format : '',a)', fileout OPEN(unit=2,file=fileout,status='unknown',form='formatted') @@ -56,10 +56,10 @@ PROGRAM casino2upf END PROGRAM casino2upf MODULE casino - + ! ! All variables read from CASINO file format - ! + ! ! trailing underscore means that a variable with the same name ! is used in module 'upf' containing variables to be written ! @@ -82,17 +82,17 @@ MODULE casino REAL(8), ALLOCATABLE:: chi_(:,:), oc_(:) END MODULE casino -! +! ! ---------------------------------------------------------- SUBROUTINE read_casino(iunps,nofiles) ! ---------------------------------------------------------- - ! + ! USE casino USE upf , ONLY : els USE kinds IMPLICIT NONE TYPE :: wavfun_list - INTEGER :: occ,eup,edwn, nquant, lquant + INTEGER :: occ,eup,edwn, nquant, lquant CHARACTER(len=2) :: label #ifdef __STD_F95 REAL*8, POINTER :: wavefunc(:) @@ -112,7 +112,7 @@ SUBROUTINE read_casino(iunps,nofiles) LOGICAL :: groundstate, found CHARACTER(len=1), DIMENSION(0:3) :: convel=(/'s','p','d','f'/) CHARACTER(len=2) :: label, rellab - REAL(DP), parameter :: r_exp=20._dp/1500._dp + REAL(DP), PARAMETER :: r_exp=20._dp/1500._dp INTEGER :: l, i, ir, nb, gsorbs, j,k,m,tmp, lquant, orbs, nquant INTEGER, ALLOCATABLE :: gs(:,:) @@ -122,31 +122,31 @@ SUBROUTINE read_casino(iunps,nofiles) nlc = 0 !These two values are always 0 for numeric pps nnl = 0 !so lets just hard code them - nlcc_ = .FALSE. !Again these two are alwas false for CASINO pps - bhstype = .FALSE. + nlcc_ = .false. !Again these two are alwas false for CASINO pps + bhstype = .false. - - READ(iunps,'(a2,35x,a2)') rellab, psd_ + + READ(iunps,'(a2,35x,a2)') rellab, psd_ READ(iunps,*) - IF ( rellab .EQ. 'DF' ) THEN + IF ( rellab == 'DF' ) THEN rel_=1 ELSE rel_=0 ENDIF - + READ(iunps,*) zmesh,zp_ !Here we are reading zmesh (atomic #) and DO i=1,3 !zp_ (pseudo charge) READ(iunps,*) ENDDO READ(iunps,*) lmax_ !reading in lmax - IF ( zp_.LE.0d0 ) & + IF ( zp_<=0d0 ) & CALL errore( 'read_casino','Wrong zp ',1 ) - IF ( lmax_.GT.3.OR.lmax_.LT.0 ) & + IF ( lmax_>3.or.lmax_<0 ) & CALL errore( 'read_casino','Wrong lmax ',1 ) lloc=lmax_ !think lloc shoudl always = lmax for this case yes/no ?? - + ! ! compute the radial mesh ! @@ -165,7 +165,7 @@ SUBROUTINE read_casino(iunps,nofiles) ENDDO DO ir = 1, mesh_ rab_(ir) = r_exp * r_(ir) !hardcoded at the moment - END DO + ENDDO ALLOCATE(vnl(mesh_,0:lmax_)) @@ -179,7 +179,7 @@ SUBROUTINE read_casino(iunps,nofiles) DO ir = 1, mesh_ vnl(ir,l) = vnl(ir,l)/r_(ir) !Removing the factor of r CASINO has ENDDO - vnl(1,l) = 0 !correcting for the divide by zero + vnl(1,l) = 0 !correcting for the divide by zero ENDDO ALLOCATE(rho_atc_(mesh_)) @@ -192,7 +192,7 @@ SUBROUTINE read_casino(iunps,nofiles) ! DO l = 0, lmax_ - IF ( l.NE.lloc ) vnl(:,l) = vnl(:,l) - vnl(:,lloc) + IF ( l/=lloc ) vnl(:,l) = vnl(:,l) - vnl(:,lloc) ENDDO @@ -202,7 +202,7 @@ SUBROUTINE read_casino(iunps,nofiles) mptr => mhead NULLIFY(mtail%p) - groundstate=.TRUE. + groundstate=.true. DO j=1,nofiles DO i=1,4 @@ -212,13 +212,13 @@ SUBROUTINE read_casino(iunps,nofiles) READ(j,*) orbs - IF ( groundstate ) THEN + IF ( groundstate ) THEN ALLOCATE( gs(orbs,3) ) gs = 0 gsorbs = orbs - END IF + ENDIF DO i=1,2 READ(j,*) @@ -231,20 +231,20 @@ SUBROUTINE read_casino(iunps,nofiles) READ(j,*) tmp, nquant, lquant IF ( groundstate ) THEN - found = .TRUE. + found = .true. DO m=1,orbs - IF ( (nquant==gs(m,1) .AND. lquant==gs(m,2)) ) THEN + IF ( (nquant==gs(m,1) .and. lquant==gs(m,2)) ) THEN gs(m,3) = gs(m,3) + 1 - EXIT - END IF + exit + ENDIF - found = .FALSE. + found = .false. ENDDO - IF (.NOT. found ) THEN + IF (.not. found ) THEN DO m=1,orbs @@ -253,7 +253,7 @@ SUBROUTINE read_casino(iunps,nofiles) gs(m,2) = lquant gs(m,3) = 1 - EXIT + exit ENDIF ENDDO @@ -265,9 +265,9 @@ SUBROUTINE read_casino(iunps,nofiles) ENDDO READ(j,*) - READ(j,*) + READ(j,*) - DO i=1,mesh_ + DO i=1,mesh_ READ(j,*) ENDDO @@ -275,34 +275,34 @@ SUBROUTINE read_casino(iunps,nofiles) READ(j,'(13x,a2)', err=300) label READ(j,*) tmp, nquant, lquant - IF ( .NOT. groundstate ) THEN - found = .FALSE. + IF ( .not. groundstate ) THEN + found = .false. DO m = 1,gsorbs - IF ( nquant == gs(m,1) .AND. lquant == gs(m,2) ) THEN - found = .TRUE. - EXIT - END IF - END DO + IF ( nquant == gs(m,1) .and. lquant == gs(m,2) ) THEN + found = .true. + exit + ENDIF + ENDDO mptr => mhead DO - IF ( .NOT. ASSOCIATED(mptr) )EXIT - IF ( nquant == mptr%nquant .AND. lquant == mptr%lquant ) found = .TRUE. + IF ( .not. associated(mptr) )exit + IF ( nquant == mptr%nquant .and. lquant == mptr%lquant ) found = .true. mptr =>mptr%p - END DO + ENDDO IF ( found ) THEN DO i=1,mesh_ READ(j,*) ENDDO CYCLE - END IF - END IF + ENDIF + ENDIF #ifdef __STD_F95 - IF ( ASSOCIATED(mtail%wavefunc) ) THEN + IF ( associated(mtail%wavefunc) ) THEN #else - IF ( ALLOCATED(mtail%wavefunc) ) THEN + IF ( allocated(mtail%wavefunc) ) THEN #endif ALLOCATE(mtail%p) mtail=>mtail%p @@ -310,7 +310,7 @@ SUBROUTINE read_casino(iunps,nofiles) ALLOCATE( mtail%wavefunc(mesh_) ) ELSE ALLOCATE( mtail%wavefunc(mesh_) ) - END IF + ENDIF mtail%label = label mtail%nquant = nquant mtail%lquant = lquant @@ -318,17 +318,17 @@ SUBROUTINE read_casino(iunps,nofiles) READ(j, *, err=300) (mtail%wavefunc(ir),ir=1,mesh_) ENDDO - groundstate = .FALSE. + groundstate = .false. ENDDO nchi =0 mptr => mhead DO - IF ( .NOT. ASSOCIATED(mptr) )EXIT + IF ( .not. associated(mptr) )exit nchi=nchi+1 mptr =>mptr%p - END DO + ENDDO ALLOCATE(lchi_(nchi), els(nchi), nns_(nchi)) ALLOCATE(oc_(nchi)) @@ -336,16 +336,16 @@ SUBROUTINE read_casino(iunps,nofiles) oc_ = 0 !Sort out the occupation numbers - DO i=1,gsorbs + DO i=1,gsorbs oc_(i)=gs(i,3) ENDDO - deallocate( gs ) + DEALLOCATE( gs ) i=1 mptr => mhead DO - IF ( .NOT. ASSOCIATED(mptr) )EXIT - nns_(i) = mptr%nquant + IF ( .not. associated(mptr) )exit + nns_(i) = mptr%nquant lchi_(i) = mptr%lquant els(i) = mptr%label @@ -353,18 +353,18 @@ SUBROUTINE read_casino(iunps,nofiles) chi_(ir:,i) = mptr%wavefunc(ir) ENDDO - deallocate( mptr%wavefunc ) + DEALLOCATE( mptr%wavefunc ) mptr =>mptr%p i=i+1 - END DO + ENDDO !Clean up the linked list (deallocate it) DO - IF ( .NOT. ASSOCIATED(mhead) )EXIT + IF ( .not. associated(mhead) )exit mptr => mhead mhead => mhead%p - deallocate( mptr ) - END DO + DEALLOCATE( mptr ) + ENDDO ! @@ -373,9 +373,9 @@ SUBROUTINE read_casino(iunps,nofiles) ALLOCATE(rho_at_(mesh_)) rho_at_(:)=0.d0 DO nb = 1, nchi - IF( oc_(nb).NE.0.d0) & + IF( oc_(nb)/=0.d0) & & rho_at_(:) = rho_at_(:) + oc_(nb)*chi_(:,nb)**2 - END DO + ENDDO ! ---------------------------------------------------------- WRITE (6,'(a)') 'Pseudopotential successfully read' ! ---------------------------------------------------------- @@ -405,7 +405,7 @@ SUBROUTINE convert_casino rcloc = 0.0d0 - nwfs = nchi + nwfs = nchi ALLOCATE( oc(nwfs), epseu(nwfs)) ALLOCATE(lchi(nwfs), nns(nwfs) ) ALLOCATE(rcut (nwfs), rcutus (nwfs)) @@ -416,7 +416,7 @@ SUBROUTINE convert_casino rcutus(i)= 0.0d0 oc (i) = oc_(i) epseu(i) = 0.0d0 - END DO + ENDDO DEALLOCATE (lchi_, oc_, nns_) psd = psd_ @@ -430,8 +430,8 @@ SUBROUTINE convert_casino lmax = lmax_-1 ELSE lmax = lmax_ - END IF - nbeta= lmax_ + ENDIF + nbeta= lmax_ mesh = mesh_ ntwfc= nchi ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) @@ -439,7 +439,7 @@ SUBROUTINE convert_casino lchiw(i) = lchi(i) ocw(i) = oc(i) elsw(i) = els(i) - END DO + ENDDO CALL set_dft_from_name(dft_) iexch = get_iexch() icorr = get_icorr() @@ -465,12 +465,12 @@ SUBROUTINE convert_casino DO ir = 1,mesh IF ( r(ir) > rmax ) THEN kkbeta=ir - EXIT - END IF - END DO + exit + ENDIF + ENDDO ! make sure kkbeta is odd as required for simpson - IF(MOD(kkbeta,2) == 0) kkbeta=kkbeta-1 + IF(mod(kkbeta,2) == 0) kkbeta=kkbeta-1 ikk2(:) = kkbeta ALLOCATE(aux(kkbeta)) ALLOCATE(betar(mesh,nbeta)) @@ -483,18 +483,18 @@ SUBROUTINE convert_casino iv=0 DO i=1,nchi l=lchi(i) - IF (l.NE.lloc) THEN + IF (l/=lloc) THEN iv=iv+1 lll(iv)=l DO ir=1,kkbeta betar(ir,iv)=chi_(ir,i)*vnl(ir,l) aux(ir) = chi_(ir,i)**2*vnl(ir,l) - END DO + ENDDO CALL simpson(kkbeta,aux,rab,vll) dion(iv,iv) = 1.0d0/vll - END IF - IF(iv >= nbeta) EXIT ! skip additional pseudo wfns + ENDIF + IF(iv >= nbeta) exit ! skip additional pseudo wfns ENDDO @@ -506,13 +506,13 @@ SUBROUTINE convert_casino DO iv=1,nbeta ikk2(iv)=kkbeta DO ir = kkbeta,1,-1 - IF ( ABS(betar(ir,iv)) > 1.d-12 ) THEN + IF ( abs(betar(ir,iv)) > 1.d-12 ) THEN ikk2(iv)=ir - EXIT - END IF - END DO - END DO - END IF + exit + ENDIF + ENDDO + ENDDO + ENDIF ALLOCATE (rho_at(mesh)) rho_at = rho_at_ DEALLOCATE (rho_at_) diff --git a/upftools/cpmd2upf.f90 b/upftools/cpmd2upf.f90 index 855809dea..c7781e892 100644 --- a/upftools/cpmd2upf.f90 +++ b/upftools/cpmd2upf.f90 @@ -7,236 +7,236 @@ ! ! !--------------------------------------------------------------------- -program cpmd2upf +PROGRAM cpmd2upf !--------------------------------------------------------------------- ! ! Convert a pseudopotential written in the CPMD format ! (TYPE=NORMCONSERVING NUMERIC only, single radial grid) ! to unified pseudopotential format ! - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open (unit = 1, file = filein, status = 'old', form = 'formatted') - call read_cpmd(1) - close (1) + CALL get_file ( filein ) + OPEN (unit = 1, file = filein, status = 'old', form = 'formatted') + CALL read_cpmd(1) + CLOSE (1) ! convert variables read from CPMD format into those needed ! by the upf format - add missing quantities - call convert_cpmd + CALL convert_cpmd fileout=trim(filein)//'.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) -stop -20 call errore ('cpmd2upf', 'Reading pseudo file name ', 1) +STOP +20 CALL errore ('cpmd2upf', 'Reading pseudo file name ', 1) -end program cpmd2upf +END PROGRAM cpmd2upf -module cpmd +MODULE cpmd ! ! All variables read from CPMD file format ! - character (len=80) title + CHARACTER (len=80) title ! - integer :: ixc + INTEGER :: ixc real(8) :: alphaxc - integer :: z, zv + INTEGER :: z, zv ! - integer :: mesh_ + INTEGER :: mesh_ real(8) :: amesh, amesh_ - real(8), allocatable :: r_(:) + real(8), ALLOCATABLE :: r_(:) ! - integer ::lmax_ - real(8), allocatable :: vnl(:,:) - real(8), allocatable :: chi_(:,:) + INTEGER ::lmax_ + real(8), ALLOCATABLE :: vnl(:,:) + real(8), ALLOCATABLE :: chi_(:,:) ! - logical :: nlcc_ - real(8), allocatable :: rho_atc_(:) + LOGICAL :: nlcc_ + real(8), ALLOCATABLE :: rho_atc_(:) ! - integer :: maxinfo_, info_lines_ - parameter (maxinfo_ = 100) - character (len=80), allocatable :: info_sect_(:) + INTEGER :: maxinfo_, info_lines_ + PARAMETER (maxinfo_ = 100) + CHARACTER (len=80), ALLOCATABLE :: info_sect_(:) !------------------------------ -end module cpmd -! +END MODULE cpmd +! ! ---------------------------------------------------------- -subroutine read_cpmd(iunps) +SUBROUTINE read_cpmd(iunps) ! ---------------------------------------------------------- - ! - use cpmd - implicit none - integer :: iunps - ! - integer :: found = 0, closed = 0, unknown = 0 - integer :: i, l, ios - character (len=80) line - character (len=4) token + ! + USE cpmd + IMPLICIT NONE + INTEGER :: iunps + ! + INTEGER :: found = 0, closed = 0, unknown = 0 + INTEGER :: i, l, ios + CHARACTER (len=80) line + CHARACTER (len=4) token real (8) :: vnl0(0:3) - logical, external :: matches - integer, external :: locate + LOGICAL, EXTERNAL :: matches + INTEGER, EXTERNAL :: locate ! nlcc_ = .false. info_lines_ = 0 -10 read (iunps,'(A)',end=20,err=20) line - if (matches ("&ATOM", trim(line)) ) then +10 READ (iunps,'(A)',end=20,err=20) line + IF (matches ("&ATOM", trim(line)) ) THEN found = found + 1 ! Z - read (iunps,'(a)',end=200,err=200) line + READ (iunps,'(a)',end=200,err=200) line l = len_trim(line) i = locate('=',line) - read (line(i+1:l),*) z + READ (line(i+1:l),*) z ! ZV - read (iunps,'(a)',end=200,err=200) line + READ (iunps,'(a)',end=200,err=200) line l = len_trim(line) i = locate('=',line) - read (line(i+1:l),*) zv + READ (line(i+1:l),*) zv ! XC - read (iunps,'(a)',end=200,err=200) line + READ (iunps,'(a)',end=200,err=200) line l = len_trim(line) i = locate('=',line) - read (line(i+1:l),*) ixc, alphaxc + READ (line(i+1:l),*) ixc, alphaxc ! TYPE - read (iunps,'(a)',end=200,err=200) line - if (.not. matches("NORMCONSERVING",line) .or. & + READ (iunps,'(a)',end=200,err=200) line + IF (.not. matches("NORMCONSERVING",line) .or. & .not. matches("NUMERIC",line) ) & - call errore('read_cpmd','unknown type: '//line,1) - else if (matches ("&INFO", trim(line)) ) then + CALL errore('read_cpmd','unknown type: '//line,1) + ELSEIF (matches ("&INFO", trim(line)) ) THEN found = found + 1 ! read (iunps,'(a)') title ! store info section for later perusal (FIXME: not yet implemented. 2004/10/12, AK) - allocate (info_sect_(maxinfo_)) - do i=1,maxinfo_ - read (iunps,'(a)',end=20,err=20) title - if (matches ("&END", trim(title)) ) then + ALLOCATE (info_sect_(maxinfo_)) + DO i=1,maxinfo_ + READ (iunps,'(a)',end=20,err=20) title + IF (matches ("&END", trim(title)) ) THEN closed = closed + 1 - goto 10 - else + GOTO 10 + ELSE info_sect_(i) = trim(title) info_lines_ = i - end if - enddo - else if (matches ("&POTENTIAL", trim(line)) ) then + ENDIF + ENDDO + ELSEIF (matches ("&POTENTIAL", trim(line)) ) THEN found = found + 1 - read (iunps,'(a)') line - read (line,*,iostat=ios) mesh_, amesh - if ( ios /= 0) then - read (line,*,iostat=ios) mesh_ + READ (iunps,'(a)') line + READ (line,*,iostat=ios) mesh_, amesh + IF ( ios /= 0) THEN + READ (line,*,iostat=ios) mesh_ amesh = -1.0d0 - end if - allocate (r_(mesh_)) + ENDIF + ALLOCATE (r_(mesh_)) ! ! determine the number of angular momenta ! - read (iunps, '(a)') line + READ (iunps, '(a)') line ios = 1 lmax_=4 - do while (ios /= 0) + DO WHILE (ios /= 0) lmax_ = lmax_ - 1 - read(line,*,iostat=ios) r_(1),(vnl0(l),l=0,lmax_) - end do - allocate (vnl(mesh_,0:lmax_)) + READ(line,*,iostat=ios) r_(1),(vnl0(l),l=0,lmax_) + ENDDO + ALLOCATE (vnl(mesh_,0:lmax_)) vnl(1,0:lmax_) = vnl0(0:lmax_) - do i=2,mesh_ - read(iunps, *) r_(i),(vnl(i,l),l=0,lmax_) - end do + DO i=2,mesh_ + READ(iunps, *) r_(i),(vnl(i,l),l=0,lmax_) + ENDDO ! get amesh if not available directly, check its value otherwise - print "('Radial grid r(i) has ',i4,' points')", mesh_ - print "('Assuming log radial grid: r(i)=exp[(i-1)*amesh]*r(1), with:')" - if (amesh < 0.0d0) then + PRINT "('Radial grid r(i) has ',i4,' points')", mesh_ + PRINT "('Assuming log radial grid: r(i)=exp[(i-1)*amesh]*r(1), with:')" + IF (amesh < 0.0d0) THEN amesh = log (r_(mesh_)/r_(1))/(mesh_-1) - print "('amesh = log (r(mesh)/r(1))/(mesh-1) = ',f10.6)",amesh - else + PRINT "('amesh = log (r(mesh)/r(1))/(mesh-1) = ',f10.6)",amesh + ELSE ! not clear whether the value of amesh read from file ! matches the above definition, or if it is exp(amesh) ... amesh_ = log (r_(mesh_)/r_(1))/(mesh_-1) - if ( abs ( amesh - amesh_ ) > 1.0d-5 ) then - if ( abs ( amesh - exp(amesh_) ) < 1.0d-5 ) then + IF ( abs ( amesh - amesh_ ) > 1.0d-5 ) THEN + IF ( abs ( amesh - exp(amesh_) ) < 1.0d-5 ) THEN amesh = log(amesh) - print "('amesh = log (value read from file) = ',f10.6)",amesh - else - call errore ('cpmd2upf', 'unknown real-space grid',2) - end if - else - print "('amesh = value read from file = ',f10.6)",amesh - end if - end if + PRINT "('amesh = log (value read from file) = ',f10.6)",amesh + ELSE + CALL errore ('cpmd2upf', 'unknown real-space grid',2) + ENDIF + ELSE + PRINT "('amesh = value read from file = ',f10.6)",amesh + ENDIF + ENDIF ! check if the grid is what we expect - do i=2,mesh_ - if ( abs(r_(i) - exp((i-1)*amesh)*r_(1)) > 1.0d-5) then - print "('grid point ',i4,': found ',f10.6,', expected ',f10.6)",& + DO i=2,mesh_ + IF ( abs(r_(i) - exp((i-1)*amesh)*r_(1)) > 1.0d-5) THEN + PRINT "('grid point ',i4,': found ',f10.6,', expected ',f10.6)",& i, r_(i), exp((i-1)*amesh)*r_(1) - call errore ('cpmd2upf', 'unknown real-space grid',1) - end if - end do - else if (matches ("&WAVEFUNCTION", trim(line)) ) then + CALL errore ('cpmd2upf', 'unknown real-space grid',1) + ENDIF + ENDDO + ELSEIF (matches ("&WAVEFUNCTION", trim(line)) ) THEN found = found + 1 ! read (iunps,*) mesh_, amesh - read (iunps,'(a)') line - read (line,*,iostat=ios) mesh_ - allocate(chi_(mesh_,lmax_+1)) - do i=1,mesh_ - read(iunps, *) r_(i),(chi_(i,l+1),l=0,lmax_) - end do - else if (matches ("&NLCC", trim(line)) ) then + READ (iunps,'(a)') line + READ (line,*,iostat=ios) mesh_ + ALLOCATE(chi_(mesh_,lmax_+1)) + DO i=1,mesh_ + READ(iunps, *) r_(i),(chi_(i,l+1),l=0,lmax_) + ENDDO + ELSEIF (matches ("&NLCC", trim(line)) ) THEN found = found + 1 nlcc_ = .true. - read (iunps, '(a)') line - if (.not. matches ("NUMERIC", trim(line)) ) & - call errore('read_cpmd',' only NUMERIC core-correction supported',1) - read(iunps, *) mesh_ - allocate (rho_atc_(mesh_)) - read(iunps, * ) (r_(i), rho_atc_(i), i=1,mesh_) - else if (matches ("&ATDENS", trim(line)) ) then + READ (iunps, '(a)') line + IF (.not. matches ("NUMERIC", trim(line)) ) & + CALL errore('read_cpmd',' only NUMERIC core-correction supported',1) + READ(iunps, *) mesh_ + ALLOCATE (rho_atc_(mesh_)) + READ(iunps, * ) (r_(i), rho_atc_(i), i=1,mesh_) + ELSEIF (matches ("&ATDENS", trim(line)) ) THEN ! skip over &ATDENS section, add others here, if there are more. - do while(.not. matches("&END", trim(line))) - read (iunps,'(a)') line - end do - else if (matches ("&END", trim(line)) ) then + DO WHILE(.not. matches("&END", trim(line))) + READ (iunps,'(a)') line + ENDDO + ELSEIF (matches ("&END", trim(line)) ) THEN closed = closed + 1 - else - print*, 'line ignored: ', line + ELSE + PRINT*, 'line ignored: ', line unknown = unknown + 1 - end if - go to 10 + ENDIF + GOTO 10 -20 continue - if (nlcc_ .and. found /= 5 .or. .not.nlcc_ .and. found /= 4) & - call errore('read_cpmd','some &FIELD card missing',found) - if (closed /= found) & - call errore('read_cpmd','some &END card missing',closed) - if (unknown /= 0 ) print '("WARNING: ",i3," cards not read")', unknown +20 CONTINUE + IF (nlcc_ .and. found /= 5 .or. .not.nlcc_ .and. found /= 4) & + CALL errore('read_cpmd','some &FIELD card missing',found) + IF (closed /= found) & + CALL errore('read_cpmd','some &END card missing',closed) + IF (unknown /= 0 ) PRINT '("WARNING: ",i3," cards not read")', unknown - return -200 call errore('read_cpmd','error in reading file',1) + RETURN +200 CALL errore('read_cpmd','error in reading file',1) -end subroutine read_cpmd +END SUBROUTINE read_cpmd ! ---------------------------------------------------------- -subroutine convert_cpmd +SUBROUTINE convert_cpmd ! ---------------------------------------------------------- ! - use cpmd - use upf - implicit none - real(8), parameter :: rmax = 10.0d0 - real(8), allocatable :: aux(:) + USE cpmd + USE upf + IMPLICIT NONE + real(8), PARAMETER :: rmax = 10.0d0 + real(8), ALLOCATABLE :: aux(:) real(8) :: vll - character (len=20):: dft - character (len=2), external :: atom_name - integer :: lloc, kkbeta, my_lmax - integer :: l, i, ir, iv + CHARACTER (len=20):: dft + CHARACTER (len=2), EXTERNAL :: atom_name + INTEGER :: lloc, kkbeta, my_lmax + INTEGER :: l, i, ir, iv ! - write(generated, '("Generated using unknown code")') - write(date_author,'("Author: unknown Generation date: as well")') + WRITE(generated, '("Generated using unknown code")') + WRITE(date_author,'("Author: unknown Generation date: as well")') comment = 'Info: automatically converted from CPMD format' ! NOTE: many CPMD pseudopotentials created with the 'Hamann' code @@ -246,34 +246,34 @@ subroutine convert_cpmd ! we need to be able to ignore that part or the resulting UPF file ! will be useless. so we first print the info section and ask ! for the LMAX to really use. AK 2005/03/30. - do i=1,info_lines_ - print '(A)', info_sect_(i) - enddo - print '("lmax to use. (max.",I2,") > ",$)', lmax_ - read (5,*) my_lmax - if ((my_lmax <= lmax_) .and. (my_lmax >= 0)) lmax_ = my_lmax - print '("l local (max.",I2,") > ",$)', lmax_ - read (5,*) lloc + DO i=1,info_lines_ + PRINT '(A)', info_sect_(i) + ENDDO + PRINT '("lmax to use. (max.",I2,") > ",$)', lmax_ + READ (5,*) my_lmax + IF ((my_lmax <= lmax_) .and. (my_lmax >= 0)) lmax_ = my_lmax + PRINT '("l local (max.",I2,") > ",$)', lmax_ + READ (5,*) lloc ! reasonable assumption - if (z > 18) then + IF (z > 18) THEN rel = 1 - else + ELSE rel = 0 - end if + ENDIF rcloc = 0.0d0 nwfs = lmax_+1 - allocate( els(nwfs), oc(nwfs), epseu(nwfs)) - allocate(lchi(nwfs), nns(nwfs) ) - allocate(rcut (nwfs), rcutus (nwfs)) - do i=1, nwfs - print '("Wavefunction # ",i1,": label, occupancy > ",$)', i - read (5,*) els(i), oc(i) + ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs)) + ALLOCATE(lchi(nwfs), nns(nwfs) ) + ALLOCATE(rcut (nwfs), rcutus (nwfs)) + DO i=1, nwfs + PRINT '("Wavefunction # ",i1,": label, occupancy > ",$)', i + READ (5,*) els(i), oc(i) nns (i) = 0 lchi(i) = i-1 rcut(i) = 0.0d0 rcutus(i)= 0.0d0 epseu(i) = 0.0d0 - end do + ENDDO psd = atom_name (z) pseudotype = 'NC' nlcc = nlcc_ @@ -281,20 +281,20 @@ subroutine convert_cpmd etotps =0.0d0 ecutrho=0.0d0 ecutwfc=0.0d0 - if ( lmax_ == lloc) then + IF ( lmax_ == lloc) THEN lmax = lmax_-1 - else + ELSE lmax = lmax_ - end if + ENDIF nbeta= lmax_ mesh = mesh_ ntwfc= nwfs - allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) - do i=1, nwfs + ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) + DO i=1, nwfs lchiw(i) = lchi(i) ocw(i) = oc(i) elsw(i) = els(i) - end do + ENDDO iexch = ixc/1000 icorr = (ixc-1000*iexch)/100 igcx = (ixc-1000*iexch-100*icorr)/10 @@ -302,91 +302,91 @@ subroutine convert_cpmd ! ! We have igcc=2 (PW91) and 3 (LYP) exchanged wrt CPMD conventions ! - if (igcc.eq.3) then + IF (igcc==3) THEN igcc=2 - else if (igcc.eq.2) then + ELSEIF (igcc==2) THEN igcc=3 - end if + ENDIF - allocate(rab(mesh)) - allocate( r(mesh)) + ALLOCATE(rab(mesh)) + ALLOCATE( r(mesh)) r = r_ rab = r * amesh - allocate (rho_atc(mesh)) - if (nlcc) rho_atc = rho_atc_ + ALLOCATE (rho_atc(mesh)) + IF (nlcc) rho_atc = rho_atc_ - allocate (vloc0(mesh)) + ALLOCATE (vloc0(mesh)) ! the factor 2 converts from Hartree to Rydberg vloc0(:) = vnl(:,lloc)*2.d0 - if (nbeta > 0) then + IF (nbeta > 0) THEN - allocate(ikk2(nbeta), lll(nbeta)) + ALLOCATE(ikk2(nbeta), lll(nbeta)) kkbeta=mesh - do ir = 1,mesh - if ( r(ir) > rmax ) then + DO ir = 1,mesh + IF ( r(ir) > rmax ) THEN kkbeta=ir exit - end if - end do + ENDIF + ENDDO ikk2(:) = kkbeta - allocate(aux(kkbeta)) - allocate(betar(mesh,nbeta)) - allocate(qfunc(mesh,nbeta,nbeta)) - allocate(dion(nbeta,nbeta)) - allocate(qqq (nbeta,nbeta)) + ALLOCATE(aux(kkbeta)) + ALLOCATE(betar(mesh,nbeta)) + ALLOCATE(qfunc(mesh,nbeta,nbeta)) + ALLOCATE(dion(nbeta,nbeta)) + ALLOCATE(qqq (nbeta,nbeta)) qfunc(:,:,:)=0.0d0 dion(:,:) =0.d0 qqq(:,:) =0.d0 iv=0 - do i=1,nwfs + DO i=1,nwfs l=lchi(i) - if (l.ne.lloc) then + IF (l/=lloc) THEN iv=iv+1 lll(iv)=l - do ir=1,kkbeta + DO ir=1,kkbeta ! the factor 2 converts from Hartree to Rydberg betar(ir,iv) = 2.d0 * chi_(ir,l+1) * & ( vnl(ir,l) - vnl(ir,lloc) ) aux(ir) = chi_(ir,l+1) * betar(ir,iv) - end do - call simpson(kkbeta,aux,rab,vll) + ENDDO + CALL simpson(kkbeta,aux,rab,vll) dion(iv,iv) = 1.0d0/vll - end if - enddo + ENDIF + ENDDO - end if + ENDIF - allocate (rho_at(mesh)) + ALLOCATE (rho_at(mesh)) rho_at = 0.d0 - do i=1,nwfs + DO i=1,nwfs rho_at(:) = rho_at(:) + ocw(i) * chi_(:,i) ** 2 - end do - - allocate (chi(mesh,ntwfc)) + ENDDO + + ALLOCATE (chi(mesh,ntwfc)) chi = chi_ ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully converted' + WRITE (6,'(a)') 'Pseudopotential successfully converted' ! ---------------------------------------------------------- - return -end subroutine convert_cpmd + RETURN +END SUBROUTINE convert_cpmd ! ! ------------------------------------------------------------------ -integer function locate(onechar,string) +INTEGER FUNCTION locate(onechar,string) ! ------------------------------------------------------------------ ! - character(len=1) :: onechar - character(len=*) :: string + CHARACTER(len=1) :: onechar + CHARACTER(len=*) :: string ! - integer:: i + INTEGER:: i ! - do i=1,len_trim(string) - if (string(i:i) .eq. "=") then + DO i=1,len_trim(string) + IF (string(i:i) == "=") THEN locate = i - return - end if - end do + RETURN + ENDIF + ENDDO locate = 0 - return -end function locate + RETURN +END FUNCTION locate diff --git a/upftools/fhi2upf.f90 b/upftools/fhi2upf.f90 index f9e10beb5..06dd3dc8f 100644 --- a/upftools/fhi2upf.f90 +++ b/upftools/fhi2upf.f90 @@ -7,246 +7,246 @@ ! ! !--------------------------------------------------------------------- -program fhi2upf +PROGRAM fhi2upf !--------------------------------------------------------------------- ! ! Convert a pseudopotential file in Fritz-Haber numerical format ! either ".cpi" (fhi88pp) or ".fhi" (abinit) ! to unified pseudopotential format ! May or may not work: carefully check what you get - ! Adapted from the converter written by Andrea Ferretti + ! Adapted from the converter written by Andrea Ferretti ! - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open (unit = 1, file = filein, status = 'old', form = 'formatted') - call read_fhi(1) - close (1) + CALL get_file ( filein ) + OPEN (unit = 1, file = filein, status = 'old', form = 'formatted') + CALL read_fhi(1) + CLOSE (1) ! convert variables read from FHI format into those needed ! by the upf format - add missing quantities - call convert_fhi + CALL convert_fhi fileout=trim(filein)//'.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) -stop -20 write (6,'("fhi2upf: error reading pseudopotential file name")') - stop -end program fhi2upf +STOP +20 WRITE (6,'("fhi2upf: error reading pseudopotential file name")') + STOP +END PROGRAM fhi2upf -module fhi +MODULE fhi ! ! All variables read from FHI file format ! - type angular_comp - real(8), pointer :: pot(:) - real(8), pointer :: wfc(:) - real(8), pointer :: grid(:) + TYPE angular_comp + real(8), POINTER :: pot(:) + real(8), POINTER :: wfc(:) + real(8), POINTER :: grid(:) real(8) :: amesh - integer :: nmesh - integer :: lcomp - end type angular_comp + INTEGER :: nmesh + INTEGER :: lcomp + END TYPE angular_comp !------------------------------ real(8) :: Zval ! valence charge - integer :: lmax_ ! max l-component used + INTEGER :: lmax_ ! max l-component used - logical :: nlcc_ - real(8), allocatable :: rho_atc_(:) ! core charge + LOGICAL :: nlcc_ + real(8), ALLOCATABLE :: rho_atc_(:) ! core charge - type (angular_comp), pointer :: comp(:) ! PP numerical info + TYPE (angular_comp), POINTER :: comp(:) ! PP numerical info ! (wfc, grid, potentials...) !------------------------------ ! variables for the abinit header real(8) :: Zatom, Zion, r2well, rchrg, fchrg, qchrg - integer :: pspdat = 0, pspcod = 0 , pspxc = 0, lloc_ = -1, mmax = 0 - character(len=256) :: info + INTEGER :: pspdat = 0, pspcod = 0 , pspxc = 0, lloc_ = -1, mmax = 0 + CHARACTER(len=256) :: info -end module fhi -! +END MODULE fhi +! ! ---------------------------------------------------------- -subroutine read_fhi(iunps) +SUBROUTINE read_fhi(iunps) ! ---------------------------------------------------------- - ! - use fhi - implicit none - integer, parameter :: Nl=7 ! max number of l-components - integer :: iunps + ! + USE fhi + IMPLICIT NONE + INTEGER, PARAMETER :: Nl=7 ! max number of l-components + INTEGER :: iunps real(8) :: r, rhoc, drhoc, d2rhoc ! - - integer :: l, i, idum, mesh + + INTEGER :: l, i, idum, mesh ! Start reading file - read(iunps,'(a)') info - read(info,*,iostat=i) Zval, l - if ( i /= 0 .or. zval <= 0.0 .or. zval > 100.0 ) then - write (6,'("read_fhi: assuming abinit format")') - read(iunps,*) Zatom, Zion, pspdat - read(iunps,*) pspcod, pspxc, lmax_,lloc_, mmax, r2well - if (pspcod /= 6) then - write (6,'("read_fhi: unknown PP type ",i1,"...stopping")') pspcod - stop - end if - read(iunps,*) rchrg, fchrg, qchrg + READ(iunps,'(a)') info + READ(info,*,iostat=i) Zval, l + IF ( i /= 0 .or. zval <= 0.0 .or. zval > 100.0 ) THEN + WRITE (6,'("read_fhi: assuming abinit format")') + READ(iunps,*) Zatom, Zion, pspdat + READ(iunps,*) pspcod, pspxc, lmax_,lloc_, mmax, r2well + IF (pspcod /= 6) THEN + WRITE (6,'("read_fhi: unknown PP type ",i1,"...stopping")') pspcod + STOP + ENDIF + READ(iunps,*) rchrg, fchrg, qchrg ! - read(iunps,*) - read(iunps,*) - read(iunps,*) + READ(iunps,*) + READ(iunps,*) + READ(iunps,*) ! - read(iunps,*) Zval, l - if (abs(Zion-Zval) > 1.0d-8) then - write (6,'("read_fhi: Zval/Zion mismatch...stopping")') - stop - end if - if (l-1 /= lmax_) then - write (6,'("read_fhi: lmax mismatch...stopping")') - stop - end if - else + READ(iunps,*) Zval, l + IF (abs(Zion-Zval) > 1.0d-8) THEN + WRITE (6,'("read_fhi: Zval/Zion mismatch...stopping")') + STOP + ENDIF + IF (l-1 /= lmax_) THEN + WRITE (6,'("read_fhi: lmax mismatch...stopping")') + STOP + ENDIF + ELSE info = ' ' - end if + ENDIF lmax_ = l - 1 - if (lmax_+1 > Nl) then - write (6,'("read_fhi: too many l-components...stopping")') - stop - end if + IF (lmax_+1 > Nl) THEN + WRITE (6,'("read_fhi: too many l-components...stopping")') + STOP + ENDIF - do i=1,10 - read(iunps,*) ! skipping 11 lines - end do + DO i=1,10 + READ(iunps,*) ! skipping 11 lines + ENDDO - allocate( comp(0:lmax_) ) + ALLOCATE( comp(0:lmax_) ) - do l=0,lmax_ + DO l=0,lmax_ comp(l)%lcomp = l - read(iunps,*) comp(l)%nmesh, comp(l)%amesh - if (mmax > 0 .and. mmax /= comp(l)%nmesh) then - write (6,'("read_fhi: mismatched number of grid points...stopping")') - stop - end if - if ( l > 0) then - if (comp(l)%nmesh /= comp(0)%nmesh .or. & - comp(l)%amesh /= comp(0)%amesh ) then - write(6,'("read_fhi: different radial grids not allowed...stopping")') - stop - end if - end if + READ(iunps,*) comp(l)%nmesh, comp(l)%amesh + IF (mmax > 0 .and. mmax /= comp(l)%nmesh) THEN + WRITE (6,'("read_fhi: mismatched number of grid points...stopping")') + STOP + ENDIF + IF ( l > 0) THEN + IF (comp(l)%nmesh /= comp(0)%nmesh .or. & + comp(l)%amesh /= comp(0)%amesh ) THEN + WRITE(6,'("read_fhi: different radial grids not allowed...stopping")') + STOP + ENDIF + ENDIF mesh = comp(l)%nmesh - allocate( comp(l)%wfc(mesh), & ! wave-functions + ALLOCATE( comp(l)%wfc(mesh), & ! wave-functions comp(l)%pot(mesh), & ! potentials comp(l)%grid(mesh) ) ! real space radial grid ! read the above quantities - do i=1,mesh - read(iunps,*) idum, comp(l)%grid(i), & + DO i=1,mesh + READ(iunps,*) idum, comp(l)%grid(i), & comp(l)%wfc(i), & - comp(l)%pot(i) - end do - end do + comp(l)%pot(i) + ENDDO + ENDDO nlcc_ =.false. - allocate(rho_atc_(comp(0)%nmesh)) + ALLOCATE(rho_atc_(comp(0)%nmesh)) mesh = comp(0)%nmesh - do i=1,mesh - read(iunps,*,end=10, err=20) r, rho_atc_(i), drhoc, d2rhoc - if ( abs( r - comp(0)%grid(i) ) > 1.d-6 ) then - write(6,'("read_fhi: radial grid for core charge? stopping")') - stop - end if - end do + DO i=1,mesh + READ(iunps,*,end=10, err=20) r, rho_atc_(i), drhoc, d2rhoc + IF ( abs( r - comp(0)%grid(i) ) > 1.d-6 ) THEN + WRITE(6,'("read_fhi: radial grid for core charge? stopping")') + STOP + ENDIF + ENDDO nlcc_ = .true. ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential with NLCC successfully read' + WRITE (6,'(a)') 'Pseudopotential with NLCC successfully read' ! ---------------------------------------------------------- - return -10 continue + RETURN +10 CONTINUE ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential without NLCC successfully read' + WRITE (6,'(a)') 'Pseudopotential without NLCC successfully read' ! ---------------------------------------------------------- - return + RETURN ! -20 write(6,'("read_fhi: error reading core charge")') - stop +20 WRITE(6,'("read_fhi: error reading core charge")') + STOP ! -100 write(6,'("read_fhi: error reading pseudopotential file")') - stop +100 WRITE(6,'("read_fhi: error reading pseudopotential file")') + STOP -end subroutine read_fhi +END SUBROUTINE read_fhi ! ---------------------------------------------------------- -subroutine convert_fhi +SUBROUTINE convert_fhi ! ---------------------------------------------------------- ! - use fhi - use upf - use funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc - use constants, ONLY : fpi - implicit none - real(8), parameter :: rmax = 10.0d0 - real(8), allocatable :: aux(:) + USE fhi + USE upf + USE funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc + USE constants, ONLY : fpi + IMPLICIT NONE + real(8), PARAMETER :: rmax = 10.0d0 + real(8), ALLOCATABLE :: aux(:) real(8) :: vll - character (len=20):: dft - character (len=2), external:: atom_name - integer :: lloc, kkbeta - integer :: l, i, ir, iv + CHARACTER (len=20):: dft + CHARACTER (len=2), EXTERNAL:: atom_name + INTEGER :: lloc, kkbeta + INTEGER :: l, i, ir, iv ! - if (nint(Zatom) > 0) then + IF (nint(Zatom) > 0) THEN psd = atom_name(nint(Zatom)) - else - print '("Atom name > ",$)' - read (5,'(a)') psd - end if - if ( lloc_ < 0 ) then - print '("l local (max: ",i1,") > ",$)', lmax_ - read (5,*) lloc - else + ELSE + PRINT '("Atom name > ",$)' + READ (5,'(a)') psd + ENDIF + IF ( lloc_ < 0 ) THEN + PRINT '("l local (max: ",i1,") > ",$)', lmax_ + READ (5,*) lloc + ELSE lloc = lloc_ - end if - if (pspxc == 7) then + ENDIF + IF (pspxc == 7) THEN dft = 'PW' - else - if (pspxc > 0) then - print '("DFT read from abinit file: ",i1)', pspxc - end if - print '("DFT > ",$)' - read (5,'(a)') dft - end if - write(generated, '("Generated using Fritz-Haber code")') - write(date_author,'("Author: unknown Generation date: as well")') - if (trim(info) /= ' ') then + ELSE + IF (pspxc > 0) THEN + PRINT '("DFT read from abinit file: ",i1)', pspxc + ENDIF + PRINT '("DFT > ",$)' + READ (5,'(a)') dft + ENDIF + WRITE(generated, '("Generated using Fritz-Haber code")') + WRITE(date_author,'("Author: unknown Generation date: as well")') + IF (trim(info) /= ' ') THEN comment = trim(info) - else + ELSE comment = 'Info: automatically converted from FHI format' - end if + ENDIF ! reasonable assumption rel = 1 rcloc = 0.0d0 nwfs = lmax_+1 - allocate( els(nwfs), oc(nwfs), epseu(nwfs)) - allocate(lchi(nwfs), nns(nwfs) ) - allocate(rcut (nwfs), rcutus (nwfs)) - do i=1, nwfs - print '("Wavefunction # ",i1,": label, occupancy > ",$)', i - read (5,*) els(i), oc(i) + ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs)) + ALLOCATE(lchi(nwfs), nns(nwfs) ) + ALLOCATE(rcut (nwfs), rcutus (nwfs)) + DO i=1, nwfs + PRINT '("Wavefunction # ",i1,": label, occupancy > ",$)', i + READ (5,*) els(i), oc(i) nns (i) = 0 lchi(i) = i-1 rcut(i) = 0.0d0 rcutus(i)= 0.0d0 epseu(i) = 0.0d0 - end do + ENDDO pseudotype = 'NC' nlcc = nlcc_ @@ -254,91 +254,91 @@ subroutine convert_fhi etotps = 0.0d0 ecutrho=0.0d0 ecutwfc=0.0d0 - if ( lmax_ == lloc) then + IF ( lmax_ == lloc) THEN lmax = lmax_-1 - else + ELSE lmax = lmax_ - end if + ENDIF nbeta= lmax_ mesh = comp(0)%nmesh ntwfc= nwfs - allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) - do i=1, nwfs + ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) + DO i=1, nwfs lchiw(i) = lchi(i) ocw(i) = oc(i) elsw(i) = els(i) - end do - call set_dft_from_name(dft) + ENDDO + CALL set_dft_from_name(dft) iexch = get_iexch() icorr = get_icorr() igcx = get_igcx() igcc = get_igcc() - allocate(rab(mesh)) - allocate( r(mesh)) + ALLOCATE(rab(mesh)) + ALLOCATE( r(mesh)) r = comp(0)%grid rab = r * log( comp(0)%amesh ) - if (nlcc) then - allocate (rho_atc(mesh)) + IF (nlcc) THEN + ALLOCATE (rho_atc(mesh)) rho_atc(:) = rho_atc_(:) / fpi - end if + ENDIF - allocate (vloc0(mesh)) + ALLOCATE (vloc0(mesh)) ! the factor 2 converts from Hartree to Rydberg vloc0 = 2.d0*comp(lloc)%pot - if (nbeta > 0) then + IF (nbeta > 0) THEN - allocate(ikk2(nbeta), lll(nbeta)) + ALLOCATE(ikk2(nbeta), lll(nbeta)) kkbeta=mesh - do ir = 1,mesh - if ( r(ir) > rmax ) then + DO ir = 1,mesh + IF ( r(ir) > rmax ) THEN kkbeta=ir exit - end if - end do + ENDIF + ENDDO ikk2(:) = kkbeta - allocate(aux(kkbeta)) - allocate(betar(mesh,nbeta)) - allocate(qfunc(mesh,nbeta,nbeta)) - allocate(dion(nbeta,nbeta)) - allocate(qqq (nbeta,nbeta)) + ALLOCATE(aux(kkbeta)) + ALLOCATE(betar(mesh,nbeta)) + ALLOCATE(qfunc(mesh,nbeta,nbeta)) + ALLOCATE(dion(nbeta,nbeta)) + ALLOCATE(qqq (nbeta,nbeta)) qfunc(:,:,:)=0.0d0 dion(:,:) =0.d0 qqq(:,:) =0.d0 iv=0 - do i=1,nwfs + DO i=1,nwfs l=lchi(i) - if (l.ne.lloc) then + IF (l/=lloc) THEN iv=iv+1 lll(iv)=l - do ir=1,kkbeta - ! FHI potentials are in Hartree + DO ir=1,kkbeta + ! FHI potentials are in Hartree betar(ir,iv) = 2.d0 * comp(l)%wfc(ir) * & ( comp(l)%pot(ir) - comp(lloc)%pot(ir) ) aux(ir) = comp(l)%wfc(ir) * betar(ir,iv) - end do - call simpson(kkbeta,aux,rab,vll) + ENDDO + CALL simpson(kkbeta,aux,rab,vll) dion(iv,iv) = 1.0d0/vll - end if - enddo + ENDIF + ENDDO - end if + ENDIF - allocate (rho_at(mesh)) + ALLOCATE (rho_at(mesh)) rho_at = 0.d0 - do i=1,nwfs + DO i=1,nwfs l=lchi(i) rho_at = rho_at + ocw(i) * comp(l)%wfc ** 2 - end do - - allocate (chi(mesh,ntwfc)) - do i=1,ntwfc + ENDDO + + ALLOCATE (chi(mesh,ntwfc)) + DO i=1,ntwfc chi(:,i) = comp(i-1)%wfc(:) - end do + ENDDO ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully converted' + WRITE (6,'(a)') 'Pseudopotential successfully converted' ! ---------------------------------------------------------- - return -end subroutine convert_fhi + RETURN +END SUBROUTINE convert_fhi diff --git a/upftools/fpmd2upf.f90 b/upftools/fpmd2upf.f90 index 0658c16a3..a64c5f317 100644 --- a/upftools/fpmd2upf.f90 +++ b/upftools/fpmd2upf.f90 @@ -29,20 +29,20 @@ ! Example: -module fpmd2upf_module +MODULE fpmd2upf_module USE kinds, ONLY: DP USE parameters - use radial_grids, ONLY: ndmx + USE radial_grids, ONLY: ndmx - implicit none - save + IMPLICIT NONE + SAVE REAL(DP), PRIVATE :: TOLMESH = 1.d-5 TYPE pseudo_ncpp - CHARACTER(LEN=4) :: psd ! Element label - CHARACTER(LEN=20) :: pottyp ! Potential type + CHARACTER(len=4) :: psd ! Element label + CHARACTER(len=20) :: pottyp ! Potential type LOGICAL :: tmix LOGICAL :: tnlcc INTEGER :: igau @@ -72,53 +72,53 @@ module fpmd2upf_module END TYPE pseudo_ncpp -contains +CONTAINS - subroutine read_pseudo_fpmd( ap, psfile ) - type(pseudo_ncpp) :: ap - character(len=256) :: psfile - character(len=80) :: error_msg - integer :: info, iunit + SUBROUTINE read_pseudo_fpmd( ap, psfile ) + TYPE(pseudo_ncpp) :: ap + CHARACTER(len=256) :: psfile + CHARACTER(len=80) :: error_msg + INTEGER :: info, iunit iunit = 11 OPEN(UNIT=iunit,FILE=psfile,STATUS='OLD') REWIND( iunit ) CALL read_head_pp( iunit, ap, error_msg, info) - IF( info /= 0 ) GO TO 200 - + IF( info /= 0 ) GOTO 200 + IF( ap%pottyp == 'GIANNOZ' ) THEN CALL read_giannoz(iunit, ap, info) - IF( info /= 0 ) GO TO 200 + IF( info /= 0 ) GOTO 200 - ELSE IF( ap%pottyp == 'NUMERIC' ) THEN + ELSEIF( ap%pottyp == 'NUMERIC' ) THEN CALL read_numeric_pp( iunit, ap, error_msg, info) - IF( info /= 0 ) GO TO 200 + IF( info /= 0 ) GOTO 200 - ELSE IF( ap%pottyp == 'ANALYTIC' ) THEN + ELSEIF( ap%pottyp == 'ANALYTIC' ) THEN CALL read_analytic_pp( iunit, ap, error_msg, info) - IF( info /= 0 ) GO TO 200 + IF( info /= 0 ) GOTO 200 ELSE info = 1 - error_msg = ' Pseudopotential type '//TRIM(ap%pottyp)//' not implemented ' - GO TO 200 + error_msg = ' Pseudopotential type '//trim(ap%pottyp)//' not implemented ' + GOTO 200 - END IF + ENDIF 200 CONTINUE IF( info /= 0 ) THEN - CALL errore(' readpseudo ', error_msg, ABS(info) ) - END IF + CALL errore(' readpseudo ', error_msg, abs(info) ) + ENDIF CLOSE(iunit) - return - end subroutine read_pseudo_fpmd + RETURN + END SUBROUTINE read_pseudo_fpmd !=----------------------------------------------------------------------------=! @@ -127,28 +127,28 @@ contains ! ... This sub. check if a given fortran unit 'iunit' contains a UPF pseudopot. - INTEGER, INTENT(IN) :: iunit - INTEGER, INTENT(OUT) :: info - CHARACTER(LEN=80) :: dummy + INTEGER, INTENT(in) :: iunit + INTEGER, INTENT(out) :: info + CHARACTER(len=80) :: dummy INTEGER :: ios LOGICAL, EXTERNAL :: matches info = 0 ios = 0 - header_loop: do while (ios == 0) - read (iunit, *, iostat = ios, err = 200) dummy - if (matches ("", dummy) ) then + header_loop: DO WHILE (ios == 0) + READ (iunit, *, iostat = ios, err = 200) dummy + IF (matches ("", dummy) ) THEN info = 1 exit header_loop - endif - enddo header_loop -200 continue + ENDIF + ENDDO header_loop +200 CONTINUE RETURN END SUBROUTINE check_file_type !=----------------------------------------------------------------------------=! SUBROUTINE analytic_to_numeric(ap) - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap + TYPE (pseudo_ncpp), INTENT(inout) :: ap INTEGER :: ir, mesh, lmax, l, n, il, ib, ll REAL(DP) :: xmin, zmesh, dx, x ! REAL(DP) :: pi = 3.14159265358979323846_DP @@ -158,36 +158,36 @@ contains IF( ap%mesh == 0 ) THEN ! ... Local pseudopotential, define a logaritmic grid - mesh = SIZE( ap%rw ) + mesh = size( ap%rw ) xmin = -5.0d0 zmesh = 6.0d0 dx = 0.025d0 DO ir = 1, mesh - x = xmin + DBLE(ir-1) * dx - ap%rw(ir) = EXP(x) / zmesh - IF( ap%rw(ir) > 1000.0d0 ) EXIT - END DO + x = xmin + dble(ir-1) * dx + ap%rw(ir) = exp(x) / zmesh + IF( ap%rw(ir) > 1000.0d0 ) exit + ENDDO ap%mesh = mesh ap%dx = dx ap%rab = ap%dx * ap%rw - END IF + ENDIF ap%vnl = 0.0d0 ap%vloc = 0.0d0 ap%vrps = 0.0d0 - do l = 1, 3 - do ir = 1, ap%mesh - ap%vnl(ir,l)= - ( ap%wrc(1) * qe_erf(SQRT(ap%rc(1))*ap%rw(ir)) + & - ap%wrc(2) * qe_erf(SQRT(ap%rc(2))*ap%rw(ir)) ) *& + DO l = 1, 3 + DO ir = 1, ap%mesh + ap%vnl(ir,l)= - ( ap%wrc(1) * qe_erf(sqrt(ap%rc(1))*ap%rw(ir)) + & + ap%wrc(2) * qe_erf(sqrt(ap%rc(2))*ap%rw(ir)) ) *& ap%zv / ap%rw(ir) - end do - do ir = 1, ap%mesh - do n = 1, ap%igau + ENDDO + DO ir = 1, ap%mesh + DO n = 1, ap%igau ap%vnl(ir,l)= ap%vnl(ir,l)+ (ap%al(n,l)+ ap%bl(n,l)*ap%rw(ir)**2 )* & - EXP(-ap%rcl(n,l)*ap%rw(ir)**2) - end do - end do - end do + exp(-ap%rcl(n,l)*ap%rw(ir)**2) + ENDDO + ENDDO + ENDDO ! ... Copy local component to a separate array ap%vloc(:) = ap%vnl(:,ap%lloc) @@ -195,7 +195,7 @@ contains ll=ap%lll(l) + 1 ! find out the angular momentum (ll-1) of the component stored ! in position l ap%vrps(:,l) = ( ap%vnl(:,ll) - ap%vloc(:) ) * ap%rps(:,ll) - END DO + ENDDO RETURN END SUBROUTINE analytic_to_numeric @@ -205,24 +205,24 @@ contains SUBROUTINE read_giannoz(uni, ap, ierr) ! USE constants, ONLY : fpi IMPLICIT NONE - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap - INTEGER, INTENT(IN) :: uni - INTEGER, INTENT(OUT) :: ierr - REAL(DP) :: chi( SIZE(ap%rps, 1), SIZE(ap%rps, 2) ) - REAL(DP) :: vnl( SIZE(ap%vnl, 1), SIZE(ap%vnl, 2) ) - REAL(DP) :: rho_core( SIZE(ap%rhoc, 1) ) + TYPE (pseudo_ncpp), INTENT(inout) :: ap + INTEGER, INTENT(in) :: uni + INTEGER, INTENT(out) :: ierr + REAL(DP) :: chi( size(ap%rps, 1), size(ap%rps, 2) ) + REAL(DP) :: vnl( size(ap%vnl, 1), size(ap%vnl, 2) ) + REAL(DP) :: rho_core( size(ap%rhoc, 1) ) REAL(DP) :: r, ra, rb, fac - REAL(DP) :: oc( SIZE(ap%rps, 2) ) - REAL(DP) :: enl( SIZE(ap%rps, 2) ) + REAL(DP) :: oc( size(ap%rps, 2) ) + REAL(DP) :: enl( size(ap%rps, 2) ) REAL(DP) :: zmesh, xmin, dx, etot REAL(DP) :: zval - INTEGER :: nn(SIZE(ap%rps, 2)), ll(SIZE(ap%rps, 2)) + INTEGER :: nn(size(ap%rps, 2)), ll(size(ap%rps, 2)) INTEGER :: nwf, mesh, i, j, in1, in2, in3, in4, m INTEGER :: lmax, nlc, nnl, lloc, l, il LOGICAL :: nlcc CHARACTER(len=80) :: dft CHARACTER(len=4) :: atom - CHARACTER(len=2) :: el( SIZE(ap%rps, 2) ) + CHARACTER(len=2) :: el( size(ap%rps, 2) ) CHARACTER(len=80) :: ppinfo CHARACTER(len=80) :: strdum CHARACTER(len=2) :: sdum1, sdum2 @@ -236,38 +236,38 @@ contains ! WRITE(6,*) ' DEBUG ', atom, zval,lmax, nlc, nnl, nlcc, lloc, ppinfo - IF( (lmax+1) > SIZE(ap%vnl, 2) ) THEN + IF( (lmax+1) > size(ap%vnl, 2) ) THEN ierr = 1 RETURN - END IF - IF( (nlcc .AND. .NOT.ap%tnlcc) .OR. (.NOT.nlcc .AND. ap%tnlcc) ) THEN + ENDIF + IF( (nlcc .and. .not.ap%tnlcc) .or. (.not.nlcc .and. ap%tnlcc) ) THEN ierr = 2 RETURN - END IF + ENDIF READ(uni,fmt='(f8.2,f8.4,f10.6,2i6)') zmesh, xmin, dx, mesh, nwf - IF( mesh > SIZE(ap%rps, 1) ) THEN + IF( mesh > size(ap%rps, 1) ) THEN ierr = 3 RETURN - END IF - IF( nwf > SIZE(ap%rps, 2) ) THEN + ENDIF + IF( nwf > size(ap%rps, 2) ) THEN ierr = 4 RETURN - END IF + ENDIF DO j = 0, lmax READ(uni,fmt="(A16,i1)") strdum, l READ(uni,'(4e16.8)') (vnl(i,j+1), i=1,mesh) - END DO + ENDDO IF (nlcc) THEN READ(uni,fmt='(4e16.8)') (rho_core(i), i=1,mesh) - END IF + ENDIF DO j = 1, nwf READ(uni,fmt="(A16,a2)") strdum,el(j) READ(uni,fmt='(i5,f6.2)') ll(j),oc(j) READ(uni,fmt='(4e16.8)') (chi(i,j), i=1,mesh) - END DO + ENDDO ap%zv = zval ap%nchan = lmax+1 @@ -280,16 +280,16 @@ contains ! WRITE(6,*) ' DEBUG ', ap%lloc, ap%numeric, ap%nbeta, ap%raggio, ap%zv DO i = 1, mesh - r = EXP(xmin+DBLE(i-1)*dx)/zmesh + r = exp(xmin+dble(i-1)*dx)/zmesh ap%rw(i) = r DO j = 1, lmax+1 ap%vnl(i,j) = vnl(i,j) * fac - END DO - END DO - IF( MINVAL( ap%rw(1:mesh) ) <= 0.0d0 ) THEN + ENDDO + ENDDO + IF( minval( ap%rw(1:mesh) ) <= 0.0d0 ) THEN ierr = 5 RETURN - END IF + ENDIF ap%dx = dx ap%rab = ap%dx * ap%rw ap%vloc(:) = ap%vnl(:,ap%lloc) @@ -302,27 +302,27 @@ contains ! fac = 1.0d0/SQRT(fpi) fac = 1.0d0 DO i = 1, mesh - r = EXP(xmin+DBLE(i-1)*dx)/zmesh + r = exp(xmin+dble(i-1)*dx)/zmesh DO j = 1, nwf ap%rps(i,j) = chi(i,j) * fac - END DO - END DO + ENDDO + ENDDO DO l = 1, ap%nbeta il=ap%lll(l) + 1 ! find out the angular momentum (il-1) of the component stored ! in position l DO i = 1, mesh ap%vrps(i,l) = ( ap%vnl(i,il) - ap%vloc(i) ) * ap%rps(i,il) - END DO - END DO + ENDDO + ENDDO IF( nlcc ) THEN ap%rhoc = 0.0d0 DO i = 1, mesh - r = EXP(xmin+DBLE(i-1)*dx)/zmesh + r = exp(xmin+dble(i-1)*dx)/zmesh ap%rhoc(i) = rho_core(i) - END DO - END IF + ENDDO + ENDIF RETURN END SUBROUTINE read_giannoz @@ -331,24 +331,24 @@ contains SUBROUTINE ap_info( ap ) - TYPE (pseudo_ncpp), INTENT(IN) :: ap + TYPE (pseudo_ncpp), INTENT(in) :: ap INTEGER :: in1, in2, in3, in4, m, il, ib, l, i IF (ap%nbeta > 0) THEN WRITE(6,10) ap%pottyp IF (ap%tmix) THEN - WRITE(6,107) + WRITE(6,107) WRITE(6,106) (ap%lll(l),l=1,ap%nbeta) WRITE(6,105) (ap%wgv(l),l=1,ap%nbeta) ELSE WRITE(6,50) ap%lloc - END IF + ENDIF WRITE(6,60) (ap%lll(l),l=1,ap%nbeta) ELSE ! ... A local pseudopotential has been read. WRITE(6,11) ap%pottyp WRITE(6,50) ap%lloc - END IF + ENDIF 10 FORMAT( 3X,'Type is ',A10,' and NONLOCAL. ') 107 FORMAT( 3X,'Mixed reference potential:') @@ -391,15 +391,15 @@ contains DO il=1,3 DO ib=1,ap%igau WRITE(6,103) ap%rcl(ib,il),ap%al(ib,il),ap%bl(ib,il) - END DO - END DO + ENDDO + ENDDO 40 FORMAT( 3X,'Hsc radii and coeff. A and B :') 103 FORMAT(3X,F8.4,2(3X,F15.7)) - END IF + ENDIF - IF( ap%nrps > 0 .AND. ap%mesh > 0 ) THEN + IF( ap%nrps > 0 .and. ap%mesh > 0 ) THEN WRITE(6,141) ap%nrps, ap%mesh, ap%dx in1=1 in2=ap%mesh/4 @@ -411,7 +411,7 @@ contains WRITE(6,120) in2,ap%rw(in2),(ap%rps(in2,m),m=1,ap%nrps) WRITE(6,120) in3,ap%rw(in3),(ap%rps(in3,m),m=1,ap%nrps) WRITE(6,120) in4,ap%rw(in4),(ap%rps(in4,m),m=1,ap%nrps) - END IF + ENDIF 141 FORMAT(/, 3X,'Atomic wavefunction Grid : Channels = ',I2,& ', Mesh = ',I5,/,30X,'dx = ',F16.14) @@ -429,7 +429,7 @@ contains WRITE(6,120) in2,ap%rw(in2),ap%rhoc(in2) WRITE(6,120) in3,ap%rw(in3),ap%rhoc(in3) WRITE(6,120) in4,ap%rw(in4),ap%rhoc(in4) - END IF + ENDIF 151 FORMAT(/, 3X,'Core correction Grid : Mesh = ',I5, & ', dx = ',F16.14) @@ -441,15 +441,15 @@ contains !=----------------------------------------------------------------------------=! REAL(DP) FUNCTION calculate_dx( a, m ) - REAL(DP), INTENT(IN) :: a(:) - INTEGER, INTENT(IN) :: m + REAL(DP), INTENT(in) :: a(:) + INTEGER, INTENT(in) :: m INTEGER :: n - REAL(DP) :: ra, rb - n = MIN( SIZE( a ), m ) + REAL(DP) :: ra, rb + n = min( size( a ), m ) ra = a(1) rb = a(n) - calculate_dx = LOG( rb / ra ) / DBLE( n - 1 ) - write(6,*) 'amesh (dx) = ', calculate_dx + calculate_dx = log( rb / ra ) / dble( n - 1 ) + WRITE(6,*) 'amesh (dx) = ', calculate_dx RETURN END FUNCTION calculate_dx @@ -457,12 +457,12 @@ contains SUBROUTINE read_atomic_wf( iunit, ap, err_msg, ierr) USE parser, ONLY: field_count IMPLICIT NONE - INTEGER, INTENT(IN) :: iunit - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap - CHARACTER(LEN=*) :: err_msg - INTEGER, INTENT(OUT) :: ierr + INTEGER, INTENT(in) :: iunit + TYPE (pseudo_ncpp), INTENT(inout) :: ap + CHARACTER(len=*) :: err_msg + INTEGER, INTENT(out) :: ierr ! - CHARACTER(LEN=80) :: input_line + CHARACTER(len=80) :: input_line INTEGER :: i, j, m, strlen, info, nf, mesh REAL(DP) :: rdum @@ -479,7 +479,7 @@ SUBROUTINE read_atomic_wf( iunit, ap, err_msg, ierr) ! this is for local pseudopotentials IF( ap%nbeta == 0 ) RETURN - + READ(iunit,'(A80)',end=100) input_line CALL field_count(nf, input_line) @@ -488,40 +488,40 @@ SUBROUTINE read_atomic_wf( iunit, ap, err_msg, ierr) IF( nf == 2 ) THEN READ(input_line(1:strlen),*,IOSTAT=ierr) mesh, ap%nrps ELSE - READ(input_line(1:strlen),*,IOSTAT=ierr) mesh, ap%nrps, ( ap%oc(j), j=1, MIN(ap%nrps,SIZE(ap%oc)) ) - END IF - IF( ap%nrps > SIZE(ap%rps,2) ) THEN - ierr = 2 + READ(input_line(1:strlen),*,IOSTAT=ierr) mesh, ap%nrps, ( ap%oc(j), j=1, min(ap%nrps,size(ap%oc)) ) + ENDIF + IF( ap%nrps > size(ap%rps,2) ) THEN + ierr = 2 WRITE( 6, * ) ' nchan = (wf) ', ap%nrps err_msg = ' NCHAN NOT PROGRAMMED ' - GO TO 110 - END IF - IF( mesh > SIZE(ap%rw) .OR. mesh < 0) THEN + GOTO 110 + ENDIF + IF( mesh > size(ap%rw) .or. mesh < 0) THEN ierr = 4 err_msg = ' WAVMESH OUT OF RANGE ' - GO TO 110 - END IF + GOTO 110 + ENDIF DO j = 1, mesh READ(iunit,*,IOSTAT=ierr) rdum, (ap%rps(j,m),m=1,ap%nrps) IF( ap%mesh == 0 ) ap%rw(j) = rdum - IF( ABS(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN + IF( abs(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN ierr = 5 err_msg = ' radial meshes do not match ' - GO TO 110 - END IF - END DO + GOTO 110 + ENDIF + ENDDO IF( ap%mesh == 0 ) THEN ap%mesh = mesh ap%dx = calculate_dx( ap%rw, ap%mesh ) ap%rab = ap%dx * ap%rw - END IF + ENDIF GOTO 110 100 ierr = 1 110 CONTINUE - + RETURN END SUBROUTINE read_atomic_wf @@ -529,12 +529,12 @@ END SUBROUTINE read_atomic_wf SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr) IMPLICIT NONE - INTEGER, INTENT(IN) :: iunit - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap - CHARACTER(LEN=*) :: err_msg - INTEGER, INTENT(OUT) :: ierr + INTEGER, INTENT(in) :: iunit + TYPE (pseudo_ncpp), INTENT(inout) :: ap + CHARACTER(len=*) :: err_msg + INTEGER, INTENT(out) :: ierr ! - CHARACTER(LEN=80) :: input_line + CHARACTER(len=80) :: input_line INTEGER :: i, j, m, strlen, info, nf, l, ll ! ... read numeric atomic pseudopotential @@ -545,22 +545,22 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr) IF(ap%tmix) THEN READ(iunit,*) (ap%wgv(l),l=1,ap%nbeta) - END IF + ENDIF READ(iunit,*,IOSTAT=ierr) ap%zv READ(iunit,*,IOSTAT=ierr) ap%mesh, ap%nchan - IF((ap%nchan > SIZE(ap%vnl,2) ) .OR. (ap%nchan < 1)) THEN + IF((ap%nchan > size(ap%vnl,2) ) .or. (ap%nchan < 1)) THEN ierr = 1 WRITE( 6, * ) ' nchan (pp) = ', ap%nchan err_msg = ' NCHAN NOT PROGRAMMED ' - GO TO 110 - END IF - IF((ap%mesh > SIZE(ap%rw) ) .OR. (ap%mesh < 0)) THEN + GOTO 110 + ENDIF + IF((ap%mesh > size(ap%rw) ) .or. (ap%mesh < 0)) THEN info = 2 err_msg = ' NPOTMESH OUT OF RANGE ' - GO TO 110 - END IF + GOTO 110 + ENDIF ap%rw = 0.0d0 ap%vnl = 0.0d0 @@ -568,13 +568,13 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr) ap%vrps = 0.0d0 DO j = 1, ap%mesh READ(iunit,*,IOSTAT=ierr) ap%rw(j), (ap%vnl(j,l),l=1,ap%nchan) - END DO + ENDDO - IF( MINVAL( ap%rw(1:ap%mesh) ) <= 0.0d0 ) THEN + IF( minval( ap%rw(1:ap%mesh) ) <= 0.0d0 ) THEN info = 30 err_msg = ' ap rw too small ' - GO TO 110 - END IF + GOTO 110 + ENDIF ! ... mixed reference potential is in vr(lloc) IF(ap%tmix) THEN @@ -583,31 +583,31 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr) DO l=1,ap%nchan IF(l /= ap%lloc) THEN ap%vnl(j,ap%lloc)= ap%vnl(j,ap%lloc) + ap%wgv(l) * ap%vnl(j,l) - END IF - END DO - END DO - END IF + ENDIF + ENDDO + ENDDO + ENDIF ap%vloc(:) = ap%vnl(:,ap%lloc) ap%dx = calculate_dx( ap%rw, ap%mesh ) ap%rab = ap%dx * ap%rw CALL read_atomic_wf( iunit, ap, err_msg, ierr) - IF( ierr /= 0 ) GO TO 110 + IF( ierr /= 0 ) GOTO 110 DO l = 1, ap%nbeta - ll=ap%lll(l) + 1 + ll=ap%lll(l) + 1 ap%vrps(:,l) = ( ap%vnl(:,ll) - ap%vloc(:) ) * ap%rps(:,ll) - END DO + ENDDO IF(ap%tnlcc) THEN CALL read_atomic_cc( iunit, ap, err_msg, ierr) - IF( ierr /= 0 ) GO TO 110 - END IF + IF( ierr /= 0 ) GOTO 110 + ENDIF GOTO 110 100 ierr = 1 110 CONTINUE - + RETURN END SUBROUTINE read_numeric_pp @@ -615,10 +615,10 @@ END SUBROUTINE read_numeric_pp SUBROUTINE read_head_pp( iunit, ap, err_msg, ierr) IMPLICIT NONE - INTEGER, INTENT(IN) :: iunit - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap - CHARACTER(LEN=*) :: err_msg - INTEGER, INTENT(OUT) :: ierr + INTEGER, INTENT(in) :: iunit + TYPE (pseudo_ncpp), INTENT(inout) :: ap + CHARACTER(len=*) :: err_msg + INTEGER, INTENT(out) :: ierr ! INTEGER :: i, l @@ -629,44 +629,44 @@ SUBROUTINE read_head_pp( iunit, ap, err_msg, ierr) ap%lll = 0 READ(iunit, *) ap%tnlcc, ap%tmix - READ(iunit, *) ap%pottyp, ap%lloc, ap%nbeta, (ap%lll(l), l = 1, MIN(ap%nbeta, SIZE(ap%lll)) ) + READ(iunit, *) ap%pottyp, ap%lloc, ap%nbeta, (ap%lll(l), l = 1, min(ap%nbeta, size(ap%lll)) ) ap%lll = ap%lll - 1 - IF( ap%nbeta > SIZE(ap%lll) .OR. ap%nbeta < 0 ) THEN + IF( ap%nbeta > size(ap%lll) .or. ap%nbeta < 0 ) THEN ierr = 1 err_msg = 'LNL out of range' - GO TO 110 - END IF - IF( ap%lloc < 0 .OR. ap%lloc > SIZE(ap%vnl,2) ) THEN + GOTO 110 + ENDIF + IF( ap%lloc < 0 .or. ap%lloc > size(ap%vnl,2) ) THEN ierr = 3 err_msg = 'LLOC out of range' - GO TO 110 - END IF - IF( ap%tmix .AND. ap%pottyp /= 'NUMERIC' ) THEN + GOTO 110 + ENDIF + IF( ap%tmix .and. ap%pottyp /= 'NUMERIC' ) THEN ierr = 4 err_msg = 'tmix not implemented for pseudo ' // ap%pottyp - GO TO 110 - END IF + GOTO 110 + ENDIF DO l = 2, ap%nbeta IF( ap%lll(l) <= ap%lll(l-1)) THEN ierr = 5 err_msg =' NONLOCAL COMPONENTS MUST BE GIVEN IN ASCENDING ORDER' - GO TO 110 - END IF - END DO + GOTO 110 + ENDIF + ENDDO DO l = 1, ap%nbeta IF( ap%lll(l)+1 == ap%lloc) THEN ierr = 6 - err_msg = ' LLOC.EQ.L NON LOCAL!!' - GO TO 110 - END IF - END DO + err_msg = ' LLOC.EQ.L NON LOCAL!!' + GOTO 110 + ENDIF + ENDDO GOTO 110 100 ierr = 1 110 CONTINUE - + RETURN END SUBROUTINE read_head_pp @@ -674,10 +674,10 @@ END SUBROUTINE read_head_pp SUBROUTINE read_analytic_pp( iunit, ap, err_msg, ierr) IMPLICIT NONE - INTEGER, INTENT(IN) :: iunit - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap - CHARACTER(LEN=*) :: err_msg - INTEGER, INTENT(OUT) :: ierr + INTEGER, INTENT(in) :: iunit + TYPE (pseudo_ncpp), INTENT(inout) :: ap + CHARACTER(len=*) :: err_msg + INTEGER, INTENT(out) :: ierr ! INTEGER :: i, l @@ -688,8 +688,8 @@ SUBROUTINE read_analytic_pp( iunit, ap, err_msg, ierr) READ(iunit,*,IOSTAT=ierr) ap%zv, ap%igau - ap%mesh = 0 - ap%nchan = 0 + ap%mesh = 0 + ap%nchan = 0 ap%dx = 0.0d0 ap%rab = 0.0d0 ap%rw = 0.0d0 @@ -708,22 +708,22 @@ SUBROUTINE read_analytic_pp( iunit, ap, err_msg, ierr) CASE DEFAULT ierr = 1 err_msg = ' IGAU NOT PROGRAMMED ' - GO TO 110 + GOTO 110 END SELECT DO l=1,3 DO i=1,ap%igau READ(iunit,*,IOSTAT=ierr) ap%rcl(i,l), ap%al(i,l), ap%bl(i,l) - END DO - END DO + ENDDO + ENDDO CALL read_atomic_wf( iunit, ap, err_msg, ierr) - IF( ierr /= 0 ) GO TO 110 + IF( ierr /= 0 ) GOTO 110 IF(ap%tnlcc) THEN CALL read_atomic_cc( iunit, ap, err_msg, ierr) - IF( ierr /= 0 ) GO TO 110 - END IF + IF( ierr /= 0 ) GOTO 110 + ENDIF ! ... Analytic pseudo are not supported anymore, conversion ! ... to numeric form is forced @@ -732,7 +732,7 @@ SUBROUTINE read_analytic_pp( iunit, ap, err_msg, ierr) GOTO 110 100 ierr = 1 110 CONTINUE - + RETURN END SUBROUTINE read_analytic_pp @@ -741,12 +741,12 @@ END SUBROUTINE read_analytic_pp SUBROUTINE read_atomic_cc( iunit, ap, err_msg, ierr) IMPLICIT NONE - INTEGER, INTENT(IN) :: iunit - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap - CHARACTER(LEN=*) :: err_msg - INTEGER, INTENT(OUT) :: ierr + INTEGER, INTENT(in) :: iunit + TYPE (pseudo_ncpp), INTENT(inout) :: ap + CHARACTER(len=*) :: err_msg + INTEGER, INTENT(out) :: ierr ! - CHARACTER(LEN=80) :: input_line + CHARACTER(len=80) :: input_line INTEGER :: j, mesh REAL(DP) :: rdum @@ -758,41 +758,41 @@ SUBROUTINE read_atomic_cc( iunit, ap, err_msg, ierr) ap%rhoc = 0.0d0 READ(iunit,*,IOSTAT=ierr) mesh - IF(mesh > SIZE(ap%rw) .OR. mesh < 0 ) THEN + IF(mesh > size(ap%rw) .or. mesh < 0 ) THEN ierr = 17 err_msg = ' CORE CORRECTION MESH OUT OF RANGE ' - GO TO 110 - END IF + GOTO 110 + ENDIF DO j = 1, mesh READ(iunit,*,IOSTAT=ierr) rdum, ap%rhoc(j) IF( ap%mesh == 0 ) ap%rw(j) = rdum - IF( ABS(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN + IF( abs(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN ierr = 5 err_msg = ' core cor. radial mesh does not match ' - GO TO 110 - END IF - END DO + GOTO 110 + ENDIF + ENDDO IF( ap%mesh == 0 ) THEN ap%mesh = mesh ap%dx = calculate_dx( ap%rw, ap%mesh ) ap%rab = ap%dx * ap%rw - END IF + ENDIF GOTO 110 100 ierr = 1 110 CONTINUE - + RETURN END SUBROUTINE read_atomic_cc -end module fpmd2upf_module +END MODULE fpmd2upf_module -program fpmd2upf +PROGRAM fpmd2upf ! ! Convert a pseudopotential written in the FPMD format @@ -807,39 +807,39 @@ program fpmd2upf IMPLICIT NONE TYPE (pseudo_ncpp) :: ap - CHARACTER(LEN=256) :: psfile - CHARACTER(LEN=2) :: wfl( 10 ) + CHARACTER(len=256) :: psfile + CHARACTER(len=2) :: wfl( 10 ) REAL(8) :: wfoc( 10 ) INTEGER :: nsp, nspnl, i, lloc, l, ir, iv, kkbeta REAL(8) :: rmax = 10 REAL(8) :: vll - REAL(8), allocatable :: aux(:) + REAL(8), ALLOCATABLE :: aux(:) - namelist / fpmd_pseudo / psfile, nwfs, wfl, wfoc, psd, & + NAMELIST / fpmd_pseudo / psfile, nwfs, wfl, wfoc, psd, & iexch, icorr, igcx, igcc, zp ! ... end of declarations - call input_from_file() + CALL input_from_file() - read( 5, fpmd_pseudo ) + READ( 5, fpmd_pseudo ) nsp = 1 CALL read_pseudo_fpmd(ap, psfile) - write(generated, '("Generated using unknown code")') - write(date_author,'("Author: unknown Generation date: as well")') + WRITE(generated, '("Generated using unknown code")') + WRITE(date_author,'("Author: unknown Generation date: as well")') comment = 'Info: automatically converted from CPMD format' rcloc = 0.0d0 - allocate( els(nwfs), oc(nwfs), epseu(nwfs) ) - allocate( lchi(nwfs), nns(nwfs) ) - allocate( rcut (nwfs), rcutus (nwfs) ) + ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs) ) + ALLOCATE( lchi(nwfs), nns(nwfs) ) + ALLOCATE( rcut (nwfs), rcutus (nwfs) ) - els = '?' + els = '?' oc = 0.0d0 - do i = 1, nwfs + DO i = 1, nwfs els(i) = wfl(i) oc(i) = wfoc(i) lchi(i) = i - 1 @@ -847,110 +847,110 @@ program fpmd2upf rcut(i) = 0.0d0 rcutus(i)= 0.0d0 epseu(i) = 0.0d0 - end do + ENDDO pseudotype = 'NC' nlcc = ap%tnlcc - if( ap%zv > 0.0d0 ) zp = ap%zv + IF( ap%zv > 0.0d0 ) zp = ap%zv etotps = 0.0d0 - lloc = ap%lloc - lmax = MAX( MAXVAL( ap%lll( 1:ap%nbeta ) ), ap%lloc - 1 ) + lloc = ap%lloc + lmax = max( maxval( ap%lll( 1:ap%nbeta ) ), ap%lloc - 1 ) nbeta = ap%nbeta mesh = ap%mesh ntwfc = nwfs - allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) - do i = 1, nwfs + ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) + DO i = 1, nwfs lchiw(i) = lchi(i) ocw(i) = oc(i) elsw(i) = els(i) - end do + ENDDO - allocate(rab(mesh)) - allocate( r(mesh)) + ALLOCATE(rab(mesh)) + ALLOCATE( r(mesh)) r = ap%rw ap%dx = calculate_dx( ap%rw, ap%mesh ) rab = ap%rw * ap%dx - write(6,*) ap%lloc, ap%lll( 1:ap%nbeta ) , ap%nbeta, ap%dx + WRITE(6,*) ap%lloc, ap%lll( 1:ap%nbeta ) , ap%nbeta, ap%dx - allocate (rho_atc(mesh)) - if (nlcc) rho_atc = ap%rhoc + ALLOCATE (rho_atc(mesh)) + IF (nlcc) rho_atc = ap%rhoc - allocate (vloc0(mesh)) + ALLOCATE (vloc0(mesh)) ! the factor 2 converts from Hartree to Rydberg vloc0(:) = ap%vloc * 2.0d0 - if (nbeta > 0) then + IF (nbeta > 0) THEN - allocate(ikk2(nbeta), lll(nbeta)) + ALLOCATE(ikk2(nbeta), lll(nbeta)) kkbeta = mesh - do ir = 1,mesh - if ( r(ir) > rmax ) then + DO ir = 1,mesh + IF ( r(ir) > rmax ) THEN kkbeta=ir exit - end if - end do + ENDIF + ENDDO ikk2(:) = kkbeta - allocate(aux(kkbeta)) - allocate(betar(mesh,nbeta)) - allocate(qfunc(mesh,nbeta,nbeta)) - allocate(dion(nbeta,nbeta)) - allocate(qqq (nbeta,nbeta)) + ALLOCATE(aux(kkbeta)) + ALLOCATE(betar(mesh,nbeta)) + ALLOCATE(qfunc(mesh,nbeta,nbeta)) + ALLOCATE(dion(nbeta,nbeta)) + ALLOCATE(qqq (nbeta,nbeta)) qfunc(:,:,:)=0.0d0 dion(:,:) =0.d0 qqq(:,:) =0.d0 iv = 0 - do i = 1, nwfs + DO i = 1, nwfs l = lchi(i) - if ( l .ne. (lloc-1) ) then + IF ( l /= (lloc-1) ) THEN iv = iv + 1 lll( iv ) = l - do ir = 1, kkbeta + DO ir = 1, kkbeta ! the factor 2 converts from Hartree to Rydberg betar(ir, iv) = 2.d0 * ap%vrps( ir, iv ) aux(ir) = ap%rps(ir, (l+1) ) * betar(ir, iv) - end do - call simpson2(kkbeta, aux(1), rab(1), vll) + ENDDO + CALL simpson2(kkbeta, aux(1), rab(1), vll) dion(iv,iv) = 1.0d0/vll - write(6,*) aux(2), rab(2), kkbeta, vll - end if - enddo + WRITE(6,*) aux(2), rab(2), kkbeta, vll + ENDIF + ENDDO - end if + ENDIF - allocate (rho_at(mesh)) + ALLOCATE (rho_at(mesh)) rho_at = 0.d0 - do i = 1, nwfs + DO i = 1, nwfs rho_at(:) = rho_at(:) + ocw(i) * ap%rps(:, i) ** 2 - end do + ENDDO - allocate (chi(mesh,ntwfc)) + ALLOCATE (chi(mesh,ntwfc)) chi = ap%rps ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully converted' + WRITE (6,'(a)') 'Pseudopotential successfully converted' ! ---------------------------------------------------------- - call write_upf( 10 ) + CALL write_upf( 10 ) -100 continue +100 CONTINUE -end program fpmd2upf +END PROGRAM fpmd2upf !---------------------------------------------------------------------- -subroutine simpson2(mesh,func,rab,asum) +SUBROUTINE simpson2(mesh,func,rab,asum) !----------------------------------------------------------------------- ! ! simpson's rule integrator for function stored on the ! radial logarithmic mesh ! - implicit none + IMPLICIT NONE - integer :: i, mesh + INTEGER :: i, mesh real(8) :: rab(mesh), func(mesh), f1, f2, f3, r12, asum ! routine assumes that mesh is an odd number so run check @@ -964,15 +964,15 @@ subroutine simpson2(mesh,func,rab,asum) r12 = 1.0d0 / 12.0d0 f3 = func(1) * rab(1) * r12 - do i = 2,mesh-1,2 + DO i = 2,mesh-1,2 f1 = f3 f2 = func(i) * rab(i) * r12 f3 = func(i+1) * rab(i+1) * r12 asum = asum + 4.0d0*f1 + 16.0d0*f2 + 4.0d0*f3 - enddo + ENDDO - return -end subroutine simpson2 + RETURN +END SUBROUTINE simpson2 diff --git a/upftools/ncpp2upf.f90 b/upftools/ncpp2upf.f90 index a2e40ceda..2c50b2ebd 100644 --- a/upftools/ncpp2upf.f90 +++ b/upftools/ncpp2upf.f90 @@ -7,281 +7,281 @@ ! ! !--------------------------------------------------------------------- -program ncpp2upf +PROGRAM ncpp2upf !--------------------------------------------------------------------- ! ! Convert a pseudopotential written in PWSCF format ! (norm-conserving) to unified pseudopotential format - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open(unit=1,file=filein,status='old',form='formatted') - call read_ncpp(1) - close (unit=1) + CALL get_file ( filein ) + OPEN(unit=1,file=filein,status='old',form='formatted') + CALL read_ncpp(1) + CLOSE (unit=1) ! convert variables read from NCPP format into those needed ! by the upf format - add missing quantities - call convert_ncpp + CALL convert_ncpp fileout=trim(filein)//'.UPF' - print '(''Output PP file in US format : '',a)', fileout + PRINT '(''Output PP file in US format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) - stop -20 call errore ('ncpp2upf', 'Reading pseudo file name ', 1) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) + STOP +20 CALL errore ('ncpp2upf', 'Reading pseudo file name ', 1) -end program ncpp2upf +END PROGRAM ncpp2upf -module ncpp +MODULE ncpp ! ! All variables read from NCPP file format - ! + ! ! trailing underscore means that a variable with the same name ! is used in module 'upf' containing variables to be written ! - character(len=20) :: dft_ - character(len=2) :: psd_ + CHARACTER(len=20) :: dft_ + CHARACTER(len=2) :: psd_ real(8) :: zp_ - integer nlc, nnl, lmax_, lloc, nchi - logical :: numeric, bhstype, nlcc_ + INTEGER nlc, nnl, lmax_, lloc, nchi + LOGICAL :: numeric, bhstype, nlcc_ real(8) :: alpc(2), cc(2), alps(3,0:3), aps(6,0:3) real(8) :: a_nlcc, b_nlcc, alpha_nlcc real(8) :: zmesh, xmin, dx - real(8), allocatable:: r_(:), rab_(:) - integer :: mesh_ + real(8), ALLOCATABLE:: r_(:), rab_(:) + INTEGER :: mesh_ - real(8), allocatable:: vnl(:,:), rho_atc_(:), rho_at_(:) - integer, allocatable:: lchi_(:) - real(8), allocatable:: chi_(:,:), oc_(:) + real(8), ALLOCATABLE:: vnl(:,:), rho_atc_(:), rho_at_(:) + INTEGER, ALLOCATABLE:: lchi_(:) + real(8), ALLOCATABLE:: chi_(:,:), oc_(:) -end module ncpp -! +END MODULE ncpp +! ! ---------------------------------------------------------- -subroutine read_ncpp(iunps) +SUBROUTINE read_ncpp(iunps) ! ---------------------------------------------------------- - ! - use ncpp - use upf , only : els - implicit none - integer :: iunps ! - character(len=1), dimension(0:3) :: convel=(/'S','P','D','F'/) - character(len=2) :: label + USE ncpp + USE upf , ONLY : els + IMPLICIT NONE + INTEGER :: iunps + ! + CHARACTER(len=1), DIMENSION(0:3) :: convel=(/'S','P','D','F'/) + CHARACTER(len=2) :: label real (8) :: x, qe_erf - integer :: l, i, ir, nb, n - character (len=255) line - external qe_erf + INTEGER :: l, i, ir, nb, n + CHARACTER (len=255) line + EXTERNAL qe_erf - read(iunps, '(a)', end=300, err=300 ) dft_ - if (dft_(1:2).eq.'**') dft_ = 'PZ' + READ(iunps, '(a)', end=300, err=300 ) dft_ + IF (dft_(1:2)=='**') dft_ = 'PZ' - read (iunps, *, err=300) psd_, zp_, lmax_, nlc, nnl, nlcc_, & + READ (iunps, *, err=300) psd_, zp_, lmax_, nlc, nnl, nlcc_, & lloc, bhstype - if ( nlc.gt.2 .or. nnl.gt.3) & - call errore( 'read_ncpp','Wrong nlc or nnl',1 ) - if ( nlc* nnl .lt. 0 ) & - call errore( 'read_ncpp','nlc*nnl < 0 ? ',1 ) - if ( zp_.le.0d0 ) & - call errore( 'read_ncpp','Wrong zp ',1 ) - if ( lmax_.gt.3.or.lmax_.lt.0 ) & - call errore( 'read_ncpp','Wrong lmax ',1 ) + IF ( nlc>2 .or. nnl>3) & + CALL errore( 'read_ncpp','Wrong nlc or nnl',1 ) + IF ( nlc* nnl < 0 ) & + CALL errore( 'read_ncpp','nlc*nnl < 0 ? ',1 ) + IF ( zp_<=0d0 ) & + CALL errore( 'read_ncpp','Wrong zp ',1 ) + IF ( lmax_>3.or.lmax_<0 ) & + CALL errore( 'read_ncpp','Wrong lmax ',1 ) - if (lloc.eq.-1000) lloc=lmax_ + IF (lloc==-1000) lloc=lmax_ ! ! In numeric pseudopotentials both nlc and nnl are zero. ! - numeric = nlc.le.0 .and. nnl.le.0 + numeric = nlc<=0 .and. nnl<=0 - if (.not.numeric) then + IF (.not.numeric) THEN ! ! read pseudopotentials in analytic form ! - read(iunps, *, err=300) & + READ(iunps, *, err=300) & ( alpc(i), i=1, 2 ), ( cc(i), i=1,2 ) - if ( abs(cc(1)+cc(2)-1.d0).gt.1.0d-6) & - call errore ('read_ncpp','wrong pseudopotential coefficients',1) - do l = 0, lmax_ - read (iunps, *, err=300) & + IF ( abs(cc(1)+cc(2)-1.d0)>1.0d-6) & + CALL errore ('read_ncpp','wrong pseudopotential coefficients',1) + DO l = 0, lmax_ + READ (iunps, *, err=300) & ( alps(i,l),i=1,3 ), (aps(i,l),i=1,6) - enddo + ENDDO - if (nlcc_) then - read(iunps, *, err=300) & + IF (nlcc_) THEN + READ(iunps, *, err=300) & a_nlcc, b_nlcc, alpha_nlcc - if (alpha_nlcc.le.0.d0) & - call errore('read_ncpp','nlcc but alpha=0',1) - end if + IF (alpha_nlcc<=0.d0) & + CALL errore('read_ncpp','nlcc but alpha=0',1) + ENDIF - if (bhstype) call bachel(alps,aps,1,lmax_) - end if + IF (bhstype) CALL bachel(alps,aps,1,lmax_) + ENDIF - read(iunps, *, err=300) zmesh, xmin, dx, mesh_, nchi + READ(iunps, *, err=300) zmesh, xmin, dx, mesh_, nchi - if ( mesh_.le.0) call errore( 'read_ncpp', 'mesh too small', 1) - if ( (nchi.lt.lmax_ .and. lloc.eq.lmax_).or. & - (nchi.lt.lmax_+1 .and. lloc.ne.lmax_) ) & - call errore( 'read_ncpp', 'wrong no. of wfcts', 1 ) + IF ( mesh_<=0) CALL errore( 'read_ncpp', 'mesh too small', 1) + IF ( (nchilmax_ .or. lchi_(nb)<0) & + CALL errore('read_ncpp','wrong lchi',nb) + IF ( oc_(nb)<0.d0 .or. & + oc_(nb)>2.d0*(2*lchi_(nb)+1)) & + CALL errore('read_ncpp','wrong oc',nb) ! ! parse and check wavefunction label - read(line,'(14x,a2)', err=222, end=222) label - if (label(2:2).ne.convel(lchi_(nb))) goto 222 - do l = 0, lmax_ - if (label(2:2).eq.convel(l)) then + READ(line,'(14x,a2)', err=222, end=222) label + IF (label(2:2)/=convel(lchi_(nb))) GOTO 222 + DO l = 0, lmax_ + IF (label(2:2)==convel(l)) THEN els(nb) = label(1:2) - goto 223 - endif - end do -222 continue - els(nb) = '*'//convel(lchi_(nb)) -223 continue + GOTO 223 + ENDIF + ENDDO +222 CONTINUE + els(nb) = '*'//convel(lchi_(nb)) +223 CONTINUE ! ! finally read the wavefunction - read(iunps, *, err=300) (chi_(ir,nb),ir=1,mesh_) - enddo + READ(iunps, *, err=300) (chi_(ir,nb),ir=1,mesh_) + ENDDO ! ! compute the atomic charges ! - allocate(rho_at_(mesh_)) + ALLOCATE(rho_at_(mesh_)) rho_at_(:)=0.d0 - do nb = 1, nchi - if( oc_(nb).ne.0.d0) & + DO nb = 1, nchi + IF( oc_(nb)/=0.d0) & rho_at_(:) = rho_at_(:) + oc_(nb)*chi_(:,nb)**2 - end do + ENDDO ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully read' + WRITE (6,'(a)') 'Pseudopotential successfully read' ! ---------------------------------------------------------- - return + RETURN -300 call errore('read_ncpp','pseudo file is empty or wrong',1) +300 CALL errore('read_ncpp','pseudo file is empty or wrong',1) -end subroutine read_ncpp +END SUBROUTINE read_ncpp ! ---------------------------------------------------------- -subroutine convert_ncpp +SUBROUTINE convert_ncpp ! ---------------------------------------------------------- - use ncpp - use upf - use funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc - implicit none - real(8), parameter :: rmax = 10.0d0 - real(8), allocatable :: aux(:) + USE ncpp + USE upf + USE funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc + IMPLICIT NONE + real(8), PARAMETER :: rmax = 10.0d0 + real(8), ALLOCATABLE :: aux(:) real(8) :: vll - integer :: kkbeta, l, iv, ir, i + INTEGER :: kkbeta, l, iv, ir, i - write(generated, '("Generated using ld1 code (maybe, or maybe not)")') - write(date_author,'("Author: unknown Generation date: as well")') + WRITE(generated, '("Generated using ld1 code (maybe, or maybe not)")') + WRITE(date_author,'("Author: unknown Generation date: as well")') comment = 'Info: automatically converted from PWSCF format' ! reasonable assumption - if (zmesh > 18) then + IF (zmesh > 18) THEN rel = 1 - else + ELSE rel = 0 - end if + ENDIF rcloc = 0.0d0 - nwfs = nchi - allocate( oc(nwfs), epseu(nwfs)) - allocate(lchi(nwfs), nns(nwfs) ) - allocate(rcut (nwfs), rcutus (nwfs)) - do i=1, nwfs + nwfs = nchi + ALLOCATE( oc(nwfs), epseu(nwfs)) + ALLOCATE(lchi(nwfs), nns(nwfs) ) + ALLOCATE(rcut (nwfs), rcutus (nwfs)) + DO i=1, nwfs nns (i) = 0 lchi(i) = lchi_(i) rcut(i) = 0.0d0 rcutus(i)= 0.0d0 oc (i) = oc_(i) epseu(i) = 0.0d0 - end do - deallocate (lchi_, oc_) + ENDDO + DEALLOCATE (lchi_, oc_) psd = psd_ pseudotype = 'NC' @@ -290,96 +290,96 @@ subroutine convert_ncpp etotps = 0.0d0 ecutrho=0.0d0 ecutwfc=0.0d0 - if ( lmax_ == lloc) then + IF ( lmax_ == lloc) THEN lmax = lmax_-1 - else + ELSE lmax = lmax_ - end if + ENDIF nbeta= lmax_ mesh = mesh_ ntwfc= nchi - allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) - do i=1, nchi + ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) + DO i=1, nchi lchiw(i) = lchi(i) ocw(i) = oc(i) elsw(i) = els(i) - end do - call set_dft_from_name(dft_) + ENDDO + CALL set_dft_from_name(dft_) iexch = get_iexch() icorr = get_icorr() igcx = get_igcx() igcc = get_igcc() - allocate(rab(mesh)) - allocate( r(mesh)) + ALLOCATE(rab(mesh)) + ALLOCATE( r(mesh)) rab = rab_ r = r_ - allocate (rho_atc(mesh)) + ALLOCATE (rho_atc(mesh)) rho_atc = rho_atc_ - deallocate (rho_atc_) + DEALLOCATE (rho_atc_) - allocate (vloc0(mesh)) + ALLOCATE (vloc0(mesh)) vloc0(:) = vnl(:,lloc) - if (nbeta > 0) then + IF (nbeta > 0) THEN - allocate(ikk2(nbeta), lll(nbeta)) + ALLOCATE(ikk2(nbeta), lll(nbeta)) kkbeta=mesh - do ir = 1,mesh - if ( r(ir) > rmax ) then + DO ir = 1,mesh + IF ( r(ir) > rmax ) THEN kkbeta=ir exit - end if - end do + ENDIF + ENDDO ! make sure kkbeta is odd as required for simpson - if(mod(kkbeta,2) == 0) kkbeta=kkbeta-1 + IF(mod(kkbeta,2) == 0) kkbeta=kkbeta-1 ikk2(:) = kkbeta - allocate(aux(kkbeta)) - allocate(betar(mesh,nbeta)) - allocate(qfunc(mesh,nbeta,nbeta)) - allocate(dion(nbeta,nbeta)) - allocate(qqq (nbeta,nbeta)) + ALLOCATE(aux(kkbeta)) + ALLOCATE(betar(mesh,nbeta)) + ALLOCATE(qfunc(mesh,nbeta,nbeta)) + ALLOCATE(dion(nbeta,nbeta)) + ALLOCATE(qqq (nbeta,nbeta)) qfunc(:,:,:)=0.0d0 dion(:,:) =0.d0 qqq(:,:) =0.d0 iv=0 - do i=1,nchi + DO i=1,nchi l=lchi(i) - if (l.ne.lloc) then + IF (l/=lloc) THEN iv=iv+1 lll(iv)=l - do ir=1,kkbeta + DO ir=1,kkbeta betar(ir,iv)=chi_(ir,i)*vnl(ir,l) aux(ir) = chi_(ir,i)**2*vnl(ir,l) - end do - call simpson(kkbeta,aux,rab,vll) + ENDDO + CALL simpson(kkbeta,aux,rab,vll) dion(iv,iv) = 1.0d0/vll - end if - if(iv >= nbeta) exit ! skip additional pseudo wfns - enddo - deallocate (vnl, aux) + ENDIF + IF(iv >= nbeta) exit ! skip additional pseudo wfns + ENDDO + DEALLOCATE (vnl, aux) ! ! redetermine ikk2 ! - do iv=1,nbeta + DO iv=1,nbeta ikk2(iv)=kkbeta - do ir = kkbeta,1,-1 - if ( abs(betar(ir,iv)) > 1.d-12 ) then + DO ir = kkbeta,1,-1 + IF ( abs(betar(ir,iv)) > 1.d-12 ) THEN ikk2(iv)=ir exit - end if - end do - end do - end if - allocate (rho_at(mesh)) + ENDIF + ENDDO + ENDDO + ENDIF + ALLOCATE (rho_at(mesh)) rho_at = rho_at_ - deallocate (rho_at_) - - allocate (chi(mesh,ntwfc)) - chi = chi_ - deallocate (chi_) + DEALLOCATE (rho_at_) - return -end subroutine convert_ncpp + ALLOCATE (chi(mesh,ntwfc)) + chi = chi_ + DEALLOCATE (chi_) + + RETURN +END SUBROUTINE convert_ncpp diff --git a/upftools/oldcp2upf.f90 b/upftools/oldcp2upf.f90 index 391190e2f..0b8809632 100644 --- a/upftools/oldcp2upf.f90 +++ b/upftools/oldcp2upf.f90 @@ -7,139 +7,139 @@ ! ! !--------------------------------------------------------------------- -program oldcp2upf +PROGRAM oldcp2upf !--------------------------------------------------------------------- ! ! Convert a pseudopotential written in the old CP90 format ! (without core correction) to unified pseudopotential format ! - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open (unit = 1, file = filein, status = 'old', form = 'formatted') - call read_oldcp(1) - close (1) + CALL get_file ( filein ) + OPEN (unit = 1, file = filein, status = 'old', form = 'formatted') + CALL read_oldcp(1) + CLOSE (1) ! convert variables read from old CP90 format into those needed ! by the upf format - add missing quantities - call convert_oldcp + CALL convert_oldcp fileout=trim(filein)//'.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) -stop -20 call errore ('oldcp2upf', 'Reading pseudo file name ', 1) +STOP +20 CALL errore ('oldcp2upf', 'Reading pseudo file name ', 1) -end program oldcp2upf +END PROGRAM oldcp2upf -module oldcp +MODULE oldcp ! ! All variables read from old CP90 file format ! real(8) :: amesh, z, zv - integer :: exfact, lloc, nbeta_, mesh_ + INTEGER :: exfact, lloc, nbeta_, mesh_ real(8) :: wrc1, rc1, wrc2, rc2, rcl(3,3), al(3,3), bl(3,3) - real(8), allocatable :: r_(:), vnl(:,:), chi_(:,:) + real(8), ALLOCATABLE :: r_(:), vnl(:,:), chi_(:,:) ! !------------------------------ -end module oldcp -! +END MODULE oldcp +! ! ---------------------------------------------------------- -subroutine read_oldcp(iunps) +SUBROUTINE read_oldcp(iunps) ! ---------------------------------------------------------- - ! - use oldcp - implicit none - integer :: iunps ! - real(8), external :: qe_erf - integer :: i, l, j, jj + USE oldcp + IMPLICIT NONE + INTEGER :: iunps ! - read(iunps,*, end=10, err=10) z, zv, nbeta_, lloc, exfact - if (z < 1 .or. z > 100 .or. zv < 1 .or. zv > 25 ) & - call errore ('read_oldcp','wrong potential read',1) - read(iunps,*, end=10, err=10) wrc1, rc1, wrc2, rc2 - read(iunps,*, end=10, err=10) ( ( rcl(i,l), al(i,l), & + real(8), EXTERNAL :: qe_erf + INTEGER :: i, l, j, jj + ! + READ(iunps,*, end=10, err=10) z, zv, nbeta_, lloc, exfact + IF (z < 1 .or. z > 100 .or. zv < 1 .or. zv > 25 ) & + CALL errore ('read_oldcp','wrong potential read',1) + READ(iunps,*, end=10, err=10) wrc1, rc1, wrc2, rc2 + READ(iunps,*, end=10, err=10) ( ( rcl(i,l), al(i,l), & bl(i,l), i = 1, 3), l = 1, 3) - read(iunps,*, end=10, err=10) mesh_, amesh - allocate(r_(mesh_)) - allocate (chi_(mesh_,nbeta_)) - do l = 1, nbeta_ - if (l > 1) read(iunps,*, end=10, err=10) mesh_, amesh - do j = 1, mesh_ - read(iunps,*, end=10, err=10) jj, r_(j), chi_(j,l) - end do - end do + READ(iunps,*, end=10, err=10) mesh_, amesh + ALLOCATE(r_(mesh_)) + ALLOCATE (chi_(mesh_,nbeta_)) + DO l = 1, nbeta_ + IF (l > 1) READ(iunps,*, end=10, err=10) mesh_, amesh + DO j = 1, mesh_ + READ(iunps,*, end=10, err=10) jj, r_(j), chi_(j,l) + ENDDO + ENDDO ! ! convert analytic to numeric form ! - allocate (vnl(mesh_,0:nbeta_)) - do l=0,nbeta_ + ALLOCATE (vnl(mesh_,0:nbeta_)) + DO l=0,nbeta_ ! ! DO NOT USE f90 ARRAY SYNTAX: qe_erf IS NOT AN INTRINSIC FUNCTION!!! ! - do j=1, mesh_ + DO j=1, mesh_ vnl(j,l)= - (wrc1*qe_erf(sqrt(rc1)*r_(j)) + & wrc2*qe_erf(sqrt(rc2)*r_(j)) ) * zv/r_(j) - end do + ENDDO ! - do i=1,3 + DO i=1,3 vnl(:,l)= vnl(:,l)+ (al(i,l+1)+ bl(i,l+1)*r_(:)**2) * & exp(-rcl(i,l+1)*r_(:)**2) - end do - end do + ENDDO + ENDDO - return -10 call errore('read_oldcp','error in reading file',1) + RETURN +10 CALL errore('read_oldcp','error in reading file',1) -end subroutine read_oldcp +END SUBROUTINE read_oldcp ! ---------------------------------------------------------- -subroutine convert_oldcp +SUBROUTINE convert_oldcp ! ---------------------------------------------------------- ! - use oldcp - use upf - implicit none - real(8), parameter :: rmax = 10.0d0 - real(8), allocatable :: aux(:) + USE oldcp + USE upf + IMPLICIT NONE + real(8), PARAMETER :: rmax = 10.0d0 + real(8), ALLOCATABLE :: aux(:) real(8) :: vll - character (len=20):: dft - character (len=2), external :: atom_name - integer :: kkbeta - integer :: l, i, ir, iv + CHARACTER (len=20):: dft + CHARACTER (len=2), EXTERNAL :: atom_name + INTEGER :: kkbeta + INTEGER :: l, i, ir, iv ! - write(generated, '("Generated using unknown code")') - write(date_author,'("Author: unknown Generation date: as well")') + WRITE(generated, '("Generated using unknown code")') + WRITE(date_author,'("Author: unknown Generation date: as well")') comment = 'Info: automatically converted from old CP90 format' ! reasonable assumption - if (z > 18) then + IF (z > 18) THEN rel = 1 - else + ELSE rel = 0 - end if + ENDIF rcloc = 0.0d0 nwfs = nbeta_ - allocate( els(nwfs), oc(nwfs), epseu(nwfs)) - allocate(lchi(nwfs), nns(nwfs) ) - allocate(rcut (nwfs), rcutus (nwfs)) - do i=1, nwfs - print '("Wavefunction # ",i1,": label, occupancy > ",$)', i - read (5,*) els(i), oc(i) + ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs)) + ALLOCATE(lchi(nwfs), nns(nwfs) ) + ALLOCATE(rcut (nwfs), rcutus (nwfs)) + DO i=1, nwfs + PRINT '("Wavefunction # ",i1,": label, occupancy > ",$)', i + READ (5,*) els(i), oc(i) nns (i) = 0 lchi(i) = i-1 rcut(i) = 0.0d0 rcutus(i)= 0.0d0 epseu(i) = 0.0d0 - end do + ENDDO psd = atom_name (nint(z)) pseudotype = 'NC' nlcc = .false. @@ -151,90 +151,90 @@ subroutine convert_oldcp nbeta = nbeta_ mesh = mesh_ ntwfc = nwfs - allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) - do i=1, nwfs + ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) + DO i=1, nwfs lchiw(i) = lchi(i) ocw(i) = oc(i) elsw(i) = els(i) - end do + ENDDO ! - if ( exfact.eq.0) then + IF ( exfact==0) THEN iexch=1; icorr=1; igcx=0; igcc=0 ! Perdew-Zunger - else if ( exfact.eq.1) then + ELSEIF ( exfact==1) THEN iexch=1; icorr=3; igcx=1; igcc=3 ! Becke-Lee-Yang-Parr - else if ( exfact.eq.2) then + ELSEIF ( exfact==2) THEN iexch=1; icorr=1; igcx=1; igcc=0 ! Becke88 exchange - else if (exfact.eq.-5.or.exfact.eq.3) then + ELSEIF (exfact==-5.or.exfact==3) THEN iexch=1; icorr=1; igcx=1; igcc=1 ! Becke88-Perdew 86 - else if (exfact.eq.-6.or.exfact.eq.4) then + ELSEIF (exfact==-6.or.exfact==4) THEN iexch=1; icorr=4; igcx=2; igcc=2 ! Perdew-Wang 91 - else if (exfact.eq. 5) then + ELSEIF (exfact== 5) THEN iexch=1; icorr=4; igcx=3; igcc=4 ! Perdew-Becke-Erkerhof - else - call errore('convert','Wrong xc in pseudopotential',1) - end if + ELSE + CALL errore('convert','Wrong xc in pseudopotential',1) + ENDIF - allocate(rab(mesh)) - allocate( r(mesh)) + ALLOCATE(rab(mesh)) + ALLOCATE( r(mesh)) r = r_ rab = r * log( amesh ) ! ! convert analytic to numeric form ! ! - allocate (vloc0(mesh)) + ALLOCATE (vloc0(mesh)) ! the factor 2 converts from Hartree to Rydberg vloc0(:) = vnl(:,lloc)*2.d0 - if (nbeta > 0) then + IF (nbeta > 0) THEN - allocate(ikk2(nbeta), lll(nbeta)) + ALLOCATE(ikk2(nbeta), lll(nbeta)) kkbeta=mesh - do ir = 1,mesh - if ( r(ir) > rmax ) then + DO ir = 1,mesh + IF ( r(ir) > rmax ) THEN kkbeta=ir exit - end if - end do + ENDIF + ENDDO ikk2(:) = kkbeta - allocate(aux(kkbeta)) - allocate(betar(mesh,nbeta)) - allocate(qfunc(mesh,nbeta,nbeta)) - allocate(dion(nbeta,nbeta)) - allocate(qqq (nbeta,nbeta)) + ALLOCATE(aux(kkbeta)) + ALLOCATE(betar(mesh,nbeta)) + ALLOCATE(qfunc(mesh,nbeta,nbeta)) + ALLOCATE(dion(nbeta,nbeta)) + ALLOCATE(qqq (nbeta,nbeta)) qfunc(:,:,:)=0.0d0 dion(:,:) =0.d0 qqq(:,:) =0.d0 iv=0 - do i=1,nwfs + DO i=1,nwfs l=lchi(i) - if (l.ne.lloc) then + IF (l/=lloc) THEN iv=iv+1 lll(iv)=l - do ir=1,kkbeta + DO ir=1,kkbeta ! the factor 2 converts from Hartree to Rydberg betar(ir,iv) = 2.d0 * chi_(ir,l+1) * & ( vnl(ir,l) - vnl(ir,lloc) ) aux(ir) = chi_(ir,l+1) * betar(ir,iv) - end do - call simpson(kkbeta,aux,rab,vll) + ENDDO + CALL simpson(kkbeta,aux,rab,vll) dion(iv,iv) = 1.0d0/vll - end if - enddo + ENDIF + ENDDO - end if + ENDIF - allocate (rho_at(mesh)) + ALLOCATE (rho_at(mesh)) rho_at = 0.d0 - do i=1,nwfs + DO i=1,nwfs rho_at(:) = rho_at(:) + ocw(i) * chi_(:,i) ** 2 - end do - - allocate (chi(mesh,ntwfc)) + ENDDO + + ALLOCATE (chi(mesh,ntwfc)) chi = chi_ ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully converted' + WRITE (6,'(a)') 'Pseudopotential successfully converted' ! ---------------------------------------------------------- - return -end subroutine convert_oldcp + RETURN +END SUBROUTINE convert_oldcp diff --git a/upftools/read_upf.f90 b/upftools/read_upf.f90 index fc8c78295..4819c8bf6 100644 --- a/upftools/read_upf.f90 +++ b/upftools/read_upf.f90 @@ -5,32 +5,32 @@ ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! -module pseudo +MODULE pseudo ! ! All variables to be read from the UPF file ! (UPF = unified pseudopotential format) ! - integer ,parameter :: npsx = 6 + INTEGER ,PARAMETER :: npsx = 6 ! npsx : maximum number of different pseudopotentials - integer, parameter :: lmaxx = 3, nchix = 6, ndm = 2000 - ! lmaxx : maximum non local angular momentum in PP + INTEGER, PARAMETER :: lmaxx = 3, nchix = 6, ndm = 2000 + ! lmaxx : maximum non local angular momentum in PP ! nchix : maximum number of atomic wavefunctions per PP ! ndm : maximum number of points in the radial mesh - integer, parameter :: nbrx = 8, lqmax = 5, nqfx = 8 - ! nbrx : maximum number of beta functions - ! lqmax : maximum number of angular momentum of Q + INTEGER, PARAMETER :: nbrx = 8, lqmax = 5, nqfx = 8 + ! nbrx : maximum number of beta functions + ! lqmax : maximum number of angular momentum of Q ! nqfx : maximum number of coefficients in Q smoothing ! ! pp_header - character (len=80):: generated, date_author, comment - character (len=2) :: psd(npsx), pseudotype - character (len=20):: dft(npsx) - integer :: lmax(npsx), mesh(npsx), nbeta(npsx), ntwfc(npsx) - logical :: nlcc(npsx), isus(npsx) + CHARACTER (len=80):: generated, date_author, comment + CHARACTER (len=2) :: psd(npsx), pseudotype + CHARACTER (len=20):: dft(npsx) + INTEGER :: lmax(npsx), mesh(npsx), nbeta(npsx), ntwfc(npsx) + LOGICAL :: nlcc(npsx), isus(npsx) real(8) :: zp(npsx), ecutrho, ecutwfc, etotps real(8) :: oc(nchix,npsx) - character(len=2) :: els(nchix,npsx) - integer :: lchi(nchix,npsx) + CHARACTER(len=2) :: els(nchix,npsx) + INTEGER :: lchi(nchix,npsx) ! ! pp_mesh real(8) :: r(ndm,npsx), rab(ndm,npsx) @@ -43,11 +43,11 @@ module pseudo ! pp_nonlocal ! pp_beta real(8) :: betar(ndm, nbrx, npsx) - integer :: lll(nbrx,npsx), ikk2(nbrx,npsx) + INTEGER :: lll(nbrx,npsx), ikk2(nbrx,npsx) ! pp_dij real(8) :: dion(nbrx,nbrx,npsx) ! pp_qij - integer :: nqf(npsx), nqlc(npsx) + INTEGER :: nqf(npsx), nqlc(npsx) real(8) :: rinner(lqmax,npsx), qqq(nbrx,nbrx,npsx), & qfunc(ndm,nbrx,nbrx,npsx) ! pp_qfcoef @@ -58,339 +58,339 @@ module pseudo ! ! pp_rhoatom real(8) :: rho_at(ndm,npsx) -end module pseudo +END MODULE pseudo ! !--------------------------------------------------------------------- -program read_ps +PROGRAM read_ps !--------------------------------------------------------------------- ! ! Read pseudopotentials in the Unified Pseudopotential Format (UPF) ! - implicit none - integer :: is, ios, iunps = 4 - character (len=256) :: filein + IMPLICIT NONE + INTEGER :: is, ios, iunps = 4 + CHARACTER (len=256) :: filein ! is = 0 -10 print '('' Input PP file # '',i2,'' in UPF format > '',$)', is+1 - read (5, '(a)', end = 20, err = 20) filein - open(unit=iunps,file=filein,status='old',form='formatted',iostat=ios) - if (ios.ne.0) stop +10 PRINT '('' Input PP file # '',i2,'' in UPF format > '',$)', is+1 + READ (5, '(a)', end = 20, err = 20) filein + OPEN(unit=iunps,file=filein,status='old',form='formatted',iostat=ios) + IF (ios/=0) STOP is = is + 1 - call read_pseudo(is, iunps) - close (unit=iunps) - go to 10 -20 stop -end program read_ps + CALL read_pseudo(is, iunps) + CLOSE (unit=iunps) + GOTO 10 +20 STOP +END PROGRAM read_ps ! !--------------------------------------------------------------------- -subroutine read_pseudo (is, iunps) +SUBROUTINE read_pseudo (is, iunps) !--------------------------------------------------------------------- ! ! Read pseudopotential in the Unified Pseudopotential Format (UPF) ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! is : index of this pseudopotential ! iunps: unit connected with pseudopotential file ! - if (is < 0 .or. is > npsx ) call errore ('read_pseudo', 'Wrong is number', 1) - write ( *, * ) " Reading pseudopotential file in UPF format..." + IF (is < 0 .or. is > npsx ) CALL errore ('read_pseudo', 'Wrong is number', 1) + WRITE ( *, * ) " Reading pseudopotential file in UPF format..." !------->Search for Header - call scan_begin (iunps, "HEADER", .true.) - call read_pseudo_header (is, iunps) - call scan_end (iunps, "HEADER") + CALL scan_begin (iunps, "HEADER", .true.) + CALL read_pseudo_header (is, iunps) + CALL scan_end (iunps, "HEADER") !-------->Search for mesh information - call scan_begin (iunps, "MESH", .true.) - call read_pseudo_mesh (is, iunps) - call scan_end (iunps, "MESH") + CALL scan_begin (iunps, "MESH", .true.) + CALL read_pseudo_mesh (is, iunps) + CALL scan_end (iunps, "MESH") !-------->If present, search for nlcc - if (nlcc (is) ) then - call scan_begin (iunps, "NLCC", .true.) - call read_pseudo_nlcc (is, iunps) - call scan_end (iunps, "NLCC") - endif + IF (nlcc (is) ) THEN + CALL scan_begin (iunps, "NLCC", .true.) + CALL read_pseudo_nlcc (is, iunps) + CALL scan_end (iunps, "NLCC") + ENDIF !-------->Search for Local potential - call scan_begin (iunps, "LOCAL", .true.) - call read_pseudo_local (is, iunps) - call scan_end (iunps, "LOCAL") + CALL scan_begin (iunps, "LOCAL", .true.) + CALL read_pseudo_local (is, iunps) + CALL scan_end (iunps, "LOCAL") !-------->Search for Nonlocal potential - call scan_begin (iunps, "NONLOCAL", .true.) - call read_pseudo_nl (is, iunps) - call scan_end (iunps, "NONLOCAL") + CALL scan_begin (iunps, "NONLOCAL", .true.) + CALL read_pseudo_nl (is, iunps) + CALL scan_end (iunps, "NONLOCAL") !-------->Search for atomic wavefunctions - call scan_begin (iunps, "PSWFC", .true.) - call read_pseudo_pswfc (is, iunps) - call scan_end (iunps, "PSWFC") + CALL scan_begin (iunps, "PSWFC", .true.) + CALL read_pseudo_pswfc (is, iunps) + CALL scan_end (iunps, "PSWFC") !-------->Search for atomic charge - call scan_begin (iunps, "RHOATOM", .true.) - call read_pseudo_rhoatom (is, iunps) - call scan_end (iunps, "RHOATOM") + CALL scan_begin (iunps, "RHOATOM", .true.) + CALL read_pseudo_rhoatom (is, iunps) + CALL scan_end (iunps, "RHOATOM") ! - write ( *, * ) " ...done" - return -end subroutine read_pseudo + WRITE ( *, * ) " ...done" + RETURN +END SUBROUTINE read_pseudo !--------------------------------------------------------------------- -subroutine scan_begin (iunps, string, rew) +SUBROUTINE scan_begin (iunps, string, rew) !--------------------------------------------------------------------- ! - implicit none + IMPLICIT NONE ! Unit of the input file - integer :: iunps + INTEGER :: iunps ! Label to be matched - character (len=*) :: string - logical :: rew + CHARACTER (len=*) :: string + LOGICAL :: rew ! Flag: if .true. rewind the file - character (len=80) :: rstring + CHARACTER (len=80) :: rstring ! String read from file - integer :: ios - logical, external :: matches + INTEGER :: ios + LOGICAL, EXTERNAL :: matches ios = 0 - if (rew) rewind (iunps) - do while (ios.eq.0) - read (iunps, *, iostat = ios, err = 300) rstring - if (matches ("", rstring) ) return - enddo -300 call errore ('scan_begin', 'No '//string//' block', abs (ios) ) + IF (rew) REWIND (iunps) + DO WHILE (ios==0) + READ (iunps, *, iostat = ios, err = 300) rstring + IF (matches ("", rstring) ) RETURN + ENDDO +300 CALL errore ('scan_begin', 'No '//string//' block', abs (ios) ) -end subroutine scan_begin +END SUBROUTINE scan_begin !--------------------------------------------------------------------- -subroutine scan_end (iunps, string) +SUBROUTINE scan_end (iunps, string) !--------------------------------------------------------------------- - implicit none + IMPLICIT NONE ! Unit of the input file - integer :: iunps + INTEGER :: iunps ! Label to be matched - character (len=*) :: string + CHARACTER (len=*) :: string ! String read from file - character (len=80) :: rstring - integer :: ios - logical, external :: matches + CHARACTER (len=80) :: rstring + INTEGER :: ios + LOGICAL, EXTERNAL :: matches - read (iunps, '(a)', iostat = ios, err = 300) rstring - if (matches ("", rstring) ) return -300 call errore ('scan_end', & + READ (iunps, '(a)', iostat = ios, err = 300) rstring + IF (matches ("", rstring) ) RETURN +300 CALL errore ('scan_end', & 'No '//string//' block end statement, possibly corrupted file', - 1) -end subroutine scan_end +END SUBROUTINE scan_end ! !--------------------------------------------------------------------- -subroutine read_pseudo_header (is, iunps) +SUBROUTINE read_pseudo_header (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: nv, ios, nw - character (len=75) :: dummy - logical, external :: matches + INTEGER :: nv, ios, nw + CHARACTER (len=75) :: dummy + LOGICAL, EXTERNAL :: matches - read (iunps, *, err = 100, iostat = ios) nv, dummy - read (iunps, *, err = 100, iostat = ios) psd (is), dummy - read (iunps, *, err = 100, iostat = ios) pseudotype - if (matches (pseudotype, "US") ) isus (is) = .true. - read (iunps, *, err = 100, iostat = ios) nlcc (is), dummy - read (iunps, '(a20,t24,a)', err = 100, iostat = ios) dft(is), dummy - read (iunps, * ) zp (is), dummy - read (iunps, * ) etotps, dummy - read (iunps, * ) ecutwfc, ecutrho - read (iunps, * ) lmax (is), dummy - read (iunps, *, err = 100, iostat = ios) mesh (is), dummy - read (iunps, *, err = 100, iostat = ios) ntwfc(is), nbeta (is), dummy - read (iunps, '(a)', err = 100, iostat = ios) dummy - do nw = 1, ntwfc(is) - read (iunps, * ) els (nw,is), lchi (nw, is), oc (nw, is) - enddo - return -100 call errore ('read_pseudo_header', 'Reading pseudo file', abs (ios)) -end subroutine read_pseudo_header + READ (iunps, *, err = 100, iostat = ios) nv, dummy + READ (iunps, *, err = 100, iostat = ios) psd (is), dummy + READ (iunps, *, err = 100, iostat = ios) pseudotype + IF (matches (pseudotype, "US") ) isus (is) = .true. + READ (iunps, *, err = 100, iostat = ios) nlcc (is), dummy + READ (iunps, '(a20,t24,a)', err = 100, iostat = ios) dft(is), dummy + READ (iunps, * ) zp (is), dummy + READ (iunps, * ) etotps, dummy + READ (iunps, * ) ecutwfc, ecutrho + READ (iunps, * ) lmax (is), dummy + READ (iunps, *, err = 100, iostat = ios) mesh (is), dummy + READ (iunps, *, err = 100, iostat = ios) ntwfc(is), nbeta (is), dummy + READ (iunps, '(a)', err = 100, iostat = ios) dummy + DO nw = 1, ntwfc(is) + READ (iunps, * ) els (nw,is), lchi (nw, is), oc (nw, is) + ENDDO + RETURN +100 CALL errore ('read_pseudo_header', 'Reading pseudo file', abs (ios)) +END SUBROUTINE read_pseudo_header ! !--------------------------------------------------------------------- -subroutine read_pseudo_local (is, iunps) +SUBROUTINE read_pseudo_local (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios ! - read (iunps, *, err=100, iostat=ios) (vloc0(ir,is) , ir=1,mesh(is)) + READ (iunps, *, err=100, iostat=ios) (vloc0(ir,is) , ir=1,mesh(is)) -100 call errore ('read_pseudo_local','Reading pseudo file', abs(ios) ) +100 CALL errore ('read_pseudo_local','Reading pseudo file', abs(ios) ) - return -end subroutine read_pseudo_local + RETURN +END SUBROUTINE read_pseudo_local ! !--------------------------------------------------------------------- -subroutine read_pseudo_mesh (is, iunps) +SUBROUTINE read_pseudo_mesh (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios ! - call scan_begin (iunps, "R", .false.) - read (iunps, *, err = 100, iostat = ios) (r(ir,is), ir=1,mesh(is) ) - call scan_end (iunps, "R") - call scan_begin (iunps, "RAB", .false.) - read (iunps, *, err = 100, iostat = ios) (rab(ir,is), ir=1,mesh(is) ) - call scan_end (iunps, "RAB") + CALL scan_begin (iunps, "R", .false.) + READ (iunps, *, err = 100, iostat = ios) (r(ir,is), ir=1,mesh(is) ) + CALL scan_end (iunps, "R") + CALL scan_begin (iunps, "RAB", .false.) + READ (iunps, *, err = 100, iostat = ios) (rab(ir,is), ir=1,mesh(is) ) + CALL scan_end (iunps, "RAB") - return + RETURN -100 call errore ('read_pseudo_mesh', 'Reading pseudo file', abs (ios) ) -end subroutine read_pseudo_mesh +100 CALL errore ('read_pseudo_mesh', 'Reading pseudo file', abs (ios) ) +END SUBROUTINE read_pseudo_mesh ! !--------------------------------------------------------------------- -subroutine read_pseudo_nl (is, iunps) +SUBROUTINE read_pseudo_nl (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: nb, mb, n, ir, nd, ios, idum, ldum, icon, lp, i + INTEGER :: nb, mb, n, ir, nd, ios, idum, ldum, icon, lp, i ! counters - character (len=75) :: dummy + CHARACTER (len=75) :: dummy ! - do nb = 1, nbeta (is) - call scan_begin (iunps, "BETA", .false.) - read (iunps, *, err = 100, iostat = ios) idum, lll(nb,is), dummy - read (iunps, '(i6)', err = 100, iostat = ios) ikk2(nb,is) - read (iunps, *, err = 100, iostat = ios) & + DO nb = 1, nbeta (is) + CALL scan_begin (iunps, "BETA", .false.) + READ (iunps, *, err = 100, iostat = ios) idum, lll(nb,is), dummy + READ (iunps, '(i6)', err = 100, iostat = ios) ikk2(nb,is) + READ (iunps, *, err = 100, iostat = ios) & (betar(ir,nb,is), ir=1,ikk2(nb,is)) - do ir = ikk2(nb,is) + 1, mesh (is) - betar (ir, nb, is) = 0.d0 - enddo - call scan_end (iunps, "BETA") - enddo + DO ir = ikk2(nb,is) + 1, mesh (is) + betar (ir, nb, is) = 0.d0 + ENDDO + CALL scan_end (iunps, "BETA") + ENDDO - call scan_begin (iunps, "DIJ", .false.) - read (iunps, *, err = 100, iostat = ios) nd, dummy + CALL scan_begin (iunps, "DIJ", .false.) + READ (iunps, *, err = 100, iostat = ios) nd, dummy dion (:,:,is) = 0.d0 - do icon = 1, nd - read (iunps, *, err = 100, iostat = ios) nb, mb, dion(nb,mb,is) - dion (mb,nb,is) = dion (nb,mb,is) - enddo - call scan_end (iunps, "DIJ") + DO icon = 1, nd + READ (iunps, *, err = 100, iostat = ios) nb, mb, dion(nb,mb,is) + dion (mb,nb,is) = dion (nb,mb,is) + ENDDO + CALL scan_end (iunps, "DIJ") - if (isus (is) ) then - call scan_begin (iunps, "QIJ", .false.) - read (iunps, *, err = 100, iostat = ios) nqf(is) + IF (isus (is) ) THEN + CALL scan_begin (iunps, "QIJ", .false.) + READ (iunps, *, err = 100, iostat = ios) nqf(is) nqlc (is)= 2 * lmax (is) + 1 - if (nqlc(is).gt.lqmax .or. nqlc(is).lt.0) & - call errore (' read_pseudo_nl', 'Wrong nqlc', nqlc (is) ) - if (nqf(is).ne.0) then - call scan_begin (iunps, "RINNER", .false.) - read (iunps,*,err=100,iostat=ios) & + IF (nqlc(is)>lqmax .or. nqlc(is)<0) & + CALL errore (' read_pseudo_nl', 'Wrong nqlc', nqlc (is) ) + IF (nqf(is)/=0) THEN + CALL scan_begin (iunps, "RINNER", .false.) + READ (iunps,*,err=100,iostat=ios) & (idum,rinner(i,is),i=1,nqlc(is)) - call scan_end (iunps, "RINNER") - end if - do nb = 1, nbeta(is) - do mb = nb, nbeta(is) + CALL scan_end (iunps, "RINNER") + ENDIF + DO nb = 1, nbeta(is) + DO mb = nb, nbeta(is) - read (iunps,*,err=100,iostat=ios) idum, idum, ldum, dummy + READ (iunps,*,err=100,iostat=ios) idum, idum, ldum, dummy !" i j (l)" - if (ldum.ne.lll(mb,is) ) call errore ('read_pseudo_nl', & + IF (ldum/=lll(mb,is) ) CALL errore ('read_pseudo_nl', & 'inconsistent angular momentum for Q_ij', 1) - read (iunps,*,err=100,iostat=ios) qqq(nb,mb,is), dummy + READ (iunps,*,err=100,iostat=ios) qqq(nb,mb,is), dummy ! "Q_int" - qqq(mb,nb,is) = qqq(nb,mb,is) + qqq(mb,nb,is) = qqq(nb,mb,is) - read (iunps,*,err=100,iostat=ios) & + READ (iunps,*,err=100,iostat=ios) & (qfunc(n,nb,mb,is), n=1,mesh(is)) - do n = 0, mesh (is) - qfunc(n,mb,nb,is) = qfunc(n,nb,mb,is) - enddo + DO n = 0, mesh (is) + qfunc(n,mb,nb,is) = qfunc(n,nb,mb,is) + ENDDO - if (nqf(is).gt.0) then - call scan_begin (iunps, "QFCOEF", .false.) - read (iunps,*,err=100,iostat=ios) & + IF (nqf(is)>0) THEN + CALL scan_begin (iunps, "QFCOEF", .false.) + READ (iunps,*,err=100,iostat=ios) & ((qfcoef(i,lp,nb,mb,is),i=1,nqf(is)),lp=1,nqlc(is)) - call scan_end (iunps, "QFCOEF") - end if + CALL scan_end (iunps, "QFCOEF") + ENDIF - enddo - enddo - call scan_end (iunps, "QIJ") - else + ENDDO + ENDDO + CALL scan_end (iunps, "QIJ") + ELSE qqq (:,:,is) = 0.d0 qfunc(:,:,:,is) =0.d0 - endif + ENDIF -100 call errore ('read_pseudo_nl', 'Reading pseudo file', abs (ios) ) - return -end subroutine read_pseudo_nl +100 CALL errore ('read_pseudo_nl', 'Reading pseudo file', abs (ios) ) + RETURN +END SUBROUTINE read_pseudo_nl ! !--------------------------------------------------------------------- -subroutine read_pseudo_nlcc (is, iunps) +SUBROUTINE read_pseudo_nlcc (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios - read (iunps, *, err = 100, iostat = ios) (rho_atc(ir,is), ir=1,mesh(is) ) + READ (iunps, *, err = 100, iostat = ios) (rho_atc(ir,is), ir=1,mesh(is) ) ! -100 call errore ('read_pseudo_nlcc', 'Reading pseudo file', abs (ios) ) - return -end subroutine read_pseudo_nlcc +100 CALL errore ('read_pseudo_nlcc', 'Reading pseudo file', abs (ios) ) + RETURN +END SUBROUTINE read_pseudo_nlcc ! !--------------------------------------------------------------------- -subroutine read_pseudo_pswfc (is, iunps) +SUBROUTINE read_pseudo_pswfc (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - character (len=75) :: dummy - integer :: nb, ir, ios + CHARACTER (len=75) :: dummy + INTEGER :: nb, ir, ios ! - do nb = 1, ntwfc(is) - read (iunps,*,err=100,iostat=ios) dummy !Wavefunction labels - read (iunps,*,err=100,iostat=ios) (chi(ir,nb,is), ir=1,mesh(is)) - enddo -100 call errore ('read_pseudo_pswfc', 'Reading pseudo file', abs(ios)) - return + DO nb = 1, ntwfc(is) + READ (iunps,*,err=100,iostat=ios) dummy !Wavefunction labels + READ (iunps,*,err=100,iostat=ios) (chi(ir,nb,is), ir=1,mesh(is)) + ENDDO +100 CALL errore ('read_pseudo_pswfc', 'Reading pseudo file', abs(ios)) + RETURN -end subroutine read_pseudo_pswfc +END SUBROUTINE read_pseudo_pswfc ! !--------------------------------------------------------------------- -subroutine read_pseudo_rhoatom (is, iunps) +SUBROUTINE read_pseudo_rhoatom (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios - read (iunps,*,err=100,iostat=ios) (rho_at(ir,is), ir=1,mesh(is)) - return + READ (iunps,*,err=100,iostat=ios) (rho_at(ir,is), ir=1,mesh(is)) + RETURN -100 call errore ('read_pseudo_rhoatom','Reading pseudo file',abs(ios)) +100 CALL errore ('read_pseudo_rhoatom','Reading pseudo file',abs(ios)) -end subroutine read_pseudo_rhoatom +END SUBROUTINE read_pseudo_rhoatom diff --git a/upftools/read_upf_tofile.f90 b/upftools/read_upf_tofile.f90 index 9287d0c65..4a4b6bbfe 100644 --- a/upftools/read_upf_tofile.f90 +++ b/upftools/read_upf_tofile.f90 @@ -6,10 +6,10 @@ ! or http://www.gnu.org/copyleft/gpl.txt . ! !--------------------------------------------------------------------- -PROGRAM read_upf_tofile +PROGRAM read_upf_tofile !--------------------------------------------------------------------- ! - ! This small program reads the pseudopotential in the Unified + ! This small program reads the pseudopotential in the Unified ! Pseudopotential Format and writes three files ! in a format which can be plotted. The files are: ! @@ -21,10 +21,10 @@ PROGRAM read_upf_tofile ! PWSCF modules ! ! - USE constants, only : fpi + USE constants, ONLY : fpi USE pseudo_types USE upf_module - USE radial_grids, only : radial_grid_type, nullify_radial_grid + USE radial_grids, ONLY : radial_grid_type, nullify_radial_grid ! IMPLICIT NONE ! @@ -38,11 +38,11 @@ PROGRAM read_upf_tofile TYPE (pseudo_upf) :: upf TYPE (radial_grid_type) :: grid ! - WRITE(6,'("Name of the upf file > ", $)') - READ(5,'(a)') file_pseudo + WRITE(6,'("Name of the upf file > ", $)') + READ(5,'(a)') file_pseudo ! nullify objects as soon as they are instantiated - + CALL nullify_pseudo_upf( upf ) CALL nullify_radial_grid( grid ) @@ -53,39 +53,39 @@ PROGRAM read_upf_tofile CALL read_upf(upf, grid, ierr, unit=iunps) ! - IF (ierr .NE. 0) & - CALL errore('read_upf_tofile','reading pseudo upf', ABS(ierr)) + IF (ierr /= 0) & + CALL errore('read_upf_tofile','reading pseudo upf', abs(ierr)) ! CLOSE(iunps) ! OPEN(UNIT=iunps,FILE='filewfc',STATUS='unknown',FORM='formatted', & ERR=200, IOSTAT=ios) -200 CALL errore('read_upf_tofile','open error on file filewfc',ABS(ios)) +200 CALL errore('read_upf_tofile','open error on file filewfc',abs(ios)) DO n=1,upf%mesh WRITE(iunps,'(30f12.6)') upf%r(n), (upf%chi(n,j), j=1,upf%nwfc) - END DO + ENDDO CLOSE(iunps) OPEN(UNIT=iunps,FILE='filebeta',STATUS='unknown',FORM='formatted', & ERR=300, IOSTAT=ios) -300 CALL errore('read_upf_tofile','open error on file filebeta',ABS(ios)) +300 CALL errore('read_upf_tofile','open error on file filebeta',abs(ios)) DO n=1,upf%mesh WRITE(iunps,'(30f12.6)') upf%r(n), (upf%beta(n,j), j=1,upf%nbeta) - END DO + ENDDO CLOSE(iunps) OPEN(UNIT=iunps,FILE='filepot',STATUS='unknown',FORM='formatted', & ERR=400, IOSTAT=ios) -400 CALL errore('read_upf_tofile','open error on file filepot',ABS(ios)) +400 CALL errore('read_upf_tofile','open error on file filepot',abs(ios)) DO n=1,upf%mesh WRITE(iunps,'(4f12.6)') upf%r(n), upf%vloc(n), & - upf%rho_at(n), upf%rho_atc(n)*fpi*upf%r(n)**2 - END DO + upf%rho_at(n), upf%rho_atc(n)*fpi*upf%r(n)**2 + ENDDO CLOSE(iunps) diff --git a/upftools/rrkj2upf.f90 b/upftools/rrkj2upf.f90 index 1f25114de..c54806287 100644 --- a/upftools/rrkj2upf.f90 +++ b/upftools/rrkj2upf.f90 @@ -7,178 +7,178 @@ ! ! !--------------------------------------------------------------------- -program rrkj2upf +PROGRAM rrkj2upf !--------------------------------------------------------------------- ! ! Convert a pseudopotential written in "rrkj3" format ! (Rabe-Rappe-Kaxiras-Joannopoulos with 3 Bessel functions) ! to unified pseudopotential format ! - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open (unit = 1, file = filein, status = 'old', form = 'formatted') - call read_rrkj(1) - close (1) + CALL get_file ( filein ) + OPEN (unit = 1, file = filein, status = 'old', form = 'formatted') + CALL read_rrkj(1) + CLOSE (1) ! convert variables read from rrkj3 format into those needed ! by the upf format - add missing quantities - call convert_rrkj + CALL convert_rrkj fileout=trim(filein)//'.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) -stop -20 write (6,'("rrkj2upf: error reading pseudopotential file name")') - stop +STOP +20 WRITE (6,'("rrkj2upf: error reading pseudopotential file name")') + STOP -end program rrkj2upf +END PROGRAM rrkj2upf -module rrkj3 +MODULE rrkj3 ! ! All variables read from RRKJ3 file format - ! + ! ! trailing underscore means that a variable with the same name ! is used in module 'upf' containing variables to be written ! - character(len=75):: titleps - character (len=2), allocatable :: els_(:) - integer :: pseudotype_, iexch_, icorr_, igcx_, igcc_, mesh_, & + CHARACTER(len=75):: titleps + CHARACTER (len=2), ALLOCATABLE :: els_(:) + INTEGER :: pseudotype_, iexch_, icorr_, igcx_, igcc_, mesh_, & nwfs_, nbeta_, lmax_ - logical :: rel_, nlcc_ + LOGICAL :: rel_, nlcc_ real (8) :: zp_, etotps_, xmin, rmax, zmesh, dx, rcloc_ - integer, allocatable:: lchi_(:), nns_(:), ikk2_(:) - real (8), allocatable :: rcut_(:), rcutus_(:), oc_(:), & + INTEGER, ALLOCATABLE:: lchi_(:), nns_(:), ikk2_(:) + real (8), ALLOCATABLE :: rcut_(:), rcutus_(:), oc_(:), & beta(:,:), dion_(:,:), qqq_(:,:), ddd(:,:), qfunc_(:,:,:), & rho_atc_(:), rho_at_(:), chi_(:,:), vloc_(:) -end module rrkj3 -! +END MODULE rrkj3 +! ! ---------------------------------------------------------- -subroutine read_rrkj(iunps) +SUBROUTINE read_rrkj(iunps) ! ---------------------------------------------------------- - ! - use rrkj3 - implicit none - integer :: iunps - integer :: nb, mb, n, ir, ios + ! + USE rrkj3 + IMPLICIT NONE + INTEGER :: iunps + INTEGER :: nb, mb, n, ir, ios !--- > Start the header reading - read (iunps, '(a75)', err = 100) titleps - read (iunps, *, err = 100) pseudotype_ - read (iunps, *, err = 100) rel_, nlcc_ - read (iunps, *, err=100) iexch_, icorr_, igcx_, igcc_ - read (iunps, '(2e17.11,i5)') zp_, etotps_, lmax_ - read (iunps, '(4e17.11,i5)', err=100) xmin, rmax, zmesh, dx, mesh_ - read (iunps, *, err=100) nwfs_, nbeta_ + READ (iunps, '(a75)', err = 100) titleps + READ (iunps, *, err = 100) pseudotype_ + READ (iunps, *, err = 100) rel_, nlcc_ + READ (iunps, *, err=100) iexch_, icorr_, igcx_, igcc_ + READ (iunps, '(2e17.11,i5)') zp_, etotps_, lmax_ + READ (iunps, '(4e17.11,i5)', err=100) xmin, rmax, zmesh, dx, mesh_ + READ (iunps, *, err=100) nwfs_, nbeta_ - allocate(rcut_(nwfs_), rcutus_(nwfs_)) - read (iunps, *, err=100) (rcut_(nb), nb=1,nwfs_) - read (iunps, *, err=100) (rcutus_(nb), nb=1,nwfs_) + ALLOCATE(rcut_(nwfs_), rcutus_(nwfs_)) + READ (iunps, *, err=100) (rcut_(nb), nb=1,nwfs_) + READ (iunps, *, err=100) (rcutus_(nb), nb=1,nwfs_) - allocate(els_(nwfs_), nns_(nwfs_), lchi_(nwfs_), oc_(nwfs_)) - do nb = 1, nwfs_ - read (iunps, '(a2,2i3,f6.2)', err = 100) els_(nb), & + ALLOCATE(els_(nwfs_), nns_(nwfs_), lchi_(nwfs_), oc_(nwfs_)) + DO nb = 1, nwfs_ + READ (iunps, '(a2,2i3,f6.2)', err = 100) els_(nb), & nns_(nb), lchi_(nb) , oc_(nb) - enddo + ENDDO - allocate(ikk2_(nbeta_)) - allocate(beta( mesh_,nbeta_)) - allocate(dion_(nbeta_,nbeta_)) - allocate(ddd (nbeta_,nbeta_)) - allocate(qqq_(nbeta_,nbeta_)) - allocate(qfunc_(mesh_,nbeta_,nbeta_)) + ALLOCATE(ikk2_(nbeta_)) + ALLOCATE(beta( mesh_,nbeta_)) + ALLOCATE(dion_(nbeta_,nbeta_)) + ALLOCATE(ddd (nbeta_,nbeta_)) + ALLOCATE(qqq_(nbeta_,nbeta_)) + ALLOCATE(qfunc_(mesh_,nbeta_,nbeta_)) - do nb = 1, nbeta_ - read (iunps, *, err = 100) ikk2_(nb) - read (iunps, *, err = 100) (beta (ir, nb) , ir = 1,ikk2_(nb) ) - do ir = ikk2_(nb) + 1, mesh_ - beta (ir, nb) = 0.d0 - enddo - do mb = 1, nb - read (iunps, *, err = 100) dion_(nb, mb) - dion_(mb, nb) = dion_(nb, mb) - if (pseudotype_.eq.3) then - read (iunps, *, err = 100) qqq_(nb, mb) - qqq_(mb, nb) = qqq_(nb, mb) - read (iunps, *, err = 100) (qfunc_(n,nb, mb), n = 1, mesh_) - do n = 1, mesh_ - qfunc_(n, mb, nb) = qfunc_(n, nb, mb) - enddo - else - qqq_(nb, mb) = 0.d0 - qqq_(mb, nb) = 0.d0 - do n = 1, mesh_ - qfunc_(n, nb, mb) = 0.d0 - qfunc_(n, mb, nb) = 0.d0 - enddo - endif - enddo - enddo + DO nb = 1, nbeta_ + READ (iunps, *, err = 100) ikk2_(nb) + READ (iunps, *, err = 100) (beta (ir, nb) , ir = 1,ikk2_(nb) ) + DO ir = ikk2_(nb) + 1, mesh_ + beta (ir, nb) = 0.d0 + ENDDO + DO mb = 1, nb + READ (iunps, *, err = 100) dion_(nb, mb) + dion_(mb, nb) = dion_(nb, mb) + IF (pseudotype_==3) THEN + READ (iunps, *, err = 100) qqq_(nb, mb) + qqq_(mb, nb) = qqq_(nb, mb) + READ (iunps, *, err = 100) (qfunc_(n,nb, mb), n = 1, mesh_) + DO n = 1, mesh_ + qfunc_(n, mb, nb) = qfunc_(n, nb, mb) + ENDDO + ELSE + qqq_(nb, mb) = 0.d0 + qqq_(mb, nb) = 0.d0 + DO n = 1, mesh_ + qfunc_(n, nb, mb) = 0.d0 + qfunc_(n, mb, nb) = 0.d0 + ENDDO + ENDIF + ENDDO + ENDDO ! ! read the local potential ! - allocate(vloc_(mesh_)) - read (iunps, *, err = 100) rcloc_, (vloc_(ir ) , ir = 1, mesh_ ) + ALLOCATE(vloc_(mesh_)) + READ (iunps, *, err = 100) rcloc_, (vloc_(ir ) , ir = 1, mesh_ ) ! ! read the atomic charge ! - allocate(rho_at_(mesh_)) - read (iunps, *, err=100) (rho_at_(ir), ir=1,mesh_) + ALLOCATE(rho_at_(mesh_)) + READ (iunps, *, err=100) (rho_at_(ir), ir=1,mesh_) ! ! if present read the core charge ! - allocate(rho_atc_(mesh_)) - if (nlcc_) then - read (iunps, *, err=100) (rho_atc_(ir), ir=1, mesh_) - endif + ALLOCATE(rho_atc_(mesh_)) + IF (nlcc_) THEN + READ (iunps, *, err=100) (rho_atc_(ir), ir=1, mesh_) + ENDIF ! ! read the pseudo wavefunctions of the atom ! - allocate(chi_(mesh_,nwfs_)) - read (iunps, *, err=100) ( (chi_(ir,nb), ir = 1,mesh_) , nb = 1, nwfs_) + ALLOCATE(chi_(mesh_,nwfs_)) + READ (iunps, *, err=100) ( (chi_(ir,nb), ir = 1,mesh_) , nb = 1, nwfs_) ! ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully read' + WRITE (6,'(a)') 'Pseudopotential successfully read' ! ---------------------------------------------------------- ! - return -100 write (6,'("read_rrkj: error reading pseudopotential file")') - stop + RETURN +100 WRITE (6,'("read_rrkj: error reading pseudopotential file")') + STOP -end subroutine read_rrkj +END SUBROUTINE read_rrkj -subroutine convert_rrkj +SUBROUTINE convert_rrkj ! ---------------------------------------------------------- ! - use rrkj3 - use upf - use constants, only : fpi - implicit none - integer i, n + USE rrkj3 + USE upf + USE constants, ONLY : fpi + IMPLICIT NONE + INTEGER i, n real(8) :: x - write(generated, '("Generated using Andrea Dal Corso code (rrkj3)")') - write(date_author,'("Author: Andrea Dal Corso Generation date: unknown")') + WRITE(generated, '("Generated using Andrea Dal Corso code (rrkj3)")') + WRITE(date_author,'("Author: Andrea Dal Corso Generation date: unknown")') comment = 'Info:'//titleps - if (rel_) then + IF (rel_) THEN rel = 1 - else + ELSE rel = 0 - end if + ENDIF rcloc = rcloc_ nwfs = nwfs_ - allocate( els(nwfs), oc(nwfs), epseu(nwfs)) - allocate(lchi(nwfs), nns(nwfs) ) - allocate(rcut (nwfs), rcutus (nwfs)) - do i=1, nwfs + ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs)) + ALLOCATE(lchi(nwfs), nns(nwfs) ) + ALLOCATE(rcut (nwfs), rcutus (nwfs)) + DO i=1, nwfs nns (i) = nns_(i) lchi(i) = lchi_(i) rcut(i) = rcut_(i) @@ -186,15 +186,15 @@ subroutine convert_rrkj oc (i) = oc_(i) els(i) = els_(i) epseu(i) = 0.0d0 - end do - deallocate (els_, oc_, rcutus_, rcut_, nns_) + ENDDO + DEALLOCATE (els_, oc_, rcutus_, rcut_, nns_) - psd = titleps (7:8) - if (pseudotype_.eq.3) then + psd = titleps (7:8) + IF (pseudotype_==3) THEN pseudotype = 'US' - else + ELSE pseudotype = 'NC' - endif + ENDIF nlcc = nlcc_ zp = zp_ etotps = etotps_ @@ -204,85 +204,85 @@ subroutine convert_rrkj mesh = mesh_ nbeta = nbeta_ ntwfc = 0 - do i=1, nwfs - if (oc(i) .gt. 1.0d-12) ntwfc = ntwfc + 1 - end do - allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) + DO i=1, nwfs + IF (oc(i) > 1.0d-12) ntwfc = ntwfc + 1 + ENDDO + ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) n = 0 - do i=1, nwfs - if (oc(i) .gt. 1.0d-12) then + DO i=1, nwfs + IF (oc(i) > 1.0d-12) THEN n = n + 1 elsw(n) = els(i) ocw (n) = oc (i) lchiw(n)=lchi(i) - end if - end do + ENDIF + ENDDO iexch = iexch_ icorr = icorr_ igcx = igcx_ igcc = igcc_ - allocate(rab(mesh)) - allocate( r(mesh)) + ALLOCATE(rab(mesh)) + ALLOCATE( r(mesh)) ! define logarithmic mesh - do i = 1, mesh - x = xmin + DBLE(i-1) * dx + DO i = 1, mesh + x = xmin + dble(i-1) * dx r (i) = exp(x) / zmesh rab(i) = dx * r(i) - end do + ENDDO - allocate (rho_atc(mesh)) + ALLOCATE (rho_atc(mesh)) ! rrkj rho_core(r) = 4pi*r^2*rho_core(r) UPF rho_atc (:) = rho_atc_(:) / fpi / r(:)**2 - deallocate (rho_atc_) + DEALLOCATE (rho_atc_) - allocate (vloc0(mesh)) + ALLOCATE (vloc0(mesh)) vloc0 = vloc_ - deallocate (vloc_) + DEALLOCATE (vloc_) - allocate(ikk2(nbeta), lll(nbeta)) + ALLOCATE(ikk2(nbeta), lll(nbeta)) ikk2 = ikk2_ lll = lchi_ - deallocate (ikk2_, lchi_) -! kkbeta = 0 + DEALLOCATE (ikk2_, lchi_) +! kkbeta = 0 ! do nb=1,nbeta -! kkbeta = max (kkbeta , ikk2(nb) ) +! kkbeta = max (kkbeta , ikk2(nb) ) ! end do - allocate(betar(mesh,nbeta)) + ALLOCATE(betar(mesh,nbeta)) betar = 0.0d0 - do i=1, nbeta + DO i=1, nbeta betar(1:ikk2(i),i) = beta(1:ikk2(i),i) - end do - deallocate (beta) + ENDDO + DEALLOCATE (beta) - allocate(dion(nbeta,nbeta)) + ALLOCATE(dion(nbeta,nbeta)) dion = dion_ - deallocate (dion_) + DEALLOCATE (dion_) - allocate(qqq(nbeta,nbeta)) + ALLOCATE(qqq(nbeta,nbeta)) qqq = qqq_ - deallocate (qqq_) + DEALLOCATE (qqq_) - allocate(qfunc(mesh,nbeta,nbeta)) + ALLOCATE(qfunc(mesh,nbeta,nbeta)) qfunc = qfunc_ nqf = 0 nqlc= 0 - allocate (rho_at(mesh)) + ALLOCATE (rho_at(mesh)) rho_at = rho_at_ - deallocate (rho_at_) + DEALLOCATE (rho_at_) - allocate (chi(mesh,ntwfc)) + ALLOCATE (chi(mesh,ntwfc)) n = 0 - do i=1, nwfs - if (oc(i) .gt. 1.0d-12) then + DO i=1, nwfs + IF (oc(i) > 1.0d-12) THEN n = n + 1 chi(:,n) = chi_(:,i) - end if - end do - deallocate (chi_) + ENDIF + ENDDO + DEALLOCATE (chi_) - return -end subroutine convert_rrkj + RETURN +END SUBROUTINE convert_rrkj diff --git a/upftools/uspp2upf.f90 b/upftools/uspp2upf.f90 index 05616f1b0..22a8d03c0 100644 --- a/upftools/uspp2upf.f90 +++ b/upftools/uspp2upf.f90 @@ -7,35 +7,35 @@ ! ! !--------------------------------------------------------------------- -program uspp2upf +PROGRAM uspp2upf !--------------------------------------------------------------------- ! ! Convert a pseudopotential written in Vanderbilt format ! (unformatted) to unified pseudopotential format ! - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open(unit=1,file=filein,status='old',form='unformatted') - call read_uspp(1) - close (unit=1) + CALL get_file ( filein ) + OPEN(unit=1,file=filein,status='old',form='unformatted') + CALL read_uspp(1) + CLOSE (unit=1) ! convert variables read from Vanderbilt format into those needed ! by the upf format - add missing quantities - call convert_uspp + CALL convert_uspp fileout=trim(filein)//'.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) - stop -20 write (6,'("uspp2upf: error reading pseudopotential file name")') - stop -end program uspp2upf + STOP +20 WRITE (6,'("uspp2upf: error reading pseudopotential file name")') + STOP +END PROGRAM uspp2upf diff --git a/upftools/vanderbilt.f90 b/upftools/vanderbilt.f90 index f4b9e164b..506751810 100644 --- a/upftools/vanderbilt.f90 +++ b/upftools/vanderbilt.f90 @@ -6,273 +6,273 @@ ! or http://www.gnu.org/copyleft/gpl.txt . ! ! -module Vanderbilt +MODULE Vanderbilt ! ! All variables read from Vanderbilt's file format - ! + ! ! trailing underscore means that a variable with the same name ! is used in module 'upf' containing variables to be written ! - integer :: nvalps, nang, nbeta_, kkbeta, nchi, ifpcor, keyps, & + INTEGER :: nvalps, nang, nbeta_, kkbeta, nchi, ifpcor, keyps, & mesh_, iver(3), idmy(3), nnlz, ifqopt, nqf_, irel, npf, & nlc, lloc real(8) :: z_, zp_, exfact, etot, eloc, rcloc_, rpcor, & qtryc, ptryc, rinner1_ - real(8), allocatable:: wwnlps(:), eeps(:), rinner_(:), rc(:), & + real(8), ALLOCATABLE:: wwnlps(:), eeps(:), rinner_(:), rc(:), & beta(:,:), ddd0(:,:), ddd(:,:), qqq_(:,:), eee(:), rho_atc_(:), & r_(:), rab_(:), rho_at_(:), qfunc_(:,:,:), vloc(:), vloc_(:), & wf(:,:), qfcoef_(:,:,:,:) - integer, allocatable :: lll_(:), nnlzps(:), iptype(:) - Character(len=20):: title -end module Vanderbilt -! + INTEGER, ALLOCATABLE :: lll_(:), nnlzps(:), iptype(:) + CHARACTER(len=20):: title +END MODULE Vanderbilt +! ! ---------------------------------------------------------- -subroutine read_uspp(iunit) +SUBROUTINE read_uspp(iunit) ! ---------------------------------------------------------- ! - use Vanderbilt - implicit none - integer :: iunit + USE Vanderbilt + IMPLICIT NONE + INTEGER :: iunit ! - integer :: i, j, k, lp + INTEGER :: i, j, k, lp real(8) :: rinner1 ! ! - read (iunit) (iver(i),i=1,3),(idmy(i),i=1,3) - read (iunit) title, z_, zp_, exfact, nvalps, mesh_, etot + READ (iunit) (iver(i),i=1,3),(idmy(i),i=1,3) + READ (iunit) title, z_, zp_, exfact, nvalps, mesh_, etot - allocate(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps)) - read (iunit) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps) + ALLOCATE(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps)) + READ (iunit) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps) - read (iunit) keyps, ifpcor, rinner1 + READ (iunit) keyps, ifpcor, rinner1 - if ( iver(1) .eq. 1 ) then + IF ( iver(1) == 1 ) THEN nang = nvalps nqf_ = 3 nlc = 5 - elseif ( iver(1) .eq. 2 ) then + ELSEIF ( iver(1) == 2 ) THEN nang = nvalps nqf_ = 3 nlc = 2 * nvalps - 1 - else if ( iver(1) .ge. 3 ) then - read (iunit) nang, lloc, eloc, ifqopt, nqf_, qtryc + ELSEIF ( iver(1) >= 3 ) THEN + READ (iunit) nang, lloc, eloc, ifqopt, nqf_, qtryc nlc = 2 * nang - 1 - endif + ENDIF - allocate(rinner_(2*nang-1)) + ALLOCATE(rinner_(2*nang-1)) rinner_(1) = rinner1 rinner1_ = rinner1 - if (10*iver(1)+iver(2).ge.51) & - read (iunit) (rinner_(i),i=1,nang*2-1) + IF (10*iver(1)+iver(2)>=51) & + READ (iunit) (rinner_(i),i=1,nang*2-1) - if ( iver(1) .ge. 4 ) then - read (iunit) irel - else + IF ( iver(1) >= 4 ) THEN + READ (iunit) irel + ELSE irel = 0 - end if + ENDIF - allocate(rc(nang)) - read (iunit) (rc(i),i=1,nang) + ALLOCATE(rc(nang)) + READ (iunit) (rc(i),i=1,nang) - read (iunit) nbeta_,kkbeta + READ (iunit) nbeta_,kkbeta ! - allocate(beta(kkbeta,nbeta_)) - allocate(qfunc_(kkbeta,nbeta_,nbeta_)) - allocate(ddd0(nbeta_,nbeta_)) - allocate(ddd (nbeta_,nbeta_)) - allocate(qqq_(nbeta_,nbeta_)) - allocate(lll_(nbeta_)) - allocate(eee(nbeta_)) - allocate(qfcoef_(nqf_,nlc,nbeta_,nbeta_)) + ALLOCATE(beta(kkbeta,nbeta_)) + ALLOCATE(qfunc_(kkbeta,nbeta_,nbeta_)) + ALLOCATE(ddd0(nbeta_,nbeta_)) + ALLOCATE(ddd (nbeta_,nbeta_)) + ALLOCATE(qqq_(nbeta_,nbeta_)) + ALLOCATE(lll_(nbeta_)) + ALLOCATE(eee(nbeta_)) + ALLOCATE(qfcoef_(nqf_,nlc,nbeta_,nbeta_)) ! - do j=1,nbeta_ - read (iunit) lll_(j),eee(j),(beta(i,j),i=1,kkbeta) - do k=j,nbeta_ - read (iunit) ddd0(j,k),ddd(j,k),qqq_(j,k), & + DO j=1,nbeta_ + READ (iunit) lll_(j),eee(j),(beta(i,j),i=1,kkbeta) + DO k=j,nbeta_ + READ (iunit) ddd0(j,k),ddd(j,k),qqq_(j,k), & (qfunc_(i,j,k),i=1,kkbeta), & ((qfcoef_(i,lp,j,k),i=1,nqf_),lp=1,2*nang-1) - end do - end do + ENDDO + ENDDO ! - allocate(iptype(nbeta_)) - if (10*iver(1)+iver(2).ge.72) & - read (iunit) (iptype(j),j=1,nbeta_),npf,ptryc + ALLOCATE(iptype(nbeta_)) + IF (10*iver(1)+iver(2)>=72) & + READ (iunit) (iptype(j),j=1,nbeta_),npf,ptryc ! - allocate(vloc_(mesh_)) - read (iunit) rcloc_,(vloc_(i),i=1,mesh_) + ALLOCATE(vloc_(mesh_)) + READ (iunit) rcloc_,(vloc_(i),i=1,mesh_) ! - allocate(rho_atc_(mesh_)) - if (ifpcor.gt.0) then - read (iunit) rpcor - read (iunit) (rho_atc_(i),i=1,mesh_) - end if + ALLOCATE(rho_atc_(mesh_)) + IF (ifpcor>0) THEN + READ (iunit) rpcor + READ (iunit) (rho_atc_(i),i=1,mesh_) + ENDIF ! - allocate(rho_at_(mesh_), vloc(mesh_)) - read (iunit) (vloc(i),i=1,mesh_) - read (iunit) (rho_at_(i),i=1,mesh_) + ALLOCATE(rho_at_(mesh_), vloc(mesh_)) + READ (iunit) (vloc(i),i=1,mesh_) + READ (iunit) (rho_at_(i),i=1,mesh_) - allocate(r_(mesh_), rab_(mesh_)) - read (iunit) (r_(i),i=1,mesh_) - read (iunit) (rab_(i),i=1,mesh_) - if (iver(1) .ge. 6) then + ALLOCATE(r_(mesh_), rab_(mesh_)) + READ (iunit) (r_(i),i=1,mesh_) + READ (iunit) (rab_(i),i=1,mesh_) + IF (iver(1) >= 6) THEN nchi = nvalps - if (iver(1) .ge. 7) read (iunit) nchi - allocate(wf(mesh_,nchi)) - read (iunit) ((wf(i,j), i=1,mesh_),j=1,nchi) - end if + IF (iver(1) >= 7) READ (iunit) nchi + ALLOCATE(wf(mesh_,nchi)) + READ (iunit) ((wf(i,j), i=1,mesh_),j=1,nchi) + ENDIF ! ! ---------------------------------------------------------- - write (6,'(a)') 'Pseudopotential successfully read' + WRITE (6,'(a)') 'Pseudopotential successfully read' ! ---------------------------------------------------------- ! -end subroutine read_uspp +END SUBROUTINE read_uspp ! ---------------------------------------------------------- ! ---------------------------------------------------------- -subroutine read_vdb(iunit) +SUBROUTINE read_vdb(iunit) ! ---------------------------------------------------------- ! - use Vanderbilt - implicit none - integer :: iunit + USE Vanderbilt + IMPLICIT NONE + INTEGER :: iunit ! - integer :: i, j, k, lp + INTEGER :: i, j, k, lp real(8) :: rinner1 ! ! - read(iunit, *) (iver(i),i=1,3),(idmy(i),i=1,3) - read(iunit,'(a20,3f15.9)' ) title, z_, zp_, exfact - read(iunit, *) nvalps, mesh_, etot + READ(iunit, *) (iver(i),i=1,3),(idmy(i),i=1,3) + READ(iunit,'(a20,3f15.9)' ) title, z_, zp_, exfact + READ(iunit, *) nvalps, mesh_, etot - allocate(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps)) - do i = 1,nvalps - read(iunit, *) nnlzps(i), wwnlps(i), eeps(i) - end do + ALLOCATE(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps)) + DO i = 1,nvalps + READ(iunit, *) nnlzps(i), wwnlps(i), eeps(i) + ENDDO - read(iunit, *) keyps, ifpcor, rinner1 + READ(iunit, *) keyps, ifpcor, rinner1 - if ( iver(1) .eq. 1 ) then + IF ( iver(1) == 1 ) THEN nang = nvalps nqf_ = 3 nlc = 5 - elseif ( iver(1) .eq. 2 ) then + ELSEIF ( iver(1) == 2 ) THEN nang = nvalps nqf_ = 3 nlc = 2 * nvalps - 1 - else if ( iver(1) .ge. 3 ) then - read(iunit, *) nang, lloc, eloc, ifqopt, nqf_, qtryc + ELSEIF ( iver(1) >= 3 ) THEN + READ(iunit, *) nang, lloc, eloc, ifqopt, nqf_, qtryc nlc = 2 * nang - 1 - endif + ENDIF - allocate(rinner_(2*nang-1)) + ALLOCATE(rinner_(2*nang-1)) rinner_(1) = rinner1 - if (10*iver(1)+iver(2).ge.51) & - read (iunit, *) (rinner_(i),i=1,nang*2-1) - if ( iver(1) .ge. 4 ) then - read (iunit, *) irel - else + IF (10*iver(1)+iver(2)>=51) & + READ (iunit, *) (rinner_(i),i=1,nang*2-1) + IF ( iver(1) >= 4 ) THEN + READ (iunit, *) irel + ELSE irel = 0 - end if + ENDIF - allocate(rc(nang)) - read(iunit, *) ( rc(i), i=1,nang) + ALLOCATE(rc(nang)) + READ(iunit, *) ( rc(i), i=1,nang) - read (iunit,* ) nbeta_, kkbeta + READ (iunit,* ) nbeta_, kkbeta - allocate(beta(kkbeta,nbeta_)) - allocate(qfunc_(kkbeta,nbeta_,nbeta_)) - allocate(ddd0(nbeta_,nbeta_)) - allocate(ddd (nbeta_,nbeta_)) - allocate(qqq_(nbeta_,nbeta_)) - allocate(lll_(nbeta_)) - allocate(eee (nbeta_)) - allocate(qfcoef_(nqf_,nlc,nbeta_,nbeta_)) + ALLOCATE(beta(kkbeta,nbeta_)) + ALLOCATE(qfunc_(kkbeta,nbeta_,nbeta_)) + ALLOCATE(ddd0(nbeta_,nbeta_)) + ALLOCATE(ddd (nbeta_,nbeta_)) + ALLOCATE(qqq_(nbeta_,nbeta_)) + ALLOCATE(lll_(nbeta_)) + ALLOCATE(eee (nbeta_)) + ALLOCATE(qfcoef_(nqf_,nlc,nbeta_,nbeta_)) - do j=1,nbeta_ - read ( iunit, *) lll_(j) - read ( iunit, *) eee(j), ( beta(i,j), i=1,kkbeta ) - do k=j,nbeta_ - read( iunit, *) ddd0(j,k), ddd(j,k), qqq_(j,k), & + DO j=1,nbeta_ + READ ( iunit, *) lll_(j) + READ ( iunit, *) eee(j), ( beta(i,j), i=1,kkbeta ) + DO k=j,nbeta_ + READ( iunit, *) ddd0(j,k), ddd(j,k), qqq_(j,k), & (qfunc_(i,j,k),i=1,kkbeta),& ((qfcoef_(i,lp,j,k),i=1,nqf_),lp=1,2*nang-1) - enddo - enddo + ENDDO + ENDDO - allocate(iptype(nbeta_)) - if (10*iver(1)+iver(2).ge.72) then - read ( iunit, * ) (iptype(i), i=1,nbeta_) - read ( iunit, * ) npf, ptryc - end if + ALLOCATE(iptype(nbeta_)) + IF (10*iver(1)+iver(2)>=72) THEN + READ ( iunit, * ) (iptype(i), i=1,nbeta_) + READ ( iunit, * ) npf, ptryc + ENDIF - allocate(vloc_(mesh_)) - read(iunit, *) rcloc_, ( vloc_(i), i=1,mesh_) - - allocate(rho_atc_(mesh_)) - if ( ifpcor.gt.0 ) then - read(iunit, *) rpcor - read(iunit, *) ( rho_atc_(i), i=1,mesh_) - endif + ALLOCATE(vloc_(mesh_)) + READ(iunit, *) rcloc_, ( vloc_(i), i=1,mesh_) - allocate(rho_at_(mesh_), vloc(mesh_)) - read(iunit, *) (vloc(i), i=1,mesh_) - read(iunit, *) (rho_at_(i), i=1,mesh_) + ALLOCATE(rho_atc_(mesh_)) + IF ( ifpcor>0 ) THEN + READ(iunit, *) rpcor + READ(iunit, *) ( rho_atc_(i), i=1,mesh_) + ENDIF - allocate(r_(mesh_),rab_(mesh_)) - read(iunit, *) (r_(i), i=1,mesh_) - read(iunit, *) (rab_(i),i=1,mesh_) + ALLOCATE(rho_at_(mesh_), vloc(mesh_)) + READ(iunit, *) (vloc(i), i=1,mesh_) + READ(iunit, *) (rho_at_(i), i=1,mesh_) - if (iver(1) .ge. 6) then + ALLOCATE(r_(mesh_),rab_(mesh_)) + READ(iunit, *) (r_(i), i=1,mesh_) + READ(iunit, *) (rab_(i),i=1,mesh_) + + IF (iver(1) >= 6) THEN nchi = nvalps - if (iver(1) .ge. 7) read (iunit, *) nchi - allocate(wf(mesh_,nchi)) - read (iunit, *) ((wf(i,j), i=1,mesh_),j=1,nchi) - end if + IF (iver(1) >= 7) READ (iunit, *) nchi + ALLOCATE(wf(mesh_,nchi)) + READ (iunit, *) ((wf(i,j), i=1,mesh_),j=1,nchi) + ENDIF - return -end subroutine read_vdb + RETURN +END SUBROUTINE read_vdb -subroutine convert_uspp +SUBROUTINE convert_uspp ! ---------------------------------------------------------- ! - use Vanderbilt - use constants, only : fpi - use upf - implicit none - integer i - character(len=1), dimension(0:3) :: convel=(/'S','P','D','F'/) + USE Vanderbilt + USE constants, ONLY : fpi + USE upf + IMPLICIT NONE + INTEGER i + CHARACTER(len=1), DIMENSION(0:3) :: convel=(/'S','P','D','F'/) - write(generated, '("Generated using Vanderbilt code, version ",3i3)') iver - write(date_author,'("Author: unknown Generation date:",3i5)') idmy - write(comment,'("Automatically converted from original format")') - if (irel == 0) then + WRITE(generated, '("Generated using Vanderbilt code, version ",3i3)') iver + WRITE(date_author,'("Author: unknown Generation date:",3i5)') idmy + WRITE(comment,'("Automatically converted from original format")') + IF (irel == 0) THEN rel = 0 - else if (irel == 1) then + ELSEIF (irel == 1) THEN rel = 2 - else if (irel == 2) then + ELSEIF (irel == 2) THEN rel = 1 - end if + ENDIF rcloc = rcloc_ nwfs = nvalps - allocate( els(nwfs), oc(nwfs), epseu(nwfs)) - allocate(lchi(nwfs), nns(nwfs) ) - allocate(rcut (nwfs), rcutus (nwfs)) - do i=1, nwfs + ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs)) + ALLOCATE(lchi(nwfs), nns(nwfs) ) + ALLOCATE(rcut (nwfs), rcutus (nwfs)) + DO i=1, nwfs nns (i) = nnlzps(i)/100 lchi(i) = mod (nnlzps(i)/10,10) rcut(i) = rinner1_ - rcutus(i)= rc(lchi(i)+1) + rcutus(i)= rc(lchi(i)+1) oc (i) = wwnlps(i) - write(els(i),'(i1,a1)') nns(i), convel(lchi(i)) + WRITE(els(i),'(i1,a1)') nns(i), convel(lchi(i)) epseu(i) = eeps(i) - end do - deallocate (nnlzps, rc, wwnlps, eeps) + ENDDO + DEALLOCATE (nnlzps, rc, wwnlps, eeps) psd = title - if (keyps.le.2) then + IF (keyps<=2) THEN pseudotype = 'NC' - else + ELSE pseudotype = 'US' - end if - nlcc = ifpcor.gt.0 + ENDIF + nlcc = ifpcor>0 zp = zp_ etotps = etot ecutrho=0.0d0 @@ -280,87 +280,87 @@ subroutine convert_uspp lmax = nang - 1 mesh = mesh_ nbeta = nbeta_ - if (nvalps .ne. nchi) then - print *, 'WARNING: verify info on atomic wavefunctions' - end if + IF (nvalps /= nchi) THEN + PRINT *, 'WARNING: verify info on atomic wavefunctions' + ENDIF ntwfc = nchi - allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) - do i=1, min(ntwfc,nwfs) + ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) ) + DO i=1, min(ntwfc,nwfs) elsw(i) = els(i) ocw(i) = oc (i) lchiw(i)=lchi(i) - end do - if ( exfact.eq.0) then + ENDDO + IF ( exfact==0) THEN iexch=1; icorr=1; igcx=0; igcc=0 ! Perdew-Zunger - else if ( exfact.eq.1) then + ELSEIF ( exfact==1) THEN iexch=1; icorr=3; igcx=1; igcc=3 ! Becke-Lee-Yang-Parr - else if ( exfact.eq.2) then + ELSEIF ( exfact==2) THEN iexch=1; icorr=1; igcx=1; igcc=0 ! Becke88 exchange - else if (exfact.eq.-5.or.exfact.eq.3) then + ELSEIF (exfact==-5.or.exfact==3) THEN iexch=1; icorr=1; igcx=1; igcc=1 ! Becke88-Perdew 86 - else if (exfact.eq.-6.or.exfact.eq.4) then + ELSEIF (exfact==-6.or.exfact==4) THEN iexch=1; icorr=4; igcx=2; igcc=2 ! Perdew-Wang 91 - else if (exfact.eq. 5) then + ELSEIF (exfact== 5) THEN iexch=1; icorr=4; igcx=3; igcc=4 ! Perdew-Becke-Erkerhof - else - write (6,'("convert: wrong xc in pseudopotential ",f12.6)') exfact - stop - end if + ELSE + WRITE (6,'("convert: wrong xc in pseudopotential ",f12.6)') exfact + STOP + ENDIF - allocate (r(mesh), rab(mesh)) + ALLOCATE (r(mesh), rab(mesh)) r = r_ rab=rab_ - deallocate (r_, rab_) + DEALLOCATE (r_, rab_) - allocate (rho_atc(mesh)) + ALLOCATE (rho_atc(mesh)) ! Vanderbilt rho_core(r) = 4pi*r^2*rho_core(r) UPF rho_atc (1) = 0.d0 rho_atc (2:mesh) = rho_atc_(2:mesh) / fpi / r(2:mesh)**2 - deallocate (rho_atc_) + DEALLOCATE (rho_atc_) - allocate (vloc0(mesh)) + ALLOCATE (vloc0(mesh)) vloc0(2:mesh) = vloc_(2:mesh)/r(2:mesh) vloc0(1) = vloc0(2) - deallocate (vloc_) + DEALLOCATE (vloc_) - allocate(ikk2(nbeta), lll(nbeta)) + ALLOCATE(ikk2(nbeta), lll(nbeta)) ikk2 = kkbeta lll = lll_ - deallocate (lll_) - allocate(betar(kkbeta,nbeta)) + DEALLOCATE (lll_) + ALLOCATE(betar(kkbeta,nbeta)) betar = beta - deallocate (beta) + DEALLOCATE (beta) - allocate(dion(nbeta,nbeta)) + ALLOCATE(dion(nbeta,nbeta)) dion = ddd0 - deallocate (ddd0) + DEALLOCATE (ddd0) - allocate(qqq(nbeta,nbeta)) + ALLOCATE(qqq(nbeta,nbeta)) qqq = qqq_ - deallocate (qqq_) + DEALLOCATE (qqq_) - allocate(qfunc(mesh,nbeta,nbeta)) + ALLOCATE(qfunc(mesh,nbeta,nbeta)) qfunc(1:kkbeta,:,:) = qfunc_(1:kkbeta,:,:) qfunc(kkbeta+1:mesh,:,:) = 0.d0 - deallocate (qfunc_) + DEALLOCATE (qfunc_) nqf = nqf_ nqlc= nlc - allocate(rinner(nqlc)) + ALLOCATE(rinner(nqlc)) rinner = rinner_ - deallocate(rinner_) - allocate(qfcoef(nqf,nqlc,nbeta,nbeta)) + DEALLOCATE(rinner_) + ALLOCATE(qfcoef(nqf,nqlc,nbeta,nbeta)) qfcoef = qfcoef_ - deallocate (qfcoef_) + DEALLOCATE (qfcoef_) - allocate (rho_at(mesh)) + ALLOCATE (rho_at(mesh)) rho_at = rho_at_ - deallocate (rho_at_) + DEALLOCATE (rho_at_) - allocate (chi(mesh,ntwfc)) + ALLOCATE (chi(mesh,ntwfc)) chi = wf - deallocate (wf) + DEALLOCATE (wf) - return -end subroutine convert_uspp + RETURN +END SUBROUTINE convert_uspp diff --git a/upftools/vdb2upf.f90 b/upftools/vdb2upf.f90 index dea3b0908..58480423a 100644 --- a/upftools/vdb2upf.f90 +++ b/upftools/vdb2upf.f90 @@ -7,32 +7,32 @@ ! ! !--------------------------------------------------------------------- -program vdb2upf +PROGRAM vdb2upf !--------------------------------------------------------------------- ! ! Convert a pseudopotential written in Vanderbilt format ! (formatted) to unified pseudopotential format ! - implicit none - character(len=256) filein, fileout + IMPLICIT NONE + CHARACTER(len=256) filein, fileout ! ! - call get_file ( filein ) - open(unit=1,file=filein,status='old',form='formatted') - call read_vdb(1) - close (unit=1) + CALL get_file ( filein ) + OPEN(unit=1,file=filein,status='old',form='formatted') + CALL read_vdb(1) + CLOSE (unit=1) ! convert variables read from Vanderbilt format into those needed ! by the upf format - add missing quantities - call convert_uspp + CALL convert_uspp fileout=trim(filein)//'.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) - stop -end program vdb2upf + STOP +END PROGRAM vdb2upf diff --git a/upftools/virtual.f90 b/upftools/virtual.f90 index ffe4e8ec8..c56589d2c 100644 --- a/upftools/virtual.f90 +++ b/upftools/virtual.f90 @@ -6,32 +6,32 @@ ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! -module pseudo +MODULE pseudo ! ! All variables to be read from the UPF file ! (UPF = unified pseudopotential format) ! - integer ,parameter :: npsx = 2 + INTEGER ,PARAMETER :: npsx = 2 ! npsx : maximum number of different pseudopotentials - integer, parameter :: lmaxx = 3, nchix = 6, ndm = 2000 - ! lmaxx : maximum non local angular momentum in PP + INTEGER, PARAMETER :: lmaxx = 3, nchix = 6, ndm = 2000 + ! lmaxx : maximum non local angular momentum in PP ! nchix : maximum number of atomic wavefunctions per PP ! ndm : maximum number of points in the radial mesh - integer, parameter :: nbrx = 8, lqmax = 5, nqfx = 8 - ! nbrx : maximum number of beta functions - ! lqmax : maximum number of angular momentum of Q + INTEGER, PARAMETER :: nbrx = 8, lqmax = 5, nqfx = 8 + ! nbrx : maximum number of beta functions + ! lqmax : maximum number of angular momentum of Q ! nqfx : maximum number of coefficients in Q smoothing ! ! pp_header - character (len=80):: generated, date_author, comment - character (len=2) :: psd(npsx), pseudotype - character (len=20):: dft(npsx) - integer :: lmax(npsx), mesh(npsx), nbeta(npsx), ntwfc(npsx) - logical :: nlcc(npsx), isus(npsx) + CHARACTER (len=80):: generated, date_author, comment + CHARACTER (len=2) :: psd(npsx), pseudotype + CHARACTER (len=20):: dft(npsx) + INTEGER :: lmax(npsx), mesh(npsx), nbeta(npsx), ntwfc(npsx) + LOGICAL :: nlcc(npsx), isus(npsx) real(8) :: zp(npsx), ecutrho, ecutwfc, etotps real(8) :: oc(nchix,npsx) - character(len=2) :: els(nchix,npsx) - integer :: lchi(nchix,npsx) + CHARACTER(len=2) :: els(nchix,npsx) + INTEGER :: lchi(nchix,npsx) ! ! pp_mesh real(8) :: r(ndm,npsx), rab(ndm,npsx) @@ -44,11 +44,11 @@ module pseudo ! pp_nonlocal ! pp_beta real(8) :: betar(ndm, nbrx, npsx) - integer :: lll(nbrx,npsx), ikk2(nbrx,npsx) + INTEGER :: lll(nbrx,npsx), ikk2(nbrx,npsx) ! pp_dij real(8) :: dion(nbrx,nbrx,npsx) ! pp_qij - integer :: nqf(npsx), nqlc(npsx) + INTEGER :: nqf(npsx), nqlc(npsx) real(8) :: rinner(lqmax,npsx), qqq(nbrx,nbrx,npsx), & qfunc(ndm,nbrx,nbrx,npsx) ! pp_qfcoef @@ -59,54 +59,54 @@ module pseudo ! ! pp_rhoatom real(8) :: rho_at(ndm,npsx) -end module pseudo +END MODULE pseudo ! -program virtual +PROGRAM virtual !--------------------------------------------------------------------- ! ! Read pseudopotentials in the Unified Pseudopotential Format (UPF) ! - implicit none - integer :: is, ios, iunps = 4 + IMPLICIT NONE + INTEGER :: is, ios, iunps = 4 real (8) :: x - character (len=256) :: filein(2), fileout - print '('' '')' - print '('' Generate the UPF pseudopotential for a virtual atom '')' - print '('' combining two pseudopootentials in UPF format '')' - print '('' '')' + CHARACTER (len=256) :: filein(2), fileout + PRINT '('' '')' + PRINT '('' Generate the UPF pseudopotential for a virtual atom '')' + PRINT '('' combining two pseudopootentials in UPF format '')' + PRINT '('' '')' ! - do is=1,2 - print '('' Input PP file # '',i2,'' in UPF format > '',$)', is - read (5, '(a)', end = 20, err = 20) filein(is) - open(unit=iunps,file=filein(is),status='old',form='formatted',iostat=ios) - if (ios.ne.0) stop - write (*,*) " IOS= ", ios, is, iunps - call read_pseudo(is, iunps) - close (unit=iunps) - print '('' '')' - end do - print '('' New Pseudo = x '',a,'' + (1-x) '',a)', (trim(filein(is)), is=1,2) -10 continue - print '('' mixing parameter x [01) go to 10 + DO is=1,2 + PRINT '('' Input PP file # '',i2,'' in UPF format > '',$)', is + READ (5, '(a)', end = 20, err = 20) filein(is) + OPEN(unit=iunps,file=filein(is),status='old',form='formatted',iostat=ios) + IF (ios/=0) STOP + WRITE (*,*) " IOS= ", ios, is, iunps + CALL read_pseudo(is, iunps) + CLOSE (unit=iunps) + PRINT '('' '')' + ENDDO + PRINT '('' New Pseudo = x '',a,'' + (1-x) '',a)', (trim(filein(is)), is=1,2) +10 CONTINUE + PRINT '('' mixing parameter x [01) GOTO 10 - call compute_virtual(x,filein) + CALL compute_virtual(x,filein) fileout='NewPseudo.UPF' - print '(''Output PP file in UPF format : '',a)', fileout + PRINT '(''Output PP file in UPF format : '',a)', fileout - open(unit=2,file=fileout,status='unknown',form='formatted') - call write_upf(2) - close (unit=2) + OPEN(unit=2,file=fileout,status='unknown',form='formatted') + CALL write_upf(2) + CLOSE (unit=2) -20 stop -end program virtual +20 STOP +END PROGRAM virtual ! !--------------------------------------------------------------------- -subroutine compute_virtual(x,filein) - use pseudo - use upf, ONLY : & +SUBROUTINE compute_virtual(x,filein) + USE pseudo + USE upf, ONLY : & upf_rel => rel, upf_rcloc => rcloc, upf_nwfs => nwfs, & upf_oc => oc, upf_rcut => rcut, upf_rcutus => rcutus, & upf_epseu => epseu, upf_els => els, & @@ -133,15 +133,15 @@ subroutine compute_virtual(x,filein) upf_qfcoef => qfcoef, & upf_chi => chi, & upf_rho_at => rho_at - use splinelib - use funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc - implicit none - integer :: i, j, ib - character (len=256) :: filein(2) - character (len=5) :: xlabel + USE splinelib + USE funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc + IMPLICIT NONE + INTEGER :: i, j, ib + CHARACTER (len=256) :: filein(2) + CHARACTER (len=5) :: xlabel real (8) :: x, capel - real (8), allocatable :: aux1(:,:), aux2(:,:) - logical :: interpolate + real (8), ALLOCATABLE :: aux1(:,:), aux2(:,:) + LOGICAL :: interpolate interpolate = .false. ! !pp_info @@ -152,29 +152,29 @@ subroutine compute_virtual(x,filein) upf_generated = 'Generated using virtual.x code ' upf_date_author= 'Author and generation date: unknown. '//& 'Refer to original pseudopotential files' - write( xlabel, '(f5.3)' ) x + WRITE( xlabel, '(f5.3)' ) x upf_comment = 'Pseudo = x '//trim(filein(1))//& ' + (1-x) '//trim(filein(2))//', with x='//xlabel upf_psd = "Xx" upf_pseudotype = "NC" - if (isus(1) .or. isus(2)) upf_pseudotype = "US" - call set_dft_from_name(dft(1)) + IF (isus(1) .or. isus(2)) upf_pseudotype = "US" + CALL set_dft_from_name(dft(1)) upf_iexch = get_iexch() upf_icorr = get_icorr() upf_igcx = get_igcx() upf_igcc = get_igcc() - call set_dft_from_name(dft(2)) - if (get_iexch().ne.upf_iexch .or. get_icorr().ne.upf_icorr .or. & - get_igcx().ne.upf_igcx .or. get_igcc().ne.upf_igcc) & - call errore ('virtual','conflicting DFT functionals',1) + CALL set_dft_from_name(dft(2)) + IF (get_iexch()/=upf_iexch .or. get_icorr()/=upf_icorr .or. & + get_igcx()/=upf_igcx .or. get_igcc()/=upf_igcc) & + CALL errore ('virtual','conflicting DFT functionals',1) upf_lmax = max(lmax(1), lmax(2)) - if (mesh(1).ne.mesh(2) ) then - write (*,*) " pseudopotentials have different mesh " - write (*,*) mesh(1),mesh(2) - write (*,*) r(1,1), r(1,2) - write (*,*) r(mesh(1),1),r(mesh(2),2) + IF (mesh(1)/=mesh(2) ) THEN + WRITE (*,*) " pseudopotentials have different mesh " + WRITE (*,*) mesh(1),mesh(2) + WRITE (*,*) r(1,1), r(1,2) + WRITE (*,*) r(mesh(1),1),r(mesh(2),2) interpolate = .true. - end if + ENDIF upf_mesh = mesh(1) upf_nbeta = nbeta(1)+nbeta(2) upf_ntwfc = ntwfc(1) @@ -182,7 +182,7 @@ subroutine compute_virtual(x,filein) upf_ecutrho = ecutrho upf_ecutwfc = ecutwfc upf_etotps = etotps - allocate( upf_ocw(upf_ntwfc), upf_elsw(upf_ntwfc), upf_lchiw(upf_ntwfc) ) + ALLOCATE( upf_ocw(upf_ntwfc), upf_elsw(upf_ntwfc), upf_lchiw(upf_ntwfc) ) upf_ocw(1:upf_ntwfc) = oc(1:upf_ntwfc,1) upf_elsw(1:upf_ntwfc) = els(1:upf_ntwfc,1) upf_lchiw(1:upf_ntwfc) = lchi(1:upf_ntwfc,1) @@ -190,513 +190,513 @@ subroutine compute_virtual(x,filein) ! !pp_mesh capel = 0.d0 - do i=1,upf_mesh + DO i=1,upf_mesh capel = capel + abs(r(i,1)-r(i,2)) + abs(rab(i,1)-rab(i,2)) - end do - if (capel.gt.1.d-6) then - write (*,*) " pseudopotentials have different mesh " + ENDDO + IF (capel>1.d-6) THEN + WRITE (*,*) " pseudopotentials have different mesh " interpolate = .true. - end if - write (*,*) "INTERPOLATE =", interpolate + ENDIF + WRITE (*,*) "INTERPOLATE =", interpolate !if (interpolate) call errore ("virtual", & ! "grid interpolation is not working yet",1) - if (interpolate) allocate ( aux1(1,mesh(1)), aux2(1,mesh(2)) ) + IF (interpolate) ALLOCATE ( aux1(1,mesh(1)), aux2(1,mesh(2)) ) - allocate( upf_r(upf_mesh), upf_rab(upf_mesh) ) + ALLOCATE( upf_r(upf_mesh), upf_rab(upf_mesh) ) upf_r(1:upf_mesh) = r(1:upf_mesh,1) upf_rab(1:upf_mesh) = rab(1:upf_mesh,1) ! !pp_nlcc - allocate( upf_rho_atc(upf_mesh) ) - if (interpolate) then - write (*,*) "interpolate rho_atc" + ALLOCATE( upf_rho_atc(upf_mesh) ) + IF (interpolate) THEN + WRITE (*,*) "interpolate rho_atc" aux2(1,1:mesh(2)) = rho_atc(1:mesh(2),2) - call dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) + CALL dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) rho_atc(1:upf_mesh,2) = aux1(1,1:upf_mesh) - write (*,*) " done" - end if + WRITE (*,*) " done" + ENDIF upf_rho_atc(1:upf_mesh) = x * rho_atc(1:upf_mesh,1) + & (1.d0-x) * rho_atc(1:upf_mesh,2) ! !pp_local - allocate( upf_vloc0(upf_mesh) ) - if (interpolate) then - write (*,*) " interpolate vloc0" + ALLOCATE( upf_vloc0(upf_mesh) ) + IF (interpolate) THEN + WRITE (*,*) " interpolate vloc0" aux2(1,1:mesh(2)) = vloc0(1:mesh(2),2) - - call dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) - + + CALL dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) + vloc0(1:upf_mesh,2) = aux1(1,1:upf_mesh) - - ! Jivtesh - if the mesh of the first atom extends to a larger radius + + ! Jivtesh - if the mesh of the first atom extends to a larger radius ! than the mesh of the second atom, then, for those radii that are ! greater than the maximum radius of the second atom, the local potential - ! of the second atom is calculated using the expression - ! v_local = (-2)*Z/r instead of using the extrapolated value. - ! This is because, typically extrapolation leads to positive potentials. + ! of the second atom is calculated using the expression + ! v_local = (-2)*Z/r instead of using the extrapolated value. + ! This is because, typically extrapolation leads to positive potentials. ! This is implemented in lines 240-242 - - do i=1,mesh(1) - if(r(i,1).GT.r(mesh(2),2)) vloc0(i,2) = -(2.0*zp(2))/r(i,1) - end do - - end if + + DO i=1,mesh(1) + IF(r(i,1)>r(mesh(2),2)) vloc0(i,2) = -(2.0*zp(2))/r(i,1) + ENDDO + + ENDIF upf_vloc0(1:upf_mesh) = x * vloc0(1:upf_mesh,1) + & (1.d0-x) * vloc0(1:upf_mesh,2) ! !pp_nonlocal !pp_beta - allocate( upf_betar(upf_mesh,upf_nbeta), & + ALLOCATE( upf_betar(upf_mesh,upf_nbeta), & upf_lll(upf_nbeta), upf_ikk2(upf_nbeta) ) ib = 0 - do i=1,nbeta(1) + DO i=1,nbeta(1) ib = ib + 1 upf_betar(1:upf_mesh,ib) = betar(1:upf_mesh,i,1) upf_lll(ib) = lll(i,1) upf_ikk2(ib) = ikk2(i,1) - end do - do i=1,nbeta(2) + ENDDO + DO i=1,nbeta(2) ib = ib + 1 - if (interpolate) then - write (*,*) " interpolate betar" + IF (interpolate) THEN + WRITE (*,*) " interpolate betar" aux2(1,1:mesh(2)) = betar(1:mesh(2),i,2) - call dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) + CALL dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) betar(1:upf_mesh,i,2) = aux1(1,1:upf_mesh) - end if + ENDIF upf_betar(1:upf_mesh,ib) = betar(1:upf_mesh,i,2) upf_lll(ib) = lll(i,2) - ! SdG - when the meshes of the two pseudo are different the ikk2 limits + ! SdG - when the meshes of the two pseudo are different the ikk2 limits ! for the beta functions of the second one must be set properly ! This is done in lines 273-277 - if (interpolate) then + IF (interpolate) THEN j = 1 - do while ( upf_r(j) < r( ikk2(i,2), 2) ) + DO WHILE ( upf_r(j) < r( ikk2(i,2), 2) ) j = j + 1 - end do + ENDDO upf_ikk2(ib) = j - else + ELSE upf_ikk2(ib) = ikk2(i,2) - end if - end do + ENDIF + ENDDO ! !pp_dij - allocate( upf_dion(upf_nbeta, upf_nbeta) ) + ALLOCATE( upf_dion(upf_nbeta, upf_nbeta) ) upf_dion(:,:) = 0.d0 - do i=1,nbeta(1) - do j=1,nbeta(1) + DO i=1,nbeta(1) + DO j=1,nbeta(1) upf_dion(i,j) = x * dion(i,j,1) - end do - end do - do i=1,nbeta(2) - do j=1,nbeta(2) + ENDDO + ENDDO + DO i=1,nbeta(2) + DO j=1,nbeta(2) upf_dion(nbeta(1)+i,nbeta(1)+j) = (1.d0-x) * dion(i,j,2) - end do - end do + ENDDO + ENDDO ! !pp_qij - if (nqf(1).ne.nqf(2)) & - call errore ("Virtual","different nqf are not implemented (yet)", 1) - if (nqlc(1).ne.nqlc(2)) & - call errore ("Virtual","different nqlc are not implemented (yet)", 1) + IF (nqf(1)/=nqf(2)) & + CALL errore ("Virtual","different nqf are not implemented (yet)", 1) + IF (nqlc(1)/=nqlc(2)) & + CALL errore ("Virtual","different nqlc are not implemented (yet)", 1) upf_nqf = nqf(1) upf_nqlc = nqlc(1) - allocate( upf_rinner(upf_nqlc), upf_qqq(upf_nbeta,upf_nbeta), & + ALLOCATE( upf_rinner(upf_nqlc), upf_qqq(upf_nbeta,upf_nbeta), & upf_qfunc(upf_mesh,upf_nbeta,upf_nbeta) ) - do i=1,upf_nqlc - if(rinner(i,1).ne.rinner(i,2)) & - call errore("Virtual","different rinner are not implemented (yet)",i) - end do + DO i=1,upf_nqlc + IF(rinner(i,1)/=rinner(i,2)) & + CALL errore("Virtual","different rinner are not implemented (yet)",i) + ENDDO upf_rinner(1:upf_nqlc) = rinner(1:upf_nqlc,1) - + upf_qqq(:,:) = 0.d0 upf_qfunc(:,:,:) = 0.d0 - do i=1,nbeta(1) - do j=1,nbeta(1) + DO i=1,nbeta(1) + DO j=1,nbeta(1) upf_qqq(i,j) = x * qqq(i, j,1) upf_qfunc(1:upf_mesh,i,j) = x * qfunc(1:upf_mesh,i,j,1) - end do - end do - do i=1,nbeta(2) - do j=1,nbeta(2) + ENDDO + ENDDO + DO i=1,nbeta(2) + DO j=1,nbeta(2) upf_qqq(nbeta(1)+i,nbeta(1)+j) = (1.d0-x) * qqq(i, j, 2) - if (interpolate) then - write (*,*) " interpolate qfunc" + IF (interpolate) THEN + WRITE (*,*) " interpolate qfunc" aux2(1,1:mesh(2) ) = qfunc(1:mesh(2),i,j,2) - call dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) + CALL dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) qfunc(1:upf_mesh,i,j,2) = aux1(1,1:upf_mesh) - write (*,*) " done" - end if + WRITE (*,*) " done" + ENDIF upf_qfunc(1:upf_mesh,nbeta(1)+i,nbeta(1)+j) = (1.d0-x) * qfunc(1:upf_mesh,i,j,2) - end do - end do + ENDDO + ENDDO ! !pp_qfcoef - allocate( upf_qfcoef(upf_nqf,upf_nqlc,upf_nbeta,upf_nbeta) ) + ALLOCATE( upf_qfcoef(upf_nqf,upf_nqlc,upf_nbeta,upf_nbeta) ) upf_qfcoef(:,:,:,:) = 0.d0 - do i=1,nbeta(1) - do j=1,nbeta(1) + DO i=1,nbeta(1) + DO j=1,nbeta(1) upf_qfcoef(1:upf_nqf,1:upf_nqlc,i,j) = & x * qfcoef(1:upf_nqf,1:upf_nqlc,i,j, 1) - end do - end do - do i=1,nbeta(2) - do j=1,nbeta(2) + ENDDO + ENDDO + DO i=1,nbeta(2) + DO j=1,nbeta(2) upf_qfcoef(1:upf_nqf,1:upf_nqlc,nbeta(1)+i,nbeta(1)+j) = & (1.d0-x) * qfcoef(1:upf_nqf,1:upf_nqlc,i,j, 2) - end do - end do + ENDDO + ENDDO ! !pp_pswfc - - allocate (upf_chi(upf_mesh,upf_ntwfc) ) - - if (ntwfc(1).EQ.ntwfc(2)) then - - do i=1,ntwfc(2) - if (interpolate) then - write (*,*) " interpolate chi" + + ALLOCATE (upf_chi(upf_mesh,upf_ntwfc) ) + + IF (ntwfc(1)==ntwfc(2)) THEN + + DO i=1,ntwfc(2) + IF (interpolate) THEN + WRITE (*,*) " interpolate chi" aux2(1,1:mesh(2)) = chi(1:mesh(2),i,2) - call dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) + CALL dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) chi(1:upf_mesh,i,2) = aux1(1,1:upf_mesh) - write (*,*) " done" - end if - ! Jivtesh - The wavefunctions are calcuated to be the average of the + WRITE (*,*) " done" + ENDIF + ! Jivtesh - The wavefunctions are calcuated to be the average of the ! wavefunctions of the two atoms - lines 365-366 upf_chi(1:upf_mesh,i) = x * chi(1:upf_mesh,i,1) + & - (1.d0-x) * chi(1:upf_mesh,i,2) - enddo - else - write (*,*) "Number of wavefunctions not the same for the two pseudopotentials" - endif - !upf_chi(1:upf_mesh,1:upf_ntwfc) = chi(1:upf_mesh,1:upf_ntwfc,1) + (1.d0-x) * chi(1:upf_mesh,i,2) + ENDDO + ELSE + WRITE (*,*) "Number of wavefunctions not the same for the two pseudopotentials" + ENDIF + !upf_chi(1:upf_mesh,1:upf_ntwfc) = chi(1:upf_mesh,1:upf_ntwfc,1) ! !pp_rhoatm - - allocate (upf_rho_at(upf_mesh) ) - if (interpolate) then - write (*,*) " interpolate rho_at" - aux2(1,1:mesh(2)) = rho_at(1:mesh(2),2) - call dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) - rho_at(1:upf_mesh,2) = aux1(1,1:upf_mesh) - write (*,*) " done" - end if - upf_rho_at(1:upf_mesh) = x * rho_at(1:upf_mesh,1) + & - (1.d0-x) * rho_at(1:upf_mesh,2) -end subroutine compute_virtual + ALLOCATE (upf_rho_at(upf_mesh) ) + IF (interpolate) THEN + WRITE (*,*) " interpolate rho_at" + aux2(1,1:mesh(2)) = rho_at(1:mesh(2),2) + CALL dosplineint( r(1:mesh(2),2), aux2, upf_r(1:upf_mesh), aux1 ) + rho_at(1:upf_mesh,2) = aux1(1,1:upf_mesh) + WRITE (*,*) " done" + ENDIF + upf_rho_at(1:upf_mesh) = x * rho_at(1:upf_mesh,1) + & + (1.d0-x) * rho_at(1:upf_mesh,2) + +END SUBROUTINE compute_virtual ! !--------------------------------------------------------------------- -subroutine read_pseudo (is, iunps) +SUBROUTINE read_pseudo (is, iunps) !--------------------------------------------------------------------- ! ! Read pseudopotential in the Unified Pseudopotential Format (UPF) ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! is : index of this pseudopotential ! iunps: unit connected with pseudopotential file ! - if (is < 0 .or. is > npsx ) call errore ('read_pseudo', 'Wrong is number', 1) - write ( *, * ) " Reading pseudopotential file in UPF format..." + IF (is < 0 .or. is > npsx ) CALL errore ('read_pseudo', 'Wrong is number', 1) + WRITE ( *, * ) " Reading pseudopotential file in UPF format..." !------->Search for Header - call scan_begin (iunps, "HEADER", .true.) - call read_pseudo_header (is, iunps) - call scan_end (iunps, "HEADER") + CALL scan_begin (iunps, "HEADER", .true.) + CALL read_pseudo_header (is, iunps) + CALL scan_end (iunps, "HEADER") !-------->Search for mesh information - call scan_begin (iunps, "MESH", .true.) - call read_pseudo_mesh (is, iunps) - call scan_end (iunps, "MESH") + CALL scan_begin (iunps, "MESH", .true.) + CALL read_pseudo_mesh (is, iunps) + CALL scan_end (iunps, "MESH") !-------->If present, search for nlcc - if (nlcc (is) ) then - call scan_begin (iunps, "NLCC", .true.) - call read_pseudo_nlcc (is, iunps) - call scan_end (iunps, "NLCC") - endif + IF (nlcc (is) ) THEN + CALL scan_begin (iunps, "NLCC", .true.) + CALL read_pseudo_nlcc (is, iunps) + CALL scan_end (iunps, "NLCC") + ENDIF !-------->Search for Local potential - call scan_begin (iunps, "LOCAL", .true.) - call read_pseudo_local (is, iunps) - call scan_end (iunps, "LOCAL") + CALL scan_begin (iunps, "LOCAL", .true.) + CALL read_pseudo_local (is, iunps) + CALL scan_end (iunps, "LOCAL") !-------->Search for Nonlocal potential - call scan_begin (iunps, "NONLOCAL", .true.) - call read_pseudo_nl (is, iunps) - call scan_end (iunps, "NONLOCAL") + CALL scan_begin (iunps, "NONLOCAL", .true.) + CALL read_pseudo_nl (is, iunps) + CALL scan_end (iunps, "NONLOCAL") !-------->Search for atomic wavefunctions - call scan_begin (iunps, "PSWFC", .true.) - call read_pseudo_pswfc (is, iunps) - call scan_end (iunps, "PSWFC") + CALL scan_begin (iunps, "PSWFC", .true.) + CALL read_pseudo_pswfc (is, iunps) + CALL scan_end (iunps, "PSWFC") !-------->Search for atomic charge - call scan_begin (iunps, "RHOATOM", .true.) - call read_pseudo_rhoatom (is, iunps) - call scan_end (iunps, "RHOATOM") + CALL scan_begin (iunps, "RHOATOM", .true.) + CALL read_pseudo_rhoatom (is, iunps) + CALL scan_end (iunps, "RHOATOM") ! - write ( *, * ) " ...done" - return -end subroutine read_pseudo + WRITE ( *, * ) " ...done" + RETURN +END SUBROUTINE read_pseudo !--------------------------------------------------------------------- -subroutine scan_begin (iunps, string, rew) +SUBROUTINE scan_begin (iunps, string, rew) !--------------------------------------------------------------------- ! - implicit none + IMPLICIT NONE ! Unit of the input file - integer :: iunps + INTEGER :: iunps ! Label to be matched - character (len=*) :: string - logical :: rew + CHARACTER (len=*) :: string + LOGICAL :: rew ! Flag: if .true. rewind the file - character (len=80) :: rstring + CHARACTER (len=80) :: rstring ! String read from file - integer :: ios - logical, external :: matches + INTEGER :: ios + LOGICAL, EXTERNAL :: matches ios = 0 - if (rew) rewind (iunps) - do while (ios.eq.0) - read (iunps, *, iostat = ios, err = 300) rstring - if (matches ("", rstring) ) return - enddo -300 call errore ('scan_begin', 'No '//string//' block', abs (ios) ) + IF (rew) REWIND (iunps) + DO WHILE (ios==0) + READ (iunps, *, iostat = ios, err = 300) rstring + IF (matches ("", rstring) ) RETURN + ENDDO +300 CALL errore ('scan_begin', 'No '//string//' block', abs (ios) ) -end subroutine scan_begin +END SUBROUTINE scan_begin !--------------------------------------------------------------------- -subroutine scan_end (iunps, string) +SUBROUTINE scan_end (iunps, string) !--------------------------------------------------------------------- - implicit none + IMPLICIT NONE ! Unit of the input file - integer :: iunps + INTEGER :: iunps ! Label to be matched - character (len=*) :: string + CHARACTER (len=*) :: string ! String read from file - character (len=80) :: rstring - integer :: ios - logical, external :: matches + CHARACTER (len=80) :: rstring + INTEGER :: ios + LOGICAL, EXTERNAL :: matches - read (iunps, '(a)', iostat = ios, err = 300) rstring - if (matches ("", rstring) ) return -300 call errore ('scan_end', & + READ (iunps, '(a)', iostat = ios, err = 300) rstring + IF (matches ("", rstring) ) RETURN +300 CALL errore ('scan_end', & 'No '//string//' block end statement, possibly corrupted file', - 1) -end subroutine scan_end +END SUBROUTINE scan_end ! !--------------------------------------------------------------------- -subroutine read_pseudo_header (is, iunps) +SUBROUTINE read_pseudo_header (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: nv, ios, nw - character (len=75) :: dummy - logical, external :: matches + INTEGER :: nv, ios, nw + CHARACTER (len=75) :: dummy + LOGICAL, EXTERNAL :: matches - read (iunps, *, err = 100, iostat = ios) nv, dummy - read (iunps, *, err = 100, iostat = ios) psd (is), dummy - read (iunps, *, err = 100, iostat = ios) pseudotype - if (matches (pseudotype, "US") ) isus (is) = .true. - read (iunps, *, err = 100, iostat = ios) nlcc (is), dummy - read (iunps, '(a20,t24,a)', err = 100, iostat = ios) dft(is), dummy - read (iunps, * ) zp (is), dummy - read (iunps, * ) etotps, dummy - read (iunps, * ) ecutwfc, ecutrho - read (iunps, * ) lmax (is), dummy - read (iunps, *, err = 100, iostat = ios) mesh (is), dummy - read (iunps, *, err = 100, iostat = ios) ntwfc(is), nbeta (is), dummy - read (iunps, '(a)', err = 100, iostat = ios) dummy - do nw = 1, ntwfc(is) - read (iunps, * ) els (nw,is), lchi (nw, is), oc (nw, is) - enddo - return -100 call errore ('read_pseudo_header', 'Reading pseudo file', abs (ios)) -end subroutine read_pseudo_header + READ (iunps, *, err = 100, iostat = ios) nv, dummy + READ (iunps, *, err = 100, iostat = ios) psd (is), dummy + READ (iunps, *, err = 100, iostat = ios) pseudotype + IF (matches (pseudotype, "US") ) isus (is) = .true. + READ (iunps, *, err = 100, iostat = ios) nlcc (is), dummy + READ (iunps, '(a20,t24,a)', err = 100, iostat = ios) dft(is), dummy + READ (iunps, * ) zp (is), dummy + READ (iunps, * ) etotps, dummy + READ (iunps, * ) ecutwfc, ecutrho + READ (iunps, * ) lmax (is), dummy + READ (iunps, *, err = 100, iostat = ios) mesh (is), dummy + READ (iunps, *, err = 100, iostat = ios) ntwfc(is), nbeta (is), dummy + READ (iunps, '(a)', err = 100, iostat = ios) dummy + DO nw = 1, ntwfc(is) + READ (iunps, * ) els (nw,is), lchi (nw, is), oc (nw, is) + ENDDO + RETURN +100 CALL errore ('read_pseudo_header', 'Reading pseudo file', abs (ios)) +END SUBROUTINE read_pseudo_header ! !--------------------------------------------------------------------- -subroutine read_pseudo_local (is, iunps) +SUBROUTINE read_pseudo_local (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios ! - read (iunps, *, err=100, iostat=ios) (vloc0(ir,is) , ir=1,mesh(is)) + READ (iunps, *, err=100, iostat=ios) (vloc0(ir,is) , ir=1,mesh(is)) -100 call errore ('read_pseudo_local','Reading pseudo file', abs(ios) ) +100 CALL errore ('read_pseudo_local','Reading pseudo file', abs(ios) ) - return -end subroutine read_pseudo_local + RETURN +END SUBROUTINE read_pseudo_local ! !--------------------------------------------------------------------- -subroutine read_pseudo_mesh (is, iunps) +SUBROUTINE read_pseudo_mesh (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios ! - call scan_begin (iunps, "R", .false.) - read (iunps, *, err = 100, iostat = ios) (r(ir,is), ir=1,mesh(is) ) - call scan_end (iunps, "R") - call scan_begin (iunps, "RAB", .false.) - read (iunps, *, err = 100, iostat = ios) (rab(ir,is), ir=1,mesh(is) ) - call scan_end (iunps, "RAB") + CALL scan_begin (iunps, "R", .false.) + READ (iunps, *, err = 100, iostat = ios) (r(ir,is), ir=1,mesh(is) ) + CALL scan_end (iunps, "R") + CALL scan_begin (iunps, "RAB", .false.) + READ (iunps, *, err = 100, iostat = ios) (rab(ir,is), ir=1,mesh(is) ) + CALL scan_end (iunps, "RAB") - return + RETURN -100 call errore ('read_pseudo_mesh', 'Reading pseudo file', abs (ios) ) -end subroutine read_pseudo_mesh +100 CALL errore ('read_pseudo_mesh', 'Reading pseudo file', abs (ios) ) +END SUBROUTINE read_pseudo_mesh ! !--------------------------------------------------------------------- -subroutine read_pseudo_nl (is, iunps) +SUBROUTINE read_pseudo_nl (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: nb, mb, n, ir, nd, ios, idum, ldum, icon, lp, i + INTEGER :: nb, mb, n, ir, nd, ios, idum, ldum, icon, lp, i ! counters - character (len=75) :: dummy + CHARACTER (len=75) :: dummy ! - do nb = 1, nbeta (is) - call scan_begin (iunps, "BETA", .false.) - read (iunps, *, err = 100, iostat = ios) idum, lll(nb,is), dummy - read (iunps, '(i6)', err = 100, iostat = ios) ikk2(nb,is) - read (iunps, *, err = 100, iostat = ios) & + DO nb = 1, nbeta (is) + CALL scan_begin (iunps, "BETA", .false.) + READ (iunps, *, err = 100, iostat = ios) idum, lll(nb,is), dummy + READ (iunps, '(i6)', err = 100, iostat = ios) ikk2(nb,is) + READ (iunps, *, err = 100, iostat = ios) & (betar(ir,nb,is), ir=1,ikk2(nb,is)) - do ir = ikk2(nb,is) + 1, mesh (is) - betar (ir, nb, is) = 0.d0 - enddo - call scan_end (iunps, "BETA") - enddo + DO ir = ikk2(nb,is) + 1, mesh (is) + betar (ir, nb, is) = 0.d0 + ENDDO + CALL scan_end (iunps, "BETA") + ENDDO WRITE(*,*)'ikk2',ikk2 - call scan_begin (iunps, "DIJ", .false.) - read (iunps, *, err = 100, iostat = ios) nd, dummy + CALL scan_begin (iunps, "DIJ", .false.) + READ (iunps, *, err = 100, iostat = ios) nd, dummy dion (:,:,is) = 0.d0 - do icon = 1, nd - read (iunps, *, err = 100, iostat = ios) nb, mb, dion(nb,mb,is) - dion (mb,nb,is) = dion (nb,mb,is) - enddo - call scan_end (iunps, "DIJ") + DO icon = 1, nd + READ (iunps, *, err = 100, iostat = ios) nb, mb, dion(nb,mb,is) + dion (mb,nb,is) = dion (nb,mb,is) + ENDDO + CALL scan_end (iunps, "DIJ") - if (isus (is) ) then - call scan_begin (iunps, "QIJ", .false.) - read (iunps, *, err = 100, iostat = ios) nqf(is) + IF (isus (is) ) THEN + CALL scan_begin (iunps, "QIJ", .false.) + READ (iunps, *, err = 100, iostat = ios) nqf(is) nqlc (is)= 2 * lmax (is) + 1 - if (nqlc(is).gt.lqmax .or. nqlc(is).lt.0) & - call errore (' read_pseudo_nl', 'Wrong nqlc', nqlc (is) ) - if (nqf(is).ne.0) then - call scan_begin (iunps, "RINNER", .false.) - read (iunps,*,err=100,iostat=ios) & + IF (nqlc(is)>lqmax .or. nqlc(is)<0) & + CALL errore (' read_pseudo_nl', 'Wrong nqlc', nqlc (is) ) + IF (nqf(is)/=0) THEN + CALL scan_begin (iunps, "RINNER", .false.) + READ (iunps,*,err=100,iostat=ios) & (idum,rinner(i,is),i=1,nqlc(is)) - call scan_end (iunps, "RINNER") - end if - do nb = 1, nbeta(is) - do mb = nb, nbeta(is) + CALL scan_end (iunps, "RINNER") + ENDIF + DO nb = 1, nbeta(is) + DO mb = nb, nbeta(is) - read (iunps,*,err=100,iostat=ios) idum, idum, ldum, dummy + READ (iunps,*,err=100,iostat=ios) idum, idum, ldum, dummy !" i j (l)" - if (ldum.ne.lll(mb,is) ) call errore ('read_pseudo_nl', & + IF (ldum/=lll(mb,is) ) CALL errore ('read_pseudo_nl', & 'inconsistent angular momentum for Q_ij', 1) - read (iunps,*,err=100,iostat=ios) qqq(nb,mb,is), dummy + READ (iunps,*,err=100,iostat=ios) qqq(nb,mb,is), dummy ! "Q_int" - qqq(mb,nb,is) = qqq(nb,mb,is) + qqq(mb,nb,is) = qqq(nb,mb,is) - read (iunps,*,err=100,iostat=ios) & + READ (iunps,*,err=100,iostat=ios) & (qfunc(n,nb,mb,is), n=1,mesh(is)) - do n = 0, mesh (is) - qfunc(n,mb,nb,is) = qfunc(n,nb,mb,is) - enddo + DO n = 0, mesh (is) + qfunc(n,mb,nb,is) = qfunc(n,nb,mb,is) + ENDDO - if (nqf(is).gt.0) then - call scan_begin (iunps, "QFCOEF", .false.) - read (iunps,*,err=100,iostat=ios) & + IF (nqf(is)>0) THEN + CALL scan_begin (iunps, "QFCOEF", .false.) + READ (iunps,*,err=100,iostat=ios) & ((qfcoef(i,lp,nb,mb,is),i=1,nqf(is)),lp=1,nqlc(is)) - call scan_end (iunps, "QFCOEF") - end if + CALL scan_end (iunps, "QFCOEF") + ENDIF - enddo - enddo - call scan_end (iunps, "QIJ") - else + ENDDO + ENDDO + CALL scan_end (iunps, "QIJ") + ELSE qqq (:,:,is) = 0.d0 qfunc(:,:,:,is) =0.d0 - endif + ENDIF -100 call errore ('read_pseudo_nl', 'Reading pseudo file', abs (ios) ) - return -end subroutine read_pseudo_nl +100 CALL errore ('read_pseudo_nl', 'Reading pseudo file', abs (ios) ) + RETURN +END SUBROUTINE read_pseudo_nl ! !--------------------------------------------------------------------- -subroutine read_pseudo_nlcc (is, iunps) +SUBROUTINE read_pseudo_nlcc (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios - read (iunps, *, err = 100, iostat = ios) (rho_atc(ir,is), ir=1,mesh(is) ) + READ (iunps, *, err = 100, iostat = ios) (rho_atc(ir,is), ir=1,mesh(is) ) ! -100 call errore ('read_pseudo_nlcc', 'Reading pseudo file', abs (ios) ) - return -end subroutine read_pseudo_nlcc +100 CALL errore ('read_pseudo_nlcc', 'Reading pseudo file', abs (ios) ) + RETURN +END SUBROUTINE read_pseudo_nlcc ! !--------------------------------------------------------------------- -subroutine read_pseudo_pswfc (is, iunps) +SUBROUTINE read_pseudo_pswfc (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - character (len=75) :: dummy - integer :: nb, ir, ios + CHARACTER (len=75) :: dummy + INTEGER :: nb, ir, ios ! - do nb = 1, ntwfc(is) - read (iunps,*,err=100,iostat=ios) dummy !Wavefunction labels - read (iunps,*,err=100,iostat=ios) (chi(ir,nb,is), ir=1,mesh(is)) - enddo -100 call errore ('read_pseudo_pswfc', 'Reading pseudo file', abs(ios)) - return + DO nb = 1, ntwfc(is) + READ (iunps,*,err=100,iostat=ios) dummy !Wavefunction labels + READ (iunps,*,err=100,iostat=ios) (chi(ir,nb,is), ir=1,mesh(is)) + ENDDO +100 CALL errore ('read_pseudo_pswfc', 'Reading pseudo file', abs(ios)) + RETURN -end subroutine read_pseudo_pswfc +END SUBROUTINE read_pseudo_pswfc ! !--------------------------------------------------------------------- -subroutine read_pseudo_rhoatom (is, iunps) +SUBROUTINE read_pseudo_rhoatom (is, iunps) !--------------------------------------------------------------------- ! - use pseudo - implicit none + USE pseudo + IMPLICIT NONE ! - integer :: is, iunps + INTEGER :: is, iunps ! - integer :: ir, ios + INTEGER :: ir, ios - read (iunps,*,err=100,iostat=ios) (rho_at(ir,is), ir=1,mesh(is)) - return + READ (iunps,*,err=100,iostat=ios) (rho_at(ir,is), ir=1,mesh(is)) + RETURN -100 call errore ('read_pseudo_rhoatom','Reading pseudo file',abs(ios)) +100 CALL errore ('read_pseudo_rhoatom','Reading pseudo file',abs(ios)) -end subroutine read_pseudo_rhoatom +END SUBROUTINE read_pseudo_rhoatom diff --git a/upftools/write_upf.f90 b/upftools/write_upf.f90 index 86201b7fc..610dbb87c 100644 --- a/upftools/write_upf.f90 +++ b/upftools/write_upf.f90 @@ -5,445 +5,445 @@ ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! -module upf +MODULE upf ! ! All variables to be written into the UPF file ! (UPF = unified pseudopotential format) ! ! pp_info - integer :: rel + INTEGER :: rel real(8) :: rcloc - integer :: nwfs - real(8), allocatable :: oc(:), rcut(:), rcutus(:), epseu(:) - character(len=2), allocatable :: els(:) - integer, allocatable:: lchi (:), nns (:) + INTEGER :: nwfs + real(8), ALLOCATABLE :: oc(:), rcut(:), rcutus(:), epseu(:) + CHARACTER(len=2), ALLOCATABLE :: els(:) + INTEGER, ALLOCATABLE:: lchi (:), nns (:) ! ! pp_header - character (len=80):: generated, date_author, comment - character (len=2) :: psd, pseudotype - integer :: nv = 0 - integer :: iexch, icorr, igcx, igcc - integer :: lmax, mesh, nbeta, ntwfc - logical :: nlcc + CHARACTER (len=80):: generated, date_author, comment + CHARACTER (len=2) :: psd, pseudotype + INTEGER :: nv = 0 + INTEGER :: iexch, icorr, igcx, igcc + INTEGER :: lmax, mesh, nbeta, ntwfc + LOGICAL :: nlcc real(8) :: zp, ecutrho, ecutwfc, etotps - real(8), allocatable :: ocw(:) - character(len=2), allocatable :: elsw(:) - integer, allocatable:: lchiw(:) + real(8), ALLOCATABLE :: ocw(:) + CHARACTER(len=2), ALLOCATABLE :: elsw(:) + INTEGER, ALLOCATABLE:: lchiw(:) ! ! pp_mesh - real(8), allocatable :: r(:), rab(:) + real(8), ALLOCATABLE :: r(:), rab(:) ! ! pp_nlcc - real(8), allocatable :: rho_atc(:) + real(8), ALLOCATABLE :: rho_atc(:) ! ! pp_local - real(8), allocatable :: vloc0(:) + real(8), ALLOCATABLE :: vloc0(:) ! ! pp_nonlocal ! pp_beta - real(8), allocatable :: betar(:,:) - integer, allocatable:: lll(:), ikk2(:) + real(8), ALLOCATABLE :: betar(:,:) + INTEGER, ALLOCATABLE:: lll(:), ikk2(:) ! pp_dij - real(8), allocatable :: dion(:,:) + real(8), ALLOCATABLE :: dion(:,:) ! pp_qij - integer :: nqf, nqlc - real(8), allocatable :: rinner(:), qqq(:,:), qfunc(:,:,:) + INTEGER :: nqf, nqlc + real(8), ALLOCATABLE :: rinner(:), qqq(:,:), qfunc(:,:,:) ! pp_qfcoef - real(8), allocatable :: qfcoef(:,:,:,:) + real(8), ALLOCATABLE :: qfcoef(:,:,:,:) ! ! pp_pswfc - real(8), allocatable :: chi(:,:) + real(8), ALLOCATABLE :: chi(:,:) ! ! pp_rhoatom - real(8), allocatable :: rho_at(:) -end module upf + real(8), ALLOCATABLE :: rho_at(:) +END MODULE upf ! -subroutine write_upf(ounps) +SUBROUTINE write_upf(ounps) - use upf, only: nlcc + USE upf, ONLY: nlcc - integer :: ounps + INTEGER :: ounps - call write_pseudo_comment(ounps) - call write_pseudo_header(ounps) - call write_pseudo_mesh(ounps) - if (nlcc) call write_pseudo_nlcc(ounps) - call write_pseudo_local(ounps) - call write_pseudo_nl(ounps) - call write_pseudo_pswfc(ounps) - call write_pseudo_rhoatom(ounps) + CALL write_pseudo_comment(ounps) + CALL write_pseudo_header(ounps) + CALL write_pseudo_mesh(ounps) + IF (nlcc) CALL write_pseudo_nlcc(ounps) + CALL write_pseudo_local(ounps) + CALL write_pseudo_nl(ounps) + CALL write_pseudo_pswfc(ounps) + CALL write_pseudo_rhoatom(ounps) ! - print '("*** PLEASE TEST BEFORE USING!!! ***")' - print '("review the content of the PP_INFO fields")' + PRINT '("*** PLEASE TEST BEFORE USING!!! ***")' + PRINT '("review the content of the PP_INFO fields")' ! -end subroutine write_upf +END SUBROUTINE write_upf ! !--------------------------------------------------------------------- - subroutine write_pseudo_comment (ounps) + SUBROUTINE write_pseudo_comment (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the comments of the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps - integer :: nb, ios + INTEGER :: nb, ios - write (ounps, '(a9)', err = 100, iostat = ios) "" + WRITE (ounps, '(a9)', err = 100, iostat = ios) "" - write (ounps, '(a)', err = 100, iostat = ios) generated - write (ounps, '(a)', err = 100, iostat = ios) date_author - write (ounps, '(a)', err = 100, iostat = ios) comment - if (rel==2) then - write (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,& + WRITE (ounps, '(a)', err = 100, iostat = ios) generated + WRITE (ounps, '(a)', err = 100, iostat = ios) date_author + WRITE (ounps, '(a)', err = 100, iostat = ios) comment + IF (rel==2) THEN + WRITE (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,& &"The Pseudo was generated with a Full-Relativistic Calculation" - else if (rel==1) then - write (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,& + ELSEIF (rel==1) THEN + WRITE (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,& &"The Pseudo was generated with a Scalar-Relativistic Calculation" - else if (rel==0) then - write (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel, & + ELSEIF (rel==0) THEN + WRITE (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel, & & "The Pseudo was generated with a Non-Relativistic Calculation" - endif + ENDIF - if (rcloc > 0.d0) & - write (ounps, '(1pe19.11,t24,a)', err = 100, iostat = ios) & + IF (rcloc > 0.d0) & + WRITE (ounps, '(1pe19.11,t24,a)', err = 100, iostat = ios) & rcloc, "Local Potential cutoff radius" - if (nwfs>0) & - write (ounps, '(a2,2a3,a6,3a19)', err = 100, iostat = ios) "nl", & + IF (nwfs>0) & + WRITE (ounps, '(a2,2a3,a6,3a19)', err = 100, iostat = ios) "nl", & &" pn", "l", "occ", "Rcut", "Rcut US", "E pseu" - do nb = 1, nwfs - write (ounps, '(a2,2i3,f6.2,3f19.11)') els (nb) , nns (nb) , & + DO nb = 1, nwfs + WRITE (ounps, '(a2,2i3,f6.2,3f19.11)') els (nb) , nns (nb) , & lchi (nb) , oc (nb) , rcut (nb) , rcutus (nb) , epseu(nb) - enddo + ENDDO - write (ounps, '(a10)', err = 100, iostat = ios) "" - return -100 write(6,'("write_pseudo_comment: error writing pseudopotential file")') - stop + WRITE (ounps, '(a10)', err = 100, iostat = ios) "" + RETURN +100 WRITE(6,'("write_pseudo_comment: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_comment + END SUBROUTINE write_pseudo_comment ! !--------------------------------------------------------------------- - subroutine write_pseudo_header (ounps) + SUBROUTINE write_pseudo_header (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the header of the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps ! - character (len=4) :: shortname - character (len=20):: dft - integer :: nb, ios + CHARACTER (len=4) :: shortname + CHARACTER (len=20):: dft + INTEGER :: nb, ios ! ! - write (ounps, '(//a11)', err = 100, iostat = ios) "" + WRITE (ounps, '(//a11)', err = 100, iostat = ios) "" - write (ounps, '(t3,i2,t24,a)', err = 100, iostat = ios) nv, & + WRITE (ounps, '(t3,i2,t24,a)', err = 100, iostat = ios) nv, & "Version Number" - write (ounps, '(t3,a,t24,a)', err = 100, iostat = ios) psd , & + WRITE (ounps, '(t3,a,t24,a)', err = 100, iostat = ios) psd , & "Element" - if (pseudotype == 'NC') then - write (ounps, '(a5,t24,a)', err = 100, iostat = ios) "NC", & + IF (pseudotype == 'NC') THEN + WRITE (ounps, '(a5,t24,a)', err = 100, iostat = ios) "NC", & "Norm - Conserving pseudopotential" - else if (pseudotype == 'US') then - write (ounps, '(a5,t24,a)', err = 100, iostat = ios) "US", & + ELSEIF (pseudotype == 'US') THEN + WRITE (ounps, '(a5,t24,a)', err = 100, iostat = ios) "US", & "Ultrasoft pseudopotential" - else - write(6,'("write_pseudo_header: unknown PP type ",A)') pseudotype - stop - endif - write (ounps, '(l5,t24,a)', err = 100, iostat = ios) nlcc , & + ELSE + WRITE(6,'("write_pseudo_header: unknown PP type ",A)') pseudotype + STOP + ENDIF + WRITE (ounps, '(l5,t24,a)', err = 100, iostat = ios) nlcc , & "Nonlinear Core Correction" - call dftname (iexch, icorr, igcx, igcc, dft, shortname) - write (ounps, '(a,t24,a4,a)', err = 100, iostat = ios) & + CALL dftname (iexch, icorr, igcx, igcc, dft, shortname) + WRITE (ounps, '(a,t24,a4,a)', err = 100, iostat = ios) & dft, shortname," Exchange-Correlation functional" - write (ounps, '(f17.11,t24,a)') zp , "Z valence" - write (ounps, '(f17.11,t24,a)') etotps, "Total energy" - write (ounps, '(2f11.7,t24,a)') ecutrho, ecutwfc, & - "Suggested cutoff for wfc and rho" + WRITE (ounps, '(f17.11,t24,a)') zp , "Z valence" + WRITE (ounps, '(f17.11,t24,a)') etotps, "Total energy" + WRITE (ounps, '(2f11.7,t24,a)') ecutrho, ecutwfc, & + "Suggested cutoff for wfc and rho" - write (ounps, '(i5,t24,a)') lmax, "Max angular momentum component" - write (ounps, '(i5,t24,a)') mesh, "Number of points in mesh" - write (ounps, '(2i5,t24,a)', err = 100, iostat = ios) ntwfc, & + WRITE (ounps, '(i5,t24,a)') lmax, "Max angular momentum component" + WRITE (ounps, '(i5,t24,a)') mesh, "Number of points in mesh" + WRITE (ounps, '(2i5,t24,a)', err = 100, iostat = ios) ntwfc, & nbeta , "Number of Wavefunctions, Number of Projectors" - write (ounps, '(a,t24,a2,a3,a6)', err = 100, iostat = ios) & + WRITE (ounps, '(a,t24,a2,a3,a6)', err = 100, iostat = ios) & " Wavefunctions", "nl", "l", "occ" - do nb = 1, ntwfc - write (ounps, '(t24,a2,i3,f6.2)') elsw(nb), lchiw(nb), ocw(nb) - enddo + DO nb = 1, ntwfc + WRITE (ounps, '(t24,a2,i3,f6.2)') elsw(nb), lchiw(nb), ocw(nb) + ENDDO !---> End header writing - write (ounps, '(a12)', err = 100, iostat = ios) "" - return -100 write(6,'("write_pseudo_header: error writing pseudopotential file")') - stop + WRITE (ounps, '(a12)', err = 100, iostat = ios) "" + RETURN +100 WRITE(6,'("write_pseudo_header: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_header + END SUBROUTINE write_pseudo_header ! !--------------------------------------------------------------------- - subroutine write_pseudo_mesh (ounps) + SUBROUTINE write_pseudo_mesh (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the atomic charge density to the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps ! - integer :: ir, ios + INTEGER :: ir, ios ! - write (ounps, '(//a9)', err = 100, iostat = ios) "" + WRITE (ounps, '(//a9)', err = 100, iostat = ios) "" - write (ounps, '(t3,a6)', err = 100, iostat = ios) "" - write (ounps, '(1p4e19.11)', err=100, iostat=ios) (r(ir), ir=1,mesh ) - write (ounps, '(t3,a7)', err = 100, iostat = ios) "" - write (ounps, '(t3,a8)', err = 100, iostat = ios) "" - write (ounps, '(1p4e19.11)', err=100, iostat=ios) (rab(ir), ir=1,mesh ) - write (ounps, '(t3,a9)', err = 100, iostat = ios) "" + WRITE (ounps, '(t3,a6)', err = 100, iostat = ios) "" + WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) (r(ir), ir=1,mesh ) + WRITE (ounps, '(t3,a7)', err = 100, iostat = ios) "" + WRITE (ounps, '(t3,a8)', err = 100, iostat = ios) "" + WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) (rab(ir), ir=1,mesh ) + WRITE (ounps, '(t3,a9)', err = 100, iostat = ios) "" - write (ounps, '(a10)', err = 100, iostat = ios) "" + WRITE (ounps, '(a10)', err = 100, iostat = ios) "" - return + RETURN -100 write(6,'("write_pseudo_mesh: error writing pseudopotential file")') - stop +100 WRITE(6,'("write_pseudo_mesh: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_mesh + END SUBROUTINE write_pseudo_mesh ! !--------------------------------------------------------------------- - subroutine write_pseudo_nlcc (ounps) + SUBROUTINE write_pseudo_nlcc (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the core charge for the nonlinear core ! correction of the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps ! - integer :: ir, ios + INTEGER :: ir, ios - write (ounps, '(//a9)', err = 100, iostat = ios) "" + WRITE (ounps, '(//a9)', err = 100, iostat = ios) "" - write (ounps, '(1p4e19.11)', err=100, iostat=ios) & + WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) & ( rho_atc(ir), ir = 1, mesh ) - write (ounps, '(a10)', err = 100, iostat = ios) "" - return + WRITE (ounps, '(a10)', err = 100, iostat = ios) "" + RETURN -100 write(6,'("write_pseudo_nlcc: error writing pseudopotential file")') - stop +100 WRITE(6,'("write_pseudo_nlcc: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_nlcc + END SUBROUTINE write_pseudo_nlcc ! !--------------------------------------------------------------------- - subroutine write_pseudo_local (ounps) + SUBROUTINE write_pseudo_local (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the local part of the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps ! - integer :: ir, ios + INTEGER :: ir, ios - write (ounps, '(//a10)', err = 100, iostat = ios) "" - write (ounps, '(1p4e19.11)', err=100, iostat=ios) & + WRITE (ounps, '(//a10)', err = 100, iostat = ios) "" + WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) & ( vloc0(ir), ir = 1, mesh ) - write (ounps, '(a11)', err = 100, iostat = ios) "" - return -100 write(6,'("write_pseudo_local: error writing pseudopotential file")') - stop + WRITE (ounps, '(a11)', err = 100, iostat = ios) "" + RETURN +100 WRITE(6,'("write_pseudo_local: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_local + END SUBROUTINE write_pseudo_local ! !--------------------------------------------------------------------- - subroutine write_pseudo_nl (ounps) + SUBROUTINE write_pseudo_nl (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the non local part of the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps ! - integer :: nb, mb, n, ir, nd, i, lp, ios + INTEGER :: nb, mb, n, ir, nd, i, lp, ios - write (ounps, '(//a13)', err = 100, iostat = ios) "" - do nb = 1, nbeta - write (ounps, '(t3,a9)', err = 100, iostat = ios) "" - write (ounps, '(2i5,t24,a)', err=100, iostat=ios) & + WRITE (ounps, '(//a13)', err = 100, iostat = ios) "" + DO nb = 1, nbeta + WRITE (ounps, '(t3,a9)', err = 100, iostat = ios) "" + WRITE (ounps, '(2i5,t24,a)', err=100, iostat=ios) & nb, lll(nb), "Beta L" - write (ounps, '(i6)', err=100, iostat=ios) ikk2 (nb) - write (ounps, '(1p4e19.11)', err=100, iostat=ios) & + WRITE (ounps, '(i6)', err=100, iostat=ios) ikk2 (nb) + WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) & ( betar(ir,nb), ir=1,ikk2(nb) ) - write (ounps, '(t3,a10)', err = 100, iostat = ios) "" - enddo + WRITE (ounps, '(t3,a10)', err = 100, iostat = ios) "" + ENDDO - write (ounps, '(t3,a8)', err = 100, iostat = ios) "" - nd = 0 - do nb = 1, nbeta - do mb = nb, nbeta - if ( abs(dion(nb,mb)) .gt. 1.0d-12 ) nd = nd + 1 - enddo - enddo - write (ounps, '(1p,i5,t24,a)', err=100, iostat=ios) & + WRITE (ounps, '(t3,a8)', err = 100, iostat = ios) "" + nd = 0 + DO nb = 1, nbeta + DO mb = nb, nbeta + IF ( abs(dion(nb,mb)) > 1.0d-12 ) nd = nd + 1 + ENDDO + ENDDO + WRITE (ounps, '(1p,i5,t24,a)', err=100, iostat=ios) & nd, "Number of nonzero Dij" - do nb = 1, nbeta - do mb = nb, nbeta - if ( abs(dion(nb,mb)) .gt. 1.0d-12 ) & - write(ounps,'(1p,2i5,e19.11)', err=100, iostat=ios) & + DO nb = 1, nbeta + DO mb = nb, nbeta + IF ( abs(dion(nb,mb)) > 1.0d-12 ) & + WRITE(ounps,'(1p,2i5,e19.11)', err=100, iostat=ios) & nb, mb, dion(nb,mb) - enddo - enddo - write (ounps, '(t3,a9)', err=100, iostat=ios) "" + ENDDO + ENDDO + WRITE (ounps, '(t3,a9)', err=100, iostat=ios) "" - if (pseudotype == 'US') then - write (ounps, '(t3,a8)', err = 100, iostat = ios) "" - write (ounps, '(i5,a)',err=100, iostat=ios) nqf," nqf.& + IF (pseudotype == 'US') THEN + WRITE (ounps, '(t3,a8)', err = 100, iostat = ios) "" + WRITE (ounps, '(i5,a)',err=100, iostat=ios) nqf," nqf.& & If not zero, Qij's inside rinner are computed using qfcoef's" - if (nqf.gt.0) then - write (ounps, '(t5,a11)', err=100, iostat=ios) "" - write (ounps,'(i5,1pe19.11)', err=100, iostat=ios) & + IF (nqf>0) THEN + WRITE (ounps, '(t5,a11)', err=100, iostat=ios) "" + WRITE (ounps,'(i5,1pe19.11)', err=100, iostat=ios) & (i, rinner(i), i = 1, nqlc) - write (ounps, '(t5,a12)', err=100, iostat=ios) "" - end if - do nb = 1, nbeta - do mb = nb, nbeta - write (ounps, '(3i5,t24,a)', err=100, iostat=ios) & + WRITE (ounps, '(t5,a12)', err=100, iostat=ios) "" + ENDIF + DO nb = 1, nbeta + DO mb = nb, nbeta + WRITE (ounps, '(3i5,t24,a)', err=100, iostat=ios) & nb, mb, lll(mb) , "i j (l(j))" - write (ounps, '(1pe19.11,t24,a)', err=100, iostat=ios) & + WRITE (ounps, '(1pe19.11,t24,a)', err=100, iostat=ios) & qqq(nb,mb), "Q_int" - write (ounps, '(1p4e19.11)', err=100, iostat=ios) & + WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) & ( qfunc (n,nb,mb), n=1,mesh ) - if (nqf.gt.0) then - write (ounps, '(t5,a11)', err=100, iostat=ios) & - "" - write(ounps,'(1p4e19.11)', err=100, iostat=ios) & + IF (nqf>0) THEN + WRITE (ounps, '(t5,a11)', err=100, iostat=ios) & + "" + WRITE(ounps,'(1p4e19.11)', err=100, iostat=ios) & ((qfcoef(i,lp,nb,mb),i=1,nqf),lp=1,nqlc) - write (ounps, '(t5,a12)', err=100, iostat=ios) & + WRITE (ounps, '(t5,a12)', err=100, iostat=ios) & "" - end if - enddo - enddo - write (ounps, '(t3,a9)', err = 100, iostat = ios) "" + ENDIF + ENDDO + ENDDO + WRITE (ounps, '(t3,a9)', err = 100, iostat = ios) "" - endif - write (ounps, '(a14)', err = 100, iostat = ios) "" - return + ENDIF + WRITE (ounps, '(a14)', err = 100, iostat = ios) "" + RETURN -100 write(6,'("write_pseudo_nl: error writing pseudopotential file")') - stop +100 WRITE(6,'("write_pseudo_nl: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_nl + END SUBROUTINE write_pseudo_nl ! !--------------------------------------------------------------------- - subroutine write_pseudo_pswfc (ounps) + SUBROUTINE write_pseudo_pswfc (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the pseudo atomic functions ! of the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps ! - integer :: nb, ir, ios + INTEGER :: nb, ir, ios - write (ounps, '(//a10)', err = 100, iostat = ios) "" - do nb = 1, ntwfc - write (ounps,'(a2,i5,f6.2,t24,a)', err=100, iostat=ios) & + WRITE (ounps, '(//a10)', err = 100, iostat = ios) "" + DO nb = 1, ntwfc + WRITE (ounps,'(a2,i5,f6.2,t24,a)', err=100, iostat=ios) & elsw(nb), lchiw(nb), ocw(nb), "Wavefunction" - write (ounps, '(1p4e19.11)', err=100, iostat=ios) & + WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) & ( chi(ir,nb), ir=1,mesh ) - enddo - write (ounps, '(a11)', err = 100, iostat = ios) "" - return + ENDDO + WRITE (ounps, '(a11)', err = 100, iostat = ios) "" + RETURN -100 write(6,'("write_pseudo_pswfc: error writing pseudopotential file")') - stop +100 WRITE(6,'("write_pseudo_pswfc: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_pswfc + END SUBROUTINE write_pseudo_pswfc ! !--------------------------------------------------------------------- - subroutine write_pseudo_rhoatom (ounps) + SUBROUTINE write_pseudo_rhoatom (ounps) !--------------------------------------------------------------------- ! ! ! This routine writes the atomic charge density to the new UPF file ! - use upf - implicit none - integer :: ounps + USE upf + IMPLICIT NONE + INTEGER :: ounps ! - integer :: ir, ios + INTEGER :: ir, ios - write (ounps, '(//a12)', err = 100, iostat = ios) "" - write (ounps, '(1p4e19.11)', err = 100, iostat = ios) & + WRITE (ounps, '(//a12)', err = 100, iostat = ios) "" + WRITE (ounps, '(1p4e19.11)', err = 100, iostat = ios) & ( rho_at(ir), ir=1,mesh ) - write (ounps, '(a13)', err = 100, iostat = ios) "" - return + WRITE (ounps, '(a13)', err = 100, iostat = ios) "" + RETURN -100 write(6,'("write_pseudo_rhoatom: error writing pseudopotential file")') - stop +100 WRITE(6,'("write_pseudo_rhoatom: error writing pseudopotential file")') + STOP - end subroutine write_pseudo_rhoatom + END SUBROUTINE write_pseudo_rhoatom !--------------------------------------------------------------------- - subroutine dftname(iexch, icorr, igcx, igcc, longname, shortname) + SUBROUTINE dftname(iexch, icorr, igcx, igcc, longname, shortname) !--------------------------------------------------------------------- - implicit none - integer iexch, icorr, igcx, igcc - character (len=4) :: shortname - character (len=20):: longname + IMPLICIT NONE + INTEGER iexch, icorr, igcx, igcc + CHARACTER (len=4) :: shortname + CHARACTER (len=20):: longname ! ! The data used to convert iexch, icorr, igcx, igcc ! into a user-readable string ! - integer, parameter :: nxc = 6, ncc = 9, ngcx = 4, ngcc = 5 - character (len=20) :: exc, corr, gradx, gradc - dimension exc (0:nxc), corr (0:ncc), gradx (0:ngcx), gradc (0:ngcc) - data exc / 'NOX ', 'SLA ', 'SL1 ', 'RXC ', 'OEP ', 'HF ', 'PB0X' / + INTEGER, PARAMETER :: nxc = 6, ncc = 9, ngcx = 4, ngcc = 5 + CHARACTER (len=20) :: exc, corr, gradx, gradc + DIMENSION exc (0:nxc), corr (0:ncc), gradx (0:ngcx), gradc (0:ngcc) + data exc / 'NOX ', 'SLA ', 'SL1 ', 'RXC ', 'OEP ', 'HF ', 'PB0X' / data corr / 'NOC ', 'PZ ', 'VWN ', 'LYP ', 'PW ', 'WIG ', 'HL ',& 'OBZ ', 'OBW ', 'GL ' / - data gradx / 'NOGX', 'B88 ', 'GGX ', 'PBE ', 'TPSS' / - data gradc / 'NOGC', 'P86 ', 'GGC ', 'BLYP', 'PBE ', 'TPSS' / + data gradx / 'NOGX', 'B88 ', 'GGX ', 'PBE ', 'TPSS' / + data gradc / 'NOGC', 'P86 ', 'GGC ', 'BLYP', 'PBE ', 'TPSS' / - if (iexch==1.and.igcx==0.and.igcc==0) then + IF (iexch==1.and.igcx==0.and.igcc==0) THEN shortname = corr(icorr) - else if (iexch==1.and.icorr==3.and.igcx==1.and.igcc==3) then + ELSEIF (iexch==1.and.icorr==3.and.igcx==1.and.igcc==3) THEN shortname = 'BLYP' - else if (iexch==1.and.icorr==1.and.igcx==1.and.igcc==0) then + ELSEIF (iexch==1.and.icorr==1.and.igcx==1.and.igcc==0) THEN shortname = 'B88' - else if (iexch==1.and.icorr==1.and.igcx==1.and.igcc==1) then + ELSEIF (iexch==1.and.icorr==1.and.igcx==1.and.igcc==1) THEN shortname = 'BP' - else if (iexch==1.and.icorr==4.and.igcx==2.and.igcc==2) then + ELSEIF (iexch==1.and.icorr==4.and.igcx==2.and.igcc==2) THEN shortname = 'PW91' - else if (iexch==1.and.icorr==4.and.igcx==3.and.igcc==4) then + ELSEIF (iexch==1.and.icorr==4.and.igcx==3.and.igcc==4) THEN shortname = 'PBE' - else if (iexch==1.and.icorr==4.and.igcx==4.and.igcc==5) then + ELSEIF (iexch==1.and.icorr==4.and.igcx==4.and.igcc==5) THEN shortname = 'TPSS' - else + ELSE shortname = ' ' - end if - write(longname,'(4a5)') exc(iexch),corr(icorr),gradx(igcx),gradc(igcc) + ENDIF + WRITE(longname,'(4a5)') exc(iexch),corr(icorr),gradx(igcx),gradc(igcc) - return -end subroutine dftname + RETURN +END SUBROUTINE dftname