qfix.f90 3.67 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 15
  SUBROUTINE qfix(stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
                  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 ..
Gregor Michalicek's avatar
Gregor Michalicek committed
23 24 25 26 27 28 29 30 31
    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
32
    LOGICAL,INTENT(IN)           :: l_noco,l_printData,force_fix
Gregor Michalicek's avatar
Gregor Michalicek committed
33 34
    REAL,    INTENT (OUT)        :: fix

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

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

Daniel Wortmann's avatar
Daniel Wortmann committed
60 61 62 63 64 65 66
    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
67
          den%mt(:jm,0:lh,n,:) = fix*den%mt(:jm,0:lh,n,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
68 69
          na = na + atoms%neq(n)
       ENDDO
Gregor Michalicek's avatar
Gregor Michalicek committed
70
       den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
71
       IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
72
          den%vacz(:vacuum%nmz,:vacuum%nvac,:) = fix*den%vacz(:vacuum%nmz,:vacuum%nvac,:)
73 74
          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
75 76 77 78
       END IF
       WRITE (6,FMT=8000) zc,fix
    ELSE
       fix = (zc - qtot) / qis + 1.
Gregor Michalicek's avatar
Gregor Michalicek committed
79
       den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
80 81
       WRITE (6,FMT=8001) zc,fix
    ENDIF
82

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

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

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

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

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