qfix.f90 3.76 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3 4 5 6 7 8 9 10 11
MODULE m_qfix
  USE m_juDFT
  !Calculate total charge
  !Depending on variable input%qfix, the following will be done to fix the charge
  !Input qfix can be 1 or 2
  ! qfix=0 (no qfix in inp.xml) means we usually do not run the code
  !                       if force_fix is .true. we run the code and assume qfix=2
  !                       in the call to qfix we will always run it
  ! qfix=1 (qfix=f in inp.xml) means we fix only in INT (only done in firstcall)
  ! qfix=2 (qfix=t in inp.xml) means we fix total charge
  ! qfix file no longer supported!
12

Daniel Wortmann's avatar
Daniel Wortmann committed
13 14 15
CONTAINS
  SUBROUTINE qfix( stars,atoms,sym,vacuum,&
       sphhar,input,cell,oneD,qpw,rhtxy,rho,rht,l_printData,force_fix,fix)
16

Daniel Wortmann's avatar
Daniel Wortmann committed
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
    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
    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)
    !     ..
    !     .. Local Scalars ..
    LOGICAL :: l_qfixfile,fixtotal
    LOGICAL :: l_firstcall=.true.
    REAL    :: qtot,qis,zc
    INTEGER :: jm,lh,n,na
    !     ..
    fixtotal=.true. !this is the default
    fix=1.0
    if (l_firstcall) THEN
       INQUIRE(file='qfix',exist=l_qfixfile)
       IF (l_qfixfile) CALL judft_info("qfix file no longer supported, check the qfix option in inp.xml","INFO")
       if (input%qfix==1) fixtotal=.false.
       l_firstcall=.false.
    ELSE
       IF (input%qfix==0.AND..NOT.force_fix) RETURN
    ENDIF
    ! 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,&
         qpw,rho,rht,.FALSE., qtot,qis)
60

Daniel Wortmann's avatar
Daniel Wortmann committed
61 62 63
    !The total nucleii charge
    zc=SUM(atoms%neq(:)*atoms%zatom(:))
    zc = zc + 2*input%efield%sigma
64

Daniel Wortmann's avatar
Daniel Wortmann committed
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
   
    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)
94

Daniel Wortmann's avatar
Daniel Wortmann committed
95
    ENDIF
96

Daniel Wortmann's avatar
Daniel Wortmann committed
97 98
    IF (fix>1.1) CALL juDFT_WARN("You lost too much charge")
    IF (fix<.9) CALL juDFT_WARN("You gained too much charge")
99

Daniel Wortmann's avatar
Daniel Wortmann committed
100 101 102 103 104 105 106

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