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