(git:9111030)
Loading...
Searching...
No Matches
qs_vcd_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8MODULE qs_vcd_utils
9 USE cell_types, ONLY: cell_type
10 USE commutator_rpnl, ONLY: build_com_mom_nl
11 USE cp_control_types, ONLY: dft_control_type
12 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
13 dbcsr_create,&
14 dbcsr_init_p,&
15 dbcsr_p_type,&
16 dbcsr_set,&
17 dbcsr_type_antisymmetric,&
18 dbcsr_type_no_symmetry
19 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
20 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
21 dbcsr_deallocate_matrix_set
22 USE cp_files, ONLY: close_file,&
23 open_file
24 USE cp_fm_types, ONLY: cp_fm_create,&
25 cp_fm_get_info,&
26 cp_fm_get_submatrix,&
27 cp_fm_release,&
28 cp_fm_set_submatrix,&
29 cp_fm_type
30 USE cp_log_handling, ONLY: cp_get_default_logger,&
31 cp_logger_get_default_io_unit,&
32 cp_logger_type,&
33 cp_to_string
34 USE cp_output_handling, ONLY: cp_p_file,&
35 cp_print_key_finished_output,&
36 cp_print_key_generate_filename,&
37 cp_print_key_should_output,&
38 cp_print_key_unit_nr
39 USE cp_result_methods, ONLY: get_results
40 USE cp_result_types, ONLY: cp_result_type
41 USE input_constants, ONLY: use_mom_ref_user
42 USE input_section_types, ONLY: section_vals_get_subs_vals,&
43 section_vals_type,&
44 section_vals_val_get
45 USE kinds, ONLY: default_path_length,&
46 default_string_length,&
47 dp
48 USE message_passing, ONLY: mp_para_env_type
49 USE molecule_types, ONLY: molecule_type
50 USE moments_utils, ONLY: get_reference_point
51 USE orbital_pointers, ONLY: init_orbital_pointers
52 USE particle_types, ONLY: particle_type
53 USE qs_dcdr_utils, ONLY: dcdr_env_cleanup,&
54 dcdr_env_init
55 USE qs_environment_types, ONLY: get_qs_env,&
56 qs_environment_type
57 USE qs_kind_types, ONLY: qs_kind_type
58 USE qs_ks_types, ONLY: qs_ks_env_type
59 USE qs_linres_types, ONLY: vcd_env_type
60 USE qs_mo_types, ONLY: get_mo_set,&
61 mo_set_type
62 USE qs_moments, ONLY: build_local_moments_der_matrix
63 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
64 USE qs_operators_ao, ONLY: build_lin_mom_matrix
65 USE qs_vcd_ao, ONLY: build_com_rpnl_r,&
66 build_matrix_r_vhxc,&
67 build_rcore_matrix,&
68 build_rpnl_matrix,&
69 build_tr_matrix
70 USE string_utilities, ONLY: xstring
71#include "./base/base_uses.f90"
72
73 IMPLICIT NONE
74
75 PRIVATE
76 PUBLIC :: vcd_env_cleanup, vcd_env_init
77 PUBLIC :: vcd_read_restart, vcd_write_restart
78 PUBLIC :: vcd_print
79
80 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd_utils'
81
82 REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = reshape([ &
83 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
84 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
85 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp], [3, 3, 3])
86
87CONTAINS
88
89! *****************************************************************************
90!> \brief Initialize the vcd environment
91!> \param vcd_env ...
92!> \param qs_env ...
93!> \author Edward Ditler
94! **************************************************************************************************
95 SUBROUTINE vcd_env_init(vcd_env, qs_env)
96 TYPE(vcd_env_type), TARGET :: vcd_env
97 TYPE(qs_environment_type), POINTER :: qs_env
98
99 CHARACTER(LEN=*), PARAMETER :: routinen = 'vcd_env_init'
100
101 INTEGER :: handle, i, idir, ispin, j, natom, &
102 nspins, output_unit, reference, &
103 unit_number
104 LOGICAL :: explicit
105 REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
106 TYPE(cell_type), POINTER :: cell
107 TYPE(cp_logger_type), POINTER :: logger
108 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, my_matrix_hr_1d
109 TYPE(dft_control_type), POINTER :: dft_control
110 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
111 POINTER :: sab_all, sab_orb, sap_ppnl
112 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
113 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
114 TYPE(qs_ks_env_type), POINTER :: ks_env
115 TYPE(section_vals_type), POINTER :: lr_section, vcd_section
116
117 CALL timeset(routinen, handle)
118 vcd_env%do_mfp = .false.
119
120 ! Set up the logger
121 NULLIFY (logger, vcd_section, lr_section)
122 logger => cp_get_default_logger()
123 vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
124 vcd_env%output_unit = cp_print_key_unit_nr(logger, vcd_section, "PRINT%VCD", &
125 extension=".data", middle_name="vcd", log_filename=.false., &
126 file_position="REWIND", file_status="REPLACE")
127
128 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
129 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
130 extension=".linresLog")
131 unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
132
133 ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
134 CALL dcdr_env_init(vcd_env%dcdr_env, qs_env)
135 ! vcd_env%dcdr_env%output_unit = vcd_env%output_unit
136
137 IF (output_unit > 0) THEN
138 WRITE (output_unit, "(/,T20,A,/)") "*** Start NVPT/MFPT calculation ***"
139 END IF
140
141 ! Just to make sure. The memory requirements are tiny.
142 CALL init_orbital_pointers(12)
143
144 CALL section_vals_val_get(vcd_section, "DISTRIBUTED_ORIGIN", l_val=vcd_env%distributed_origin)
145 CALL section_vals_val_get(vcd_section, "ORIGIN_DEPENDENT_MFP", l_val=vcd_env%origin_dependent_op_mfp)
146
147 ! Reference point
148 vcd_env%magnetic_origin = 0._dp
149 vcd_env%spatial_origin = 0._dp
150 ! Get the magnetic origin from the input
151 CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN", i_val=reference)
152 CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", explicit=explicit)
153 IF (explicit) THEN
154 CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", r_vals=ref_point)
155 ELSE
156 IF (reference == use_mom_ref_user) THEN
157 cpabort("User-defined reference point should be given explicitly")
158 END IF
159 END IF
160
161 CALL get_reference_point(rpoint=vcd_env%magnetic_origin, qs_env=qs_env, &
162 reference=reference, &
163 ref_point=ref_point)
164
165 ! Get the spatial origin from the input
166 CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN", i_val=reference)
167 CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", explicit=explicit)
168 IF (explicit) THEN
169 CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", r_vals=ref_point)
170 ELSE
171 IF (reference == use_mom_ref_user) THEN
172 cpabort("User-defined reference point should be given explicitly")
173 END IF
174 END IF
175
176 CALL get_reference_point(rpoint=vcd_env%spatial_origin, qs_env=qs_env, &
177 reference=reference, &
178 ref_point=ref_point)
179
180 IF (vcd_env%distributed_origin .AND. any(vcd_env%magnetic_origin /= vcd_env%spatial_origin)) THEN
181 cpwarn("The magnetic and spatial origins don't match")
182 ! This is fine for NVP but will give unphysical results for MFP.
183 END IF
184
185 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
186 'The reference point is', vcd_env%dcdr_env%ref_point
187 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
188 'The magnetic origin is', vcd_env%magnetic_origin
189 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
190 'The velocity origin is', vcd_env%spatial_origin
191
192 vcd_env%magnetic_origin_atom = vcd_env%magnetic_origin
193 vcd_env%spatial_origin_atom = vcd_env%spatial_origin
194
195 CALL get_qs_env(qs_env=qs_env, &
196 ks_env=ks_env, &
197 dft_control=dft_control, &
198 sab_orb=sab_orb, &
199 sab_all=sab_all, &
200 sap_ppnl=sap_ppnl, &
201 particle_set=particle_set, &
202 matrix_ks=matrix_ks, &
203 cell=cell, &
204 qs_kind_set=qs_kind_set)
205
206 natom = SIZE(particle_set)
207 nspins = dft_control%nspins
208
209 ALLOCATE (vcd_env%apt_el_nvpt(3, 3, natom))
210 ALLOCATE (vcd_env%apt_nuc_nvpt(3, 3, natom))
211 ALLOCATE (vcd_env%apt_total_nvpt(3, 3, natom))
212 ALLOCATE (vcd_env%aat_atom_nvpt(3, 3, natom))
213 ALLOCATE (vcd_env%aat_atom_mfp(3, 3, natom))
214 vcd_env%apt_el_nvpt = 0._dp
215 vcd_env%apt_nuc_nvpt = 0._dp
216 vcd_env%apt_total_nvpt = 0._dp
217 vcd_env%aat_atom_nvpt = 0._dp
218 vcd_env%aat_atom_mfp = 0._dp
219
220 ALLOCATE (vcd_env%dCV(nspins))
221 ALLOCATE (vcd_env%dCV_prime(nspins))
222 ALLOCATE (vcd_env%op_dV(nspins))
223 ALLOCATE (vcd_env%op_dB(nspins))
224 DO ispin = 1, nspins
225 CALL cp_fm_create(vcd_env%dCV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.true.)
226 CALL cp_fm_create(vcd_env%dCV_prime(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.true.)
227 CALL cp_fm_create(vcd_env%op_dV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.true.)
228 CALL cp_fm_create(vcd_env%op_dB(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.true.)
229 END DO
230
231 ALLOCATE (vcd_env%dCB(3))
232 ALLOCATE (vcd_env%dCB_prime(3))
233 DO i = 1, 3
234 CALL cp_fm_create(vcd_env%dCB(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.true.)
235 CALL cp_fm_create(vcd_env%dCB_prime(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.true.)
236 END DO
237
238 ! DBCSR matrices
239 CALL dbcsr_allocate_matrix_set(vcd_env%moments_der, 9, 3)
240 CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_right, 9, 3)
241 CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_left, 9, 3)
242 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_difdip2, 3, 3)
243 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdV, 3)
244 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdB, 3)
245 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hxc_dsdv, nspins)
246
247 CALL dbcsr_allocate_matrix_set(vcd_env%hcom, 3)
248 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rcomr, 3, 3)
249 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rrcom, 3, 3)
250 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dcom, 3, 3)
251 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hr, nspins, 3)
252 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rh, nspins, 3)
253 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_drpnl, 3)
254 CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao, 3)
255 CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao_delta, 3)
256 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxrv, 3)
257 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_rxvr, 3, 3)
258 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxvr_r, 3, 3)
259 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_doublecom, 3, 3)
260
261 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp_33, 3, 3)
262 CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp2_33, 3, 3)
263 DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
264 DO idir = 1, 3 ! d/dx, d/dy, d/dz
265 CALL dbcsr_init_p(vcd_env%moments_der(i, idir)%matrix)
266 CALL dbcsr_init_p(vcd_env%moments_der_right(i, idir)%matrix)
267 CALL dbcsr_init_p(vcd_env%moments_der_left(i, idir)%matrix)
268
269 CALL dbcsr_create(vcd_env%moments_der(i, idir)%matrix, template=matrix_ks(1)%matrix, &
270 matrix_type=dbcsr_type_antisymmetric)
271 CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%moments_der(i, idir)%matrix, sab_orb)
272 CALL dbcsr_set(vcd_env%moments_der(i, idir)%matrix, 0.0_dp)
273
274 ! And the ones which will be multiplied by delta_(mu/nu)
275 CALL dbcsr_copy(vcd_env%moments_der_right(i, idir)%matrix, &
276 vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
277 CALL dbcsr_copy(vcd_env%moments_der_left(i, idir)%matrix, &
278 vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
279 END DO
280 END DO
281
282 DO i = 1, 3
283 DO j = 1, 3
284 CALL dbcsr_init_p(vcd_env%matrix_difdip2(i, j)%matrix)
285 CALL dbcsr_copy(vcd_env%matrix_difdip2(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
286 CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0.0_dp)
287
288 CALL dbcsr_init_p(vcd_env%matrix_nosym_temp_33(i, j)%matrix)
289 CALL dbcsr_create(vcd_env%matrix_nosym_temp_33(i, j)%matrix, template=matrix_ks(1)%matrix, &
290 matrix_type=dbcsr_type_no_symmetry)
291 CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp_33(i, j)%matrix, sab_all)
292 CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(i, j)%matrix, 0._dp)
293
294 CALL dbcsr_init_p(vcd_env%matrix_nosym_temp2_33(i, j)%matrix)
295 CALL dbcsr_create(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, template=matrix_ks(1)%matrix, &
296 matrix_type=dbcsr_type_no_symmetry)
297 CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, sab_all)
298 CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, 0._dp)
299
300 END DO
301 CALL dbcsr_init_p(vcd_env%matrix_dSdV(i)%matrix)
302 CALL dbcsr_copy(vcd_env%matrix_dSdV(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
303 CALL dbcsr_init_p(vcd_env%matrix_dSdB(i)%matrix)
304 CALL dbcsr_copy(vcd_env%matrix_dSdB(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
305 END DO
306
307 DO ispin = 1, nspins
308 CALL dbcsr_init_p(vcd_env%matrix_hxc_dsdv(ispin)%matrix)
309 CALL dbcsr_copy(vcd_env%matrix_hxc_dsdv(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
310 END DO
311
312 ! Things for op_dV
313 ! lin_mom
314 DO i = 1, 3
315 CALL dbcsr_init_p(vcd_env%dipvel_ao(i)%matrix)
316 CALL dbcsr_copy(vcd_env%dipvel_ao(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
317
318 CALL dbcsr_init_p(vcd_env%dipvel_ao_delta(i)%matrix)
319 CALL dbcsr_copy(vcd_env%dipvel_ao_delta(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
320 END DO
321
322 ! [V, r]
323 DO i = 1, 3
324 CALL dbcsr_init_p(vcd_env%hcom(i)%matrix)
325 CALL dbcsr_create(vcd_env%hcom(i)%matrix, template=matrix_ks(1)%matrix, &
326 matrix_type=dbcsr_type_antisymmetric)
327 CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%hcom(i)%matrix, sab_orb)
328
329 CALL dbcsr_init_p(vcd_env%matrix_rxrv(i)%matrix)
330 CALL dbcsr_create(vcd_env%matrix_rxrv(i)%matrix, template=matrix_ks(1)%matrix, &
331 matrix_type=dbcsr_type_antisymmetric)
332 CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_rxrv(i)%matrix, sab_orb)
333
334 DO j = 1, 3
335 CALL dbcsr_init_p(vcd_env%matrix_rcomr(i, j)%matrix)
336 CALL dbcsr_copy(vcd_env%matrix_rcomr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
337 CALL dbcsr_init_p(vcd_env%matrix_rrcom(i, j)%matrix)
338 CALL dbcsr_copy(vcd_env%matrix_rrcom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
339 CALL dbcsr_init_p(vcd_env%matrix_dcom(i, j)%matrix)
340 CALL dbcsr_copy(vcd_env%matrix_dcom(i, j)%matrix, matrix_ks(1)%matrix)
341 CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0._dp)
342
343 CALL dbcsr_init_p(vcd_env%matrix_r_rxvr(i, j)%matrix)
344 CALL dbcsr_copy(vcd_env%matrix_r_rxvr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
345 CALL dbcsr_set(vcd_env%matrix_r_rxvr(i, j)%matrix, 0._dp)
346
347 CALL dbcsr_init_p(vcd_env%matrix_rxvr_r(i, j)%matrix)
348 CALL dbcsr_copy(vcd_env%matrix_rxvr_r(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
349 CALL dbcsr_set(vcd_env%matrix_rxvr_r(i, j)%matrix, 0._dp)
350
351 CALL dbcsr_init_p(vcd_env%matrix_r_doublecom(i, j)%matrix)
352 CALL dbcsr_copy(vcd_env%matrix_r_doublecom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
353 CALL dbcsr_set(vcd_env%matrix_r_doublecom(i, j)%matrix, 0._dp)
354 END DO
355 END DO
356
357 ! matrix_hr: nonsymmetric dbcsr matrix
358 DO ispin = 1, nspins
359 DO i = 1, 3
360 CALL dbcsr_init_p(vcd_env%matrix_hr(ispin, i)%matrix)
361 CALL dbcsr_copy(vcd_env%matrix_hr(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
362
363 CALL dbcsr_init_p(vcd_env%matrix_rh(ispin, i)%matrix)
364 CALL dbcsr_copy(vcd_env%matrix_rh(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
365 END DO
366 END DO
367
368 ! drpnl for the operator
369 DO i = 1, 3
370 CALL dbcsr_init_p(vcd_env%matrix_drpnl(i)%matrix)
371 CALL dbcsr_copy(vcd_env%matrix_drpnl(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
372 END DO
373
374 ! NVP matrices
375 ! hr matrices
376 my_matrix_hr_1d => vcd_env%matrix_hr(1, 1:3)
377 CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
378 dft_control%qs_control%eps_ppnl, cell, [0._dp, 0._dp, 0._dp], &
379 direction_or=.true.)
380 CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
381 direction_or=.true., rc=[0._dp, 0._dp, 0._dp])
382 CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
383 CALL build_matrix_r_vhxc(vcd_env%matrix_hr, qs_env, [0._dp, 0._dp, 0._dp])
384
385 my_matrix_hr_1d => vcd_env%matrix_rh(1, 1:3)
386 CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
387 dft_control%qs_control%eps_ppnl, cell, [0._dp, 0._dp, 0._dp], &
388 direction_or=.false.)
389 CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
390 direction_or=.false., rc=[0._dp, 0._dp, 0._dp])
391 CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
392 CALL build_matrix_r_vhxc(vcd_env%matrix_rh, qs_env, [0._dp, 0._dp, 0._dp])
393
394 ! commutator terms
395 ! - [V, r]
396 CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
397 particle_set, cell=cell, matrix_rv=vcd_env%hcom)
398 ! <[V, r] * r> and <r * [V, r]>
399 CALL build_com_rpnl_r(vcd_env%matrix_rcomr, qs_kind_set, sab_all, sap_ppnl, &
400 dft_control%qs_control%eps_ppnl, particle_set, cell, .true.)
401 CALL build_com_rpnl_r(vcd_env%matrix_rrcom, qs_kind_set, sab_all, sap_ppnl, &
402 dft_control%qs_control%eps_ppnl, particle_set, cell, .false.)
403
404 ! lin_mom
405 CALL build_lin_mom_matrix(qs_env, vcd_env%dipvel_ao)
406
407 ! AAT
408 ! The moments are set to zero and then recomputed in the routine.
409 CALL build_local_moments_der_matrix(qs_env, moments_der=vcd_env%moments_der, &
410 nmoments_der=2, nmoments=0, ref_point=[0._dp, 0._dp, 0._dp])
411
412 ! PP terms
413 CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
414 particle_set, matrix_rxrv=vcd_env%matrix_rxrv, ref_point=[0._dp, 0._dp, 0._dp], &
415 cell=cell)
416
417 CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
418 particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
419 matrix_r_rxvr=vcd_env%matrix_r_rxvr)
420
421 CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
422 particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
423 matrix_rxvr_r=vcd_env%matrix_rxvr_r)
424
425 ! Done with NVP matrices
426
427 CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
428 "PRINT%PROGRAM_RUN_INFO")
429
430 CALL timestop(handle)
431
432 END SUBROUTINE vcd_env_init
433
434! *****************************************************************************
435!> \brief Deallocate the vcd environment
436!> \param qs_env ...
437!> \param vcd_env ...
438!> \author Edward Ditler
439! **************************************************************************************************
440 SUBROUTINE vcd_env_cleanup(qs_env, vcd_env)
441
442 TYPE(qs_environment_type), POINTER :: qs_env
443 TYPE(vcd_env_type) :: vcd_env
444
445 CHARACTER(LEN=*), PARAMETER :: routinen = 'vcd_env_cleanup'
446
447 INTEGER :: handle
448
449 CALL timeset(routinen, handle)
450
451 ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
452 CALL dcdr_env_cleanup(qs_env, vcd_env%dcdr_env)
453
454 DEALLOCATE (vcd_env%apt_el_nvpt)
455 DEALLOCATE (vcd_env%apt_nuc_nvpt)
456 DEALLOCATE (vcd_env%apt_total_nvpt)
457 DEALLOCATE (vcd_env%aat_atom_nvpt)
458 DEALLOCATE (vcd_env%aat_atom_mfp)
459
460 CALL cp_fm_release(vcd_env%dCV)
461 CALL cp_fm_release(vcd_env%dCV_prime)
462 CALL cp_fm_release(vcd_env%op_dV)
463 CALL cp_fm_release(vcd_env%op_dB)
464
465 CALL cp_fm_release(vcd_env%dCB)
466 CALL cp_fm_release(vcd_env%dCB_prime)
467
468 ! DBCSR matrices
469 ! Probably, the memory requirements could be reduced by quite a bit
470 ! by not storing each term in its own set of matrices.
471 ! On the other hand, the memory bottleneck is usually the numerical
472 ! integration grid.
473 CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der)
474 CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_right)
475 CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_left)
476 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_difdip2)
477 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdV)
478 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdB)
479 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hxc_dsdv)
480 CALL dbcsr_deallocate_matrix_set(vcd_env%hcom)
481 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rcomr)
482 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rrcom)
483 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dcom)
484 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hr)
485 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rh)
486 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_drpnl)
487 CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao)
488 CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao_delta)
489 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxrv)
490 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_rxvr)
491 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxvr_r)
492 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_doublecom)
493 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp_33)
494 CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp2_33)
495 CALL timestop(handle)
496
497 END SUBROUTINE vcd_env_cleanup
498
499! **************************************************************************************************
500!> \brief Copied from linres_read_restart
501!> \param qs_env ...
502!> \param linres_section ...
503!> \param vec ...
504!> \param lambda ...
505!> \param beta ...
506!> \param tag ...
507!> \author Edward Ditler
508! **************************************************************************************************
509 SUBROUTINE vcd_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
510 TYPE(qs_environment_type), POINTER :: qs_env
511 TYPE(section_vals_type), POINTER :: linres_section
512 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
513 INTEGER, INTENT(IN) :: lambda, beta
514 CHARACTER(LEN=*) :: tag
515
516 CHARACTER(LEN=*), PARAMETER :: routinen = 'vcd_read_restart'
517
518 CHARACTER(LEN=default_path_length) :: filename
519 CHARACTER(LEN=default_string_length) :: my_middle
520 INTEGER :: beta_tmp, handle, i, i_block, ia, ie, iostat, iounit, ispin, j, lambda_tmp, &
521 max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit
522 LOGICAL :: file_exists
523 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
524 TYPE(cp_fm_type), POINTER :: mo_coeff
525 TYPE(cp_logger_type), POINTER :: logger
526 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
527 TYPE(mp_para_env_type), POINTER :: para_env
528 TYPE(section_vals_type), POINTER :: print_key
529
530 file_exists = .false.
531
532 CALL timeset(routinen, handle)
533
534 NULLIFY (mos, para_env, logger, print_key, vecbuffer)
535 logger => cp_get_default_logger()
536
537 iounit = cp_print_key_unit_nr(logger, linres_section, &
538 "PRINT%PROGRAM_RUN_INFO", extension=".Log")
539
540 CALL get_qs_env(qs_env=qs_env, &
541 para_env=para_env, &
542 mos=mos)
543
544 nspins = SIZE(mos)
545
546 rst_unit = -1
547 IF (para_env%is_source()) THEN
548 CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
549 n_rep_val=n_rep_val)
550
551 CALL xstring(tag, ia, ie)
552 my_middle = "RESTART-"//tag(ia:ie)//trim("-")//trim(adjustl(cp_to_string(beta))) &
553 //trim("-")//trim(adjustl(cp_to_string(lambda)))
554
555 IF (n_rep_val > 0) THEN
556 CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
557 CALL xstring(filename, ia, ie)
558 filename = filename(ia:ie)//trim(my_middle)//".lr"
559 ELSE
560 ! try to read from the filename that is generated automatically from the printkey
561 print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
562 filename = cp_print_key_generate_filename(logger, print_key, &
563 extension=".lr", middle_name=trim(my_middle), my_local=.false.)
564 END IF
565 INQUIRE (file=filename, exist=file_exists)
566 !
567 ! open file
568 IF (file_exists) THEN
569 CALL open_file(file_name=trim(filename), &
570 file_action="READ", &
571 file_form="UNFORMATTED", &
572 file_position="REWIND", &
573 file_status="OLD", &
574 unit_number=rst_unit)
575
576 IF (iounit > 0) WRITE (iounit, "(T2,A)") &
577 "LINRES| Reading response wavefunctions from the restart file <"//trim(adjustl(filename))//">"
578 ELSE
579 IF (iounit > 0) WRITE (iounit, "(T2,A)") &
580 "LINRES| Restart file <"//trim(adjustl(filename))//"> not found"
581 END IF
582 END IF
583
584 CALL para_env%bcast(file_exists)
585
586 IF (file_exists) THEN
587
588 CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
589 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
590
591 ALLOCATE (vecbuffer(nao, max_block))
592 !
593 ! read headers
594 IF (rst_unit > 0) READ (rst_unit, iostat=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
595 CALL para_env%bcast(iostat)
596
597 CALL para_env%bcast(beta_tmp)
598 CALL para_env%bcast(lambda_tmp)
599 CALL para_env%bcast(nspins_tmp)
600 CALL para_env%bcast(nao_tmp)
601
602 ! check that the number nao, nmo and nspins are
603 ! the same as in the current mos
604 IF (nspins_tmp /= nspins) THEN
605 cpabort("nspins not consistent")
606 END IF
607 IF (nao_tmp /= nao) cpabort("nao not consistent")
608 ! check that it's the right file
609 ! the same as in the current mos
610 IF (lambda_tmp /= lambda) cpabort("lambda not consistent")
611 IF (beta_tmp /= beta) cpabort("beta not consistent")
612 !
613 DO ispin = 1, nspins
614 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
615 CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
616 !
617 IF (rst_unit > 0) READ (rst_unit) nmo_tmp
618 CALL para_env%bcast(nmo_tmp)
619 IF (nmo_tmp /= nmo) cpabort("nmo not consistent")
620 !
621 ! read the response
622 DO i = 1, nmo, max(max_block, 1)
623 i_block = min(max_block, nmo - i + 1)
624 DO j = 1, i_block
625 IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
626 END DO
627 CALL para_env%bcast(vecbuffer)
628 CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
629 END DO
630 END DO
631
632 IF (iostat /= 0) THEN
633 IF (iounit > 0) WRITE (iounit, "(T2,A)") &
634 "LINRES| Restart file <"//trim(adjustl(filename))//"> not found"
635 END IF
636
637 DEALLOCATE (vecbuffer)
638
639 END IF
640
641 IF (para_env%is_source()) THEN
642 IF (file_exists) CALL close_file(unit_number=rst_unit)
643 END IF
644
645 CALL timestop(handle)
646
647 END SUBROUTINE vcd_read_restart
648
649! **************************************************************************************************
650!> \brief Copied from linres_write_restart
651!> \param qs_env ...
652!> \param linres_section ...
653!> \param vec ...
654!> \param lambda ...
655!> \param beta ...
656!> \param tag ...
657!> \author Edward Ditler
658! **************************************************************************************************
659 SUBROUTINE vcd_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
660 TYPE(qs_environment_type), POINTER :: qs_env
661 TYPE(section_vals_type), POINTER :: linres_section
662 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
663 INTEGER, INTENT(IN) :: lambda, beta
664 CHARACTER(LEN=*) :: tag
665
666 CHARACTER(LEN=*), PARAMETER :: routinen = 'vcd_write_restart'
667
668 CHARACTER(LEN=default_path_length) :: filename
669 CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status
670 INTEGER :: handle, i, i_block, ia, ie, iounit, &
671 ispin, j, max_block, nao, nmo, nspins, &
672 rst_unit
673 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
674 TYPE(cp_fm_type), POINTER :: mo_coeff
675 TYPE(cp_logger_type), POINTER :: logger
676 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
677 TYPE(mp_para_env_type), POINTER :: para_env
678 TYPE(section_vals_type), POINTER :: print_key
679
680 NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
681
682 CALL timeset(routinen, handle)
683
684 logger => cp_get_default_logger()
685
686 IF (btest(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
687 used_print_key=print_key), &
688 cp_p_file)) THEN
689
690 iounit = cp_print_key_unit_nr(logger, linres_section, &
691 "PRINT%PROGRAM_RUN_INFO", extension=".Log")
692
693 CALL get_qs_env(qs_env=qs_env, &
694 mos=mos, &
695 para_env=para_env)
696
697 nspins = SIZE(mos)
698
699 my_status = "REPLACE"
700 my_pos = "REWIND"
701 CALL xstring(tag, ia, ie)
702 my_middle = "RESTART-"//tag(ia:ie)//trim("-")//trim(adjustl(cp_to_string(beta))) &
703 //trim("-")//trim(adjustl(cp_to_string(lambda)))
704 rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
705 extension=".lr", middle_name=trim(my_middle), file_status=trim(my_status), &
706 file_position=trim(my_pos), file_action="WRITE", file_form="UNFORMATTED")
707
708 filename = cp_print_key_generate_filename(logger, print_key, &
709 extension=".lr", middle_name=trim(my_middle), my_local=.false.)
710
711 IF (iounit > 0) THEN
712 WRITE (unit=iounit, fmt="(T2,A)") &
713 "LINRES| Writing response functions to the restart file <"//trim(adjustl(filename))//">"
714 END IF
715
716 !
717 ! write data to file
718 ! use the scalapack block size as a default for buffering columns
719 CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
720 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
721 ALLOCATE (vecbuffer(nao, max_block))
722
723 IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
724
725 DO ispin = 1, nspins
726 CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
727
728 IF (rst_unit > 0) WRITE (rst_unit) nmo
729
730 DO i = 1, nmo, max(max_block, 1)
731 i_block = min(max_block, nmo - i + 1)
732 CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
733 ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
734 ! to old ones, and in cases where max_block is different between runs, as might happen during
735 ! restarts with a different number of CPUs
736 DO j = 1, i_block
737 IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
738 END DO
739 END DO
740 END DO
741
742 DEALLOCATE (vecbuffer)
743
744 CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
745 "PRINT%RESTART")
746 END IF
747
748 CALL timestop(handle)
749
750 END SUBROUTINE vcd_write_restart
751
752! **************************************************************************************************
753!> \brief Print the APTs, AATs, and sum rules
754!> \param vcd_env ...
755!> \param qs_env ...
756!> \author Edward Ditler
757! **************************************************************************************************
758 SUBROUTINE vcd_print(vcd_env, qs_env)
759 TYPE(vcd_env_type) :: vcd_env
760 TYPE(qs_environment_type), POINTER :: qs_env
761
762 CHARACTER(len=*), PARAMETER :: routinen = 'vcd_print'
763
764 CHARACTER(LEN=default_string_length) :: description
765 INTEGER :: alpha, beta, delta, gamma, handle, i, l, &
766 lambda, natom, nsubset, output_unit
767 REAL(dp) :: mean, standard_deviation, &
768 standard_deviation_sum
769 REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_el_nvpt, apt_nuc_dcdr, &
770 apt_nuc_nvpt, apt_total_dcdr, &
771 apt_total_nvpt
772 REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center_dcdr, apt_subset_dcdr
773 REAL(kind=dp), DIMENSION(3, 3) :: sum_rule_0, sum_rule_0_second, &
774 sum_rule_1, sum_rule_2, &
775 sum_rule_2_second, sum_rule_3_mfp, &
776 sum_rule_3_second
777 TYPE(cp_logger_type), POINTER :: logger
778 TYPE(cp_result_type), POINTER :: results
779 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
780 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
781 TYPE(section_vals_type), POINTER :: vcd_section
782
783 CALL timeset(routinen, handle)
784
785 NULLIFY (logger)
786
787 logger => cp_get_default_logger()
788 output_unit = cp_logger_get_default_io_unit(logger)
789
790 vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
791
792 NULLIFY (particle_set)
793 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
794 natom = SIZE(particle_set)
795 nsubset = SIZE(molecule_set)
796
797 apt_el_dcdr => vcd_env%dcdr_env%apt_el_dcdr
798 apt_nuc_dcdr => vcd_env%dcdr_env%apt_nuc_dcdr
799 apt_total_dcdr => vcd_env%dcdr_env%apt_total_dcdr
800 apt_subset_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_subset
801 apt_center_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_center
802
803 apt_el_nvpt => vcd_env%apt_el_nvpt
804 apt_nuc_nvpt => vcd_env%apt_nuc_nvpt
805 apt_total_nvpt => vcd_env%apt_total_nvpt
806
807 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
808 'APT | Write the final APT matrix per atom (Position perturbation)'
809 DO l = 1, natom
810 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
811 'APT | Atom', l, ' - GAPT ', &
812 (apt_total_dcdr(1, 1, l) &
813 + apt_total_dcdr(2, 2, l) &
814 + apt_total_dcdr(3, 3, l))/3._dp
815 DO i = 1, 3
816 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l)
817 END DO
818 END DO
819
820 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
821 'NVP | Write the final APT matrix per atom (Velocity perturbation)'
822 DO l = 1, natom
823 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
824 'NVP | Atom', l, ' - GAPT ', &
825 (apt_total_nvpt(1, 1, l) &
826 + apt_total_nvpt(2, 2, l) &
827 + apt_total_nvpt(3, 3, l))/3._dp
828 DO i = 1, 3
829 IF (vcd_env%output_unit > 0) THEN
830 WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
831 "NVP | ", apt_total_nvpt(i, :, l)
832 END IF
833 END DO
834 END DO
835
836 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
837 'NVP | Write the final AAT matrix per atom (Velocity perturbation)'
838 DO l = 1, natom
839 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
840 'NVP | Atom', l
841 DO i = 1, 3
842 IF (vcd_env%output_unit > 0) THEN
843 WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
844 "NVP | ", vcd_env%aat_atom_nvpt(i, :, l)
845 END IF
846 END DO
847 END DO
848
849 IF (vcd_env%do_mfp) THEN
850 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
851 'MFP | Write the final AAT matrix per atom (Magnetic Field perturbation)'
852 DO l = 1, natom
853 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
854 'MFP | Atom', l
855 DO i = 1, 3
856 IF (vcd_env%output_unit > 0) THEN
857 WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
858 "MFP | ", vcd_env%aat_atom_mfp(i, :, l)
859 END IF
860 END DO
861 END DO
862 END IF
863
864 ! Get the dipole
865 CALL get_qs_env(qs_env, results=results)
866 description = "[DIPOLE]"
867 CALL get_results(results=results, description=description, values=vcd_env%dcdr_env%dipole_pos(1:3))
868
869 ! Sum rules [for all alpha, beta]
870 sum_rule_0 = 0._dp
871 sum_rule_1 = 0._dp
872 sum_rule_2 = 0._dp
873 sum_rule_0_second = 0._dp
874 sum_rule_2_second = 0._dp
875 sum_rule_3_second = 0._dp
876 sum_rule_3_mfp = 0._dp
877 standard_deviation = 0._dp
878 standard_deviation_sum = 0._dp
879
880 DO alpha = 1, 3
881 DO beta = 1, 3
882 ! 0: sum_lambda apt(alpha, beta, lambda)
883 DO lambda = 1, natom
884 sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
885 + apt_total_dcdr(alpha, beta, lambda)
886 sum_rule_0_second(alpha, beta) = sum_rule_0_second(alpha, beta) &
887 + apt_total_nvpt(alpha, beta, lambda)
888 END DO
889
890 ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
891 DO gamma = 1, 3
892 sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
893 + levi_civita(alpha, beta, gamma)*vcd_env%dcdr_env%dipole_pos(gamma)
894 END DO
895
896 ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda)
897 DO lambda = 1, natom
898 DO gamma = 1, 3
899 DO delta = 1, 3
900 sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) &
901 + levi_civita(beta, gamma, delta) &
902 *particle_set(lambda)%r(gamma) &
903 *apt_total_dcdr(delta, alpha, lambda)
904 sum_rule_2_second(alpha, beta) = sum_rule_2_second(alpha, beta) &
905 + levi_civita(beta, gamma, delta) &
906 *particle_set(lambda)%r(gamma) &
907 *apt_total_nvpt(delta, alpha, lambda)
908 END DO
909 END DO
910 END DO
911
912 ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
913 DO lambda = 1, natom
914 sum_rule_3_second(alpha, beta) = sum_rule_3_second(alpha, beta) &
915 + vcd_env%aat_atom_nvpt(alpha, beta, lambda)
916 ! + 2._dp*c_light_au*vcd_env%aat_atom_nvpt(alpha, beta, lambda)
917 END DO
918
919 IF (vcd_env%do_mfp) THEN
920 ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
921 DO lambda = 1, natom
922 sum_rule_3_mfp(alpha, beta) = sum_rule_3_mfp(alpha, beta) &
923 + vcd_env%aat_atom_mfp(alpha, beta, lambda)
924 END DO
925 END IF
926
927 END DO ! beta
928 END DO ! alpha
929
930 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") "APT | Position perturbation sum rules"
931 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") &
932 "APT |", " Total APT", "Dipole", "R * APT", "AAT"
933 standard_deviation_sum = 0._dp
934 DO alpha = 1, 3
935 DO beta = 1, 3
936 mean = (sum_rule_1(alpha, beta) + sum_rule_2(alpha, beta) + sum_rule_3_mfp(alpha, beta))/3
937 standard_deviation = &
938 sqrt((sum_rule_1(alpha, beta)**2 + sum_rule_2(alpha, beta)**2 + sum_rule_3_mfp(alpha, beta)**2)/3 &
939 - mean**2)
940 standard_deviation_sum = standard_deviation_sum + standard_deviation
941
942 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, &
943 "(A,I3,I3,F15.6,F15.6,F15.6,F15.6,F15.6)") &
944 "APT | ", &
945 alpha, beta, &
946 sum_rule_0(alpha, beta), &
947 sum_rule_1(alpha, beta), &
948 sum_rule_2(alpha, beta), &
949 sum_rule_3_mfp(alpha, beta), &
950 standard_deviation
951 END DO
952 END DO
953 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
954
955 IF (vcd_env%output_unit > 0) THEN
956 WRITE (vcd_env%output_unit, "(A)") "NVP | Velocity perturbation sum rules"
957 WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") "NVP |", " Total APT", "Dipole", "R * APT", "AAT"
958 END IF
959
960 standard_deviation_sum = 0._dp
961 DO alpha = 1, 3
962 DO beta = 1, 3
963 mean = (sum_rule_1(alpha, beta) + sum_rule_2_second(alpha, beta) + sum_rule_3_second(alpha, beta))/3
964 standard_deviation = &
965 sqrt((sum_rule_1(alpha, beta)**2 + sum_rule_2_second(alpha, beta)**2 + sum_rule_3_second(alpha, beta)**2)/3 &
966 - mean**2)
967 standard_deviation_sum = standard_deviation_sum + standard_deviation
968 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, &
969 "(A,I3,I3,F15.6,F15.6,F15.6,F15.6,F15.6)") &
970 "NVP | ", &
971 alpha, &
972 beta, &
973 sum_rule_0_second(alpha, beta), &
974 sum_rule_1(alpha, beta), &
975 sum_rule_2_second(alpha, beta), &
976 sum_rule_3_second(alpha, beta), &
977 standard_deviation
978 END DO
979 END DO
980 IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
981
982 CALL timestop(handle)
983 END SUBROUTINE vcd_print
984
985END MODULE qs_vcd_utils
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
cp_result_methods::get_results
Definition cp_result_methods.F:43
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition commutator_rpnl.F:14
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition commutator_rpnl.F:439
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:504
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:754
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition cp_result_types.F:18
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::use_mom_ref_user
integer, parameter, public use_mom_ref_user
Definition input_constants.F:510
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
moments_utils
Calculates the moment integrals <a|r^m|b>
Definition moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition moments_utils.F:61
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_dcdr_utils
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition qs_dcdr_utils.F:13
qs_dcdr_utils::dcdr_env_cleanup
subroutine, public dcdr_env_cleanup(qs_env, dcdr_env)
Deallocate the dcdr environment.
Definition qs_dcdr_utils.F:1015
qs_dcdr_utils::dcdr_env_init
subroutine, public dcdr_env_init(dcdr_env, qs_env)
Initialize the dcdr environment.
Definition qs_dcdr_utils.F:654
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:533
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_ks_types
Definition qs_ks_types.F:15
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_local_moments_der_matrix
subroutine, public build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, ref_point, moments)
Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b > Optionally sto...
Definition qs_moments.F:825
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_operators_ao
Definition qs_operators_ao.F:13
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition qs_operators_ao.F:74
qs_vcd_ao
Definition qs_vcd_ao.F:7
qs_vcd_ao::build_rpnl_matrix
subroutine, public build_rpnl_matrix(matrix_rv, qs_kind_set, particle_set, sab_all, sap_ppnl, eps_ppnl, cell, ref_point, direction_or)
Product of r with V_nl. Adapted from build_com_rpnl.
Definition qs_vcd_ao.F:188
qs_vcd_ao::build_matrix_r_vhxc
subroutine, public build_matrix_r_vhxc(matrix_rv, qs_env, rc)
Commutator of the Hartree+XC potentials with r.
Definition qs_vcd_ao.F:979
qs_vcd_ao::build_rcore_matrix
subroutine, public build_rcore_matrix(matrix_rcore, qs_env, qs_kind_set, basis_type, sab_nl, rf)
Commutator of the of the local part of the pseudopotential with r.
Definition qs_vcd_ao.F:629
qs_vcd_ao::build_com_rpnl_r
subroutine, public build_com_rpnl_r(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, direction_or)
Builds the [Vnl, r] * r from either side.
Definition qs_vcd_ao.F:1550
qs_vcd_ao::build_tr_matrix
subroutine, public build_tr_matrix(matrix_tr, qs_env, qs_kind_set, basis_type, sab_nl, direction_or, rc)
Calculation of the product Tr or rT over Cartesian Gaussian functions.
Definition qs_vcd_ao.F:389
qs_vcd_utils
Definition qs_vcd_utils.F:8
qs_vcd_utils::vcd_print
subroutine, public vcd_print(vcd_env, qs_env)
Print the APTs, AATs, and sum rules.
Definition qs_vcd_utils.F:759
qs_vcd_utils::vcd_env_cleanup
subroutine, public vcd_env_cleanup(qs_env, vcd_env)
Deallocate the vcd environment.
Definition qs_vcd_utils.F:441
qs_vcd_utils::vcd_write_restart
subroutine, public vcd_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_write_restart.
Definition qs_vcd_utils.F:660
qs_vcd_utils::vcd_env_init
subroutine, public vcd_env_init(vcd_env, qs_env)
Initialize the vcd environment.
Definition qs_vcd_utils.F:96
qs_vcd_utils::vcd_read_restart
subroutine, public vcd_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_read_restart.
Definition qs_vcd_utils.F:510
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::xstring
elemental subroutine, public xstring(string, ia, ib)
...
Definition string_utilities.F:3382
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:751
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
molecule_types::molecule_type
Definition molecule_types.F:159
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_linres_types::vcd_env_type
Definition qs_linres_types.F:305
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67