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