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