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