(git:4d70172)
Loading...
Searching...
No Matches
pint_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Methods to performs a path integral run
10!> \author fawzi
11!> \par History
12!> 02.2005 created [fawzi]
13!> 11.2006 modified so it might actually work [hforbert]
14!> 10.2015 added RPMD propagator
15!> 10.2015 added exact harmonic integrator [Felix Uhl]
16!> \note quick & dirty rewrite of my python program
17! **************************************************************************************************
18MODULE pint_methods
19
20 USE atomic_kind_list_types, ONLY: atomic_kind_list_type
21 USE atomic_kind_types, ONLY: atomic_kind_type,&
22 get_atomic_kind
23 USE bibliography, ONLY: brieuc2016,&
24 ceriotti2010,&
25 ceriotti2012,&
26 shiga2022,&
27 cite_reference
28 USE cell_types, ONLY: cell_type
29 USE constraint, ONLY: rattle_control,&
30 shake_control,&
31 shake_update_targets
32 USE constraint_util, ONLY: getold
33 USE cp_external_control, ONLY: external_control
34 USE cp_log_handling, ONLY: cp_get_default_logger,&
35 cp_logger_get_default_io_unit,&
36 cp_logger_type,&
37 cp_to_string
38 USE cp_output_handling, ONLY: cp_add_iter_level,&
39 cp_iterate,&
40 cp_p_file,&
41 cp_print_key_finished_output,&
42 cp_print_key_should_output,&
43 cp_print_key_unit_nr,&
44 cp_rm_iter_level
45 USE cp_subsys_types, ONLY: cp_subsys_get,&
46 cp_subsys_type
47 USE cp_units, ONLY: cp_unit_from_cp2k,&
48 cp_unit_to_cp2k
49 USE distribution_1d_types, ONLY: distribution_1d_type
50 USE f77_interface, ONLY: f_env_add_defaults,&
51 f_env_rm_defaults,&
52 f_env_type
53 USE force_env_types, ONLY: force_env_get
54 USE gle_system_dynamics, ONLY: gle_cholesky_stab,&
55 gle_matrix_exp,&
56 restart_gle
57 USE gle_system_types, ONLY: gle_dealloc,&
58 gle_init,&
59 gle_thermo_create
60 USE global_types, ONLY: global_environment_type
61 USE helium_interactions, ONLY: helium_intpot_scan
62 USE helium_io, ONLY: helium_write_cubefile
63 USE helium_methods, ONLY: helium_create,&
64 helium_init,&
65 helium_release
66 USE helium_sampling, ONLY: helium_do_run,&
67 helium_step
68 USE helium_types, ONLY: helium_solvent_p_type
69 USE input_constants, ONLY: integrate_exact,&
70 integrate_numeric,&
71 propagator_bcmd,&
72 propagator_cmd,&
73 propagator_rpmd,&
74 transformation_normal,&
75 transformation_stage
76 USE input_cp2k_restarts, ONLY: write_restart
77 USE input_section_types, ONLY: &
78 section_type, section_vals_add_values, section_vals_get, section_vals_get_subs_vals, &
79 section_vals_release, section_vals_retain, section_vals_type, section_vals_val_get, &
80 section_vals_val_set, section_vals_val_unset
81 USE kinds, ONLY: default_path_length,&
82 default_string_length,&
83 dp
84 USE machine, ONLY: m_flush,&
85 m_walltime
86 USE mathconstants, ONLY: twopi
87 USE mathlib, ONLY: gcd
88 USE message_passing, ONLY: mp_comm_self,&
89 mp_para_env_type
90 USE molecule_kind_list_types, ONLY: molecule_kind_list_type
91 USE molecule_kind_types, ONLY: molecule_kind_type
92 USE molecule_list_types, ONLY: molecule_list_type
93 USE molecule_types, ONLY: global_constraint_type,&
94 molecule_type
95 USE parallel_rng_types, ONLY: gaussian,&
96 rng_stream_type
97 USE particle_list_types, ONLY: particle_list_type
98 USE particle_types, ONLY: particle_type
99 USE pint_gle, ONLY: pint_calc_gle_energy,&
100 pint_gle_init,&
101 pint_gle_step
102 USE pint_io, ONLY: pint_write_action,&
103 pint_write_centroids,&
104 pint_write_com,&
105 pint_write_ener,&
106 pint_write_line,&
107 pint_write_rgyr,&
108 pint_write_step_info,&
109 pint_write_trajectory
110 USE pint_normalmode, ONLY: normalmode_calc_uf_h,&
111 normalmode_env_create,&
112 normalmode_init_masses,&
113 normalmode_release
114 USE pint_piglet, ONLY: pint_calc_piglet_energy,&
115 pint_piglet_create,&
116 pint_piglet_init,&
117 pint_piglet_release,&
118 pint_piglet_step
119 USE pint_pile, ONLY: pint_calc_pile_energy,&
120 pint_pile_init,&
121 pint_pile_release,&
122 pint_pile_step
123 USE pint_public, ONLY: pint_levy_walk
124 USE pint_qtb, ONLY: pint_calc_qtb_energy,&
125 pint_qtb_init,&
126 pint_qtb_release,&
127 pint_qtb_step
128 USE pint_staging, ONLY: staging_calc_uf_h,&
129 staging_env_create,&
130 staging_init_masses,&
131 staging_release
132 USE pint_transformations, ONLY: pint_f2uf,&
133 pint_u2x,&
134 pint_x2u
135 USE pint_types, ONLY: &
136 e_conserved_id, e_kin_thermo_id, e_kin_virial_id, e_potential_id, pint_env_type, &
137 thermostat_gle, thermostat_none, thermostat_nose, thermostat_piglet, thermostat_pile, &
138 thermostat_qtb
139 USE replica_methods, ONLY: rep_env_calc_e_f,&
140 rep_env_create
141 USE replica_types, ONLY: rep_env_release,&
142 replica_env_type
143 USE simpar_types, ONLY: create_simpar_type,&
144 release_simpar_type
145#include "../base/base_uses.f90"
146
147 IMPLICIT NONE
148 PRIVATE
149
150 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .true.
151 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_methods'
152
153 PUBLIC :: do_pint_run
154
155CONTAINS
156
157! ***************************************************************************
158!> \brief Create a path integral environment
159!> \param pint_env ...
160!> \param input ...
161!> \param input_declaration ...
162!> \param para_env ...
163!> \par History
164!> Fixed some bugs [hforbert]
165!> Added normal mode transformation [hforbert]
166!> 10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl]
167!> 10.2018 Added centroid constraints [cschran+rperez]
168!> 10.2021 Added beadwise constraints [lduran]
169!> \author fawzi
170!> \note Might return an unassociated pointer in parallel on the processors
171!> that are not needed.
172! **************************************************************************************************
173 SUBROUTINE pint_create(pint_env, input, input_declaration, para_env)
174
175 TYPE(pint_env_type), INTENT(OUT) :: pint_env
176 TYPE(section_vals_type), POINTER :: input
177 TYPE(section_type), POINTER :: input_declaration
178 TYPE(mp_para_env_type), POINTER :: para_env
179
180 CHARACTER(len=*), PARAMETER :: routineN = 'pint_create'
181
182 CHARACTER(len=2*default_string_length) :: msg
183 CHARACTER(len=default_path_length) :: output_file_name, project_name
184 INTEGER :: handle, iat, ibead, icont, idim, idir, &
185 ierr, ig, itmp, nrep, prep, stat
186 LOGICAL :: explicit, ltmp
187 REAL(kind=dp) :: dt, mass, omega
188 TYPE(cp_subsys_type), POINTER :: subsys
189 TYPE(f_env_type), POINTER :: f_env
190 TYPE(global_constraint_type), POINTER :: gci
191 TYPE(particle_list_type), POINTER :: particles
192 TYPE(replica_env_type), POINTER :: rep_env
193 TYPE(section_vals_type), POINTER :: constraint_section, gle_section, nose_section, &
194 piglet_section, pile_section, pint_section, qtb_section, transform_section
195
196 CALL timeset(routinen, handle)
197
198 NULLIFY (f_env, subsys, particles, nose_section, gle_section, gci)
199
200 cpassert(ASSOCIATED(input))
201 cpassert(input%ref_count > 0)
202 NULLIFY (rep_env)
203 pint_section => section_vals_get_subs_vals(input, "MOTION%PINT")
204 CALL section_vals_val_get(pint_section, "p", i_val=nrep)
205 CALL section_vals_val_get(pint_section, "proc_per_replica", &
206 i_val=prep)
207 ! Maybe let the user have his/her way as long as prep is
208 ! within the bounds of number of CPUs??
209 IF ((prep < 1) .OR. (prep > para_env%num_pe) .OR. &
210 (mod(prep*nrep, para_env%num_pe) /= 0)) THEN
211 prep = para_env%num_pe/gcd(para_env%num_pe, nrep)
212 IF (para_env%is_source()) THEN
213 WRITE (unit=msg, fmt=*) "PINT WARNING: Adjusting number of processors per replica to ", prep
214 cpwarn(msg)
215 END IF
216 END IF
217
218 ! replica_env modifies the global input structure - which is wrong - one
219 ! of the side effects is the inifite adding of the -r-N string to the
220 ! project name and the output file name, which corrupts restart files.
221 ! For now: save the project name and output file name and restore them
222 ! after the rep_env_create has executed - the initialization of the
223 ! replicas will run correctly anyways.
224 ! TODO: modify rep_env so that it behaves better
225 CALL section_vals_val_get(input, "GLOBAL%PROJECT_NAME", c_val=project_name)
226 CALL section_vals_val_get(input, "GLOBAL%OUTPUT_FILE_NAME", c_val=output_file_name)
227 CALL rep_env_create(rep_env, para_env=para_env, input=input, &
228 input_declaration=input_declaration, nrep=nrep, prep=prep, row_force=.true.)
229 CALL section_vals_val_set(input, "GLOBAL%PROJECT_NAME", c_val=trim(project_name))
230 IF (len_trim(output_file_name) > 0) THEN
231 CALL section_vals_val_set(input, "GLOBAL%OUTPUT_FILE_NAME", c_val=trim(output_file_name))
232 ELSE
233 CALL section_vals_val_unset(input, "GLOBAL%OUTPUT_FILE_NAME")
234 END IF
235 IF (.NOT. ASSOCIATED(rep_env)) RETURN
236
237 NULLIFY (pint_env%logger)
238 pint_env%logger => cp_get_default_logger()
239 CALL cp_add_iter_level(pint_env%logger%iter_info, "PINT")
240
241 NULLIFY (pint_env%replicas, pint_env%input, pint_env%staging_env, &
242 pint_env%normalmode_env, pint_env%propagator)
243 pint_env%p = nrep
244 pint_env%replicas => rep_env
245 pint_env%ndim = rep_env%ndim
246 pint_env%input => input
247
248 CALL section_vals_retain(pint_env%input)
249
250 ! get first step, last step, number of steps, etc
251 CALL section_vals_val_get(input, "MOTION%PINT%ITERATION", &
252 i_val=itmp)
253 pint_env%first_step = itmp
254 CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
255 explicit=explicit)
256 IF (explicit) THEN
257 CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", &
258 i_val=itmp)
259 pint_env%last_step = itmp
260 pint_env%num_steps = pint_env%last_step - pint_env%first_step
261 ELSE
262 CALL section_vals_val_get(input, "MOTION%PINT%NUM_STEPS", &
263 i_val=itmp)
264 pint_env%num_steps = itmp
265 pint_env%last_step = pint_env%first_step + pint_env%num_steps
266 END IF
267
268 CALL section_vals_val_get(pint_section, "DT", &
269 r_val=pint_env%dt)
270 pint_env%t = pint_env%first_step*pint_env%dt
271
272 CALL section_vals_val_get(pint_section, "nrespa", i_val=pint_env%nrespa)
273 CALL section_vals_val_get(pint_section, "Temp", r_val=pint_env%kT)
274 CALL section_vals_val_get(pint_section, "T_TOL", &
275 r_val=pint_env%t_tol)
276
277 CALL section_vals_val_get(pint_section, "HARM_INT", i_val=pint_env%harm_integrator)
278
279 ALLOCATE (pint_env%propagator)
280 CALL section_vals_val_get(pint_section, "propagator", &
281 i_val=pint_env%propagator%prop_kind)
282 !Initialize simulation temperature depending on the propagator
283 !As well as the scaling factor for the physical potential
284 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
285 pint_env%propagator%temp_phys2sim = real(pint_env%p, dp)
286 pint_env%propagator%physpotscale = 1.0_dp
287 ELSE
288 pint_env%propagator%temp_phys2sim = 1.0_dp
289 pint_env%propagator%physpotscale = 1.0_dp/real(pint_env%p, dp)
290 END IF
291 pint_env%propagator%temp_sim2phys = 1.0_dp/pint_env%propagator%temp_phys2sim
292 pint_env%kT = pint_env%kT*pint_env%propagator%temp_phys2sim
293
294 CALL section_vals_val_get(pint_section, "transformation", &
295 i_val=pint_env%transform)
296
297 IF ((pint_env%propagator%prop_kind == propagator_cmd) .AND. &
298 (pint_env%transform /= transformation_normal)) THEN
299 cpabort("CMD propagator without normal modes not implemented!")
300 END IF
301
302 IF ((pint_env%propagator%prop_kind == propagator_bcmd) .AND. &
303 (pint_env%transform /= transformation_normal)) THEN
304 cpabort("BCMD propagator without normal modes not implemented!")
305 END IF
306
307 NULLIFY (pint_env%tx, pint_env%tv, pint_env%tv_t, pint_env%tv_old, pint_env%tv_new, pint_env%tf)
308
309 pint_env%nnos = 0
310 pint_env%pimd_thermostat = thermostat_none
311 nose_section => section_vals_get_subs_vals(input, "MOTION%PINT%NOSE")
312 CALL section_vals_get(nose_section, explicit=explicit)
313 IF (explicit) THEN
314 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
315 cpabort("RPMD propagator with Nose-thermostat not implemented!")
316 END IF
317 IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN
318 cpabort("BCMD propagator with Nose-thermostat not implemented!")
319 END IF
320 CALL section_vals_val_get(nose_section, "nnos", i_val=pint_env%nnos)
321 IF (pint_env%nnos > 0) THEN
322 pint_env%pimd_thermostat = thermostat_nose
323 ALLOCATE ( &
324 pint_env%tx(pint_env%nnos, pint_env%p, pint_env%ndim), &
325 pint_env%tv(pint_env%nnos, pint_env%p, pint_env%ndim), &
326 pint_env%tv_t(pint_env%nnos, pint_env%p, pint_env%ndim), &
327 pint_env%tv_old(pint_env%nnos, pint_env%p, pint_env%ndim), &
328 pint_env%tv_new(pint_env%nnos, pint_env%p, pint_env%ndim), &
329 pint_env%tf(pint_env%nnos, pint_env%p, pint_env%ndim))
330 pint_env%tx = 0._dp
331 pint_env%tv = 0._dp
332 pint_env%tv_t = 0._dp
333 pint_env%tv_old = 0._dp
334 pint_env%tv_new = 0._dp
335 pint_env%tf = 0._dp
336 END IF
337 END IF
338
339 pint_env%beta = 1._dp/(pint_env%kT*pint_env%propagator%temp_sim2phys)
340!TODO
341! v_tol not in current input structure
342! should also probably be part of nose_section
343! CALL section_vals_val_get(transform_section,"v_tol_nose",r_val=pint_env%v_tol)
344!MK ... but we have to initialise v_tol
345 pint_env%v_tol = 0.0_dp ! to be fixed
346
347 pint_env%randomG = rng_stream_type( &
348 name="pint_randomG", &
349 distribution_type=gaussian, &
350 extended_precision=.true.)
351
352 ALLOCATE (pint_env%e_pot_bead(pint_env%p))
353 pint_env%e_pot_bead = 0._dp
354 pint_env%e_pot_h = 0._dp
355 pint_env%e_kin_beads = 0._dp
356 pint_env%e_pot_t = 0._dp
357 pint_env%e_gle = 0._dp
358 pint_env%e_pile = 0._dp
359 pint_env%e_piglet = 0._dp
360 pint_env%e_qtb = 0._dp
361 pint_env%e_kin_t = 0._dp
362 pint_env%energy(:) = 0.0_dp
363
364!TODO: rearrange to use standard nose hoover chain functions/data types
365
366 ALLOCATE ( &
367 pint_env%x(pint_env%p, pint_env%ndim), &
368 pint_env%v(pint_env%p, pint_env%ndim), &
369 pint_env%f(pint_env%p, pint_env%ndim), &
370 pint_env%external_f(pint_env%p, pint_env%ndim), &
371 pint_env%ux(pint_env%p, pint_env%ndim), &
372 pint_env%ux_t(pint_env%p, pint_env%ndim), &
373 pint_env%uv(pint_env%p, pint_env%ndim), &
374 pint_env%uv_t(pint_env%p, pint_env%ndim), &
375 pint_env%uv_new(pint_env%p, pint_env%ndim), &
376 pint_env%uf(pint_env%p, pint_env%ndim), &
377 pint_env%uf_h(pint_env%p, pint_env%ndim), &
378 pint_env%centroid(pint_env%ndim), &
379 pint_env%rtmp_ndim(pint_env%ndim), &
380 pint_env%rtmp_natom(pint_env%ndim/3), &
381 stat=stat)
382 cpassert(stat == 0)
383 pint_env%x = 0._dp
384 pint_env%v = 0._dp
385 pint_env%f = 0._dp
386 pint_env%external_f = 0._dp
387 pint_env%ux = 0._dp
388 pint_env%ux_t = 0._dp
389 pint_env%uv = 0._dp
390 pint_env%uv_t = 0._dp
391 pint_env%uv_new = 0._dp
392 pint_env%uf = 0._dp
393 pint_env%uf_h = 0._dp
394 pint_env%centroid(:) = 0.0_dp
395 pint_env%rtmp_ndim = 0._dp
396 pint_env%rtmp_natom = 0._dp
397 pint_env%time_per_step = 0.0_dp
398
399 IF (pint_env%transform == transformation_stage) THEN
400 transform_section => section_vals_get_subs_vals(input, &
401 "MOTION%PINT%STAGING")
402 ALLOCATE (pint_env%staging_env)
403 CALL staging_env_create(pint_env%staging_env, transform_section, &
404 p=pint_env%p, kt=pint_env%kT)
405 ELSE
406 transform_section => section_vals_get_subs_vals(input, &
407 "MOTION%PINT%NORMALMODE")
408 IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN
409 CALL section_vals_val_get(transform_section, "GAMMA", explicit=explicit)
410 IF (explicit) THEN
411 CALL section_vals_val_set(transform_section, "GAMMA", &
412 r_val=sqrt(2.0_dp/(pint_env%p*pint_env%dt*pint_env%kT)))
413 ELSE
414 CALL section_vals_val_set(transform_section, "MODEFACTOR", &
415 r_val=0.5_dp*pint_env%p*pint_env%dt*pint_env%kT)
416 END IF
417 END IF
418 ALLOCATE (pint_env%normalmode_env)
419 CALL normalmode_env_create(pint_env%normalmode_env, &
420 transform_section, p=pint_env%p, kt=pint_env%kT, propagator=pint_env%propagator%prop_kind)
421 IF (para_env%is_source()) THEN
422 IF (pint_env%harm_integrator == integrate_numeric) THEN
423 IF (10.0_dp*pint_env%dt/real(pint_env%nrespa, dp) > &
424 twopi/(pint_env%p*sqrt(maxval(pint_env%normalmode_env%lambda))* &
425 pint_env%normalmode_env%modefactor)) THEN
426 msg = "PINT WARNING| Number of RESPA steps to small "// &
427 "to integrate the harmonic springs."
428 cpwarn(msg)
429 END IF
430 END IF
431 END IF
432 END IF
433 ALLOCATE (pint_env%mass(pint_env%ndim))
434 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
435 f_env=f_env)
436 CALL force_env_get(force_env=f_env%force_env, subsys=subsys)
437 CALL cp_subsys_get(subsys, particles=particles, gci=gci)
438
439!TODO length of pint_env%mass is redundant
440 idim = 0
441 DO iat = 1, pint_env%ndim/3
442 CALL get_atomic_kind(particles%els(iat)%atomic_kind, mass=mass)
443 DO idir = 1, 3
444 idim = idim + 1
445 pint_env%mass(idim) = mass
446 END DO
447 END DO
448 CALL f_env_rm_defaults(f_env, ierr)
449 cpassert(ierr == 0)
450
451 ALLOCATE (pint_env%Q(pint_env%p), &
452 pint_env%mass_beads(pint_env%p, pint_env%ndim), &
453 pint_env%mass_fict(pint_env%p, pint_env%ndim))
454 IF (pint_env%transform == transformation_stage) THEN
455 CALL staging_init_masses(pint_env%staging_env, mass=pint_env%mass, &
456 mass_beads=pint_env%mass_beads, mass_fict=pint_env%mass_fict, &
457 q=pint_env%Q)
458 ELSE
459 CALL normalmode_init_masses(pint_env%normalmode_env, &
460 mass=pint_env%mass, mass_beads=pint_env%mass_beads, &
461 mass_fict=pint_env%mass_fict, q=pint_env%Q)
462 END IF
463
464 NULLIFY (pint_env%gle)
465 gle_section => section_vals_get_subs_vals(input, "MOTION%PINT%GLE")
466 CALL section_vals_get(gle_section, explicit=explicit)
467 IF (explicit) THEN
468 ALLOCATE (pint_env%gle)
469 CALL gle_init(pint_env%gle, dt=pint_env%dt/pint_env%nrespa, temp=pint_env%kT, &
470 section=gle_section)
471 IF (pint_env%pimd_thermostat == thermostat_none .AND. pint_env%gle%ndim > 0) THEN
472 pint_env%pimd_thermostat = thermostat_gle
473
474 ! initialize a GLE with ALL degrees of freedom on node 0,
475 ! as it seems to me that here everything but force eval is replicated
476 pint_env%gle%loc_num_gle = pint_env%p*pint_env%ndim
477 pint_env%gle%glob_num_gle = pint_env%gle%loc_num_gle
478 ALLOCATE (pint_env%gle%map_info%index(pint_env%gle%loc_num_gle))
479 cpassert(stat == 0)
480 DO itmp = 1, pint_env%gle%loc_num_gle
481 pint_env%gle%map_info%index(itmp) = itmp
482 END DO
483 CALL gle_thermo_create(pint_env%gle, pint_env%gle%loc_num_gle)
484
485 ! here we should have read a_mat and c_mat;
486 !we can therefore compute the matrices needed for the propagator
487 ! deterministic part of the propagator
488 CALL gle_matrix_exp((-pint_env%dt/pint_env%nrespa*0.5_dp)*pint_env%gle%a_mat, &
489 pint_env%gle%ndim, 15, 15, pint_env%gle%gle_t)
490 ! stochastic part
491 CALL gle_cholesky_stab(pint_env%gle%c_mat - matmul(pint_env%gle%gle_t, &
492 matmul(pint_env%gle%c_mat, transpose(pint_env%gle%gle_t))), &
493 pint_env%gle%gle_s, pint_env%gle%ndim)
494 ! and initialize the additional momenta
495 CALL pint_gle_init(pint_env)
496 END IF
497 END IF
498
499 !Setup pile thermostat
500 NULLIFY (pint_env%pile_therm)
501 pile_section => section_vals_get_subs_vals(input, "MOTION%PINT%PILE")
502 CALL section_vals_get(pile_section, explicit=explicit)
503 IF (explicit) THEN
504 CALL cite_reference(ceriotti2010)
505 CALL section_vals_val_get(pint_env%input, &
506 "MOTION%PINT%INIT%THERMOSTAT_SEED", &
507 i_val=pint_env%thermostat_rng_seed)
508 IF (pint_env%pimd_thermostat == thermostat_none) THEN
509 pint_env%pimd_thermostat = thermostat_pile
510 ALLOCATE (pint_env%pile_therm)
511 CALL pint_pile_init(pile_therm=pint_env%pile_therm, &
512 pint_env=pint_env, &
513 normalmode_env=pint_env%normalmode_env, &
514 section=pile_section)
515 ELSE
516 cpabort("PILE thermostat can't be used with another thermostat.")
517 END IF
518 END IF
519
520 !Setup PIGLET thermostat
521 NULLIFY (pint_env%piglet_therm)
522 piglet_section => section_vals_get_subs_vals(input, "MOTION%PINT%PIGLET")
523 CALL section_vals_get(piglet_section, explicit=explicit)
524 IF (explicit) THEN
525 CALL cite_reference(ceriotti2012)
526 CALL section_vals_val_get(pint_env%input, &
527 "MOTION%PINT%INIT%THERMOSTAT_SEED", &
528 i_val=pint_env%thermostat_rng_seed)
529 IF (pint_env%pimd_thermostat == thermostat_none) THEN
530 pint_env%pimd_thermostat = thermostat_piglet
531 ALLOCATE (pint_env%piglet_therm)
532 CALL pint_piglet_create(pint_env%piglet_therm, &
533 pint_env, &
534 piglet_section)
535 CALL pint_piglet_init(pint_env%piglet_therm, &
536 pint_env, &
537 piglet_section, &
538 dt=pint_env%dt, para_env=para_env)
539 ELSE
540 cpabort("PIGLET thermostat can't be used with another thermostat.")
541 END IF
542 END IF
543
544 !Setup qtb thermostat
545 NULLIFY (pint_env%qtb_therm)
546 qtb_section => section_vals_get_subs_vals(input, "MOTION%PINT%QTB")
547 CALL section_vals_get(qtb_section, explicit=explicit)
548 IF (explicit) THEN
549 CALL cite_reference(brieuc2016)
550 CALL section_vals_val_get(pint_env%input, &
551 "MOTION%PINT%INIT%THERMOSTAT_SEED", &
552 i_val=pint_env%thermostat_rng_seed)
553 IF (pint_env%pimd_thermostat == thermostat_none) THEN
554 pint_env%pimd_thermostat = thermostat_qtb
555 CALL pint_qtb_init(qtb_therm=pint_env%qtb_therm, &
556 pint_env=pint_env, &
557 normalmode_env=pint_env%normalmode_env, &
558 section=qtb_section)
559 ELSE
560 cpabort("QTB thermostat can't be used with another thermostat.")
561 END IF
562 END IF
563
564 ! force PILE for BCMD
565 IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN
566 IF (pint_env%pimd_thermostat == thermostat_none) THEN
567 !crate section
568 pile_section => section_vals_get_subs_vals(input, "MOTION%PINT%PILE")
569 CALL section_vals_get(pile_section, explicit=explicit)
570 IF (.NOT. explicit) THEN
571 CALL section_vals_add_values(pile_section)
572 END IF
573 CALL section_vals_val_get(pint_env%input, &
574 "MOTION%PINT%INIT%THERMOSTAT_SEED", &
575 i_val=pint_env%thermostat_rng_seed)
576 pint_env%pimd_thermostat = thermostat_pile
577 ALLOCATE (pint_env%pile_therm)
578 CALL pint_pile_init(pile_therm=pint_env%pile_therm, &
579 pint_env=pint_env, &
580 normalmode_env=pint_env%normalmode_env, &
581 section=pile_section)
582 END IF
583 IF (pint_env%pimd_thermostat /= thermostat_pile) THEN
584 cpabort("PILE/no thermostat currently needed for BCMD")
585 END IF
586 CALL cite_reference(shiga2022)
587 END IF
588
589 !Initialize integrator scheme
590 CALL section_vals_val_get(pint_section, "HARM_INT", i_val=pint_env%harm_integrator)
591 IF (pint_env%harm_integrator == integrate_exact) THEN
592 IF (pint_env%pimd_thermostat == thermostat_nose) THEN
593 WRITE (unit=msg, fmt=*) "PINT WARNING| Nose Thermostat only available in "// &
594 "the numeric harmonic integrator. Switching to numeric harmonic integrator."
595 cpwarn(msg)
596 pint_env%harm_integrator = integrate_numeric
597 END IF
598 IF (pint_env%pimd_thermostat == thermostat_gle) THEN
599 WRITE (unit=msg, fmt=*) "PINT WARNING| GLE Thermostat only available in "// &
600 "the numeric harmonic integrator. Switching to numeric harmonic integrator."
601 cpwarn(msg)
602 pint_env%harm_integrator = integrate_numeric
603 END IF
604 ELSE IF (pint_env%harm_integrator == integrate_numeric) THEN
605 IF (pint_env%pimd_thermostat == thermostat_pile) THEN
606 WRITE (unit=msg, fmt=*) "PINT WARNING| PILE Thermostat only available in "// &
607 "the exact harmonic integrator. Switching to exact harmonic integrator."
608 cpwarn(msg)
609 pint_env%harm_integrator = integrate_exact
610 END IF
611 IF (pint_env%pimd_thermostat == thermostat_piglet) THEN
612 WRITE (unit=msg, fmt=*) "PINT WARNING| PIGLET Thermostat only available in "// &
613 "the exact harmonic integrator. Switching to exact harmonic integrator."
614 cpwarn(msg)
615 pint_env%harm_integrator = integrate_exact
616 END IF
617 IF (pint_env%pimd_thermostat == thermostat_qtb) THEN
618 WRITE (unit=msg, fmt=*) "PINT WARNING| QTB Thermostat only available in "// &
619 "the exact harmonic integrator. Switching to exact harmonic integrator."
620 cpwarn(msg)
621 pint_env%harm_integrator = integrate_exact
622 END IF
623 IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN
624 WRITE (unit=msg, fmt=*) "PINT WARNING| BCMD needs the exact harmonic "// &
625 "integrator. Switching to exact harmonic integrator."
626 cpwarn(msg)
627 pint_env%harm_integrator = integrate_exact
628 END IF
629 END IF
630
631 IF (pint_env%harm_integrator == integrate_exact) THEN
632 IF (pint_env%nrespa /= 1) THEN
633 pint_env%nrespa = 1
634 WRITE (unit=msg, fmt=*) "PINT WARNING| Adjusting NRESPA to 1 for exact harmonic integration."
635 cpwarn(msg)
636 END IF
637 NULLIFY (pint_env%wsinex)
638 ALLOCATE (pint_env%wsinex(pint_env%p))
639 NULLIFY (pint_env%iwsinex)
640 ALLOCATE (pint_env%iwsinex(pint_env%p))
641 NULLIFY (pint_env%cosex)
642 ALLOCATE (pint_env%cosex(pint_env%p))
643 dt = pint_env%dt/real(pint_env%nrespa, kind=dp)
644 !Centroid mode shoud not be propagated
645 pint_env%wsinex(1) = 0.0_dp
646 pint_env%iwsinex(1) = dt
647 pint_env%cosex(1) = 1.0_dp
648 DO ibead = 2, pint_env%p
649 omega = sqrt(pint_env%normalmode_env%lambda(ibead))
650 pint_env%wsinex(ibead) = sin(omega*dt)*omega
651 pint_env%iwsinex(ibead) = sin(omega*dt)/omega
652 pint_env%cosex(ibead) = cos(omega*dt)
653 END DO
654 END IF
655
656 CALL section_vals_val_get(pint_section, "FIX_CENTROID_POS", &
657 l_val=ltmp)
658 IF (ltmp .AND. (pint_env%transform == transformation_normal)) THEN
659 pint_env%first_propagated_mode = 2
660 ELSE
661 pint_env%first_propagated_mode = 1
662 END IF
663
664 ! Constraint information:
665 NULLIFY (pint_env%simpar)
666 constraint_section => section_vals_get_subs_vals(pint_env%input, &
667 "MOTION%CONSTRAINT")
668 CALL section_vals_get(constraint_section, explicit=explicit)
669 CALL create_simpar_type(pint_env%simpar)
670 pint_env%simpar%constraint = explicit
671 pint_env%kTcorr = 1.0_dp
672
673 ! Determine if beadwise constraints are activated
674 pint_env%beadwise_constraints = .false.
675 CALL section_vals_val_get(constraint_section, "PIMD_BEADWISE_CONSTRAINT", &
676 l_val=pint_env%beadwise_constraints)
677 IF (pint_env%simpar%constraint) THEN
678 IF (pint_env%beadwise_constraints) THEN
679 CALL pint_write_line("Using beadwise constraints")
680 ELSE
681 CALL pint_write_line("Using centroid constraints")
682 END IF
683 END IF
684
685 IF (explicit) THEN
686 ! Staging not supported yet, since lowest mode is assumed
687 ! to be related to centroid
688 IF (pint_env%transform == transformation_stage) THEN
689 cpabort("Constraints are not supported for staging transformation")
690 END IF
691
692 ! Check thermostats that are not supported:
693 IF (pint_env%pimd_thermostat == thermostat_gle) THEN
694 WRITE (unit=msg, fmt=*) "GLE Thermostat not supported for "// &
695 "constraints. Switch to NOSE for numeric integration."
696 cpabort(msg)
697 END IF
698 ! Warn for NOSE
699 IF (pint_env%pimd_thermostat == thermostat_nose) THEN
700 !Beadwise constraints not supported
701 IF (pint_env%beadwise_constraints) THEN
702 cpabort("Beadwise constraints are not supported for NOSE Thermostat.")
703 !Centroid constraints supported
704 ELSE
705 WRITE (unit=msg, fmt=*) "PINT WARNING| Nose Thermostat set to "// &
706 "zero for constrained atoms. Careful interpretation of temperature."
707 cpwarn(msg)
708 WRITE (unit=msg, fmt=*) "PINT WARNING| Lagrange multipliers are "// &
709 "are printed every RESPA step and need to be treated carefully."
710 cpwarn(msg)
711 END IF
712 END IF
713
714 CALL section_vals_val_get(constraint_section, "SHAKE_TOLERANCE", &
715 r_val=pint_env%simpar%shake_tol)
716 pint_env%simpar%info_constraint = cp_print_key_unit_nr(pint_env%logger, &
717 constraint_section, &
718 "CONSTRAINT_INFO", &
719 extension=".shakeLog", &
720 log_filename=.false.)
721 pint_env%simpar%lagrange_multipliers = cp_print_key_unit_nr(pint_env%logger, &
722 constraint_section, &
723 "LAGRANGE_MULTIPLIERS", &
724 extension=".LagrangeMultLog", &
725 log_filename=.false.)
726 pint_env%simpar%dump_lm = &
727 btest(cp_print_key_should_output(pint_env%logger%iter_info, &
728 constraint_section, &
729 "LAGRANGE_MULTIPLIERS"), cp_p_file)
730
731 ! Determine constrained atoms:
732 pint_env%n_atoms_constraints = 0
733 DO ig = 1, gci%ncolv%ntot
734 ! Double counts, if the same atom is involved in different collective variables
735 pint_env%n_atoms_constraints = pint_env%n_atoms_constraints + SIZE(gci%colv_list(ig)%i_atoms)
736 END DO
737
738 ALLOCATE (pint_env%atoms_constraints(pint_env%n_atoms_constraints))
739 icont = 0
740 DO ig = 1, gci%ncolv%ntot
741 DO iat = 1, SIZE(gci%colv_list(ig)%i_atoms)
742 icont = icont + 1
743 pint_env%atoms_constraints(icont) = gci%colv_list(ig)%i_atoms(iat)
744 END DO
745 END DO
746
747 ! Set the correction to the temperature due to the frozen degrees of freedom in NOSE:
748 CALL section_vals_val_get(pint_section, "kT_CORRECTION", &
749 l_val=ltmp)
750 IF (ltmp) THEN
751 pint_env%kTcorr = 1.0_dp + real(3*pint_env%n_atoms_constraints, dp)/(real(pint_env%ndim, dp)*real(pint_env%p, dp))
752 END IF
753 END IF
754
755 CALL timestop(handle)
756
757 END SUBROUTINE pint_create
758
759! ***************************************************************************
760!> \brief Release a path integral environment
761!> \param pint_env the pint_env to release
762!> \par History
763!> Added normal mode transformation [hforbert]
764!> \author Fawzi Mohamed
765! **************************************************************************************************
766 SUBROUTINE pint_release(pint_env)
767 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
768
769 CALL rep_env_release(pint_env%replicas)
770 CALL section_vals_release(pint_env%input)
771 IF (ASSOCIATED(pint_env%staging_env)) THEN
772 CALL staging_release(pint_env%staging_env)
773 DEALLOCATE (pint_env%staging_env)
774 END IF
775 IF (ASSOCIATED(pint_env%normalmode_env)) THEN
776 CALL normalmode_release(pint_env%normalmode_env)
777 DEALLOCATE (pint_env%normalmode_env)
778 END IF
779
780 DEALLOCATE (pint_env%mass)
781 DEALLOCATE (pint_env%e_pot_bead)
782
783 DEALLOCATE (pint_env%x)
784 DEALLOCATE (pint_env%v)
785 DEALLOCATE (pint_env%f)
786 DEALLOCATE (pint_env%external_f)
787 DEALLOCATE (pint_env%mass_beads)
788 DEALLOCATE (pint_env%mass_fict)
789 DEALLOCATE (pint_env%ux)
790 DEALLOCATE (pint_env%ux_t)
791 DEALLOCATE (pint_env%uv)
792 DEALLOCATE (pint_env%uv_t)
793 DEALLOCATE (pint_env%uv_new)
794 DEALLOCATE (pint_env%uf)
795 DEALLOCATE (pint_env%uf_h)
796 DEALLOCATE (pint_env%centroid)
797 DEALLOCATE (pint_env%rtmp_ndim)
798 DEALLOCATE (pint_env%rtmp_natom)
799 DEALLOCATE (pint_env%propagator)
800
801 IF (pint_env%simpar%constraint) THEN
802 DEALLOCATE (pint_env%atoms_constraints)
803 END IF
804 CALL release_simpar_type(pint_env%simpar)
805
806 IF (pint_env%harm_integrator == integrate_exact) THEN
807 DEALLOCATE (pint_env%wsinex)
808 DEALLOCATE (pint_env%iwsinex)
809 DEALLOCATE (pint_env%cosex)
810 END IF
811
812 SELECT CASE (pint_env%pimd_thermostat)
813 CASE (thermostat_nose)
814 DEALLOCATE (pint_env%tx)
815 DEALLOCATE (pint_env%tv)
816 DEALLOCATE (pint_env%tv_t)
817 DEALLOCATE (pint_env%tv_old)
818 DEALLOCATE (pint_env%tv_new)
819 DEALLOCATE (pint_env%tf)
820 CASE (thermostat_gle)
821 CALL gle_dealloc(pint_env%gle)
822 CASE (thermostat_pile)
823 CALL pint_pile_release(pint_env%pile_therm)
824 DEALLOCATE (pint_env%pile_therm)
825 CASE (thermostat_piglet)
826 CALL pint_piglet_release(pint_env%piglet_therm)
827 DEALLOCATE (pint_env%piglet_therm)
828 CASE (thermostat_qtb)
829 CALL pint_qtb_release(pint_env%qtb_therm)
830 DEALLOCATE (pint_env%qtb_therm)
831 END SELECT
832
833 DEALLOCATE (pint_env%Q)
834
835 END SUBROUTINE pint_release
836
837! ***************************************************************************
838!> \brief Tests the path integral methods
839!> \param para_env parallel environment
840!> \param input the input to test
841!> \param input_declaration ...
842!> \author fawzi
843! **************************************************************************************************
844 SUBROUTINE pint_test(para_env, input, input_declaration)
845 TYPE(mp_para_env_type), POINTER :: para_env
846 TYPE(section_vals_type), POINTER :: input
847 TYPE(section_type), POINTER :: input_declaration
848
849 INTEGER :: i, ib, idim, unit_nr
850 REAL(kind=dp) :: c, e_h, err
851 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: x1
852 TYPE(pint_env_type) :: pint_env
853
854 cpassert(ASSOCIATED(para_env))
855 cpassert(ASSOCIATED(input))
856 cpassert(para_env%is_valid())
857 cpassert(input%ref_count > 0)
858 unit_nr = cp_logger_get_default_io_unit()
859 CALL pint_create(pint_env, input, input_declaration, para_env)
860 ALLOCATE (x1(pint_env%ndim, pint_env%p))
861 x1(:, :) = pint_env%x
862 CALL pint_x2u(pint_env)
863 pint_env%x = 0._dp
864 CALL pint_u2x(pint_env)
865 err = 0._dp
866 DO i = 1, pint_env%ndim
867 err = max(err, abs(x1(1, i) - pint_env%x(1, i)))
868 END DO
869 IF (unit_nr > 0) WRITE (unit_nr, *) "diff_r1="//cp_to_string(err)
870
871 CALL pint_calc_uf_h(pint_env, e_h=e_h)
872 c = -pint_env%staging_env%w_p**2
873 pint_env%f = 0._dp
874 DO idim = 1, pint_env%ndim
875 DO ib = 1, pint_env%p
876 pint_env%f(ib, idim) = pint_env%f(ib, idim) + &
877 c*(2._dp*pint_env%x(ib, idim) &
878 - pint_env%x(modulo(ib - 2, pint_env%p) + 1, idim) &
879 - pint_env%x(modulo(ib, pint_env%p) + 1, idim))
880 END DO
881 END DO
882 CALL pint_f2uf(pint_env)
883 err = 0._dp
884 DO idim = 1, pint_env%ndim
885 DO ib = 1, pint_env%p
886 err = max(err, abs(pint_env%uf(ib, idim) - pint_env%uf_h(ib, idim)))
887 END DO
888 END DO
889 IF (unit_nr > 0) WRITE (unit_nr, *) "diff_f_h="//cp_to_string(err)
890
891 END SUBROUTINE pint_test
892
893! ***************************************************************************
894!> \brief Perform a path integral simulation
895!> \param para_env parallel environment
896!> \param input the input to test
897!> \param input_declaration ...
898!> \param globenv ...
899!> \par History
900!> 2003-11 created [fawzi]
901!> 2009-12-14 globenv parameter added to handle soft exit
902!> requests [lwalewski]
903!> 2016-07-14 Modified to work with independent helium_env [cschran]
904!> \author Fawzi Mohamed
905! **************************************************************************************************
906 SUBROUTINE do_pint_run(para_env, input, input_declaration, globenv)
907 TYPE(mp_para_env_type), POINTER :: para_env
908 TYPE(section_vals_type), POINTER :: input
909 TYPE(section_type), POINTER :: input_declaration
910 TYPE(global_environment_type), POINTER :: globenv
911
912 CHARACTER(len=*), PARAMETER :: routinen = 'do_pint_run'
913 INTEGER, PARAMETER :: helium_only_mid = 1, &
914 int_pot_scan_mid = 4, &
915 solute_only_mid = 2, &
916 solute_with_helium_mid = 3
917
918 CHARACTER(len=default_string_length) :: stmp
919 INTEGER :: handle, mode
920 LOGICAL :: explicit, helium_only, int_pot_scan, &
921 solvent_present
922 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
923 TYPE(pint_env_type) :: pint_env
924 TYPE(section_vals_type), POINTER :: helium_section
925
926 CALL timeset(routinen, handle)
927
928 cpassert(ASSOCIATED(para_env))
929 cpassert(ASSOCIATED(input))
930 cpassert(para_env%is_valid())
931 cpassert(input%ref_count > 0)
932
933 ! check if helium solvent is present
934 NULLIFY (helium_section)
935 helium_section => section_vals_get_subs_vals(input, &
936 "MOTION%PINT%HELIUM")
937 CALL section_vals_get(helium_section, explicit=explicit)
938 IF (explicit) THEN
939 CALL section_vals_val_get(helium_section, "_SECTION_PARAMETERS_", &
940 l_val=solvent_present)
941 ELSE
942 solvent_present = .false.
943 END IF
944
945 ! check if there is anything but helium
946 IF (solvent_present) THEN
947 CALL section_vals_val_get(helium_section, "HELIUM_ONLY", &
948 l_val=helium_only)
949 ELSE
950 helium_only = .false.
951 END IF
952
953 ! check wheather to perform solute-helium interaction pot scan
954 IF (solvent_present) THEN
955 CALL section_vals_val_get(helium_section, "INTERACTION_POT_SCAN", &
956 l_val=int_pot_scan)
957 ELSE
958 int_pot_scan = .false.
959 END IF
960
961 ! input consistency check
962 IF (helium_only .AND. int_pot_scan) THEN
963 stmp = "Options HELIUM_ONLY and INTERACTION_POT_SCAN are exclusive"
964 cpabort(stmp)
965 END IF
966
967 ! select mode of operation
968 mode = 0
969 IF (solvent_present) THEN
970 IF (helium_only) THEN
971 mode = helium_only_mid
972 ELSE
973 IF (int_pot_scan) THEN
974 mode = int_pot_scan_mid
975 ELSE
976 mode = solute_with_helium_mid
977 END IF
978 END IF
979 ELSE
980 mode = solute_only_mid
981 END IF
982
983 ! perform the simulation according to the chosen mode
984 SELECT CASE (mode)
985
986 CASE (helium_only_mid)
987 CALL helium_create(helium_env, input)
988 CALL helium_init(helium_env, pint_env)
989 CALL helium_do_run(helium_env, globenv)
990 CALL helium_release(helium_env)
991
992 CASE (solute_only_mid)
993 CALL pint_create(pint_env, input, input_declaration, para_env)
994 CALL pint_init(pint_env)
995 CALL pint_do_run(pint_env, globenv)
996 CALL pint_release(pint_env)
997
998 CASE (int_pot_scan_mid)
999 CALL pint_create(pint_env, input, input_declaration, para_env)
1000! TODO only initialization of positions is necessary, but rep_env_calc_e_f called
1001! from within pint_init_f does something to the replica environments which can not be
1002! avoided (has something to do with f_env_add_defaults) so leaving for now..
1003 CALL pint_init(pint_env)
1004 CALL helium_create(helium_env, input, solute=pint_env)
1005 CALL pint_run_scan(pint_env, helium_env)
1006 CALL helium_release(helium_env)
1007 CALL pint_release(pint_env)
1008
1009 CASE (solute_with_helium_mid)
1010 CALL pint_create(pint_env, input, input_declaration, para_env)
1011 ! init pint without helium forces (they are not yet initialized)
1012 CALL pint_init(pint_env)
1013 ! init helium with solute's positions (they are already initialized)
1014 CALL helium_create(helium_env, input, solute=pint_env)
1015 CALL helium_init(helium_env, pint_env)
1016 ! reinit pint forces with helium forces (they are now initialized)
1017 CALL pint_init_f(pint_env, helium_env=helium_env)
1018
1019 CALL pint_do_run(pint_env, globenv, helium_env=helium_env)
1020 CALL helium_release(helium_env)
1021 CALL pint_release(pint_env)
1022
1023 CASE DEFAULT
1024 cpabort("Unknown mode ("//trim(adjustl(cp_to_string(mode)))//")")
1025 END SELECT
1026
1027 CALL timestop(handle)
1028
1029 END SUBROUTINE do_pint_run
1030
1031! ***************************************************************************
1032!> \brief Reads the restart, initializes the beads, etc.
1033!> \param pint_env ...
1034!> \par History
1035!> 11.2003 created [fawzi]
1036!> actually ASSIGN input pointer [hforbert]
1037!> 2010-12-16 turned into a wrapper routine [lwalewski]
1038!> \author Fawzi Mohamed
1039! **************************************************************************************************
1040 SUBROUTINE pint_init(pint_env)
1041
1042 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1043
1044 CALL pint_init_x(pint_env)
1045 CALL pint_init_v(pint_env)
1046 CALL pint_init_t(pint_env)
1047 CALL pint_init_f(pint_env)
1048
1049 END SUBROUTINE pint_init
1050
1051! ***************************************************************************
1052!> \brief Assign initial postions to the beads.
1053!> \param pint_env ...
1054!> \date 2010-12-15
1055!> \author Lukasz Walewski
1056!> \note Initialization is done in the following way:
1057!> 1. assign all beads with the same classical positions from
1058!> FORCE_EVAL (hot start)
1059!> 2. spread the beads around classical positions as if they were
1060!> free particles (if requested)
1061!> 3. replace positions generated in steps 1-2 with the explicit
1062!> ones if they are explicitly given in the input structure
1063!> 4. apply Gaussian noise to the positions generated so far (if
1064!> requested)
1065! **************************************************************************************************
1066 SUBROUTINE pint_init_x(pint_env)
1067
1068 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1069
1070 CHARACTER(len=5*default_string_length) :: msg, tmp
1071 INTEGER :: ia, ib, ic, idim, input_seed, n_rep_val
1072 LOGICAL :: done_init, done_levy, done_rand, &
1073 explicit, levycorr, ltmp
1074 REAL(kind=dp) :: tcorr, var
1075 REAL(kind=dp), DIMENSION(3) :: x0
1076 REAL(kind=dp), DIMENSION(3, 2) :: seed
1077 REAL(kind=dp), DIMENSION(:), POINTER :: bx, r_vals
1078 TYPE(rng_stream_type), ALLOCATABLE :: rng_gaussian
1079 TYPE(section_vals_type), POINTER :: input_section
1080
1081 DO idim = 1, pint_env%ndim
1082 DO ib = 1, pint_env%p
1083 pint_env%x(ib, idim) = pint_env%replicas%r(idim, ib)
1084 END DO
1085 END DO
1086
1087 done_levy = .false.
1088 CALL section_vals_val_get(pint_env%input, &
1089 "MOTION%PINT%INIT%LEVY_POS_SAMPLE", &
1090 l_val=ltmp)
1091 CALL section_vals_val_get(pint_env%input, &
1092 "MOTION%PINT%INIT%LEVY_TEMP_FACTOR", &
1093 r_val=tcorr)
1094 IF (ltmp) THEN
1095
1096 IF (pint_env%beadwise_constraints) THEN
1097 WRITE (unit=msg, fmt=*) "Beadwise constraints are not supported for "// &
1098 "the initialization of the beads as free particles. "// &
1099 "Please use hot start (default)."
1100 cpabort(msg)
1101 END IF
1102
1103 NULLIFY (bx)
1104 ALLOCATE (bx(3*pint_env%p))
1105 CALL section_vals_val_get(pint_env%input, &
1106 "MOTION%PINT%INIT%LEVY_SEED", i_val=input_seed)
1107 seed(:, :) = real(input_seed, kind=dp)
1108! seed(:,:) = next_rng_seed()
1109 rng_gaussian = rng_stream_type( &
1110 name="tmp_rng_gaussian", &
1111 distribution_type=gaussian, &
1112 extended_precision=.true., &
1113 seed=seed)
1114
1115 CALL section_vals_val_get(pint_env%input, &
1116 "MOTION%PINT%INIT%LEVY_CORRELATED", &
1117 l_val=levycorr)
1118
1119 IF (levycorr) THEN
1120
1121 ! correlated Levy walk - the same path for all atoms
1122 x0 = [0.0_dp, 0.0_dp, 0.0_dp]
1123 CALL pint_levy_walk(x0, pint_env%p, 1.0_dp, bx, rng_gaussian)
1124 idim = 0
1125 DO ia = 1, pint_env%ndim/3
1126 var = sqrt(1.0_dp/(pint_env%kT*tcorr*pint_env%mass(3*ia)))
1127 DO ic = 1, 3
1128 idim = idim + 1
1129 DO ib = 1, pint_env%p
1130 pint_env%x(ib, idim) = pint_env%x(ib, idim) + bx(3*(ib - 1) + ic)*var
1131 END DO
1132 END DO
1133 END DO
1134
1135 ELSE
1136
1137 ! uncorrelated bead initialization - distinct Levy walk for each atom
1138 idim = 0
1139 DO ia = 1, pint_env%ndim/3
1140 x0(1) = pint_env%x(1, 3*(ia - 1) + 1)
1141 x0(2) = pint_env%x(1, 3*(ia - 1) + 2)
1142 x0(3) = pint_env%x(1, 3*(ia - 1) + 3)
1143 var = sqrt(1.0_dp/(pint_env%kT*tcorr*pint_env%mass(3*ia)))
1144 CALL pint_levy_walk(x0, pint_env%p, var, bx, rng_gaussian)
1145 DO ic = 1, 3
1146 idim = idim + 1
1147 DO ib = 1, pint_env%p
1148 pint_env%x(ib, idim) = pint_env%x(ib, idim) + bx(3*(ib - 1) + ic)
1149 END DO
1150 END DO
1151 END DO
1152
1153 END IF
1154
1155 DEALLOCATE (bx)
1156 done_levy = .true.
1157 END IF
1158
1159 done_init = .false.
1160 NULLIFY (input_section)
1161 input_section => section_vals_get_subs_vals(pint_env%input, &
1162 "MOTION%PINT%BEADS%COORD")
1163 CALL section_vals_get(input_section, explicit=explicit)
1164 IF (explicit) THEN
1165 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1166 n_rep_val=n_rep_val)
1167 IF (n_rep_val > 0) THEN
1168 cpassert(n_rep_val == 1)
1169 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1170 r_vals=r_vals)
1171 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim) &
1172 cpabort("Invalid size of MOTION%PINT%BEADS%COORD")
1173 ic = 0
1174 DO idim = 1, pint_env%ndim
1175 DO ib = 1, pint_env%p
1176 ic = ic + 1
1177 pint_env%x(ib, idim) = r_vals(ic)
1178 END DO
1179 END DO
1180 done_init = .true.
1181 END IF
1182 END IF
1183
1184 done_rand = .false.
1185 CALL section_vals_val_get(pint_env%input, &
1186 "MOTION%PINT%INIT%RANDOMIZE_POS", &
1187 l_val=ltmp)
1188 IF (ltmp) THEN
1189
1190 IF (pint_env%beadwise_constraints) THEN
1191 WRITE (unit=msg, fmt=*) "Beadwise constraints are not supported if "// &
1192 "a random noise is applied to the initialization of the bead positions. "// &
1193 "Please use hot start (default)."
1194 cpabort(msg)
1195 END IF
1196
1197 DO idim = 1, pint_env%ndim
1198 DO ib = 1, pint_env%p
1199 pint_env%x(ib, idim) = pint_env%x(ib, idim) + &
1200 pint_env%randomG%next(variance=pint_env%beta/ &
1201 sqrt(12.0_dp*pint_env%mass(idim)))
1202 END DO
1203 END DO
1204 done_rand = .true.
1205 END IF
1206
1207 WRITE (tmp, '(A)') "Bead positions initialization:"
1208 IF (done_init) THEN
1209 WRITE (msg, '(A,A)') trim(tmp), " input structure"
1210 ELSE IF (done_levy) THEN
1211 WRITE (msg, '(A,A)') trim(tmp), " Levy random walk"
1212 ELSE
1213 WRITE (msg, '(A,A)') trim(tmp), " hot start"
1214 END IF
1215 CALL pint_write_line(msg)
1216
1217 IF (done_levy) THEN
1218 WRITE (msg, '(A,F6.3)') "Levy walk at effective temperature: ", tcorr
1219 END IF
1220
1221 IF (done_rand) THEN
1222 WRITE (msg, '(A)') "Added gaussian noise to the positions of the beads."
1223 CALL pint_write_line(msg)
1224 END IF
1225
1226 END SUBROUTINE pint_init_x
1227
1228! ***************************************************************************
1229!> \brief Initialize velocities
1230!> \param pint_env the pint env in which you should initialize the
1231!> velocity
1232!> \par History
1233!> 2010-12-16 gathered all velocity-init code here [lwalewski]
1234!> 2011-04-05 added centroid velocity initialization [lwalewski]
1235!> 2011-12-19 removed optional parameter kT, target temperature is
1236!> now determined from the input directly [lwalewski]
1237!> \author fawzi
1238!> \note Initialization is done according to the following protocol:
1239!> 1. set all the velocities to FORCE_EVAL%SUBSYS%VELOCITY if present
1240!> 2. scale the velocities according to the actual temperature
1241!> (has no effect if vels not present in 1.)
1242!> 3. draw vels for the remaining dof from MB distribution
1243!> (all or non-centroid modes only depending on 1.)
1244!> 4. add random noise to the centroid vels if CENTROID_SPEED == T
1245!> 5. set the vels for all dof to 0.0 if VELOCITY_QUENCH == T
1246!> 6. set the vels according to the explicit values from the input
1247!> if present
1248! **************************************************************************************************
1249 SUBROUTINE pint_init_v(pint_env)
1250 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1251
1252 CHARACTER(len=default_string_length) :: msg, stmp, stmp1, stmp2, unit_str
1253 INTEGER :: first_mode, i, ia, ib, ic, idim, ierr, &
1254 itmp, j, n_rep_val, nparticle, &
1255 nparticle_kind
1256 LOGICAL :: done_init, done_quench, done_scale, &
1257 done_sped, explicit, ltmp, vels_present
1258 REAL(kind=dp) :: actual_t, ek, factor, rtmp, target_t, &
1259 unit_conv
1260 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vel
1261 REAL(kind=dp), DIMENSION(:), POINTER :: r_vals
1262 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1263 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1264 TYPE(cell_type), POINTER :: cell
1265 TYPE(cp_logger_type), POINTER :: logger
1266 TYPE(cp_subsys_type), POINTER :: subsys
1267 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
1268 TYPE(f_env_type), POINTER :: f_env
1269 TYPE(global_constraint_type), POINTER :: gci
1270 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1271 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1272 TYPE(molecule_list_type), POINTER :: molecules
1273 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1274 TYPE(particle_list_type), POINTER :: particles
1275 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1276 TYPE(section_vals_type), POINTER :: input_section
1277
1278 NULLIFY (logger)
1279 logger => cp_get_default_logger()
1280
1281 ! Get constraint info, if needed
1282 ! Create a force environment which will be identical to
1283 ! the bead that is being processed by the processor.
1284 IF (pint_env%simpar%constraint) THEN
1285 NULLIFY (subsys, cell)
1286 NULLIFY (atomic_kinds, local_particles, particles)
1287 NULLIFY (local_molecules, molecules, molecule_kinds, gci)
1288 NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set)
1289
1290 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env)
1291 CALL force_env_get(force_env=f_env%force_env, subsys=subsys)
1292 CALL f_env_rm_defaults(f_env, ierr)
1293 cpassert(ierr == 0)
1294
1295 ! Get gci and more from subsys
1296 CALL cp_subsys_get(subsys=subsys, &
1297 cell=cell, &
1298 atomic_kinds=atomic_kinds, &
1299 local_particles=local_particles, &
1300 particles=particles, &
1301 local_molecules=local_molecules, &
1302 molecules=molecules, &
1303 molecule_kinds=molecule_kinds, &
1304 gci=gci)
1305
1306 nparticle_kind = atomic_kinds%n_els
1307 atomic_kind_set => atomic_kinds%els
1308 molecule_kind_set => molecule_kinds%els
1309 nparticle = particles%n_els
1310 particle_set => particles%els
1311 molecule_set => molecules%els
1312
1313 ! Allocate work storage
1314 ALLOCATE (vel(3, nparticle))
1315 vel(:, :) = 0.0_dp
1316 CALL getold(gci, local_molecules, molecule_set, &
1317 molecule_kind_set, particle_set, cell)
1318 END IF
1319
1320 ! read the velocities from the input file if they are given explicitly
1321 vels_present = .false.
1322 NULLIFY (input_section)
1323 input_section => section_vals_get_subs_vals(pint_env%input, &
1324 "FORCE_EVAL%SUBSYS%VELOCITY")
1325 CALL section_vals_get(input_section, explicit=explicit)
1326 IF (explicit) THEN
1327
1328 CALL section_vals_val_get(input_section, "PINT_UNIT", &
1329 c_val=unit_str)
1330 unit_conv = cp_unit_to_cp2k(1.0_dp, trim(unit_str))
1331
1332 ! assign all the beads with the same velocities from FORCE_EVAL%SUBSYS%VELOCITY
1333 NULLIFY (r_vals)
1334 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1335 n_rep_val=n_rep_val)
1336 stmp = ""
1337 WRITE (stmp, *) n_rep_val
1338 msg = "Invalid number of atoms in FORCE_EVAL%SUBSYS%VELOCITY ("// &
1339 trim(adjustl(stmp))//")."
1340 IF (3*n_rep_val /= pint_env%ndim) &
1341 cpabort(msg)
1342 DO ia = 1, pint_env%ndim/3
1343 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1344 i_rep_val=ia, r_vals=r_vals)
1345 itmp = SIZE(r_vals)
1346 stmp = ""
1347 WRITE (stmp, *) itmp
1348 msg = "Number of coordinates != 3 in FORCE_EVAL%SUBSYS%VELOCITY ("// &
1349 trim(adjustl(stmp))//")."
1350 IF (itmp /= 3) &
1351 cpabort(msg)
1352 DO ib = 1, pint_env%p
1353 DO ic = 1, 3
1354 idim = 3*(ia - 1) + ic
1355 pint_env%v(ib, idim) = r_vals(ic)*unit_conv
1356 END DO
1357 END DO
1358 END DO
1359
1360 vels_present = .true.
1361 END IF
1362
1363 ! set the actual temperature...
1364 IF (vels_present) THEN
1365 ! ...from the initial velocities
1366 ek = 0.0_dp
1367 DO ia = 1, pint_env%ndim/3
1368 rtmp = 0.0_dp
1369 DO ic = 1, 3
1370 idim = 3*(ia - 1) + ic
1371 rtmp = rtmp + pint_env%v(1, idim)*pint_env%v(1, idim)
1372 END DO
1373 ek = ek + 0.5_dp*pint_env%mass(idim)*rtmp
1374 END DO
1375 actual_t = 2.0_dp*ek/pint_env%ndim
1376 ELSE
1377 ! ...using the temperature value from the input
1378 actual_t = pint_env%kT
1379 END IF
1380
1381 ! set the target temperature
1382 target_t = pint_env%kT
1383 CALL section_vals_val_get(pint_env%input, &
1384 "MOTION%PINT%INIT%VELOCITY_SCALE", &
1385 l_val=done_scale)
1386 IF (vels_present) THEN
1387 IF (done_scale) THEN
1388 ! rescale the velocities to match the target temperature
1389 rtmp = sqrt(target_t/actual_t)
1390 DO ia = 1, pint_env%ndim/3
1391 DO ib = 1, pint_env%p
1392 DO ic = 1, 3
1393 idim = 3*(ia - 1) + ic
1394 pint_env%v(ib, idim) = rtmp*pint_env%v(ib, idim)
1395 END DO
1396 END DO
1397 END DO
1398 ELSE
1399 target_t = actual_t
1400 END IF
1401 END IF
1402
1403 ! draw velocities from the M-B distribution...
1404 IF (vels_present) THEN
1405 ! ...for non-centroid modes only
1406 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
1407 first_mode = 2
1408 ELSE
1409 ! ...for all the modes
1410 first_mode = 1
1411 END IF
1412 DO idim = 1, SIZE(pint_env%uv, 2)
1413 DO ib = first_mode, SIZE(pint_env%uv, 1)
1414 pint_env%uv(ib, idim) = &
1415 pint_env%randomG%next(variance=target_t/pint_env%mass_fict(ib, idim))
1416 END DO
1417 END DO
1418
1419 ! add random component to the centroid velocity if requested
1420 done_sped = .false.
1421 CALL section_vals_val_get(pint_env%input, &
1422 "MOTION%PINT%INIT%CENTROID_SPEED", &
1423 l_val=ltmp)
1424 IF (ltmp) THEN
1425 CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
1426 DO idim = 1, pint_env%ndim
1427 rtmp = pint_env%randomG%next(variance=pint_env%mass(idim)*pint_env%kT) &
1428 /pint_env%mass(idim)
1429 DO ib = 1, pint_env%p
1430 pint_env%v(ib, idim) = pint_env%v(ib, idim) + rtmp
1431 END DO
1432 END DO
1433 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
1434 done_sped = .true.
1435 END IF
1436
1437 ! quench (set to zero) velocities for all the modes if requested
1438 ! (disregard all the initialization done so far)
1439 done_quench = .false.
1440 CALL section_vals_val_get(pint_env%input, &
1441 "MOTION%PINT%INIT%VELOCITY_QUENCH", &
1442 l_val=ltmp)
1443 IF (ltmp) THEN
1444 DO idim = 1, pint_env%ndim
1445 DO ib = 1, pint_env%p
1446 pint_env%v(ib, idim) = 0.0_dp
1447 END DO
1448 END DO
1449 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
1450 done_quench = .true.
1451 END IF
1452
1453 ! set the velocities to the values from the input if they are explicit
1454 ! (disregard all the initialization done so far)
1455 done_init = .false.
1456 NULLIFY (input_section)
1457 input_section => section_vals_get_subs_vals(pint_env%input, &
1458 "MOTION%PINT%BEADS%VELOCITY")
1459 CALL section_vals_get(input_section, explicit=explicit)
1460 IF (explicit) THEN
1461 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1462 n_rep_val=n_rep_val)
1463 IF (n_rep_val > 0) THEN
1464 cpassert(n_rep_val == 1)
1465 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1466 r_vals=r_vals)
1467 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim) &
1468 cpabort("Invalid size of MOTION%PINT%BEAD%VELOCITY")
1469 itmp = 0
1470 DO idim = 1, pint_env%ndim
1471 DO ib = 1, pint_env%p
1472 itmp = itmp + 1
1473 pint_env%v(ib, idim) = r_vals(itmp)
1474 END DO
1475 END DO
1476 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
1477 done_init = .true.
1478 END IF
1479 END IF
1480
1481 unit_conv = cp_unit_from_cp2k(1.0_dp, "K")
1482 WRITE (stmp1, '(F10.2)') target_t*pint_env%propagator%temp_sim2phys*unit_conv
1483 msg = "Bead velocities initialization:"
1484 IF (done_init) THEN
1485 msg = trim(msg)//" input structure"
1486 ELSE IF (done_quench) THEN
1487 msg = trim(msg)//" quenching (set to 0.0)"
1488 ELSE
1489 IF (vels_present) THEN
1490 msg = trim(adjustl(msg))//" centroid +"
1491 END IF
1492 msg = trim(adjustl(msg))//" Maxwell-Boltzmann at "//trim(adjustl(stmp1))//" K."
1493 END IF
1494 CALL pint_write_line(msg)
1495
1496 IF (done_init .AND. done_quench) THEN
1497 msg = "WARNING: exclusive options requested (velocity restart and quenching)"
1498 cpwarn(msg)
1499 msg = "WARNING: velocity restart took precedence"
1500 cpwarn(msg)
1501 END IF
1502
1503 IF ((.NOT. done_init) .AND. (.NOT. done_quench)) THEN
1504 IF (vels_present .AND. done_scale) THEN
1505 WRITE (stmp1, '(F10.2)') actual_t*unit_conv
1506 WRITE (stmp2, '(F10.2)') target_t*unit_conv
1507 msg = "Scaled initial velocities from "//trim(adjustl(stmp1))// &
1508 " to "//trim(adjustl(stmp2))//" K as requested."
1509 cpwarn(msg)
1510 END IF
1511 IF (done_sped) THEN
1512 msg = "Added random component to the initial centroid velocities."
1513 cpwarn(msg)
1514 END IF
1515 END IF
1516
1517 ! Apply constraints to the initial velocities
1518 IF (pint_env%simpar%constraint) THEN
1519 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
1520 ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD
1521 factor = sqrt(real(pint_env%p, dp))
1522 ELSE
1523 ! lowest NM is centroid
1524 factor = 1.0_dp
1525 END IF
1526 ! Beadwise constraints
1527 IF (pint_env%beadwise_constraints) THEN
1528 IF (pint_env%logger%para_env%is_source()) THEN
1529 CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
1530 DO ib = 1, pint_env%p
1531 DO i = 1, nparticle
1532 DO j = 1, 3
1533 ! Centroid is also constrained. This has to be changed if the initialization
1534 ! of the positions of the beads is done as free particles (LEVY_POS_SAMPLE)
1535 ! or if a Gaussian noise is added (RANDOMIZE_POS)
1536 particle_set(i)%r(j) = pint_env%x(1, j + (i - 1)*3)/factor
1537 vel(j, i) = pint_env%v(ib, j + (i - 1)*3)
1538 END DO
1539 END DO
1540 ! Possibly update the target values
1541 CALL shake_update_targets(gci, local_molecules, molecule_set, &
1542 molecule_kind_set, pint_env%dt, &
1543 f_env%force_env%root_section)
1544 CALL rattle_control(gci, local_molecules, molecule_set, &
1545 molecule_kind_set, particle_set, &
1546 vel, pint_env%dt, pint_env%simpar%shake_tol, &
1547 pint_env%simpar%info_constraint, &
1548 pint_env%simpar%lagrange_multipliers, &
1549 .false., &
1550 cell, mp_comm_self, &
1551 local_particles)
1552 DO i = 1, nparticle
1553 DO j = 1, 3
1554 pint_env%v(ib, j + (i - 1)*3) = vel(j, i)
1555 END DO
1556 END DO
1557 END DO
1558 ! Transform back to normal modes:
1559 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
1560 END IF
1561 ! Broadcast updated velocities to other nodes
1562 CALL pint_env%logger%para_env%bcast(pint_env%uv)
1563 ! Centroid constraints
1564 ELSE
1565 ! Transform positions and velocities to Cartesian coordinates:
1566 IF (pint_env%logger%para_env%is_source()) THEN
1567 DO i = 1, nparticle
1568 DO j = 1, 3
1569 particle_set(i)%r(j) = pint_env%x(1, j + (i - 1)*3)/factor
1570 vel(j, i) = pint_env%uv(1, j + (i - 1)*3)/factor
1571 END DO
1572 END DO
1573 ! Possibly update the target values
1574 CALL shake_update_targets(gci, local_molecules, molecule_set, &
1575 molecule_kind_set, pint_env%dt, &
1576 f_env%force_env%root_section)
1577 CALL rattle_control(gci, local_molecules, molecule_set, &
1578 molecule_kind_set, particle_set, &
1579 vel, pint_env%dt, pint_env%simpar%shake_tol, &
1580 pint_env%simpar%info_constraint, &
1581 pint_env%simpar%lagrange_multipliers, &
1582 .false., &
1583 cell, mp_comm_self, &
1584 local_particles)
1585 END IF
1586 ! Broadcast updated velocities to other nodes
1587 CALL pint_env%logger%para_env%bcast(vel)
1588 ! Transform back to normal modes
1589 DO i = 1, nparticle
1590 DO j = 1, 3
1591 pint_env%uv(1, j + (i - 1)*3) = vel(j, i)*factor
1592 END DO
1593 END DO
1594 END IF
1595 END IF
1596
1597 END SUBROUTINE pint_init_v
1598
1599! ***************************************************************************
1600!> \brief Assign initial postions and velocities to the thermostats.
1601!> \param pint_env ...
1602!> \param kT ...
1603!> \date 2010-12-15
1604!> \author Lukasz Walewski
1605!> \note Extracted from pint_init
1606! **************************************************************************************************
1607 SUBROUTINE pint_init_t(pint_env, kT)
1608
1609 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1610 REAL(kind=dp), INTENT(in), OPTIONAL :: kt
1611
1612 INTEGER :: ib, idim, ii, inos, n_rep_val
1613 LOGICAL :: explicit, gle_restart
1614 REAL(kind=dp) :: mykt
1615 REAL(kind=dp), DIMENSION(:), POINTER :: r_vals
1616 TYPE(section_vals_type), POINTER :: input_section
1617
1618 IF (pint_env%pimd_thermostat == thermostat_nose) THEN
1619
1620 mykt = pint_env%kT
1621 IF (PRESENT(kt)) mykt = kt
1622 DO idim = 1, SIZE(pint_env%tv, 3)
1623 DO ib = 1, SIZE(pint_env%tv, 2)
1624 DO inos = 1, SIZE(pint_env%tv, 1)
1625 pint_env%tv(inos, ib, idim) = &
1626 pint_env%randomG%next(variance=mykt/pint_env%Q(ib))
1627 END DO
1628 END DO
1629 END DO
1630 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
1631 pint_env%tv(:, 1, :) = 0.0_dp
1632 END IF
1633
1634 NULLIFY (input_section)
1635 input_section => section_vals_get_subs_vals(pint_env%input, &
1636 "MOTION%PINT%NOSE%COORD")
1637 CALL section_vals_get(input_section, explicit=explicit)
1638 IF (explicit) THEN
1639 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1640 n_rep_val=n_rep_val)
1641 IF (n_rep_val > 0) THEN
1642 cpassert(n_rep_val == 1)
1643 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1644 r_vals=r_vals)
1645 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim*pint_env%nnos) &
1646 cpabort("Invalid size of MOTION%PINT%NOSE%COORD")
1647 ii = 0
1648 DO idim = 1, pint_env%ndim
1649 DO ib = 1, pint_env%p
1650 DO inos = 1, pint_env%nnos
1651 ii = ii + 1
1652 pint_env%tx(inos, ib, idim) = r_vals(ii)
1653 END DO
1654 END DO
1655 END DO
1656 END IF
1657 END IF
1658 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
1659 pint_env%tx(:, 1, :) = 0.0_dp
1660 END IF
1661
1662 NULLIFY (input_section)
1663 input_section => section_vals_get_subs_vals(pint_env%input, &
1664 "MOTION%PINT%NOSE%VELOCITY")
1665 CALL section_vals_get(input_section, explicit=explicit)
1666 IF (explicit) THEN
1667 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1668 n_rep_val=n_rep_val)
1669 IF (n_rep_val > 0) THEN
1670 cpassert(n_rep_val == 1)
1671 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", &
1672 r_vals=r_vals)
1673 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim*pint_env%nnos) &
1674 cpabort("Invalid size of MOTION%PINT%NOSE%VELOCITY")
1675 ii = 0
1676 DO idim = 1, pint_env%ndim
1677 DO ib = 1, pint_env%p
1678 DO inos = 1, pint_env%nnos
1679 ii = ii + 1
1680 pint_env%tv(inos, ib, idim) = r_vals(ii)
1681 END DO
1682 END DO
1683 END DO
1684 END IF
1685 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
1686 pint_env%tv(:, 1, :) = 0.0_dp
1687 END IF
1688 END IF
1689
1690 ELSEIF (pint_env%pimd_thermostat == thermostat_gle) THEN
1691 NULLIFY (input_section)
1692 input_section => section_vals_get_subs_vals(pint_env%input, &
1693 "MOTION%PINT%GLE")
1694 CALL section_vals_get(input_section, explicit=explicit)
1695 IF (explicit) THEN
1696 CALL restart_gle(pint_env%gle, input_section, save_mem=.false., &
1697 restart=gle_restart)
1698 END IF
1699 END IF
1700
1701 END SUBROUTINE pint_init_t
1702
1703! ***************************************************************************
1704!> \brief Prepares the forces, etc. to perform an PIMD step
1705!> \param pint_env ...
1706!> \param helium_env ...
1707!> \par History
1708!> Added nh_energy calculation [hforbert]
1709!> Bug fixes for no thermostats [hforbert]
1710!> 2016-07-14 Modified to work with independent helium_env [cschran]
1711!> \author fawzi
1712! **************************************************************************************************
1713 SUBROUTINE pint_init_f(pint_env, helium_env)
1714 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1715 TYPE(helium_solvent_p_type), DIMENSION(:), &
1716 OPTIONAL, POINTER :: helium_env
1717
1718 INTEGER :: ib, idim, inos
1719 REAL(kind=dp) :: e_h
1720 TYPE(cp_logger_type), POINTER :: logger
1721
1722 NULLIFY (logger)
1723 logger => cp_get_default_logger()
1724
1725 ! initialize iteration info
1726 CALL cp_iterate(logger%iter_info, iter_nr=pint_env%first_step)
1727 CALL cp_iterate(pint_env%logger%iter_info, iter_nr=pint_env%first_step)
1728
1729 CALL pint_x2u(pint_env)
1730 CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h)
1731 CALL pint_calc_f(pint_env)
1732
1733 ! add helium forces to the solute's internal ones
1734 ! Assume that helium has been already initialized and helium_env(1)
1735 ! contains proper forces in force_avrg array at ionode
1736 IF (PRESENT(helium_env)) THEN
1737 IF (logger%para_env%is_source()) THEN
1738 pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
1739 END IF
1740 CALL logger%para_env%bcast(pint_env%f)
1741 END IF
1742 CALL pint_f2uf(pint_env)
1743
1744 ! set the centroid forces to 0 if FIX_CENTROID_POS
1745 IF (pint_env%first_propagated_mode == 2) THEN
1746 pint_env%uf(1, :) = 0.0_dp
1747 END IF
1748
1749 CALL pint_calc_e_kin_beads_u(pint_env)
1750 CALL pint_calc_e_vir(pint_env)
1751 DO idim = 1, SIZE(pint_env%uf_h, 2)
1752 DO ib = pint_env%first_propagated_mode, SIZE(pint_env%uf_h, 1)
1753 pint_env%uf(ib, idim) = real(pint_env%nrespa, dp)*pint_env%uf(ib, idim)
1754 END DO
1755 END DO
1756
1757 IF (pint_env%nnos > 0) THEN
1758 DO idim = 1, SIZE(pint_env%uf_h, 2)
1759 DO ib = 1, SIZE(pint_env%uf_h, 1)
1760 pint_env%tf(1, ib, idim) = (pint_env%mass_fict(ib, idim)* &
1761 pint_env%uv(ib, idim)**2 - pint_env%kT)/pint_env%Q(ib)
1762 END DO
1763 END DO
1764
1765 DO idim = 1, pint_env%ndim
1766 DO ib = 1, pint_env%p
1767 DO inos = 1, pint_env%nnos - 1
1768 pint_env%tf(inos + 1, ib, idim) = pint_env%tv(inos, ib, idim)**2 - &
1769 pint_env%kT/pint_env%Q(ib)
1770 END DO
1771 DO inos = 1, pint_env%nnos - 1
1772 pint_env%tf(inos, ib, idim) = pint_env%tf(inos, ib, idim) &
1773 - pint_env%tv(inos, ib, idim)*pint_env%tv(inos + 1, ib, idim)
1774 END DO
1775 END DO
1776 END DO
1777 CALL pint_calc_nh_energy(pint_env)
1778 END IF
1779
1780 END SUBROUTINE pint_init_f
1781
1782! ***************************************************************************
1783!> \brief Perform the PIMD simulation (main MD loop)
1784!> \param pint_env ...
1785!> \param globenv ...
1786!> \param helium_env ...
1787!> \par History
1788!> 2003-11 created [fawzi]
1789!> renamed from pint_run to pint_do_run because of conflicting name
1790!> of pint_run in input_constants [hforbert]
1791!> 2009-12-14 globenv parameter added to handle soft exit
1792!> requests [lwalewski]
1793!> 2016-07-14 Modified to work with independent helium_env [cschran]
1794!> \author Fawzi Mohamed
1795!> \note Everything should be read for an md step.
1796! **************************************************************************************************
1797 SUBROUTINE pint_do_run(pint_env, globenv, helium_env)
1798 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1799 TYPE(global_environment_type), POINTER :: globenv
1800 TYPE(helium_solvent_p_type), DIMENSION(:), &
1801 OPTIONAL, POINTER :: helium_env
1802
1803 INTEGER :: k, step
1804 LOGICAL :: should_stop
1805 REAL(kind=dp) :: scal
1806 TYPE(cp_logger_type), POINTER :: logger
1807 TYPE(f_env_type), POINTER :: f_env
1808
1809 ! initialize iteration info
1810 CALL cp_iterate(pint_env%logger%iter_info, iter_nr=pint_env%first_step)
1811
1812 ! iterate replica pint counter by accessing the globally saved
1813 ! force environment error/logger variables and setting them
1814 ! explicitly to the pimd "PINT" step value
1815 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, &
1816 f_env=f_env)
1817 NULLIFY (logger)
1818 logger => cp_get_default_logger()
1819 CALL cp_iterate(logger%iter_info, &
1820 iter_nr=pint_env%first_step)
1821 CALL f_env_rm_defaults(f_env)
1822
1823 pint_env%iter = pint_env%first_step
1824
1825 IF (PRESENT(helium_env)) THEN
1826 IF (ASSOCIATED(helium_env)) THEN
1827 ! set the properties accumulated over the whole MC process to 0
1828 DO k = 1, SIZE(helium_env)
1829 helium_env(k)%helium%proarea%accu(:) = 0.0_dp
1830 helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
1831 helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
1832 helium_env(k)%helium%mominer%accu(:) = 0.0_dp
1833 IF (helium_env(k)%helium%rho_present) THEN
1834 helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
1835 END IF
1836 IF (helium_env(k)%helium%rdf_present) THEN
1837 helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp
1838 END IF
1839 END DO
1840 END IF
1841 END IF
1842
1843 ! write the properties at 0-th step
1844 CALL pint_calc_energy(pint_env)
1845 CALL pint_calc_total_action(pint_env)
1846 CALL pint_write_ener(pint_env)
1847 CALL pint_write_action(pint_env)
1848 CALL pint_write_centroids(pint_env)
1849 CALL pint_write_trajectory(pint_env)
1850 CALL pint_write_com(pint_env)
1851 CALL pint_write_rgyr(pint_env)
1852
1853 ! main PIMD loop
1854 DO step = 1, pint_env%num_steps
1855
1856 pint_env%iter = pint_env%iter + 1
1857 CALL cp_iterate(pint_env%logger%iter_info, &
1858 last=(step == pint_env%num_steps), &
1859 iter_nr=pint_env%iter)
1860 CALL cp_iterate(logger%iter_info, &
1861 last=(step == pint_env%num_steps), &
1862 iter_nr=pint_env%iter)
1863 pint_env%t = pint_env%t + pint_env%dt
1864
1865 IF (pint_env%t_tol > 0.0_dp) THEN
1866 IF (abs(2._dp*pint_env%e_kin_beads/(pint_env%p*pint_env%ndim) &
1867 - pint_env%kT) > pint_env%t_tol) THEN
1868 scal = sqrt(pint_env%kT*(pint_env%p*pint_env%ndim)/(2.0_dp*pint_env%e_kin_beads))
1869 pint_env%uv = scal*pint_env%uv
1870 CALL pint_init_f(pint_env, helium_env=helium_env)
1871 END IF
1872 END IF
1873 CALL pint_step(pint_env, helium_env=helium_env)
1874
1875 CALL pint_write_ener(pint_env)
1876 CALL pint_write_action(pint_env)
1877 CALL pint_write_centroids(pint_env)
1878 CALL pint_write_trajectory(pint_env)
1879 CALL pint_write_com(pint_env)
1880 CALL pint_write_rgyr(pint_env)
1881
1882 CALL write_restart(root_section=pint_env%input, &
1883 pint_env=pint_env, helium_env=helium_env)
1884
1885 ! exit from the main loop if soft exit has been requested
1886 CALL external_control(should_stop, "PINT", globenv=globenv)
1887 IF (should_stop) EXIT
1888
1889 END DO
1890
1891 ! remove iteration level
1892 CALL cp_rm_iter_level(pint_env%logger%iter_info, "PINT")
1893
1894 END SUBROUTINE pint_do_run
1895
1896! ***************************************************************************
1897!> \brief Performs a scan of the helium-solute interaction energy
1898!> \param pint_env ...
1899!> \param helium_env ...
1900!> \date 2013-11-26
1901!> \parm History
1902!> 2016-07-14 Modified to work with independent helium_env [cschran]
1903!> \author Lukasz Walewski
1904! **************************************************************************************************
1905 SUBROUTINE pint_run_scan(pint_env, helium_env)
1906 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1907 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
1908
1909 CHARACTER(len=default_string_length) :: comment
1910 INTEGER :: unit_nr
1911 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: DATA
1912 TYPE(section_vals_type), POINTER :: print_key
1913
1914 NULLIFY (pint_env%logger, print_key)
1915 pint_env%logger => cp_get_default_logger()
1916
1917 ! assume that ionode always has at least one helium_env
1918 IF (pint_env%logger%para_env%is_source()) THEN
1919 print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
1920 "MOTION%PINT%HELIUM%PRINT%RHO")
1921 END IF
1922
1923 ! perform the actual scan wrt the COM of the solute
1924 CALL helium_intpot_scan(pint_env, helium_env)
1925
1926 ! output the interaction potential into a cubefile
1927 ! assume that ionode always has at least one helium_env
1928 IF (pint_env%logger%para_env%is_source()) THEN
1929
1930 unit_nr = cp_print_key_unit_nr( &
1931 pint_env%logger, &
1932 print_key, &
1933 middle_name="helium-pot", &
1934 extension=".cube", &
1935 file_position="REWIND", &
1936 do_backup=.false.)
1937
1938 comment = "Solute - helium interaction potential"
1939 NULLIFY (data)
1940 DATA => helium_env(1)%helium%rho_inst(1, :, :, :)
1941 CALL helium_write_cubefile( &
1942 unit_nr, &
1943 comment, &
1944 helium_env(1)%helium%center - 0.5_dp* &
1945 (helium_env(1)%helium%rho_maxr - helium_env(1)%helium%rho_delr), &
1946 helium_env(1)%helium%rho_delr, &
1947 helium_env(1)%helium%rho_nbin, &
1948 data)
1949
1950 CALL m_flush(unit_nr)
1951 CALL cp_print_key_finished_output(unit_nr, pint_env%logger, print_key)
1952
1953 END IF
1954
1955 ! output solute positions
1956 CALL pint_write_centroids(pint_env)
1957 CALL pint_write_trajectory(pint_env)
1958
1959 END SUBROUTINE pint_run_scan
1960
1961! ***************************************************************************
1962!> \brief Does an PINT step (and nrespa harmonic evaluations)
1963!> \param pint_env ...
1964!> \param helium_env ...
1965!> \par History
1966!> various bug fixes [hforbert]
1967!> 10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl]
1968!> 04.2016 Changed to work with helium_env [cschran]
1969!> 10.2018 Added centroid constraints [cschran+rperez]
1970!> 10.2021 Added beadwise constraints [lduran]
1971!> \author fawzi
1972! **************************************************************************************************
1973 SUBROUTINE pint_step(pint_env, helium_env)
1974 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
1975 TYPE(helium_solvent_p_type), DIMENSION(:), &
1976 OPTIONAL, POINTER :: helium_env
1977
1978 CHARACTER(len=*), PARAMETER :: routinen = 'pint_step'
1979
1980 INTEGER :: handle, i, ia, ib, idim, ierr, inos, &
1981 iresp, j, k, nbeads, nparticle, &
1982 nparticle_kind
1983 REAL(kind=dp) :: dt_temp, dti, dti2, dti22, e_h, factor, &
1984 rn, tdti, time_start, time_stop, tol
1985 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pos, vel
1986 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: tmp
1987 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1988 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1989 TYPE(cell_type), POINTER :: cell
1990 TYPE(cp_subsys_type), POINTER :: subsys
1991 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
1992 TYPE(f_env_type), POINTER :: f_env
1993 TYPE(global_constraint_type), POINTER :: gci
1994 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1995 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1996 TYPE(molecule_list_type), POINTER :: molecules
1997 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1998 TYPE(particle_list_type), POINTER :: particles
1999 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2000
2001 CALL timeset(routinen, handle)
2002 time_start = m_walltime()
2003
2004 rn = real(pint_env%nrespa, dp)
2005 dti = pint_env%dt/rn
2006 dti2 = dti/2._dp
2007 tdti = 2.*dti
2008 dti22 = dti**2/2._dp
2009
2010 ! Get constraint info, if needed
2011 ! Create a force environment which will be identical to
2012 ! the bead that is being processed by the processor.
2013 IF (pint_env%simpar%constraint) THEN
2014 NULLIFY (subsys, cell)
2015 NULLIFY (atomic_kinds, local_particles, particles)
2016 NULLIFY (local_molecules, molecules, molecule_kinds, gci)
2017 NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set)
2018
2019 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env)
2020 CALL force_env_get(force_env=f_env%force_env, subsys=subsys)
2021 CALL f_env_rm_defaults(f_env, ierr)
2022 cpassert(ierr == 0)
2023
2024 ! Get gci and more from subsys
2025 CALL cp_subsys_get(subsys=subsys, &
2026 cell=cell, &
2027 atomic_kinds=atomic_kinds, &
2028 local_particles=local_particles, &
2029 particles=particles, &
2030 local_molecules=local_molecules, &
2031 molecules=molecules, &
2032 molecule_kinds=molecule_kinds, &
2033 gci=gci)
2034
2035 nparticle_kind = atomic_kinds%n_els
2036 atomic_kind_set => atomic_kinds%els
2037 molecule_kind_set => molecule_kinds%els
2038 nparticle = particles%n_els
2039 nbeads = pint_env%p
2040 particle_set => particles%els
2041 molecule_set => molecules%els
2042
2043 ! Allocate work storage
2044 ALLOCATE (pos(3, nparticle))
2045 ALLOCATE (vel(3, nparticle))
2046 pos(:, :) = 0.0_dp
2047 vel(:, :) = 0.0_dp
2048
2049 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN
2050 ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD
2051 factor = sqrt(real(pint_env%p, dp))
2052 ELSE
2053 factor = 1.0_dp
2054 END IF
2055
2056 CALL getold(gci, local_molecules, molecule_set, &
2057 molecule_kind_set, particle_set, cell)
2058 END IF
2059
2060 SELECT CASE (pint_env%harm_integrator)
2061 CASE (integrate_numeric)
2062
2063 DO iresp = 1, pint_env%nrespa
2064
2065 ! integrate bead positions, first_propagated_mode = { 1, 2 }
2066 ! Nose needs an extra step
2067 IF (pint_env%pimd_thermostat == thermostat_nose) THEN
2068
2069 !Set thermostat action of constrained DoF to zero:
2070 IF (pint_env%simpar%constraint) THEN
2071 DO k = 1, pint_env%n_atoms_constraints
2072 ia = pint_env%atoms_constraints(k)
2073 DO j = 3*(ia - 1) + 1, 3*ia
2074 pint_env%tv(:, 1, j) = 0.0_dp
2075 END DO
2076 END DO
2077 END IF
2078
2079 ! Exempt centroid from thermostat for CMD
2080 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2081 pint_env%tx(:, 1, :) = 0.0_dp
2082 pint_env%tv(:, 1, :) = 0.0_dp
2083 pint_env%tf(:, 1, :) = 0.0_dp
2084 END IF
2085
2086 DO i = pint_env%first_propagated_mode, pint_env%p
2087 pint_env%ux(i, :) = pint_env%ux(i, :) - &
2088 dti22*pint_env%uv(i, :)*pint_env%tv(1, i, :)
2089 END DO
2090 pint_env%tx = pint_env%tx + dti*pint_env%tv + dti22*pint_env%tf
2091
2092 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2093 pint_env%tx(:, 1, :) = 0.0_dp
2094 pint_env%tv(:, 1, :) = 0.0_dp
2095 pint_env%tf(:, 1, :) = 0.0_dp
2096 END IF
2097
2098 END IF
2099 !Integrate position in harmonic springs (uf_h) and physical potential
2100 !(uf)
2101 DO i = pint_env%first_propagated_mode, pint_env%p
2102 pint_env%ux_t(i, :) = pint_env%ux(i, :) + &
2103 dti*pint_env%uv(i, :) + &
2104 dti22*(pint_env%uf_h(i, :) + &
2105 pint_env%uf(i, :))
2106 END DO
2107
2108 ! apply thermostats to velocities
2109 SELECT CASE (pint_env%pimd_thermostat)
2110 CASE (thermostat_nose)
2111
2112 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2113 pint_env%tx(:, 1, :) = 0.0_dp
2114 pint_env%tv(:, 1, :) = 0.0_dp
2115 pint_env%tf(:, 1, :) = 0.0_dp
2116 END IF
2117
2118 pint_env%uv_t = pint_env%uv - dti2* &
2119 pint_env%uv*pint_env%tv(1, :, :)
2120 tmp => pint_env%tv_t
2121 pint_env%tv_t => pint_env%tv
2122 pint_env%tv => tmp
2123 pint_env%tv = pint_env%tv_old + tdti*pint_env%tf
2124 pint_env%tv_old = pint_env%tv_t
2125 pint_env%tv_t = pint_env%tv_t + dti2*pint_env%tf
2126 CASE DEFAULT
2127 pint_env%uv_t = pint_env%uv
2128 END SELECT
2129
2130 !Set thermostat action of constrained DoF to zero:
2131 IF (pint_env%simpar%constraint) THEN
2132 DO k = 1, pint_env%n_atoms_constraints
2133 ia = pint_env%atoms_constraints(k)
2134 DO j = 3*(ia - 1) + 1, 3*ia
2135 pint_env%tv(:, 1, j) = 0.0_dp
2136 pint_env%tv_t(:, 1, j) = 0.0_dp
2137 END DO
2138 END DO
2139 END IF
2140
2141 ! Exempt centroid from thermostat for CMD
2142 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2143 pint_env%tx(:, 1, :) = 0.0_dp
2144 pint_env%tv(:, 1, :) = 0.0_dp
2145 pint_env%tf(:, 1, :) = 0.0_dp
2146 END IF
2147
2148 !Integrate harmonic velocities and physical velocities
2149 pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf)
2150
2151 ! physical forces are only applied in first respa step.
2152 pint_env%uf = 0.0_dp
2153 ! calc harmonic forces at new pos
2154 pint_env%ux = pint_env%ux_t
2155
2156 ! Apply centroid constraints (SHAKE)
2157 IF (pint_env%simpar%constraint) THEN
2158 IF (pint_env%logger%para_env%is_source()) THEN
2159 DO i = 1, nparticle
2160 DO j = 1, 3
2161 pos(j, i) = pint_env%ux(1, j + (i - 1)*3)
2162 vel(j, i) = pint_env%uv_t(1, j + (i - 1)*3)
2163 END DO
2164 END DO
2165
2166 ! Possibly update the target values
2167 CALL shake_update_targets(gci, local_molecules, molecule_set, &
2168 molecule_kind_set, dti, &
2169 f_env%force_env%root_section)
2170 CALL shake_control(gci, local_molecules, molecule_set, &
2171 molecule_kind_set, particle_set, &
2172 pos, vel, dti, pint_env%simpar%shake_tol, &
2173 pint_env%simpar%info_constraint, &
2174 pint_env%simpar%lagrange_multipliers, &
2175 pint_env%simpar%dump_lm, cell, &
2176 mp_comm_self, local_particles)
2177 END IF
2178 ! Positions and velocities of centroid were constrained by SHAKE
2179 CALL pint_env%logger%para_env%bcast(pos)
2180 CALL pint_env%logger%para_env%bcast(vel)
2181 ! Transform back to normal modes:
2182 DO i = 1, nparticle
2183 DO j = 1, 3
2184 pint_env%ux(1, j + (i - 1)*3) = pos(j, i)
2185 pint_env%uv_t(1, j + (i - 1)*3) = vel(j, i)
2186 END DO
2187 END DO
2188
2189 END IF
2190 ! Exempt centroid from thermostat for CMD
2191 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2192 pint_env%tx(:, 1, :) = 0.0_dp
2193 pint_env%tv(:, 1, :) = 0.0_dp
2194 pint_env%tf(:, 1, :) = 0.0_dp
2195 END IF
2196
2197 CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h)
2198 pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf)
2199
2200 ! For last respa step include integration of physical and helium
2201 ! forces
2202 IF (iresp == pint_env%nrespa) THEN
2203 CALL pint_u2x(pint_env)
2204 CALL pint_calc_f(pint_env)
2205 ! perform helium step and add helium forces
2206 IF (PRESENT(helium_env)) THEN
2207 CALL helium_step(helium_env, pint_env)
2208 !Update force of solute in pint_env
2209 IF (pint_env%logger%para_env%is_source()) THEN
2210 pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
2211 END IF
2212 CALL pint_env%logger%para_env%bcast(pint_env%f)
2213 END IF
2214
2215 CALL pint_f2uf(pint_env)
2216 ! set the centroid forces to 0 if FIX_CENTROID_POS
2217 IF (pint_env%first_propagated_mode == 2) THEN
2218 pint_env%uf(1, :) = 0.0_dp
2219 END IF
2220 !Scale physical forces and integrate velocities with physical
2221 !forces
2222 pint_env%uf = pint_env%uf*rn
2223 pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf
2224
2225 END IF
2226
2227 ! Apply second half of thermostats
2228 SELECT CASE (pint_env%pimd_thermostat)
2229 CASE (thermostat_nose)
2230 ! Exempt centroid from thermostat for CMD
2231 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2232 pint_env%tx(:, 1, :) = 0.0_dp
2233 pint_env%tv(:, 1, :) = 0.0_dp
2234 pint_env%tf(:, 1, :) = 0.0_dp
2235 END IF
2236 DO i = 1, 6
2237 tol = 0._dp
2238 pint_env%uv_new = pint_env%uv_t/(1.+dti2*pint_env%tv(1, :, :))
2239 DO idim = 1, pint_env%ndim
2240 DO ib = 1, pint_env%p
2241 pint_env%tf(1, ib, idim) = (pint_env%mass_fict(ib, idim)* &
2242 pint_env%uv_new(ib, idim)**2 - pint_env%kT*pint_env%kTcorr)/ &
2243 pint_env%Q(ib)
2244 END DO
2245 END DO
2246
2247 !Set thermostat action of constrained DoF to zero:
2248 IF (pint_env%simpar%constraint) THEN
2249 DO k = 1, pint_env%n_atoms_constraints
2250 ia = pint_env%atoms_constraints(k)
2251 DO j = 3*(ia - 1) + 1, 3*ia
2252 pint_env%tf(:, 1, j) = 0.0_dp
2253 END DO
2254 END DO
2255 END IF
2256
2257 ! Exempt centroid from thermostat for CMD
2258 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2259 pint_env%tx(:, 1, :) = 0.0_dp
2260 pint_env%tv(:, 1, :) = 0.0_dp
2261 pint_env%tf(:, 1, :) = 0.0_dp
2262 END IF
2263
2264 DO idim = 1, pint_env%ndim
2265 DO ib = 1, pint_env%p
2266 DO inos = 1, pint_env%nnos - 1
2267 pint_env%tv_new(inos, ib, idim) = &
2268 (pint_env%tv_t(inos, ib, idim) + dti2*pint_env%tf(inos, ib, idim))/ &
2269 (1._dp + dti2*pint_env%tv(inos + 1, ib, idim))
2270 pint_env%tf(inos + 1, ib, idim) = &
2271 (pint_env%tv_new(inos, ib, idim)**2 - &
2272 pint_env%kT*pint_env%kTcorr/pint_env%Q(ib))
2273 tol = max(tol, abs(pint_env%tv(inos, ib, idim) &
2274 - pint_env%tv_new(inos, ib, idim)))
2275 END DO
2276 !Set thermostat action of constrained DoF to zero:
2277 IF (pint_env%simpar%constraint) THEN
2278 DO k = 1, pint_env%n_atoms_constraints
2279 ia = pint_env%atoms_constraints(k)
2280 DO j = 3*(ia - 1) + 1, 3*ia
2281 pint_env%tv_new(:, 1, j) = 0.0_dp
2282 pint_env%tf(:, 1, j) = 0.0_dp
2283 END DO
2284 END DO
2285 END IF
2286
2287 ! Exempt centroid from thermostat for CMD
2288 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2289 pint_env%tx(:, 1, :) = 0.0_dp
2290 pint_env%tv(:, 1, :) = 0.0_dp
2291 pint_env%tf(:, 1, :) = 0.0_dp
2292 END IF
2293
2294 pint_env%tv_new(pint_env%nnos, ib, idim) = &
2295 pint_env%tv_t(pint_env%nnos, ib, idim) + &
2296 dti2*pint_env%tf(pint_env%nnos, ib, idim)
2297 tol = max(tol, abs(pint_env%tv(pint_env%nnos, ib, idim) &
2298 - pint_env%tv_new(pint_env%nnos, ib, idim)))
2299 tol = max(tol, abs(pint_env%uv(ib, idim) &
2300 - pint_env%uv_new(ib, idim)))
2301 !Set thermostat action of constrained DoF to zero:
2302 IF (pint_env%simpar%constraint) THEN
2303 DO k = 1, pint_env%n_atoms_constraints
2304 ia = pint_env%atoms_constraints(k)
2305 DO j = 3*(ia - 1) + 1, 3*ia
2306 pint_env%tv_new(:, 1, j) = 0.0_dp
2307 END DO
2308 END DO
2309 END IF
2310 ! Exempt centroid from thermostat for CMD
2311 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2312 pint_env%tx(:, 1, :) = 0.0_dp
2313 pint_env%tv(:, 1, :) = 0.0_dp
2314 pint_env%tf(:, 1, :) = 0.0_dp
2315 END IF
2316
2317 END DO
2318 END DO
2319
2320 pint_env%uv = pint_env%uv_new
2321 pint_env%tv = pint_env%tv_new
2322 IF (tol <= pint_env%v_tol) EXIT
2323 ! Exempt centroid from thermostat for CMD
2324 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2325 pint_env%tx(:, 1, :) = 0.0_dp
2326 pint_env%tv(:, 1, :) = 0.0_dp
2327 pint_env%tf(:, 1, :) = 0.0_dp
2328 END IF
2329 END DO
2330
2331 ! Apply centroid constraints (RATTLE)
2332 IF (pint_env%simpar%constraint) THEN
2333 IF (pint_env%logger%para_env%is_source()) THEN
2334 ! Reset particle r, due to force calc:
2335 DO i = 1, nparticle
2336 DO j = 1, 3
2337 vel(j, i) = pint_env%uv(1, j + (i - 1)*3)
2338 particle_set(i)%r(j) = pint_env%ux(1, j + (i - 1)*3)
2339 END DO
2340 END DO
2341
2342 ! Small time step for all small integrations steps
2343 ! Big step for last RESPA
2344 IF (iresp == pint_env%nrespa) THEN
2345 dt_temp = dti
2346 ELSE
2347 dt_temp = dti*rn
2348 END IF
2349 CALL rattle_control(gci, local_molecules, molecule_set, &
2350 molecule_kind_set, particle_set, &
2351 vel, dt_temp, pint_env%simpar%shake_tol, &
2352 pint_env%simpar%info_constraint, &
2353 pint_env%simpar%lagrange_multipliers, &
2354 pint_env%simpar%dump_lm, cell, &
2355 mp_comm_self, local_particles)
2356 END IF
2357 ! Velocities of centroid were constrained by RATTLE
2358 ! Broadcast updated velocities to other nodes
2359 CALL pint_env%logger%para_env%bcast(vel)
2360
2361 DO i = 1, nparticle
2362 DO j = 1, 3
2363 pint_env%uv(1, j + (i - 1)*3) = vel(j, i)
2364 END DO
2365 END DO
2366 END IF
2367
2368 DO inos = 1, pint_env%nnos - 1
2369 pint_env%tf(inos, :, :) = pint_env%tf(inos, :, :) &
2370 - pint_env%tv(inos, :, :)*pint_env%tv(inos + 1, :, :)
2371 END DO
2372
2373 ! Exempt centroid from thermostat for CMD
2374 IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
2375 pint_env%tx(:, 1, :) = 0.0_dp
2376 pint_env%tv(:, 1, :) = 0.0_dp
2377 pint_env%tf(:, 1, :) = 0.0_dp
2378 END IF
2379
2380 CASE (thermostat_gle)
2381 CALL pint_gle_step(pint_env)
2382 pint_env%uv = pint_env%uv_t
2383 CASE DEFAULT
2384 pint_env%uv = pint_env%uv_t
2385 END SELECT
2386 END DO
2387
2388 CASE (integrate_exact)
2389 ! The Liouvillian splitting is as follows:
2390 ! 1. Thermostat
2391 ! 2. 0.5*physical integration
2392 ! 3. Exact harmonic integration + apply constraints (SHAKE)
2393 ! 4. 0.5*physical integration
2394 ! 5. Thermostat + apply constraints (RATTLE)
2395
2396 ! 1. Apply thermostats
2397 SELECT CASE (pint_env%pimd_thermostat)
2398 CASE (thermostat_pile)
2399 CALL pint_pile_step(vold=pint_env%uv, &
2400 vnew=pint_env%uv_t, &
2401 p=pint_env%p, &
2402 ndim=pint_env%ndim, &
2403 first_mode=pint_env%first_propagated_mode, &
2404 masses=pint_env%mass_fict, &
2405 pile_therm=pint_env%pile_therm)
2406 CASE (thermostat_piglet)
2407 CALL pint_piglet_step(vold=pint_env%uv, &
2408 vnew=pint_env%uv_t, &
2409 first_mode=pint_env%first_propagated_mode, &
2410 masses=pint_env%mass_fict, &
2411 piglet_therm=pint_env%piglet_therm)
2412 CASE (thermostat_qtb)
2413 CALL pint_qtb_step(vold=pint_env%uv, &
2414 vnew=pint_env%uv_t, &
2415 p=pint_env%p, &
2416 ndim=pint_env%ndim, &
2417 masses=pint_env%mass_fict, &
2418 qtb_therm=pint_env%qtb_therm)
2419 CASE DEFAULT
2420 pint_env%uv_t = pint_env%uv
2421 END SELECT
2422
2423 ! 2. 1/2*Physical integration
2424 pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf
2425
2426 ! 3. Exact harmonic integration
2427 IF (pint_env%first_propagated_mode == 1) THEN
2428 ! The centroid is integrated via standard velocity-verlet
2429 ! Commented out code is only there to show similarities to
2430 ! Numeric integrator
2431 pint_env%ux_t(1, :) = pint_env%ux(1, :) + &
2432 dti*pint_env%uv_t(1, :) !+ &
2433 ! dti22*pint_env%uf_h(1, :)
2434 !pint_env%uv_t(1, :) = pint_env%uv_t(1, :)+ &
2435 ! dti2*pint_env%uf_h(1, :)
2436 ELSE
2437 ! set velocities to zero for fixed centroids
2438 pint_env%ux_t(1, :) = pint_env%ux(1, :)
2439 pint_env%uv_t(1, :) = 0.0_dp
2440 END IF
2441 ! Other modes are integrated exactly
2442 DO i = 2, pint_env%p
2443 pint_env%ux_t(i, :) = pint_env%cosex(i)*pint_env%ux(i, :) &
2444 + pint_env%iwsinex(i)*pint_env%uv_t(i, :)
2445 pint_env%uv_t(i, :) = pint_env%cosex(i)*pint_env%uv_t(i, :) &
2446 - pint_env%wsinex(i)*pint_env%ux(i, :)
2447 END DO
2448
2449 ! Apply constraints (SHAKE)
2450 IF (pint_env%simpar%constraint) THEN
2451 ! Beadwise constraints
2452 IF (pint_env%beadwise_constraints) THEN
2453 IF (pint_env%logger%para_env%is_source()) THEN
2454 ! Transform positions and velocities to Cartesian coordinates:
2455 CALL pint_u2x(pint_env, ux=pint_env%ux_t, x=pint_env%x)
2456 CALL pint_u2x(pint_env, ux=pint_env%uv_t, x=pint_env%v)
2457 DO ib = 1, nbeads
2458 DO i = 1, nparticle
2459 DO j = 1, 3
2460 pos(j, i) = pint_env%x(ib, j + (i - 1)*3)
2461 vel(j, i) = pint_env%v(ib, j + (i - 1)*3)
2462 END DO
2463 END DO
2464 ! Possibly update the target values
2465 CALL shake_update_targets(gci, local_molecules, molecule_set, &
2466 molecule_kind_set, dti, &
2467 f_env%force_env%root_section)
2468 CALL shake_control(gci, local_molecules, molecule_set, &
2469 molecule_kind_set, particle_set, &
2470 pos, vel, dti, pint_env%simpar%shake_tol, &
2471 pint_env%simpar%info_constraint, &
2472 pint_env%simpar%lagrange_multipliers, &
2473 pint_env%simpar%dump_lm, cell, &
2474 mp_comm_self, local_particles)
2475 DO i = 1, nparticle
2476 DO j = 1, 3
2477 pint_env%x(ib, j + (i - 1)*3) = pos(j, i)
2478 pint_env%v(ib, j + (i - 1)*3) = vel(j, i)
2479 END DO
2480 END DO
2481 END DO
2482 ! Transform back to normal modes:
2483 CALL pint_x2u(pint_env, ux=pint_env%ux_t, x=pint_env%x)
2484 CALL pint_x2u(pint_env, ux=pint_env%uv_t, x=pint_env%v)
2485 END IF
2486 ! Broadcast positions and velocities to all nodes
2487 CALL pint_env%logger%para_env%bcast(pint_env%ux_t)
2488 CALL pint_env%logger%para_env%bcast(pint_env%uv_t)
2489 ! Centroid constraints
2490 ELSE
2491 IF (pint_env%logger%para_env%is_source()) THEN
2492 ! Transform positions and velocities to Cartesian coordinates:
2493 DO i = 1, nparticle
2494 DO j = 1, 3
2495 pos(j, i) = pint_env%ux_t(1, j + (i - 1)*3)/factor
2496 vel(j, i) = pint_env%uv_t(1, j + (i - 1)*3)/factor
2497 END DO
2498 END DO
2499 ! Possibly update the target values
2500 CALL shake_update_targets(gci, local_molecules, molecule_set, &
2501 molecule_kind_set, dti, &
2502 f_env%force_env%root_section)
2503 CALL shake_control(gci, local_molecules, molecule_set, &
2504 molecule_kind_set, particle_set, &
2505 pos, vel, dti, pint_env%simpar%shake_tol, &
2506 pint_env%simpar%info_constraint, &
2507 pint_env%simpar%lagrange_multipliers, &
2508 pint_env%simpar%dump_lm, cell, &
2509 mp_comm_self, local_particles)
2510 END IF
2511 ! Broadcast positions and velocities to all nodes
2512 CALL pint_env%logger%para_env%bcast(pos)
2513 CALL pint_env%logger%para_env%bcast(vel)
2514 ! Transform back to normal modes:
2515 DO i = 1, nparticle
2516 DO j = 1, 3
2517 pint_env%ux_t(1, j + (i - 1)*3) = pos(j, i)*factor
2518 pint_env%uv_t(1, j + (i - 1)*3) = vel(j, i)*factor
2519 END DO
2520 END DO
2521 END IF
2522 ! Positions and velocities were constrained by SHAKE
2523 END IF
2524 ! Update positions
2525 pint_env%ux = pint_env%ux_t
2526
2527 ! 4. 1/2*Physical integration
2528 pint_env%uf = 0.0_dp
2529 CALL pint_u2x(pint_env)
2530 CALL pint_calc_f(pint_env)
2531 ! perform helium step and add helium forces
2532 IF (PRESENT(helium_env)) THEN
2533 CALL helium_step(helium_env, pint_env)
2534 !Update force of solute in pint_env
2535 IF (pint_env%logger%para_env%is_source()) THEN
2536 pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :)
2537 END IF
2538 CALL pint_env%logger%para_env%bcast(pint_env%f)
2539 END IF
2540 CALL pint_f2uf(pint_env)
2541 ! set the centroid forces to 0 if FIX_CENTROID_POS
2542 IF (pint_env%first_propagated_mode == 2) THEN
2543 pint_env%uf(1, :) = 0.0_dp
2544 END IF
2545 pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf
2546
2547 ! 5. Apply thermostats
2548 SELECT CASE (pint_env%pimd_thermostat)
2549 CASE (thermostat_pile)
2550 CALL pint_pile_step(vold=pint_env%uv_t, &
2551 vnew=pint_env%uv, &
2552 p=pint_env%p, &
2553 ndim=pint_env%ndim, &
2554 first_mode=pint_env%first_propagated_mode, &
2555 masses=pint_env%mass_fict, &
2556 pile_therm=pint_env%pile_therm)
2557 CASE (thermostat_piglet)
2558 CALL pint_piglet_step(vold=pint_env%uv_t, &
2559 vnew=pint_env%uv, &
2560 first_mode=pint_env%first_propagated_mode, &
2561 masses=pint_env%mass_fict, &
2562 piglet_therm=pint_env%piglet_therm)
2563 CASE (thermostat_qtb)
2564 CALL pint_qtb_step(vold=pint_env%uv_t, &
2565 vnew=pint_env%uv, &
2566 p=pint_env%p, &
2567 ndim=pint_env%ndim, &
2568 masses=pint_env%mass_fict, &
2569 qtb_therm=pint_env%qtb_therm)
2570 CASE DEFAULT
2571 pint_env%uv = pint_env%uv_t
2572 END SELECT
2573
2574 ! Apply constraints (RATTLE)
2575 IF (pint_env%simpar%constraint) THEN
2576 ! Beadwise constraints
2577 IF (pint_env%beadwise_constraints) THEN
2578 IF (pint_env%logger%para_env%is_source()) THEN
2579 ! Transform positions and velocities to Cartesian coordinates:
2580 ! Reset particle r, due to force calc:
2581 CALL pint_u2x(pint_env, ux=pint_env%ux, x=pint_env%x)
2582 CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v)
2583 DO ib = 1, nbeads
2584 DO i = 1, nparticle
2585 DO j = 1, 3
2586 particle_set(i)%r(j) = pint_env%x(ib, j + (i - 1)*3)
2587 vel(j, i) = pint_env%v(ib, j + (i - 1)*3)
2588 END DO
2589 END DO
2590 CALL rattle_control(gci, local_molecules, &
2591 molecule_set, molecule_kind_set, &
2592 particle_set, vel, dti, &
2593 pint_env%simpar%shake_tol, &
2594 pint_env%simpar%info_constraint, &
2595 pint_env%simpar%lagrange_multipliers, &
2596 pint_env%simpar%dump_lm, cell, &
2597 mp_comm_self, local_particles)
2598 DO i = 1, nparticle
2599 DO j = 1, 3
2600 pint_env%v(ib, j + (i - 1)*3) = vel(j, i)
2601 END DO
2602 END DO
2603 END DO
2604 ! Transform back to normal modes:
2605 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v)
2606 END IF
2607 CALL pint_env%logger%para_env%bcast(pint_env%uv)
2608 ! Centroid constraints
2609 ELSE
2610 IF (pint_env%logger%para_env%is_source()) THEN
2611 ! Transform positions and velocities to Cartesian coordinates:
2612 ! Reset particle r, due to force calc:
2613 DO i = 1, nparticle
2614 DO j = 1, 3
2615 vel(j, i) = pint_env%uv(1, j + (i - 1)*3)/factor
2616 particle_set(i)%r(j) = pint_env%ux(1, j + (i - 1)*3)/factor
2617 END DO
2618 END DO
2619 CALL rattle_control(gci, local_molecules, &
2620 molecule_set, molecule_kind_set, &
2621 particle_set, vel, dti, &
2622 pint_env%simpar%shake_tol, &
2623 pint_env%simpar%info_constraint, &
2624 pint_env%simpar%lagrange_multipliers, &
2625 pint_env%simpar%dump_lm, cell, &
2626 mp_comm_self, local_particles)
2627 END IF
2628 ! Velocities of centroid were constrained by RATTLE
2629 ! Broadcast updated velocities to other nodes
2630 CALL pint_env%logger%para_env%bcast(vel)
2631
2632 ! Transform back to normal modes:
2633 ! Multiply with SQRT(n_beads) due to normal mode transformation
2634 DO i = 1, nparticle
2635 DO j = 1, 3
2636 pint_env%uv(1, j + (i - 1)*3) = vel(j, i)*factor
2637 END DO
2638 END DO
2639 END IF
2640 END IF
2641
2642 END SELECT
2643
2644 IF (pint_env%simpar%constraint) THEN
2645 DEALLOCATE (pos, vel)
2646 END IF
2647
2648 ! calculate the energy components
2649 CALL pint_calc_energy(pint_env)
2650 CALL pint_calc_total_action(pint_env)
2651
2652 ! check that the number of PINT steps matches
2653 ! the number of force evaluations done so far
2654!TODO make this check valid if we start from ITERATION != 0
2655! CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id,&
2656! f_env=f_env,new_error=new_error)
2657! NULLIFY(logger)
2658! logger => cp_error_get_logger(new_error)
2659! IF(logger%iter_info%iteration(2)/=pint_env%iter+1)&
2660! CPABORT("md & force_eval lost sychro")
2661! CALL f_env_rm_defaults(f_env,new_error,ierr)
2662
2663 time_stop = m_walltime()
2664 pint_env%time_per_step = time_stop - time_start
2665 CALL pint_write_step_info(pint_env)
2666 CALL timestop(handle)
2667
2668 END SUBROUTINE pint_step
2669
2670! ***************************************************************************
2671!> \brief Calculate the energy components (private wrapper function)
2672!> \param pint_env ...
2673!> \date 2011-01-07
2674!> \author Lukasz Walewski
2675! **************************************************************************************************
2676 SUBROUTINE pint_calc_energy(pint_env)
2677
2678 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
2679
2680 REAL(kind=dp) :: e_h
2681
2682 CALL pint_calc_e_kin_beads_u(pint_env)
2683 CALL pint_calc_e_vir(pint_env)
2684
2685 CALL pint_calc_uf_h(pint_env, e_h=e_h)
2686 pint_env%e_pot_h = e_h
2687
2688 SELECT CASE (pint_env%pimd_thermostat)
2689 CASE (thermostat_nose)
2690 CALL pint_calc_nh_energy(pint_env)
2691 CASE (thermostat_gle)
2692 CALL pint_calc_gle_energy(pint_env)
2693 CASE (thermostat_pile)
2694 CALL pint_calc_pile_energy(pint_env)
2695 CASE (thermostat_qtb)
2696 CALL pint_calc_qtb_energy(pint_env)
2697 CASE (thermostat_piglet)
2698 CALL pint_calc_piglet_energy(pint_env)
2699 END SELECT
2700
2701 pint_env%energy(e_kin_thermo_id) = &
2702 (0.5_dp*real(pint_env%p, dp)*real(pint_env%ndim, dp)*pint_env%kT - &
2703 pint_env%e_pot_h)*pint_env%propagator%temp_sim2phys
2704
2705 pint_env%energy(e_potential_id) = sum(pint_env%e_pot_bead)
2706
2707 pint_env%energy(e_conserved_id) = &
2708 pint_env%energy(e_potential_id)*pint_env%propagator%physpotscale + &
2709 pint_env%e_pot_h + &
2710 pint_env%e_kin_beads + &
2711 pint_env%e_pot_t + &
2712 pint_env%e_kin_t + &
2713 pint_env%e_gle + pint_env%e_pile + pint_env%e_piglet + pint_env%e_qtb
2714
2715 pint_env%energy(e_potential_id) = &
2716 pint_env%energy(e_potential_id)/real(pint_env%p, dp)
2717
2718 END SUBROUTINE pint_calc_energy
2719
2720! ***************************************************************************
2721!> \brief Calculate the harmonic force in the u basis
2722!> \param pint_env the path integral environment in which the harmonic
2723!> forces should be calculated
2724!> \param e_h ...
2725!> \par History
2726!> Added normal mode transformation [hforbert]
2727!> \author fawzi
2728! **************************************************************************************************
2729 SUBROUTINE pint_calc_uf_h(pint_env, e_h)
2730 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
2731 REAL(kind=dp), INTENT(OUT) :: e_h
2732
2733 IF (pint_env%transform == transformation_stage) THEN
2734 CALL staging_calc_uf_h(pint_env%staging_env, &
2735 pint_env%mass_beads, &
2736 pint_env%ux, &
2737 pint_env%uf_h, &
2738 pint_env%e_pot_h)
2739 ELSE
2740 CALL normalmode_calc_uf_h(pint_env%normalmode_env, &
2741 pint_env%mass_beads, &
2742 pint_env%ux, &
2743 pint_env%uf_h, &
2744 pint_env%e_pot_h)
2745 END IF
2746 e_h = pint_env%e_pot_h
2747 pint_env%uf_h = pint_env%uf_h/pint_env%mass_fict
2748 END SUBROUTINE pint_calc_uf_h
2749
2750! ***************************************************************************
2751!> \brief calculates the force (and energy) in each bead, returns the sum
2752!> of the potential energy
2753!> \param pint_env path integral environment on which you want to calculate
2754!> the forces
2755!> \param x positions at which you want to evaluate the forces
2756!> \param f the forces
2757!> \param e potential energy on each bead
2758!> \par History
2759!> 2009-06-15 moved helium calls out from here [lwalewski]
2760!> \author fawzi
2761! **************************************************************************************************
2762 SUBROUTINE pint_calc_f(pint_env, x, f, e)
2763 TYPE(pint_env_type), INTENT(IN) :: pint_env
2764 REAL(kind=dp), DIMENSION(:, :), INTENT(in), &
2765 OPTIONAL, TARGET :: x
2766 REAL(kind=dp), DIMENSION(:, :), INTENT(out), &
2767 OPTIONAL, TARGET :: f
2768 REAL(kind=dp), DIMENSION(:), INTENT(out), &
2769 OPTIONAL, TARGET :: e
2770
2771 INTEGER :: ib, idim
2772 REAL(kind=dp), DIMENSION(:), POINTER :: my_e
2773 REAL(kind=dp), DIMENSION(:, :), POINTER :: my_f, my_x
2774
2775 my_x => pint_env%x
2776 IF (PRESENT(x)) my_x => x
2777 my_f => pint_env%f
2778 IF (PRESENT(f)) my_f => f
2779 my_e => pint_env%e_pot_bead
2780 IF (PRESENT(e)) my_e => e
2781 DO idim = 1, pint_env%ndim
2782 DO ib = 1, pint_env%p
2783 pint_env%replicas%r(idim, ib) = my_x(ib, idim)
2784 END DO
2785 END DO
2786 CALL rep_env_calc_e_f(pint_env%replicas, calc_f=.true.)
2787 DO idim = 1, pint_env%ndim
2788 DO ib = 1, pint_env%p
2789 !ljw: is that fine ? - idim <-> ib
2790 my_f(ib, idim) = pint_env%replicas%f(idim, ib)
2791 END DO
2792 END DO
2793 my_e = pint_env%replicas%f(SIZE(pint_env%replicas%f, 1), :)
2794
2795 END SUBROUTINE pint_calc_f
2796
2797! ***************************************************************************
2798!> \brief Calculate the kinetic energy of the beads (in the u variables)
2799!> \param pint_env ...
2800!> \param uv ...
2801!> \param e_k ...
2802!> \par History
2803!> Bug fix to give my_uv a default location if not given in call [hforbert]
2804!> \author fawzi
2805! **************************************************************************************************
2806 SUBROUTINE pint_calc_e_kin_beads_u(pint_env, uv, e_k)
2807 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
2808 REAL(kind=dp), DIMENSION(:, :), INTENT(in), &
2809 OPTIONAL, TARGET :: uv
2810 REAL(kind=dp), INTENT(out), OPTIONAL :: e_k
2811
2812 INTEGER :: ib, idim
2813 REAL(kind=dp) :: res
2814 REAL(kind=dp), DIMENSION(:, :), POINTER :: my_uv
2815
2816 res = -1.0_dp
2817 my_uv => pint_env%uv
2818 IF (PRESENT(uv)) my_uv => uv
2819 res = 0._dp
2820 DO idim = 1, pint_env%ndim
2821 DO ib = 1, pint_env%p
2822 res = res + pint_env%mass_fict(ib, idim)*my_uv(ib, idim)**2
2823 END DO
2824 END DO
2825 res = res*0.5
2826 IF (.NOT. PRESENT(uv)) pint_env%e_kin_beads = res
2827 IF (PRESENT(e_k)) e_k = res
2828 END SUBROUTINE pint_calc_e_kin_beads_u
2829
2830! ***************************************************************************
2831!> \brief Calculate the virial estimator of the real (quantum) kinetic energy
2832!> \param pint_env ...
2833!> \param e_vir ...
2834!> \author hforbert
2835!> \note This subroutine modifies pint_env%energy(e_kin_virial_id) global
2836!> variable [lwalewski]
2837! **************************************************************************************************
2838 ELEMENTAL SUBROUTINE pint_calc_e_vir(pint_env, e_vir)
2839 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
2840 REAL(kind=dp), INTENT(out), OPTIONAL :: e_vir
2841
2842 INTEGER :: ib, idim
2843 REAL(kind=dp) :: res, xcentroid
2844
2845 res = -1.0_dp
2846 res = 0._dp
2847 DO idim = 1, pint_env%ndim
2848 ! calculate the centroid
2849 xcentroid = 0._dp
2850 DO ib = 1, pint_env%p
2851 xcentroid = xcentroid + pint_env%x(ib, idim)
2852 END DO
2853 xcentroid = xcentroid/real(pint_env%p, dp)
2854 DO ib = 1, pint_env%p
2855 res = res + (pint_env%x(ib, idim) - xcentroid)*pint_env%f(ib, idim)
2856 END DO
2857 END DO
2858 res = 0.5_dp*(real(pint_env%ndim, dp)* &
2859 (pint_env%kT*pint_env%propagator%temp_sim2phys) - res/real(pint_env%p, dp))
2860 pint_env%energy(e_kin_virial_id) = res
2861 IF (PRESENT(e_vir)) e_vir = res
2862 END SUBROUTINE pint_calc_e_vir
2863
2864! ***************************************************************************
2865!> \brief calculates the energy (potential and kinetic) of the Nose-Hoover
2866!> chain thermostats
2867!> \param pint_env the path integral environment
2868!> \author fawzi
2869! **************************************************************************************************
2870 ELEMENTAL SUBROUTINE pint_calc_nh_energy(pint_env)
2871 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
2872
2873 INTEGER :: ib, idim, inos
2874 REAL(kind=dp) :: ekin, epot
2875
2876 ekin = 0._dp
2877 DO idim = 1, pint_env%ndim
2878 DO ib = 1, pint_env%p
2879 DO inos = 1, pint_env%nnos
2880 ekin = ekin + pint_env%Q(ib)*pint_env%tv(inos, ib, idim)**2
2881 END DO
2882 END DO
2883 END DO
2884 pint_env%e_kin_t = 0.5_dp*ekin
2885 epot = 0._dp
2886 DO idim = 1, pint_env%ndim
2887 DO ib = 1, pint_env%p
2888 DO inos = 1, pint_env%nnos
2889 epot = epot + pint_env%tx(inos, ib, idim)
2890 END DO
2891 END DO
2892 END DO
2893 pint_env%e_pot_t = pint_env%kT*epot
2894 END SUBROUTINE pint_calc_nh_energy
2895
2896! ***************************************************************************
2897!> \brief calculates the total link action of the PI system (excluding helium)
2898!> \param pint_env the path integral environment
2899!> \return ...
2900!> \author Felix Uhl
2901! **************************************************************************************************
2902 ELEMENTAL FUNCTION pint_calc_total_link_action(pint_env) RESULT(link_action)
2903 TYPE(pint_env_type), INTENT(IN) :: pint_env
2904 REAL(kind=dp) :: link_action
2905
2906 INTEGER :: iatom, ibead, idim, indx
2907 REAL(kind=dp) :: hb2m, tau, tmp_link_action
2908 REAL(kind=dp), DIMENSION(3) :: r
2909
2910 !tau = 1/(k_B T p)
2911 tau = pint_env%beta/real(pint_env%p, dp)
2912
2913 link_action = 0.0_dp
2914 DO iatom = 1, pint_env%ndim/3
2915 ! hbar / (2.0*m)
2916 hb2m = 1.0_dp/pint_env%mass((iatom - 1)*3 + 1)
2917 tmp_link_action = 0.0_dp
2918 DO ibead = 1, pint_env%p - 1
2919 DO idim = 1, 3
2920 indx = (iatom - 1)*3 + idim
2921 r(idim) = pint_env%x(ibead, indx) - pint_env%x(ibead + 1, indx)
2922 END DO
2923 tmp_link_action = tmp_link_action + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
2924 END DO
2925 DO idim = 1, 3
2926 indx = (iatom - 1)*3 + idim
2927 r(idim) = pint_env%x(pint_env%p, indx) - pint_env%x(1, indx)
2928 END DO
2929 tmp_link_action = tmp_link_action + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
2930 link_action = link_action + tmp_link_action/hb2m
2931 END DO
2932
2933 link_action = link_action/(2.0_dp*tau)
2934
2935 END FUNCTION pint_calc_total_link_action
2936
2937! ***************************************************************************
2938!> \brief calculates the potential action of the PI system (excluding helium)
2939!> \param pint_env the path integral environment
2940!> \return ...
2941!> \author Felix Uhl
2942! **************************************************************************************************
2943 ELEMENTAL FUNCTION pint_calc_total_pot_action(pint_env) RESULT(pot_action)
2944 TYPE(pint_env_type), INTENT(IN) :: pint_env
2945 REAL(kind=dp) :: pot_action
2946
2947 REAL(kind=dp) :: tau
2948
2949 tau = pint_env%beta/real(pint_env%p, dp)
2950 pot_action = tau*sum(pint_env%e_pot_bead)
2951
2952 END FUNCTION pint_calc_total_pot_action
2953
2954! ***************************************************************************
2955!> \brief calculates the total action of the PI system (excluding helium)
2956!> \param pint_env the path integral environment
2957!> \author Felix Uhl
2958! **************************************************************************************************
2959 ELEMENTAL SUBROUTINE pint_calc_total_action(pint_env)
2960 TYPE(pint_env_type), INTENT(INOUT) :: pint_env
2961
2962 pint_env%pot_action = pint_calc_total_pot_action(pint_env)
2963 pint_env%link_action = pint_calc_total_link_action(pint_env)
2964
2965 END SUBROUTINE pint_calc_total_action
2966
2967END MODULE pint_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
atomic_kind_list_types
represent a simple array based list of the given type
Definition atomic_kind_list_types.F:15
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::shiga2022
integer, save, public shiga2022
Definition bibliography.F:36
bibliography::ceriotti2012
integer, save, public ceriotti2012
Definition bibliography.F:36
bibliography::ceriotti2010
integer, save, public ceriotti2010
Definition bibliography.F:36
bibliography::brieuc2016
integer, save, public brieuc2016
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
constraint_util
Contains routines useful for the application of constraints during MD.
Definition constraint_util.F:13
constraint_util::getold
subroutine, public getold(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, cell)
saves all of the old variables
Definition constraint_util.F:63
constraint::rattle_control
subroutine, public rattle_control(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, vel, dt, rattle_tol, log_unit, lagrange_mult, dump_lm, cell, group, local_particles)
...
Definition constraint.F:229
constraint::shake_control
subroutine, public shake_control(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, pos, vel, dt, shake_tol, log_unit, lagrange_mult, dump_lm, cell, group, local_particles)
...
Definition constraint.F:101
constraint::shake_update_targets
subroutine, public shake_update_targets(gci, local_molecules, molecule_set, molecule_kind_set, dt, root_section)
Updates the TARGET of the COLLECTIVE constraints if the growth speed is different from zero.
Definition constraint.F:843
cp_external_control
Routines to handle the external control of CP2K.
Definition cp_external_control.F:15
cp_external_control::external_control
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
Definition cp_external_control.F:90
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1067
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_iterate
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
Definition cp_output_handling.F:574
cp_output_handling::cp_rm_iter_level
subroutine, public cp_rm_iter_level(iteration_info, level_name, n_rlevel_att)
Removes an iteration level.
Definition cp_output_handling.F:705
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_output_handling::cp_add_iter_level
subroutine, public cp_add_iter_level(iteration_info, level_name, n_rlevel_new)
Adds an iteration level.
Definition cp_output_handling.F:662
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1239
cp_units::cp_unit_to_cp2k
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1210
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
f77_interface
interface to use cp2k as library
Definition f77_interface.F:20
f77_interface::f_env_add_defaults
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
Definition f77_interface.F:477
f77_interface::f_env_rm_defaults
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
Definition f77_interface.F:521
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
gle_system_dynamics
Definition gle_system_dynamics.F:12
gle_system_dynamics::restart_gle
subroutine, public restart_gle(gle, gle_section, save_mem, restart)
...
Definition gle_system_dynamics.F:417
gle_system_dynamics::gle_matrix_exp
subroutine, public gle_matrix_exp(m, n, j, k, em)
...
Definition gle_system_dynamics.F:231
gle_system_dynamics::gle_cholesky_stab
subroutine, public gle_cholesky_stab(sst, s, n)
...
Definition gle_system_dynamics.F:278
gle_system_types
Definition gle_system_types.F:13
gle_system_types::gle_dealloc
subroutine, public gle_dealloc(gle)
Deallocate type for GLE thermostat.
Definition gle_system_types.F:217
gle_system_types::gle_thermo_create
subroutine, public gle_thermo_create(gle, mal_size)
...
Definition gle_system_types.F:166
gle_system_types::gle_init
subroutine, public gle_init(gle, dt, temp, section)
...
Definition gle_system_types.F:69
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
helium_interactions
Methods that handle helium-solvent and helium-helium interactions.
Definition helium_interactions.F:13
helium_interactions::helium_intpot_scan
subroutine, public helium_intpot_scan(pint_env, helium_env)
Scan the helium-solute interaction energy within the periodic cell.
Definition helium_interactions.F:517
helium_io
I/O subroutines for helium.
Definition helium_io.F:13
helium_io::helium_write_cubefile
subroutine, public helium_write_cubefile(unit, comment, origin, deltar, ndim, data)
Write volumetric data to an orthorhombic cubefile.
Definition helium_io.F:1661
helium_methods
Methods dealing with helium_solvent_type.
Definition helium_methods.F:13
helium_methods::helium_release
subroutine, public helium_release(helium_env)
Releases helium_solvent_type.
Definition helium_methods.F:856
helium_methods::helium_create
subroutine, public helium_create(helium_env, input, solute)
Data-structure that holds all needed information about (superfluid) helium solvent.
Definition helium_methods.F:113
helium_methods::helium_init
subroutine, public helium_init(helium_env, pint_env)
Initialize helium data structures.
Definition helium_methods.F:1056
helium_sampling
Methods for sampling helium variables.
Definition helium_sampling.F:13
helium_sampling::helium_step
subroutine, public helium_step(helium_env, pint_env)
Perform MC step for helium.
Definition helium_sampling.F:420
helium_sampling::helium_do_run
subroutine, public helium_do_run(helium_env, globenv)
Performs MC simulation for helium (only)
Definition helium_sampling.F:81
helium_types
Data types representing superfluid helium.
Definition helium_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::propagator_cmd
integer, parameter, public propagator_cmd
Definition input_constants.F:314
input_constants::propagator_rpmd
integer, parameter, public propagator_rpmd
Definition input_constants.F:314
input_constants::integrate_exact
integer, parameter, public integrate_exact
Definition input_constants.F:319
input_constants::transformation_stage
integer, parameter, public transformation_stage
Definition input_constants.F:310
input_constants::integrate_numeric
integer, parameter, public integrate_numeric
Definition input_constants.F:319
input_constants::transformation_normal
integer, parameter, public transformation_normal
Definition input_constants.F:310
input_constants::propagator_bcmd
integer, parameter, public propagator_bcmd
Definition input_constants.F:314
input_cp2k_restarts
Set of routines to dump the restart file of CP2K.
Definition input_cp2k_restarts.F:13
input_cp2k_restarts::write_restart
subroutine, public write_restart(md_env, force_env, root_section, coords, vels, pint_env, helium_env)
checks if a restart needs to be written and does so, updating all necessary fields in the input file....
Definition input_cp2k_restarts.F:112
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_unset
subroutine, public section_vals_val_unset(section_vals, keyword_name, i_rep_section, i_rep_val)
unsets (removes) the requested value (if it is a keyword repetitions removes the repetition,...
Definition input_section_types.F:1354
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1223
input_section_types::section_vals_retain
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
Definition input_section_types.F:638
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
input_section_types::section_vals_add_values
subroutine, public section_vals_add_values(section_vals)
adds the place to store the values of a repetition of the section
Definition input_section_types.F:848
input_section_types::section_vals_release
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Definition input_section_types.F:653
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::gcd
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1287
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_comm_self
type(mp_comm_type), parameter, public mp_comm_self
Definition message_passing.F:792
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_list_types
represent a simple array based list of the given type
Definition molecule_list_types.F:15
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition parallel_rng_types.F:51
parallel_rng_types::gaussian
integer, parameter, public gaussian
Definition parallel_rng_types.F:73
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pint_gle
Methods to apply GLE to PI runs.
Definition pint_gle.F:16
pint_gle::pint_gle_init
subroutine, public pint_gle_init(pint_env)
...
Definition pint_gle.F:54
pint_gle::pint_gle_step
subroutine, public pint_gle_step(pint_env)
...
Definition pint_gle.F:78
pint_gle::pint_calc_gle_energy
elemental subroutine, public pint_calc_gle_energy(pint_env)
...
Definition pint_gle.F:37
pint_io
I/O subroutines for pint_env.
Definition pint_io.F:13
pint_io::pint_write_action
subroutine, public pint_write_action(pint_env)
Writes out the actions according to PINTPRINTACTION.
Definition pint_io.F:557
pint_io::pint_write_centroids
subroutine, public pint_write_centroids(pint_env)
Write out the trajectory of the centroid (positions and velocities)
Definition pint_io.F:104
pint_io::pint_write_rgyr
subroutine, public pint_write_rgyr(pint_env)
Write radii of gyration according to PINTPRINTCENTROID_GYR.
Definition pint_io.F:675
pint_io::pint_write_step_info
subroutine, public pint_write_step_info(pint_env)
Write step info to the output file.
Definition pint_io.F:620
pint_io::pint_write_ener
subroutine, public pint_write_ener(pint_env)
Writes out the energies according to PINTPRINTENERGY.
Definition pint_io.F:481
pint_io::pint_write_line
subroutine, public pint_write_line(line)
Writes out a line of text to the default output unit.
Definition pint_io.F:76
pint_io::pint_write_trajectory
subroutine, public pint_write_trajectory(pint_env)
Write out the trajectory of the beads (positions and velocities)
Definition pint_io.F:259
pint_io::pint_write_com
subroutine, public pint_write_com(pint_env)
Write center of mass (COM) position according to PINTPRINTCOM.
Definition pint_io.F:415
pint_methods
Methods to performs a path integral run.
Definition pint_methods.F:18
pint_methods::do_pint_run
subroutine, public do_pint_run(para_env, input, input_declaration, globenv)
Perform a path integral simulation.
Definition pint_methods.F:907
pint_normalmode
Data type and methods dealing with PI calcs in normal mode coords.
Definition pint_normalmode.F:18
pint_normalmode::normalmode_calc_uf_h
pure subroutine, public normalmode_calc_uf_h(normalmode_env, mass_beads, ux, uf_h, e_h)
calculates the harmonic force in the normal mode basis
Definition pint_normalmode.F:307
pint_normalmode::normalmode_release
pure subroutine, public normalmode_release(normalmode_env)
releases the normalmode environment
Definition pint_normalmode.F:190
pint_normalmode::normalmode_env_create
subroutine, public normalmode_env_create(normalmode_env, normalmode_section, p, kt, propagator)
creates the data needed for a normal mode transformation
Definition pint_normalmode.F:57
pint_normalmode::normalmode_init_masses
pure subroutine, public normalmode_init_masses(normalmode_env, mass, mass_beads, mass_fict, q)
initializes the masses and fictitious masses compatible with the normal mode information
Definition pint_normalmode.F:211
pint_piglet
Methods to apply the piglet thermostat to PI runs.
Definition pint_piglet.F:14
pint_piglet::pint_calc_piglet_energy
elemental subroutine, public pint_calc_piglet_energy(pint_env)
returns the piglet kinetic energy contribution
Definition pint_piglet.F:673
pint_piglet::pint_piglet_release
subroutine, public pint_piglet_release(piglet_therm)
releases the piglet environment
Definition pint_piglet.F:509
pint_piglet::pint_piglet_create
subroutine, public pint_piglet_create(piglet_therm, pint_env, section)
creates the data structure for a piglet thermostating in PI runs
Definition pint_piglet.F:92
pint_piglet::pint_piglet_step
subroutine, public pint_piglet_step(vold, vnew, first_mode, masses, piglet_therm)
...
Definition pint_piglet.F:411
pint_piglet::pint_piglet_init
subroutine, public pint_piglet_init(piglet_therm, pint_env, section, dt, para_env)
initializes the data for a piglet run
Definition pint_piglet.F:142
pint_pile
Methods to apply a simple Lagevin thermostat to PI runs. v_new = c1*vold + SQRT(kT/m)*c2*random.
Definition pint_pile.F:15
pint_pile::pint_pile_step
subroutine, public pint_pile_step(vold, vnew, p, ndim, first_mode, masses, pile_therm)
...
Definition pint_pile.F:149
pint_pile::pint_pile_init
subroutine, public pint_pile_init(pile_therm, pint_env, normalmode_env, section)
initializes the data for a pile run
Definition pint_pile.F:54
pint_pile::pint_pile_release
subroutine, public pint_pile_release(pile_therm)
releases the pile environment
Definition pint_pile.F:182
pint_pile::pint_calc_pile_energy
subroutine, public pint_calc_pile_energy(pint_env)
returns the pile kinetic energy contribution
Definition pint_pile.F:198
pint_public
Public path integral routines that can be called from other modules.
Definition pint_public.F:15
pint_public::pint_levy_walk
subroutine, public pint_levy_walk(x0, n, v, x, rng_gaussian)
Perform a Brownian walk of length n around x0 with the variance v.
Definition pint_public.F:204
pint_qtb
Methods to apply the QTB thermostat to PI runs. Based on the PILE implementation from Felix Uhl (pint...
Definition pint_qtb.F:15
pint_qtb::pint_qtb_step
subroutine, public pint_qtb_step(vold, vnew, p, ndim, masses, qtb_therm)
...
Definition pint_qtb.F:155
pint_qtb::pint_calc_qtb_energy
subroutine, public pint_calc_qtb_energy(pint_env)
returns the qtb kinetic energy contribution
Definition pint_qtb.F:240
pint_qtb::pint_qtb_release
subroutine, public pint_qtb_release(qtb_therm)
releases the qtb environment
Definition pint_qtb.F:219
pint_qtb::pint_qtb_init
subroutine, public pint_qtb_init(qtb_therm, pint_env, normalmode_env, section)
initializes the data for a QTB run
Definition pint_qtb.F:69
pint_staging
Data type and methods dealing with PI calcs in staging coordinates.
Definition pint_staging.F:16
pint_staging::staging_release
elemental subroutine, public staging_release(staging_env)
releases the staging environment, kept for symmetry reasons with staging_env_create
Definition pint_staging.F:69
pint_staging::staging_env_create
subroutine, public staging_env_create(staging_env, staging_section, p, kt)
creates the data needed for a staging transformation
Definition pint_staging.F:45
pint_staging::staging_calc_uf_h
pure subroutine, public staging_calc_uf_h(staging_env, mass_beads, ux, uf_h, e_h)
calculates the harmonic force in the staging basis
Definition pint_staging.F:260
pint_staging::staging_init_masses
pure subroutine, public staging_init_masses(staging_env, mass, mass_beads, mass_fict, q)
initializes the masses and fictitious masses compatibly with the staging information
Definition pint_staging.F:88
pint_transformations
Definition pint_transformations.F:8
pint_transformations::pint_x2u
subroutine, public pint_x2u(pint_env, ux, x)
Transforms from the x into the u variables (at the moment a staging transformation for the positions)
Definition pint_transformations.F:43
pint_transformations::pint_u2x
subroutine, public pint_u2x(pint_env, ux, x)
transform from the u variable to the x (inverse of x2u)
Definition pint_transformations.F:75
pint_transformations::pint_f2uf
subroutine, public pint_f2uf(pint_env, uf, f)
transformation x to u for the forces
Definition pint_transformations.F:110
pint_types
Definition pint_types.F:8
pint_types::e_kin_thermo_id
integer, parameter, public e_kin_thermo_id
Definition pint_types.F:25
pint_types::e_conserved_id
integer, parameter, public e_conserved_id
Definition pint_types.F:25
pint_types::thermostat_none
integer, parameter, public thermostat_none
Definition pint_types.F:33
pint_types::thermostat_gle
integer, parameter, public thermostat_gle
Definition pint_types.F:33
pint_types::e_potential_id
integer, parameter, public e_potential_id
Definition pint_types.F:25
pint_types::thermostat_pile
integer, parameter, public thermostat_pile
Definition pint_types.F:33
pint_types::thermostat_piglet
integer, parameter, public thermostat_piglet
Definition pint_types.F:33
pint_types::thermostat_nose
integer, parameter, public thermostat_nose
Definition pint_types.F:33
pint_types::e_kin_virial_id
integer, parameter, public e_kin_virial_id
Definition pint_types.F:25
pint_types::thermostat_qtb
integer, parameter, public thermostat_qtb
Definition pint_types.F:33
replica_methods
methods to setup replicas of the same system differing only by atom positions and velocities (as used...
Definition replica_methods.F:16
replica_methods::rep_env_create
subroutine, public rep_env_create(rep_env, para_env, input, input_declaration, nrep, prep, sync_v, keep_wf_history, row_force)
creates a replica environment together with its force environment
Definition replica_methods.F:91
replica_methods::rep_env_calc_e_f
subroutine, public rep_env_calc_e_f(rep_env, calc_f)
evaluates the forces
Definition replica_methods.F:366
replica_types
types used to handle many replica of the same system that differ only in atom positions,...
Definition replica_types.F:21
replica_types::rep_env_release
subroutine, public rep_env_release(rep_env)
releases the given replica environment
Definition replica_types.F:138
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
simpar_types::release_simpar_type
subroutine, public release_simpar_type(simpar)
Releases the simulation parameters type.
Definition simpar_types.F:118
simpar_types::create_simpar_type
subroutine, public create_simpar_type(simpar)
Creates the simulation parameters type.
Definition simpar_types.F:106
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
f77_interface::f_env_type
Definition f77_interface.F:169
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
helium_types::helium_solvent_p_type
data structure for array of solvent helium environments
Definition helium_types.F:300
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:151
molecule_list_types::molecule_list_type
represent a list of objects
Definition molecule_list_types.F:46
molecule_types::global_constraint_type
Definition molecule_types.F:134
molecule_types::molecule_type
Definition molecule_types.F:159
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pint_types::pint_env_type
environment for a path integral run
Definition pint_types.F:112
replica_types::replica_env_type
keeps replicated information about the replicas
Definition replica_types.F:96