(git:97501a3)
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-2025 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 .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
480 weight_act = 0
481 IF (ASSOCIATED(ana_env%last_elem)) THEN
482 weight_act = elem%nr - ana_env%last_elem%nr
483 END IF
484
485 IF (weight_act .GT. 0) THEN
486 ! calculates the 3 dimensional distributed density
487 IF (ASSOCIATED(ana_env%density_3d)) &
488 CALL calc_density_3d(elem=ana_env%last_elem, &
489 weight=weight_act, atoms=ana_env%atoms, &
490 ana_env=ana_env)
491 ! calculated the radial distribution function for each atom type
492 IF (ASSOCIATED(ana_env%pair_correl)) &
493 CALL calc_paircorrelation(elem=ana_env%last_elem, weight=weight_act, &
494 atoms=ana_env%atoms, ana_env=ana_env)
495 ! calculates the classical dipole moments
496 IF (ASSOCIATED(ana_env%dip_mom)) &
497 CALL calc_dipole_moment(elem=ana_env%last_elem, weight=weight_act, &
498 ana_env=ana_env)
499 ! calculates the dipole moments analysis and dielectric constant
500 IF (ASSOCIATED(ana_env%dip_ana)) THEN
501 ! in symmetric case use also the dipoles
502 ! (-x,y,z) .. .. (-x,-y,z).... (-x,-y-z) all have the same energy
503 IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) THEN
504 ! (-x,y,z)
505 ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
506 dip_tmp(:) = ana_env%last_elem%dipole(:)
507 IF (ASSOCIATED(ana_env%dip_mom)) &
508 ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
509 CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
510 ana_env=ana_env)
511 ! (-x,-y,z)
512 ana_env%last_elem%dipole(:) = dip_tmp(:)
513 ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
514 dip_tmp(:) = ana_env%last_elem%dipole(:)
515 IF (ASSOCIATED(ana_env%dip_mom)) &
516 ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
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(3) = -ana_env%last_elem%dipole(3)
522 dip_tmp(:) = ana_env%last_elem%dipole(:)
523 IF (ASSOCIATED(ana_env%dip_mom)) &
524 ana_env%dip_mom%last_dip_cl(3) = -ana_env%dip_mom%last_dip_cl(3)
525 CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
526 ana_env=ana_env)
527 ! (x,-y,-z)
528 ana_env%last_elem%dipole(:) = dip_tmp(:)
529 ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
530 dip_tmp(:) = ana_env%last_elem%dipole(:)
531 IF (ASSOCIATED(ana_env%dip_mom)) &
532 ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
533 CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
534 ana_env=ana_env)
535 ! (x,y,-z)
536 ana_env%last_elem%dipole(:) = dip_tmp(:)
537 ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
538 dip_tmp(:) = ana_env%last_elem%dipole(:)
539 IF (ASSOCIATED(ana_env%dip_mom)) &
540 ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
541 CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
542 ana_env=ana_env)
543 ! (-x,y,-z)
544 ana_env%last_elem%dipole(:) = dip_tmp(:)
545 ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
546 dip_tmp(:) = ana_env%last_elem%dipole(:)
547 IF (ASSOCIATED(ana_env%dip_mom)) &
548 ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
549 CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
550 ana_env=ana_env)
551 ! (x,-y,z)
552 ana_env%last_elem%dipole(:) = dip_tmp(:)
553 ana_env%last_elem%dipole(:) = -ana_env%last_elem%dipole(:)
554 dip_tmp(:) = ana_env%last_elem%dipole(:)
555 IF (ASSOCIATED(ana_env%dip_mom)) &
556 ana_env%dip_mom%last_dip_cl(:) = -ana_env%dip_mom%last_dip_cl(:)
557 CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
558 ana_env=ana_env)
559 ! back to (x,y,z)
560 ana_env%last_elem%dipole(:) = dip_tmp(:)
561 ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
562 dip_tmp(:) = ana_env%last_elem%dipole(:)
563 IF (ASSOCIATED(ana_env%dip_mom)) &
564 ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
565 END IF
566 CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
567 ana_env=ana_env)
568 CALL print_act_dipole_analysis(elem=ana_env%last_elem, &
569 ana_env=ana_env)
570 END IF
571
572 ! calculates the cell displacement from last cell
573 IF (ASSOCIATED(ana_env%displace)) THEN
574 CALL calc_displacement(elem=elem, ana_env=ana_env)
575 END IF
576 END IF
577 ! swap elem with last elem, to delete original last element and store the actual one
578 elem_tmp => ana_env%last_elem
579 ana_env%last_elem => elem
580 elem => elem_tmp
581 ! end the timing
582 CALL timestop(handle)
583 END SUBROUTINE do_tmc_analysis
584
585! **************************************************************************************************
586!> \brief call all the necessarry analysis printing routines
587!> \param ana_env ...
588!> \param
589!> \author Mandes 02.2013
590! **************************************************************************************************
591 SUBROUTINE finalize_tmc_analysis(ana_env)
592 TYPE(tmc_analysis_env), POINTER :: ana_env
593
594 CHARACTER(LEN=*), PARAMETER :: routinen = 'finalize_tmc_analysis'
595
596 INTEGER :: handle
597
598 cpassert(ASSOCIATED(ana_env))
599
600 ! start the timing
601 CALL timeset(routinen, handle)
602 IF (ASSOCIATED(ana_env%density_3d)) THEN
603 IF (ana_env%density_3d%conf_counter .GT. 0) &
604 CALL print_density_3d(ana_env=ana_env)
605 END IF
606 IF (ASSOCIATED(ana_env%pair_correl)) THEN
607 IF (ana_env%pair_correl%conf_counter .GT. 0) &
608 CALL print_paircorrelation(ana_env=ana_env)
609 END IF
610 IF (ASSOCIATED(ana_env%dip_mom)) THEN
611 IF (ana_env%dip_mom%conf_counter .GT. 0) &
612 CALL print_dipole_moment(ana_env)
613 END IF
614 IF (ASSOCIATED(ana_env%dip_ana)) THEN
615 IF (ana_env%dip_ana%conf_counter .GT. 0) &
616 CALL print_dipole_analysis(ana_env)
617 END IF
618 IF (ASSOCIATED(ana_env%displace)) THEN
619 IF (ana_env%displace%conf_counter .GT. 0) &
620 CALL print_average_displacement(ana_env)
621 END IF
622
623 ! end the timing
624 CALL timestop(handle)
625 END SUBROUTINE finalize_tmc_analysis
626
627! **************************************************************************************************
628!> \brief read the files and analyze the configurations
629!> \param start_id ...
630!> \param end_id ...
631!> \param dir_ind ...
632!> \param ana_env ...
633!> \param tmc_params ...
634!> \author Mandes 03.2013
635! **************************************************************************************************
636 SUBROUTINE analyze_file_configurations(start_id, end_id, dir_ind, &
637 ana_env, tmc_params)
638 INTEGER :: start_id, end_id
639 INTEGER, OPTIONAL :: dir_ind
640 TYPE(tmc_analysis_env), POINTER :: ana_env
641 TYPE(tmc_param_type), POINTER :: tmc_params
642
643 CHARACTER(LEN=*), PARAMETER :: routinen = 'analyze_file_configurations'
644
645 INTEGER :: conf_nr, handle, nr_dim, stat
646 TYPE(tree_type), POINTER :: elem
647
648 NULLIFY (elem)
649 conf_nr = -1
650 stat = tmc_status_wait_for_new_task
651 cpassert(ASSOCIATED(ana_env))
652 cpassert(ASSOCIATED(tmc_params))
653
654 ! start the timing
655 CALL timeset(routinen, handle)
656
657 ! open the files
658 CALL analyse_files_open(tmc_ana=ana_env, stat=stat, dir_ind=dir_ind)
659 ! set the existence of exact dipoles (from file)
660 IF (ana_env%id_dip .GT. 0) THEN
661 tmc_params%print_dipole = .true.
662 ELSE
663 tmc_params%print_dipole = .false.
664 END IF
665
666 ! allocate the actual element structure
667 CALL allocate_new_sub_tree_node(tmc_params=tmc_params, next_el=elem, &
668 nr_dim=ana_env%nr_dim)
669
670 IF (ASSOCIATED(ana_env%last_elem)) conf_nr = ana_env%last_elem%nr
671 nr_dim = SIZE(elem%pos)
672
673 IF (stat .EQ. tmc_status_ok) THEN
674 conf_loop: DO
675 CALL read_element_from_file(elem=elem, tmc_ana=ana_env, conf_nr=conf_nr, &
676 stat=stat)
677 IF (stat .EQ. tmc_status_wait_for_new_task) THEN
678 CALL deallocate_sub_tree_node(tree_elem=elem)
679 EXIT conf_loop
680 END IF
681 ! if we want just a certain part of the trajectory
682 IF (start_id .LT. 0 .OR. conf_nr .GE. start_id) THEN
683 IF (end_id .LT. 0 .OR. conf_nr .LE. end_id) THEN
684 ! do the analysis calculations
685 CALL do_tmc_analysis(elem=elem, ana_env=ana_env)
686 END IF
687 END IF
688
689 ! clean temporary element (already analyzed)
690 IF (ASSOCIATED(elem)) THEN
691 CALL deallocate_sub_tree_node(tree_elem=elem)
692 END IF
693 ! if there was no previous element, create a new temp element to write in
694 IF (.NOT. ASSOCIATED(elem)) &
695 CALL allocate_new_sub_tree_node(tmc_params=tmc_params, next_el=elem, &
696 nr_dim=nr_dim)
697 END DO conf_loop
698 END IF
699 ! close the files
700 CALL analyse_files_close(tmc_ana=ana_env)
701
702 IF (ASSOCIATED(elem)) &
703 CALL deallocate_sub_tree_node(tree_elem=elem)
704
705 ! end the timing
706 CALL timestop(handle)
707 END SUBROUTINE analyze_file_configurations
708
709 !============================================================================
710 ! density calculations
711 !============================================================================
712
713! **************************************************************************************************
714!> \brief calculates the density in rectantangulares
715!> defined by the number of bins in each direction
716!> \param elem ...
717!> \param weight ...
718!> \param atoms ...
719!> \param ana_env ...
720!> \param
721!> \author Mandes 02.2013
722! **************************************************************************************************
723 SUBROUTINE calc_density_3d(elem, weight, atoms, ana_env)
724 TYPE(tree_type), POINTER :: elem
725 INTEGER :: weight
726 TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
727 TYPE(tmc_analysis_env), POINTER :: ana_env
728
729 CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_density_3d'
730
731 CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp
732 INTEGER :: atom, bin_x, bin_y, bin_z, file_ptr, &
733 handle
734 LOGICAL :: flag
735 REAL(kind=dp) :: mass_total, r_tmp, vol_cell, vol_sub_box
736 REAL(kind=dp), DIMENSION(3) :: atom_pos, cell_size, interval_size
737 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mass_bin
738
739 NULLIFY (mass_bin)
740
741 cpassert(ASSOCIATED(elem))
742 cpassert(ASSOCIATED(elem%pos))
743 cpassert(weight .GT. 0)
744 cpassert(ASSOCIATED(atoms))
745 cpassert(ASSOCIATED(ana_env))
746 cpassert(ASSOCIATED(ana_env%cell))
747 cpassert(ASSOCIATED(ana_env%density_3d))
748 cpassert(ASSOCIATED(ana_env%density_3d%sum_density))
749 cpassert(ASSOCIATED(ana_env%density_3d%sum_dens2))
750
751 ! start the timing
752 CALL timeset(routinen, handle)
753
754 atom_pos(:) = 0.0_dp
755 cell_size(:) = 0.0_dp
756 interval_size(:) = 0.0_dp
757 mass_total = 0.0_dp
758
759 bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
760 bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
761 bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
762 ALLOCATE (mass_bin(bin_x, bin_y, bin_z))
763 mass_bin(:, :, :) = 0.0_dp
764
765 ! if NPT -> box_scale.NE.1.0 use the scaled cell
766 ! ATTENTION then the sub box middle points are not correct in the output
767 ! espacially if we use multiple sub boxes
768 CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
769 abc=cell_size, vol=vol_cell)
770 ! volume summed over configurations for average volume [A]
771 ana_env%density_3d%sum_vol = ana_env%density_3d%sum_vol + &
772 vol_cell*(au2a)**3*weight
773 ana_env%density_3d%sum_vol2 = ana_env%density_3d%sum_vol2 + &
774 (vol_cell*(au2a)**3)**2*weight
775
776 ana_env%density_3d%sum_box_length(:) = ana_env%density_3d%sum_box_length(:) &
777 + cell_size(:)*(au2a)*weight
778 ana_env%density_3d%sum_box_length2(:) = ana_env%density_3d%sum_box_length2(:) &
779 + (cell_size(:)*(au2a))**2*weight
780
781 ! sub interval length
782 interval_size(1) = cell_size(1)/real(bin_x, dp)
783 interval_size(2) = cell_size(2)/real(bin_y, dp)
784 interval_size(3) = cell_size(3)/real(bin_z, dp)
785
786 ! volume in [cm^3]
787 vol_cell = vol_cell*(au2a*1e-8)**3
788 vol_sub_box = interval_size(1)*interval_size(2)*interval_size(3)* &
789 (au2a*1e-8)**3
790
791 ! count every atom
792 DO atom = 1, SIZE(elem%pos), ana_env%dim_per_elem
793
794 atom_pos(:) = elem%pos(atom:atom + 2)
795 ! fold into box
796 CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
797 vec=atom_pos)
798 ! shifts the box to positive values (before 0,0,0 is the center)
799 atom_pos(:) = atom_pos(:) + 0.5_dp*cell_size(:)
800 ! calculate the index of the sub box
801 bin_x = int(atom_pos(1)/interval_size(1)) + 1
802 bin_y = int(atom_pos(2)/interval_size(2)) + 1
803 bin_z = int(atom_pos(3)/interval_size(3)) + 1
804 cpassert(bin_x .GT. 0 .AND. bin_y .GT. 0 .AND. bin_z .GT. 0)
805 cpassert(bin_x .LE. SIZE(ana_env%density_3d%sum_density(:, 1, 1)))
806 cpassert(bin_y .LE. SIZE(ana_env%density_3d%sum_density(1, :, 1)))
807 cpassert(bin_z .LE. SIZE(ana_env%density_3d%sum_density(1, 1, :)))
808
809 ! sum mass in [g] (in bins and total)
810 mass_bin(bin_x, bin_y, bin_z) = mass_bin(bin_x, bin_y, bin_z) + &
811 atoms(int(atom/real(ana_env%dim_per_elem, kind=dp)) + 1)%mass/massunit*1000*a_mass
812 mass_total = mass_total + &
813 atoms(int(atom/real(ana_env%dim_per_elem, kind=dp)) + 1)%mass/massunit*1000*a_mass
814 !mass_bin(bin_x,bin_y,bin_z) = mass_bin(bin_x,bin_y,bin_z) + &
815 ! atoms(INT(atom/REAL(ana_env%dim_per_elem,KIND=dp))+1)%mass/&
816 ! massunit/n_avogadro
817 !mass_total = mass_total + &
818 ! atoms(INT(atom/REAL(ana_env%dim_per_elem,KIND=dp))+1)%mass/&
819 ! massunit/n_avogadro
820 END DO
821 ! check total cell density
822 r_tmp = mass_total/vol_cell - sum(mass_bin(:, :, :))/vol_sub_box/SIZE(mass_bin(:, :, :))
823 cpassert(abs(r_tmp) .LT. 1e-5)
824
825 ! calculate density (mass per volume) and sum up for average value
826 ana_env%density_3d%sum_density(:, :, :) = ana_env%density_3d%sum_density(:, :, :) + &
827 weight*mass_bin(:, :, :)/vol_sub_box
828
829 ! calculate density squared ( (mass per volume)^2 ) for variance and sum up for average value
830 ana_env%density_3d%sum_dens2(:, :, :) = ana_env%density_3d%sum_dens2(:, :, :) + &
831 weight*(mass_bin(:, :, :)/vol_sub_box)**2
832
833 ana_env%density_3d%conf_counter = ana_env%density_3d%conf_counter + weight
834
835 ! print out the actual and average density in file
836 IF (ana_env%density_3d%print_dens) THEN
837 file_name_tmp = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
838 tmc_default_trajectory_file_name, &
839 ana_env%temperature)
840 file_name = trim(expand_file_name_char(file_name_tmp, &
841 "dens"))
842 INQUIRE (file=file_name, exist=flag)
843 CALL open_file(file_name=file_name, file_status="UNKNOWN", &
844 file_action="WRITE", file_position="APPEND", &
845 unit_number=file_ptr)
846 IF (.NOT. flag) &
847 WRITE (file_ptr, fmt='(A8,11A20)') "# conf_nr", "dens_act[g/cm^3]", &
848 "dens_average[g/cm^3]", "density_variance", &
849 "averages:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z", &
850 "variances:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z"
851 WRITE (file_ptr, fmt="(I8,11F20.10)") ana_env%density_3d%conf_counter + 1 - weight, &
852 sum(mass_bin(:, :, :))/vol_sub_box/SIZE(mass_bin(:, :, :)), &
853 sum(ana_env%density_3d%sum_density(:, :, :))/ &
854 SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
855 REAL(ana_env%density_3d%conf_counter, kind=dp), &
856 sum(ana_env%density_3d%sum_dens2(:, :, :))/ &
857 SIZE(ana_env%density_3d%sum_dens2(:, :, :))/ &
858 REAL(ana_env%density_3d%conf_counter, kind=dp) - &
859 (sum(ana_env%density_3d%sum_density(:, :, :))/ &
860 SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
861 REAL(ana_env%density_3d%conf_counter, kind=dp))**2, &
862 ana_env%density_3d%sum_vol/ &
863 REAL(ana_env%density_3d%conf_counter, kind=dp), &
864 ana_env%density_3d%sum_box_length(:)/ &
865 REAL(ana_env%density_3d%conf_counter, kind=dp), &
866 ana_env%density_3d%sum_vol2/ &
867 REAL(ana_env%density_3d%conf_counter, kind=dp) - &
868 (ana_env%density_3d%sum_vol/ &
869 REAL(ana_env%density_3d%conf_counter, kind=dp))**2, &
870 ana_env%density_3d%sum_box_length2(:)/ &
871 REAL(ana_env%density_3d%conf_counter, kind=dp) - &
872 (ana_env%density_3d%sum_box_length(:)/ &
873 REAL(ana_env%density_3d%conf_counter, kind=dp))**2
874 CALL close_file(unit_number=file_ptr)
875 END IF
876
877 DEALLOCATE (mass_bin)
878 ! end the timing
879 CALL timestop(handle)
880 END SUBROUTINE calc_density_3d
881
882! **************************************************************************************************
883!> \brief print the density in rectantangulares
884!> defined by the number of bins in each direction
885!> \param ana_env ...
886!> \param
887!> \author Mandes 02.2013
888! **************************************************************************************************
889 SUBROUTINE print_density_3d(ana_env)
890 TYPE(tmc_analysis_env), POINTER :: ana_env
891
892 CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA", &
893 routinen = 'print_density_3d'
894
895 CHARACTER(LEN=default_path_length) :: file_name, file_name_vari
896 INTEGER :: bin_x, bin_y, bin_z, file_ptr_dens, &
897 file_ptr_vari, handle, i, j, k
898 REAL(kind=dp), DIMENSION(3) :: cell_size, interval_size
899
900 cpassert(ASSOCIATED(ana_env))
901 cpassert(ASSOCIATED(ana_env%density_3d))
902 cpassert(ASSOCIATED(ana_env%density_3d%sum_density))
903 cpassert(ASSOCIATED(ana_env%density_3d%sum_dens2))
904
905 ! start the timing
906 CALL timeset(routinen, handle)
907
908 file_name = ""
909 file_name_vari = ""
910
911 bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
912 bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
913 bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
914 CALL get_cell(cell=ana_env%cell, abc=cell_size)
915 interval_size(1) = cell_size(1)/real(bin_x, kind=dp)*au2a
916 interval_size(2) = cell_size(2)/real(bin_y, kind=dp)*au2a
917 interval_size(3) = cell_size(3)/real(bin_z, kind=dp)*au2a
918
919 file_name = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
920 tmc_ana_density_file_name, &
921 ana_env%temperature)
922 CALL open_file(file_name=file_name, file_status="REPLACE", &
923 file_action="WRITE", file_position="APPEND", &
924 unit_number=file_ptr_dens)
925 WRITE (file_ptr_dens, fmt='(A,1X,I0,1X,A,3(I0,1X),1X,A,1X,3F10.5)') &
926 "# configurations", ana_env%density_3d%conf_counter, "bins", &
927 ana_env%density_3d%nr_bins, "interval size", interval_size(:)
928 WRITE (file_ptr_dens, fmt='(A,3A10,A20)') "#", " x [A] ", " y [A] ", " z [A] ", " density [g/cm^3] "
929
930 file_name_vari = expand_file_name_temp(expand_file_name_char( &
931 trim(ana_env%out_file_prefix)// &
932 tmc_ana_density_file_name, "vari"), &
933 ana_env%temperature)
934 CALL open_file(file_name=file_name_vari, file_status="REPLACE", &
935 file_action="WRITE", file_position="APPEND", &
936 unit_number=file_ptr_vari)
937 WRITE (file_ptr_vari, fmt='(A,1X,I0,1X,A,3(I0,1X),1X,A,1X,3F10.5)') &
938 "# configurations", ana_env%density_3d%conf_counter, "bins", &
939 ana_env%density_3d%nr_bins, "interval size", interval_size(:)
940 WRITE (file_ptr_vari, fmt='(A,3A10,A20)') "#", " x [A] ", " y [A] ", " z [A] ", " variance"
941
942 DO i = 1, SIZE(ana_env%density_3d%sum_density(:, 1, 1))
943 DO j = 1, SIZE(ana_env%density_3d%sum_density(1, :, 1))
944 DO k = 1, SIZE(ana_env%density_3d%sum_density(1, 1, :))
945 WRITE (file_ptr_dens, fmt='(3F10.2,F20.10)') &
946 (i - 0.5_dp)*interval_size(1), (j - 0.5_dp)*interval_size(2), (k - 0.5_dp)*interval_size(3), &
947 ana_env%density_3d%sum_density(i, j, k)/real(ana_env%density_3d%conf_counter, kind=dp)
948 WRITE (file_ptr_vari, fmt='(3F10.2,F20.10)') &
949 (i - 0.5_dp)*interval_size(1), (j - 0.5_dp)*interval_size(2), (k - 0.5_dp)*interval_size(3), &
950 ana_env%density_3d%sum_dens2(i, j, k)/real(ana_env%density_3d%conf_counter, kind=dp) - &
951 (ana_env%density_3d%sum_density(i, j, k)/real(ana_env%density_3d%conf_counter, kind=dp))**2
952 END DO
953 END DO
954 END DO
955 CALL close_file(unit_number=file_ptr_dens)
956 CALL close_file(unit_number=file_ptr_vari)
957
958 WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
959 WRITE (ana_env%io_unit, fmt="(T2,A,T35,A,T80,A)") "-", "density calculation", "-"
960 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
961 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "used configurations", &
962 cp_to_string(real(ana_env%density_3d%conf_counter, kind=dp))
963 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "average volume", &
964 cp_to_string(ana_env%density_3d%sum_vol/ &
965 REAL(ana_env%density_3d%conf_counter, kind=dp))
966 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "average density in the cell: ", &
967 cp_to_string(sum(ana_env%density_3d%sum_density(:, :, :))/ &
968 SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
969 REAL(ana_env%density_3d%conf_counter, kind=dp))
970 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "density variance:", &
971 cp_to_string(sum(ana_env%density_3d%sum_dens2(:, :, :))/ &
972 SIZE(ana_env%density_3d%sum_dens2(:, :, :))/ &
973 REAL(ana_env%density_3d%conf_counter, kind=dp) - &
974 (sum(ana_env%density_3d%sum_density(:, :, :))/ &
975 SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
976 REAL(ana_env%density_3d%conf_counter, kind=dp))**2)
977 WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
978 IF (ana_env%print_test_output) &
979 WRITE (ana_env%io_unit, *) "TMC|ANALYSIS_CELL_DENSITY_X= ", &
980 sum(ana_env%density_3d%sum_density(:, :, :))/ &
981 SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
982 REAL(ana_env%density_3d%conf_counter, kind=dp)
983 ! end the timing
984 CALL timestop(handle)
985 END SUBROUTINE print_density_3d
986
987 !============================================================================
988 ! radial distribution function
989 !============================================================================
990
991! **************************************************************************************************
992!> \brief init radial distribution function structures
993!> \param ana_pair_correl ...
994!> \param atoms ...
995!> \param cell ...
996!> \param
997!> \author Mandes 02.2013
998! **************************************************************************************************
999 SUBROUTINE ana_pair_correl_init(ana_pair_correl, atoms, cell)
1000 TYPE(pair_correl_type), POINTER :: ana_pair_correl
1001 TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
1002 TYPE(cell_type), POINTER :: cell
1003
1004 CHARACTER(LEN=*), PARAMETER :: routinen = 'ana_pair_correl_init'
1005
1006 INTEGER :: counter, f_n, handle, list, list_ind, s_n
1007 REAL(kind=dp), DIMENSION(3) :: cell_size
1008 TYPE(atom_pairs_type), DIMENSION(:), POINTER :: pairs_tmp
1009
1010 cpassert(ASSOCIATED(ana_pair_correl))
1011 cpassert(.NOT. ASSOCIATED(ana_pair_correl%g_r))
1012 cpassert(.NOT. ASSOCIATED(ana_pair_correl%pairs))
1013 cpassert(ASSOCIATED(atoms))
1014 cpassert(SIZE(atoms) .GT. 1)
1015 cpassert(ASSOCIATED(cell))
1016
1017 ! start the timing
1018 CALL timeset(routinen, handle)
1019
1020 CALL get_cell(cell=cell, abc=cell_size)
1021 IF (ana_pair_correl%nr_bins .LE. 0) THEN
1022 ana_pair_correl%nr_bins = ceiling(maxval(cell_size(:))/2.0_dp/(0.03/au2a))
1023 END IF
1024 ana_pair_correl%step_length = maxval(cell_size(:))/2.0_dp/ &
1025 ana_pair_correl%nr_bins
1026 ana_pair_correl%conf_counter = 0
1027
1028 counter = 1
1029 ! initialise the atom pairs
1030 ALLOCATE (pairs_tmp(SIZE(atoms)))
1031 DO f_n = 1, SIZE(atoms)
1032 DO s_n = f_n + 1, SIZE(atoms)
1033 ! search if atom pair is already selected
1034 list_ind = search_pair_in_list(pair_list=pairs_tmp, n1=atoms(f_n)%name, &
1035 n2=atoms(s_n)%name, list_end=counter - 1)
1036 ! add to list
1037 IF (list_ind .LT. 0) THEN
1038 pairs_tmp(counter)%f_n = atoms(f_n)%name
1039 pairs_tmp(counter)%s_n = atoms(s_n)%name
1040 pairs_tmp(counter)%pair_count = 1
1041 counter = counter + 1
1042 ELSE
1043 pairs_tmp(list_ind)%pair_count = pairs_tmp(list_ind)%pair_count + 1
1044 END IF
1045 END DO
1046 END DO
1047
1048 ALLOCATE (ana_pair_correl%pairs(counter - 1))
1049 DO list = 1, counter - 1
1050 ana_pair_correl%pairs(list)%f_n = pairs_tmp(list)%f_n
1051 ana_pair_correl%pairs(list)%s_n = pairs_tmp(list)%s_n
1052 ana_pair_correl%pairs(list)%pair_count = pairs_tmp(list)%pair_count
1053 END DO
1054 DEALLOCATE (pairs_tmp)
1055
1056 ALLOCATE (ana_pair_correl%g_r(SIZE(ana_pair_correl%pairs(:)), ana_pair_correl%nr_bins))
1057 ana_pair_correl%g_r = 0.0_dp
1058 ! end the timing
1059 CALL timestop(handle)
1060 END SUBROUTINE ana_pair_correl_init
1061
1062! **************************************************************************************************
1063!> \brief calculates the radial distribution function
1064!> \param elem ...
1065!> \param weight ...
1066!> \param atoms ...
1067!> \param ana_env ...
1068!> \param
1069!> \author Mandes 02.2013
1070! **************************************************************************************************
1071 SUBROUTINE calc_paircorrelation(elem, weight, atoms, ana_env)
1072 TYPE(tree_type), POINTER :: elem
1073 INTEGER :: weight
1074 TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
1075 TYPE(tmc_analysis_env), POINTER :: ana_env
1076
1077 CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_paircorrelation'
1078
1079 INTEGER :: handle, i, ind, j, pair_ind
1080 REAL(kind=dp) :: dist
1081 REAL(kind=dp), DIMENSION(3) :: cell_size
1082
1083 cpassert(ASSOCIATED(elem))
1084 cpassert(ASSOCIATED(elem%pos))
1085 cpassert(all(elem%box_scale(:) .GT. 0.0_dp))
1086 cpassert(weight .GT. 0)
1087 cpassert(ASSOCIATED(atoms))
1088 cpassert(ASSOCIATED(ana_env))
1089 cpassert(ASSOCIATED(ana_env%cell))
1090 cpassert(ASSOCIATED(ana_env%pair_correl))
1091 cpassert(ASSOCIATED(ana_env%pair_correl%g_r))
1092 cpassert(ASSOCIATED(ana_env%pair_correl%pairs))
1093
1094 ! start the timing
1095 CALL timeset(routinen, handle)
1096
1097 dist = -1.0_dp
1098
1099 first_elem_loop: DO i = 1, SIZE(elem%pos), ana_env%dim_per_elem
1100 second_elem_loop: DO j = i + 3, SIZE(elem%pos), ana_env%dim_per_elem
1101 dist = nearest_distance(x1=elem%pos(i:i + ana_env%dim_per_elem - 1), &
1102 x2=elem%pos(j:j + ana_env%dim_per_elem - 1), &
1103 cell=ana_env%cell, box_scale=elem%box_scale)
1104 ind = ceiling(dist/ana_env%pair_correl%step_length)
1105 IF (ind .LE. ana_env%pair_correl%nr_bins) THEN
1106 pair_ind = search_pair_in_list(pair_list=ana_env%pair_correl%pairs, &
1107 n1=atoms(int(i/real(ana_env%dim_per_elem, kind=dp)) + 1)%name, &
1108 n2=atoms(int(j/real(ana_env%dim_per_elem, kind=dp)) + 1)%name)
1109 cpassert(pair_ind .GT. 0)
1110 ana_env%pair_correl%g_r(pair_ind, ind) = &
1111 ana_env%pair_correl%g_r(pair_ind, ind) + weight
1112 END IF
1113 END DO second_elem_loop
1114 END DO first_elem_loop
1115 ana_env%pair_correl%conf_counter = ana_env%pair_correl%conf_counter + weight
1116 CALL get_cell(cell=ana_env%cell, abc=cell_size)
1117 ana_env%pair_correl%sum_box_scale = ana_env%pair_correl%sum_box_scale + &
1118 (elem%box_scale(:)*weight)
1119 ! end the timing
1120 CALL timestop(handle)
1121 END SUBROUTINE calc_paircorrelation
1122
1123! **************************************************************************************************
1124!> \brief print the radial distribution function for each pair of atoms
1125!> \param ana_env ...
1126!> \param
1127!> \author Mandes 02.2013
1128! **************************************************************************************************
1129 SUBROUTINE print_paircorrelation(ana_env)
1130 TYPE(tmc_analysis_env), POINTER :: ana_env
1131
1132 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_paircorrelation'
1133
1134 CHARACTER(LEN=default_path_length) :: file_name
1135 INTEGER :: bin, file_ptr, handle, pair
1136 REAL(kind=dp) :: aver_box_scale(3), vol, voldr
1137 REAL(kind=dp), DIMENSION(3) :: cell_size
1138
1139 cpassert(ASSOCIATED(ana_env))
1140 cpassert(ASSOCIATED(ana_env%pair_correl))
1141
1142 ! start the timing
1143 CALL timeset(routinen, handle)
1144
1145 CALL get_cell(cell=ana_env%cell, abc=cell_size)
1146 aver_box_scale(:) = ana_env%pair_correl%sum_box_scale(:)/ana_env%pair_correl%conf_counter
1147 vol = (cell_size(1)*aver_box_scale(1))* &
1148 (cell_size(2)*aver_box_scale(2))* &
1149 (cell_size(3)*aver_box_scale(3))
1150
1151 DO pair = 1, SIZE(ana_env%pair_correl%pairs)
1152 file_name = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1153 tmc_ana_pair_correl_file_name, &
1154 ana_env%temperature)
1155 CALL open_file(file_name=expand_file_name_char( &
1156 expand_file_name_char(file_name, &
1157 ana_env%pair_correl%pairs(pair)%f_n), &
1158 ana_env%pair_correl%pairs(pair)%s_n), &
1159 file_status="REPLACE", &
1160 file_action="WRITE", file_position="APPEND", &
1161 unit_number=file_ptr)
1162 WRITE (file_ptr, *) "# radial distribution function of "// &
1163 trim(ana_env%pair_correl%pairs(pair)%f_n)//" and "// &
1164 trim(ana_env%pair_correl%pairs(pair)%s_n)//" of ", &
1165 ana_env%pair_correl%conf_counter, " configurations"
1166 WRITE (file_ptr, *) "# using a bin size of ", &
1167 ana_env%pair_correl%step_length*au2a, &
1168 "[A] (for Vol changes: referring to the reference cell)"
1169 DO bin = 1, ana_env%pair_correl%nr_bins
1170 voldr = 4.0/3.0*pi*ana_env%pair_correl%step_length**3* &
1171 (real(bin, kind=dp)**3 - real(bin - 1, kind=dp)**3)
1172 WRITE (file_ptr, *) (bin - 0.5)*ana_env%pair_correl%step_length*au2a, &
1173 (ana_env%pair_correl%g_r(pair, bin)/ana_env%pair_correl%conf_counter)/ &
1174 (voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
1175 END DO
1176 CALL close_file(unit_number=file_ptr)
1177
1178 IF (ana_env%print_test_output) &
1179 WRITE (*, *) "TMC|ANALYSIS_G_R_"// &
1180 trim(ana_env%pair_correl%pairs(pair)%f_n)//"_"// &
1181 trim(ana_env%pair_correl%pairs(pair)%s_n)//"_X= ", &
1182 sum(ana_env%pair_correl%g_r(pair, :)/ana_env%pair_correl%conf_counter/ &
1183 voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
1184 END DO
1185
1186 ! end the timing
1187 CALL timestop(handle)
1188 END SUBROUTINE print_paircorrelation
1189
1190 !============================================================================
1191 ! classical cell dipole moment
1192 !============================================================================
1193
1194! **************************************************************************************************
1195!> \brief init radial distribution function structures
1196!> \param ana_dip_mom ...
1197!> \param atoms ...
1198!> \param
1199!> \author Mandes 02.2013
1200! **************************************************************************************************
1201 SUBROUTINE ana_dipole_moment_init(ana_dip_mom, atoms)
1202 TYPE(dipole_moment_type), POINTER :: ana_dip_mom
1203 TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
1204
1205 CHARACTER(LEN=*), PARAMETER :: routinen = 'ana_dipole_moment_init'
1206
1207 INTEGER :: atom, charge, handle
1208
1209 cpassert(ASSOCIATED(ana_dip_mom))
1210 cpassert(ASSOCIATED(ana_dip_mom%charges_inp))
1211 cpassert(ASSOCIATED(atoms))
1212
1213 ! start the timing
1214 CALL timeset(routinen, handle)
1215
1216 ALLOCATE (ana_dip_mom%charges(SIZE(atoms)))
1217 ana_dip_mom%charges = 0.0_dp
1218 ! for every atom searcht the correct charge
1219 DO atom = 1, SIZE(atoms)
1220 charge_loop: DO charge = 1, SIZE(ana_dip_mom%charges_inp)
1221 IF (atoms(atom)%name .EQ. ana_dip_mom%charges_inp(charge)%name) THEN
1222 ana_dip_mom%charges(atom) = ana_dip_mom%charges_inp(charge)%mass
1223 EXIT charge_loop
1224 END IF
1225 END DO charge_loop
1226 END DO
1227
1228 DEALLOCATE (ana_dip_mom%charges_inp)
1229 ! end the timing
1230 CALL timestop(handle)
1231 END SUBROUTINE ana_dipole_moment_init
1232
1233! **************************************************************************************************
1234!> \brief calculates the classical cell dipole moment
1235!> \param elem ...
1236!> \param weight ...
1237!> \param ana_env ...
1238!> \param
1239!> \author Mandes 02.2013
1240! **************************************************************************************************
1241 SUBROUTINE calc_dipole_moment(elem, weight, ana_env)
1242 TYPE(tree_type), POINTER :: elem
1243 INTEGER :: weight
1244 TYPE(tmc_analysis_env), POINTER :: ana_env
1245
1246 CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_dipole_moment'
1247
1248 CHARACTER(LEN=default_path_length) :: file_name
1249 INTEGER :: handle, i
1250 REAL(kind=dp), DIMENSION(:), POINTER :: dip_cl
1251
1252 cpassert(ASSOCIATED(elem))
1253 cpassert(ASSOCIATED(elem%pos))
1254 cpassert(ASSOCIATED(ana_env))
1255 cpassert(ASSOCIATED(ana_env%dip_mom))
1256 cpassert(ASSOCIATED(ana_env%dip_mom%charges))
1257
1258 ! start the timing
1259 CALL timeset(routinen, handle)
1260
1261 ALLOCATE (dip_cl(ana_env%dim_per_elem))
1262 dip_cl(:) = 0.0_dp
1263
1264 DO i = 1, SIZE(elem%pos, 1), ana_env%dim_per_elem
1265 dip_cl(:) = dip_cl(:) + elem%pos(i:i + ana_env%dim_per_elem - 1)* &
1266 ana_env%dip_mom%charges(int(i/real(ana_env%dim_per_elem, kind=dp)) + 1)
1267 END DO
1268
1269 ! if there are no exact dipoles save these ones in element structure
1270 IF (.NOT. ASSOCIATED(elem%dipole)) THEN
1271 ALLOCATE (elem%dipole(ana_env%dim_per_elem))
1272 elem%dipole(:) = dip_cl(:)
1273 END IF
1274
1275 IF (ana_env%dip_mom%print_cl_dip) THEN
1276 file_name = expand_file_name_temp(tmc_default_trajectory_file_name, &
1277 ana_env%temperature)
1278 CALL write_dipoles_in_file(file_name=file_name, &
1279 conf_nr=ana_env%dip_mom%conf_counter + 1, dip=dip_cl, &
1280 file_ext="dip_cl")
1281 END IF
1282 ana_env%dip_mom%conf_counter = ana_env%dip_mom%conf_counter + weight
1283 ana_env%dip_mom%last_dip_cl(:) = dip_cl
1284
1285 DEALLOCATE (dip_cl)
1286
1287 ! end the timing
1288 CALL timestop(handle)
1289 END SUBROUTINE calc_dipole_moment
1290
1291! **************************************************************************************************
1292!> \brief prints final values for classical cell dipole moment calculation
1293!> \param ana_env ...
1294!> \param
1295!> \author Mandes 02.2013
1296! **************************************************************************************************
1297 SUBROUTINE print_dipole_moment(ana_env)
1298 TYPE(tmc_analysis_env), POINTER :: ana_env
1299
1300 IF (ana_env%print_test_output) &
1301 WRITE (*, *) "TMC|ANALYSIS_FINAL_CLASS_CELL_DIPOLE_MOMENT_X= ", &
1302 ana_env%dip_mom%last_dip_cl(:)
1303 END SUBROUTINE print_dipole_moment
1304
1305! **************************************************************************************************
1306!> \brief calculates the dipole moment analysis
1307!> \param elem ...
1308!> \param weight ...
1309!> \param ana_env ...
1310!> \param
1311!> \author Mandes 03.2013
1312! **************************************************************************************************
1313 SUBROUTINE calc_dipole_analysis(elem, weight, ana_env)
1314 TYPE(tree_type), POINTER :: elem
1315 INTEGER :: weight
1316 TYPE(tmc_analysis_env), POINTER :: ana_env
1317
1318 REAL(kind=dp) :: vol, weight_act
1319 REAL(kind=dp), DIMENSION(3, 3) :: tmp_dip
1320 TYPE(cell_type), POINTER :: scaled_cell
1321
1322 NULLIFY (scaled_cell)
1323
1324 cpassert(ASSOCIATED(elem))
1325 cpassert(ASSOCIATED(elem%dipole))
1326 cpassert(ASSOCIATED(ana_env))
1327 cpassert(ASSOCIATED(ana_env%dip_ana))
1328
1329 weight_act = weight
1330 IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
1331 weight_act = weight_act/real(8.0, kind=dp)
1332
1333 ! get the volume
1334 ALLOCATE (scaled_cell)
1335 CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, vol=vol, &
1336 scaled_cell=scaled_cell)
1337
1338 ! fold exact dipole moments using the classical ones
1339 IF (ASSOCIATED(ana_env%dip_mom)) THEN
1340 IF (all(ana_env%dip_mom%last_dip_cl .NE. elem%dipole)) THEN
1341 elem%dipole = pbc(r=elem%dipole(:) - ana_env%dip_mom%last_dip_cl, &
1342 cell=scaled_cell) + ana_env%dip_mom%last_dip_cl
1343 END IF
1344 END IF
1345
1346 ana_env%dip_ana%conf_counter = ana_env%dip_ana%conf_counter + weight_act
1347
1348 ! dipole sqared absolut value summed and weight_acted with volume and conf weight_act
1349 ana_env%dip_ana%mu2_pv_s = ana_env%dip_ana%mu2_pv_s + &
1350 dot_product(elem%dipole(:), elem%dipole(:))/vol*weight_act
1351
1352 tmp_dip(:, :) = 0.0_dp
1353 tmp_dip(:, 1) = elem%dipole(:)
1354
1355 ! dipole sum, weight_acted with volume and conf weight_act
1356 ana_env%dip_ana%mu_pv(:) = ana_env%dip_ana%mu_pv(:) + &
1357 tmp_dip(:, 1)/vol*weight_act
1358
1359 ! dipole sum, weight_acted with square root of volume and conf weight_act
1360 ana_env%dip_ana%mu_psv(:) = ana_env%dip_ana%mu_psv(:) + &
1361 tmp_dip(:, 1)/sqrt(vol)*weight_act
1362
1363 ! dipole squared sum, weight_acted with volume and conf weight_act
1364 ana_env%dip_ana%mu2_pv(:) = ana_env%dip_ana%mu2_pv(:) + &
1365 tmp_dip(:, 1)**2/vol*weight_act
1366
1367 ! calculate the directional average with componentwise correlation per volume
1368 tmp_dip(:, :) = matmul(tmp_dip(:, :), transpose(tmp_dip(:, :)))
1369 ana_env%dip_ana%mu2_pv_mat(:, :) = ana_env%dip_ana%mu2_pv_mat(:, :) + &
1370 tmp_dip(:, :)/vol*weight_act
1371
1372 END SUBROUTINE calc_dipole_analysis
1373
1374! **************************************************************************************************
1375!> \brief prints the actual dipole moment analysis (trajectories)
1376!> \param elem ...
1377!> \param ana_env ...
1378!> \param
1379!> \author Mandes 03.2013
1380! **************************************************************************************************
1381 SUBROUTINE print_act_dipole_analysis(elem, ana_env)
1382 TYPE(tree_type), POINTER :: elem
1383 TYPE(tmc_analysis_env), POINTER :: ana_env
1384
1385 CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp
1386 INTEGER :: counter_tmp, file_ptr
1387 LOGICAL :: flag
1388 REAL(kind=dp) :: diel_const, diel_const_norm, &
1389 diel_const_sym, e0, kb
1390 REAL(kind=dp), DIMENSION(3, 3) :: tmp_dip
1391
1392 kb = boltzmann/joule
1393 counter_tmp = int(ana_env%dip_ana%conf_counter)
1394
1395 ! TODO get correct constant using physcon
1396 e0 = 0.07957747154594767_dp !e^2*a0*me*hbar^-2
1397 diel_const_norm = 1/(3.0_dp*e0*kb*ana_env%temperature)
1398
1399 file_name = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1400 tmc_default_trajectory_file_name, &
1401 ana_env%temperature)
1402 CALL write_dipoles_in_file(file_name=file_name, &
1403 conf_nr=int(ana_env%dip_ana%conf_counter) + 1, dip=elem%dipole, &
1404 file_ext="dip_folded")
1405
1406 ! set output file name
1407 file_name_tmp = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1408 tmc_default_trajectory_file_name, &
1409 ana_env%temperature)
1410
1411 SELECT CASE (ana_env%dip_ana%ana_type)
1412 CASE (ana_type_default)
1413 file_name = trim(expand_file_name_char(file_name_tmp, &
1414 "diel_const"))
1415 file_name_tmp = trim(expand_file_name_char(file_name_tmp, &
1416 "diel_const_tensor"))
1417 CASE (ana_type_sym_xyz)
1418 file_name = trim(expand_file_name_char(file_name_tmp, &
1419 "diel_const_sym"))
1420 file_name_tmp = trim(expand_file_name_char(file_name_tmp, &
1421 "diel_const_tensor_sym"))
1422 CASE DEFAULT
1423 cpwarn('unknown analysis type "'//cp_to_string(ana_env%dip_ana%ana_type)//'" used.')
1424 END SELECT
1425
1426 ! calc the dielectric constant
1427 ! 1+( <M^2> - <M>^2 ) / (3*e_0*V*k*T)
1428 diel_const = 1.0_dp + (ana_env%dip_ana%mu2_pv_s/(ana_env%dip_ana%conf_counter) - &
1429 dot_product(ana_env%dip_ana%mu_psv(:)/(ana_env%dip_ana%conf_counter), &
1430 ana_env%dip_ana%mu_psv(:)/(ana_env%dip_ana%conf_counter)))* &
1431 diel_const_norm
1432 ! symmetrized dielctric constant
1433 ! 1+( <M^2> ) / (3*e_0*V*k*T)
1434 diel_const_sym = 1.0_dp + ana_env%dip_ana%mu2_pv_s/(ana_env%dip_ana%conf_counter)* &
1435 diel_const_norm
1436 ! print dielectric constant trajectory
1437 ! if szmetry used print only every 8th configuration, hence every different (not mirrowed)
1438 INQUIRE (file=file_name, exist=flag)
1439 CALL open_file(file_name=file_name, file_status="UNKNOWN", &
1440 file_action="WRITE", file_position="APPEND", &
1441 unit_number=file_ptr)
1442 IF (.NOT. flag) THEN
1443 WRITE (file_ptr, fmt='(A8,5A20)') "# conf", "diel_const", &
1444 "diel_const_sym", "diel_const_sym_x", &
1445 "diel_const_sym_y", "diel_const_sym_z"
1446 END IF
1447 WRITE (file_ptr, fmt="(I8,10F20.10)") counter_tmp, diel_const, &
1448 diel_const_sym, &
1449 4.0_dp*pi/(kb*ana_env%temperature)* &
1450 ana_env%dip_ana%mu2_pv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1451 CALL close_file(unit_number=file_ptr)
1452
1453 ! print dielectric constant tensor trajectory
1454 INQUIRE (file=file_name_tmp, exist=flag)
1455 CALL open_file(file_name=file_name_tmp, file_status="UNKNOWN", &
1456 file_action="WRITE", file_position="APPEND", &
1457 unit_number=file_ptr)
1458 IF (.NOT. flag) THEN
1459 WRITE (file_ptr, fmt='(A8,9A20)') "# conf", "xx", "xy", "xz", &
1460 "yx", "yy", "yz", &
1461 "zx", "zy", "zz"
1462 END IF
1463 tmp_dip(:, :) = 0.0_dp
1464 tmp_dip(:, 1) = ana_env%dip_ana%mu_psv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1465
1466 WRITE (file_ptr, fmt="(I8,10F20.10)") counter_tmp, &
1467 4.0_dp*pi/(kb*ana_env%temperature)* &
1468 (ana_env%dip_ana%mu2_pv_mat(:, :)/real(ana_env%dip_ana%conf_counter, kind=dp) - &
1469 matmul(tmp_dip(:, :), transpose(tmp_dip(:, :))))
1470 CALL close_file(unit_number=file_ptr)
1471 END SUBROUTINE print_act_dipole_analysis
1472
1473! **************************************************************************************************
1474!> \brief prints the dipole moment analysis
1475!> \param ana_env ...
1476!> \param
1477!> \author Mandes 03.2013
1478! **************************************************************************************************
1479 SUBROUTINE print_dipole_analysis(ana_env)
1480 TYPE(tmc_analysis_env), POINTER :: ana_env
1481
1482 CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA"
1483
1484 INTEGER :: i
1485 REAL(kind=dp) :: diel_const_scalar, kb
1486 REAL(kind=dp), DIMENSION(3) :: diel_const_sym, dielec_ev
1487 REAL(kind=dp), DIMENSION(3, 3) :: diel_const, tmp_dip, tmp_ev
1488
1489 kb = boltzmann/joule
1490
1491 cpassert(ASSOCIATED(ana_env))
1492 cpassert(ASSOCIATED(ana_env%dip_ana))
1493
1494 tmp_dip(:, :) = 0.0_dp
1495 diel_const(:, :) = 0.0_dp
1496 diel_const_scalar = 0.0_dp
1497 diel_const_sym = 0.0_dp
1498
1499 !dielectric constant
1500 tmp_dip(:, 1) = ana_env%dip_ana%mu_psv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1501 diel_const(:, :) = 4.0_dp*pi/(kb*ana_env%temperature)* &
1502 (ana_env%dip_ana%mu2_pv_mat(:, :)/real(ana_env%dip_ana%conf_counter, kind=dp) - &
1503 matmul(tmp_dip(:, :), transpose(tmp_dip(:, :))))
1504
1505 !dielectric constant for symmetric case
1506 diel_const_sym(:) = 4.0_dp*pi/(kb*ana_env%temperature)* &
1507 ana_env%dip_ana%mu2_pv(:)/real(ana_env%dip_ana%conf_counter, kind=dp)
1508
1509 DO i = 1, 3
1510 diel_const(i, i) = diel_const(i, i) + 1.0_dp ! +1 for unpolarizable models, 1.592 for polarizable
1511 diel_const_scalar = diel_const_scalar + diel_const(i, i)
1512 END DO
1513 diel_const_scalar = diel_const_scalar/real(3, kind=dp)
1514
1515 tmp_dip(:, :) = diel_const
1516 CALL diag(3, tmp_dip, dielec_ev, tmp_ev)
1517
1518 ! print out results
1519 WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
1520 WRITE (ana_env%io_unit, fmt="(T2,A,T35,A,T80,A)") "-", "average dipoles", "-"
1521 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
1522 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "used configurations ", &
1523 cp_to_string(real(ana_env%dip_ana%conf_counter, kind=dp))
1524 IF (ana_env%dip_ana%ana_type .EQ. ana_type_ice) &
1525 WRITE (ana_env%io_unit, fmt='(T2,A,"| ",A)') plabel, &
1526 "ice analysis with directions of hexagonal structure"
1527 IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
1528 WRITE (ana_env%io_unit, fmt='(T2,A,"| ",A)') plabel, &
1529 "ice analysis with symmetrized dipoles in each direction."
1530
1531 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "for product of 2 directions(per vol):"
1532 DO i = 1, 3
1533 WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", ana_env%dip_ana%mu2_pv_mat(i, :)/ &
1534 REAL(ana_env%dip_ana%conf_counter, kind=dp), " |"
1535 END DO
1536
1537 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric constant tensor:"
1538 DO i = 1, 3
1539 WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", diel_const(i, :), " |"
1540 END DO
1541
1542 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric tensor eigenvalues", &
1543 cp_to_string(dielec_ev(1))//" "// &
1544 cp_to_string(dielec_ev(2))//" "// &
1545 cp_to_string(dielec_ev(3))
1546 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric constant symm ", &
1547 cp_to_string(diel_const_sym(1))//" | "// &
1548 cp_to_string(diel_const_sym(2))//" | "// &
1549 cp_to_string(diel_const_sym(3))
1550 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "dielectric constant ", &
1551 cp_to_string(diel_const_scalar)
1552 WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
1553
1554 END SUBROUTINE print_dipole_analysis
1555
1556 !============================================================================
1557 ! particle displacement in cell (from one configuration to the next)
1558 !============================================================================
1559
1560! **************************************************************************************************
1561!> \brief calculates the mean square displacement
1562!> \param elem ...
1563!> \param ana_env ...
1564!> \param
1565!> \author Mandes 02.2013
1566! **************************************************************************************************
1567 SUBROUTINE calc_displacement(elem, ana_env)
1568 TYPE(tree_type), POINTER :: elem
1569 TYPE(tmc_analysis_env), POINTER :: ana_env
1570
1571 CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_displacement'
1572
1573 CHARACTER(LEN=default_path_length) :: file_name, file_name_tmp
1574 INTEGER :: file_ptr, handle, ind
1575 LOGICAL :: flag
1576 REAL(kind=dp) :: disp
1577 REAL(kind=dp), DIMENSION(3) :: atom_disp
1578
1579 disp = 0.0_dp
1580
1581 cpassert(ASSOCIATED(elem))
1582 cpassert(ASSOCIATED(elem%pos))
1583 cpassert(ASSOCIATED(ana_env))
1584 cpassert(ASSOCIATED(ana_env%displace))
1585 cpassert(ASSOCIATED(ana_env%last_elem))
1586
1587 ! start the timing
1588 CALL timeset(routinen, handle)
1589
1590 DO ind = 1, SIZE(elem%pos), ana_env%dim_per_elem
1591 ! fold into box
1592 atom_disp(:) = elem%pos(ind:ind + 2) - ana_env%last_elem%pos(ind:ind + 2)
1593 CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
1594 vec=atom_disp)
1595 disp = disp + sum((atom_disp(:)*au2a)**2)
1596 END DO
1597 ana_env%displace%disp = ana_env%displace%disp + disp
1598 ana_env%displace%conf_counter = ana_env%displace%conf_counter + 1
1599
1600 IF (ana_env%displace%print_disp) THEN
1601 file_name_tmp = expand_file_name_temp(trim(ana_env%out_file_prefix)// &
1602 tmc_default_trajectory_file_name, &
1603 ana_env%temperature)
1604 file_name = trim(expand_file_name_char(file_name_tmp, &
1605 "devi"))
1606 INQUIRE (file=file_name, exist=flag)
1607 CALL open_file(file_name=file_name, file_status="UNKNOWN", &
1608 file_action="WRITE", file_position="APPEND", &
1609 unit_number=file_ptr)
1610 IF (.NOT. flag) &
1611 WRITE (file_ptr, *) "# conf squared deviation of the cell"
1612 WRITE (file_ptr, *) elem%nr, disp
1613 CALL close_file(unit_number=file_ptr)
1614 END IF
1615
1616 ! end the timing
1617 CALL timestop(handle)
1618
1619 END SUBROUTINE calc_displacement
1620
1621! **************************************************************************************************
1622!> \brief prints final values for the displacement calculations
1623!> \param ana_env ...
1624!> \param
1625!> \author Mandes 02.2013
1626! **************************************************************************************************
1627 SUBROUTINE print_average_displacement(ana_env)
1628 TYPE(tmc_analysis_env), POINTER :: ana_env
1629
1630 CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA"
1631
1632 WRITE (ana_env%io_unit, fmt="(/,T2,A)") repeat("-", 79)
1633 WRITE (ana_env%io_unit, fmt="(T2,A,T35,A,T80,A)") "-", "average displacement", "-"
1634 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "temperature ", &
1635 cp_to_string(ana_env%temperature)
1636 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "used configurations ", &
1637 cp_to_string(real(ana_env%displace%conf_counter, kind=dp))
1638 WRITE (ana_env%io_unit, fmt=fmt_my) plabel, "cell root mean square deviation: ", &
1639 cp_to_string(sqrt(ana_env%displace%disp/ &
1640 REAL(ana_env%displace%conf_counter, kind=dp)))
1641 IF (ana_env%print_test_output) &
1642 WRITE (*, *) "TMC|ANALYSIS_AVERAGE_CELL_DISPLACEMENT_X= ", &
1643 sqrt(ana_env%displace%disp/ &
1644 REAL(ana_env%displace%conf_counter, kind=dp))
1645 END SUBROUTINE print_average_displacement
1646END MODULE tmc_analysis
cell_types::pbc
Definition cell_types.F:92
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: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:1934
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:731
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:704
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:1047
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:1496
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:592
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:638
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:1045
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:574
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:620
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:945
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:55
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