qfix.f90 3.63 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,fix_pw_only)
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
    REAL,    INTENT (OUT)        :: fix
35
    LOGICAL,INTENT(IN),OPTIONAL  :: fix_pw_only
Daniel Wortmann's avatar
Daniel Wortmann committed
36 37 38 39 40 41 42
    !     .. 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
43
    IF (PRESENT(fix_pw_only)) fixtotal=.NOT.fix_pw_only
Daniel Wortmann's avatar
Daniel Wortmann committed
44
    fix=1.0
45
    IF (l_firstcall) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
46 47
       l_firstcall=.false.
    ELSE
48
       IF (MOD(input%qfix,2)==0.AND..NOT.force_fix) RETURN
Daniel Wortmann's avatar
Daniel Wortmann committed
49 50 51 52
    ENDIF
    ! qfix==0 means no qfix was given in inp.xml. 
    ! In this case do nothing except when forced to fix!
    
Matthias Redies's avatar
Matthias Redies committed
53
    CALL cdntot(stars,atoms,sym,vacuum,input,cell,oneD,den,.TRUE.,qtot,qis)
54

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

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

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

93
    CALL openXMLElementNoAttributes('fixedCharges')
Matthias Redies's avatar
Matthias Redies committed
94
    CALL cdntot(stars,atoms,sym,vacuum,input,cell,oneD,den,l_printData,qtot,qis)
95 96
    CALL closeXMLElement('fixedCharges')

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