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