(git:e5fdd81)
tmc_analysis.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 ! **************************************************************************************************
9 !> \brief module analyses element of the TMC tree element structure
10 !> e.g. density, radial distribution function, dipole correlation,...
11 !> \par History
12 !> 02.2013 created [Mandes Schoenherr]
13 !> \author Mandes
14 ! **************************************************************************************************
15 
16 MODULE tmc_analysis
17  USE cell_types, ONLY: cell_type,&
18  get_cell,&
19  pbc
20  USE cp_files, ONLY: close_file,&
21  open_file
22  USE cp_log_handling, ONLY: cp_to_string
23  USE force_fields_input, ONLY: read_chrg_section
24  USE input_section_types, ONLY: section_vals_get,&
25  section_vals_get_subs_vals,&
26  section_vals_type,&
27  section_vals_val_get
28  USE kinds, ONLY: default_path_length,&
29  default_string_length,&
30  dp
31  USE mathconstants, ONLY: pi
32  USE mathlib, ONLY: diag
33  USE physcon, ONLY: a_mass,&
34  au2a => angstrom,&
35  boltzmann,&
36  joule,&
37  massunit
38  USE tmc_analysis_types, ONLY: &
39  ana_type_default, ana_type_ice, ana_type_sym_xyz, atom_pairs_type, dipole_moment_type, &
40  pair_correl_type, search_pair_in_list, tmc_ana_density_create, tmc_ana_density_file_name, &
41  tmc_ana_dipole_analysis_create, tmc_ana_dipole_moment_create, tmc_ana_displacement_create, &
42  tmc_ana_env_create, tmc_ana_pair_correl_create, tmc_ana_pair_correl_file_name, &
43  tmc_analysis_env
44  USE tmc_calculations, ONLY: get_scaled_cell,&
45  nearest_distance
46  USE tmc_file_io, ONLY: analyse_files_close,&
47  analyse_files_open,&
48  expand_file_name_char,&
49  expand_file_name_temp,&
50  read_element_from_file,&
51  write_dipoles_in_file
52  USE tmc_stati, ONLY: tmc_status_ok,&
53  tmc_status_wait_for_new_task,&
54  tmc_default_restart_in_file_name,&
55  tmc_default_restart_out_file_name,&
56  tmc_default_trajectory_file_name,&
57  tmc_default_unspecified_name
58  USE tmc_tree_build, ONLY: allocate_new_sub_tree_node,&
59  deallocate_sub_tree_node
60  USE tmc_tree_types, ONLY: read_subtree_elem_unformated,&
61  tree_type,&
62  write_subtree_elem_unformated
63  USE tmc_types, ONLY: tmc_atom_type,&
64  tmc_param_type
65 #include "../base/base_uses.f90"
66 
67  IMPLICIT NONE
68 
69  PRIVATE
70 
71  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_analysis'
72 
73  PUBLIC :: tmc_read_ana_input
74  PUBLIC :: analysis_init, do_tmc_analysis, analyze_file_configurations, finalize_tmc_analysis
75  PUBLIC :: analysis_restart_print, analysis_restart_read
76 
77 CONTAINS
78 
79 ! **************************************************************************************************
80 !> \brief creates a new para environment for tmc analysis
81 !> \param tmc_ana_section ...
82 !> \param tmc_ana TMC analysis environment
83 !> \author Mandes 02.2013
84 ! **************************************************************************************************
85  SUBROUTINE tmc_read_ana_input(tmc_ana_section, tmc_ana)
86  TYPE(section_vals_type), POINTER :: tmc_ana_section
87  TYPE(tmc_analysis_env), POINTER :: tmc_ana
88 
89  CHARACTER(LEN=default_path_length) :: c_tmp
90  CHARACTER(LEN=default_string_length), POINTER :: charge_atm(:)
91  INTEGER :: i_tmp, ntot
92  INTEGER, DIMENSION(3) :: nr_bins
93  INTEGER, DIMENSION(:), POINTER :: i_arr_tmp
94  LOGICAL :: explicit, explicit_key, flag
95  REAL(kind=dp), POINTER :: charge(:)
96  TYPE(section_vals_type), POINTER :: tmp_section
97 
98  NULLIFY (tmp_section, charge_atm, i_arr_tmp, charge)
99 
100  cpassert(ASSOCIATED(tmc_ana_section))
101  cpassert(.NOT. ASSOCIATED(tmc_ana))
102 
103  CALL section_vals_get(tmc_ana_section, explicit=explicit)
104  IF (explicit) THEN
105  CALL tmc_ana_env_create(tmc_ana=tmc_ana)
106  ! restarting
107  CALL section_vals_val_get(tmc_ana_section, "RESTART", l_val=tmc_ana%restart)
108  ! file name prefix
109  CALL section_vals_val_get(tmc_ana_section, "PREFIX_ANA_FILES", &
110  c_val=tmc_ana%out_file_prefix)
111  IF (tmc_ana%out_file_prefix .NE. "") THEN
112  tmc_ana%out_file_prefix = trim(tmc_ana%out_file_prefix)//"_"
113  END IF
114 
115  ! density calculation
116  CALL section_vals_val_get(tmc_ana_section, "DENSITY", explicit=explicit_key)
117  IF (explicit_key) THEN
118  CALL section_vals_val_get(tmc_ana_section, "DENSITY", i_vals=i_arr_tmp)
119 
120  IF (SIZE(i_arr_tmp(:)) .EQ. 3) THEN
121  IF (any(i_arr_tmp(:) .LE. 0)) &
122  CALL cp_abort(__location__, "The amount of intervals in each "// &
123  "direction has to be greater than 0.")
124  nr_bins(:) = i_arr_tmp(:)
125  ELSE IF (SIZE(i_arr_tmp(:)) .EQ. 1) THEN
126  IF (any(i_arr_tmp(:) .LE. 0)) &
127  cpabort("The amount of intervals has to be greater than 0.")
128  nr_bins(:) = i_arr_tmp(1)
129  ELSE IF (SIZE(i_arr_tmp(:)) .EQ. 0) THEN
130  nr_bins(:) = 1
131  ELSE
132  cpabort("unknown amount of dimensions for the binning.")
133  END IF
134  CALL tmc_ana_density_create(tmc_ana%density_3d, nr_bins)
135  END IF
136 
137  ! radial distribution function calculation
138  CALL section_vals_val_get(tmc_ana_section, "G_R", explicit=explicit_key)
139  IF (explicit_key) THEN
140  CALL section_vals_val_get(tmc_ana_section, "G_R", i_val=i_tmp)
141  CALL tmc_ana_pair_correl_create(ana_pair_correl=tmc_ana%pair_correl, &
142  nr_bins=i_tmp)
143  END IF
144 
145  ! radial distribution function calculation
146  CALL section_vals_val_get(tmc_ana_section, "CLASSICAL_DIPOLE_MOMENTS", explicit=explicit_key)
147  IF (explicit_key) THEN
148  ! charges for dipoles needed
149  tmp_section => section_vals_get_subs_vals(tmc_ana_section, "CHARGE")
150  CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=i_tmp)
151  IF (explicit) THEN
152  ntot = 0
153  ALLOCATE (charge_atm(i_tmp))
154  ALLOCATE (charge(i_tmp))
155  CALL read_chrg_section(charge_atm, charge, tmp_section, ntot)
156  ELSE
157  CALL cp_abort(__location__, &
158  "to calculate the classical cell dipole moment "// &
159  "the charges has to be specified")
160  END IF
161 
162  CALL tmc_ana_dipole_moment_create(tmc_ana%dip_mom, charge_atm, charge, &
163  tmc_ana%dim_per_elem)
164 
165  IF (ASSOCIATED(charge_atm)) DEALLOCATE (charge_atm)
166  IF (ASSOCIATED(charge)) DEALLOCATE (charge)
167  END IF
168 
169  ! dipole moment analysis
170  CALL section_vals_val_get(tmc_ana_section, "DIPOLE_ANALYSIS", explicit=explicit_key)
171  IF (explicit_key) THEN
172  CALL tmc_ana_dipole_analysis_create(tmc_ana%dip_ana)
173  CALL section_vals_val_get(tmc_ana_section, "DIPOLE_ANALYSIS", c_val=c_tmp)
174  SELECT CASE (trim(c_tmp))
175  CASE (trim(tmc_default_unspecified_name))
176  tmc_ana%dip_ana%ana_type = ana_type_default
177  CASE ("ICE")
178  tmc_ana%dip_ana%ana_type = ana_type_ice
179  CASE ("SYM_XYZ")
180  tmc_ana%dip_ana%ana_type = ana_type_sym_xyz
181  CASE DEFAULT
182  cpwarn('unknown analysis type "'//trim(c_tmp)//'" specified. Set to default.')
183  tmc_ana%dip_ana%ana_type = ana_type_default
184  END SELECT
185  END IF
186 
187  END IF
188 
189  ! cell displacement (deviation)
190  CALL section_vals_val_get(tmc_ana_section, "DEVIATION", l_val=flag)
191  IF (flag) THEN
192  CALL tmc_ana_displacement_create(ana_disp=tmc_ana%displace, &
193  dim_per_elem=tmc_ana%dim_per_elem)
194  END IF
195  END SUBROUTINE tmc_read_ana_input
196 
197 ! **************************************************************************************************
198 !> \brief initialize all the necessarry analysis structures
199 !> \param ana_env ...
200 !> \param nr_dim dimension of the pos, frc etc. array
201 !> \author Mandes 02.2013
202 ! **************************************************************************************************
203  SUBROUTINE analysis_init(ana_env, nr_dim)
204  TYPE(tmc_analysis_env), POINTER :: ana_env
205  INTEGER :: nr_dim
206 
207  CHARACTER(LEN=default_path_length) :: tmp_cell_file, tmp_dip_file, tmp_pos_file
208 
209  cpassert(ASSOCIATED(ana_env))
210  cpassert(nr_dim > 0)
211 
212  ana_env%nr_dim = nr_dim
213 
214  ! save file names
215  tmp_pos_file = ana_env%costum_pos_file_name
216  tmp_cell_file = ana_env%costum_cell_file_name
217  tmp_dip_file = ana_env%costum_dip_file_name
218 
219  ! unset all filenames
220  ana_env%costum_pos_file_name = tmc_default_unspecified_name
221  ana_env%costum_cell_file_name = tmc_default_unspecified_name
222  ana_env%costum_dip_file_name = tmc_default_unspecified_name
223 
224  ! set the necessary files for ...
225  ! density
226  IF (ASSOCIATED(ana_env%density_3d)) THEN
227  ana_env%costum_pos_file_name = tmp_pos_file
228  ana_env%costum_cell_file_name = tmp_cell_file
229  END IF
230  ! pair correlation
231  IF (ASSOCIATED(ana_env%pair_correl)) THEN
232  ana_env%costum_pos_file_name = tmp_pos_file
233  ana_env%costum_cell_file_name = tmp_cell_file
234  END IF
235  ! dipole moment
236  IF (ASSOCIATED(ana_env%dip_mom)) THEN
237  ana_env%costum_pos_file_name = tmp_pos_file
238  ana_env%costum_cell_file_name = tmp_cell_file
239  END IF
240  ! dipole analysis
241  IF (ASSOCIATED(ana_env%dip_ana)) THEN
242  ana_env%costum_pos_file_name = tmp_pos_file
243  ana_env%costum_cell_file_name = tmp_cell_file
244  ana_env%costum_dip_file_name = tmp_dip_file
245  END IF
246  ! deviation / displacement
247  IF (ASSOCIATED(ana_env%displace)) THEN
248  ana_env%costum_pos_file_name = tmp_pos_file
249  ana_env%costum_cell_file_name = tmp_cell_file
250  END IF
251 
252  ! init radial distribution function
253  IF (ASSOCIATED(ana_env%pair_correl)) &
254  CALL ana_pair_correl_init(ana_pair_correl=ana_env%pair_correl, &
255  atoms=ana_env%atoms, cell=ana_env%cell)
256  ! init classical dipole moment calculations
257  IF (ASSOCIATED(ana_env%dip_mom)) &
258  CALL ana_dipole_moment_init(ana_dip_mom=ana_env%dip_mom, &
259  atoms=ana_env%atoms)
260  END SUBROUTINE analysis_init
261 
262 ! **************************************************************************************************
263 !> \brief print analysis restart file
264 !> \param ana_env ...
265 !> \param
266 !> \author Mandes 02.2013
267 ! **************************************************************************************************
268  SUBROUTINE analysis_restart_print(ana_env)
269  TYPE(tmc_analysis_env), POINTER :: ana_env
270 
271  CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp, &
272  restart_file_name
273  INTEGER :: file_ptr
274  LOGICAL :: l_tmp
275 
276  cpassert(ASSOCIATED(ana_env))
277  cpassert(ASSOCIATED(ana_env%last_elem))
278  IF (.NOT. ana_env%restart) RETURN
279 
280  WRITE (file_name, fmt='(I9.9)') ana_env%last_elem%nr
281  file_name_tmp = trim(expand_file_name_temp(expand_file_name_char( &
282  trim(ana_env%out_file_prefix)// &
283  tmc_default_restart_out_file_name, &
284  "ana"), ana_env%temperature))
285  restart_file_name = expand_file_name_char(file_name_tmp, &
286  file_name)
287  CALL open_file(file_name=restart_file_name, file_status="REPLACE", &
288  file_action="WRITE", file_form="UNFORMATTED", &
289  unit_number=file_ptr)
290  WRITE (file_ptr) ana_env%temperature
291  CALL write_subtree_elem_unformated(ana_env%last_elem, file_ptr)
292 
293  ! first mention the different kind of anlysis types initialized
294  ! then the variables for each calculation type
295  l_tmp = ASSOCIATED(ana_env%density_3d)
296  WRITE (file_ptr) l_tmp
297  IF (l_tmp) THEN
298  WRITE (file_ptr) ana_env%density_3d%conf_counter, &
299  ana_env%density_3d%nr_bins, &
300  ana_env%density_3d%sum_vol, &
301  ana_env%density_3d%sum_vol2, &
302  ana_env%density_3d%sum_box_length, &
303  ana_env%density_3d%sum_box_length2, &
304  ana_env%density_3d%sum_density, &
305  ana_env%density_3d%sum_dens2
306  END IF
307 
308  l_tmp = ASSOCIATED(ana_env%pair_correl)
309  WRITE (file_ptr) l_tmp
310  IF (l_tmp) THEN
311  WRITE (file_ptr) ana_env%pair_correl%conf_counter, &
312  ana_env%pair_correl%nr_bins, &
313  ana_env%pair_correl%step_length, &
314  ana_env%pair_correl%pairs, &
315  ana_env%pair_correl%g_r
316  END IF
317 
318  l_tmp = ASSOCIATED(ana_env%dip_mom)
319  WRITE (file_ptr) l_tmp
320  IF (l_tmp) THEN
321  WRITE (file_ptr) ana_env%dip_mom%conf_counter, &
322  ana_env%dip_mom%charges, &
323  ana_env%dip_mom%last_dip_cl
324  END IF
325 
326  l_tmp = ASSOCIATED(ana_env%dip_ana)
327  WRITE (file_ptr) l_tmp
328  IF (l_tmp) THEN
329  WRITE (file_ptr) ana_env%dip_ana%conf_counter, &
330  ana_env%dip_ana%ana_type, &
331  ana_env%dip_ana%mu2_pv_s, &
332  ana_env%dip_ana%mu_psv, &
333  ana_env%dip_ana%mu_pv, &
334  ana_env%dip_ana%mu2_pv_mat, &
335  ana_env%dip_ana%mu2_pv_mat
336  END IF
337 
338  l_tmp = ASSOCIATED(ana_env%displace)
339  WRITE (file_ptr) l_tmp
340  IF (l_tmp) THEN
341  WRITE (file_ptr) ana_env%displace%conf_counter, &
342  ana_env%displace%disp
343  END IF
344 
345  CALL close_file(unit_number=file_ptr)
346 
347  file_name_tmp = expand_file_name_char(trim(ana_env%out_file_prefix)// &
348  tmc_default_restart_in_file_name, "ana")
349  file_name = expand_file_name_temp(file_name_tmp, &
350  ana_env%temperature)
351  CALL open_file(file_name=file_name, &
352  file_action="WRITE", file_status="REPLACE", &
353  unit_number=file_ptr)
354  WRITE (file_ptr, *) trim(restart_file_name)
355  CALL close_file(unit_number=file_ptr)
356  END SUBROUTINE analysis_restart_print
357 
358 ! **************************************************************************************************
359 !> \brief read analysis restart file
360 !> \param ana_env ...
361 !> \param elem ...
362 !> \param
363 !> \author Mandes 02.2013
364 ! **************************************************************************************************
365  SUBROUTINE analysis_restart_read(ana_env, elem)
366  TYPE(tmc_analysis_env), POINTER :: ana_env
367  TYPE(tree_type), POINTER :: elem
368 
369  CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp
370  INTEGER :: file_ptr
371  LOGICAL :: l_tmp
372  REAL(kind=dp) :: temp
373 
374  cpassert(ASSOCIATED(ana_env))
375  cpassert(ASSOCIATED(elem))
376  IF (.NOT. ana_env%restart) RETURN
377 
378  file_name_tmp = expand_file_name_char(trim(ana_env%out_file_prefix)// &
379  tmc_default_restart_in_file_name, "ana")
380  file_name = expand_file_name_temp(file_name_tmp, &
381  ana_env%temperature)
382  INQUIRE (file=file_name, exist=l_tmp)
383  IF (l_tmp) THEN
384  CALL open_file(file_name=file_name, file_status="OLD", &
385  file_action="READ", unit_number=file_ptr)
386  READ (file_ptr, *) file_name_tmp
387  CALL close_file(unit_number=file_ptr)
388 
389  CALL open_file(file_name=file_name_tmp, file_status="OLD", file_form="UNFORMATTED", &
390  file_action="READ", unit_number=file_ptr)
391  READ (file_ptr) temp
392  cpassert(ana_env%temperature .EQ. temp)
393  ana_env%last_elem => elem
394  CALL read_subtree_elem_unformated(elem, file_ptr)
395 
396  ! first mention the different kind of anlysis types initialized
397  ! then the variables for each calculation type
398  READ (file_ptr) l_tmp
399  cpassert(ASSOCIATED(ana_env%density_3d) .EQV. l_tmp)
400  IF (l_tmp) THEN
401  READ (file_ptr) ana_env%density_3d%conf_counter, &
402  ana_env%density_3d%nr_bins, &
403  ana_env%density_3d%sum_vol, &
404  ana_env%density_3d%sum_vol2, &
405  ana_env%density_3d%sum_box_length, &
406  ana_env%density_3d%sum_box_length2, &
407  ana_env%density_3d%sum_density, &
408  ana_env%density_3d%sum_dens2
409  END IF
410 
411  READ (file_ptr) l_tmp
412  cpassert(ASSOCIATED(ana_env%pair_correl) .EQV. l_tmp)
413  IF (l_tmp) THEN
414  READ (file_ptr) ana_env%pair_correl%conf_counter, &
415  ana_env%pair_correl%nr_bins, &
416  ana_env%pair_correl%step_length, &
417  ana_env%pair_correl%pairs, &
418  ana_env%pair_correl%g_r
419  END IF
420 
421  READ (file_ptr) l_tmp
422  cpassert(ASSOCIATED(ana_env%dip_mom) .EQV. l_tmp)
423  IF (l_tmp) THEN
424  READ (file_ptr) ana_env%dip_mom%conf_counter, &
425  ana_env%dip_mom%charges, &
426  ana_env%dip_mom%last_dip_cl
427  END IF
428 
429  READ (file_ptr) l_tmp
430  cpassert(ASSOCIATED(ana_env%dip_ana) .EQV. l_tmp)
431  IF (l_tmp) THEN
432  READ (file_ptr) ana_env%dip_ana%conf_counter, &
433  ana_env%dip_ana%ana_type, &
434  ana_env%dip_ana%mu2_pv_s, &
435  ana_env%dip_ana%mu_psv, &
436  ana_env%dip_ana%mu_pv, &
437  ana_env%dip_ana%mu2_pv_mat, &
438  ana_env%dip_ana%mu2_pv_mat
439  END IF
440 
441  READ (file_ptr) l_tmp
442  cpassert(ASSOCIATED(ana_env%displace) .EQV. l_tmp)
443  IF (l_tmp) THEN
444  READ (file_ptr) ana_env%displace%conf_counter, &
445  ana_env%displace%disp
446  END IF
447 
448  CALL close_file(unit_number=file_ptr)
449  elem => null()
450  END IF
451  END SUBROUTINE analysis_restart_read
452 
453 ! **************************************************************************************************
454 !> \brief call all the necessarry analysis routines
455 !> analysis the previous element with the weight of the different
456 !> configuration numbers
457 !> and stores the actual in the structur % last_elem
458 !> afterwards the previous configuration can be deallocated (outside)
459 !> \param elem ...
460 !> \param ana_env ...
461 !> \param
462 !> \author Mandes 02.2013
463 ! **************************************************************************************************
464  SUBROUTINE do_tmc_analysis(elem, ana_env)
465  TYPE(tree_type), POINTER :: elem
466  TYPE(tmc_analysis_env), POINTER :: ana_env
467 
468  CHARACTER(LEN=*), PARAMETER :: routinen = 'do_tmc_analysis'
469 
470  INTEGER :: handle, weight_act
471  REAL(kind=dp), DIMENSION(3) :: dip_tmp
472  TYPE(tree_type), POINTER :: elem_tmp
473 
474  cpassert(ASSOCIATED(elem))
475  cpassert(ASSOCIATED(ana_env))
476 
477  ! start the timing
478  CALL timeset(routinen, handle)
479  IF (ASSOCIATED(ana_env%last_elem) .AND. &
480  (ana_env%last_elem%nr .LT. elem%nr)) THEN
481  weight_act = elem%nr - ana_env%last_elem%nr
482  ! calculates the 3 dimensional distributed density
483  IF (ASSOCIATED(ana_env%density_3d)) &
484  CALL calc_density_3d(elem=ana_env%last_elem, &
485  weight=weight_act, atoms=ana_env%atoms, &
486  ana_env=ana_env)
487  ! calculated the radial distribution function for each atom type
488  IF (ASSOCIATED(ana_env%pair_correl)) &
489  CALL calc_paircorrelation(elem=ana_env%last_elem, weight=weight_act, &
490  atoms=ana_env%atoms, ana_env=ana_env)
491  ! calculates the classical dipole moments
492  IF (ASSOCIATED(ana_env%dip_mom)) &
493  CALL calc_dipole_moment(elem=ana_env%last_elem, weight=weight_act, &
494  ana_env=ana_env)
495  ! calculates the dipole moments analysis and dielectric constant
496  IF (ASSOCIATED(ana_env%dip_ana)) THEN
497  ! in symmetric case use also the dipoles
498  ! (-x,y,z) .. .. (-x,-y,z).... (-x,-y-z) all have the same energy
499  IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) THEN
500  ! (-x,y,z)
501  ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
502  dip_tmp(:) = ana_env%last_elem%dipole(:)
503  IF (ASSOCIATED(ana_env%dip_mom)) &
504  ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
505  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
506  ana_env=ana_env)
507  ! (-x,-y,z)
508  ana_env%last_elem%dipole(:) = dip_tmp(:)
509  ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
510  dip_tmp(:) = ana_env%last_elem%dipole(:)
511  IF (ASSOCIATED(ana_env%dip_mom)) &
512  ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
513  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
514  ana_env=ana_env)
515  ! (-x,-y,-z)
516  ana_env%last_elem%dipole(:) = dip_tmp(:)
517  ana_env%last_elem%dipole(3) = -ana_env%last_elem%dipole(3)
518  dip_tmp(:) = ana_env%last_elem%dipole(:)
519  IF (ASSOCIATED(ana_env%dip_mom)) &
520  ana_env%dip_mom%last_dip_cl(3) = -ana_env%dip_mom%last_dip_cl(3)
521  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
522  ana_env=ana_env)
523  ! (x,-y,-z)
524  ana_env%last_elem%dipole(:) = dip_tmp(:)
525  ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
526  dip_tmp(:) = ana_env%last_elem%dipole(:)
527  IF (ASSOCIATED(ana_env%dip_mom)) &
528  ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
529  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
530  ana_env=ana_env)
531  ! (x,y,-z)
532  ana_env%last_elem%dipole(:) = dip_tmp(:)
533  ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
534  dip_tmp(:) = ana_env%last_elem%dipole(:)
535  IF (ASSOCIATED(ana_env%dip_mom)) &
536  ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
537  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
538  ana_env=ana_env)
539  ! (-x,y,-z)
540  ana_env%last_elem%dipole(:) = dip_tmp(:)
541  ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
542  dip_tmp(:) = ana_env%last_elem%dipole(:)
543  IF (ASSOCIATED(ana_env%dip_mom)) &
544  ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
545  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
546  ana_env=ana_env)
547  ! (x,-y,z)
548  ana_env%last_elem%dipole(:) = dip_tmp(:)
549  ana_env%last_elem%dipole(:) = -ana_env%last_elem%dipole(:)
550  dip_tmp(:) = ana_env%last_elem%dipole(:)
551  IF (ASSOCIATED(ana_env%dip_mom)) &
552  ana_env%dip_mom%last_dip_cl(:) = -ana_env%dip_mom%last_dip_cl(:)
553  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
554  ana_env=ana_env)
555  ! back to (x,y,z)
556  ana_env%last_elem%dipole(:) = dip_tmp(:)
557  ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
558  dip_tmp(:) = ana_env%last_elem%dipole(:)
559  IF (ASSOCIATED(ana_env%dip_mom)) &
560  ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
561  END IF
562  CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
563  ana_env=ana_env)
564  CALL print_act_dipole_analysis(elem=ana_env%last_elem, &
565  ana_env=ana_env)
566  END IF
567 
568  ! calculates the cell displacement from last cell
569  IF (ASSOCIATED(ana_env%displace) .AND. weight_act .GT. 0) &
570  CALL calc_displacement(elem=elem, ana_env=ana_env)
571  END IF
572  ! swap elem with last elem, to delete original last element and store the actual one
573  elem_tmp => ana_env%last_elem
574  ana_env%last_elem => elem
575  elem => elem_tmp
576  ! end the timing
577  CALL timestop(handle)
578  END SUBROUTINE do_tmc_analysis
579 
580 ! **************************************************************************************************
581 !> \brief call all the necessarry analysis printing routines
582 !> \param ana_env ...
583 !> \param
584 !> \author Mandes 02.2013
585 ! **************************************************************************************************
586  SUBROUTINE finalize_tmc_analysis(ana_env)
587  TYPE(tmc_analysis_env), POINTER :: ana_env
588 
589  CHARACTER(LEN=*), PARAMETER :: routinen = 'finalize_tmc_analysis'
590 
591  INTEGER :: handle
592 
593  cpassert(ASSOCIATED(ana_env))
594 
595  ! start the timing
596  CALL timeset(routinen, handle)
597  IF (ASSOCIATED(ana_env%density_3d)) THEN
598  IF (ana_env%density_3d%conf_counter .GT. 0) &
599  CALL print_density_3d(ana_env=ana_env)
600  END IF
601  IF (ASSOCIATED(ana_env%pair_correl)) THEN
602  IF (ana_env%pair_correl%conf_counter .GT. 0) &
603  CALL print_paircorrelation(ana_env=ana_env)
604  END IF
605  IF (ASSOCIATED(ana_env%dip_mom)) THEN
606  IF (ana_env%dip_mom%conf_counter .GT. 0) &
607  CALL print_dipole_moment(ana_env)
608  END IF
609  IF (ASSOCIATED(ana_env%dip_ana)) THEN
610  IF (ana_env%dip_ana%conf_counter .GT. 0) &
611  CALL print_dipole_analysis(ana_env)
612  END IF
613  IF (ASSOCIATED(ana_env%displace)) THEN
614  IF (ana_env%displace%conf_counter .GT. 0) &
615  CALL print_average_displacement(ana_env)
616  END IF
617 
618  ! end the timing
619  CALL timestop(handle)
620  END SUBROUTINE finalize_tmc_analysis
621 
622 ! **************************************************************************************************
623 !> \brief read the files and analyze the configurations
624 !> \param start_id ...
625 !> \param end_id ...
626 !> \param dir_ind ...
627 !> \param ana_env ...
628 !> \param tmc_params ...
629 !> \author Mandes 03.2013
630 ! **************************************************************************************************
631  SUBROUTINE analyze_file_configurations(start_id, end_id, dir_ind, &
632  ana_env, tmc_params)
633  INTEGER :: start_id, end_id
634  INTEGER, OPTIONAL :: dir_ind
635  TYPE(tmc_analysis_env), POINTER :: ana_env
636  TYPE(tmc_param_type), POINTER :: tmc_params
637 
638  CHARACTER(LEN=*), PARAMETER :: routinen = 'analyze_file_configurations'
639 
640  INTEGER :: conf_nr, handle, nr_dim, stat
641  TYPE(tree_type), POINTER :: elem
642 
643  NULLIFY (elem)
644  conf_nr = -1
645  stat = tmc_status_wait_for_new_task
646  cpassert(ASSOCIATED(ana_env))
647  cpassert(ASSOCIATED(tmc_params))
648 
649  ! start the timing
650  CALL timeset(routinen, handle)
651 
652  ! open the files
653  CALL analyse_files_open(tmc_ana=ana_env, stat=stat, dir_ind=dir_ind)
654  ! set the existence of exact dipoles (from file)
655  IF (ana_env%id_dip .GT. 0) THEN
656  tmc_params%print_dipole = .true.
657  ELSE
658  tmc_params%print_dipole = .false.
659  END IF
660 
661  ! allocate the actual element structure
662  CALL allocate_new_sub_tree_node(tmc_params=tmc_params, next_el=elem, &
663  nr_dim=ana_env%nr_dim)
664 
665  IF (ASSOCIATED(ana_env%last_elem)) conf_nr = ana_env%last_elem%nr
666  nr_dim = SIZE(elem%pos)
667 
668  IF (stat .EQ. tmc_status_ok) THEN
669  conf_loop: DO
670  CALL read_element_from_file(elem=elem, tmc_ana=ana_env, conf_nr=conf_nr, &
671  stat=stat)
672  IF (stat .EQ. tmc_status_wait_for_new_task) THEN
673  CALL deallocate_sub_tree_node(tree_elem=elem)
674  EXIT conf_loop
675  END IF
676  ! if we want just a certain part of the trajectory
677  IF (start_id .LT. 0 .OR. conf_nr .GE. start_id) THEN
678  IF (end_id .LT. 0 .OR. conf_nr .LE. end_id) THEN
679  ! do the analysis calculations
680  CALL do_tmc_analysis(elem=elem, ana_env=ana_env)
681  END IF
682  END IF
683 
684  ! clean temporary element (already analyzed)
685  IF (ASSOCIATED(elem)) THEN
686  CALL deallocate_sub_tree_node(tree_elem=elem)
687  END IF
688  ! if there was no previous element, create a new temp element to write in
689  IF (.NOT. ASSOCIATED(elem)) &
690  CALL allocate_new_sub_tree_node(tmc_params=tmc_params, next_el=elem, &
691  nr_dim=nr_dim)
692  END DO conf_loop
693  END IF
694  ! close the files
695  CALL analyse_files_close(tmc_ana=ana_env)
696 
697  IF (ASSOCIATED(elem)) &
698  CALL deallocate_sub_tree_node(tree_elem=elem)
699 
700  ! end the timing
701  CALL timestop(handle)
702  END SUBROUTINE analyze_file_configurations
703 
704  !============================================================================
705  ! density calculations
706  !============================================================================
707 
708 ! **************************************************************************************************
709 !> \brief calculates the density in rectantangulares
710 !> defined by the number of bins in each direction
711 !> \param elem ...
712 !> \param weight ...
713 !> \param atoms ...
714 !> \param ana_env ...
715 !> \param
716 !> \author Mandes 02.2013
717 ! **************************************************************************************************
718  SUBROUTINE calc_density_3d(elem, weight, atoms, ana_env)
719  TYPE(tree_type), POINTER :: elem
720  INTEGER :: weight
721  TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
722  TYPE(tmc_analysis_env), POINTER :: ana_env
723 
724  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_density_3d'
725 
726  CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp
727  INTEGER :: atom, bin_x, bin_y, bin_z, file_ptr, &
728  handle
729  LOGICAL :: flag
730  REAL(kind=dp) :: mass_total, r_tmp, vol_cell, vol_sub_box
731  REAL(kind=dp), DIMENSION(3) :: atom_pos, cell_size, interval_size
732  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mass_bin
733 
734  NULLIFY (mass_bin)
735 
736  cpassert(ASSOCIATED(elem))
737  cpassert(ASSOCIATED(elem%pos))
738  cpassert(weight .GT. 0)
739  cpassert(ASSOCIATED(atoms))
740  cpassert(ASSOCIATED(ana_env))
741  cpassert(ASSOCIATED(ana_env%cell))
742  cpassert(ASSOCIATED(ana_env%density_3d))
743  cpassert(ASSOCIATED(ana_env%density_3d%sum_density))
744  cpassert(ASSOCIATED(ana_env%density_3d%sum_dens2))
745 
746  ! start the timing
747  CALL timeset(routinen, handle)
748 
749  atom_pos(:) = 0.0_dp
750  cell_size(:) = 0.0_dp
751  interval_size(:) = 0.0_dp
752  mass_total = 0.0_dp
753 
754  bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
755  bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
756  bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
757  ALLOCATE (mass_bin(bin_x, bin_y, bin_z))
758  mass_bin(:, :, :) = 0.0_dp
759 
760  ! if NPT -> box_scale.NE.1.0 use the scaled cell
761  ! ATTENTION then the sub box middle points are not correct in the output
762  ! espacially if we use multiple sub boxes
763  CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
764  abc=cell_size, vol=vol_cell)
765  ! volume summed over configurations for average volume [A]
766  ana_env%density_3d%sum_vol = ana_env%density_3d%sum_vol + &
767  vol_cell*(au2a)**3*weight
768  ana_env%density_3d%sum_vol2 = ana_env%density_3d%sum_vol2 + &
769  (vol_cell*(au2a)**3)**2*weight
770 
771  ana_env%density_3d%sum_box_length(:) = ana_env%density_3d%sum_box_length(:) &
772  + cell_size(:)*(au2a)*weight
773  ana_env%density_3d%sum_box_length2(:) = ana_env%density_3d%sum_box_length2(:) &
774  + (cell_size(:)*(au2a))**2*weight
775 
776  ! sub interval length
777  interval_size(1) = cell_size(1)/real(bin_x, dp)
778  interval_size(2) = cell_size(2)/real(bin_y, dp)
779  interval_size(3) = cell_size(3)/real(bin_z, dp)
780 
781  ! volume in [cm^3]
782  vol_cell = vol_cell*(au2a*1e-8)**3
783  vol_sub_box = interval_size(1)*interval_size(2)*interval_size(3)* &
784  (au2a*1e-8)**3
785 
786  ! count every atom
787  DO atom = 1, SIZE(elem%pos), ana_env%dim_per_elem
788 
789  atom_pos(:) = elem%pos(atom:atom + 2)
790  ! fold into box
791  CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
792  vec=atom_pos)
793  ! shifts the box to positive values (before 0,0,0 is the center)
794  atom_pos(:) = atom_pos(:) + 0.5_dp*cell_size(:)
795  ! calculate the index of the sub box
796  bin_x = int(atom_pos(1)/interval_size(1)) + 1
797  bin_y = int(atom_pos(2)/interval_size(2)) + 1
798  bin_z = int(atom_pos(3)/interval_size(3)) + 1
799  cpassert(bin_x .GT. 0 .AND. bin_y .GT. 0 .AND. bin_z .GT. 0)
800  cpassert(bin_x .LE. SIZE(ana_env%density_3d%sum_density(:, 1, 1)))
801  cpassert(bin_y .LE. SIZE(ana_env%density_3d%sum_density(1, :, 1)))
802  cpassert(bin_z .LE. SIZE(ana_env%density_3d%sum_density(1, 1, :)))
803 
804  ! sum mass in [g] (in bins and total)
805  mass_bin(bin_x, bin_y, bin_z) = mass_bin(bin_x, bin_y, bin_z) + &
806  atoms(int(atom/real(ana_env%dim_per_elem, kind=dp)) + 1)%mass/massunit*1000*a_mass
807  mass_total = mass_total + &
808  atoms(int(atom/real(ana_env%dim_per_elem, kind=dp)) + 1)%mass/massunit*1000*a_mass
809  !mass_bin(bin_x,bin_y,bin_z) = mass_bin(bin_x,bin_y,bin_z) + &
810  ! atoms(INT(atom/REAL(ana_env%dim_per_elem,KIND=dp))+1)%mass/&
811  ! massunit/n_avogadro
812  !mass_total = mass_total + &
813  ! atoms(INT(atom/REAL(ana_env%dim_per_elem,KIND=dp))+1)%mass/&
814  ! massunit/n_avogadro
815  END DO
816  ! check total cell density
817  r_tmp = mass_total/vol_cell - sum(mass_bin(:, :, :))/vol_sub_box/SIZE(mass_bin(:, :, :))
818  cpassert(abs(r_tmp) .LT. 1e-5)
819 
820  ! calculate density (mass per volume) and sum up for average value
821  ana_env%density_3d%sum_density(:, :, :) = ana_env%density_3d%sum_density(:, :, :) + &
822  weight*mass_bin(:, :, :)/vol_sub_box
823 
824  ! calculate density squared ( (mass per volume)^2 ) for variance and sum up for average value
825  ana_env%density_3d%sum_dens2(:, :, :) = ana_env%density_3d%sum_dens2(:, :, :) + &
826  weight*(mass_bin(:, :, :)/vol_sub_box)**2
827 
828  ana_env%density_3d%conf_counter = ana_env%density_3d%conf_counter + weight
829 
830  ! print out the actual and average density in file
831  IF (ana_env%density_3d%print_dens) THEN
832  file_name_tmp = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
833  tmc_default_trajectory_file_name, &
834  ana_env%temperature)
835  file_name = trim(expand_file_name_char(file_name_tmp, &
836  "dens"))
837  INQUIRE (file=file_name, exist=flag)
838  CALL open_file(file_name=file_name, file_status="UNKNOWN", &
839  file_action="WRITE", file_position="APPEND", &
840  unit_number=file_ptr)
841  IF (.NOT. flag) &
842  WRITE (file_ptr, fmt='(A8,11A20)') "# conf_nr", "dens_act[g/cm^3]", &
843  "dens_average[g/cm^3]", "density_variance", &
844  "averages:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z", &
845  "variances:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z"
846  WRITE (file_ptr, fmt="(I8,11F20.10)") ana_env%density_3d%conf_counter + 1 - weight, &
847  sum(mass_bin(:, :, :))/vol_sub_box/SIZE(mass_bin(:, :, :)), &
848  sum(ana_env%density_3d%sum_density(:, :, :))/ &
849  SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
850  REAL(ana_env%density_3d%conf_counter, kind=dp), &
851  sum(ana_env%density_3d%sum_dens2(:, :, :))/ &
852  SIZE(ana_env%density_3d%sum_dens2(:, :, :))/ &
853  REAL(ana_env%density_3d%conf_counter, kind=dp) - &
854  (sum(ana_env%density_3d%sum_density(:, :, :))/ &
855  SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
856  REAL(ana_env%density_3d%conf_counter, kind=dp))**2, &
857  ana_env%density_3d%sum_vol/ &
858  REAL(ana_env%density_3d%conf_counter, kind=dp), &
859  ana_env%density_3d%sum_box_length(:)/ &
860  REAL(ana_env%density_3d%conf_counter, kind=dp), &
861  ana_env%density_3d%sum_vol2/ &
862  REAL(ana_env%density_3d%conf_counter, kind=dp) - &
863  (ana_env%density_3d%sum_vol/ &
864  REAL(ana_env%density_3d%conf_counter, kind=dp))**2, &
865  ana_env%density_3d%sum_box_length2(:)/ &
866  REAL(ana_env%density_3d%conf_counter, kind=dp) - &
867  (ana_env%density_3d%sum_box_length(:)/ &
868  REAL(ana_env%density_3d%conf_counter, kind=dp))**2
869  CALL close_file(unit_number=file_ptr)
870  END IF
871 
872  DEALLOCATE (mass_bin)
873  ! end the timing
874  CALL timestop(handle)
875  END SUBROUTINE calc_density_3d
876 
877 ! **************************************************************************************************
878 !> \brief print the density in rectantangulares
879 !> defined by the number of bins in each direction
880 !> \param ana_env ...
881 !> \param
882 !> \author Mandes 02.2013
883 ! **************************************************************************************************
884  SUBROUTINE print_density_3d(ana_env)
885  TYPE(tmc_analysis_env), POINTER :: ana_env
886 
887  CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA", &
888  routinen = 'print_density_3d'
889 
890  CHARACTER(LEN=default_path_length) :: file_name, file_name_vari
891  INTEGER :: bin_x, bin_y, bin_z, file_ptr_dens, &
892  file_ptr_vari, handle, i, j, k
893  REAL(kind=dp), DIMENSION(3) :: cell_size, interval_size
894 
895  cpassert(ASSOCIATED(ana_env))
896  cpassert(ASSOCIATED(ana_env%density_3d))
897  cpassert(ASSOCIATED(ana_env%density_3d%sum_density))
898  cpassert(ASSOCIATED(ana_env%density_3d%sum_dens2))
899 
900  ! start the timing
901  CALL timeset(routinen, handle)
902 
903  file_name = ""
904  file_name_vari = ""
905 
906  bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
907  bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
908  bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
909  CALL get_cell(cell=ana_env%cell, abc=cell_size)
910  interval_size(1) = cell_size(1)/real(bin_x, kind=dp)*au2a
911  interval_size(2) = cell_size(2)/real(bin_y, kind=dp)*au2a
912  interval_size(3) = cell_size(3)/real(bin_z, kind=dp)*au2a
913 
914  file_name = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
915  tmc_ana_density_file_name, &
916  ana_env%temperature)
917  CALL open_file(file_name=file_name, file_status="REPLACE", &
918  file_action="WRITE", file_position="APPEND", &
919  unit_number=file_ptr_dens)
920  WRITE (file_ptr_dens, fmt='(A,1X,I0,1X,A,3(I0,1X),1X,A,1X,3F10.5)') &
921  "# configurations", ana_env%density_3d%conf_counter, "bins", &
922  ana_env%density_3d%nr_bins, "interval size", interval_size(:)
923  WRITE (file_ptr_dens, fmt='(A,3A10,A20)') "#", " x [A] ", " y [A] ", " z [A] ", " density [g/cm^3] "
924 
925  file_name_vari = expand_file_name_temp(expand_file_name_char( &
926  trim(ana_env%out_file_prefix)// &
927  tmc_ana_density_file_name, "vari"), &
928  ana_env%temperature)
929  CALL open_file(file_name=file_name_vari, file_status="REPLACE", &
930  file_action="WRITE", file_position="APPEND", &
931  unit_number=file_ptr_vari)
932  WRITE (file_ptr_vari, fmt='(A,1X,I0,1X,A,3(I0,1X),1X,A,1X,3F10.5)') &
933  "# configurations", ana_env%density_3d%conf_counter, "bins", &
934  ana_env%density_3d%nr_bins, "interval size", interval_size(:)
935  WRITE (file_ptr_vari, fmt='(A,3A10,A20)') "#", " x [A] ", " y [A] ", " z [A] ", " variance"
936 
937  DO i = 1, SIZE(ana_env%density_3d%sum_density(:, 1, 1))
938  DO j = 1, SIZE(ana_env%density_3d%sum_density(1, :, 1))
939  DO k = 1, SIZE(ana_env%density_3d%sum_density(1, 1, :))
940  WRITE (file_ptr_dens, fmt='(3F10.2,F20.10)') &
941  (i - 0.5_dp)*interval_size(1), (j - 0.5_dp)*interval_size(2), (k - 0.5_dp)*interval_size(3), &
942  ana_env%density_3d%sum_density(i, j, k)/real(ana_env%density_3d%conf_counter, kind=dp)
943  WRITE (file_ptr_vari, fmt='(3F10.2,F20.10)') &
944  (i - 0.5_dp)*interval_size(1), (j - 0.5_dp)*interval_size(2), (k - 0.5_dp)*interval_size(3), &
945  ana_env%density_3d%sum_dens2(i, j, k)/real(ana_env%density_3d%conf_counter, kind=dp) - &
946  (ana_env%density_3d%sum_density(i, j, k)/real(ana_env%density_3d%conf_counter, kind=dp))**2
947  END DO
948  END DO
949  END DO
950  CALL close_file(unit_number=file_ptr_dens)
951  CALL close_file(unit_number=file_ptr_vari)
952 
953  WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
954  WRITE (ana_env%io_unit, fmt="(T2,A,T35,A,T80,A)") "-", "density calculation", "-"
955  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
956  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "used configurations", &
957  cp_to_string(real(ana_env%density_3d%conf_counter, kind=dp))
958  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "average volume", &
959  cp_to_string(ana_env%density_3d%sum_vol/ &
960  REAL(ana_env%density_3d%conf_counter, kind=dp))
961  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "average density in the cell: ", &
962  cp_to_string(sum(ana_env%density_3d%sum_density(:, :, :))/ &
963  SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
964  REAL(ana_env%density_3d%conf_counter, kind=dp))
965  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "density variance:", &
966  cp_to_string(sum(ana_env%density_3d%sum_dens2(:, :, :))/ &
967  SIZE(ana_env%density_3d%sum_dens2(:, :, :))/ &
968  REAL(ana_env%density_3d%conf_counter, kind=dp) - &
969  (sum(ana_env%density_3d%sum_density(:, :, :))/ &
970  SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
971  REAL(ana_env%density_3d%conf_counter, kind=dp))**2)
972  WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
973  IF (ana_env%print_test_output) &
974  WRITE (ana_env%io_unit, *) "TMC|ANALYSIS_CELL_DENSITY_X= ", &
975  sum(ana_env%density_3d%sum_density(:, :, :))/ &
976  SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
977  REAL(ana_env%density_3d%conf_counter, kind=dp)
978  ! end the timing
979  CALL timestop(handle)
980  END SUBROUTINE print_density_3d
981 
982  !============================================================================
983  ! radial distribution function
984  !============================================================================
985 
986 ! **************************************************************************************************
987 !> \brief init radial distribution function structures
988 !> \param ana_pair_correl ...
989 !> \param atoms ...
990 !> \param cell ...
991 !> \param
992 !> \author Mandes 02.2013
993 ! **************************************************************************************************
994  SUBROUTINE ana_pair_correl_init(ana_pair_correl, atoms, cell)
995  TYPE(pair_correl_type), POINTER :: ana_pair_correl
996  TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
997  TYPE(cell_type), POINTER :: cell
998 
999  CHARACTER(LEN=*), PARAMETER :: routinen = 'ana_pair_correl_init'
1000 
1001  INTEGER :: counter, f_n, handle, list, list_ind, s_n
1002  REAL(kind=dp), DIMENSION(3) :: cell_size
1003  TYPE(atom_pairs_type), DIMENSION(:), POINTER :: pairs_tmp
1004 
1005  cpassert(ASSOCIATED(ana_pair_correl))
1006  cpassert(.NOT. ASSOCIATED(ana_pair_correl%g_r))
1007  cpassert(.NOT. ASSOCIATED(ana_pair_correl%pairs))
1008  cpassert(ASSOCIATED(atoms))
1009  cpassert(SIZE(atoms) .GT. 1)
1010  cpassert(ASSOCIATED(cell))
1011 
1012  ! start the timing
1013  CALL timeset(routinen, handle)
1014 
1015  CALL get_cell(cell=cell, abc=cell_size)
1016  IF (ana_pair_correl%nr_bins .LE. 0) THEN
1017  ana_pair_correl%nr_bins = ceiling(maxval(cell_size(:))/2.0_dp/(0.03/au2a))
1018  END IF
1019  ana_pair_correl%step_length = maxval(cell_size(:))/2.0_dp/ &
1020  ana_pair_correl%nr_bins
1021  ana_pair_correl%conf_counter = 0
1022 
1023  counter = 1
1024  ! initialise the atom pairs
1025  ALLOCATE (pairs_tmp(SIZE(atoms)))
1026  DO f_n = 1, SIZE(atoms)
1027  DO s_n = f_n + 1, SIZE(atoms)
1028  ! search if atom pair is already selected
1029  list_ind = search_pair_in_list(pair_list=pairs_tmp, n1=atoms(f_n)%name, &
1030  n2=atoms(s_n)%name, list_end=counter - 1)
1031  ! add to list
1032  IF (list_ind .LT. 0) THEN
1033  pairs_tmp(counter)%f_n = atoms(f_n)%name
1034  pairs_tmp(counter)%s_n = atoms(s_n)%name
1035  pairs_tmp(counter)%pair_count = 1
1036  counter = counter + 1
1037  ELSE
1038  pairs_tmp(list_ind)%pair_count = pairs_tmp(list_ind)%pair_count + 1
1039  END IF
1040  END DO
1041  END DO
1042 
1043  ALLOCATE (ana_pair_correl%pairs(counter - 1))
1044  DO list = 1, counter - 1
1045  ana_pair_correl%pairs(list)%f_n = pairs_tmp(list)%f_n
1046  ana_pair_correl%pairs(list)%s_n = pairs_tmp(list)%s_n
1047  ana_pair_correl%pairs(list)%pair_count = pairs_tmp(list)%pair_count
1048  END DO
1049  DEALLOCATE (pairs_tmp)
1050 
1051  ALLOCATE (ana_pair_correl%g_r(SIZE(ana_pair_correl%pairs(:)), ana_pair_correl%nr_bins))
1052  ana_pair_correl%g_r = 0.0_dp
1053  ! end the timing
1054  CALL timestop(handle)
1055  END SUBROUTINE ana_pair_correl_init
1056 
1057 ! **************************************************************************************************
1058 !> \brief calculates the radial distribution function
1059 !> \param elem ...
1060 !> \param weight ...
1061 !> \param atoms ...
1062 !> \param ana_env ...
1063 !> \param
1064 !> \author Mandes 02.2013
1065 ! **************************************************************************************************
1066  SUBROUTINE calc_paircorrelation(elem, weight, atoms, ana_env)
1067  TYPE(tree_type), POINTER :: elem
1068  INTEGER :: weight
1069  TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
1070  TYPE(tmc_analysis_env), POINTER :: ana_env
1071 
1072  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_paircorrelation'
1073 
1074  INTEGER :: handle, i, ind, j, pair_ind
1075  REAL(kind=dp) :: dist
1076  REAL(kind=dp), DIMENSION(3) :: cell_size
1077 
1078  cpassert(ASSOCIATED(elem))
1079  cpassert(ASSOCIATED(elem%pos))
1080  cpassert(all(elem%box_scale(:) .GT. 0.0_dp))
1081  cpassert(weight .GT. 0)
1082  cpassert(ASSOCIATED(atoms))
1083  cpassert(ASSOCIATED(ana_env))
1084  cpassert(ASSOCIATED(ana_env%cell))
1085  cpassert(ASSOCIATED(ana_env%pair_correl))
1086  cpassert(ASSOCIATED(ana_env%pair_correl%g_r))
1087  cpassert(ASSOCIATED(ana_env%pair_correl%pairs))
1088 
1089  ! start the timing
1090  CALL timeset(routinen, handle)
1091 
1092  dist = -1.0_dp
1093 
1094  first_elem_loop: DO i = 1, SIZE(elem%pos), ana_env%dim_per_elem
1095  second_elem_loop: DO j = i + 3, SIZE(elem%pos), ana_env%dim_per_elem
1096  dist = nearest_distance(x1=elem%pos(i:i + ana_env%dim_per_elem - 1), &
1097  x2=elem%pos(j:j + ana_env%dim_per_elem - 1), &
1098  cell=ana_env%cell, box_scale=elem%box_scale)
1099  ind = ceiling(dist/ana_env%pair_correl%step_length)
1100  IF (ind .LE. ana_env%pair_correl%nr_bins) THEN
1101  pair_ind = search_pair_in_list(pair_list=ana_env%pair_correl%pairs, &
1102  n1=atoms(int(i/real(ana_env%dim_per_elem, kind=dp)) + 1)%name, &
1103  n2=atoms(int(j/real(ana_env%dim_per_elem, kind=dp)) + 1)%name)
1104  cpassert(pair_ind .GT. 0)
1105  ana_env%pair_correl%g_r(pair_ind, ind) = &
1106  ana_env%pair_correl%g_r(pair_ind, ind) + weight
1107  END IF
1108  END DO second_elem_loop
1109  END DO first_elem_loop
1110  ana_env%pair_correl%conf_counter = ana_env%pair_correl%conf_counter + weight
1111  CALL get_cell(cell=ana_env%cell, abc=cell_size)
1112  ana_env%pair_correl%sum_box_scale = ana_env%pair_correl%sum_box_scale + &
1113  (elem%box_scale(:)*weight)
1114  ! end the timing
1115  CALL timestop(handle)
1116  END SUBROUTINE calc_paircorrelation
1117 
1118 ! **************************************************************************************************
1119 !> \brief print the radial distribution function for each pair of atoms
1120 !> \param ana_env ...
1121 !> \param
1122 !> \author Mandes 02.2013
1123 ! **************************************************************************************************
1124  SUBROUTINE print_paircorrelation(ana_env)
1125  TYPE(tmc_analysis_env), POINTER :: ana_env
1126 
1127  CHARACTER(LEN=*), PARAMETER :: routinen = 'print_paircorrelation'
1128 
1129  CHARACTER(LEN=default_path_length) :: file_name
1130  INTEGER :: bin, file_ptr, handle, pair
1131  REAL(kind=dp) :: aver_box_scale(3), vol, voldr
1132  REAL(kind=dp), DIMENSION(3) :: cell_size
1133 
1134  cpassert(ASSOCIATED(ana_env))
1135  cpassert(ASSOCIATED(ana_env%pair_correl))
1136 
1137  ! start the timing
1138  CALL timeset(routinen, handle)
1139 
1140  CALL get_cell(cell=ana_env%cell, abc=cell_size)
1141  aver_box_scale(:) = ana_env%pair_correl%sum_box_scale(:)/ana_env%pair_correl%conf_counter
1142  vol = (cell_size(1)*aver_box_scale(1))* &
1143  (cell_size(2)*aver_box_scale(2))* &
1144  (cell_size(3)*aver_box_scale(3))
1145 
1146  DO pair = 1, SIZE(ana_env%pair_correl%pairs)
1147  file_name = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1148  tmc_ana_pair_correl_file_name, &
1149  ana_env%temperature)
1150  CALL open_file(file_name=expand_file_name_char( &
1151  expand_file_name_char(file_name, &
1152  ana_env%pair_correl%pairs(pair)%f_n), &
1153  ana_env%pair_correl%pairs(pair)%s_n), &
1154  file_status="REPLACE", &
1155  file_action="WRITE", file_position="APPEND", &
1156  unit_number=file_ptr)
1157  WRITE (file_ptr, *) "# radial distribution function of "// &
1158  trim(ana_env%pair_correl%pairs(pair)%f_n)//" and "// &
1159  trim(ana_env%pair_correl%pairs(pair)%s_n)//" of ", &
1160  ana_env%pair_correl%conf_counter, " configurations"
1161  WRITE (file_ptr, *) "# using a bin size of ", &
1162  ana_env%pair_correl%step_length*au2a, &
1163  "[A] (for Vol changes: referring to the reference cell)"
1164  DO bin = 1, ana_env%pair_correl%nr_bins
1165  voldr = 4.0/3.0*pi*ana_env%pair_correl%step_length**3* &
1166  (real(bin, kind=dp)**3 - real(bin - 1, kind=dp)**3)
1167  WRITE (file_ptr, *) (bin - 0.5)*ana_env%pair_correl%step_length*au2a, &
1168  (ana_env%pair_correl%g_r(pair, bin)/ana_env%pair_correl%conf_counter)/ &
1169  (voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
1170  END DO
1171  CALL close_file(unit_number=file_ptr)
1172 
1173  IF (ana_env%print_test_output) &
1174  WRITE (*, *) "TMC|ANALYSIS_G_R_"// &
1175  trim(ana_env%pair_correl%pairs(pair)%f_n)//"_"// &
1176  trim(ana_env%pair_correl%pairs(pair)%s_n)//"_X= ", &
1177  sum(ana_env%pair_correl%g_r(pair, :)/ana_env%pair_correl%conf_counter/ &
1178  voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
1179  END DO
1180 
1181  ! end the timing
1182  CALL timestop(handle)
1183  END SUBROUTINE print_paircorrelation
1184 
1185  !============================================================================
1186  ! classical cell dipole moment
1187  !============================================================================
1188 
1189 ! **************************************************************************************************
1190 !> \brief init radial distribution function structures
1191 !> \param ana_dip_mom ...
1192 !> \param atoms ...
1193 !> \param
1194 !> \author Mandes 02.2013
1195 ! **************************************************************************************************
1196  SUBROUTINE ana_dipole_moment_init(ana_dip_mom, atoms)
1197  TYPE(dipole_moment_type), POINTER :: ana_dip_mom
1198  TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
1199 
1200  CHARACTER(LEN=*), PARAMETER :: routinen = 'ana_dipole_moment_init'
1201 
1202  INTEGER :: atom, charge, handle
1203 
1204  cpassert(ASSOCIATED(ana_dip_mom))
1205  cpassert(ASSOCIATED(ana_dip_mom%charges_inp))
1206  cpassert(ASSOCIATED(atoms))
1207 
1208  ! start the timing
1209  CALL timeset(routinen, handle)
1210 
1211  ALLOCATE (ana_dip_mom%charges(SIZE(atoms)))
1212  ana_dip_mom%charges = 0.0_dp
1213  ! for every atom searcht the correct charge
1214  DO atom = 1, SIZE(atoms)
1215  charge_loop: DO charge = 1, SIZE(ana_dip_mom%charges_inp)
1216  IF (atoms(atom)%name .EQ. ana_dip_mom%charges_inp(charge)%name) THEN
1217  ana_dip_mom%charges(atom) = ana_dip_mom%charges_inp(charge)%mass
1218  EXIT charge_loop
1219  END IF
1220  END DO charge_loop
1221  END DO
1222 
1223  DEALLOCATE (ana_dip_mom%charges_inp)
1224  ! end the timing
1225  CALL timestop(handle)
1226  END SUBROUTINE ana_dipole_moment_init
1227 
1228 ! **************************************************************************************************
1229 !> \brief calculates the classical cell dipole moment
1230 !> \param elem ...
1231 !> \param weight ...
1232 !> \param ana_env ...
1233 !> \param
1234 !> \author Mandes 02.2013
1235 ! **************************************************************************************************
1236  SUBROUTINE calc_dipole_moment(elem, weight, ana_env)
1237  TYPE(tree_type), POINTER :: elem
1238  INTEGER :: weight
1239  TYPE(tmc_analysis_env), POINTER :: ana_env
1240 
1241  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_dipole_moment'
1242 
1243  CHARACTER(LEN=default_path_length) :: file_name
1244  INTEGER :: handle, i
1245  REAL(kind=dp), DIMENSION(:), POINTER :: dip_cl
1246 
1247  cpassert(ASSOCIATED(elem))
1248  cpassert(ASSOCIATED(elem%pos))
1249  cpassert(ASSOCIATED(ana_env))
1250  cpassert(ASSOCIATED(ana_env%dip_mom))
1251  cpassert(ASSOCIATED(ana_env%dip_mom%charges))
1252 
1253  ! start the timing
1254  CALL timeset(routinen, handle)
1255 
1256  ALLOCATE (dip_cl(ana_env%dim_per_elem))
1257  dip_cl(:) = 0.0_dp
1258 
1259  DO i = 1, SIZE(elem%pos, 1), ana_env%dim_per_elem
1260  dip_cl(:) = dip_cl(:) + elem%pos(i:i + ana_env%dim_per_elem - 1)* &
1261  ana_env%dip_mom%charges(int(i/real(ana_env%dim_per_elem, kind=dp)) + 1)
1262  END DO
1263 
1264  ! if there are no exact dipoles save these ones in element structure
1265  IF (.NOT. ASSOCIATED(elem%dipole)) THEN
1266  ALLOCATE (elem%dipole(ana_env%dim_per_elem))
1267  elem%dipole(:) = dip_cl(:)
1268  END IF
1269 
1270  IF (ana_env%dip_mom%print_cl_dip) THEN
1271  file_name = expand_file_name_temp(tmc_default_trajectory_file_name, &
1272  ana_env%temperature)
1273  CALL write_dipoles_in_file(file_name=file_name, &
1274  conf_nr=ana_env%dip_mom%conf_counter + 1, dip=dip_cl, &
1275  file_ext="dip_cl")
1276  END IF
1277  ana_env%dip_mom%conf_counter = ana_env%dip_mom%conf_counter + weight
1278  ana_env%dip_mom%last_dip_cl(:) = dip_cl
1279 
1280  DEALLOCATE (dip_cl)
1281 
1282  ! end the timing
1283  CALL timestop(handle)
1284  END SUBROUTINE calc_dipole_moment
1285 
1286 ! **************************************************************************************************
1287 !> \brief prints final values for classical cell dipole moment calculation
1288 !> \param ana_env ...
1289 !> \param
1290 !> \author Mandes 02.2013
1291 ! **************************************************************************************************
1292  SUBROUTINE print_dipole_moment(ana_env)
1293  TYPE(tmc_analysis_env), POINTER :: ana_env
1294 
1295  IF (ana_env%print_test_output) &
1296  WRITE (*, *) "TMC|ANALYSIS_FINAL_CLASS_CELL_DIPOLE_MOMENT_X= ", &
1297  ana_env%dip_mom%last_dip_cl(:)
1298  END SUBROUTINE print_dipole_moment
1299 
1300 ! **************************************************************************************************
1301 !> \brief calculates the dipole moment analysis
1302 !> \param elem ...
1303 !> \param weight ...
1304 !> \param ana_env ...
1305 !> \param
1306 !> \author Mandes 03.2013
1307 ! **************************************************************************************************
1308  SUBROUTINE calc_dipole_analysis(elem, weight, ana_env)
1309  TYPE(tree_type), POINTER :: elem
1310  INTEGER :: weight
1311  TYPE(tmc_analysis_env), POINTER :: ana_env
1312 
1313  REAL(kind=dp) :: vol, weight_act
1314  REAL(kind=dp), DIMENSION(3, 3) :: tmp_dip
1315  TYPE(cell_type), POINTER :: scaled_cell
1316 
1317  NULLIFY (scaled_cell)
1318 
1319  cpassert(ASSOCIATED(elem))
1320  cpassert(ASSOCIATED(elem%dipole))
1321  cpassert(ASSOCIATED(ana_env))
1322  cpassert(ASSOCIATED(ana_env%dip_ana))
1323 
1324  weight_act = weight
1325  IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
1326  weight_act = weight_act/real(8.0, kind=dp)
1327 
1328  ! get the volume
1329  ALLOCATE (scaled_cell)
1330  CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, vol=vol, &
1331  scaled_cell=scaled_cell)
1332 
1333  ! fold exact dipole moments using the classical ones
1334  IF (ASSOCIATED(ana_env%dip_mom)) THEN
1335  IF (all(ana_env%dip_mom%last_dip_cl .NE. elem%dipole)) THEN
1336  elem%dipole = pbc(r=elem%dipole(:) - ana_env%dip_mom%last_dip_cl, &
1337  cell=scaled_cell) + ana_env%dip_mom%last_dip_cl
1338  END IF
1339  END IF
1340 
1341  ana_env%dip_ana%conf_counter = ana_env%dip_ana%conf_counter + weight_act
1342 
1343  ! dipole sqared absolut value summed and weight_acted with volume and conf weight_act
1344  ana_env%dip_ana%mu2_pv_s = ana_env%dip_ana%mu2_pv_s + &
1345  dot_product(elem%dipole(:), elem%dipole(:))/vol*weight_act
1346 
1347  tmp_dip(:, :) = 0.0_dp
1348  tmp_dip(:, 1) = elem%dipole(:)
1349 
1350  ! dipole sum, weight_acted with volume and conf weight_act
1351  ana_env%dip_ana%mu_pv(:) = ana_env%dip_ana%mu_pv(:) + &
1352  tmp_dip(:, 1)/vol*weight_act
1353 
1354  ! dipole sum, weight_acted with square root of volume and conf weight_act
1355  ana_env%dip_ana%mu_psv(:) = ana_env%dip_ana%mu_psv(:) + &
1356  tmp_dip(:, 1)/sqrt(vol)*weight_act
1357 
1358  ! dipole squared sum, weight_acted with volume and conf weight_act
1359  ana_env%dip_ana%mu2_pv(:) = ana_env%dip_ana%mu2_pv(:) + &
1360  tmp_dip(:, 1)**2/vol*weight_act
1361 
1362  ! calculate the directional average with componentwise correlation per volume
1363  tmp_dip(:, :) = matmul(tmp_dip(:, :), transpose(tmp_dip(:, :)))
1364  ana_env%dip_ana%mu2_pv_mat(:, :) = ana_env%dip_ana%mu2_pv_mat(:, :) + &
1365  tmp_dip(:, :)/vol*weight_act
1366 
1367  END SUBROUTINE calc_dipole_analysis
1368 
1369 ! **************************************************************************************************
1370 !> \brief prints the actual dipole moment analysis (trajectories)
1371 !> \param elem ...
1372 !> \param ana_env ...
1373 !> \param
1374 !> \author Mandes 03.2013
1375 ! **************************************************************************************************
1376  SUBROUTINE print_act_dipole_analysis(elem, ana_env)
1377  TYPE(tree_type), POINTER :: elem
1378  TYPE(tmc_analysis_env), POINTER :: ana_env
1379 
1380  CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp
1381  INTEGER :: counter_tmp, file_ptr
1382  LOGICAL :: flag
1383  REAL(kind=dp) :: diel_const, diel_const_norm, &
1384  diel_const_sym, e0, kb
1385  REAL(kind=dp), DIMENSION(3, 3) :: tmp_dip
1386 
1387  kb = boltzmann/joule
1388  counter_tmp = int(ana_env%dip_ana%conf_counter)
1389 
1390  ! TODO get correct constant using physcon
1391  e0 = 0.07957747154594767_dp !e^2*a0*me*hbar^-2
1392  diel_const_norm = 1/(3.0_dp*e0*kb*ana_env%temperature)
1393 
1394  file_name = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1395  tmc_default_trajectory_file_name, &
1396  ana_env%temperature)
1397  CALL write_dipoles_in_file(file_name=file_name, &
1398  conf_nr=int(ana_env%dip_ana%conf_counter) + 1, dip=elem%dipole, &
1399  file_ext="dip_folded")
1400 
1401  ! set output file name
1402  file_name_tmp = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1403  tmc_default_trajectory_file_name, &
1404  ana_env%temperature)
1405 
1406  SELECT CASE (ana_env%dip_ana%ana_type)
1407  CASE (ana_type_default)
1408  file_name = trim(expand_file_name_char(file_name_tmp, &
1409  "diel_const"))
1410  file_name_tmp = trim(expand_file_name_char(file_name_tmp, &
1411  "diel_const_tensor"))
1412  CASE (ana_type_sym_xyz)
1413  file_name = trim(expand_file_name_char(file_name_tmp, &
1414  "diel_const_sym"))
1415  file_name_tmp = trim(expand_file_name_char(file_name_tmp, &
1416  "diel_const_tensor_sym"))
1417  CASE DEFAULT
1418  cpwarn('unknown analysis type "'//cp_to_string(ana_env%dip_ana%ana_type)//'" used.')
1419  END SELECT
1420 
1421  ! calc the dielectric constant
1422  ! 1+( <M^2> - <M>^2 ) / (3*e_0*V*k*T)
1423  diel_const = 1.0_dp + (ana_env%dip_ana%mu2_pv_s/(ana_env%dip_ana%conf_counter) - &
1424  dot_product(ana_env%dip_ana%mu_psv(:)/(ana_env%dip_ana%conf_counter), &
1425  ana_env%dip_ana%mu_psv(:)/(ana_env%dip_ana%conf_counter)))* &
1426  diel_const_norm
1427  ! symmetrized dielctric constant
1428  ! 1+( <M^2> ) / (3*e_0*V*k*T)
1429  diel_const_sym = 1.0_dp + ana_env%dip_ana%mu2_pv_s/(ana_env%dip_ana%conf_counter)* &
1430  diel_const_norm
1431  ! print dielectric constant trajectory
1432  ! if szmetry used print only every 8th configuration, hence every different (not mirrowed)
1433  INQUIRE (file=file_name, exist=flag)
1434  CALL open_file(file_name=file_name, file_status="UNKNOWN", &
1435  file_action="WRITE", file_position="APPEND", &
1436  unit_number=file_ptr)
1437  IF (.NOT. flag) THEN
1438  WRITE (file_ptr, fmt='(A8,5A20)') "# conf", "diel_const", &
1439  "diel_const_sym", "diel_const_sym_x", &
1440  "diel_const_sym_y", "diel_const_sym_z"
1441  END IF
1442  WRITE (file_ptr, fmt="(I8,10F20.10)") counter_tmp, diel_const, &
1443  diel_const_sym, &
1444  4.0_dp*pi/(kb*ana_env%temperature)* &
1445  ana_env%dip_ana%mu2_pv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1446  CALL close_file(unit_number=file_ptr)
1447 
1448  ! print dielectric constant tensor trajectory
1449  INQUIRE (file=file_name_tmp, exist=flag)
1450  CALL open_file(file_name=file_name_tmp, file_status="UNKNOWN", &
1451  file_action="WRITE", file_position="APPEND", &
1452  unit_number=file_ptr)
1453  IF (.NOT. flag) THEN
1454  WRITE (file_ptr, fmt='(A8,9A20)') "# conf", "xx", "xy", "xz", &
1455  "yx", "yy", "yz", &
1456  "zx", "zy", "zz"
1457  END IF
1458  tmp_dip(:, :) = 0.0_dp
1459  tmp_dip(:, 1) = ana_env%dip_ana%mu_psv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1460 
1461  WRITE (file_ptr, fmt="(I8,10F20.10)") counter_tmp, &
1462  4.0_dp*pi/(kb*ana_env%temperature)* &
1463  (ana_env%dip_ana%mu2_pv_mat(:, :)/real(ana_env%dip_ana%conf_counter, kind=dp) - &
1464  matmul(tmp_dip(:, :), transpose(tmp_dip(:, :))))
1465  CALL close_file(unit_number=file_ptr)
1466  END SUBROUTINE print_act_dipole_analysis
1467 
1468 ! **************************************************************************************************
1469 !> \brief prints the dipole moment analysis
1470 !> \param ana_env ...
1471 !> \param
1472 !> \author Mandes 03.2013
1473 ! **************************************************************************************************
1474  SUBROUTINE print_dipole_analysis(ana_env)
1475  TYPE(tmc_analysis_env), POINTER :: ana_env
1476 
1477  CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA"
1478 
1479  INTEGER :: i
1480  REAL(kind=dp) :: diel_const_scalar, kb
1481  REAL(kind=dp), DIMENSION(3) :: diel_const_sym, dielec_ev
1482  REAL(kind=dp), DIMENSION(3, 3) :: diel_const, tmp_dip, tmp_ev
1483 
1484  kb = boltzmann/joule
1485 
1486  cpassert(ASSOCIATED(ana_env))
1487  cpassert(ASSOCIATED(ana_env%dip_ana))
1488 
1489  tmp_dip(:, :) = 0.0_dp
1490  diel_const(:, :) = 0.0_dp
1491  diel_const_scalar = 0.0_dp
1492  diel_const_sym = 0.0_dp
1493 
1494  !dielectric constant
1495  tmp_dip(:, 1) = ana_env%dip_ana%mu_psv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1496  diel_const(:, :) = 4.0_dp*pi/(kb*ana_env%temperature)* &
1497  (ana_env%dip_ana%mu2_pv_mat(:, :)/real(ana_env%dip_ana%conf_counter, kind=dp) - &
1498  matmul(tmp_dip(:, :), transpose(tmp_dip(:, :))))
1499 
1500  !dielectric constant for symmetric case
1501  diel_const_sym(:) = 4.0_dp*pi/(kb*ana_env%temperature)* &
1502  ana_env%dip_ana%mu2_pv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1503 
1504  DO i = 1, 3
1505  diel_const(i, i) = diel_const(i, i) + 1.0_dp ! +1 for unpolarizable models, 1.592 for polarizable
1506  diel_const_scalar = diel_const_scalar + diel_const(i, i)
1507  END DO
1508  diel_const_scalar = diel_const_scalar/real(3, kind=dp)
1509 
1510  tmp_dip(:, :) = diel_const
1511  CALL diag(3, tmp_dip, dielec_ev, tmp_ev)
1512 
1513  ! print out results
1514  WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
1515  WRITE (ana_env%io_unit, fmt="(T2,A,T35,A,T80,A)") "-", "average dipoles", "-"
1516  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
1517  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "used configurations ", &
1518  cp_to_string(real(ana_env%dip_ana%conf_counter, kind=dp))
1519  IF (ana_env%dip_ana%ana_type .EQ. ana_type_ice) &
1520  WRITE (ana_env%io_unit, fmt='(T2,A,"| ",A)') plabel, &
1521  "ice analysis with directions of hexagonal structure"
1522  IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
1523  WRITE (ana_env%io_unit, fmt='(T2,A,"| ",A)') plabel, &
1524  "ice analysis with symmetrized dipoles in each direction."
1525 
1526  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "for product of 2 directions(per vol):"
1527  DO i = 1, 3
1528  WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", ana_env%dip_ana%mu2_pv_mat(i, :)/ &
1529  REAL(ana_env%dip_ana%conf_counter, kind=dp), " |"
1530  END DO
1531 
1532  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric constant tensor:"
1533  DO i = 1, 3
1534  WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", diel_const(i, :), " |"
1535  END DO
1536 
1537  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric tensor eigenvalues", &
1538  cp_to_string(dielec_ev(1))//" "// &
1539  cp_to_string(dielec_ev(2))//" "// &
1540  cp_to_string(dielec_ev(3))
1541  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric constant symm ", &
1542  cp_to_string(diel_const_sym(1))//" | "// &
1543  cp_to_string(diel_const_sym(2))//" | "// &
1544  cp_to_string(diel_const_sym(3))
1545  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric constant ", &
1546  cp_to_string(diel_const_scalar)
1547  WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
1548 
1549  END SUBROUTINE print_dipole_analysis
1550 
1551  !============================================================================
1552  ! particle displacement in cell (from one configuration to the next)
1553  !============================================================================
1554 
1555 ! **************************************************************************************************
1556 !> \brief calculates the mean square displacement
1557 !> \param elem ...
1558 !> \param ana_env ...
1559 !> \param
1560 !> \author Mandes 02.2013
1561 ! **************************************************************************************************
1562  SUBROUTINE calc_displacement(elem, ana_env)
1563  TYPE(tree_type), POINTER :: elem
1564  TYPE(tmc_analysis_env), POINTER :: ana_env
1565 
1566  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_displacement'
1567 
1568  CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp
1569  INTEGER :: file_ptr, handle, ind
1570  LOGICAL :: flag
1571  REAL(kind=dp) :: disp
1572  REAL(kind=dp), DIMENSION(3) :: atom_disp
1573 
1574  disp = 0.0_dp
1575 
1576  cpassert(ASSOCIATED(elem))
1577  cpassert(ASSOCIATED(elem%pos))
1578  cpassert(ASSOCIATED(ana_env))
1579  cpassert(ASSOCIATED(ana_env%displace))
1580  cpassert(ASSOCIATED(ana_env%last_elem))
1581 
1582  ! start the timing
1583  CALL timeset(routinen, handle)
1584 
1585  DO ind = 1, SIZE(elem%pos), ana_env%dim_per_elem
1586  ! fold into box
1587  atom_disp(:) = elem%pos(ind:ind + 2) - ana_env%last_elem%pos(ind:ind + 2)
1588  CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
1589  vec=atom_disp)
1590  disp = disp + sum((atom_disp(:)*au2a)**2)
1591  END DO
1592  ana_env%displace%disp = ana_env%displace%disp + disp
1593  ana_env%displace%conf_counter = ana_env%displace%conf_counter + 1
1594 
1595  IF (ana_env%displace%print_disp) THEN
1596  file_name_tmp = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1597  tmc_default_trajectory_file_name, &
1598  ana_env%temperature)
1599  file_name = trim(expand_file_name_char(file_name_tmp, &
1600  "devi"))
1601  INQUIRE (file=file_name, exist=flag)
1602  CALL open_file(file_name=file_name, file_status="UNKNOWN", &
1603  file_action="WRITE", file_position="APPEND", &
1604  unit_number=file_ptr)
1605  IF (.NOT. flag) &
1606  WRITE (file_ptr, *) "# conf squared deviation of the cell"
1607  WRITE (file_ptr, *) elem%nr, disp
1608  CALL close_file(unit_number=file_ptr)
1609  END IF
1610 
1611  ! end the timing
1612  CALL timestop(handle)
1613 
1614  END SUBROUTINE calc_displacement
1615 
1616 ! **************************************************************************************************
1617 !> \brief prints final values for the displacement calculations
1618 !> \param ana_env ...
1619 !> \param
1620 !> \author Mandes 02.2013
1621 ! **************************************************************************************************
1622  SUBROUTINE print_average_displacement(ana_env)
1623  TYPE(tmc_analysis_env), POINTER :: ana_env
1624 
1625  CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA"
1626 
1627  WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
1628  WRITE (ana_env%io_unit, fmt="(T2,A,T35,A,T80,A)") "-", "average displacement", "-"
1629  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "temperature ", &
1630  cp_to_string(ana_env%temperature)
1631  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "used configurations ", &
1632  cp_to_string(real(ana_env%displace%conf_counter, kind=dp))
1633  WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "cell root mean square deviation: ", &
1634  cp_to_string(sqrt(ana_env%displace%disp/ &
1635  REAL(ana_env%displace%conf_counter, kind=dp)))
1636  IF (ana_env%print_test_output) &
1637  WRITE (*, *) "TMC|ANALYSIS_AVERAGE_CELL_DISPLACEMENT_X= ", &
1638  sqrt(ana_env%displace%disp/ &
1639  REAL(ana_env%displace%conf_counter, kind=dp))
1640  END SUBROUTINE print_average_displacement
1641 END MODULE tmc_analysis
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
atom
Definition: atom.F:9
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
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_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
force_fields_input
Definition: force_fields_input.F:20
force_fields_input::read_chrg_section
subroutine, public read_chrg_section(charge_atm, charge, section, start)
Reads the CHARGE section.
Definition: force_fields_input.F:1842
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_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
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
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::diag
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition: mathlib.F:1515
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::boltzmann
real(kind=dp), parameter, public boltzmann
Definition: physcon.F:129
physcon::a_mass
real(kind=dp), parameter, public a_mass
Definition: physcon.F:132
physcon::joule
real(kind=dp), parameter, public joule
Definition: physcon.F:159
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
physcon::massunit
real(kind=dp), parameter, public massunit
Definition: physcon.F:141
tmc_analysis_types
module provides variables for the TMC analysis tool
Definition: tmc_analysis_types.F:15
tmc_analysis_types::search_pair_in_list
integer function, public search_pair_in_list(pair_list, n1, n2, list_end)
search the pair of two atom types in list
Definition: tmc_analysis_types.F:268
tmc_analysis_types::tmc_ana_displacement_create
subroutine, public tmc_ana_displacement_create(ana_disp, dim_per_elem)
creates a new structure environment for TMC analysis
Definition: tmc_analysis_types.F:402
tmc_analysis_types::tmc_ana_dipole_analysis_create
subroutine, public tmc_ana_dipole_analysis_create(ana_dip_ana)
creates a new structure environment for TMC analysis
Definition: tmc_analysis_types.F:356
tmc_analysis_types::tmc_ana_dipole_moment_create
subroutine, public tmc_ana_dipole_moment_create(ana_dip_mom, charge_atm, charge, dim_per_elem)
creates a new structure environment for TMC analysis
Definition: tmc_analysis_types.F:311
tmc_analysis_types::tmc_ana_env_create
subroutine, public tmc_ana_env_create(tmc_ana)
creates a new structure environment for TMC analysis
Definition: tmc_analysis_types.F:141
tmc_analysis_types::ana_type_default
integer, parameter, public ana_type_default
Definition: tmc_analysis_types.F:44
tmc_analysis_types::ana_type_ice
integer, parameter, public ana_type_ice
Definition: tmc_analysis_types.F:45
tmc_analysis_types::tmc_ana_pair_correl_file_name
character(len=default_path_length), parameter, public tmc_ana_pair_correl_file_name
Definition: tmc_analysis_types.F:41
tmc_analysis_types::tmc_ana_density_file_name
character(len=default_path_length), parameter, public tmc_ana_density_file_name
Definition: tmc_analysis_types.F:39
tmc_analysis_types::ana_type_sym_xyz
integer, parameter, public ana_type_sym_xyz
Definition: tmc_analysis_types.F:46
tmc_analysis_types::tmc_ana_density_create
subroutine, public tmc_ana_density_create(ana_dens, nr_bins)
creates a new structure environment for TMC analysis
Definition: tmc_analysis_types.F:191
tmc_analysis_types::tmc_ana_pair_correl_create
subroutine, public tmc_ana_pair_correl_create(ana_pair_correl, nr_bins)
creates a new structure environment for TMC analysis
Definition: tmc_analysis_types.F:233
tmc_analysis
module analyses element of the TMC tree element structure e.g. density, radial distribution function,...
Definition: tmc_analysis.F:16
tmc_analysis::analysis_restart_read
subroutine, public analysis_restart_read(ana_env, elem)
read analysis restart file
Definition: tmc_analysis.F:366
tmc_analysis::finalize_tmc_analysis
subroutine, public finalize_tmc_analysis(ana_env)
call all the necessarry analysis printing routines
Definition: tmc_analysis.F:587
tmc_analysis::tmc_read_ana_input
subroutine, public tmc_read_ana_input(tmc_ana_section, tmc_ana)
creates a new para environment for tmc analysis
Definition: tmc_analysis.F:86
tmc_analysis::analyze_file_configurations
subroutine, public analyze_file_configurations(start_id, end_id, dir_ind, ana_env, tmc_params)
read the files and analyze the configurations
Definition: tmc_analysis.F:633
tmc_analysis::analysis_init
subroutine, public analysis_init(ana_env, nr_dim)
initialize all the necessarry analysis structures
Definition: tmc_analysis.F:204
tmc_analysis::do_tmc_analysis
subroutine, public do_tmc_analysis(elem, ana_env)
call all the necessarry analysis routines analysis the previous element with the weight of the differ...
Definition: tmc_analysis.F:465
tmc_analysis::analysis_restart_print
subroutine, public analysis_restart_print(ana_env)
print analysis restart file
Definition: tmc_analysis.F:269
tmc_calculations
calculation section for TreeMonteCarlo
Definition: tmc_calculations.F:15
tmc_calculations::get_scaled_cell
subroutine, public get_scaled_cell(cell, box_scale, scaled_hmat, scaled_cell, vol, abc, vec)
handles properties and calculations of a scaled cell
Definition: tmc_calculations.F:153
tmc_calculations::nearest_distance
real(kind=dp) function, public nearest_distance(x1, x2, cell, box_scale)
neares distance of atoms within the periodic boundary condition
Definition: tmc_calculations.F:233
tmc_file_io
writing and printing the files, trajectory (pos, cell, dipoles) as well as restart files
Definition: tmc_file_io.F:20
tmc_file_io::analyse_files_close
subroutine, public analyse_files_close(tmc_ana)
close the files for reading configurations data to analyze
Definition: tmc_file_io.F:1043
tmc_file_io::write_dipoles_in_file
subroutine, public write_dipoles_in_file(file_name, conf_nr, dip, file_ext)
writes the cell dipoles in dipole trajectory file
Definition: tmc_file_io.F:573
tmc_file_io::read_element_from_file
subroutine, public read_element_from_file(elem, tmc_ana, conf_nr, stat)
read the trajectory element from a file from sub tree element
Definition: tmc_file_io.F:619
tmc_file_io::analyse_files_open
subroutine, public analyse_files_open(tmc_ana, stat, dir_ind)
opens the files for reading configurations data to analyze
Definition: tmc_file_io.F:943
tmc_file_io::expand_file_name_char
character(len=default_path_length) function, public expand_file_name_char(file_name, extra)
placing a character string at the end of a file name (before the file extension)
Definition: tmc_file_io.F:102
tmc_file_io::expand_file_name_temp
character(len=default_path_length) function, public expand_file_name_temp(file_name, rvalue)
placing the temperature at the end of a file name (before the file extension)
Definition: tmc_file_io.F:129
tmc_stati
tree nodes creation, searching, deallocation, references etc.
Definition: tmc_stati.F:15
tmc_stati::tmc_default_trajectory_file_name
character(len= *), parameter, public tmc_default_trajectory_file_name
Definition: tmc_stati.F:24
tmc_stati::tmc_default_unspecified_name
character(len= *), parameter, public tmc_default_unspecified_name
Definition: tmc_stati.F:40
tmc_stati::tmc_status_wait_for_new_task
integer, parameter, public tmc_status_wait_for_new_task
Definition: tmc_stati.F:52
tmc_stati::tmc_default_restart_in_file_name
character(len= *), parameter, public tmc_default_restart_in_file_name
Definition: tmc_stati.F:28
tmc_stati::tmc_default_restart_out_file_name
character(len= *), parameter, public tmc_default_restart_out_file_name
Definition: tmc_stati.F:26
tmc_stati::tmc_status_ok
integer, parameter, public tmc_status_ok
Definition: tmc_stati.F:51
tmc_tree_build
tree nodes creation, deallocation, references etc.
Definition: tmc_tree_build.F:35
tmc_tree_build::deallocate_sub_tree_node
subroutine, public deallocate_sub_tree_node(tree_elem)
deallocates an elements of the subtree element structure
Definition: tmc_tree_build.F:220
tmc_tree_build::allocate_new_sub_tree_node
subroutine, public allocate_new_sub_tree_node(tmc_params, next_el, nr_dim)
allocates an elements of the subtree element structure
Definition: tmc_tree_build.F:160
tmc_tree_types
module handles definition of the tree nodes for the global and the subtrees binary tree parent elemen...
Definition: tmc_tree_types.F:32
tmc_tree_types::read_subtree_elem_unformated
subroutine, public read_subtree_elem_unformated(elem, io_unit)
reads the TMC sub tree structure element unformated in file
Definition: tmc_tree_types.F:274
tmc_tree_types::write_subtree_elem_unformated
subroutine, public write_subtree_elem_unformated(elem, io_unit)
prints out the TMC sub tree structure element unformated in file
Definition: tmc_tree_types.F:244
tmc_types
module handles definition of the tree nodes for the global and the subtrees binary tree parent elemen...
Definition: tmc_types.F:32