Line data    Source code

       1             : MODULE m_qfix
       2             :   USE m_juDFT
       3             :   !Calculate total charge
       4             :   !Depending on variable input%qfix, the following will be done to fix the charge
       5             :   !Input qfix can be 1 or 2
       6             :   ! qfix=0 (no qfix in inp.xml) means we usually do not run the code
       7             :   !                       if force_fix is .true. we run the code and assume qfix=2
       8             :   !                       in the call to qfix we will always run it
       9             :   ! qfix=1 (qfix=f in inp.xml) means we fix only in INT (only done in firstcall)
      10             :   ! qfix=2 (qfix=t in inp.xml) means we fix total charge
      11             :   ! qfix file no longer supported!
      12             :   ! If l_par=.TRUE., MPI parallelization in the cdntot will be used.
      13             :   ! Be carefull not to set it to .TRUE. if you are calling only from one MPI
      14             :   ! rank.
      15             : 
      16             : CONTAINS
      17        2126 :   SUBROUTINE qfix(fmpi,stars,nococonv,atoms,sym,vacuum,sphhar,input,cell ,&
      18             :                   den,l_noco,l_printData,l_par,force_fix,fix,fix_pw_only)
      19             : 
      20             :     USE m_types
      21             :     USE m_constants
      22             :     USE m_cdntot
      23             :     USE m_xmlOutput
      24             : 
      25             :     IMPLICIT NONE
      26             : 
      27             :     !     .. Scalar Arguments ..
      28             :     TYPE(t_mpi),INTENT(IN)       :: fmpi
      29             :     TYPE(t_nococonv),INTENT(IN)  :: nococonv
      30             :     TYPE(t_stars),INTENT(IN)     :: stars
      31             :     TYPE(t_atoms),INTENT(IN)     :: atoms
      32             :     TYPE(t_sym),INTENT(IN)       :: sym
      33             :     TYPE(t_vacuum),INTENT(IN)    :: vacuum
      34             :     TYPE(t_sphhar),INTENT(IN)    :: sphhar
      35             :     TYPE(t_input),INTENT(IN)     :: input
      36             : 
      37             :     TYPE(t_cell),INTENT(IN)      :: cell
      38             :     TYPE(t_potden),INTENT(INOUT) :: den
      39             :     LOGICAL,INTENT(IN)           :: l_noco,l_printData,force_fix,l_par
      40             :     REAL,    INTENT (OUT)        :: fix
      41             :     LOGICAL,INTENT(IN),OPTIONAL  :: fix_pw_only
      42             :     !     .. Local Scalars ..
      43             :     LOGICAL :: l_qfixfile,fixtotal
      44             :     LOGICAL :: l_firstcall=.true.
      45             :     REAL    :: qtot,qis,zc
      46             :     INTEGER :: jm,lh,n,na
      47             :     !     ..
      48        1295 :     fixtotal=.true. !this is the default
      49        1295 :     IF (PRESENT(fix_pw_only)) fixtotal=.NOT.fix_pw_only
      50        1295 :     fix=1.0
      51        1295 :     IF (l_firstcall) THEN
      52         158 :        l_firstcall=.false.
      53             :     ELSE
      54        1536 :        IF (MOD(input%qfix,2)==0.AND..NOT.force_fix) RETURN
      55             :     ENDIF
      56             :     ! qfix==0 means no qfix was given in inp.xml.
      57             :     ! In this case do nothing except when forced to fix!
      58             : 
      59         831 :     CALL cdntot(stars,nococonv,atoms,sym,vacuum,input,cell ,den,l_printData,qtot,qis,fmpi,l_par)
      60             : 
      61         831 :     IF (fmpi%irank.EQ.0) THEN
      62             :        !The total nucleii charge
      63        1158 :        zc=SUM(atoms%neq(:)*atoms%zatom(:))
      64             :        !zc = zc + 2*input%sigma !TODO : reactivate fields
      65             : 
      66         417 :        IF (fixtotal) THEN
      67             :           !-roa
      68         417 :           fix = zc/qtot
      69        1158 :           DO n = 1,atoms%ntype
      70         741 :              na = atoms%firstAtom(n)
      71         741 :              lh = sphhar%nlh(sym%ntypsy(na))
      72         741 :              jm = atoms%jri(n)
      73    28934282 :              den%mt(:jm,0:lh,n,:) = fix*den%mt(:jm,0:lh,n,:)
      74             :           ENDDO
      75     2035979 :           den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
      76         417 :           IF (input%film) THEN
      77     6747106 :              den%vac(:vacuum%nmz,:stars%ng2,:vacuum%nvac,:) = fix*den%vac(:vacuum%nmz,:stars%ng2,:vacuum%nvac,:)
      78             :           END IF
      79         417 :           WRITE (oUnit,FMT=8000) zc,fix
      80             :        ELSE
      81           0 :           fix = (zc - qtot) / qis + 1.
      82           0 :           den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
      83           0 :           WRITE (oUnit,FMT=8001) zc,fix
      84             :        ENDIF
      85             : 
      86             :        ! TODO: This looks spooky.
      87             :        ! a) All noco quantities are already included in the fix above.
      88             :        ! b) Below, MT is missing.
      89             :        !IF (l_noco) THEN
      90             :           !fix also the off-diagonal part of the density matrix
      91             :           !den%pw(:stars%ng3,3) = fix*den%pw(:stars%ng3,3)
      92             :           !IF (input%film.AND.fixtotal) THEN
      93             :              !den%vacz(:,:,3:4) = fix*den%vacz(:,:,3:4)
      94             :              !den%vacxy(:,:,:,3) = fix*den%vacxy(:,:,:,3)
      95             :           !END IF
      96             :        !END IF
      97             : 
      98         417 :        IF (ABS(fix-1.0)<1.E-6) RETURN !no second calculation of cdntot as nothing was fixed
      99             : 
     100          18 :        IF(l_printData) CALL openXMLElementNoAttributes('fixedCharges')
     101          18 :        CALL cdntot(stars,nococonv,atoms,sym,vacuum,input,cell ,den,l_printData,qtot,qis,fmpi,.FALSE.)
     102          18 :        IF(l_printData) CALL closeXMLElement('fixedCharges')
     103             : 
     104          18 :        IF (fix>1.1) CALL juDFT_WARN("You lost too much charge")
     105          18 :        IF (fix<.9) CALL juDFT_WARN("You gained too much charge")
     106             : 
     107             : 8000   FORMAT (/,10x,'zc= ',f12.6,5x,'qfix=  ',f10.6)
     108             : 8001   FORMAT (/,' > broy only qis: ','zc= ',f12.6,5x,'qfix=  ',f10.6)
     109             :        !-roa
     110             :     ENDIF
     111             : 
     112             :   END SUBROUTINE qfix
     113             : END MODULE m_qfix