LCOV - code coverage report
Current view: top level - rdmft - rdmft.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 365 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2018 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             : 
       7             : MODULE m_rdmft
       8             : 
       9             : CONTAINS
      10             : 
      11           0 : SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars,vacuum,dimension,&
      12             :                  sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,hybrid,results,coreSpecInput,archiveType,outDen)
      13             : 
      14             :    USE m_types
      15             :    USE m_juDFT
      16             :    USE m_constants
      17             :    USE m_eig66_io
      18             : #ifndef CPP_OLDINTEL
      19             :    USE m_cdnval
      20             :    USE m_cdngen
      21             :    USE m_cdn_io
      22             :    USE m_cdncore
      23             :    USE m_qfix
      24             :    USE m_vgen_coulomb
      25             :    USE m_convol
      26             :    USE m_intnv
      27             : 
      28             :    USE m_mixedbasis
      29             :    USE m_coulombmatrix
      30             :    USE m_hf_init
      31             :    USE m_hf_setup
      32             :    USE m_io_hybrid
      33             :    USE m_symm_hf
      34             :    USE m_exchange_valence_hf
      35             :    USE m_exchange_core
      36             :    USE m_symmetrizeh
      37             :    USE m_bfgs_b2
      38             : 
      39             : #ifdef CPP_MPI
      40             :    USE m_mpi_bc_potden
      41             : #endif
      42             : #endif
      43             : 
      44             :    IMPLICIT NONE
      45             : 
      46             :    TYPE(t_mpi),           INTENT(IN)    :: mpi
      47             :    TYPE(t_input),         INTENT(IN)    :: input
      48             :    TYPE(t_kpts),          INTENT(IN)    :: kpts
      49             :    TYPE(t_banddos),       INTENT(IN)    :: banddos
      50             :    TYPE(t_sliceplot),     INTENT(IN)    :: sliceplot
      51             :    TYPE(t_cell),          INTENT(IN)    :: cell
      52             :    TYPE(t_atoms),         INTENT(IN)    :: atoms
      53             :    TYPE(t_enpara),        INTENT(INOUT) :: enpara
      54             :    TYPE(t_stars),         INTENT(IN)    :: stars
      55             :    TYPE(t_vacuum),        INTENT(IN)    :: vacuum
      56             :    TYPE(t_dimension),     INTENT(IN)    :: dimension
      57             :    TYPE(t_sphhar),        INTENT(IN)    :: sphhar
      58             :    TYPE(t_sym),           INTENT(IN)    :: sym
      59             :    TYPE(t_field),         INTENT(INOUT) :: field
      60             :    TYPE(t_potden),        INTENT(INOUT) :: vTot
      61             :    TYPE(t_potden),        INTENT(INOUT) :: vCoul
      62             :    TYPE(t_oneD),          INTENT(IN)    :: oneD
      63             :    TYPE(t_noco),          INTENT(INOUT) :: noco
      64             :    TYPE(t_xcpot_inbuild), INTENT(INOUT) :: xcpot
      65             :    TYPE(t_hybrid),        INTENT(INOUT) :: hybrid
      66             :    TYPE(t_results),       INTENT(INOUT) :: results
      67             :    TYPE(t_coreSpecInput), INTENT(IN)    :: coreSpecInput
      68             :    TYPE(t_potden),        INTENT(INOUT) :: outDen
      69             : 
      70             :    INTEGER,               INTENT(IN)    :: eig_id
      71             :    INTEGER,               INTENT(IN)    :: archiveType
      72           0 :    TYPE(t_potden)                       :: EnergyDen
      73             : 
      74             : #ifndef CPP_OLDINTEL
      75           0 :    TYPE(t_cdnvalJob)                    :: cdnvalJob
      76           0 :    TYPE(t_potden)                       :: singleStateDen, overallDen, overallVCoul, vTotTemp
      77           0 :    TYPE(t_regionCharges)                :: regCharges
      78           0 :    TYPE(t_dos)                          :: dos
      79           0 :    TYPE(t_moments)                      :: moments
      80           0 :    TYPE(t_mat)                          :: exMat, zMat, olap, trafo, invtrafo, tmpMat, exMatLAPW
      81           0 :    TYPE(t_lapw)                         :: lapw
      82           0 :    TYPE(t_hybdat)                       :: hybdat
      83             :    INTEGER                              :: ikpt,ikpt_i,iband_i, iBand, jkpt, jBand, iAtom, i, na, itype, lh, j
      84             :    INTEGER                              :: jspin, jspmax, jsp, isp, ispin, nbasfcn, nbands
      85             :    INTEGER                              :: nsymop, nkpt_EIBZ, ikptf, iterHF, mnobd
      86             :    INTEGER                              :: iState, iStep, numStates, maxHistoryLength, numRelevantStates
      87             :    REAL                                 :: fix, potDenInt, fermiEnergyTemp, spinDegenFac
      88             :    REAL                                 :: rdmftFunctionalValue, occStateI, gradSum
      89             :    REAL                                 :: exchangeTerm, lagrangeMultiplier, equalityCriterion
      90             :    REAL                                 :: mixParam, rdmftEnergy
      91             :    REAL                                 :: sumOcc, tempOcc, addCharge, subCharge, addChargeWeight, subChargeWeight
      92             :    REAL, PARAMETER                      :: degenEps = 0.00001
      93             :    REAL, PARAMETER                      :: convCrit = 1.0e-6
      94             :    REAL, PARAMETER                      :: minOcc = 1.0e-8
      95             :    LOGICAL                              :: converged, l_qfix, l_restart, l_zref
      96             :    CHARACTER(LEN=20)                    :: filename
      97             : 
      98           0 :    INTEGER                              :: nsest(dimension%neigd) ! probably too large
      99           0 :    INTEGER                              :: indx_sest(dimension%neigd,dimension%neigd) ! probably too large
     100           0 :    INTEGER                              :: rrot(3,3,sym%nsym)
     101           0 :    INTEGER                              :: psym(sym%nsym) ! Note: psym is only filled up to index nsymop
     102           0 :    INTEGER                              :: lowestState(kpts%nkpt,input%jspins)
     103           0 :    INTEGER                              :: highestState(kpts%nkpt,input%jspins)
     104           0 :    INTEGER                              :: neigTemp(kpts%nkpt,input%jspins)
     105             : 
     106           0 :    REAL                                 :: wl_iks(dimension%neigd,kpts%nkptf)
     107             : 
     108           0 :    REAL, ALLOCATABLE                    :: overallVCoulSSDen(:,:,:)
     109           0 :    REAL, ALLOCATABLE                    :: vTotSSDen(:,:,:)
     110           0 :    REAL, ALLOCATABLE                    :: dEdOcc(:,:,:)
     111             : 
     112           0 :    REAL, ALLOCATABLE                    :: exDiag(:,:,:)
     113           0 :    REAL, ALLOCATABLE                    :: eig_irr(:,:)
     114             : 
     115           0 :    REAL, ALLOCATABLE                    :: gradient(:), lastGradient(:)
     116           0 :    REAL, ALLOCATABLE                    :: parameters(:), lastParameters(:)
     117           0 :    REAL, ALLOCATABLE                    :: minConstraints(:)
     118           0 :    REAL, ALLOCATABLE                    :: maxConstraints(:)
     119           0 :    REAL, ALLOCATABLE                    :: gradientCorrections(:,:)
     120           0 :    REAL, ALLOCATABLE                    :: paramCorrections(:,:)
     121           0 :    REAL, ALLOCATABLE                    :: equalityLinCombi(:)
     122             : 
     123           0 :    REAL, ALLOCATABLE                    :: occupationVec(:)
     124             : 
     125           0 :    INTEGER, ALLOCATABLE                 :: parent(:)
     126           0 :    INTEGER, ALLOCATABLE                 :: pointer_EIBZ(:)
     127           0 :    INTEGER, ALLOCATABLE                 :: n_q(:)
     128             : 
     129           0 :    LOGICAL, ALLOCATABLE                 :: enabledConstraints(:)
     130             : 
     131             : #endif
     132             : 
     133             : #ifndef CPP_OLDINTEL
     134             : 
     135           0 :    WRITE(*,*) 'entered RDMFT subroutine'
     136             : 
     137             :    ! General initializations
     138           0 :    mixParam = 0.0001
     139           0 :    lagrangeMultiplier = 0.1 !results%ef
     140           0 :    spinDegenFac = 2.0 / input%jspins ! This factor is used to compensate the missing second spin in non-spinpolarized calculations
     141             : 
     142           0 :    neigTemp(:,:) = results%neig(:,:)
     143             : 
     144             :    ! Determine which states have to be considered
     145           0 :    lowestState(:,:) = -1
     146           0 :    highestState(:,:) = -1
     147           0 :    DO jspin = 1, input%jspins
     148           0 :       jsp = MERGE(1,jspin,noco%l_noco)
     149           0 :       DO ikpt = 1, kpts%nkpt
     150             :          ! determine lowest state
     151           0 :          DO iBand = 1, results%neig(ikpt,jsp)
     152           0 :             IF((results%w_iks(iBand,ikpt,jspin) / kpts%wtkpt(ikpt)).LT.(1.0-input%rdmftOccEps)) THEN
     153           0 :                lowestState(ikpt,jspin) = iBand
     154           0 :                EXIT
     155             :             END IF
     156             :          END DO
     157           0 :          lowestState(ikpt,jspin) = lowestState(ikpt,jspin) - input%rdmftStatesBelow
     158           0 :          DO iBand = lowestState(ikpt,jspin)-1, 1, -1
     159           0 :             IF((results%eig(iBand+1,ikpt,jsp) - results%eig(iBand,ikpt,jsp)).GT.degenEps) THEN
     160             :                EXIT
     161             :             END IF
     162           0 :             lowestState(ikpt,jspin) = iBand
     163             :          END DO
     164           0 :          lowestState(ikpt,jspin) = MAX(lowestState(ikpt,jspin),1)
     165             : 
     166             :          ! determine highest state
     167           0 :          DO iBand = results%neig(ikpt,jsp)-1, 1, -1
     168           0 :             IF((results%w_iks(iBand,ikpt,jspin) / kpts%wtkpt(ikpt)).GT.(0.0+input%rdmftOccEps)) THEN
     169           0 :                highestState(ikpt,jspin) = iBand
     170           0 :                EXIT
     171             :             END IF
     172             :          END DO
     173           0 :          highestState(ikpt,jspin) = highestState(ikpt,jspin) + input%rdmftStatesAbove
     174           0 :          IF((results%neig(ikpt,jsp)-1).LT.highestState(ikpt,jspin)) THEN
     175           0 :             WRITE(6,*) 'Error: Not enough states calculated:'
     176           0 :             WRITE(6,*) 'ikpt, jsp: ', ikpt, jsp
     177           0 :             WRITE(6,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
     178           0 :             WRITE(6,*) 'results%neig(ikpt,jsp): ', results%neig(ikpt,jsp)
     179           0 :             CALL juDFT_error('Not enough states calculated', calledby = 'rdmft')
     180             :          END IF
     181           0 :          DO iBand = highestState(ikpt,jspin)+1, results%neig(ikpt,jsp)
     182           0 :             IF((results%eig(iBand,ikpt,jsp) - results%eig(iBand-1,ikpt,jsp)).GT.degenEps) THEN
     183             :                EXIT
     184             :             END IF
     185           0 :             highestState(ikpt,jspin) = iBand
     186             :          END DO
     187           0 :          IF(highestState(ikpt,jspin).EQ.results%neig(ikpt,jsp)) THEN
     188           0 :             WRITE(6,*) 'Error: Highest state is degenerate:'
     189           0 :             WRITE(6,*) 'ikpt, jsp: ', ikpt, jsp
     190           0 :             WRITE(6,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
     191           0 :             WRITE(6,*) 'results%eig(highestState(ikpt,jspin),ikpt,jsp): ', results%eig(highestState(ikpt,jspin),ikpt,jsp)
     192           0 :             WRITE(6,*) 'results%eig(highestState(ikpt,jspin)-1,ikpt,jsp): ', results%eig(highestState(ikpt,jspin)-1,ikpt,jsp)
     193           0 :             CALL juDFT_error('Highest state is degenerate', calledby = 'rdmft')
     194             :          END IF
     195             :       END DO
     196             :    END DO
     197             : 
     198             :    ! Move occupations of relevant states well into allowed region
     199           0 :    numRelevantStates = SUM(highestState(:,:)) - SUM(lowestState(:,:)) + input%jspins*kpts%nkpt
     200           0 :    ALLOCATE(occupationVec(numRelevantStates))
     201           0 :    occupationVec(:) = 0.0
     202           0 :    sumOcc = 0.0
     203           0 :    tempOcc = 0.0
     204           0 :    addCharge = 0.0
     205           0 :    subCharge = 0.0
     206           0 :    addChargeWeight = 0.0
     207           0 :    subChargeWeight = 0.0
     208           0 :    iState = 0
     209           0 :    DO jspin = 1, input%jspins
     210           0 :       jsp = MERGE(1,jspin,noco%l_noco)
     211           0 :       DO ikpt = 1, kpts%nkpt
     212           0 :          DO iBand = lowestState(ikpt,jspin), highestState(ikpt,jspin)
     213           0 :             iState = iState + 1
     214           0 :             occupationVec(iState) = results%w_iks(iBand,ikpt,jsp) / (kpts%wtkpt(ikpt))
     215           0 :             sumOcc = sumOcc + results%w_iks(iBand,ikpt,jsp)
     216           0 :             IF(occupationVec(iState).LT.0.01) THEN
     217           0 :                addCharge = addCharge + (0.01-occupationVec(iState))*kpts%wtkpt(ikpt)
     218           0 :                addChargeWeight = addChargeWeight + kpts%wtkpt(ikpt)
     219             :             END IF
     220           0 :             IF(occupationVec(iState).GT.0.99) THEN
     221           0 :                subCharge = subCharge + (occupationVec(iState)-0.99)*kpts%wtkpt(ikpt)
     222           0 :                subChargeWeight = subChargeWeight + kpts%wtkpt(ikpt)
     223             :             END IF
     224             :          END DO
     225             :       END DO
     226             :    END DO
     227           0 :    iState = 0
     228           0 :    DO jspin = 1, input%jspins
     229           0 :       jsp = MERGE(1,jspin,noco%l_noco)
     230           0 :       DO ikpt = 1, kpts%nkpt
     231           0 :          DO iBand = lowestState(ikpt,jspin), highestState(ikpt,jspin)
     232           0 :             iState = iState + 1
     233           0 :             IF(occupationVec(iState).LT.0.01) THEN
     234           0 :                occupationVec(iState) = occupationVec(iState) + 0.5*(subCharge+addCharge)*(kpts%wtkpt(ikpt)/addChargeWeight)
     235             :             END IF
     236           0 :             IF(occupationVec(iState).GT.0.99) THEN
     237           0 :                occupationVec(iState) = occupationVec(iState) - 0.5*(subCharge+addCharge)*(kpts%wtkpt(ikpt)/subChargeWeight)
     238             :             END IF
     239           0 :             results%w_iks(iBand,ikpt,jsp) = occupationVec(iState) * kpts%wtkpt(ikpt)
     240           0 :             tempOcc = tempOcc + occupationVec(iState) * kpts%wtkpt(ikpt)
     241             :          END DO
     242             :       END DO
     243             :    END DO
     244           0 :    DEALLOCATE(occupationVec)
     245             : 
     246             :    ! Some more initializations
     247             : 
     248           0 :    results%neig(:,:) = highestState(:,:) + 1
     249             : 
     250           0 :    ALLOCATE(overallVCoulSSDen(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
     251           0 :    ALLOCATE(vTotSSDen(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
     252           0 :    ALLOCATE(dEdOcc(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
     253             : 
     254           0 :    CALL regCharges%init(input,atoms)
     255           0 :    CALL dos%init(input,atoms,dimension,kpts,vacuum)
     256           0 :    CALL moments%init(input,atoms)
     257           0 :    CALL overallDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
     258           0 :    CALL overallVCoul%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTCOUL)
     259           0 :    IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
     260           0 :    ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),input%jspins))
     261           0 :    DO jspin = 1, input%jspins
     262           0 :       CALL convol(stars,vTot%pw_w(:,jspin),vTot%pw(:,jspin),stars%ufft)
     263             :    END DO
     264             : 
     265           0 :    CALL vTotTemp%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTTOT)
     266           0 :    vTotTemp = vTot
     267             : 
     268           0 :    DO jsp = 1,SIZE(vTot%mt,4)
     269           0 :       DO iAtom = 1,atoms%ntype
     270           0 :          vTotTemp%mt(:atoms%jri(iAtom),0,iAtom,jsp)  = sfp_const * vTot%mt(:atoms%jri(iAtom),0,iAtom,jsp) / atoms%rmsh(:atoms%jri(iAtom),iAtom)
     271             :       END DO
     272             :    END DO
     273             : 
     274           0 :    vCoul%pw_w = CMPLX(0.0,0.0)
     275           0 :    DO jspin = 1, input%jspins
     276           0 :       CALL convol(stars,vCoul%pw_w(:,jspin),vCoul%pw(:,jspin),stars%ufft)
     277             :    END DO
     278             : 
     279           0 :    vTotSSDen = 0.0
     280             : 
     281             :    ! Calculate all single state densities
     282           0 :    CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
     283             : 
     284           0 :    numStates = 0
     285           0 :    DO jspin = 1, input%jspins
     286           0 :       jsp = MERGE(1,jspin,noco%l_noco)
     287           0 :       DO ikpt_i = 1, SIZE(mpi%k_list)
     288           0 :          ikpt= mpi%k_list(ikpt_i)
     289           0 :          DO iBand_i = 1,size(cdnvalJOB%ev_list)
     290           0 :             iband=mpi%ev_list(iband_i)
     291           0 :             IF (iband>highestState(ikpt,jsp)) CYCLE
     292           0 :             numStates = numStates + 1
     293             :             ! Construct cdnvalJob object for this state
     294             :             ! (Reasonable parallelization is not yet done - should be placed over the loops enclosing this section)
     295           0 :             cdnvalJob%k_list=[ikpt]
     296           0 :             cdnvalJob%ev_list=[iband]
     297           0 :             cdnvalJob%weights(iBand,ikpt) = spinDegenFac
     298             :      
     299             :             ! Call cdnval to construct density
     300           0 :             WRITE(*,*) 'Note: some optional flags may have to be reset in rdmft before the cdnval call'
     301           0 :             WRITE(*,*) 'This is not yet implemented!'
     302           0 :             CALL singleStateDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
     303             :             CALL cdnval(eig_id,mpi,kpts,jsp,noco,input,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
     304           0 :                         sphhar,sym,vTot,oneD,cdnvalJob,singleStateDen,regCharges,dos,results,moments)
     305             : 
     306             :             ! Store the density on disc (These are probably way too many densities to keep them in memory)
     307           0 :             filename = ''
     308           0 :             WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
     309           0 :             IF (mpi%irank.EQ.0) THEN
     310             :                CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     311           0 :                                  0,-1.0,0.0,.FALSE.,singleStateDen,TRIM(ADJUSTL(filename)))
     312             :             END IF
     313             : #ifdef CPP_MPI
     314           0 :             CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,singleStateDen)
     315             : #endif
     316             :             ! For each state calculate Integral over KS effective potential times single state density
     317           0 :             potDenInt = 0.0
     318           0 :             CALL int_nv(jsp,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vTotTemp,singleStateDen,potDenInt)
     319           0 :             vTotSSDen(iBand,ikpt,jsp) = potDenInt
     320             :          END DO
     321             :       END DO
     322             :    END DO
     323             : 
     324             :    ! Initializations for exchange contributions
     325             : 
     326           0 :    WRITE(*,*) 'RDMFT: HF initializations start'
     327             : 
     328           0 :    IF(ALLOCATED(hybrid%ne_eig)) DEALLOCATE(hybrid%ne_eig)
     329           0 :    IF(ALLOCATED(hybrid%nbands)) DEALLOCATE(hybrid%nbands)
     330           0 :    IF(ALLOCATED(hybrid%nobd)) DEALLOCATE(hybrid%nobd)
     331           0 :    IF(ALLOCATED(hybrid%nbasm)) DEALLOCATE(hybrid%nbasm)
     332           0 :    IF(ALLOCATED(hybrid%div_vv)) DEALLOCATE(hybrid%div_vv)
     333           0 :    ALLOCATE(hybrid%ne_eig(kpts%nkpt),hybrid%nbands(kpts%nkpt),hybrid%nobd(kpts%nkptf))
     334           0 :    ALLOCATE(hybrid%nbasm(kpts%nkptf))
     335           0 :    ALLOCATE(hybrid%div_vv(DIMENSION%neigd,kpts%nkpt,input%jspins))
     336           0 :    l_restart = .FALSE.
     337           0 :    l_zref = (sym%zrfs.AND.(SUM(ABS(kpts%bk(3,:kpts%nkpt))).LT.1e-9).AND..NOT.noco%l_noco)
     338           0 :    iterHF = 0
     339           0 :    hybrid%l_calhf = .TRUE.
     340             : 
     341             : !   CALL open_hybrid_io1(DIMENSION,sym%invs)
     342             : 
     343           0 :    CALL mixedbasis(atoms,kpts,dimension,input,cell,sym,xcpot,hybrid,enpara,mpi,vTot,l_restart)
     344             : 
     345           0 :    CALL open_hybrid_io2(hybrid,DIMENSION,atoms,sym%invs)
     346             : 
     347           0 :    CALL coulombmatrix(mpi,atoms,kpts,cell,sym,hybrid,xcpot,l_restart)
     348             : 
     349           0 :    CALL hf_init(hybrid,kpts,atoms,input,DIMENSION,hybdat,sym%invs)
     350             : 
     351           0 :    WRITE(*,*) 'RDMFT: HF initializations end'
     352             : 
     353           0 :    ALLOCATE(parent(kpts%nkptf))
     354           0 :    ALLOCATE(exDiag(dimension%neigd,ikpt,input%jspins))
     355           0 :    ALLOCATE(lastGradient(numStates+1))
     356           0 :    ALLOCATE(lastParameters(numStates+1))
     357           0 :    lastGradient = 0.0
     358           0 :    lastParameters = 0.0
     359           0 :    maxHistoryLength = 17!7
     360           0 :    ALLOCATE(gradientCorrections(numStates+1,maxHistoryLength))
     361           0 :    ALLOCATE(paramCorrections(numStates+1,maxHistoryLength))
     362           0 :    gradientCorrections = 0.0
     363           0 :    paramCorrections = 0.0
     364           0 :    istep = 0
     365             : 
     366             :    ! Occupation number optimization loop
     367             : 
     368           0 :    converged = .FALSE.
     369           0 :    DO WHILE (.NOT.converged)
     370             : 
     371           0 :       WRITE(*,*) 'RDMFT: convergence loop start'
     372             : 
     373           0 :       DO jspin = 1, input%jspins
     374           0 :          DO ikpt = 1,kpts%nkpt
     375           0 :             WRITE(*,*) 'jspin, ikpt: ', jspin, ikpt
     376           0 :             WRITE(*,'(8f10.5)') results%w_iks(1:10,ikpt,jspin)
     377             :          END DO
     378             :       END DO
     379             : 
     380             :       ! Calculate overall density with current occupation numbers (don't forget core electron density)
     381           0 :       CALL overallDen%resetPotDen()
     382           0 :       jspmax = input%jspins
     383           0 :       IF (noco%l_mperp) jspmax = 1
     384           0 :       DO jspin = 1,jspmax
     385           0 :          CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
     386             :          CALL cdnval(eig_id,mpi,kpts,jsp,noco,input,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
     387           0 :                      sphhar,sym,vTot,oneD,cdnvalJob,overallDen,regCharges,dos,results,moments)
     388             :       END DO
     389             : 
     390             :       CALL cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
     391           0 :                    stars,cell,sphhar,atoms,vTot,overallDen,moments,results)
     392           0 :       IF (mpi%irank.EQ.0) THEN
     393           0 :          CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,overallDen,noco%l_noco,.TRUE.,.true.,fix)
     394             :       END IF
     395             : #ifdef CPP_MPI
     396           0 :       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallDen)
     397             : #endif
     398             : 
     399             :       ! Calculate Coulomb potential for overall density (+including external potential)
     400           0 :       CALL overallDen%sum_both_spin()!workden)
     401           0 :       CALL overallVCoul%resetPotDen()
     402           0 :       ALLOCATE(overallVCoul%pw_w(size(overallVCoul%pw,1),size(overallVCoul%pw,2)))
     403           0 :       overallVCoul%pw_w(:,:) = 0.0
     404           0 :       CALL vgen_coulomb(1,mpi,DIMENSION,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,overallDen,overallVCoul)
     405           0 :       CALL convol(stars,overallVCoul%pw_w(:,1),overallVCoul%pw(:,1),stars%ufft)   ! Is there a problem with a second spin?!
     406             : #ifdef CPP_MPI
     407           0 :       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallVCoul)
     408             : #endif
     409             : 
     410           0 :       overallVCoulSSDen = 0.0
     411           0 :       DO jspin = 1, input%jspins
     412           0 :          jsp = MERGE(1,jspin,noco%l_noco)
     413           0 :          DO ikpt = 1, kpts%nkpt
     414           0 :             DO iBand = 1, highestState(ikpt,jsp)
     415             :                ! Read the single-state density from disc
     416           0 :                filename = ''
     417           0 :                WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
     418           0 :                IF (mpi%irank.EQ.0) THEN
     419             :                   CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,&
     420           0 :                                    CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,singleStateDen,TRIM(ADJUSTL(filename)))
     421           0 :                   CALL singleStateDen%sum_both_spin()!workden)
     422             :                END IF
     423             : #ifdef CPP_MPI
     424           0 :                CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,singleStateDen)
     425             : #endif
     426             : 
     427             :                ! For each state calculate integral over Coulomb potential times single state density
     428           0 :                potDenInt = 0.0
     429           0 :                CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
     430           0 :                overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
     431             :             END DO
     432             :          END DO
     433             :       END DO
     434             : 
     435             :       ! Construct exchange matrix diagonal
     436             : 
     437           0 :       exDiag = 0.0
     438           0 :       DO jspin = 1, input%jspins
     439           0 :          jsp = MERGE(1,jspin,noco%l_noco)
     440             :          ! remove weights(wtkpt) in w_iks
     441           0 :          wl_iks = 0.0
     442           0 :          DO ikptf=1,kpts%nkptf
     443           0 :             ikpt = kpts%bkp(ikptf)
     444           0 :             DO iBand=1, results%neig(ikpt,jsp)
     445           0 :                wl_iks(iBand,ikptf) = results%w_iks(iBand,ikpt,jspin) / (kpts%wtkpt(ikpt))!*kpts%nkptf) Last term to be included after applying functional
     446           0 :                wl_iks(iBand,ikptf) = sqrt(wl_iks(iBand,ikptf)) ! For Müller functional
     447           0 :                wl_iks(iBand,ikptf) = wl_iks(iBand,ikptf) / kpts%nkptf ! This is the last part of two lines above
     448             :             END DO
     449             :          END DO
     450             : 
     451           0 :          IF(ALLOCATED(eig_irr)) DEALLOCATE (eig_irr)
     452           0 :          IF(ALLOCATED(hybdat%kveclo_eig)) DEALLOCATE (hybdat%kveclo_eig)
     453           0 :          IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE (hybdat%pntgptd)
     454           0 :          IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE (hybdat%pntgpt)
     455           0 :          IF(ALLOCATED(hybdat%prodm)) DEALLOCATE (hybdat%prodm)
     456           0 :          IF(ALLOCATED(hybdat%prod)) DEALLOCATE (hybdat%prod)
     457           0 :          IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE (hybdat%nindxp1)
     458             : 
     459           0 :          results%neig(:,:) = neigTemp(:,:)
     460             : 
     461             :          CALL HF_setup(hybrid,input,sym,kpts,dimension,atoms,mpi,noco,cell,oneD,results,jspin,enpara,eig_id,&
     462           0 :                        hybdat,iterHF,sym%invs,vTot%mt(:,0,:,:),eig_irr)
     463             : 
     464           0 :          results%neig(:,:) = highestState(:,:) + 1
     465             : 
     466           0 :          mnobd = MAXVAL(hybrid%nobd)
     467             : 
     468           0 :          DO ikpt = 1,kpts%nkpt
     469             : 
     470           0 :             CALL lapw%init(input,noco,kpts,atoms,sym,ikpt,cell,l_zref)
     471             : 
     472           0 :             parent = 0
     473           0 :             CALL symm_hf_init(sym,kpts,ikpt,nsymop,rrot,psym)
     474             :             CALL symm_hf(kpts,ikpt,sym,dimension,hybdat,eig_irr,atoms,hybrid,cell,lapw,jspin,mpi,&
     475           0 :                          rrot,nsymop,psym,nkpt_EIBZ,n_q,parent,pointer_EIBZ,nsest,indx_sest)
     476             : 
     477           0 :             exMat%l_real=sym%invs
     478             :             CALL exchange_valence_hf(ikpt,kpts,nkpt_EIBZ, sym,atoms,hybrid,cell,dimension,input,jspin,hybdat,mnobd,lapw,&
     479             :                                      eig_irr,results,parent,pointer_EIBZ,n_q,wl_iks,iterHF,xcpot,noco,nsest,indx_sest,&
     480           0 :                                      mpi,exMat)
     481           0 :             CALL exchange_vccv1(ikpt,atoms,hybrid,hybdat,dimension,jspin,lapw,nsymop,nsest,indx_sest,mpi,1.0,results,exMat)
     482             : 
     483             :             !Start of workaround for increased functionality of symmetrizeh (call it))
     484             : 
     485           0 :             nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
     486             : 
     487           0 :             CALL olap%alloc(sym%invs,nbasfcn)
     488           0 :             CALL read_olap(olap, kpts%nkpt*(jspin-1)+ikpt)
     489           0 :             IF (olap%l_real) THEN
     490           0 :                DO i = 1, nbasfcn
     491           0 :                   DO j = 1, i
     492           0 :                      olap%data_r(i,j) = olap%data_r(j,i)
     493             :                   END DO
     494             :                END DO
     495             :             ELSE
     496           0 :                DO i = 1, nbasfcn
     497           0 :                   DO j = 1, i
     498           0 :                      olap%data_c(i,j) = CONJG(olap%data_c(j,i))
     499             :                   END DO
     500             :                END DO
     501           0 :                olap%data_c = conjg(olap%data_c)
     502             :             END IF
     503             : 
     504           0 :             CALL zMat%init(olap%l_real,nbasfcn,dimension%neigd)
     505             : 
     506           0 :             CALL read_eig(eig_id,ikpt,jspin,list=[(i,i=1,hybrid%nbands(ikpt))],neig=nbands,zmat=zMat)
     507             : 
     508             : !            CALL read_z(zMat,kpts%nkpt*(jspin-1)+ikpt)
     509           0 :             zMat%matsize2 = hybrid%nbands(ikpt) ! reduce "visible matsize" for the following computations
     510             : 
     511           0 :             CALL olap%multiply(zMat,trafo)
     512             : 
     513           0 :             CALL invtrafo%alloc(olap%l_real,hybrid%nbands(ikpt),nbasfcn)
     514           0 :             CALL trafo%TRANSPOSE(invtrafo)
     515           0 :             IF(.NOT.invtrafo%l_real) invtrafo%data_c = CONJG(invtrafo%data_c)
     516             : 
     517           0 :             DO i = 1, hybrid%nbands(ikpt)
     518           0 :                DO j = 1, i-1
     519           0 :                   IF (exMat%l_real) THEN
     520           0 :                      exMat%data_r(i,j)=exMat%data_r(j,i)
     521             :                   ELSE
     522           0 :                      exMat%data_c(i,j)=conjg(exMat%data_c(j,i))
     523             :                   END IF
     524             :                END DO
     525             :             END DO
     526             : 
     527           0 :             CALL exMat%multiply(invtrafo,tmpMat)
     528           0 :             CALL trafo%multiply(tmpMat,exMatLAPW)
     529             : 
     530           0 :             CALL symmetrizeh(atoms,kpts%bkf(:,ikpt),dimension,jspin,lapw,sym,hybdat%kveclo_eig,cell,nsymop,psym,exMatLAPW)
     531             : 
     532           0 :             IF (.NOT.exMatLAPW%l_real) exMatLAPW%data_c=conjg(exMatLAPW%data_c) 
     533           0 :             zMat%matsize1=MIN(zMat%matsize1,exMatLAPW%matsize2)
     534             : 
     535           0 :             CALL exMatLAPW%multiply(zMat,tmpMat)
     536             : 
     537           0 :             DO iBand = 1, highestState(ikpt,jsp)
     538           0 :                IF (zMat%l_real) THEN
     539           0 :                   exDiag(iBand,ikpt,jspin) = dot_product(zMat%data_r(:zMat%matsize1,iband),tmpMat%data_r(:,iband))
     540             :                ELSE
     541           0 :                   exDiag(iBand,ikpt,jspin) = REAL(dot_product(zMat%data_c(:zMat%matsize1,iband),tmpMat%data_c(:,iband)))
     542             :                END IF
     543             :             END DO
     544             : 
     545             :             !End of workaround for increased functionality of symmetrizeh (call it))
     546             : 
     547             : !            DO iBand = 1, highestState(ikpt,jsp)
     548             : !               IF (exMat%l_real) THEN
     549             : !                  exDiag(iBand,ikpt,jspin) = exMat%data_r(iBand,iBand)
     550             : !               ELSE
     551             : !                  exDiag(iBand,ikpt,jspin) = REAL(exMat%data_c(iBand,iBand))
     552             : !               END IF
     553             : !            END DO
     554             :          END DO
     555             :       END DO
     556             : 
     557             :       ! Calculate total energy derivative with respect to occupations (dEdOcc)
     558             : 
     559           0 :       gradSum = 0.0
     560           0 :       DO ispin = 1, input%jspins
     561           0 :          isp = MERGE(1,ispin,noco%l_noco)
     562             : !         CALL cdnvalJob%init(mpi,input,kpts,noco,results,isp,banddos=banddos)
     563           0 :          DO ikpt = 1, kpts%nkpt
     564           0 :             DO iBand = 1, highestState(ikpt,isp)
     565           0 :                occStateI = results%w_iks(iBand,ikpt,isp) / (kpts%wtkpt(ikpt))!*kpts%nkptf)
     566           0 :                occStateI = MAX(occStateI,minOcc)
     567             : !               IF(occStateI.LT.1.0e-7) occStateI = 5.0e-4 ! This is preliminary. I have to discuss what do do here.
     568             : !               occStateI = cdnvalJob%weights(iBand,ikpt)
     569           0 :                rdmftFunctionalValue = 0.5*0.5*SQRT(1.0/occStateI) ! for Müller functional derivative
     570             : 
     571           0 :                exchangeTerm = - rdmftFunctionalValue * exDiag(iBand,ikpt,isp) * kpts%wtkpt(ikpt) * spinDegenFac !*kpts%nkptf
     572             : 
     573             :                dEdOcc(iBand,ikpt,isp) = +((spinDegenFac * results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
     574           0 :                                               overallVCoulSSDen(iBand,ikpt,isp) + exchangeTerm)
     575             : 
     576           0 :                WRITE(*,*) 'ENERGY GRADIENT CONTRIBUTIONS'
     577           0 :                WRITE(*,*) 'ispin, ikpt, iBand', ispin, ikpt, iBand
     578           0 :                WRITE(*,*) 'results%eig(iBand,ikpt,isp)', results%eig(iBand,ikpt,isp)
     579           0 :                WRITE(*,*) 'vTotSSDen(iBand,ikpt,isp)', vTotSSDen(iBand,ikpt,isp)
     580           0 :                WRITE(*,*) 'overallVCoulSSDen(iBand,ikpt,isp)', overallVCoulSSDen(iBand,ikpt,isp)
     581           0 :                WRITE(*,*) 'exchangeTerm', exchangeTerm
     582           0 :                WRITE(*,*) 'exDiag(iBand,ikpt,isp)', exDiag(iBand,ikpt,isp)
     583           0 :                WRITE(*,*) 'rdmftFunctionalValue', rdmftFunctionalValue
     584             : 
     585             : 
     586           0 :                gradSum = gradSum + dEdOcc(iBand,ikpt,isp) ! * results%w_iks(iBand,ikpt,isp)
     587             :             END DO
     588             :          END DO
     589             :       END DO
     590           0 :       lagrangeMultiplier = -gradSum / numStates !(input%zelec/(2.0/REAL(input%jspins)))
     591             : 
     592           0 :    WRITE(*,*) 'lagrangeMultiplier: ', lagrangeMultiplier
     593             : 
     594             :       ! Optimize occupation numbers
     595             : 
     596           0 :       ALLOCATE (gradient(numStates+1))
     597           0 :       ALLOCATE (parameters(numStates+1))
     598           0 :       ALLOCATE (equalityLinCombi(numStates+1))
     599           0 :       ALLOCATE (minConstraints(numStates+1))
     600           0 :       ALLOCATE (maxConstraints(numStates+1))
     601           0 :       ALLOCATE (enabledConstraints(numStates+1))
     602           0 :       enabledConstraints = .FALSE.
     603           0 :       enabledConstraints(1:numStates) = .TRUE.
     604           0 :       minConstraints(:) = 0.0
     605           0 :       maxConstraints(:) = 1.0
     606           0 :       equalityLinCombi = 0.0
     607           0 :       gradient = 0.0
     608           0 :       parameters = 0.0
     609           0 :       iState = 0
     610           0 :       DO ispin = 1, input%jspins
     611           0 :          isp = MERGE(1,ispin,noco%l_noco)
     612           0 :          DO ikpt = 1, kpts%nkpt
     613           0 :             DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
     614           0 :                iState = iState + 1
     615           0 :                occStateI = results%w_iks(iBand,ikpt,isp) / kpts%wtkpt(ikpt)
     616           0 :                occStateI = MAX(occStateI,minOcc)
     617           0 :                equalityLinCombi(iState) = kpts%wtkpt(ikpt)
     618           0 :                gradient(iState) = dEdOcc(iBand,ikpt,isp) + lagrangeMultiplier
     619           0 :                gradient(numStates+1) = gradient(numStates+1) + occStateI * kpts%wtkpt(ikpt)
     620           0 :                parameters(iState) = occStateI
     621             :             END DO
     622             :          END DO
     623             :       END DO
     624           0 :       equalityCriterion = input%zelec/(2.0/REAL(input%jspins))
     625           0 :       gradient(numStates+1) = gradient(numStates+1) - equalityCriterion ! This should actually always be 0.0
     626           0 :       parameters(numStates+1) = lagrangeMultiplier
     627             : 
     628           0 :       mixParam = 0.01 / MAXVAL(ABS(gradient(:numStates)))
     629           0 :       WRITE(*,*) 'mixParam: ', mixParam
     630             : 
     631             :       CALL bfgs_b2(numStates+1,gradient,lastGradient,minConstraints,maxConstraints,enabledConstraints,parameters,&
     632             :                    lastParameters,equalityLinCombi,equalityCriterion,maxHistoryLength,paramCorrections,&
     633           0 :                    gradientCorrections,iStep,mixParam,converged,convCrit)
     634             : 
     635           0 :       iState = 0
     636           0 :       DO ispin = 1, input%jspins
     637           0 :          isp = MERGE(1,ispin,noco%l_noco)
     638           0 :          DO ikpt = 1, kpts%nkpt
     639           0 :             DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
     640           0 :                iState = iState + 1
     641           0 :                results%w_iks(iBand,ikpt,isp) = MERGE(parameters(iState) * kpts%wtkpt(ikpt),0.0,parameters(iState).GT.minOcc)
     642             :             END DO
     643             :          END DO
     644             :       END DO
     645             : 
     646           0 :       DEALLOCATE (enabledConstraints,maxConstraints,minConstraints)
     647           0 :       DEALLOCATE (parameters,gradient,equalityLinCombi)
     648             : 
     649             : 
     650             :    END DO ! WHILE (.NOT.converged)
     651             : 
     652           0 :    WRITE(*,*) 'RDMFT: convergence loop end'
     653             : 
     654           0 :    hybrid%l_calhf = .FALSE.
     655             : 
     656             :    ! Calculate final overall density
     657             : 
     658             :    !I think we need most of cdngen at this place so I just use cdngen
     659           0 :    CALL outDen%resetPotDen()
     660             :    CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,DIMENSION,kpts,atoms,sphhar,stars,sym,&
     661           0 :                enpara,cell,noco,vTot,results,oneD,coreSpecInput,archiveType,xcpot,outDen, EnergyDen)
     662             : 
     663             :    ! Calculate RDMFT energy
     664           0 :    rdmftEnergy = 0.0
     665           0 :    DO ispin = 1, input%jspins
     666           0 :       isp = MERGE(1,ispin,noco%l_noco)
     667             : !      CALL cdnvalJob%init(mpi,input,kpts,noco,results,isp,banddos=banddos)
     668           0 :       DO ikpt = 1, kpts%nkpt
     669           0 :          DO iBand = 1, highestState(ikpt,isp)
     670           0 :             occStateI = results%w_iks(iBand,ikpt,isp) / (kpts%wtkpt(ikpt))!*kpts%nkptf)
     671           0 :             rdmftFunctionalValue = 0.5*SQRT(occStateI) ! for Müller functional
     672             : 
     673           0 :             exchangeTerm = -rdmftFunctionalValue * exDiag(iBand,ikpt,isp) * kpts%wtkpt(ikpt) * spinDegenFac !*kpts%nkptf
     674             : 
     675             :             rdmftEnergy = rdmftEnergy + exchangeTerm + &
     676             :                           occStateI * ((spinDegenFac*results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
     677           0 :                                        overallVCoulSSDen(iBand,ikpt,isp))
     678             :          END DO
     679             :       END DO
     680             :    END DO
     681             : 
     682           0 :    results%neig(:,:) = neigTemp(:,:)
     683             : 
     684           0 :    WRITE(6,'(a,f20.10,a)') 'RDMFT energy: ', rdmftEnergy, ' Htr'
     685             : 
     686             : #endif
     687           0 : END SUBROUTINE rdmft
     688             : 
     689             : END MODULE m_rdmft

Generated by: LCOV version 1.13