Commit a93a5b2a authored by Gregor Michalicek's avatar Gregor Michalicek

Some more simplifications in main/cdngen.F90

parent 36be5e32
......@@ -8,7 +8,7 @@ MODULE m_magMoms
CONTAINS
SUBROUTINE magMoms(dimension,input,atoms,noco,vTot,chmom,qa21)
SUBROUTINE magMoms(dimension,input,atoms,noco,vTot,stdn,svdn,chmom,qa21)
USE m_types
USE m_xmlOutput
......@@ -24,11 +24,28 @@ SUBROUTINE magMoms(dimension,input,atoms,noco,vTot,chmom,qa21)
REAL, INTENT(INOUT) :: chmom(atoms%ntype,dimension%jspd)
COMPLEX, INTENT(IN) :: qa21(atoms%ntype)
REAL, INTENT(IN) :: stdn(atoms%ntype,dimension%jspd)
REAL, INTENT(IN) :: svdn(atoms%ntype,dimension%jspd)
INTEGER :: iType, j, iRepAtom
REAL :: smom
REAL :: sval,stot,scor,smom
CHARACTER(LEN=20) :: attributes(4)
WRITE (6,FMT=8000)
WRITE (16,FMT=8000)
DO iType = 1,atoms%ntype
sval = svdn(iType,1) - svdn(iType,input%jspins)
stot = stdn(iType,1) - stdn(iType,input%jspins)
scor = stot - sval
WRITE (6,FMT=8010) iType,stot,sval,scor,svdn(iType,1),stdn(iType,1)
WRITE (16,FMT=8010) iType,stot,sval,scor,svdn(iType,1),stdn(iType,1)
END DO
8000 FORMAT (/,/,10x,'spin density at the nucleus:',/,10x,'type',t25,&
'input%total',t42,'valence',t65,'core',t90,&
'majority valence and input%total density',/)
8010 FORMAT (i13,2x,3e20.8,5x,2e20.8)
WRITE (6,FMT=8020)
WRITE (16,FMT=8020)
......
......@@ -11,8 +11,8 @@ MODULE m_qfix
! qfix file no longer supported!
CONTAINS
SUBROUTINE qfix( stars,atoms,sym,vacuum,&
sphhar,input,cell,oneD,den,l_printData,force_fix,fix)
SUBROUTINE qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
den,l_noco,l_printData,force_fix,fix)
USE m_types
USE m_cdntot
......@@ -29,7 +29,7 @@ CONTAINS
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_potden),INTENT(INOUT) :: den
LOGICAL,INTENT(IN) :: l_printData,force_fix
LOGICAL,INTENT(IN) :: l_noco,l_printData,force_fix
REAL, INTENT (OUT) :: fix
! .. Local Scalars ..
......@@ -51,14 +51,13 @@ CONTAINS
! qfix==0 means no qfix was given in inp.xml.
! In this case do nothing except when forced to fix!
CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
den%pw,den%mt,den%vacz,.TRUE., qtot,qis)
CALL cdntot(stars,atoms,sym, vacuum,input,cell,oneD,&
den%pw,den%mt,den%vacz,.TRUE., qtot,qis)
!The total nucleii charge
zc=SUM(atoms%neq(:)*atoms%zatom(:))
zc = zc + 2*input%efield%sigma
IF (fixtotal) THEN
!-roa
fix = zc/qtot
......@@ -76,23 +75,27 @@ CONTAINS
den%vacxy(:vacuum%nmzxy,:oneD%odi%nq2-1,:vacuum%nvac,:)
END IF
WRITE (6,FMT=8000) zc,fix
IF (ABS(fix-1.0)<1.E-6) RETURN !no second calculation of cdntot as nothing was fixed
CALL openXMLElementNoAttributes('fixedCharges')
CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
den%pw,den%mt,den%vacz,l_printData, qtot,qis)
CALL closeXMLElement('fixedCharges')
!+roa
ELSE
fix = (zc - qtot) / qis + 1.
den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
WRITE (6,FMT=8001) zc,fix
IF (ABS(fix-1.0)<1.E-6) RETURN !no second calculation of cdntot as nothing was fixed
CALL openXMLElementNoAttributes('fixedCharges')
CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
den%pw,den%mt,den%vacz,l_printData, qtot,qis)
CALL closeXMLElement('fixedCharges')
ENDIF
IF (l_noco) THEN
!fix also the off-diagonal part of the density matrix
den%pw(:stars%ng3,3) = fix*den%pw(:stars%ng3,3)
IF (input%film.AND.fixtotal) THEN
den%vacz(:,:,3:4) = fix*den%vacz(:,:,3:4)
den%vacxy(:,:,:,3) = fix*den%vacxy(:,:,:,3)
END IF
END IF
IF (ABS(fix-1.0)<1.E-6) RETURN !no second calculation of cdntot as nothing was fixed
CALL openXMLElementNoAttributes('fixedCharges')
CALL cdntot(stars,atoms,sym, vacuum,input,cell,oneD,&
den%pw,den%mt,den%vacz,l_printData, qtot,qis)
CALL closeXMLElement('fixedCharges')
IF (fix>1.1) CALL juDFT_WARN("You lost too much charge")
IF (fix<.9) CALL juDFT_WARN("You gained too much charge")
......
......@@ -70,7 +70,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
TYPE(t_noco) :: noco_new
!Local Scalars
REAL fix,qtot,scor,stot,sval,dummy
REAL fix,qtot,dummy
INTEGER ivac,j,jspin,jspmax,k,iType
LOGICAL l_enpara
......@@ -82,7 +82,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
!pk non-collinear (start)
INTEGER igq2_fft(0:stars%kq1_fft*stars%kq2_fft-1)
COMPLEX,ALLOCATABLE :: qa21(:), cdomvz(:,:)
COMPLEX,ALLOCATABLE :: qa21(:)
!pk non-collinear (end)
IF (mpi%irank.EQ.0) THEN
......@@ -133,55 +133,6 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
stars,cell,sphhar,atoms,vTot,outDen,stdn,svdn)
IF (mpi%irank.EQ.0) THEN
!block 2 unnecessary for slicing: begin
IF (.NOT.sliceplot%slice) THEN
CALL openXMLElementNoAttributes('allElectronCharges')
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,outDen,.TRUE.,.true.,fix)
CALL closeXMLElement('allElectronCharges')
!pk non-collinear (start)
IF (noco%l_noco) THEN
!fix also the off-diagonal part of the density matrix
outDen%pw(:stars%ng3,3) = fix*outDen%pw(:stars%ng3,3)
IF (input%film) THEN
outDen%vacz(:,:,3:4) = fix*outDen%vacz(:,:,3:4)
outDen%vacxy(:,:,:,3) = fix*outDen%vacxy(:,:,:,3)
END IF
END IF
!pk non-collinear (end)
!spin densities at the nucleus
!and magnetic moment in the spheres
IF (input%jspins.EQ.2) THEN
WRITE (6,FMT=8000)
WRITE (16,FMT=8000)
DO iType = 1,atoms%ntype
sval = svdn(iType,1) - svdn(iType,input%jspins)
stot = stdn(iType,1) - stdn(iType,input%jspins)
scor = stot - sval
WRITE (6,FMT=8010) iType,stot,sval,scor,svdn(iType,1),stdn(iType,1)
WRITE (16,FMT=8010) iType,stot,sval,scor,svdn(iType,1),stdn(iType,1)
END DO
noco_new = noco
CALL magMoms(dimension,input,atoms,noco_new,vTot,chmom,qa21)
!Generate and save the new nocoinp file if the directions of the local
!moments are relaxed or a constraint B-field is calculated.
IF (ANY(noco%l_relax(:atoms%ntype)).OR.noco%l_constr) THEN
CALL genNewNocoInp(input,atoms,jij,noco,noco_new)
END IF
IF (noco%l_soc) CALL orbMagMoms(dimension,atoms,noco,clmom)
END IF
!block 2 unnecessary for slicing: end
END IF ! .NOT.sliceplot%slice
8000 FORMAT (/,/,10x,'spin density at the nucleus:',/,10x,'type',t25,&
'input%total',t42,'valence',t65,'core',t90,&
'majority valence and input%total density',/)
8010 FORMAT (i13,2x,3e20.8,5x,2e20.8)
IF(vacuum%nvac.EQ.1) THEN
outDen%vacz(:,2,:) = outDen%vacz(:,1,:)
IF (sym%invs) THEN
......@@ -195,11 +146,31 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,-1.0,0.0,.FALSE.,outDen,'cdn_slice')
END IF
ENDIF ! mpi%irank.EQ.0
END IF
IF (sliceplot%slice) CALL juDFT_end("slice OK",mpi%irank)
IF (mpi%irank.EQ.0) THEN
CALL openXMLElementNoAttributes('allElectronCharges')
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,outDen,noco%l_noco,.TRUE.,.true.,fix)
CALL closeXMLElement('allElectronCharges')
IF (input%jspins.EQ.2) THEN
noco_new = noco
!Calculate and write out spin densities at the nucleus and magnetic moments in the spheres
CALL magMoms(dimension,input,atoms,noco_new,vTot,stdn,svdn,chmom,qa21)
!Generate and save the new nocoinp file if the directions of the local
!moments are relaxed or a constraint B-field is calculated.
IF (ANY(noco%l_relax(:atoms%ntype)).OR.noco%l_constr) THEN
CALL genNewNocoInp(input,atoms,jij,noco,noco_new)
END IF
IF (noco%l_soc) CALL orbMagMoms(dimension,atoms,noco,clmom)
END IF
END IF ! mpi%irank.EQ.0
#ifdef CPP_MPI
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,outDen)
#endif
......@@ -207,8 +178,6 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
DEALLOCATE (qvac,qvlay,qa21)
DEALLOCATE (igq_fft)
RETURN
END SUBROUTINE cdngen
END MODULE m_cdngen
......@@ -184,7 +184,7 @@ CONTAINS
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,fermiEnergyTemp,l_qfix,inDen)
CALL timestart("Qfix")
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,.FALSE.,.false.,fix)
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false.,fix)
CALL timestop("Qfix")
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,-1.0,0.0,.FALSE.,inDen)
......
......@@ -243,7 +243,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
CALL closeXMLElement('densityConvergence')
!fix charge of the new density
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,.FALSE.,.false., fix)
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false., fix)
IF(atoms%n_u.NE.n_u_keep) THEN
inDen%mmpMat = n_mmpTemp
......
......@@ -153,7 +153,7 @@ SUBROUTINE pldngen(sym,stars,atoms,sphhar,vacuum,&
den%vacz(:,:,4) = AIMAG(cdomvz(:,:))
den%vacxy(:,:,:,3) = cdomvxy
END IF
CALL qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,.FALSE.,.true.,fix)
CALL qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,noco%l_noco,.FALSE.,.true.,fix)
rho(:,0:,1:,:input%jspins) = den%mt(:,0:,1:,:input%jspins)
qpw(1:,:input%jspins) = den%pw(1:,:input%jspins)
rht(1:,1:,:input%jspins) = den%vacz(1:,1:,:input%jspins)
......
......@@ -204,7 +204,7 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,&
IF (mpi%irank == 0) THEN
! Check the normalization of total density
CALL qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,.FALSE.,.true.,fix)
CALL qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,.FALSE.,.FALSE.,.true.,fix)
z=SUM(atoms%neq(:)*atoms%zatom(:))
IF (ABS(fix*z-z)>0.5) THEN
CALL judft_warn("Starting density not charge neutral",hint= &
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment