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