38 #include "../base/base_uses.f90"
44 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'tmc_calculations'
67 TYPE(tree_type),
POINTER :: conf
69 LOGICAL :: exact_approx_pot
70 TYPE(tmc_env_type),
POINTER :: tmc_env
74 REAL(kind=
dp) :: e_pot, rnd
75 TYPE(cell_type),
POINTER :: tmp_cell
79 cpassert(
ASSOCIATED(conf))
80 cpassert(env_id .GT. 0)
81 cpassert(
ASSOCIATED(tmc_env))
83 SELECT CASE (tmc_env%params%task_type)
87 IF (tmc_env%params%pressure .GE. 0.0_dp)
THEN
89 CALL get_scaled_cell(cell=tmc_env%params%cell, box_scale=conf%box_scale, &
91 CALL set_cell(env_id=env_id, new_cell=tmp_cell%hmat, ierr=ierr)
98 IF (flag .EQV. .true.)
THEN
99 IF (tmc_env%params%print_forces .OR. &
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 n_el_force=
SIZE(conf%frc), ierr=ierr)
108 CALL calc_energy(env_id=env_id, pos=conf%pos, n_el=
SIZE(conf%pos), e_pot=e_pot, ierr=ierr)
116 CALL cp_abort(__location__, &
117 "worker task typ is unknown "// &
118 cp_to_string(tmc_env%params%task_type))
122 rnd = tmc_env%rng_stream%next()
133 IF (exact_approx_pot)
THEN
134 conf%potential = e_pot
136 conf%e_pot_approx = e_pot
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
161 LOGICAL :: new_scaled_cell
162 TYPE(cell_type),
POINTER :: tmp_cell
164 cpassert(
ASSOCIATED(cell))
165 cpassert(
ASSOCIATED(box_scale))
167 new_scaled_cell = .false.
169 IF (.NOT.
PRESENT(scaled_cell))
THEN
171 new_scaled_cell = .true.
173 tmp_cell => scaled_cell
175 CALL cell_copy(cell_in=cell, cell_out=tmp_cell)
176 tmp_cell%hmat(:, 1) = tmp_cell%hmat(:, 1)*box_scale(1)
177 tmp_cell%hmat(:, 2) = tmp_cell%hmat(:, 2)*box_scale(2)
178 tmp_cell%hmat(:, 3) = tmp_cell%hmat(:, 3)*box_scale(3)
181 IF (
PRESENT(scaled_hmat)) &
182 scaled_hmat(:, :) = tmp_cell%hmat
184 IF (
PRESENT(vec))
THEN
185 vec =
pbc(r=vec, cell=tmp_cell)
188 IF (
PRESENT(vol))
CALL get_cell(cell=tmp_cell, deth=vol)
189 IF (
PRESENT(abc))
CALL get_cell(cell=tmp_cell, abc=abc)
190 IF (new_scaled_cell)
DEALLOCATE (tmp_cell)
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
206 REAL(kind=
dp),
DIMENSION(3) :: abc_new, abc_orig
207 TYPE(cell_type),
POINTER :: tmp_cell
209 cpassert(
ASSOCIATED(cell))
212 CALL cell_copy(cell_in=cell, cell_out=tmp_cell)
213 tmp_cell%hmat(:, :) = scaled_hmat(:, :)
215 CALL get_cell(cell=cell, abc=abc_orig)
216 CALL get_cell(cell=tmp_cell, abc=abc_new)
218 box_scale(:) = abc_new(:)/abc_orig(:)
220 DEALLOCATE (tmp_cell)
233 REAL(kind=
dp),
DIMENSION(:) :: x1, x2
234 TYPE(cell_type),
POINTER :: cell
235 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: box_scale
238 REAL(kind=
dp),
DIMENSION(3) :: dist_vec
239 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tmp_box_scale
241 NULLIFY (tmp_box_scale)
243 cpassert(
ASSOCIATED(cell))
244 cpassert(
SIZE(x1) .EQ. 3)
245 cpassert(
SIZE(x2) .EQ. 3)
247 dist_vec(:) = x2(:) - x1(:)
248 ALLOCATE (tmp_box_scale(3))
249 IF (
PRESENT(box_scale))
THEN
250 cpassert(
SIZE(box_scale) .EQ. 3)
251 tmp_box_scale(:) = box_scale
253 tmp_box_scale(:) = 1.0_dp
256 res = sqrt(sum(dist_vec(:)*dist_vec(:)))
257 DEALLOCATE (tmp_box_scale)
268 REAL(kind=
dp),
DIMENSION(:) :: pos
269 REAL(kind=
dp),
DIMENSION(:),
POINTER :: center
271 CHARACTER(LEN=*),
PARAMETER :: routinen =
'geometrical_center'
275 cpassert(
ASSOCIATED(center))
276 cpassert(
SIZE(pos) .GE.
SIZE(center))
279 CALL timeset(routinen, handle)
282 DO i = 1,
SIZE(pos),
SIZE(center)
283 center(:) = center(:) + &
284 pos(i:i +
SIZE(center) - 1)/(
SIZE(pos)/real(
SIZE(center), kind=
dp))
287 CALL timestop(handle)
301 REAL(kind=
dp),
DIMENSION(:) :: pos
302 TYPE(tmc_atom_type),
DIMENSION(:),
OPTIONAL :: atoms
303 REAL(kind=
dp),
DIMENSION(:),
POINTER :: center
305 CHARACTER(LEN=*),
PARAMETER :: routinen =
'center_of_mass'
308 REAL(kind=
dp) :: mass_sum, mass_tmp
310 cpassert(
ASSOCIATED(center))
311 cpassert(
SIZE(pos) .GE.
SIZE(center))
314 CALL timeset(routinen, handle)
318 DO i = 1,
SIZE(pos),
SIZE(center)
319 IF (
PRESENT(atoms))
THEN
320 cpassert(
SIZE(atoms) .EQ.
SIZE(pos)/
SIZE(center))
321 mass_tmp = atoms(int(i/real(
SIZE(center), kind=
dp)) + 1)%mass
322 center(:) = center(:) + pos(i:i +
SIZE(center) - 1)/ &
323 (
SIZE(pos)/real(
SIZE(center), kind=
dp))*mass_tmp
324 mass_sum = mass_sum + mass_tmp
326 cpwarn(
"try to calculate center of mass without any mass.")
327 center(:) = center(:) + pos(i:i +
SIZE(center) - 1)/ &
328 (
SIZE(pos)/real(
SIZE(center), kind=
dp))
332 center(:) = center(:)/mass_sum
334 CALL timestop(handle)
347 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
355 REAL(kind=
dp) :: kb, mass_tmp, rnd1, rnd2
359 cpassert(
ASSOCIATED(vel))
360 cpassert(
ASSOCIATED(atoms))
362 CALL rng_stream%set(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
364 rnd1 = rng_stream%next()
365 rnd2 = rng_stream%next()
367 mass_tmp = atoms(int(i/real(3, kind=
dp)) + 1)%mass
369 vel(i) = sqrt(-2.0_dp*log(rnd1))*cos(2.0_dp*
pi*rnd2)* &
370 sqrt(kb*temerature/mass_tmp)
372 CALL rng_stream%get(bg=rnd_seed(:, :, 1), cg=rnd_seed(:, :, 2), ig=rnd_seed(:, :, 3))
385 REAL(kind=
dp),
DIMENSION(:),
POINTER :: vel
386 TYPE(tmc_atom_type),
DIMENSION(:),
POINTER :: atoms
387 REAL(kind=
dp) :: ekin
390 REAL(kind=
dp) :: mass_tmp
392 cpassert(
ASSOCIATED(vel))
393 cpassert(
ASSOCIATED(atoms))
397 mass_tmp = atoms(int(i/real(3, kind=
dp)) + 1)%mass
398 ekin = ekin + 0.5_dp*mass_tmp*vel(i)*vel(i)
412 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
416 REAL(kind=
dp) :: e1, e2, e3
417 REAL(kind=
dp) :: a, b, c, d12, d23, ddd
419 extrapolate = huge(extrapolate)
429 e1 = v1; e2 = v2; e3 = v3
438 IF (d12 == 0 .OR. d23 == 0 .OR. abs(ddd) == 0)
THEN
444 b = (d12**3/(d23*ddd))
445 c = e2 - (d12*d23)/ddd
448 extrapolate = a**7*b + c
449 res_err = e3 - extrapolate
451 cpassert(extrapolate .NE. huge(extrapolate))
458 SUBROUTINE swap(x1, x2)
459 REAL(kind=
dp) :: x1, x2
469 END SUBROUTINE three_point_extrapolate
489 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
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 (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 prob = min(1.0_dp - epsilon(1.0_dp), max(epsilon(1.0_dp), prob))
505 END FUNCTION compute_prob
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
527 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_estimated_prob'
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
533 cpassert(
ASSOCIATED(elem_old))
534 cpassert(
ASSOCIATED(elem_new))
535 cpassert(rnd_nr .GT. 0.0_dp)
538 CALL timeset(routinen, handle)
541 IF ((elem_new%scf_energies_count .GE. 3) .AND. &
542 (elem_old%scf_energies_count .GE. 3) .AND. &
543 tmc_params%prior_NMC_acc%counter .GE. 10)
THEN
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 extrapolate=e_mu_tmp, res_err=e_sigma_tmp)
552 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 extrapolate=e_n_mu, res_err=e_n_sigma)
557 e_n_sigma = max(e_n_sigma, abs(e_n_mu - e_mu_tmp))
559 e_n_sigma = e_sigma_tmp
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 extrapolate=e_mu_tmp, res_err=e_sigma_tmp)
568 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 extrapolate=e_o_mu, res_err=e_o_sigma)
573 e_o_sigma = max(e_o_sigma, abs(e_o_mu - e_mu_tmp))
575 e_o_sigma = e_sigma_tmp
580 prior_sigma = sqrt(abs(tmc_params%prior_NMC_acc%aver_2 &
581 - tmc_params%prior_NMC_acc%aver**2))
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, &
591 CALL timestop(handle)
602 TYPE(tmc_env_type),
POINTER :: tmc_env
603 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eff
607 cpassert(
ASSOCIATED(tmc_env))
608 cpassert(
ASSOCIATED(tmc_env%params))
609 cpassert(
ASSOCIATED(tmc_env%m_env))
613 DO i = 1, tmc_env%params%nr_temp
614 IF (tmc_env%m_env%tree_node_count(i) > 0) &
615 eff(i) = tmc_env%params%move_types%mv_count(0, i)/ &
616 (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 (sum(tmc_env%m_env%tree_node_count(1:))*1.0_dp)
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Handles all functions related to the CELL.
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
Handles all functions related to the CELL.
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
subroutine, public cell_copy(cell_in, cell_out, tag)
Copy cell variable.
various routines to log and control the output. The idea is that decisions about where to log should ...
interface to use cp2k as library
recursive subroutine, public calc_energy(env_id, pos, n_el, e_pot, ierr)
returns the energy of the configuration given by the positions passed as argument
subroutine, public set_cell(env_id, new_cell, ierr)
sets a new cell
recursive subroutine, public calc_force(env_id, pos, n_el_pos, e_pot, force, n_el_force, ierr)
returns the energy of the configuration given by the positions passed as argument
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition of physical constants:
real(kind=dp), parameter, public boltzmann
real(kind=dp), parameter, public joule
calculation section for TreeMonteCarlo
subroutine, public init_vel(vel, atoms, temerature, rng_stream, rnd_seed)
routine sets initial velocity, using the Box-Muller Method for Normal (Gaussian) Deviates
subroutine, public geometrical_center(pos, center)
calculate the geometrical center of an amount of atoms array size should be multiple of dim_per_elem
subroutine, public get_scaled_cell(cell, box_scale, scaled_hmat, scaled_cell, vol, abc, vec)
handles properties and calculations of a scaled cell
real(kind=dp) function, public compute_estimated_prob(elem_old, elem_new, E_classical_diff, rnd_nr, beta, tmc_params)
extimates the probability of acceptance considering the intermetiate step energies
subroutine, public get_subtree_efficiency(tmc_env, eff)
calculated the rate of used tree elements to created tree elements for every temperature
subroutine, public get_cell_scaling(cell, scaled_hmat, box_scale)
handles properties and calculations of a scaled cell
real(kind=dp) function, public calc_e_kin(vel, atoms)
routine calculates the kinetic energy, using the velocities and atom mass, both in atomic units
subroutine, public center_of_mass(pos, atoms, center)
calculate the center of mass of an amount of atoms array size should be multiple of dim_per_elem
subroutine, public calc_potential_energy(conf, env_id, exact_approx_pot, tmc_env)
start the calculation of the energy (distinguish between exact and approximate)
real(kind=dp) function, public nearest_distance(x1, x2, cell, box_scale)
neares distance of atoms within the periodic boundary condition
tree nodes creation, searching, deallocation, references etc.
integer, parameter, public mv_type_md
tree nodes creation, searching, deallocation, references etc.
integer, parameter, public task_type_gaussian_adaptation
integer, parameter, public task_type_mc
integer, parameter, public task_type_ideal_gas
module handles definition of the tree nodes for the global and the subtrees binary tree parent elemen...
module handles definition of the tree nodes for the global and the subtrees binary tree parent elemen...