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
10 : ! **************************************************************************************************
11 : MODULE qs_loc_dipole
12 : USE atomic_kind_types, ONLY: get_atomic_kind
13 : USE cell_types, ONLY: cell_type,&
14 : pbc
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE cp_log_handling, ONLY: cp_logger_type
17 : USE cp_output_handling, ONLY: cp_iter_string,&
18 : cp_p_file,&
19 : cp_print_key_finished_output,&
20 : cp_print_key_should_output,&
21 : cp_print_key_unit_nr
22 : USE cp_result_methods, ONLY: cp_results_erase,&
23 : get_results,&
24 : put_results
25 : USE cp_result_types, ONLY: cp_result_type
26 : USE input_section_types, ONLY: section_get_ival,&
27 : section_vals_get_subs_vals,&
28 : section_vals_type,&
29 : section_vals_val_get
30 : USE kinds, ONLY: default_string_length,&
31 : dp
32 : USE mathconstants, ONLY: twopi,&
33 : z_one
34 : USE moments_utils, ONLY: get_reference_point
35 : USE particle_types, ONLY: particle_type
36 : USE physcon, ONLY: debye
37 : USE qs_environment_types, ONLY: get_qs_env,&
38 : qs_environment_type
39 : USE qs_kind_types, ONLY: get_qs_kind,&
40 : qs_kind_type
41 : USE qs_loc_types, ONLY: qs_loc_env_type
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 : PRIVATE
46 :
47 : ! Global parameters
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_dipole'
49 : PUBLIC :: loc_dipole
50 :
51 : ! **************************************************************************************************
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Computes and prints the Dipole (using localized charges)
57 : !> \param input ...
58 : !> \param dft_control ...
59 : !> \param qs_loc_env ...
60 : !> \param logger ...
61 : !> \param qs_env the qs_env in which the qs_env lives
62 : ! **************************************************************************************************
63 80 : SUBROUTINE loc_dipole(input, dft_control, qs_loc_env, logger, qs_env)
64 : TYPE(section_vals_type), POINTER :: input
65 : TYPE(dft_control_type), POINTER :: dft_control
66 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
67 : TYPE(cp_logger_type), POINTER :: logger
68 : TYPE(qs_environment_type), POINTER :: qs_env
69 :
70 : CHARACTER(len=*), PARAMETER :: routineN = 'loc_dipole'
71 :
72 : CHARACTER(LEN=default_string_length) :: description, descriptionThisDip, iter
73 : COMPLEX(KIND=dp) :: zeta
74 : COMPLEX(KIND=dp), DIMENSION(3) :: ggamma, zphase
75 : INTEGER :: handle, i, ikind, ispins, j, n_rep, &
76 : reference, unit_nr
77 : LOGICAL :: do_berry, first_time, floating, ghost
78 : REAL(KIND=dp) :: charge_tot, theta, zeff, zwfc
79 : REAL(KIND=dp), DIMENSION(3) :: ci, dipole, dipole_old, gvec, rcc, ria
80 80 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
81 : TYPE(cell_type), POINTER :: cell
82 : TYPE(cp_result_type), POINTER :: results
83 80 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
84 80 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
85 : TYPE(section_vals_type), POINTER :: print_key
86 :
87 80 : CALL timeset(routineN, handle)
88 :
89 80 : print_key => section_vals_get_subs_vals(input, "DFT%LOCALIZE%PRINT%TOTAL_DIPOLE")
90 80 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, first_time=first_time) &
91 : , cp_p_file)) THEN
92 16 : NULLIFY (cell, particle_set, qs_kind_set, ref_point, results)
93 : CALL get_qs_env(qs_env=qs_env, &
94 : cell=cell, &
95 : particle_set=particle_set, &
96 : qs_kind_set=qs_kind_set, &
97 16 : results=results)
98 :
99 16 : reference = section_get_ival(print_key, keyword_name="REFERENCE")
100 16 : CALL section_vals_val_get(print_key, "REF_POINT", r_vals=ref_point)
101 16 : CALL section_vals_val_get(print_key, "PERIODIC", l_val=do_berry)
102 16 : description = '[DIPOLE]'
103 16 : descriptionThisDip = '[TOTAL_DIPOLE]'
104 16 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
105 :
106 16 : dipole = 0.0_dp
107 16 : IF (do_berry) THEN
108 64 : rcc = pbc(rcc, cell)
109 16 : charge_tot = REAL(dft_control%charge, KIND=dp)
110 256 : ria = twopi*MATMUL(cell%h_inv, rcc)
111 64 : zphase = CMPLX(COS(ria), SIN(ria), KIND=dp)**charge_tot
112 64 : ggamma = z_one
113 :
114 : ! Nuclear charges
115 102 : DO i = 1, SIZE(particle_set)
116 86 : CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
117 86 : CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost, floating=floating)
118 188 : IF (.NOT. ghost .AND. .NOT. floating) THEN
119 86 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=zeff)
120 86 : ria = pbc(particle_set(i)%r, cell)
121 344 : DO j = 1, 3
122 1032 : gvec = twopi*cell%h_inv(j, :)
123 1032 : theta = SUM(ria(:)*gvec(:))
124 258 : zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(zeff)
125 344 : ggamma(j) = ggamma(j)*zeta
126 : END DO
127 : END IF
128 : END DO
129 :
130 : ! Charges of the wfc involved
131 : ! Warning, this assumes the same occupation for all states
132 16 : zwfc = 3.0_dp - REAL(dft_control%nspins, dp)
133 :
134 32 : DO ispins = 1, dft_control%nspins
135 148 : DO i = 1, SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array, 2)
136 116 : ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3, i), cell)
137 480 : DO j = 1, 3
138 1392 : gvec = twopi*cell%h_inv(j, :)
139 1392 : theta = SUM(ria(:)*gvec(:))
140 348 : zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-zwfc)
141 464 : ggamma(j) = ggamma(j)*zeta
142 : END DO
143 : END DO
144 : END DO
145 64 : ggamma = ggamma*zphase
146 64 : ci = AIMAG(LOG(ggamma))/twopi
147 208 : dipole = MATMUL(cell%hmat, ci)
148 : ELSE
149 : ! Charges of the atoms involved
150 0 : DO i = 1, SIZE(particle_set)
151 0 : CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
152 0 : CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost, floating=floating)
153 0 : IF (.NOT. ghost .AND. .NOT. floating) THEN
154 0 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=zeff)
155 0 : ria = pbc(particle_set(i)%r, cell)
156 0 : dipole = dipole + zeff*(ria - rcc)
157 : END IF
158 : END DO
159 :
160 : ! Charges of the wfc involved
161 : ! Warning, this assumes the same occupation for all states
162 0 : zwfc = 3.0_dp - REAL(dft_control%nspins, dp)
163 :
164 0 : DO ispins = 1, dft_control%nspins
165 0 : DO i = 1, SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array, 2)
166 0 : ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3, i), cell)
167 0 : dipole = dipole - zwfc*(ria - rcc)
168 : END DO
169 : END DO
170 : END IF
171 :
172 : ! Print and possibly store results
173 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".Dipole", &
174 16 : middle_name="TOTAL_DIPOLE")
175 16 : IF (unit_nr > 0) THEN
176 8 : IF (first_time) THEN
177 : WRITE (unit=unit_nr, fmt="(A,T31,A,T88,A,T136,A)") &
178 8 : "# iter_level", "dipole(x,y,z)[atomic units]", &
179 8 : "dipole(x,y,z)[debye]", &
180 16 : "delta_dipole(x,y,z)[atomic units]"
181 : END IF
182 8 : iter = cp_iter_string(logger%iter_info)
183 8 : CALL get_results(results, descriptionThisDip, n_rep=n_rep)
184 8 : IF (n_rep == 0) THEN
185 3 : dipole_old = 0._dp
186 : ELSE
187 5 : CALL get_results(results, descriptionThisDip, dipole_old, nval=n_rep)
188 : END IF
189 8 : IF (do_berry) THEN
190 : WRITE (unit=unit_nr, fmt="(a,9(es18.8))") &
191 88 : iter(1:15), dipole, dipole*debye, pbc(dipole - dipole_old, cell)
192 : ELSE
193 : WRITE (unit=unit_nr, fmt="(a,9(es18.8))") &
194 0 : iter(1:15), dipole, dipole*debye, (dipole - dipole_old)
195 : END IF
196 : END IF
197 16 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
198 16 : CALL cp_results_erase(results, description)
199 16 : CALL put_results(results, description, dipole)
200 16 : CALL cp_results_erase(results, descriptionThisDip)
201 16 : CALL put_results(results, descriptionThisDip, dipole)
202 : END IF
203 :
204 80 : CALL timestop(handle)
205 :
206 80 : END SUBROUTINE loc_dipole
207 :
208 : END MODULE qs_loc_dipole
|
