(git:34ef472)
helium_methods.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 Methods dealing with helium_solvent_type
10 !> \author Lukasz Walewski
11 !> \date 2009-06-10
12 ! **************************************************************************************************
13 MODULE helium_methods
14 
15  USE atomic_kind_types, ONLY: get_atomic_kind
16  USE bibliography, ONLY: walewski2014,&
17  cite_reference
18  USE cell_types, ONLY: cell_type,&
19  get_cell
20  USE cp_files, ONLY: close_file,&
21  open_file
22  USE cp_log_handling, ONLY: cp_add_default_logger,&
23  cp_get_default_logger,&
24  cp_logger_create,&
25  cp_logger_release,&
26  cp_logger_type,&
27  cp_rm_default_logger
28  USE cp_output_handling, ONLY: cp_printkey_is_on
29  USE cp_subsys_types, ONLY: cp_subsys_get,&
30  cp_subsys_type
31  USE f77_interface, ONLY: f_env_add_defaults,&
32  f_env_rm_defaults,&
33  f_env_type
34  USE force_env_types, ONLY: force_env_get
35  USE helium_common, ONLY: helium_com,&
36  helium_pbc
37  USE helium_interactions, ONLY: helium_vij
38  USE helium_io, ONLY: helium_write_line,&
39  helium_write_setup
40  USE helium_nnp, ONLY: helium_init_nnp
41  USE helium_sampling, ONLY: helium_sample
42  USE helium_types, ONLY: helium_solvent_p_type,&
43  helium_solvent_type,&
44  rho_atom_number,&
45  rho_moment_of_inertia,&
46  rho_num,&
47  rho_projected_area,&
48  rho_winding_cycle,&
49  rho_winding_number
50  USE input_constants, ONLY: helium_cell_shape_cube,&
51  helium_cell_shape_octahedron,&
52  helium_sampling_ceperley,&
53  helium_sampling_worm,&
54  helium_solute_intpot_nnp,&
55  helium_solute_intpot_none
56  USE input_section_types, ONLY: section_vals_get,&
57  section_vals_get_subs_vals,&
58  section_vals_type,&
59  section_vals_val_get,&
60  section_vals_val_set
61  USE kinds, ONLY: default_path_length,&
62  default_string_length,&
63  dp,&
64  max_line_length
65  USE mathconstants, ONLY: pi,&
66  twopi
67  USE message_passing, ONLY: mp_para_env_type
68  USE nnp_environment_types, ONLY: nnp_env_release
69  USE parallel_rng_types, ONLY: gaussian,&
70  uniform,&
71  rng_stream_p_type,&
72  rng_stream_type
73  USE particle_list_types, ONLY: particle_list_type
74  USE physcon, ONLY: a_mass,&
75  angstrom,&
76  boltzmann,&
77  h_bar,&
78  kelvin,&
79  massunit
80  USE pint_public, ONLY: pint_com_pos
81  USE pint_types, ONLY: pint_env_type
82  USE splines_methods, ONLY: init_spline,&
83  init_splinexy
84  USE splines_types, ONLY: spline_data_create,&
85  spline_data_release
86 #include "../base/base_uses.f90"
87 
88  IMPLICIT NONE
89 
90  PRIVATE
91 
92  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .true.
93  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_methods'
94 
95  PUBLIC :: helium_create
96  PUBLIC :: helium_init
97  PUBLIC :: helium_release
98 
99 CONTAINS
100 
101 ! ***************************************************************************
102 !> \brief Data-structure that holds all needed information about
103 !> (superfluid) helium solvent
104 !> \param helium_env ...
105 !> \param input ...
106 !> \param solute ...
107 !> \par History
108 !> 2016-07-14 Modified to work with independent helium_env [cschran]
109 !> 2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
110 !> \author hforbert
111 ! **************************************************************************************************
112  SUBROUTINE helium_create(helium_env, input, solute)
113  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
114  TYPE(section_vals_type), POINTER :: input
115  TYPE(pint_env_type), INTENT(IN), OPTIONAL :: solute
116 
117  CHARACTER(len=*), PARAMETER :: routinen = 'helium_create'
118 
119  CHARACTER(len=default_path_length) :: msg_str, potential_file_name
120  INTEGER :: color_sub, handle, i, input_unit, isize, &
121  itmp, j, k, mepos, nlines, ntab, &
122  num_env, pdx
123  INTEGER, DIMENSION(:), POINTER :: env_all
124  LOGICAL :: expl_cell, expl_dens, expl_nats, &
125  expl_pot, explicit, ltmp
126  REAL(kind=dp) :: cgeof, dx, he_mass, mhe, rtmp, t, tau, &
127  tcheck, x1, x_spline
128  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pot_transfer
129  TYPE(cp_logger_type), POINTER :: logger, tmplogger
130  TYPE(mp_para_env_type), POINTER :: new_comm
131  TYPE(section_vals_type), POINTER :: helium_section, input_worm, nnp_section
132 
133  CALL timeset(routinen, handle)
134 
135  CALL cite_reference(walewski2014)
136  NULLIFY (logger)
137  logger => cp_get_default_logger()
138 
139  NULLIFY (helium_section)
140  helium_section => section_vals_get_subs_vals(input, &
141  "MOTION%PINT%HELIUM")
142  CALL section_vals_get(helium_section, explicit=explicit)
143  cpassert(explicit)
144 
145  ! get number of environments
146  CALL section_vals_val_get(helium_section, "NUM_ENV", &
147  explicit=explicit)
148  IF (explicit) THEN
149  CALL section_vals_val_get(helium_section, "NUM_ENV", &
150  i_val=num_env)
151  ELSE
152  num_env = logger%para_env%num_pe
153  END IF
154  cpassert(num_env >= 0)
155  IF (num_env .NE. logger%para_env%num_pe) THEN
156  msg_str = "NUM_ENV not equal to number of processors"
157  cpwarn(msg_str)
158  END IF
159  ! calculate number of tasks for each processor
160  mepos = num_env/logger%para_env%num_pe &
161  + min(mod(num_env, logger%para_env%num_pe)/(logger%para_env%mepos + 1), 1)
162  ! gather result
163  NULLIFY (env_all)
164  ALLOCATE (env_all(logger%para_env%num_pe))
165  env_all(:) = 0
166  CALL logger%para_env%allgather(mepos, env_all)
167 
168  ! create new communicator for processors with helium_env
169  IF (mepos == 0) THEN
170  color_sub = 0
171  ELSE
172  color_sub = 1
173  END IF
174  ALLOCATE (new_comm)
175  CALL new_comm%from_split(logger%para_env, color_sub)
176  ! release new_comm for processors without helium_env
177  IF (mepos == 0) THEN
178  CALL new_comm%free()
179  DEALLOCATE (new_comm)
180  NULLIFY (new_comm)
181  END IF
182 
183  NULLIFY (helium_env)
184  IF (mepos .GT. 0) THEN
185  ALLOCATE (helium_env(mepos))
186  DO k = 1, mepos
187  helium_env(k)%comm => new_comm
188  NULLIFY (helium_env(k)%env_all)
189  helium_env(k)%env_all => env_all
190  ALLOCATE (helium_env(k)%helium)
191  NULLIFY (helium_env(k)%helium%input)
192  helium_env(k)%helium%input => input
193  helium_env(k)%helium%num_env = num_env
194  END DO
195  ! RNG state create & init
196  CALL helium_rng_init(helium_env)
197 
198  DO k = 1, mepos
199  NULLIFY (helium_env(k)%helium%ptable, &
200  helium_env(k)%helium%permutation, &
201  helium_env(k)%helium%savepermutation, &
202  helium_env(k)%helium%iperm, &
203  helium_env(k)%helium%saveiperm, &
204  helium_env(k)%helium%itmp_atoms_1d, &
205  helium_env(k)%helium%ltmp_atoms_1d, &
206  helium_env(k)%helium%itmp_atoms_np_1d, &
207  helium_env(k)%helium%pos, &
208  helium_env(k)%helium%savepos, &
209  helium_env(k)%helium%work, &
210  helium_env(k)%helium%force_avrg, &
211  helium_env(k)%helium%force_inst, &
212  helium_env(k)%helium%rtmp_3_np_1d, &
213  helium_env(k)%helium%rtmp_p_ndim_1d, &
214  helium_env(k)%helium%rtmp_p_ndim_np_1d, &
215  helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
216  helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
217  helium_env(k)%helium%rtmp_p_ndim_2d, &
218  helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
219  helium_env(k)%helium%tmatrix, helium_env(k)%helium%pmatrix, &
220  helium_env(k)%helium%nmatrix, helium_env(k)%helium%ipmatrix, &
221  helium_env(k)%helium%u0, helium_env(k)%helium%e0, &
222  helium_env(k)%helium%uoffdiag, helium_env(k)%helium%eoffdiag, &
223  helium_env(k)%helium%vij, &
224  helium_env(k)%helium%rdf_inst, &
225  helium_env(k)%helium%plength_avrg, &
226  helium_env(k)%helium%plength_inst, &
227  helium_env(k)%helium%atom_plength, &
228  helium_env(k)%helium%ename &
229  )
230 
231  helium_env(k)%helium%accepts = 0
232  helium_env(k)%helium%relrot = 0
233 
234  ! check if solute is present in our simulation
235  helium_env(k)%helium%solute_present = .false.
236  helium_env(k)%helium%solute_atoms = 0
237  helium_env(k)%helium%solute_beads = 0
238  CALL section_vals_val_get( &
239  helium_section, &
240  "HELIUM_ONLY", &
241  l_val=ltmp)
242  IF (.NOT. ltmp) THEN
243  IF (PRESENT(solute)) THEN
244  helium_env(k)%helium%solute_present = .true.
245  helium_env(k)%helium%solute_atoms = solute%ndim/3
246  helium_env(k)%helium%solute_beads = solute%p
247  END IF
248  END IF
249 
250  CALL section_vals_val_get(helium_section, "NBEADS", &
251  i_val=helium_env(k)%helium%beads)
252  CALL section_vals_val_get(helium_section, "INOROT", &
253  i_val=helium_env(k)%helium%iter_norot)
254  CALL section_vals_val_get(helium_section, "IROT", &
255  i_val=helium_env(k)%helium%iter_rot)
256 
257  ! get number of steps and current step number from PINT
258  CALL section_vals_val_get(input, "MOTION%PINT%ITERATION", &
259  i_val=itmp)
260  helium_env(k)%helium%first_step = itmp
261  CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
262  explicit=explicit)
263  IF (explicit) THEN
264  CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
265  i_val=itmp)
266  helium_env(k)%helium%last_step = itmp
267  helium_env(k)%helium%num_steps = helium_env(k)%helium%last_step &
268  - helium_env(k)%helium%first_step
269  ELSE
270  CALL section_vals_val_get(input, "MOTION%PINT%NUM_STEPS", &
271  i_val=itmp)
272  helium_env(k)%helium%num_steps = itmp
273  helium_env(k)%helium%last_step = helium_env(k)%helium%first_step &
274  + helium_env(k)%helium%num_steps
275  END IF
276 
277  ! boundary conditions
278  CALL section_vals_val_get(helium_section, "PERIODIC", &
279  l_val=helium_env(k)%helium%periodic)
280  CALL section_vals_val_get(helium_section, "CELL_SHAPE", &
281  i_val=helium_env(k)%helium%cell_shape)
282 
283  CALL section_vals_val_get(helium_section, "DROPLET_RADIUS", &
284  r_val=helium_env(k)%helium%droplet_radius)
285 
286  ! Set density Rho, number of atoms N and volume V ( Rho = N / V ).
287  ! Allow only 2 out of 3 values to be defined at the same time, calculate
288  ! the third.
289  ! Note, that DENSITY and NATOMS keywords have default values, while
290  ! CELL_SIZE does not. Thus if CELL_SIZE is given explicitly then one and
291  ! only one of the two remaining options must be give explicitly as well.
292  ! If CELL_SIZE is not given explicitly then all four combinations of the
293  ! two other options are valid.
294  CALL section_vals_val_get(helium_section, "DENSITY", &
295  explicit=expl_dens, r_val=helium_env(k)%helium%density)
296  CALL section_vals_val_get(helium_section, "NATOMS", &
297  explicit=expl_nats, i_val=helium_env(k)%helium%atoms)
298  CALL section_vals_val_get(helium_section, "CELL_SIZE", &
299  explicit=expl_cell)
300  cgeof = 1.0_dp
301  IF (helium_env(k)%helium%periodic) THEN
302  IF (helium_env(k)%helium%cell_shape .EQ. helium_cell_shape_octahedron) cgeof = 2.0_dp
303  END IF
304  rtmp = (cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%density)**(1.0_dp/3.0_dp)
305  IF (.NOT. expl_cell) THEN
306  helium_env(k)%helium%cell_size = rtmp
307  ELSE
308  CALL section_vals_val_get(helium_section, "CELL_SIZE", &
309  r_val=helium_env(k)%helium%cell_size)
310  ! only more work if not all three values are consistent:
311  IF (abs(helium_env(k)%helium%cell_size - rtmp) .GT. 100.0_dp*epsilon(0.0_dp)* &
312  (abs(helium_env(k)%helium%cell_size) + rtmp)) THEN
313  IF (expl_dens .AND. expl_nats) THEN
314  msg_str = "DENSITY, NATOMS and CELL_SIZE options "// &
315  "contradict each other"
316  cpwarn(msg_str)
317  END IF
318  !ok we have enough freedom to resolve the conflict:
319  IF (.NOT. expl_dens) THEN
320  helium_env(k)%helium%density = cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%cell_size**3.0_dp
321  IF (.NOT. expl_nats) THEN
322  msg_str = "CELL_SIZE defined but neither "// &
323  "NATOMS nor DENSITY given, using default NATOMS."
324  cpwarn(msg_str)
325  END IF
326  ELSE ! ( expl_dens .AND. .NOT. expl_nats )
327  ! calculate the nearest number of atoms for given conditions
328  helium_env(k)%helium%atoms = nint(helium_env(k)%helium%density* &
329  helium_env(k)%helium%cell_size**3.0_dp/cgeof)
330  ! adjust cell size to maintain correct density
331  ! (should be a small correction)
332  rtmp = (cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%density &
333  )**(1.0_dp/3.0_dp)
334  IF (abs(helium_env(k)%helium%cell_size - rtmp) .GT. 100.0_dp*epsilon(0.0_dp) &
335  *(abs(helium_env(k)%helium%cell_size) + rtmp)) THEN
336  msg_str = "Adjusting actual cell size "// &
337  "to maintain correct density."
338  cpwarn(msg_str)
339  helium_env(k)%helium%cell_size = rtmp
340  END IF
341  END IF
342  END IF
343  END IF
344  helium_env(k)%helium%cell_size_inv = 1.0_dp/helium_env(k)%helium%cell_size
345  ! From now on helium%density, helium%atoms and helium%cell_size are
346  ! correctly defined.
347 
348  ! set the M matrix for winding number calculations
349  SELECT CASE (helium_env(k)%helium%cell_shape)
350 
351  CASE (helium_cell_shape_octahedron)
352  helium_env(k)%helium%cell_m(1, 1) = -0.5_dp*helium_env(k)%helium%cell_size
353  helium_env(k)%helium%cell_m(2, 1) = 0.5_dp*helium_env(k)%helium%cell_size
354  helium_env(k)%helium%cell_m(3, 1) = 0.5_dp*helium_env(k)%helium%cell_size
355  helium_env(k)%helium%cell_m(1, 2) = 0.5_dp*helium_env(k)%helium%cell_size
356  helium_env(k)%helium%cell_m(2, 2) = -0.5_dp*helium_env(k)%helium%cell_size
357  helium_env(k)%helium%cell_m(3, 2) = 0.5_dp*helium_env(k)%helium%cell_size
358  helium_env(k)%helium%cell_m(1, 3) = 0.5_dp*helium_env(k)%helium%cell_size
359  helium_env(k)%helium%cell_m(2, 3) = 0.5_dp*helium_env(k)%helium%cell_size
360  helium_env(k)%helium%cell_m(3, 3) = -0.5_dp*helium_env(k)%helium%cell_size
361  helium_env(k)%helium%cell_m_inv(1, 1) = 0.0_dp
362  helium_env(k)%helium%cell_m_inv(2, 1) = 1.0_dp/helium_env(k)%helium%cell_size
363  helium_env(k)%helium%cell_m_inv(3, 1) = 1.0_dp/helium_env(k)%helium%cell_size
364  helium_env(k)%helium%cell_m_inv(1, 2) = 1.0_dp/helium_env(k)%helium%cell_size
365  helium_env(k)%helium%cell_m_inv(2, 2) = 0.0_dp
366  helium_env(k)%helium%cell_m_inv(3, 2) = 1.0_dp/helium_env(k)%helium%cell_size
367  helium_env(k)%helium%cell_m_inv(1, 3) = 1.0_dp/helium_env(k)%helium%cell_size
368  helium_env(k)%helium%cell_m_inv(2, 3) = 1.0_dp/helium_env(k)%helium%cell_size
369  helium_env(k)%helium%cell_m_inv(3, 3) = 0.0_dp
370  CASE (helium_cell_shape_cube)
371 
372  helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size
373  helium_env(k)%helium%cell_m(2, 1) = 0.0_dp
374  helium_env(k)%helium%cell_m(3, 1) = 0.0_dp
375  helium_env(k)%helium%cell_m(1, 2) = 0.0_dp
376  helium_env(k)%helium%cell_m(2, 2) = helium_env(k)%helium%cell_size
377  helium_env(k)%helium%cell_m(3, 2) = 0.0_dp
378  helium_env(k)%helium%cell_m(1, 3) = 0.0_dp
379  helium_env(k)%helium%cell_m(2, 3) = 0.0_dp
380  helium_env(k)%helium%cell_m(3, 3) = helium_env(k)%helium%cell_size
381  helium_env(k)%helium%cell_m_inv(1, 1) = 1.0_dp/helium_env(k)%helium%cell_size
382  helium_env(k)%helium%cell_m_inv(2, 1) = 0.0_dp
383  helium_env(k)%helium%cell_m_inv(3, 1) = 0.0_dp
384  helium_env(k)%helium%cell_m_inv(1, 2) = 0.0_dp
385  helium_env(k)%helium%cell_m_inv(2, 2) = 1.0_dp/helium_env(k)%helium%cell_size
386  helium_env(k)%helium%cell_m_inv(3, 2) = 0.0_dp
387  helium_env(k)%helium%cell_m_inv(1, 3) = 0.0_dp
388  helium_env(k)%helium%cell_m_inv(2, 3) = 0.0_dp
389  helium_env(k)%helium%cell_m_inv(3, 3) = 1.0_dp/helium_env(k)%helium%cell_size
390  CASE DEFAULT
391  helium_env(k)%helium%cell_m(:, :) = 0.0_dp
392  helium_env(k)%helium%cell_m_inv(:, :) = 0.0_dp
393 
394  END SELECT
395 
396  END DO ! k
397 
398  IF (logger%para_env%is_source()) THEN
399  CALL section_vals_val_get(helium_section, "POTENTIAL_FILE_NAME", &
400  c_val=potential_file_name)
401  CALL open_file(file_name=trim(potential_file_name), &
402  file_action="READ", file_status="OLD", unit_number=input_unit)
403  READ (input_unit, "(A)") msg_str
404  READ (msg_str, *, iostat=i) nlines, pdx, tau, &
405  x_spline, dx, he_mass
406  IF (i /= 0) THEN
407  ! old style potential file, use default mass and potential
408  he_mass = 4.00263037059764_dp !< 4He mass in [u]
409  expl_pot = .false.
410  READ (msg_str, *, iostat=i) nlines, pdx, tau, &
411  x_spline, dx
412  IF (i /= 0) THEN
413  msg_str = "Format/Read Error from Solvent POTENTIAL_FILE"
414  cpabort(msg_str)
415  END IF
416  ELSE
417  expl_pot = .true.
418  ! in file really hb2m in kelvin angstrom**2
419  he_mass = angstrom**2*kelvin/massunit/he_mass
420  ! tau might be negative to get older versions of cp2k,
421  ! that cannot handle the new potential file format to
422  ! crash and not run a calculation with wrong mass/potential
423  tau = abs(tau)
424  END IF
425  tau = kelvin/tau
426  x_spline = x_spline/angstrom
427  dx = dx/angstrom
428  END IF
429 
430  CALL helium_env(1)%comm%bcast(nlines, logger%para_env%source)
431  CALL helium_env(1)%comm%bcast(pdx, logger%para_env%source)
432  CALL helium_env(1)%comm%bcast(tau, logger%para_env%source)
433  CALL helium_env(1)%comm%bcast(x_spline, logger%para_env%source)
434  CALL helium_env(1)%comm%bcast(dx, logger%para_env%source)
435  CALL helium_env(1)%comm%bcast(he_mass, logger%para_env%source)
436  isize = (pdx + 1)*(pdx + 2) + 1
437  ALLOCATE (pot_transfer(nlines, isize))
438  IF (logger%para_env%is_source()) THEN
439  IF (expl_pot) THEN
440  DO i = 1, nlines
441  READ (input_unit, *) pot_transfer(i, :)
442  END DO
443  ELSE
444  DO i = 1, nlines
445  READ (input_unit, *) pot_transfer(i, 1:isize - 1)
446  ! potential implicit, calculate it here now:
447  pot_transfer(i, isize) = helium_vij(x_spline + (i - 1)*dx)*kelvin
448  END DO
449  END IF
450  CALL close_file(unit_number=input_unit)
451  END IF
452  CALL helium_env(1)%comm%bcast(pot_transfer, logger%para_env%source)
453 
454  CALL spline_data_create(helium_env(1)%helium%vij)
455  CALL init_splinexy(helium_env(1)%helium%vij, nlines)
456  helium_env(1)%helium%vij%x1 = x_spline
457 
458  CALL spline_data_create(helium_env(1)%helium%u0)
459  CALL init_splinexy(helium_env(1)%helium%u0, nlines)
460  helium_env(1)%helium%u0%x1 = x_spline
461 
462  CALL spline_data_create(helium_env(1)%helium%e0)
463  CALL init_splinexy(helium_env(1)%helium%e0, nlines)
464  helium_env(1)%helium%e0%x1 = x_spline
465 
466  isize = pdx + 1
467  ntab = ((isize + 1)*isize)/2 - 1 ! -1 because endpoint approx treated separately
468  ALLOCATE (helium_env(1)%helium%uoffdiag(ntab, 2, nlines))
469  ALLOCATE (helium_env(1)%helium%eoffdiag(ntab, 2, nlines))
470  DO j = 1, isize
471  DO i = j, isize
472  IF (i + j == 2) cycle ! endpoint approx later separately
473  k = ((i - 1)*i)/2 + j
474  helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2)
475  CALL init_spline(helium_env(1)%helium%vij, dx=dx)
476  helium_env(1)%helium%uoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:)
477  helium_env(1)%helium%uoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:)
478  k = k + ((isize + 1)*isize)/2
479  helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2)/kelvin
480  CALL init_spline(helium_env(1)%helium%vij, dx=dx)
481  helium_env(1)%helium%eoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:)
482  helium_env(1)%helium%eoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:)
483  ntab = ntab - 1
484  END DO
485  END DO
486 
487  ntab = SIZE(pot_transfer, 2)
488  helium_env(1)%helium%vij%y(:) = pot_transfer(:, ntab)/kelvin
489  CALL init_spline(helium_env(1)%helium%vij, dx=dx)
490 
491  helium_env(1)%helium%u0%y(:) = pot_transfer(:, 1)
492  CALL init_spline(helium_env(1)%helium%u0, dx=dx)
493  k = ((isize + 1)*isize)/2 + 1
494  helium_env(1)%helium%e0%y(:) = pot_transfer(:, k)/kelvin
495  CALL init_spline(helium_env(1)%helium%e0, dx=dx)
496 
497  DO k = 2, mepos
498  helium_env(k)%helium%vij => helium_env(1)%helium%vij
499  helium_env(k)%helium%u0 => helium_env(1)%helium%u0
500  helium_env(k)%helium%e0 => helium_env(1)%helium%e0
501  helium_env(k)%helium%uoffdiag => helium_env(1)%helium%uoffdiag
502  helium_env(k)%helium%eoffdiag => helium_env(1)%helium%eoffdiag
503  END DO
504 
505  DO k = 1, mepos
506 
507  helium_env(k)%helium%pdx = pdx
508  helium_env(k)%helium%tau = tau
509 
510  ! boltzmann : Boltzmann constant [J/K]
511  ! h_bar : Planck constant [J*s]
512  ! J = kg*m^2/s^2
513  ! 4He mass in [kg]
514  mhe = he_mass*a_mass
515  ! physical temperature [K]
516  t = kelvin/helium_env(k)%helium%tau/helium_env(k)%helium%beads
517  ! prefactors for calculating superfluid fractions [Angstrom^-2]
518  helium_env(k)%helium%wpref = (((1e-20/h_bar)*mhe)/h_bar)*boltzmann*t
519  helium_env(k)%helium%apref = (((4e-20/h_bar)*mhe)/h_bar)*boltzmann*t
520 
521  helium_env(k)%helium%he_mass_au = he_mass*massunit
522  helium_env(k)%helium%hb2m = 1.0_dp/helium_env(k)%helium%he_mass_au
523  helium_env(k)%helium%pweight = 0.0_dp
524 
525  ! Default in case sampling_method is not helium_sampling_worm.
526  helium_env(k)%helium%worm_max_open_cycles = 0
527 
528  ! Choose sampling method:
529  CALL section_vals_val_get(helium_section, "SAMPLING_METHOD", &
530  i_val=helium_env(k)%helium%sampling_method)
531 
532  SELECT CASE (helium_env(k)%helium%sampling_method)
533  CASE (helium_sampling_ceperley)
534  ! check value of maxcycle
535  CALL section_vals_val_get(helium_section, "CEPERLEY%MAX_PERM_CYCLE", &
536  i_val=helium_env(k)%helium%maxcycle)
537  i = helium_env(k)%helium%maxcycle
538  cpassert(i >= 0)
539  i = helium_env(k)%helium%atoms - helium_env(k)%helium%maxcycle
540  cpassert(i >= 0)
541 
542  ! set m-distribution parameters
543  CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%DISTRIBUTION-TYPE", &
544  i_val=i)
545  cpassert(i >= 1)
546  cpassert(i <= 6)
547  helium_env(k)%helium%m_dist_type = i
548  CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-VALUE", &
549  i_val=i)
550  cpassert(i >= 1)
551  cpassert(i <= helium_env(k)%helium%maxcycle)
552  helium_env(k)%helium%m_value = i
553  CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-RATIO", &
554  r_val=rtmp)
555  cpassert(rtmp > 0.0_dp)
556  cpassert(rtmp <= 1.0_dp)
557  helium_env(k)%helium%m_ratio = rtmp
558 
559  CALL section_vals_val_get(helium_section, "CEPERLEY%BISECTION", &
560  i_val=helium_env(k)%helium%bisection)
561  ! precheck bisection value (not all invalids are filtered out here yet)
562  i = helium_env(k)%helium%bisection
563  cpassert(i > 1)
564  i = helium_env(k)%helium%beads - helium_env(k)%helium%bisection
565  cpassert(i > 0)
566  !
567  itmp = helium_env(k)%helium%bisection
568  rtmp = 2.0_dp**(anint(log(real(itmp, dp))/log(2.0_dp)))
569  tcheck = abs(real(itmp, kind=dp) - rtmp)
570  IF (tcheck > 100.0_dp*epsilon(0.0_dp)) THEN
571  msg_str = "BISECTION should be integer power of 2."
572  cpabort(msg_str)
573  END IF
574  helium_env(k)%helium%bisctlog2 = nint(log(real(itmp, dp))/log(2.0_dp))
575 
576  CASE (helium_sampling_worm)
577  NULLIFY (input_worm)
578  input_worm => section_vals_get_subs_vals(helium_env(k)%helium%input, &
579  "MOTION%PINT%HELIUM%WORM")
580  CALL section_vals_val_get(helium_section, "WORM%CENTROID_DRMAX", &
581  r_val=helium_env(k)%helium%worm_centroid_drmax)
582 
583  CALL section_vals_val_get(helium_section, "WORM%STAGING_L", &
584  i_val=helium_env(k)%helium%worm_staging_l)
585 
586  CALL section_vals_val_get(helium_section, "WORM%OPEN_CLOSE_SCALE", &
587  r_val=helium_env(k)%helium%worm_open_close_scale)
588 
589  CALL section_vals_val_get(helium_section, "WORM%ALLOW_OPEN", &
590  l_val=helium_env(k)%helium%worm_allow_open)
591  IF (helium_env(k)%helium%atoms == 1) THEN
592  IF (helium_env(k)%helium%worm_allow_open) THEN
593  msg_str = "Default enabled open state sampling "// &
594  "for only 1 He might be inefficient."
595  cpwarn(msg_str)
596  END IF
597  END IF
598 
599  CALL section_vals_val_get(helium_section, "WORM%MAX_OPEN_CYCLES", &
600  i_val=helium_env(k)%helium%worm_max_open_cycles)
601 
602  IF (helium_env(k)%helium%worm_staging_l + 1 >= helium_env(k)%helium%beads) THEN
603  msg_str = "STAGING_L for worm sampling is too large"
604  cpabort(msg_str)
605  ELSE IF (helium_env(k)%helium%worm_staging_l < 1) THEN
606  msg_str = "STAGING_L must be positive integer"
607  cpabort(msg_str)
608  END IF
609 
610  CALL section_vals_val_get(helium_section, "WORM%SHOW_STATISTICS", &
611  l_val=helium_env(k)%helium%worm_show_statistics)
612 
613  ! precompute an expensive scaling for the open and close acceptance probability
614  ! tau is not included here, as It will be first defined in a few lines
615  rtmp = 2.0_dp*pi*helium_env(k)%helium%hb2m
616  rtmp = sqrt(rtmp)
617  rtmp = rtmp**3
618  rtmp = rtmp*helium_env(k)%helium%worm_open_close_scale
619  IF (helium_env(k)%helium%periodic) THEN
620  rtmp = rtmp*helium_env(k)%helium%density
621  ELSE
622  rtmp = rtmp*helium_env(k)%helium%atoms/ &
623  (4.0_dp/3.0_dp*pi*helium_env(k)%helium%droplet_radius**3)
624  END IF
625  helium_env(k)%helium%worm_ln_openclose_scale = log(rtmp)
626 
627  ! deal with acceptance statistics without changing the ceperley stuff
628  helium_env(k)%helium%maxcycle = 1
629  helium_env(k)%helium%bisctlog2 = 0
630 
631  ! get the absolute weights of the individual moves
632  helium_env(k)%helium%worm_all_limit = 0
633  CALL section_vals_val_get(helium_section, "WORM%CENTROID_WEIGHT", &
634  i_val=itmp)
635  helium_env(k)%helium%worm_centroid_min = 1
636  helium_env(k)%helium%worm_centroid_max = itmp
637  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
638  CALL section_vals_val_get(helium_section, "WORM%STAGING_WEIGHT", &
639  i_val=itmp)
640  helium_env(k)%helium%worm_staging_min = helium_env(k)%helium%worm_centroid_max + 1
641  helium_env(k)%helium%worm_staging_max = helium_env(k)%helium%worm_centroid_max + itmp
642  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
643  IF (helium_env(k)%helium%worm_allow_open) THEN
644  CALL section_vals_val_get(helium_section, "WORM%CRAWL_WEIGHT", &
645  i_val=itmp)
646  helium_env(k)%helium%worm_fcrawl_min = helium_env(k)%helium%worm_staging_max + 1
647  helium_env(k)%helium%worm_fcrawl_max = helium_env(k)%helium%worm_staging_max + itmp
648  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
649  helium_env(k)%helium%worm_bcrawl_min = helium_env(k)%helium%worm_fcrawl_max + 1
650  helium_env(k)%helium%worm_bcrawl_max = helium_env(k)%helium%worm_fcrawl_max + itmp
651  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
652  CALL section_vals_val_get(helium_section, "WORM%HEAD_TAIL_WEIGHT", &
653  i_val=itmp)
654  helium_env(k)%helium%worm_head_min = helium_env(k)%helium%worm_bcrawl_max + 1
655  helium_env(k)%helium%worm_head_max = helium_env(k)%helium%worm_bcrawl_max + itmp
656  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
657  helium_env(k)%helium%worm_tail_min = helium_env(k)%helium%worm_head_max + 1
658  helium_env(k)%helium%worm_tail_max = helium_env(k)%helium%worm_head_max + itmp
659  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
660  CALL section_vals_val_get(helium_section, "WORM%SWAP_WEIGHT", &
661  i_val=itmp)
662  helium_env(k)%helium%worm_swap_min = helium_env(k)%helium%worm_tail_max + 1
663  helium_env(k)%helium%worm_swap_max = helium_env(k)%helium%worm_tail_max + itmp
664  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
665  CALL section_vals_val_get(helium_section, "WORM%OPEN_CLOSE_WEIGHT", &
666  i_val=itmp)
667  helium_env(k)%helium%worm_open_close_min = helium_env(k)%helium%worm_swap_max + 1
668  helium_env(k)%helium%worm_open_close_max = helium_env(k)%helium%worm_swap_max + itmp
669  helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
670  CALL section_vals_val_get(helium_section, "WORM%CRAWL_REPETITION", &
671  i_val=helium_env(k)%helium%worm_repeat_crawl)
672  END IF
673 
674  !CPPostcondition(i<helium_env(k)%helium%beads,cp_failure_level,routineP,failure)
675  ! end of worm
676  CASE DEFAULT
677  msg_str = "Unknown helium sampling method!"
678  cpabort(msg_str)
679  END SELECT
680 
681  ! ALLOCATE helium-related arrays
682  i = helium_env(k)%helium%atoms
683  j = helium_env(k)%helium%beads
684  ALLOCATE (helium_env(k)%helium%pos(3, i, j))
685  helium_env(k)%helium%pos = 0.0_dp
686  ALLOCATE (helium_env(k)%helium%work(3, i, j))
687  ALLOCATE (helium_env(k)%helium%ptable(helium_env(k)%helium%maxcycle + 1))
688  ALLOCATE (helium_env(k)%helium%permutation(i))
689  ALLOCATE (helium_env(k)%helium%iperm(i))
690  ALLOCATE (helium_env(k)%helium%tmatrix(i, i))
691  ALLOCATE (helium_env(k)%helium%nmatrix(i, 2*i))
692  ALLOCATE (helium_env(k)%helium%pmatrix(i, i))
693  ALLOCATE (helium_env(k)%helium%ipmatrix(i, i))
694  itmp = helium_env(k)%helium%bisctlog2 + 2
695  ALLOCATE (helium_env(k)%helium%num_accepted(itmp, helium_env(k)%helium%maxcycle))
696  ALLOCATE (helium_env(k)%helium%plength_avrg(helium_env(k)%helium%atoms))
697  ALLOCATE (helium_env(k)%helium%plength_inst(helium_env(k)%helium%atoms))
698  ALLOCATE (helium_env(k)%helium%atom_plength(helium_env(k)%helium%atoms))
699  IF (helium_env(k)%helium%worm_max_open_cycles > 0) THEN
700  ALLOCATE (helium_env(k)%helium%savepermutation(i))
701  ALLOCATE (helium_env(k)%helium%saveiperm(i))
702  ALLOCATE (helium_env(k)%helium%savepos(3, i, j))
703  END IF
704 
705  ! check whether rdfs should be calculated and printed
706  helium_env(k)%helium%rdf_present = helium_property_active(helium_env(k)%helium, "RDF")
707  IF (helium_env(k)%helium%rdf_present) THEN
708  ! allocate & initialize rdf related data structures
709  CALL helium_rdf_init(helium_env(k)%helium)
710  END IF
711 
712  ! check whether densities should be calculated and printed
713  helium_env(k)%helium%rho_present = helium_property_active(helium_env(k)%helium, "RHO")
714  IF (helium_env(k)%helium%rho_present) THEN
715  ! allocate & initialize density related data structures
716  NULLIFY (helium_env(k)%helium%rho_property)
717  CALL helium_rho_init(helium_env(k)%helium)
718  END IF
719 
720  END DO
721 
722  ! restore averages calculated in previous runs
723  CALL helium_averages_restore(helium_env)
724 
725  DO k = 1, mepos
726  ! fill in the solute-related data structures
727  helium_env(k)%helium%e_corr = 0.0_dp
728  IF (helium_env(k)%helium%solute_present) THEN
729  IF (helium_env(k)%helium%solute_beads > helium_env(k)%helium%beads) THEN
730  ! Imaginary time striding for solute:
731  helium_env(k)%helium%bead_ratio = helium_env(k)%helium%solute_beads/ &
732  helium_env(k)%helium%beads
733  ! check if bead numbers are commensurate:
734  i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%beads - helium_env(k)%helium%solute_beads
735  IF (i /= 0) THEN
736  msg_str = "Adjust number of solute beads to multiple of solvent beads."
737  cpabort(msg_str)
738  END IF
739  msg_str = "Using multiple-time stepping in imaginary time for solute to couple "// &
740  "to fewer solvent beads, e.g. for factor 3: "// &
741  "he_1 - 3*sol_1; he_2 - 3*sol_4... "// &
742  "Avoid too large coupling factors."
743  cpwarn(msg_str)
744  ELSE IF (helium_env(k)%helium%solute_beads < helium_env(k)%helium%beads) THEN
745  ! Imaginary time striding for solvent:
746  helium_env(k)%helium%bead_ratio = helium_env(k)%helium%beads/ &
747  helium_env(k)%helium%solute_beads
748  ! check if bead numbers are commensurate:
749  i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%solute_beads - helium_env(k)%helium%beads
750  IF (i /= 0) THEN
751  msg_str = "Adjust number of solvent beads to multiple of solute beads."
752  cpabort(msg_str)
753  END IF
754  msg_str = "Coupling solvent beads to fewer solute beads via "// &
755  "direct coupling, e.g. for factor 3: "// &
756  "sol_1 - he_1,2,3; sol_2 - he_4,5,6..."
757  cpwarn(msg_str)
758  END IF
759 !TODO Adjust helium bead number if not comm. and if coords not given expl.
760 
761  ! check if tau, temperature and bead number are consistent:
762  tcheck = abs((helium_env(k)%helium%tau*helium_env(k)%helium%beads - solute%beta)/solute%beta)
763  IF (tcheck > 1.0e-14_dp) THEN
764  msg_str = "Tau, temperature and bead number are inconsistent."
765  cpabort(msg_str)
766  END IF
767 
768  CALL helium_set_solute_indices(helium_env(k)%helium, solute)
769  CALL helium_set_solute_cell(helium_env(k)%helium, solute)
770 
771  ! set the interaction potential type
772  CALL section_vals_val_get(helium_section, "SOLUTE_INTERACTION", &
773  i_val=helium_env(k)%helium%solute_interaction)
774  IF (helium_env(k)%helium%solute_interaction .EQ. helium_solute_intpot_nnp) THEN
775  IF (k == 1) THEN
776  NULLIFY (nnp_section)
777  nnp_section => section_vals_get_subs_vals(helium_section, "NNP")
778  CALL section_vals_get(nnp_section, explicit=explicit)
779  msg_str = "NNP section not explicitly stated. Using default file names."
780  IF (.NOT. explicit) cpwarn(msg_str)
781  END IF
782  ALLOCATE (helium_env(k)%helium%nnp)
783  CALL cp_logger_create(tmplogger, para_env=helium_env(k)%comm, template_logger=logger)
784  CALL cp_add_default_logger(tmplogger)
785  CALL helium_init_nnp(helium_env(k)%helium, helium_env(k)%helium%nnp, nnp_section)
786  CALL cp_rm_default_logger()
787  CALL cp_logger_release(tmplogger)
788  END IF
789  IF (helium_env(k)%helium%solute_interaction .EQ. helium_solute_intpot_none) THEN
790  WRITE (msg_str, '(A,I0,A)') &
791  "Solute found but no helium-solute interaction selected "// &
792  "(see SOLUTE_INTERACTION keyword)"
793  cpabort(msg_str)
794  END IF
795 
796  ! ALLOCATE solute-related arrays
797  ALLOCATE (helium_env(k)%helium%force_avrg(helium_env(k)%helium%solute_beads, &
798  helium_env(k)%helium%solute_atoms*3))
799  ALLOCATE (helium_env(k)%helium%force_inst(helium_env(k)%helium%solute_beads, &
800  helium_env(k)%helium%solute_atoms*3))
801 
802  ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_1d(solute%p*solute%ndim))
803  ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_np_1d(solute%p*solute%ndim*helium_env(k)%helium%num_env))
804  ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_2d(solute%p, solute%ndim))
805 
806  ELSE
807  helium_env(k)%helium%bead_ratio = 0
808  IF (helium_env(k)%helium%periodic) THEN
809  ! this assumes a specific potential (and its ugly):
810  x1 = angstrom*0.5_dp*helium_env(k)%helium%cell_size
811  ! 10.8 is in Kelvin, x1 needs to be in Angstrom,
812  ! since 2.9673 is in Angstrom
813  helium_env(k)%helium%e_corr = (twopi*helium_env(k)%helium%density/angstrom**3*10.8_dp* &
814  (544850.4_dp*exp(-13.353384_dp*x1/2.9673_dp)*(2.9673_dp/13.353384_dp)**3* &
815  (2.0_dp + 2.0_dp*13.353384_dp*x1/2.9673_dp + (13.353384_dp*x1/2.9673_dp)**2) - &
816  (((0.1781_dp/7.0_dp*(2.9673_dp/x1)**2 + 0.4253785_dp/5.0_dp)*(2.9673_dp/x1)**2 + &
817  1.3732412_dp/3.0_dp)*(2.9673_dp/x1)**3)*2.9673_dp**3))/kelvin
818  END IF
819  END IF
820 
821  ! ALLOCATE temporary arrays
822  ALLOCATE (helium_env(k)%helium%itmp_atoms_1d(helium_env(k)%helium%atoms))
823  ALLOCATE (helium_env(k)%helium%ltmp_atoms_1d(helium_env(k)%helium%atoms))
824  ALLOCATE (helium_env(k)%helium%itmp_atoms_np_1d(helium_env(k)%helium%atoms*helium_env(k)%helium%num_env))
825  ALLOCATE (helium_env(k)%helium%rtmp_3_np_1d(3*helium_env(k)%helium%num_env))
826  ALLOCATE (helium_env(k)%helium%rtmp_3_atoms_beads_1d(3*helium_env(k)%helium%atoms* &
827  helium_env(k)%helium%beads))
828  ALLOCATE (helium_env(k)%helium%rtmp_3_atoms_beads_np_1d(3*helium_env(k)%helium%atoms* &
829  helium_env(k)%helium%beads*helium_env(k)%helium%num_env))
830  ALLOCATE (helium_env(k)%helium%ltmp_3_atoms_beads_3d(3, helium_env(k)%helium%atoms, &
831  helium_env(k)%helium%beads))
832  IF (k .EQ. 1) THEN
833  CALL helium_write_setup(helium_env(k)%helium)
834  END IF
835  END DO
836  DEALLOCATE (pot_transfer)
837  ELSE
838  ! Deallocate env_all on processors without helium_env
839  DEALLOCATE (env_all)
840  END IF ! mepos .GT. 0
841 
842  NULLIFY (env_all)
843  CALL timestop(handle)
844 
845  RETURN
846  END SUBROUTINE helium_create
847 
848 ! ***************************************************************************
849 !> \brief Releases helium_solvent_type
850 !> \param helium_env ...
851 !> \par History
852 !> 2016-07-14 Modified to work with independent helium_env [cschran]
853 !> \author hforbert
854 ! **************************************************************************************************
855  SUBROUTINE helium_release(helium_env)
856  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
857 
858  INTEGER :: k
859 
860  IF (ASSOCIATED(helium_env)) THEN
861  DO k = 1, SIZE(helium_env)
862  IF (k .EQ. 1) THEN
863  CALL helium_env(k)%comm%free()
864  DEALLOCATE (helium_env(k)%comm)
865  DEALLOCATE (helium_env(k)%env_all)
866  END IF
867  NULLIFY (helium_env(k)%env_all)
868 
869  ! DEALLOCATE temporary arrays
870  DEALLOCATE ( &
871  helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
872  helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
873  helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
874  helium_env(k)%helium%rtmp_3_np_1d, &
875  helium_env(k)%helium%itmp_atoms_np_1d, &
876  helium_env(k)%helium%ltmp_atoms_1d, &
877  helium_env(k)%helium%itmp_atoms_1d)
878 
879  NULLIFY ( &
880  helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
881  helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
882  helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
883  helium_env(k)%helium%rtmp_3_np_1d, &
884  helium_env(k)%helium%itmp_atoms_np_1d, &
885  helium_env(k)%helium%ltmp_atoms_1d, &
886  helium_env(k)%helium%itmp_atoms_1d &
887  )
888 
889  IF (helium_env(k)%helium%solute_present) THEN
890  ! DEALLOCATE solute-related arrays
891  DEALLOCATE ( &
892  helium_env(k)%helium%rtmp_p_ndim_2d, &
893  helium_env(k)%helium%rtmp_p_ndim_np_1d, &
894  helium_env(k)%helium%rtmp_p_ndim_1d, &
895  helium_env(k)%helium%force_inst, &
896  helium_env(k)%helium%force_avrg)
897  NULLIFY ( &
898  helium_env(k)%helium%rtmp_p_ndim_2d, &
899  helium_env(k)%helium%rtmp_p_ndim_np_1d, &
900  helium_env(k)%helium%rtmp_p_ndim_1d, &
901  helium_env(k)%helium%force_inst, &
902  helium_env(k)%helium%force_avrg)
903  END IF
904 
905  IF (helium_env(k)%helium%rho_present) THEN
906  DEALLOCATE ( &
907  helium_env(k)%helium%rho_rstr, &
908  helium_env(k)%helium%rho_accu, &
909  helium_env(k)%helium%rho_inst, &
910  helium_env(k)%helium%rho_incr)
911  NULLIFY ( &
912  helium_env(k)%helium%rho_rstr, &
913  helium_env(k)%helium%rho_accu, &
914  helium_env(k)%helium%rho_inst, &
915  helium_env(k)%helium%rho_incr)
916  ! DEALLOCATE everything
917  DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix)
918  DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_name)
919  DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_index)
920  NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix)
921  NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_name)
922  NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_index)
923  DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix)
924  DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_name)
925  DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_index)
926  NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix)
927  NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_name)
928  NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_index)
929  DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix)
930  DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name)
931  DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index)
932  NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix)
933  NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name)
934  NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index)
935  DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix)
936  DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_name)
937  DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_index)
938  NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix)
939  NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_name)
940  NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_index)
941  DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix)
942  DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name)
943  DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index)
944  NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix)
945  NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name)
946  NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index)
947  DEALLOCATE (helium_env(k)%helium%rho_property)
948  NULLIFY (helium_env(k)%helium%rho_property)
949  END IF
950 
951  CALL helium_rdf_release(helium_env(k)%helium)
952 
953  ! DEALLOCATE helium-related arrays
954  DEALLOCATE ( &
955  helium_env(k)%helium%atom_plength, &
956  helium_env(k)%helium%plength_inst, &
957  helium_env(k)%helium%plength_avrg, &
958  helium_env(k)%helium%num_accepted, &
959  helium_env(k)%helium%ipmatrix, &
960  helium_env(k)%helium%pmatrix, &
961  helium_env(k)%helium%nmatrix, &
962  helium_env(k)%helium%tmatrix, &
963  helium_env(k)%helium%iperm, &
964  helium_env(k)%helium%permutation, &
965  helium_env(k)%helium%ptable, &
966  helium_env(k)%helium%work, &
967  helium_env(k)%helium%pos)
968  IF (helium_env(k)%helium%worm_max_open_cycles > 0) THEN
969  DEALLOCATE (helium_env(k)%helium%saveiperm, &
970  helium_env(k)%helium%savepermutation, &
971  helium_env(k)%helium%savepos)
972  END IF
973  NULLIFY ( &
974  helium_env(k)%helium%atom_plength, &
975  helium_env(k)%helium%plength_inst, &
976  helium_env(k)%helium%plength_avrg, &
977  helium_env(k)%helium%num_accepted, &
978  helium_env(k)%helium%ipmatrix, &
979  helium_env(k)%helium%pmatrix, &
980  helium_env(k)%helium%nmatrix, &
981  helium_env(k)%helium%tmatrix, &
982  helium_env(k)%helium%iperm, &
983  helium_env(k)%helium%saveiperm, &
984  helium_env(k)%helium%permutation, &
985  helium_env(k)%helium%savepermutation, &
986  helium_env(k)%helium%ptable, &
987  helium_env(k)%helium%work, &
988  helium_env(k)%helium%pos, &
989  helium_env(k)%helium%savepos &
990  )
991 
992  IF (k == 1) THEN
993  CALL spline_data_release(helium_env(k)%helium%vij)
994  CALL spline_data_release(helium_env(k)%helium%u0)
995  CALL spline_data_release(helium_env(k)%helium%e0)
996  DEALLOCATE (helium_env(k)%helium%uoffdiag)
997  DEALLOCATE (helium_env(k)%helium%eoffdiag)
998  END IF
999  NULLIFY (helium_env(k)%helium%uoffdiag, &
1000  helium_env(k)%helium%eoffdiag, &
1001  helium_env(k)%helium%vij, &
1002  helium_env(k)%helium%u0, &
1003  helium_env(k)%helium%e0)
1004 
1005  DEALLOCATE (helium_env(k)%helium%rng_stream_uniform)
1006  DEALLOCATE (helium_env(k)%helium%rng_stream_gaussian)
1007 
1008  ! deallocate solute-related arrays
1009  IF (helium_env(k)%helium%solute_present) THEN
1010  DEALLOCATE (helium_env(k)%helium%solute_element)
1011  NULLIFY (helium_env(k)%helium%solute_element)
1012  END IF
1013 
1014  ! Deallocate everything from the helium_set_solute_indices
1015  IF (ASSOCIATED(helium_env(k)%helium%ename)) THEN
1016  DEALLOCATE (helium_env(k)%helium%ename)
1017  NULLIFY (helium_env(k)%helium%ename)
1018  END IF
1019 
1020  ! NNP interaction
1021  IF (ASSOCIATED(helium_env(k)%helium%nnp)) THEN
1022  CALL nnp_env_release(helium_env(k)%helium%nnp)
1023  DEALLOCATE (helium_env(k)%helium%nnp)
1024  NULLIFY (helium_env(k)%helium%nnp)
1025  END IF
1026  IF (ASSOCIATED(helium_env(k)%helium%nnp_sr_cut)) THEN
1027  DEALLOCATE (helium_env(k)%helium%nnp_sr_cut)
1028  NULLIFY (helium_env(k)%helium%nnp_sr_cut)
1029  END IF
1030 
1031  DEALLOCATE (helium_env(k)%helium)
1032 
1033  END DO
1034 
1035  DEALLOCATE (helium_env)
1036  NULLIFY (helium_env)
1037  END IF
1038  RETURN
1039  END SUBROUTINE helium_release
1040 
1041 ! ***************************************************************************
1042 !> \brief Initialize helium data structures.
1043 !> \param helium_env ...
1044 !> \param pint_env ...
1045 !> \par History
1046 !> removed references to pint_env_type data structure [lwalewski]
1047 !> 2009-11-10 init/restore coords, perm, RNG and forces [lwalewski]
1048 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1049 !> \author hforbert
1050 !> \note Initializes helium coordinates either as random positions or from
1051 !> HELIUM%COORD section if it's present in the input file.
1052 !> Initializes helium permutation state as identity permutation or
1053 !> from HELIUM%PERM section if it's present in the input file.
1054 ! **************************************************************************************************
1055  SUBROUTINE helium_init(helium_env, pint_env)
1056 
1057  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1058  TYPE(pint_env_type), INTENT(IN) :: pint_env
1059 
1060  CHARACTER(len=*), PARAMETER :: routinen = 'helium_init'
1061 
1062  INTEGER :: handle, k
1063  LOGICAL :: coords_presampled, explicit, presample
1064  REAL(kind=dp) :: initkt, solute_radius
1065  TYPE(cp_logger_type), POINTER :: logger
1066  TYPE(section_vals_type), POINTER :: helium_section, sec
1067 
1068  CALL timeset(routinen, handle)
1069 
1070  NULLIFY (logger)
1071  logger => cp_get_default_logger()
1072 
1073  IF (ASSOCIATED(helium_env)) THEN
1074 
1075  NULLIFY (helium_section)
1076  helium_section => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1077  "MOTION%PINT%HELIUM")
1078 
1079  ! restore RNG state
1080  NULLIFY (sec)
1081  sec => section_vals_get_subs_vals(helium_section, "RNG_STATE")
1082  CALL section_vals_get(sec, explicit=explicit)
1083  IF (explicit) THEN
1084  CALL helium_rng_restore(helium_env)
1085  ELSE
1086  CALL helium_write_line("RNG state initialized as new.")
1087  END IF
1088 
1089  ! init/restore permutation state
1090  NULLIFY (sec)
1091  sec => section_vals_get_subs_vals(helium_section, "PERM")
1092  CALL section_vals_get(sec, explicit=explicit)
1093  IF (explicit) THEN
1094  CALL helium_perm_restore(helium_env)
1095  ELSE
1096  CALL helium_perm_init(helium_env)
1097  CALL helium_write_line("Permutation state initialized as identity.")
1098  END IF
1099 
1100  ! Specify if forces should be obtained as AVG or LAST
1101  DO k = 1, SIZE(helium_env)
1102  CALL section_vals_val_get(helium_section, "GET_FORCES", &
1103  i_val=helium_env(k)%helium%get_helium_forces)
1104  END DO
1105 
1106  DO k = 1, SIZE(helium_env)
1107  ! init center of mass
1108  IF (helium_env(k)%helium%solute_present) THEN
1109  helium_env(k)%helium%center(:) = pint_com_pos(pint_env)
1110  ELSE
1111  helium_env(k)%helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
1112  END IF
1113  END DO
1114 
1115  ! init/restore coordinates
1116  NULLIFY (sec)
1117  sec => section_vals_get_subs_vals(helium_section, "COORD")
1118  CALL section_vals_get(sec, explicit=explicit)
1119  IF (explicit) THEN
1120  CALL helium_coord_restore(helium_env)
1121  CALL helium_write_line("Coordinates restarted.")
1122  ELSE
1123  CALL section_vals_val_get(helium_section, "COORD_INIT_TEMP", r_val=initkt)
1124  CALL section_vals_val_get(helium_section, "SOLUTE_RADIUS", r_val=solute_radius)
1125  CALL helium_coord_init(helium_env, initkt, solute_radius)
1126  IF (initkt > 0.0_dp) THEN
1127  CALL helium_write_line("Coordinates initialized with thermal gaussian.")
1128  ELSE
1129  CALL helium_write_line("Coordinates initialized as point particles.")
1130  END IF
1131  END IF
1132 
1133  DO k = 1, SIZE(helium_env)
1134 
1135  helium_env(k)%helium%worm_is_closed = .true.
1136  helium_env(k)%helium%worm_atom_idx = 0
1137  helium_env(k)%helium%worm_bead_idx = 0
1138 
1139  helium_env(k)%helium%work(:, :, :) = helium_env(k)%helium%pos(:, :, :)
1140 
1141  ! init center of mass
1142  IF (helium_env(k)%helium%solute_present) THEN
1143  helium_env(k)%helium%center(:) = pint_com_pos(pint_env)
1144  ELSE
1145  IF (helium_env(k)%helium%periodic) THEN
1146  helium_env(k)%helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
1147  ELSE
1148  helium_env(k)%helium%center(:) = helium_com(helium_env(k)%helium)
1149  END IF
1150  END IF
1151  END DO
1152 
1153  ! Optional helium coordinate presampling:
1154  ! Assume IONODE to have always at least one helium_env
1155  CALL section_vals_val_get(helium_section, "PRESAMPLE", &
1156  l_val=presample)
1157  coords_presampled = .false.
1158  IF (presample) THEN
1159  DO k = 1, SIZE(helium_env)
1160  helium_env(k)%helium%current_step = 0
1161  END DO
1162  CALL helium_sample(helium_env, pint_env)
1163  DO k = 1, SIZE(helium_env)
1164  IF (helium_env(k)%helium%solute_present) helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
1165  helium_env(k)%helium%energy_avrg(:) = 0.0_dp
1166  helium_env(k)%helium%plength_avrg(:) = 0.0_dp
1167  helium_env(k)%helium%num_accepted(:, :) = 0.0_dp
1168  ! Reset properties accumulated over presample:
1169  helium_env(k)%helium%proarea%accu(:) = 0.0_dp
1170  helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
1171  helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
1172  helium_env(k)%helium%mominer%accu(:) = 0.0_dp
1173  IF (helium_env(k)%helium%rho_present) THEN
1174  helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
1175  END IF
1176  IF (helium_env(k)%helium%rdf_present) THEN
1177  helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp
1178  END IF
1179  END DO
1180  coords_presampled = .true.
1181  CALL helium_write_line("Bead coordinates pre-sampled.")
1182  END IF
1183 
1184  IF (helium_env(1)%helium%solute_present) THEN
1185  ! restore helium forces
1186  NULLIFY (sec)
1187  sec => section_vals_get_subs_vals(helium_section, "FORCE")
1188  CALL section_vals_get(sec, explicit=explicit)
1189  IF (explicit) THEN
1190  IF (.NOT. coords_presampled) THEN
1191  CALL helium_force_restore(helium_env)
1192  END IF
1193  ELSE
1194  IF (.NOT. coords_presampled) THEN
1195  CALL helium_force_init(helium_env)
1196  CALL helium_write_line("Forces on the solute initialized as zero.")
1197  END IF
1198  END IF
1199  !! Update pint_env force, assume always one helium_env at IONODE
1200  !IF (pint_env%logger%para_env%is_source()) THEN
1201  ! pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
1202  !END IF
1203  !CALL pint_env%logger%para_env%bcast(pint_env%f,&
1204  ! pint_env%logger%para_env%source)
1205 
1206  END IF
1207  END IF
1208 
1209  CALL timestop(handle)
1210 
1211  RETURN
1212  END SUBROUTINE helium_init
1213 
1214 ! ***************************************************************************
1215 ! Data transfer functions.
1216 !
1217 ! These functions manipulate and transfer data between the runtime
1218 ! environment and the input structure.
1219 ! ***************************************************************************
1220 
1221 ! ***************************************************************************
1222 !> \brief Initialize helium coordinates with random positions.
1223 !> \param helium_env ...
1224 !> \param initkT ...
1225 !> \param solute_radius ...
1226 !> \date 2009-11-09
1227 !> \par History
1228 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1229 !> 2018-04-30 Useful initialization for droplets [fuhl]
1230 !> \author Lukasz Walewski
1231 ! **************************************************************************************************
1232  SUBROUTINE helium_coord_init(helium_env, initkT, solute_radius)
1233  TYPE(helium_solvent_p_type), DIMENSION(:), &
1234  INTENT(INOUT), POINTER :: helium_env
1235  REAL(kind=dp), INTENT(IN) :: initkt, solute_radius
1236 
1237  REAL(kind=dp), PARAMETER :: minhehedst = 5.669177966_dp
1238 
1239  INTEGER :: ia, ib, ic, id, iter, k
1240  LOGICAL :: invalidpos
1241  REAL(kind=dp) :: minhehedsttmp, r1, r2
1242  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: centroids
1243  REAL(kind=dp), DIMENSION(3) :: cvek, rvek, tvek
1244 
1245  !corresponds to three angstrom (roughly first maximum of He-He-rdf)
1246  minhehedsttmp = minhehedst
1247 
1248  DO k = 1, SIZE(helium_env)
1249  IF (helium_env(k)%helium%solute_present) THEN
1250  cvek(:) = helium_env(k)%helium%center(:)
1251  CALL helium_pbc(helium_env(k)%helium, cvek)
1252  END IF
1253 
1254  ALLOCATE (centroids(3, helium_env(k)%helium%atoms))
1255  IF (helium_env(k)%helium%periodic) THEN
1256  DO ia = 1, helium_env(k)%helium%atoms
1257  invalidpos = .true.
1258  iter = 0
1259  DO WHILE (invalidpos)
1260  iter = iter + 1
1261  invalidpos = .false.
1262  ! if sampling fails to often, reduce he he criterion
1263  !CS TODO:
1264  !minHeHedsttmp = 0.90_dp**(iter/100)*minHeHedst
1265  minhehedsttmp = 0.90_dp**min(0, iter - 2)*minhehedst
1266  DO ic = 1, 3
1267  r1 = helium_env(k)%helium%rng_stream_uniform%next()
1268  r1 = 2.0_dp*r1 - 1.0_dp
1269  r1 = r1*helium_env(k)%helium%cell_size
1270  centroids(ic, ia) = r1
1271  END DO
1272  ! check if helium is outside of cell
1273  tvek(:) = centroids(:, ia)
1274  CALL helium_pbc(helium_env(k)%helium, tvek(:))
1275  rvek(:) = tvek(:) - centroids(:, ia)
1276  r2 = dot_product(rvek, rvek)
1277  IF (r2 > 1.0_dp*10.0_dp**(-6)) THEN
1278  invalidpos = .true.
1279  ELSE
1280  ! check for helium-helium collision
1281  DO id = 1, ia - 1
1282  rvek = centroids(:, ia) - centroids(:, id)
1283  CALL helium_pbc(helium_env(k)%helium, rvek)
1284  r2 = dot_product(rvek, rvek)
1285  IF (r2 < minhehedsttmp**2) THEN
1286  invalidpos = .true.
1287  EXIT
1288  END IF
1289  END DO
1290  END IF
1291  IF (.NOT. invalidpos) THEN
1292  ! check if centroid collides with molecule
1293  IF (helium_env(k)%helium%solute_present) THEN
1294  rvek(:) = (cvek(:) - centroids(:, ia))
1295  r2 = dot_product(rvek, rvek)
1296  IF (r2 <= solute_radius**2) invalidpos = .true.
1297  END IF
1298  END IF
1299  END DO
1300  END DO
1301  ! do thermal gaussian delocalization of hot start
1302  IF (initkt > 0.0_dp) THEN
1303  CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkt)
1304  ELSE
1305  DO ia = 1, helium_env(k)%helium%atoms
1306  DO ib = 1, helium_env(k)%helium%beads
1307  helium_env(k)%helium%pos(:, ia, ib) = centroids(:, ia)
1308  END DO
1309  END DO
1310  END IF
1311  ! apply PBC to bead coords
1312  DO ia = 1, helium_env(k)%helium%atoms
1313  DO ib = 1, helium_env(k)%helium%beads
1314  CALL helium_pbc(helium_env(k)%helium, helium_env(k)%helium%pos(:, ia, ib))
1315  ! check if bead collides with molecule
1316  IF (helium_env(k)%helium%solute_present) THEN
1317  rvek(:) = (cvek(:) - helium_env(k)%helium%pos(:, ia, ib))
1318  r2 = dot_product(rvek, rvek)
1319  IF (r2 <= solute_radius**2) THEN
1320  r1 = sqrt(r2)
1321  helium_env(k)%helium%pos(:, ia, ib) = &
1322  cvek(:) + solute_radius/r1*rvek(:)
1323  END IF
1324  END IF
1325  END DO
1326  END DO
1327  ELSE
1328  DO ia = 1, helium_env(k)%helium%atoms
1329  iter = 0
1330  invalidpos = .true.
1331  DO WHILE (invalidpos)
1332  invalidpos = .false.
1333  iter = iter + 1
1334  ! if sampling fails to often, reduce he he criterion
1335  minhehedsttmp = 0.90_dp**min(0, iter - 2)*minhehedst
1336  DO ic = 1, 3
1337  rvek(ic) = helium_env(k)%helium%rng_stream_uniform%next()
1338  rvek(ic) = 2.0_dp*rvek(ic) - 1.0_dp
1339  rvek(ic) = rvek(ic)*helium_env(k)%helium%droplet_radius
1340  END DO
1341  centroids(:, ia) = rvek(:)
1342  ! check if helium is outside of the droplet
1343  r2 = dot_product(rvek, rvek)
1344  IF (r2 > helium_env(k)%helium%droplet_radius**2) THEN
1345  invalidpos = .true.
1346  ELSE
1347  ! check for helium-helium collision
1348  DO id = 1, ia - 1
1349  rvek = centroids(:, ia) - centroids(:, id)
1350  r2 = dot_product(rvek, rvek)
1351  IF (r2 < minhehedsttmp**2) THEN
1352  invalidpos = .true.
1353  EXIT
1354  END IF
1355  END DO
1356  END IF
1357  IF (.NOT. invalidpos) THEN
1358  ! make sure the helium does not collide with the solute
1359  IF (helium_env(k)%helium%solute_present) THEN
1360  rvek(:) = (cvek(:) - centroids(:, ia))
1361  r2 = dot_product(rvek, rvek)
1362  IF (r2 <= solute_radius**2) invalidpos = .true.
1363  END IF
1364  END IF
1365  END DO
1366  END DO
1367  ! do thermal gaussian delocalization of hot start
1368  IF (initkt > 0.0_dp) THEN
1369  CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkt)
1370  ELSE
1371  DO ia = 1, helium_env(k)%helium%atoms
1372  DO ib = 1, helium_env(k)%helium%beads
1373  helium_env(k)%helium%pos(:, ia, ib) = centroids(:, ia)
1374  END DO
1375  END DO
1376  END IF
1377  DO ia = 1, helium_env(k)%helium%atoms
1378  DO ib = 1, helium_env(k)%helium%beads
1379  ! Make sure, that nothing lies outside the droplet radius
1380  r1 = dot_product(helium_env(k)%helium%pos(:, ia, ib), &
1381  helium_env(k)%helium%pos(:, ia, ib))
1382  IF (r1 > helium_env(k)%helium%droplet_radius**2) THEN
1383  r1 = sqrt(r1)
1384  helium_env(k)%helium%pos(:, ia, ib) = &
1385  helium_env(k)%helium%droplet_radius/r1* &
1386  helium_env(k)%helium%pos(:, ia, ib)
1387  ELSE IF (helium_env(k)%helium%solute_present) THEN
1388  IF (r1 < solute_radius**2) THEN
1389  !make sure that nothing lies within the molecule
1390  r1 = sqrt(r1)
1391  helium_env(k)%helium%pos(:, ia, ib) = &
1392  solute_radius/r1* &
1393  helium_env(k)%helium%pos(:, ia, ib)
1394  END IF
1395  END IF
1396  ! transfer to position around actual center of droplet
1397  helium_env(k)%helium%pos(:, ia, ib) = &
1398  helium_env(k)%helium%pos(:, ia, ib) + &
1399  helium_env(k)%helium%center(:)
1400  END DO
1401  END DO
1402  END IF
1403  helium_env(k)%helium%work = helium_env(k)%helium%pos
1404  DEALLOCATE (centroids)
1405  END DO
1406 
1407  RETURN
1408  END SUBROUTINE helium_coord_init
1409 
1410 ! **************************************************************************************************
1411 !> \brief ...
1412 !> \param helium_env ...
1413 !> \param centroids ...
1414 !> \param kbT ...
1415 ! **************************************************************************************************
1416  SUBROUTINE helium_thermal_gaussian_beads_init(helium_env, centroids, kbT)
1417 
1418  TYPE(helium_solvent_type), POINTER :: helium_env
1419  REAL(kind=dp), DIMENSION(3, helium_env%atoms), &
1420  INTENT(IN) :: centroids
1421  REAL(kind=dp), INTENT(IN) :: kbt
1422 
1423  INTEGER :: i, iatom, idim, imode, j, p
1424  REAL(kind=dp) :: invsqrtp, omega, pip, rand, sqrt2p, &
1425  sqrtp, twopip, variance
1426  REAL(kind=dp), &
1427  DIMENSION(helium_env%beads, helium_env%beads) :: u2x
1428  REAL(kind=dp), DIMENSION(helium_env%beads) :: nmhecoords
1429 
1430  p = helium_env%beads
1431 
1432  sqrt2p = sqrt(2.0_dp/real(p, dp))
1433  twopip = twopi/real(p, dp)
1434  pip = pi/real(p, dp)
1435  sqrtp = sqrt(real(p, dp))
1436  invsqrtp = 1.0_dp/sqrt(real(p, dp))
1437 
1438  ! set up normal mode backtransform matrix
1439  u2x(:, :) = 0.0_dp
1440  u2x(:, 1) = invsqrtp
1441  DO i = 2, p/2 + 1
1442  DO j = 1, p
1443  u2x(j, i) = sqrt2p*cos(twopip*(i - 1)*(j - 1))
1444  END DO
1445  END DO
1446  DO i = p/2 + 2, p
1447  DO j = 1, p
1448  u2x(j, i) = sqrt2p*sin(twopip*(i - 1)*(j - 1))
1449  END DO
1450  END DO
1451  IF (mod(p, 2) == 0) THEN
1452  DO i = 1, p - 1, 2
1453  u2x(i, p/2 + 1) = invsqrtp
1454  u2x(i + 1, p/2 + 1) = -1.0_dp*invsqrtp
1455  END DO
1456  END IF
1457 
1458  DO iatom = 1, helium_env%atoms
1459  DO idim = 1, 3
1460  nmhecoords(1) = sqrtp*centroids(idim, iatom)
1461  DO imode = 2, p
1462  omega = 2.0_dp*p*kbt*sin((imode - 1)*pip)
1463  variance = kbt*p/(helium_env%he_mass_au*omega**2)
1464  rand = helium_env%rng_stream_gaussian%next()
1465  nmhecoords(imode) = rand*sqrt(variance)
1466  END DO
1467  helium_env%pos(idim, iatom, 1:p) = matmul(u2x, nmhecoords)
1468  END DO
1469  END DO
1470 
1471  END SUBROUTINE helium_thermal_gaussian_beads_init
1472 
1473 ! ***************************************************************************
1474 !> \brief Restore coordinates from the input structure.
1475 !> \param helium_env ...
1476 !> \date 2009-11-09
1477 !> \par History
1478 !> 2010-07-22 accommodate additional cpus in the runtime wrt the
1479 !> restart [lwalewski]
1480 !> 2016-07-14 Modified to work with independent helium_env
1481 !> [cschran]
1482 !> \author Lukasz Walewski
1483 ! **************************************************************************************************
1484  SUBROUTINE helium_coord_restore(helium_env)
1485  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1486 
1487  CHARACTER(len=default_string_length) :: err_str, stmp
1488  INTEGER :: actlen, i, k, msglen, num_env_restart, &
1489  off, offset
1490  LOGICAL, DIMENSION(:, :, :), POINTER :: m
1491  REAL(kind=dp), DIMENSION(:), POINTER :: message
1492  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: f
1493  TYPE(cp_logger_type), POINTER :: logger
1494 
1495  NULLIFY (logger)
1496  logger => cp_get_default_logger()
1497 
1498  ! assign the pointer to the memory location of the input structure, where
1499  ! the coordinates are stored
1500  NULLIFY (message)
1501  CALL section_vals_val_get(helium_env(1)%helium%input, &
1502  "MOTION%PINT%HELIUM%COORD%_DEFAULT_KEYWORD_", &
1503  r_vals=message)
1504 
1505  ! check that the number of values in the input match the current runtime
1506  actlen = SIZE(message)
1507  num_env_restart = actlen/helium_env(1)%helium%atoms/helium_env(1)%helium%beads/3
1508 
1509  IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
1510  err_str = "Reading bead coordinates from the input file."
1511  CALL helium_write_line(err_str)
1512  err_str = "Number of environments in the restart...: '"
1513  stmp = ""
1514  WRITE (stmp, *) num_env_restart
1515  err_str = trim(adjustl(err_str))//trim(adjustl(stmp))//"'."
1516  CALL helium_write_line(err_str)
1517  err_str = "Number of current run time environments.: '"
1518  stmp = ""
1519  WRITE (stmp, *) helium_env(1)%helium%num_env
1520  err_str = trim(adjustl(err_str))//trim(adjustl(stmp))//"'."
1521  CALL helium_write_line(err_str)
1522  err_str = "Missmatch between number of bead coord. in input file and helium environments."
1523  cpabort(err_str)
1524  ELSE
1525  CALL helium_write_line("Bead coordinates read from the input file.")
1526 
1527  offset = 0
1528  DO i = 1, logger%para_env%mepos
1529  offset = offset + helium_env(1)%env_all(i)
1530  END DO
1531 
1532  ! distribute coordinates over processors (no message passing)
1533  DO k = 1, SIZE(helium_env)
1534  msglen = helium_env(k)%helium%atoms*helium_env(k)%helium%beads*3
1535  off = msglen*mod(offset + k - 1, num_env_restart)
1536  NULLIFY (m, f)
1537  ALLOCATE (m(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads))
1538  ALLOCATE (f(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads))
1539  m(:, :, :) = .true.
1540  f(:, :, :) = 0.0_dp
1541  helium_env(k)%helium%pos(:, :, 1:helium_env(k)%helium%beads) = unpack(message(off + 1:off + msglen), mask=m, field=f)
1542  DEALLOCATE (f, m)
1543  END DO
1544 
1545  END IF
1546 
1547  NULLIFY (message)
1548 
1549  RETURN
1550  END SUBROUTINE helium_coord_restore
1551 
1552 ! ***************************************************************************
1553 !> \brief Initialize forces exerted on the solute
1554 !> \param helium_env ...
1555 !> \date 2009-11-10
1556 !> \par History
1557 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1558 !> \author Lukasz Walewski
1559 ! **************************************************************************************************
1560  SUBROUTINE helium_force_init(helium_env)
1561 
1562  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1563 
1564  INTEGER :: k
1565 
1566  DO k = 1, SIZE(helium_env)
1567  IF (helium_env(k)%helium%solute_present) THEN
1568  helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
1569  helium_env(k)%helium%force_inst(:, :) = 0.0_dp
1570  END IF
1571  END DO
1572 
1573  RETURN
1574  END SUBROUTINE helium_force_init
1575 
1576 ! ***************************************************************************
1577 !> \brief Restore forces from the input structure to the runtime environment.
1578 !> \param helium_env ...
1579 !> \date 2009-11-10
1580 !> \par History
1581 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1582 !> \author Lukasz Walewski
1583 ! **************************************************************************************************
1584  SUBROUTINE helium_force_restore(helium_env)
1585  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1586 
1587  CHARACTER(len=default_string_length) :: err_str, stmp
1588  INTEGER :: actlen, k, msglen
1589  LOGICAL, DIMENSION(:, :), POINTER :: m
1590  REAL(kind=dp), DIMENSION(:), POINTER :: message
1591  REAL(kind=dp), DIMENSION(:, :), POINTER :: f
1592 
1593 ! assign the pointer to the memory location of the input structure, where
1594 ! the forces are stored
1595 
1596  NULLIFY (message)
1597  CALL section_vals_val_get(helium_env(1)%helium%input, &
1598  "MOTION%PINT%HELIUM%FORCE%_DEFAULT_KEYWORD_", &
1599  r_vals=message)
1600 
1601  ! check if the destination array has correct size
1602  msglen = helium_env(1)%helium%solute_atoms*helium_env(1)%helium%solute_beads*3
1603  actlen = SIZE(helium_env(1)%helium%force_avrg)
1604  err_str = "Invalid size of helium%force_avrg array: actual '"
1605  stmp = ""
1606  WRITE (stmp, *) actlen
1607  err_str = trim(adjustl(err_str))// &
1608  trim(adjustl(stmp))//"' but expected '"
1609  stmp = ""
1610  WRITE (stmp, *) msglen
1611  IF (actlen /= msglen) THEN
1612  err_str = trim(adjustl(err_str))// &
1613  trim(adjustl(stmp))//"'."
1614  cpabort(err_str)
1615  END IF
1616 
1617  ! restore forces on all processors (no message passing)
1618  NULLIFY (m, f)
1619  ALLOCATE (m(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3))
1620  ALLOCATE (f(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3))
1621  m(:, :) = .true.
1622  f(:, :) = 0.0_dp
1623  DO k = 1, SIZE(helium_env)
1624  helium_env(k)%helium%force_avrg(:, :) = unpack(message(1:msglen), mask=m, field=f)
1625  helium_env(k)%helium%force_inst(:, :) = 0.0_dp
1626  END DO
1627  DEALLOCATE (f, m)
1628 
1629  CALL helium_write_line("Forces on the solute read from the input file.")
1630 
1631  NULLIFY (message)
1632 
1633  RETURN
1634  END SUBROUTINE helium_force_restore
1635 
1636 ! ***************************************************************************
1637 !> \brief Initialize the permutation state.
1638 !> \param helium_env ...
1639 !> \date 2009-11-05
1640 !> \par History
1641 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1642 !> \author Lukasz Walewski
1643 !> \note Assign the identity permutation at each processor. Inverse
1644 !> permutation array gets assigned as well.
1645 ! **************************************************************************************************
1646  SUBROUTINE helium_perm_init(helium_env)
1647  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1648 
1649  INTEGER :: ia, k
1650 
1651  DO k = 1, SIZE(helium_env)
1652  DO ia = 1, helium_env(k)%helium%atoms
1653  helium_env(k)%helium%permutation(ia) = ia
1654  helium_env(k)%helium%iperm(ia) = ia
1655  END DO
1656  END DO
1657 
1658  RETURN
1659  END SUBROUTINE helium_perm_init
1660 
1661 ! ***************************************************************************
1662 !> \brief Restore permutation state from the input structure.
1663 !> \param helium_env ...
1664 !> \date 2009-11-05
1665 !> \par History
1666 !> 2010-07-22 accommodate additional cpus in the runtime wrt the
1667 !> restart [lwalewski]
1668 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1669 !> \author Lukasz Walewski
1670 !> \note Transfer permutation state from the input tree to the runtime
1671 !> data structures on each processor. Inverse permutation array is
1672 !> recalculated according to the restored permutation state.
1673 ! **************************************************************************************************
1674  SUBROUTINE helium_perm_restore(helium_env)
1675  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1676 
1677  CHARACTER(len=default_string_length) :: err_str, stmp
1678  INTEGER :: actlen, i, ia, ic, k, msglen, &
1679  num_env_restart, off, offset
1680  INTEGER, DIMENSION(:), POINTER :: message
1681  TYPE(cp_logger_type), POINTER :: logger
1682 
1683  NULLIFY (logger)
1684  logger => cp_get_default_logger()
1685 
1686  ! assign the pointer to the memory location of the input structure, where
1687  ! the permutation state is stored
1688  NULLIFY (message)
1689  CALL section_vals_val_get(helium_env(1)%helium%input, &
1690  "MOTION%PINT%HELIUM%PERM%_DEFAULT_KEYWORD_", &
1691  i_vals=message)
1692 
1693  ! check the number of environments presumably stored in the restart
1694  actlen = SIZE(message)
1695  num_env_restart = actlen/helium_env(1)%helium%atoms
1696 
1697 !TODO maybe add some sanity checks here:
1698 ! is num_env_restart integer ?
1699 
1700  IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
1701  err_str = "Reading permutation state from the input file."
1702  CALL helium_write_line(err_str)
1703  err_str = "Number of environments in the restart...: '"
1704  stmp = ""
1705  WRITE (stmp, *) num_env_restart
1706  err_str = trim(adjustl(err_str))//trim(adjustl(stmp))//"'."
1707  CALL helium_write_line(err_str)
1708  err_str = "Number of current run time environments.: '"
1709  stmp = ""
1710  WRITE (stmp, *) helium_env(1)%helium%num_env
1711  err_str = trim(adjustl(err_str))//trim(adjustl(stmp))//"'."
1712  CALL helium_write_line(err_str)
1713  err_str = "Missmatch between number of perm. states in input file and helium environments."
1714  cpabort(err_str)
1715  ELSE
1716  CALL helium_write_line("Permutation state read from the input file.")
1717 
1718  ! distribute permutation state over processors
1719  offset = 0
1720  DO i = 1, logger%para_env%mepos
1721  offset = offset + helium_env(1)%env_all(i)
1722  END DO
1723 
1724  DO k = 1, SIZE(helium_env)
1725  msglen = helium_env(k)%helium%atoms
1726  off = msglen*mod(k - 1 + offset, num_env_restart)
1727  helium_env(k)%helium%permutation(:) = message(off + 1:off + msglen)
1728  END DO
1729  END IF
1730 
1731  ! recalculate the inverse permutation array
1732  DO k = 1, SIZE(helium_env)
1733  helium_env(k)%helium%iperm(:) = 0
1734  ic = 0
1735  DO ia = 1, msglen
1736  IF ((helium_env(k)%helium%permutation(ia) > 0) .AND. (helium_env(k)%helium%permutation(ia) <= msglen)) THEN
1737  helium_env(k)%helium%iperm(helium_env(k)%helium%permutation(ia)) = ia
1738  ic = ic + 1
1739  END IF
1740  END DO
1741  err_str = "Invalid HELIUM%PERM state: some numbers not within (1,"
1742  stmp = ""
1743  WRITE (stmp, *) msglen
1744  IF (ic /= msglen) THEN
1745  err_str = trim(adjustl(err_str))// &
1746  trim(adjustl(stmp))//")."
1747  cpabort(err_str)
1748  END IF
1749  END DO
1750  NULLIFY (message)
1751 
1752  RETURN
1753  END SUBROUTINE helium_perm_restore
1754 
1755 ! ***************************************************************************
1756 !> \brief Restore averages from the input structure
1757 !> \param helium_env ...
1758 !> \date 2014-06-25
1759 !> \par History
1760 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1761 !> \author Lukasz Walewski
1762 ! **************************************************************************************************
1763  SUBROUTINE helium_averages_restore(helium_env)
1764 
1765  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1766 
1767  INTEGER :: i, k, msglen, num_env_restart, off, &
1768  offset
1769  LOGICAL :: explicit
1770  REAL(kind=dp), DIMENSION(:), POINTER :: message
1771  TYPE(cp_logger_type), POINTER :: logger
1772 
1773  NULLIFY (logger)
1774  logger => cp_get_default_logger()
1775 
1776  offset = 0
1777  DO i = 1, logger%para_env%mepos
1778  offset = offset + helium_env(1)%env_all(i)
1779  END DO
1780 
1781  ! restore projected area
1782  CALL section_vals_val_get(helium_env(1)%helium%input, &
1783  "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", &
1784  explicit=explicit)
1785  IF (explicit) THEN
1786  NULLIFY (message)
1787  CALL section_vals_val_get(helium_env(1)%helium%input, &
1788  "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", &
1789  r_vals=message)
1790  num_env_restart = SIZE(message)/3 ! apparent number of environments
1791  msglen = 3
1792  DO k = 1, SIZE(helium_env)
1793  off = msglen*mod(offset + k - 1, num_env_restart)
1794  helium_env(k)%helium%proarea%rstr(:) = message(off + 1:off + msglen)
1795  END DO
1796  ELSE
1797  DO k = 1, SIZE(helium_env)
1798  helium_env(k)%helium%proarea%rstr(:) = 0.0_dp
1799  END DO
1800  END IF
1801 
1802  ! restore projected area squared
1803  CALL section_vals_val_get(helium_env(1)%helium%input, &
1804  "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", &
1805  explicit=explicit)
1806  IF (explicit) THEN
1807  NULLIFY (message)
1808  CALL section_vals_val_get(helium_env(1)%helium%input, &
1809  "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", &
1810  r_vals=message)
1811  num_env_restart = SIZE(message)/3 ! apparent number of environments
1812  msglen = 3
1813  DO k = 1, SIZE(helium_env)
1814  off = msglen*mod(offset + k - 1, num_env_restart)
1815  helium_env(k)%helium%prarea2%rstr(:) = message(off + 1:off + msglen)
1816  END DO
1817  ELSE
1818  DO k = 1, SIZE(helium_env)
1819  helium_env(k)%helium%prarea2%rstr(:) = 0.0_dp
1820  END DO
1821  END IF
1822 
1823  ! restore winding number squared
1824  CALL section_vals_val_get(helium_env(1)%helium%input, &
1825  "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", &
1826  explicit=explicit)
1827  IF (explicit) THEN
1828  NULLIFY (message)
1829  CALL section_vals_val_get(helium_env(1)%helium%input, &
1830  "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", &
1831  r_vals=message)
1832  num_env_restart = SIZE(message)/3 ! apparent number of environments
1833  msglen = 3
1834  DO k = 1, SIZE(helium_env)
1835  off = msglen*mod(offset + k - 1, num_env_restart)
1836  helium_env(k)%helium%wnmber2%rstr(:) = message(off + 1:off + msglen)
1837  END DO
1838  ELSE
1839  DO k = 1, SIZE(helium_env)
1840  helium_env(k)%helium%wnmber2%rstr(:) = 0.0_dp
1841  END DO
1842  END IF
1843 
1844  ! restore moment of inertia
1845  CALL section_vals_val_get(helium_env(1)%helium%input, &
1846  "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", &
1847  explicit=explicit)
1848  IF (explicit) THEN
1849  NULLIFY (message)
1850  CALL section_vals_val_get(helium_env(1)%helium%input, &
1851  "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", &
1852  r_vals=message)
1853  num_env_restart = SIZE(message)/3 ! apparent number of environments
1854  msglen = 3
1855  DO k = 1, SIZE(helium_env)
1856  off = msglen*mod(offset + k - 1, num_env_restart)
1857  helium_env(k)%helium%mominer%rstr(:) = message(off + 1:off + msglen)
1858  END DO
1859  ELSE
1860  DO k = 1, SIZE(helium_env)
1861  helium_env(k)%helium%mominer%rstr(:) = 0.0_dp
1862  END DO
1863  END IF
1864 
1865  IF (helium_env(1)%helium%rdf_present) THEN
1866  CALL helium_rdf_restore(helium_env)
1867  END IF
1868 
1869  IF (helium_env(1)%helium%rho_present) THEN
1870  ! restore densities
1871  CALL helium_rho_restore(helium_env)
1872  END IF
1873 
1874  ! get the weighting factor
1875  DO k = 1, SIZE(helium_env)
1876  CALL section_vals_val_get( &
1877  helium_env(k)%helium%input, &
1878  "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
1879  i_val=helium_env(k)%helium%averages_iweight)
1880 
1881  ! set the flag indicating whether the averages have been restarted
1882  CALL section_vals_val_get( &
1883  helium_env(k)%helium%input, &
1884  "EXT_RESTART%RESTART_HELIUM_AVERAGES", &
1885  l_val=helium_env(k)%helium%averages_restarted)
1886  END DO
1887 
1888  RETURN
1889  END SUBROUTINE helium_averages_restore
1890 
1891 ! ***************************************************************************
1892 !> \brief Create RNG streams and initialize their state.
1893 !> \param helium_env ...
1894 !> \date 2009-11-04
1895 !> \par History
1896 !> 2016-07-14 Modified to work with independent helium_env [cschran]
1897 !> \author Lukasz Walewski
1898 !> \note TODO: This function shouldn't create (allocate) objects! Only
1899 !> initialization, i.e. setting the seed values etc, should be done
1900 !> here, allocation should be moved to helium_create
1901 ! **************************************************************************************************
1902  SUBROUTINE helium_rng_init(helium_env)
1903  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1904 
1905  INTEGER :: helium_seed, i, offset
1906  REAL(kind=dp), DIMENSION(3, 2) :: initial_seed
1907  TYPE(cp_logger_type), POINTER :: logger
1908  TYPE(rng_stream_p_type), DIMENSION(:), POINTER :: gaussian_array, uniform_array
1909 
1910  NULLIFY (logger)
1911  logger => cp_get_default_logger()
1912  IF (logger%para_env%is_source()) THEN
1913  CALL section_vals_val_get(helium_env(1)%helium%input, &
1914  "MOTION%PINT%HELIUM%RNG_SEED", &
1915  i_val=helium_seed)
1916  END IF
1917  CALL helium_env(1)%comm%bcast(helium_seed, &
1918  logger%para_env%source)
1919  initial_seed(:, :) = real(helium_seed, dp)
1920 
1921  ALLOCATE (uniform_array(helium_env(1)%helium%num_env), &
1922  gaussian_array(helium_env(1)%helium%num_env))
1923  DO i = 1, helium_env(1)%helium%num_env
1924  ALLOCATE (uniform_array(i)%stream, gaussian_array(i)%stream)
1925  END DO
1926 
1927  ! Create num_env RNG streams on processor all processors
1928  ! and distribute them so, that each processor gets unique
1929  ! RN sequences for his helium environments
1930  ! COMMENT: rng_stream can not be used with mp_bcast
1931 
1932  uniform_array(1)%stream = rng_stream_type(name="helium_rns_uniform", &
1933  distribution_type=uniform, &
1934  extended_precision=.true., &
1935  seed=initial_seed)
1936 
1937  gaussian_array(1)%stream = rng_stream_type(name="helium_rns_gaussian", &
1938  distribution_type=gaussian, &
1939  extended_precision=.true., &
1940  last_rng_stream=uniform_array(1)%stream)
1941  DO i = 2, helium_env(1)%helium%num_env
1942  uniform_array(i)%stream = rng_stream_type(name="helium_rns_uniform", &
1943  distribution_type=uniform, &
1944  extended_precision=.true., &
1945  last_rng_stream=gaussian_array(i - 1)%stream)
1946 
1947  gaussian_array(i)%stream = rng_stream_type(name="helium_rns_uniform", &
1948  distribution_type=gaussian, &
1949  extended_precision=.true., &
1950  last_rng_stream=uniform_array(i)%stream)
1951  END DO
1952 
1953  offset = 0
1954  DO i = 1, logger%para_env%mepos
1955  offset = offset + helium_env(1)%env_all(i)
1956  END DO
1957 
1958  DO i = 1, SIZE(helium_env)
1959  NULLIFY (helium_env(i)%helium%rng_stream_uniform, &
1960  helium_env(i)%helium%rng_stream_gaussian)
1961  helium_env(i)%helium%rng_stream_uniform => uniform_array(offset + i)%stream
1962  helium_env(i)%helium%rng_stream_gaussian => gaussian_array(offset + i)%stream
1963  END DO
1964 
1965  DO i = 1, helium_env(1)%helium%num_env
1966  IF (i .LE. offset .OR. i .GT. offset + SIZE(helium_env)) THEN
1967  ! only deallocate pointers here which were not passed on to helium_env(*)%helium
1968  DEALLOCATE (uniform_array(i)%stream)
1969  DEALLOCATE (gaussian_array(i)%stream)
1970  END IF
1971  NULLIFY (uniform_array(i)%stream)
1972  NULLIFY (gaussian_array(i)%stream)
1973  END DO
1974 
1975  DEALLOCATE (uniform_array)
1976  DEALLOCATE (gaussian_array)
1977  END SUBROUTINE helium_rng_init
1978 
1979 ! ***************************************************************************
1980 !> \brief Restore RNG state from the input structure.
1981 !> \param helium_env ...
1982 !> \date 2009-11-04
1983 !> \par History
1984 !> 2010-07-22 Create new rng streams if more cpus available in the
1985 !> runtime than in the restart [lwalewski]
1986 !> 2016-04-18 Modified for independet number of helium_env [cschran]
1987 !> \author Lukasz Walewski
1988 ! **************************************************************************************************
1989  SUBROUTINE helium_rng_restore(helium_env)
1990  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1991 
1992  CHARACTER(len=default_string_length) :: err_str, stmp
1993  INTEGER :: actlen, i, k, msglen, num_env_restart, &
1994  off, offset
1995  LOGICAL :: lbf
1996  LOGICAL, DIMENSION(3, 2) :: m
1997  REAL(kind=dp) :: bf, bu
1998  REAL(kind=dp), DIMENSION(3, 2) :: bg, cg, f, ig
1999  REAL(kind=dp), DIMENSION(:), POINTER :: message
2000  TYPE(cp_logger_type), POINTER :: logger
2001 
2002  NULLIFY (logger)
2003  logger => cp_get_default_logger()
2004 
2005  ! assign the pointer to the memory location of the input structure
2006  ! where the RNG state is stored
2007  NULLIFY (message)
2008  CALL section_vals_val_get(helium_env(1)%helium%input, &
2009  "MOTION%PINT%HELIUM%RNG_STATE%_DEFAULT_KEYWORD_", &
2010  r_vals=message)
2011 
2012  ! check the number of environments presumably stored in the restart
2013  actlen = SIZE(message)
2014  num_env_restart = actlen/40
2015 
2016  ! check, if RNG restart has the same dimension as helium%num_env
2017  IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
2018  err_str = "Reading RNG state from the input file."
2019  CALL helium_write_line(err_str)
2020  err_str = "Number of environments in the restart...: '"
2021  stmp = ""
2022  WRITE (stmp, *) num_env_restart
2023  err_str = trim(adjustl(err_str))//trim(adjustl(stmp))//"'."
2024  CALL helium_write_line(err_str)
2025  err_str = "Number of current run time environments.: '"
2026  stmp = ""
2027  WRITE (stmp, *) helium_env(1)%helium%num_env
2028  err_str = trim(adjustl(err_str))//trim(adjustl(stmp))//"'."
2029  CALL helium_write_line(err_str)
2030  err_str = "Missmatch between number of RNG states in input file and helium environments."
2031  cpabort(err_str)
2032  ELSE
2033  CALL helium_write_line("RNG state read from the input file.")
2034 
2035  ! unpack the buffer at each processor, set RNG state
2036  offset = 0
2037  DO i = 1, logger%para_env%mepos
2038  offset = offset + helium_env(1)%env_all(i)
2039  END DO
2040 
2041  DO k = 1, SIZE(helium_env)
2042  msglen = 40
2043  off = msglen*(offset + k - 1)
2044  m(:, :) = .true.
2045  f(:, :) = 0.0_dp
2046  bg(:, :) = unpack(message(off + 1:off + 6), mask=m, field=f)
2047  cg(:, :) = unpack(message(off + 7:off + 12), mask=m, field=f)
2048  ig(:, :) = unpack(message(off + 13:off + 18), mask=m, field=f)
2049  bf = message(off + 19)
2050  bu = message(off + 20)
2051  IF (bf .GT. 0) THEN
2052  lbf = .true.
2053  ELSE
2054  lbf = .false.
2055  END IF
2056  CALL helium_env(k)%helium%rng_stream_uniform%set(bg=bg, cg=cg, ig=ig, &
2057  buffer=bu, buffer_filled=lbf)
2058  bg(:, :) = unpack(message(off + 21:off + 26), mask=m, field=f)
2059  cg(:, :) = unpack(message(off + 27:off + 32), mask=m, field=f)
2060  ig(:, :) = unpack(message(off + 33:off + 38), mask=m, field=f)
2061  bf = message(off + 39)
2062  bu = message(off + 40)
2063  IF (bf .GT. 0) THEN
2064  lbf = .true.
2065  ELSE
2066  lbf = .false.
2067  END IF
2068  CALL helium_env(k)%helium%rng_stream_gaussian%set(bg=bg, cg=cg, ig=ig, &
2069  buffer=bu, buffer_filled=lbf)
2070  END DO
2071  END IF
2072 
2073  NULLIFY (message)
2074 
2075  RETURN
2076  END SUBROUTINE helium_rng_restore
2077 
2078 ! ***************************************************************************
2079 !> \brief Create the RDF related data structures and initialize
2080 !> \param helium ...
2081 !> \date 2014-02-25
2082 !> \par History
2083 !> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
2084 !> \author Lukasz Walewski
2085 ! **************************************************************************************************
2086  SUBROUTINE helium_rdf_init(helium)
2087 
2088  TYPE(helium_solvent_type), POINTER :: helium
2089 
2090  CHARACTER(len=2*default_string_length) :: err_str, stmp
2091  INTEGER :: ii, ij
2092  LOGICAL :: explicit
2093  TYPE(cp_logger_type), POINTER :: logger
2094 
2095  ! read parameters
2096  NULLIFY (logger)
2097  logger => cp_get_default_logger()
2098  CALL section_vals_val_get( &
2099  helium%input, &
2100  "MOTION%PINT%HELIUM%RDF%SOLUTE_HE", &
2101  l_val=helium%rdf_sol_he)
2102  CALL section_vals_val_get( &
2103  helium%input, &
2104  "MOTION%PINT%HELIUM%RDF%HE_HE", &
2105  l_val=helium%rdf_he_he)
2106 
2107  ! Prevent He_He Rdfs for a single he atom:
2108  IF (helium%atoms <= 1) THEN
2109  helium%rdf_he_he = .false.
2110  END IF
2111 
2112  helium%rdf_num = 0
2113  IF (helium%rdf_sol_he) THEN
2114  IF (helium%solute_present) THEN
2115  ! get number of centers from solute to store solute positions
2116  ALLOCATE (helium%rdf_centers(helium%beads, helium%solute_atoms*3))
2117  helium%rdf_centers(:, :) = 0.0_dp
2118  helium%rdf_num = helium%solute_atoms
2119  ELSE
2120  helium%rdf_sol_he = .false.
2121  END IF
2122  END IF
2123 
2124  IF (helium%rdf_he_he) THEN
2125  helium%rdf_num = helium%rdf_num + 1
2126  END IF
2127 
2128  ! set the flag for RDF and either proceed or return
2129  IF (helium%rdf_num > 0) THEN
2130  helium%rdf_present = .true.
2131  ELSE
2132  helium%rdf_present = .false.
2133  err_str = "HELIUM RDFs requested, but no valid choice of "// &
2134  "parameters specified. No RDF will be computed."
2135  cpwarn(err_str)
2136  RETURN
2137  END IF
2138 
2139  ! set the maximum RDF range
2140  CALL section_vals_val_get( &
2141  helium%input, &
2142  "MOTION%PINT%HELIUM%RDF%MAXR", &
2143  explicit=explicit)
2144  IF (explicit) THEN
2145  ! use the value explicitly set in the input
2146  CALL section_vals_val_get( &
2147  helium%input, &
2148  "MOTION%PINT%HELIUM%RDF%MAXR", &
2149  r_val=helium%rdf_maxr)
2150  ELSE
2151  ! use the default value
2152  CALL section_vals_val_get( &
2153  helium%input, &
2154  "MOTION%PINT%HELIUM%DROPLET_RADIUS", &
2155  explicit=explicit)
2156  IF (explicit) THEN
2157  ! use the droplet radius
2158  IF (helium%solute_present .AND. .NOT. helium%periodic) THEN
2159  ! COM of solute is used as center of the box.
2160  ! Therefore distances became larger then droplet_radius
2161  ! when parts of the solute are on opposite site of
2162  ! COM and helium.
2163  ! Use two times droplet_radius for security:
2164  helium%rdf_maxr = helium%droplet_radius*2.0_dp
2165  ELSE
2166  helium%rdf_maxr = helium%droplet_radius
2167  END IF
2168  ELSE
2169  ! use cell_size and cell_shape
2170  ! (they are set regardless of us being periodic or not)
2171  SELECT CASE (helium%cell_shape)
2172  CASE (helium_cell_shape_cube)
2173  helium%rdf_maxr = helium%cell_size/2.0_dp
2174  CASE (helium_cell_shape_octahedron)
2175  helium%rdf_maxr = helium%cell_size*sqrt(3.0_dp)/4.0_dp
2176  CASE DEFAULT
2177  helium%rdf_maxr = 0.0_dp
2178  WRITE (stmp, *) helium%cell_shape
2179  err_str = "cell shape unknown ("//trim(adjustl(stmp))//")"
2180  cpabort(err_str)
2181  END SELECT
2182  END IF
2183  END IF
2184 
2185  ! get number of bins and set bin spacing
2186  CALL section_vals_val_get( &
2187  helium%input, &
2188  "MOTION%PINT%HELIUM%RDF%NBIN", &
2189  i_val=helium%rdf_nbin)
2190  helium%rdf_delr = helium%rdf_maxr/real(helium%rdf_nbin, dp)
2191 
2192  ! get the weighting factor
2193  CALL section_vals_val_get( &
2194  helium%input, &
2195  "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
2196  i_val=helium%rdf_iweight)
2197 
2198  ! allocate and initialize memory for RDF storage
2199  ii = helium%rdf_num
2200  ij = helium%rdf_nbin
2201  ALLOCATE (helium%rdf_inst(ii, ij))
2202  ALLOCATE (helium%rdf_accu(ii, ij))
2203  ALLOCATE (helium%rdf_rstr(ii, ij))
2204  helium%rdf_inst(:, :) = 0.0_dp
2205  helium%rdf_accu(:, :) = 0.0_dp
2206  helium%rdf_rstr(:, :) = 0.0_dp
2207 
2208  RETURN
2209  END SUBROUTINE helium_rdf_init
2210 
2211 ! ***************************************************************************
2212 !> \brief Restore the RDFs from the input structure
2213 !> \param helium_env ...
2214 !> \date 2011-06-22
2215 !> \par History
2216 !> 2016-07-14 Modified to work with independent helium_env [cschran]
2217 !> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
2218 !> \author Lukasz Walewski
2219 ! **************************************************************************************************
2220  SUBROUTINE helium_rdf_restore(helium_env)
2221 
2222  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
2223 
2224  CHARACTER(len=default_string_length) :: stmp1, stmp2
2225  CHARACTER(len=max_line_length) :: err_str
2226  INTEGER :: ii, ij, itmp, k, msglen
2227  LOGICAL :: explicit, ltmp
2228  LOGICAL, DIMENSION(:, :), POINTER :: m
2229  REAL(kind=dp), DIMENSION(:), POINTER :: message
2230  REAL(kind=dp), DIMENSION(:, :), POINTER :: f
2231 
2232  CALL section_vals_val_get(helium_env(1)%helium%input, &
2233  "MOTION%PINT%HELIUM%AVERAGES%RDF", &
2234  explicit=explicit)
2235  IF (explicit) THEN
2236  NULLIFY (message)
2237  CALL section_vals_val_get(helium_env(1)%helium%input, &
2238  "MOTION%PINT%HELIUM%AVERAGES%RDF", &
2239  r_vals=message)
2240  msglen = SIZE(message)
2241  itmp = SIZE(helium_env(1)%helium%rdf_rstr)
2242  ltmp = (msglen .EQ. itmp)
2243  IF (.NOT. ltmp) THEN
2244  stmp1 = ""
2245  WRITE (stmp1, *) msglen
2246  stmp2 = ""
2247  WRITE (stmp2, *) itmp
2248  err_str = "Size of the RDF array in the input ("// &
2249  trim(adjustl(stmp1))// &
2250  .NE.") that in the runtime ("// &
2251  trim(adjustl(stmp2))//")."
2252  cpabort(err_str)
2253  END IF
2254  ELSE
2255  RETURN
2256  END IF
2257 
2258  ii = helium_env(1)%helium%rdf_num
2259  ij = helium_env(1)%helium%rdf_nbin
2260  NULLIFY (m, f)
2261  ALLOCATE (m(ii, ij))
2262  ALLOCATE (f(ii, ij))
2263  m(:, :) = .true.
2264  f(:, :) = 0.0_dp
2265 
2266  DO k = 1, SIZE(helium_env)
2267  helium_env(k)%helium%rdf_rstr(:, :) = unpack(message(1:msglen), mask=m, field=f)
2268  CALL section_vals_val_get(helium_env(k)%helium%input, &
2269  "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
2270  i_val=helium_env(k)%helium%rdf_iweight)
2271  END DO
2272 
2273  DEALLOCATE (f, m)
2274  NULLIFY (message)
2275 
2276  RETURN
2277  END SUBROUTINE helium_rdf_restore
2278 
2279 ! ***************************************************************************
2280 !> \brief Release/deallocate RDF related data structures
2281 !> \param helium ...
2282 !> \date 2014-02-25
2283 !> \author Lukasz Walewski
2284 ! **************************************************************************************************
2285  SUBROUTINE helium_rdf_release(helium)
2286 
2287  TYPE(helium_solvent_type), POINTER :: helium
2288 
2289  IF (helium%rdf_present) THEN
2290 
2291  DEALLOCATE ( &
2292  helium%rdf_rstr, &
2293  helium%rdf_accu, &
2294  helium%rdf_inst)
2295 
2296  NULLIFY ( &
2297  helium%rdf_rstr, &
2298  helium%rdf_accu, &
2299  helium%rdf_inst)
2300 
2301  IF (helium%rdf_sol_he) THEN
2302  DEALLOCATE (helium%rdf_centers)
2303  NULLIFY (helium%rdf_centers)
2304  END IF
2305 
2306  END IF
2307 
2308  RETURN
2309  END SUBROUTINE helium_rdf_release
2310 
2311 ! ***************************************************************************
2312 !> \brief Check whether property <prop> is activated in the input structure
2313 !> \param helium ...
2314 !> \param prop ...
2315 !> \return ...
2316 !> \date 2014-06-26
2317 !> \author Lukasz Walewski
2318 !> \note The property is controlled by two items in the input structure,
2319 !> the printkey and the control section. Two settings result in
2320 !> the property being considered active:
2321 !> 1. printkey is on at the given print level
2322 !> 2. control section is explicit and on
2323 !> If the property is considered active it should be calculated
2324 !> and printed through out the run.
2325 ! **************************************************************************************************
2326  FUNCTION helium_property_active(helium, prop) RESULT(is_active)
2327 
2328  TYPE(helium_solvent_type), POINTER :: helium
2329  CHARACTER(len=*) :: prop
2330  LOGICAL :: is_active
2331 
2332  CHARACTER(len=default_string_length) :: input_path
2333  INTEGER :: print_level
2334  LOGICAL :: explicit, is_on
2335  TYPE(cp_logger_type), POINTER :: logger
2336  TYPE(section_vals_type), POINTER :: print_key, section
2337 
2338  is_active = .false.
2339  NULLIFY (logger)
2340  logger => cp_get_default_logger()
2341 
2342  ! if the printkey is on at this runlevel consider prop to be active
2343  NULLIFY (print_key)
2344  input_path = "MOTION%PINT%HELIUM%PRINT%"//trim(adjustl(prop))
2345  print_key => section_vals_get_subs_vals( &
2346  helium%input, &
2347  input_path)
2348  is_on = cp_printkey_is_on( &
2349  iteration_info=logger%iter_info, &
2350  print_key=print_key)
2351  IF (is_on) THEN
2352  is_active = .true.
2353  RETURN
2354  END IF
2355 
2356  ! if the control section is explicit and on consider prop to be active
2357  ! and activate the printkey
2358  is_active = .false.
2359  NULLIFY (section)
2360  input_path = "MOTION%PINT%HELIUM%"//trim(adjustl(prop))
2361  section => section_vals_get_subs_vals( &
2362  helium%input, &
2363  input_path)
2364  CALL section_vals_get(section, explicit=explicit)
2365  IF (explicit) THEN
2366  ! control section explicitly present, continue checking
2367  CALL section_vals_val_get( &
2368  section, &
2369  "_SECTION_PARAMETERS_", &
2370  l_val=is_on)
2371  IF (is_on) THEN
2372  ! control section is explicit and on, activate the property
2373  is_active = .true.
2374  ! activate the corresponding print_level as well
2375  print_level = logger%iter_info%print_level
2376  CALL section_vals_val_set( &
2377  print_key, &
2378  "_SECTION_PARAMETERS_", &
2379  i_val=print_level)
2380  END IF
2381  END IF
2382 
2383  RETURN
2384  END FUNCTION helium_property_active
2385 
2386 ! ***************************************************************************
2387 !> \brief Create the density related data structures and initialize
2388 !> \param helium ...
2389 !> \date 2014-02-25
2390 !> \author Lukasz Walewski
2391 ! **************************************************************************************************
2392  SUBROUTINE helium_rho_property_init(helium)
2393 
2394  TYPE(helium_solvent_type), POINTER :: helium
2395 
2396  INTEGER :: nc
2397 
2398  ALLOCATE (helium%rho_property(rho_num))
2399 
2400  helium%rho_property(rho_atom_number)%name = 'Atom number density'
2401  nc = 1
2402  helium%rho_property(rho_atom_number)%num_components = nc
2403  ALLOCATE (helium%rho_property(rho_atom_number)%filename_suffix(nc))
2404  ALLOCATE (helium%rho_property(rho_atom_number)%component_name(nc))
2405  ALLOCATE (helium%rho_property(rho_atom_number)%component_index(nc))
2406  helium%rho_property(rho_atom_number)%filename_suffix(1) = 'an'
2407  helium%rho_property(rho_atom_number)%component_name(1) = ''
2408  helium%rho_property(rho_atom_number)%component_index(:) = 0
2409 
2410  helium%rho_property(rho_projected_area)%name = 'Projected area squared density, A*A(r)'
2411  nc = 3
2412  helium%rho_property(rho_projected_area)%num_components = nc
2413  ALLOCATE (helium%rho_property(rho_projected_area)%filename_suffix(nc))
2414  ALLOCATE (helium%rho_property(rho_projected_area)%component_name(nc))
2415  ALLOCATE (helium%rho_property(rho_projected_area)%component_index(nc))
2416  helium%rho_property(rho_projected_area)%filename_suffix(1) = 'pa_x'
2417  helium%rho_property(rho_projected_area)%filename_suffix(2) = 'pa_y'
2418  helium%rho_property(rho_projected_area)%filename_suffix(3) = 'pa_z'
2419  helium%rho_property(rho_projected_area)%component_name(1) = 'component x'
2420  helium%rho_property(rho_projected_area)%component_name(2) = 'component y'
2421  helium%rho_property(rho_projected_area)%component_name(3) = 'component z'
2422  helium%rho_property(rho_projected_area)%component_index(:) = 0
2423 
2424  helium%rho_property(rho_winding_number)%name = 'Winding number squared density, W*W(r)'
2425  nc = 3
2426  helium%rho_property(rho_winding_number)%num_components = nc
2427  ALLOCATE (helium%rho_property(rho_winding_number)%filename_suffix(nc))
2428  ALLOCATE (helium%rho_property(rho_winding_number)%component_name(nc))
2429  ALLOCATE (helium%rho_property(rho_winding_number)%component_index(nc))
2430  helium%rho_property(rho_winding_number)%filename_suffix(1) = 'wn_x'
2431  helium%rho_property(rho_winding_number)%filename_suffix(2) = 'wn_y'
2432  helium%rho_property(rho_winding_number)%filename_suffix(3) = 'wn_z'
2433  helium%rho_property(rho_winding_number)%component_name(1) = 'component x'
2434  helium%rho_property(rho_winding_number)%component_name(2) = 'component y'
2435  helium%rho_property(rho_winding_number)%component_name(3) = 'component z'
2436  helium%rho_property(rho_winding_number)%component_index(:) = 0
2437 
2438  helium%rho_property(rho_winding_cycle)%name = 'Winding number squared density, W^2(r)'
2439  nc = 3
2440  helium%rho_property(rho_winding_cycle)%num_components = nc
2441  ALLOCATE (helium%rho_property(rho_winding_cycle)%filename_suffix(nc))
2442  ALLOCATE (helium%rho_property(rho_winding_cycle)%component_name(nc))
2443  ALLOCATE (helium%rho_property(rho_winding_cycle)%component_index(nc))
2444  helium%rho_property(rho_winding_cycle)%filename_suffix(1) = 'wc_x'
2445  helium%rho_property(rho_winding_cycle)%filename_suffix(2) = 'wc_y'
2446  helium%rho_property(rho_winding_cycle)%filename_suffix(3) = 'wc_z'
2447  helium%rho_property(rho_winding_cycle)%component_name(1) = 'component x'
2448  helium%rho_property(rho_winding_cycle)%component_name(2) = 'component y'
2449  helium%rho_property(rho_winding_cycle)%component_name(3) = 'component z'
2450  helium%rho_property(rho_winding_cycle)%component_index(:) = 0
2451 
2452  helium%rho_property(rho_moment_of_inertia)%name = 'Moment of inertia'
2453  nc = 3
2454  helium%rho_property(rho_moment_of_inertia)%num_components = nc
2455  ALLOCATE (helium%rho_property(rho_moment_of_inertia)%filename_suffix(nc))
2456  ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_name(nc))
2457  ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_index(nc))
2458  helium%rho_property(rho_moment_of_inertia)%filename_suffix(1) = 'mi_x'
2459  helium%rho_property(rho_moment_of_inertia)%filename_suffix(2) = 'mi_y'
2460  helium%rho_property(rho_moment_of_inertia)%filename_suffix(3) = 'mi_z'
2461  helium%rho_property(rho_moment_of_inertia)%component_name(1) = 'component x'
2462  helium%rho_property(rho_moment_of_inertia)%component_name(2) = 'component y'
2463  helium%rho_property(rho_moment_of_inertia)%component_name(3) = 'component z'
2464  helium%rho_property(rho_moment_of_inertia)%component_index(:) = 0
2465 
2466  helium%rho_property(:)%is_calculated = .false.
2467 
2468  RETURN
2469  END SUBROUTINE helium_rho_property_init
2470 
2471 ! ***************************************************************************
2472 !> \brief Create the density related data structures and initialize
2473 !> \param helium ...
2474 !> \date 2014-02-25
2475 !> \author Lukasz Walewski
2476 ! **************************************************************************************************
2477  SUBROUTINE helium_rho_init(helium)
2478 
2479  TYPE(helium_solvent_type), POINTER :: helium
2480 
2481  INTEGER :: ii, itmp, jtmp
2482  LOGICAL :: explicit, ltmp
2483 
2484  CALL helium_rho_property_init(helium)
2485 
2486  helium%rho_num_act = 0
2487 
2488  ! check for atom number density
2489  CALL section_vals_val_get( &
2490  helium%input, &
2491  "MOTION%PINT%HELIUM%RHO%ATOM_NUMBER", &
2492  l_val=ltmp)
2493  IF (ltmp) THEN
2494  helium%rho_property(rho_atom_number)%is_calculated = .true.
2495  helium%rho_num_act = helium%rho_num_act + 1
2496  helium%rho_property(rho_atom_number)%component_index(1) = helium%rho_num_act
2497  END IF
2498 
2499  ! check for projected area density
2500  CALL section_vals_val_get( &
2501  helium%input, &
2502  "MOTION%PINT%HELIUM%RHO%PROJECTED_AREA_2", &
2503  l_val=ltmp)
2504  IF (ltmp) THEN
2505  helium%rho_property(rho_projected_area)%is_calculated = .true.
2506  DO ii = 1, helium%rho_property(rho_projected_area)%num_components
2507  helium%rho_num_act = helium%rho_num_act + 1
2508  helium%rho_property(rho_projected_area)%component_index(ii) = helium%rho_num_act
2509  END DO
2510  END IF
2511 
2512  ! check for winding number density
2513  CALL section_vals_val_get( &
2514  helium%input, &
2515  "MOTION%PINT%HELIUM%RHO%WINDING_NUMBER_2", &
2516  l_val=ltmp)
2517  IF (ltmp) THEN
2518  helium%rho_property(rho_winding_number)%is_calculated = .true.
2519  DO ii = 1, helium%rho_property(rho_winding_number)%num_components
2520  helium%rho_num_act = helium%rho_num_act + 1
2521  helium%rho_property(rho_winding_number)%component_index(ii) = helium%rho_num_act
2522  END DO
2523  END IF
2524 
2525  ! check for winding cycle density
2526  CALL section_vals_val_get( &
2527  helium%input, &
2528  "MOTION%PINT%HELIUM%RHO%WINDING_CYCLE_2", &
2529  l_val=ltmp)
2530  IF (ltmp) THEN
2531  helium%rho_property(rho_winding_cycle)%is_calculated = .true.
2532  DO ii = 1, helium%rho_property(rho_winding_cycle)%num_components
2533  helium%rho_num_act = helium%rho_num_act + 1
2534  helium%rho_property(rho_winding_cycle)%component_index(ii) = helium%rho_num_act
2535  END DO
2536  END IF
2537 
2538  ! check for moment of inertia density
2539  CALL section_vals_val_get( &
2540  helium%input, &
2541  "MOTION%PINT%HELIUM%RHO%MOMENT_OF_INERTIA", &
2542  l_val=ltmp)
2543  IF (ltmp) THEN
2544  helium%rho_property(rho_moment_of_inertia)%is_calculated = .true.
2545  DO ii = 1, helium%rho_property(rho_moment_of_inertia)%num_components
2546  helium%rho_num_act = helium%rho_num_act + 1
2547  helium%rho_property(rho_moment_of_inertia)%component_index(ii) = helium%rho_num_act
2548  END DO
2549  END IF
2550 
2551  ! set the cube dimensions, etc (common to all estimators)
2552  helium%rho_maxr = helium%cell_size
2553  CALL section_vals_val_get( &
2554  helium%input, &
2555  "MOTION%PINT%HELIUM%RHO%NBIN", &
2556  i_val=helium%rho_nbin)
2557  helium%rho_delr = helium%rho_maxr/real(helium%rho_nbin, dp)
2558 
2559  ! check for optional estimators based on winding paths
2560  helium%rho_num_min_len_wdg = 0
2561  CALL section_vals_val_get( &
2562  helium%input, &
2563  "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_WDG", &
2564  explicit=explicit)
2565  IF (explicit) THEN
2566  NULLIFY (helium%rho_min_len_wdg_vals)
2567  CALL section_vals_val_get( &
2568  helium%input, &
2569  "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_WDG", &
2570  i_vals=helium%rho_min_len_wdg_vals)
2571  itmp = SIZE(helium%rho_min_len_wdg_vals)
2572  IF (itmp .GT. 0) THEN
2573  helium%rho_num_min_len_wdg = itmp
2574  helium%rho_num_act = helium%rho_num_act + itmp
2575  END IF
2576  END IF
2577 
2578  ! check for optional estimators based on non-winding paths
2579  helium%rho_num_min_len_non = 0
2580  CALL section_vals_val_get( &
2581  helium%input, &
2582  "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_NON", &
2583  explicit=explicit)
2584  IF (explicit) THEN
2585  NULLIFY (helium%rho_min_len_non_vals)
2586  CALL section_vals_val_get( &
2587  helium%input, &
2588  "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_NON", &
2589  i_vals=helium%rho_min_len_non_vals)
2590  itmp = SIZE(helium%rho_min_len_non_vals)
2591  IF (itmp .GT. 0) THEN
2592  helium%rho_num_min_len_non = itmp
2593  helium%rho_num_act = helium%rho_num_act + itmp
2594  END IF
2595  END IF
2596 
2597  ! check for optional estimators based on all paths
2598  helium%rho_num_min_len_all = 0
2599  CALL section_vals_val_get( &
2600  helium%input, &
2601  "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_ALL", &
2602  explicit=explicit)
2603  IF (explicit) THEN
2604  NULLIFY (helium%rho_min_len_all_vals)
2605  CALL section_vals_val_get( &
2606  helium%input, &
2607  "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_ALL", &
2608  i_vals=helium%rho_min_len_all_vals)
2609  itmp = SIZE(helium%rho_min_len_all_vals)
2610  IF (itmp .GT. 0) THEN
2611  helium%rho_num_min_len_all = itmp
2612  helium%rho_num_act = helium%rho_num_act + itmp
2613  END IF
2614  END IF
2615 
2616  ! get the weighting factor
2617  CALL section_vals_val_get( &
2618  helium%input, &
2619  "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
2620  i_val=helium%rho_iweight)
2621 
2622  ! allocate and initialize memory for density storage
2623  itmp = helium%rho_nbin
2624  jtmp = helium%rho_num_act
2625  ALLOCATE (helium%rho_inst(jtmp, itmp, itmp, itmp))
2626  ALLOCATE (helium%rho_accu(jtmp, itmp, itmp, itmp))
2627  ALLOCATE (helium%rho_rstr(jtmp, itmp, itmp, itmp))
2628  ALLOCATE (helium%rho_incr(jtmp, helium%atoms, helium%beads))
2629 
2630  helium%rho_incr(:, :, :) = 0.0_dp
2631  helium%rho_inst(:, :, :, :) = 0.0_dp
2632  helium%rho_accu(:, :, :, :) = 0.0_dp
2633  helium%rho_rstr(:, :, :, :) = 0.0_dp
2634 
2635  RETURN
2636  END SUBROUTINE helium_rho_init
2637 
2638 ! ***************************************************************************
2639 !> \brief Restore the densities from the input structure.
2640 !> \param helium_env ...
2641 !> \date 2011-06-22
2642 !> \par History
2643 !> 2016-07-14 Modified to work with independent helium_env [cschran]
2644 !> \author Lukasz Walewski
2645 ! **************************************************************************************************
2646  SUBROUTINE helium_rho_restore(helium_env)
2647 
2648  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
2649 
2650  CHARACTER(len=default_string_length) :: stmp1, stmp2
2651  CHARACTER(len=max_line_length) :: err_str
2652  INTEGER :: itmp, k, msglen
2653  LOGICAL :: explicit, ltmp
2654  LOGICAL, DIMENSION(:, :, :, :), POINTER :: m
2655  REAL(kind=dp), DIMENSION(:), POINTER :: message
2656  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: f
2657 
2658  CALL section_vals_val_get(helium_env(1)%helium%input, &
2659  "MOTION%PINT%HELIUM%AVERAGES%RHO", &
2660  explicit=explicit)
2661  IF (explicit) THEN
2662  NULLIFY (message)
2663  CALL section_vals_val_get(helium_env(1)%helium%input, &
2664  "MOTION%PINT%HELIUM%AVERAGES%RHO", &
2665  r_vals=message)
2666  msglen = SIZE(message)
2667  itmp = SIZE(helium_env(1)%helium%rho_rstr)
2668  ltmp = (msglen .EQ. itmp)
2669  IF (.NOT. ltmp) THEN
2670  stmp1 = ""
2671  WRITE (stmp1, *) msglen
2672  stmp2 = ""
2673  WRITE (stmp2, *) itmp
2674  err_str = "Size of the S density array in the input ("// &
2675  trim(adjustl(stmp1))// &
2676  .NE.") that in the runtime ("// &
2677  trim(adjustl(stmp2))//")."
2678  cpabort(err_str)
2679  END IF
2680  ELSE
2681  RETURN
2682  END IF
2683 
2684  itmp = helium_env(1)%helium%rho_nbin
2685  NULLIFY (m, f)
2686  ALLOCATE (m(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp))
2687  ALLOCATE (f(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp))
2688  m(:, :, :, :) = .true.
2689  f(:, :, :, :) = 0.0_dp
2690 
2691  DO k = 1, SIZE(helium_env)
2692  helium_env(k)%helium%rho_rstr(:, :, :, :) = unpack(message(1:msglen), mask=m, field=f)
2693  END DO
2694 
2695  DEALLOCATE (f, m)
2696  NULLIFY (message)
2697 
2698  RETURN
2699  END SUBROUTINE helium_rho_restore
2700 
2701 ! ***************************************************************************
2702 !> \brief Count atoms of different types and store their global indices.
2703 !> \param helium ...
2704 !> \param pint_env ...
2705 !> \author Lukasz Walewski
2706 !> \note Arrays ALLOCATEd here are (should be) DEALLOCATEd in
2707 !> helium_release.
2708 ! **************************************************************************************************
2709  SUBROUTINE helium_set_solute_indices(helium, pint_env)
2710  TYPE(helium_solvent_type), POINTER :: helium
2711  TYPE(pint_env_type), INTENT(IN) :: pint_env
2712 
2713  INTEGER :: iatom, natoms
2714  REAL(kind=dp) :: mass
2715  TYPE(cp_subsys_type), POINTER :: my_subsys
2716  TYPE(f_env_type), POINTER :: my_f_env
2717  TYPE(particle_list_type), POINTER :: my_particles
2718 
2719 ! set up my_particles structure
2720 
2721  NULLIFY (my_f_env, my_subsys, my_particles)
2722  CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
2723  f_env=my_f_env)
2724  CALL force_env_get(force_env=my_f_env%force_env, subsys=my_subsys)
2725  CALL cp_subsys_get(my_subsys, particles=my_particles)
2726  CALL f_env_rm_defaults(my_f_env)
2727 
2728  natoms = helium%solute_atoms
2729  NULLIFY (helium%solute_element)
2730  ALLOCATE (helium%solute_element(natoms))
2731 
2732  ! find out how many different atomic types are there
2733  helium%enum = 0
2734  DO iatom = 1, natoms
2735  CALL get_atomic_kind(my_particles%els(iatom)%atomic_kind, &
2736  mass=mass, &
2737  element_symbol=helium%solute_element(iatom))
2738  END DO
2739 
2740  RETURN
2741  END SUBROUTINE helium_set_solute_indices
2742 
2743 ! ***************************************************************************
2744 !> \brief Sets helium%solute_cell based on the solute's force_env.
2745 !> \param helium ...
2746 !> \param pint_env ...
2747 !> \author Lukasz Walewski
2748 !> \note The simulation cell for the solvated molecule is taken from force_env
2749 !> which should assure that we get proper cell dimensions regardless of
2750 !> the method used for the solute (QS, FIST). Helium solvent needs the
2751 !> solute's cell dimensions to calculate the solute-solvent distances
2752 !> correctly.
2753 !> \note At the moment only orthorhombic cells are supported.
2754 ! **************************************************************************************************
2755  SUBROUTINE helium_set_solute_cell(helium, pint_env)
2756  TYPE(helium_solvent_type), POINTER :: helium
2757  TYPE(pint_env_type), INTENT(IN) :: pint_env
2758 
2759  LOGICAL :: my_orthorhombic
2760  TYPE(cell_type), POINTER :: my_cell
2761  TYPE(f_env_type), POINTER :: my_f_env
2762 
2763 ! get the cell structure from pint_env
2764 
2765  NULLIFY (my_f_env, my_cell)
2766  CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
2767  f_env=my_f_env)
2768  CALL force_env_get(force_env=my_f_env%force_env, cell=my_cell)
2769  CALL f_env_rm_defaults(my_f_env)
2770 
2771  CALL get_cell(my_cell, orthorhombic=my_orthorhombic)
2772  IF (.NOT. my_orthorhombic) THEN
2773  cpabort("Helium solvent not implemented for non-orthorhombic cells.")
2774  ELSE
2775  helium%solute_cell => my_cell
2776  END IF
2777 
2778  RETURN
2779  END SUBROUTINE helium_set_solute_cell
2780 
2781 END MODULE helium_methods
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::walewski2014
integer, save, public walewski2014
Definition: bibliography.F:43
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_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
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_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition: cp_log_handling.F:213
cp_log_handling::cp_logger_release
subroutine, public cp_logger_release(logger)
releases this logger
Definition: cp_log_handling.F:423
cp_log_handling::cp_logger_create
subroutine, public cp_logger_create(logger, para_env, print_level, default_global_unit_nr, default_local_unit_nr, global_filename, local_filename, close_global_unit_on_dealloc, iter_info, close_local_unit_on_dealloc, suffix, template_logger)
initializes a logger
Definition: cp_log_handling.F:285
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition: cp_log_handling.F:188
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_printkey_is_on
logical function, public cp_printkey_is_on(iteration_info, print_key)
returns true if the printlevel activates this printkey does not look if this iteration it should be p...
Definition: cp_output_handling.F:429
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
f77_interface
interface to use cp2k as library
Definition: f77_interface.F:20
f77_interface::f_env_add_defaults
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
Definition: f77_interface.F:450
f77_interface::f_env_rm_defaults
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
Definition: f77_interface.F:494
force_env_types
Interface for the force calculations.
Definition: force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
Definition: force_env_types.F:329
helium_common
Independent helium subroutines shared with other modules.
Definition: helium_common.F:14
helium_common::helium_com
pure real(kind=dp) function, dimension(3), public helium_com(helium)
Calculate center of mass.
Definition: helium_common.F:2387
helium_common::helium_pbc
subroutine, public helium_pbc(helium, r, enforce)
General PBC routine for helium.
Definition: helium_common.F:81
helium_interactions
Methods that handle helium-solvent and helium-helium interactions.
Definition: helium_interactions.F:13
helium_interactions::helium_vij
elemental real(kind=dp) function, public helium_vij(r)
Helium-helium pair interaction potential.
Definition: helium_interactions.F:834
helium_io
I/O subroutines for helium.
Definition: helium_io.F:13
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_methods
Methods dealing with helium_solvent_type.
Definition: helium_methods.F:13
helium_methods::helium_release
subroutine, public helium_release(helium_env)
Releases helium_solvent_type.
Definition: helium_methods.F:856
helium_methods::helium_create
subroutine, public helium_create(helium_env, input, solute)
Data-structure that holds all needed information about (superfluid) helium solvent.
Definition: helium_methods.F:113
helium_methods::helium_init
subroutine, public helium_init(helium_env, pint_env)
Initialize helium data structures.
Definition: helium_methods.F:1056
helium_nnp
Methods dealing with Neural Network interaction potential.
Definition: helium_nnp.F:13
helium_nnp::helium_init_nnp
subroutine, public helium_init_nnp(helium, nnp, input)
Read and initialize all the information for neural network potentials.
Definition: helium_nnp.F:68
helium_sampling
Methods for sampling helium variables.
Definition: helium_sampling.F:13
helium_sampling::helium_sample
subroutine, public helium_sample(helium_env, pint_env)
Sample the helium phase space.
Definition: helium_sampling.F:183
helium_types
Data types representing superfluid helium.
Definition: helium_types.F:15
helium_types::rho_moment_of_inertia
integer, parameter, public rho_moment_of_inertia
Definition: helium_types.F:53
helium_types::rho_winding_number
integer, parameter, public rho_winding_number
Definition: helium_types.F:53
helium_types::rho_atom_number
integer, parameter, public rho_atom_number
density function identifier names
Definition: helium_types.F:53
helium_types::rho_winding_cycle
integer, parameter, public rho_winding_cycle
Definition: helium_types.F:53
helium_types::rho_projected_area
integer, parameter, public rho_projected_area
Definition: helium_types.F:53
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_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::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::helium_solute_intpot_nnp
integer, parameter, public helium_solute_intpot_nnp
Definition: input_constants.F:322
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_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition: input_section_types.F:1216
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_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
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::max_line_length
integer, parameter, public max_line_length
Definition: kinds.F:59
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
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
nnp_environment_types
Data types for neural network potentials.
Definition: nnp_environment_types.F:13
nnp_environment_types::nnp_env_release
subroutine, public nnp_env_release(nnp_env)
Release data structure that holds all the information for neural network potentials.
Definition: nnp_environment_types.F:276
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition: parallel_rng_types.F:51
parallel_rng_types::uniform
integer, parameter, public uniform
Definition: parallel_rng_types.F:73
parallel_rng_types::gaussian
integer, parameter, public gaussian
Definition: parallel_rng_types.F:73
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::boltzmann
real(kind=dp), parameter, public boltzmann
Definition: physcon.F:129
physcon::a_mass
real(kind=dp), parameter, public a_mass
Definition: physcon.F:132
physcon::kelvin
real(kind=dp), parameter, public kelvin
Definition: physcon.F:165
physcon::h_bar
real(kind=dp), parameter, public h_bar
Definition: physcon.F:103
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_public
Public path integral routines that can be called from other modules.
Definition: pint_public.F:15
pint_public::pint_com_pos
pure real(kind=dp) function, dimension(3), public pint_com_pos(pint_env)
Return the center of mass of the PI system.
Definition: pint_public.F:45
pint_types
Definition: pint_types.F:8
splines_methods
routines for handling splines
Definition: splines_methods.F:14
splines_methods::init_splinexy
pure subroutine, public init_splinexy(spl, nn)
allocates storage for function table to be interpolated both x and y are allocated
Definition: splines_methods.F:48
splines_methods::init_spline
pure subroutine, public init_spline(spl, dx, y1a, y1b)
allocates storage for y2 table calculates y2 table and other spline parameters
Definition: splines_methods.F:100
splines_types
routines for handling splines_types
Definition: splines_types.F:14
splines_types::spline_data_create
subroutine, public spline_data_create(spline_data)
Data-structure that constains spline table.
Definition: splines_types.F:287
splines_types::spline_data_release
subroutine, public spline_data_release(spline_data)
releases spline_data
Definition: splines_types.F:106