qfix.f90 3.97 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
    USE m_types
    USE m_cdntot
19
    USE m_xmlOutput
Daniel Wortmann's avatar
Daniel Wortmann committed
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
    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,&
60
         qpw,rho,rht,.TRUE., qtot,qis)
61

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

Daniel Wortmann's avatar
Daniel Wortmann committed
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
   
    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
85
       CALL openXMLElementNoAttributes('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
86 87
       CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
            qpw,rho,rht,l_printData, qtot,qis)
88
       CALL closeXMLElement('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
89 90 91 92 93 94
       !+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
95
       CALL openXMLElementNoAttributes('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
96 97
       CALL cdntot( stars,atoms,sym, vacuum,input,cell,oneD,&
            qpw,rho,rht,l_printData, qtot,qis)
98
       CALL closeXMLElement('fixedCharges')
Daniel Wortmann's avatar
Daniel Wortmann committed
99
    ENDIF
100

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

Daniel Wortmann's avatar
Daniel Wortmann committed
104 105 106 107 108 109 110

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