Use potden type in qfix

...in optional/pldngen.f90 this has to be strongly polished.

global/qfix.f90

@@ -12,31 +12,26 @@ MODULE m_qfix
 
 CONTAINS
   SUBROUTINE qfix( stars,atoms,sym,vacuum,&
-       sphhar,input,cell,oneD,qpw,rhtxy,rho,rht,l_printData,force_fix,fix)
+       sphhar,input,cell,oneD,den,l_printData,force_fix,fix)
 
     USE m_types
     USE m_cdntot
     USE m_xmlOutput
     IMPLICIT NONE
-    !     ..
+
     !     .. 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,force_fix
-    REAL,    INTENT (OUT)    :: fix
-    !     ..
-    !     .. 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)
-    !     ..
+    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
+    TYPE(t_potden),INTENT(INOUT) :: den
+    LOGICAL,INTENT(IN)           :: l_printData,force_fix
+    REAL,    INTENT (OUT)        :: fix
+
     !     .. Local Scalars ..
     LOGICAL :: l_qfixfile,fixtotal
     LOGICAL :: l_firstcall=.true.
@@ -57,7 +52,7 @@ CONTAINS
     ! In this case do nothing except when forced to fix!
     
     CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
-         qpw,rho,rht,.TRUE., qtot,qis)
+         den%pw,den%mt,den%vacz,.TRUE., qtot,qis)
 
     !The total nucleii charge
     zc=SUM(atoms%neq(:)*atoms%zatom(:))
@@ -71,30 +66,30 @@ CONTAINS
        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,:)
+          den%mt(:jm,0:lh,n,:) = fix*den%mt(:jm,0:lh,n,:)
           na = na + atoms%neq(n)
        ENDDO
-       qpw(:stars%ng3,:) = fix*qpw(:stars%ng3,:)
+       den%pw(:stars%ng3,:) = fix*den%pw(: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,:)
+          den%vacz(:vacuum%nmz,:vacuum%nvac,:) = fix*den%vacz(:vacuum%nmz,:vacuum%nvac,:)
+          den%vacxy(:vacuum%nmzxy,:oneD%odi%nq2-1,:vacuum%nvac,:) = fix*&
+             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,&
-            qpw,rho,rht,l_printData, qtot,qis)
+            den%pw,den%mt,den%vacz,l_printData, qtot,qis)
        CALL closeXMLElement('fixedCharges')
        !+roa 
     ELSE
        fix = (zc - qtot) / qis + 1.
-       qpw(:stars%ng3,:) = fix*qpw(:stars%ng3,:)
+       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,&
-            qpw,rho,rht,l_printData, qtot,qis)
+            den%pw,den%mt,den%vacz,l_printData, qtot,qis)
        CALL closeXMLElement('fixedCharges')
     ENDIF
 

main/cdngen.F90

@@ -146,8 +146,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
       !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%pw,outDen%vacxy,outDen%mt,outDen%vacz,.TRUE.,.true.,fix)
+         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

main/fleur.F90

@@ -185,8 +185,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%pw,inDen%vacxy,inDen%mt,inDen%vacz,&
-                 .FALSE.,.false.,fix)
+       CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,.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)

main/mix.F90

@@ -241,8 +241,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%pw,inDen%vacxy,inDen%mt,inDen%vacz,.FALSE.,.false., fix)
+   CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,.FALSE.,.false., fix)
 
    IF(atoms%n_u.NE.n_u_keep) THEN
       inDen%mmpMat = n_mmpTemp

optional/pldngen.f90

@@ -147,8 +147,28 @@ SUBROUTINE pldngen(sym,stars,atoms,sphhar,vacuum,&
    WRITE(6,*)'Assuming collinear magnetization.'
    150 CLOSE (nrhomfile)
    IF (.NOT. sliceplot%slice) THEN 
-      CALL qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
-                qpw,rhtxy,rho,rht,.FALSE.,.true.,fix)
+      CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,inp%jspins,noco%l_noco,POTDEN_TYPE_DEN)
+      den%iter = iter
+      den%mt(:,0:,1:,:input%jspins) = rho(:,0:,1:,:input%jspins)
+      den%pw(1:,:input%jspins) = qpw(1:,:input%jspins)
+      den%vacz(1:,1:,:input%jspins) = rht(1:,1:,:input%jspins)
+      den%vacxy(1:,1:,1:,:input%jspins) = rhtxy(1:,1:,1:,:input%jspins)
+      IF(noco%l_noco) THEN
+         den%pw(:,3) = cdom
+         den%vacz(:,:,3) = REAL(cdomvz(:,:))
+         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)
+      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)
+      rhtxy(1:,1:,1:,:input%jspins) = den%vacxy(1:,1:,1:,:input%jspins)
+      IF(noco%l_noco) THEN
+         cdom = den%pw(:,3)
+         cdomvz(:,:) = CMPLX(den%vacz(:,:,3),den%vacz(:,:,4))
+         cdomvxy = den%vacxy(:,:,:,3)
+      END IF
    END IF
 
    !---> for testing: read offdiag. output density matrix

optional/stden.f90

@@ -204,8 +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%pw,den%vacxy,den%mt,den%vacz,.FALSE.,.true.,fix)
+      CALL qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,.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= &