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