Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : MODULE qs_charge_mixing
      10             : 
      11             :    USE kinds,                           ONLY: dp
      12             :    USE mathlib,                         ONLY: get_pseudo_inverse_svd
      13             :    USE message_passing,                 ONLY: mp_para_env_type
      14             :    USE qs_density_mixing_types,         ONLY: broyden_mixing_nr,&
      15             :                                               gspace_mixing_nr,&
      16             :                                               mixing_storage_type,&
      17             :                                               multisecant_mixing_nr,&
      18             :                                               pulay_mixing_nr
      19             : #include "./base/base_uses.f90"
      20             : 
      21             :    IMPLICIT NONE
      22             : 
      23             :    PRIVATE
      24             : 
      25             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_charge_mixing'
      26             : 
      27             :    PUBLIC :: charge_mixing
      28             : 
      29             : CONTAINS
      30             : 
      31             : ! **************************************************************************************************
      32             : !> \brief  Driver for the charge mixing, calls the proper routine given the requested method
      33             : !> \param mixing_method ...
      34             : !> \param mixing_store ...
      35             : !> \param charges ...
      36             : !> \param para_env ...
      37             : !> \param iter_count ...
      38             : !> \par History
      39             : !> \author JGH
      40             : ! **************************************************************************************************
      41       20960 :    SUBROUTINE charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count)
      42             :       INTEGER, INTENT(IN)                                :: mixing_method
      43             :       TYPE(mixing_storage_type), POINTER                 :: mixing_store
      44             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: charges
      45             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      46             :       INTEGER, INTENT(IN)                                :: iter_count
      47             : 
      48             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'charge_mixing'
      49             : 
      50             :       INTEGER                                            :: handle, ia, ii, imin, inow, nbuffer, ns, &
      51             :                                                             nvec
      52             :       REAL(dp)                                           :: alpha
      53             : 
      54       20960 :       CALL timeset(routineN, handle)
      55             : 
      56       20960 :       IF (mixing_method >= gspace_mixing_nr) THEN
      57         556 :          CPASSERT(ASSOCIATED(mixing_store))
      58         556 :          mixing_store%ncall = mixing_store%ncall + 1
      59         556 :          ns = SIZE(charges, 2)
      60         556 :          ns = MIN(ns, mixing_store%max_shell)
      61         556 :          alpha = mixing_store%alpha
      62         556 :          nbuffer = mixing_store%nbuffer
      63         556 :          inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
      64         556 :          imin = inow - 1
      65         556 :          IF (imin == 0) imin = nbuffer
      66         556 :          IF (mixing_store%ncall > nbuffer) THEN
      67         334 :             nvec = nbuffer
      68             :          ELSE
      69         222 :             nvec = mixing_store%ncall - 1
      70             :          END IF
      71         556 :          IF (mixing_store%ncall > 1) THEN
      72             :             ! store in/out charge difference
      73        2322 :             DO ia = 1, mixing_store%nat_local
      74        1802 :                ii = mixing_store%atlist(ia)
      75       11332 :                mixing_store%dacharge(ia, 1:ns, imin) = mixing_store%acharge(ia, 1:ns, imin) - charges(ii, 1:ns)
      76             :             END DO
      77             :          END IF
      78         556 :          IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
      79             :             ! skip mixing
      80          36 :             mixing_store%iter_method = "NoMix"
      81         520 :          ELSEIF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
      82          36 :             CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
      83          36 :             mixing_store%iter_method = "Mixing"
      84         484 :          ELSEIF (mixing_method == gspace_mixing_nr) THEN
      85           0 :             CPABORT("Kerker method not available for Charge Mixing")
      86         484 :          ELSEIF (mixing_method == pulay_mixing_nr) THEN
      87           0 :             CPABORT("Pulay method not available for Charge Mixing")
      88         484 :          ELSEIF (mixing_method == broyden_mixing_nr) THEN
      89         484 :             CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
      90         484 :             mixing_store%iter_method = "Broy."
      91           0 :          ELSEIF (mixing_method == multisecant_mixing_nr) THEN
      92           0 :             CPABORT("Multisecant_mixing method not available for Charge Mixing")
      93             :          END IF
      94             : 
      95             :          ! store new 'input' charges
      96        2452 :          DO ia = 1, mixing_store%nat_local
      97        1896 :             ii = mixing_store%atlist(ia)
      98       11932 :             mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
      99             :          END DO
     100             : 
     101             :       END IF
     102             : 
     103       20960 :       CALL timestop(handle)
     104             : 
     105       20960 :    END SUBROUTINE charge_mixing
     106             : 
     107             : ! **************************************************************************************************
     108             : !> \brief Simple charge mixing
     109             : !> \param mixing_store ...
     110             : !> \param charges ...
     111             : !> \param alpha ...
     112             : !> \param imin ...
     113             : !> \param ns ...
     114             : !> \param para_env ...
     115             : !> \author JGH
     116             : ! **************************************************************************************************
     117          36 :    SUBROUTINE mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
     118             :       TYPE(mixing_storage_type), POINTER                 :: mixing_store
     119             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: charges
     120             :       REAL(KIND=dp), INTENT(IN)                          :: alpha
     121             :       INTEGER, INTENT(IN)                                :: imin, ns
     122             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     123             : 
     124             :       INTEGER                                            :: ia, ii
     125             : 
     126        1156 :       charges = 0.0_dp
     127             : 
     128         130 :       DO ia = 1, mixing_store%nat_local
     129          94 :          ii = mixing_store%atlist(ia)
     130         600 :          charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
     131             :       END DO
     132             : 
     133        2276 :       CALL para_env%sum(charges)
     134             : 
     135          36 :    END SUBROUTINE mix_charges_only
     136             : 
     137             : ! **************************************************************************************************
     138             : !> \brief Broyden charge mixing
     139             : !> \param mixing_store ...
     140             : !> \param charges ...
     141             : !> \param inow ...
     142             : !> \param nvec ...
     143             : !> \param ns ...
     144             : !> \param para_env ...
     145             : !> \author JGH
     146             : ! **************************************************************************************************
     147         484 :    SUBROUTINE broyden_mixing(mixing_store, charges, inow, nvec, ns, para_env)
     148             :       TYPE(mixing_storage_type), POINTER                 :: mixing_store
     149             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: charges
     150             :       INTEGER, INTENT(IN)                                :: inow, nvec, ns
     151             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     152             : 
     153             :       INTEGER                                            :: i, ia, ii, imin, j, nbuffer, nv
     154             :       REAL(KIND=dp)                                      :: alpha, maxw, minw, omega0, rskip, wdf, &
     155             :                                                             wfac
     156         484 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cvec, gammab
     157         484 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: amat, beta
     158         484 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: dq_last, dq_now, q_last, q_now
     159             : 
     160           0 :       CPASSERT(nvec > 1)
     161             : 
     162         484 :       nbuffer = mixing_store%nbuffer
     163         484 :       alpha = mixing_store%alpha
     164         484 :       imin = inow - 1
     165         484 :       IF (imin == 0) imin = nvec
     166         484 :       nv = nvec - 1
     167             : 
     168             :       ! charge vectors
     169         484 :       q_now => mixing_store%acharge(:, :, inow)
     170         484 :       q_last => mixing_store%acharge(:, :, imin)
     171         484 :       dq_now => mixing_store%dacharge(:, :, inow)
     172         484 :       dq_last => mixing_store%dacharge(:, :, imin)
     173             : 
     174         484 :       IF (nvec == nbuffer) THEN
     175             :          ! reshuffel Broyden storage n->n-1
     176        1194 :          DO i = 1, nv - 1
     177         860 :             mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
     178       21080 :             mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
     179       21414 :             mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
     180             :          END DO
     181        1194 :          DO i = 1, nv - 1
     182        4450 :             DO j = 1, nv - 1
     183        4116 :                mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
     184             :             END DO
     185             :          END DO
     186             :       END IF
     187             : 
     188         484 :       omega0 = 0.01_dp
     189         484 :       minw = 1.0_dp
     190         484 :       maxw = 100000.0_dp
     191         484 :       wfac = 0.01_dp
     192             : 
     193       11444 :       mixing_store%wbroy(nv) = SUM(dq_now(:, :)**2)
     194         484 :       CALL para_env%sum(mixing_store%wbroy(nv))
     195         484 :       mixing_store%wbroy(nv) = SQRT(mixing_store%wbroy(nv))
     196         484 :       IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN
     197         472 :          mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv)
     198             :       ELSE
     199          12 :          mixing_store%wbroy(nv) = maxw
     200             :       END IF
     201         484 :       IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw
     202             : 
     203             :       ! dfbroy
     204       22404 :       mixing_store%dfbroy(:, :, nv) = dq_now(:, :) - dq_last(:, :)
     205       11444 :       wdf = SUM(mixing_store%dfbroy(:, :, nv)**2)
     206         484 :       CALL para_env%sum(wdf)
     207         484 :       wdf = 1.0_dp/SQRT(wdf)
     208       11444 :       mixing_store%dfbroy(:, :, nv) = wdf*mixing_store%dfbroy(:, :, nv)
     209             : 
     210             :       ! abroy matrix
     211        2350 :       DO i = 1, nv
     212       45846 :          wfac = SUM(mixing_store%dfbroy(:, :, i)*mixing_store%dfbroy(:, :, nv))
     213        1866 :          CALL para_env%sum(wfac)
     214        1866 :          mixing_store%abroy(i, nv) = wfac
     215        2350 :          mixing_store%abroy(nv, i) = wfac
     216             :       END DO
     217             : 
     218             :       ! broyden matrices
     219        4356 :       ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
     220        2350 :       DO i = 1, nv
     221       45846 :          wfac = SUM(mixing_store%dfbroy(:, :, i)*dq_now(:, :))
     222        1866 :          CALL para_env%sum(wfac)
     223        2350 :          cvec(i) = mixing_store%wbroy(i)*wfac
     224             :       END DO
     225             : 
     226        2350 :       DO i = 1, nv
     227       12508 :          DO j = 1, nv
     228       12508 :             beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
     229             :          END DO
     230        2350 :          beta(i, i) = beta(i, i) + omega0*omega0
     231             :       END DO
     232             : 
     233         484 :       rskip = 1.e-12_dp
     234         484 :       CALL get_pseudo_inverse_svd(beta, amat, rskip)
     235       14858 :       gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
     236             : 
     237             :       ! build ubroy
     238       22404 :       mixing_store%ubroy(:, :, nv) = alpha*mixing_store%dfbroy(:, :, nv) + wdf*(q_now(:, :) - q_last(:, :))
     239             : 
     240       19984 :       charges = 0.0_dp
     241        2192 :       DO ia = 1, mixing_store%nat_local
     242        1708 :          ii = mixing_store%atlist(ia)
     243       10732 :          charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
     244             :       END DO
     245        2350 :       DO i = 1, nv
     246        9280 :          DO ia = 1, mixing_store%nat_local
     247        6930 :             ii = mixing_store%atlist(ia)
     248       43446 :             charges(ii, 1:ns) = charges(ii, 1:ns) - mixing_store%wbroy(i)*gammab(i)*mixing_store%ubroy(ia, 1:ns, i)
     249             :          END DO
     250             :       END DO
     251       39484 :       CALL para_env%sum(charges)
     252             : 
     253         484 :       DEALLOCATE (amat, beta, cvec, gammab)
     254             : 
     255         484 :    END SUBROUTINE broyden_mixing
     256             : 
     257             : END MODULE qs_charge_mixing