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