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