qfix.f90 3.75 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
CONTAINS
  SUBROUTINE qfix( stars,atoms,sym,vacuum,&
Gregor Michalicek's avatar
Gregor Michalicek committed
15
       sphhar,input,cell,oneD,den,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 32 33 34
    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

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 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!
    
    CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
Gregor Michalicek's avatar
Gregor Michalicek committed
55
         den%pw,den%mt,den%vacz,.TRUE., qtot,qis)
56

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

Daniel Wortmann's avatar
Daniel Wortmann committed
61 62 63 64 65 66 67 68
   
    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
69
          den%mt(:jm,0:lh,n,:) = fix*den%mt(:jm,0:lh,n,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
70 71
          na = na + atoms%neq(n)
       ENDDO
Gregor Michalicek's avatar
Gregor Michalicek committed
72
       den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
73
       IF (input%film) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
74 75 76
          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,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
77 78 79
       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
80
       CALL openXMLElementNoAttributes('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
81
       CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
Gregor Michalicek's avatar
Gregor Michalicek committed
82
            den%pw,den%mt,den%vacz,l_printData, qtot,qis)
83
       CALL closeXMLElement('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
84 85 86
       !+roa 
    ELSE
       fix = (zc - qtot) / qis + 1.
Gregor Michalicek's avatar
Gregor Michalicek committed
87
       den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
88 89
       WRITE (6,FMT=8001) zc,fix
       IF (ABS(fix-1.0)<1.E-6) RETURN !no second calculation of cdntot as nothing was fixed
90
       CALL openXMLElementNoAttributes('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
91
       CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
Gregor Michalicek's avatar
Gregor Michalicek committed
92
            den%pw,den%mt,den%vacz,l_printData, qtot,qis)
93
       CALL closeXMLElement('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
94
    ENDIF
95

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

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

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