From 776215457376ce4c1e1f965fdd3ee2d3f0cefb52 Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Sat, 29 Oct 2022 12:22:34 +0900 Subject: [PATCH] Implement all things for fractional calculus --- MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.c | 45 ++++ MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.f90 | 192 ++++++++++++++++++ MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.h | 9 + .../COMPILATION_METHOD.txt | 18 +- MultiWFN/Multiwfn_3.8_dev_src_Linux/Makefile | 33 ++- .../Multiwfn_3.8_dev_src_Linux/Multiwfn.f90 | 4 + MultiWFN/Multiwfn_3.8_dev_src_Linux/no2F2.c | 11 + 7 files changed, 307 insertions(+), 5 deletions(-) create mode 100644 MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.c create mode 100644 MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.f90 create mode 100644 MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.h create mode 100644 MultiWFN/Multiwfn_3.8_dev_src_Linux/no2F2.c diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.c b/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.c new file mode 100644 index 0000000..f8c99aa --- /dev/null +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.c @@ -0,0 +1,45 @@ +#include "acb_hypgeom.h" + +#include "2F2.h" + +double hyp2F2_prec(double a1, double a2, double b1, double b2, double z, + slong prec) { + const slong p = 2, q = 2; + const int regularized = 0; + arb_t res_r; + acb_t res_ri, _z; + acb_ptr a, b; + arf_t resf; + a = _acb_vec_init(p); + b = _acb_vec_init(q); + arb_init(res_r); + acb_init(res_ri); + acb_init(_z); + arf_init(resf); + acb_set_d(a + 0, a1); + acb_set_d(a + 1, a2); + acb_set_d(b + 0, b1); + acb_set_d(b + 1, b2); + acb_set_d(_z, z); + acb_hypgeom_pfq(res_ri, a, p, b, q, _z, regularized, prec); + if (acb_is_real(res_ri)) { + acb_get_real(res_r, res_ri); + } else { + return 0.0; + } + arb_get_abs_ubound_arf(resf, res_r, prec); + double hyp2F2 = arf_get_d(resf, ARF_RND_DOWN); + return hyp2F2; +} + +double hyp2F2(double a1, double a2, double b1, double b2, double z) { + slong prec = 1000; + double hyp2F2 = 0.0; + // for our purpose, 0 < hyp2F2 <= 1 + while (hyp2F2 == 0.0 || hyp2F2 > 1.0) { + hyp2F2 = hyp2F2_prec(a1, a2, b1, b2, z, prec); + if(hyp2F2 == 0.0 && prec > 100000) break; + prec *= 10; + } + return hyp2F2; +} diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.f90 b/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.f90 new file mode 100644 index 0000000..155dae2 --- /dev/null +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.f90 @@ -0,0 +1,192 @@ +module mod_2F2 + use util, only: to_lower => struc2lc + implicit none + private + integer, parameter :: GENERAL_EXACT_EVALUATION = 1, & + INTERPOLATION_EVALUATION = 2, & + SPECIAL_EXACT_EVALUATION = 3 + ! normally, it should be 0 (fractional integrals) or 1 (fractional derivatives) + integer :: p = 2000 + ! order of fractional derivative + real(8) :: alpha_ + ! used algorithm + integer :: algorithm_ = -1 + ! interface to C function; see 2F2.h + ! note, this routine returns 0 < hyp2F2 <= 1 + ! z should be <= 0 + interface + real(8) function hyp2F2(a, b, c, d, z) bind(C, name="hyp2F2") + real(8), intent(in), value :: a, b, c, d, z + end function hyp2F2 + end interface + ! internal global arrays for interpolation algorithm + integer, parameter :: n_2F2 = 101, k_2F2 = 10 + real(8) eprr_2F2(0:n_2F2 - 1) + real(8) f_2F2(0:n_2F2 - 1, 0:8), Bs_2F2(0:n_2F2 - 1, 0:8) + real(8) t_2F2(n_2F2 + k_2F2) + ! + public :: set_alpha_level, GTO_fractional_integral + public :: p +contains + subroutine set_alpha_level() + character(len=80) :: tempstr + do while (.true.) + write (*, *) "Type `get` for getting information about current fractional derivative parameters" + write (*, *) "Type `set_p` for updating fractional derivative parameters with specific integer order of internal derivatives" + write (*, *) "Type `set_auto` for updating fractional derivative parameters with automatic integer order of internal derivatives" + read (*, *) tempstr + call to_lower(tempstr) + select case(trim(tempstr)) + case("get") + call get_alpha_level() + return + case("set_p") + write (*, *) "Input fractional order of derivative (alpha); it should be <= 1" + read (*, *) alpha_ + write (*, *) "Input integer order of internal derivative (p); it can be 0 or 1" + read (*, *) p + if (alpha_ > p) then + write (*, *) "alpha can not be larger than p" + write (*, *) "" + cycle + end if + if (p /= 0 .and. p /= 1) then + write (*, *) "p can be only 0 or 1" + write (*, *) "" + cycle + end if + case("set_auto") + write (*, *) "Input fractional order of derivative (alpha); it should be <= 1" + read (*, *) alpha_ + p = ceiling(alpha_) + if (p < 0) p = 0 + write (*, '(A,I0)') "Automatically chosen integer order of internal derivative (p) is ", p + case default + write (*, *) "Try again..." + cycle + end select + exit + end do + write (*, *) "Type `exact` for enabling evaluation of general analytic form of fractional derivative" + write (*, *) "Type `approx` for interpolation of fractional derivative" + write (*, *) "Type `special` for enabling evaluation of specific analytic form of fractional derivative" + write (*, *) "For special, only several pairs are available:" + write (*, *) " p alpha" + write (*, *) " 1 1.0" + write (*, *) " 0 0.0" + do while (.true.) + read (*, *) tempstr + call to_lower(tempstr) + select case(trim(tempstr)) + case("exact") + algorithm_ = GENERAL_EXACT_EVALUATION + case("approx") + algorithm_ = INTERPOLATION_EVALUATION + call prepare_interpolation() + case("special") + algorithm_ = SPECIAL_EXACT_EVALUATION + case default + write (*, *) "Try again..." + cycle + end select + exit + end do + call get_alpha_level() + end subroutine set_alpha_level + subroutine get_alpha_level() + write (*, *) "Fractional derivative parameters:" + write (*, *) "Alpha: ", alpha_ + write (*, '(A,I0,A,I0)') "Order of derivative is ", p, "; should be ", ceiling(alpha_) + if (algorithm_ == GENERAL_EXACT_EVALUATION) then + write (*, *) "Evaluation of fractional derivative using general analytic form" + else if (algorithm_ == INTERPOLATION_EVALUATION) then + write (*, *) "Evaluation of fractional derivative using interpolation" + else if (algorithm_ == SPECIAL_EXACT_EVALUATION) then + write (*, *) "Evaluation of fractional derivative using specific analytic form" + else + write (*, *) "Unknown algorithm" + end if + write (*, *) + end subroutine get_alpha_level + ! This routine evaluates the following expression: + ! \frac{1}{\Gamma(p-\alpha)} \int\limits_0^1 (1-t)^{p-\alpha-1} t^n exp(x t^2) dt + ! using general, approximational or special algorithms + ! n is integer, principal quantum number + ! x is real(8) <= 0 + real(8) function GTO_fractional_integral(n, x) + use bspline_sub_module, only: db1val + real(8), parameter :: one = 1e0_8, zero = 0e0_8 + integer, intent(in) :: n + real(8), intent(in) :: x + real(8) :: top, down + real(8) :: gtop, gdown + real(8) :: val2F2 + integer :: iout, inbvx + if (x > zero) then + write (*, *) "x is ", x + error stop "x greater than zero!" + end if + val2F2 = -one + GTO_fractional_integral = 0e0_8 + if (algorithm_ == GENERAL_EXACT_EVALUATION) then + top = real(n + 1, kind=8) + down = top + real(p, kind=8) - alpha_ + gtop = gamma(top) + gdown = gamma(down) + val2F2 = hyp2F2(top, top + one, down, down + one, x) + GTO_fractional_integral = gtop / gdown * val2F2 + else if (algorithm_ == INTERPOLATION_EVALUATION) then + top = real(n + 1, kind=8) + down = top + real(p, kind=8) - alpha_ + gtop = gamma(top) + gdown = gamma(down) + inbvx = 1 + ! normally, x shall be passed, but here some logic is broken since Bspline expects increasing values of x + call db1val(-x, 0, t_2F2, n_2F2, k_2F2, Bs_2F2(:, n), val2F2, iout, inbvx) + GTO_fractional_integral = gtop / gdown * val2F2 + else if (algorithm_ == SPECIAL_EXACT_EVALUATION) then + select case(p) + case(1) + if (alpha_ == one) then + GTO_fractional_integral = exp(x) + end if + case(0) + if (alpha_ == zero) then + GTO_fractional_integral = exp(x) + end if + end select + if(GTO_fractional_integral /= -one) return + write (*, *) "alpha: ", alpha_, " p: ", p + error stop "No specialization for given alpha and p" + else + write (*, '(A,I0)') "Algorithm internal value is ", algorithm_ + error stop "Algorithm is not implemented!" + end if + end function GTO_fractional_integral + ! this routine fills arrays for interpolation evaluation of 2F2 function + subroutine prepare_interpolation() + use bspline_sub_module, only: db1ink + real(8), parameter :: one = 1e0_8, zero = 0e0_8 + real(8) :: top, down + integer :: i, j + integer :: iout + do i = 0, n_2F2 - 1 + eprr_2F2(i) = exp(0.1_8 * real(i, 8)) - one + do j = 0, 8 + if (i > 0 .and. f_2F2(max(i - 1, 0), j) < 1e-40_8) then + f_2F2(i, j) = zero + else + top = real(j + 1, kind=8) + down = top + real(p, kind=8) - alpha_ + ! we evaluate negative eprr + f_2F2(i, j) = hyp2F2(top, top + one, down, down + one, -eprr_2F2(i)) + end if + end do + write(*,*) i, eprr_2F2(i), f_2F2(i, 0:4) + end do + do j = 0, 8 + ! here, eprr_2F2 is from 0 to +inf + call db1ink(eprr_2F2, n_2F2, f_2F2(:, j), k_2F2, 0, t_2F2, Bs_2F2(:, j), iout) + end do + end subroutine prepare_interpolation +end module mod_2F2 diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.h b/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.h new file mode 100644 index 0000000..bc2736b --- /dev/null +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/2F2.h @@ -0,0 +1,9 @@ +#ifdef __cplusplus +extern "C" { +#endif + +double hyp2F2(double a1, double a2, double b1, double b2, double z); + +#ifdef __cplusplus +} +#endif diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/COMPILATION_METHOD.txt b/MultiWFN/Multiwfn_3.8_dev_src_Linux/COMPILATION_METHOD.txt index 1b58c9f..b8cecf3 100644 --- a/MultiWFN/Multiwfn_3.8_dev_src_Linux/COMPILATION_METHOD.txt +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/COMPILATION_METHOD.txt @@ -27,4 +27,20 @@ If this is the first time you use Linux version of Multiwfn, do not forget to re ============ Compile Multiwfn without GUI and plotting functionalities ============ -If you encountered difficulties when compiling or running Multiwfn due to missing or incompatibility of some graphical library files, and in the meantime you do not need visualization functions of Multiwfn, you can compile Multiwfn without GUI supported, all functions irrelevant to GUI and plotting will still work normally. The compilation steps is simply running "make noGUI -j" instead of the "make -j" mentioned above. \ No newline at end of file +If you encountered difficulties when compiling or running Multiwfn due to missing or incompatibility of some graphical library files, and in the meantime you do not need visualization functions of Multiwfn, you can compile Multiwfn without GUI supported, all functions irrelevant to GUI and plotting will still work normally. The compilation steps is simply running "make noGUI -j" instead of the "make -j" mentioned above. + + +============ Compile Multiwfn with fractional derivatives support ============ + +1) Be sure that you have installed arb (https://arblib.org/) and flint (http://flintlib.org/) libraries + for Fedora (I suppose, for RHEL/CentOS the same command will work): + dnf install arb-devel + or + yum install arb-devel + for Ubuntu: + apt install libflint-arb-dev + +2) provide `WITH_FD=1` to your `make` command and `OS=Ubuntu` (For Ubuntu) or `OS=RHEL` (for Fedora, CentOS, RHEL) + So, on my PC with Fedora 35: `make noGUI WITH_FD=1 OS=RHEL -j 6` + on Ubuntu 22.04: `make noGUI WITH_FD=1 OS=Ubuntu -j 24` +3) If there is some problem with compilation, please, revise `Makefile` about libraries and include directories below line containing `ifeq ($(WITH_FD),1)` diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/Makefile b/MultiWFN/Multiwfn_3.8_dev_src_Linux/Makefile index 91c2378..0076eec 100644 --- a/MultiWFN/Multiwfn_3.8_dev_src_Linux/Makefile +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/Makefile @@ -4,9 +4,12 @@ OPT1 = -O2 -fopenmp -ffree-line-length-none -cpp #Options in the next line is for debugging purpose #OPTDBG = -O0 -qopenmp -diag-disable 8290,8291,6371 -threads -qopenmp-link=static -debug all -g -traceback -check all -fstack-protector -fpp -mkl -static-intel -LIB_GUI = ./dislin_d-11.0.a -lXm -lXt -lX11 -lGL #At least works for CentOS 7.x -LIB_noGUI = +LIB_base = +LIB_GUI = $(LIB_base) ./dislin_d-11.0.a -lXm -lXt -lX11 -lGL #At least works for CentOS 7.x +LIB_noGUI = $(LIB_base) +INCLUDE = FC = ifort +CC = icc EXE = Multiwfn EXE_noGUI = Multiwfn_noGUI LIBRETAPATH = ./libreta_hybrid @@ -16,7 +19,21 @@ DFTxclib.o edflib.o fparser.o fileIO.o spectrum.o DOS.o Multiwfn.o 0123dim.o LSB population.o orbcomp.o bondorder.o topology.o excittrans.o otherfunc.o \ otherfunc2.o otherfunc3.o O1.o surfana.o procgriddata.o AdNDP.o fuzzy.o CDA.o basin.o \ orbloc.o visweak.o EDA.o CDFT.o ETS_NOCV.o atmraddens.o NAONBO.o grid.o PBC.o hyper_polar.o deloc_aromat.o cp2kmate.o\ -minpack.o blockhrr_012345.o ean.o hrr_012345.o eanvrr_012345.o boysfunc.o naiveeri.o ryspoly.o +minpack.o blockhrr_012345.o ean.o hrr_012345.o eanvrr_012345.o boysfunc.o naiveeri.o ryspoly.o 2F2.f90.o + +ifeq ($(WITH_FD),1) + objects += 2F2.c.o +# for Ubuntu +ifeq ($(OS),Ubuntu) + LIB_base += -lflint -lflint-arb +# for Fedora, CentOS, RHEL +else ifeq ($(OS),RHEL) + INCLUDE += -I/usr/include/arb + LIB_base += -lflint -larb +endif +else + objects += no2F2.c.o +endif objects_noGUI = noGUI/dislin_d_empty.o @@ -74,7 +91,10 @@ plot.o : plot.f90 define.o util.o GUI.o : GUI.f90 define.o plot.o function.o $(FC) $(OPT) -c GUI.f90 -modules = define.o util.o function.o plot.o GUI.o libreta.o +2F2.f90.o : 2F2.f90 util.o Bspline.o + $(FC) $(OPT) -c 2F2.f90 -o 2F2.f90.o + +modules = define.o util.o function.o plot.o GUI.o libreta.o 2F2.f90.o #Library or adpated third-part codes @@ -100,6 +120,11 @@ minpack.o : minpack.f90 fparser.o : fparser.f90 $(FC) $(OPT) -c fparser.f90 +2F2.c.o : 2F2.c + $(CC) $(INCLUDE) -c 2F2.c -o 2F2.c.o + +no2F2.c.o : no2F2.c + $(CC) -c no2F2.c -o no2F2.c.o #Others diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/Multiwfn.f90 b/MultiWFN/Multiwfn_3.8_dev_src_Linux/Multiwfn.f90 index ea49251..9c8b635 100644 --- a/MultiWFN/Multiwfn_3.8_dev_src_Linux/Multiwfn.f90 +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/Multiwfn.f90 @@ -2,6 +2,7 @@ program multiwfn use defvar use util use GUI +use mod_2F2, only: set_alpha_level !use libreta !use function implicit real*8 (a-h,o-z) @@ -586,6 +587,7 @@ do while(.true.) !Main loop write(*,*) "15 Make orbitals equivalent to basis functions" write(*,*) "16 Define one or two fragments for special purpose" write(*,*) "17 Generate promolecular wavefunction by calculating and combining atomic ones" + write(*,*) "26 Setup fractional calculus" write(*,*) "90 Calculate nuclear attractive energy between a fragment and an orbital" write(*,*) "91 Exchange orbital energies and occupations" write(*,*) "92 Calculate result of various kinetic energy functionals" @@ -736,6 +738,8 @@ do while(.true.) !Main loop end do else if (i==17) then call generate_promolwfn + else if (i==26) then + call set_alpha_level() else if (i==90) then call attene_orb_fragnuc else if (i==91) then diff --git a/MultiWFN/Multiwfn_3.8_dev_src_Linux/no2F2.c b/MultiWFN/Multiwfn_3.8_dev_src_Linux/no2F2.c new file mode 100644 index 0000000..f7d8605 --- /dev/null +++ b/MultiWFN/Multiwfn_3.8_dev_src_Linux/no2F2.c @@ -0,0 +1,11 @@ +#include +#include + +#include "2F2.h" + +double hyp2F2(double a1, double a2, double b1, double b2, double z) { + fprintf(stderr, "Fractional derivatives are not supported!"); + fprintf(stderr, "Please, recompile Multiwfn with its support!"); + exit(1); + return 0.0; +} -- 2.34.1

map