(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
13MODULE 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
90CONTAINS
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
1944END MODULE helium_io
cp_parser_methods::parser_get_object
Definition cp_parser_methods.F:54
memory_utilities::reallocate
Definition memory_utilities.F:27
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_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_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_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_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_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_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
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_parser_types::cp_parser_type
represent a parser
Definition cp_parser_types.F:77
helium_types::helium_solvent_p_type
data structure for array of solvent helium environments
Definition helium_types.F:300
helium_types::helium_solvent_type
data structure for solvent helium
Definition helium_types.F:93
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pint_types::pint_env_type
environment for a path integral run
Definition pint_types.F:112