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