(git:e7e05ae)
helium_io.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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  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  CALL timeset(routinen, handle)
115 
116  ! check if the file can be accessed
117  INQUIRE (file=file_name, exist=found, iostat=istat)
118  IF (istat /= 0) THEN
119  WRITE (unit=strtmp, fmt="(A,I0,A)") &
120  "An error occurred inquiring the file <"// &
121  trim(file_name)//">"
122  cpabort(trim(strtmp))
123  END IF
124  IF (.NOT. found) THEN
125  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  parse_white_lines=.true.)
133 
134  natom = 0
135  frame = 0
136  CALL parser_get_next_line(parser, 1)
137  frames: DO
138  ! Atom number
139  CALL parser_get_object(parser, natom)
140  frame = frame + 1
141  IF (frame == 1) THEN
142  ALLOCATE (coords(3*natom))
143  ELSE
144  strtmp = "Warning: more than one frame on file <"//trim(file_name)//">"
145  CALL helium_write_line(strtmp)
146  strtmp = "Warning: using the first frame only!"
147  CALL helium_write_line(strtmp)
148  EXIT frames
149  END IF
150  ! Dummy line
151  CALL parser_get_next_line(parser, 2)
152  DO j = 1, natom
153  ! Atom coordinates
154  READ (parser%input_line, *) strtmp, &
155  coords(3*(j - 1) + 1), &
156  coords(3*(j - 1) + 2), &
157  coords(3*(j - 1) + 3)
158  coords(3*(j - 1) + 1) = cp_unit_to_cp2k(coords(3*(j - 1) + 1), "angstrom")
159  coords(3*(j - 1) + 2) = cp_unit_to_cp2k(coords(3*(j - 1) + 2), "angstrom")
160  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  CALL parser_get_next_line(parser, 1, at_end=my_end)
163  my_end = my_end .OR. (len_trim(parser%input_line) == 0)
164  IF (my_end) THEN
165  IF (j /= natom) THEN
166  cpabort("Error in XYZ format.")
167  END IF
168  EXIT frames
169  END IF
170  END DO
171  END DO frames
172  CALL parser_release(parser)
173 
174  CALL timestop(handle)
175 
176  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  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  NULLIFY (logger)
199  logger => cp_get_default_logger()
200  my_label = "HELIUM| "
201 
202  IF (logger%para_env%is_source()) THEN
203  unit_nr = cp_logger_get_default_unit_nr(logger)
204 
205  WRITE (unit_nr, *)
206  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
207  " Number of helium environments: ", helium%num_env
208 
209  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
210  " Number of solvent atoms: ", helium%atoms
211  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
212  " Number of solvent beads: ", helium%beads
213  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
214  " Total number of solvent particles: ", helium%atoms*helium%beads
215 
216  unit_str = "angstrom^-3"
217  rtmp = cp_unit_from_cp2k(helium%density, &
218  unit_str)
219  WRITE (unit_nr, '(T2,A,F12.6)') trim(my_label)//" Density ["// &
220  trim(unit_str)//"]:", rtmp
221 
222  unit_str = "a.m.u."
223  rtmp = helium%he_mass_au/massunit
224  WRITE (unit_nr, '(T2,A,F12.6)') trim(my_label)//" He Mass ["// &
225  trim(unit_str)//"]: ", rtmp
226 
227  unit_str = "angstrom"
228  rtmp = cp_unit_from_cp2k(helium%cell_size, &
229  unit_str)
230  WRITE (unit_nr, '(T2,A,F12.6)') trim(my_label)//" Cell size ["// &
231  trim(unit_str)//"]: ", rtmp
232 
233  IF (helium%periodic) THEN
234  SELECT CASE (helium%cell_shape)
235  CASE (helium_cell_shape_cube)
236  CALL helium_write_line("PBC cell shape: CUBE.")
237  CASE (helium_cell_shape_octahedron)
238  CALL helium_write_line("PBC cell shape: TRUNCATED OCTAHEDRON.")
239  CASE DEFAULT
240  CALL helium_write_line("*** Warning: unknown cell shape.")
241  END SELECT
242  ELSE
243  CALL helium_write_line("PBC turned off.")
244  END IF
245 
246  IF (helium%droplet_radius .LT. huge(1.0_dp)) THEN
247  WRITE (stmp, *) helium%droplet_radius*angstrom
248  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  " First MC step :", helium%first_step + 1
254  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
255  " Last MC step :", helium%last_step
256  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
257  " Total number of MC steps :", helium%num_steps
258  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
259  " Number of outer MC trials per step :", helium%iter_rot
260  i8tmp = int(helium%iter_rot, int_8)
261  IF (helium%sampling_method == helium_sampling_ceperley) THEN
262  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
263  " Number of inner MC trials per step :", helium%iter_norot
264  i8tmp = i8tmp*int(helium%iter_norot, int_8)
265  END IF
266  stmp = ""
267  WRITE (stmp, *) i8tmp
268  WRITE (unit_nr, '(T2,A)') trim(my_label)// &
269  " Total number of MC trials per step : "//trim(adjustl(stmp))
270  i8tmp = int(helium%num_steps, int_8)
271  i8tmp = i8tmp*int(helium%iter_rot, int_8)
272  IF (helium%sampling_method == helium_sampling_ceperley) THEN
273  i8tmp = i8tmp*int(helium%iter_norot, int_8)
274  END IF
275  stmp = ""
276  WRITE (stmp, *) i8tmp
277  WRITE (unit_nr, '(T2,A)') trim(my_label)// &
278  " Total number of MC trials : "//trim(adjustl(stmp))
279 
280  SELECT CASE (helium%sampling_method)
281 
282  CASE (helium_sampling_ceperley)
283 
284  ! permutation cycle length sampling
285  stmp = ""
286  CALL helium_write_line(stmp)
287  WRITE (stmp, *) helium%maxcycle
288  stmp2 = ""
289  WRITE (stmp2, *) "Using maximum permutation cycle length: "// &
290  trim(adjustl(stmp))
291  CALL helium_write_line(stmp2)
292  stmp = ""
293  WRITE (stmp, *) "Permutation cycle length distribution: "// &
294  trim(adjustl(m_dist_name(helium%m_dist_type)))
295  CALL helium_write_line(stmp)
296  stmp = ""
297  stmp1 = ""
298  WRITE (stmp1, *) helium%m_ratio
299  stmp2 = ""
300  WRITE (stmp2, *) helium%m_value
301  WRITE (stmp, *) "Using ratio "//trim(adjustl(stmp1))// &
302  " for M = "//trim(adjustl(stmp2))
303  CALL helium_write_line(stmp)
304  stmp = ""
305  CALL helium_write_line(stmp)
306 
307  CASE (helium_sampling_worm)
308 
309  stmp1 = ""
310  stmp2 = ""
311  CALL helium_write_line(stmp1)
312  WRITE (stmp1, *) helium%worm_centroid_drmax
313  WRITE (stmp2, *) "WORM| Centroid move max. displacement: "// &
314  trim(adjustl(stmp1))
315  CALL helium_write_line(stmp2)
316  stmp1 = ""
317  stmp2 = ""
318  WRITE (stmp1, *) helium%worm_staging_l
319  WRITE (stmp2, *) "WORM| Maximal staging length: "//trim(adjustl(stmp1))
320  CALL helium_write_line(stmp2)
321  stmp1 = ""
322  stmp2 = ""
323  WRITE (stmp1, *) helium%worm_open_close_scale
324  WRITE (stmp2, *) "WORM| Open/Close scaling: "//trim(adjustl(stmp1))
325  CALL helium_write_line(stmp2)
326  stmp1 = ""
327  stmp2 = ""
328  WRITE (stmp1, *) helium%worm_allow_open
329  WRITE (stmp2, *) "WORM| Open worm sector: "//trim(adjustl(stmp1))
330  CALL helium_write_line(stmp2)
331  stmp1 = ""
332  stmp2 = ""
333  WRITE (stmp1, *) helium%worm_show_statistics
334  WRITE (stmp2, *) "WORM| Print statistics: "//trim(adjustl(stmp1))
335  CALL helium_write_line(stmp2)
336  IF (helium%worm_max_open_cycles > 0 .AND. helium%worm_allow_open) THEN
337  stmp1 = ""
338  stmp2 = ""
339  WRITE (stmp1, *) helium%worm_max_open_cycles
340  WRITE (stmp2, *) "WORM| Max. nb of MC cycles in open sector: "//trim(adjustl(stmp1))
341  CALL helium_write_line(stmp2)
342  END IF
343  stmp1 = ""
344  CALL helium_write_line(stmp1)
345 
346  CASE DEFAULT
347  WRITE (msg_str, *) helium%sampling_method
348  msg_str = "Sampling method ("//trim(adjustl(msg_str))//") not supported."
349  cpabort(msg_str)
350 
351  END SELECT
352 
353  ! solute related data
354  IF (helium%solute_present) THEN
355  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
356  " Number of solute atoms: ", helium%solute_atoms
357  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
358  " Number of solute beads: ", helium%solute_beads
359  WRITE (unit_nr, '(T2,A,1X,I0)') trim(my_label)// &
360  " Total number of solute particles: ", helium%solute_atoms* &
361  helium%solute_beads
362 
363  stmp1 = ""
364  SELECT CASE (helium%solute_interaction)
365  CASE (helium_solute_intpot_none)
366  WRITE (stmp1, *) "NONE"
367  CASE (helium_solute_intpot_mwater)
368  WRITE (stmp1, *) "MWATER"
369  CASE (helium_solute_intpot_nnp)
370  WRITE (stmp1, *) "NNP"
371  CASE DEFAULT
372  WRITE (stmp1, *) "***UNKNOWN***"
373  END SELECT
374  WRITE (unit_nr, '(T2,A,A,A)') &
375  trim(my_label), &
376  " Solute interaction type: ", &
377  trim(adjustl(stmp1))
378 
379  CALL get_cell(helium%solute_cell, abc=my_abc)
380  unit_str = "angstrom"
381  v1 = cp_unit_from_cp2k(my_abc(1), unit_str)
382  v2 = cp_unit_from_cp2k(my_abc(2), unit_str)
383  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  trim(unit_str)//"]: ", v1, v2, v3
387  ELSE
388  WRITE (unit_nr, '(T2,A)') trim(my_label)//" Solute is NOT present"
389  END IF
390  END IF
391  CALL helium_write_line("")
392 
393  ! RDF related settings
394  IF (helium%rdf_present) THEN
395  WRITE (stmp, '(I6)') helium%rdf_num
396  CALL helium_write_line("RDF| number of centers: "//trim(stmp))
397  rtmp = cp_unit_from_cp2k(helium%rdf_delr, "angstrom")
398  WRITE (stmp, '(1X,F12.6)') rtmp
399  CALL helium_write_line("RDF| delr [angstrom] : "//trim(stmp))
400  rtmp = cp_unit_from_cp2k(helium%rdf_maxr, "angstrom")
401  WRITE (stmp, '(1X,F12.6)') rtmp
402  CALL helium_write_line("RDF| maxr [angstrom] : "//trim(stmp))
403  itmp = helium%rdf_nbin
404  WRITE (stmp, '(I6)') itmp
405  CALL helium_write_line("RDF| nbin : "//trim(stmp))
406  CALL helium_write_line("")
407  ELSE
408  CALL helium_write_line("RDF| radial distributions will NOT be calculated.")
409  CALL helium_write_line("")
410  END IF
411 
412  ! RHO related settings
413  IF (helium%rho_present) THEN
414  CALL helium_write_line('RHO| The following densities will be calculated:')
415  DO i = 1, rho_num
416  IF (helium%rho_property(i)%is_calculated) THEN
417  WRITE (stmp, '(A)') 'RHO| '//trim(helium%rho_property(i)%name)
418  CALL helium_write_line(trim(stmp))
419  END IF
420  END DO
421  CALL helium_write_line('RHO|')
422  rtmp = cp_unit_from_cp2k(helium%rho_delr, "angstrom")
423  WRITE (stmp, '(1X,F12.6)') rtmp
424  CALL helium_write_line("RHO| delr [angstrom] : "//trim(stmp))
425  rtmp = cp_unit_from_cp2k(helium%rho_maxr, "angstrom")
426  WRITE (stmp, '(1X,F12.6)') rtmp
427  CALL helium_write_line("RHO| maxr [angstrom] : "//trim(stmp))
428  itmp = helium%rho_nbin
429  WRITE (stmp, '(I6)') itmp
430  CALL helium_write_line("RHO| nbin : "//trim(stmp))
431  CALL helium_write_line("")
432  ELSE
433  CALL helium_write_line("RHO| Density distributions will NOT be calculated.")
434  CALL helium_write_line("")
435  END IF
436 
437  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  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  NULLIFY (logger)
455  logger => cp_get_default_logger()
456  my_label = "HELIUM|"
457 
458  IF (logger%para_env%is_source()) THEN
459  unit_nr = cp_logger_get_default_unit_nr(logger)
460  WRITE (unit_nr, '(T2,A)') trim(my_label)//" "//trim(line)
461  END IF
462 
463  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: my_energy
484  TYPE(cp_logger_type), POINTER :: logger
485  TYPE(section_vals_type), POINTER :: print_key
486 
487  CALL timeset(routinen, handle)
488 
489  NULLIFY (logger, print_key)
490  logger => cp_get_default_logger()
491 
492  IF (logger%para_env%is_source()) THEN
493  print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
494  "MOTION%PINT%HELIUM%PRINT%ENERGY")
495  should_output = btest(cp_print_key_should_output( &
496  iteration_info=logger%iter_info, &
497  basis_section=print_key), cp_p_file)
498  cepsample = helium_env(1)%helium%sampling_method == helium_sampling_ceperley
499  END IF
500  CALL helium_env(1)%comm%bcast(should_output, logger%para_env%source)
501  CALL helium_env(1)%comm%bcast(cepsample, logger%para_env%source)
502 
503  IF (should_output) THEN
504 
505  naccptd = 0.0_dp
506  IF (cepsample) THEN
507  ntot = 0.0_dp
508  DO k = 1, SIZE(helium_env)
509  ntot = ntot + helium_env(1)%helium%iter_norot*helium_env(1)%helium%iter_rot
510  DO m = 1, helium_env(k)%helium%maxcycle
511  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  ntot = 0.0_dp
516  DO k = 1, SIZE(helium_env)
517  naccptd = naccptd + helium_env(k)%helium%num_accepted(1, 1)
518  ntot = ntot + helium_env(k)%helium%num_accepted(2, 1)
519  END DO
520  END IF
521  CALL helium_env(1)%comm%sum(naccptd)
522  CALL helium_env(1)%comm%sum(ntot)
523 
524  IF (logger%para_env%is_source()) THEN
525  my_energy => helium_env(1)%helium%energy_avrg
526  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  is_new_file=file_is_new)
534 
535  IF (file_is_new) THEN
536  WRITE (unit_nr, '(A9,7(1X,A20))') &
537  "# Step", &
538  " AccRatio", &
539  " E_pot", &
540  " E_kin", &
541  " E_thermo", &
542  " E_virial", &
543  " E_inter", &
544  " E_tot"
545  END IF
546 
547  WRITE (unit_nr, "(I9,7(1X,F20.9))") &
548  iteration, &
549  naccptd/ntot, &
550  my_energy(e_id_potential), &
551  my_energy(e_id_kinetic), &
552  my_energy(e_id_thermo), &
553  my_energy(e_id_virial), &
554  my_energy(e_id_interact), &
555  my_energy(e_id_total)
556 
557  CALL m_flush(unit_nr)
558  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
559 
560  END IF
561  END IF
562 
563  CALL timestop(handle)
564 
565  RETURN
566  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: data_p
603  TYPE(cp_logger_type), POINTER :: logger
604  TYPE(section_vals_type), POINTER :: psection
605 
606  CALL timeset(routinen, handle)
607 
608  IF (PRESENT(avg)) THEN
609  my_avg = avg
610  ELSE
611  my_avg = .false.
612  END IF
613 
614  IF (PRESENT(fpos)) THEN
615  my_fpos = fpos
616  ELSE
617  my_fpos = "APPEND"
618  END IF
619 
620  NULLIFY (logger, psection)
621  logger => cp_get_default_logger()
622 
623  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  basis_section=psection), cp_p_file)
627 
628  IF (.NOT. should_output) THEN
629  CALL timestop(handle)
630  RETURN
631  END IF
632 
633  IF (my_avg) THEN
634  ! average data over all processors and environments
635  avg_data(:) = 0.0_dp
636  msglen = SIZE(avg_data)
637  DO i = 0, SIZE(helium_env) - 1
638  avg_data(:) = avg_data(:) + DATA(msglen*i + 1:msglen*(i + 1))
639  END DO
640  CALL helium_env(1)%comm%sum(avg_data)
641  nenv = helium_env(1)%helium%num_env
642  avg_data(:) = avg_data(:)/real(nenv, dp)
643  ELSE
644  ! gather data from all processors
645  offset = 0
646  DO i = 1, logger%para_env%mepos
647  offset = offset + helium_env(1)%env_all(i)
648  END DO
649 
650  helium_env(1)%helium%rtmp_3_np_1d = 0.0_dp
651  msglen = SIZE(avg_data)
652  DO i = 0, SIZE(helium_env) - 1
653  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  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  is_new_file=is_new)
665 
666  IF (logger%para_env%is_source()) THEN
667 
668  IF (is_new) THEN
669  ! comment
670  IF (cmmnt .NE. "") THEN
671  WRITE (unit_nr, '(A)') "# "//cmmnt
672  END IF
673  ! column labels
674  WRITE (unit_nr, '(A2,A18,1X,A20,1X,A20)') &
675  "# ", &
676  col_label(1), &
677  col_label(2), &
678  col_label(3)
679  END IF
680 
681  IF (my_avg) THEN
682  DO i = 1, 3
683  WRITE (unit_nr, '(E27.20)', advance='NO') uconv*avg_data(i)
684  IF (i .LT. 3) THEN
685  WRITE (unit_nr, '(1X)', advance='NO')
686  END IF
687  END DO
688  WRITE (unit_nr, '(A)') ""
689  ELSE
690  ! iterate over processors/helium environments
691  DO irank = 1, helium_env(1)%helium%num_env
692  ! unpack data (actually point to the right fragment only)
693  msglen = SIZE(avg_data)
694  offset = (irank - 1)*msglen
695  data_p => helium_env(1)%helium%rtmp_3_np_1d(offset + 1:offset + msglen)
696  ! write out the data
697  DO i = 1, 3
698  WRITE (unit_nr, '(E27.20)', advance='NO') uconv*data_p(i)
699  IF (i .LT. 3) THEN
700  WRITE (unit_nr, '(1X)', advance='NO')
701  END IF
702  END DO
703  WRITE (unit_nr, '(A)') ""
704  END DO
705  END IF
706 
707  CALL m_flush(unit_nr)
708  CALL cp_print_key_finished_output(unit_nr, logger, psection)
709 
710  END IF
711 
712  CALL timestop(handle)
713 
714  RETURN
715  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  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  CALL timeset(routinen, handle)
737 
738  NULLIFY (logger, print_key)
739  logger => cp_get_default_logger()
740 
741  IF (logger%para_env%is_source()) THEN
742  print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
743  "MOTION%PINT%HELIUM%PRINT%ACCEPTS")
744  should_output = btest(cp_print_key_should_output( &
745  iteration_info=logger%iter_info, &
746  basis_section=print_key), cp_p_file)
747 
748  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  is_new_file=file_is_new)
755 
756  IF (file_is_new) THEN
757  WRITE (unit_nr, '(A8,1X,A15,1X,A20)', advance='NO') &
758  "# Length", &
759  " Trials", &
760  " Selected"
761  DO j = 1, helium_env(1)%helium%bisctlog2
762  WRITE (unit_nr, '(A17,1X,I3)', advance='NO') " Level", j
763  END DO
764  WRITE (unit_nr, '(A)') ""
765  END IF
766 
767  DO i = 1, helium_env(1)%helium%maxcycle
768  WRITE (unit_nr, '(I3)', advance='NO') i
769  DO j = 1, helium_env(1)%helium%bisctlog2 + 2
770  WRITE (unit_nr, '(1X,F20.2)', advance='NO') helium_env(1)%helium%num_accepted(j, i)
771  END DO
772  WRITE (unit_nr, '(A)') ""
773  END DO
774  WRITE (unit_nr, '(A)') "&"
775 
776  CALL m_flush(unit_nr)
777  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
778 
779  END IF
780  END IF
781 
782  CALL timestop(handle)
783  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  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  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  CALL timeset(routinen, handle)
810 
811  NULLIFY (logger, print_key)
812  NULLIFY (used_indices)
813  logger => cp_get_default_logger()
814 
815  ! determine offset for arrays
816  offset = 0
817  DO i = 1, logger%para_env%mepos
818  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  "MOTION%PINT%HELIUM%PRINT%PERM")
823  should_output = btest(cp_print_key_should_output( &
824  iteration_info=logger%iter_info, &
825  basis_section=print_key), cp_p_file)
826 
827  IF (.NOT. should_output) THEN
828  CALL timestop(handle)
829  RETURN
830  END IF
831 
832  ! get the output type
833  CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
834  IF (outformat .EQ. perm_cycle) THEN
835  ! gather positions from all processors
836  helium_env(1)%helium%rtmp_3_atoms_beads_np_1d = 0.0_dp
837  j = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
838  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  pack(helium_env(i)%helium%pos(:, :, 1:helium_env(1)%helium%beads), .true.)
841  END DO
842  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  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  helium_env(1)%helium%itmp_atoms_np_1d(:) = 0
849  msglen = SIZE(helium_env(1)%helium%permutation)
850  DO i = 1, SIZE(helium_env)
851  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  CALL helium_env(1)%comm%sum(helium_env(1)%helium%itmp_atoms_np_1d)
855 
856  IF (logger%para_env%is_source()) THEN
857 
858  ! iterate over helium environments
859  DO irank = 1, helium_env(1)%helium%num_env
860 
861  ! generate one file per environment
862  stmp = ""
863  WRITE (stmp, *) irank
864  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  extension=".dat")
870 
871  ! unpack permutation state (actually point to the right section only)
872  lb = (irank - 1)*helium_env(1)%helium%atoms + 1
873  ub = irank*helium_env(1)%helium%atoms
874  my_perm => helium_env(1)%helium%itmp_atoms_np_1d(lb:ub)
875 
876  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  msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
883  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  mask=helium_env(1)%helium%ltmp_3_atoms_beads_3d, field=0.0_dp)
887 
888  curat = 1
889  nused = 0
890  first_cycle = .true.
891  DO WHILE (curat .LE. helium_env(1)%helium%atoms)
892 
893  ! get the permutation cycle the current atom belongs to
894  my_cycle => helium_cycle_of(curat, my_perm)
895 
896  ! include the current cycle in the pool of "used" indices
897  nused = nused + SIZE(my_cycle)
898  CALL reallocate(used_indices, 1, nused)
899  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  IF (helium_is_winding(helium_env(1)%helium, curat, helium_env(1)%helium%work, my_perm)) THEN
903  left_delim = "["
904  right_delim = "]"
905  ELSE
906  left_delim = "("
907  right_delim = ")"
908  END IF
909 
910  ! cycle delimiter
911  IF (first_cycle) THEN
912  first_cycle = .false.
913  ELSE
914  WRITE (unit_nr, '(1X)', advance='NO')
915  END IF
916 
917  ! write out the current cycle
918  WRITE (unit_nr, '(A1)', advance='NO') left_delim
919  first = .true.
920  DO i = 1, SIZE(my_cycle)
921  IF (first) THEN
922  first = .false.
923  ELSE
924  WRITE (unit_nr, '(1X)', advance='NO')
925  END IF
926  WRITE (unit_nr, '(I0)', advance='NO') my_cycle(i)
927  END DO
928  WRITE (unit_nr, '(A1)', advance='NO') right_delim
929 
930  ! clean up
931  DEALLOCATE (my_cycle)
932  NULLIFY (my_cycle)
933 
934  ! try to increment the current atom index
935  DO WHILE (any(used_indices .EQ. curat))
936  curat = curat + 1
937  END DO
938 
939  END DO
940  WRITE (unit_nr, *)
941 
942  DEALLOCATE (used_indices)
943  NULLIFY (used_indices)
944 
945  CASE (perm_plain)
946  ! write the plain permutation
947 
948  first = .true.
949  DO i = 1, helium_env(1)%helium%atoms
950  IF (first) THEN
951  first = .false.
952  ELSE
953  WRITE (unit_nr, '(1X)', advance='NO')
954  END IF
955  WRITE (unit_nr, '(I0)', advance='NO') my_perm(i)
956  END DO
957  WRITE (unit_nr, '(A)') ""
958 
959  CASE default
960 
961  msg_str = "Format for permutation output unknown."
962  cpabort(msg_str)
963 
964  END SELECT
965 
966  CALL m_flush(unit_nr)
967  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
968 
969  END DO
970 
971  END IF
972 
973  CALL timestop(handle)
974 
975  RETURN
976  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: tmp_inter_action, tmp_link_action, &
999  tmp_pair_action
1000  TYPE(cp_logger_type), POINTER :: logger
1001  TYPE(section_vals_type), POINTER :: print_key
1002 
1003  CALL timeset(routinen, handle)
1004 
1005  NULLIFY (logger, print_key)
1006  logger => cp_get_default_logger()
1007 
1008  ! determine offset for arrays
1009  offset = 0
1010  DO i = 1, logger%para_env%mepos
1011  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  "MOTION%PINT%HELIUM%PRINT%ACTION")
1016  should_output = btest(cp_print_key_should_output( &
1017  iteration_info=logger%iter_info, &
1018  basis_section=print_key), cp_p_file)
1019 
1020  IF (.NOT. should_output) THEN
1021  CALL timestop(handle)
1022  RETURN
1023  END IF
1024 
1025  DO k = 1, SIZE(helium_env)
1026  helium_env(k)%helium%link_action = helium_total_link_action(helium_env(k)%helium)
1027  helium_env(k)%helium%pair_action = helium_total_pair_action(helium_env(k)%helium)
1028  helium_env(k)%helium%inter_action = helium_total_inter_action(pint_env, helium_env(k)%helium)
1029  END DO
1030 
1031  NULLIFY (tmp_link_action)
1032  NULLIFY (tmp_pair_action)
1033  NULLIFY (tmp_inter_action)
1034  ALLOCATE (tmp_link_action(helium_env(1)%helium%num_env))
1035  ALLOCATE (tmp_pair_action(helium_env(1)%helium%num_env))
1036  ALLOCATE (tmp_inter_action(helium_env(1)%helium%num_env))
1037  tmp_link_action(:) = 0.0_dp
1038  tmp_pair_action(:) = 0.0_dp
1039  tmp_inter_action(:) = 0.0_dp
1040  ! gather Action from all processors to logger%para_env%source
1041  DO k = 1, SIZE(helium_env)
1042  tmp_link_action(offset + k) = helium_env(k)%helium%link_action
1043  tmp_pair_action(offset + k) = helium_env(k)%helium%pair_action
1044  tmp_inter_action(offset + k) = helium_env(k)%helium%inter_action
1045  END DO
1046  CALL helium_env(1)%comm%sum(tmp_link_action)
1047  CALL helium_env(1)%comm%sum(tmp_pair_action)
1048  CALL helium_env(1)%comm%sum(tmp_inter_action)
1049 
1050  IF (logger%para_env%is_source()) THEN
1051  iteration = logger%iter_info%iteration(2)
1052  ! iterate over processors/helium environments
1053  DO irank = 1, helium_env(1)%helium%num_env
1054 
1055  ! generate one file per helium_env
1056  stmp = ""
1057  WRITE (stmp, *) irank
1058  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  is_new_file=file_is_new)
1065 
1066  IF (file_is_new) THEN
1067  WRITE (unit_nr, '(A9,3(1X,A25))') &
1068  "# Step", &
1069  " He_Total_Link_Action", &
1070  " He_Total_Pair_Action", &
1071  " He_Total_Interaction"
1072  END IF
1073 
1074  WRITE (unit_nr, "(I9,3(1X,F25.14))") &
1075  iteration, &
1076  tmp_link_action(irank), &
1077  tmp_pair_action(irank), &
1078  tmp_inter_action(irank)
1079 
1080  CALL m_flush(unit_nr)
1081  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1082 
1083  END DO
1084  END IF
1085 
1086  DEALLOCATE (tmp_link_action)
1087  DEALLOCATE (tmp_pair_action)
1088  DEALLOCATE (tmp_inter_action)
1089  NULLIFY (tmp_link_action)
1090  NULLIFY (tmp_pair_action)
1091  NULLIFY (tmp_inter_action)
1092 
1093  CALL timestop(handle)
1094 
1095  RETURN
1096  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  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  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  CALL timeset(routinen, handle)
1126 
1127  NULLIFY (logger, print_key)
1128  logger => cp_get_default_logger()
1129 
1130  ! determine offset for arrays
1131  offset = 0
1132  DO i = 1, logger%para_env%mepos
1133  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  "MOTION%PINT%HELIUM%PRINT%COORDINATES")
1138  should_output = btest(cp_print_key_should_output( &
1139  iteration_info=logger%iter_info, &
1140  basis_section=print_key), cp_p_file)
1141 
1142  IF (.NOT. should_output) THEN
1143  CALL timestop(handle)
1144  RETURN
1145  END IF
1146 
1147  ! prepare the coordinates for output (use unit cell centered around r0)
1148  DO k = 1, SIZE(helium_env)
1149  r0(:) = helium_env(k)%helium%center(:)
1150  DO ia = 1, helium_env(k)%helium%atoms
1151  DO ib = 1, helium_env(k)%helium%beads
1152  r1(:) = helium_env(k)%helium%pos(:, ia, ib) - r0(:)
1153  r2(:) = helium_env(k)%helium%pos(:, ia, ib) - r0(:)
1154  CALL helium_pbc(helium_env(k)%helium, r2)
1155  ltmp = .false.
1156  DO ic = 1, 3
1157  IF (abs(r1(ic) - r2(ic)) .GT. 100.0_dp*epsilon(0.0_dp)) THEN
1158  ltmp = .true.
1159  cycle
1160  END IF
1161  END DO
1162  IF (ltmp) THEN
1163  helium_env(k)%helium%work(:, ia, ib) = r0(:) + r2(:)
1164  ELSE
1165  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  helium_env(1)%helium%rtmp_3_atoms_beads_np_1d = 0.0_dp
1173  j = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
1174  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  pack(helium_env(i)%helium%pos(:, :, 1:helium_env(1)%helium%beads), .true.)
1177  END DO
1178  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  helium_env(1)%helium%itmp_atoms_np_1d(:) = 0
1182  j = SIZE(helium_env(1)%helium%permutation)
1183  DO i = 1, SIZE(helium_env)
1184  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  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  helium_env(1)%helium%ltmp_3_atoms_beads_3d(:, :, :) = .true.
1191 
1192  IF (logger%para_env%is_source()) THEN
1193 
1194  CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
1195 
1196  ! iterate over helium environments
1197  DO irank = 1, helium_env(1)%helium%num_env
1198  IF (outformat == fmt_id_pdb) THEN
1199  ! generate one file per environment
1200  stmp = ""
1201  WRITE (stmp, *) irank
1202  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  extension=".pdb")
1208 
1209  ! write out the unit cell parameters
1210  fmt_string = "(A6,3F9.3,3F7.2,1X,A11,1X,I3)"
1211  xtmp = helium_env(1)%helium%cell_size
1212  xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
1213  SELECT CASE (helium_env(1)%helium%cell_shape)
1214  CASE (helium_cell_shape_cube)
1215  stmp = "C "
1216  CASE (helium_cell_shape_octahedron)
1217  stmp = "O "
1218  CASE DEFAULT
1219  stmp = "UNKNOWN "
1220  END SELECT
1221  WRITE (unit_nr, fmt_string) "CRYST1", &
1222  xtmp, xtmp, xtmp, &
1223  90.0_dp, 90.0_dp, 90.0_dp, &
1224  stmp, helium_env(1)%helium%beads
1225 
1226  ! unpack coordinates
1227  msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
1228  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  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  msglen = SIZE(helium_env(1)%helium%permutation)
1235  offset = (irank - 1)*msglen
1236  my_perm => helium_env(1)%helium%itmp_atoms_np_1d(offset + 1:offset + msglen)
1237 
1238  ! write out coordinates
1239  fmt_string = &
1240  "(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2)"
1241  DO ia = 1, helium_env(1)%helium%atoms
1242  icycle = helium_cycle_number(helium_env(1)%helium, ia, my_perm)
1243  is_winding = helium_is_winding(helium_env(1)%helium, ia, helium_env(1)%helium%work, my_perm)
1244  IF (is_winding) THEN
1245  resname = "SPR"
1246  ELSE
1247  resname = "NRM"
1248  END IF
1249  DO ib = 1, helium_env(1)%helium%beads
1250  xtmp = helium_env(1)%helium%work(1, ia, ib)
1251  xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
1252  ytmp = helium_env(1)%helium%work(2, ia, ib)
1253  ytmp = cp_unit_from_cp2k(ytmp, "angstrom")
1254  ztmp = helium_env(1)%helium%work(3, ia, ib)
1255  ztmp = cp_unit_from_cp2k(ztmp, "angstrom")
1256  WRITE (unit_nr, fmt_string) "ATOM ", &
1257  (ia - 1)*helium_env(1)%helium%beads + ib, &
1258  " He ", " ", resname, "X", &
1259  icycle, &
1260  " ", &
1261  xtmp, ytmp, ztmp, &
1262  1.0_dp, 0.0_dp, "HE", " "
1263  END DO
1264  END DO
1265 
1266  ! write out the bead connectivity information
1267  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  IF (helium_path_length(helium_env(1)%helium, ia, my_perm) .LE. 1) THEN
1272  cycle
1273  END IF
1274 
1275  DO ib = 1, helium_env(1)%helium%beads - 1
1276  ! check whether the consecutive beads belong to the same box
1277  r1(:) = helium_env(1)%helium%work(:, ia, ib) - helium_env(1)%helium%work(:, ia, ib + 1)
1278  r2(:) = r1(:)
1279  CALL helium_pbc(helium_env(1)%helium, r2)
1280  are_connected = .true.
1281  DO ic = 1, 3
1282  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  are_connected = .false.
1286  cycle
1287  END IF
1288  END DO
1289  IF (are_connected) THEN
1290  tmp1 = (ia - 1)*helium_env(1)%helium%beads + ib
1291  tmp2 = (ia - 1)*helium_env(1)%helium%beads + ib + 1
1292  ! smaller value has to go first
1293  IF (tmp1 .LT. tmp2) THEN
1294  ib1 = tmp1
1295  ib2 = tmp2
1296  ELSE
1297  ib1 = tmp2
1298  ib2 = tmp1
1299  END IF
1300  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  r1(:) = helium_env(1)%helium%work(:, ia, helium_env(1)%helium%beads) - helium_env(1)%helium%work(:, my_perm(ia), 1)
1307  r2(:) = r1(:)
1308  CALL helium_pbc(helium_env(1)%helium, r2)
1309  are_connected = .true.
1310  DO ic = 1, 3
1311  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  are_connected = .false.
1315  cycle
1316  END IF
1317  END DO
1318  IF (are_connected) THEN
1319  tmp1 = ia*helium_env(1)%helium%beads
1320  tmp2 = (my_perm(ia) - 1)*helium_env(1)%helium%beads + 1
1321  IF (tmp1 .LT. tmp2) THEN
1322  ib1 = tmp1
1323  ib2 = tmp2
1324  ELSE
1325  ib1 = tmp2
1326  ib2 = tmp1
1327  END IF
1328  WRITE (unit_nr, '(A6,2I5)') "CONECT", ib1, ib2
1329  END IF
1330  END DO
1331  WRITE (unit_nr, '(A)') "END"
1332 
1333  CALL m_flush(unit_nr)
1334  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1335  ELSE IF (outformat == fmt_id_xyz) THEN
1336  ! generate one file per environment and bead
1337  DO ib = 1, helium_env(1)%helium%beads
1338  stmp = ""
1339  WRITE (stmp, *) irank
1340  my_middle_name = "helium-pos-"//trim(adjustl(stmp))
1341  WRITE (stmp, *) ib
1342  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  extension=".xyz")
1348  ! write out xyz header
1349  WRITE (unit_nr, *) helium_env(1)%helium%atoms
1350  stmp = ""
1351  WRITE (stmp, *) logger%iter_info%n_rlevel
1352  fmt_string = "(A6,"//trim(adjustl(stmp))//"I12)"
1353  WRITE (unit_nr, fmt_string) "iter= ", logger%iter_info%iteration(:)
1354  fmt_string = "(A6,3F9.3,3F7.2,1X,A11,1X,I3)"
1355 
1356  ! unpack coordinates
1357  msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
1358  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  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  msglen = SIZE(helium_env(1)%helium%permutation)
1365  offset = (irank - 1)*msglen
1366  my_perm => helium_env(1)%helium%itmp_atoms_np_1d(offset + 1:offset + msglen)
1367 
1368  ! write out coordinates
1369  fmt_string = "(A2,3(1X,F20.10))"
1370  DO ia = 1, helium_env(1)%helium%atoms
1371  xtmp = helium_env(1)%helium%work(1, ia, ib)
1372  xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
1373  ytmp = helium_env(1)%helium%work(2, ia, ib)
1374  ytmp = cp_unit_from_cp2k(ytmp, "angstrom")
1375  ztmp = helium_env(1)%helium%work(3, ia, ib)
1376  ztmp = cp_unit_from_cp2k(ztmp, "angstrom")
1377  WRITE (unit_nr, fmt_string) "He", xtmp, ytmp, ztmp
1378  END DO
1379  CALL m_flush(unit_nr)
1380  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1381  END DO
1382  ELSE
1383  cpabort("")
1384  END IF
1385  END DO
1386 
1387  END IF
1388 
1389  CALL timestop(handle)
1390 
1391  RETURN
1392  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  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  REAL(kind=dp), DIMENSION(:, :), POINTER :: message
1414  TYPE(cp_logger_type), POINTER :: logger
1415  TYPE(section_vals_type), POINTER :: print_key
1416 
1417  CALL timeset(routinen, handle)
1418 
1419  NULLIFY (logger, print_key)
1420  logger => cp_get_default_logger()
1421 
1422  IF (logger%para_env%is_source()) THEN
1423  print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1424  "MOTION%PINT%HELIUM%PRINT%RDF")
1425  should_output = btest(cp_print_key_should_output( &
1426  iteration_info=logger%iter_info, &
1427  basis_section=print_key), cp_p_file)
1428  END IF
1429  CALL helium_env(1)%comm%bcast(should_output, logger%para_env%source)
1430 
1431  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  helium_env(1)%helium%rdf_inst(:, :) = 0.0_dp
1435  DO k = 1, SIZE(helium_env)
1436  helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :) + &
1437  helium_env(k)%helium%rdf_accu(:, :)
1438  END DO
1439 
1440  ! average over processors
1441  NULLIFY (message)
1442  message => helium_env(1)%helium%rdf_inst(:, :)
1443  CALL helium_env(1)%comm%sum(message)
1444  itmp = helium_env(1)%helium%num_env
1445  inv_norm = 1.0_dp/real(itmp, dp)
1446  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  nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1450  inv_norm = 1.0_dp/real(nsteps, dp)
1451  helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*inv_norm
1452 
1453  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  iweight*helium_env(1)%helium%rdf_rstr(:, :)
1457  helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/real(nsteps + iweight, dp)
1458 
1459  IF (logger%para_env%is_source()) THEN
1460 
1461  ia = 0
1462  rtmp = cp_unit_from_cp2k(1.0_dp, "angstrom")
1463  rtmp2 = 1.0_dp
1464  IF (.NOT. helium_env(1)%helium%periodic) THEN
1465  ! RDF in non-periodic case has unit 1/bohr**3, convert to Angstrom:
1466  rtmp2 = rtmp**(-3)
1467  END IF
1468 
1469  IF (helium_env(1)%helium%rdf_he_he) THEN
1470  ! overwrite RDF file each time it is written
1471  ia = 1
1472  stmp = ""
1473  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  do_backup=.false.)
1481 
1482  DO ic = 1, helium_env(1)%helium%rdf_nbin
1483  WRITE (unit_nr, '(F20.10)', advance='NO') (real(ic, dp) - 0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
1484  WRITE (unit_nr, '(F20.10)', advance='NO') helium_env(1)%helium%rdf_inst(ia, ic)*rtmp2
1485  WRITE (unit_nr, *)
1486  END DO
1487 
1488  CALL m_flush(unit_nr)
1489  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1490  END IF
1491 
1492  IF (helium_env(1)%helium%rdf_sol_he) THEN
1493  ! overwrite RDF file each time it is written
1494  stmp = ""
1495  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  do_backup=.false.)
1503 
1504  DO ic = 1, helium_env(1)%helium%rdf_nbin
1505  WRITE (unit_nr, '(F20.10)', advance='NO') (real(ic, dp) - 0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
1506  DO id = 1 + ia, helium_env(1)%helium%rdf_num
1507  WRITE (unit_nr, '(F20.10)', advance='NO') helium_env(1)%helium%rdf_inst(id, ic)*rtmp2
1508  END DO
1509  WRITE (unit_nr, *)
1510  END DO
1511 
1512  CALL m_flush(unit_nr)
1513  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1514  END IF
1515 
1516  END IF
1517  END IF
1518 
1519  CALL timestop(handle)
1520 
1521  RETURN
1522  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  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  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: cubdata
1546  TYPE(cp_logger_type), POINTER :: logger
1547  TYPE(section_vals_type), POINTER :: print_key
1548 
1549  CALL timeset(routinen, handle)
1550 
1551  NULLIFY (logger, print_key)
1552  logger => cp_get_default_logger()
1553 
1554  print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1555  "MOTION%PINT%HELIUM%PRINT%RHO")
1556  should_output = btest(cp_print_key_should_output( &
1557  iteration_info=logger%iter_info, &
1558  basis_section=print_key), cp_p_file)
1559 
1560  IF (.NOT. should_output) THEN
1561  CALL timestop(handle)
1562  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  helium_env(1)%helium%rho_inst(:, :, :, :) = 0.0_dp
1568  DO k = 1, SIZE(helium_env)
1569  helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :) + &
1570  helium_env(k)%helium%rho_accu(:, :, :, :)
1571  END DO
1572 
1573  ! average over processors
1574  CALL helium_env(1)%comm%sum(helium_env(1)%helium%rho_inst)
1575  itmp = helium_env(1)%helium%num_env
1576  invproc = 1.0_dp/real(itmp, dp)
1577  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  nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
1581  inv_norm = 1.0_dp/real(nsteps, dp)
1582  helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)*inv_norm
1583 
1584  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  iweight*helium_env(1)%helium%rho_rstr(:, :, :, :)
1588  helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)/real(nsteps + iweight, dp)
1589 
1590  ! set center of the cubefile
1591  IF (helium_env(1)%helium%solute_present) THEN
1592  ! should be set to solute's COM
1593  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  center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
1599  END IF
1600 
1601  DO id = 1, rho_num ! loop over densities ---
1602 
1603  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  DO ic = 1, helium_env(1)%helium%rho_property(id)%num_components ! loop over components
1609 
1610  WRITE (fname, '(A)') "helium-rho-"// &
1611  trim(adjustl(helium_env(1)%helium%rho_property(id)%filename_suffix(ic)))
1612  IF (helium_env(1)%helium%rho_property(id)%component_name(ic) .EQ. "") THEN
1613  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  trim(helium_env(1)%helium%rho_property(id)%component_name(ic))
1618  END IF
1619  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  do_backup=.false.)
1628 
1629  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  cubdata)
1637  CALL m_flush(unit_nr)
1638  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  CALL timestop(handle)
1646 
1647  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  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  CALL timeset(routinen, handle)
1677 
1678  ! convert cubefile data to the proper units of measure
1679  delr = angstrom*deltar
1680  orig(:) = angstrom*origin(:)
1681 
1682  ! output cube file header
1683  WRITE (unit, '(A)') comment
1684  WRITE (unit, '(A)') "Volumetric data in cubefile format generated by CP2K"
1685  WRITE (unit, '(I5,3(1X,F12.8))') 0, orig(1), orig(2), orig(3)
1686  WRITE (unit, '(I5,3(1X,F12.8))') ndim, delr, 0.0_dp, 0.0_dp
1687  WRITE (unit, '(I5,3(1X,F12.8))') ndim, 0.0_dp, delr, 0.0_dp
1688  WRITE (unit, '(I5,3(1X,F12.8))') ndim, 0.0_dp, 0.0_dp, delr
1689 
1690  ! output cube file data
1691  nw = 0
1692  inva3 = 1.0_dp/(angstrom*angstrom*angstrom)
1693  DO ix = 1, ndim
1694  DO iy = 1, ndim
1695  DO iz = 1, ndim
1696  WRITE (unit, '(1X,E13.5)', advance='NO') inva3*DATA(ix, iy, iz)
1697  nw = nw + 1
1698  IF (mod(nw, 6) .EQ. 0) THEN
1699  nw = 0
1700  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  IF (nw .NE. 0) THEN
1710  WRITE (unit, *)
1711  END IF
1712 
1713  CALL timestop(handle)
1714 
1715  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  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  CALL timeset(routinen, handle)
1737 
1738  NULLIFY (logger, print_key)
1739  logger => cp_get_default_logger()
1740 
1741  IF (logger%para_env%is_source()) THEN
1742  print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1743  "MOTION%PINT%HELIUM%PRINT%PLENGTH")
1744  should_output = btest(cp_print_key_should_output( &
1745  iteration_info=logger%iter_info, &
1746  basis_section=print_key), cp_p_file)
1747 
1748  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  is_new_file=is_new)
1756 
1757  DO i = 1, helium_env(1)%helium%atoms
1758  WRITE (unit_nr, '(F20.10)', advance='NO') helium_env(1)%helium%plength_avrg(i)
1759  IF (i .LT. helium_env(1)%helium%atoms) THEN
1760  WRITE (unit_nr, '(1X)', advance='NO')
1761  END IF
1762  END DO
1763  WRITE (unit_nr, '(A)') ""
1764 
1765  CALL m_flush(unit_nr)
1766  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1767 
1768  END IF
1769  END IF
1770 
1771  CALL timestop(handle)
1772 
1773  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  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  CALL timeset(routinen, handle)
1797 
1798  NULLIFY (logger, print_key)
1799  logger => cp_get_default_logger()
1800 
1801  print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1802  "MOTION%PINT%HELIUM%PRINT%FORCES")
1803  should_output = btest(cp_print_key_should_output( &
1804  logger%iter_info, &
1805  basis_section=print_key), cp_p_file)
1806 
1807  IF (.NOT. should_output) THEN
1808  CALL timestop(handle)
1809  RETURN
1810  END IF
1811 
1812  ! check if there is anything to be printed out
1813  IF (.NOT. helium_env(1)%helium%solute_present) THEN
1814  msgstr = "Warning: force printout requested but there is no solute!"
1815  CALL helium_write_line(msgstr)
1816  CALL timestop(handle)
1817  RETURN
1818  END IF
1819 
1820  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  extension=".dat")
1827 
1828  ! print all force components in one line
1829  DO ib = 1, helium_env(1)%helium%solute_beads
1830  idim = 0
1831  DO ia = 1, helium_env(1)%helium%solute_atoms
1832  DO ic = 1, 3
1833  idim = idim + 1
1834  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  WRITE (unit_nr, *)
1839 
1840  CALL m_flush(unit_nr)
1841  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
1842 
1843  END IF
1844 
1845  CALL timestop(handle)
1846 
1847  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
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition: cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Definition: cp_parser_methods.F:421
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition: cp_parser_types.F:109
cp_parser_types::parser_create
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Definition: cp_parser_types.F:147
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition: cp_units.F:1179
cp_units::cp_unit_to_cp2k
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition: cp_units.F:1150
helium_common
Independent helium subroutines shared with other modules.
Definition: helium_common.F:14
helium_common::helium_cycle_number
integer function, public helium_cycle_number(helium, atom_number, permutation)
Given the atom number and permutation state return the cycle number the atom belongs to.
Definition: helium_common.F:1485
helium_common::helium_pbc
subroutine, public helium_pbc(helium, r, enforce)
General PBC routine for helium.
Definition: helium_common.F:81
helium_common::helium_cycle_of
pure integer function, dimension(:), pointer, public helium_cycle_of(element, permutation)
Given an element of a permutation return the cycle it belongs to.
Definition: helium_common.F:1633
helium_common::helium_path_length
pure integer function, public helium_path_length(helium, atom_number, permutation)
Given the atom number and permutation state return the length of the path this atom belongs to.
Definition: helium_common.F:1595
helium_common::helium_is_winding
logical function, public helium_is_winding(helium, atmidx, pos, permutation)
Given the atom index and permutation state returns .TRUE. if the atom belongs to a winding path,...
Definition: helium_common.F:2496
helium_interactions
Methods that handle helium-solvent and helium-helium interactions.
Definition: helium_interactions.F:13
helium_interactions::helium_total_pair_action
real(kind=dp) function, public helium_total_pair_action(helium)
Computes the total pair action of the helium.
Definition: helium_interactions.F:215
helium_interactions::helium_total_inter_action
real(kind=dp) function, public helium_total_inter_action(pint_env, helium)
Computes the total interaction of the helium with the solute.
Definition: helium_interactions.F:296
helium_interactions::helium_total_link_action
real(kind=dp) function, public helium_total_link_action(helium)
Computes the total harmonic link action of the helium.
Definition: helium_interactions.F:178
helium_io
I/O subroutines for helium.
Definition: helium_io.F:13
helium_io::helium_print_vector
subroutine, public helium_print_vector(helium_env, pkey, DATA, uconv, col_label, cmmnt, fname, fpos, avg)
Print a 3D real vector according to printkey <pkey>
Definition: helium_io.F:585
helium_io::helium_print_action
subroutine, public helium_print_action(pint_env, helium_env)
Print helium action file to HELIUMPRINTACTION.
Definition: helium_io.F:988
helium_io::helium_print_rdf
subroutine, public helium_print_rdf(helium_env)
Print radial distribution functions according to HELIUMPRINTRDF.
Definition: helium_io.F:1403
helium_io::helium_print_force
subroutine, public helium_print_force(helium_env)
Print helium force according to HELIUMPRINTFORCE.
Definition: helium_io.F:1785
helium_io::helium_print_accepts
subroutine, public helium_print_accepts(helium_env)
Print acceptance counts according to HELIUMPRINTACCEPTS.
Definition: helium_io.F:726
helium_io::helium_print_energy
subroutine, public helium_print_energy(helium_env)
Print energies according to HELIUMPRINTENERGY.
Definition: helium_io.F:475
helium_io::helium_print_perm
subroutine, public helium_print_perm(helium_env)
Print permutation state according to HELIUMPRINTPERM.
Definition: helium_io.F:795
helium_io::helium_read_xyz
subroutine, public helium_read_xyz(coords, file_name, para_env)
Read XYZ coordinates from file.
Definition: helium_io.F:102
helium_io::helium_print_rho
subroutine, public helium_print_rho(helium_env)
Print densities according to HELIUMPRINTRHO.
Definition: helium_io.F:1534
helium_io::helium_print_plength
subroutine, public helium_print_plength(helium_env)
Print permutation length according to HELIUMPRINTPLENGTH.
Definition: helium_io.F:1726
helium_io::helium_write_setup
subroutine, public helium_write_setup(helium)
Write helium parameters to the output unit.
Definition: helium_io.F:187
helium_io::helium_write_line
subroutine, public helium_write_line(line)
Writes out a line of text to the default output unit.
Definition: helium_io.F:447
helium_io::helium_write_cubefile
subroutine, public helium_write_cubefile(unit, comment, origin, deltar, ndim, DATA)
Write volumetric data to an orthorhombic cubefile.
Definition: helium_io.F:1661
helium_io::helium_print_coordinates
subroutine, public helium_print_coordinates(helium_env)
Print coordinates according to HELIUMPRINTCOORDINATES.
Definition: helium_io.F:1107
helium_types
Data types representing superfluid helium.
Definition: helium_types.F:15
helium_types::e_id_potential
integer, parameter, public e_id_potential
Definition: helium_types.F:38
helium_types::e_id_thermo
integer, parameter, public e_id_thermo
Definition: helium_types.F:38
helium_types::e_id_virial
integer, parameter, public e_id_virial
Definition: helium_types.F:38
helium_types::e_id_interact
integer, parameter, public e_id_interact
Definition: helium_types.F:38
helium_types::e_id_kinetic
integer, parameter, public e_id_kinetic
Definition: helium_types.F:38
helium_types::e_id_total
integer, parameter, public e_id_total
Energy contributions - symbolic names for indexing energy arrays.
Definition: helium_types.F:38
helium_types::rho_num
integer, parameter, public rho_num
number of density function identifiers
Definition: helium_types.F:50
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::helium_solute_intpot_mwater
integer, parameter, public helium_solute_intpot_mwater
Definition: input_constants.F:322
input_constants::helium_sampling_ceperley
integer, parameter, public helium_sampling_ceperley
Definition: input_constants.F:347
input_constants::helium_cell_shape_octahedron
integer, parameter, public helium_cell_shape_octahedron
Definition: input_constants.F:317
input_constants::helium_solute_intpot_none
integer, parameter, public helium_solute_intpot_none
Definition: input_constants.F:322
input_constants::perm_cycle
integer, parameter, public perm_cycle
Definition: input_constants.F:358
input_constants::perm_plain
integer, parameter, public perm_plain
Definition: input_constants.F:358
input_constants::helium_sampling_worm
integer, parameter, public helium_sampling_worm
Definition: input_constants.F:347
input_constants::helium_cell_shape_cube
integer, parameter, public helium_cell_shape_cube
Definition: input_constants.F:317
input_constants::fmt_id_pdb
integer, parameter, public fmt_id_pdb
Definition: input_constants.F:381
input_constants::helium_solute_intpot_nnp
integer, parameter, public helium_solute_intpot_nnp
Definition: input_constants.F:322
input_constants::fmt_id_xyz
integer, parameter, public fmt_id_xyz
Definition: input_constants.F:381
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition: kinds.F:54
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
physcon::massunit
real(kind=dp), parameter, public massunit
Definition: physcon.F:141
pint_types
Definition: pint_types.F:8