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