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