29 use,
intrinsic :: iso_fortran_env
50 real(real64),
public :: A
51 real(real64),
public :: alpha
52 real(real64),
public :: beta
53 real(real64),
allocatable :: pos(:, :)
55 type(root_solver_t) :: rs
68 procedure curv_gygi_constructor
72 class(curv_gygi_t),
pointer :: gygi_p
73 type(curv_gygi_t),
target :: gygi_global
74 real(real64),
allocatable :: chi_p(:)
80 type(namespace_t),
intent(in) :: namespace
81 integer,
intent(in) :: dim
82 integer,
intent(in) :: npos
83 real(real64),
intent(in) :: pos(1:dim,1:npos)
84 class(curv_gygi_t),
pointer :: gygi
91 gygi%local_basis = .
true.
92 gygi%orthogonal = .
true.
95 safe_allocate(gygi%pos(1:dim, 1:gygi%npos))
96 gygi%pos(:, :) = pos(:, :)
142 gygi%min_mesh_scaling_product = (
m_one / (
m_one + gygi%A))**gygi%dim
157 this_out%A = this_in%A
158 this_out%alpha = this_in%alpha
159 this_out%beta = this_in%beta
160 safe_allocate_source_a(this_out%pos, this_in%pos)
161 this_out%pos = this_in%pos
162 this_out%npos = this_in%npos
163 this_out%dim = this_in%dim
164 this_out%local_basis = this_in%local_basis
165 this_out%orthogonal = this_in%orthogonal
167 maxiter = 500, abs_tolerance = 1.0e-10_real64)
178 safe_deallocate_a(this%pos)
184 real(real64)
pure function
gygi_f(this, r)
186 real(real64),
intent(in) :: r
190 if (r < 1.0e-4_real64)
then
196 gygi_f = this%A * this%alpha/r *
tanh(r/this%alpha) *
exp(-(r/this%beta)**2)
204 real(real64),
intent(in) :: r
208 if (r < 3.0e-3_real64)
then
211 gygi_dfdr_over_r = this%A * ((-m_two/(m_three*this%alpha**2) - m_two/this%beta**2) + &
212 (8._real64/(15._real64*this%alpha**4) + m_four/(m_three*this%alpha**2*this%beta**2) &
213 + m_two/this%beta**4) * r**2)
214 else if (-(r/this%beta)**2 <= m_min_exp_arg)
then
218 this%alpha*(this%beta**2+m_two*r**2) *
tanh(r/this%alpha)) &
219 / (this%beta**2*r**3) *
exp(-(r/this%beta)**2)
228 real(real64),
intent(in) :: r
232 if (r < 3.0e-3_real64)
then
236 m_two*(8._real64/(15._real64*this%alpha**4) + m_four/(m_three*this%alpha**2*this%beta**2) &
237 + m_two/this%beta**4) &
238 - m_two*(34._real64/(105._real64*this%alpha**6) + m_four/(m_five*this%alpha**4*this%beta**2) &
239 + m_one/(this%alpha**2*this%beta**4) + m_one/this%beta**6) * r**2 &
241 else if (-(r/this%beta)**2 <= m_min_exp_arg)
then
245 -m_two*
tanh(r/this%alpha)/(
cosh(r/this%alpha)**2*this%alpha*r) &
246 +
tanh(r/this%alpha)*(m_three*this%alpha/r**3+m_four*this%alpha/(this%beta**2*r)+m_four*this%alpha*r/this%beta**4) &
247 + m_one/
cosh(r/this%alpha)**2*(-m_three/r**2-m_four/this%beta**2))/r**2
255 real(real64),
intent(in) :: chi(:)
256 real(real64) :: xx(1:this%dim), xx_start(1:this%dim)
259 integer :: i_conv, n_conv
264 safe_allocate(
chi_p(1:this%dim))
267 call droot_solver_run(this%rs,
getf, xx, conv, startval = chi)
275 do i_conv = 1, n_conv
277 call droot_solver_run(
gygi_p%rs,
getf, xx, conv, startval = chi)
285 do i_conv = n_conv, 1, -1
288 call droot_solver_run(
gygi_p%rs,
getf, xx, conv, startval = xx_start)
298 message(1) =
"During the construction of the adaptive grid, the Newton-Raphson"
299 message(2) =
"method did not converge for point:"
300 write(message(3),
'(9f14.6)') xx(1:this%dim)
301 message(4) =
"Try varying the Gygi parameters -- usually reducing CurvGygiA or"
302 message(5) =
"CurvGygiAlpha (or both) solves the problem."
303 call messages_fatal(5)
311 real(real64),
intent(in) :: xx(:)
312 real(real64) :: chi(1:this%dim)
315 real(real64) :: diff(this%dim),
f
319 chi(1:this%dim) = xx(1:this%dim)
321 diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
324 chi(i) = chi(i) + diff(i) *
f
333 integer,
optional,
intent(in) :: iunit
334 type(namespace_t),
optional,
intent(in) :: namespace
338 write(message(1),
'(a)')
' Curvilinear Method = gygi'
339 write(message(2),
'(a)')
' Gygi Parameters:'
340 write(message(3),
'(4x,a,f6.3)')
'A = ', this%a
341 write(message(4),
'(4x,3a,f6.3)')
'alpha [', trim(units_abbrev(units_out%length)),
'] = ', &
342 units_from_atomic(units_out%length, this%alpha)
343 write(message(5),
'(4x,3a,f6.3)')
'beta [', trim(units_abbrev(units_out%length)),
'] = ', &
344 units_from_atomic(units_out%length, this%beta)
345 call messages_info(5, iunit=iunit, namespace=namespace)
353 integer,
intent(in) :: idir
356 message(1) =
'Surface element with gygi curvilinear coordinates not implemented'
357 call messages_fatal(1)
364 real(real64),
intent(in) :: chi(:)
365 real(real64) :: jacobian(1:this%dim, 1:this%dim)
367 jacobian = this%jacobian_inverse(chi)
368 call lalg_inverse(this%dim, jacobian,
"dir")
374 real(real64),
intent(in) :: chi(:)
375 real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
377 real(real64) :: xx(1:this%dim)
381 xx(:) = this%to_cartesian(chi)
389 real(real64),
intent(in) :: xx(:)
390 real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
393 real(real64) :: r, diff(1:this%dim), dfdr_over_r,
f
397 jacobian_inverse(1:this%dim, 1:this%dim) = diagonal_matrix(this%dim, m_one)
400 diff = xx(1:this%dim) - this%pos(1:this%dim, i)
407 jacobian_inverse(ix, ix) = jacobian_inverse(ix, ix) +
f
409 jacobian_inverse(ix, iy) = jacobian_inverse(ix, iy) + diff(ix)*diff(iy)*dfdr_over_r
420 real(real64),
intent(in) :: xx(:)
421 real(real64),
intent(out) :: hessian(:, :, :)
422 integer,
optional,
intent(in) :: natoms
424 integer :: ia, i, j, k, natoms_
425 real(real64) :: r, diff(this%dim), dfdr_over_r, d2fdr2_combination
429 hessian(1:this%dim, 1:this%dim, 1:this%dim) = m_zero
432 if (
present(natoms)) natoms_ = natoms
435 diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
445 hessian(i, j, k) = hessian(i, j, k) + diff(k)*dfdr_over_r
448 hessian(i, j, k) = hessian(i, j, k) + diff(j)*dfdr_over_r
451 hessian(i, j, k) = hessian(i, j, k) + diff(i)*dfdr_over_r
453 hessian(i, j, k) = hessian(i, j, k) + diff(i)*diff(j)*diff(k)*d2fdr2_combination
464 real(real64),
intent(in) :: chi(:)
465 real(real64) :: trace_hessian(1:this%dim)
468 real(real64) :: r, xx(this%dim), diff(this%dim), dfdr_over_r, d2fdr2_combination
472 xx(:) = this%to_cartesian(chi)
473 trace_hessian(:) = m_zero
476 diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
485 trace_hessian(i) = trace_hessian(i) + m_two*diff(k)*dfdr_over_r
487 trace_hessian(i) = trace_hessian(i) + diff(i)*dfdr_over_r
488 trace_hessian(i) = trace_hessian(i) + diff(i)*diff(k)*diff(k)*d2fdr2_combination
496 subroutine getf(y, f, jf)
497 real(real64),
intent(in) :: y(:)
498 real(real64),
intent(out) ::
f(:), jf(:, :)
double exp(double __x) __attribute__((__nothrow__
double tanh(double __x) __attribute__((__nothrow__
double cosh(double __x) __attribute__((__nothrow__
This module implements the curvilinear coordinates given in F. Gygi and G. Galli, PRB 52 R2229 (1996)...
type(curv_gygi_t), target gygi_global
subroutine curv_gygi_write_info(this, iunit, namespace)
real(real64), dimension(:), allocatable chi_p
real(real64) pure function gygi_f(this, r)
real(real64) pure function gygi_d2fdr2_combination(this, r)
real(real64) pure function gygi_dfdr_over_r(this, r)
pure subroutine, public curv_gygi_hessian(this, xx, hessian, natoms)
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_gygi_jacobian_inverse_cartesian(this, xx)
real(real64) function, dimension(1:this%dim) curv_gygi_trace_hessian(this, chi)
pure real(real64) function, dimension(1:this%dim) curv_gygi_from_cartesian(this, xx)
subroutine, public curv_gygi_copy(this_out, this_in)
real(real64) function, dimension(1:this%dim) curv_gygi_to_cartesian(this, chi)
subroutine getf(y, f, jf)
real(real64) function curv_gygi_surface_element(this, idir)
class(curv_gygi_t) function, pointer curv_gygi_constructor(namespace, dim, npos, pos)
class(curv_gygi_t), pointer gygi_p
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_gygi_jacobian(this, chi)
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_gygi_jacobian_inverse(this, chi)
subroutine curv_gygi_finalize(this)
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_min_exp_arg
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_three
This module is intended to contain "only mathematical" functions and procedures.
subroutine, public messages_input_error(namespace, var, details, row, column)
type(namespace_t), public global_namespace
integer, parameter, public root_newton
subroutine, public root_solver_init(rs, namespace, dimensionality, solver_type, maxiter, rel_tolerance, abs_tolerance)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
This module defines the unit system, used for input and output.
type(unit_system_t), public units_inp
the units systems for reading and writing
static double f(double w, void *p)
abstract class to describe coordinate systems