cp2k/tools/Fun2D/functions.f90
Joost VandeVondele 3ff9d0a356 also added routines for the yukawa potential, and cleaned up to code in order to make
adding new functions easier.


svn-origin-rev: 8076
2008-12-17 11:44:07 +00:00

85 lines
2.5 KiB
Fortran

!
! wrap all possible functions together, they need unique ids
!
MODULE functions
USE mpfr
USE mpfr_ops
USE mpfr_cutoff_gamma
USE mpfr_yukawa
USE function_types
IMPLICIT NONE
INTEGER, SAVE :: functionid=-1
INTEGER, SAVE :: should_output=-1
CONTAINS
! evaluate the function, + nderiv other functions in the point X1,X2 in a reduced domain
! as a first step, this transforms X1,X2 into the natural variables of the original domain (e.g. R,T)
! write additionally to an open unit (should_output>0) connected to a file, for later checking the quality of the interpolation
SUBROUTINE f(res,nderiv,x1,x2)
INTEGER :: nderiv
TYPE(mpfr_type) :: x1,x2
TYPE(mpfr_type) :: res(0:nderiv)
TYPE(mpfr_type) :: T,r,upper,lower,zero
TYPE(mpfr_type) :: dummy(0:21)
INTEGER :: i
zero=.CONVERT."0"
SELECT CASE(functionid)
CASE(fun_trunc_coulomb_farfield,fun_trunc_coulomb_nearfield)
! 1 = farfield (R>11)
! 2 = nearfield (R<11)
SELECT CASE(functionid)
CASE (fun_trunc_coulomb_farfield)
! if R=+Infinity the result is zero
IF (mpfr_cmp(x2,zero)<=0) THEN
res=.CONVERT."0"
return
ENDIF
r=11/x2
upper=r**2 + 11*r + 50
lower=r**2 - 11*r
CASE (fun_trunc_coulomb_nearfield)
R=x2*11
upper=r**2 + 11*r + 50
lower=.CONVERT."0"
CASE DEFAULT
STOP "Function ID not implemented"
END SELECT
t=lower+x1*(upper-lower)
!t is zero, return the limiting expansion
IF (mpfr_cmp(t,zero)<=0) THEN
CALL cutoff_gamma_T0(21,R,dummy)
ELSE
CALL cutoff_gamma(21,t,r,dummy)
END IF
res(0:nderiv)=dummy(0:nderiv)
CASE(fun_yukawa)
! deal with infinite T,R locally
IF((X1==zero) .OR. (X2==zero)) THEN
res=zero
return
ELSE
T=(1/X1-1)**2
R=(1/X2-1)**2
CALL yukawa_gn_all(nderiv,T,R,dummy)
ENDIF
res(0:nderiv)=dummy(0:nderiv)
CASE DEFAULT
STOP "Function ID not implemented"
END SELECT
IF (should_output>0) THEN
!$OMP CRITICAL
WRITE(should_output,*) REAL(R),REAL(T),REAL(res(0:nderiv))
!$OMP END CRITICAL
ENDIF
END SUBROUTINE f
END MODULE functions