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