Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods dealing with helium_solvent_type
10 : !> \author Lukasz Walewski
11 : !> \date 2009-06-10
12 : ! **************************************************************************************************
13 : MODULE helium_methods
14 :
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE bibliography, ONLY: Walewski2014,&
17 : cite_reference
18 : USE cell_types, ONLY: cell_type,&
19 : get_cell
20 : USE cp_files, ONLY: close_file,&
21 : open_file
22 : USE cp_log_handling, ONLY: cp_add_default_logger,&
23 : cp_get_default_logger,&
24 : cp_logger_create,&
25 : cp_logger_release,&
26 : cp_logger_type,&
27 : cp_rm_default_logger
28 : USE cp_output_handling, ONLY: cp_printkey_is_on
29 : USE cp_subsys_types, ONLY: cp_subsys_get,&
30 : cp_subsys_type
31 : USE f77_interface, ONLY: f_env_add_defaults,&
32 : f_env_rm_defaults,&
33 : f_env_type
34 : USE force_env_types, ONLY: force_env_get
35 : USE helium_common, ONLY: helium_com,&
36 : helium_pbc
37 : USE helium_interactions, ONLY: helium_vij
38 : USE helium_io, ONLY: helium_write_line,&
39 : helium_write_setup
40 : USE helium_nnp, ONLY: helium_init_nnp
41 : USE helium_sampling, ONLY: helium_sample
42 : USE helium_types, ONLY: helium_solvent_p_type,&
43 : helium_solvent_type,&
44 : rho_atom_number,&
45 : rho_moment_of_inertia,&
46 : rho_num,&
47 : rho_projected_area,&
48 : rho_winding_cycle,&
49 : rho_winding_number
50 : USE input_constants, ONLY: helium_cell_shape_cube,&
51 : helium_cell_shape_octahedron,&
52 : helium_sampling_ceperley,&
53 : helium_sampling_worm,&
54 : helium_solute_intpot_nnp,&
55 : helium_solute_intpot_none
56 : USE input_section_types, ONLY: section_vals_get,&
57 : section_vals_get_subs_vals,&
58 : section_vals_type,&
59 : section_vals_val_get,&
60 : section_vals_val_set
61 : USE kinds, ONLY: default_path_length,&
62 : default_string_length,&
63 : dp,&
64 : max_line_length
65 : USE mathconstants, ONLY: pi,&
66 : twopi
67 : USE message_passing, ONLY: mp_para_env_type
68 : USE nnp_environment_types, ONLY: nnp_env_release
69 : USE parallel_rng_types, ONLY: GAUSSIAN,&
70 : UNIFORM,&
71 : rng_stream_p_type,&
72 : rng_stream_type
73 : USE particle_list_types, ONLY: particle_list_type
74 : USE physcon, ONLY: a_mass,&
75 : angstrom,&
76 : boltzmann,&
77 : h_bar,&
78 : kelvin,&
79 : massunit
80 : USE pint_public, ONLY: pint_com_pos
81 : USE pint_types, ONLY: pint_env_type
82 : USE splines_methods, ONLY: init_spline,&
83 : init_splinexy
84 : USE splines_types, ONLY: spline_data_create,&
85 : spline_data_release
86 : #include "../base/base_uses.f90"
87 :
88 : IMPLICIT NONE
89 :
90 : PRIVATE
91 :
92 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
93 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_methods'
94 :
95 : PUBLIC :: helium_create
96 : PUBLIC :: helium_init
97 : PUBLIC :: helium_release
98 :
99 : CONTAINS
100 :
101 : ! ***************************************************************************
102 : !> \brief Data-structure that holds all needed information about
103 : !> (superfluid) helium solvent
104 : !> \param helium_env ...
105 : !> \param input ...
106 : !> \param solute ...
107 : !> \par History
108 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
109 : !> 2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
110 : !> \author hforbert
111 : ! **************************************************************************************************
112 26 : 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 26 : 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 26 : 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 26 : CALL timeset(routineN, handle)
134 :
135 26 : CALL cite_reference(Walewski2014)
136 26 : NULLIFY (logger)
137 26 : logger => cp_get_default_logger()
138 :
139 26 : NULLIFY (helium_section)
140 : helium_section => section_vals_get_subs_vals(input, &
141 26 : "MOTION%PINT%HELIUM")
142 26 : CALL section_vals_get(helium_section, explicit=explicit)
143 26 : CPASSERT(explicit)
144 :
145 : ! get number of environments
146 : CALL section_vals_val_get(helium_section, "NUM_ENV", &
147 26 : explicit=explicit)
148 26 : IF (explicit) THEN
149 : CALL section_vals_val_get(helium_section, "NUM_ENV", &
150 26 : i_val=num_env)
151 : ELSE
152 0 : num_env = logger%para_env%num_pe
153 : END IF
154 26 : CPASSERT(num_env >= 0)
155 26 : IF (num_env .NE. logger%para_env%num_pe) THEN
156 6 : msg_str = "NUM_ENV not equal to number of processors"
157 6 : CPWARN(msg_str)
158 : END IF
159 : ! calculate number of tasks for each processor
160 : mepos = num_env/logger%para_env%num_pe &
161 26 : + MIN(MOD(num_env, logger%para_env%num_pe)/(logger%para_env%mepos + 1), 1)
162 : ! gather result
163 26 : NULLIFY (env_all)
164 78 : ALLOCATE (env_all(logger%para_env%num_pe))
165 78 : env_all(:) = 0
166 78 : CALL logger%para_env%allgather(mepos, env_all)
167 :
168 : ! create new communicator for processors with helium_env
169 26 : IF (mepos == 0) THEN
170 1 : color_sub = 0
171 : ELSE
172 25 : color_sub = 1
173 : END IF
174 26 : ALLOCATE (new_comm)
175 26 : CALL new_comm%from_split(logger%para_env, color_sub)
176 : ! release new_comm for processors without helium_env
177 26 : IF (mepos == 0) THEN
178 1 : CALL new_comm%free()
179 1 : DEALLOCATE (new_comm)
180 1 : NULLIFY (new_comm)
181 : END IF
182 :
183 26 : NULLIFY (helium_env)
184 26 : IF (mepos .GT. 0) THEN
185 106 : ALLOCATE (helium_env(mepos))
186 56 : DO k = 1, mepos
187 31 : helium_env(k)%comm => new_comm
188 31 : NULLIFY (helium_env(k)%env_all)
189 31 : helium_env(k)%env_all => env_all
190 3627 : ALLOCATE (helium_env(k)%helium)
191 31 : NULLIFY (helium_env(k)%helium%input)
192 31 : helium_env(k)%helium%input => input
193 56 : helium_env(k)%helium%num_env = num_env
194 : END DO
195 : ! RNG state create & init
196 25 : CALL helium_rng_init(helium_env)
197 :
198 56 : DO k = 1, mepos
199 : NULLIFY (helium_env(k)%helium%ptable, &
200 31 : helium_env(k)%helium%permutation, &
201 31 : helium_env(k)%helium%savepermutation, &
202 31 : helium_env(k)%helium%iperm, &
203 31 : helium_env(k)%helium%saveiperm, &
204 31 : helium_env(k)%helium%itmp_atoms_1d, &
205 31 : helium_env(k)%helium%ltmp_atoms_1d, &
206 31 : helium_env(k)%helium%itmp_atoms_np_1d, &
207 31 : helium_env(k)%helium%pos, &
208 31 : helium_env(k)%helium%savepos, &
209 31 : helium_env(k)%helium%work, &
210 31 : helium_env(k)%helium%force_avrg, &
211 31 : helium_env(k)%helium%force_inst, &
212 31 : helium_env(k)%helium%rtmp_3_np_1d, &
213 31 : helium_env(k)%helium%rtmp_p_ndim_1d, &
214 31 : helium_env(k)%helium%rtmp_p_ndim_np_1d, &
215 31 : helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
216 31 : helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
217 31 : helium_env(k)%helium%rtmp_p_ndim_2d, &
218 31 : helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
219 31 : helium_env(k)%helium%tmatrix, helium_env(k)%helium%pmatrix, &
220 31 : helium_env(k)%helium%nmatrix, helium_env(k)%helium%ipmatrix, &
221 31 : helium_env(k)%helium%u0, helium_env(k)%helium%e0, &
222 31 : helium_env(k)%helium%uoffdiag, helium_env(k)%helium%eoffdiag, &
223 31 : helium_env(k)%helium%vij, &
224 31 : helium_env(k)%helium%rdf_inst, &
225 31 : helium_env(k)%helium%plength_avrg, &
226 31 : helium_env(k)%helium%plength_inst, &
227 31 : helium_env(k)%helium%atom_plength, &
228 31 : helium_env(k)%helium%ename &
229 31 : )
230 :
231 31 : helium_env(k)%helium%accepts = 0
232 31 : helium_env(k)%helium%relrot = 0
233 :
234 : ! check if solute is present in our simulation
235 31 : helium_env(k)%helium%solute_present = .FALSE.
236 31 : helium_env(k)%helium%solute_atoms = 0
237 31 : helium_env(k)%helium%solute_beads = 0
238 : CALL section_vals_val_get( &
239 : helium_section, &
240 : "HELIUM_ONLY", &
241 31 : l_val=ltmp)
242 31 : IF (.NOT. ltmp) THEN
243 21 : IF (PRESENT(solute)) THEN
244 21 : helium_env(k)%helium%solute_present = .TRUE.
245 21 : helium_env(k)%helium%solute_atoms = solute%ndim/3
246 21 : 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 31 : i_val=helium_env(k)%helium%beads)
252 : CALL section_vals_val_get(helium_section, "INOROT", &
253 31 : i_val=helium_env(k)%helium%iter_norot)
254 : CALL section_vals_val_get(helium_section, "IROT", &
255 31 : 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 31 : i_val=itmp)
260 31 : helium_env(k)%helium%first_step = itmp
261 : CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
262 31 : explicit=explicit)
263 31 : IF (explicit) THEN
264 : CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
265 0 : i_val=itmp)
266 0 : helium_env(k)%helium%last_step = itmp
267 : helium_env(k)%helium%num_steps = helium_env(k)%helium%last_step &
268 0 : - helium_env(k)%helium%first_step
269 : ELSE
270 : CALL section_vals_val_get(input, "MOTION%PINT%NUM_STEPS", &
271 31 : i_val=itmp)
272 31 : helium_env(k)%helium%num_steps = itmp
273 : helium_env(k)%helium%last_step = helium_env(k)%helium%first_step &
274 31 : + helium_env(k)%helium%num_steps
275 : END IF
276 :
277 : ! boundary conditions
278 : CALL section_vals_val_get(helium_section, "PERIODIC", &
279 31 : l_val=helium_env(k)%helium%periodic)
280 : CALL section_vals_val_get(helium_section, "CELL_SHAPE", &
281 31 : i_val=helium_env(k)%helium%cell_shape)
282 :
283 : CALL section_vals_val_get(helium_section, "DROPLET_RADIUS", &
284 31 : 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 31 : explicit=expl_dens, r_val=helium_env(k)%helium%density)
296 : CALL section_vals_val_get(helium_section, "NATOMS", &
297 31 : explicit=expl_nats, i_val=helium_env(k)%helium%atoms)
298 : CALL section_vals_val_get(helium_section, "CELL_SIZE", &
299 31 : explicit=expl_cell)
300 31 : cgeof = 1.0_dp
301 31 : IF (helium_env(k)%helium%periodic) THEN
302 16 : IF (helium_env(k)%helium%cell_shape .EQ. helium_cell_shape_octahedron) cgeof = 2.0_dp
303 : END IF
304 31 : rtmp = (cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%density)**(1.0_dp/3.0_dp)
305 31 : IF (.NOT. expl_cell) THEN
306 31 : helium_env(k)%helium%cell_size = rtmp
307 : ELSE
308 : CALL section_vals_val_get(helium_section, "CELL_SIZE", &
309 0 : r_val=helium_env(k)%helium%cell_size)
310 : ! only more work if not all three values are consistent:
311 0 : 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 0 : IF (expl_dens .AND. expl_nats) THEN
314 : msg_str = "DENSITY, NATOMS and CELL_SIZE options "// &
315 0 : "contradict each other"
316 0 : CPWARN(msg_str)
317 : END IF
318 : !ok we have enough freedom to resolve the conflict:
319 0 : IF (.NOT. expl_dens) THEN
320 0 : helium_env(k)%helium%density = cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%cell_size**3.0_dp
321 0 : IF (.NOT. expl_nats) THEN
322 : msg_str = "CELL_SIZE defined but neither "// &
323 0 : "NATOMS nor DENSITY given, using default NATOMS."
324 0 : 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 0 : 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 0 : )**(1.0_dp/3.0_dp)
334 0 : 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 0 : "to maintain correct density."
338 0 : CPWARN(msg_str)
339 0 : helium_env(k)%helium%cell_size = rtmp
340 : END IF
341 : END IF
342 : END IF
343 : END IF
344 31 : 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 118 : SELECT CASE (helium_env(k)%helium%cell_shape)
350 :
351 : CASE (helium_cell_shape_octahedron)
352 21 : helium_env(k)%helium%cell_m(1, 1) = -0.5_dp*helium_env(k)%helium%cell_size
353 21 : helium_env(k)%helium%cell_m(2, 1) = 0.5_dp*helium_env(k)%helium%cell_size
354 21 : helium_env(k)%helium%cell_m(3, 1) = 0.5_dp*helium_env(k)%helium%cell_size
355 21 : helium_env(k)%helium%cell_m(1, 2) = 0.5_dp*helium_env(k)%helium%cell_size
356 21 : helium_env(k)%helium%cell_m(2, 2) = -0.5_dp*helium_env(k)%helium%cell_size
357 21 : helium_env(k)%helium%cell_m(3, 2) = 0.5_dp*helium_env(k)%helium%cell_size
358 21 : helium_env(k)%helium%cell_m(1, 3) = 0.5_dp*helium_env(k)%helium%cell_size
359 21 : helium_env(k)%helium%cell_m(2, 3) = 0.5_dp*helium_env(k)%helium%cell_size
360 21 : helium_env(k)%helium%cell_m(3, 3) = -0.5_dp*helium_env(k)%helium%cell_size
361 21 : helium_env(k)%helium%cell_m_inv(1, 1) = 0.0_dp
362 21 : helium_env(k)%helium%cell_m_inv(2, 1) = 1.0_dp/helium_env(k)%helium%cell_size
363 21 : helium_env(k)%helium%cell_m_inv(3, 1) = 1.0_dp/helium_env(k)%helium%cell_size
364 21 : helium_env(k)%helium%cell_m_inv(1, 2) = 1.0_dp/helium_env(k)%helium%cell_size
365 21 : helium_env(k)%helium%cell_m_inv(2, 2) = 0.0_dp
366 21 : helium_env(k)%helium%cell_m_inv(3, 2) = 1.0_dp/helium_env(k)%helium%cell_size
367 21 : helium_env(k)%helium%cell_m_inv(1, 3) = 1.0_dp/helium_env(k)%helium%cell_size
368 21 : helium_env(k)%helium%cell_m_inv(2, 3) = 1.0_dp/helium_env(k)%helium%cell_size
369 21 : helium_env(k)%helium%cell_m_inv(3, 3) = 0.0_dp
370 : CASE (helium_cell_shape_cube)
371 :
372 10 : helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size
373 10 : helium_env(k)%helium%cell_m(2, 1) = 0.0_dp
374 10 : helium_env(k)%helium%cell_m(3, 1) = 0.0_dp
375 10 : helium_env(k)%helium%cell_m(1, 2) = 0.0_dp
376 10 : helium_env(k)%helium%cell_m(2, 2) = helium_env(k)%helium%cell_size
377 10 : helium_env(k)%helium%cell_m(3, 2) = 0.0_dp
378 10 : helium_env(k)%helium%cell_m(1, 3) = 0.0_dp
379 10 : helium_env(k)%helium%cell_m(2, 3) = 0.0_dp
380 10 : helium_env(k)%helium%cell_m(3, 3) = helium_env(k)%helium%cell_size
381 10 : helium_env(k)%helium%cell_m_inv(1, 1) = 1.0_dp/helium_env(k)%helium%cell_size
382 10 : helium_env(k)%helium%cell_m_inv(2, 1) = 0.0_dp
383 10 : helium_env(k)%helium%cell_m_inv(3, 1) = 0.0_dp
384 10 : helium_env(k)%helium%cell_m_inv(1, 2) = 0.0_dp
385 10 : helium_env(k)%helium%cell_m_inv(2, 2) = 1.0_dp/helium_env(k)%helium%cell_size
386 10 : helium_env(k)%helium%cell_m_inv(3, 2) = 0.0_dp
387 10 : helium_env(k)%helium%cell_m_inv(1, 3) = 0.0_dp
388 10 : helium_env(k)%helium%cell_m_inv(2, 3) = 0.0_dp
389 10 : helium_env(k)%helium%cell_m_inv(3, 3) = 1.0_dp/helium_env(k)%helium%cell_size
390 : CASE DEFAULT
391 0 : helium_env(k)%helium%cell_m(:, :) = 0.0_dp
392 31 : helium_env(k)%helium%cell_m_inv(:, :) = 0.0_dp
393 :
394 : END SELECT
395 :
396 : END DO ! k
397 :
398 25 : IF (logger%para_env%is_source()) THEN
399 : CALL section_vals_val_get(helium_section, "POTENTIAL_FILE_NAME", &
400 13 : c_val=potential_file_name)
401 : CALL open_file(file_name=TRIM(potential_file_name), &
402 13 : file_action="READ", file_status="OLD", unit_number=input_unit)
403 13 : READ (input_unit, "(A)") msg_str
404 13 : READ (msg_str, *, IOSTAT=i) nlines, pdx, tau, &
405 26 : x_spline, dx, he_mass
406 13 : IF (i /= 0) THEN
407 : ! old style potential file, use default mass and potential
408 13 : he_mass = 4.00263037059764_dp !< 4He mass in [u]
409 13 : expl_pot = .FALSE.
410 13 : READ (msg_str, *, IOSTAT=i) nlines, pdx, tau, &
411 26 : x_spline, dx
412 13 : IF (i /= 0) THEN
413 0 : msg_str = "Format/Read Error from Solvent POTENTIAL_FILE"
414 0 : CPABORT(msg_str)
415 : END IF
416 : ELSE
417 0 : expl_pot = .TRUE.
418 : ! in file really hb2m in kelvin angstrom**2
419 0 : 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 0 : tau = ABS(tau)
424 : END IF
425 13 : tau = kelvin/tau
426 13 : x_spline = x_spline/angstrom
427 13 : dx = dx/angstrom
428 : END IF
429 :
430 25 : CALL helium_env(1)%comm%bcast(nlines, logger%para_env%source)
431 25 : CALL helium_env(1)%comm%bcast(pdx, logger%para_env%source)
432 25 : CALL helium_env(1)%comm%bcast(tau, logger%para_env%source)
433 25 : CALL helium_env(1)%comm%bcast(x_spline, logger%para_env%source)
434 25 : CALL helium_env(1)%comm%bcast(dx, logger%para_env%source)
435 25 : CALL helium_env(1)%comm%bcast(he_mass, logger%para_env%source)
436 25 : isize = (pdx + 1)*(pdx + 2) + 1
437 100 : ALLOCATE (pot_transfer(nlines, isize))
438 25 : IF (logger%para_env%is_source()) THEN
439 13 : IF (expl_pot) THEN
440 0 : DO i = 1, nlines
441 0 : READ (input_unit, *) pot_transfer(i, :)
442 : END DO
443 : ELSE
444 1313 : DO i = 1, nlines
445 1300 : READ (input_unit, *) pot_transfer(i, 1:isize - 1)
446 : ! potential implicit, calculate it here now:
447 1313 : pot_transfer(i, isize) = helium_vij(x_spline + (i - 1)*dx)*kelvin
448 : END DO
449 : END IF
450 13 : CALL close_file(unit_number=input_unit)
451 : END IF
452 25 : CALL helium_env(1)%comm%bcast(pot_transfer, logger%para_env%source)
453 :
454 25 : CALL spline_data_create(helium_env(1)%helium%vij)
455 25 : CALL init_splinexy(helium_env(1)%helium%vij, nlines)
456 25 : helium_env(1)%helium%vij%x1 = x_spline
457 :
458 25 : CALL spline_data_create(helium_env(1)%helium%u0)
459 25 : CALL init_splinexy(helium_env(1)%helium%u0, nlines)
460 25 : helium_env(1)%helium%u0%x1 = x_spline
461 :
462 25 : CALL spline_data_create(helium_env(1)%helium%e0)
463 25 : CALL init_splinexy(helium_env(1)%helium%e0, nlines)
464 25 : helium_env(1)%helium%e0%x1 = x_spline
465 :
466 25 : isize = pdx + 1
467 25 : ntab = ((isize + 1)*isize)/2 - 1 ! -1 because endpoint approx treated separately
468 100 : ALLOCATE (helium_env(1)%helium%uoffdiag(ntab, 2, nlines))
469 100 : ALLOCATE (helium_env(1)%helium%eoffdiag(ntab, 2, nlines))
470 100 : DO j = 1, isize
471 250 : DO i = j, isize
472 150 : IF (i + j == 2) CYCLE ! endpoint approx later separately
473 125 : k = ((i - 1)*i)/2 + j
474 12625 : helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2)
475 125 : CALL init_spline(helium_env(1)%helium%vij, dx=dx)
476 12625 : helium_env(1)%helium%uoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:)
477 12625 : helium_env(1)%helium%uoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:)
478 125 : k = k + ((isize + 1)*isize)/2
479 12625 : helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2)/kelvin
480 125 : CALL init_spline(helium_env(1)%helium%vij, dx=dx)
481 12625 : helium_env(1)%helium%eoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:)
482 12625 : helium_env(1)%helium%eoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:)
483 225 : ntab = ntab - 1
484 : END DO
485 : END DO
486 :
487 25 : ntab = SIZE(pot_transfer, 2)
488 2525 : helium_env(1)%helium%vij%y(:) = pot_transfer(:, ntab)/kelvin
489 25 : CALL init_spline(helium_env(1)%helium%vij, dx=dx)
490 :
491 2525 : helium_env(1)%helium%u0%y(:) = pot_transfer(:, 1)
492 25 : CALL init_spline(helium_env(1)%helium%u0, dx=dx)
493 25 : k = ((isize + 1)*isize)/2 + 1
494 2525 : helium_env(1)%helium%e0%y(:) = pot_transfer(:, k)/kelvin
495 25 : CALL init_spline(helium_env(1)%helium%e0, dx=dx)
496 :
497 31 : DO k = 2, mepos
498 6 : helium_env(k)%helium%vij => helium_env(1)%helium%vij
499 6 : helium_env(k)%helium%u0 => helium_env(1)%helium%u0
500 6 : helium_env(k)%helium%e0 => helium_env(1)%helium%e0
501 6 : helium_env(k)%helium%uoffdiag => helium_env(1)%helium%uoffdiag
502 31 : helium_env(k)%helium%eoffdiag => helium_env(1)%helium%eoffdiag
503 : END DO
504 :
505 56 : DO k = 1, mepos
506 :
507 31 : helium_env(k)%helium%pdx = pdx
508 31 : 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 31 : mHe = he_mass*a_mass
515 : ! physical temperature [K]
516 31 : T = kelvin/helium_env(k)%helium%tau/helium_env(k)%helium%beads
517 : ! prefactors for calculating superfluid fractions [Angstrom^-2]
518 31 : helium_env(k)%helium%wpref = (((1e-20/h_bar)*mHe)/h_bar)*boltzmann*T
519 31 : helium_env(k)%helium%apref = (((4e-20/h_bar)*mHe)/h_bar)*boltzmann*T
520 :
521 31 : helium_env(k)%helium%he_mass_au = he_mass*massunit
522 31 : helium_env(k)%helium%hb2m = 1.0_dp/helium_env(k)%helium%he_mass_au
523 31 : helium_env(k)%helium%pweight = 0.0_dp
524 :
525 : ! Default in case sampling_method is not helium_sampling_worm.
526 31 : 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 31 : i_val=helium_env(k)%helium%sampling_method)
531 :
532 53 : 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 22 : i_val=helium_env(k)%helium%maxcycle)
537 22 : i = helium_env(k)%helium%maxcycle
538 22 : CPASSERT(i >= 0)
539 22 : i = helium_env(k)%helium%atoms - helium_env(k)%helium%maxcycle
540 22 : CPASSERT(i >= 0)
541 :
542 : ! set m-distribution parameters
543 : CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%DISTRIBUTION-TYPE", &
544 22 : i_val=i)
545 22 : CPASSERT(i >= 1)
546 22 : CPASSERT(i <= 6)
547 22 : helium_env(k)%helium%m_dist_type = i
548 : CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-VALUE", &
549 22 : i_val=i)
550 22 : CPASSERT(i >= 1)
551 22 : CPASSERT(i <= helium_env(k)%helium%maxcycle)
552 22 : helium_env(k)%helium%m_value = i
553 : CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-RATIO", &
554 22 : r_val=rtmp)
555 22 : CPASSERT(rtmp > 0.0_dp)
556 22 : CPASSERT(rtmp <= 1.0_dp)
557 22 : helium_env(k)%helium%m_ratio = rtmp
558 :
559 : CALL section_vals_val_get(helium_section, "CEPERLEY%BISECTION", &
560 22 : i_val=helium_env(k)%helium%bisection)
561 : ! precheck bisection value (not all invalids are filtered out here yet)
562 22 : i = helium_env(k)%helium%bisection
563 22 : CPASSERT(i > 1)
564 22 : i = helium_env(k)%helium%beads - helium_env(k)%helium%bisection
565 22 : CPASSERT(i > 0)
566 : !
567 22 : itmp = helium_env(k)%helium%bisection
568 22 : rtmp = 2.0_dp**(ANINT(LOG(REAL(itmp, dp))/LOG(2.0_dp)))
569 22 : tcheck = ABS(REAL(itmp, KIND=dp) - rtmp)
570 22 : IF (tcheck > 100.0_dp*EPSILON(0.0_dp)) THEN
571 0 : msg_str = "BISECTION should be integer power of 2."
572 0 : CPABORT(msg_str)
573 : END IF
574 22 : helium_env(k)%helium%bisctlog2 = NINT(LOG(REAL(itmp, dp))/LOG(2.0_dp))
575 :
576 : CASE (helium_sampling_worm)
577 9 : NULLIFY (input_worm)
578 : input_worm => section_vals_get_subs_vals(helium_env(k)%helium%input, &
579 9 : "MOTION%PINT%HELIUM%WORM")
580 : CALL section_vals_val_get(helium_section, "WORM%CENTROID_DRMAX", &
581 9 : r_val=helium_env(k)%helium%worm_centroid_drmax)
582 :
583 : CALL section_vals_val_get(helium_section, "WORM%STAGING_L", &
584 9 : i_val=helium_env(k)%helium%worm_staging_l)
585 :
586 : CALL section_vals_val_get(helium_section, "WORM%OPEN_CLOSE_SCALE", &
587 9 : r_val=helium_env(k)%helium%worm_open_close_scale)
588 :
589 : CALL section_vals_val_get(helium_section, "WORM%ALLOW_OPEN", &
590 9 : l_val=helium_env(k)%helium%worm_allow_open)
591 9 : IF (helium_env(k)%helium%atoms == 1) THEN
592 0 : IF (helium_env(k)%helium%worm_allow_open) THEN
593 : msg_str = "Default enabled open state sampling "// &
594 0 : "for only 1 He might be inefficient."
595 0 : CPWARN(msg_str)
596 : END IF
597 : END IF
598 :
599 : CALL section_vals_val_get(helium_section, "WORM%MAX_OPEN_CYCLES", &
600 9 : i_val=helium_env(k)%helium%worm_max_open_cycles)
601 :
602 9 : IF (helium_env(k)%helium%worm_staging_l + 1 >= helium_env(k)%helium%beads) THEN
603 0 : msg_str = "STAGING_L for worm sampling is too large"
604 0 : CPABORT(msg_str)
605 9 : ELSE IF (helium_env(k)%helium%worm_staging_l < 1) THEN
606 0 : msg_str = "STAGING_L must be positive integer"
607 0 : CPABORT(msg_str)
608 : END IF
609 :
610 : CALL section_vals_val_get(helium_section, "WORM%SHOW_STATISTICS", &
611 9 : 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 9 : rtmp = 2.0_dp*pi*helium_env(k)%helium%hb2m
616 9 : rtmp = SQRT(rtmp)
617 9 : rtmp = rtmp**3
618 9 : rtmp = rtmp*helium_env(k)%helium%worm_open_close_scale
619 9 : IF (helium_env(k)%helium%periodic) THEN
620 4 : rtmp = rtmp*helium_env(k)%helium%density
621 : ELSE
622 : rtmp = rtmp*helium_env(k)%helium%atoms/ &
623 5 : (4.0_dp/3.0_dp*pi*helium_env(k)%helium%droplet_radius**3)
624 : END IF
625 9 : helium_env(k)%helium%worm_ln_openclose_scale = LOG(rtmp)
626 :
627 : ! deal with acceptance statistics without changing the ceperley stuff
628 9 : helium_env(k)%helium%maxcycle = 1
629 9 : helium_env(k)%helium%bisctlog2 = 0
630 :
631 : ! get the absolute weights of the individual moves
632 9 : helium_env(k)%helium%worm_all_limit = 0
633 : CALL section_vals_val_get(helium_section, "WORM%CENTROID_WEIGHT", &
634 9 : i_val=itmp)
635 9 : helium_env(k)%helium%worm_centroid_min = 1
636 9 : helium_env(k)%helium%worm_centroid_max = itmp
637 9 : 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 9 : i_val=itmp)
640 9 : helium_env(k)%helium%worm_staging_min = helium_env(k)%helium%worm_centroid_max + 1
641 9 : helium_env(k)%helium%worm_staging_max = helium_env(k)%helium%worm_centroid_max + itmp
642 9 : helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
643 9 : IF (helium_env(k)%helium%worm_allow_open) THEN
644 : CALL section_vals_val_get(helium_section, "WORM%CRAWL_WEIGHT", &
645 9 : i_val=itmp)
646 9 : helium_env(k)%helium%worm_fcrawl_min = helium_env(k)%helium%worm_staging_max + 1
647 9 : helium_env(k)%helium%worm_fcrawl_max = helium_env(k)%helium%worm_staging_max + itmp
648 9 : helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
649 9 : helium_env(k)%helium%worm_bcrawl_min = helium_env(k)%helium%worm_fcrawl_max + 1
650 9 : helium_env(k)%helium%worm_bcrawl_max = helium_env(k)%helium%worm_fcrawl_max + itmp
651 9 : 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 9 : i_val=itmp)
654 9 : helium_env(k)%helium%worm_head_min = helium_env(k)%helium%worm_bcrawl_max + 1
655 9 : helium_env(k)%helium%worm_head_max = helium_env(k)%helium%worm_bcrawl_max + itmp
656 9 : helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp
657 9 : helium_env(k)%helium%worm_tail_min = helium_env(k)%helium%worm_head_max + 1
658 9 : helium_env(k)%helium%worm_tail_max = helium_env(k)%helium%worm_head_max + itmp
659 9 : 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 9 : i_val=itmp)
662 9 : helium_env(k)%helium%worm_swap_min = helium_env(k)%helium%worm_tail_max + 1
663 9 : helium_env(k)%helium%worm_swap_max = helium_env(k)%helium%worm_tail_max + itmp
664 9 : 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 9 : i_val=itmp)
667 9 : helium_env(k)%helium%worm_open_close_min = helium_env(k)%helium%worm_swap_max + 1
668 9 : helium_env(k)%helium%worm_open_close_max = helium_env(k)%helium%worm_swap_max + itmp
669 9 : 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 9 : 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 0 : msg_str = "Unknown helium sampling method!"
678 106 : CPABORT(msg_str)
679 : END SELECT
680 :
681 : ! ALLOCATE helium-related arrays
682 31 : i = helium_env(k)%helium%atoms
683 31 : j = helium_env(k)%helium%beads
684 124 : ALLOCATE (helium_env(k)%helium%pos(3, i, j))
685 57577 : helium_env(k)%helium%pos = 0.0_dp
686 124 : ALLOCATE (helium_env(k)%helium%work(3, i, j))
687 93 : ALLOCATE (helium_env(k)%helium%ptable(helium_env(k)%helium%maxcycle + 1))
688 93 : ALLOCATE (helium_env(k)%helium%permutation(i))
689 93 : ALLOCATE (helium_env(k)%helium%iperm(i))
690 124 : ALLOCATE (helium_env(k)%helium%tmatrix(i, i))
691 124 : ALLOCATE (helium_env(k)%helium%nmatrix(i, 2*i))
692 124 : ALLOCATE (helium_env(k)%helium%pmatrix(i, i))
693 124 : ALLOCATE (helium_env(k)%helium%ipmatrix(i, i))
694 31 : itmp = helium_env(k)%helium%bisctlog2 + 2
695 124 : ALLOCATE (helium_env(k)%helium%num_accepted(itmp, helium_env(k)%helium%maxcycle))
696 93 : ALLOCATE (helium_env(k)%helium%plength_avrg(helium_env(k)%helium%atoms))
697 93 : ALLOCATE (helium_env(k)%helium%plength_inst(helium_env(k)%helium%atoms))
698 93 : ALLOCATE (helium_env(k)%helium%atom_plength(helium_env(k)%helium%atoms))
699 31 : IF (helium_env(k)%helium%worm_max_open_cycles > 0) THEN
700 0 : ALLOCATE (helium_env(k)%helium%savepermutation(i))
701 0 : ALLOCATE (helium_env(k)%helium%saveiperm(i))
702 0 : ALLOCATE (helium_env(k)%helium%savepos(3, i, j))
703 : END IF
704 :
705 : ! check whether rdfs should be calculated and printed
706 31 : helium_env(k)%helium%rdf_present = helium_property_active(helium_env(k)%helium, "RDF")
707 31 : IF (helium_env(k)%helium%rdf_present) THEN
708 : ! allocate & initialize rdf related data structures
709 2 : CALL helium_rdf_init(helium_env(k)%helium)
710 : END IF
711 :
712 : ! check whether densities should be calculated and printed
713 31 : helium_env(k)%helium%rho_present = helium_property_active(helium_env(k)%helium, "RHO")
714 56 : IF (helium_env(k)%helium%rho_present) THEN
715 : ! allocate & initialize density related data structures
716 2 : NULLIFY (helium_env(k)%helium%rho_property)
717 2 : CALL helium_rho_init(helium_env(k)%helium)
718 : END IF
719 :
720 : END DO
721 :
722 : ! restore averages calculated in previous runs
723 25 : CALL helium_averages_restore(helium_env)
724 :
725 56 : DO k = 1, mepos
726 : ! fill in the solute-related data structures
727 31 : helium_env(k)%helium%e_corr = 0.0_dp
728 31 : IF (helium_env(k)%helium%solute_present) THEN
729 21 : 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 2 : helium_env(k)%helium%beads
733 : ! check if bead numbers are commensurate:
734 2 : i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%beads - helium_env(k)%helium%solute_beads
735 2 : IF (i /= 0) THEN
736 0 : msg_str = "Adjust number of solute beads to multiple of solvent beads."
737 0 : 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 2 : "Avoid too large coupling factors."
743 2 : CPWARN(msg_str)
744 19 : 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 18 : helium_env(k)%helium%solute_beads
748 : ! check if bead numbers are commensurate:
749 18 : i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%solute_beads - helium_env(k)%helium%beads
750 18 : IF (i /= 0) THEN
751 0 : msg_str = "Adjust number of solvent beads to multiple of solute beads."
752 0 : 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 18 : "sol_1 - he_1,2,3; sol_2 - he_4,5,6..."
757 18 : 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 21 : tcheck = ABS((helium_env(k)%helium%tau*helium_env(k)%helium%beads - solute%beta)/solute%beta)
763 21 : IF (tcheck > 1.0e-14_dp) THEN
764 0 : msg_str = "Tau, temperature and bead number are inconsistent."
765 0 : CPABORT(msg_str)
766 : END IF
767 :
768 21 : CALL helium_set_solute_indices(helium_env(k)%helium, solute)
769 21 : 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 21 : i_val=helium_env(k)%helium%solute_interaction)
774 21 : IF (helium_env(k)%helium%solute_interaction .EQ. helium_solute_intpot_nnp) THEN
775 1 : IF (k == 1) THEN
776 1 : NULLIFY (nnp_section)
777 1 : nnp_section => section_vals_get_subs_vals(helium_section, "NNP")
778 1 : CALL section_vals_get(nnp_section, explicit=explicit)
779 1 : msg_str = "NNP section not explicitly stated. Using default file names."
780 1 : IF (.NOT. explicit) CPWARN(msg_str)
781 : END IF
782 1 : ALLOCATE (helium_env(k)%helium%nnp)
783 1 : CALL cp_logger_create(tmplogger, para_env=helium_env(k)%comm, template_logger=logger)
784 1 : CALL cp_add_default_logger(tmplogger)
785 1 : CALL helium_init_nnp(helium_env(k)%helium, helium_env(k)%helium%nnp, nnp_section)
786 1 : CALL cp_rm_default_logger()
787 1 : CALL cp_logger_release(tmplogger)
788 : END IF
789 21 : 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 0 : "(see SOLUTE_INTERACTION keyword)"
793 0 : 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 84 : helium_env(k)%helium%solute_atoms*3))
799 : ALLOCATE (helium_env(k)%helium%force_inst(helium_env(k)%helium%solute_beads, &
800 84 : helium_env(k)%helium%solute_atoms*3))
801 :
802 63 : ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_1d(solute%p*solute%ndim))
803 63 : ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_np_1d(solute%p*solute%ndim*helium_env(k)%helium%num_env))
804 84 : ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_2d(solute%p, solute%ndim))
805 :
806 : ELSE
807 10 : helium_env(k)%helium%bead_ratio = 0
808 10 : IF (helium_env(k)%helium%periodic) THEN
809 : ! this assumes a specific potential (and its ugly):
810 10 : 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 10 : 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 93 : ALLOCATE (helium_env(k)%helium%itmp_atoms_1d(helium_env(k)%helium%atoms))
823 93 : ALLOCATE (helium_env(k)%helium%ltmp_atoms_1d(helium_env(k)%helium%atoms))
824 93 : ALLOCATE (helium_env(k)%helium%itmp_atoms_np_1d(helium_env(k)%helium%atoms*helium_env(k)%helium%num_env))
825 93 : 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 93 : helium_env(k)%helium%beads))
828 : ALLOCATE (helium_env(k)%helium%rtmp_3_atoms_beads_np_1d(3*helium_env(k)%helium%atoms* &
829 93 : 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 124 : helium_env(k)%helium%beads))
832 56 : IF (k .EQ. 1) THEN
833 25 : CALL helium_write_setup(helium_env(k)%helium)
834 : END IF
835 : END DO
836 25 : DEALLOCATE (pot_transfer)
837 : ELSE
838 : ! Deallocate env_all on processors without helium_env
839 1 : DEALLOCATE (env_all)
840 : END IF ! mepos .GT. 0
841 :
842 26 : NULLIFY (env_all)
843 26 : CALL timestop(handle)
844 :
845 26 : RETURN
846 104 : 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 26 : SUBROUTINE helium_release(helium_env)
856 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
857 :
858 : INTEGER :: k
859 :
860 26 : IF (ASSOCIATED(helium_env)) THEN
861 56 : DO k = 1, SIZE(helium_env)
862 31 : IF (k .EQ. 1) THEN
863 25 : CALL helium_env(k)%comm%free()
864 25 : DEALLOCATE (helium_env(k)%comm)
865 25 : DEALLOCATE (helium_env(k)%env_all)
866 : END IF
867 31 : NULLIFY (helium_env(k)%env_all)
868 :
869 : ! DEALLOCATE temporary arrays
870 0 : DEALLOCATE ( &
871 : helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
872 0 : helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
873 0 : helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
874 0 : helium_env(k)%helium%rtmp_3_np_1d, &
875 0 : helium_env(k)%helium%itmp_atoms_np_1d, &
876 0 : helium_env(k)%helium%ltmp_atoms_1d, &
877 31 : helium_env(k)%helium%itmp_atoms_1d)
878 :
879 : NULLIFY ( &
880 : helium_env(k)%helium%ltmp_3_atoms_beads_3d, &
881 31 : helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, &
882 31 : helium_env(k)%helium%rtmp_3_atoms_beads_1d, &
883 31 : helium_env(k)%helium%rtmp_3_np_1d, &
884 31 : helium_env(k)%helium%itmp_atoms_np_1d, &
885 31 : helium_env(k)%helium%ltmp_atoms_1d, &
886 31 : helium_env(k)%helium%itmp_atoms_1d &
887 31 : )
888 :
889 31 : IF (helium_env(k)%helium%solute_present) THEN
890 : ! DEALLOCATE solute-related arrays
891 0 : DEALLOCATE ( &
892 : helium_env(k)%helium%rtmp_p_ndim_2d, &
893 0 : helium_env(k)%helium%rtmp_p_ndim_np_1d, &
894 0 : helium_env(k)%helium%rtmp_p_ndim_1d, &
895 0 : helium_env(k)%helium%force_inst, &
896 21 : helium_env(k)%helium%force_avrg)
897 : NULLIFY ( &
898 : helium_env(k)%helium%rtmp_p_ndim_2d, &
899 21 : helium_env(k)%helium%rtmp_p_ndim_np_1d, &
900 21 : helium_env(k)%helium%rtmp_p_ndim_1d, &
901 21 : helium_env(k)%helium%force_inst, &
902 21 : helium_env(k)%helium%force_avrg)
903 : END IF
904 :
905 31 : IF (helium_env(k)%helium%rho_present) THEN
906 0 : DEALLOCATE ( &
907 : helium_env(k)%helium%rho_rstr, &
908 0 : helium_env(k)%helium%rho_accu, &
909 0 : helium_env(k)%helium%rho_inst, &
910 2 : helium_env(k)%helium%rho_incr)
911 : NULLIFY ( &
912 : helium_env(k)%helium%rho_rstr, &
913 2 : helium_env(k)%helium%rho_accu, &
914 2 : helium_env(k)%helium%rho_inst, &
915 2 : helium_env(k)%helium%rho_incr)
916 : ! DEALLOCATE everything
917 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix)
918 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_name)
919 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_index)
920 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix)
921 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_name)
922 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_index)
923 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix)
924 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_name)
925 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_index)
926 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix)
927 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_name)
928 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_index)
929 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix)
930 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name)
931 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index)
932 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix)
933 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name)
934 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index)
935 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix)
936 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_name)
937 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_index)
938 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix)
939 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_name)
940 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_index)
941 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix)
942 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name)
943 2 : DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index)
944 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix)
945 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name)
946 2 : NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index)
947 2 : DEALLOCATE (helium_env(k)%helium%rho_property)
948 2 : NULLIFY (helium_env(k)%helium%rho_property)
949 : END IF
950 :
951 31 : CALL helium_rdf_release(helium_env(k)%helium)
952 :
953 : ! DEALLOCATE helium-related arrays
954 0 : DEALLOCATE ( &
955 : helium_env(k)%helium%atom_plength, &
956 0 : helium_env(k)%helium%plength_inst, &
957 0 : helium_env(k)%helium%plength_avrg, &
958 0 : helium_env(k)%helium%num_accepted, &
959 0 : helium_env(k)%helium%ipmatrix, &
960 0 : helium_env(k)%helium%pmatrix, &
961 0 : helium_env(k)%helium%nmatrix, &
962 0 : helium_env(k)%helium%tmatrix, &
963 0 : helium_env(k)%helium%iperm, &
964 0 : helium_env(k)%helium%permutation, &
965 0 : helium_env(k)%helium%ptable, &
966 0 : helium_env(k)%helium%work, &
967 31 : helium_env(k)%helium%pos)
968 31 : IF (helium_env(k)%helium%worm_max_open_cycles > 0) THEN
969 0 : DEALLOCATE (helium_env(k)%helium%saveiperm, &
970 0 : helium_env(k)%helium%savepermutation, &
971 0 : helium_env(k)%helium%savepos)
972 : END IF
973 : NULLIFY ( &
974 : helium_env(k)%helium%atom_plength, &
975 31 : helium_env(k)%helium%plength_inst, &
976 31 : helium_env(k)%helium%plength_avrg, &
977 31 : helium_env(k)%helium%num_accepted, &
978 31 : helium_env(k)%helium%ipmatrix, &
979 31 : helium_env(k)%helium%pmatrix, &
980 31 : helium_env(k)%helium%nmatrix, &
981 31 : helium_env(k)%helium%tmatrix, &
982 31 : helium_env(k)%helium%iperm, &
983 31 : helium_env(k)%helium%saveiperm, &
984 31 : helium_env(k)%helium%permutation, &
985 31 : helium_env(k)%helium%savepermutation, &
986 31 : helium_env(k)%helium%ptable, &
987 31 : helium_env(k)%helium%work, &
988 31 : helium_env(k)%helium%pos, &
989 31 : helium_env(k)%helium%savepos &
990 31 : )
991 :
992 31 : IF (k == 1) THEN
993 25 : CALL spline_data_release(helium_env(k)%helium%vij)
994 25 : CALL spline_data_release(helium_env(k)%helium%u0)
995 25 : CALL spline_data_release(helium_env(k)%helium%e0)
996 25 : DEALLOCATE (helium_env(k)%helium%uoffdiag)
997 25 : DEALLOCATE (helium_env(k)%helium%eoffdiag)
998 : END IF
999 : NULLIFY (helium_env(k)%helium%uoffdiag, &
1000 31 : helium_env(k)%helium%eoffdiag, &
1001 31 : helium_env(k)%helium%vij, &
1002 31 : helium_env(k)%helium%u0, &
1003 31 : helium_env(k)%helium%e0)
1004 :
1005 31 : DEALLOCATE (helium_env(k)%helium%rng_stream_uniform)
1006 31 : DEALLOCATE (helium_env(k)%helium%rng_stream_gaussian)
1007 :
1008 : ! deallocate solute-related arrays
1009 31 : IF (helium_env(k)%helium%solute_present) THEN
1010 21 : DEALLOCATE (helium_env(k)%helium%solute_element)
1011 21 : NULLIFY (helium_env(k)%helium%solute_element)
1012 : END IF
1013 :
1014 : ! Deallocate everything from the helium_set_solute_indices
1015 31 : IF (ASSOCIATED(helium_env(k)%helium%ename)) THEN
1016 0 : DEALLOCATE (helium_env(k)%helium%ename)
1017 0 : NULLIFY (helium_env(k)%helium%ename)
1018 : END IF
1019 :
1020 : ! NNP interaction
1021 31 : IF (ASSOCIATED(helium_env(k)%helium%nnp)) THEN
1022 1 : CALL nnp_env_release(helium_env(k)%helium%nnp)
1023 1 : DEALLOCATE (helium_env(k)%helium%nnp)
1024 1 : NULLIFY (helium_env(k)%helium%nnp)
1025 : END IF
1026 31 : IF (ASSOCIATED(helium_env(k)%helium%nnp_sr_cut)) THEN
1027 1 : DEALLOCATE (helium_env(k)%helium%nnp_sr_cut)
1028 1 : NULLIFY (helium_env(k)%helium%nnp_sr_cut)
1029 : END IF
1030 :
1031 56 : DEALLOCATE (helium_env(k)%helium)
1032 :
1033 : END DO
1034 :
1035 25 : DEALLOCATE (helium_env)
1036 : NULLIFY (helium_env)
1037 : END IF
1038 26 : 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 26 : 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 26 : CALL timeset(routineN, handle)
1069 :
1070 26 : NULLIFY (logger)
1071 26 : logger => cp_get_default_logger()
1072 :
1073 26 : IF (ASSOCIATED(helium_env)) THEN
1074 :
1075 25 : NULLIFY (helium_section)
1076 : helium_section => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1077 25 : "MOTION%PINT%HELIUM")
1078 :
1079 : ! restore RNG state
1080 25 : NULLIFY (sec)
1081 25 : sec => section_vals_get_subs_vals(helium_section, "RNG_STATE")
1082 25 : CALL section_vals_get(sec, explicit=explicit)
1083 25 : IF (explicit) THEN
1084 6 : CALL helium_rng_restore(helium_env)
1085 : ELSE
1086 19 : CALL helium_write_line("RNG state initialized as new.")
1087 : END IF
1088 :
1089 : ! init/restore permutation state
1090 25 : NULLIFY (sec)
1091 25 : sec => section_vals_get_subs_vals(helium_section, "PERM")
1092 25 : CALL section_vals_get(sec, explicit=explicit)
1093 25 : IF (explicit) THEN
1094 6 : CALL helium_perm_restore(helium_env)
1095 : ELSE
1096 19 : CALL helium_perm_init(helium_env)
1097 19 : 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 56 : DO k = 1, SIZE(helium_env)
1102 : CALL section_vals_val_get(helium_section, "GET_FORCES", &
1103 56 : i_val=helium_env(k)%helium%get_helium_forces)
1104 : END DO
1105 :
1106 56 : DO k = 1, SIZE(helium_env)
1107 : ! init center of mass
1108 56 : IF (helium_env(k)%helium%solute_present) THEN
1109 84 : helium_env(k)%helium%center(:) = pint_com_pos(pint_env)
1110 : ELSE
1111 40 : 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 25 : NULLIFY (sec)
1117 25 : sec => section_vals_get_subs_vals(helium_section, "COORD")
1118 25 : CALL section_vals_get(sec, explicit=explicit)
1119 25 : IF (explicit) THEN
1120 6 : CALL helium_coord_restore(helium_env)
1121 6 : CALL helium_write_line("Coordinates restarted.")
1122 : ELSE
1123 19 : CALL section_vals_val_get(helium_section, "COORD_INIT_TEMP", r_val=initkT)
1124 19 : CALL section_vals_val_get(helium_section, "SOLUTE_RADIUS", r_val=solute_radius)
1125 19 : CALL helium_coord_init(helium_env, initkT, solute_radius)
1126 19 : IF (initkT > 0.0_dp) THEN
1127 12 : CALL helium_write_line("Coordinates initialized with thermal gaussian.")
1128 : ELSE
1129 7 : CALL helium_write_line("Coordinates initialized as point particles.")
1130 : END IF
1131 : END IF
1132 :
1133 56 : DO k = 1, SIZE(helium_env)
1134 :
1135 31 : helium_env(k)%helium%worm_is_closed = .TRUE.
1136 31 : helium_env(k)%helium%worm_atom_idx = 0
1137 31 : helium_env(k)%helium%worm_bead_idx = 0
1138 :
1139 57577 : helium_env(k)%helium%work(:, :, :) = helium_env(k)%helium%pos(:, :, :)
1140 :
1141 : ! init center of mass
1142 56 : IF (helium_env(k)%helium%solute_present) THEN
1143 84 : helium_env(k)%helium%center(:) = pint_com_pos(pint_env)
1144 : ELSE
1145 10 : IF (helium_env(k)%helium%periodic) THEN
1146 40 : helium_env(k)%helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
1147 : ELSE
1148 0 : 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 25 : l_val=presample)
1157 25 : coords_presampled = .FALSE.
1158 25 : IF (presample) THEN
1159 30 : DO k = 1, SIZE(helium_env)
1160 30 : helium_env(k)%helium%current_step = 0
1161 : END DO
1162 12 : CALL helium_sample(helium_env, pint_env)
1163 30 : DO k = 1, SIZE(helium_env)
1164 558 : IF (helium_env(k)%helium%solute_present) helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
1165 198 : helium_env(k)%helium%energy_avrg(:) = 0.0_dp
1166 324 : helium_env(k)%helium%plength_avrg(:) = 0.0_dp
1167 228 : helium_env(k)%helium%num_accepted(:, :) = 0.0_dp
1168 : ! Reset properties accumulated over presample:
1169 72 : helium_env(k)%helium%proarea%accu(:) = 0.0_dp
1170 72 : helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
1171 72 : helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
1172 72 : helium_env(k)%helium%mominer%accu(:) = 0.0_dp
1173 18 : IF (helium_env(k)%helium%rho_present) THEN
1174 4222 : helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
1175 : END IF
1176 30 : IF (helium_env(k)%helium%rdf_present) THEN
1177 1002 : helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp
1178 : END IF
1179 : END DO
1180 12 : coords_presampled = .TRUE.
1181 12 : CALL helium_write_line("Bead coordinates pre-sampled.")
1182 : END IF
1183 :
1184 25 : IF (helium_env(1)%helium%solute_present) THEN
1185 : ! restore helium forces
1186 15 : NULLIFY (sec)
1187 15 : sec => section_vals_get_subs_vals(helium_section, "FORCE")
1188 15 : CALL section_vals_get(sec, explicit=explicit)
1189 15 : IF (explicit) THEN
1190 2 : IF (.NOT. coords_presampled) THEN
1191 2 : CALL helium_force_restore(helium_env)
1192 : END IF
1193 : ELSE
1194 13 : IF (.NOT. coords_presampled) THEN
1195 7 : CALL helium_force_init(helium_env)
1196 7 : 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 26 : CALL timestop(handle)
1210 :
1211 26 : 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 19 : 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 19 : 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 19 : minHeHedsttmp = minHeHedst
1247 :
1248 44 : DO k = 1, SIZE(helium_env)
1249 25 : IF (helium_env(k)%helium%solute_present) THEN
1250 76 : cvek(:) = helium_env(k)%helium%center(:)
1251 19 : CALL helium_pbc(helium_env(k)%helium, cvek)
1252 : END IF
1253 :
1254 75 : ALLOCATE (centroids(3, helium_env(k)%helium%atoms))
1255 25 : IF (helium_env(k)%helium%periodic) THEN
1256 330 : DO ia = 1, helium_env(k)%helium%atoms
1257 : invalidpos = .TRUE.
1258 : iter = 0
1259 343843 : DO WHILE (invalidpos)
1260 343513 : iter = iter + 1
1261 343513 : invalidpos = .FALSE.
1262 : ! if sampling fails to often, reduce he he criterion
1263 : !CS TODO:
1264 : !minHeHedsttmp = 0.90_dp**(iter/100)*minHeHedst
1265 343513 : minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst
1266 1374052 : DO ic = 1, 3
1267 1030539 : r1 = helium_env(k)%helium%rng_stream_uniform%next()
1268 1030539 : r1 = 2.0_dp*r1 - 1.0_dp
1269 1030539 : r1 = r1*helium_env(k)%helium%cell_size
1270 1374052 : centroids(ic, ia) = r1
1271 : END DO
1272 : ! check if helium is outside of cell
1273 1374052 : tvek(:) = centroids(:, ia)
1274 343513 : CALL helium_pbc(helium_env(k)%helium, tvek(:))
1275 1374052 : rvek(:) = tvek(:) - centroids(:, ia)
1276 1374052 : r2 = DOT_PRODUCT(rvek, rvek)
1277 343513 : IF (r2 > 1.0_dp*10.0_dp**(-6)) THEN
1278 : invalidpos = .TRUE.
1279 : ELSE
1280 : ! check for helium-helium collision
1281 201908 : DO id = 1, ia - 1
1282 806348 : rvek = centroids(:, ia) - centroids(:, id)
1283 201587 : CALL helium_pbc(helium_env(k)%helium, rvek)
1284 806348 : r2 = DOT_PRODUCT(rvek, rvek)
1285 201908 : IF (r2 < minHeHedsttmp**2) THEN
1286 : invalidpos = .TRUE.
1287 : EXIT
1288 : END IF
1289 : END DO
1290 : END IF
1291 21775 : IF (.NOT. invalidpos) THEN
1292 : ! check if centroid collides with molecule
1293 321 : IF (helium_env(k)%helium%solute_present) THEN
1294 516 : rvek(:) = (cvek(:) - centroids(:, ia))
1295 516 : r2 = DOT_PRODUCT(rvek, rvek)
1296 129 : 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 10 : IF (initkT > 0.0_dp) THEN
1303 10 : CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkT)
1304 : ELSE
1305 0 : DO ia = 1, helium_env(k)%helium%atoms
1306 0 : DO ib = 1, helium_env(k)%helium%beads
1307 0 : 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 330 : DO ia = 1, helium_env(k)%helium%atoms
1313 7178 : DO ib = 1, helium_env(k)%helium%beads
1314 6848 : CALL helium_pbc(helium_env(k)%helium, helium_env(k)%helium%pos(:, ia, ib))
1315 : ! check if bead collides with molecule
1316 7168 : IF (helium_env(k)%helium%solute_present) THEN
1317 8192 : rvek(:) = (cvek(:) - helium_env(k)%helium%pos(:, ia, ib))
1318 8192 : r2 = DOT_PRODUCT(rvek, rvek)
1319 2048 : IF (r2 <= solute_radius**2) THEN
1320 2 : r1 = SQRT(r2)
1321 : helium_env(k)%helium%pos(:, ia, ib) = &
1322 8 : cvek(:) + solute_radius/r1*rvek(:)
1323 : END IF
1324 : END IF
1325 : END DO
1326 : END DO
1327 : ELSE
1328 213 : DO ia = 1, helium_env(k)%helium%atoms
1329 : iter = 0
1330 : invalidpos = .TRUE.
1331 601 : DO WHILE (invalidpos)
1332 388 : invalidpos = .FALSE.
1333 388 : iter = iter + 1
1334 : ! if sampling fails to often, reduce he he criterion
1335 388 : minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst
1336 1552 : DO ic = 1, 3
1337 1164 : rvek(ic) = helium_env(k)%helium%rng_stream_uniform%next()
1338 1164 : rvek(ic) = 2.0_dp*rvek(ic) - 1.0_dp
1339 1552 : rvek(ic) = rvek(ic)*helium_env(k)%helium%droplet_radius
1340 : END DO
1341 1552 : centroids(:, ia) = rvek(:)
1342 : ! check if helium is outside of the droplet
1343 1552 : r2 = DOT_PRODUCT(rvek, rvek)
1344 388 : IF (r2 > helium_env(k)%helium%droplet_radius**2) THEN
1345 : invalidpos = .TRUE.
1346 : ELSE
1347 : ! check for helium-helium collision
1348 2550 : DO id = 1, ia - 1
1349 9408 : rvek = centroids(:, ia) - centroids(:, id)
1350 9408 : r2 = DOT_PRODUCT(rvek, rvek)
1351 2550 : IF (r2 < minHeHedsttmp**2) THEN
1352 : invalidpos = .TRUE.
1353 : EXIT
1354 : END IF
1355 : END DO
1356 : END IF
1357 404 : IF (.NOT. invalidpos) THEN
1358 : ! make sure the helium does not collide with the solute
1359 198 : IF (helium_env(k)%helium%solute_present) THEN
1360 792 : rvek(:) = (cvek(:) - centroids(:, ia))
1361 792 : r2 = DOT_PRODUCT(rvek, rvek)
1362 198 : 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 15 : IF (initkT > 0.0_dp) THEN
1369 5 : CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkT)
1370 : ELSE
1371 183 : DO ia = 1, helium_env(k)%helium%atoms
1372 2951 : DO ib = 1, helium_env(k)%helium%beads
1373 11245 : helium_env(k)%helium%pos(:, ia, ib) = centroids(:, ia)
1374 : END DO
1375 : END DO
1376 : END IF
1377 213 : DO ia = 1, helium_env(k)%helium%atoms
1378 3381 : 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 12672 : helium_env(k)%helium%pos(:, ia, ib))
1382 3168 : IF (r1 > helium_env(k)%helium%droplet_radius**2) THEN
1383 16 : r1 = SQRT(r1)
1384 : helium_env(k)%helium%pos(:, ia, ib) = &
1385 : helium_env(k)%helium%droplet_radius/r1* &
1386 64 : helium_env(k)%helium%pos(:, ia, ib)
1387 3152 : ELSE IF (helium_env(k)%helium%solute_present) THEN
1388 3152 : IF (r1 < solute_radius**2) THEN
1389 : !make sure that nothing lies within the molecule
1390 16 : r1 = SQRT(r1)
1391 : helium_env(k)%helium%pos(:, ia, ib) = &
1392 : solute_radius/r1* &
1393 64 : 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 12870 : helium_env(k)%helium%center(:)
1400 : END DO
1401 : END DO
1402 : END IF
1403 40543 : helium_env(k)%helium%work = helium_env(k)%helium%pos
1404 44 : DEALLOCATE (centroids)
1405 : END DO
1406 :
1407 19 : RETURN
1408 19 : END SUBROUTINE helium_coord_init
1409 :
1410 : ! **************************************************************************************************
1411 : !> \brief ...
1412 : !> \param helium_env ...
1413 : !> \param centroids ...
1414 : !> \param kbT ...
1415 : ! **************************************************************************************************
1416 15 : 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 30 : DIMENSION(helium_env%beads, helium_env%beads) :: u2x
1428 30 : REAL(KIND=dp), DIMENSION(helium_env%beads) :: nmhecoords
1429 :
1430 15 : p = helium_env%beads
1431 :
1432 15 : sqrt2p = SQRT(2.0_dp/REAL(p, dp))
1433 15 : twopip = twopi/REAL(p, dp)
1434 15 : pip = pi/REAL(p, dp)
1435 15 : sqrtp = SQRT(REAL(p, dp))
1436 15 : invsqrtp = 1.0_dp/SQRT(REAL(p, dp))
1437 :
1438 : ! set up normal mode backtransform matrix
1439 6363 : u2x(:, :) = 0.0_dp
1440 309 : u2x(:, 1) = invsqrtp
1441 159 : DO i = 2, p/2 + 1
1442 3111 : DO j = 1, p
1443 3096 : u2x(j, i) = sqrt2p*COS(twopip*(i - 1)*(j - 1))
1444 : END DO
1445 : END DO
1446 150 : DO i = p/2 + 2, p
1447 2958 : DO j = 1, p
1448 2943 : u2x(j, i) = sqrt2p*SIN(twopip*(i - 1)*(j - 1))
1449 : END DO
1450 : END DO
1451 15 : IF (MOD(p, 2) == 0) THEN
1452 9 : DO i = 1, p - 1, 2
1453 72 : u2x(i, p/2 + 1) = invsqrtp
1454 72 : u2x(i + 1, p/2 + 1) = -1.0_dp*invsqrtp
1455 : END DO
1456 : END IF
1457 :
1458 360 : DO iatom = 1, helium_env%atoms
1459 1395 : DO idim = 1, 3
1460 1035 : nmhecoords(1) = sqrtp*centroids(idim, iatom)
1461 21744 : DO imode = 2, p
1462 20709 : omega = 2.0_dp*p*kbT*SIN((imode - 1)*pip)
1463 20709 : variance = kbT*p/(helium_env%he_mass_au*omega**2)
1464 20709 : rand = helium_env%rng_stream_gaussian%next()
1465 21744 : nmhecoords(imode) = rand*SQRT(variance)
1466 : END DO
1467 522372 : helium_env%pos(idim, iatom, 1:p) = MATMUL(u2x, nmhecoords)
1468 : END DO
1469 : END DO
1470 :
1471 15 : 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 6 : 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 6 : LOGICAL, DIMENSION(:, :, :), POINTER :: m
1491 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: message
1492 6 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: f
1493 : TYPE(cp_logger_type), POINTER :: logger
1494 :
1495 6 : NULLIFY (logger)
1496 12 : 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 6 : NULLIFY (message)
1501 : CALL section_vals_val_get(helium_env(1)%helium%input, &
1502 : "MOTION%PINT%HELIUM%COORD%_DEFAULT_KEYWORD_", &
1503 6 : r_vals=message)
1504 :
1505 : ! check that the number of values in the input match the current runtime
1506 6 : actlen = SIZE(message)
1507 6 : num_env_restart = actlen/helium_env(1)%helium%atoms/helium_env(1)%helium%beads/3
1508 :
1509 6 : IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
1510 0 : err_str = "Reading bead coordinates from the input file."
1511 0 : CALL helium_write_line(err_str)
1512 0 : err_str = "Number of environments in the restart...: '"
1513 0 : stmp = ""
1514 0 : WRITE (stmp, *) num_env_restart
1515 0 : err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
1516 0 : CALL helium_write_line(err_str)
1517 0 : err_str = "Number of current run time environments.: '"
1518 0 : stmp = ""
1519 0 : WRITE (stmp, *) helium_env(1)%helium%num_env
1520 0 : err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
1521 0 : CALL helium_write_line(err_str)
1522 0 : err_str = "Missmatch between number of bead coord. in input file and helium environments."
1523 0 : CPABORT(err_str)
1524 : ELSE
1525 6 : CALL helium_write_line("Bead coordinates read from the input file.")
1526 :
1527 6 : offset = 0
1528 9 : DO i = 1, logger%para_env%mepos
1529 9 : offset = offset + helium_env(1)%env_all(i)
1530 : END DO
1531 :
1532 : ! distribute coordinates over processors (no message passing)
1533 12 : DO k = 1, SIZE(helium_env)
1534 6 : msglen = helium_env(k)%helium%atoms*helium_env(k)%helium%beads*3
1535 6 : off = msglen*MOD(offset + k - 1, num_env_restart)
1536 : NULLIFY (m, f)
1537 24 : ALLOCATE (m(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads))
1538 24 : ALLOCATE (f(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads))
1539 17034 : m(:, :, :) = .TRUE.
1540 17034 : f(:, :, :) = 0.0_dp
1541 17034 : helium_env(k)%helium%pos(:, :, 1:helium_env(k)%helium%beads) = UNPACK(message(off + 1:off + msglen), MASK=m, FIELD=f)
1542 12 : DEALLOCATE (f, m)
1543 : END DO
1544 :
1545 : END IF
1546 :
1547 : NULLIFY (message)
1548 :
1549 6 : RETURN
1550 6 : 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 7 : SUBROUTINE helium_force_init(helium_env)
1561 :
1562 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1563 :
1564 : INTEGER :: k
1565 :
1566 14 : DO k = 1, SIZE(helium_env)
1567 14 : IF (helium_env(k)%helium%solute_present) THEN
1568 934 : helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
1569 934 : helium_env(k)%helium%force_inst(:, :) = 0.0_dp
1570 : END IF
1571 : END DO
1572 :
1573 7 : 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 2 : 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 2 : LOGICAL, DIMENSION(:, :), POINTER :: m
1590 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: message
1591 2 : 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 2 : NULLIFY (message)
1597 : CALL section_vals_val_get(helium_env(1)%helium%input, &
1598 : "MOTION%PINT%HELIUM%FORCE%_DEFAULT_KEYWORD_", &
1599 2 : r_vals=message)
1600 :
1601 : ! check if the destination array has correct size
1602 2 : msglen = helium_env(1)%helium%solute_atoms*helium_env(1)%helium%solute_beads*3
1603 2 : actlen = SIZE(helium_env(1)%helium%force_avrg)
1604 2 : err_str = "Invalid size of helium%force_avrg array: actual '"
1605 2 : stmp = ""
1606 2 : WRITE (stmp, *) actlen
1607 : err_str = TRIM(ADJUSTL(err_str))// &
1608 2 : TRIM(ADJUSTL(stmp))//"' but expected '"
1609 2 : stmp = ""
1610 2 : WRITE (stmp, *) msglen
1611 2 : IF (actlen /= msglen) THEN
1612 : err_str = TRIM(ADJUSTL(err_str))// &
1613 0 : TRIM(ADJUSTL(stmp))//"'."
1614 0 : CPABORT(err_str)
1615 : END IF
1616 :
1617 : ! restore forces on all processors (no message passing)
1618 : NULLIFY (m, f)
1619 8 : ALLOCATE (m(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3))
1620 8 : ALLOCATE (f(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3))
1621 92 : m(:, :) = .TRUE.
1622 92 : f(:, :) = 0.0_dp
1623 4 : DO k = 1, SIZE(helium_env)
1624 92 : helium_env(k)%helium%force_avrg(:, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f)
1625 94 : helium_env(k)%helium%force_inst(:, :) = 0.0_dp
1626 : END DO
1627 2 : DEALLOCATE (f, m)
1628 :
1629 2 : CALL helium_write_line("Forces on the solute read from the input file.")
1630 :
1631 : NULLIFY (message)
1632 :
1633 2 : RETURN
1634 2 : 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 19 : SUBROUTINE helium_perm_init(helium_env)
1647 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1648 :
1649 : INTEGER :: ia, k
1650 :
1651 44 : DO k = 1, SIZE(helium_env)
1652 562 : DO ia = 1, helium_env(k)%helium%atoms
1653 518 : helium_env(k)%helium%permutation(ia) = ia
1654 543 : helium_env(k)%helium%iperm(ia) = ia
1655 : END DO
1656 : END DO
1657 :
1658 19 : 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 6 : 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 6 : INTEGER, DIMENSION(:), POINTER :: message
1681 : TYPE(cp_logger_type), POINTER :: logger
1682 :
1683 6 : NULLIFY (logger)
1684 12 : 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 6 : NULLIFY (message)
1689 : CALL section_vals_val_get(helium_env(1)%helium%input, &
1690 : "MOTION%PINT%HELIUM%PERM%_DEFAULT_KEYWORD_", &
1691 6 : i_vals=message)
1692 :
1693 : ! check the number of environments presumably stored in the restart
1694 6 : actlen = SIZE(message)
1695 6 : 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 6 : IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
1701 0 : err_str = "Reading permutation state from the input file."
1702 0 : CALL helium_write_line(err_str)
1703 0 : err_str = "Number of environments in the restart...: '"
1704 0 : stmp = ""
1705 0 : WRITE (stmp, *) num_env_restart
1706 0 : err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
1707 0 : CALL helium_write_line(err_str)
1708 0 : err_str = "Number of current run time environments.: '"
1709 0 : stmp = ""
1710 0 : WRITE (stmp, *) helium_env(1)%helium%num_env
1711 0 : err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
1712 0 : CALL helium_write_line(err_str)
1713 0 : err_str = "Missmatch between number of perm. states in input file and helium environments."
1714 0 : CPABORT(err_str)
1715 : ELSE
1716 6 : CALL helium_write_line("Permutation state read from the input file.")
1717 :
1718 : ! distribute permutation state over processors
1719 6 : offset = 0
1720 9 : DO i = 1, logger%para_env%mepos
1721 9 : offset = offset + helium_env(1)%env_all(i)
1722 : END DO
1723 :
1724 12 : DO k = 1, SIZE(helium_env)
1725 6 : msglen = helium_env(k)%helium%atoms
1726 6 : off = msglen*MOD(k - 1 + offset, num_env_restart)
1727 402 : helium_env(k)%helium%permutation(:) = message(off + 1:off + msglen)
1728 : END DO
1729 : END IF
1730 :
1731 : ! recalculate the inverse permutation array
1732 12 : DO k = 1, SIZE(helium_env)
1733 198 : helium_env(k)%helium%iperm(:) = 0
1734 6 : ic = 0
1735 198 : DO ia = 1, msglen
1736 198 : IF ((helium_env(k)%helium%permutation(ia) > 0) .AND. (helium_env(k)%helium%permutation(ia) <= msglen)) THEN
1737 192 : helium_env(k)%helium%iperm(helium_env(k)%helium%permutation(ia)) = ia
1738 192 : ic = ic + 1
1739 : END IF
1740 : END DO
1741 6 : err_str = "Invalid HELIUM%PERM state: some numbers not within (1,"
1742 6 : stmp = ""
1743 6 : WRITE (stmp, *) msglen
1744 12 : IF (ic /= msglen) THEN
1745 : err_str = TRIM(ADJUSTL(err_str))// &
1746 0 : TRIM(ADJUSTL(stmp))//")."
1747 0 : CPABORT(err_str)
1748 : END IF
1749 : END DO
1750 : NULLIFY (message)
1751 :
1752 6 : RETURN
1753 6 : 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 125 : 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 25 : REAL(KIND=dp), DIMENSION(:), POINTER :: message
1771 : TYPE(cp_logger_type), POINTER :: logger
1772 :
1773 25 : NULLIFY (logger)
1774 50 : logger => cp_get_default_logger()
1775 :
1776 25 : offset = 0
1777 37 : DO i = 1, logger%para_env%mepos
1778 37 : 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 25 : explicit=explicit)
1785 25 : IF (explicit) THEN
1786 0 : NULLIFY (message)
1787 : CALL section_vals_val_get(helium_env(1)%helium%input, &
1788 : "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", &
1789 0 : r_vals=message)
1790 0 : num_env_restart = SIZE(message)/3 ! apparent number of environments
1791 0 : msglen = 3
1792 0 : DO k = 1, SIZE(helium_env)
1793 0 : off = msglen*MOD(offset + k - 1, num_env_restart)
1794 0 : helium_env(k)%helium%proarea%rstr(:) = message(off + 1:off + msglen)
1795 : END DO
1796 : ELSE
1797 56 : DO k = 1, SIZE(helium_env)
1798 149 : 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 25 : explicit=explicit)
1806 25 : IF (explicit) THEN
1807 0 : NULLIFY (message)
1808 : CALL section_vals_val_get(helium_env(1)%helium%input, &
1809 : "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", &
1810 0 : r_vals=message)
1811 0 : num_env_restart = SIZE(message)/3 ! apparent number of environments
1812 0 : msglen = 3
1813 0 : DO k = 1, SIZE(helium_env)
1814 0 : off = msglen*MOD(offset + k - 1, num_env_restart)
1815 0 : helium_env(k)%helium%prarea2%rstr(:) = message(off + 1:off + msglen)
1816 : END DO
1817 : ELSE
1818 56 : DO k = 1, SIZE(helium_env)
1819 149 : 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 25 : explicit=explicit)
1827 25 : IF (explicit) THEN
1828 0 : NULLIFY (message)
1829 : CALL section_vals_val_get(helium_env(1)%helium%input, &
1830 : "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", &
1831 0 : r_vals=message)
1832 0 : num_env_restart = SIZE(message)/3 ! apparent number of environments
1833 0 : msglen = 3
1834 0 : DO k = 1, SIZE(helium_env)
1835 0 : off = msglen*MOD(offset + k - 1, num_env_restart)
1836 0 : helium_env(k)%helium%wnmber2%rstr(:) = message(off + 1:off + msglen)
1837 : END DO
1838 : ELSE
1839 56 : DO k = 1, SIZE(helium_env)
1840 149 : 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 25 : explicit=explicit)
1848 25 : IF (explicit) THEN
1849 0 : NULLIFY (message)
1850 : CALL section_vals_val_get(helium_env(1)%helium%input, &
1851 : "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", &
1852 0 : r_vals=message)
1853 0 : num_env_restart = SIZE(message)/3 ! apparent number of environments
1854 0 : msglen = 3
1855 0 : DO k = 1, SIZE(helium_env)
1856 0 : off = msglen*MOD(offset + k - 1, num_env_restart)
1857 0 : helium_env(k)%helium%mominer%rstr(:) = message(off + 1:off + msglen)
1858 : END DO
1859 : ELSE
1860 56 : DO k = 1, SIZE(helium_env)
1861 149 : helium_env(k)%helium%mominer%rstr(:) = 0.0_dp
1862 : END DO
1863 : END IF
1864 :
1865 25 : IF (helium_env(1)%helium%rdf_present) THEN
1866 2 : CALL helium_rdf_restore(helium_env)
1867 : END IF
1868 :
1869 25 : IF (helium_env(1)%helium%rho_present) THEN
1870 : ! restore densities
1871 2 : CALL helium_rho_restore(helium_env)
1872 : END IF
1873 :
1874 : ! get the weighting factor
1875 56 : 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 31 : 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 56 : l_val=helium_env(k)%helium%averages_restarted)
1886 : END DO
1887 :
1888 25 : RETURN
1889 25 : 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 25 : 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 25 : TYPE(rng_stream_p_type), DIMENSION(:), POINTER :: gaussian_array, uniform_array
1909 :
1910 25 : NULLIFY (logger)
1911 50 : logger => cp_get_default_logger()
1912 25 : 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 13 : i_val=helium_seed)
1916 : END IF
1917 : CALL helium_env(1)%comm%bcast(helium_seed, &
1918 25 : logger%para_env%source)
1919 225 : initial_seed(:, :) = REAL(helium_seed, dp)
1920 :
1921 : ALLOCATE (uniform_array(helium_env(1)%helium%num_env), &
1922 247 : gaussian_array(helium_env(1)%helium%num_env))
1923 86 : DO i = 1, helium_env(1)%helium%num_env
1924 3014 : 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 25 : 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 25 : last_rng_stream=uniform_array(1)%stream)
1941 61 : 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 36 : 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 61 : last_rng_stream=uniform_array(i)%stream)
1951 : END DO
1952 :
1953 25 : offset = 0
1954 37 : DO i = 1, logger%para_env%mepos
1955 37 : offset = offset + helium_env(1)%env_all(i)
1956 : END DO
1957 :
1958 56 : DO i = 1, SIZE(helium_env)
1959 : NULLIFY (helium_env(i)%helium%rng_stream_uniform, &
1960 31 : helium_env(i)%helium%rng_stream_gaussian)
1961 31 : helium_env(i)%helium%rng_stream_uniform => uniform_array(offset + i)%stream
1962 56 : helium_env(i)%helium%rng_stream_gaussian => gaussian_array(offset + i)%stream
1963 : END DO
1964 :
1965 86 : DO i = 1, helium_env(1)%helium%num_env
1966 61 : 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 30 : DEALLOCATE (uniform_array(i)%stream)
1969 30 : DEALLOCATE (gaussian_array(i)%stream)
1970 : END IF
1971 61 : NULLIFY (uniform_array(i)%stream)
1972 86 : NULLIFY (gaussian_array(i)%stream)
1973 : END DO
1974 :
1975 25 : DEALLOCATE (uniform_array)
1976 25 : DEALLOCATE (gaussian_array)
1977 25 : 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 6 : 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 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: message
2000 : TYPE(cp_logger_type), POINTER :: logger
2001 :
2002 6 : NULLIFY (logger)
2003 12 : 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 6 : NULLIFY (message)
2008 : CALL section_vals_val_get(helium_env(1)%helium%input, &
2009 : "MOTION%PINT%HELIUM%RNG_STATE%_DEFAULT_KEYWORD_", &
2010 6 : r_vals=message)
2011 :
2012 : ! check the number of environments presumably stored in the restart
2013 6 : actlen = SIZE(message)
2014 6 : num_env_restart = actlen/40
2015 :
2016 : ! check, if RNG restart has the same dimension as helium%num_env
2017 6 : IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN
2018 0 : err_str = "Reading RNG state from the input file."
2019 0 : CALL helium_write_line(err_str)
2020 0 : err_str = "Number of environments in the restart...: '"
2021 0 : stmp = ""
2022 0 : WRITE (stmp, *) num_env_restart
2023 0 : err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
2024 0 : CALL helium_write_line(err_str)
2025 0 : err_str = "Number of current run time environments.: '"
2026 0 : stmp = ""
2027 0 : WRITE (stmp, *) helium_env(1)%helium%num_env
2028 0 : err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'."
2029 0 : CALL helium_write_line(err_str)
2030 0 : err_str = "Missmatch between number of RNG states in input file and helium environments."
2031 0 : CPABORT(err_str)
2032 : ELSE
2033 6 : CALL helium_write_line("RNG state read from the input file.")
2034 :
2035 : ! unpack the buffer at each processor, set RNG state
2036 6 : offset = 0
2037 9 : DO i = 1, logger%para_env%mepos
2038 9 : offset = offset + helium_env(1)%env_all(i)
2039 : END DO
2040 :
2041 12 : DO k = 1, SIZE(helium_env)
2042 6 : msglen = 40
2043 6 : off = msglen*(offset + k - 1)
2044 54 : m(:, :) = .TRUE.
2045 6 : f(:, :) = 0.0_dp
2046 6 : bg(:, :) = UNPACK(message(off + 1:off + 6), MASK=m, FIELD=f)
2047 6 : cg(:, :) = UNPACK(message(off + 7:off + 12), MASK=m, FIELD=f)
2048 6 : ig(:, :) = UNPACK(message(off + 13:off + 18), MASK=m, FIELD=f)
2049 6 : bf = message(off + 19)
2050 6 : bu = message(off + 20)
2051 6 : IF (bf .GT. 0) THEN
2052 0 : lbf = .TRUE.
2053 : ELSE
2054 6 : lbf = .FALSE.
2055 : END IF
2056 : CALL helium_env(k)%helium%rng_stream_uniform%set(bg=bg, cg=cg, ig=ig, &
2057 6 : buffer=bu, buffer_filled=lbf)
2058 6 : bg(:, :) = UNPACK(message(off + 21:off + 26), MASK=m, FIELD=f)
2059 6 : cg(:, :) = UNPACK(message(off + 27:off + 32), MASK=m, FIELD=f)
2060 6 : ig(:, :) = UNPACK(message(off + 33:off + 38), MASK=m, FIELD=f)
2061 6 : bf = message(off + 39)
2062 6 : bu = message(off + 40)
2063 6 : IF (bf .GT. 0) THEN
2064 4 : lbf = .TRUE.
2065 : ELSE
2066 2 : lbf = .FALSE.
2067 : END IF
2068 : CALL helium_env(k)%helium%rng_stream_gaussian%set(bg=bg, cg=cg, ig=ig, &
2069 12 : buffer=bu, buffer_filled=lbf)
2070 : END DO
2071 : END IF
2072 :
2073 : NULLIFY (message)
2074 :
2075 6 : RETURN
2076 6 : 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 2 : 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 2 : NULLIFY (logger)
2097 2 : logger => cp_get_default_logger()
2098 : CALL section_vals_val_get( &
2099 : helium%input, &
2100 : "MOTION%PINT%HELIUM%RDF%SOLUTE_HE", &
2101 2 : l_val=helium%rdf_sol_he)
2102 : CALL section_vals_val_get( &
2103 : helium%input, &
2104 : "MOTION%PINT%HELIUM%RDF%HE_HE", &
2105 2 : l_val=helium%rdf_he_he)
2106 :
2107 : ! Prevent He_He Rdfs for a single he atom:
2108 2 : IF (helium%atoms <= 1) THEN
2109 0 : helium%rdf_he_he = .FALSE.
2110 : END IF
2111 :
2112 2 : helium%rdf_num = 0
2113 2 : IF (helium%rdf_sol_he) THEN
2114 2 : IF (helium%solute_present) THEN
2115 : ! get number of centers from solute to store solute positions
2116 0 : ALLOCATE (helium%rdf_centers(helium%beads, helium%solute_atoms*3))
2117 0 : helium%rdf_centers(:, :) = 0.0_dp
2118 0 : helium%rdf_num = helium%solute_atoms
2119 : ELSE
2120 2 : helium%rdf_sol_he = .FALSE.
2121 : END IF
2122 : END IF
2123 :
2124 2 : IF (helium%rdf_he_he) THEN
2125 2 : helium%rdf_num = helium%rdf_num + 1
2126 : END IF
2127 :
2128 : ! set the flag for RDF and either proceed or return
2129 2 : IF (helium%rdf_num > 0) THEN
2130 2 : helium%rdf_present = .TRUE.
2131 : ELSE
2132 0 : helium%rdf_present = .FALSE.
2133 : err_str = "HELIUM RDFs requested, but no valid choice of "// &
2134 0 : "parameters specified. No RDF will be computed."
2135 0 : CPWARN(err_str)
2136 0 : 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 2 : explicit=explicit)
2144 2 : 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 0 : 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 2 : explicit=explicit)
2156 2 : IF (explicit) THEN
2157 : ! use the droplet radius
2158 0 : 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 0 : helium%rdf_maxr = helium%droplet_radius*2.0_dp
2165 : ELSE
2166 0 : 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 2 : SELECT CASE (helium%cell_shape)
2172 : CASE (helium_cell_shape_cube)
2173 0 : helium%rdf_maxr = helium%cell_size/2.0_dp
2174 : CASE (helium_cell_shape_octahedron)
2175 2 : helium%rdf_maxr = helium%cell_size*SQRT(3.0_dp)/4.0_dp
2176 : CASE DEFAULT
2177 0 : helium%rdf_maxr = 0.0_dp
2178 0 : WRITE (stmp, *) helium%cell_shape
2179 0 : err_str = "cell shape unknown ("//TRIM(ADJUSTL(stmp))//")"
2180 2 : 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 2 : i_val=helium%rdf_nbin)
2190 2 : 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 2 : i_val=helium%rdf_iweight)
2197 :
2198 : ! allocate and initialize memory for RDF storage
2199 2 : ii = helium%rdf_num
2200 2 : ij = helium%rdf_nbin
2201 8 : ALLOCATE (helium%rdf_inst(ii, ij))
2202 8 : ALLOCATE (helium%rdf_accu(ii, ij))
2203 8 : ALLOCATE (helium%rdf_rstr(ii, ij))
2204 1002 : helium%rdf_inst(:, :) = 0.0_dp
2205 1002 : helium%rdf_accu(:, :) = 0.0_dp
2206 1002 : helium%rdf_rstr(:, :) = 0.0_dp
2207 :
2208 : RETURN
2209 2 : 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 2 : 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 2 : LOGICAL, DIMENSION(:, :), POINTER :: m
2229 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: message
2230 2 : 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 2 : explicit=explicit)
2235 2 : IF (explicit) THEN
2236 0 : NULLIFY (message)
2237 : CALL section_vals_val_get(helium_env(1)%helium%input, &
2238 : "MOTION%PINT%HELIUM%AVERAGES%RDF", &
2239 0 : r_vals=message)
2240 0 : msglen = SIZE(message)
2241 0 : itmp = SIZE(helium_env(1)%helium%rdf_rstr)
2242 0 : ltmp = (msglen .EQ. itmp)
2243 0 : IF (.NOT. ltmp) THEN
2244 0 : stmp1 = ""
2245 0 : WRITE (stmp1, *) msglen
2246 0 : stmp2 = ""
2247 0 : 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 0 : TRIM(ADJUSTL(stmp2))//")."
2252 0 : CPABORT(err_str)
2253 : END IF
2254 : ELSE
2255 : RETURN
2256 : END IF
2257 :
2258 0 : ii = helium_env(1)%helium%rdf_num
2259 0 : ij = helium_env(1)%helium%rdf_nbin
2260 : NULLIFY (m, f)
2261 0 : ALLOCATE (m(ii, ij))
2262 0 : ALLOCATE (f(ii, ij))
2263 0 : m(:, :) = .TRUE.
2264 0 : f(:, :) = 0.0_dp
2265 :
2266 0 : DO k = 1, SIZE(helium_env)
2267 0 : 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 0 : i_val=helium_env(k)%helium%rdf_iweight)
2271 : END DO
2272 :
2273 0 : DEALLOCATE (f, m)
2274 : NULLIFY (message)
2275 :
2276 0 : RETURN
2277 2 : 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 31 : SUBROUTINE helium_rdf_release(helium)
2286 :
2287 : TYPE(helium_solvent_type), POINTER :: helium
2288 :
2289 31 : IF (helium%rdf_present) THEN
2290 :
2291 0 : DEALLOCATE ( &
2292 : helium%rdf_rstr, &
2293 0 : helium%rdf_accu, &
2294 2 : helium%rdf_inst)
2295 :
2296 : NULLIFY ( &
2297 : helium%rdf_rstr, &
2298 2 : helium%rdf_accu, &
2299 2 : helium%rdf_inst)
2300 :
2301 2 : IF (helium%rdf_sol_he) THEN
2302 0 : DEALLOCATE (helium%rdf_centers)
2303 0 : NULLIFY (helium%rdf_centers)
2304 : END IF
2305 :
2306 : END IF
2307 :
2308 31 : 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 124 : 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 62 : is_active = .FALSE.
2339 62 : NULLIFY (logger)
2340 124 : logger => cp_get_default_logger()
2341 :
2342 : ! if the printkey is on at this runlevel consider prop to be active
2343 62 : NULLIFY (print_key)
2344 62 : input_path = "MOTION%PINT%HELIUM%PRINT%"//TRIM(ADJUSTL(prop))
2345 : print_key => section_vals_get_subs_vals( &
2346 : helium%input, &
2347 62 : input_path)
2348 : is_on = cp_printkey_is_on( &
2349 : iteration_info=logger%iter_info, &
2350 62 : print_key=print_key)
2351 62 : IF (is_on) THEN
2352 62 : 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 62 : is_active = .FALSE.
2359 62 : NULLIFY (section)
2360 62 : input_path = "MOTION%PINT%HELIUM%"//TRIM(ADJUSTL(prop))
2361 : section => section_vals_get_subs_vals( &
2362 : helium%input, &
2363 62 : input_path)
2364 62 : CALL section_vals_get(section, explicit=explicit)
2365 62 : IF (explicit) THEN
2366 : ! control section explicitly present, continue checking
2367 : CALL section_vals_val_get( &
2368 : section, &
2369 : "_SECTION_PARAMETERS_", &
2370 4 : l_val=is_on)
2371 4 : IF (is_on) THEN
2372 : ! control section is explicit and on, activate the property
2373 4 : is_active = .TRUE.
2374 : ! activate the corresponding print_level as well
2375 4 : print_level = logger%iter_info%print_level
2376 : CALL section_vals_val_set( &
2377 : print_key, &
2378 : "_SECTION_PARAMETERS_", &
2379 4 : i_val=print_level)
2380 : END IF
2381 : END IF
2382 :
2383 : RETURN
2384 62 : 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 2 : SUBROUTINE helium_rho_property_init(helium)
2393 :
2394 : TYPE(helium_solvent_type), POINTER :: helium
2395 :
2396 : INTEGER :: nc
2397 :
2398 12 : ALLOCATE (helium%rho_property(rho_num))
2399 :
2400 2 : helium%rho_property(rho_atom_number)%name = 'Atom number density'
2401 2 : nc = 1
2402 2 : helium%rho_property(rho_atom_number)%num_components = nc
2403 2 : ALLOCATE (helium%rho_property(rho_atom_number)%filename_suffix(nc))
2404 2 : ALLOCATE (helium%rho_property(rho_atom_number)%component_name(nc))
2405 2 : ALLOCATE (helium%rho_property(rho_atom_number)%component_index(nc))
2406 2 : helium%rho_property(rho_atom_number)%filename_suffix(1) = 'an'
2407 2 : helium%rho_property(rho_atom_number)%component_name(1) = ''
2408 4 : helium%rho_property(rho_atom_number)%component_index(:) = 0
2409 :
2410 2 : helium%rho_property(rho_projected_area)%name = 'Projected area squared density, A*A(r)'
2411 2 : nc = 3
2412 2 : helium%rho_property(rho_projected_area)%num_components = nc
2413 2 : ALLOCATE (helium%rho_property(rho_projected_area)%filename_suffix(nc))
2414 2 : ALLOCATE (helium%rho_property(rho_projected_area)%component_name(nc))
2415 2 : ALLOCATE (helium%rho_property(rho_projected_area)%component_index(nc))
2416 2 : helium%rho_property(rho_projected_area)%filename_suffix(1) = 'pa_x'
2417 2 : helium%rho_property(rho_projected_area)%filename_suffix(2) = 'pa_y'
2418 2 : helium%rho_property(rho_projected_area)%filename_suffix(3) = 'pa_z'
2419 2 : helium%rho_property(rho_projected_area)%component_name(1) = 'component x'
2420 2 : helium%rho_property(rho_projected_area)%component_name(2) = 'component y'
2421 2 : helium%rho_property(rho_projected_area)%component_name(3) = 'component z'
2422 8 : helium%rho_property(rho_projected_area)%component_index(:) = 0
2423 :
2424 2 : helium%rho_property(rho_winding_number)%name = 'Winding number squared density, W*W(r)'
2425 2 : nc = 3
2426 2 : helium%rho_property(rho_winding_number)%num_components = nc
2427 2 : ALLOCATE (helium%rho_property(rho_winding_number)%filename_suffix(nc))
2428 2 : ALLOCATE (helium%rho_property(rho_winding_number)%component_name(nc))
2429 2 : ALLOCATE (helium%rho_property(rho_winding_number)%component_index(nc))
2430 2 : helium%rho_property(rho_winding_number)%filename_suffix(1) = 'wn_x'
2431 2 : helium%rho_property(rho_winding_number)%filename_suffix(2) = 'wn_y'
2432 2 : helium%rho_property(rho_winding_number)%filename_suffix(3) = 'wn_z'
2433 2 : helium%rho_property(rho_winding_number)%component_name(1) = 'component x'
2434 2 : helium%rho_property(rho_winding_number)%component_name(2) = 'component y'
2435 2 : helium%rho_property(rho_winding_number)%component_name(3) = 'component z'
2436 8 : helium%rho_property(rho_winding_number)%component_index(:) = 0
2437 :
2438 2 : helium%rho_property(rho_winding_cycle)%name = 'Winding number squared density, W^2(r)'
2439 2 : nc = 3
2440 2 : helium%rho_property(rho_winding_cycle)%num_components = nc
2441 2 : ALLOCATE (helium%rho_property(rho_winding_cycle)%filename_suffix(nc))
2442 2 : ALLOCATE (helium%rho_property(rho_winding_cycle)%component_name(nc))
2443 2 : ALLOCATE (helium%rho_property(rho_winding_cycle)%component_index(nc))
2444 2 : helium%rho_property(rho_winding_cycle)%filename_suffix(1) = 'wc_x'
2445 2 : helium%rho_property(rho_winding_cycle)%filename_suffix(2) = 'wc_y'
2446 2 : helium%rho_property(rho_winding_cycle)%filename_suffix(3) = 'wc_z'
2447 2 : helium%rho_property(rho_winding_cycle)%component_name(1) = 'component x'
2448 2 : helium%rho_property(rho_winding_cycle)%component_name(2) = 'component y'
2449 2 : helium%rho_property(rho_winding_cycle)%component_name(3) = 'component z'
2450 8 : helium%rho_property(rho_winding_cycle)%component_index(:) = 0
2451 :
2452 2 : helium%rho_property(rho_moment_of_inertia)%name = 'Moment of inertia'
2453 2 : nc = 3
2454 2 : helium%rho_property(rho_moment_of_inertia)%num_components = nc
2455 2 : ALLOCATE (helium%rho_property(rho_moment_of_inertia)%filename_suffix(nc))
2456 2 : ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_name(nc))
2457 2 : ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_index(nc))
2458 2 : helium%rho_property(rho_moment_of_inertia)%filename_suffix(1) = 'mi_x'
2459 2 : helium%rho_property(rho_moment_of_inertia)%filename_suffix(2) = 'mi_y'
2460 2 : helium%rho_property(rho_moment_of_inertia)%filename_suffix(3) = 'mi_z'
2461 2 : helium%rho_property(rho_moment_of_inertia)%component_name(1) = 'component x'
2462 2 : helium%rho_property(rho_moment_of_inertia)%component_name(2) = 'component y'
2463 2 : helium%rho_property(rho_moment_of_inertia)%component_name(3) = 'component z'
2464 8 : helium%rho_property(rho_moment_of_inertia)%component_index(:) = 0
2465 :
2466 12 : helium%rho_property(:)%is_calculated = .FALSE.
2467 :
2468 2 : 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 2 : 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 2 : CALL helium_rho_property_init(helium)
2485 :
2486 2 : 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 2 : l_val=ltmp)
2493 2 : IF (ltmp) THEN
2494 2 : helium%rho_property(rho_atom_number)%is_calculated = .TRUE.
2495 2 : helium%rho_num_act = helium%rho_num_act + 1
2496 2 : 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 2 : l_val=ltmp)
2504 2 : IF (ltmp) THEN
2505 0 : helium%rho_property(rho_projected_area)%is_calculated = .TRUE.
2506 0 : DO ii = 1, helium%rho_property(rho_projected_area)%num_components
2507 0 : helium%rho_num_act = helium%rho_num_act + 1
2508 0 : 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 2 : l_val=ltmp)
2517 2 : IF (ltmp) THEN
2518 0 : helium%rho_property(rho_winding_number)%is_calculated = .TRUE.
2519 0 : DO ii = 1, helium%rho_property(rho_winding_number)%num_components
2520 0 : helium%rho_num_act = helium%rho_num_act + 1
2521 0 : 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 2 : l_val=ltmp)
2530 2 : IF (ltmp) THEN
2531 0 : helium%rho_property(rho_winding_cycle)%is_calculated = .TRUE.
2532 0 : DO ii = 1, helium%rho_property(rho_winding_cycle)%num_components
2533 0 : helium%rho_num_act = helium%rho_num_act + 1
2534 0 : 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 2 : l_val=ltmp)
2543 2 : IF (ltmp) THEN
2544 0 : helium%rho_property(rho_moment_of_inertia)%is_calculated = .TRUE.
2545 0 : DO ii = 1, helium%rho_property(rho_moment_of_inertia)%num_components
2546 0 : helium%rho_num_act = helium%rho_num_act + 1
2547 0 : 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 2 : helium%rho_maxr = helium%cell_size
2553 : CALL section_vals_val_get( &
2554 : helium%input, &
2555 : "MOTION%PINT%HELIUM%RHO%NBIN", &
2556 2 : i_val=helium%rho_nbin)
2557 2 : helium%rho_delr = helium%rho_maxr/REAL(helium%rho_nbin, dp)
2558 :
2559 : ! check for optional estimators based on winding paths
2560 2 : 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 2 : explicit=explicit)
2565 2 : IF (explicit) THEN
2566 0 : 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 0 : i_vals=helium%rho_min_len_wdg_vals)
2571 0 : itmp = SIZE(helium%rho_min_len_wdg_vals)
2572 0 : IF (itmp .GT. 0) THEN
2573 0 : helium%rho_num_min_len_wdg = itmp
2574 0 : 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 2 : 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 2 : explicit=explicit)
2584 2 : IF (explicit) THEN
2585 0 : 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 0 : i_vals=helium%rho_min_len_non_vals)
2590 0 : itmp = SIZE(helium%rho_min_len_non_vals)
2591 0 : IF (itmp .GT. 0) THEN
2592 0 : helium%rho_num_min_len_non = itmp
2593 0 : 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 2 : 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 2 : explicit=explicit)
2603 2 : IF (explicit) THEN
2604 0 : 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 0 : i_vals=helium%rho_min_len_all_vals)
2609 0 : itmp = SIZE(helium%rho_min_len_all_vals)
2610 0 : IF (itmp .GT. 0) THEN
2611 0 : helium%rho_num_min_len_all = itmp
2612 0 : 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 2 : i_val=helium%rho_iweight)
2621 :
2622 : ! allocate and initialize memory for density storage
2623 2 : itmp = helium%rho_nbin
2624 2 : jtmp = helium%rho_num_act
2625 12 : ALLOCATE (helium%rho_inst(jtmp, itmp, itmp, itmp))
2626 12 : ALLOCATE (helium%rho_accu(jtmp, itmp, itmp, itmp))
2627 12 : ALLOCATE (helium%rho_rstr(jtmp, itmp, itmp, itmp))
2628 10 : ALLOCATE (helium%rho_incr(jtmp, helium%atoms, helium%beads))
2629 :
2630 3252 : helium%rho_incr(:, :, :) = 0.0_dp
2631 4222 : helium%rho_inst(:, :, :, :) = 0.0_dp
2632 4222 : helium%rho_accu(:, :, :, :) = 0.0_dp
2633 4222 : helium%rho_rstr(:, :, :, :) = 0.0_dp
2634 :
2635 2 : RETURN
2636 16 : 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 2 : 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 2 : LOGICAL, DIMENSION(:, :, :, :), POINTER :: m
2655 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: message
2656 2 : 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 2 : explicit=explicit)
2661 2 : IF (explicit) THEN
2662 0 : NULLIFY (message)
2663 : CALL section_vals_val_get(helium_env(1)%helium%input, &
2664 : "MOTION%PINT%HELIUM%AVERAGES%RHO", &
2665 0 : r_vals=message)
2666 0 : msglen = SIZE(message)
2667 0 : itmp = SIZE(helium_env(1)%helium%rho_rstr)
2668 0 : ltmp = (msglen .EQ. itmp)
2669 0 : IF (.NOT. ltmp) THEN
2670 0 : stmp1 = ""
2671 0 : WRITE (stmp1, *) msglen
2672 0 : stmp2 = ""
2673 0 : 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 0 : TRIM(ADJUSTL(stmp2))//")."
2678 0 : CPABORT(err_str)
2679 : END IF
2680 : ELSE
2681 : RETURN
2682 : END IF
2683 :
2684 0 : itmp = helium_env(1)%helium%rho_nbin
2685 : NULLIFY (m, f)
2686 0 : ALLOCATE (m(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp))
2687 0 : ALLOCATE (f(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp))
2688 0 : m(:, :, :, :) = .TRUE.
2689 0 : f(:, :, :, :) = 0.0_dp
2690 :
2691 0 : DO k = 1, SIZE(helium_env)
2692 0 : helium_env(k)%helium%rho_rstr(:, :, :, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f)
2693 : END DO
2694 :
2695 0 : DEALLOCATE (f, m)
2696 : NULLIFY (message)
2697 :
2698 0 : RETURN
2699 2 : 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 21 : 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 21 : NULLIFY (my_f_env, my_subsys, my_particles)
2722 : CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
2723 21 : f_env=my_f_env)
2724 21 : CALL force_env_get(force_env=my_f_env%force_env, subsys=my_subsys)
2725 21 : CALL cp_subsys_get(my_subsys, particles=my_particles)
2726 21 : CALL f_env_rm_defaults(my_f_env)
2727 :
2728 21 : natoms = helium%solute_atoms
2729 21 : NULLIFY (helium%solute_element)
2730 42 : ALLOCATE (helium%solute_element(natoms))
2731 :
2732 : ! find out how many different atomic types are there
2733 21 : helium%enum = 0
2734 84 : DO iatom = 1, natoms
2735 : CALL get_atomic_kind(my_particles%els(iatom)%atomic_kind, &
2736 : mass=mass, &
2737 84 : element_symbol=helium%solute_element(iatom))
2738 : END DO
2739 :
2740 21 : 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 42 : 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 21 : NULLIFY (my_f_env, my_cell)
2766 : CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
2767 21 : f_env=my_f_env)
2768 21 : CALL force_env_get(force_env=my_f_env%force_env, cell=my_cell)
2769 21 : CALL f_env_rm_defaults(my_f_env)
2770 :
2771 21 : CALL get_cell(my_cell, orthorhombic=my_orthorhombic)
2772 21 : IF (.NOT. my_orthorhombic) THEN
2773 0 : CPABORT("Helium solvent not implemented for non-orthorhombic cells.")
2774 : ELSE
2775 21 : helium%solute_cell => my_cell
2776 : END IF
2777 :
2778 21 : RETURN
2779 : END SUBROUTINE helium_set_solute_cell
2780 :
2781 : END MODULE helium_methods
|
