LCOV - code coverage report
Current view: top level - global - qfix.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 32 37 86.5 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          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             : 
      13             : CONTAINS
      14         387 :   SUBROUTINE qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
      15             :                   den,l_noco,l_printData,force_fix,fix,fix_pw_only)
      16             : 
      17             :     USE m_types
      18             :     USE m_cdntot
      19             :     USE m_xmlOutput
      20             :     IMPLICIT NONE
      21             : 
      22             :     !     .. Scalar Arguments ..
      23             :     TYPE(t_mpi),INTENT(IN)       :: mpi
      24             :     TYPE(t_stars),INTENT(IN)     :: stars
      25             :     TYPE(t_atoms),INTENT(IN)     :: atoms
      26             :     TYPE(t_sym),INTENT(IN)       :: sym
      27             :     TYPE(t_vacuum),INTENT(IN)    :: vacuum
      28             :     TYPE(t_sphhar),INTENT(IN)    :: sphhar
      29             :     TYPE(t_input),INTENT(IN)     :: input
      30             :     TYPE(t_oneD),INTENT(IN)      :: oneD
      31             :     TYPE(t_cell),INTENT(IN)      :: cell
      32             :     TYPE(t_potden),INTENT(INOUT) :: den
      33             :     LOGICAL,INTENT(IN)           :: l_noco,l_printData,force_fix
      34             :     REAL,    INTENT (OUT)        :: fix
      35             :     LOGICAL,INTENT(IN),OPTIONAL  :: fix_pw_only
      36             :     !     .. Local Scalars ..
      37             :     LOGICAL :: l_qfixfile,fixtotal
      38             :     LOGICAL :: l_firstcall=.true.
      39             :     REAL    :: qtot,qis,zc
      40             :     INTEGER :: jm,lh,n,na
      41             :     !     ..
      42         387 :     fixtotal=.true. !this is the default
      43         387 :     IF (PRESENT(fix_pw_only)) fixtotal=.NOT.fix_pw_only
      44         387 :     fix=1.0
      45         387 :     IF (l_firstcall) THEN
      46          36 :        l_firstcall=.false.
      47             :     ELSE
      48         735 :        IF (MOD(input%qfix,2)==0.AND..NOT.force_fix) RETURN
      49             :     ENDIF
      50             :     ! qfix==0 means no qfix was given in inp.xml. 
      51             :     ! In this case do nothing except when forced to fix!
      52             :     
      53         201 :     CALL cdntot(stars,atoms,sym,vacuum,input,cell,oneD,den,.TRUE.,qtot,qis)
      54             : 
      55             :     !The total nucleii charge
      56         201 :     zc=SUM(atoms%neq(:)*atoms%zatom(:))
      57         201 :     zc = zc + 2*input%sigma
      58             : 
      59         201 :     IF (fixtotal) THEN
      60             :        !-roa
      61         201 :        fix = zc/qtot
      62         201 :        na = 1
      63         609 :        DO n = 1,atoms%ntype
      64         408 :           lh = sphhar%nlh(atoms%ntypsy(na))
      65         408 :           jm = atoms%jri(n)
      66         408 :           den%mt(:jm,0:lh,n,:) = fix*den%mt(:jm,0:lh,n,:)
      67         609 :           na = na + atoms%neq(n)
      68             :        ENDDO
      69         201 :        den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
      70         201 :        IF (input%film) THEN
      71           9 :           den%vacz(:vacuum%nmz,:vacuum%nvac,:) = fix*den%vacz(:vacuum%nmz,:vacuum%nvac,:)
      72             :           den%vacxy(:vacuum%nmzxy,:stars%ng2-1,:vacuum%nvac,:) = fix*&
      73           9 :              den%vacxy(:vacuum%nmzxy,:stars%ng2-1,:vacuum%nvac,:)
      74             :        END IF
      75         201 :        WRITE (6,FMT=8000) zc,fix
      76             :     ELSE
      77           0 :        fix = (zc - qtot) / qis + 1.
      78           0 :        den%pw(:stars%ng3,:) = fix*den%pw(:stars%ng3,:)
      79           0 :        WRITE (6,FMT=8001) zc,fix
      80             :     ENDIF
      81             : 
      82         201 :     IF (l_noco) THEN
      83             :        !fix also the off-diagonal part of the density matrix
      84         125 :        den%pw(:stars%ng3,3) = fix*den%pw(:stars%ng3,3)
      85         125 :        IF (input%film.AND.fixtotal) THEN
      86           0 :           den%vacz(:,:,3:4) = fix*den%vacz(:,:,3:4)
      87           0 :           den%vacxy(:,:,:,3) = fix*den%vacxy(:,:,:,3)
      88             :        END IF
      89             :     END IF
      90             : 
      91         201 :     IF (ABS(fix-1.0)<1.E-6) RETURN !no second calculation of cdntot as nothing was fixed
      92             : 
      93           3 :     CALL openXMLElementNoAttributes('fixedCharges')
      94           3 :     CALL cdntot(stars,atoms,sym,vacuum,input,cell,oneD,den,l_printData,qtot,qis)
      95           3 :     CALL closeXMLElement('fixedCharges')
      96             : 
      97           3 :     IF (fix>1.1) CALL juDFT_WARN("You lost too much charge")
      98           3 :     IF (fix<.9) CALL juDFT_WARN("You gained too much charge")
      99             : 
     100             : 
     101             : 8000 FORMAT (/,10x,'zc= ',f12.6,5x,'qfix=  ',f10.6)
     102             : 8001 FORMAT (/,' > broy only qis: ','zc= ',f12.6,5x,'qfix=  ',f10.6)
     103             :     !-roa
     104             : 
     105             :   END SUBROUTINE qfix
     106             : END MODULE m_qfix

Generated by: LCOV version 1.13