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