Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 4771 : 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 4771 : rnd = 0.0_dp
78 :
79 4771 : CPASSERT(ASSOCIATED(conf))
80 4771 : CPASSERT(env_id .GT. 0)
81 4771 : CPASSERT(ASSOCIATED(tmc_env))
82 :
83 9542 : SELECT CASE (tmc_env%params%task_type)
84 : CASE (task_type_gaussian_adaptation)
85 : !CALL gaussian_adaptation_energy(, )
86 : CASE (task_type_MC)
87 4771 : IF (tmc_env%params%pressure .GE. 0.0_dp) THEN
88 23340 : 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 .EQ. 0)
93 778 : DEALLOCATE (tmp_cell)
94 : END IF
95 :
96 : ! TODO check for minimal distances
97 4771 : flag = .TRUE.
98 0 : IF (flag .EQV. .TRUE.) THEN
99 4771 : IF (tmc_env%params%print_forces .OR. &
100 : conf%move_type .EQ. 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 4173 : 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 4771 : cp_to_string(tmc_env%params%task_type))
119 : END SELECT
120 :
121 : ! --- wait a bit
122 4771 : 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 4771 : IF (exact_approx_pot) THEN
134 4574 : conf%potential = e_pot
135 : ELSE
136 197 : conf%e_pot_approx = e_pot
137 : END IF
138 4771 : 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 6527760 : 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)) &
182 5096 : scaled_hmat(:, :) = tmp_cell%hmat
183 :
184 218382 : IF (PRESENT(vec)) THEN
185 862448 : vec = pbc(r=vec, cell=tmp_cell)
186 : END IF
187 :
188 218382 : IF (PRESENT(vol)) CALL get_cell(cell=tmp_cell, deth=vol)
189 218382 : IF (PRESENT(abc)) CALL get_cell(cell=tmp_cell, abc=abc)
190 218382 : IF (new_scaled_cell) DEALLOCATE (tmp_cell)
191 :
192 218382 : END SUBROUTINE get_scaled_cell
193 :
194 : ! **************************************************************************************************
195 : !> \brief handles properties and calculations of a scaled cell
196 : !> \param cell original cell
197 : !> \param scaled_hmat returns the scaled h matrix (matrix of cell vectors)
198 : !> \param box_scale scaling factors for each direction
199 : !> \author Mandes 11.2012
200 : ! **************************************************************************************************
201 1206 : SUBROUTINE get_cell_scaling(cell, scaled_hmat, box_scale)
202 : TYPE(cell_type), INTENT(IN), POINTER :: cell
203 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: scaled_hmat
204 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: box_scale
205 :
206 : REAL(KIND=dp), DIMENSION(3) :: abc_new, abc_orig
207 : TYPE(cell_type), POINTER :: tmp_cell
208 :
209 1206 : CPASSERT(ASSOCIATED(cell))
210 :
211 36180 : ALLOCATE (tmp_cell)
212 1206 : CALL cell_copy(cell_in=cell, cell_out=tmp_cell)
213 15678 : tmp_cell%hmat(:, :) = scaled_hmat(:, :)
214 1206 : CALL init_cell(cell=tmp_cell)
215 1206 : CALL get_cell(cell=cell, abc=abc_orig)
216 1206 : CALL get_cell(cell=tmp_cell, abc=abc_new)
217 :
218 4824 : box_scale(:) = abc_new(:)/abc_orig(:)
219 :
220 1206 : DEALLOCATE (tmp_cell)
221 1206 : END SUBROUTINE get_cell_scaling
222 :
223 : ! **************************************************************************************************
224 : !> \brief neares distance of atoms within the periodic boundary condition
225 : !> \param x1 ...
226 : !> \param x2 ...
227 : !> \param cell ...
228 : !> \param box_scale ...
229 : !> \return ...
230 : !> \author Mandes 11.2012
231 : ! **************************************************************************************************
232 187940 : FUNCTION nearest_distance(x1, x2, cell, box_scale) RESULT(res)
233 : REAL(KIND=dp), DIMENSION(:) :: x1, x2
234 : TYPE(cell_type), POINTER :: cell
235 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: box_scale
236 : REAL(KIND=dp) :: res
237 :
238 : REAL(KIND=dp), DIMENSION(3) :: dist_vec
239 187940 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_box_scale
240 :
241 187940 : NULLIFY (tmp_box_scale)
242 :
243 0 : CPASSERT(ASSOCIATED(cell))
244 187940 : CPASSERT(SIZE(x1) .EQ. 3)
245 187940 : CPASSERT(SIZE(x2) .EQ. 3)
246 :
247 751760 : dist_vec(:) = x2(:) - x1(:) ! distance vector between atoms
248 187940 : ALLOCATE (tmp_box_scale(3))
249 187940 : IF (PRESENT(box_scale)) THEN
250 187940 : CPASSERT(SIZE(box_scale) .EQ. 3)
251 1315580 : tmp_box_scale(:) = box_scale
252 : ELSE
253 0 : tmp_box_scale(:) = 1.0_dp
254 : END IF
255 187940 : CALL get_scaled_cell(cell=cell, box_scale=box_scale, vec=dist_vec)
256 751760 : res = SQRT(SUM(dist_vec(:)*dist_vec(:)))
257 187940 : DEALLOCATE (tmp_box_scale)
258 187940 : END FUNCTION nearest_distance
259 :
260 : ! **************************************************************************************************
261 : !> \brief calculate the geometrical center of an amount of atoms
262 : !> array size should be multiple of dim_per_elem
263 : !> \param pos list of atoms
264 : !> \param center return value, the geometrical center
265 : !> \author Mandes 11.2012
266 : ! **************************************************************************************************
267 7760 : SUBROUTINE geometrical_center(pos, center)
268 : REAL(KIND=dp), DIMENSION(:) :: pos
269 : REAL(KIND=dp), DIMENSION(:), POINTER :: center
270 :
271 : CHARACTER(LEN=*), PARAMETER :: routineN = 'geometrical_center'
272 :
273 : INTEGER :: handle, i
274 :
275 7760 : CPASSERT(ASSOCIATED(center))
276 7760 : CPASSERT(SIZE(pos) .GE. SIZE(center))
277 :
278 : ! start the timing
279 7760 : CALL timeset(routineN, handle)
280 :
281 31040 : center = 0.0_dp
282 31138 : DO i = 1, SIZE(pos), SIZE(center)
283 : center(:) = center(:) + &
284 101272 : pos(i:i + SIZE(center) - 1)/(SIZE(pos)/REAL(SIZE(center), KIND=dp))
285 : END DO
286 : ! end the timing
287 7760 : CALL timestop(handle)
288 7760 : END SUBROUTINE geometrical_center
289 :
290 : ! **************************************************************************************************
291 : !> \brief calculate the center of mass of an amount of atoms
292 : !> array size should be multiple of dim_per_elem
293 : !> \param pos ...
294 : !> \param atoms ...
295 : !> \param center ...
296 : !> \param
297 : !> \param
298 : !> \author Mandes 11.2012
299 : ! **************************************************************************************************
300 0 : SUBROUTINE center_of_mass(pos, atoms, center)
301 : REAL(KIND=dp), DIMENSION(:) :: pos
302 : TYPE(tmc_atom_type), DIMENSION(:), OPTIONAL :: atoms
303 : REAL(KIND=dp), DIMENSION(:), POINTER :: center
304 :
305 : CHARACTER(LEN=*), PARAMETER :: routineN = 'center_of_mass'
306 :
307 : INTEGER :: handle, i
308 : REAL(KIND=dp) :: mass_sum, mass_tmp
309 :
310 0 : CPASSERT(ASSOCIATED(center))
311 0 : CPASSERT(SIZE(pos) .GE. SIZE(center))
312 :
313 : ! start the timing
314 0 : CALL timeset(routineN, handle)
315 :
316 0 : center = 0.0_dp
317 0 : mass_sum = 0.0_dp
318 0 : DO i = 1, SIZE(pos), SIZE(center)
319 0 : IF (PRESENT(atoms)) THEN
320 0 : CPASSERT(SIZE(atoms) .EQ. SIZE(pos)/SIZE(center))
321 0 : mass_tmp = atoms(INT(i/REAL(SIZE(center), KIND=dp)) + 1)%mass
322 : center(:) = center(:) + pos(i:i + SIZE(center) - 1)/ &
323 0 : (SIZE(pos)/REAL(SIZE(center), KIND=dp))*mass_tmp
324 0 : mass_sum = mass_sum + mass_tmp
325 : ELSE
326 0 : CPWARN("try to calculate center of mass without any mass.")
327 : center(:) = center(:) + pos(i:i + SIZE(center) - 1)/ &
328 0 : (SIZE(pos)/REAL(SIZE(center), KIND=dp))
329 0 : mass_sum = 1.0_dp
330 : END IF
331 : END DO
332 0 : center(:) = center(:)/mass_sum
333 : ! end the timing
334 0 : CALL timestop(handle)
335 0 : END SUBROUTINE center_of_mass
336 :
337 : ! **************************************************************************************************
338 : !> \brief routine sets initial velocity, using the Box-Muller Method for Normal
339 : !> (Gaussian) Deviates
340 : !> \param vel ...
341 : !> \param atoms ...
342 : !> \param temerature ...
343 : !> \param rng_stream ...
344 : !> \param rnd_seed ...
345 : !> \author Mandes 11.2012
346 : ! **************************************************************************************************
347 0 : SUBROUTINE init_vel(vel, atoms, temerature, rng_stream, rnd_seed)
348 : REAL(KIND=dp), DIMENSION(:), POINTER :: vel
349 : TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
350 : REAL(KIND=dp) :: temerature
351 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
352 : REAL(KIND=dp), DIMENSION(3, 2, 3) :: rnd_seed
353 :
354 : INTEGER :: i
355 : REAL(KIND=dp) :: kB, mass_tmp, rnd1, rnd2
356 :
357 0 : kB = boltzmann/joule
358 :
359 0 : CPASSERT(ASSOCIATED(vel))
360 0 : CPASSERT(ASSOCIATED(atoms))
361 :
362 0 : CALL rng_stream%set(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
363 0 : DO i = 1, SIZE(vel)
364 0 : rnd1 = rng_stream%next()
365 0 : rnd2 = rng_stream%next()
366 :
367 0 : mass_tmp = atoms(INT(i/REAL(3, KIND=dp)) + 1)%mass
368 :
369 : vel(i) = SQRT(-2.0_dp*LOG(rnd1))*COS(2.0_dp*PI*rnd2)* &
370 0 : SQRT(kB*temerature/mass_tmp)
371 : END DO
372 0 : CALL rng_stream%get(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
373 :
374 0 : END SUBROUTINE init_vel
375 :
376 : ! **************************************************************************************************
377 : !> \brief routine calculates the kinetic energy, using the velocities
378 : !> and atom mass, both in atomic units
379 : !> \param vel ...
380 : !> \param atoms ...
381 : !> \return ...
382 : !> \author Mandes 11.2012
383 : ! **************************************************************************************************
384 0 : FUNCTION calc_e_kin(vel, atoms) RESULT(ekin)
385 : REAL(KIND=dp), DIMENSION(:), POINTER :: vel
386 : TYPE(tmc_atom_type), DIMENSION(:), POINTER :: atoms
387 : REAL(KIND=dp) :: ekin
388 :
389 : INTEGER :: i
390 : REAL(KIND=dp) :: mass_tmp
391 :
392 0 : CPASSERT(ASSOCIATED(vel))
393 0 : CPASSERT(ASSOCIATED(atoms))
394 0 : ekin = 0.0_dp
395 :
396 0 : DO i = 1, SIZE(vel)
397 0 : mass_tmp = atoms(INT(i/REAL(3, KIND=dp)) + 1)%mass
398 0 : ekin = ekin + 0.5_dp*mass_tmp*vel(i)*vel(i)
399 : END DO
400 0 : END FUNCTION calc_e_kin
401 :
402 : ! **************************************************************************************************
403 : !> \brief assuming an (exponential) decreasing function, this function
404 : !> extrapolate the converged value
405 : !> \param v1 function values
406 : !> \param v2 function values
407 : !> \param v3 function values
408 : !> \param extrapolate extrapolated final value (result)
409 : !> \param res_err error of the result
410 : !> \author Mandes 12.2012
411 : ! **************************************************************************************************
412 0 : SUBROUTINE three_point_extrapolate(v1, v2, v3, extrapolate, res_err)
413 : REAL(KIND=dp) :: v1, v2, v3
414 : REAL(KIND=dp), INTENT(OUT) :: extrapolate, res_err
415 :
416 : REAL(KIND=dp) :: e1, e2, e3
417 : REAL(KIND=dp) :: a, b, c, d12, d23, ddd
418 :
419 : extrapolate = HUGE(extrapolate)
420 :
421 : !> solve({exp(a+b)+c = e1, exp(2*a+b)+c = e2, exp(3*a+b)+c = e3}, [a, b, c])
422 : !> solve({a*b+c = e1, a^2*b+c = e2, a^3*b+c = e3}, [a, b, c]);
423 : ! [[ 3 2 ]]
424 : ! [[ -e3 + e2 (e1 - e2) -e2 + e1 e3 ]]
425 : ! [[a = --------, b = ---------------------------, c = --------------]]
426 : ! [[ e1 - e2 (-e3 + e2) (e3 - 2 e2 + e1) e3 - 2 e2 + e1]]
427 :
428 : ! sort so that e1>=e2>=e3
429 0 : e1 = v1; e2 = v2; e3 = v3
430 0 : CALL swap(e1, e2)
431 0 : CALL swap(e1, e3)
432 0 : CALL swap(e2, e3)
433 : ! we need extra care if some of the difference e1-e2, e3-e2 are nearly zero,
434 : ! since the formulae suffer from sever loss of precision
435 0 : d12 = e1 - e2
436 0 : d23 = e2 - e3
437 0 : ddd = d12 - d23
438 0 : IF (d12 == 0 .OR. d23 == 0 .OR. ABS(ddd) == 0) THEN
439 : ! a degenerate case, we do no extrapolation
440 0 : extrapolate = e3
441 0 : res_err = e1 - e3
442 : ELSE
443 0 : a = d23/d12
444 0 : b = (d12**3/(d23*ddd))
445 0 : c = e2 - (d12*d23)/ddd
446 : ! extrapolation, let's only look 4 iterations ahead, more is presumably anyway not accurate
447 : ! fewer is maybe more stable
448 0 : extrapolate = a**7*b + c
449 0 : res_err = e3 - extrapolate
450 : END IF
451 0 : CPASSERT(extrapolate .NE. HUGE(extrapolate))
452 : CONTAINS
453 : ! **************************************************************************************************
454 : !> \brief ...
455 : !> \param x1 ...
456 : !> \param x2 ...
457 : ! **************************************************************************************************
458 0 : SUBROUTINE swap(x1, x2)
459 : REAL(KIND=dp) :: x1, x2
460 :
461 : REAL(KIND=dp) :: tmp
462 :
463 0 : IF (x2 > x1) THEN
464 0 : tmp = x2
465 0 : x2 = x1
466 0 : x1 = tmp
467 : END IF
468 0 : END SUBROUTINE swap
469 : END SUBROUTINE three_point_extrapolate
470 :
471 : ! **************************************************************************************************
472 : !> \brief calculates the probability of acceptance for given intervals of the
473 : !> exact energy
474 : !> \param E_n_mu energy distribution of new configuration
475 : !> \param E_n_sigma energy distribution of new configuration
476 : !> \param E_o_mu energy distribution of old configuration
477 : !> \param E_o_sigma energy distribution of old configuration
478 : !> \param E_classical_diff the difference in approximated energies for the
479 : !> old and new configuration (E_o-E_n)
480 : !> \param prior_mu energy distribution of the already converged
481 : !> energies
482 : !> \param prior_sigma energy distribution of the already converged
483 : !> energies
484 : !> \param p the random number, the criteria has to be smaller than this
485 : !> \param beta ...
486 : !> \return return probability of acceptance
487 : !> \author Mandes 12.2012
488 : ! **************************************************************************************************
489 0 : FUNCTION compute_prob(E_n_mu, E_n_sigma, E_o_mu, E_o_sigma, E_classical_diff, &
490 : prior_mu, prior_sigma, p, beta) RESULT(prob)
491 : REAL(KIND=dp) :: E_n_mu, E_n_sigma, E_o_mu, E_o_sigma, &
492 : E_classical_diff, prior_mu, &
493 : prior_sigma, p, beta, prob
494 :
495 : ! INTEGER :: io,in
496 : ! REAL(KIND=dp) :: diff,E_n,E_o,surface,lower_bound,upper_bound,delta
497 :
498 : prob = 0.5_dp*ERFC(-0.5_dp*SQRT(2.0_dp)*( &
499 : (-prior_sigma**2 - E_o_sigma**2 - E_n_sigma**2)*LOG(p) + &
500 : ((E_classical_diff - E_n_mu + E_o_mu)*prior_sigma**2 - prior_mu*(E_n_sigma**2 + E_o_sigma**2))*beta)/ &
501 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))
502 :
503 0 : prob = MIN(1.0_dp - EPSILON(1.0_dp), MAX(EPSILON(1.0_dp), prob))
504 :
505 0 : END FUNCTION compute_prob
506 :
507 : ! **************************************************************************************************
508 : !> \brief extimates the probability of acceptance considering the intermetiate
509 : !> step energies
510 : !> \param elem_old old/parent sub tree element
511 : !> \param elem_new new/actual sub tree element, which schould be checked
512 : !> \param E_classical_diff difference in the classical energy of the old and
513 : !> new configuration
514 : !> \param rnd_nr random number acceptance check will be done with
515 : !> \param beta 1/(kB*T) can differ for different acceptance checks
516 : !> \param tmc_params TMC environment parameters
517 : !> \return estimated acceptance probability
518 : !> \author Mandes 12.2012
519 : ! **************************************************************************************************
520 0 : FUNCTION compute_estimated_prob(elem_old, elem_new, E_classical_diff, &
521 : rnd_nr, beta, tmc_params) RESULT(prob)
522 : TYPE(tree_type), POINTER :: elem_old, elem_new
523 : REAL(KIND=dp) :: E_classical_diff, rnd_nr, beta
524 : TYPE(tmc_param_type), POINTER :: tmc_params
525 : REAL(KIND=dp) :: prob
526 :
527 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_estimated_prob'
528 :
529 : INTEGER :: handle
530 : REAL(KIND=dp) :: E_mu_tmp, E_n_mu, E_n_sigma, E_o_mu, &
531 : E_o_sigma, E_sigma_tmp, prior_sigma
532 :
533 0 : CPASSERT(ASSOCIATED(elem_old))
534 0 : CPASSERT(ASSOCIATED(elem_new))
535 0 : CPASSERT(rnd_nr .GT. 0.0_dp)
536 :
537 : ! start the timing
538 0 : CALL timeset(routineN, handle)
539 :
540 0 : prob = -1.0_dp
541 : IF ((elem_new%scf_energies_count .GE. 3) .AND. &
542 0 : (elem_old%scf_energies_count .GE. 3) .AND. &
543 : tmc_params%prior_NMC_acc%counter .GE. 10) THEN
544 : !-- first the new element energy estimation
545 : ! using 3 point extrapolation of two different intervals -> more stable estimation
546 : ! the energies are sorted in the three_point_extrapolate routine !
547 : ! But with array of length 4 we have to select the 3 connected ones
548 : CALL three_point_extrapolate(v1=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 3, 4) + 1), &
549 : v2=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 2, 4) + 1), &
550 : v3=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 1, 4) + 1), &
551 0 : extrapolate=E_mu_tmp, res_err=E_sigma_tmp)
552 0 : IF ((elem_new%scf_energies_count .GT. 3)) THEN
553 : CALL three_point_extrapolate(v1=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 4, 4) + 1), &
554 : v2=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 3, 4) + 1), &
555 : v3=elem_new%scf_energies(MOD(elem_new%scf_energies_count - 2, 4) + 1), &
556 0 : extrapolate=E_n_mu, res_err=E_n_sigma)
557 0 : E_n_sigma = MAX(E_n_sigma, ABS(E_n_mu - E_mu_tmp))
558 : ELSE
559 0 : E_n_sigma = E_sigma_tmp
560 0 : E_n_mu = E_mu_tmp
561 : END IF
562 :
563 : !-- the old/parent element energy estimation
564 : CALL three_point_extrapolate(v1=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 3, 4) + 1), &
565 : v2=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 2, 4) + 1), &
566 : v3=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 1, 4) + 1), &
567 0 : extrapolate=E_mu_tmp, res_err=E_sigma_tmp)
568 0 : IF ((elem_old%scf_energies_count .GT. 3)) THEN
569 : CALL three_point_extrapolate(v1=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 4, 4) + 1), &
570 : v2=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 3, 4) + 1), &
571 : v3=elem_old%scf_energies(MOD(elem_old%scf_energies_count - 2, 4) + 1), &
572 0 : extrapolate=E_o_mu, res_err=E_o_sigma)
573 0 : E_o_sigma = MAX(E_o_sigma, ABS(E_o_mu - E_mu_tmp))
574 : ELSE
575 0 : E_o_sigma = E_sigma_tmp
576 0 : E_o_mu = E_mu_tmp
577 : END IF
578 :
579 : ! calculate the estimation for the average of the trajectory elements
580 : prior_sigma = SQRT(ABS(tmc_params%prior_NMC_acc%aver_2 &
581 0 : - tmc_params%prior_NMC_acc%aver**2))
582 :
583 : ! calculate the probability of acceptance for those two elements with their energy
584 : ! swap and 2 potential moves are distinguished using the difference in classical energy and different betas
585 : 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, &
586 : E_classical_diff=E_classical_diff, &
587 : prior_mu=tmc_params%prior_NMC_acc%aver, prior_sigma=prior_sigma, &
588 0 : p=rnd_nr, beta=beta)
589 : END IF
590 : ! end the timing
591 0 : CALL timestop(handle)
592 0 : END FUNCTION compute_estimated_prob
593 :
594 : ! **************************************************************************************************
595 : !> \brief calculated the rate of used tree elements to created tree elements
596 : !> for every temperature
597 : !> \param tmc_env TMC environment variables
598 : !> \param eff result efficiency
599 : !> \author Mandes 01.2013
600 : ! **************************************************************************************************
601 14 : SUBROUTINE get_subtree_efficiency(tmc_env, eff)
602 : TYPE(tmc_env_type), POINTER :: tmc_env
603 : REAL(KIND=dp), DIMENSION(:), POINTER :: eff
604 :
605 : INTEGER :: i
606 :
607 14 : CPASSERT(ASSOCIATED(tmc_env))
608 14 : CPASSERT(ASSOCIATED(tmc_env%params))
609 14 : CPASSERT(ASSOCIATED(tmc_env%m_env))
610 :
611 54 : eff(:) = 0.0_dp
612 :
613 40 : DO i = 1, tmc_env%params%nr_temp
614 26 : IF (tmc_env%m_env%tree_node_count(i) > 0) &
615 : eff(i) = tmc_env%params%move_types%mv_count(0, i)/ &
616 24 : (tmc_env%m_env%tree_node_count(i)*1.0_dp)
617 : eff(0) = eff(0) + tmc_env%params%move_types%mv_count(0, i)/ &
618 102 : (SUM(tmc_env%m_env%tree_node_count(1:))*1.0_dp)
619 : END DO
620 14 : END SUBROUTINE get_subtree_efficiency
621 : END MODULE tmc_calculations
|
