1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Set of routines handling the localization for molecular properties 8! ************************************************************************************************** 9MODULE molecular_moments 10 USE atomic_kind_types, ONLY: atomic_kind_type,& 11 get_atomic_kind 12 USE cell_types, ONLY: cell_type,& 13 pbc 14 USE cp_control_types, ONLY: dft_control_type 15 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply 16 USE cp_fm_basic_linalg, ONLY: cp_fm_schur_product 17 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 18 cp_fm_struct_release,& 19 cp_fm_struct_type 20 USE cp_fm_types, ONLY: cp_fm_create,& 21 cp_fm_get_info,& 22 cp_fm_p_type,& 23 cp_fm_release,& 24 cp_fm_to_fm,& 25 cp_fm_type 26 USE cp_log_handling, ONLY: cp_get_default_logger,& 27 cp_logger_type 28 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 29 cp_print_key_unit_nr 30 USE cp_para_types, ONLY: cp_para_env_type 31 USE dbcsr_api, ONLY: dbcsr_copy,& 32 dbcsr_deallocate_matrix,& 33 dbcsr_p_type,& 34 dbcsr_set 35 USE distribution_1d_types, ONLY: distribution_1d_type 36 USE input_constants, ONLY: use_mom_ref_com 37 USE input_section_types, ONLY: section_vals_type,& 38 section_vals_val_get 39 USE kinds, ONLY: dp 40 USE message_passing, ONLY: mp_bcast,& 41 mp_maxloc,& 42 mp_sum 43 USE molecule_kind_types, ONLY: get_molecule_kind,& 44 molecule_kind_type 45 USE molecule_types, ONLY: molecule_type 46 USE moments_utils, ONLY: get_reference_point 47 USE orbital_pointers, ONLY: indco,& 48 ncoset 49 USE particle_types, ONLY: particle_type 50 USE qs_environment_types, ONLY: get_qs_env,& 51 qs_environment_type 52 USE qs_kind_types, ONLY: get_qs_kind,& 53 qs_kind_type 54 USE qs_loc_types, ONLY: qs_loc_env_new_type 55 USE qs_moments, ONLY: build_local_moment_matrix 56#include "./base/base_uses.f90" 57 58 IMPLICIT NONE 59 60 PRIVATE 61 62 ! *** Public *** 63 PUBLIC :: calculate_molecular_moments 64 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecular_moments' 66 67CONTAINS 68 69! ************************************************************************************************** 70!> \brief Calculates electrical molecular moments using local operators r-r_ref 71!> r_ref: center of mass of the molecule 72!> Output is in atomic units 73!> \param qs_env the qs_env in which the qs_env lives 74!> \param qs_loc_env ... 75!> \param mo_local ... 76!> \param loc_print_key ... 77!> \param molecule_set ... 78! ************************************************************************************************** 79 SUBROUTINE calculate_molecular_moments(qs_env, qs_loc_env, mo_local, loc_print_key, molecule_set) 80 TYPE(qs_environment_type), POINTER :: qs_env 81 TYPE(qs_loc_env_new_type), INTENT(IN) :: qs_loc_env 82 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_local 83 TYPE(section_vals_type), POINTER :: loc_print_key 84 TYPE(molecule_type), POINTER :: molecule_set(:) 85 86 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_molecular_moments', & 87 routineP = moduleN//':'//routineN 88 89 INTEGER :: akind, first_atom, i, iatom, ikind, imol, imol_now, iounit, iproc, ispin, j, lx, & 90 ly, lz, molkind, n, n1, n2, natom, ncol_global, nm, nmol, nmols, norder, nrow_global, ns, & 91 nspins 92 INTEGER, ALLOCATABLE, DIMENSION(:) :: states 93 INTEGER, DIMENSION(2) :: nstates 94 LOGICAL :: floating, ghost 95 REAL(KIND=dp) :: zeff, zwfc 96 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charge_set 97 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: moment_set 98 REAL(KIND=dp), DIMENSION(3) :: rcc, ria 99 REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point 100 TYPE(atomic_kind_type), POINTER :: atomic_kind 101 TYPE(cell_type), POINTER :: cell 102 TYPE(cp_fm_struct_type), POINTER :: fm_struct 103 TYPE(cp_fm_type), POINTER :: mo_localized, momv, mvector, omvector 104 TYPE(cp_logger_type), POINTER :: logger 105 TYPE(cp_para_env_type), POINTER :: para_env 106 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments 107 TYPE(dft_control_type), POINTER :: dft_control 108 TYPE(distribution_1d_type), POINTER :: local_molecules 109 TYPE(molecule_kind_type), POINTER :: molecule_kind 110 TYPE(particle_type), POINTER :: particle_set(:) 111 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 112 113 logger => cp_get_default_logger() 114 115 CALL get_qs_env(qs_env, dft_control=dft_control) 116 nspins = dft_control%nspins 117 zwfc = 3.0_dp - REAL(nspins, KIND=dp) 118 119 CALL section_vals_val_get(loc_print_key, "MOLECULAR_MOMENTS%ORDER", i_val=norder) 120 CPASSERT(norder >= 0) 121 nm = ncoset(norder) - 1 122 123 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell) 124 particle_set => qs_loc_env%particle_set 125 para_env => qs_loc_env%para_env 126 local_molecules => qs_loc_env%local_molecules 127 molkind = SIZE(local_molecules%n_el) 128 nmols = SIZE(molecule_set) 129 ALLOCATE (charge_set(nmols), moment_set(nm, nmols)) 130 charge_set = 0.0_dp 131 moment_set = 0.0_dp 132 133 IF (norder > 0) THEN 134 CALL get_qs_env(qs_env, matrix_s=matrix_s) 135 DO imol = 1, SIZE(molecule_set) 136 molecule_kind => molecule_set(imol)%molecule_kind 137 first_atom = molecule_set(imol)%first_atom 138 CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom) 139 ! Get reference point for this molecule 140 CALL get_reference_point(rcc, qs_env=qs_env, reference=use_mom_ref_com, & 141 ref_point=ref_point, ifirst=first_atom, & 142 ilast=first_atom + natom - 1) 143 ALLOCATE (moments(nm)) 144 DO i = 1, nm 145 ALLOCATE (moments(i)%matrix) 146 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, 'MOM MAT') 147 CALL dbcsr_set(moments(i)%matrix, 0.0_dp) 148 END DO 149 ! 150 CALL build_local_moment_matrix(qs_env, moments, norder, rcc) 151 ! 152 DO ispin = 1, nspins 153 IF (ASSOCIATED(molecule_set(imol)%lmi)) THEN 154 nstates(1) = molecule_set(imol)%lmi(ispin)%nstates 155 ELSE 156 nstates(1) = 0 157 ENDIF 158 nstates(2) = para_env%mepos 159 CALL mp_maxloc(nstates, para_env%group) 160 IF (nstates(1) == 0) CYCLE 161 ns = nstates(1) 162 iproc = nstates(2) 163 ALLOCATE (states(ns)) 164 IF (iproc == para_env%mepos) THEN 165 states(:) = molecule_set(imol)%lmi(ispin)%states(:) 166 ELSE 167 states(:) = 0 168 END IF 169 CALL mp_bcast(states, iproc, para_env%group) 170 ! assemble local states for this molecule 171 mo_localized => mo_local(ispin)%matrix 172 CALL cp_fm_get_info(mo_localized, ncol_global=ncol_global, nrow_global=nrow_global) 173 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_global, ncol_global=ns, & 174 para_env=mo_localized%matrix_struct%para_env, & 175 context=mo_localized%matrix_struct%context) 176 CALL cp_fm_create(mvector, fm_struct, name="mvector") 177 CALL cp_fm_create(omvector, fm_struct, name="omvector") 178 CALL cp_fm_create(momv, fm_struct, name="omvector") 179 CALL cp_fm_struct_release(fm_struct) 180 ! 181 DO i = 1, ns 182 CALL cp_fm_to_fm(mo_localized, mvector, 1, states(i), i) 183 END DO 184 DO i = 1, nm 185 CALL cp_dbcsr_sm_fm_multiply(moments(i)%matrix, mvector, omvector, ns) 186 CALL cp_fm_schur_product(mvector, omvector, momv) 187 moment_set(i, imol) = moment_set(i, imol) - zwfc*SUM(momv%local_data) 188 END DO 189 ! 190 CALL cp_fm_release(mvector) 191 CALL cp_fm_release(omvector) 192 CALL cp_fm_release(momv) 193 DEALLOCATE (states) 194 END DO 195 DO i = 1, nm 196 CALL dbcsr_deallocate_matrix(moments(i)%matrix) 197 END DO 198 DEALLOCATE (moments) 199 END DO 200 END IF 201 ! 202 DO ikind = 1, molkind ! loop over different molecules 203 nmol = SIZE(local_molecules%list(ikind)%array) 204 DO imol = 1, nmol ! all the molecules of the kind 205 imol_now = local_molecules%list(ikind)%array(imol) ! index in the global array 206 molecule_kind => molecule_set(imol_now)%molecule_kind 207 first_atom = molecule_set(imol_now)%first_atom 208 CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom) 209 ! Get reference point for this molecule 210 CALL get_reference_point(rcc, qs_env=qs_env, reference=use_mom_ref_com, & 211 ref_point=ref_point, ifirst=first_atom, & 212 ilast=first_atom + natom - 1) 213 ! charge 214 DO iatom = 1, natom 215 i = first_atom + iatom - 1 216 atomic_kind => particle_set(i)%atomic_kind 217 CALL get_atomic_kind(atomic_kind, kind_number=akind) 218 CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating) 219 IF (.NOT. ghost .AND. .NOT. floating) THEN 220 CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff) 221 charge_set(imol_now) = charge_set(imol_now) + zeff 222 END IF 223 END DO 224 DO ispin = 1, nspins 225 IF (ASSOCIATED(molecule_set(imol_now)%lmi(ispin)%states)) THEN 226 ns = SIZE(molecule_set(imol_now)%lmi(ispin)%states) 227 charge_set(imol_now) = charge_set(imol_now) - zwfc*ns 228 END IF 229 ENDDO 230 ! 231 IF (norder > 0) THEN 232 ! nuclear contribution 233 DO i = 1, nm 234 lx = indco(1, i + 1) 235 ly = indco(2, i + 1) 236 lz = indco(3, i + 1) 237 DO iatom = 1, natom 238 j = first_atom + iatom - 1 239 atomic_kind => particle_set(j)%atomic_kind 240 CALL get_atomic_kind(atomic_kind, kind_number=akind) 241 CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating) 242 IF (.NOT. ghost .AND. .NOT. floating) THEN 243 CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff) 244 ria = particle_set(j)%r - rcc 245 ria = pbc(ria, cell) 246 moment_set(i, imol_now) = moment_set(i, imol_now) + & 247 zeff*ria(1)**lx*ria(2)**ly*ria(3)**lz 248 END IF 249 END DO 250 END DO 251 END IF 252 END DO 253 END DO 254 CALL mp_sum(moment_set, para_env%group) 255 CALL mp_sum(charge_set, para_env%group) 256 257 iounit = cp_print_key_unit_nr(logger, loc_print_key, "MOLECULAR_MOMENTS", & 258 extension=".MolMom", middle_name="MOLECULAR_MOMENTS") 259 IF (iounit > 0) THEN 260 DO i = 1, SIZE(charge_set) 261 WRITE (UNIT=iounit, FMT='(A,I6,A,F12.6)') " # molecule nr:", i, " Charge:", charge_set(I) 262 DO n = 1, norder 263 n1 = ncoset(n - 1) 264 n2 = ncoset(n) - 1 265 WRITE (UNIT=iounit, FMT='(T4,A,I2,10(T16,6F12.6))') "Order:", n, moment_set(n1:n2, i) 266 END DO 267 ENDDO 268 ENDIF 269 CALL cp_print_key_finished_output(iounit, logger, loc_print_key, & 270 "MOLECULAR_MOMENTS") 271 272 DEALLOCATE (charge_set, moment_set) 273 274 END SUBROUTINE calculate_molecular_moments 275 !------------------------------------------------------------------------------ 276 277END MODULE molecular_moments 278 279