Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief calculation section for TreeMonteCarlo
10 : !> \par History
11 : !> 11.2012 created [Mandes Schoenherr]
12 : !> \author Mandes
13 : ! **************************************************************************************************
14 :
15 : MODULE tmc_calculations
16 : USE cell_methods, ONLY: init_cell
17 : USE cell_types, ONLY: cell_copy,&
18 : cell_type,&
19 : get_cell,&
20 : pbc
21 : USE cp_log_handling, ONLY: cp_to_string
22 : USE f77_interface, ONLY: calc_energy,&
23 : calc_force,&
24 : set_cell
25 : USE kinds, ONLY: dp
26 : USE mathconstants, ONLY: pi
27 : USE parallel_rng_types, ONLY: rng_stream_type
28 : USE physcon, ONLY: boltzmann,&
29 : joule
30 : USE tmc_move_types, ONLY: mv_type_MD
31 : USE tmc_stati, ONLY: task_type_MC,&
32 : task_type_gaussian_adaptation,&
33 : task_type_ideal_gas
34 : USE tmc_tree_types, ONLY: tree_type
35 : USE tmc_types, ONLY: tmc_atom_type,&
36 : tmc_env_type,&
37 : tmc_param_type
38 : #include "../base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_calculations'
45 :
46 : PUBLIC :: calc_potential_energy
47 : PUBLIC :: get_scaled_cell, get_cell_scaling
48 : PUBLIC :: nearest_distance
49 : PUBLIC :: geometrical_center, center_of_mass
50 : PUBLIC :: init_vel, calc_e_kin
51 : PUBLIC :: compute_estimated_prob
52 : PUBLIC :: get_subtree_efficiency
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief start the calculation of the energy
57 : !> (distinguish between exact and approximate)
58 : !> \param conf actual configurations to calculate potential energy
59 : !> \param env_id f77_interface env id
60 : !> \param exact_approx_pot flag if result should be stores in exact or approx
61 : !> energy variable
62 : !> \param tmc_env TMC environment parameters
63 : !> \author Mandes 01.2013
64 : ! **************************************************************************************************
65 4521 : SUBROUTINE calc_potential_energy(conf, env_id, exact_approx_pot, &
66 : tmc_env)
67 : TYPE(tree_type), POINTER :: conf
68 : INTEGER :: env_id
69 : LOGICAL :: exact_approx_pot
70 : TYPE(tmc_env_type), POINTER :: tmc_env
71 :
72 : INTEGER :: ierr
73 : LOGICAL :: flag
74 : REAL(KIND=dp) :: e_pot, rnd
75 : TYPE(cell_type), POINTER :: tmp_cell
76 :
77 4521 : rnd = 0.0_dp
78 :
79 4521 : CPASSERT(ASSOCIATED(conf))
80 4521 : CPASSERT(env_id > 0)
81 4521 : CPASSERT(ASSOCIATED(tmc_env))
82 :
83 9042 : SELECT CASE (tmc_env%params%task_type)
84 : CASE (task_type_gaussian_adaptation)
85 : !CALL gaussian_adaptation_energy(, )
86 : CASE (task_type_MC)
87 4521 : IF (tmc_env%params%pressure >= 0.0_dp) THEN
88 53682 : ALLOCATE (tmp_cell)
89 : CALL get_scaled_cell(cell=tmc_env%params%cell, box_scale=conf%box_scale, &
90 778 : scaled_cell=tmp_cell)
91 778 : CALL set_cell(env_id=env_id, new_cell=tmp_cell%hmat, ierr=ierr)
92 778 : CPASSERT(ierr == 0)
93 778 : DEALLOCATE (tmp_cell)
94 : END IF
95 :
96 : ! TODO check for minimal distances
97 4521 : flag = .TRUE.
98 0 : IF (flag .EQV. .TRUE.) THEN
99 4521 : IF (tmc_env%params%print_forces .OR. &
100 : conf%move_type == mv_type_MD) THEN
101 : e_pot = 0.0_dp
102 38272 : conf%frc(:) = 0.0_dp
103 : CALL calc_force(env_id=env_id, pos=conf%pos, n_el_pos=SIZE(conf%pos), &
104 : e_pot=e_pot, force=conf%frc, &
105 598 : n_el_force=SIZE(conf%frc), ierr=ierr)
106 : ELSE
107 : e_pot = 0.0_dp
108 3923 : CALL calc_energy(env_id=env_id, pos=conf%pos, n_el=SIZE(conf%pos), e_pot=e_pot, ierr=ierr)
109 : END IF
110 : ELSE
111 : e_pot = HUGE(e_pot)
112 : END IF
113 : CASE (task_type_ideal_gas)
114 0 : e_pot = 0.0_dp
115 : CASE DEFAULT
116 : CALL cp_abort(__LOCATION__, &
117 : "worker task typ is unknown "// &
118 4521 : cp_to_string(tmc_env%params%task_type))
119 : END SELECT
120 :
121 : ! --- wait a bit
122 4521 : rnd = tmc_env%rng_stream%next()
123 : !rnd = 0.5
124 : !TODO IF(worker_random_wait.AND.exact_approx_pot)THEN
125 : ! CALL SYSTEM_CLOCK(time0, time_rate, time_max)
126 : ! wait_end=time0+(1.0+rnd)*worker_wait_msec*time_rate/1000.0
127 : ! !wait_end=time0+((worker_wait_msec*time_rate+999)/1000)
128 : ! time_wait: DO
129 : ! CALL SYSTEM_CLOCK(time1, time_rate, time_max)
130 : ! IF(time1<time0.OR.time1>wait_end) exit time_wait
131 : ! END DO time_wait
132 : ! END IF
133 4521 : IF (exact_approx_pot) THEN
134 4324 : conf%potential = e_pot
135 : ELSE
136 197 : conf%e_pot_approx = e_pot
137 : END IF
138 4521 : END SUBROUTINE calc_potential_energy
139 :
140 : ! **************************************************************************************************
141 : !> \brief handles properties and calculations of a scaled cell
142 : !> \param cell original cell
143 : !> \param box_scale scaling factors for each direction
144 : !> \param scaled_hmat returns the scaled h matrix (matrix of cell vectors)
145 : !> \param scaled_cell ...
146 : !> \param vol returns the cell volume
147 : !> \param abc ...
148 : !> \param vec a vector, which will be folded (pbc) in the cell
149 : !> \author Mandes 11.2012
150 : ! **************************************************************************************************
151 218382 : SUBROUTINE get_scaled_cell(cell, box_scale, scaled_hmat, scaled_cell, vol, &
152 : abc, vec)
153 : TYPE(cell_type), INTENT(IN), POINTER :: cell
154 : REAL(KIND=dp), DIMENSION(:), POINTER :: box_scale
155 : REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: scaled_hmat
156 : TYPE(cell_type), OPTIONAL, POINTER :: scaled_cell
157 : REAL(KIND=dp), OPTIONAL :: vol
158 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: abc
159 : REAL(KIND=dp), DIMENSION(3), OPTIONAL :: vec
160 :
161 : LOGICAL :: new_scaled_cell
162 : TYPE(cell_type), POINTER :: tmp_cell
163 :
164 218382 : CPASSERT(ASSOCIATED(cell))
165 218382 : CPASSERT(ASSOCIATED(box_scale))
166 :
167 218382 : new_scaled_cell = .FALSE.
168 :
169 218382 : IF (.NOT. PRESENT(scaled_cell)) THEN
170 15013848 : ALLOCATE (tmp_cell)
171 217592 : new_scaled_cell = .TRUE.
172 : ELSE
173 790 : tmp_cell => scaled_cell
174 : END IF
175 218382 : CALL cell_copy(cell_in=cell, cell_out=tmp_cell)
176 873528 : tmp_cell%hmat(:, 1) = tmp_cell%hmat(:, 1)*box_scale(1)
177 873528 : tmp_cell%hmat(:, 2) = tmp_cell%hmat(:, 2)*box_scale(2)
178 873528 : tmp_cell%hmat(:, 3) = tmp_cell%hmat(:, 3)*box_scale(3)
179 218382 : CALL init_cell(cell=tmp_cell)
180 :
181 218382 : IF (PRESENT(scaled_hmat)) THEN
182 5096 : scaled_hmat(:, :) = tmp_cell%hmat
183 : END IF
184 :
185 218382 : IF (PRESENT(vec)) THEN
186 862448 : vec = pbc(r=vec, cell=tmp_cell)
187 : END IF
188 :
189 218382 : IF (PRESENT(vol)) CALL get_cell(cell=tmp_cell, deth=vol)
190 218382 : IF (PRESENT(abc)) CALL get_cell(cell=tmp_cell, abc=abc)
191 218382 : IF (new_scaled_cell) DEALLOCATE (tmp_cell)
192 :
193 218382 : END SUBROUTINE get_scaled_cell
194 :
195 : ! **************************************************************************************************
196 : !> \brief handles properties and calculations of a scaled cell
197 : !> \param cell original cell
198 : !> \param scaled_hmat returns the scaled h matrix (matrix of cell vectors)
199 : !> \param box_scale scaling factors for each direction
200 : !> \author Mandes 11.2012
201 : ! **************************************************************************************************
202 1206 : SUBROUTINE get_cell_scaling(cell, scaled_hmat, box_scale)
203 : TYPE(cell_type), INTENT(IN), POINTER :: cell
204 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: scaled_hmat
205 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: box_scale
206 :
207 : REAL(KIND=dp), DIMENSION(3) :: abc_new, abc_orig
208 : TYPE(cell_type), POINTER :: tmp_cell
209 :
210 1206 : CPASSERT(ASSOCIATED(cell))
211 :
212 83214 : ALLOCATE (tmp_cell)
213 1206 : CALL cell_copy(cell_in=cell, cell_out=tmp_cell)
214 15678 : tmp_cell%hmat(:, :) = scaled_hmat(:, :)
215 1206 : CALL init_cell(cell=tmp_cell)
216 1206 : CALL get_cell(cell=cell, abc=abc_orig)
217 1206 : CALL get_cell(cell=tmp_cell, abc=abc_new)
218 :
219 4824 : box_scale(:) = abc_new(:)/abc_orig(:)
220 :
221 1206 : DEALLOCATE (tmp_cell)
222 1206 : END SUBROUTINE get_cell_scaling
223 :
224 : ! **************************************************************************************************
225 : !> \brief neares distance of atoms within the periodic boundary condition
226 : !> \param x1 ...
227 : !> \param x2 ...
228 : !> \param cell ...
229 : !> \param box_scale ...
230 : !> \return ...
231 : !> \author Mandes 11.2012
232 : ! **************************************************************************************************
233 187940 : FUNCTION nearest_distance(x1, x2, cell, box_scale) RESULT(res)
234 : REAL(KIND=dp), DIMENSION(:) :: x1, x2
235 : TYPE(cell_type), POINTER :: cell
236 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: box_scale
237 : REAL(KIND=dp) :: res
238 :
239 : REAL(KIND=dp), DIMENSION(3) :: dist_vec
240 187940 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_box_scale
241 :
242 187940 : NULLIFY (tmp_box_scale)
243 :
244 0 : CPASSERT(ASSOCIATED(cell))
245 187940 : CPASSERT(SIZE(x1) == 3)
246 187940 : CPASSERT(SIZE(x2) == 3)
247 :
248 751760 : dist_vec(:) = x2(:) - x1(:) ! distance vector between atoms
249 187940 : ALLOCATE (tmp_box_scale(3))
250 187940 : IF (PRESENT(box_scale)) THEN
251 187940 : CPASSERT(SIZE(box_scale) == 3)
252 1315580 : tmp_box_scale(:) = box_scale
253 : ELSE
254 0 : tmp_box_scale(:) = 1.0_dp
255 : END IF
256 187940 : CALL get_scaled_cell(cell=cell, box_scale=box_scale, vec=dist_vec)
257 751760 : res = SQRT(SUM(dist_vec(:)*dist_vec(:)))
258 187940 : DEALLOCATE (tmp_box_scale)
259 187940 : END FUNCTION nearest_distance
260 :
261 : ! **************************************************************************************************
262 : !> \brief calculate the geometrical center of an amount of atoms
263 : !> array size should be multiple of dim_per_elem
264 : !> \param pos list of atoms
265 : !> \param center return value, the geometrical center
266 : !> \author Mandes 11.2012
267 : ! **************************************************************************************************
268 7760 : SUBROUTINE geometrical_center(pos, center)
269 : REAL(KIND=dp), DIMENSION(:) :: pos
270 : REAL(KIND=dp), DIMENSION(:), POINTER :: center
271 :
272 : CHARACTER(LEN=*), PARAMETER :: routineN = 'geometrical_center'
273 :
274 : INTEGER :: handle, i
275 :
276 7760 : CPASSERT(ASSOCIATED(center))
277 7760 : CPASSERT(SIZE(pos) >= SIZE(center))
278 :
279 : ! start the timing
280 7760 : CALL timeset(routineN, handle)
281 :
282 31040 : center = 0.0_dp
283 31138 : DO i = 1, SIZE(pos), SIZE(center)
284 : center(:) = center(:) + &
285 101272 : pos(i:i + SIZE(center) - 1)/(SIZE(pos)/REAL(SIZE(center), KIND=dp))
286 : END DO
287 : ! end the timing
288 7760 : CALL timestop(handle)
289 7760 : END SUBROUTINE geometrical_center
290 :
291 : ! **************************************************************************************************
292 : !> \brief calculate the center of mass of an amount of atoms
293 : !> array size should be multiple of dim_per_elem
294 : !> \param pos ...
295 : !> \param atoms ...
296 : !> \param center ...
297 : !> \param
298 : !> \param
299 : !> \author Mandes 11.2012
300 : ! **************************************************************************************************
301 0 : SUBROUTINE center_of_mass(pos, atoms, center)
302 : REAL(KIND=dp), DIMENSION(:) :: pos
303 : TYPE(tmc_atom_type), DIMENSION(:), OPTIONAL :: atoms
304 : REAL(KIND=dp), DIMENSION(:), POINTER :: center
305 :
306 : CHARACTER(LEN=*), PARAMETER :: routineN = 'center_of_mass'
307 :
308 : INTEGER :: handle, i
309 : REAL(KIND=dp) :: mass_sum, mass_tmp
310 :
311 0 : CPASSERT(ASSOCIATED(center))
312 0 : CPASSERT(SIZE(pos) >= SIZE(center))
313 :
314 : ! start the timing
315 0 : CALL timeset(routineN, handle)
316 :
317 0 : center = 0.0_dp
318 0 : mass_sum = 0.0_dp
319 0 : DO i = 1, SIZE(pos), SIZE(center)
320 0 : IF (PRESENT(atoms)) THEN
321 0 : CPASSERT(SIZE(atoms) == SIZE(pos)/SIZE(center))
322 0 : mass_tmp = atoms(INT(i/REAL(SIZE(center), KIND=dp)) + 1)%mass
323 : center(:) = center(:) + pos(i:i + SIZE(center) - 1)/ &
324 0 : (SIZE(pos)/REAL(SIZE(center), KIND=dp))*mass_tmp
325 0 : mass_sum = mass_sum + mass_tmp
326 : ELSE
327 0 : CPWARN("try to calculate center of mass without any mass.")
328 : center(:) = center(:) + pos(i:i + SIZE(center) - 1)/ &
329 0 : (SIZE(pos)/REAL(SIZE(center), KIND=dp))
330 0 : mass_sum = 1.0_dp
331 : END IF
332 : END DO
333 0 : center(:) = center(:)/mass_sum
334 : ! end the timing
335 0 : CALL timestop(handle)
336 0 : END SUBROUTINE center_of_mass
337 :
338 : ! **************************************************************************************************
339 : !> \brief routine sets initial velocity, using the Box-Muller Method for Normal
340 : !> (Gaussian) Deviates
341 : !> \param vel ...
342 : !> \param atoms ...
343 : !> \param temerature ...
344 : !> \param rng_stream ...
345 : !> \param rnd_seed ...
346 : !> \author Mandes 11.2012
347 : ! **************************************************************************************************
348 0 : SUBROUTINE init_vel(vel, atoms, temerature, rng_stream, rnd_seed)
349 : REAL(KIND=dp), DIMENSION(:), POINTER :: vel
350 : TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
351 : REAL(KIND=dp) :: temerature
352 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
353 : REAL(KIND=dp), DIMENSION(3, 2, 3) :: rnd_seed
354 :
355 : INTEGER :: i
356 : REAL(KIND=dp) :: kB, mass_tmp, rnd1, rnd2
357 :
358 0 : kB = boltzmann/joule
359 :
360 0 : CPASSERT(ASSOCIATED(vel))
361 0 : CPASSERT(ASSOCIATED(atoms))
362 :
363 0 : CALL rng_stream%set(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
364 0 : DO i = 1, SIZE(vel)
365 0 : rnd1 = rng_stream%next()
366 0 : rnd2 = rng_stream%next()
367 :
368 0 : mass_tmp = atoms(INT(i/REAL(3, KIND=dp)) + 1)%mass
369 :
370 : vel(i) = SQRT(-2.0_dp*LOG(rnd1))*COS(2.0_dp*PI*rnd2)* &
371 0 : SQRT(kB*temerature/mass_tmp)
372 : END DO
373 0 : CALL rng_stream%get(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
374 :
375 0 : END SUBROUTINE init_vel
376 :
377 : ! **************************************************************************************************
378 : !> \brief routine calculates the kinetic energy, using the velocities
379 : !> and atom mass, both in atomic units
380 : !> \param vel ...
381 : !> \param atoms ...
382 : !> \return ...
383 : !> \author Mandes 11.2012
384 : ! **************************************************************************************************
385 0 : FUNCTION calc_e_kin(vel, atoms) RESULT(ekin)
386 : REAL(KIND=dp), DIMENSION(:), POINTER :: vel
387 : TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
388 : REAL(KIND=dp) :: ekin
389 :
390 : INTEGER :: i
391 : REAL(KIND=dp) :: mass_tmp
392 :
393 0 : CPASSERT(ASSOCIATED(vel))
394 0 : CPASSERT(ASSOCIATED(atoms))
395 0 : ekin = 0.0_dp
396 :
397 0 : DO i = 1, SIZE(vel)
398 0 : mass_tmp = atoms(INT(i/REAL(3, KIND=dp)) + 1)%mass
399 0 : ekin = ekin + 0.5_dp*mass_tmp*vel(i)*vel(i)
400 : END DO
401 0 : END FUNCTION calc_e_kin
402 :
403 : ! **************************************************************************************************
404 : !> \brief assuming an (exponential) decreasing function, this function
405 : !> extrapolate the converged value
406 : !> \param v1 function values
407 : !> \param v2 function values
408 : !> \param v3 function values
409 : !> \param extrapolate extrapolated final value (result)
410 : !> \param res_err error of the result
411 : !> \author Mandes 12.2012
412 : ! **************************************************************************************************
413 0 : SUBROUTINE three_point_extrapolate(v1, v2, v3, extrapolate, res_err)
414 : REAL(KIND=dp) :: v1, v2, v3
415 : REAL(KIND=dp), INTENT(OUT) :: extrapolate, res_err
416 :
417 : REAL(KIND=dp) :: e1, e2, e3
418 : REAL(KIND=dp) :: a, b, c, d12, d23, ddd
419 :
420 : extrapolate = HUGE(extrapolate)
421 :
422 : !> solve({exp(a+b)+c = e1, exp(2*a+b)+c = e2, exp(3*a+b)+c = e3}, [a, b, c])
423 : !> solve({a*b+c = e1, a^2*b+c = e2, a^3*b+c = e3}, [a, b, c]);
424 : ! [[ 3 2 ]]
425 : ! [[ -e3 + e2 (e1 - e2) -e2 + e1 e3 ]]
426 : ! [[a = --------, b = ---------------------------, c = --------------]]
427 : ! [[ e1 - e2 (-e3 + e2) (e3 - 2 e2 + e1) e3 - 2 e2 + e1]]
428 :
429 : ! sort so that e1>=e2>=e3
430 0 : e1 = v1; e2 = v2; e3 = v3
431 0 : CALL swap(e1, e2)
432 0 : CALL swap(e1, e3)
433 0 : CALL swap(e2, e3)
434 : ! we need extra care if some of the difference e1-e2, e3-e2 are nearly zero,
435 : ! since the formulae suffer from sever loss of precision
436 0 : d12 = e1 - e2
437 0 : d23 = e2 - e3
438 0 : ddd = d12 - d23
439 0 : IF (d12 == 0 .OR. d23 == 0 .OR. ABS(ddd) == 0) THEN
440 : ! a degenerate case, we do no extrapolation
441 0 : extrapolate = e3
442 0 : res_err = e1 - e3
443 : ELSE
444 0 : a = d23/d12
445 0 : b = (d12**3/(d23*ddd))
446 0 : c = e2 - (d12*d23)/ddd
447 : ! extrapolation, let's only look 4 iterations ahead, more is presumably anyway not accurate
448 : ! fewer is maybe more stable
449 0 : extrapolate = a**7*b + c
450 0 : res_err = e3 - extrapolate
451 : END IF
452 0 : CPASSERT(extrapolate /= HUGE(extrapolate))
453 : CONTAINS
454 : ! **************************************************************************************************
455 : !> \brief ...
456 : !> \param x1 ...
457 : !> \param x2 ...
458 : ! **************************************************************************************************
459 0 : SUBROUTINE swap(x1, x2)
460 : REAL(KIND=dp) :: x1, x2
461 :
462 : REAL(KIND=dp) :: tmp
463 :
464 0 : IF (x2 > x1) THEN
465 0 : tmp = x2
466 0 : x2 = x1
467 0 : x1 = tmp
468 : END IF
469 0 : END SUBROUTINE swap
470 : END SUBROUTINE three_point_extrapolate
471 :
472 : ! **************************************************************************************************
473 : !> \brief calculates the probability of acceptance for given intervals of the
474 : !> exact energy
475 : !> \param E_n_mu energy distribution of new configuration
476 : !> \param E_n_sigma energy distribution of new configuration
477 : !> \param E_o_mu energy distribution of old configuration
478 : !> \param E_o_sigma energy distribution of old configuration
479 : !> \param E_classical_diff the difference in approximated energies for the
480 : !> old and new configuration (E_o-E_n)
481 : !> \param prior_mu energy distribution of the already converged
482 : !> energies
483 : !> \param prior_sigma energy distribution of the already converged
484 : !> energies
485 : !> \param p the random number, the criteria has to be smaller than this
486 : !> \param beta ...
487 : !> \return return probability of acceptance
488 : !> \author Mandes 12.2012
489 : ! **************************************************************************************************
490 0 : FUNCTION compute_prob(E_n_mu, E_n_sigma, E_o_mu, E_o_sigma, E_classical_diff, &
491 : prior_mu, prior_sigma, p, beta) RESULT(prob)
492 : REAL(KIND=dp) :: E_n_mu, E_n_sigma, E_o_mu, E_o_sigma, &
493 : E_classical_diff, prior_mu, &
494 : prior_sigma, p, beta, prob
495 :
496 : ! INTEGER :: io,in
497 : ! REAL(KIND=dp) :: diff,E_n,E_o,surface,lower_bound,upper_bound,delta
498 :
499 : prob = 0.5_dp*ERFC(-0.5_dp*SQRT(2.0_dp)*( &
500 : (-prior_sigma**2 - E_o_sigma**2 - E_n_sigma**2)*LOG(p) + &
501 : ((E_classical_diff - E_n_mu + E_o_mu)*prior_sigma**2 - prior_mu*(E_n_sigma**2 + E_o_sigma**2))*beta)/ &
502 0 : (SQRT(E_o_sigma**2 + E_n_sigma**2)*SQRT(prior_sigma**2 + E_o_sigma**2 + E_n_sigma**2)*prior_sigma*beta))
503 :
504 0 : prob = MIN(1.0_dp - EPSILON(1.0_dp), MAX(EPSILON(1.0_dp), prob))
505 :
506 0 : END FUNCTION compute_prob
507 :
508 : ! **************************************************************************************************
509 : !> \brief extimates the probability of acceptance considering the intermetiate
510 : !> step energies
511 : !> \param elem_old old/parent sub tree element
512 : !> \param elem_new new/actual sub tree element, which schould be checked
513 : !> \param E_classical_diff difference in the classical energy of the old and
514 : !> new configuration
515 : !> \param rnd_nr random number acceptance check will be done with
516 : !> \param beta 1/(kB*T) can differ for different acceptance checks
517 : !> \param tmc_params TMC environment parameters
518 : !> \return estimated acceptance probability
519 : !> \author Mandes 12.2012
520 : ! **************************************************************************************************
521 0 : FUNCTION compute_estimated_prob(elem_old, elem_new, E_classical_diff, &
522 : rnd_nr, beta, tmc_params) RESULT(prob)
523 : TYPE(tree_type), POINTER :: elem_old, elem_new
524 : REAL(KIND=dp) :: E_classical_diff, rnd_nr, beta
525 : TYPE(tmc_param_type), POINTER :: tmc_params
526 : REAL(KIND=dp) :: prob
527 :
528 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_estimated_prob'
529 :
530 : INTEGER :: handle
531 : REAL(KIND=dp) :: E_mu_tmp, E_n_mu, E_n_sigma, E_o_mu, &
532 : E_o_sigma, E_sigma_tmp, prior_sigma
533 :
534 0 : CPASSERT(ASSOCIATED(elem_old))
535 0 : CPASSERT(ASSOCIATED(elem_new))
536 0 : CPASSERT(rnd_nr > 0.0_dp)
537 :
538 : ! start the timing
539 0 : CALL timeset(routineN, handle)
540 :
541 0 : prob = -1.0_dp
542 : IF ((elem_new%scf_energies_count >= 3) .AND. &
543 0 : (elem_old%scf_energies_count >= 3) .AND. &
544 : tmc_params%prior_NMC_acc%counter >= 10) THEN
545 : !-- first the new element energy estimation
546 : ! using 3 point extrapolation of two different intervals -> more stable estimation
547 : ! the energies are sorted in the three_point_extrapolate routine !
548 : ! But with array of length 4 we have to select the 3 connected ones
549 : CALL three_point_extrapolate(v1=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 3, 4) + 1), &
550 : v2=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 2, 4) + 1), &
551 : v3=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 1, 4) + 1), &
552 0 : extrapolate=E_mu_tmp, res_err=E_sigma_tmp)
553 0 : IF ((elem_new%scf_energies_count > 3)) THEN
554 : CALL three_point_extrapolate(v1=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 4, 4) + 1), &
555 : v2=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 3, 4) + 1), &
556 : v3=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 2, 4) + 1), &
557 0 : extrapolate=E_n_mu, res_err=E_n_sigma)
558 0 : E_n_sigma = MAX(E_n_sigma, ABS(E_n_mu - E_mu_tmp))
559 : ELSE
560 0 : E_n_sigma = E_sigma_tmp
561 0 : E_n_mu = E_mu_tmp
562 : END IF
563 :
564 : !-- the old/parent element energy estimation
565 : CALL three_point_extrapolate(v1=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 3, 4) + 1), &
566 : v2=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 2, 4) + 1), &
567 : v3=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 1, 4) + 1), &
568 0 : extrapolate=E_mu_tmp, res_err=E_sigma_tmp)
569 0 : IF ((elem_old%scf_energies_count > 3)) THEN
570 : CALL three_point_extrapolate(v1=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 4, 4) + 1), &
571 : v2=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 3, 4) + 1), &
572 : v3=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 2, 4) + 1), &
573 0 : extrapolate=E_o_mu, res_err=E_o_sigma)
574 0 : E_o_sigma = MAX(E_o_sigma, ABS(E_o_mu - E_mu_tmp))
575 : ELSE
576 0 : E_o_sigma = E_sigma_tmp
577 0 : E_o_mu = E_mu_tmp
578 : END IF
579 :
580 : ! calculate the estimation for the average of the trajectory elements
581 : prior_sigma = SQRT(ABS(tmc_params%prior_NMC_acc%aver_2 &
582 0 : - tmc_params%prior_NMC_acc%aver**2))
583 :
584 : ! calculate the probability of acceptance for those two elements with their energy
585 : ! swap and 2 potential moves are distinguished using the difference in classical energy and different betas
586 : prob = compute_prob(E_n_mu=E_n_mu, E_n_sigma=E_n_sigma, E_o_mu=E_o_mu, E_o_sigma=E_o_sigma, &
587 : E_classical_diff=E_classical_diff, &
588 : prior_mu=tmc_params%prior_NMC_acc%aver, prior_sigma=prior_sigma, &
589 0 : p=rnd_nr, beta=beta)
590 : END IF
591 : ! end the timing
592 0 : CALL timestop(handle)
593 0 : END FUNCTION compute_estimated_prob
594 :
595 : ! **************************************************************************************************
596 : !> \brief calculated the rate of used tree elements to created tree elements
597 : !> for every temperature
598 : !> \param tmc_env TMC environment variables
599 : !> \param eff result efficiency
600 : !> \author Mandes 01.2013
601 : ! **************************************************************************************************
602 14 : SUBROUTINE get_subtree_efficiency(tmc_env, eff)
603 : TYPE(tmc_env_type), POINTER :: tmc_env
604 : REAL(KIND=dp), DIMENSION(:), POINTER :: eff
605 :
606 : INTEGER :: i
607 :
608 14 : CPASSERT(ASSOCIATED(tmc_env))
609 14 : CPASSERT(ASSOCIATED(tmc_env%params))
610 14 : CPASSERT(ASSOCIATED(tmc_env%m_env))
611 :
612 54 : eff(:) = 0.0_dp
613 :
614 40 : DO i = 1, tmc_env%params%nr_temp
615 26 : IF (tmc_env%m_env%tree_node_count(i) > 0) THEN
616 : eff(i) = tmc_env%params%move_types%mv_count(0, i)/ &
617 24 : (tmc_env%m_env%tree_node_count(i)*1.0_dp)
618 : END IF
619 : eff(0) = eff(0) + tmc_env%params%move_types%mv_count(0, i)/ &
620 102 : (SUM(tmc_env%m_env%tree_node_count(1:))*1.0_dp)
621 : END DO
622 14 : END SUBROUTINE get_subtree_efficiency
623 : END MODULE tmc_calculations
|
