Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief I/O subroutines for helium
10 : !> \author Lukasz Walewski
11 : !> \date 2009-06-08
12 : ! **************************************************************************************************
13 : MODULE helium_io
14 :
15 : USE cell_types, ONLY: get_cell
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_get_default_unit_nr,&
18 : cp_logger_type
19 : USE cp_output_handling, ONLY: cp_p_file,&
20 : cp_print_key_finished_output,&
21 : cp_print_key_should_output,&
22 : cp_print_key_unit_nr
23 : USE cp_parser_methods, ONLY: parser_get_next_line,&
24 : parser_get_object
25 : USE cp_parser_types, ONLY: cp_parser_type,&
26 : parser_create,&
27 : parser_release
28 : USE cp_units, ONLY: cp_unit_from_cp2k,&
29 : cp_unit_to_cp2k
30 : USE helium_common, ONLY: helium_cycle_number,&
31 : helium_cycle_of,&
32 : helium_is_winding,&
33 : helium_path_length,&
34 : helium_pbc
35 : USE helium_interactions, ONLY: helium_total_inter_action,&
36 : helium_total_link_action,&
37 : helium_total_pair_action
38 : USE helium_types, ONLY: &
39 : e_id_interact, e_id_kinetic, e_id_potential, e_id_thermo, e_id_total, e_id_virial, &
40 : helium_solvent_p_type, helium_solvent_type, rho_num
41 : USE input_constants, ONLY: &
42 : fmt_id_pdb, fmt_id_xyz, helium_cell_shape_cube, helium_cell_shape_octahedron, &
43 : helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, &
44 : helium_solute_intpot_nnp, helium_solute_intpot_none, perm_cycle, perm_plain
45 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
46 : section_vals_type,&
47 : section_vals_val_get
48 : USE kinds, ONLY: default_path_length,&
49 : default_string_length,&
50 : dp,&
51 : int_8
52 : USE machine, ONLY: m_flush
53 : USE memory_utilities, ONLY: reallocate
54 : USE message_passing, ONLY: mp_para_env_type
55 : USE physcon, ONLY: angstrom,&
56 : massunit
57 : USE pint_types, ONLY: pint_env_type
58 : #include "../base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 :
62 : PRIVATE
63 :
64 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_io'
66 :
67 : CHARACTER(len=*), DIMENSION(6), PARAMETER, PRIVATE :: m_dist_name = (/ &
68 : "SINGLEV ", &
69 : "UNIFORM ", &
70 : "LINEAR ", &
71 : "QUADRATIC ", &
72 : "EXPONENTIAL", &
73 : "GAUSSIAN "/)
74 :
75 : PUBLIC :: helium_read_xyz
76 : PUBLIC :: helium_print_rdf
77 : PUBLIC :: helium_print_rho
78 : PUBLIC :: helium_write_line
79 : PUBLIC :: helium_write_setup
80 : PUBLIC :: helium_print_energy
81 : PUBLIC :: helium_print_vector
82 : PUBLIC :: helium_print_plength
83 : PUBLIC :: helium_print_coordinates
84 : PUBLIC :: helium_print_force
85 : PUBLIC :: helium_print_accepts
86 : PUBLIC :: helium_print_perm
87 : PUBLIC :: helium_print_action
88 : PUBLIC :: helium_write_cubefile
89 :
90 : CONTAINS
91 :
92 : ! ***************************************************************************
93 : !> \brief Read XYZ coordinates from file
94 : !> \param coords ...
95 : !> \param file_name ...
96 : !> \param para_env ...
97 : !> \date 2009-06-03
98 : !> \author Lukasz Walewski
99 : !> \note This is a parallel routine, all ranks get coords defined
100 : ! **************************************************************************************************
101 0 : SUBROUTINE helium_read_xyz(coords, file_name, para_env)
102 :
103 : REAL(KIND=dp), DIMENSION(:), POINTER :: coords
104 : CHARACTER(LEN=default_path_length) :: file_name
105 : TYPE(mp_para_env_type), POINTER :: para_env
106 :
107 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_read_xyz'
108 :
109 : CHARACTER(LEN=default_string_length) :: strtmp
110 : INTEGER :: frame, handle, istat, j, natom
111 : LOGICAL :: found, my_end
112 : TYPE(cp_parser_type) :: parser
113 :
114 0 : CALL timeset(routineN, handle)
115 :
116 : ! check if the file can be accessed
117 0 : INQUIRE (FILE=file_name, EXIST=found, IOSTAT=istat)
118 0 : IF (istat /= 0) THEN
119 : WRITE (UNIT=strtmp, FMT="(A,I0,A)") &
120 : "An error occurred inquiring the file <"// &
121 0 : TRIM(file_name)//">"
122 0 : CPABORT(TRIM(strtmp))
123 : END IF
124 0 : IF (.NOT. found) THEN
125 0 : CPABORT("Coordinate file <"//TRIM(file_name)//"> not found.")
126 : END IF
127 :
128 : CALL parser_create( &
129 : parser, &
130 : file_name, &
131 : para_env=para_env, &
132 0 : parse_white_lines=.TRUE.)
133 :
134 : natom = 0
135 0 : frame = 0
136 0 : CALL parser_get_next_line(parser, 1)
137 0 : Frames: DO
138 : ! Atom number
139 0 : CALL parser_get_object(parser, natom)
140 0 : frame = frame + 1
141 0 : IF (frame == 1) THEN
142 0 : ALLOCATE (coords(3*natom))
143 : ELSE
144 0 : strtmp = "Warning: more than one frame on file <"//TRIM(file_name)//">"
145 0 : CALL helium_write_line(strtmp)
146 0 : strtmp = "Warning: using the first frame only!"
147 0 : CALL helium_write_line(strtmp)
148 0 : EXIT Frames
149 : END IF
150 : ! Dummy line
151 0 : CALL parser_get_next_line(parser, 2)
152 0 : DO j = 1, natom
153 : ! Atom coordinates
154 0 : READ (parser%input_line, *) strtmp, &
155 0 : coords(3*(j - 1) + 1), &
156 0 : coords(3*(j - 1) + 2), &
157 0 : coords(3*(j - 1) + 3)
158 0 : coords(3*(j - 1) + 1) = cp_unit_to_cp2k(coords(3*(j - 1) + 1), "angstrom")
159 0 : coords(3*(j - 1) + 2) = cp_unit_to_cp2k(coords(3*(j - 1) + 2), "angstrom")
160 0 : coords(3*(j - 1) + 3) = cp_unit_to_cp2k(coords(3*(j - 1) + 3), "angstrom")
161 : ! If there's a white line or end of file exit.. otherwise go on
162 0 : CALL parser_get_next_line(parser, 1, at_end=my_end)
163 0 : my_end = my_end .OR. (LEN_TRIM(parser%input_line) == 0)
164 0 : IF (my_end) THEN
165 0 : IF (j /= natom) THEN
166 0 : CPABORT("Error in XYZ format.")
167 : END IF
168 : EXIT Frames
169 : END IF
170 : END DO
171 : END DO Frames
172 0 : CALL parser_release(parser)
173 :
174 0 : CALL timestop(handle)
175 :
176 0 : END SUBROUTINE helium_read_xyz
177 :
178 : ! ***************************************************************************
179 : !> \brief Write helium parameters to the output unit
180 : !> \param helium ...
181 : !> \date 2009-06-03
182 : !> \par History
183 : !> 2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
184 : !> \author Lukasz Walewski
185 : ! **************************************************************************************************
186 25 : SUBROUTINE helium_write_setup(helium)
187 :
188 : TYPE(helium_solvent_type), POINTER :: helium
189 :
190 : CHARACTER(len=default_string_length) :: msg_str, my_label, stmp, stmp1, stmp2, &
191 : unit_str
192 : INTEGER :: i, itmp, unit_nr
193 : INTEGER(KIND=int_8) :: i8tmp
194 : REAL(KIND=dp) :: rtmp, v1, v2, v3
195 : REAL(KIND=dp), DIMENSION(3) :: my_abc
196 : TYPE(cp_logger_type), POINTER :: logger
197 :
198 25 : NULLIFY (logger)
199 25 : logger => cp_get_default_logger()
200 25 : my_label = "HELIUM| "
201 :
202 25 : IF (logger%para_env%is_source()) THEN
203 13 : unit_nr = cp_logger_get_default_unit_nr(logger)
204 :
205 13 : WRITE (unit_nr, *)
206 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
207 13 : " Number of helium environments: ", helium%num_env
208 :
209 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
210 13 : " Number of solvent atoms: ", helium%atoms
211 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
212 13 : " Number of solvent beads: ", helium%beads
213 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
214 13 : " Total number of solvent particles: ", helium%atoms*helium%beads
215 :
216 13 : unit_str = "angstrom^-3"
217 : rtmp = cp_unit_from_cp2k(helium%density, &
218 13 : unit_str)
219 : WRITE (unit_nr, '(T2,A,F12.6)') TRIM(my_label)//" Density ["// &
220 13 : TRIM(unit_str)//"]:", rtmp
221 :
222 13 : unit_str = "a.m.u."
223 13 : rtmp = helium%he_mass_au/massunit
224 : WRITE (unit_nr, '(T2,A,F12.6)') TRIM(my_label)//" He Mass ["// &
225 13 : TRIM(unit_str)//"]: ", rtmp
226 :
227 13 : unit_str = "angstrom"
228 : rtmp = cp_unit_from_cp2k(helium%cell_size, &
229 13 : unit_str)
230 : WRITE (unit_nr, '(T2,A,F12.6)') TRIM(my_label)//" Cell size ["// &
231 13 : TRIM(unit_str)//"]: ", rtmp
232 :
233 13 : IF (helium%periodic) THEN
234 8 : SELECT CASE (helium%cell_shape)
235 : CASE (helium_cell_shape_cube)
236 0 : CALL helium_write_line("PBC cell shape: CUBE.")
237 : CASE (helium_cell_shape_octahedron)
238 8 : CALL helium_write_line("PBC cell shape: TRUNCATED OCTAHEDRON.")
239 : CASE DEFAULT
240 8 : CALL helium_write_line("*** Warning: unknown cell shape.")
241 : END SELECT
242 : ELSE
243 5 : CALL helium_write_line("PBC turned off.")
244 : END IF
245 :
246 13 : IF (helium%droplet_radius .LT. HUGE(1.0_dp)) THEN
247 5 : WRITE (stmp, *) helium%droplet_radius*angstrom
248 5 : CALL helium_write_line("Droplet radius: "//TRIM(ADJUSTL(stmp))//" [angstrom]")
249 : END IF
250 :
251 : ! first step gets incremented during first iteration
252 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
253 13 : " First MC step :", helium%first_step + 1
254 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
255 13 : " Last MC step :", helium%last_step
256 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
257 13 : " Total number of MC steps :", helium%num_steps
258 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
259 13 : " Number of outer MC trials per step :", helium%iter_rot
260 13 : i8tmp = INT(helium%iter_rot, int_8)
261 13 : IF (helium%sampling_method == helium_sampling_ceperley) THEN
262 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
263 8 : " Number of inner MC trials per step :", helium%iter_norot
264 8 : i8tmp = i8tmp*INT(helium%iter_norot, int_8)
265 : END IF
266 13 : stmp = ""
267 13 : WRITE (stmp, *) i8tmp
268 : WRITE (unit_nr, '(T2,A)') TRIM(my_label)// &
269 13 : " Total number of MC trials per step : "//TRIM(ADJUSTL(stmp))
270 13 : i8tmp = INT(helium%num_steps, int_8)
271 13 : i8tmp = i8tmp*INT(helium%iter_rot, int_8)
272 13 : IF (helium%sampling_method == helium_sampling_ceperley) THEN
273 8 : i8tmp = i8tmp*INT(helium%iter_norot, int_8)
274 : END IF
275 13 : stmp = ""
276 13 : WRITE (stmp, *) i8tmp
277 : WRITE (unit_nr, '(T2,A)') TRIM(my_label)// &
278 13 : " Total number of MC trials : "//TRIM(ADJUSTL(stmp))
279 :
280 21 : SELECT CASE (helium%sampling_method)
281 :
282 : CASE (helium_sampling_ceperley)
283 :
284 : ! permutation cycle length sampling
285 8 : stmp = ""
286 8 : CALL helium_write_line(stmp)
287 8 : WRITE (stmp, *) helium%maxcycle
288 8 : stmp2 = ""
289 : WRITE (stmp2, *) "Using maximum permutation cycle length: "// &
290 8 : TRIM(ADJUSTL(stmp))
291 8 : CALL helium_write_line(stmp2)
292 8 : stmp = ""
293 : WRITE (stmp, *) "Permutation cycle length distribution: "// &
294 8 : TRIM(ADJUSTL(m_dist_name(helium%m_dist_type)))
295 8 : CALL helium_write_line(stmp)
296 8 : stmp = ""
297 8 : stmp1 = ""
298 8 : WRITE (stmp1, *) helium%m_ratio
299 8 : stmp2 = ""
300 8 : WRITE (stmp2, *) helium%m_value
301 : WRITE (stmp, *) "Using ratio "//TRIM(ADJUSTL(stmp1))// &
302 8 : " for M = "//TRIM(ADJUSTL(stmp2))
303 8 : CALL helium_write_line(stmp)
304 8 : stmp = ""
305 8 : CALL helium_write_line(stmp)
306 :
307 : CASE (helium_sampling_worm)
308 :
309 5 : stmp1 = ""
310 5 : stmp2 = ""
311 5 : CALL helium_write_line(stmp1)
312 5 : WRITE (stmp1, *) helium%worm_centroid_drmax
313 : WRITE (stmp2, *) "WORM| Centroid move max. displacement: "// &
314 5 : TRIM(ADJUSTL(stmp1))
315 5 : CALL helium_write_line(stmp2)
316 5 : stmp1 = ""
317 5 : stmp2 = ""
318 5 : WRITE (stmp1, *) helium%worm_staging_l
319 5 : WRITE (stmp2, *) "WORM| Maximal staging length: "//TRIM(ADJUSTL(stmp1))
320 5 : CALL helium_write_line(stmp2)
321 5 : stmp1 = ""
322 5 : stmp2 = ""
323 5 : WRITE (stmp1, *) helium%worm_open_close_scale
324 5 : WRITE (stmp2, *) "WORM| Open/Close scaling: "//TRIM(ADJUSTL(stmp1))
325 5 : CALL helium_write_line(stmp2)
326 5 : stmp1 = ""
327 5 : stmp2 = ""
328 5 : WRITE (stmp1, *) helium%worm_allow_open
329 5 : WRITE (stmp2, *) "WORM| Open worm sector: "//TRIM(ADJUSTL(stmp1))
330 5 : CALL helium_write_line(stmp2)
331 5 : stmp1 = ""
332 5 : stmp2 = ""
333 5 : WRITE (stmp1, *) helium%worm_show_statistics
334 5 : WRITE (stmp2, *) "WORM| Print statistics: "//TRIM(ADJUSTL(stmp1))
335 5 : CALL helium_write_line(stmp2)
336 5 : IF (helium%worm_max_open_cycles > 0 .AND. helium%worm_allow_open) THEN
337 0 : stmp1 = ""
338 0 : stmp2 = ""
339 0 : WRITE (stmp1, *) helium%worm_max_open_cycles
340 0 : WRITE (stmp2, *) "WORM| Max. nb of MC cycles in open sector: "//TRIM(ADJUSTL(stmp1))
341 0 : CALL helium_write_line(stmp2)
342 : END IF
343 5 : stmp1 = ""
344 5 : CALL helium_write_line(stmp1)
345 :
346 : CASE DEFAULT
347 0 : WRITE (msg_str, *) helium%sampling_method
348 0 : msg_str = "Sampling method ("//TRIM(ADJUSTL(msg_str))//") not supported."
349 13 : CPABORT(msg_str)
350 :
351 : END SELECT
352 :
353 : ! solute related data
354 13 : IF (helium%solute_present) THEN
355 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
356 8 : " Number of solute atoms: ", helium%solute_atoms
357 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
358 8 : " Number of solute beads: ", helium%solute_beads
359 : WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
360 8 : " Total number of solute particles: ", helium%solute_atoms* &
361 16 : helium%solute_beads
362 :
363 8 : stmp1 = ""
364 8 : SELECT CASE (helium%solute_interaction)
365 : CASE (helium_solute_intpot_none)
366 0 : WRITE (stmp1, *) "NONE"
367 : CASE (helium_solute_intpot_mwater)
368 7 : WRITE (stmp1, *) "MWATER"
369 : CASE (helium_solute_intpot_nnp)
370 1 : WRITE (stmp1, *) "NNP"
371 : CASE DEFAULT
372 8 : WRITE (stmp1, *) "***UNKNOWN***"
373 : END SELECT
374 : WRITE (unit_nr, '(T2,A,A,A)') &
375 8 : TRIM(my_label), &
376 8 : " Solute interaction type: ", &
377 16 : TRIM(ADJUSTL(stmp1))
378 :
379 8 : CALL get_cell(helium%solute_cell, abc=my_abc)
380 8 : unit_str = "angstrom"
381 8 : v1 = cp_unit_from_cp2k(my_abc(1), unit_str)
382 8 : v2 = cp_unit_from_cp2k(my_abc(2), unit_str)
383 8 : v3 = cp_unit_from_cp2k(my_abc(3), unit_str)
384 : WRITE (unit_nr, '(T2,A,F12.6,1X,F12.6,1X,F12.6)') &
385 : TRIM(my_label)//" Solute cell size ["// &
386 8 : TRIM(unit_str)//"]: ", v1, v2, v3
387 : ELSE
388 5 : WRITE (unit_nr, '(T2,A)') TRIM(my_label)//" Solute is NOT present"
389 : END IF
390 : END IF
391 25 : CALL helium_write_line("")
392 :
393 : ! RDF related settings
394 25 : IF (helium%rdf_present) THEN
395 2 : WRITE (stmp, '(I6)') helium%rdf_num
396 2 : CALL helium_write_line("RDF| number of centers: "//TRIM(stmp))
397 2 : rtmp = cp_unit_from_cp2k(helium%rdf_delr, "angstrom")
398 2 : WRITE (stmp, '(1X,F12.6)') rtmp
399 2 : CALL helium_write_line("RDF| delr [angstrom] : "//TRIM(stmp))
400 2 : rtmp = cp_unit_from_cp2k(helium%rdf_maxr, "angstrom")
401 2 : WRITE (stmp, '(1X,F12.6)') rtmp
402 2 : CALL helium_write_line("RDF| maxr [angstrom] : "//TRIM(stmp))
403 2 : itmp = helium%rdf_nbin
404 2 : WRITE (stmp, '(I6)') itmp
405 2 : CALL helium_write_line("RDF| nbin : "//TRIM(stmp))
406 2 : CALL helium_write_line("")
407 : ELSE
408 23 : CALL helium_write_line("RDF| radial distributions will NOT be calculated.")
409 23 : CALL helium_write_line("")
410 : END IF
411 :
412 : ! RHO related settings
413 25 : IF (helium%rho_present) THEN
414 2 : CALL helium_write_line('RHO| The following densities will be calculated:')
415 12 : DO i = 1, rho_num
416 12 : IF (helium%rho_property(i)%is_calculated) THEN
417 2 : WRITE (stmp, '(A)') 'RHO| '//TRIM(helium%rho_property(i)%name)
418 2 : CALL helium_write_line(TRIM(stmp))
419 : END IF
420 : END DO
421 2 : CALL helium_write_line('RHO|')
422 2 : rtmp = cp_unit_from_cp2k(helium%rho_delr, "angstrom")
423 2 : WRITE (stmp, '(1X,F12.6)') rtmp
424 2 : CALL helium_write_line("RHO| delr [angstrom] : "//TRIM(stmp))
425 2 : rtmp = cp_unit_from_cp2k(helium%rho_maxr, "angstrom")
426 2 : WRITE (stmp, '(1X,F12.6)') rtmp
427 2 : CALL helium_write_line("RHO| maxr [angstrom] : "//TRIM(stmp))
428 2 : itmp = helium%rho_nbin
429 2 : WRITE (stmp, '(I6)') itmp
430 2 : CALL helium_write_line("RHO| nbin : "//TRIM(stmp))
431 2 : CALL helium_write_line("")
432 : ELSE
433 23 : CALL helium_write_line("RHO| Density distributions will NOT be calculated.")
434 23 : CALL helium_write_line("")
435 : END IF
436 :
437 25 : RETURN
438 : END SUBROUTINE helium_write_setup
439 :
440 : ! ***************************************************************************
441 : !> \brief Writes out a line of text to the default output unit
442 : !> \param line ...
443 : !> \date 2009-07-10
444 : !> \author Lukasz Walewski
445 : ! **************************************************************************************************
446 897 : SUBROUTINE helium_write_line(line)
447 :
448 : CHARACTER(len=*), INTENT(IN) :: line
449 :
450 : CHARACTER(len=default_string_length) :: my_label
451 : INTEGER :: unit_nr
452 : TYPE(cp_logger_type), POINTER :: logger
453 :
454 897 : NULLIFY (logger)
455 897 : logger => cp_get_default_logger()
456 897 : my_label = "HELIUM|"
457 :
458 897 : IF (logger%para_env%is_source()) THEN
459 537 : unit_nr = cp_logger_get_default_unit_nr(logger)
460 537 : WRITE (unit_nr, '(T2,A)') TRIM(my_label)//" "//TRIM(line)
461 : END IF
462 :
463 897 : RETURN
464 : END SUBROUTINE helium_write_line
465 :
466 : ! ***************************************************************************
467 : !> \brief Print energies according to HELIUM%PRINT%ENERGY
468 : !> \param helium_env ...
469 : !> \date 2009-06-08
470 : !> \par History
471 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
472 : !> \author Lukasz Walewski
473 : ! **************************************************************************************************
474 105 : SUBROUTINE helium_print_energy(helium_env)
475 :
476 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
477 :
478 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_energy'
479 :
480 : INTEGER :: handle, iteration, k, m, unit_nr
481 : LOGICAL :: cepsample, file_is_new, should_output
482 : REAL(KIND=dp) :: naccptd, ntot
483 105 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_energy
484 : TYPE(cp_logger_type), POINTER :: logger
485 : TYPE(section_vals_type), POINTER :: print_key
486 :
487 105 : CALL timeset(routineN, handle)
488 :
489 105 : NULLIFY (logger, print_key)
490 105 : logger => cp_get_default_logger()
491 :
492 105 : IF (logger%para_env%is_source()) THEN
493 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
494 55 : "MOTION%PINT%HELIUM%PRINT%ENERGY")
495 : should_output = BTEST(cp_print_key_should_output( &
496 : iteration_info=logger%iter_info, &
497 55 : basis_section=print_key), cp_p_file)
498 55 : cepsample = helium_env(1)%helium%sampling_method == helium_sampling_ceperley
499 : END IF
500 105 : CALL helium_env(1)%comm%bcast(should_output, logger%para_env%source)
501 105 : CALL helium_env(1)%comm%bcast(cepsample, logger%para_env%source)
502 :
503 105 : IF (should_output) THEN
504 :
505 105 : naccptd = 0.0_dp
506 105 : IF (cepsample) THEN
507 80 : ntot = 0.0_dp
508 190 : DO k = 1, SIZE(helium_env)
509 110 : ntot = ntot + helium_env(1)%helium%iter_norot*helium_env(1)%helium%iter_rot
510 480 : DO m = 1, helium_env(k)%helium%maxcycle
511 400 : naccptd = naccptd + helium_env(k)%helium%num_accepted(helium_env(k)%helium%bisctlog2 + 2, m)
512 : END DO
513 : END DO
514 : ELSE !(wormsample)
515 25 : ntot = 0.0_dp
516 50 : DO k = 1, SIZE(helium_env)
517 25 : naccptd = naccptd + helium_env(k)%helium%num_accepted(1, 1)
518 50 : ntot = ntot + helium_env(k)%helium%num_accepted(2, 1)
519 : END DO
520 : END IF
521 105 : CALL helium_env(1)%comm%sum(naccptd)
522 105 : CALL helium_env(1)%comm%sum(ntot)
523 :
524 105 : IF (logger%para_env%is_source()) THEN
525 55 : my_energy => helium_env(1)%helium%energy_avrg
526 55 : iteration = logger%iter_info%iteration(2)
527 :
528 : unit_nr = cp_print_key_unit_nr( &
529 : logger, &
530 : print_key, &
531 : middle_name="helium-energy", &
532 : extension=".dat", &
533 55 : is_new_file=file_is_new)
534 :
535 55 : IF (file_is_new) THEN
536 : WRITE (unit_nr, '(A9,7(1X,A20))') &
537 11 : "# Step", &
538 11 : " AccRatio", &
539 11 : " E_pot", &
540 11 : " E_kin", &
541 11 : " E_thermo", &
542 11 : " E_virial", &
543 11 : " E_inter", &
544 22 : " E_tot"
545 : END IF
546 :
547 : WRITE (unit_nr, "(I9,7(1X,F20.9))") &
548 55 : iteration, &
549 55 : naccptd/ntot, &
550 55 : my_energy(e_id_potential), &
551 55 : my_energy(e_id_kinetic), &
552 55 : my_energy(e_id_thermo), &
553 55 : my_energy(e_id_virial), &
554 55 : my_energy(e_id_interact), &
555 110 : my_energy(e_id_total)
556 :
557 55 : CALL m_flush(unit_nr)
558 55 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
559 :
560 : END IF
561 : END IF
562 :
563 105 : CALL timestop(handle)
564 :
565 105 : RETURN
566 105 : END SUBROUTINE helium_print_energy
567 :
568 : ! ***************************************************************************
569 : !> \brief Print a 3D real vector according to printkey <pkey>
570 : !> \param helium_env ...
571 : !> \param pkey ...
572 : !> \param DATA ...
573 : !> \param uconv ...
574 : !> \param col_label ...
575 : !> \param cmmnt ...
576 : !> \param fname ...
577 : !> \param fpos ...
578 : !> \param avg ...
579 : !> \date 2014-09-09
580 : !> \par History
581 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
582 : !> \author Lukasz Walewski
583 : ! **************************************************************************************************
584 630 : SUBROUTINE helium_print_vector(helium_env, pkey, DATA, uconv, col_label, cmmnt, fname, fpos, avg)
585 :
586 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
587 : CHARACTER(len=*) :: pkey
588 : REAL(KIND=dp), DIMENSION(:), POINTER :: DATA
589 : REAL(KIND=dp) :: uconv
590 : CHARACTER(len=*), DIMENSION(3) :: col_label
591 : CHARACTER(len=*) :: cmmnt, fname
592 : CHARACTER(len=*), OPTIONAL :: fpos
593 : LOGICAL, OPTIONAL :: avg
594 :
595 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_vector'
596 :
597 : CHARACTER(len=default_string_length) :: my_fpos
598 : INTEGER :: handle, i, irank, msglen, nenv, offset, &
599 : unit_nr
600 : LOGICAL :: is_new, my_avg, should_output
601 : REAL(KIND=dp), DIMENSION(3) :: avg_data
602 630 : REAL(KIND=dp), DIMENSION(:), POINTER :: data_p
603 : TYPE(cp_logger_type), POINTER :: logger
604 : TYPE(section_vals_type), POINTER :: psection
605 :
606 630 : CALL timeset(routineN, handle)
607 :
608 630 : IF (PRESENT(avg)) THEN
609 315 : my_avg = avg
610 : ELSE
611 : my_avg = .FALSE.
612 : END IF
613 :
614 630 : IF (PRESENT(fpos)) THEN
615 315 : my_fpos = fpos
616 : ELSE
617 315 : my_fpos = "APPEND"
618 : END IF
619 :
620 630 : NULLIFY (logger, psection)
621 630 : logger => cp_get_default_logger()
622 :
623 630 : psection => section_vals_get_subs_vals(helium_env(1)%helium%input, pkey)
624 : should_output = BTEST(cp_print_key_should_output( &
625 : iteration_info=logger%iter_info, &
626 630 : basis_section=psection), cp_p_file)
627 :
628 630 : IF (.NOT. should_output) THEN
629 315 : CALL timestop(handle)
630 315 : RETURN
631 : END IF
632 :
633 315 : IF (my_avg) THEN
634 : ! average data over all processors and environments
635 315 : avg_data(:) = 0.0_dp
636 315 : msglen = SIZE(avg_data)
637 720 : DO i = 0, SIZE(helium_env) - 1
638 1935 : avg_data(:) = avg_data(:) + DATA(msglen*i + 1:msglen*(i + 1))
639 : END DO
640 315 : CALL helium_env(1)%comm%sum(avg_data)
641 315 : nenv = helium_env(1)%helium%num_env
642 1260 : avg_data(:) = avg_data(:)/REAL(nenv, dp)
643 : ELSE
644 : ! gather data from all processors
645 0 : offset = 0
646 0 : DO i = 1, logger%para_env%mepos
647 0 : offset = offset + helium_env(1)%env_all(i)
648 : END DO
649 :
650 0 : helium_env(1)%helium%rtmp_3_np_1d = 0.0_dp
651 0 : msglen = SIZE(avg_data)
652 0 : DO i = 0, SIZE(helium_env) - 1
653 0 : helium_env(1)%helium%rtmp_3_np_1d(msglen*(offset + i) + 1:msglen*(offset + i + 1)) = DATA(msglen*i + 1:msglen*(i + 1))
654 : END DO
655 0 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rtmp_3_np_1d)
656 : END IF
657 :
658 : unit_nr = cp_print_key_unit_nr( &
659 : logger, &
660 : psection, &
661 : middle_name=fname, &
662 : extension=".dat", &
663 : file_position=my_fpos, &
664 315 : is_new_file=is_new)
665 :
666 315 : IF (logger%para_env%is_source()) THEN
667 :
668 165 : IF (is_new) THEN
669 : ! comment
670 165 : IF (cmmnt .NE. "") THEN
671 165 : WRITE (unit_nr, '(A)') "# "//cmmnt
672 : END IF
673 : ! column labels
674 : WRITE (unit_nr, '(A2,A18,1X,A20,1X,A20)') &
675 165 : "# ", &
676 165 : col_label(1), &
677 165 : col_label(2), &
678 330 : col_label(3)
679 : END IF
680 :
681 165 : IF (my_avg) THEN
682 660 : DO i = 1, 3
683 495 : WRITE (unit_nr, '(E27.20)', ADVANCE='NO') uconv*avg_data(i)
684 660 : IF (i .LT. 3) THEN
685 330 : WRITE (unit_nr, '(1X)', ADVANCE='NO')
686 : END IF
687 : END DO
688 165 : WRITE (unit_nr, '(A)') ""
689 : ELSE
690 : ! iterate over processors/helium environments
691 0 : DO irank = 1, helium_env(1)%helium%num_env
692 : ! unpack data (actually point to the right fragment only)
693 0 : msglen = SIZE(avg_data)
694 0 : offset = (irank - 1)*msglen
695 0 : data_p => helium_env(1)%helium%rtmp_3_np_1d(offset + 1:offset + msglen)
696 : ! write out the data
697 0 : DO i = 1, 3
698 0 : WRITE (unit_nr, '(E27.20)', ADVANCE='NO') uconv*data_p(i)
699 0 : IF (i .LT. 3) THEN
700 0 : WRITE (unit_nr, '(1X)', ADVANCE='NO')
701 : END IF
702 : END DO
703 0 : WRITE (unit_nr, '(A)') ""
704 : END DO
705 : END IF
706 :
707 165 : CALL m_flush(unit_nr)
708 165 : CALL cp_print_key_finished_output(unit_nr, logger, psection)
709 :
710 : END IF
711 :
712 315 : CALL timestop(handle)
713 :
714 315 : RETURN
715 630 : END SUBROUTINE helium_print_vector
716 :
717 : ! ***************************************************************************
718 : !> \brief Print acceptance counts according to HELIUM%PRINT%ACCEPTS
719 : !> \param helium_env ...
720 : !> \date 2010-05-27
721 : !> \par History
722 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
723 : !> \author Lukasz Walewski
724 : ! **************************************************************************************************
725 105 : SUBROUTINE helium_print_accepts(helium_env)
726 :
727 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
728 :
729 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_accepts'
730 :
731 : INTEGER :: handle, i, j, unit_nr
732 : LOGICAL :: file_is_new, should_output
733 : TYPE(cp_logger_type), POINTER :: logger
734 : TYPE(section_vals_type), POINTER :: print_key
735 :
736 105 : CALL timeset(routineN, handle)
737 :
738 105 : NULLIFY (logger, print_key)
739 105 : logger => cp_get_default_logger()
740 :
741 105 : IF (logger%para_env%is_source()) THEN
742 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
743 55 : "MOTION%PINT%HELIUM%PRINT%ACCEPTS")
744 : should_output = BTEST(cp_print_key_should_output( &
745 : iteration_info=logger%iter_info, &
746 55 : basis_section=print_key), cp_p_file)
747 :
748 55 : IF (should_output) THEN
749 : unit_nr = cp_print_key_unit_nr( &
750 : logger, &
751 : print_key, &
752 : middle_name="helium-accepts", &
753 : extension=".dat", &
754 0 : is_new_file=file_is_new)
755 :
756 0 : IF (file_is_new) THEN
757 : WRITE (unit_nr, '(A8,1X,A15,1X,A20)', ADVANCE='NO') &
758 0 : "# Length", &
759 0 : " Trials", &
760 0 : " Selected"
761 0 : DO j = 1, helium_env(1)%helium%bisctlog2
762 0 : WRITE (unit_nr, '(A17,1X,I3)', ADVANCE='NO') " Level", j
763 : END DO
764 0 : WRITE (unit_nr, '(A)') ""
765 : END IF
766 :
767 0 : DO i = 1, helium_env(1)%helium%maxcycle
768 0 : WRITE (unit_nr, '(I3)', ADVANCE='NO') i
769 0 : DO j = 1, helium_env(1)%helium%bisctlog2 + 2
770 0 : WRITE (unit_nr, '(1X,F20.2)', ADVANCE='NO') helium_env(1)%helium%num_accepted(j, i)
771 : END DO
772 0 : WRITE (unit_nr, '(A)') ""
773 : END DO
774 0 : WRITE (unit_nr, '(A)') "&"
775 :
776 0 : CALL m_flush(unit_nr)
777 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
778 :
779 : END IF
780 : END IF
781 :
782 105 : CALL timestop(handle)
783 105 : RETURN
784 : END SUBROUTINE helium_print_accepts
785 :
786 : ! ***************************************************************************
787 : !> \brief Print permutation state according to HELIUM%PRINT%PERM
788 : !> \param helium_env ...
789 : !> \date 2010-06-07
790 : !> \par History
791 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
792 : !> \author Lukasz Walewski
793 : ! **************************************************************************************************
794 105 : SUBROUTINE helium_print_perm(helium_env)
795 :
796 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
797 :
798 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_perm'
799 :
800 : CHARACTER :: left_delim, right_delim
801 : CHARACTER(len=default_string_length) :: msg_str, my_middle_name, stmp
802 : INTEGER :: curat, handle, i, irank, j, lb, msglen, &
803 : nused, offset, outformat, ub, unit_nr
804 105 : INTEGER, DIMENSION(:), POINTER :: my_cycle, my_perm, used_indices
805 : LOGICAL :: first, first_cycle, should_output
806 : TYPE(cp_logger_type), POINTER :: logger
807 : TYPE(section_vals_type), POINTER :: print_key
808 :
809 105 : CALL timeset(routineN, handle)
810 :
811 105 : NULLIFY (logger, print_key)
812 105 : NULLIFY (used_indices)
813 105 : logger => cp_get_default_logger()
814 :
815 : ! determine offset for arrays
816 105 : offset = 0
817 155 : DO i = 1, logger%para_env%mepos
818 155 : offset = offset + helium_env(1)%env_all(i)
819 : END DO
820 :
821 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
822 105 : "MOTION%PINT%HELIUM%PRINT%PERM")
823 : should_output = BTEST(cp_print_key_should_output( &
824 : iteration_info=logger%iter_info, &
825 105 : basis_section=print_key), cp_p_file)
826 :
827 105 : IF (.NOT. should_output) THEN
828 105 : CALL timestop(handle)
829 105 : RETURN
830 : END IF
831 :
832 : ! get the output type
833 0 : CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
834 0 : IF (outformat .EQ. perm_cycle) THEN
835 : ! gather positions from all processors
836 0 : helium_env(1)%helium%rtmp_3_atoms_beads_np_1d = 0.0_dp
837 0 : j = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
838 0 : DO i = 1, SIZE(helium_env)
839 : helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(j*(offset + i - 1) + 1:j*(offset + i)) = &
840 0 : PACK(helium_env(i)%helium%pos(:, :, 1:helium_env(1)%helium%beads), .TRUE.)
841 : END DO
842 0 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d)
843 : ! set logical mask for unpacking coordinates gathered from other ranks
844 0 : helium_env(1)%helium%ltmp_3_atoms_beads_3d(:, :, :) = .TRUE.
845 : END IF
846 :
847 : ! gather permutation state from all processors to logger%para_env%source
848 0 : helium_env(1)%helium%itmp_atoms_np_1d(:) = 0
849 0 : msglen = SIZE(helium_env(1)%helium%permutation)
850 0 : DO i = 1, SIZE(helium_env)
851 0 : helium_env(1)%helium%itmp_atoms_np_1d(msglen*(offset + i - 1) + 1:msglen*(offset + i)) = helium_env(i)%helium%permutation
852 : END DO
853 :
854 0 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%itmp_atoms_np_1d)
855 :
856 0 : IF (logger%para_env%is_source()) THEN
857 :
858 : ! iterate over helium environments
859 0 : DO irank = 1, helium_env(1)%helium%num_env
860 :
861 : ! generate one file per environment
862 0 : stmp = ""
863 0 : WRITE (stmp, *) irank
864 0 : my_middle_name = "helium-perm-"//TRIM(ADJUSTL(stmp))
865 : unit_nr = cp_print_key_unit_nr( &
866 : logger, &
867 : print_key, &
868 : middle_name=TRIM(my_middle_name), &
869 0 : extension=".dat")
870 :
871 : ! unpack permutation state (actually point to the right section only)
872 0 : lb = (irank - 1)*helium_env(1)%helium%atoms + 1
873 0 : ub = irank*helium_env(1)%helium%atoms
874 0 : my_perm => helium_env(1)%helium%itmp_atoms_np_1d(lb:ub)
875 :
876 0 : SELECT CASE (outformat)
877 :
878 : CASE (perm_cycle)
879 : ! write the permutation grouped into cycles
880 :
881 : ! unpack coordinates (necessary only for winding path delimiters)
882 0 : msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
883 0 : offset = (irank - 1)*msglen
884 : helium_env(1)%helium%work(:, :, 1:helium_env(1)%helium%beads) = &
885 : UNPACK(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(offset + 1:offset + msglen), &
886 0 : MASK=helium_env(1)%helium%ltmp_3_atoms_beads_3d, FIELD=0.0_dp)
887 :
888 0 : curat = 1
889 0 : nused = 0
890 0 : first_cycle = .TRUE.
891 0 : DO WHILE (curat .LE. helium_env(1)%helium%atoms)
892 :
893 : ! get the permutation cycle the current atom belongs to
894 0 : my_cycle => helium_cycle_of(curat, my_perm)
895 :
896 : ! include the current cycle in the pool of "used" indices
897 0 : nused = nused + SIZE(my_cycle)
898 0 : CALL reallocate(used_indices, 1, nused)
899 0 : used_indices(nused - SIZE(my_cycle) + 1:nused) = my_cycle(1:SIZE(my_cycle))
900 :
901 : ! select delimiters according to the cycle's winding state
902 0 : IF (helium_is_winding(helium_env(1)%helium, curat, helium_env(1)%helium%work, my_perm)) THEN
903 0 : left_delim = "["
904 0 : right_delim = "]"
905 : ELSE
906 0 : left_delim = "("
907 0 : right_delim = ")"
908 : END IF
909 :
910 : ! cycle delimiter
911 0 : IF (first_cycle) THEN
912 : first_cycle = .FALSE.
913 : ELSE
914 0 : WRITE (unit_nr, '(1X)', ADVANCE='NO')
915 : END IF
916 :
917 : ! write out the current cycle
918 0 : WRITE (unit_nr, '(A1)', ADVANCE='NO') left_delim
919 0 : first = .TRUE.
920 0 : DO i = 1, SIZE(my_cycle)
921 0 : IF (first) THEN
922 : first = .FALSE.
923 : ELSE
924 0 : WRITE (unit_nr, '(1X)', ADVANCE='NO')
925 : END IF
926 0 : WRITE (unit_nr, '(I0)', ADVANCE='NO') my_cycle(i)
927 : END DO
928 0 : WRITE (unit_nr, '(A1)', ADVANCE='NO') right_delim
929 :
930 : ! clean up
931 0 : DEALLOCATE (my_cycle)
932 0 : NULLIFY (my_cycle)
933 :
934 : ! try to increment the current atom index
935 0 : DO WHILE (ANY(used_indices .EQ. curat))
936 0 : curat = curat + 1
937 : END DO
938 :
939 : END DO
940 0 : WRITE (unit_nr, *)
941 :
942 0 : DEALLOCATE (used_indices)
943 0 : NULLIFY (used_indices)
944 :
945 : CASE (perm_plain)
946 : ! write the plain permutation
947 :
948 0 : first = .TRUE.
949 0 : DO i = 1, helium_env(1)%helium%atoms
950 0 : IF (first) THEN
951 : first = .FALSE.
952 : ELSE
953 0 : WRITE (unit_nr, '(1X)', ADVANCE='NO')
954 : END IF
955 0 : WRITE (unit_nr, '(I0)', ADVANCE='NO') my_perm(i)
956 : END DO
957 0 : WRITE (unit_nr, '(A)') ""
958 :
959 : CASE default
960 :
961 0 : msg_str = "Format for permutation output unknown."
962 0 : CPABORT(msg_str)
963 :
964 : END SELECT
965 :
966 0 : CALL m_flush(unit_nr)
967 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
968 :
969 : END DO
970 :
971 : END IF
972 :
973 0 : CALL timestop(handle)
974 :
975 0 : RETURN
976 105 : END SUBROUTINE helium_print_perm
977 :
978 : ! **************************************************************************************************
979 : !> \brief Print helium action file to HELIUM%PRINT%ACTION
980 : !> \param pint_env ...
981 : !> \param helium_env ...
982 : !> \date 2016-06-07
983 : !> \par History
984 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
985 : !> \author Felix Uhl
986 : ! **************************************************************************************************
987 105 : SUBROUTINE helium_print_action(pint_env, helium_env)
988 :
989 : TYPE(pint_env_type), INTENT(IN) :: pint_env
990 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
991 :
992 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_action'
993 :
994 : CHARACTER(len=default_string_length) :: my_middle_name, stmp
995 : INTEGER :: handle, i, irank, iteration, k, offset, &
996 : unit_nr
997 : LOGICAL :: file_is_new, should_output
998 105 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_inter_action, tmp_link_action, &
999 105 : tmp_pair_action
1000 : TYPE(cp_logger_type), POINTER :: logger
1001 : TYPE(section_vals_type), POINTER :: print_key
1002 :
1003 105 : CALL timeset(routineN, handle)
1004 :
1005 105 : NULLIFY (logger, print_key)
1006 105 : logger => cp_get_default_logger()
1007 :
1008 : ! determine offset for arrays
1009 105 : offset = 0
1010 155 : DO i = 1, logger%para_env%mepos
1011 155 : offset = offset + helium_env(1)%env_all(i)
1012 : END DO
1013 :
1014 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1015 105 : "MOTION%PINT%HELIUM%PRINT%ACTION")
1016 : should_output = BTEST(cp_print_key_should_output( &
1017 : iteration_info=logger%iter_info, &
1018 105 : basis_section=print_key), cp_p_file)
1019 :
1020 105 : IF (.NOT. should_output) THEN
1021 85 : CALL timestop(handle)
1022 85 : RETURN
1023 : END IF
1024 :
1025 70 : DO k = 1, SIZE(helium_env)
1026 50 : helium_env(k)%helium%link_action = helium_total_link_action(helium_env(k)%helium)
1027 50 : helium_env(k)%helium%pair_action = helium_total_pair_action(helium_env(k)%helium)
1028 70 : helium_env(k)%helium%inter_action = helium_total_inter_action(pint_env, helium_env(k)%helium)
1029 : END DO
1030 :
1031 20 : NULLIFY (tmp_link_action)
1032 20 : NULLIFY (tmp_pair_action)
1033 20 : NULLIFY (tmp_inter_action)
1034 60 : ALLOCATE (tmp_link_action(helium_env(1)%helium%num_env))
1035 40 : ALLOCATE (tmp_pair_action(helium_env(1)%helium%num_env))
1036 40 : ALLOCATE (tmp_inter_action(helium_env(1)%helium%num_env))
1037 120 : tmp_link_action(:) = 0.0_dp
1038 120 : tmp_pair_action(:) = 0.0_dp
1039 120 : tmp_inter_action(:) = 0.0_dp
1040 : ! gather Action from all processors to logger%para_env%source
1041 70 : DO k = 1, SIZE(helium_env)
1042 50 : tmp_link_action(offset + k) = helium_env(k)%helium%link_action
1043 50 : tmp_pair_action(offset + k) = helium_env(k)%helium%pair_action
1044 70 : tmp_inter_action(offset + k) = helium_env(k)%helium%inter_action
1045 : END DO
1046 220 : CALL helium_env(1)%comm%sum(tmp_link_action)
1047 220 : CALL helium_env(1)%comm%sum(tmp_pair_action)
1048 220 : CALL helium_env(1)%comm%sum(tmp_inter_action)
1049 :
1050 20 : IF (logger%para_env%is_source()) THEN
1051 10 : iteration = logger%iter_info%iteration(2)
1052 : ! iterate over processors/helium environments
1053 60 : DO irank = 1, helium_env(1)%helium%num_env
1054 :
1055 : ! generate one file per helium_env
1056 50 : stmp = ""
1057 50 : WRITE (stmp, *) irank
1058 50 : my_middle_name = "helium-action-"//TRIM(ADJUSTL(stmp))
1059 : unit_nr = cp_print_key_unit_nr( &
1060 : logger, &
1061 : print_key, &
1062 : middle_name=TRIM(my_middle_name), &
1063 : extension=".dat", &
1064 50 : is_new_file=file_is_new)
1065 :
1066 50 : IF (file_is_new) THEN
1067 : WRITE (unit_nr, '(A9,3(1X,A25))') &
1068 5 : "# Step", &
1069 5 : " He_Total_Link_Action", &
1070 5 : " He_Total_Pair_Action", &
1071 10 : " He_Total_Interaction"
1072 : END IF
1073 :
1074 : WRITE (unit_nr, "(I9,3(1X,F25.14))") &
1075 50 : iteration, &
1076 50 : tmp_link_action(irank), &
1077 50 : tmp_pair_action(irank), &
1078 100 : tmp_inter_action(irank)
1079 :
1080 50 : CALL m_flush(unit_nr)
1081 60 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1082 :
1083 : END DO
1084 : END IF
1085 :
1086 20 : DEALLOCATE (tmp_link_action)
1087 20 : DEALLOCATE (tmp_pair_action)
1088 20 : DEALLOCATE (tmp_inter_action)
1089 20 : NULLIFY (tmp_link_action)
1090 20 : NULLIFY (tmp_pair_action)
1091 20 : NULLIFY (tmp_inter_action)
1092 :
1093 20 : CALL timestop(handle)
1094 :
1095 20 : RETURN
1096 105 : END SUBROUTINE helium_print_action
1097 :
1098 : ! ***************************************************************************
1099 : !> \brief Print coordinates according to HELIUM%PRINT%COORDINATES
1100 : !> \param helium_env ...
1101 : !> \date 2009-07-16
1102 : !> \par History
1103 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
1104 : !> \author Lukasz Walewski
1105 : ! **************************************************************************************************
1106 105 : SUBROUTINE helium_print_coordinates(helium_env)
1107 :
1108 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1109 :
1110 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_coordinates'
1111 :
1112 : CHARACTER(3) :: resName
1113 : CHARACTER(len=default_string_length) :: fmt_string, my_middle_name, stmp
1114 : INTEGER :: handle, i, ia, ib, ib1, ib2, ic, icycle, &
1115 : irank, j, k, msglen, offset, &
1116 : outformat, tmp1, tmp2, unit_nr
1117 105 : INTEGER, DIMENSION(:), POINTER :: my_perm
1118 : LOGICAL :: are_connected, is_winding, ltmp, &
1119 : should_output
1120 : REAL(KIND=dp) :: xtmp, ytmp, ztmp
1121 : REAL(KIND=dp), DIMENSION(3) :: r0, r1, r2
1122 : TYPE(cp_logger_type), POINTER :: logger
1123 : TYPE(section_vals_type), POINTER :: print_key
1124 :
1125 105 : CALL timeset(routineN, handle)
1126 :
1127 105 : NULLIFY (logger, print_key)
1128 105 : logger => cp_get_default_logger()
1129 :
1130 : ! determine offset for arrays
1131 105 : offset = 0
1132 155 : DO i = 1, logger%para_env%mepos
1133 155 : offset = offset + helium_env(1)%env_all(i)
1134 : END DO
1135 :
1136 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1137 105 : "MOTION%PINT%HELIUM%PRINT%COORDINATES")
1138 : should_output = BTEST(cp_print_key_should_output( &
1139 : iteration_info=logger%iter_info, &
1140 105 : basis_section=print_key), cp_p_file)
1141 :
1142 105 : IF (.NOT. should_output) THEN
1143 85 : CALL timestop(handle)
1144 85 : RETURN
1145 : END IF
1146 :
1147 : ! prepare the coordinates for output (use unit cell centered around r0)
1148 70 : DO k = 1, SIZE(helium_env)
1149 200 : r0(:) = helium_env(k)%helium%center(:)
1150 320 : DO ia = 1, helium_env(k)%helium%atoms
1151 4300 : DO ib = 1, helium_env(k)%helium%beads
1152 16000 : r1(:) = helium_env(k)%helium%pos(:, ia, ib) - r0(:)
1153 16000 : r2(:) = helium_env(k)%helium%pos(:, ia, ib) - r0(:)
1154 4000 : CALL helium_pbc(helium_env(k)%helium, r2)
1155 4000 : ltmp = .FALSE.
1156 16000 : DO ic = 1, 3
1157 16000 : IF (ABS(r1(ic) - r2(ic)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
1158 0 : ltmp = .TRUE.
1159 0 : CYCLE
1160 : END IF
1161 : END DO
1162 4250 : IF (ltmp) THEN
1163 0 : helium_env(k)%helium%work(:, ia, ib) = r0(:) + r2(:)
1164 : ELSE
1165 16000 : helium_env(k)%helium%work(:, ia, ib) = helium_env(k)%helium%pos(:, ia, ib)
1166 : END IF
1167 : END DO
1168 : END DO
1169 : END DO
1170 :
1171 : ! gather positions from all processors to logger%para_env%source
1172 24020 : helium_env(1)%helium%rtmp_3_atoms_beads_np_1d = 0.0_dp
1173 20 : j = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
1174 70 : DO i = 1, SIZE(helium_env)
1175 : helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(j*(offset + i - 1) + 1:j*(offset + i)) = &
1176 12070 : PACK(helium_env(i)%helium%pos(:, :, 1:helium_env(1)%helium%beads), .TRUE.)
1177 : END DO
1178 48020 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d)
1179 :
1180 : ! gather permutation state from all processors to logger%para_env%source
1181 520 : helium_env(1)%helium%itmp_atoms_np_1d(:) = 0
1182 20 : j = SIZE(helium_env(1)%helium%permutation)
1183 70 : DO i = 1, SIZE(helium_env)
1184 320 : helium_env(1)%helium%itmp_atoms_np_1d(j*(offset + i - 1) + 1:j*(offset + i)) = helium_env(i)%helium%permutation
1185 : END DO
1186 :
1187 1020 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%itmp_atoms_np_1d)
1188 :
1189 : ! set logical mask for unpacking coordinates gathered from other ranks
1190 6740 : helium_env(1)%helium%ltmp_3_atoms_beads_3d(:, :, :) = .TRUE.
1191 :
1192 20 : IF (logger%para_env%is_source()) THEN
1193 :
1194 10 : CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
1195 :
1196 : ! iterate over helium environments
1197 60 : DO irank = 1, helium_env(1)%helium%num_env
1198 60 : IF (outformat == fmt_id_pdb) THEN
1199 : ! generate one file per environment
1200 50 : stmp = ""
1201 50 : WRITE (stmp, *) irank
1202 50 : my_middle_name = "helium-pos-"//TRIM(ADJUSTL(stmp))
1203 : unit_nr = cp_print_key_unit_nr( &
1204 : logger, &
1205 : print_key, &
1206 : middle_name=TRIM(my_middle_name), &
1207 50 : extension=".pdb")
1208 :
1209 : ! write out the unit cell parameters
1210 50 : fmt_string = "(A6,3F9.3,3F7.2,1X,A11,1X,I3)"
1211 50 : xtmp = helium_env(1)%helium%cell_size
1212 50 : xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
1213 100 : SELECT CASE (helium_env(1)%helium%cell_shape)
1214 : CASE (helium_cell_shape_cube)
1215 50 : stmp = "C "
1216 : CASE (helium_cell_shape_octahedron)
1217 0 : stmp = "O "
1218 : CASE DEFAULT
1219 50 : stmp = "UNKNOWN "
1220 : END SELECT
1221 50 : WRITE (unit_nr, fmt_string) "CRYST1", &
1222 50 : xtmp, xtmp, xtmp, &
1223 50 : 90.0_dp, 90.0_dp, 90.0_dp, &
1224 100 : stmp, helium_env(1)%helium%beads
1225 :
1226 : ! unpack coordinates
1227 50 : msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
1228 50 : offset = (irank - 1)*msglen
1229 : helium_env(1)%helium%work(:, :, 1:helium_env(1)%helium%beads) = &
1230 : UNPACK(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(offset + 1:offset + msglen), &
1231 16850 : MASK=helium_env(1)%helium%ltmp_3_atoms_beads_3d, FIELD=0.0_dp)
1232 :
1233 : ! unpack permutation state (actually point to the right section only)
1234 50 : msglen = SIZE(helium_env(1)%helium%permutation)
1235 50 : offset = (irank - 1)*msglen
1236 50 : my_perm => helium_env(1)%helium%itmp_atoms_np_1d(offset + 1:offset + msglen)
1237 :
1238 : ! write out coordinates
1239 : fmt_string = &
1240 50 : "(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2)"
1241 300 : DO ia = 1, helium_env(1)%helium%atoms
1242 250 : icycle = helium_cycle_number(helium_env(1)%helium, ia, my_perm)
1243 250 : is_winding = helium_is_winding(helium_env(1)%helium, ia, helium_env(1)%helium%work, my_perm)
1244 250 : IF (is_winding) THEN
1245 0 : resName = "SPR"
1246 : ELSE
1247 250 : resName = "NRM"
1248 : END IF
1249 4300 : DO ib = 1, helium_env(1)%helium%beads
1250 4000 : xtmp = helium_env(1)%helium%work(1, ia, ib)
1251 4000 : xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
1252 4000 : ytmp = helium_env(1)%helium%work(2, ia, ib)
1253 4000 : ytmp = cp_unit_from_cp2k(ytmp, "angstrom")
1254 4000 : ztmp = helium_env(1)%helium%work(3, ia, ib)
1255 4000 : ztmp = cp_unit_from_cp2k(ztmp, "angstrom")
1256 4000 : WRITE (unit_nr, fmt_string) "ATOM ", &
1257 4000 : (ia - 1)*helium_env(1)%helium%beads + ib, &
1258 4000 : " He ", " ", resName, "X", &
1259 4000 : icycle, &
1260 4000 : " ", &
1261 4000 : xtmp, ytmp, ztmp, &
1262 8250 : 1.0_dp, 0.0_dp, "HE", " "
1263 : END DO
1264 : END DO
1265 :
1266 : ! write out the bead connectivity information
1267 300 : DO ia = 1, helium_env(1)%helium%atoms
1268 :
1269 : ! write connectivity records for this atom only if the path
1270 : ! it belongs to is longer than 1.
1271 250 : IF (helium_path_length(helium_env(1)%helium, ia, my_perm) .LE. 1) THEN
1272 : CYCLE
1273 : END IF
1274 :
1275 0 : DO ib = 1, helium_env(1)%helium%beads - 1
1276 : ! check whether the consecutive beads belong to the same box
1277 0 : r1(:) = helium_env(1)%helium%work(:, ia, ib) - helium_env(1)%helium%work(:, ia, ib + 1)
1278 0 : r2(:) = r1(:)
1279 0 : CALL helium_pbc(helium_env(1)%helium, r2)
1280 0 : are_connected = .TRUE.
1281 0 : DO ic = 1, 3
1282 0 : IF (ABS(r1(ic) - r2(ic)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
1283 : ! if the distance betw ib and ib+1 changes upon applying
1284 : ! PBC do not connect them
1285 0 : are_connected = .FALSE.
1286 0 : CYCLE
1287 : END IF
1288 : END DO
1289 0 : IF (are_connected) THEN
1290 0 : tmp1 = (ia - 1)*helium_env(1)%helium%beads + ib
1291 0 : tmp2 = (ia - 1)*helium_env(1)%helium%beads + ib + 1
1292 : ! smaller value has to go first
1293 : IF (tmp1 .LT. tmp2) THEN
1294 0 : ib1 = tmp1
1295 0 : ib2 = tmp2
1296 : ELSE
1297 : ib1 = tmp2
1298 : ib2 = tmp1
1299 : END IF
1300 0 : WRITE (unit_nr, '(A6,2I5)') "CONECT", ib1, ib2
1301 : END IF
1302 : END DO
1303 :
1304 : ! last bead of atom <ia> connects to the first bead
1305 : ! of the next atom in the permutation cycle
1306 0 : r1(:) = helium_env(1)%helium%work(:, ia, helium_env(1)%helium%beads) - helium_env(1)%helium%work(:, my_perm(ia), 1)
1307 0 : r2(:) = r1(:)
1308 0 : CALL helium_pbc(helium_env(1)%helium, r2)
1309 0 : are_connected = .TRUE.
1310 0 : DO ic = 1, 3
1311 0 : IF (ABS(r1(ic) - r2(ic)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
1312 : ! if the distance betw ib and ib+1 changes upon applying
1313 : ! PBC do not connect them
1314 0 : are_connected = .FALSE.
1315 0 : CYCLE
1316 : END IF
1317 : END DO
1318 50 : IF (are_connected) THEN
1319 0 : tmp1 = ia*helium_env(1)%helium%beads
1320 0 : tmp2 = (my_perm(ia) - 1)*helium_env(1)%helium%beads + 1
1321 0 : IF (tmp1 .LT. tmp2) THEN
1322 0 : ib1 = tmp1
1323 0 : ib2 = tmp2
1324 : ELSE
1325 0 : ib1 = tmp2
1326 0 : ib2 = tmp1
1327 : END IF
1328 0 : WRITE (unit_nr, '(A6,2I5)') "CONECT", ib1, ib2
1329 : END IF
1330 : END DO
1331 50 : WRITE (unit_nr, '(A)') "END"
1332 :
1333 50 : CALL m_flush(unit_nr)
1334 50 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1335 0 : ELSE IF (outformat == fmt_id_xyz) THEN
1336 : ! generate one file per environment and bead
1337 0 : DO ib = 1, helium_env(1)%helium%beads
1338 0 : stmp = ""
1339 0 : WRITE (stmp, *) irank
1340 0 : my_middle_name = "helium-pos-"//TRIM(ADJUSTL(stmp))
1341 0 : WRITE (stmp, *) ib
1342 0 : my_middle_name = TRIM(my_middle_name)//"-"//TRIM(ADJUSTL(stmp))
1343 : unit_nr = cp_print_key_unit_nr( &
1344 : logger, &
1345 : print_key, &
1346 : middle_name=TRIM(my_middle_name), &
1347 0 : extension=".xyz")
1348 : ! write out xyz header
1349 0 : WRITE (unit_nr, *) helium_env(1)%helium%atoms
1350 0 : stmp = ""
1351 0 : WRITE (stmp, *) logger%iter_info%n_rlevel
1352 0 : fmt_string = "(A6,"//TRIM(ADJUSTL(stmp))//"I12)"
1353 0 : WRITE (unit_nr, fmt_string) "iter= ", logger%iter_info%iteration(:)
1354 0 : fmt_string = "(A6,3F9.3,3F7.2,1X,A11,1X,I3)"
1355 :
1356 : ! unpack coordinates
1357 0 : msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
1358 0 : offset = (irank - 1)*msglen
1359 : helium_env(1)%helium%work(:, :, 1:helium_env(1)%helium%beads) = &
1360 : UNPACK(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(offset + 1:offset + msglen), &
1361 0 : MASK=helium_env(1)%helium%ltmp_3_atoms_beads_3d, FIELD=0.0_dp)
1362 :
1363 : ! unpack permutation state (actually point to the right section only)
1364 0 : msglen = SIZE(helium_env(1)%helium%permutation)
1365 0 : offset = (irank - 1)*msglen
1366 0 : my_perm => helium_env(1)%helium%itmp_atoms_np_1d(offset + 1:offset + msglen)
1367 :
1368 : ! write out coordinates
1369 0 : fmt_string = "(A2,3(1X,F20.10))"
1370 0 : DO ia = 1, helium_env(1)%helium%atoms
1371 0 : xtmp = helium_env(1)%helium%work(1, ia, ib)
1372 0 : xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
1373 0 : ytmp = helium_env(1)%helium%work(2, ia, ib)
1374 0 : ytmp = cp_unit_from_cp2k(ytmp, "angstrom")
1375 0 : ztmp = helium_env(1)%helium%work(3, ia, ib)
1376 0 : ztmp = cp_unit_from_cp2k(ztmp, "angstrom")
1377 0 : WRITE (unit_nr, fmt_string) "He", xtmp, ytmp, ztmp
1378 : END DO
1379 0 : CALL m_flush(unit_nr)
1380 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1381 : END DO
1382 : ELSE
1383 0 : CPABORT("")
1384 : END IF
1385 : END DO
1386 :
1387 : END IF
1388 :
1389 20 : CALL timestop(handle)
1390 :
1391 20 : RETURN
1392 105 : END SUBROUTINE helium_print_coordinates
1393 :
1394 : ! ***************************************************************************
1395 : !> \brief Print radial distribution functions according to HELIUM%PRINT%RDF
1396 : !> \param helium_env ...
1397 : !> \date 2009-07-23
1398 : !> \par History
1399 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
1400 : !> \author Lukasz Walewski
1401 : ! **************************************************************************************************
1402 10 : SUBROUTINE helium_print_rdf(helium_env)
1403 :
1404 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1405 :
1406 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_rdf'
1407 :
1408 : CHARACTER(len=default_string_length) :: stmp
1409 : INTEGER :: handle, ia, ic, id, itmp, iweight, k, &
1410 : nsteps, unit_nr
1411 : LOGICAL :: should_output
1412 : REAL(KIND=dp) :: inv_norm, rtmp, rtmp2
1413 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: message
1414 : TYPE(cp_logger_type), POINTER :: logger
1415 : TYPE(section_vals_type), POINTER :: print_key
1416 :
1417 10 : CALL timeset(routineN, handle)
1418 :
1419 10 : NULLIFY (logger, print_key)
1420 10 : logger => cp_get_default_logger()
1421 :
1422 10 : IF (logger%para_env%is_source()) THEN
1423 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1424 5 : "MOTION%PINT%HELIUM%PRINT%RDF")
1425 : should_output = BTEST(cp_print_key_should_output( &
1426 : iteration_info=logger%iter_info, &
1427 5 : basis_section=print_key), cp_p_file)
1428 : END IF
1429 10 : CALL helium_env(1)%comm%bcast(should_output, logger%para_env%source)
1430 :
1431 10 : IF (should_output) THEN
1432 : ! work on the temporary array so that accumulated data remains intact
1433 : ! save accumulated data of different env on same core in first temp
1434 5010 : helium_env(1)%helium%rdf_inst(:, :) = 0.0_dp
1435 20 : DO k = 1, SIZE(helium_env)
1436 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :) + &
1437 5020 : helium_env(k)%helium%rdf_accu(:, :)
1438 : END DO
1439 :
1440 : ! average over processors
1441 10 : NULLIFY (message)
1442 10 : message => helium_env(1)%helium%rdf_inst(:, :)
1443 10010 : CALL helium_env(1)%comm%sum(message)
1444 10 : itmp = helium_env(1)%helium%num_env
1445 10 : inv_norm = 1.0_dp/REAL(itmp, dp)
1446 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*inv_norm
1447 :
1448 : ! average over steps performed so far in this run
1449 10 : nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1450 10 : inv_norm = 1.0_dp/REAL(nsteps, dp)
1451 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*inv_norm
1452 :
1453 10 : iweight = helium_env(1)%helium%rdf_iweight
1454 : ! average over the old and the current density (observe the weights!)
1455 : helium_env(1)%helium%rdf_inst(:, :) = nsteps*helium_env(1)%helium%rdf_inst(:, :) + &
1456 5010 : iweight*helium_env(1)%helium%rdf_rstr(:, :)
1457 5010 : helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/REAL(nsteps + iweight, dp)
1458 :
1459 10 : IF (logger%para_env%is_source()) THEN
1460 :
1461 5 : ia = 0
1462 5 : rtmp = cp_unit_from_cp2k(1.0_dp, "angstrom")
1463 5 : rtmp2 = 1.0_dp
1464 5 : IF (.NOT. helium_env(1)%helium%periodic) THEN
1465 : ! RDF in non-periodic case has unit 1/bohr**3, convert to Angstrom:
1466 0 : rtmp2 = rtmp**(-3)
1467 : END IF
1468 :
1469 5 : IF (helium_env(1)%helium%rdf_he_he) THEN
1470 : ! overwrite RDF file each time it is written
1471 5 : ia = 1
1472 5 : stmp = ""
1473 5 : WRITE (stmp, *) "He-He"
1474 : unit_nr = cp_print_key_unit_nr( &
1475 : logger, &
1476 : print_key, &
1477 : middle_name="helium-rdf-"//TRIM(ADJUSTL(stmp)), &
1478 : extension=".dat", &
1479 : file_position="REWIND", &
1480 5 : do_backup=.FALSE.)
1481 :
1482 1255 : DO ic = 1, helium_env(1)%helium%rdf_nbin
1483 1250 : WRITE (unit_nr, '(F20.10)', ADVANCE='NO') (REAL(ic, dp) - 0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
1484 1250 : WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%rdf_inst(ia, ic)*rtmp2
1485 1255 : WRITE (unit_nr, *)
1486 : END DO
1487 :
1488 5 : CALL m_flush(unit_nr)
1489 5 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1490 : END IF
1491 :
1492 5 : IF (helium_env(1)%helium%rdf_sol_he) THEN
1493 : ! overwrite RDF file each time it is written
1494 0 : stmp = ""
1495 0 : WRITE (stmp, *) "Solute-He"
1496 : unit_nr = cp_print_key_unit_nr( &
1497 : logger, &
1498 : print_key, &
1499 : middle_name="helium-rdf-"//TRIM(ADJUSTL(stmp)), &
1500 : extension=".dat", &
1501 : file_position="REWIND", &
1502 0 : do_backup=.FALSE.)
1503 :
1504 0 : DO ic = 1, helium_env(1)%helium%rdf_nbin
1505 0 : WRITE (unit_nr, '(F20.10)', ADVANCE='NO') (REAL(ic, dp) - 0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
1506 0 : DO id = 1 + ia, helium_env(1)%helium%rdf_num
1507 0 : WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%rdf_inst(id, ic)*rtmp2
1508 : END DO
1509 0 : WRITE (unit_nr, *)
1510 : END DO
1511 :
1512 0 : CALL m_flush(unit_nr)
1513 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1514 : END IF
1515 :
1516 : END IF
1517 : END IF
1518 :
1519 10 : CALL timestop(handle)
1520 :
1521 10 : RETURN
1522 10 : END SUBROUTINE helium_print_rdf
1523 :
1524 : ! ***************************************************************************
1525 : !> \brief Print densities according to HELIUM%PRINT%RHO
1526 : !> \param helium_env ...
1527 : !> \date 2011-06-21
1528 : !> \par History
1529 : !> 08.2015 cleaned code from unneeded arrays
1530 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
1531 : !> \author Lukasz Walewski
1532 : ! **************************************************************************************************
1533 10 : SUBROUTINE helium_print_rho(helium_env)
1534 :
1535 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1536 :
1537 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_rho'
1538 :
1539 : CHARACTER(len=default_string_length) :: comment, fname
1540 : INTEGER :: handle, ic, id, itmp, iweight, k, &
1541 : nsteps, unit_nr
1542 : LOGICAL :: should_output
1543 : REAL(KIND=dp) :: inv_norm, invproc
1544 : REAL(KIND=dp), DIMENSION(3) :: center
1545 10 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cubdata
1546 : TYPE(cp_logger_type), POINTER :: logger
1547 : TYPE(section_vals_type), POINTER :: print_key
1548 :
1549 10 : CALL timeset(routineN, handle)
1550 :
1551 10 : NULLIFY (logger, print_key)
1552 10 : logger => cp_get_default_logger()
1553 :
1554 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1555 10 : "MOTION%PINT%HELIUM%PRINT%RHO")
1556 : should_output = BTEST(cp_print_key_should_output( &
1557 : iteration_info=logger%iter_info, &
1558 10 : basis_section=print_key), cp_p_file)
1559 :
1560 10 : IF (.NOT. should_output) THEN
1561 8 : CALL timestop(handle)
1562 8 : RETURN
1563 : END IF
1564 :
1565 : ! work on temporary array so that the average remains intact
1566 : ! save accumulated data of different env on same core in first temp
1567 4222 : helium_env(1)%helium%rho_inst(:, :, :, :) = 0.0_dp
1568 4 : DO k = 1, SIZE(helium_env)
1569 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :) + &
1570 4224 : helium_env(k)%helium%rho_accu(:, :, :, :)
1571 : END DO
1572 :
1573 : ! average over processors
1574 8442 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%rho_inst)
1575 2 : itmp = helium_env(1)%helium%num_env
1576 2 : invproc = 1.0_dp/REAL(itmp, dp)
1577 4222 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)*invproc
1578 :
1579 : ! average over steps performed so far in this run
1580 2 : nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1581 2 : inv_norm = 1.0_dp/REAL(nsteps, dp)
1582 4222 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)*inv_norm
1583 :
1584 2 : iweight = helium_env(1)%helium%rho_iweight
1585 : ! average over the old and the current density (observe the weights!)
1586 : helium_env(1)%helium%rho_inst(:, :, :, :) = nsteps*helium_env(1)%helium%rho_inst(:, :, :, :) + &
1587 4222 : iweight*helium_env(1)%helium%rho_rstr(:, :, :, :)
1588 4222 : helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)/REAL(nsteps + iweight, dp)
1589 :
1590 : ! set center of the cubefile
1591 2 : IF (helium_env(1)%helium%solute_present) THEN
1592 : ! should be set to solute's COM
1593 0 : center(:) = helium_env(1)%helium%center(:)
1594 : ELSE
1595 : ! regardless of whether we are periodic or not use the origin, since
1596 : ! pure cluster's COM might drift, but we want the cube fixed (note that
1597 : ! the densities are correctly calculated wrt to COM in such case)
1598 2 : center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
1599 : END IF
1600 :
1601 12 : DO id = 1, rho_num ! loop over densities ---
1602 :
1603 10 : IF (.NOT. helium_env(1)%helium%rho_property(id)%is_calculated) THEN
1604 : ! skip densities that are not requested by the user
1605 : CYCLE
1606 : END IF
1607 :
1608 6 : DO ic = 1, helium_env(1)%helium%rho_property(id)%num_components ! loop over components
1609 :
1610 : WRITE (fname, '(A)') "helium-rho-"// &
1611 2 : TRIM(ADJUSTL(helium_env(1)%helium%rho_property(id)%filename_suffix(ic)))
1612 2 : IF (helium_env(1)%helium%rho_property(id)%component_name(ic) .EQ. "") THEN
1613 2 : WRITE (comment, '(A)') TRIM(helium_env(1)%helium%rho_property(id)%name)
1614 : ELSE
1615 : WRITE (comment, '(A)') TRIM(helium_env(1)%helium%rho_property(id)%name)// &
1616 : ", "// &
1617 0 : TRIM(helium_env(1)%helium%rho_property(id)%component_name(ic))
1618 : END IF
1619 2 : cubdata => helium_env(1)%helium%rho_inst(helium_env(1)%helium%rho_property(id)%component_index(ic), :, :, :)
1620 :
1621 : unit_nr = cp_print_key_unit_nr( &
1622 : logger, &
1623 : print_key, &
1624 : middle_name=TRIM(ADJUSTL(fname)), &
1625 : extension=".cube", &
1626 : file_position="REWIND", &
1627 2 : do_backup=.FALSE.)
1628 :
1629 12 : IF (logger%para_env%is_source()) THEN
1630 : CALL helium_write_cubefile( &
1631 : unit_nr, &
1632 : comment, &
1633 : center - 0.5_dp*(helium_env(1)%helium%rho_maxr - helium_env(1)%helium%rho_delr), &
1634 : helium_env(1)%helium%rho_delr, &
1635 : helium_env(1)%helium%rho_nbin, &
1636 4 : cubdata)
1637 1 : CALL m_flush(unit_nr)
1638 1 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1639 : END IF
1640 :
1641 : END DO ! loop over components
1642 :
1643 : END DO ! loop over densities ---
1644 :
1645 2 : CALL timestop(handle)
1646 :
1647 10 : END SUBROUTINE helium_print_rho
1648 :
1649 : ! ***************************************************************************
1650 : !> \brief Write volumetric data to an orthorhombic cubefile
1651 : !> \param unit unit number to which output will be sent
1652 : !> \param comment description of the data stored in the cubefile
1653 : !> \param origin position of the cubefile origin
1654 : !> \param deltar voxel side length
1655 : !> \param ndim number of voxels in each dimension
1656 : !> \param DATA array (ndim x ndim x ndim) with the data for output
1657 : !> \date 2013-11-25
1658 : !> \author Lukasz Walewski
1659 : ! **************************************************************************************************
1660 1 : SUBROUTINE helium_write_cubefile(unit, comment, origin, deltar, ndim, DATA)
1661 :
1662 : INTEGER, INTENT(IN) :: unit
1663 : CHARACTER(len=default_string_length), INTENT(IN) :: comment
1664 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: origin
1665 : REAL(KIND=dp), INTENT(IN) :: deltar
1666 : INTEGER, INTENT(IN) :: ndim
1667 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
1668 : POINTER :: DATA
1669 :
1670 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_write_cubefile'
1671 :
1672 : INTEGER :: handle, ix, iy, iz, nw
1673 : REAL(KIND=dp) :: delr, inva3
1674 : REAL(KIND=dp), DIMENSION(3) :: orig
1675 :
1676 1 : CALL timeset(routineN, handle)
1677 :
1678 : ! convert cubefile data to the proper units of measure
1679 1 : delr = angstrom*deltar
1680 4 : orig(:) = angstrom*origin(:)
1681 :
1682 : ! output cube file header
1683 1 : WRITE (unit, '(A)') comment
1684 1 : WRITE (unit, '(A)') "Volumetric data in cubefile format generated by CP2K"
1685 1 : WRITE (unit, '(I5,3(1X,F12.8))') 0, orig(1), orig(2), orig(3)
1686 1 : WRITE (unit, '(I5,3(1X,F12.8))') ndim, delr, 0.0_dp, 0.0_dp
1687 1 : WRITE (unit, '(I5,3(1X,F12.8))') ndim, 0.0_dp, delr, 0.0_dp
1688 1 : WRITE (unit, '(I5,3(1X,F12.8))') ndim, 0.0_dp, 0.0_dp, delr
1689 :
1690 : ! output cube file data
1691 1 : nw = 0
1692 1 : inva3 = 1.0_dp/(angstrom*angstrom*angstrom)
1693 11 : DO ix = 1, ndim
1694 111 : DO iy = 1, ndim
1695 1110 : DO iz = 1, ndim
1696 1000 : WRITE (unit, '(1X,E13.5)', ADVANCE='NO') inva3*DATA(ix, iy, iz)
1697 1000 : nw = nw + 1
1698 1100 : IF (MOD(nw, 6) .EQ. 0) THEN
1699 166 : nw = 0
1700 166 : WRITE (unit, *)
1701 : END IF
1702 : END DO
1703 : END DO
1704 : END DO
1705 : ! some compilers write over the first entry on a line losing all previous
1706 : ! values written on that line unless line terminator is written at the end
1707 : ! so make sure that the last WRITE statement runs without ADVANCE='NO'
1708 : ! (observed for ifort 14.0.1 and 14.0.2 but not for gfortran 4.8.2)
1709 1 : IF (nw .NE. 0) THEN
1710 1 : WRITE (unit, *)
1711 : END IF
1712 :
1713 1 : CALL timestop(handle)
1714 :
1715 1 : END SUBROUTINE helium_write_cubefile
1716 :
1717 : ! ***************************************************************************
1718 : !> \brief Print permutation length according to HELIUM%PRINT%PLENGTH
1719 : !> \param helium_env ...
1720 : !> \date 2010-06-07
1721 : !> \par History
1722 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
1723 : !> \author Lukasz Walewski
1724 : ! **************************************************************************************************
1725 105 : SUBROUTINE helium_print_plength(helium_env)
1726 :
1727 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1728 :
1729 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_plength'
1730 :
1731 : INTEGER :: handle, i, unit_nr
1732 : LOGICAL :: is_new, should_output
1733 : TYPE(cp_logger_type), POINTER :: logger
1734 : TYPE(section_vals_type), POINTER :: print_key
1735 :
1736 105 : CALL timeset(routineN, handle)
1737 :
1738 105 : NULLIFY (logger, print_key)
1739 105 : logger => cp_get_default_logger()
1740 :
1741 105 : IF (logger%para_env%is_source()) THEN
1742 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1743 55 : "MOTION%PINT%HELIUM%PRINT%PLENGTH")
1744 : should_output = BTEST(cp_print_key_should_output( &
1745 : iteration_info=logger%iter_info, &
1746 55 : basis_section=print_key), cp_p_file)
1747 :
1748 55 : IF (should_output) THEN
1749 :
1750 : unit_nr = cp_print_key_unit_nr( &
1751 : logger, &
1752 : print_key, &
1753 : middle_name="helium-plength", &
1754 : extension=".dat", &
1755 0 : is_new_file=is_new)
1756 :
1757 0 : DO i = 1, helium_env(1)%helium%atoms
1758 0 : WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%plength_avrg(i)
1759 0 : IF (i .LT. helium_env(1)%helium%atoms) THEN
1760 0 : WRITE (unit_nr, '(1X)', ADVANCE='NO')
1761 : END IF
1762 : END DO
1763 0 : WRITE (unit_nr, '(A)') ""
1764 :
1765 0 : CALL m_flush(unit_nr)
1766 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1767 :
1768 : END IF
1769 : END IF
1770 :
1771 105 : CALL timestop(handle)
1772 :
1773 105 : RETURN
1774 : END SUBROUTINE helium_print_plength
1775 :
1776 : ! ***************************************************************************
1777 : !> \brief Print helium force according to HELIUM%PRINT%FORCE
1778 : !> \param helium_env ...
1779 : !> \date 2010-01-27
1780 : !> \par History
1781 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
1782 : !> \author Lukasz Walewski
1783 : ! **************************************************************************************************
1784 105 : SUBROUTINE helium_print_force(helium_env)
1785 :
1786 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1787 :
1788 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_force'
1789 :
1790 : CHARACTER(len=default_string_length) :: msgstr
1791 : INTEGER :: handle, ia, ib, ic, idim, unit_nr
1792 : LOGICAL :: should_output
1793 : TYPE(cp_logger_type), POINTER :: logger
1794 : TYPE(section_vals_type), POINTER :: print_key
1795 :
1796 105 : CALL timeset(routineN, handle)
1797 :
1798 105 : NULLIFY (logger, print_key)
1799 105 : logger => cp_get_default_logger()
1800 :
1801 : print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1802 105 : "MOTION%PINT%HELIUM%PRINT%FORCES")
1803 : should_output = BTEST(cp_print_key_should_output( &
1804 : logger%iter_info, &
1805 105 : basis_section=print_key), cp_p_file)
1806 :
1807 105 : IF (.NOT. should_output) THEN
1808 105 : CALL timestop(handle)
1809 105 : RETURN
1810 : END IF
1811 :
1812 : ! check if there is anything to be printed out
1813 0 : IF (.NOT. helium_env(1)%helium%solute_present) THEN
1814 0 : msgstr = "Warning: force printout requested but there is no solute!"
1815 0 : CALL helium_write_line(msgstr)
1816 0 : CALL timestop(handle)
1817 0 : RETURN
1818 : END IF
1819 :
1820 0 : IF (logger%para_env%is_source()) THEN
1821 :
1822 : unit_nr = cp_print_key_unit_nr( &
1823 : logger, &
1824 : print_key, &
1825 : middle_name="helium-force", &
1826 0 : extension=".dat")
1827 :
1828 : ! print all force components in one line
1829 0 : DO ib = 1, helium_env(1)%helium%solute_beads
1830 0 : idim = 0
1831 0 : DO ia = 1, helium_env(1)%helium%solute_atoms
1832 0 : DO ic = 1, 3
1833 0 : idim = idim + 1
1834 0 : WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%force_avrg(ib, idim)
1835 : END DO
1836 : END DO
1837 : END DO
1838 0 : WRITE (unit_nr, *)
1839 :
1840 0 : CALL m_flush(unit_nr)
1841 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1842 :
1843 : END IF
1844 :
1845 0 : CALL timestop(handle)
1846 :
1847 0 : RETURN
1848 : END SUBROUTINE helium_print_force
1849 :
1850 : #if 0
1851 :
1852 : ! ***************************************************************************
1853 : !> \brief Print instantaneous force according to HELIUM%PRINT%FORCES_INST
1854 : !> \param helium ...
1855 : !> \date 2010-01-29
1856 : !> \author Lukasz Walewski
1857 : !> \note Collects instantaneous helium forces from all processors on
1858 : !> logger%para_env%source and writes them to files - one file per processor.
1859 : !> !TODO: Generalize for helium_env
1860 : ! **************************************************************************************************
1861 : SUBROUTINE helium_print_force_inst(helium)
1862 :
1863 : TYPE(helium_solvent_type), POINTER :: helium
1864 :
1865 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_force_inst', &
1866 : routineP = moduleN//':'//routineN
1867 :
1868 : CHARACTER(len=default_string_length) :: my_middle_name, stmp
1869 : INTEGER :: handle, ia, ib, ic, idim, irank, offset, &
1870 : unit_nr
1871 : LOGICAL :: should_output
1872 : TYPE(cp_logger_type), POINTER :: logger
1873 : TYPE(section_vals_type), POINTER :: print_key
1874 :
1875 : CALL timeset(routineN, handle)
1876 :
1877 : NULLIFY (logger, print_key)
1878 : logger => cp_get_default_logger()
1879 : print_key => section_vals_get_subs_vals(helium%input, &
1880 : "MOTION%PINT%HELIUM%PRINT%FORCES_INST")
1881 : should_output = BTEST(cp_print_key_should_output( &
1882 : logger%iter_info, &
1883 : basis_section=print_key), cp_p_file)
1884 :
1885 : IF (should_output) THEN
1886 :
1887 : ! check if there is anything to be printed out
1888 : IF (.NOT. helium%solute_present) THEN
1889 : stmp = "Warning: force printout requested but there is no solute!"
1890 : CALL helium_write_line(stmp)
1891 : CALL timestop(handle)
1892 : RETURN
1893 : END IF
1894 :
1895 : ! fill the tmp buffer with instantaneous helium forces at each proc
1896 : helium%rtmp_p_ndim_1d(:) = PACK(helium%force_inst, .TRUE.)
1897 :
1898 : ! pass the message from all processors to logger%para_env%source
1899 : helium%rtmp_p_ndim_np_1d(:) = 0.0_dp
1900 : CALL logger%para_env%gather(helium%rtmp_p_ndim_1d, helium%rtmp_p_ndim_np_1d)
1901 :
1902 : IF (logger%para_env%is_source()) THEN
1903 :
1904 : ! iterate over processors/helium environments
1905 : DO irank = 1, helium%num_env
1906 :
1907 : ! generate one file per processor
1908 : stmp = ""
1909 : WRITE (stmp, *) irank
1910 : my_middle_name = "helium-force-inst-"//TRIM(ADJUSTL(stmp))
1911 : unit_nr = cp_print_key_unit_nr( &
1912 : logger, &
1913 : print_key, &
1914 : middle_name=TRIM(my_middle_name), &
1915 : extension=".dat")
1916 :
1917 : ! unpack and actually print the forces - all components in one line
1918 : offset = (irank - 1)*SIZE(helium%rtmp_p_ndim_1d)
1919 : idim = 0
1920 : DO ib = 1, helium%solute_beads
1921 : DO ia = 1, helium%solute_atoms
1922 : DO ic = 1, 3
1923 : idim = idim + 1
1924 : WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium%rtmp_p_ndim_np_1d(offset + idim)
1925 : END DO
1926 : END DO
1927 : END DO
1928 : WRITE (unit_nr, *)
1929 :
1930 : CALL m_flush(unit_nr)
1931 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1932 :
1933 : END DO
1934 :
1935 : END IF
1936 : END IF
1937 :
1938 : CALL timestop(handle)
1939 :
1940 : END SUBROUTINE helium_print_force_inst
1941 :
1942 : #endif
1943 :
1944 : END MODULE helium_io
|
