Commit 1cde15db authored by Daniel Wortmann's avatar Daniel Wortmann

Bugfixes

parent ff0cec65
...@@ -47,8 +47,8 @@ ...@@ -47,8 +47,8 @@
TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_atoms),INTENT(IN) :: atoms
REAL, INTENT(IN) :: efermiarg, bandgap REAL, INTENT(IN) :: efermiarg, bandgap
LOGICAL, INTENT(IN) :: l_mcd LOGICAL, INTENT(IN) :: l_mcd
! locals ! locals
INTEGER, PARAMETER :: lmax= 4, ned = 1301 INTEGER, PARAMETER :: lmax= 4, ned = 1301
INTEGER i,s,v,index,jspin,k,l,l1,l2,ln,n,nl,ntb,ntria,ntetra INTEGER i,s,v,index,jspin,k,l,l1,l2,ln,n,nl,ntb,ntria,ntetra
...@@ -73,7 +73,7 @@ ...@@ -73,7 +73,7 @@
qdim = lmax*atoms%ntype+3 qdim = lmax*atoms%ntype+3
l_orbcomp = banddos%l_orb l_orbcomp = banddos%l_orb
IF (banddos%ndir.EQ.-3) THEN IF (banddos%ndir.EQ.-3) THEN
qdim = 2*slab%nsld qdim = 2*slab%nsld
n_orb = 0 n_orb = 0
IF (banddos%l_orb) THEN IF (banddos%l_orb) THEN
n_orb = banddos%orbCompAtom n_orb = banddos%orbCompAtom
...@@ -84,7 +84,7 @@ ...@@ -84,7 +84,7 @@
ALLOCATE( qal(qdim,dimension%neigd,kpts%nkpt),& ALLOCATE( qal(qdim,dimension%neigd,kpts%nkpt),&
& qval(vacuum%nstars*vacuum%layers*vacuum%nvac,dimension%neigd,kpts%nkpt),& & qval(vacuum%nstars*vacuum%layers*vacuum%nvac,dimension%neigd,kpts%nkpt),&
& qlay(dimension%neigd,vacuum%layerd,2)) & qlay(dimension%neigd,vacuum%layerd,2))
IF (l_mcd) THEN IF (l_mcd) THEN
ALLOCATE(mcd_local(3*atoms%ntype*ncored,dimension%neigd,kpts%nkpt) ) ALLOCATE(mcd_local(3*atoms%ntype*ncored,dimension%neigd,kpts%nkpt) )
ELSE ELSE
ALLOCATE(mcd_local(0,0,0)) ALLOCATE(mcd_local(0,0,0))
...@@ -95,11 +95,11 @@ ...@@ -95,11 +95,11 @@
emin =min(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const) emin =min(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const)
emax =max(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const) emax =max(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const)
efermi = efermiarg*hartree_to_ev_const efermi = efermiarg*hartree_to_ev_const
WRITE (6,'(a)') 'DOS-Output is generated!' WRITE (6,'(a)') 'DOS-Output is generated!'
IF ( NINT((emax - emin)/sigma) > ned ) THEN IF ( NINT((emax - emin)/sigma) > ned ) THEN
WRITE(6,*) 'sig_dos too small for DOS smoothing:' WRITE(6,*) 'sig_dos too small for DOS smoothing:'
WRITE(6,*) 'Reduce energy window or enlarge banddos%sig_dos!' WRITE(6,*) 'Reduce energy window or enlarge banddos%sig_dos!'
WRITE(6,*) 'For now: setting sigma to zero !' WRITE(6,*) 'For now: setting sigma to zero !'
sigma = 0.0 sigma = 0.0
...@@ -117,10 +117,10 @@ ...@@ -117,10 +117,10 @@
DO i=1,ned DO i=1,ned
e(i) = emin + (i-1)*de e(i) = emin + (i-1)*de
ENDDO ENDDO
IF ( l_mcd ) THEN ! create an energy grid for mcd-spectra IF ( l_mcd ) THEN ! create an energy grid for mcd-spectra
e_lo = 9.9*10.0**9 e_lo = 9.9*10.0**9
e_up = -9.9*10.0**9 e_up = -9.9*10.0**9
DO jspin = 1,input%jspins DO jspin = 1,input%jspins
DO n = 1,atoms%ntype DO n = 1,atoms%ntype
DO icore = 1 , mcd%ncore(n) DO icore = 1 , mcd%ncore(n)
...@@ -129,7 +129,7 @@ ...@@ -129,7 +129,7 @@
ENDDO ENDDO
ENDDO ENDDO
ENDDO ENDDO
e_lo = e_lo*hartree_to_ev_const - efermi - emax e_lo = e_lo*hartree_to_ev_const - efermi - emax
e_up = e_up*hartree_to_ev_const - efermi e_up = e_up*hartree_to_ev_const - efermi
de = (e_up-e_lo)/(ned-1) de = (e_up-e_lo)/(ned-1)
DO i=1,ned DO i=1,ned
...@@ -207,7 +207,7 @@ ...@@ -207,7 +207,7 @@
qal(lmax*atoms%ntype+1,i,k) = qal(lmax*atoms%ntype+1,i,k) - qmt qal(lmax*atoms%ntype+1,i,k) = qal(lmax*atoms%ntype+1,i,k) - qmt
ENDDO ENDDO
qal(lmax*atoms%ntype+1,i,k) = qal(lmax*atoms%ntype+1,i,k)& qal(lmax*atoms%ntype+1,i,k) = qal(lmax*atoms%ntype+1,i,k)&
-qal(lmax*atoms%ntype+2,i,k)*(3-vacuum%nvac) -qal(lmax*atoms%ntype+3,i,k)*(vacuum%nvac-1) -qal(lmax*atoms%ntype+2,i,k)*(3-vacuum%nvac) -qal(lmax*atoms%ntype+3,i,k)*(vacuum%nvac-1)
ENDDO ENDDO
ENDIF ENDIF
! !
...@@ -237,27 +237,17 @@ ...@@ -237,27 +237,17 @@
write(*,*) as,sym%nop2,l_tria write(*,*) as,sym%nop2,l_tria
! l_tria=.true. ! l_tria=.true.
ELSE ELSE
IF (input%l_inpXML) THEN IF (input%tria) THEN
IF (input%tria) THEN ntetra = kpts%ntet
ntetra = kpts%ntet DO i = 1, ntetra
DO i = 1, ntetra itetra(1:4,i) = kpts%ntetra(1:4,i)
itetra(1:4,i) = kpts%ntetra(1:4,i) voltet(i) = kpts%voltet(i) / ntetra
voltet(i) = kpts%voltet(i) / ntetra END DO
END DO l_tria = input%tria
l_tria = input%tria GOTO 67
GOTO 67 ELSE
ELSE GOTO 66
GOTO 66
END IF
END IF END IF
OPEN (41,file='kpts',FORM='formatted',STATUS='old')
DO i = 1, kpts%nkpt+1
READ (41,*,END=66,ERR=66)
ENDDO
READ (41,'(i5)',END=66,ERR=66) ntetra
READ (41,'(4(4i6,4x))') ((itetra(i,k),i=1,4),k=1,ntetra)
READ (41,'(4f20.13)') (voltet(k),k=1,ntetra)
CLOSE(41)
voltet(1:ntetra) = voltet(1:ntetra) / ntetra voltet(1:ntetra) = voltet(1:ntetra) / ntetra
l_tria=.true. l_tria=.true.
GOTO 67 GOTO 67
...@@ -265,7 +255,7 @@ ...@@ -265,7 +255,7 @@
CALL triang(kpts%bk,kpts%nkpt,itria,ntria,atr,as,l_tria) CALL triang(kpts%bk,kpts%nkpt,itria,ntria,atr,as,l_tria)
l_tria=.true. l_tria=.true.
! YM: tetrahedrons is not the way in 1D ! YM: tetrahedrons is not the way in 1D
IF (oneD%odi%d1) as = 0.0 IF (oneD%odi%d1) as = 0.0
IF (sym%invs) THEN IF (sym%invs) THEN
IF (abs(sym%nop2*as-1.0).GT.0.000001) l_tria=.false. IF (abs(sym%nop2*as-1.0).GT.0.000001) l_tria=.false.
ELSE ELSE
...@@ -328,9 +318,9 @@ ...@@ -328,9 +318,9 @@
CALL smooth(e,g(1,ln),sigma,ned) CALL smooth(e,g(1,ln),sigma,ned)
ENDDO ENDDO
ENDIF ENDIF
!*** sum up for all atoms !*** sum up for all atoms
IF (banddos%ndir.NE.-3) THEN IF (banddos%ndir.NE.-3) THEN
DO l = 1 , atoms%ntype DO l = 1 , atoms%ntype
l1 = lmax*(l-1) + 1 l1 = lmax*(l-1) + 1
...@@ -350,7 +340,7 @@ ...@@ -350,7 +340,7 @@
ENDDO ENDDO
ENDDO ENDDO
ENDIF ENDIF
!**** write out DOS !**** write out DOS
OPEN (18,FILE='DOS'//spin12(jspin)) OPEN (18,FILE='DOS'//spin12(jspin))
...@@ -409,14 +399,14 @@ ...@@ -409,14 +399,14 @@
ENDDO ENDDO
ENDDO ENDDO
ENDDO ENDDO
CLOSE (18) CLOSE (18)
ENDIF ENDIF
DEALLOCATE (g) DEALLOCATE (g)
! !
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
! now calculate the VACOS ! now calculate the VACOS
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
IF ( banddos%vacdos .and. input%film ) THEN IF ( banddos%vacdos .and. input%film ) THEN
ALLOCATE(g(ned,vacuum%nstars*vacuum%layers*vacuum%nvac)) ALLOCATE(g(ned,vacuum%nstars*vacuum%layers*vacuum%nvac))
! CALL ptdos( ! CALL ptdos(
...@@ -425,16 +415,16 @@ ...@@ -425,16 +415,16 @@
! < g) ! < g)
CALL ptdos(emin,emax,input%jspins,ned,vacuum%nstars*vacuum%nvac*vacuum%layers,ntb,ntria& CALL ptdos(emin,emax,input%jspins,ned,vacuum%nstars*vacuum%nvac*vacuum%layers,ntb,ntria&
,as,atr,2*kpts%nkpt,itria,kpts%nkpt,ev(1:ntb,1:kpts%nkpt), qval(:,1:ntb,1:kpts%nkpt),e,g) ,as,atr,2*kpts%nkpt,itria,kpts%nkpt,ev(1:ntb,1:kpts%nkpt), qval(:,1:ntb,1:kpts%nkpt),e,g)
!---- > smoothening !---- > smoothening
IF ( sigma.GT.0.0 ) THEN IF ( sigma.GT.0.0 ) THEN
DO ln = 1 , vacuum%nstars*vacuum%nvac*vacuum%layers DO ln = 1 , vacuum%nstars*vacuum%nvac*vacuum%layers
CALL smooth(e,g(1,ln),sigma,ned) CALL smooth(e,g(1,ln),sigma,ned)
ENDDO ENDDO
ENDIF ENDIF
! write VACDOS ! write VACDOS
OPEN (18,FILE='VACDOS'//spin12(jspin)) OPEN (18,FILE='VACDOS'//spin12(jspin))
! WRITE (18,'(i2,25(2x,i3))') Layers , (Zlay(l),l=1,Layers) ! WRITE (18,'(i2,25(2x,i3))') Layers , (Zlay(l),l=1,Layers)
DO i = 1 , ned DO i = 1 , ned
...@@ -464,11 +454,11 @@ ...@@ -464,11 +454,11 @@
END IF END IF
OPEN (18,FILE='bands'//spin12(jspin)) OPEN (18,FILE='bands'//spin12(jspin))
ntb = minval(results%neig(:,jspin)) ntb = minval(results%neig(:,jspin))
kx(1) = 0.0 kx(1) = 0.0
vkr(:,1)=matmul(kpts%bk(:,1),cell%bmat) vkr(:,1)=matmul(kpts%bk(:,1),cell%bmat)
DO k = 2, kpts%nkpt DO k = 2, kpts%nkpt
vkr(:,k)=matmul(kpts%bk(:,k),cell%bmat) vkr(:,k)=matmul(kpts%bk(:,k),cell%bmat)
dk = (vkr(1,k)-vkr(1,k-1))**2 + (vkr(2,k)-vkr(2,k-1) )**2 + & dk = (vkr(1,k)-vkr(1,k-1))**2 + (vkr(2,k)-vkr(2,k-1) )**2 + &
(vkr(3,k)-vkr(3,k-1))**2 (vkr(3,k)-vkr(3,k-1))**2
...@@ -483,7 +473,7 @@ ...@@ -483,7 +473,7 @@
ENDIF ENDIF
ENDDO ENDDO
! !
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
! for MCD calculations ... ! for MCD calculations ...
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
...@@ -512,4 +502,4 @@ ...@@ -512,4 +502,4 @@
99001 FORMAT (f10.5,110(1x,e10.3)) 99001 FORMAT (f10.5,110(1x,e10.3))
END SUBROUTINE evaldos END SUBROUTINE evaldos
END MODULE m_evaldos END MODULE m_evaldos
...@@ -80,11 +80,11 @@ CONTAINS ...@@ -80,11 +80,11 @@ CONTAINS
!*********************************************************************** !***********************************************************************
! ABBREVIATIONS ! ABBREVIATIONS
! !
! eig : array of eigenvalues ! eig : array of eigenvalues
! wtkpt : list of the weights of each k-point (from inp-file) ! wtkpt : list of the weights of each k-point (from inp-file)
! e : linear list of the eigenvalues ! e : linear list of the eigenvalues
! we : list of weights of the eigenvalues in e ! we : list of weights of the eigenvalues in e
! zelec : number of electrons ! zelec : number of electrons
! spindg : spindegeneracy (2 in nonmagnetic calculations) ! spindg : spindegeneracy (2 in nonmagnetic calculations)
! seigv : weighted sum of the occupied valence eigenvalues ! seigv : weighted sum of the occupied valence eigenvalues
! seigsc : weighted sum of the semi-core eigenvalues ! seigsc : weighted sum of the semi-core eigenvalues
...@@ -125,7 +125,7 @@ CONTAINS ...@@ -125,7 +125,7 @@ CONTAINS
IF (mpi%irank == 0) THEN IF (mpi%irank == 0) THEN
CALL read_eig(eig_id,k,jsp,neig=ne(k,jsp),eig=eig(:,k,jsp)) CALL read_eig(eig_id,k,jsp,neig=ne(k,jsp),eig=eig(:,k,jsp))
WRITE (6,'(a2,3f10.5,f12.6)') 'at',kpts%bk(:,k),kpts%wtkpt(k) WRITE (6,'(a2,3f10.5,f12.6)') 'at',kpts%bk(:,k),kpts%wtkpt(k)
WRITE (6,'(i5,a14)') ne(k,jsp),' eigenvalues :' WRITE (6,'(i5,a14)') ne(k,jsp),' eigenvalues :'
WRITE (6,'(8f12.6)') (eig(i,k,jsp),i=1,ne(k,jsp)) WRITE (6,'(8f12.6)') (eig(i,k,jsp),i=1,ne(k,jsp))
IF(.NOT.judft_was_argument("-minimalOutput")) THEN IF(.NOT.judft_was_argument("-minimalOutput")) THEN
attributes = '' attributes = ''
...@@ -144,7 +144,7 @@ CONTAINS ...@@ -144,7 +144,7 @@ CONTAINS
ENDDO ENDDO
!finished reading of eigenvalues !finished reading of eigenvalues
IF (mpi%irank == 0) CALL closeXMLElement('eigenvalues') IF (mpi%irank == 0) CALL closeXMLElement('eigenvalues')
IF (mpi%irank == 0) THEN IF (mpi%irank == 0) THEN
IF (ABS(input%fixed_moment)<1E-6) THEN IF (ABS(input%fixed_moment)<1E-6) THEN
!this is a standard calculation !this is a standard calculation
...@@ -167,7 +167,7 @@ CONTAINS ...@@ -167,7 +167,7 @@ CONTAINS
!Generate a list of energies !Generate a list of energies
DO k = 1,kpts%nkpt DO k = 1,kpts%nkpt
! !
!---> STORE EIGENVALUES AND WEIGHTS IN A LINEAR LIST. AND MEMORIZE !---> STORE EIGENVALUES AND WEIGHTS IN A LINEAR LIST. AND MEMORIZE
!---> CONECTION TO THE ORIGINAL ARRAYS !---> CONECTION TO THE ORIGINAL ARRAYS
! !
DO j = 1,ne(k,jsp) DO j = 1,ne(k,jsp)
...@@ -208,7 +208,7 @@ CONTAINS ...@@ -208,7 +208,7 @@ CONTAINS
IF ( mpi%irank == 0 ) THEN IF ( mpi%irank == 0 ) THEN
WRITE (6,FMT=8010) n,ws,weight WRITE (6,FMT=8010) n,ws,weight
END IF END IF
CALL juDFT_error("Not enough eavefunctions",calledby="fermie") CALL juDFT_error("Not enough weavefunctions",calledby="fermie")
8010 FORMAT (/,10x,'error: not enough wavefunctions.',i10,2d20.10) 8010 FORMAT (/,10x,'error: not enough wavefunctions.',i10,2d20.10)
END IF END IF
ws = ws + we(INDEX(l)) ws = ws + we(INDEX(l))
......
...@@ -86,38 +86,14 @@ c ...@@ -86,38 +86,14 @@ c
c---> write results of triang c---> write results of triang
IF (.not.film) THEN IF (.not.film) THEN
IF(input%l_inpXML) THEN ntetra = kpts%ntet
ntetra = kpts%ntet DO j = 1, ntetra
DO j = 1, ntetra itetra(1:4,j) = kpts%ntetra(1:4,j)
itetra(1:4,j) = kpts%ntetra(1:4,j) voltet(j) = kpts%voltet(j) / ntetra
voltet(j) = kpts%voltet(j) / ntetra END DO
END DO lb = MINVAL(eig(:,:,:)) - 0.01
ELSE ub = ef + 0.2
IF ( irank == 0 ) THEN CALL tetra_ef(
WRITE (6,*) 'reading tetrahedrons from file kpts'
END IF
OPEN (41,file='kpts',FORM='formatted',STATUS='old')
DO i = 1, nkpt+1
READ (41,*)
ENDDO
READ (41,'(i5)',ERR=66,END=66) ntetra
IF (ntetra>6*nkpt) CALL juDFT_error("ntetra > 6 nkpt"
+ ,calledby ="fertri")
READ (41,'(4(4i6,4x))') ((itetra(i,j),i=1,4),j=1,ntetra)
READ (41,'(4f20.13)') (voltet(j),j=1,ntetra)
voltet(1:ntetra) = voltet(1:ntetra) / ntetra
GOTO 67
66 CONTINUE ! no tetrahedron-information of file
CALL make_tetra(
> nkpt,bk,ntria,itria,atr,
< ntetra,itetra,voltet)!keep
67 CONTINUE ! tetrahedron-information read or created
CLOSE(41)
END IF
lb = MINVAL(eig(:,:,:)) - 0.01
ub = ef + 0.2
CALL tetra_ef(
> jspins,nkpt, > jspins,nkpt,
> lb,ub,eig,zc,sfac, > lb,ub,eig,zc,sfac,
> ntetra,itetra,voltet, > ntetra,itetra,voltet,
......
...@@ -278,9 +278,10 @@ MODULE m_types_atoms ...@@ -278,9 +278,10 @@ MODULE m_types_atoms
!force type !force type
xpath='' xpath=''
IF(xml%getNumberOfNodes(TRIM(ADJUSTL(xPaths))//'/force')==1) xpath=xpaths IF(xml%getNumberOfNodes(TRIM(ADJUSTL(xPaths))//'/force')==1) xpath=xpaths
IF(xml%getNumberOfNodes(TRIM(ADJUSTL(xPathg))//'/force')==1) xpath=xpathg
IF (xpath.NE.'') THEN IF (xpath.NE.'') THEN
this%l_geo(n) = evaluateFirstBoolOnly(xml%getAttributeValue(TRIM(ADJUSTL(xPathg))//'/force/@calculate')) this%l_geo(n) = evaluateFirstBoolOnly(xml%getAttributeValue(TRIM(ADJUSTL(xPath))//'/force/@calculate'))
valueString = xml%getAttributeValue(TRIM(ADJUSTL(xPathg))//'force/@relaxXYZ') valueString = xml%getAttributeValue(TRIM(ADJUSTL(xPath))//'/force/@relaxXYZ')
READ(valueString,'(3l1)') relaxX, relaxY, relaxZ READ(valueString,'(3l1)') relaxX, relaxY, relaxZ
IF (relaxX) this%relax(1,n) = 1 IF (relaxX) this%relax(1,n) = 1
IF (relaxY) this%relax(2,n) = 1 IF (relaxY) this%relax(2,n) = 1
......
...@@ -9,7 +9,7 @@ MODULE m_types_econfig ...@@ -9,7 +9,7 @@ MODULE m_types_econfig
PRIVATE PRIVATE
!This is used by t_atoms and does not extend t_fleurinput_base by itself !This is used by t_atoms and does not extend t_fleurinput_base by itself
TYPE:: t_econfig TYPE:: t_econfig
CHARACTER(len=100) :: coreconfig CHARACTER(len=200) :: coreconfig
CHARACTER(len=100) :: valenceconfig CHARACTER(len=100) :: valenceconfig
INTEGER :: num_core_states INTEGER :: num_core_states
INTEGER :: num_states INTEGER :: num_states
...@@ -29,7 +29,7 @@ MODULE m_types_econfig ...@@ -29,7 +29,7 @@ MODULE m_types_econfig
END TYPE t_econfig END TYPE t_econfig
PUBLIC :: t_econfig PUBLIC :: t_econfig
CONTAINS CONTAINS
SUBROUTINE get_core(econf,nst,nprnc,kappa,occupation,l_valence) SUBROUTINE get_core(econf,nst,nprnc,kappa,occupation,l_valence)
CLASS(t_econfig),INTENT(IN) :: econf CLASS(t_econfig),INTENT(IN) :: econf
INTEGER ,INTENT(out) :: nst INTEGER ,INTENT(out) :: nst
...@@ -46,9 +46,9 @@ CONTAINS ...@@ -46,9 +46,9 @@ CONTAINS
occupation(:nst,:)=econf%occupation(:nst,:) occupation(:nst,:)=econf%occupation(:nst,:)
if (size(occupation,2)==1) occupation=occupation*2 if (size(occupation,2)==1) occupation=occupation*2
END SUBROUTINE get_core END SUBROUTINE get_core
FUNCTION get_state_string(econf,i)RESULT(str) FUNCTION get_state_string(econf,i)RESULT(str)
CLASS(t_econfig),INTENT(IN):: econf CLASS(t_econfig),INTENT(IN):: econf
INTEGER,INTENT(in) :: i INTEGER,INTENT(in) :: i
...@@ -74,17 +74,17 @@ CONTAINS ...@@ -74,17 +74,17 @@ CONTAINS
CASE default CASE default
call judft_error("Invalid reqest for string with kappa") call judft_error("Invalid reqest for string with kappa")
END SELECT END SELECT
WRITE(str,"(a1,i1,a)") "(",econf%nprnc(i),s WRITE(str,"(a1,i1,a)") "(",econf%nprnc(i),s
END FUNCTION get_state_string END FUNCTION get_state_string
SUBROUTINE broadcast(econf,mpi_comm) SUBROUTINE broadcast(econf,mpi_comm)
USE m_mpi_bc_tool USE m_mpi_bc_tool
CLASS(t_econfig),INTENT(INOUT):: econf CLASS(t_econfig),INTENT(INOUT):: econf
INTEGER,INTENT(in) :: mpi_comm INTEGER,INTENT(in) :: mpi_comm
#ifdef CPP_MPI #ifdef CPP_MPI
INCLUDE 'mpif.h' INCLUDE 'mpif.h'
INTEGER :: ierr,irank INTEGER :: ierr,irank
...@@ -97,7 +97,7 @@ CONTAINS ...@@ -97,7 +97,7 @@ CONTAINS
CALL mpi_bc(econf%valence_electrons,0,mpi_comm) CALL mpi_bc(econf%valence_electrons,0,mpi_comm)
#endif #endif
END SUBROUTINE broadcast END SUBROUTINE broadcast
SUBROUTINE init_num(econf,nc,nz) SUBROUTINE init_num(econf,nc,nz)
USE m_constants USE m_constants
...@@ -114,8 +114,8 @@ CONTAINS ...@@ -114,8 +114,8 @@ CONTAINS
CALL econf%init(core,nz) CALL econf%init(core,nz)
END SUBROUTINE init_num END SUBROUTINE init_num
SUBROUTINE init_simple(econf,str) SUBROUTINE init_simple(econf,str)
CLASS(t_econfig),INTENT(OUT):: econf CLASS(t_econfig),INTENT(OUT):: econf
CHARACTER(len=*),INTENT(IN) :: str CHARACTER(len=*),INTENT(IN) :: str
...@@ -123,10 +123,10 @@ CONTAINS ...@@ -123,10 +123,10 @@ CONTAINS
CHARACTER(len=200)::core,valence CHARACTER(len=200)::core,valence
INTEGER :: n INTEGER :: n
REAL,allocatable :: core_occ(:),valence_occ(:) REAL,allocatable :: core_occ(:),valence_occ(:)
n=INDEX(str,"|") n=INDEX(str,"|")
IF (n==0) CALL judft_error(("Invalid econfig:"//TRIM(str))) IF (n==0) CALL judft_error(("Invalid econfig:"//TRIM(str)))
IF (INDEX(str,"|")==1) THEN IF (INDEX(str,"|")==1) THEN
! No core ! No core
core="" core=""
...@@ -157,15 +157,15 @@ CONTAINS ...@@ -157,15 +157,15 @@ CONTAINS
INTEGER :: np(40),kap(40) INTEGER :: np(40),kap(40)
econf%coreconfig=core econf%coreconfig=core
econf%valenceconfig=valence econf%valenceconfig=valence
CALL expand_noble_gas(core) !extend noble gas config CALL expand_noble_gas(core) !extend noble gas config
IF (VERIFY(core,"(1234567spdf/) ")>0) call judft_error(("Invalid econfig:"//TRIM(core))) IF (VERIFY(core,"(1234567spdf/) ")>0) call judft_error(("Invalid econfig:"//TRIM(core)))
IF (VERIFY(valence,"(1234567spdf/) ")>0) CALL judft_error(("Invalid econfig:"//TRIM(valence))) IF (VERIFY(valence,"(1234567spdf/) ")>0) CALL judft_error(("Invalid econfig:"//TRIM(valence)))
econf%num_core_states=0 econf%num_core_states=0
DO WHILE (len_TRIM(core)>1) DO WHILE (len_TRIM(core)>1)
CALL extract_next(core,np(econf%num_core_states+1),kap(econf%num_core_states+1)) CALL extract_next(core,np(econf%num_core_states+1),kap(econf%num_core_states+1))
...@@ -180,7 +180,7 @@ CONTAINS ...@@ -180,7 +180,7 @@ CONTAINS
econf%nprnc=np(:econf%num_states) econf%nprnc=np(:econf%num_states)
econf%kappa=kap(:econf%num_states) econf%kappa=kap(:econf%num_states)
ALLOCATE(econf%occupation(econf%num_states,2)) ALLOCATE(econf%occupation(econf%num_states,2))
CALL econf%set_occupation("(1s1/2)",-1.,-1.) CALL econf%set_occupation("(1s1/2)",-1.,-1.)
END SUBROUTINE init_all END SUBROUTINE init_all
...@@ -195,34 +195,34 @@ CONTAINS ...@@ -195,34 +195,34 @@ CONTAINS
CHARACTER(len=7)::str CHARACTER(len=7)::str
IF (nz>54) CALL judft_warn("Specifying no explicit valence config for systems with f-states might lead to broken configs") IF (nz>54) CALL judft_warn("Specifying no explicit valence config for systems with f-states might lead to broken configs")
econf%coreconfig=core econf%coreconfig=core
CALL expand_noble_gas(core) !extend noble gas config CALL expand_noble_gas(core) !extend noble gas config
IF (VERIFY(core,"(1234567spdf/) ")>0) call judft_error(("Invalid econfig:"//TRIM(core))) IF (VERIFY(core,"(1234567spdf/) ")>0) call judft_error(("Invalid econfig:"//TRIM(core)))
econf%num_core_states=0 econf%num_core_states=0
DO WHILE (len_TRIM(core)>1) DO WHILE (len_TRIM(core)>1)
CALL extract_next(core,np(econf%num_core_states+1),kap(econf%num_core_states+1)) CALL extract_next(core,np(econf%num_core_states+1),kap(econf%num_core_states+1))
econf%num_core_states=econf%num_core_states+1 econf%num_core_states=econf%num_core_states+1
ENDDO ENDDO
econf%num_states=econf%num_core_states econf%num_states=econf%num_core_states
!valence charge !valence charge
charge=nz-SUM(ABS(kap(:econf%num_core_states)))*2 charge=nz-SUM(ABS(kap(:econf%num_core_states))**2)*2
DO WHILE(charge>0) !Add valence DO WHILE(charge>0) !Add valence
str=coreStateList_const(econf%num_states+1) str=coreStateList_const(econf%num_states+1)
PRINT*,econf%num_states,str,charge PRINT*,econf%num_states,str,charge
CALL extract_next(str,np(econf%num_states+1),kap(econf%num_states+1)) CALL extract_next(str,np(econf%num_states+1),kap(econf%num_states+1))
econf%num_states=econf%num_states+1 econf%num_states=econf%num_states+1
charge=charge-ABS(kap(econf%num_states)) charge=charge-ABS(kap(econf%num_states)**2)*2
ENDDO ENDDO
ALLOCATE(econf%nprnc(econf