Line data Source code
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 :
16 : MODULE tmc_analysis
17 : USE cell_types, ONLY: cell_type,&
18 : get_cell,&
19 : pbc
20 : USE cp_files, ONLY: close_file,&
21 : open_file
22 : USE cp_log_handling, ONLY: cp_to_string
23 : USE force_fields_input, ONLY: read_chrg_section
24 : USE input_section_types, ONLY: section_vals_get,&
25 : section_vals_get_subs_vals,&
26 : section_vals_type,&
27 : section_vals_val_get
28 : USE kinds, ONLY: default_path_length,&
29 : default_string_length,&
30 : dp
31 : USE mathconstants, ONLY: pi
32 : USE mathlib, ONLY: diag
33 : USE physcon, ONLY: a_mass,&
34 : au2a => angstrom,&
35 : boltzmann,&
36 : joule,&
37 : massunit
38 : USE tmc_analysis_types, ONLY: &
39 : ana_type_default, ana_type_ice, ana_type_sym_xyz, atom_pairs_type, dipole_moment_type, &
40 : pair_correl_type, search_pair_in_list, tmc_ana_density_create, tmc_ana_density_file_name, &
41 : tmc_ana_dipole_analysis_create, tmc_ana_dipole_moment_create, tmc_ana_displacement_create, &
42 : tmc_ana_env_create, tmc_ana_pair_correl_create, tmc_ana_pair_correl_file_name, &
43 : tmc_analysis_env
44 : USE tmc_calculations, ONLY: get_scaled_cell,&
45 : nearest_distance
46 : USE tmc_file_io, ONLY: analyse_files_close,&
47 : analyse_files_open,&
48 : expand_file_name_char,&
49 : expand_file_name_temp,&
50 : read_element_from_file,&
51 : write_dipoles_in_file
52 : USE tmc_stati, ONLY: TMC_STATUS_OK,&
53 : TMC_STATUS_WAIT_FOR_NEW_TASK,&
54 : tmc_default_restart_in_file_name,&
55 : tmc_default_restart_out_file_name,&
56 : tmc_default_trajectory_file_name,&
57 : tmc_default_unspecified_name
58 : USE tmc_tree_build, ONLY: allocate_new_sub_tree_node,&
59 : deallocate_sub_tree_node
60 : USE tmc_tree_types, ONLY: read_subtree_elem_unformated,&
61 : tree_type,&
62 : write_subtree_elem_unformated
63 : USE tmc_types, ONLY: tmc_atom_type,&
64 : tmc_param_type
65 : #include "../base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 :
69 : PRIVATE
70 :
71 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_analysis'
72 :
73 : PUBLIC :: tmc_read_ana_input
74 : PUBLIC :: analysis_init, do_tmc_analysis, analyze_file_configurations, finalize_tmc_analysis
75 : PUBLIC :: analysis_restart_print, analysis_restart_read
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief creates a new para environment for tmc analysis
81 : !> \param tmc_ana_section ...
82 : !> \param tmc_ana TMC analysis environment
83 : !> \author Mandes 02.2013
84 : ! **************************************************************************************************
85 54 : 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 18 : CHARACTER(LEN=default_string_length), POINTER :: charge_atm(:)
91 : INTEGER :: i_tmp, ntot
92 : INTEGER, DIMENSION(3) :: nr_bins
93 18 : INTEGER, DIMENSION(:), POINTER :: i_arr_tmp
94 : LOGICAL :: explicit, explicit_key, flag
95 18 : REAL(KIND=dp), POINTER :: charge(:)
96 : TYPE(section_vals_type), POINTER :: tmp_section
97 :
98 18 : NULLIFY (tmp_section, charge_atm, i_arr_tmp, charge)
99 :
100 0 : CPASSERT(ASSOCIATED(tmc_ana_section))
101 18 : CPASSERT(.NOT. ASSOCIATED(tmc_ana))
102 :
103 18 : CALL section_vals_get(tmc_ana_section, explicit=explicit)
104 18 : IF (explicit) THEN
105 18 : CALL tmc_ana_env_create(tmc_ana=tmc_ana)
106 : ! restarting
107 18 : 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 18 : c_val=tmc_ana%out_file_prefix)
111 18 : IF (tmc_ana%out_file_prefix .NE. "") THEN
112 0 : tmc_ana%out_file_prefix = TRIM(tmc_ana%out_file_prefix)//"_"
113 : END IF
114 :
115 : ! density calculation
116 18 : CALL section_vals_val_get(tmc_ana_section, "DENSITY", explicit=explicit_key)
117 18 : IF (explicit_key) THEN
118 9 : CALL section_vals_val_get(tmc_ana_section, "DENSITY", i_vals=i_arr_tmp)
119 :
120 9 : IF (SIZE(i_arr_tmp(:)) .EQ. 3) THEN
121 36 : IF (ANY(i_arr_tmp(:) .LE. 0)) &
122 : CALL cp_abort(__LOCATION__, "The amount of intervals in each "// &
123 0 : "direction has to be greater than 0.")
124 36 : nr_bins(:) = i_arr_tmp(:)
125 0 : ELSE IF (SIZE(i_arr_tmp(:)) .EQ. 1) THEN
126 0 : IF (ANY(i_arr_tmp(:) .LE. 0)) &
127 0 : CPABORT("The amount of intervals has to be greater than 0.")
128 0 : nr_bins(:) = i_arr_tmp(1)
129 0 : ELSE IF (SIZE(i_arr_tmp(:)) .EQ. 0) THEN
130 0 : nr_bins(:) = 1
131 : ELSE
132 0 : CPABORT("unknown amount of dimensions for the binning.")
133 : END IF
134 9 : CALL tmc_ana_density_create(tmc_ana%density_3d, nr_bins)
135 : END IF
136 :
137 : ! radial distribution function calculation
138 18 : CALL section_vals_val_get(tmc_ana_section, "G_R", explicit=explicit_key)
139 18 : IF (explicit_key) THEN
140 9 : 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 9 : nr_bins=i_tmp)
143 : END IF
144 :
145 : ! radial distribution function calculation
146 18 : CALL section_vals_val_get(tmc_ana_section, "CLASSICAL_DIPOLE_MOMENTS", explicit=explicit_key)
147 18 : IF (explicit_key) THEN
148 : ! charges for dipoles needed
149 9 : tmp_section => section_vals_get_subs_vals(tmc_ana_section, "CHARGE")
150 9 : CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=i_tmp)
151 9 : IF (explicit) THEN
152 9 : ntot = 0
153 27 : ALLOCATE (charge_atm(i_tmp))
154 27 : ALLOCATE (charge(i_tmp))
155 9 : 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 0 : "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 9 : tmc_ana%dim_per_elem)
164 :
165 9 : IF (ASSOCIATED(charge_atm)) DEALLOCATE (charge_atm)
166 9 : IF (ASSOCIATED(charge)) DEALLOCATE (charge)
167 : END IF
168 :
169 : ! dipole moment analysis
170 18 : CALL section_vals_val_get(tmc_ana_section, "DIPOLE_ANALYSIS", explicit=explicit_key)
171 18 : IF (explicit_key) THEN
172 0 : CALL tmc_ana_dipole_analysis_create(tmc_ana%dip_ana)
173 0 : CALL section_vals_val_get(tmc_ana_section, "DIPOLE_ANALYSIS", c_val=c_tmp)
174 0 : SELECT CASE (TRIM(c_tmp))
175 : CASE (TRIM(tmc_default_unspecified_name))
176 0 : tmc_ana%dip_ana%ana_type = ana_type_default
177 : CASE ("ICE")
178 0 : tmc_ana%dip_ana%ana_type = ana_type_ice
179 : CASE ("SYM_XYZ")
180 0 : tmc_ana%dip_ana%ana_type = ana_type_sym_xyz
181 : CASE DEFAULT
182 0 : CPWARN('unknown analysis type "'//TRIM(c_tmp)//'" specified. Set to default.')
183 0 : 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 18 : CALL section_vals_val_get(tmc_ana_section, "DEVIATION", l_val=flag)
191 18 : IF (flag) THEN
192 : CALL tmc_ana_displacement_create(ana_disp=tmc_ana%displace, &
193 9 : dim_per_elem=tmc_ana%dim_per_elem)
194 : END IF
195 18 : 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 18 : 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 18 : CPASSERT(ASSOCIATED(ana_env))
210 18 : CPASSERT(nr_dim > 0)
211 :
212 18 : ana_env%nr_dim = nr_dim
213 :
214 : ! save file names
215 18 : tmp_pos_file = ana_env%costum_pos_file_name
216 18 : tmp_cell_file = ana_env%costum_cell_file_name
217 18 : tmp_dip_file = ana_env%costum_dip_file_name
218 :
219 : ! unset all filenames
220 18 : ana_env%costum_pos_file_name = tmc_default_unspecified_name
221 18 : ana_env%costum_cell_file_name = tmc_default_unspecified_name
222 18 : ana_env%costum_dip_file_name = tmc_default_unspecified_name
223 :
224 : ! set the necessary files for ...
225 : ! density
226 18 : IF (ASSOCIATED(ana_env%density_3d)) THEN
227 9 : ana_env%costum_pos_file_name = tmp_pos_file
228 9 : ana_env%costum_cell_file_name = tmp_cell_file
229 : END IF
230 : ! pair correlation
231 18 : IF (ASSOCIATED(ana_env%pair_correl)) THEN
232 9 : ana_env%costum_pos_file_name = tmp_pos_file
233 9 : ana_env%costum_cell_file_name = tmp_cell_file
234 : END IF
235 : ! dipole moment
236 18 : IF (ASSOCIATED(ana_env%dip_mom)) THEN
237 9 : ana_env%costum_pos_file_name = tmp_pos_file
238 9 : ana_env%costum_cell_file_name = tmp_cell_file
239 : END IF
240 : ! dipole analysis
241 18 : IF (ASSOCIATED(ana_env%dip_ana)) THEN
242 0 : ana_env%costum_pos_file_name = tmp_pos_file
243 0 : ana_env%costum_cell_file_name = tmp_cell_file
244 0 : ana_env%costum_dip_file_name = tmp_dip_file
245 : END IF
246 : ! deviation / displacement
247 18 : IF (ASSOCIATED(ana_env%displace)) THEN
248 9 : ana_env%costum_pos_file_name = tmp_pos_file
249 9 : ana_env%costum_cell_file_name = tmp_cell_file
250 : END IF
251 :
252 : ! init radial distribution function
253 18 : IF (ASSOCIATED(ana_env%pair_correl)) &
254 : CALL ana_pair_correl_init(ana_pair_correl=ana_env%pair_correl, &
255 9 : atoms=ana_env%atoms, cell=ana_env%cell)
256 : ! init classical dipole moment calculations
257 18 : IF (ASSOCIATED(ana_env%dip_mom)) &
258 : CALL ana_dipole_moment_init(ana_dip_mom=ana_env%dip_mom, &
259 9 : atoms=ana_env%atoms)
260 18 : 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 18 : 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 18 : CPASSERT(ASSOCIATED(ana_env))
277 18 : CPASSERT(ASSOCIATED(ana_env%last_elem))
278 18 : IF (.NOT. ana_env%restart) RETURN
279 :
280 6 : 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 6 : "ana"), ana_env%temperature))
285 : restart_file_name = expand_file_name_char(file_name_tmp, &
286 6 : file_name)
287 : CALL open_file(file_name=restart_file_name, file_status="REPLACE", &
288 : file_action="WRITE", file_form="UNFORMATTED", &
289 6 : unit_number=file_ptr)
290 6 : WRITE (file_ptr) ana_env%temperature
291 6 : 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 6 : l_tmp = ASSOCIATED(ana_env%density_3d)
296 6 : WRITE (file_ptr) l_tmp
297 6 : IF (l_tmp) THEN
298 6 : WRITE (file_ptr) ana_env%density_3d%conf_counter, &
299 6 : ana_env%density_3d%nr_bins, &
300 6 : ana_env%density_3d%sum_vol, &
301 6 : ana_env%density_3d%sum_vol2, &
302 6 : ana_env%density_3d%sum_box_length, &
303 6 : ana_env%density_3d%sum_box_length2, &
304 30 : ana_env%density_3d%sum_density, &
305 36 : ana_env%density_3d%sum_dens2
306 : END IF
307 :
308 6 : l_tmp = ASSOCIATED(ana_env%pair_correl)
309 6 : WRITE (file_ptr) l_tmp
310 6 : IF (l_tmp) THEN
311 6 : WRITE (file_ptr) ana_env%pair_correl%conf_counter, &
312 6 : ana_env%pair_correl%nr_bins, &
313 6 : ana_env%pair_correl%step_length, &
314 24 : ana_env%pair_correl%pairs, &
315 5436 : ana_env%pair_correl%g_r
316 : END IF
317 :
318 6 : l_tmp = ASSOCIATED(ana_env%dip_mom)
319 6 : WRITE (file_ptr) l_tmp
320 6 : IF (l_tmp) THEN
321 6 : WRITE (file_ptr) ana_env%dip_mom%conf_counter, &
322 132 : ana_env%dip_mom%charges, &
323 30 : ana_env%dip_mom%last_dip_cl
324 : END IF
325 :
326 6 : l_tmp = ASSOCIATED(ana_env%dip_ana)
327 6 : WRITE (file_ptr) l_tmp
328 6 : IF (l_tmp) THEN
329 0 : WRITE (file_ptr) ana_env%dip_ana%conf_counter, &
330 0 : ana_env%dip_ana%ana_type, &
331 0 : ana_env%dip_ana%mu2_pv_s, &
332 0 : ana_env%dip_ana%mu_psv, &
333 0 : ana_env%dip_ana%mu_pv, &
334 0 : ana_env%dip_ana%mu2_pv_mat, &
335 0 : ana_env%dip_ana%mu2_pv_mat
336 : END IF
337 :
338 6 : l_tmp = ASSOCIATED(ana_env%displace)
339 6 : WRITE (file_ptr) l_tmp
340 6 : IF (l_tmp) THEN
341 6 : WRITE (file_ptr) ana_env%displace%conf_counter, &
342 12 : ana_env%displace%disp
343 : END IF
344 :
345 6 : CALL close_file(unit_number=file_ptr)
346 :
347 : file_name_tmp = expand_file_name_char(TRIM(ana_env%out_file_prefix)// &
348 6 : tmc_default_restart_in_file_name, "ana")
349 : file_name = expand_file_name_temp(file_name_tmp, &
350 6 : ana_env%temperature)
351 : CALL open_file(file_name=file_name, &
352 : file_action="WRITE", file_status="REPLACE", &
353 6 : unit_number=file_ptr)
354 6 : WRITE (file_ptr, *) TRIM(restart_file_name)
355 6 : 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 18 : 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 18 : CPASSERT(ASSOCIATED(ana_env))
375 18 : CPASSERT(ASSOCIATED(elem))
376 18 : IF (.NOT. ana_env%restart) RETURN
377 :
378 : file_name_tmp = expand_file_name_char(TRIM(ana_env%out_file_prefix)// &
379 6 : tmc_default_restart_in_file_name, "ana")
380 : file_name = expand_file_name_temp(file_name_tmp, &
381 6 : ana_env%temperature)
382 6 : INQUIRE (FILE=file_name, EXIST=l_tmp)
383 6 : IF (l_tmp) THEN
384 : CALL open_file(file_name=file_name, file_status="OLD", &
385 3 : file_action="READ", unit_number=file_ptr)
386 3 : READ (file_ptr, *) file_name_tmp
387 3 : CALL close_file(unit_number=file_ptr)
388 :
389 : CALL open_file(file_name=file_name_tmp, file_status="OLD", file_form="UNFORMATTED", &
390 3 : file_action="READ", unit_number=file_ptr)
391 3 : READ (file_ptr) temp
392 3 : CPASSERT(ana_env%temperature .EQ. temp)
393 3 : ana_env%last_elem => elem
394 3 : 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 3 : READ (file_ptr) l_tmp
399 3 : CPASSERT(ASSOCIATED(ana_env%density_3d) .EQV. l_tmp)
400 3 : IF (l_tmp) THEN
401 3 : READ (file_ptr) ana_env%density_3d%conf_counter, &
402 3 : ana_env%density_3d%nr_bins, &
403 3 : ana_env%density_3d%sum_vol, &
404 3 : ana_env%density_3d%sum_vol2, &
405 3 : ana_env%density_3d%sum_box_length, &
406 3 : ana_env%density_3d%sum_box_length2, &
407 15 : ana_env%density_3d%sum_density, &
408 18 : ana_env%density_3d%sum_dens2
409 : END IF
410 :
411 3 : READ (file_ptr) l_tmp
412 3 : CPASSERT(ASSOCIATED(ana_env%pair_correl) .EQV. l_tmp)
413 3 : IF (l_tmp) THEN
414 3 : READ (file_ptr) ana_env%pair_correl%conf_counter, &
415 3 : ana_env%pair_correl%nr_bins, &
416 3 : ana_env%pair_correl%step_length, &
417 12 : ana_env%pair_correl%pairs, &
418 2718 : ana_env%pair_correl%g_r
419 : END IF
420 :
421 3 : READ (file_ptr) l_tmp
422 3 : CPASSERT(ASSOCIATED(ana_env%dip_mom) .EQV. l_tmp)
423 3 : IF (l_tmp) THEN
424 3 : READ (file_ptr) ana_env%dip_mom%conf_counter, &
425 66 : ana_env%dip_mom%charges, &
426 15 : ana_env%dip_mom%last_dip_cl
427 : END IF
428 :
429 3 : READ (file_ptr) l_tmp
430 3 : CPASSERT(ASSOCIATED(ana_env%dip_ana) .EQV. l_tmp)
431 3 : IF (l_tmp) THEN
432 0 : READ (file_ptr) ana_env%dip_ana%conf_counter, &
433 0 : ana_env%dip_ana%ana_type, &
434 0 : ana_env%dip_ana%mu2_pv_s, &
435 0 : ana_env%dip_ana%mu_psv, &
436 0 : ana_env%dip_ana%mu_pv, &
437 0 : ana_env%dip_ana%mu2_pv_mat, &
438 0 : ana_env%dip_ana%mu2_pv_mat
439 : END IF
440 :
441 3 : READ (file_ptr) l_tmp
442 3 : CPASSERT(ASSOCIATED(ana_env%displace) .EQV. l_tmp)
443 3 : IF (l_tmp) THEN
444 3 : READ (file_ptr) ana_env%displace%conf_counter, &
445 6 : ana_env%displace%disp
446 : END IF
447 :
448 3 : CALL close_file(unit_number=file_ptr)
449 3 : 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 2062 : 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 1031 : CPASSERT(ASSOCIATED(elem))
475 1031 : CPASSERT(ASSOCIATED(ana_env))
476 :
477 : ! start the timing
478 1031 : CALL timeset(routineN, handle)
479 :
480 1031 : weight_act = 0
481 1031 : IF (ASSOCIATED(ana_env%last_elem)) THEN
482 1016 : weight_act = elem%nr - ana_env%last_elem%nr
483 : END IF
484 :
485 1031 : IF (weight_act .GT. 0) THEN
486 : ! calculates the 3 dimensional distributed density
487 1016 : 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 500 : ana_env=ana_env)
491 : ! calculated the radial distribution function for each atom type
492 1016 : IF (ASSOCIATED(ana_env%pair_correl)) &
493 : CALL calc_paircorrelation(elem=ana_env%last_elem, weight=weight_act, &
494 500 : atoms=ana_env%atoms, ana_env=ana_env)
495 : ! calculates the classical dipole moments
496 1016 : IF (ASSOCIATED(ana_env%dip_mom)) &
497 : CALL calc_dipole_moment(elem=ana_env%last_elem, weight=weight_act, &
498 500 : ana_env=ana_env)
499 : ! calculates the dipole moments analysis and dielectric constant
500 1016 : 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 0 : IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) THEN
504 : ! (-x,y,z)
505 0 : ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
506 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
507 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
508 0 : 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 0 : ana_env=ana_env)
511 : ! (-x,-y,z)
512 0 : ana_env%last_elem%dipole(:) = dip_tmp(:)
513 0 : ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
514 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
515 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
516 0 : 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 0 : ana_env=ana_env)
519 : ! (-x,-y,-z)
520 0 : ana_env%last_elem%dipole(:) = dip_tmp(:)
521 0 : ana_env%last_elem%dipole(3) = -ana_env%last_elem%dipole(3)
522 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
523 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
524 0 : 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 0 : ana_env=ana_env)
527 : ! (x,-y,-z)
528 0 : ana_env%last_elem%dipole(:) = dip_tmp(:)
529 0 : ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
530 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
531 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
532 0 : 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 0 : ana_env=ana_env)
535 : ! (x,y,-z)
536 0 : ana_env%last_elem%dipole(:) = dip_tmp(:)
537 0 : ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
538 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
539 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
540 0 : 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 0 : ana_env=ana_env)
543 : ! (-x,y,-z)
544 0 : ana_env%last_elem%dipole(:) = dip_tmp(:)
545 0 : ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
546 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
547 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
548 0 : 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 0 : ana_env=ana_env)
551 : ! (x,-y,z)
552 0 : ana_env%last_elem%dipole(:) = dip_tmp(:)
553 0 : ana_env%last_elem%dipole(:) = -ana_env%last_elem%dipole(:)
554 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
555 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
556 0 : 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 0 : ana_env=ana_env)
559 : ! back to (x,y,z)
560 0 : ana_env%last_elem%dipole(:) = dip_tmp(:)
561 0 : ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
562 0 : dip_tmp(:) = ana_env%last_elem%dipole(:)
563 0 : IF (ASSOCIATED(ana_env%dip_mom)) &
564 0 : 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 0 : ana_env=ana_env)
568 : CALL print_act_dipole_analysis(elem=ana_env%last_elem, &
569 0 : ana_env=ana_env)
570 : END IF
571 :
572 : ! calculates the cell displacement from last cell
573 1016 : IF (ASSOCIATED(ana_env%displace)) THEN
574 500 : 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 1031 : elem_tmp => ana_env%last_elem
579 1031 : ana_env%last_elem => elem
580 1031 : elem => elem_tmp
581 : ! end the timing
582 1031 : CALL timestop(handle)
583 1031 : 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 36 : 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 18 : CPASSERT(ASSOCIATED(ana_env))
599 :
600 : ! start the timing
601 18 : CALL timeset(routineN, handle)
602 18 : IF (ASSOCIATED(ana_env%density_3d)) THEN
603 9 : IF (ana_env%density_3d%conf_counter .GT. 0) &
604 9 : CALL print_density_3d(ana_env=ana_env)
605 : END IF
606 18 : IF (ASSOCIATED(ana_env%pair_correl)) THEN
607 9 : IF (ana_env%pair_correl%conf_counter .GT. 0) &
608 9 : CALL print_paircorrelation(ana_env=ana_env)
609 : END IF
610 18 : IF (ASSOCIATED(ana_env%dip_mom)) THEN
611 9 : IF (ana_env%dip_mom%conf_counter .GT. 0) &
612 9 : CALL print_dipole_moment(ana_env)
613 : END IF
614 18 : IF (ASSOCIATED(ana_env%dip_ana)) THEN
615 0 : IF (ana_env%dip_ana%conf_counter .GT. 0) &
616 0 : CALL print_dipole_analysis(ana_env)
617 : END IF
618 18 : IF (ASSOCIATED(ana_env%displace)) THEN
619 9 : IF (ana_env%displace%conf_counter .GT. 0) &
620 9 : CALL print_average_displacement(ana_env)
621 : END IF
622 :
623 : ! end the timing
624 18 : CALL timestop(handle)
625 18 : 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 36 : 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 18 : NULLIFY (elem)
649 18 : conf_nr = -1
650 18 : stat = TMC_STATUS_WAIT_FOR_NEW_TASK
651 18 : CPASSERT(ASSOCIATED(ana_env))
652 18 : CPASSERT(ASSOCIATED(tmc_params))
653 :
654 : ! start the timing
655 18 : CALL timeset(routineN, handle)
656 :
657 : ! open the files
658 18 : CALL analyse_files_open(tmc_ana=ana_env, stat=stat, dir_ind=dir_ind)
659 : ! set the existence of exact dipoles (from file)
660 18 : IF (ana_env%id_dip .GT. 0) THEN
661 0 : tmc_params%print_dipole = .TRUE.
662 : ELSE
663 18 : 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 18 : nr_dim=ana_env%nr_dim)
669 :
670 18 : IF (ASSOCIATED(ana_env%last_elem)) conf_nr = ana_env%last_elem%nr
671 18 : nr_dim = SIZE(elem%pos)
672 :
673 18 : 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 1049 : stat=stat)
677 1049 : IF (stat .EQ. TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
678 18 : 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 1031 : IF (start_id .LT. 0 .OR. conf_nr .GE. start_id) THEN
683 1031 : IF (end_id .LT. 0 .OR. conf_nr .LE. end_id) THEN
684 : ! do the analysis calculations
685 1031 : CALL do_tmc_analysis(elem=elem, ana_env=ana_env)
686 : END IF
687 : END IF
688 :
689 : ! clean temporary element (already analyzed)
690 1031 : IF (ASSOCIATED(elem)) THEN
691 1016 : 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 1031 : IF (.NOT. ASSOCIATED(elem)) &
695 : CALL allocate_new_sub_tree_node(tmc_params=tmc_params, next_el=elem, &
696 1031 : nr_dim=nr_dim)
697 : END DO conf_loop
698 : END IF
699 : ! close the files
700 18 : CALL analyse_files_close(tmc_ana=ana_env)
701 :
702 18 : IF (ASSOCIATED(elem)) &
703 0 : CALL deallocate_sub_tree_node(tree_elem=elem)
704 :
705 : ! end the timing
706 18 : CALL timestop(handle)
707 18 : 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 500 : 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 500 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mass_bin
738 :
739 500 : NULLIFY (mass_bin)
740 :
741 0 : CPASSERT(ASSOCIATED(elem))
742 500 : CPASSERT(ASSOCIATED(elem%pos))
743 500 : CPASSERT(weight .GT. 0)
744 500 : CPASSERT(ASSOCIATED(atoms))
745 500 : CPASSERT(ASSOCIATED(ana_env))
746 500 : CPASSERT(ASSOCIATED(ana_env%cell))
747 500 : CPASSERT(ASSOCIATED(ana_env%density_3d))
748 500 : CPASSERT(ASSOCIATED(ana_env%density_3d%sum_density))
749 500 : CPASSERT(ASSOCIATED(ana_env%density_3d%sum_dens2))
750 :
751 : ! start the timing
752 500 : CALL timeset(routineN, handle)
753 :
754 500 : atom_pos(:) = 0.0_dp
755 500 : cell_size(:) = 0.0_dp
756 500 : interval_size(:) = 0.0_dp
757 500 : mass_total = 0.0_dp
758 :
759 500 : bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
760 500 : bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
761 500 : bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
762 2500 : ALLOCATE (mass_bin(bin_x, bin_y, bin_z))
763 2500 : 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 500 : 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 500 : vol_cell*(au2a)**3*weight
773 : ana_env%density_3d%sum_vol2 = ana_env%density_3d%sum_vol2 + &
774 500 : (vol_cell*(au2a)**3)**2*weight
775 :
776 : ana_env%density_3d%sum_box_length(:) = ana_env%density_3d%sum_box_length(:) &
777 2000 : + cell_size(:)*(au2a)*weight
778 : ana_env%density_3d%sum_box_length2(:) = ana_env%density_3d%sum_box_length2(:) &
779 2000 : + (cell_size(:)*(au2a))**2*weight
780 :
781 : ! sub interval length
782 500 : interval_size(1) = cell_size(1)/REAL(bin_x, dp)
783 500 : interval_size(2) = cell_size(2)/REAL(bin_y, dp)
784 500 : interval_size(3) = cell_size(3)/REAL(bin_z, dp)
785 :
786 : ! volume in [cm^3]
787 500 : vol_cell = vol_cell*(au2a*1E-8)**3
788 : vol_sub_box = interval_size(1)*interval_size(2)*interval_size(3)* &
789 500 : (au2a*1E-8)**3
790 :
791 : ! count every atom
792 500 : DO atom = 1, SIZE(elem%pos), ana_env%dim_per_elem
793 :
794 42000 : 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 10500 : vec=atom_pos)
798 : ! shifts the box to positive values (before 0,0,0 is the center)
799 42000 : atom_pos(:) = atom_pos(:) + 0.5_dp*cell_size(:)
800 : ! calculate the index of the sub box
801 10500 : bin_x = INT(atom_pos(1)/interval_size(1)) + 1
802 10500 : bin_y = INT(atom_pos(2)/interval_size(2)) + 1
803 10500 : bin_z = INT(atom_pos(3)/interval_size(3)) + 1
804 10500 : CPASSERT(bin_x .GT. 0 .AND. bin_y .GT. 0 .AND. bin_z .GT. 0)
805 10500 : CPASSERT(bin_x .LE. SIZE(ana_env%density_3d%sum_density(:, 1, 1)))
806 10500 : CPASSERT(bin_y .LE. SIZE(ana_env%density_3d%sum_density(1, :, 1)))
807 10500 : 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 10500 : atoms(INT(atom/REAL(ana_env%dim_per_elem, KIND=dp)) + 1)%mass/massunit*1000*a_mass
812 : mass_total = mass_total + &
813 10500 : 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 4000 : r_tmp = mass_total/vol_cell - SUM(mass_bin(:, :, :))/vol_sub_box/SIZE(mass_bin(:, :, :))
823 500 : 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 4500 : 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 4500 : weight*(mass_bin(:, :, :)/vol_sub_box)**2
832 :
833 500 : 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 500 : 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 500 : ana_env%temperature)
840 : file_name = TRIM(expand_file_name_char(file_name_tmp, &
841 500 : "dens"))
842 500 : 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 500 : unit_number=file_ptr)
846 500 : IF (.NOT. flag) &
847 3 : WRITE (file_ptr, FMT='(A8,11A20)') "# conf_nr", "dens_act[g/cm^3]", &
848 3 : "dens_average[g/cm^3]", "density_variance", &
849 3 : "averages:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z", &
850 6 : "variances:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z"
851 500 : WRITE (file_ptr, FMT="(I8,11F20.10)") ana_env%density_3d%conf_counter + 1 - weight, &
852 4000 : 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 4000 : 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 7500 : REAL(ana_env%density_3d%conf_counter, KIND=dp))**2, &
862 : ana_env%density_3d%sum_vol/ &
863 500 : REAL(ana_env%density_3d%conf_counter, KIND=dp), &
864 : ana_env%density_3d%sum_box_length(:)/ &
865 2000 : 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 500 : 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 2500 : REAL(ana_env%density_3d%conf_counter, KIND=dp))**2
874 500 : CALL close_file(unit_number=file_ptr)
875 : END IF
876 :
877 500 : DEALLOCATE (mass_bin)
878 : ! end the timing
879 500 : CALL timestop(handle)
880 1000 : 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 18 : 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 9 : CPASSERT(ASSOCIATED(ana_env))
901 9 : CPASSERT(ASSOCIATED(ana_env%density_3d))
902 9 : CPASSERT(ASSOCIATED(ana_env%density_3d%sum_density))
903 9 : CPASSERT(ASSOCIATED(ana_env%density_3d%sum_dens2))
904 :
905 : ! start the timing
906 9 : CALL timeset(routineN, handle)
907 :
908 : file_name = ""
909 9 : file_name_vari = ""
910 :
911 9 : bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
912 9 : bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
913 9 : bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
914 9 : CALL get_cell(cell=ana_env%cell, abc=cell_size)
915 9 : interval_size(1) = cell_size(1)/REAL(bin_x, KIND=dp)*au2a
916 9 : interval_size(2) = cell_size(2)/REAL(bin_y, KIND=dp)*au2a
917 9 : 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 9 : ana_env%temperature)
922 : CALL open_file(file_name=file_name, file_status="REPLACE", &
923 : file_action="WRITE", file_position="APPEND", &
924 9 : 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 9 : "# configurations", ana_env%density_3d%conf_counter, "bins", &
927 18 : ana_env%density_3d%nr_bins, "interval size", interval_size(:)
928 9 : 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 9 : ana_env%temperature)
934 : CALL open_file(file_name=file_name_vari, file_status="REPLACE", &
935 : file_action="WRITE", file_position="APPEND", &
936 9 : 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 9 : "# configurations", ana_env%density_3d%conf_counter, "bins", &
939 18 : ana_env%density_3d%nr_bins, "interval size", interval_size(:)
940 9 : WRITE (file_ptr_vari, FMT='(A,3A10,A20)') "#", " x [A] ", " y [A] ", " z [A] ", " variance"
941 :
942 27 : DO i = 1, SIZE(ana_env%density_3d%sum_density(:, 1, 1))
943 45 : DO j = 1, SIZE(ana_env%density_3d%sum_density(1, :, 1))
944 54 : DO k = 1, SIZE(ana_env%density_3d%sum_density(1, 1, :))
945 : WRITE (file_ptr_dens, FMT='(3F10.2,F20.10)') &
946 18 : (i - 0.5_dp)*interval_size(1), (j - 0.5_dp)*interval_size(2), (k - 0.5_dp)*interval_size(3), &
947 36 : 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 18 : (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 54 : (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 9 : CALL close_file(unit_number=file_ptr_dens)
956 9 : CALL close_file(unit_number=file_ptr_vari)
957 :
958 9 : WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
959 9 : WRITE (ana_env%io_unit, FMT="(T2,A,T35,A,T80,A)") "-", "density calculation", "-"
960 9 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
961 9 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "used configurations", &
962 18 : cp_to_string(REAL(ana_env%density_3d%conf_counter, KIND=dp))
963 9 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "average volume", &
964 : cp_to_string(ana_env%density_3d%sum_vol/ &
965 18 : REAL(ana_env%density_3d%conf_counter, KIND=dp))
966 9 : 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 81 : REAL(ana_env%density_3d%conf_counter, KIND=dp))
970 9 : 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 144 : REAL(ana_env%density_3d%conf_counter, KIND=dp))**2)
977 9 : WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
978 9 : IF (ana_env%print_test_output) &
979 9 : 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 81 : REAL(ana_env%density_3d%conf_counter, KIND=dp)
983 : ! end the timing
984 9 : CALL timestop(handle)
985 9 : 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 9 : 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 9 : TYPE(atom_pairs_type), DIMENSION(:), POINTER :: pairs_tmp
1009 :
1010 0 : CPASSERT(ASSOCIATED(ana_pair_correl))
1011 9 : CPASSERT(.NOT. ASSOCIATED(ana_pair_correl%g_r))
1012 9 : CPASSERT(.NOT. ASSOCIATED(ana_pair_correl%pairs))
1013 9 : CPASSERT(ASSOCIATED(atoms))
1014 9 : CPASSERT(SIZE(atoms) .GT. 1)
1015 9 : CPASSERT(ASSOCIATED(cell))
1016 :
1017 : ! start the timing
1018 9 : CALL timeset(routineN, handle)
1019 :
1020 9 : CALL get_cell(cell=cell, abc=cell_size)
1021 9 : IF (ana_pair_correl%nr_bins .LE. 0) THEN
1022 45 : 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 45 : ana_pair_correl%nr_bins
1026 9 : ana_pair_correl%conf_counter = 0
1027 :
1028 9 : counter = 1
1029 : ! initialise the atom pairs
1030 216 : ALLOCATE (pairs_tmp(SIZE(atoms)))
1031 198 : DO f_n = 1, SIZE(atoms)
1032 2088 : 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 1890 : n2=atoms(s_n)%name, list_end=counter - 1)
1036 : ! add to list
1037 2079 : IF (list_ind .LT. 0) THEN
1038 27 : pairs_tmp(counter)%f_n = atoms(f_n)%name
1039 27 : pairs_tmp(counter)%s_n = atoms(s_n)%name
1040 27 : pairs_tmp(counter)%pair_count = 1
1041 27 : counter = counter + 1
1042 : ELSE
1043 1863 : pairs_tmp(list_ind)%pair_count = pairs_tmp(list_ind)%pair_count + 1
1044 : END IF
1045 : END DO
1046 : END DO
1047 :
1048 54 : ALLOCATE (ana_pair_correl%pairs(counter - 1))
1049 36 : DO list = 1, counter - 1
1050 27 : ana_pair_correl%pairs(list)%f_n = pairs_tmp(list)%f_n
1051 27 : ana_pair_correl%pairs(list)%s_n = pairs_tmp(list)%s_n
1052 36 : ana_pair_correl%pairs(list)%pair_count = pairs_tmp(list)%pair_count
1053 : END DO
1054 9 : DEALLOCATE (pairs_tmp)
1055 :
1056 36 : ALLOCATE (ana_pair_correl%g_r(SIZE(ana_pair_correl%pairs(:)), ana_pair_correl%nr_bins))
1057 8145 : ana_pair_correl%g_r = 0.0_dp
1058 : ! end the timing
1059 9 : CALL timestop(handle)
1060 18 : 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 1000 : 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 500 : CPASSERT(ASSOCIATED(elem))
1084 500 : CPASSERT(ASSOCIATED(elem%pos))
1085 2000 : CPASSERT(ALL(elem%box_scale(:) .GT. 0.0_dp))
1086 500 : CPASSERT(weight .GT. 0)
1087 500 : CPASSERT(ASSOCIATED(atoms))
1088 500 : CPASSERT(ASSOCIATED(ana_env))
1089 500 : CPASSERT(ASSOCIATED(ana_env%cell))
1090 500 : CPASSERT(ASSOCIATED(ana_env%pair_correl))
1091 500 : CPASSERT(ASSOCIATED(ana_env%pair_correl%g_r))
1092 500 : CPASSERT(ASSOCIATED(ana_env%pair_correl%pairs))
1093 :
1094 : ! start the timing
1095 500 : CALL timeset(routineN, handle)
1096 :
1097 500 : dist = -1.0_dp
1098 :
1099 11000 : first_elem_loop: DO i = 1, SIZE(elem%pos), ana_env%dim_per_elem
1100 116000 : 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 105000 : cell=ana_env%cell, box_scale=elem%box_scale)
1104 105000 : ind = CEILING(dist/ana_env%pair_correl%step_length)
1105 115500 : 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 86745 : n2=atoms(INT(j/REAL(ana_env%dim_per_elem, KIND=dp)) + 1)%name)
1109 86745 : CPASSERT(pair_ind .GT. 0)
1110 : ana_env%pair_correl%g_r(pair_ind, ind) = &
1111 86745 : 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 500 : ana_env%pair_correl%conf_counter = ana_env%pair_correl%conf_counter + weight
1116 500 : 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 4000 : (elem%box_scale(:)*weight)
1119 : ! end the timing
1120 500 : CALL timestop(handle)
1121 500 : 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 18 : 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 9 : CPASSERT(ASSOCIATED(ana_env))
1140 9 : CPASSERT(ASSOCIATED(ana_env%pair_correl))
1141 :
1142 : ! start the timing
1143 9 : CALL timeset(routineN, handle)
1144 :
1145 9 : CALL get_cell(cell=ana_env%cell, abc=cell_size)
1146 36 : 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 9 : (cell_size(3)*aver_box_scale(3))
1150 :
1151 36 : 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 27 : 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 27 : unit_number=file_ptr)
1162 : WRITE (file_ptr, *) "# radial distribution function of "// &
1163 : TRIM(ana_env%pair_correl%pairs(pair)%f_n)//" and "// &
1164 27 : TRIM(ana_env%pair_correl%pairs(pair)%s_n)//" of ", &
1165 54 : ana_env%pair_correl%conf_counter, " configurations"
1166 27 : WRITE (file_ptr, *) "# using a bin size of ", &
1167 27 : ana_env%pair_correl%step_length*au2a, &
1168 54 : "[A] (for Vol changes: referring to the reference cell)"
1169 6129 : DO bin = 1, ana_env%pair_correl%nr_bins
1170 : voldr = 4.0/3.0*PI*ana_env%pair_correl%step_length**3* &
1171 6102 : (REAL(bin, KIND=dp)**3 - REAL(bin - 1, KIND=dp)**3)
1172 6102 : 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 12231 : (voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
1175 : END DO
1176 27 : CALL close_file(unit_number=file_ptr)
1177 :
1178 27 : IF (ana_env%print_test_output) &
1179 : WRITE (*, *) "TMC|ANALYSIS_G_R_"// &
1180 : TRIM(ana_env%pair_correl%pairs(pair)%f_n)//"_"// &
1181 27 : 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 6165 : voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
1184 : END DO
1185 :
1186 : ! end the timing
1187 9 : CALL timestop(handle)
1188 9 : 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 9 : 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 9 : CPASSERT(ASSOCIATED(ana_dip_mom))
1210 9 : CPASSERT(ASSOCIATED(ana_dip_mom%charges_inp))
1211 9 : CPASSERT(ASSOCIATED(atoms))
1212 :
1213 : ! start the timing
1214 9 : CALL timeset(routineN, handle)
1215 :
1216 27 : ALLOCATE (ana_dip_mom%charges(SIZE(atoms)))
1217 198 : ana_dip_mom%charges = 0.0_dp
1218 : ! for every atom searcht the correct charge
1219 198 : DO atom = 1, SIZE(atoms)
1220 324 : charge_loop: DO charge = 1, SIZE(ana_dip_mom%charges_inp)
1221 315 : IF (atoms(atom)%name .EQ. ana_dip_mom%charges_inp(charge)%name) THEN
1222 189 : ana_dip_mom%charges(atom) = ana_dip_mom%charges_inp(charge)%mass
1223 189 : EXIT charge_loop
1224 : END IF
1225 : END DO charge_loop
1226 : END DO
1227 :
1228 9 : DEALLOCATE (ana_dip_mom%charges_inp)
1229 : ! end the timing
1230 9 : CALL timestop(handle)
1231 9 : 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 500 : 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 500 : REAL(KIND=dp), DIMENSION(:), POINTER :: dip_cl
1251 :
1252 0 : CPASSERT(ASSOCIATED(elem))
1253 500 : CPASSERT(ASSOCIATED(elem%pos))
1254 500 : CPASSERT(ASSOCIATED(ana_env))
1255 500 : CPASSERT(ASSOCIATED(ana_env%dip_mom))
1256 500 : CPASSERT(ASSOCIATED(ana_env%dip_mom%charges))
1257 :
1258 : ! start the timing
1259 500 : CALL timeset(routineN, handle)
1260 :
1261 1500 : ALLOCATE (dip_cl(ana_env%dim_per_elem))
1262 2000 : dip_cl(:) = 0.0_dp
1263 :
1264 500 : 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 84000 : 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 500 : IF (.NOT. ASSOCIATED(elem%dipole)) THEN
1271 1500 : ALLOCATE (elem%dipole(ana_env%dim_per_elem))
1272 4000 : elem%dipole(:) = dip_cl(:)
1273 : END IF
1274 :
1275 500 : IF (ana_env%dip_mom%print_cl_dip) THEN
1276 : file_name = expand_file_name_temp(tmc_default_trajectory_file_name, &
1277 500 : 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 500 : file_ext="dip_cl")
1281 : END IF
1282 500 : ana_env%dip_mom%conf_counter = ana_env%dip_mom%conf_counter + weight
1283 4000 : ana_env%dip_mom%last_dip_cl(:) = dip_cl
1284 :
1285 500 : DEALLOCATE (dip_cl)
1286 :
1287 : ! end the timing
1288 500 : CALL timestop(handle)
1289 1000 : 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 9 : SUBROUTINE print_dipole_moment(ana_env)
1298 : TYPE(tmc_analysis_env), POINTER :: ana_env
1299 :
1300 9 : IF (ana_env%print_test_output) &
1301 9 : WRITE (*, *) "TMC|ANALYSIS_FINAL_CLASS_CELL_DIPOLE_MOMENT_X= ", &
1302 18 : ana_env%dip_mom%last_dip_cl(:)
1303 9 : 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 0 : 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 0 : NULLIFY (scaled_cell)
1323 :
1324 0 : CPASSERT(ASSOCIATED(elem))
1325 0 : CPASSERT(ASSOCIATED(elem%dipole))
1326 0 : CPASSERT(ASSOCIATED(ana_env))
1327 0 : CPASSERT(ASSOCIATED(ana_env%dip_ana))
1328 :
1329 0 : weight_act = weight
1330 0 : IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
1331 0 : weight_act = weight_act/REAL(8.0, KIND=dp)
1332 :
1333 : ! get the volume
1334 0 : ALLOCATE (scaled_cell)
1335 : CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, vol=vol, &
1336 0 : scaled_cell=scaled_cell)
1337 :
1338 : ! fold exact dipole moments using the classical ones
1339 0 : IF (ASSOCIATED(ana_env%dip_mom)) THEN
1340 0 : 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 0 : cell=scaled_cell) + ana_env%dip_mom%last_dip_cl
1343 : END IF
1344 : END IF
1345 :
1346 0 : 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 0 : DOT_PRODUCT(elem%dipole(:), elem%dipole(:))/vol*weight_act
1351 :
1352 0 : tmp_dip(:, :) = 0.0_dp
1353 0 : 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 0 : 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 0 : 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 0 : tmp_dip(:, 1)**2/vol*weight_act
1366 :
1367 : ! calculate the directional average with componentwise correlation per volume
1368 0 : tmp_dip(:, :) = MATMUL(tmp_dip(:, :), TRANSPOSE(tmp_dip(:, :)))
1369 : ana_env%dip_ana%mu2_pv_mat(:, :) = ana_env%dip_ana%mu2_pv_mat(:, :) + &
1370 0 : tmp_dip(:, :)/vol*weight_act
1371 :
1372 0 : 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 0 : 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 0 : kB = boltzmann/joule
1393 0 : counter_tmp = INT(ana_env%dip_ana%conf_counter)
1394 :
1395 : ! TODO get correct constant using physcon
1396 0 : e0 = 0.07957747154594767_dp !e^2*a0*me*hbar^-2
1397 0 : 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 0 : 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 0 : 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 0 : ana_env%temperature)
1410 :
1411 0 : 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 0 : "diel_const"))
1415 : file_name_tmp = TRIM(expand_file_name_char(file_name_tmp, &
1416 0 : "diel_const_tensor"))
1417 : CASE (ana_type_sym_xyz)
1418 : file_name = TRIM(expand_file_name_char(file_name_tmp, &
1419 0 : "diel_const_sym"))
1420 : file_name_tmp = TRIM(expand_file_name_char(file_name_tmp, &
1421 0 : "diel_const_tensor_sym"))
1422 : CASE DEFAULT
1423 0 : 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 0 : 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 0 : diel_const_norm
1436 : ! print dielectric constant trajectory
1437 : ! if szmetry used print only every 8th configuration, hence every different (not mirrowed)
1438 0 : 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 0 : unit_number=file_ptr)
1442 0 : IF (.NOT. flag) THEN
1443 0 : WRITE (file_ptr, FMT='(A8,5A20)') "# conf", "diel_const", &
1444 0 : "diel_const_sym", "diel_const_sym_x", &
1445 0 : "diel_const_sym_y", "diel_const_sym_z"
1446 : END IF
1447 0 : WRITE (file_ptr, FMT="(I8,10F20.10)") counter_tmp, diel_const, &
1448 0 : diel_const_sym, &
1449 : 4.0_dp*PI/(kB*ana_env%temperature)* &
1450 0 : ana_env%dip_ana%mu2_pv(:)/REAL(ana_env%dip_ana%conf_counter, KIND=dp)
1451 0 : CALL close_file(unit_number=file_ptr)
1452 :
1453 : ! print dielectric constant tensor trajectory
1454 0 : 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 0 : unit_number=file_ptr)
1458 0 : IF (.NOT. flag) THEN
1459 0 : WRITE (file_ptr, FMT='(A8,9A20)') "# conf", "xx", "xy", "xz", &
1460 0 : "yx", "yy", "yz", &
1461 0 : "zx", "zy", "zz"
1462 : END IF
1463 0 : tmp_dip(:, :) = 0.0_dp
1464 0 : tmp_dip(:, 1) = ana_env%dip_ana%mu_psv(:)/REAL(ana_env%dip_ana%conf_counter, KIND=dp)
1465 :
1466 0 : 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 0 : MATMUL(tmp_dip(:, :), TRANSPOSE(tmp_dip(:, :))))
1470 0 : CALL close_file(unit_number=file_ptr)
1471 0 : 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 0 : 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 0 : kB = boltzmann/joule
1490 :
1491 0 : CPASSERT(ASSOCIATED(ana_env))
1492 0 : CPASSERT(ASSOCIATED(ana_env%dip_ana))
1493 :
1494 0 : tmp_dip(:, :) = 0.0_dp
1495 : diel_const(:, :) = 0.0_dp
1496 0 : diel_const_scalar = 0.0_dp
1497 0 : diel_const_sym = 0.0_dp
1498 :
1499 : !dielectric constant
1500 0 : 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 0 : 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 0 : ana_env%dip_ana%mu2_pv(:)/REAL(ana_env%dip_ana%conf_counter, KIND=dp)
1508 :
1509 0 : DO i = 1, 3
1510 0 : diel_const(i, i) = diel_const(i, i) + 1.0_dp ! +1 for unpolarizable models, 1.592 for polarizable
1511 0 : diel_const_scalar = diel_const_scalar + diel_const(i, i)
1512 : END DO
1513 0 : diel_const_scalar = diel_const_scalar/REAL(3, KIND=dp)
1514 :
1515 0 : tmp_dip(:, :) = diel_const
1516 0 : CALL diag(3, tmp_dip, dielec_ev, tmp_ev)
1517 :
1518 : ! print out results
1519 0 : WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
1520 0 : WRITE (ana_env%io_unit, FMT="(T2,A,T35,A,T80,A)") "-", "average dipoles", "-"
1521 0 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
1522 0 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "used configurations ", &
1523 0 : cp_to_string(REAL(ana_env%dip_ana%conf_counter, KIND=dp))
1524 0 : IF (ana_env%dip_ana%ana_type .EQ. ana_type_ice) &
1525 0 : WRITE (ana_env%io_unit, FMT='(T2,A,"| ",A)') plabel, &
1526 0 : "ice analysis with directions of hexagonal structure"
1527 0 : IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
1528 0 : WRITE (ana_env%io_unit, FMT='(T2,A,"| ",A)') plabel, &
1529 0 : "ice analysis with symmetrized dipoles in each direction."
1530 :
1531 0 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "for product of 2 directions(per vol):"
1532 0 : DO i = 1, 3
1533 0 : WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", ana_env%dip_ana%mu2_pv_mat(i, :)/ &
1534 0 : REAL(ana_env%dip_ana%conf_counter, KIND=dp), " |"
1535 : END DO
1536 :
1537 0 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "dielectric constant tensor:"
1538 0 : DO i = 1, 3
1539 0 : WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", diel_const(i, :), " |"
1540 : END DO
1541 :
1542 0 : 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 0 : cp_to_string(dielec_ev(3))
1546 0 : 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 0 : cp_to_string(diel_const_sym(3))
1550 0 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "dielectric constant ", &
1551 0 : cp_to_string(diel_const_scalar)
1552 0 : WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
1553 :
1554 0 : 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 1000 : 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 500 : disp = 0.0_dp
1580 :
1581 500 : CPASSERT(ASSOCIATED(elem))
1582 500 : CPASSERT(ASSOCIATED(elem%pos))
1583 500 : CPASSERT(ASSOCIATED(ana_env))
1584 500 : CPASSERT(ASSOCIATED(ana_env%displace))
1585 500 : CPASSERT(ASSOCIATED(ana_env%last_elem))
1586 :
1587 : ! start the timing
1588 500 : CALL timeset(routineN, handle)
1589 :
1590 500 : DO ind = 1, SIZE(elem%pos), ana_env%dim_per_elem
1591 : ! fold into box
1592 42000 : 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 10500 : vec=atom_disp)
1595 42000 : disp = disp + SUM((atom_disp(:)*au2a)**2)
1596 : END DO
1597 500 : ana_env%displace%disp = ana_env%displace%disp + disp
1598 500 : ana_env%displace%conf_counter = ana_env%displace%conf_counter + 1
1599 :
1600 500 : 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 500 : ana_env%temperature)
1604 : file_name = TRIM(expand_file_name_char(file_name_tmp, &
1605 500 : "devi"))
1606 500 : 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 500 : unit_number=file_ptr)
1610 500 : IF (.NOT. flag) &
1611 3 : WRITE (file_ptr, *) "# conf squared deviation of the cell"
1612 500 : WRITE (file_ptr, *) elem%nr, disp
1613 500 : CALL close_file(unit_number=file_ptr)
1614 : END IF
1615 :
1616 : ! end the timing
1617 500 : CALL timestop(handle)
1618 :
1619 500 : 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 9 : 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 9 : WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
1633 9 : WRITE (ana_env%io_unit, FMT="(T2,A,T35,A,T80,A)") "-", "average displacement", "-"
1634 9 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "temperature ", &
1635 18 : cp_to_string(ana_env%temperature)
1636 9 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "used configurations ", &
1637 18 : cp_to_string(REAL(ana_env%displace%conf_counter, KIND=dp))
1638 9 : WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "cell root mean square deviation: ", &
1639 : cp_to_string(SQRT(ana_env%displace%disp/ &
1640 18 : REAL(ana_env%displace%conf_counter, KIND=dp)))
1641 9 : IF (ana_env%print_test_output) &
1642 9 : WRITE (*, *) "TMC|ANALYSIS_AVERAGE_CELL_DISPLACEMENT_X= ", &
1643 : SQRT(ana_env%displace%disp/ &
1644 18 : REAL(ana_env%displace%conf_counter, KIND=dp))
1645 9 : END SUBROUTINE print_average_displacement
1646 : END MODULE tmc_analysis
|