qfix.f90 4.16 KB
Newer Older
1 2 3 4 5 6 7 8 9
      MODULE m_qfix
      USE m_juDFT
!     *******************************************************
!     check total charge and renormalize        c,l.fu
!     *******************************************************
      CONTAINS
      SUBROUTINE qfix(&
     &                stars,atoms,sym,vacuum,&
     &                sphhar,input,cell,oneD,&
10
     &                qpw,rhtxy,rho,rht,l_printData,&
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
     &                fix)

      USE m_types
      USE m_cdntot
      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
26
      LOGICAL,INTENT(IN)       :: l_printData
27 28 29 30 31
      REAL,    INTENT (OUT) :: fix
!-odim
!+odim
!     ..
!     .. Array Arguments ..
32
      COMPLEX,INTENT (INOUT) :: qpw(stars%ng3,input%jspins)
33
      COMPLEX,INTENT (INOUT) :: rhtxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,input%jspins)
Daniel Wortmann's avatar
Daniel Wortmann committed
34
      REAL,   INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
35 36 37 38 39 40 41 42 43 44
      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,&
45
     &            qpw,rho,rht,.FALSE.,&
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
     &            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)
      fixtot=.TRUE.
      INQUIRE(file='qfix',exist=l99)
      IF (l99) then
        OPEN (99,file='qfix',status='old',form='formatted')
        READ (99,'(1x,l1)',end=1199) fixtot
        IF (.NOT.fixtot ) THEN
           REWIND (99)
           WRITE (99,'(1x,l1,70a)') .TRUE.,&
     &      ' (1x,l1) F..fix interstitial T..fix total charge '
        ENDIF
 1199   CLOSE (99)
      ENDIF
      zc = zc + 2*input%efield%sigma
      IF ( fixtot ) THEN
!-roa
         fix = zc/qtot
         DO 100 jspin = 1,input%jspins
            na = 1
            DO 40 n = 1,atoms%ntype
               ns = atoms%ntypsy(na)
               lh = sphhar%nlh(ns)
               jm = atoms%jri(n)
               DO 30 nl = 0,lh
                  DO 20 j = 1,jm
                     rho(j,nl,n,jspin) = fix*rho(j,nl,n,jspin)
   20             CONTINUE
   30          CONTINUE
               na = na + atoms%neq(n)
   40       CONTINUE
            DO 50 k = 1,stars%ng3
               qpw(k,jspin) = fix*qpw(k,jspin)
   50       CONTINUE
            IF (input%film) THEN
               DO 90 ivac = 1,vacuum%nvac
                  DO 60 n = 1,vacuum%nmz
                     rht(n,ivac,jspin) = fix*rht(n,ivac,jspin)
   60             CONTINUE
                  DO 80 n = 1,vacuum%nmzxy
                     DO 70 k = 2,oneD%odi%nq2
                        rhtxy(n,k-1,ivac,jspin) = fix*&
     &                    rhtxy(n,k-1,ivac,jspin)
   70                CONTINUE
   80             CONTINUE
   90          CONTINUE
            END IF
  100    CONTINUE
         WRITE (6,FMT=8000) zc,fix
         CALL cdntot(&
     &               stars,atoms,sym,&
     &               vacuum,input,cell,oneD,&
103
     &               qpw,rho,rht,l_printData,&
104 105 106 107 108 109 110 111 112 113 114 115 116
     &               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,&
117
     &               qpw,rho,rht,l_printData,&
118 119 120 121 122 123 124 125 126 127 128 129 130 131
     &               qtot,qis)

      ENDIF

      IF (fix>1.1) CALL juDFT_WARN("You lost too much charge")
      IF (fix<.9) CALL juDFT_WARN("You gained too much charge")


 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