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