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_dipoles 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_log_handling, ONLY: cp_get_default_logger,& 16 cp_logger_type 17 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 18 cp_print_key_unit_nr 19 USE cp_para_types, ONLY: cp_para_env_type 20 USE distribution_1d_types, ONLY: distribution_1d_type 21 USE input_section_types, ONLY: section_get_ival,& 22 section_vals_type,& 23 section_vals_val_get 24 USE kinds, ONLY: dp 25 USE mathconstants, ONLY: twopi 26 USE message_passing, ONLY: mp_sum 27 USE molecule_kind_types, ONLY: get_molecule_kind,& 28 molecule_kind_type 29 USE molecule_types, ONLY: molecule_type 30 USE moments_utils, ONLY: get_reference_point 31 USE particle_types, ONLY: particle_type 32 USE physcon, ONLY: debye 33 USE qs_environment_types, ONLY: get_qs_env,& 34 qs_environment_type 35 USE qs_kind_types, ONLY: get_qs_kind,& 36 qs_kind_type 37 USE qs_loc_types, ONLY: qs_loc_env_new_type 38#include "./base/base_uses.f90" 39 40 IMPLICIT NONE 41 42 PRIVATE 43 44 ! *** Public *** 45 PUBLIC :: calculate_molecular_dipole 46 47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecular_dipoles' 48 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief maps wfc's to molecules and also prints molecular dipoles 53!> \param qs_env the qs_env in which the qs_env lives 54!> \param qs_loc_env ... 55!> \param loc_print_key ... 56!> \param molecule_set ... 57! ************************************************************************************************** 58 SUBROUTINE calculate_molecular_dipole(qs_env, qs_loc_env, loc_print_key, molecule_set) 59 TYPE(qs_environment_type), POINTER :: qs_env 60 TYPE(qs_loc_env_new_type), INTENT(IN) :: qs_loc_env 61 TYPE(section_vals_type), POINTER :: loc_print_key 62 TYPE(molecule_type), POINTER :: molecule_set(:) 63 64 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_molecular_dipole', & 65 routineP = moduleN//':'//routineN 66 67 COMPLEX(KIND=dp) :: zeta 68 COMPLEX(KIND=dp), DIMENSION(3) :: ggamma, zphase 69 INTEGER :: akind, first_atom, i, iatom, ikind, & 70 imol, imol_now, iounit, ispin, istate, & 71 j, natom, nkind, nmol, nspins, nstate, & 72 reference 73 LOGICAL :: do_berry, floating, ghost 74 REAL(KIND=dp) :: charge_tot, dipole(3), ria(3), theta, & 75 zeff, zwfc 76 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charge_set 77 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dipole_set 78 REAL(KIND=dp), DIMENSION(3) :: ci, gvec, rcc 79 REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point 80 REAL(KIND=dp), DIMENSION(:, :), POINTER :: center(:, :) 81 TYPE(atomic_kind_type), POINTER :: atomic_kind 82 TYPE(cell_type), POINTER :: cell 83 TYPE(cp_logger_type), POINTER :: logger 84 TYPE(cp_para_env_type), POINTER :: para_env 85 TYPE(dft_control_type), POINTER :: dft_control 86 TYPE(distribution_1d_type), POINTER :: local_molecules 87 TYPE(molecule_kind_type), POINTER :: molecule_kind 88 TYPE(particle_type), POINTER :: particle_set(:) 89 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 90 91 logger => cp_get_default_logger() 92 93 CALL get_qs_env(qs_env, dft_control=dft_control) 94 nspins = dft_control%nspins 95 96 ! Setup reference point 97 reference = section_get_ival(loc_print_key, keyword_name="MOLECULAR_DIPOLES%REFERENCE") 98 CALL section_vals_val_get(loc_print_key, "MOLECULAR_DIPOLES%REF_POINT", r_vals=ref_point) 99 CALL section_vals_val_get(loc_print_key, "MOLECULAR_DIPOLES%PERIODIC", l_val=do_berry) 100 101 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell) 102 particle_set => qs_loc_env%particle_set 103 para_env => qs_loc_env%para_env 104 local_molecules => qs_loc_env%local_molecules 105 nkind = SIZE(local_molecules%n_el) 106 zwfc = 3.0_dp - REAL(nspins, KIND=dp) 107 108 ALLOCATE (dipole_set(3, SIZE(molecule_set))) 109 ALLOCATE (charge_set(SIZE(molecule_set))) 110 dipole_set = 0.0_dp 111 charge_set = 0.0_dp 112 113 DO ispin = 1, nspins 114 center => qs_loc_env%localized_wfn_control%centers_set(ispin)%array 115 nstate = SIZE(center, 2) 116 DO ikind = 1, nkind ! loop over different molecules 117 nmol = SIZE(local_molecules%list(ikind)%array) 118 DO imol = 1, nmol ! all the molecules of the kind 119 imol_now = local_molecules%list(ikind)%array(imol) ! index in the global array 120 IF (.NOT. ASSOCIATED(molecule_set(imol_now)%lmi(ispin)%states)) CYCLE 121 molecule_kind => molecule_set(imol_now)%molecule_kind 122 first_atom = molecule_set(imol_now)%first_atom 123 CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom) 124 125 ! Get reference point for this molecule 126 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, & 127 ref_point=ref_point, ifirst=first_atom, & 128 ilast=first_atom + natom - 1) 129 130 dipole = 0.0_dp 131 IF (do_berry) THEN 132 rcc = pbc(rcc, cell) 133 ! Find out the total charge of the molecule 134 DO iatom = 1, natom 135 i = first_atom + iatom - 1 136 atomic_kind => particle_set(i)%atomic_kind 137 CALL get_atomic_kind(atomic_kind, kind_number=akind) 138 CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating) 139 IF (.NOT. ghost .AND. .NOT. floating) THEN 140 CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff) 141 charge_set(imol_now) = charge_set(imol_now) + zeff 142 END IF 143 END DO 144 ! Charges of the wfc involved 145 DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states) 146 charge_set(imol_now) = charge_set(imol_now) - zwfc 147 ENDDO 148 149 charge_tot = charge_set(imol_now) 150 ria = twopi*MATMUL(cell%h_inv, rcc) 151 zphase = CMPLX(COS(ria), SIN(ria), KIND=dp)**charge_tot 152 ggamma = CMPLX(1.0_dp, 0.0_dp, KIND=dp) 153 154 ! Nuclear charges 155 IF (ispin == 1) THEN 156 DO iatom = 1, natom 157 i = first_atom + iatom - 1 158 atomic_kind => particle_set(i)%atomic_kind 159 CALL get_atomic_kind(atomic_kind, kind_number=akind) 160 CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating) 161 IF (.NOT. ghost .AND. .NOT. floating) THEN 162 CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff) 163 ria = pbc(particle_set(i)%r, cell) 164 DO j = 1, 3 165 gvec = twopi*cell%h_inv(j, :) 166 theta = SUM(ria(:)*gvec(:)) 167 zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(zeff) 168 ggamma(j) = ggamma(j)*zeta 169 END DO 170 END IF 171 END DO 172 END IF 173 174 ! Charges of the wfc involved 175 DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states) 176 i = molecule_set(imol_now)%lmi(ispin)%states(istate) 177 ria = pbc(center(1:3, i), cell) 178 DO j = 1, 3 179 gvec = twopi*cell%h_inv(j, :) 180 theta = SUM(ria(:)*gvec(:)) 181 zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-zwfc) 182 ggamma(j) = ggamma(j)*zeta 183 END DO 184 ENDDO 185 186 ggamma = ggamma*zphase 187 ci = AIMAG(LOG(ggamma))/twopi 188 dipole = MATMUL(cell%hmat, ci) 189 ELSE 190 IF (ispin == 1) THEN 191 ! Nuclear charges 192 DO iatom = 1, natom 193 i = first_atom + iatom - 1 194 atomic_kind => particle_set(i)%atomic_kind 195 CALL get_atomic_kind(atomic_kind, kind_number=akind) 196 CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating) 197 IF (.NOT. ghost .AND. .NOT. floating) THEN 198 CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff) 199 ria = pbc(particle_set(i)%r, cell) - rcc 200 dipole = dipole + zeff*(ria - rcc) 201 charge_set(imol_now) = charge_set(imol_now) + zeff 202 END IF 203 END DO 204 END IF 205 ! Charges of the wfc involved 206 DO istate = 1, SIZE(molecule_set(imol_now)%lmi(ispin)%states) 207 i = molecule_set(imol_now)%lmi(ispin)%states(istate) 208 ria = pbc(center(1:3, i), cell) 209 dipole = dipole - zwfc*(ria - rcc) 210 charge_set(imol_now) = charge_set(imol_now) - zwfc 211 ENDDO 212 END IF 213 dipole_set(:, imol_now) = dipole_set(:, imol_now) + dipole ! a.u. 214 ENDDO 215 ENDDO 216 END DO 217 CALL mp_sum(dipole_set, para_env%group) 218 CALL mp_sum(charge_set, para_env%group) 219 220 iounit = cp_print_key_unit_nr(logger, loc_print_key, "MOLECULAR_DIPOLES", & 221 extension=".MolDip", middle_name="MOLECULAR_DIPOLES") 222 IF (iounit > 0) THEN 223 WRITE (UNIT=iounit, FMT='(A80)') & 224 "# molecule nr, charge, dipole vector, dipole[Debye]" 225 dipole_set(:, :) = dipole_set(:, :)*debye ! Debye 226 DO I = 1, SIZE(dipole_set, 2) 227 WRITE (UNIT=iounit, FMT='(T8,I6,T21,5F12.6)') I, charge_set(I), dipole_set(1:3, I), & 228 SQRT(DOT_PRODUCT(dipole_set(1:3, I), dipole_set(1:3, I))) 229 ENDDO 230 WRITE (UNIT=iounit, FMT="(T2,A,T61,E20.12)") ' DIPOLE : CheckSum =', SUM(dipole_set) 231 ENDIF 232 CALL cp_print_key_finished_output(iounit, logger, loc_print_key, & 233 "MOLECULAR_DIPOLES") 234 235 DEALLOCATE (dipole_set, charge_set) 236 237 END SUBROUTINE calculate_molecular_dipole 238 !------------------------------------------------------------------------------ 239 240END MODULE molecular_dipoles 241 242