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