Line data Source code
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 :
8 : MODULE 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 :
87 : CONTAINS
88 :
89 : ! *****************************************************************************
90 : !> \brief Initialize the vcd environment
91 : !> \param vcd_env ...
92 : !> \param qs_env ...
93 : !> \author Edward Ditler
94 : ! **************************************************************************************************
95 2 : 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 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
106 : TYPE(cell_type), POINTER :: cell
107 : TYPE(cp_logger_type), POINTER :: logger
108 2 : 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 2 : POINTER :: sab_all, sab_orb, sap_ppnl
112 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
113 2 : 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 2 : CALL timeset(routineN, handle)
118 2 : vcd_env%do_mfp = .FALSE.
119 :
120 : ! Set up the logger
121 2 : NULLIFY (logger, vcd_section, lr_section)
122 2 : logger => cp_get_default_logger()
123 2 : 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 2 : file_position="REWIND", file_status="REPLACE")
127 :
128 2 : 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 2 : extension=".linresLog")
131 2 : 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 2 : CALL dcdr_env_init(vcd_env%dcdr_env, qs_env)
135 : ! vcd_env%dcdr_env%output_unit = vcd_env%output_unit
136 :
137 2 : IF (output_unit > 0) THEN
138 1 : WRITE (output_unit, "(/,T20,A,/)") "*** Start NVPT/MFPT calculation ***"
139 : END IF
140 :
141 : ! Just to make sure. The memory requirements are tiny.
142 2 : CALL init_orbital_pointers(12)
143 :
144 2 : CALL section_vals_val_get(vcd_section, "DISTRIBUTED_ORIGIN", l_val=vcd_env%distributed_origin)
145 2 : CALL section_vals_val_get(vcd_section, "ORIGIN_DEPENDENT_MFP", l_val=vcd_env%origin_dependent_op_mfp)
146 :
147 : ! Reference point
148 8 : vcd_env%magnetic_origin = 0._dp
149 8 : vcd_env%spatial_origin = 0._dp
150 : ! Get the magnetic origin from the input
151 2 : CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN", i_val=reference)
152 2 : CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", explicit=explicit)
153 2 : IF (explicit) THEN
154 0 : CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", r_vals=ref_point)
155 : ELSE
156 2 : IF (reference == use_mom_ref_user) THEN
157 0 : 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 2 : ref_point=ref_point)
164 :
165 : ! Get the spatial origin from the input
166 2 : CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN", i_val=reference)
167 2 : CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", explicit=explicit)
168 2 : IF (explicit) THEN
169 0 : CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", r_vals=ref_point)
170 : ELSE
171 2 : IF (reference == use_mom_ref_user) THEN
172 0 : 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 2 : ref_point=ref_point)
179 :
180 8 : IF (vcd_env%distributed_origin .AND. ANY(vcd_env%magnetic_origin /= vcd_env%spatial_origin)) THEN
181 0 : 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 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
186 1 : 'The reference point is', vcd_env%dcdr_env%ref_point
187 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
188 1 : 'The magnetic origin is', vcd_env%magnetic_origin
189 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
190 1 : 'The velocity origin is', vcd_env%spatial_origin
191 :
192 8 : vcd_env%magnetic_origin_atom = vcd_env%magnetic_origin
193 8 : 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 2 : qs_kind_set=qs_kind_set)
205 :
206 2 : natom = SIZE(particle_set)
207 2 : nspins = dft_control%nspins
208 :
209 6 : ALLOCATE (vcd_env%apt_el_nvpt(3, 3, natom))
210 4 : ALLOCATE (vcd_env%apt_nuc_nvpt(3, 3, natom))
211 4 : ALLOCATE (vcd_env%apt_total_nvpt(3, 3, natom))
212 4 : ALLOCATE (vcd_env%aat_atom_nvpt(3, 3, natom))
213 4 : ALLOCATE (vcd_env%aat_atom_mfp(3, 3, natom))
214 80 : vcd_env%apt_el_nvpt = 0._dp
215 80 : vcd_env%apt_nuc_nvpt = 0._dp
216 80 : vcd_env%apt_total_nvpt = 0._dp
217 80 : vcd_env%aat_atom_nvpt = 0._dp
218 80 : vcd_env%aat_atom_mfp = 0._dp
219 :
220 8 : ALLOCATE (vcd_env%dCV(nspins))
221 6 : ALLOCATE (vcd_env%dCV_prime(nspins))
222 6 : ALLOCATE (vcd_env%op_dV(nspins))
223 6 : ALLOCATE (vcd_env%op_dB(nspins))
224 4 : DO ispin = 1, nspins
225 2 : CALL cp_fm_create(vcd_env%dCV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.TRUE.)
226 2 : CALL cp_fm_create(vcd_env%dCV_prime(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.TRUE.)
227 2 : CALL cp_fm_create(vcd_env%op_dV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.TRUE.)
228 4 : 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 8 : ALLOCATE (vcd_env%dCB(3))
232 8 : ALLOCATE (vcd_env%dCB_prime(3))
233 8 : DO i = 1, 3
234 6 : CALL cp_fm_create(vcd_env%dCB(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct, set_zero=.TRUE.)
235 8 : 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 2 : CALL dbcsr_allocate_matrix_set(vcd_env%moments_der, 9, 3)
240 2 : CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_right, 9, 3)
241 2 : CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_left, 9, 3)
242 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_difdip2, 3, 3)
243 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdV, 3)
244 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdB, 3)
245 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hxc_dsdv, nspins)
246 :
247 2 : CALL dbcsr_allocate_matrix_set(vcd_env%hcom, 3)
248 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rcomr, 3, 3)
249 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rrcom, 3, 3)
250 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dcom, 3, 3)
251 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hr, nspins, 3)
252 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rh, nspins, 3)
253 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_drpnl, 3)
254 2 : CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao, 3)
255 2 : CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao_delta, 3)
256 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxrv, 3)
257 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_rxvr, 3, 3)
258 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxvr_r, 3, 3)
259 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_doublecom, 3, 3)
260 :
261 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp_33, 3, 3)
262 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp2_33, 3, 3)
263 20 : DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
264 74 : DO idir = 1, 3 ! d/dx, d/dy, d/dz
265 54 : CALL dbcsr_init_p(vcd_env%moments_der(i, idir)%matrix)
266 54 : CALL dbcsr_init_p(vcd_env%moments_der_right(i, idir)%matrix)
267 54 : 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 54 : matrix_type=dbcsr_type_antisymmetric)
271 54 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%moments_der(i, idir)%matrix, sab_orb)
272 54 : 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 54 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
277 : CALL dbcsr_copy(vcd_env%moments_der_left(i, idir)%matrix, &
278 72 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
279 : END DO
280 : END DO
281 :
282 8 : DO i = 1, 3
283 24 : DO j = 1, 3
284 18 : CALL dbcsr_init_p(vcd_env%matrix_difdip2(i, j)%matrix)
285 18 : CALL dbcsr_copy(vcd_env%matrix_difdip2(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
286 18 : CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0.0_dp)
287 :
288 18 : 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 18 : matrix_type=dbcsr_type_no_symmetry)
291 18 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp_33(i, j)%matrix, sab_all)
292 18 : CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(i, j)%matrix, 0._dp)
293 :
294 18 : 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 18 : matrix_type=dbcsr_type_no_symmetry)
297 18 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, sab_all)
298 24 : CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, 0._dp)
299 :
300 : END DO
301 6 : CALL dbcsr_init_p(vcd_env%matrix_dSdV(i)%matrix)
302 6 : CALL dbcsr_copy(vcd_env%matrix_dSdV(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
303 6 : CALL dbcsr_init_p(vcd_env%matrix_dSdB(i)%matrix)
304 8 : CALL dbcsr_copy(vcd_env%matrix_dSdB(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
305 : END DO
306 :
307 4 : DO ispin = 1, nspins
308 2 : CALL dbcsr_init_p(vcd_env%matrix_hxc_dsdv(ispin)%matrix)
309 4 : 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 8 : DO i = 1, 3
315 6 : CALL dbcsr_init_p(vcd_env%dipvel_ao(i)%matrix)
316 6 : CALL dbcsr_copy(vcd_env%dipvel_ao(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
317 :
318 6 : CALL dbcsr_init_p(vcd_env%dipvel_ao_delta(i)%matrix)
319 8 : 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 8 : DO i = 1, 3
324 6 : CALL dbcsr_init_p(vcd_env%hcom(i)%matrix)
325 : CALL dbcsr_create(vcd_env%hcom(i)%matrix, template=matrix_ks(1)%matrix, &
326 6 : matrix_type=dbcsr_type_antisymmetric)
327 6 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%hcom(i)%matrix, sab_orb)
328 :
329 6 : 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 6 : matrix_type=dbcsr_type_antisymmetric)
332 6 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_rxrv(i)%matrix, sab_orb)
333 :
334 26 : DO j = 1, 3
335 18 : CALL dbcsr_init_p(vcd_env%matrix_rcomr(i, j)%matrix)
336 18 : CALL dbcsr_copy(vcd_env%matrix_rcomr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
337 18 : CALL dbcsr_init_p(vcd_env%matrix_rrcom(i, j)%matrix)
338 18 : CALL dbcsr_copy(vcd_env%matrix_rrcom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
339 18 : CALL dbcsr_init_p(vcd_env%matrix_dcom(i, j)%matrix)
340 18 : CALL dbcsr_copy(vcd_env%matrix_dcom(i, j)%matrix, matrix_ks(1)%matrix)
341 18 : CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0._dp)
342 :
343 18 : CALL dbcsr_init_p(vcd_env%matrix_r_rxvr(i, j)%matrix)
344 18 : CALL dbcsr_copy(vcd_env%matrix_r_rxvr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
345 18 : CALL dbcsr_set(vcd_env%matrix_r_rxvr(i, j)%matrix, 0._dp)
346 :
347 18 : CALL dbcsr_init_p(vcd_env%matrix_rxvr_r(i, j)%matrix)
348 18 : CALL dbcsr_copy(vcd_env%matrix_rxvr_r(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
349 18 : CALL dbcsr_set(vcd_env%matrix_rxvr_r(i, j)%matrix, 0._dp)
350 :
351 18 : CALL dbcsr_init_p(vcd_env%matrix_r_doublecom(i, j)%matrix)
352 18 : CALL dbcsr_copy(vcd_env%matrix_r_doublecom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
353 24 : 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 4 : DO ispin = 1, nspins
359 10 : DO i = 1, 3
360 6 : CALL dbcsr_init_p(vcd_env%matrix_hr(ispin, i)%matrix)
361 6 : CALL dbcsr_copy(vcd_env%matrix_hr(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
362 :
363 6 : CALL dbcsr_init_p(vcd_env%matrix_rh(ispin, i)%matrix)
364 8 : 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 8 : DO i = 1, 3
370 6 : CALL dbcsr_init_p(vcd_env%matrix_drpnl(i)%matrix)
371 8 : 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 2 : 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 2 : direction_Or=.TRUE.)
380 : CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
381 2 : direction_Or=.TRUE., rc=[0._dp, 0._dp, 0._dp])
382 2 : CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
383 2 : CALL build_matrix_r_vhxc(vcd_env%matrix_hr, qs_env, [0._dp, 0._dp, 0._dp])
384 :
385 2 : 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 2 : direction_Or=.FALSE.)
389 : CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
390 2 : direction_Or=.FALSE., rc=[0._dp, 0._dp, 0._dp])
391 2 : CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
392 2 : 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 2 : 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 2 : 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 2 : dft_control%qs_control%eps_ppnl, particle_set, cell, .FALSE.)
403 :
404 : ! lin_mom
405 2 : 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 2 : 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 2 : 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 2 : 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 2 : 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 2 : "PRINT%PROGRAM_RUN_INFO")
429 :
430 2 : CALL timestop(handle)
431 :
432 6 : 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 2 : 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 2 : CALL timeset(routineN, handle)
450 :
451 : ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
452 2 : CALL dcdr_env_cleanup(qs_env, vcd_env%dcdr_env)
453 :
454 2 : DEALLOCATE (vcd_env%apt_el_nvpt)
455 2 : DEALLOCATE (vcd_env%apt_nuc_nvpt)
456 2 : DEALLOCATE (vcd_env%apt_total_nvpt)
457 2 : DEALLOCATE (vcd_env%aat_atom_nvpt)
458 2 : DEALLOCATE (vcd_env%aat_atom_mfp)
459 :
460 2 : CALL cp_fm_release(vcd_env%dCV)
461 2 : CALL cp_fm_release(vcd_env%dCV_prime)
462 2 : CALL cp_fm_release(vcd_env%op_dV)
463 2 : CALL cp_fm_release(vcd_env%op_dB)
464 :
465 2 : CALL cp_fm_release(vcd_env%dCB)
466 2 : 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 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der)
474 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_right)
475 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_left)
476 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_difdip2)
477 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdV)
478 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdB)
479 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hxc_dsdv)
480 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%hcom)
481 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rcomr)
482 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rrcom)
483 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dcom)
484 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hr)
485 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rh)
486 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_drpnl)
487 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao)
488 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao_delta)
489 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxrv)
490 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_rxvr)
491 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxvr_r)
492 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_doublecom)
493 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp_33)
494 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp2_33)
495 2 : CALL timestop(handle)
496 :
497 2 : 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 18 : 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 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
524 : TYPE(cp_fm_type), POINTER :: mo_coeff
525 : TYPE(cp_logger_type), POINTER :: logger
526 18 : 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 18 : file_exists = .FALSE.
531 :
532 18 : CALL timeset(routineN, handle)
533 :
534 18 : NULLIFY (mos, para_env, logger, print_key, vecbuffer)
535 18 : logger => cp_get_default_logger()
536 :
537 : iounit = cp_print_key_unit_nr(logger, linres_section, &
538 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
539 :
540 : CALL get_qs_env(qs_env=qs_env, &
541 : para_env=para_env, &
542 18 : mos=mos)
543 :
544 18 : nspins = SIZE(mos)
545 :
546 18 : rst_unit = -1
547 18 : IF (para_env%is_source()) THEN
548 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
549 9 : n_rep_val=n_rep_val)
550 :
551 9 : CALL XSTRING(tag, ia, ie)
552 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
553 9 : //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
554 :
555 9 : IF (n_rep_val > 0) THEN
556 0 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
557 0 : CALL xstring(filename, ia, ie)
558 0 : 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 9 : print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
562 : filename = cp_print_key_generate_filename(logger, print_key, &
563 9 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
564 : END IF
565 9 : INQUIRE (FILE=filename, exist=file_exists)
566 : !
567 : ! open file
568 9 : 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 0 : unit_number=rst_unit)
575 :
576 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
577 0 : "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">"
578 : ELSE
579 9 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
580 9 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
581 : END IF
582 : END IF
583 :
584 18 : CALL para_env%bcast(file_exists)
585 :
586 18 : IF (file_exists) THEN
587 :
588 0 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
589 0 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
590 :
591 0 : ALLOCATE (vecbuffer(nao, max_block))
592 : !
593 : ! read headers
594 0 : IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
595 0 : CALL para_env%bcast(iostat)
596 :
597 0 : CALL para_env%bcast(beta_tmp)
598 0 : CALL para_env%bcast(lambda_tmp)
599 0 : CALL para_env%bcast(nspins_tmp)
600 0 : 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 0 : IF (nspins_tmp /= nspins) THEN
605 0 : CPABORT("nspins not consistent")
606 : END IF
607 0 : 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 0 : IF (lambda_tmp /= lambda) CPABORT("lambda not consistent")
611 0 : IF (beta_tmp /= beta) CPABORT("beta not consistent")
612 : !
613 0 : DO ispin = 1, nspins
614 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
615 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
616 : !
617 0 : IF (rst_unit > 0) READ (rst_unit) nmo_tmp
618 0 : CALL para_env%bcast(nmo_tmp)
619 0 : IF (nmo_tmp /= nmo) CPABORT("nmo not consistent")
620 : !
621 : ! read the response
622 0 : DO i = 1, nmo, MAX(max_block, 1)
623 0 : i_block = MIN(max_block, nmo - i + 1)
624 0 : DO j = 1, i_block
625 0 : IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
626 : END DO
627 0 : CALL para_env%bcast(vecbuffer)
628 0 : CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
629 : END DO
630 : END DO
631 :
632 0 : IF (iostat /= 0) THEN
633 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
634 0 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
635 : END IF
636 :
637 0 : DEALLOCATE (vecbuffer)
638 :
639 : END IF
640 :
641 18 : IF (para_env%is_source()) THEN
642 9 : IF (file_exists) CALL close_file(unit_number=rst_unit)
643 : END IF
644 :
645 18 : CALL timestop(handle)
646 :
647 18 : 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 18 : 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 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
674 : TYPE(cp_fm_type), POINTER :: mo_coeff
675 : TYPE(cp_logger_type), POINTER :: logger
676 18 : 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 18 : NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
681 :
682 18 : CALL timeset(routineN, handle)
683 :
684 18 : logger => cp_get_default_logger()
685 :
686 18 : 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 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
692 :
693 : CALL get_qs_env(qs_env=qs_env, &
694 : mos=mos, &
695 18 : para_env=para_env)
696 :
697 18 : nspins = SIZE(mos)
698 :
699 18 : my_status = "REPLACE"
700 18 : my_pos = "REWIND"
701 18 : CALL XSTRING(tag, ia, ie)
702 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
703 18 : //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 18 : file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED")
707 :
708 : filename = cp_print_key_generate_filename(logger, print_key, &
709 18 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
710 :
711 18 : IF (iounit > 0) THEN
712 : WRITE (UNIT=iounit, FMT="(T2,A)") &
713 9 : "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 18 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
720 18 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
721 72 : ALLOCATE (vecbuffer(nao, max_block))
722 :
723 18 : IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
724 :
725 36 : DO ispin = 1, nspins
726 18 : CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
727 :
728 18 : IF (rst_unit > 0) WRITE (rst_unit) nmo
729 :
730 72 : DO i = 1, nmo, MAX(max_block, 1)
731 18 : i_block = MIN(max_block, nmo - i + 1)
732 18 : 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 108 : DO j = 1, i_block
737 306 : IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
738 : END DO
739 : END DO
740 : END DO
741 :
742 18 : DEALLOCATE (vecbuffer)
743 :
744 : CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
745 36 : "PRINT%RESTART")
746 : END IF
747 :
748 18 : CALL timestop(handle)
749 :
750 18 : 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 2 : 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 2 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_el_nvpt, apt_nuc_dcdr, &
770 2 : apt_nuc_nvpt, apt_total_dcdr, &
771 2 : apt_total_nvpt
772 2 : 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 2 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
780 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
781 : TYPE(section_vals_type), POINTER :: vcd_section
782 :
783 2 : CALL timeset(routineN, handle)
784 :
785 2 : NULLIFY (logger)
786 :
787 2 : logger => cp_get_default_logger()
788 2 : output_unit = cp_logger_get_default_io_unit(logger)
789 :
790 2 : vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
791 :
792 2 : NULLIFY (particle_set)
793 2 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
794 2 : natom = SIZE(particle_set)
795 2 : nsubset = SIZE(molecule_set)
796 :
797 2 : apt_el_dcdr => vcd_env%dcdr_env%apt_el_dcdr
798 2 : apt_nuc_dcdr => vcd_env%dcdr_env%apt_nuc_dcdr
799 2 : apt_total_dcdr => vcd_env%dcdr_env%apt_total_dcdr
800 2 : apt_subset_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_subset
801 2 : apt_center_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_center
802 :
803 2 : apt_el_nvpt => vcd_env%apt_el_nvpt
804 2 : apt_nuc_nvpt => vcd_env%apt_nuc_nvpt
805 2 : apt_total_nvpt => vcd_env%apt_total_nvpt
806 :
807 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
808 1 : 'APT | Write the final APT matrix per atom (Position perturbation)'
809 8 : DO l = 1, natom
810 6 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
811 3 : 'APT | Atom', l, ' - GAPT ', &
812 : (apt_total_dcdr(1, 1, l) &
813 : + apt_total_dcdr(2, 2, l) &
814 6 : + apt_total_dcdr(3, 3, l))/3._dp
815 26 : DO i = 1, 3
816 51 : 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 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
821 1 : 'NVP | Write the final APT matrix per atom (Velocity perturbation)'
822 8 : DO l = 1, natom
823 6 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
824 3 : 'NVP | Atom', l, ' - GAPT ', &
825 : (apt_total_nvpt(1, 1, l) &
826 : + apt_total_nvpt(2, 2, l) &
827 6 : + apt_total_nvpt(3, 3, l))/3._dp
828 26 : DO i = 1, 3
829 24 : IF (vcd_env%output_unit > 0) THEN
830 : WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
831 36 : "NVP | ", apt_total_nvpt(i, :, l)
832 : END IF
833 : END DO
834 : END DO
835 :
836 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
837 1 : 'NVP | Write the final AAT matrix per atom (Velocity perturbation)'
838 8 : DO l = 1, natom
839 6 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
840 3 : 'NVP | Atom', l
841 26 : DO i = 1, 3
842 24 : IF (vcd_env%output_unit > 0) THEN
843 : WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
844 9 : "NVP | ", vcd_env%aat_atom_nvpt(i, :, l)
845 : END IF
846 : END DO
847 : END DO
848 :
849 2 : IF (vcd_env%do_mfp) THEN
850 0 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
851 0 : 'MFP | Write the final AAT matrix per atom (Magnetic Field perturbation)'
852 0 : DO l = 1, natom
853 0 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
854 0 : 'MFP | Atom', l
855 0 : DO i = 1, 3
856 0 : IF (vcd_env%output_unit > 0) THEN
857 : WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
858 0 : "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 2 : CALL get_qs_env(qs_env, results=results)
866 2 : description = "[DIPOLE]"
867 2 : 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 2 : sum_rule_0 = 0._dp
871 2 : sum_rule_1 = 0._dp
872 2 : sum_rule_2 = 0._dp
873 2 : sum_rule_0_second = 0._dp
874 2 : sum_rule_2_second = 0._dp
875 2 : sum_rule_3_second = 0._dp
876 2 : sum_rule_3_mfp = 0._dp
877 : standard_deviation = 0._dp
878 2 : standard_deviation_sum = 0._dp
879 :
880 8 : DO alpha = 1, 3
881 26 : DO beta = 1, 3
882 : ! 0: sum_lambda apt(alpha, beta, lambda)
883 72 : DO lambda = 1, natom
884 : sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
885 54 : + apt_total_dcdr(alpha, beta, lambda)
886 : sum_rule_0_second(alpha, beta) = sum_rule_0_second(alpha, beta) &
887 72 : + apt_total_nvpt(alpha, beta, lambda)
888 : END DO
889 :
890 : ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
891 72 : DO gamma = 1, 3
892 : sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
893 72 : + 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 72 : DO lambda = 1, natom
898 234 : DO gamma = 1, 3
899 702 : 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 486 : *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 648 : *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 72 : DO lambda = 1, natom
914 : sum_rule_3_second(alpha, beta) = sum_rule_3_second(alpha, beta) &
915 72 : + 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 24 : IF (vcd_env%do_mfp) THEN
920 : ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
921 0 : DO lambda = 1, natom
922 : sum_rule_3_mfp(alpha, beta) = sum_rule_3_mfp(alpha, beta) &
923 0 : + 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 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") "APT | Position perturbation sum rules"
931 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") &
932 1 : "APT |", " Total APT", "Dipole", "R * APT", "AAT"
933 : standard_deviation_sum = 0._dp
934 8 : DO alpha = 1, 3
935 26 : DO beta = 1, 3
936 18 : 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 18 : - mean**2)
940 18 : standard_deviation_sum = standard_deviation_sum + standard_deviation
941 :
942 18 : 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 9 : "APT | ", &
945 9 : alpha, beta, &
946 9 : sum_rule_0(alpha, beta), &
947 9 : sum_rule_1(alpha, beta), &
948 9 : sum_rule_2(alpha, beta), &
949 9 : sum_rule_3_mfp(alpha, beta), &
950 24 : standard_deviation
951 : END DO
952 : END DO
953 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
954 :
955 2 : IF (vcd_env%output_unit > 0) THEN
956 1 : WRITE (vcd_env%output_unit, "(A)") "NVP | Velocity perturbation sum rules"
957 1 : 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 2 : standard_deviation_sum = 0._dp
961 8 : DO alpha = 1, 3
962 26 : DO beta = 1, 3
963 18 : 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 18 : - mean**2)
967 18 : standard_deviation_sum = standard_deviation_sum + standard_deviation
968 18 : 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 9 : "NVP | ", &
971 9 : alpha, &
972 9 : beta, &
973 9 : sum_rule_0_second(alpha, beta), &
974 9 : sum_rule_1(alpha, beta), &
975 9 : sum_rule_2_second(alpha, beta), &
976 9 : sum_rule_3_second(alpha, beta), &
977 24 : standard_deviation
978 : END DO
979 : END DO
980 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
981 :
982 2 : CALL timestop(handle)
983 2 : END SUBROUTINE vcd_print
984 :
985 : END MODULE qs_vcd_utils
|
