(git:3add494)
qs_vcd_utils.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
13  USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
14  dbcsr_deallocate_matrix_set
15  USE cp_files, ONLY: close_file,&
16  open_file
17  USE cp_fm_types, ONLY: cp_fm_create,&
18  cp_fm_get_info,&
19  cp_fm_get_submatrix,&
20  cp_fm_release,&
21  cp_fm_set_submatrix,&
22  cp_fm_type
23  USE cp_log_handling, ONLY: cp_get_default_logger,&
24  cp_logger_get_default_io_unit,&
25  cp_logger_type,&
26  cp_to_string
27  USE cp_output_handling, ONLY: cp_p_file,&
28  cp_print_key_finished_output,&
29  cp_print_key_generate_filename,&
30  cp_print_key_should_output,&
31  cp_print_key_unit_nr
32  USE cp_result_methods, ONLY: get_results
33  USE cp_result_types, ONLY: cp_result_type
34  USE dbcsr_api, ONLY: dbcsr_copy,&
35  dbcsr_create,&
36  dbcsr_init_p,&
37  dbcsr_p_type,&
38  dbcsr_set,&
39  dbcsr_type_antisymmetric,&
40  dbcsr_type_no_symmetry
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  SUBROUTINE vcd_env_init(vcd_env, qs_env)
96  TYPE(vcd_env_type), TARGET :: vcd_env
97  TYPE(qs_environment_type), POINTER :: qs_env
98 
99  CHARACTER(LEN=*), PARAMETER :: routinen = 'vcd_env_init'
100 
101  INTEGER :: handle, i, idir, ispin, j, natom, &
102  nspins, output_unit, reference, &
103  unit_number
104  LOGICAL :: explicit
105  REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
106  TYPE(cell_type), POINTER :: cell
107  TYPE(cp_logger_type), POINTER :: logger
108  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, my_matrix_hr_1d
109  TYPE(dft_control_type), POINTER :: dft_control
110  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
111  POINTER :: sab_all, sab_orb, sap_ppnl
112  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
113  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
114  TYPE(qs_ks_env_type), POINTER :: ks_env
115  TYPE(section_vals_type), POINTER :: lr_section, vcd_section
116 
117  CALL timeset(routinen, handle)
118  vcd_env%do_mfp = .false.
119 
120  ! Set up the logger
121  NULLIFY (logger, vcd_section, lr_section)
122  logger => cp_get_default_logger()
123  vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
124  vcd_env%output_unit = cp_print_key_unit_nr(logger, vcd_section, "PRINT%VCD", &
125  extension=".data", middle_name="vcd", log_filename=.false., &
126  file_position="REWIND", file_status="REPLACE")
127 
128  lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
129  output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
130  extension=".linresLog")
131  unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
132 
133  ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
134  CALL dcdr_env_init(vcd_env%dcdr_env, qs_env)
135  ! vcd_env%dcdr_env%output_unit = vcd_env%output_unit
136 
137  IF (output_unit > 0) THEN
138  WRITE (output_unit, "(/,T20,A,/)") "*** Start NVPT/MFPT calculation ***"
139  END IF
140 
141  ! Just to make sure. The memory requirements are tiny.
142  CALL init_orbital_pointers(12)
143 
144  CALL section_vals_val_get(vcd_section, "DISTRIBUTED_ORIGIN", l_val=vcd_env%distributed_origin)
145  CALL section_vals_val_get(vcd_section, "ORIGIN_DEPENDENT_MFP", l_val=vcd_env%origin_dependent_op_mfp)
146 
147  ! Reference point
148  vcd_env%magnetic_origin = 0._dp
149  vcd_env%spatial_origin = 0._dp
150  ! Get the magnetic origin from the input
151  CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN", i_val=reference)
152  CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", explicit=explicit)
153  IF (explicit) THEN
154  CALL section_vals_val_get(vcd_section, "MAGNETIC_ORIGIN_REFERENCE", r_vals=ref_point)
155  ELSE
156  IF (reference == use_mom_ref_user) &
157  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  ref_point=ref_point)
163 
164  ! Get the spatial origin from the input
165  CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN", i_val=reference)
166  CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", explicit=explicit)
167  IF (explicit) THEN
168  CALL section_vals_val_get(vcd_section, "SPATIAL_ORIGIN_REFERENCE", r_vals=ref_point)
169  ELSE
170  IF (reference == use_mom_ref_user) &
171  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  ref_point=ref_point)
177 
178  IF (vcd_env%distributed_origin .AND. any(vcd_env%magnetic_origin /= vcd_env%spatial_origin)) THEN
179  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  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
184  'The reference point is', vcd_env%dcdr_env%ref_point
185  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
186  'The magnetic origin is', vcd_env%magnetic_origin
187  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,3F10.6)") &
188  'The velocity origin is', vcd_env%spatial_origin
189 
190  vcd_env%magnetic_origin_atom = vcd_env%magnetic_origin
191  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  qs_kind_set=qs_kind_set)
203 
204  natom = SIZE(particle_set)
205  nspins = dft_control%nspins
206 
207  ALLOCATE (vcd_env%apt_el_nvpt(3, 3, natom))
208  ALLOCATE (vcd_env%apt_nuc_nvpt(3, 3, natom))
209  ALLOCATE (vcd_env%apt_total_nvpt(3, 3, natom))
210  ALLOCATE (vcd_env%aat_atom_nvpt(3, 3, natom))
211  ALLOCATE (vcd_env%aat_atom_mfp(3, 3, natom))
212  vcd_env%apt_el_nvpt = 0._dp
213  vcd_env%apt_nuc_nvpt = 0._dp
214  vcd_env%apt_total_nvpt = 0._dp
215  vcd_env%aat_atom_nvpt = 0._dp
216  vcd_env%aat_atom_mfp = 0._dp
217 
218  ALLOCATE (vcd_env%dCV(nspins))
219  ALLOCATE (vcd_env%dCV_prime(nspins))
220  ALLOCATE (vcd_env%op_dV(nspins))
221  ALLOCATE (vcd_env%op_dB(nspins))
222  DO ispin = 1, nspins
223  CALL cp_fm_create(vcd_env%dCV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
224  CALL cp_fm_create(vcd_env%dCV_prime(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
225  CALL cp_fm_create(vcd_env%op_dV(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
226  CALL cp_fm_create(vcd_env%op_dB(ispin), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
227  END DO
228 
229  ALLOCATE (vcd_env%dCB(3))
230  ALLOCATE (vcd_env%dCB_prime(3))
231  DO i = 1, 3
232  CALL cp_fm_create(vcd_env%dCB(i), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
233  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  CALL dbcsr_allocate_matrix_set(vcd_env%moments_der, 9, 3)
238  CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_right, 9, 3)
239  CALL dbcsr_allocate_matrix_set(vcd_env%moments_der_left, 9, 3)
240  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_difdip2, 3, 3)
241  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdV, 3)
242  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dSdB, 3)
243  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hxc_dsdv, nspins)
244 
245  CALL dbcsr_allocate_matrix_set(vcd_env%hcom, 3)
246  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rcomr, 3, 3)
247  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rrcom, 3, 3)
248  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_dcom, 3, 3)
249  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_hr, nspins, 3)
250  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rh, nspins, 3)
251  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_drpnl, 3)
252  CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao, 3)
253  CALL dbcsr_allocate_matrix_set(vcd_env%dipvel_ao_delta, 3)
254  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxrv, 3)
255  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_rxvr, 3, 3)
256  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_rxvr_r, 3, 3)
257  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_r_doublecom, 3, 3)
258 
259  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp_33, 3, 3)
260  CALL dbcsr_allocate_matrix_set(vcd_env%matrix_nosym_temp2_33, 3, 3)
261  DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
262  DO idir = 1, 3 ! d/dx, d/dy, d/dz
263  CALL dbcsr_init_p(vcd_env%moments_der(i, idir)%matrix)
264  CALL dbcsr_init_p(vcd_env%moments_der_right(i, idir)%matrix)
265  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  matrix_type=dbcsr_type_antisymmetric)
269  CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%moments_der(i, idir)%matrix, sab_orb)
270  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  vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
275  CALL dbcsr_copy(vcd_env%moments_der_left(i, idir)%matrix, &
276  vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
277  END DO
278  END DO
279 
280  DO i = 1, 3
281  DO j = 1, 3
282  CALL dbcsr_init_p(vcd_env%matrix_difdip2(i, j)%matrix)
283  CALL dbcsr_copy(vcd_env%matrix_difdip2(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
284  CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0.0_dp)
285 
286  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  matrix_type=dbcsr_type_no_symmetry)
289  CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp_33(i, j)%matrix, sab_all)
290  CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(i, j)%matrix, 0._dp)
291 
292  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  matrix_type=dbcsr_type_no_symmetry)
295  CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, sab_all)
296  CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(i, j)%matrix, 0._dp)
297 
298  END DO
299  CALL dbcsr_init_p(vcd_env%matrix_dSdV(i)%matrix)
300  CALL dbcsr_copy(vcd_env%matrix_dSdV(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
301  CALL dbcsr_init_p(vcd_env%matrix_dSdB(i)%matrix)
302  CALL dbcsr_copy(vcd_env%matrix_dSdB(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
303  END DO
304 
305  DO ispin = 1, nspins
306  CALL dbcsr_init_p(vcd_env%matrix_hxc_dsdv(ispin)%matrix)
307  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  DO i = 1, 3
313  CALL dbcsr_init_p(vcd_env%dipvel_ao(i)%matrix)
314  CALL dbcsr_copy(vcd_env%dipvel_ao(i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
315 
316  CALL dbcsr_init_p(vcd_env%dipvel_ao_delta(i)%matrix)
317  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  DO i = 1, 3
322  CALL dbcsr_init_p(vcd_env%hcom(i)%matrix)
323  CALL dbcsr_create(vcd_env%hcom(i)%matrix, template=matrix_ks(1)%matrix, &
324  matrix_type=dbcsr_type_antisymmetric)
325  CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%hcom(i)%matrix, sab_orb)
326 
327  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  matrix_type=dbcsr_type_antisymmetric)
330  CALL cp_dbcsr_alloc_block_from_nbl(vcd_env%matrix_rxrv(i)%matrix, sab_orb)
331 
332  DO j = 1, 3
333  CALL dbcsr_init_p(vcd_env%matrix_rcomr(i, j)%matrix)
334  CALL dbcsr_copy(vcd_env%matrix_rcomr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
335  CALL dbcsr_init_p(vcd_env%matrix_rrcom(i, j)%matrix)
336  CALL dbcsr_copy(vcd_env%matrix_rrcom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
337  CALL dbcsr_init_p(vcd_env%matrix_dcom(i, j)%matrix)
338  CALL dbcsr_copy(vcd_env%matrix_dcom(i, j)%matrix, matrix_ks(1)%matrix)
339  CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0._dp)
340 
341  CALL dbcsr_init_p(vcd_env%matrix_r_rxvr(i, j)%matrix)
342  CALL dbcsr_copy(vcd_env%matrix_r_rxvr(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
343  CALL dbcsr_set(vcd_env%matrix_r_rxvr(i, j)%matrix, 0._dp)
344 
345  CALL dbcsr_init_p(vcd_env%matrix_rxvr_r(i, j)%matrix)
346  CALL dbcsr_copy(vcd_env%matrix_rxvr_r(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
347  CALL dbcsr_set(vcd_env%matrix_rxvr_r(i, j)%matrix, 0._dp)
348 
349  CALL dbcsr_init_p(vcd_env%matrix_r_doublecom(i, j)%matrix)
350  CALL dbcsr_copy(vcd_env%matrix_r_doublecom(i, j)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
351  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  DO ispin = 1, nspins
357  DO i = 1, 3
358  CALL dbcsr_init_p(vcd_env%matrix_hr(ispin, i)%matrix)
359  CALL dbcsr_copy(vcd_env%matrix_hr(ispin, i)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
360 
361  CALL dbcsr_init_p(vcd_env%matrix_rh(ispin, i)%matrix)
362  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  DO i = 1, 3
368  CALL dbcsr_init_p(vcd_env%matrix_drpnl(i)%matrix)
369  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  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  direction_or=.true.)
378  CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
379  direction_or=.true., rc=[0._dp, 0._dp, 0._dp])
380  CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
381  CALL build_matrix_r_vhxc(vcd_env%matrix_hr, qs_env, [0._dp, 0._dp, 0._dp])
382 
383  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  direction_or=.false.)
387  CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
388  direction_or=.false., rc=[0._dp, 0._dp, 0._dp])
389  CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, [0._dp, 0._dp, 0._dp])
390  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  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  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  dft_control%qs_control%eps_ppnl, particle_set, cell, .false.)
401 
402  ! lin_mom
403  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  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  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  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  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  "PRINT%PROGRAM_RUN_INFO")
427 
428  CALL timestop(handle)
429 
430  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  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  CALL timeset(routinen, handle)
448 
449  ! We can't run a NVPT/MFPT calculation without the coefficients dC/dR.
450  CALL dcdr_env_cleanup(qs_env, vcd_env%dcdr_env)
451 
452  DEALLOCATE (vcd_env%apt_el_nvpt)
453  DEALLOCATE (vcd_env%apt_nuc_nvpt)
454  DEALLOCATE (vcd_env%apt_total_nvpt)
455  DEALLOCATE (vcd_env%aat_atom_nvpt)
456  DEALLOCATE (vcd_env%aat_atom_mfp)
457 
458  CALL cp_fm_release(vcd_env%dCV)
459  CALL cp_fm_release(vcd_env%dCV_prime)
460  CALL cp_fm_release(vcd_env%op_dV)
461  CALL cp_fm_release(vcd_env%op_dB)
462 
463  CALL cp_fm_release(vcd_env%dCB)
464  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  CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der)
472  CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_right)
473  CALL dbcsr_deallocate_matrix_set(vcd_env%moments_der_left)
474  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_difdip2)
475  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdV)
476  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dSdB)
477  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hxc_dsdv)
478  CALL dbcsr_deallocate_matrix_set(vcd_env%hcom)
479  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rcomr)
480  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rrcom)
481  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_dcom)
482  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_hr)
483  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rh)
484  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_drpnl)
485  CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao)
486  CALL dbcsr_deallocate_matrix_set(vcd_env%dipvel_ao_delta)
487  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxrv)
488  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_rxvr)
489  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_rxvr_r)
490  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_r_doublecom)
491  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp_33)
492  CALL dbcsr_deallocate_matrix_set(vcd_env%matrix_nosym_temp2_33)
493  CALL timestop(handle)
494 
495  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  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  REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
522  TYPE(cp_fm_type), POINTER :: mo_coeff
523  TYPE(cp_logger_type), POINTER :: logger
524  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  file_exists = .false.
529 
530  CALL timeset(routinen, handle)
531 
532  NULLIFY (mos, para_env, logger, print_key, vecbuffer)
533  logger => cp_get_default_logger()
534 
535  iounit = cp_print_key_unit_nr(logger, linres_section, &
536  "PRINT%PROGRAM_RUN_INFO", extension=".Log")
537 
538  CALL get_qs_env(qs_env=qs_env, &
539  para_env=para_env, &
540  mos=mos)
541 
542  nspins = SIZE(mos)
543 
544  rst_unit = -1
545  IF (para_env%is_source()) THEN
546  CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
547  n_rep_val=n_rep_val)
548 
549  CALL xstring(tag, ia, ie)
550  my_middle = "RESTART-"//tag(ia:ie)//trim("-")//trim(adjustl(cp_to_string(beta))) &
551  //trim("-")//trim(adjustl(cp_to_string(lambda)))
552 
553  IF (n_rep_val > 0) THEN
554  CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
555  CALL xstring(filename, ia, ie)
556  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  print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
560  filename = cp_print_key_generate_filename(logger, print_key, &
561  extension=".lr", middle_name=trim(my_middle), my_local=.false.)
562  END IF
563  INQUIRE (file=filename, exist=file_exists)
564  !
565  ! open file
566  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  unit_number=rst_unit)
573 
574  IF (iounit > 0) WRITE (iounit, "(T2,A)") &
575  "LINRES| Reading response wavefunctions from the restart file <"//trim(adjustl(filename))//">"
576  ELSE
577  IF (iounit > 0) WRITE (iounit, "(T2,A)") &
578  "LINRES| Restart file <"//trim(adjustl(filename))//"> not found"
579  END IF
580  END IF
581 
582  CALL para_env%bcast(file_exists)
583 
584  IF (file_exists) THEN
585 
586  CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
587  CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
588 
589  ALLOCATE (vecbuffer(nao, max_block))
590  !
591  ! read headers
592  IF (rst_unit > 0) READ (rst_unit, iostat=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
593  CALL para_env%bcast(iostat)
594 
595  CALL para_env%bcast(beta_tmp)
596  CALL para_env%bcast(lambda_tmp)
597  CALL para_env%bcast(nspins_tmp)
598  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  IF (nspins_tmp .NE. nspins) THEN
603  cpabort("nspins not consistent")
604  END IF
605  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  IF (lambda_tmp .NE. lambda) cpabort("lambda not consistent")
609  IF (beta_tmp .NE. beta) cpabort("beta not consistent")
610  !
611  DO ispin = 1, nspins
612  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
613  CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
614  !
615  IF (rst_unit > 0) READ (rst_unit) nmo_tmp
616  CALL para_env%bcast(nmo_tmp)
617  IF (nmo_tmp .NE. nmo) cpabort("nmo not consistent")
618  !
619  ! read the response
620  DO i = 1, nmo, max(max_block, 1)
621  i_block = min(max_block, nmo - i + 1)
622  DO j = 1, i_block
623  IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
624  END DO
625  CALL para_env%bcast(vecbuffer)
626  CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
627  END DO
628  END DO
629 
630  IF (iostat /= 0) THEN
631  IF (iounit > 0) WRITE (iounit, "(T2,A)") &
632  "LINRES| Restart file <"//trim(adjustl(filename))//"> not found"
633  END IF
634 
635  DEALLOCATE (vecbuffer)
636 
637  END IF
638 
639  IF (para_env%is_source()) THEN
640  IF (file_exists) CALL close_file(unit_number=rst_unit)
641  END IF
642 
643  CALL timestop(handle)
644 
645  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  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  REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
672  TYPE(cp_fm_type), POINTER :: mo_coeff
673  TYPE(cp_logger_type), POINTER :: logger
674  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  NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
679 
680  CALL timeset(routinen, handle)
681 
682  logger => cp_get_default_logger()
683 
684  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  "PRINT%PROGRAM_RUN_INFO", extension=".Log")
690 
691  CALL get_qs_env(qs_env=qs_env, &
692  mos=mos, &
693  para_env=para_env)
694 
695  nspins = SIZE(mos)
696 
697  my_status = "REPLACE"
698  my_pos = "REWIND"
699  CALL xstring(tag, ia, ie)
700  my_middle = "RESTART-"//tag(ia:ie)//trim("-")//trim(adjustl(cp_to_string(beta))) &
701  //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  file_position=trim(my_pos), file_action="WRITE", file_form="UNFORMATTED")
705 
706  filename = cp_print_key_generate_filename(logger, print_key, &
707  extension=".lr", middle_name=trim(my_middle), my_local=.false.)
708 
709  IF (iounit > 0) THEN
710  WRITE (unit=iounit, fmt="(T2,A)") &
711  "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  CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
718  CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
719  ALLOCATE (vecbuffer(nao, max_block))
720 
721  IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
722 
723  DO ispin = 1, nspins
724  CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
725 
726  IF (rst_unit > 0) WRITE (rst_unit) nmo
727 
728  DO i = 1, nmo, max(max_block, 1)
729  i_block = min(max_block, nmo - i + 1)
730  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  DO j = 1, i_block
735  IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
736  END DO
737  END DO
738  END DO
739 
740  DEALLOCATE (vecbuffer)
741 
742  CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
743  "PRINT%RESTART")
744  END IF
745 
746  CALL timestop(handle)
747 
748  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  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  REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_el_nvpt, apt_nuc_dcdr, &
768  apt_nuc_nvpt, apt_total_dcdr, &
769  apt_total_nvpt
770  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  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
778  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
779  TYPE(section_vals_type), POINTER :: vcd_section
780 
781  CALL timeset(routinen, handle)
782 
783  NULLIFY (logger)
784 
785  logger => cp_get_default_logger()
786  output_unit = cp_logger_get_default_io_unit(logger)
787 
788  vcd_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%VCD")
789 
790  NULLIFY (particle_set)
791  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
792  natom = SIZE(particle_set)
793  nsubset = SIZE(molecule_set)
794 
795  apt_el_dcdr => vcd_env%dcdr_env%apt_el_dcdr
796  apt_nuc_dcdr => vcd_env%dcdr_env%apt_nuc_dcdr
797  apt_total_dcdr => vcd_env%dcdr_env%apt_total_dcdr
798  apt_subset_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_subset
799  apt_center_dcdr => vcd_env%dcdr_env%apt_el_dcdr_per_center
800 
801  apt_el_nvpt => vcd_env%apt_el_nvpt
802  apt_nuc_nvpt => vcd_env%apt_nuc_nvpt
803  apt_total_nvpt => vcd_env%apt_total_nvpt
804 
805  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
806  'APT | Write the final APT matrix per atom (Position perturbation)'
807  DO l = 1, natom
808  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
809  'APT | Atom', l, ' - GAPT ', &
810  (apt_total_dcdr(1, 1, l) &
811  + apt_total_dcdr(2, 2, l) &
812  + apt_total_dcdr(3, 3, l))/3._dp
813  DO i = 1, 3
814  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  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
819  'NVP | Write the final APT matrix per atom (Velocity perturbation)'
820  DO l = 1, natom
821  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3,A,F15.6)") &
822  'NVP | Atom', l, ' - GAPT ', &
823  (apt_total_nvpt(1, 1, l) &
824  + apt_total_nvpt(2, 2, l) &
825  + apt_total_nvpt(3, 3, l))/3._dp
826  DO i = 1, 3
827  IF (vcd_env%output_unit > 0) &
828  WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
829  "NVP | ", apt_total_nvpt(i, :, l)
830  END DO
831  END DO
832 
833  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
834  'NVP | Write the final AAT matrix per atom (Velocity perturbation)'
835  DO l = 1, natom
836  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
837  'NVP | Atom', l
838  DO i = 1, 3
839  IF (vcd_env%output_unit > 0) &
840  WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
841  "NVP | ", vcd_env%aat_atom_nvpt(i, :, l)
842  END DO
843  END DO
844 
845  IF (vcd_env%do_mfp) THEN
846  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") &
847  'MFP | Write the final AAT matrix per atom (Magnetic Field perturbation)'
848  DO l = 1, natom
849  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,I3)") &
850  'MFP | Atom', l
851  DO i = 1, 3
852  IF (vcd_env%output_unit > 0) &
853  WRITE (vcd_env%output_unit, "(A,F15.6,F15.6,F15.6)") &
854  "MFP | ", vcd_env%aat_atom_mfp(i, :, l)
855  END DO
856  END DO
857  END IF
858 
859  ! Get the dipole
860  CALL get_qs_env(qs_env, results=results)
861  description = "[DIPOLE]"
862  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  sum_rule_0 = 0._dp
866  sum_rule_1 = 0._dp
867  sum_rule_2 = 0._dp
868  sum_rule_0_second = 0._dp
869  sum_rule_2_second = 0._dp
870  sum_rule_3_second = 0._dp
871  sum_rule_3_mfp = 0._dp
872  standard_deviation = 0._dp
873  standard_deviation_sum = 0._dp
874 
875  DO alpha = 1, 3
876  DO beta = 1, 3
877  ! 0: sum_lambda apt(alpha, beta, lambda)
878  DO lambda = 1, natom
879  sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
880  + apt_total_dcdr(alpha, beta, lambda)
881  sum_rule_0_second(alpha, beta) = sum_rule_0_second(alpha, beta) &
882  + apt_total_nvpt(alpha, beta, lambda)
883  END DO
884 
885  ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
886  DO gamma = 1, 3
887  sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
888  + 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  DO lambda = 1, natom
893  DO gamma = 1, 3
894  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  *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  *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  DO lambda = 1, natom
909  sum_rule_3_second(alpha, beta) = sum_rule_3_second(alpha, beta) &
910  + 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  IF (vcd_env%do_mfp) THEN
915  ! 3: 2c * sum_lambda aat(alpha, beta, lambda)
916  DO lambda = 1, natom
917  sum_rule_3_mfp(alpha, beta) = sum_rule_3_mfp(alpha, beta) &
918  + 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  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A)") "APT | Position perturbation sum rules"
926  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(A,T19,A,T35,A,T50,A,T65,A)") &
927  "APT |", " Total APT", "Dipole", "R * APT", "AAT"
928  standard_deviation_sum = 0._dp
929  DO alpha = 1, 3
930  DO beta = 1, 3
931  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  - mean**2)
935  standard_deviation_sum = standard_deviation_sum + standard_deviation
936 
937  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  "APT | ", &
940  alpha, beta, &
941  sum_rule_0(alpha, beta), &
942  sum_rule_1(alpha, beta), &
943  sum_rule_2(alpha, beta), &
944  sum_rule_3_mfp(alpha, beta), &
945  standard_deviation
946  END DO
947  END DO
948  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
949 
950  IF (vcd_env%output_unit > 0) THEN
951  WRITE (vcd_env%output_unit, "(A)") "NVP | Velocity perturbation sum rules"
952  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  standard_deviation_sum = 0._dp
956  DO alpha = 1, 3
957  DO beta = 1, 3
958  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  - mean**2)
962  standard_deviation_sum = standard_deviation_sum + standard_deviation
963  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  "NVP | ", &
966  alpha, &
967  beta, &
968  sum_rule_0_second(alpha, beta), &
969  sum_rule_1(alpha, beta), &
970  sum_rule_2_second(alpha, beta), &
971  sum_rule_3_second(alpha, beta), &
972  standard_deviation
973  END DO
974  END DO
975  IF (vcd_env%output_unit > 0) WRITE (vcd_env%output_unit, "(T73,F15.6)") standard_deviation_sum
976 
977  CALL timestop(handle)
978  END SUBROUTINE vcd_print
979 
980 END MODULE qs_vcd_utils
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: commutator_rpnl.F:14
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition: commutator_rpnl.F:439
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition: cp_fm_types.F:768
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition: cp_fm_types.F:901
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition: cp_output_handling.F:704
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_methods.F:17
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_types.F:18
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::use_mom_ref_user
integer, parameter, public use_mom_ref_user
Definition: input_constants.F:496
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
moments_utils
Calculates the moment integrals <a|r^m|b>
Definition: moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition: moments_utils.F:61
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition: orbital_pointers.F:253
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_dcdr_utils
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition: qs_dcdr_utils.F:13
qs_dcdr_utils::dcdr_env_cleanup
subroutine, public dcdr_env_cleanup(qs_env, dcdr_env)
Deallocate the dcdr environment.
Definition: qs_dcdr_utils.F:990
qs_dcdr_utils::dcdr_env_init
subroutine, public dcdr_env_init(dcdr_env, qs_env)
Initialize the dcdr environment.
Definition: qs_dcdr_utils.F:652
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_ks_types
Definition: qs_ks_types.F:15
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition: qs_moments.F:14
qs_moments::build_local_moments_der_matrix
subroutine, public build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, ref_point, moments)
Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b > Optionally sto...
Definition: qs_moments.F:790
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_operators_ao
Definition: qs_operators_ao.F:13
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition: qs_operators_ao.F:74
qs_vcd_ao
Definition: qs_vcd_ao.F:7
qs_vcd_ao::build_tr_matrix
subroutine, public build_tr_matrix(matrix_tr, qs_env, qs_kind_set, basis_type, sab_nl, direction_Or, rc)
Calculation of the product Tr or rT over Cartesian Gaussian functions.
Definition: qs_vcd_ao.F:389
qs_vcd_ao::build_matrix_r_vhxc
subroutine, public build_matrix_r_vhxc(matrix_rv, qs_env, rc)
Commutator of the Hartree+XC potentials with r.
Definition: qs_vcd_ao.F:979
qs_vcd_ao::build_rcore_matrix
subroutine, public build_rcore_matrix(matrix_rcore, qs_env, qs_kind_set, basis_type, sab_nl, rf)
Commutator of the of the local part of the pseudopotential with r.
Definition: qs_vcd_ao.F:629
qs_vcd_ao::build_rpnl_matrix
subroutine, public build_rpnl_matrix(matrix_rv, qs_kind_set, particle_set, sab_all, sap_ppnl, eps_ppnl, cell, ref_point, direction_Or)
Product of r with V_nl. Adapted from build_com_rpnl.
Definition: qs_vcd_ao.F:188
qs_vcd_ao::build_com_rpnl_r
subroutine, public build_com_rpnl_r(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, direction_Or)
Builds the [Vnl, r] * r from either side.
Definition: qs_vcd_ao.F:1550
qs_vcd_utils
Definition: qs_vcd_utils.F:8
qs_vcd_utils::vcd_print
subroutine, public vcd_print(vcd_env, qs_env)
Print the APTs, AATs, and sum rules.
Definition: qs_vcd_utils.F:757
qs_vcd_utils::vcd_env_cleanup
subroutine, public vcd_env_cleanup(qs_env, vcd_env)
Deallocate the vcd environment.
Definition: qs_vcd_utils.F:439
qs_vcd_utils::vcd_write_restart
subroutine, public vcd_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_write_restart.
Definition: qs_vcd_utils.F:658
qs_vcd_utils::vcd_env_init
subroutine, public vcd_env_init(vcd_env, qs_env)
Initialize the vcd environment.
Definition: qs_vcd_utils.F:96
qs_vcd_utils::vcd_read_restart
subroutine, public vcd_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_read_restart.
Definition: qs_vcd_utils.F:508
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::xstring
elemental subroutine, public xstring(string, ia, ib)
...
Definition: string_utilities.F:3381