Octopus
curv_gygi.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2025 S. Ohlmann
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
24
25module curv_gygi_oct_m
27 use debug_oct_m
28 use global_oct_m
29 use, intrinsic :: iso_fortran_env
31 use math_oct_m
34 use parser_oct_m
37 use unit_oct_m
39
40 implicit none
41
42 private
43 public :: &
47
48 type, extends(coordinate_system_t) :: curv_gygi_t
49 private
50 real(real64), public :: A
51 real(real64), public :: alpha
52 real(real64), public :: beta
53 real(real64), allocatable :: pos(:, :)
54 integer :: npos
55 type(root_solver_t) :: rs
56 contains
57 procedure :: to_cartesian => curv_gygi_to_cartesian
58 procedure :: from_cartesian => curv_gygi_from_cartesian
59 procedure :: write_info => curv_gygi_write_info
60 procedure :: surface_element => curv_gygi_surface_element
61 procedure :: jacobian => curv_gygi_jacobian
62 procedure :: jacobian_inverse => curv_gygi_jacobian_inverse
63 procedure :: trace_hessian => curv_gygi_trace_hessian
64 final :: curv_gygi_finalize
65 end type curv_gygi_t
66
67 interface curv_gygi_t
68 procedure curv_gygi_constructor
69 end interface curv_gygi_t
70
71 ! Auxiliary variables for the root solver.
72 class(curv_gygi_t), pointer :: gygi_p
73 type(curv_gygi_t), target :: gygi_global
74 real(real64), allocatable :: chi_p(:)
75
76contains
77
78 ! ---------------------------------------------------------
79 function curv_gygi_constructor(namespace, dim, npos, pos) result(gygi)
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
85
86 push_sub(curv_gygi_constructor)
87
88 safe_allocate(gygi)
89
90 gygi%dim = dim
91 gygi%local_basis = .true.
92 gygi%orthogonal = .true. ! This needs to be checked.
93
94 gygi%npos = npos
95 safe_allocate(gygi%pos(1:dim, 1:gygi%npos))
96 gygi%pos(:, :) = pos(:, :)
97
98 !%Variable CurvGygiA
99 !%Type float
100 !%Default 0.5
101 !%Section Mesh::Curvilinear::Gygi
102 !%Description
103 !% The grid spacing is reduced locally around each atom, and the reduction is
104 !% given by 1/(1+<i>A</i>), where <i>A</i> is specified by this variable. So, if
105 !% <i>A</i>=1/2 (the default), the grid spacing is reduced to two thirds = 1/(1+1/2).
106 !% [This is the <math>A_{\alpha}</math> variable in Eq. 2 of F. Gygi and G. Galli, <i>Phys.
107 !% Rev. B</i> <b>52</b>, R2229 (1995)]. It must be larger than zero.
108 !% Values larger than two lead to errors because the coordinate transformation is not
109 !% bijective anymore For larger values of A.
110 !%End
111 call parse_variable(namespace, 'CurvGygiA', m_half, gygi%A)
112
113 !%Variable CurvGygiAlpha
114 !%Type float
115 !%Default 2.0 a.u.
116 !%Section Mesh::Curvilinear::Gygi
117 !%Description
118 !% This number determines the region over which the grid is enhanced (range of
119 !% enhancement of the resolution). That is, the grid is enhanced on a sphere
120 !% around each atom, whose radius is given by this variable. [This is the <math>a_{\alpha}</math>
121 !% variable in Eq. 2 of F. Gygi and G. Galli, <i>Phys. Rev. B</i> <b>52</b>, R2229 (1995)].
122 !% It must be larger than zero.
123 !%End
124
125 call parse_variable(namespace, 'CurvGygiAlpha', m_two, gygi%alpha, units_inp%length)
126 !%Variable CurvGygiBeta
127 !%Type float
128 !%Default 4.0 a.u.
129 !%Section Mesh::Curvilinear::Gygi
130 !%Description
131 !% This number determines the distance over which Euclidean coordinates are
132 !% recovered. [This is the <math>b_{\alpha}</math> variable in Eq. 2 of F. Gygi and G. Galli,
133 !% <i>Phys. Rev. B</i> <b>52</b>, R2229 (1995)]. It must be larger than zero.
134 !%End
135 call parse_variable(namespace, 'CurvGygiBeta', m_four, gygi%beta, units_inp%length)
136
137 if (gygi%a <= m_zero) call messages_input_error(namespace, 'CurvGygiA')
138 if (gygi%a > m_two) call messages_input_error(namespace, 'CurvGygiA')
139 if (gygi%alpha <= m_zero) call messages_input_error(namespace, 'CurvGygiAlpha')
140 if (gygi%beta <= m_zero) call messages_input_error(namespace, 'CurvGygiBeta')
142 gygi%min_mesh_scaling_product = (m_one / (m_one + gygi%A))**gygi%dim
144 ! initialize root solver
145 call root_solver_init(gygi%rs, namespace, dim, solver_type = root_newton, maxiter = 500, abs_tolerance = 1.0e-10_real64)
149
150 ! ---------------------------------------------------------
151 subroutine curv_gygi_copy(this_out, this_in)
152 type(curv_gygi_t), intent(inout) :: this_out
153 type(curv_gygi_t), intent(in) :: this_in
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
166 call root_solver_init(this_out%rs, global_namespace, this_in%dim, solver_type = root_newton, &
167 maxiter = 500, abs_tolerance = 1.0e-10_real64)
168
169 pop_sub(curv_gygi_copy)
170 end subroutine curv_gygi_copy
171
172 ! ---------------------------------------------------------
173 subroutine curv_gygi_finalize(this)
174 type(curv_gygi_t), intent(inout) :: this
175
176 push_sub(curv_gygi_finalize)
177
178 safe_deallocate_a(this%pos)
179
180 pop_sub(curv_gygi_finalize)
181 end subroutine curv_gygi_finalize
182
183 ! ---------------------------------------------------------
184 real(real64) pure function gygi_f(this, r)
185 class(curv_gygi_t), intent(in) :: this
186 real(real64), intent(in) :: r
187
188 ! no PUSH_SUB, called too often
189
190 if (r < 1.0e-4_real64) then
191 ! below 1e-4, the relative deviation of this Taylor expansion from the function is less than 1e-15
192 gygi_f = this%A * (m_one + (-m_one/(m_three*this%alpha**2) - m_one/this%beta**2)*r**2)
193 else if (-(r/this%beta)**2 <= m_min_exp_arg) then
194 gygi_f = m_zero
195 else
196 gygi_f = this%A * this%alpha/r * tanh(r/this%alpha) * exp(-(r/this%beta)**2)
197 end if
198 end function gygi_f
199
200 ! ---------------------------------------------------------
201 ! absorb division by r into the function to get the correct expansion
202 real(real64) pure function gygi_dfdr_over_r(this, r)
203 class(curv_gygi_t), intent(in) :: this
204 real(real64), intent(in) :: r
205
206 ! no PUSH_SUB, called too often
207
208 if (r < 3.0e-3_real64) then
209 ! below 3e-3, the expansion is closer than 1e-13 to the function and does not suffer from
210 ! rounding errors
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
215 gygi_dfdr_over_r = m_zero
216 else
217 gygi_dfdr_over_r = this%A * (this%beta**2*r/cosh(r/this%alpha)**2 - &
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)
220 end if
221 end function gygi_dfdr_over_r
222
223 ! ---------------------------------------------------------
224 ! only a combination of the second derivative is needed, namely
225 ! (d2f/d2r - 1/r df/dr)/r**2
226 real(real64) pure function gygi_d2fdr2_combination(this, r)
227 class(curv_gygi_t), intent(in) :: this
228 real(real64), intent(in) :: r
229
230 ! no PUSH_SUB, called too often
231
232 if (r < 3.0e-3_real64) then
233 ! below 3e-3, the expansion is closer than 1e-5 to the function and does not suffer from
234 ! rounding errors
235 gygi_d2fdr2_combination = this%A * ( &
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 &
240 )
241 else if (-(r/this%beta)**2 <= m_min_exp_arg) then
243 else
244 gygi_d2fdr2_combination = this%A * exp(-(r/this%beta)**2)*( &
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
248 end if
249 end function gygi_d2fdr2_combination
250
251
252 ! ---------------------------------------------------------
253 function curv_gygi_to_cartesian(this, chi) result(xx)
254 class(curv_gygi_t), target, intent(in) :: this
255 real(real64), intent(in) :: chi(:)
256 real(real64) :: xx(1:this%dim), xx_start(1:this%dim)
257
258 logical :: conv
259 integer :: i_conv, n_conv
260
261 ! no PUSH_SUB, called too often
262
263 gygi_p => this
264 safe_allocate(chi_p(1:this%dim))
265 chi_p(:) = chi(:)
267 call droot_solver_run(this%rs, getf, xx, conv, startval = chi)
268 nullify(gygi_p)
269
270 if (.not. conv) then
271 call curv_gygi_copy(gygi_global, this)
273 ! try to converge with decreasing A
274 n_conv = 8
275 do i_conv = 1, n_conv
276 gygi_p%A = gygi_p%A / m_two
277 call droot_solver_run(gygi_p%rs, getf, xx, conv, startval = chi)
278 if (conv) then
279 exit
280 end if
281 end do
282 if (conv) then
283 ! increase A again and take previous xx as starting point
284 n_conv = i_conv
285 do i_conv = n_conv, 1, -1
286 gygi_p%A = gygi_p%A * m_two
287 xx_start = xx
288 call droot_solver_run(gygi_p%rs, getf, xx, conv, startval = xx_start)
289 end do
290 end if
291 nullify(gygi_p)
293 end if
294
295 safe_deallocate_a(chi_p)
296
297 if (.not. conv) then
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)
304 end if
305
306 end function curv_gygi_to_cartesian
307
308 ! ---------------------------------------------------------
309 pure function curv_gygi_from_cartesian(this, xx) result(chi)
310 class(curv_gygi_t), target, intent(in) :: this
311 real(real64), intent(in) :: xx(:)
312 real(real64) :: chi(1:this%dim)
313
314 integer :: i, ia
315 real(real64) :: diff(this%dim), f
316
317 ! no PUSH_SUB, called too often
318
319 chi(1:this%dim) = xx(1:this%dim)
320 do ia = 1, this%npos
321 diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
322 f = gygi_f(this, norm2(diff))
323 do i = 1, this%dim
324 chi(i) = chi(i) + diff(i) * f
325 end do
326 end do
327
328 end function curv_gygi_from_cartesian
329
330 ! ---------------------------------------------------------
331 subroutine curv_gygi_write_info(this, iunit, namespace)
332 class(curv_gygi_t), intent(in) :: this
333 integer, optional, intent(in) :: iunit
334 type(namespace_t), optional, intent(in) :: namespace
335
336 push_sub(curv_gygi_write_info)
337
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)
347 pop_sub(curv_gygi_write_info)
348 end subroutine curv_gygi_write_info
349
350 ! ---------------------------------------------------------
351 real(real64) function curv_gygi_surface_element(this, idir) result(ds)
352 class(curv_gygi_t), intent(in) :: this
353 integer, intent(in) :: idir
354
355 ds = m_zero
356 message(1) = 'Surface element with gygi curvilinear coordinates not implemented'
357 call messages_fatal(1)
358
359 end function curv_gygi_surface_element
360
361 ! ---------------------------------------------------------
362 function curv_gygi_jacobian(this, chi) result(jacobian)
363 class(curv_gygi_t), intent(in) :: this
364 real(real64), intent(in) :: chi(:)
365 real(real64) :: jacobian(1:this%dim, 1:this%dim)
366
367 jacobian = this%jacobian_inverse(chi)
368 call lalg_inverse(this%dim, jacobian, "dir")
369 end function curv_gygi_jacobian
370
371 ! ---------------------------------------------------------
372 function curv_gygi_jacobian_inverse(this, chi) result(jacobian_inverse)
373 class(curv_gygi_t), intent(in) :: this
374 real(real64), intent(in) :: chi(:)
375 real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
376
377 real(real64) :: xx(1:this%dim)
378
379 ! no PUSH_SUB, called too often
380
381 xx(:) = this%to_cartesian(chi)
382 jacobian_inverse = curv_gygi_jacobian_inverse_cartesian(this, xx)
383
384 end function curv_gygi_jacobian_inverse
385
386 ! ---------------------------------------------------------
387 function curv_gygi_jacobian_inverse_cartesian(this, xx) result(jacobian_inverse)
388 class(curv_gygi_t), intent(in) :: this
389 real(real64), intent(in) :: xx(:)
390 real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
391
392 integer :: i, ix, iy
393 real(real64) :: r, diff(1:this%dim), dfdr_over_r, f
394
395 ! no PUSH_SUB, called too often
396
397 jacobian_inverse(1:this%dim, 1:this%dim) = diagonal_matrix(this%dim, m_one)
398
399 do i = 1, this%npos
400 diff = xx(1:this%dim) - this%pos(1:this%dim, i)
401 r = norm2(diff)
403 f = gygi_f(this, r)
404 dfdr_over_r = gygi_dfdr_over_r(this, r)
405
406 do ix = 1, this%dim
407 jacobian_inverse(ix, ix) = jacobian_inverse(ix, ix) + f
408 do iy = 1, this%dim
409 jacobian_inverse(ix, iy) = jacobian_inverse(ix, iy) + diff(ix)*diff(iy)*dfdr_over_r
410 end do
411 end do
412 end do
413
415
416
417 ! ---------------------------------------------------------
418 pure subroutine curv_gygi_hessian(this, xx, hessian, natoms)
419 class(curv_gygi_t), intent(in) :: this
420 real(real64), intent(in) :: xx(:)
421 real(real64), intent(out) :: hessian(:, :, :)
422 integer, optional, intent(in) :: natoms
423
424 integer :: ia, i, j, k, natoms_
425 real(real64) :: r, diff(this%dim), dfdr_over_r, d2fdr2_combination
426
427 ! no PUSH_SUB, called too often
428
429 hessian(1:this%dim, 1:this%dim, 1:this%dim) = m_zero
430
431 natoms_ = this%npos
432 if (present(natoms)) natoms_ = natoms
433
434 do ia = 1, natoms_
435 diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
436 r = norm2(diff)
437
438 dfdr_over_r = gygi_dfdr_over_r(this, r)
439 d2fdr2_combination = gygi_d2fdr2_combination(this, r)
440
441 do i = 1, this%dim
442 do j = 1, this%dim
443 do k = 1, this%dim
444 if (i == j) then
445 hessian(i, j, k) = hessian(i, j, k) + diff(k)*dfdr_over_r
446 end if
447 if (i == k) then
448 hessian(i, j, k) = hessian(i, j, k) + diff(j)*dfdr_over_r
449 end if
450 if (j == k) then
451 hessian(i, j, k) = hessian(i, j, k) + diff(i)*dfdr_over_r
452 end if
453 hessian(i, j, k) = hessian(i, j, k) + diff(i)*diff(j)*diff(k)*d2fdr2_combination
454 end do
455 end do
456 end do
457 end do
458
459 end subroutine curv_gygi_hessian
460
461 ! ---------------------------------------------------------
462 function curv_gygi_trace_hessian(this, chi) result(trace_hessian)
463 class(curv_gygi_t), intent(in) :: this
464 real(real64), intent(in) :: chi(:)
465 real(real64) :: trace_hessian(1:this%dim)
466
467 integer :: ia, i, k
468 real(real64) :: r, xx(this%dim), diff(this%dim), dfdr_over_r, d2fdr2_combination
469
470 ! no PUSH_SUB, called too often
471
472 xx(:) = this%to_cartesian(chi)
473 trace_hessian(:) = m_zero
474
475 do ia = 1, this%npos
476 diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
477 r = norm2(diff)
478
479 dfdr_over_r = gygi_dfdr_over_r(this, r)
480 d2fdr2_combination = gygi_d2fdr2_combination(this, r)
481
482 do i = 1, this%dim
483 do k = 1, this%dim
484 if (i == k) then
485 trace_hessian(i) = trace_hessian(i) + m_two*diff(k)*dfdr_over_r
486 end if
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
489 end do
490 end do
491 end do
492
493 end function curv_gygi_trace_hessian
494
495 ! ---------------------------------------------------------
496 subroutine getf(y, f, jf)
497 real(real64), intent(in) :: y(:)
498 real(real64), intent(out) :: f(:), jf(:, :)
499
500 ! no PUSH_SUB, called too often
501
503 f = gygi_p%from_cartesian(y)
504 f(1:gygi_p%dim) = f(1:gygi_p%dim) - chi_p(1:gygi_p%dim)
505
506 end subroutine getf
507
508end module curv_gygi_oct_m
509
510!! Local Variables:
511!! mode: f90
512!! coding: utf-8
513!! End:
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)...
Definition: curv_gygi.F90:118
type(curv_gygi_t), target gygi_global
Definition: curv_gygi.F90:166
subroutine curv_gygi_write_info(this, iunit, namespace)
Definition: curv_gygi.F90:425
real(real64), dimension(:), allocatable chi_p
Definition: curv_gygi.F90:167
real(real64) pure function gygi_f(this, r)
Definition: curv_gygi.F90:278
real(real64) pure function gygi_d2fdr2_combination(this, r)
Definition: curv_gygi.F90:320
real(real64) pure function gygi_dfdr_over_r(this, r)
Definition: curv_gygi.F90:296
pure subroutine, public curv_gygi_hessian(this, xx, hessian, natoms)
Definition: curv_gygi.F90:512
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_gygi_jacobian_inverse_cartesian(this, xx)
Definition: curv_gygi.F90:481
real(real64) function, dimension(1:this%dim) curv_gygi_trace_hessian(this, chi)
Definition: curv_gygi.F90:556
pure real(real64) function, dimension(1:this%dim) curv_gygi_from_cartesian(this, xx)
Definition: curv_gygi.F90:403
subroutine, public curv_gygi_copy(this_out, this_in)
Definition: curv_gygi.F90:245
real(real64) function, dimension(1:this%dim) curv_gygi_to_cartesian(this, chi)
Definition: curv_gygi.F90:347
subroutine getf(y, f, jf)
Definition: curv_gygi.F90:590
real(real64) function curv_gygi_surface_element(this, idir)
Definition: curv_gygi.F90:445
class(curv_gygi_t) function, pointer curv_gygi_constructor(namespace, dim, npos, pos)
Definition: curv_gygi.F90:173
class(curv_gygi_t), pointer gygi_p
Definition: curv_gygi.F90:165
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_gygi_jacobian(this, chi)
Definition: curv_gygi.F90:456
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_gygi_jacobian_inverse(this, chi)
Definition: curv_gygi.F90:466
subroutine curv_gygi_finalize(this)
Definition: curv_gygi.F90:267
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_min_exp_arg
Definition: global.F90:207
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
type(namespace_t), public global_namespace
Definition: namespace.F90:132
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.
Definition: unit.F90:132
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
int true(void)