LCOV - code coverage report
Current view: top level - main - totale.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 71 80 88.8 %
Date: 2019-09-08 04:53:50 Functions: 1 3 33.3 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : MODULE m_totale
       7             : CONTAINS
       8         324 :   SUBROUTINE totale(mpi,atoms,sphhar,stars,vacuum,DIMENSION, &
       9             :        sym,input,noco,cell,oneD, xcpot,hybrid,vTot,vCoul,it,den,results)
      10             :     !
      11             :     !     ***************************************************
      12             :     !     subroutine calculates the total energy 
      13             :     !     ***************************************************
      14             :     !     single particle energies
      15             :     !     SEIGC  sum of the eigenvalues of the core states
      16             :     !            calculated in cdngen.f
      17             :     !     SEIGSCV  sum of the eigenvalues of the semicore and valence states
      18             :     !              calculated in fermie.f 
      19             :     !     TS         : entropy contribution to the free energy
      20             :     !     SEIGC,SEIGSCV, TS are calculated in fermie.f
      21             :     !     ***************************************************
      22             :     !     TE_VCOUL  :   charge density-coulomb potential integral
      23             :     !     TE_VEFF:   charge density-effective potential integral
      24             :     !     TE_EXC :   charge density-ex-corr.energy density integral
      25             :     !                 exchange-correlation energy
      26             :     !     TE_VCOUL,TE_VEFF,TE_EXC are calculated in vgen.f
      27             :     !     VMD :   Madelung term
      28             :     !     ***************************************************
      29             :     !     TOTE    :   total energy due to all electrons
      30             :     !     TOTE = SEIGC + SEIGSCV + TE_VCOUL/2 -TE_VEFF + TE_EXC + VMD
      31             :     !
      32             :     !     if HF calculation/hybrid-functional calculation :
      33             :     !     TOTE = SEIGC + SEIGSCV + TE_VCOUL/2 -TE_VEFF + TE_EXC_loc + VMD - 1/2 E_FOCK
      34             :     !
      35             :     !     E_FOCK: sum of diagonal elements of fock matrix
      36             :     !
      37             :     !     ***************************************************
      38             :     !     FREE ENRGY: F = TOTE - TS
      39             :     !     total electron energy extrapolated for T->0
      40             :     !     E0 = TOTE - TS/2
      41             :     !     ***************************************************
      42             :     !
      43             :     USE m_intgr    , ONLY : intgr3 
      44             :     USE m_constants
      45             :     USE m_force_a4
      46             :     USE m_force_a3
      47             :     USE m_forcew
      48             :     USE m_cdn_io
      49             :     USE m_types
      50             :     USE m_xmlOutput
      51             :     IMPLICIT NONE
      52             :     TYPE(t_mpi),INTENT(IN)          :: mpi
      53             :     TYPE(t_results),INTENT(INOUT)   :: results
      54             :     CLASS(t_xcpot),INTENT(IN)       :: xcpot
      55             :     TYPE(t_oneD),INTENT(IN)         :: oneD
      56             :     TYPE(t_hybrid),INTENT(IN)       :: hybrid
      57             :     TYPE(t_input),INTENT(IN)        :: input
      58             :     TYPE(t_vacuum),INTENT(IN)       :: vacuum
      59             :     TYPE(t_noco),INTENT(IN)         :: noco
      60             :     TYPE(t_sym),INTENT(IN)          :: sym
      61             :     TYPE(t_stars),INTENT(IN)        :: stars
      62             :     TYPE(t_cell),INTENT(IN)         :: cell
      63             :     TYPE(t_sphhar),INTENT(IN)       :: sphhar
      64             :     TYPE(t_atoms),INTENT(IN)        :: atoms
      65             :     TYPE(t_dimension),INTENT(IN)    :: dimension
      66             :     TYPE(t_potden),INTENT(IN)       :: vTot,vCoul
      67             :     TYPE(t_potden),INTENT(IN)       :: den
      68             :     !     ..
      69             :     !     .. Scalar Arguments ..
      70             :     INTEGER,INTENT (IN) :: it      
      71             : 
      72             :     ! Local type instances
      73             : 
      74             :     !     .. Local Scalars ..
      75             :     REAL rhs,totz, eigSum, fermiEnergyTemp
      76             :     INTEGER n,j,nt,i, archiveType
      77             :     LOGICAL l_qfix
      78             : 
      79             :     !     .. Local Arrays ..
      80         646 :     REAL vmd(atoms%ntype),zintn_r(atoms%ntype)
      81         646 :     REAL dpj(atoms%jmtd),mt(atoms%jmtd,atoms%ntype)
      82             :     CHARACTER(LEN=20) :: attributes(3)
      83             : 
      84             :     !CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,.FALSE.,POTDEN_TYPE_DEN)
      85         324 :     IF (mpi%irank==0) THEN
      86         162 :        WRITE (6,FMT=8000)
      87             : 8000   FORMAT (/,/,/,5x,'t o t a l  e n e r g y')
      88             :        !
      89             :        !      ---> sum of eigenvalues (core, semicore and valence states)
      90             :        !
      91         162 :        eigSum = results%seigscv + results%seigc
      92         162 :        results%tote = eigSum
      93         162 :        WRITE (6,FMT=8010) results%tote
      94             : 8010   FORMAT (/,10x,'sum of eigenvalues =',t40,f20.10)
      95             :        !
      96             :        !      ---> add contribution of coulomb potential
      97             :        !
      98         162 :        results%tote = results%tote + 0.5e0*results%te_vcoul
      99         162 :        WRITE (6,FMT=8020) results%te_vcoul
     100             : 8020   FORMAT (/,10x,'density-coulomb potential integral =',t40,f20.10)
     101             :        !
     102             :        !      ---> add contribution of effective potential
     103             :        !
     104         162 :        results%tote = results%tote - results%te_veff
     105         162 :        WRITE (6,FMT=8030) results%te_veff
     106             : 8030   FORMAT (/,10x,'density-effective potential integral =',t40,f20.10)
     107             :        !
     108             :        !      ---> add contribution of exchange-correlation energy
     109             :        !
     110         162 :        results%tote = results%tote + results%te_exc
     111         162 :        WRITE (6,FMT=8040) results%te_exc
     112             : 8040   FORMAT (/,10x,'charge density-ex.-corr.energy density integral=', t40,f20.10)
     113             :        !
     114             :        !      ---> Fock exchange contribution 
     115             :        !
     116         162 :        IF (xcpot%is_hybrid()) THEN
     117             :           !IF (xcpot%is_name("exx")) THEN
     118             :           !   results%tote = results%tote + 0.5e0*results%te_hfex%valence
     119             :           !ELSE
     120           0 :           results%tote = results%tote - 0.5e0*results%te_hfex%valence + 0.5e0*results%te_hfex%core
     121             :           !END IF
     122           0 :           WRITE (6,FMT=8100)  0.5e0*results%te_hfex%valence
     123           0 :           WRITE (6,FMT=8101)  0.5e0*results%te_hfex%core
     124             : 8100      FORMAT (/,10x,'Fock-exchange energy (valence)=',t40,f20.10)
     125             : 8101      FORMAT (10x,'Fock-exchange energy (core)=',t40,f20.10)
     126             :        ENDIF
     127             : 
     128             : 
     129             :        !     ----> VM terms
     130             :        !     ---> reload the density
     131             :        !
     132             :        !archiveType = CDN_ARCHIVE_TYPE_CDN1_const
     133             :        !IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_CDN_const
     134             : 
     135             :        !CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
     136             :        !                 CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den)
     137             : 
     138             : 
     139             :        ! CLASSICAL HELLMAN-FEYNMAN FORCE
     140         162 :        CALL force_a3(atoms,sphhar, input, den%mt,vCoul%mt, results%force)
     141             : 
     142         162 :        IF (input%l_f) THEN
     143             :           ! core contribution to force: needs TOTAL POTENTIAL and core charge
     144           2 :           CALL force_a4(atoms,sphhar,input,dimension, vTot%mt, results%force)
     145             : 
     146             :        ENDIF
     147             : 
     148             :        !-for
     149             :        !     ---> add spin-up and spin-down charge density for lh=0
     150             :        !
     151         474 :        mt=0.0
     152         474 :        DO  n = 1,atoms%ntype
     153      187332 :           DO  i = 1,atoms%jri(n)
     154      187170 :              mt(i,n) = den%mt(i,0,n,1) + den%mt(i,0,n,input%jspins)
     155             :           ENDDO
     156             :        ENDDO
     157         162 :        IF (input%jspins.EQ.1) mt=mt/2 !we just added the same value twice
     158             : 
     159             :        !
     160             :        ! ----> coulomb interaction between electrons and nuclei of different m.t.s
     161             :        !
     162         474 :        DO  n = 1,atoms%ntype
     163      187170 :           DO  j = 1,atoms%jri(n)
     164      187170 :              dpj(j) = mt(j,n)/atoms%rmsh(j,n)
     165             :           ENDDO
     166         312 :           CALL intgr3(dpj,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),rhs)
     167             :           !
     168         312 :           zintn_r(n) = atoms%neq(n)*atoms%zatom(n)*sfp_const*rhs/2.
     169         312 :           results%tote = results%tote - zintn_r(n)
     170             :           !
     171         312 :           WRITE (6,FMT=8045) zintn_r(n)
     172         312 :           CALL intgr3(mt(1,n),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),totz)
     173         312 :           vmd(n) = atoms%rmt(n)*vCoul%mt(atoms%jri(n),0,n,1)/sfp_const + atoms%zatom(n) - totz*sfp_const
     174         312 :           vmd(n) = -atoms%neq(n)*atoms%zatom(n)*vmd(n)/ (2.*atoms%rmt(n))
     175         312 :           WRITE (6,FMT=8050) n,vmd(n)
     176         474 :           results%tote = results%tote + vmd(n)
     177             :        ENDDO
     178         162 :        IF (atoms%n_u.GT.0) THEN
     179          15 :           WRITE ( 6,FMT=8090) results%e_ldau
     180          15 :           results%tote = results%tote - results%e_ldau             ! gu test
     181             :        ENDIF
     182             :        ! print 'HF' before total energy to make it grepable
     183         162 :        IF ( .NOT. hybrid%l_calhf ) THEN
     184         162 :           WRITE ( 6,FMT=8060) results%tote
     185             :        ELSE
     186           0 :           WRITE ( 6,FMT=8061) results%tote
     187             :        END IF
     188             :        !
     189             :        !     ---> calculate the free energy and the ground state energy,
     190             :        !          extrapolated for T->0
     191             :        !
     192             :        ! print 'HF' before all energies to make them grepable
     193         162 :        IF ( .NOT. hybrid%l_calhf ) THEN
     194         162 :           WRITE ( 6,FMT=8065) results%ts
     195         162 :           WRITE ( 6,FMT=8070) results%tote-results%ts
     196         162 :           WRITE ( 6,FMT=8080) results%tote-0.5e0*results%ts
     197             :        ELSE
     198           0 :           WRITE ( 6,FMT=8066) results%ts
     199           0 :           WRITE ( 6,FMT=8071) results%tote-results%ts
     200           0 :           WRITE ( 6,FMT=8081) results%tote-0.5e0*results%ts
     201             :        END IF
     202             : 
     203         162 :        WRITE(attributes(1),'(f20.10)') results%tote
     204         162 :        WRITE(attributes(2),'(a)') 'Htr'
     205         162 :        WRITE(attributes(3),'(a)') 'HF'
     206         162 :        IF (hybrid%l_calhf) THEN
     207           0 :           CALL openXMLElementForm('totalEnergy',(/'value  ','units  ','comment'/),attributes,reshape((/40,20/),(/1,2/)))
     208             :        ELSE
     209         162 :           CALL openXMLElementForm('totalEnergy',(/'value','units'/),attributes(1:2),reshape((/40,20/),(/1,2/)))
     210             :        END IF
     211         162 :        CALL openXMLElementFormPoly('sumOfEigenvalues',(/'value'/),(/eigSum/),reshape((/32,20/),(/1,2/)))
     212         162 :        CALL writeXMLElementFormPoly('coreElectrons',(/'value'/),(/results%seigc/),reshape((/32,20/),(/1,2/)))
     213         162 :        CALL writeXMLElementFormPoly('valenceElectrons',(/'value'/),(/results%seigscv/),reshape((/29,20/),(/1,2/)))
     214         162 :        CALL closeXMLElement('sumOfEigenvalues')
     215         162 :        CALL writeXMLElementFormPoly('densityCoulombPotentialIntegral',(/'value'/),(/results%te_vcoul/),reshape((/17,20/),(/1,2/)))
     216         162 :        CALL writeXMLElementFormPoly('densityEffectivePotentialIntegral',(/'value'/),(/results%te_veff/),reshape((/15,20/),(/1,2/)))
     217         162 :        CALL writeXMLElementFormPoly('chargeDenXCDenIntegral',(/'value'/),(/results%te_exc/),reshape((/26,20/),(/1,2/)))
     218         162 :        CALL writeXMLElementFormPoly('FockExchangeEnergyValence',(/'value'/),(/0.5e0*results%te_hfex%valence/),reshape((/23,20/),(/1,2/)))
     219         162 :        CALL writeXMLElementFormPoly('FockExchangeEnergyCore',(/'value'/),(/0.5e0*results%te_hfex%core/),reshape((/26,20/),(/1,2/)))
     220         474 :        DO  n = 1,atoms%ntype
     221         312 :           CALL openXMLElementPoly('atomTypeDependentContributions',(/'atomType'/),(/n/))
     222         312 :           CALL writeXMLElementFormPoly('electronNucleiInteractionDifferentMTs',(/'value'/),(/zintn_r(n)/),reshape((/8,20/),(/1,2/)))
     223         312 :           CALL writeXMLElementFormPoly('MadelungTerm',(/'value'/),(/vmd(n)/),reshape((/33,20/),(/1,2/)))
     224         474 :           CALL closeXMLElement('atomTypeDependentContributions')
     225             :        END DO
     226         162 :        IF (atoms%n_u.GT.0) THEN
     227          15 :           CALL writeXMLElementFormPoly('dftUCorrection',(/'value'/),(/results%e_ldau/),reshape((/34,20/),(/1,2/)))
     228             :        END IF
     229         162 :        CALL writeXMLElementFormPoly('tkbTimesEntropy',(/'value'/),(/results%ts/),reshape((/33,20/),(/1,2/)))
     230         162 :        CALL writeXMLElementFormPoly('freeEnergy',(/'value'/),(/results%tote-results%ts/),reshape((/38,20/),(/1,2/)))
     231         162 :        CALL writeXMLElementFormPoly('extrapolationTo0K',(/'value'/),(/results%tote-0.5e0*results%ts/),reshape((/31,20/),(/1,2/)))
     232         162 :        CALL closeXMLElement('totalEnergy')
     233             : 8060   FORMAT (/,/,' ---->    total energy=',t40,f20.10,' htr')
     234             : 8061   FORMAT (/,/,' ----> HF total energy=',t40,f20.10,' htr')
     235             : 8050   FORMAT (/,10x,'Madelung term for atom type:',i3,t40,f20.10)
     236             : 8045   FORMAT (/,10x,'el.-nucl. inter. diff. m.t.',t40,f20.10)
     237             : 8065   FORMAT (/,/,' ---->    (input%tkb*entropy) TS=',t40,f20.10,' htr')
     238             : 8066   FORMAT (/,/,' ----> HF (input%tkb*entropy) TS=',t40,f20.10,' htr')
     239             : 8070   FORMAT (/,/,' ---->    free energy=',t40,f20.10,' htr')
     240             : 8071   FORMAT (/,/,' ----> HF free energy=',t40,f20.10,' htr')
     241             : 8080   FORMAT (/,/,'      extrapolation for T->0',&
     242             :             /,' ---->    total electron energy=',t40,f20.10,' htr')
     243             : 8081   FORMAT (/,/,'      extrapolation for T->0',&
     244             :             /,' ----> HF total electron energy=',t40,f20.10,' htr')
     245             : 8090   FORMAT (/,/,' ---->    correction for lda+U =',t40,f20.10,' htr')
     246             :     ENDIF
     247         324 :     CALL force_w(mpi,input,atoms,sym,results,cell,oneD,vacuum)
     248             : 
     249         322 :   END SUBROUTINE totale
     250           0 : END MODULE m_totale

Generated by: LCOV version 1.13