1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate MAO's and analyze wavefunctions 8!> \par History 9!> 03.2016 created [JGH] 10!> 12.2016 split into four modules [JGH] 11!> \author JGH 12! ************************************************************************************************** 13MODULE mao_optimizer 14 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& 15 dbcsr_deallocate_matrix_set 16 USE dbcsr_api, ONLY: & 17 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_info, dbcsr_norm, & 18 dbcsr_norm_maxabsnorm, dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_type, & 19 dbcsr_type_symmetric 20 USE kinds, ONLY: dp 21 USE mao_methods, ONLY: mao_function,& 22 mao_function_gradient,& 23 mao_initialization,& 24 mao_orthogonalization,& 25 mao_scalar_product 26#include "./base/base_uses.f90" 27 28 IMPLICIT NONE 29 PRIVATE 30 31 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_optimizer' 32 33 PUBLIC :: mao_optimize 34 35! ************************************************************************************************** 36 37CONTAINS 38 39! ************************************************************************************************** 40!> \brief ... 41!> \param mao_coef ... 42!> \param matrix_q ... 43!> \param matrix_smm ... 44!> \param electra ... 45!> \param max_iter ... 46!> \param eps_grad ... 47!> \param eps1 ... 48!> \param iw ... 49! ************************************************************************************************** 50 SUBROUTINE mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, iw) 51 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, matrix_q, matrix_smm 52 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: electra 53 INTEGER, INTENT(IN) :: max_iter 54 REAL(KIND=dp), INTENT(IN) :: eps_grad, eps1 55 INTEGER, INTENT(IN) :: iw 56 57 CHARACTER(len=*), PARAMETER :: routineN = 'mao_optimize', routineP = moduleN//':'//routineN 58 59 INTEGER :: handle, i, ispin, iter, nspin 60 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_sizes_a, mao_sizes_b 61 LOGICAL :: spin_warning 62 REAL(KIND=dp) :: a1, a2, alpha, an, beta, eps_fun, fa1, & 63 fa2, fnnew, fnold, fval, grad_norm 64 TYPE(dbcsr_distribution_type) :: dbcsr_dist 65 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_grad 66 TYPE(dbcsr_type) :: amat, binv, cgmat 67 68 CALL timeset(routineN, handle) 69 70 eps_fun = 10._dp*eps_grad 71 72 nspin = SIZE(mao_coef, 1) 73 74 IF (iw > 0) THEN 75 WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!' 76 WRITE (UNIT=iw, FMT="(T32,A)") "MAO BASIS GENERATION" 77 WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!' 78 END IF 79 80 ! initialize MAO coeficients from diagonal blocks of the Q matrix 81 DO ispin = 1, nspin 82 CALL mao_initialization(mao_coef(ispin)%matrix, & 83 matrix_q(ispin)%matrix, matrix_smm(1)%matrix, eps1) 84 END DO 85 ! check for mao sizes 86 spin_warning = .FALSE. 87 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=mao_sizes_a, distribution=dbcsr_dist) 88 IF (nspin == 2) THEN 89 CALL dbcsr_get_info(mao_coef(2)%matrix, col_blk_size=mao_sizes_b, distribution=dbcsr_dist) 90 DO i = 1, SIZE(mao_sizes_a) 91 IF (mao_sizes_a(i) /= mao_sizes_b(i)) spin_warning = .TRUE. 92 END DO 93 END IF 94 IF (iw > 0 .AND. spin_warning) THEN 95 WRITE (UNIT=iw, FMT="(/,A)") ">>>> WARNING: Different number of MAOs for different spins" 96 END IF 97 IF (iw > 0) THEN 98 DO ispin = 1, nspin 99 IF (ispin == 1) THEN 100 WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_a) 101 WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_a(:) 102 ELSE IF (ispin == 2) THEN 103 WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_b) 104 WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_b(:) 105 END IF 106 END DO 107 WRITE (UNIT=iw, FMT="(A)") 108 END IF 109 110 IF (max_iter < 1) THEN 111 ! projection only 112 DO ispin = 1, nspin 113 CALL dbcsr_get_info(mao_coef(ispin)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist) 114 CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, & 115 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0) 116 CALL mao_function(mao_coef(ispin)%matrix, fval, matrix_q(ispin)%matrix, & 117 matrix_smm(1)%matrix, binv, .FALSE.) 118 IF (iw > 0) THEN 119 WRITE (UNIT=iw, FMT="(T2,A,T31,A,I2,T56,A,F12.8)") & 120 "MAO Projection", "Spin:", ispin, "Completeness:", fval/electra(ispin) 121 END IF 122 CALL dbcsr_release(binv) 123 END DO 124 ELSE 125 ! optimize MAOs 126 NULLIFY (mao_grad) 127 CALL dbcsr_allocate_matrix_set(mao_grad, nspin) 128 DO ispin = 1, nspin 129 ALLOCATE (mao_grad(ispin)%matrix) 130 CALL dbcsr_create(matrix=mao_grad(ispin)%matrix, name="MAO_GRAD", template=mao_coef(1)%matrix) 131 CALL dbcsr_reserve_diag_blocks(matrix=mao_grad(ispin)%matrix) 132 END DO 133 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist) 134 CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, & 135 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0) 136 CALL dbcsr_create(cgmat, template=mao_grad(1)%matrix) 137 CALL dbcsr_create(amat, template=mao_coef(1)%matrix) 138 DO ispin = 1, nspin 139 alpha = 0.25_dp 140 beta = 0.0_dp 141 CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, & 142 matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .FALSE.) 143 CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix) 144 CALL dbcsr_norm(mao_grad(ispin)%matrix, dbcsr_norm_maxabsnorm, norm_scalar=grad_norm) 145 fnold = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix) 146 IF (iw > 0) THEN 147 WRITE (UNIT=iw, FMT="(/,T2,A,T73,A,I2)") "MAO OPTIMIZATION", "Spin =", ispin 148 WRITE (UNIT=iw, FMT="(T2,A,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") & 149 "Initialization", "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha 150 END IF 151 DO iter = 1, max_iter 152 IF (grad_norm < eps_grad) EXIT 153 IF ((1.0_dp - fval/electra(ispin)) < eps_fun) EXIT 154 CALL dbcsr_add(mao_coef(ispin)%matrix, cgmat, 1.0_dp, alpha) 155 CALL mao_orthogonalization(mao_coef(ispin)%matrix, matrix_smm(1)%matrix) 156 CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, & 157 matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.) 158 CALL dbcsr_norm(mao_grad(ispin)%matrix, dbcsr_norm_maxabsnorm, norm_scalar=grad_norm) 159 IF (iw > 0) THEN 160 WRITE (UNIT=iw, FMT="(T2,A,i8,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") & 161 "iter=", iter, "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha 162 END IF 163 fnnew = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix) 164 beta = fnnew/fnold 165 CALL dbcsr_add(cgmat, mao_grad(ispin)%matrix, beta, 1.0_dp) 166 fnold = fnnew 167 ! line search, update alpha 168 CALL dbcsr_copy(amat, mao_coef(ispin)%matrix) 169 CALL dbcsr_add(amat, cgmat, 1.0_dp, 0.5_dp*alpha) 170 CALL mao_orthogonalization(amat, matrix_smm(1)%matrix) 171 CALL mao_function(amat, fa1, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.) 172 CALL dbcsr_copy(amat, mao_coef(ispin)%matrix) 173 CALL dbcsr_add(amat, cgmat, 1.0_dp, alpha) 174 CALL mao_orthogonalization(amat, matrix_smm(1)%matrix) 175 CALL mao_function(amat, fa2, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.) 176 a2 = (4._dp*fa1 - fa2 - 3._dp*fval)/alpha 177 a1 = (fa2 - fval - a2*alpha)/(alpha*alpha) 178 IF (ABS(a1) > 1.e-14_dp) THEN 179 an = -a2/(2._dp*a1) 180 an = MIN(an, 2.0_dp*alpha) 181 ELSE 182 an = 2.0_dp*alpha 183 END IF 184 IF (an < 0.05_dp .OR. a1 > 0.0_dp) THEN 185 CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix) 186 alpha = 0.25_dp 187 ELSE 188 alpha = an 189 END IF 190 END DO 191 END DO 192 CALL dbcsr_release(binv) 193 CALL dbcsr_release(cgmat) 194 CALL dbcsr_release(amat) 195 IF (ASSOCIATED(mao_grad)) CALL dbcsr_deallocate_matrix_set(mao_grad) 196 END IF 197 198 CALL timestop(handle) 199 200 END SUBROUTINE mao_optimize 201 202END MODULE mao_optimizer 203