(git:ed6f26b)
Loading...
Searching...
No Matches
helium_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Methods dealing with helium_solvent_type
10!> \author Lukasz Walewski
11!> \date 2009-06-10
12! **************************************************************************************************
13MODULE 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
99CONTAINS
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 cpwarn_if(.NOT. explicit, 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
2781END 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:21
bibliography::walewski2014
integer, save, public walewski2014
Definition bibliography.F:36
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:452
f77_interface::get_cell
subroutine, public get_cell(env_id, cell, per, ierr)
gets a cell
Definition f77_interface.F:1053
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:496
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, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
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:348
input_constants::helium_cell_shape_octahedron
integer, parameter, public helium_cell_shape_octahedron
Definition input_constants.F:318
input_constants::helium_solute_intpot_none
integer, parameter, public helium_solute_intpot_none
Definition input_constants.F:323
input_constants::helium_sampling_worm
integer, parameter, public helium_sampling_worm
Definition input_constants.F:348
input_constants::helium_cell_shape_cube
integer, parameter, public helium_cell_shape_cube
Definition input_constants.F:318
input_constants::helium_solute_intpot_nnp
integer, parameter, public helium_solute_intpot_nnp
Definition input_constants.F:323
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:1223
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:731
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:704
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:1047
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
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
f77_interface::f_env_type
Definition f77_interface.F:163
helium_types::helium_solvent_p_type
data structure for array of solvent helium environments
Definition helium_types.F:300
helium_types::helium_solvent_type
data structure for solvent helium
Definition helium_types.F:93
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
parallel_rng_types::rng_stream_p_type
Definition parallel_rng_types.F:195
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
pint_types::pint_env_type
environment for a path integral run
Definition pint_types.F:112