Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : 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) &
157 0 : CPABORT("User-defined reference point should be given explicitly")
158 : END IF
159 :
160 : CALL get_reference_point(rpoint=vcd_env%magnetic_origin, qs_env=qs_env, &
161 : reference=reference, &
162 2 : ref_point=ref_point)
163 :
164 : ! Get the spatial origin from the input
165 2 : CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN", i_val=reference)
166 2 : CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", explicit=explicit)
167 2 : IF (explicit) THEN
168 0 : CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", r_vals=ref_point)
169 : ELSE
170 2 : IF (reference == use_mom_ref_user) &
171 0 : CPABORT("User-defined reference point should be given explicitly")
172 : END IF
173 :
174 : CALL get_reference_point(rpoint=vcd_env%spatial_origin, qs_env=qs_env, &
175 : reference=reference, &
176 2 : ref_point=ref_point)
177 :
178 8 : IF (vcd_env%distributed_origin .AND. ANY(vcd_env%magnetic_origin /= vcd_env%spatial_origin)) THEN
179 0 : CPWARN("The magnetic and spatial origins don't match")
180 : ! This is fine for NVP but will give unphysical results for MFP.
181 : END IF
182 :
183 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
184 1 : 'The reference point is', vcd_env%dcdr_env%ref_point
185 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
186 1 : 'The magnetic origin is', vcd_env%magnetic_origin
187 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
188 1 : 'The velocity origin is', vcd_env%spatial_origin
189 :
190 8 : vcd_env%magnetic_origin_atom = vcd_env%magnetic_origin
191 8 : vcd_env%spatial_origin_atom = vcd_env%spatial_origin
192 :
193 : CALL get_qs_env(qs_env=qs_env, &
194 : ks_env=ks_env, &
195 : dft_control=dft_control, &
196 : sab_orb=sab_orb, &
197 : sab_all=sab_all, &
198 : sap_ppnl=sap_ppnl, &
199 : particle_set=particle_set, &
200 : matrix_ks=matrix_ks, &
201 : cell=cell, &
202 2 : qs_kind_set=qs_kind_set)
203 :
204 2 : natom = SIZE(particle_set)
205 2 : nspins = dft_control%nspins
206 :
207 6 : ALLOCATE (vcd_env%apt_el_nvpt(3, 3, natom))
208 4 : ALLOCATE (vcd_env%apt_nuc_nvpt(3, 3, natom))
209 4 : ALLOCATE (vcd_env%apt_total_nvpt(3, 3, natom))
210 4 : ALLOCATE (vcd_env%aat_atom_nvpt(3, 3, natom))
211 4 : ALLOCATE (vcd_env%aat_atom_mfp(3, 3, natom))
212 80 : vcd_env%apt_el_nvpt = 0._dp
213 80 : vcd_env%apt_nuc_nvpt = 0._dp
214 80 : vcd_env%apt_total_nvpt = 0._dp
215 80 : vcd_env%aat_atom_nvpt = 0._dp
216 80 : vcd_env%aat_atom_mfp = 0._dp
217 :
218 8 : ALLOCATE (vcd_env%dCV(nspins))
219 6 : ALLOCATE (vcd_env%dCV_prime(nspins))
220 6 : ALLOCATE (vcd_env%op_dV(nspins))
221 6 : ALLOCATE (vcd_env%op_dB(nspins))
222 4 : DO ispin = 1, nspins
223 2 : CALL cp_fm_create(vcd_env%dCV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
224 2 : CALL cp_fm_create(vcd_env%dCV_prime(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
225 2 : CALL cp_fm_create(vcd_env%op_dV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
226 4 : CALL cp_fm_create(vcd_env%op_dB(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
227 : END DO
228 :
229 8 : ALLOCATE (vcd_env%dCB(3))
230 8 : ALLOCATE (vcd_env%dCB_prime(3))
231 8 : DO i = 1, 3
232 6 : CALL cp_fm_create(vcd_env%dCB(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
233 8 : CALL cp_fm_create(vcd_env%dCB_prime(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
234 : END DO
235 :
236 : ! DBCSR matrices
237 2 : CALL dbcsr_allocate_matrix_set(vcd_env%moments_der, 9, 3)
238 2 : CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_right, 9, 3)
239 2 : CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_left, 9, 3)
240 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_difdip2, 3, 3)
241 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdV, 3)
242 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdB, 3)
243 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hxc_dsdv, nspins)
244 :
245 2 : CALL dbcsr_allocate_matrix_set(vcd_env%hcom, 3)
246 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rcomr, 3, 3)
247 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rrcom, 3, 3)
248 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dcom, 3, 3)
249 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hr, nspins, 3)
250 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rh, nspins, 3)
251 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_drpnl, 3)
252 2 : CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao, 3)
253 2 : CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao_delta, 3)
254 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxrv, 3)
255 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_rxvr, 3, 3)
256 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxvr_r, 3, 3)
257 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_doublecom, 3, 3)
258 :
259 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp_33, 3, 3)
260 2 : CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp2_33, 3, 3)
261 20 : DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
262 74 : DO idir = 1, 3 ! d/dx, d/dy, d/dz
263 54 : CALL dbcsr_init_p(vcd_env%moments_der(i, idir)%matrix)
264 54 : CALL dbcsr_init_p(vcd_env%moments_der_right(i, idir)%matrix)
265 54 : CALL dbcsr_init_p(vcd_env%moments_der_left(i, idir)%matrix)
266 :
267 : CALL dbcsr_create(vcd_env%moments_der(i, idir)%matrix, template=matrix_ks(1)%matrix, &
268 54 : matrix_type=dbcsr_type_antisymmetric)
269 54 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%moments_der(i, idir)%matrix, sab_orb)
270 54 : CALL dbcsr_set(vcd_env%moments_der(i, idir)%matrix, 0.0_dp)
271 :
272 : ! And the ones which will be multiplied by delta_(mu/nu)
273 : CALL dbcsr_copy(vcd_env%moments_der_right(i, idir)%matrix, &
274 54 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
275 : CALL dbcsr_copy(vcd_env%moments_der_left(i, idir)%matrix, &
276 72 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
277 : END DO
278 : END DO
279 :
280 8 : DO i = 1, 3
281 24 : DO j = 1, 3
282 18 : CALL dbcsr_init_p(vcd_env%matrix_difdip2(i, j)%matrix)
283 18 : CALL dbcsr_copy(vcd_env%matrix_difdip2(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
284 18 : CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0.0_dp)
285 :
286 18 : CALL dbcsr_init_p(vcd_env%matrix_nosym_temp_33(i, j)%matrix)
287 : CALL dbcsr_create(vcd_env%matrix_nosym_temp_33(i, j)%matrix, template=matrix_ks(1)%matrix, &
288 18 : matrix_type=dbcsr_type_no_symmetry)
289 18 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp_33(i, j)%matrix, sab_all)
290 18 : CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(i, j)%matrix, 0._dp)
291 :
292 18 : CALL dbcsr_init_p(vcd_env%matrix_nosym_temp2_33(i, j)%matrix)
293 : CALL dbcsr_create(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, template=matrix_ks(1)%matrix, &
294 18 : matrix_type=dbcsr_type_no_symmetry)
295 18 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, sab_all)
296 24 : CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, 0._dp)
297 :
298 : END DO
299 6 : CALL dbcsr_init_p(vcd_env%matrix_dSdV(i)%matrix)
300 6 : CALL dbcsr_copy(vcd_env%matrix_dSdV(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
301 6 : CALL dbcsr_init_p(vcd_env%matrix_dSdB(i)%matrix)
302 8 : CALL dbcsr_copy(vcd_env%matrix_dSdB(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
303 : END DO
304 :
305 4 : DO ispin = 1, nspins
306 2 : CALL dbcsr_init_p(vcd_env%matrix_hxc_dsdv(ispin)%matrix)
307 4 : CALL dbcsr_copy(vcd_env%matrix_hxc_dsdv(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
308 : END DO
309 :
310 : ! Things for op_dV
311 : ! lin_mom
312 8 : DO i = 1, 3
313 6 : CALL dbcsr_init_p(vcd_env%dipvel_ao(i)%matrix)
314 6 : CALL dbcsr_copy(vcd_env%dipvel_ao(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
315 :
316 6 : CALL dbcsr_init_p(vcd_env%dipvel_ao_delta(i)%matrix)
317 8 : CALL dbcsr_copy(vcd_env%dipvel_ao_delta(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
318 : END DO
319 :
320 : ! [V, r]
321 8 : DO i = 1, 3
322 6 : CALL dbcsr_init_p(vcd_env%hcom(i)%matrix)
323 : CALL dbcsr_create(vcd_env%hcom(i)%matrix, template=matrix_ks(1)%matrix, &
324 6 : matrix_type=dbcsr_type_antisymmetric)
325 6 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%hcom(i)%matrix, sab_orb)
326 :
327 6 : CALL dbcsr_init_p(vcd_env%matrix_rxrv(i)%matrix)
328 : CALL dbcsr_create(vcd_env%matrix_rxrv(i)%matrix, template=matrix_ks(1)%matrix, &
329 6 : matrix_type=dbcsr_type_antisymmetric)
330 6 : CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_rxrv(i)%matrix, sab_orb)
331 :
332 26 : DO j = 1, 3
333 18 : CALL dbcsr_init_p(vcd_env%matrix_rcomr(i, j)%matrix)
334 18 : CALL dbcsr_copy(vcd_env%matrix_rcomr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
335 18 : CALL dbcsr_init_p(vcd_env%matrix_rrcom(i, j)%matrix)
336 18 : CALL dbcsr_copy(vcd_env%matrix_rrcom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
337 18 : CALL dbcsr_init_p(vcd_env%matrix_dcom(i, j)%matrix)
338 18 : CALL dbcsr_copy(vcd_env%matrix_dcom(i, j)%matrix, matrix_ks(1)%matrix)
339 18 : CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0._dp)
340 :
341 18 : CALL dbcsr_init_p(vcd_env%matrix_r_rxvr(i, j)%matrix)
342 18 : CALL dbcsr_copy(vcd_env%matrix_r_rxvr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
343 18 : CALL dbcsr_set(vcd_env%matrix_r_rxvr(i, j)%matrix, 0._dp)
344 :
345 18 : CALL dbcsr_init_p(vcd_env%matrix_rxvr_r(i, j)%matrix)
346 18 : CALL dbcsr_copy(vcd_env%matrix_rxvr_r(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
347 18 : CALL dbcsr_set(vcd_env%matrix_rxvr_r(i, j)%matrix, 0._dp)
348 :
349 18 : CALL dbcsr_init_p(vcd_env%matrix_r_doublecom(i, j)%matrix)
350 18 : CALL dbcsr_copy(vcd_env%matrix_r_doublecom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
351 24 : CALL dbcsr_set(vcd_env%matrix_r_doublecom(i, j)%matrix, 0._dp)
352 : END DO
353 : END DO
354 :
355 : ! matrix_hr: nonsymmetric dbcsr matrix
356 4 : DO ispin = 1, nspins
357 10 : DO i = 1, 3
358 6 : CALL dbcsr_init_p(vcd_env%matrix_hr(ispin, i)%matrix)
359 6 : CALL dbcsr_copy(vcd_env%matrix_hr(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
360 :
361 6 : CALL dbcsr_init_p(vcd_env%matrix_rh(ispin, i)%matrix)
362 8 : CALL dbcsr_copy(vcd_env%matrix_rh(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
363 : END DO
364 : END DO
365 :
366 : ! drpnl for the operator
367 8 : DO i = 1, 3
368 6 : CALL dbcsr_init_p(vcd_env%matrix_drpnl(i)%matrix)
369 8 : CALL dbcsr_copy(vcd_env%matrix_drpnl(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
370 : END DO
371 :
372 : ! NVP matrices
373 : ! hr matrices
374 2 : my_matrix_hr_1d => vcd_env%matrix_hr(1, 1:3)
375 : CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
376 : dft_control%qs_control%eps_ppnl, cell, [0._dp, 0._dp, 0._dp], &
377 2 : direction_Or=.TRUE.)
378 : CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
379 2 : direction_Or=.TRUE., rc=[0._dp, 0._dp, 0._dp])
380 2 : CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
381 2 : CALL build_matrix_r_vhxc(vcd_env%matrix_hr, qs_env, [0._dp, 0._dp, 0._dp])
382 :
383 2 : my_matrix_hr_1d => vcd_env%matrix_rh(1, 1:3)
384 : CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
385 : dft_control%qs_control%eps_ppnl, cell, [0._dp, 0._dp, 0._dp], &
386 2 : direction_Or=.FALSE.)
387 : CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
388 2 : direction_Or=.FALSE., rc=[0._dp, 0._dp, 0._dp])
389 2 : CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
390 2 : CALL build_matrix_r_vhxc(vcd_env%matrix_rh, qs_env, [0._dp, 0._dp, 0._dp])
391 :
392 : ! commutator terms
393 : ! - [V, r]
394 : CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
395 2 : particle_set, cell=cell, matrix_rv=vcd_env%hcom)
396 : ! <[V, r] * r> and <r * [V, r]>
397 : CALL build_com_rpnl_r(vcd_env%matrix_rcomr, qs_kind_set, sab_all, sap_ppnl, &
398 2 : dft_control%qs_control%eps_ppnl, particle_set, cell, .TRUE.)
399 : CALL build_com_rpnl_r(vcd_env%matrix_rrcom, qs_kind_set, sab_all, sap_ppnl, &
400 2 : dft_control%qs_control%eps_ppnl, particle_set, cell, .FALSE.)
401 :
402 : ! lin_mom
403 2 : CALL build_lin_mom_matrix(qs_env, vcd_env%dipvel_ao)
404 :
405 : ! AAT
406 : ! The moments are set to zero and then recomputed in the routine.
407 : CALL build_local_moments_der_matrix(qs_env, moments_der=vcd_env%moments_der, &
408 2 : nmoments_der=2, nmoments=0, ref_point=[0._dp, 0._dp, 0._dp])
409 :
410 : ! PP terms
411 : CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
412 : particle_set, matrix_rxrv=vcd_env%matrix_rxrv, ref_point=[0._dp, 0._dp, 0._dp], &
413 2 : cell=cell)
414 :
415 : CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
416 : particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
417 2 : matrix_r_rxvr=vcd_env%matrix_r_rxvr)
418 :
419 : CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
420 : particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
421 2 : matrix_rxvr_r=vcd_env%matrix_rxvr_r)
422 :
423 : ! Done with NVP matrices
424 :
425 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
426 2 : "PRINT%PROGRAM_RUN_INFO")
427 :
428 2 : CALL timestop(handle)
429 :
430 6 : END SUBROUTINE vcd_env_init
431 :
432 : ! *****************************************************************************
433 : !> \brief Deallocate the vcd environment
434 : !> \param qs_env ...
435 : !> \param vcd_env ...
436 : !> \author Edward Ditler
437 : ! **************************************************************************************************
438 2 : SUBROUTINE vcd_env_cleanup(qs_env, vcd_env)
439 :
440 : TYPE(qs_environment_type), POINTER :: qs_env
441 : TYPE(vcd_env_type) :: vcd_env
442 :
443 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vcd_env_cleanup'
444 :
445 : INTEGER :: handle
446 :
447 2 : CALL timeset(routineN, handle)
448 :
449 : ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
450 2 : CALL dcdr_env_cleanup(qs_env, vcd_env%dcdr_env)
451 :
452 2 : DEALLOCATE (vcd_env%apt_el_nvpt)
453 2 : DEALLOCATE (vcd_env%apt_nuc_nvpt)
454 2 : DEALLOCATE (vcd_env%apt_total_nvpt)
455 2 : DEALLOCATE (vcd_env%aat_atom_nvpt)
456 2 : DEALLOCATE (vcd_env%aat_atom_mfp)
457 :
458 2 : CALL cp_fm_release(vcd_env%dCV)
459 2 : CALL cp_fm_release(vcd_env%dCV_prime)
460 2 : CALL cp_fm_release(vcd_env%op_dV)
461 2 : CALL cp_fm_release(vcd_env%op_dB)
462 :
463 2 : CALL cp_fm_release(vcd_env%dCB)
464 2 : CALL cp_fm_release(vcd_env%dCB_prime)
465 :
466 : ! DBCSR matrices
467 : ! Probably, the memory requirements could be reduced by quite a bit
468 : ! by not storing each term in its own set of matrices.
469 : ! On the other hand, the memory bottleneck is usually the numerical
470 : ! integration grid.
471 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der)
472 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_right)
473 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_left)
474 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_difdip2)
475 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdV)
476 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdB)
477 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hxc_dsdv)
478 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%hcom)
479 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rcomr)
480 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rrcom)
481 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dcom)
482 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hr)
483 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rh)
484 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_drpnl)
485 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao)
486 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao_delta)
487 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxrv)
488 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_rxvr)
489 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxvr_r)
490 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_doublecom)
491 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp_33)
492 2 : CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp2_33)
493 2 : CALL timestop(handle)
494 :
495 2 : END SUBROUTINE vcd_env_cleanup
496 :
497 : ! **************************************************************************************************
498 : !> \brief Copied from linres_read_restart
499 : !> \param qs_env ...
500 : !> \param linres_section ...
501 : !> \param vec ...
502 : !> \param lambda ...
503 : !> \param beta ...
504 : !> \param tag ...
505 : !> \author Edward Ditler
506 : ! **************************************************************************************************
507 18 : SUBROUTINE vcd_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
508 : TYPE(qs_environment_type), POINTER :: qs_env
509 : TYPE(section_vals_type), POINTER :: linres_section
510 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
511 : INTEGER, INTENT(IN) :: lambda, beta
512 : CHARACTER(LEN=*) :: tag
513 :
514 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vcd_read_restart'
515 :
516 : CHARACTER(LEN=default_path_length) :: filename
517 : CHARACTER(LEN=default_string_length) :: my_middle
518 : INTEGER :: beta_tmp, handle, i, i_block, ia, ie, iostat, iounit, ispin, j, lambda_tmp, &
519 : max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit
520 : LOGICAL :: file_exists
521 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
522 : TYPE(cp_fm_type), POINTER :: mo_coeff
523 : TYPE(cp_logger_type), POINTER :: logger
524 18 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
525 : TYPE(mp_para_env_type), POINTER :: para_env
526 : TYPE(section_vals_type), POINTER :: print_key
527 :
528 18 : file_exists = .FALSE.
529 :
530 18 : CALL timeset(routineN, handle)
531 :
532 18 : NULLIFY (mos, para_env, logger, print_key, vecbuffer)
533 18 : logger => cp_get_default_logger()
534 :
535 : iounit = cp_print_key_unit_nr(logger, linres_section, &
536 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
537 :
538 : CALL get_qs_env(qs_env=qs_env, &
539 : para_env=para_env, &
540 18 : mos=mos)
541 :
542 18 : nspins = SIZE(mos)
543 :
544 18 : rst_unit = -1
545 18 : IF (para_env%is_source()) THEN
546 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
547 9 : n_rep_val=n_rep_val)
548 :
549 9 : CALL XSTRING(tag, ia, ie)
550 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
551 9 : //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
552 :
553 9 : IF (n_rep_val > 0) THEN
554 0 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
555 0 : CALL xstring(filename, ia, ie)
556 0 : filename = filename(ia:ie)//TRIM(my_middle)//".lr"
557 : ELSE
558 : ! try to read from the filename that is generated automatically from the printkey
559 9 : print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
560 : filename = cp_print_key_generate_filename(logger, print_key, &
561 9 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
562 : END IF
563 9 : INQUIRE (FILE=filename, exist=file_exists)
564 : !
565 : ! open file
566 9 : IF (file_exists) THEN
567 : CALL open_file(file_name=TRIM(filename), &
568 : file_action="READ", &
569 : file_form="UNFORMATTED", &
570 : file_position="REWIND", &
571 : file_status="OLD", &
572 0 : unit_number=rst_unit)
573 :
574 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
575 0 : "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">"
576 : ELSE
577 9 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
578 9 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
579 : END IF
580 : END IF
581 :
582 18 : CALL para_env%bcast(file_exists)
583 :
584 18 : IF (file_exists) THEN
585 :
586 0 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
587 0 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
588 :
589 0 : ALLOCATE (vecbuffer(nao, max_block))
590 : !
591 : ! read headers
592 0 : IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
593 0 : CALL para_env%bcast(iostat)
594 :
595 0 : CALL para_env%bcast(beta_tmp)
596 0 : CALL para_env%bcast(lambda_tmp)
597 0 : CALL para_env%bcast(nspins_tmp)
598 0 : CALL para_env%bcast(nao_tmp)
599 :
600 : ! check that the number nao, nmo and nspins are
601 : ! the same as in the current mos
602 0 : IF (nspins_tmp .NE. nspins) THEN
603 0 : CPABORT("nspins not consistent")
604 : END IF
605 0 : IF (nao_tmp .NE. nao) CPABORT("nao not consistent")
606 : ! check that it's the right file
607 : ! the same as in the current mos
608 0 : IF (lambda_tmp .NE. lambda) CPABORT("lambda not consistent")
609 0 : IF (beta_tmp .NE. beta) CPABORT("beta not consistent")
610 : !
611 0 : DO ispin = 1, nspins
612 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
613 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
614 : !
615 0 : IF (rst_unit > 0) READ (rst_unit) nmo_tmp
616 0 : CALL para_env%bcast(nmo_tmp)
617 0 : IF (nmo_tmp .NE. nmo) CPABORT("nmo not consistent")
618 : !
619 : ! read the response
620 0 : DO i = 1, nmo, MAX(max_block, 1)
621 0 : i_block = MIN(max_block, nmo - i + 1)
622 0 : DO j = 1, i_block
623 0 : IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
624 : END DO
625 0 : CALL para_env%bcast(vecbuffer)
626 0 : CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
627 : END DO
628 : END DO
629 :
630 0 : IF (iostat /= 0) THEN
631 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
632 0 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
633 : END IF
634 :
635 0 : DEALLOCATE (vecbuffer)
636 :
637 : END IF
638 :
639 18 : IF (para_env%is_source()) THEN
640 9 : IF (file_exists) CALL close_file(unit_number=rst_unit)
641 : END IF
642 :
643 18 : CALL timestop(handle)
644 :
645 18 : END SUBROUTINE vcd_read_restart
646 :
647 : ! **************************************************************************************************
648 : !> \brief Copied from linres_write_restart
649 : !> \param qs_env ...
650 : !> \param linres_section ...
651 : !> \param vec ...
652 : !> \param lambda ...
653 : !> \param beta ...
654 : !> \param tag ...
655 : !> \author Edward Ditler
656 : ! **************************************************************************************************
657 18 : SUBROUTINE vcd_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
658 : TYPE(qs_environment_type), POINTER :: qs_env
659 : TYPE(section_vals_type), POINTER :: linres_section
660 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
661 : INTEGER, INTENT(IN) :: lambda, beta
662 : CHARACTER(LEN=*) :: tag
663 :
664 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vcd_write_restart'
665 :
666 : CHARACTER(LEN=default_path_length) :: filename
667 : CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status
668 : INTEGER :: handle, i, i_block, ia, ie, iounit, &
669 : ispin, j, max_block, nao, nmo, nspins, &
670 : rst_unit
671 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
672 : TYPE(cp_fm_type), POINTER :: mo_coeff
673 : TYPE(cp_logger_type), POINTER :: logger
674 18 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
675 : TYPE(mp_para_env_type), POINTER :: para_env
676 : TYPE(section_vals_type), POINTER :: print_key
677 :
678 18 : NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
679 :
680 18 : CALL timeset(routineN, handle)
681 :
682 18 : logger => cp_get_default_logger()
683 :
684 18 : IF (BTEST(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
685 : used_print_key=print_key), &
686 : cp_p_file)) THEN
687 :
688 : iounit = cp_print_key_unit_nr(logger, linres_section, &
689 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
690 :
691 : CALL get_qs_env(qs_env=qs_env, &
692 : mos=mos, &
693 18 : para_env=para_env)
694 :
695 18 : nspins = SIZE(mos)
696 :
697 18 : my_status = "REPLACE"
698 18 : my_pos = "REWIND"
699 18 : CALL XSTRING(tag, ia, ie)
700 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
701 18 : //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
702 : rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
703 : extension=".lr", middle_name=TRIM(my_middle), file_status=TRIM(my_status), &
704 18 : file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED")
705 :
706 : filename = cp_print_key_generate_filename(logger, print_key, &
707 18 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
708 :
709 18 : IF (iounit > 0) THEN
710 : WRITE (UNIT=iounit, FMT="(T2,A)") &
711 9 : "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">"
712 : END IF
713 :
714 : !
715 : ! write data to file
716 : ! use the scalapack block size as a default for buffering columns
717 18 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
718 18 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
719 72 : ALLOCATE (vecbuffer(nao, max_block))
720 :
721 18 : IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
722 :
723 36 : DO ispin = 1, nspins
724 18 : CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
725 :
726 18 : IF (rst_unit > 0) WRITE (rst_unit) nmo
727 :
728 72 : DO i = 1, nmo, MAX(max_block, 1)
729 18 : i_block = MIN(max_block, nmo - i + 1)
730 18 : CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
731 : ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
732 : ! to old ones, and in cases where max_block is different between runs, as might happen during
733 : ! restarts with a different number of CPUs
734 108 : DO j = 1, i_block
735 90 : IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
736 : END DO
737 : END DO
738 : END DO
739 :
740 18 : DEALLOCATE (vecbuffer)
741 :
742 : CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
743 36 : "PRINT%RESTART")
744 : END IF
745 :
746 18 : CALL timestop(handle)
747 :
748 18 : END SUBROUTINE vcd_write_restart
749 :
750 : ! **************************************************************************************************
751 : !> \brief Print the APTs, AATs, and sum rules
752 : !> \param vcd_env ...
753 : !> \param qs_env ...
754 : !> \author Edward Ditler
755 : ! **************************************************************************************************
756 2 : SUBROUTINE vcd_print(vcd_env, qs_env)
757 : TYPE(vcd_env_type) :: vcd_env
758 : TYPE(qs_environment_type), POINTER :: qs_env
759 :
760 : CHARACTER(len=*), PARAMETER :: routineN = 'vcd_print'
761 :
762 : CHARACTER(LEN=default_string_length) :: description
763 : INTEGER :: alpha, beta, delta, gamma, handle, i, l, &
764 : lambda, natom, nsubset, output_unit
765 : REAL(dp) :: mean, standard_deviation, &
766 : standard_deviation_sum
767 2 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_el_nvpt, apt_nuc_dcdr, &
768 2 : apt_nuc_nvpt, apt_total_dcdr, &
769 2 : apt_total_nvpt
770 2 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center_dcdr, apt_subset_dcdr
771 : REAL(kind=dp), DIMENSION(3, 3) :: sum_rule_0, sum_rule_0_second, &
772 : sum_rule_1, sum_rule_2, &
773 : sum_rule_2_second, sum_rule_3_mfp, &
774 : sum_rule_3_second
775 : TYPE(cp_logger_type), POINTER :: logger
776 : TYPE(cp_result_type), POINTER :: results
777 2 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
778 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
779 : TYPE(section_vals_type), POINTER :: vcd_section
780 :
781 2 : CALL timeset(routineN, handle)
782 :
783 2 : NULLIFY (logger)
784 :
785 2 : logger => cp_get_default_logger()
786 2 : output_unit = cp_logger_get_default_io_unit(logger)
787 :
788 2 : vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
789 :
790 2 : NULLIFY (particle_set)
791 2 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
792 2 : natom = SIZE(particle_set)
793 2 : nsubset = SIZE(molecule_set)
794 :
795 2 : apt_el_dcdr => vcd_env%dcdr_env%apt_el_dcdr
796 2 : apt_nuc_dcdr => vcd_env%dcdr_env%apt_nuc_dcdr
797 2 : apt_total_dcdr => vcd_env%dcdr_env%apt_total_dcdr
798 2 : apt_subset_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_subset
799 2 : apt_center_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_center
800 :
801 2 : apt_el_nvpt => vcd_env%apt_el_nvpt
802 2 : apt_nuc_nvpt => vcd_env%apt_nuc_nvpt
803 2 : apt_total_nvpt => vcd_env%apt_total_nvpt
804 :
805 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
806 1 : 'APT | Write the final APT matrix per atom (Position perturbation)'
807 8 : DO l = 1, natom
808 6 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
809 3 : 'APT | Atom', l, ' - GAPT ', &
810 : (apt_total_dcdr(1, 1, l) &
811 : + apt_total_dcdr(2, 2, l) &
812 6 : + apt_total_dcdr(3, 3, l))/3._dp
813 26 : DO i = 1, 3
814 24 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l)
815 : END DO
816 : END DO
817 :
818 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
819 1 : 'NVP | Write the final APT matrix per atom (Velocity perturbation)'
820 8 : DO l = 1, natom
821 6 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
822 3 : 'NVP | Atom', l, ' - GAPT ', &
823 : (apt_total_nvpt(1, 1, l) &
824 : + apt_total_nvpt(2, 2, l) &
825 6 : + apt_total_nvpt(3, 3, l))/3._dp
826 26 : DO i = 1, 3
827 18 : IF (vcd_env%output_unit > 0) &
828 : WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
829 15 : "NVP | ", apt_total_nvpt(i, :, l)
830 : END DO
831 : END DO
832 :
833 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
834 1 : 'NVP | Write the final AAT matrix per atom (Velocity perturbation)'
835 8 : DO l = 1, natom
836 6 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
837 3 : 'NVP | Atom', l
838 26 : DO i = 1, 3
839 18 : IF (vcd_env%output_unit > 0) &
840 : WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
841 15 : "NVP | ", vcd_env%aat_atom_nvpt(i, :, l)
842 : END DO
843 : END DO
844 :
845 2 : IF (vcd_env%do_mfp) THEN
846 0 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
847 0 : 'MFP | Write the final AAT matrix per atom (Magnetic Field perturbation)'
848 0 : DO l = 1, natom
849 0 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
850 0 : 'MFP | Atom', l
851 0 : DO i = 1, 3
852 0 : IF (vcd_env%output_unit > 0) &
853 : WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
854 0 : "MFP | ", vcd_env%aat_atom_mfp(i, :, l)
855 : END DO
856 : END DO
857 : END IF
858 :
859 : ! Get the dipole
860 2 : CALL get_qs_env(qs_env, results=results)
861 2 : description = "[DIPOLE]"
862 2 : CALL get_results(results=results, description=description, values=vcd_env%dcdr_env%dipole_pos(1:3))
863 :
864 : ! Sum rules [for all alpha, beta]
865 2 : sum_rule_0 = 0._dp
866 2 : sum_rule_1 = 0._dp
867 2 : sum_rule_2 = 0._dp
868 2 : sum_rule_0_second = 0._dp
869 2 : sum_rule_2_second = 0._dp
870 2 : sum_rule_3_second = 0._dp
871 2 : sum_rule_3_mfp = 0._dp
872 2 : standard_deviation = 0._dp
873 2 : standard_deviation_sum = 0._dp
874 :
875 8 : DO alpha = 1, 3
876 26 : DO beta = 1, 3
877 : ! 0: sum_lambda apt(alpha, beta, lambda)
878 72 : DO lambda = 1, natom
879 : sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
880 54 : + apt_total_dcdr(alpha, beta, lambda)
881 : sum_rule_0_second(alpha, beta) = sum_rule_0_second(alpha, beta) &
882 72 : + apt_total_nvpt(alpha, beta, lambda)
883 : END DO
884 :
885 : ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
886 72 : DO gamma = 1, 3
887 : sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
888 72 : + Levi_Civita(alpha, beta, gamma)*vcd_env%dcdr_env%dipole_pos(gamma)
889 : END DO
890 :
891 : ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda)
892 72 : DO lambda = 1, natom
893 234 : DO gamma = 1, 3
894 702 : DO delta = 1, 3
895 : sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) &
896 : + Levi_Civita(beta, gamma, delta) &
897 : *particle_set(lambda)%r(gamma) &
898 486 : *apt_total_dcdr(delta, alpha, lambda)
899 : sum_rule_2_second(alpha, beta) = sum_rule_2_second(alpha, beta) &
900 : + Levi_Civita(beta, gamma, delta) &
901 : *particle_set(lambda)%r(gamma) &
902 648 : *apt_total_nvpt(delta, alpha, lambda)
903 : END DO
904 : END DO
905 : END DO
906 :
907 : ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
908 72 : DO lambda = 1, natom
909 : sum_rule_3_second(alpha, beta) = sum_rule_3_second(alpha, beta) &
910 72 : + vcd_env%aat_atom_nvpt(alpha, beta, lambda)
911 : ! + 2._dp*c_light_au*vcd_env%aat_atom_nvpt(alpha, beta, lambda)
912 : END DO
913 :
914 24 : IF (vcd_env%do_mfp) THEN
915 : ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
916 0 : DO lambda = 1, natom
917 : sum_rule_3_mfp(alpha, beta) = sum_rule_3_mfp(alpha, beta) &
918 0 : + vcd_env%aat_atom_mfp(alpha, beta, lambda)
919 : END DO
920 : END IF
921 :
922 : END DO ! beta
923 : END DO ! alpha
924 :
925 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") "APT | Position perturbation sum rules"
926 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") &
927 1 : "APT |", " Total APT", "Dipole", "R * APT", "AAT"
928 : standard_deviation_sum = 0._dp
929 8 : DO alpha = 1, 3
930 26 : DO beta = 1, 3
931 18 : mean = (sum_rule_1(alpha, beta) + sum_rule_2(alpha, beta) + sum_rule_3_mfp(alpha, beta))/3
932 : standard_deviation = &
933 : SQRT((sum_rule_1(alpha, beta)**2 + sum_rule_2(alpha, beta)**2 + sum_rule_3_mfp(alpha, beta)**2)/3 &
934 18 : - mean**2)
935 18 : standard_deviation_sum = standard_deviation_sum + standard_deviation
936 :
937 18 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, &
938 : "(A,I3,I3,F15.6,F15.6,F15.6,F15.6,F15.6)") &
939 9 : "APT | ", &
940 9 : alpha, beta, &
941 9 : sum_rule_0(alpha, beta), &
942 9 : sum_rule_1(alpha, beta), &
943 9 : sum_rule_2(alpha, beta), &
944 9 : sum_rule_3_mfp(alpha, beta), &
945 24 : standard_deviation
946 : END DO
947 : END DO
948 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
949 :
950 2 : IF (vcd_env%output_unit > 0) THEN
951 1 : WRITE (vcd_env%output_unit, "(A)") "NVP | Velocity perturbation sum rules"
952 1 : WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") "NVP |", " Total APT", "Dipole", "R * APT", "AAT"
953 : END IF
954 :
955 2 : standard_deviation_sum = 0._dp
956 8 : DO alpha = 1, 3
957 26 : DO beta = 1, 3
958 18 : mean = (sum_rule_1(alpha, beta) + sum_rule_2_second(alpha, beta) + sum_rule_3_second(alpha, beta))/3
959 : standard_deviation = &
960 : SQRT((sum_rule_1(alpha, beta)**2 + sum_rule_2_second(alpha, beta)**2 + sum_rule_3_second(alpha, beta)**2)/3 &
961 18 : - mean**2)
962 18 : standard_deviation_sum = standard_deviation_sum + standard_deviation
963 18 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, &
964 : "(A,I3,I3,F15.6,F15.6,F15.6,F15.6,F15.6)") &
965 9 : "NVP | ", &
966 9 : alpha, &
967 9 : beta, &
968 9 : sum_rule_0_second(alpha, beta), &
969 9 : sum_rule_1(alpha, beta), &
970 9 : sum_rule_2_second(alpha, beta), &
971 9 : sum_rule_3_second(alpha, beta), &
972 24 : standard_deviation
973 : END DO
974 : END DO
975 2 : IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
976 :
977 2 : CALL timestop(handle)
978 2 : END SUBROUTINE vcd_print
979 :
980 : END MODULE qs_vcd_utils
|
