qfix.f90 3.72 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
CONTAINS
14
  SUBROUTINE qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
15
                  den,l_noco,l_printData,force_fix,fix)
16

Daniel Wortmann's avatar
Daniel Wortmann committed
17 18
    USE m_types
    USE m_cdntot
19
    USE m_xmlOutput
Daniel Wortmann's avatar
Daniel Wortmann committed
20
    IMPLICIT NONE
Gregor Michalicek's avatar
Gregor Michalicek committed
21

Daniel Wortmann's avatar
Daniel Wortmann committed
22
    !     .. Scalar Arguments ..
23
    TYPE(t_mpi),INTENT(IN)       :: mpi
Gregor Michalicek's avatar
Gregor Michalicek committed
24 25 26 27 28 29 30 31 32
    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
33
    LOGICAL,INTENT(IN)           :: l_noco,l_printData,force_fix
Gregor Michalicek's avatar
Gregor Michalicek committed
34 35
    REAL,    INTENT (OUT)        :: fix

Daniel Wortmann's avatar
Daniel Wortmann committed
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
    !     .. 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!
    
55
    CALL cdntot(mpi,stars,atoms,sym,vacuum,input,cell,oneD,den,.TRUE.,qtot,qis)
56

Daniel Wortmann's avatar
Daniel Wortmann committed
57 58
    !The total nucleii charge
    zc=SUM(atoms%neq(:)*atoms%zatom(:))
59
    zc = zc + 2*input%sigma
60

Daniel Wortmann's avatar
Daniel Wortmann committed
61 62 63 64 65 66 67
    IF (fixtotal) THEN
       !-roa
       fix = zc/qtot
       na = 1
       DO n = 1,atoms%ntype
          lh = sphhar%nlh(atoms%ntypsy(na))
          jm = atoms%jri(n)
Gregor Michalicek's avatar
Gregor Michalicek committed
68
          den%mt(:jm,0:lh,n,:) = fix*den%mt(:jm,0:lh,n,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
69 70
          na = na + atoms%neq(n)
       ENDDO
Gregor Michalicek's avatar
Gregor Michalicek committed
71
       den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
72
       IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
73
          den%vacz(:vacuum%nmz,:vacuum%nvac,:) = fix*den%vacz(:vacuum%nmz,:vacuum%nvac,:)
74 75
          den%vacxy(:vacuum%nmzxy,:stars%ng2-1,:vacuum%nvac,:) = fix*&
             den%vacxy(:vacuum%nmzxy,:stars%ng2-1,:vacuum%nvac,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
76 77 78 79
       END IF
       WRITE (6,FMT=8000) zc,fix
    ELSE
       fix = (zc - qtot) / qis + 1.
Gregor Michalicek's avatar
Gregor Michalicek committed
80
       den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
81 82
       WRITE (6,FMT=8001) zc,fix
    ENDIF
83

84 85 86 87 88 89 90 91 92 93
    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
94

95
    CALL openXMLElementNoAttributes('fixedCharges')
96
    CALL cdntot(mpi,stars,atoms,sym,vacuum,input,cell,oneD,den,l_printData,qtot,qis)
97 98
    CALL closeXMLElement('fixedCharges')

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

Daniel Wortmann's avatar
Daniel Wortmann committed
102 103 104 105 106 107 108

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