Commit a991918f authored by Daniel Wortmann's avatar Daniel Wortmann

Changed qfix behaviour:

- do not use qfix file anymore
- do not run qfix if it is not specified in inp.xml
Fixes #96
parent b0b3033f
MODULE m_qfix MODULE m_qfix
USE m_juDFT USE m_juDFT
! ******************************************************* !Calculate total charge
! check total charge and renormalize c,l.fu !Depending on variable input%qfix, the following will be done to fix the charge
! ******************************************************* !Input qfix can be 1 or 2
CONTAINS ! qfix=0 (no qfix in inp.xml) means we usually do not run the code
SUBROUTINE qfix(& ! if force_fix is .true. we run the code and assume qfix=2
& stars,atoms,sym,vacuum,& ! in the call to qfix we will always run it
& sphhar,input,cell,oneD,& ! qfix=1 (qfix=f in inp.xml) means we fix only in INT (only done in firstcall)
& qpw,rhtxy,rho,rht,l_printData,& ! qfix=2 (qfix=t in inp.xml) means we fix total charge
& fix) ! qfix file no longer supported!
USE m_types CONTAINS
USE m_cdntot SUBROUTINE qfix( stars,atoms,sym,vacuum,&
IMPLICIT NONE sphhar,input,cell,oneD,qpw,rhtxy,rho,rht,l_printData,force_fix,fix)
! ..
! .. Scalar Arguments ..
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_vacuum),INTENT(IN):: vacuum
TYPE(t_sphhar),INTENT(IN):: sphhar
TYPE(t_input),INTENT(IN) :: input
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
LOGICAL,INTENT(IN) :: l_printData
REAL, INTENT (OUT) :: fix
!-odim
!+odim
! ..
! .. Array Arguments ..
COMPLEX,INTENT (INOUT) :: qpw(stars%ng3,input%jspins)
COMPLEX,INTENT (INOUT) :: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,input%jspins)
REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
REAL, INTENT (INOUT) :: rht(vacuum%nmzd,2,input%jspins)
! ..
! .. Local Scalars ..
LOGICAL fixtot,l99
REAL qtot,qis,zc
INTEGER ivac,j,jm,jspin,k,lh,n,nl,ns,na
! ..
CALL cdntot(&
& stars,atoms,sym,&
& vacuum,input,cell,oneD,&
& qpw,rho,rht,.FALSE.,&
& qtot,qis)
zc = 0.
DO 10 n = 1,atoms%ntype
zc = zc + atoms%neq(n)*atoms%zatom(n)
10 CONTINUE
!+roa (check via qfix file if total charge or only interstitial to fix) USE m_types
fixtot=.TRUE. USE m_cdntot
INQUIRE(file='qfix',exist=l99) IMPLICIT NONE
IF (l99) then ! ..
OPEN (99,file='qfix',status='old',form='formatted') ! .. Scalar Arguments ..
READ (99,'(1x,l1)',end=1199) fixtot TYPE(t_stars),INTENT(IN) :: stars
IF (.NOT.fixtot ) THEN TYPE(t_atoms),INTENT(IN) :: atoms
REWIND (99) TYPE(t_sym),INTENT(IN) :: sym
WRITE (99,'(1x,l1,70a)') .TRUE.,& TYPE(t_vacuum),INTENT(IN):: vacuum
& ' (1x,l1) F..fix interstitial T..fix total charge ' TYPE(t_sphhar),INTENT(IN):: sphhar
ENDIF TYPE(t_input),INTENT(IN) :: input
1199 CLOSE (99) TYPE(t_oneD),INTENT(IN) :: oneD
ENDIF TYPE(t_cell),INTENT(IN) :: cell
zc = zc + 2*input%efield%sigma LOGICAL,INTENT(IN) :: l_printData,force_fix
IF ( fixtot ) THEN REAL, INTENT (OUT) :: fix
!-roa ! ..
fix = zc/qtot ! .. Array Arguments ..
DO 100 jspin = 1,input%jspins COMPLEX,INTENT (INOUT) :: qpw(stars%ng3,input%jspins)
na = 1 COMPLEX,INTENT (INOUT) :: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,input%jspins)
DO 40 n = 1,atoms%ntype REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
ns = atoms%ntypsy(na) REAL, INTENT (INOUT) :: rht(vacuum%nmzd,2,input%jspins)
lh = sphhar%nlh(ns) ! ..
jm = atoms%jri(n) ! .. Local Scalars ..
DO 30 nl = 0,lh LOGICAL :: l_qfixfile,fixtotal
DO 20 j = 1,jm LOGICAL :: l_firstcall=.true.
rho(j,nl,n,jspin) = fix*rho(j,nl,n,jspin) REAL :: qtot,qis,zc
20 CONTINUE INTEGER :: jm,lh,n,na
30 CONTINUE ! ..
na = na + atoms%neq(n) fixtotal=.true. !this is the default
40 CONTINUE fix=1.0
DO 50 k = 1,stars%ng3 if (l_firstcall) THEN
qpw(k,jspin) = fix*qpw(k,jspin) INQUIRE(file='qfix',exist=l_qfixfile)
50 CONTINUE IF (l_qfixfile) CALL judft_info("qfix file no longer supported, check the qfix option in inp.xml","INFO")
IF (input%film) THEN if (input%qfix==1) fixtotal=.false.
DO 90 ivac = 1,vacuum%nvac l_firstcall=.false.
DO 60 n = 1,vacuum%nmz ELSE
rht(n,ivac,jspin) = fix*rht(n,ivac,jspin) IF (input%qfix==0.AND..NOT.force_fix) RETURN
60 CONTINUE ENDIF
DO 80 n = 1,vacuum%nmzxy ! qfix==0 means no qfix was given in inp.xml.
DO 70 k = 2,oneD%odi%nq2 ! In this case do nothing except when forced to fix!
rhtxy(n,k-1,ivac,jspin) = fix*&
& rhtxy(n,k-1,ivac,jspin) CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
70 CONTINUE qpw,rho,rht,.FALSE., qtot,qis)
80 CONTINUE
90 CONTINUE
END IF
100 CONTINUE
WRITE (6,FMT=8000) zc,fix
CALL cdntot(&
& stars,atoms,sym,&
& vacuum,input,cell,oneD,&
& qpw,rho,rht,l_printData,&
& qtot,qis)
!+roa
ELSE
fix = (zc - qtot) / qis + 1.
DO jspin = 1,input%jspins
DO k = 1,stars%ng3
qpw(k,jspin) = fix*qpw(k,jspin)
ENDDO
ENDDO
WRITE (6,FMT=8001) zc,fix
CALL cdntot(&
& stars,atoms,sym,&
& vacuum,input,cell,oneD,&
& qpw,rho,rht,l_printData,&
& qtot,qis)
ENDIF !The total nucleii charge
zc=SUM(atoms%neq(:)*atoms%zatom(:))
zc = zc + 2*input%efield%sigma
IF (fix>1.1) CALL juDFT_WARN("You lost too much charge")
IF (fix<.9) CALL juDFT_WARN("You gained too much charge") IF (fixtotal) THEN
!-roa
fix = zc/qtot
na = 1
DO n = 1,atoms%ntype
lh = sphhar%nlh(atoms%ntypsy(na))
jm = atoms%jri(n)
rho(:jm,0:lh,n,:) = fix*rho(:jm,0:lh,n,:)
na = na + atoms%neq(n)
ENDDO
qpw(:stars%ng3,:) = fix*qpw(:stars%ng3,:)
IF (input%film) THEN
rht(:vacuum%nmz,:vacuum%nvac,:) = fix*rht(:vacuum%nmz,:vacuum%nvac,:)
rhtxy(:vacuum%nmzxy,:oneD%odi%nq2-1,:vacuum%nvac,:) = fix*&
rhtxy(: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 cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
qpw,rho,rht,l_printData, qtot,qis)
!+roa
ELSE
fix = (zc - qtot) / qis + 1.
qpw(:stars%ng3,:) = fix*qpw(: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 cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
qpw,rho,rht,l_printData, qtot,qis)
ENDIF
8000 FORMAT (/,10x,'zc= ',f12.6,5x,'qfix= ',f10.6) IF (fix>1.1) CALL juDFT_WARN("You lost too much charge")
8001 FORMAT (/,' > broy only qis: ','zc= ',f12.6,5x,'qfix= ',f10.6) IF (fix<.9) CALL juDFT_WARN("You gained too much charge")
!-roa
END SUBROUTINE qfix
END MODULE m_qfix 8000 FORMAT (/,10x,'zc= ',f12.6,5x,'qfix= ',f10.6)
8001 FORMAT (/,' > broy only qis: ','zc= ',f12.6,5x,'qfix= ',f10.6)
!-roa
END SUBROUTINE qfix
END MODULE m_qfix
...@@ -97,7 +97,9 @@ ...@@ -97,7 +97,9 @@
IF (ch_rw.eq.'r') THEN IF (ch_rw.eq.'r') THEN
!-------------------------------------------------------------------- !--------------------------------------------------------------------
OPEN (5,file='inp',form='formatted',status='old') OPEN (5,file='inp',form='formatted',status='old')
!default not read in in old inp-file
input%qfix=2
! !
a1(:) = 0 a1(:) = 0
a2(:) = 0 a2(:) = 0
......
...@@ -345,7 +345,7 @@ ...@@ -345,7 +345,7 @@
! block 2 unnecessary for slicing: begin ! block 2 unnecessary for slicing: begin
IF (.NOT.sliceplot%slice) THEN IF (.NOT.sliceplot%slice) THEN
CALL openXMLElementNoAttributes('allElectronCharges') CALL openXMLElementNoAttributes('allElectronCharges')
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD, qpw,rhtxy,rho,rht,.TRUE., fix) CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD, qpw,rhtxy,rho,rht,.TRUE.,.true., fix)
CALL closeXMLElement('allElectronCharges') CALL closeXMLElement('allElectronCharges')
!---> pk non-collinear !---> pk non-collinear
IF (noco%l_noco) THEN IF (noco%l_noco) THEN
......
...@@ -302,7 +302,7 @@ CONTAINS ...@@ -302,7 +302,7 @@ CONTAINS
! !
! ---> fix the new density ! ---> fix the new density
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,& CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,&
qpw,rhtxy,rho,rht,.FALSE., fix) qpw,rhtxy,rho,rht,.FALSE.,.false., fix)
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,& CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,results%last_distance,results%ef,.TRUE.,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy) 1,results%last_distance,results%ef,.TRUE.,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
......
...@@ -164,7 +164,7 @@ CONTAINS ...@@ -164,7 +164,7 @@ CONTAINS
IF (.NOT.l_xyav) THEN IF (.NOT.l_xyav) THEN
CALL timestart("Qfix") CALL timestart("Qfix")
CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD, qpw,rhtxy,rho,rht,.FALSE., fix) CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD, qpw,rhtxy,rho,rht,.FALSE.,.false., fix)
CALL timestop("Qfix") CALL timestop("Qfix")
ENDIF ENDIF
......
...@@ -155,7 +155,7 @@ ...@@ -155,7 +155,7 @@
CALL qfix(& CALL qfix(&
& stars,atoms,sym,vacuum,& & stars,atoms,sym,vacuum,&
& sphhar,input,cell,oneD,& & sphhar,input,cell,oneD,&
& qpw,rhtxy,rho,rht,.FALSE.,& & qpw,rhtxy,rho,rht,.FALSE.,.true.,&
& fix) & fix)
ENDIF ENDIF
......
...@@ -237,7 +237,7 @@ ...@@ -237,7 +237,7 @@
CALL qfix(& CALL qfix(&
& stars,atoms,sym,vacuum,& & stars,atoms,sym,vacuum,&
& sphhar,input,cell,oneD,& & sphhar,input,cell,oneD,&
& qpw,rhtxy,rho,rht,.FALSE.,& & qpw,rhtxy,rho,rht,.FALSE.,.true.,&
& fix) & fix)
! !
! Write superposed density onto density file ! Write superposed density onto density file
......
...@@ -122,7 +122,7 @@ CONTAINS ...@@ -122,7 +122,7 @@ CONTAINS
CALL qfix(& CALL qfix(&
& stars,atoms,sym,vacuum,& & stars,atoms,sym,vacuum,&
& sphhar,input,cell,oneD,& & sphhar,input,cell,oneD,&
& qpw,rhtxy,rho,rht,.FALSE.,& & qpw,rhtxy,rho,rht,.FALSE.,.FALSE.,&
& fix) & fix)
!---> fouriertransform the diagonal part of the density matrix !---> fouriertransform the diagonal part of the density matrix
......
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