Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for property calculations of excited states
10 : !> \par History
11 : !> 02.2020 Adapted from ec_properties
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE ex_property_calculation
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE cell_types, ONLY: cell_type,&
18 : pbc
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
21 : dbcsr_copy,&
22 : dbcsr_create,&
23 : dbcsr_dot,&
24 : dbcsr_p_type,&
25 : dbcsr_release,&
26 : dbcsr_set,&
27 : dbcsr_type
28 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
29 : dbcsr_deallocate_matrix_set
30 : USE cp_log_handling, ONLY: cp_get_default_logger,&
31 : cp_logger_get_default_unit_nr,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: cp_p_file,&
34 : cp_print_key_finished_output,&
35 : cp_print_key_should_output,&
36 : cp_print_key_unit_nr
37 : USE cp_result_methods, ONLY: cp_results_erase,&
38 : put_results
39 : USE cp_result_types, ONLY: cp_result_type
40 : USE distribution_1d_types, ONLY: distribution_1d_type
41 : USE input_section_types, ONLY: section_get_ival,&
42 : section_get_lval,&
43 : section_vals_get_subs_vals,&
44 : section_vals_type,&
45 : section_vals_val_get
46 : USE kinds, ONLY: default_path_length,&
47 : default_string_length,&
48 : dp
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE moments_utils, ONLY: get_reference_point
51 : USE mulliken, ONLY: mulliken_charges
52 : USE particle_types, ONLY: particle_type
53 : USE physcon, ONLY: debye
54 : USE qs_environment_types, ONLY: get_qs_env,&
55 : qs_environment_type
56 : USE qs_kind_types, ONLY: get_qs_kind,&
57 : qs_kind_type
58 : USE qs_moments, ONLY: build_local_moment_matrix
59 : USE qs_p_env_types, ONLY: qs_p_env_type
60 : USE qs_rho_types, ONLY: qs_rho_get,&
61 : qs_rho_type
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : ! Global parameters
69 :
70 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ex_property_calculation'
71 :
72 : PUBLIC :: ex_properties
73 :
74 : CONTAINS
75 :
76 : ! **************************************************************************************************
77 : !> \brief ...
78 : !> \param qs_env ...
79 : !> \param matrix_pe ...
80 : !> \param p_env ...
81 : ! **************************************************************************************************
82 550 : SUBROUTINE ex_properties(qs_env, matrix_pe, p_env)
83 : TYPE(qs_environment_type), POINTER :: qs_env
84 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe
85 : TYPE(qs_p_env_type) :: p_env
86 :
87 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ex_properties'
88 :
89 : CHARACTER(LEN=8), DIMENSION(3) :: rlab
90 : CHARACTER(LEN=default_path_length) :: filename
91 : CHARACTER(LEN=default_string_length) :: description
92 : INTEGER :: akind, handle, i, ia, iatom, idir, &
93 : ikind, iounit, ispin, maxmom, natom, &
94 : nspins, reference, unit_nr
95 : LOGICAL :: magnetic, periodic, tb
96 : REAL(KIND=dp) :: charge, dd, q, tmp
97 550 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
98 : REAL(KIND=dp), DIMENSION(3) :: cdip, pdip, pedip, rcc, rdip, ria, tdip
99 550 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
100 : TYPE(atomic_kind_type), POINTER :: atomic_kind
101 : TYPE(cell_type), POINTER :: cell
102 : TYPE(cp_logger_type), POINTER :: logger
103 : TYPE(cp_result_type), POINTER :: results
104 550 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s, moments
105 : TYPE(dbcsr_type), POINTER :: matrix_pall
106 : TYPE(dft_control_type), POINTER :: dft_control
107 : TYPE(distribution_1d_type), POINTER :: local_particles
108 : TYPE(mp_para_env_type), POINTER :: para_env
109 550 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 550 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
111 : TYPE(qs_rho_type), POINTER :: rho
112 : TYPE(section_vals_type), POINTER :: print_key
113 :
114 550 : CALL timeset(routineN, handle)
115 :
116 550 : rlab(1) = "X"
117 550 : rlab(2) = "Y"
118 550 : rlab(3) = "Z"
119 :
120 550 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
121 550 : tb = (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)
122 :
123 550 : logger => cp_get_default_logger()
124 550 : IF (logger%para_env%is_source()) THEN
125 275 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
126 : ELSE
127 : iounit = -1
128 : END IF
129 :
130 : print_key => section_vals_get_subs_vals(section_vals=qs_env%input, &
131 550 : subsection_name="DFT%PRINT%MOMENTS")
132 :
133 550 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
134 :
135 : maxmom = section_get_ival(section_vals=qs_env%input, &
136 154 : keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT")
137 : periodic = section_get_lval(section_vals=qs_env%input, &
138 154 : keyword_name="DFT%PRINT%MOMENTS%PERIODIC")
139 : reference = section_get_ival(section_vals=qs_env%input, &
140 154 : keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
141 : magnetic = section_get_lval(section_vals=qs_env%input, &
142 154 : keyword_name="DFT%PRINT%MOMENTS%MAGNETIC")
143 154 : NULLIFY (ref_point)
144 154 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
145 : unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=qs_env%input, &
146 : print_key_path="DFT%PRINT%MOMENTS", extension=".dat", &
147 154 : middle_name="moments", log_filename=.FALSE.)
148 :
149 154 : IF (iounit > 0) THEN
150 77 : IF (unit_nr /= iounit .AND. unit_nr > 0) THEN
151 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
152 : WRITE (UNIT=iounit, FMT="(/,T2,A,2(/,T3,A),/)") &
153 0 : "MOMENTS", "The electric/magnetic moments are written to file:", &
154 0 : TRIM(filename)
155 : ELSE
156 77 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS"
157 : END IF
158 : END IF
159 :
160 154 : IF (periodic) THEN
161 0 : CPABORT("Periodic moments not implemented with TDDFT")
162 : ELSE
163 154 : CPASSERT(maxmom < 2)
164 154 : CPASSERT(.NOT. magnetic)
165 154 : IF (maxmom == 1) THEN
166 154 : CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env)
167 : ! reference point
168 154 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
169 : ! nuclear contribution
170 154 : cdip = 0.0_dp
171 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
172 154 : qs_kind_set=qs_kind_set, local_particles=local_particles)
173 484 : DO ikind = 1, SIZE(local_particles%n_el)
174 725 : DO ia = 1, local_particles%n_el(ikind)
175 241 : iatom = local_particles%list(ikind)%array(ia)
176 : ! fold atomic positions back into unit cell
177 1928 : ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
178 964 : ria = ria - rcc
179 241 : atomic_kind => particle_set(iatom)%atomic_kind
180 241 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
181 241 : CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
182 1294 : cdip(1:3) = cdip(1:3) - charge*ria(1:3)
183 : END DO
184 : END DO
185 154 : CALL para_env%sum(cdip)
186 : !
187 : ! electronic contribution
188 154 : CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s=matrix_s)
189 154 : CALL qs_rho_get(rho, rho_ao=matrix_p)
190 154 : nspins = SIZE(matrix_p, 1)
191 154 : IF (tb) THEN
192 4 : ALLOCATE (matrix_pall)
193 4 : CALL dbcsr_create(matrix_pall, template=matrix_s(1)%matrix)
194 4 : CALL dbcsr_copy(matrix_pall, matrix_s(1)%matrix, "Moments")
195 4 : CALL dbcsr_set(matrix_pall, 0.0_dp)
196 8 : DO ispin = 1, nspins
197 4 : CALL dbcsr_add(matrix_pall, matrix_p(ispin)%matrix, 1.0_dp, 1.0_dp)
198 4 : CALL dbcsr_add(matrix_pall, matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
199 8 : CALL dbcsr_add(matrix_pall, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
200 : END DO
201 4 : CALL get_qs_env(qs_env=qs_env, natom=natom)
202 : ! Mulliken charges
203 12 : ALLOCATE (mcharge(natom))
204 : !
205 4 : CALL mulliken_charges(matrix_pall, matrix_s(1)%matrix, para_env, mcharge)
206 : !
207 4 : rdip = 0.0_dp
208 4 : pdip = 0.0_dp
209 4 : pedip = 0.0_dp
210 20 : DO i = 1, SIZE(particle_set)
211 128 : ria = pbc(particle_set(i)%r - rcc, cell) + rcc
212 64 : ria = ria - rcc
213 16 : q = mcharge(i)
214 68 : rdip = rdip + q*ria
215 : END DO
216 4 : CALL dbcsr_release(matrix_pall)
217 4 : DEALLOCATE (matrix_pall)
218 4 : DEALLOCATE (mcharge)
219 : ELSE
220 : ! KS-DFT
221 150 : NULLIFY (moments)
222 150 : CALL dbcsr_allocate_matrix_set(moments, 4)
223 750 : DO i = 1, 4
224 600 : ALLOCATE (moments(i)%matrix)
225 600 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
226 750 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
227 : END DO
228 150 : CALL build_local_moment_matrix(qs_env, moments, 1, ref_point=rcc)
229 : !
230 150 : rdip = 0.0_dp
231 150 : pdip = 0.0_dp
232 150 : pedip = 0.0_dp
233 316 : DO ispin = 1, nspins
234 814 : DO idir = 1, 3
235 498 : CALL dbcsr_dot(matrix_pe(ispin)%matrix, moments(idir)%matrix, tmp)
236 498 : pedip(idir) = pedip(idir) + tmp
237 498 : CALL dbcsr_dot(matrix_p(ispin)%matrix, moments(idir)%matrix, tmp)
238 498 : pdip(idir) = pdip(idir) + tmp
239 498 : CALL dbcsr_dot(p_env%p1(ispin)%matrix, moments(idir)%matrix, tmp)
240 664 : rdip(idir) = rdip(idir) + tmp
241 : END DO
242 : END DO
243 150 : CALL dbcsr_deallocate_matrix_set(moments)
244 : END IF
245 :
246 616 : tdip = -(rdip + pedip + pdip + cdip)
247 :
248 154 : IF (unit_nr > 0) THEN
249 77 : WRITE (unit_nr, "(T3,A)") "Dipoles are based on the traditional operator."
250 308 : dd = SQRT(SUM(tdip(1:3)**2))*debye
251 77 : WRITE (unit_nr, "(T3,A)") "Dipole moment [Debye]"
252 : WRITE (unit_nr, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
253 308 : (TRIM(rlab(i)), "=", tdip(i)*debye, i=1, 3), "Total=", dd
254 308 : WRITE (unit_nr, FMT="(T2,A,T61,E20.12)") ' DIPOLE : CheckSum =', SUM(ABS(tdip))
255 : END IF
256 : END IF
257 : END IF
258 :
259 154 : CALL get_qs_env(qs_env=qs_env, results=results)
260 154 : description = "[DIPOLE]"
261 154 : CALL cp_results_erase(results=results, description=description)
262 154 : CALL put_results(results=results, description=description, values=tdip(1:3))
263 :
264 : CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
265 154 : basis_section=qs_env%input, print_key_path="DFT%PRINT%MOMENTS")
266 : END IF
267 :
268 550 : CALL timestop(handle)
269 :
270 1100 : END SUBROUTINE ex_properties
271 :
272 : END MODULE ex_property_calculation
|
