Update the keys
This commit is contained in:
jean paul nshuti
parent
d9f4329a4c
commit
eaa8c4fb52
BIN
bin/genetic
BIN
bin/genetic
Binary file not shown.
BIN
obj/adia.o
BIN
obj/adia.o
Binary file not shown.
BIN
obj/adia_mod.mod
BIN
obj/adia_mod.mod
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
obj/diab_mod.mod
BIN
obj/diab_mod.mod
Binary file not shown.
BIN
obj/keys.o
BIN
obj/keys.o
Binary file not shown.
BIN
obj/keys_mod.mod
BIN
obj/keys_mod.mod
Binary file not shown.
Binary file not shown.
BIN
obj/model.o
BIN
obj/model.o
Binary file not shown.
BIN
obj/write.o
BIN
obj/write.o
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
|
@ -1,673 +0,0 @@
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
! % SUBROUTINE CTRANS(...)
|
|
||||||
! %
|
|
||||||
! % M. Vossel 21.03.2023
|
|
||||||
! %
|
|
||||||
! % Routine to transform symmetryinput coordinates to symmetrized
|
|
||||||
! % coordinates. Distances Are discribet by Morse coordinates or
|
|
||||||
! % TMC depending on Set Parameters in the Genetic Input.
|
|
||||||
! %
|
|
||||||
! % input variables
|
|
||||||
! % q:
|
|
||||||
! % q(1): H1x
|
|
||||||
! % q(2): y
|
|
||||||
! % q(3): z
|
|
||||||
! % q(4): H2x
|
|
||||||
! % q(5): y
|
|
||||||
! % q(6): z
|
|
||||||
! % q(7): H3x
|
|
||||||
! % q(8): y
|
|
||||||
! % q(9): z
|
|
||||||
!
|
|
||||||
!
|
|
||||||
!
|
|
||||||
! % Internal variables:
|
|
||||||
! % t: primitive coordinates (double[qn])
|
|
||||||
! % t(1):
|
|
||||||
! % t(2):
|
|
||||||
! % t(3):
|
|
||||||
! % t(4):
|
|
||||||
! % t(5):
|
|
||||||
! % t(6):
|
|
||||||
! % t(7):
|
|
||||||
! % t(8):
|
|
||||||
! % t(9):
|
|
||||||
! % t: dummy (double[qn])
|
|
||||||
! % p: parameter vector
|
|
||||||
! % npar: length of parameter vector
|
|
||||||
! %
|
|
||||||
! % Output variables
|
|
||||||
! % s: symmetrized coordinates (double[qn])
|
|
||||||
! % s(1): CH-symetric streatch
|
|
||||||
! % s(2): CH-asymetric streatch-ex
|
|
||||||
! % s(3): CH-asymetric streatch-ey
|
|
||||||
! % s(4): CH-bend-ex
|
|
||||||
! % s(5): CH-bend-ey
|
|
||||||
! % s(6): CH-umbrella
|
|
||||||
! % s(7): CH-umbrella**2
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
module ctrans_mod
|
|
||||||
use accuracy_constants, only: dp, idp
|
|
||||||
implicit none
|
|
||||||
! precalculate pi, 2*pi and angle to radian conversion
|
|
||||||
real(dp), parameter :: pi = 4.0_dp*datan(1.0_dp)
|
|
||||||
real(dp), parameter :: pi2 = 2.0_dp*pi
|
|
||||||
real(dp), parameter :: ang2rad = pi/180.0_dp
|
|
||||||
! precalculate roots
|
|
||||||
real(dp), parameter:: sq2 = 1.0_dp/dsqrt(2.0_dp)
|
|
||||||
real(dp), parameter:: sq3 = 1.0_dp/dsqrt(3.0_dp)
|
|
||||||
real(dp), parameter:: sq6 = 1.0_dp/dsqrt(6.0_dp)
|
|
||||||
! change distances for equilibrium
|
|
||||||
!real(dp), parameter :: dchequi = 1.02289024_dp
|
|
||||||
real(dp), parameter :: dchequi = 2.344419_dp ! NO3
|
|
||||||
!real(dp), parameter :: dchequi = 2.34451900_dp
|
|
||||||
|
|
||||||
! see changes
|
|
||||||
|
|
||||||
contains
|
|
||||||
subroutine ctrans(q, x1,x2, invariants)
|
|
||||||
use dim_parameter, only: qn
|
|
||||||
integer(idp) k !running indices
|
|
||||||
real(dp), intent(in) :: q(qn) !given coordinates
|
|
||||||
real(dp), intent(out) :: x1(qn) !output coordinates symmetry adapted and scaled
|
|
||||||
real(dp), intent(out) :: x2(qn) !output coordinates symmetry adapted but not scaled
|
|
||||||
! ANN Variables
|
|
||||||
real(dp), optional, intent(out) :: invariants(:)
|
|
||||||
real(dp) :: s(qn),t(qn)
|
|
||||||
! kartesian coordianates copy from MeF+ so substitute c by n and removed f
|
|
||||||
real(dp) ch1(3), ch2(3), ch3(3), c_atom(3)
|
|
||||||
real(dp) nh1(3), nh2(3), nh3(3)
|
|
||||||
real(dp) zaxis(3), xaxis(3), yaxis(3)
|
|
||||||
real(dp) ph1(3), ph2(3), ph3(3)
|
|
||||||
! primitive coordinates
|
|
||||||
real(dp) dch1, dch2, dch3 !nh-distances
|
|
||||||
real(dp) umb !Umbrella Angle from xy-plane
|
|
||||||
|
|
||||||
! Symmetry coordinates
|
|
||||||
real(dp) aR !a1-modes H-Dist.,
|
|
||||||
real(dp) exR, exAng !ex components H-Dist., H-Ang.
|
|
||||||
real(dp) eyR, eyAng !ey components H-Dist., H-Ang.
|
|
||||||
! debugging
|
|
||||||
logical, parameter :: dbg = .false.
|
|
||||||
|
|
||||||
! initialize coordinate vectors
|
|
||||||
s = 0.0_dp
|
|
||||||
t = 0.0_dp
|
|
||||||
|
|
||||||
! write kartesian coords for readability
|
|
||||||
c_atom(1:3) = q(1:3)
|
|
||||||
do k = 1, 3
|
|
||||||
ch1(k) = q(k + 3)
|
|
||||||
ch2(k) = q(k + 6)
|
|
||||||
ch3(k) = q(k + 9)
|
|
||||||
end do
|
|
||||||
|
|
||||||
! construct z-axis
|
|
||||||
nh1 = normalized(ch1)
|
|
||||||
nh2 = normalized(ch2)
|
|
||||||
nh3 = normalized(ch3)
|
|
||||||
zaxis = create_plane(nh1, nh2, nh3)
|
|
||||||
|
|
||||||
! calculate bonding distance
|
|
||||||
dch1 = norm(ch1)
|
|
||||||
dch2 = norm(ch2)
|
|
||||||
dch3 = norm(ch3)
|
|
||||||
|
|
||||||
! construct symmertic and antisymmetric strech
|
|
||||||
aR = symmetrize(dch1 - dchequi, dch2 - dchequi, dch3 - dchequi, 'a')
|
|
||||||
exR = symmetrize(dch1, dch2, dch3, 'x')
|
|
||||||
eyR = symmetrize(dch1, dch2, dch3, 'y')
|
|
||||||
|
|
||||||
! construc x-axis and y axis
|
|
||||||
ph1 = normalized(project_point_into_plane(nh1, zaxis, c_atom))
|
|
||||||
xaxis = normalized(ph1)
|
|
||||||
yaxis = xproduct(zaxis, xaxis) ! right hand side koordinates
|
|
||||||
|
|
||||||
! project H atoms into C plane
|
|
||||||
ph2 = normalized(project_point_into_plane(nh2, zaxis, c_atom))
|
|
||||||
ph3 = normalized(project_point_into_plane(nh3, zaxis, c_atom))
|
|
||||||
|
|
||||||
call construct_HBend(exAng, eyAng, ph1, ph2, ph3, xaxis, yaxis)
|
|
||||||
umb = construct_umbrella(nh1, nh2, nh3, zaxis)
|
|
||||||
|
|
||||||
! set symmetry coordinates and even powers of umbrella
|
|
||||||
!s(1) = dch1- dchequi !aR
|
|
||||||
!s(2) = dch2 - dchequi !exR
|
|
||||||
!s(3) = dch3 - dchequi !eyR
|
|
||||||
s(1) = aR
|
|
||||||
s(2) = exR
|
|
||||||
s(3) = eyR
|
|
||||||
s(4) = exAng
|
|
||||||
s(5) = eyAng
|
|
||||||
s(6) = umb
|
|
||||||
s(7) = umb**2
|
|
||||||
s(8) = 0
|
|
||||||
s(9) = 0
|
|
||||||
! pairwise distances as second coordinate set
|
|
||||||
t = 0._dp
|
|
||||||
call pair_distance(q, t(1:6))
|
|
||||||
|
|
||||||
if (dbg) write (6, '("sym coords s=",9f16.8)') s(1:qn)
|
|
||||||
if (dbg) write (6, '("sym coords t=",9f16.8)') t(1:qn)
|
|
||||||
if (present(invariants)) then
|
|
||||||
call get_invariants(s, invariants)
|
|
||||||
end if
|
|
||||||
|
|
||||||
! transform s and t to x1 and x2
|
|
||||||
x1(1:qn)=s(1:qn)
|
|
||||||
x1(5)=x1(5)
|
|
||||||
! set other x coordinate to zero other than strech
|
|
||||||
!X1(4:qn)=0.0d0
|
|
||||||
x2(1:qn)=t(1:qn)
|
|
||||||
end subroutine ctrans
|
|
||||||
|
|
||||||
subroutine pair_distance(q, r)
|
|
||||||
real(dp), intent(in) :: q(9)
|
|
||||||
real(dp), intent(out) :: r(6)
|
|
||||||
real(dp) :: atom(3, 4)
|
|
||||||
integer :: n, k, count
|
|
||||||
|
|
||||||
!atom order: H1 H2 H3 N
|
|
||||||
atom(:, 1:3) = reshape(q, [3, 3])
|
|
||||||
atom(:, 4) = (/0.0_dp, 0.0_dp, 0.0_dp/)
|
|
||||||
|
|
||||||
! disntace order 12 13 14 23 24 34
|
|
||||||
count = 0
|
|
||||||
do n = 1, size(atom, 2)
|
|
||||||
do k = n + 1, size(atom, 2)
|
|
||||||
count = count + 1
|
|
||||||
r(count) = sqrt(sum((atom(:, k) - atom(:, n))**2))
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end subroutine pair_distance
|
|
||||||
|
|
||||||
function morse_and_symmetrize(x,p,pst) result(s)
|
|
||||||
real(dp), intent(in),dimension(3) :: x
|
|
||||||
real(dp), intent(in),dimension(11) :: p
|
|
||||||
integer, intent(in),dimension(2) :: pst
|
|
||||||
integer :: k
|
|
||||||
real(dp), dimension(3) :: s
|
|
||||||
real(dp), dimension(3) :: t
|
|
||||||
|
|
||||||
! Morse transform
|
|
||||||
do k=1,3
|
|
||||||
t(k) = morse_transform(x(k), p, pst)
|
|
||||||
end do
|
|
||||||
s(1) = symmetrize(t(1), t(2), t(3), 'a')
|
|
||||||
s(2) = symmetrize(t(1), t(2), t(3), 'x')
|
|
||||||
s(3) = symmetrize(t(1), t(2), t(3), 'y')
|
|
||||||
end function morse_and_symmetrize
|
|
||||||
|
|
||||||
subroutine get_invariants(s, inv_out)
|
|
||||||
use dim_parameter, only: qn
|
|
||||||
use select_monom_mod, only: v_e_monom, v_ee_monom
|
|
||||||
real(dp), intent(in) :: s(qn)
|
|
||||||
real(dp), intent(out) :: inv_out(:)
|
|
||||||
! real(dp), parameter :: ck = 1.0_dp, dk = 1.0_dp/ck ! scaling for higher order invariants
|
|
||||||
real(dp) inv(24)
|
|
||||||
integer, parameter :: inv_order(12) = & ! the order in which the invariants are selected
|
|
||||||
& [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
|
|
||||||
real(dp) Rch, umb, xR, yR, xAng, yAng
|
|
||||||
! for readability
|
|
||||||
Rch = s(1)
|
|
||||||
xR = s(2)
|
|
||||||
yR = s(3)
|
|
||||||
xAng = s(4)
|
|
||||||
yAng = s(5)
|
|
||||||
umb = s(6)**2
|
|
||||||
! invarianten
|
|
||||||
! a moden
|
|
||||||
inv(1) = Rch
|
|
||||||
inv(2) = umb
|
|
||||||
! invariante e pairs
|
|
||||||
inv(3) = v_e_monom(xR, yR, 1)
|
|
||||||
inv(4) = v_e_monom(xAng, yAng, 1)
|
|
||||||
! third order e pairs
|
|
||||||
inv(5) = v_e_monom(xR, yR, 2)
|
|
||||||
inv(6) = v_e_monom(xAng, yAng, 2)
|
|
||||||
! invariant ee coupling
|
|
||||||
inv(7) = v_ee_monom(xR, yR, xAng, yAng, 1)
|
|
||||||
! mode combinations
|
|
||||||
inv(8) = Rch*umb
|
|
||||||
|
|
||||||
inv(9) = Rch*v_e_monom(xR, yR, 1)
|
|
||||||
inv(10) = umb*v_e_monom(xR, yR, 1)
|
|
||||||
|
|
||||||
inv(11) = Rch*v_e_monom(xAng, yAng, 1)
|
|
||||||
inv(12) = umb*v_e_monom(xAng, yAng, 1)
|
|
||||||
|
|
||||||
! damp coordinates because of second order and higher invariants
|
|
||||||
inv(3) = sign(sqrt(abs(inv(3))), inv(3))
|
|
||||||
inv(4) = sign(sqrt(abs(inv(4))), inv(4))
|
|
||||||
inv(5) = sign((abs(inv(5))**(1./3.)), inv(5))
|
|
||||||
inv(6) = sign((abs(inv(6))**(1./3.)), inv(6))
|
|
||||||
inv(7) = sign((abs(inv(7))**(1./3.)), inv(7))
|
|
||||||
inv(8) = sign(sqrt(abs(inv(8))), inv(8))
|
|
||||||
inv(9) = sign((abs(inv(9))**(1./3.)), inv(9))
|
|
||||||
inv(10) = sign((abs(inv(10))**(1./3.)), inv(10))
|
|
||||||
inv(11) = sign((abs(inv(11))**(1./3.)), inv(11))
|
|
||||||
inv(12) = sign((abs(inv(12))**(1./3.)), inv(12))
|
|
||||||
|
|
||||||
inv_out(:) = inv(inv_order(1:size(inv_out, 1)))
|
|
||||||
|
|
||||||
end subroutine get_invariants
|
|
||||||
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
! % real part of spherical harmonics
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
! Ylm shifted to 0 for theta=0
|
|
||||||
real(dp) function ylm(theta, phi, l, m)
|
|
||||||
implicit none
|
|
||||||
real(dp) theta, phi
|
|
||||||
integer(idp) l, m
|
|
||||||
ylm = plm2(dcos(theta), l, m)*cos(m*phi) - plm2(1.0_dp, l, m)
|
|
||||||
end function ylm
|
|
||||||
!----------------------------------------------------------
|
|
||||||
real(dp) function plm2(x, l, n)
|
|
||||||
implicit none
|
|
||||||
real(dp) x
|
|
||||||
integer(idp) l, m, n
|
|
||||||
|
|
||||||
real(dp) pmm, p_mp1m, pllm
|
|
||||||
integer(idp) ll
|
|
||||||
|
|
||||||
! negative m und bereich von x abfangen
|
|
||||||
if ((l .lt. 0)&
|
|
||||||
&.or. (abs(n) .gt. abs(l))&
|
|
||||||
&.or. (abs(x) .gt. 1.)) then
|
|
||||||
write (6, '(''bad arguments in legendre'')')
|
|
||||||
stop
|
|
||||||
end if
|
|
||||||
|
|
||||||
! fix sign of m to compute the positiv m
|
|
||||||
m = abs(n)
|
|
||||||
|
|
||||||
pmm = (-1)**m*dsqrt(fac(2*m))*1./((2**m)*fac(m))& !compute P(m,m) not P(l,l)
|
|
||||||
&*(dsqrt(1.-x**2))**m
|
|
||||||
|
|
||||||
if (l .eq. m) then
|
|
||||||
plm2 = pmm !P(l,m)=P(m,m)
|
|
||||||
else
|
|
||||||
p_mp1m = x*dsqrt(dble(2*m + 1))*pmm !compute P(m+1,m)
|
|
||||||
if (l .eq. m + 1) then
|
|
||||||
plm2 = p_mp1m !P(l,m)=P(m+1,m)
|
|
||||||
else
|
|
||||||
do ll = m + 2, l
|
|
||||||
pllm = x*(2*l - 1)/dsqrt(dble(l**2 - m**2))*p_mp1m& ! compute P(m+2,m) up to P(l,m) recursively
|
|
||||||
&- dsqrt(dble((l - 1)**2 - m**2))&
|
|
||||||
&/dsqrt(dble(l**2 - m**2))*pmm
|
|
||||||
! schreibe m+2 und m+1 jeweils fuer die naechste iteration
|
|
||||||
pmm = p_mp1m !P(m,m) = P(m+1,m)
|
|
||||||
p_mp1m = pllm !P(m+1,m) = P(m+2,m)
|
|
||||||
end do
|
|
||||||
plm2 = pllm !P(l,m)=P(m+k,m), k element N
|
|
||||||
end if
|
|
||||||
end if
|
|
||||||
|
|
||||||
! sets the phase of -m term right (ignored to gurantee Ylm=(Yl-m)* for JT terms
|
|
||||||
! if(n.lt.0) then
|
|
||||||
! plm2 = (-1)**m * plm2 !* fac(l-m)/fac(l+m)
|
|
||||||
! endif
|
|
||||||
|
|
||||||
end function
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
real(dp) function fac(i)
|
|
||||||
integer(idp) i
|
|
||||||
select case (i)
|
|
||||||
case (0)
|
|
||||||
fac = 1.0_dp
|
|
||||||
case (1)
|
|
||||||
fac = 1.0_dp
|
|
||||||
case (2)
|
|
||||||
fac = 2.0_dp
|
|
||||||
case (3)
|
|
||||||
fac = 6.0_dp
|
|
||||||
case (4)
|
|
||||||
fac = 24.0_dp
|
|
||||||
case (5)
|
|
||||||
fac = 120.0_dp
|
|
||||||
case (6)
|
|
||||||
fac = 720.0_dp
|
|
||||||
case (7)
|
|
||||||
fac = 5040.0_dp
|
|
||||||
case (8)
|
|
||||||
fac = 40320.0_dp
|
|
||||||
case (9)
|
|
||||||
fac = 362880.0_dp
|
|
||||||
case (10)
|
|
||||||
fac = 3628800.0_dp
|
|
||||||
case (11)
|
|
||||||
fac = 39916800.0_dp
|
|
||||||
case (12)
|
|
||||||
fac = 479001600.0_dp
|
|
||||||
case default
|
|
||||||
write (*, *) 'ERROR: no case for given faculty, Max is 12!'
|
|
||||||
stop
|
|
||||||
end select
|
|
||||||
end function fac
|
|
||||||
|
|
||||||
! Does the simplest morse transform possible
|
|
||||||
! one skaling factor + shift
|
|
||||||
function morse_transform(x, p, pst) result(t)
|
|
||||||
real(dp), intent(in) :: x
|
|
||||||
real(dp), intent(in) :: p(11)
|
|
||||||
integer, intent(in) :: pst(2)
|
|
||||||
real(dp) :: t
|
|
||||||
if (pst(2) == 11) then
|
|
||||||
t = 1.0_dp - exp(-abs(p(2))*(x - p(1)))
|
|
||||||
else
|
|
||||||
error stop 'in morse_transform key required or wrong number of parameters'
|
|
||||||
end if
|
|
||||||
end function morse_transform
|
|
||||||
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
! % FUNCTION F(...) ! MAIK DEPRICATING OVER THE TOP MORSE FUNCTION FOR MYSELF
|
|
||||||
! %
|
|
||||||
! % Returns exponent of tunable Morse coordinate
|
|
||||||
! % exponent is polynomial * gaussian (skewed)
|
|
||||||
! % ilabel = 1 or 2 selects the parameters a and sfac to be used
|
|
||||||
! %
|
|
||||||
! % Background: better representation of the prefector in the
|
|
||||||
! % exponend of the morse function.
|
|
||||||
! % Formular: f(r) = lest no3 paper
|
|
||||||
! %
|
|
||||||
! % Variables:
|
|
||||||
! % x: distance of atoms (double)
|
|
||||||
! % p: parameter vector (double[20])
|
|
||||||
! % ii: 1 for CCl and 2 for CCH (int)
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
pure function f(x, p, ii)
|
|
||||||
integer(idp), intent(in) :: ii !1 for CCL and 2 for CCH
|
|
||||||
real(dp), intent(in) :: x !coordinate
|
|
||||||
real(dp), intent(in) :: p(11) !parameter-vector
|
|
||||||
|
|
||||||
integer(idp) i !running index
|
|
||||||
|
|
||||||
real(dp) r !equilibrium distance
|
|
||||||
real(dp) gaus !gaus part of f
|
|
||||||
real(dp) poly !polynom part of f
|
|
||||||
real(dp) skew !tanh part of f
|
|
||||||
|
|
||||||
real(dp) f !prefactor of exponent and returned value
|
|
||||||
|
|
||||||
integer(idp) npoly(2) !order of polynom
|
|
||||||
|
|
||||||
! Maximum polynom order
|
|
||||||
npoly(1) = 5
|
|
||||||
npoly(2) = 5
|
|
||||||
|
|
||||||
! p(1): position of equilibrium
|
|
||||||
! p(2): constant of exponent
|
|
||||||
! p(3): constant for skewing the gaussian
|
|
||||||
! p(4): tuning for skewing the gaussian
|
|
||||||
! p(5): Gaussian exponent
|
|
||||||
! p(6): Shift of Gaussian maximum
|
|
||||||
! p(7)...: polynomial coefficients
|
|
||||||
! p(8+n)...: coefficients of Morse Power series
|
|
||||||
|
|
||||||
! 1-exp{[p(2)+exp{-p(5)[x-p(6)]^2}[Taylor{p(7+n)}(x-p(6))]][x-p(1)]}
|
|
||||||
|
|
||||||
! Tunable Morse function
|
|
||||||
! Power series in Tunable Morse coordinates of order m
|
|
||||||
! exponent is polynomial of order npoly * gaussian + switching function
|
|
||||||
|
|
||||||
! set r r-r_e
|
|
||||||
r = x
|
|
||||||
r = r - p(1)
|
|
||||||
|
|
||||||
! set up skewing function:
|
|
||||||
skew = 0.5_dp*p(3)*(dtanh(dabs(p(4))*(r - p(6))) + 1.0_dp)
|
|
||||||
|
|
||||||
! set up gaussian function:
|
|
||||||
gaus = dexp(-dabs(p(5))*(r - p(6))**2)
|
|
||||||
|
|
||||||
! set up power series:
|
|
||||||
poly = 0.0_dp
|
|
||||||
do i = 0, npoly(ii) - 1
|
|
||||||
poly = poly + p(7 + i)*(r - p(6))**i
|
|
||||||
end do
|
|
||||||
! set up full exponent function:
|
|
||||||
f = dabs(p(2)) + skew + gaus*poly
|
|
||||||
|
|
||||||
end function
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
pure function xproduct(a, b) result(axb)
|
|
||||||
real(dp), intent(in) :: a(3), b(3)
|
|
||||||
real(dp) :: axb(3) !crossproduct a x b
|
|
||||||
axb(1) = a(2)*b(3) - a(3)*b(2)
|
|
||||||
axb(2) = a(3)*b(1) - a(1)*b(3)
|
|
||||||
axb(3) = a(1)*b(2) - a(2)*b(1)
|
|
||||||
end function xproduct
|
|
||||||
|
|
||||||
pure function normalized(v) result(r)
|
|
||||||
real(dp), intent(in) :: v(:)
|
|
||||||
real(dp) :: r(size(v))
|
|
||||||
r = v/norm(v)
|
|
||||||
end function normalized
|
|
||||||
|
|
||||||
pure function norm(v) result(n)
|
|
||||||
real(dp), intent(in) :: v(:)
|
|
||||||
real(dp) n
|
|
||||||
n = dsqrt(sum(v(:)**2))
|
|
||||||
end function norm
|
|
||||||
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
! % FUNCTION Project_Point_Into_Plane(x,n,r0) result(p)
|
|
||||||
! % return the to n orthogonal part of a vector x-r0
|
|
||||||
! % p: projected point in plane
|
|
||||||
! % x: point being projected
|
|
||||||
! % n: normalvector of plane
|
|
||||||
! % r0: Point in plane
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
pure function project_point_into_plane(x, n, r0) result(p)
|
|
||||||
real(dp), intent(in) :: x(:), n(:), r0(:)
|
|
||||||
real(dp) :: p(size(x)), xs(size(x))
|
|
||||||
xs = x - r0
|
|
||||||
p = xs - plane_to_point(x, n, r0)
|
|
||||||
end function project_point_into_plane
|
|
||||||
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
! % Function Plane_To_Point(x,n,r0) result(p)
|
|
||||||
! % p: part of n in x
|
|
||||||
! % x: point being projected
|
|
||||||
! % n: normalvector of plane
|
|
||||||
! % r0: Point in plane
|
|
||||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
||||||
pure function plane_to_point(x, n, r0) result(p)
|
|
||||||
real(dp), intent(in) :: x(:), n(:), r0(:)
|
|
||||||
real(dp) p(size(x)), xs(size(x)), nn(size(n))
|
|
||||||
nn = normalized(n)
|
|
||||||
xs = x - r0
|
|
||||||
p = dot_product(nn, xs)*nn
|
|
||||||
end function plane_to_point
|
|
||||||
|
|
||||||
subroutine check_coordinates(q)
|
|
||||||
! check for faulty kartesain coordinates
|
|
||||||
real(dp), intent(in) :: q(:)
|
|
||||||
integer(idp) :: i
|
|
||||||
if (all(abs(q) <= epsilon(0.0_dp))) then
|
|
||||||
stop 'Error (ctrans): all kartesian coordinates are<=1d-8'
|
|
||||||
end if
|
|
||||||
do i = 1, 9, 3
|
|
||||||
if (all(abs(q(i:i + 2)) <= epsilon(0.0_dp))) then
|
|
||||||
write (*, *) q
|
|
||||||
stop 'Error(ctrans):kartesian coordinates zero for one atom'
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
end subroutine
|
|
||||||
|
|
||||||
pure function rotor_a_to_z(a, z) result(r)
|
|
||||||
real(dp), intent(in) :: a(3), z(3)
|
|
||||||
real(dp) :: r(3, 3)
|
|
||||||
real(dp) :: alpha
|
|
||||||
real(dp) :: s1(3), s(3, 3), rotor(3, 3)
|
|
||||||
s1 = xproduct(normalized(a), normalized(z))
|
|
||||||
alpha = asin(norm(s1))
|
|
||||||
s(:, 1) = normalized(s1)
|
|
||||||
s(:, 2) = normalized(z)
|
|
||||||
s(:, 3) = xproduct(s1, z)
|
|
||||||
rotor = init_rotor(alpha, 0.0_dp, 0.0_dp)
|
|
||||||
r = matmul(s, matmul(rotor, transpose(s)))
|
|
||||||
end function
|
|
||||||
|
|
||||||
! function returning Rz(gamma) * Ry(beta) * Rx(alpha) for basis order xyz
|
|
||||||
pure function init_rotor(alpha, beta, gamma) result(rotor)
|
|
||||||
real(dp), intent(in) :: alpha, beta, gamma
|
|
||||||
real(dp) :: rotor(3, 3)
|
|
||||||
rotor = 0.0_dp
|
|
||||||
rotor(1, 1) = dcos(beta)*dcos(gamma)
|
|
||||||
rotor(1, 2) = dsin(alpha)*dsin(beta)*dcos(gamma)&
|
|
||||||
&- dcos(alpha)*dsin(gamma)
|
|
||||||
rotor(1, 3) = dcos(alpha)*dsin(beta)*dcos(gamma)&
|
|
||||||
&+ dsin(alpha)*dsin(gamma)
|
|
||||||
|
|
||||||
rotor(2, 1) = dcos(beta)*dsin(gamma)
|
|
||||||
rotor(2, 2) = dsin(alpha)*dsin(beta)*dsin(gamma)&
|
|
||||||
&+ dcos(alpha)*dcos(gamma)
|
|
||||||
rotor(2, 3) = dcos(alpha)*dsin(beta)*dsin(gamma)&
|
|
||||||
&- dsin(alpha)*dcos(gamma)
|
|
||||||
|
|
||||||
rotor(3, 1) = -dsin(beta)
|
|
||||||
rotor(3, 2) = dsin(alpha)*dcos(beta)
|
|
||||||
rotor(3, 3) = dcos(alpha)*dcos(beta)
|
|
||||||
end function init_rotor
|
|
||||||
|
|
||||||
pure function create_plane(a, b, c) result(n)
|
|
||||||
real(dp), intent(in) :: a(3), b(3), c(3)
|
|
||||||
real(dp) :: n(3)
|
|
||||||
real(dp) :: axb(3), bxc(3), cxa(3)
|
|
||||||
axb = xproduct(a, b)
|
|
||||||
bxc = xproduct(b, c)
|
|
||||||
cxa = xproduct(c, a)
|
|
||||||
n = normalized(axb + bxc + cxa)
|
|
||||||
end function create_plane
|
|
||||||
|
|
||||||
function symmetrize(q1, q2, q3, sym) result(s)
|
|
||||||
real(dp), intent(in) :: q1, q2, q3
|
|
||||||
character, intent(in) :: sym
|
|
||||||
real(dp) :: s
|
|
||||||
select case (sym)
|
|
||||||
case ('a')
|
|
||||||
s = (q1 + q2 + q3)*sq3
|
|
||||||
case ('x')
|
|
||||||
s = sq6*(2.0_dp*q1 - q2 - q3)
|
|
||||||
case ('y')
|
|
||||||
s = sq2*(q2 - q3)
|
|
||||||
case default
|
|
||||||
write (*, *) 'ERROR: no rule for symmetrize with sym=', sym
|
|
||||||
stop
|
|
||||||
end select
|
|
||||||
end function symmetrize
|
|
||||||
|
|
||||||
subroutine construct_HBend(ex, ey, ph1, ph2, ph3, x_axis, y_axis)
|
|
||||||
real(dp), intent(in) :: ph1(3), ph2(3), ph3(3)
|
|
||||||
real(dp), intent(in) :: x_axis(3), y_axis(3)
|
|
||||||
real(dp), intent(out) :: ex, ey
|
|
||||||
real(dp) :: x1, y1, alpha1
|
|
||||||
real(dp) :: x2, y2, alpha2
|
|
||||||
real(dp) :: x3, y3, alpha3
|
|
||||||
! get x and y components of projected points
|
|
||||||
x1 = dot_product(ph1, x_axis)
|
|
||||||
y1 = dot_product(ph1, y_axis)
|
|
||||||
x2 = dot_product(ph2, x_axis)
|
|
||||||
y2 = dot_product(ph2, y_axis)
|
|
||||||
x3 = dot_product(ph3, x_axis)
|
|
||||||
y3 = dot_product(ph3, y_axis)
|
|
||||||
! -> calculate H deformation angles
|
|
||||||
alpha3 = datan2(y2, x2)
|
|
||||||
alpha2 = -datan2(y3, x3) !-120*ang2rad
|
|
||||||
! write(*,*)' atan2'
|
|
||||||
! write(*,*) 'alpha2:' , alpha2/ang2rad
|
|
||||||
! write(*,*) 'alpha3:' , alpha3/ang2rad
|
|
||||||
if (alpha2 .lt. 0) alpha2 = alpha2 + pi2
|
|
||||||
if (alpha3 .lt. 0) alpha3 = alpha3 + pi2
|
|
||||||
alpha1 = (pi2 - alpha2 - alpha3)
|
|
||||||
! write(*,*)' fixed break line'
|
|
||||||
! write(*,*) 'alpha1:' , alpha1/ang2rad
|
|
||||||
! write(*,*) 'alpha2:' , alpha2/ang2rad
|
|
||||||
! write(*,*) 'alpha3:' , alpha3/ang2rad
|
|
||||||
alpha1 = alpha1 !- 120.0_dp*ang2rad
|
|
||||||
alpha2 = alpha2 !- 120.0_dp*ang2rad
|
|
||||||
alpha3 = alpha3 !- 120.0_dp*ang2rad
|
|
||||||
! write(*,*)' delta alpha'
|
|
||||||
! write(*,*) 'alpha1:' , alpha1/ang2rad
|
|
||||||
! write(*,*) 'alpha2:' , alpha2/ang2rad
|
|
||||||
! write(*,*) 'alpha3:' , alpha3/ang2rad
|
|
||||||
! write(*,*)
|
|
||||||
|
|
||||||
! construct symmetric and antisymmetric H angles
|
|
||||||
ex = symmetrize(alpha1, alpha2, alpha3, 'x')
|
|
||||||
ey = symmetrize(alpha1, alpha2, alpha3, 'y')
|
|
||||||
end subroutine construct_HBend
|
|
||||||
|
|
||||||
pure function construct_umbrella(nh1, nh2, nh3, n)&
|
|
||||||
&result(umb)
|
|
||||||
real(dp), intent(in) :: nh1(3), nh2(3), nh3(3)
|
|
||||||
real(dp), intent(in) :: n(3)
|
|
||||||
real(dp) :: umb
|
|
||||||
real(dp) :: theta(3)
|
|
||||||
! calculate projections for umberella angle
|
|
||||||
theta(1) = dacos(dot_product(n, nh1))
|
|
||||||
theta(2) = dacos(dot_product(n, nh2))
|
|
||||||
theta(3) = dacos(dot_product(n, nh3))
|
|
||||||
! construct umberella angle
|
|
||||||
umb = sum(theta(1:3))/3.0_dp - 90.0_dp*ang2rad
|
|
||||||
end function construct_umbrella
|
|
||||||
|
|
||||||
pure subroutine construct_sphericals&
|
|
||||||
&(theta, phi, cf, xaxis, yaxis, zaxis)
|
|
||||||
real(dp), intent(in) :: cf(3), xaxis(3), yaxis(3), zaxis(3)
|
|
||||||
real(dp), intent(out) :: theta, phi
|
|
||||||
real(dp) :: x, y, z, v(3)
|
|
||||||
v = normalized(cf)
|
|
||||||
x = dot_product(v, normalized(xaxis))
|
|
||||||
y = dot_product(v, normalized(yaxis))
|
|
||||||
z = dot_product(v, normalized(zaxis))
|
|
||||||
theta = dacos(z)
|
|
||||||
phi = -datan2(y, x)
|
|
||||||
end subroutine construct_sphericals
|
|
||||||
|
|
||||||
subroutine int2kart(internal, kart)
|
|
||||||
real(dp), intent(in) :: internal(6)
|
|
||||||
real(dp), intent(out) :: kart(9)
|
|
||||||
real(dp) :: h1x, h1y, h1z
|
|
||||||
real(dp) :: h2x, h2y, h2z
|
|
||||||
real(dp) :: h3x, h3y, h3z
|
|
||||||
real(dp) :: dch0, dch1, dch2, dch3
|
|
||||||
real(dp) :: a1, a2, a3, wci
|
|
||||||
|
|
||||||
kart = 0.0_dp
|
|
||||||
dch1 = dchequi + sq3*internal(1) + 2*sq6*internal(2)
|
|
||||||
dch2 = dchequi + sq3*internal(1) - sq6*internal(2) + sq2*internal(3)
|
|
||||||
dch3 = dchequi + sq3*internal(1) - sq6*internal(2) - sq2*internal(3)
|
|
||||||
a1 = 2*sq6*internal(4)
|
|
||||||
a2 = -sq6*internal(4) + sq2*internal(5)
|
|
||||||
a3 = -sq6*internal(4) - sq2*internal(5)
|
|
||||||
wci = internal(6)
|
|
||||||
|
|
||||||
! Berechnung kartesische Koordinaten
|
|
||||||
! -----------------------
|
|
||||||
h1x = dch1*cos(wci*ang2rad)
|
|
||||||
h1y = 0.0
|
|
||||||
h1z = -dch1*sin(wci*ang2rad)
|
|
||||||
|
|
||||||
h3x = dch2*cos((a2 + 120)*ang2rad)*cos(wci*ang2rad)
|
|
||||||
h3y = dch2*sin((a2 + 120)*ang2rad)*cos(wci*ang2rad)
|
|
||||||
h3z = -dch2*sin(wci*ang2rad)
|
|
||||||
|
|
||||||
h2x = dch3*cos((-a3 - 120)*ang2rad)*cos(wci*ang2rad)
|
|
||||||
h2y = dch3*sin((-a3 - 120)*ang2rad)*cos(wci*ang2rad)
|
|
||||||
h2z = -dch3*sin(wci*ang2rad)
|
|
||||||
|
|
||||||
kart(1) = h1x
|
|
||||||
kart(2) = h1y
|
|
||||||
kart(3) = h1z
|
|
||||||
kart(4) = h2x
|
|
||||||
kart(5) = h2y
|
|
||||||
kart(6) = h2z
|
|
||||||
kart(7) = h3x
|
|
||||||
kart(8) = h3y
|
|
||||||
kart(9) = h3z
|
|
||||||
end subroutine int2kart
|
|
||||||
|
|
||||||
end module ctrans_mod
|
|
||||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,575 +0,0 @@
|
||||||
! Module contains the spherical harmonics up to l=5 m=-l,..,0,..,l listed on https://en.wikipedia.org/wiki/Table_of_spherical_harmonics from 19.07.2022
|
|
||||||
! the functions are implementde by calling switch case function for given m or l value and return the corresdpondig value for given theta and phi
|
|
||||||
! the functions are split for diffrent l values and are named by P_lm.
|
|
||||||
! example for l=1 and m=-1 the realpart of the spherical harmonic for given theta and phi
|
|
||||||
! is returned by calling Re_Y_lm(1,-1,theta,phi) which itself calls the corresponding function P_1m(m,theta) and multilpies it by cos(m*phi) to account for the real part of exp(m*phi*i)
|
|
||||||
! Attention the legendre polynoms are shifted to account for the missing zero order term in spherical harmonic expansions
|
|
||||||
module sphericalharmonics_mod
|
|
||||||
use accuracy_constants, only: dp, idp
|
|
||||||
implicit none
|
|
||||||
real(kind=dp), parameter :: PI = 4.0_dp * atan( 1.0_dp )
|
|
||||||
contains
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Y_lm( l , m , theta , phi ) result( y )
|
|
||||||
integer(kind=idp), intent( in ) :: l , m
|
|
||||||
real(kind=dp), intent( in ) :: theta , phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( l )
|
|
||||||
case (1)
|
|
||||||
y = Y_1m( m , theta , phi )
|
|
||||||
case (2)
|
|
||||||
y = Y_2m( m , theta , phi )
|
|
||||||
case (3)
|
|
||||||
y = Y_3m( m , theta , phi )
|
|
||||||
case (4)
|
|
||||||
y = Y_4m( m , theta , phi )
|
|
||||||
case (5)
|
|
||||||
y = Y_5m( m , theta , phi )
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')&
|
|
||||||
&'order of spherical harmonics not implemented', l
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
end select
|
|
||||||
end function Y_lm
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Re_Y_lm( l , m , theta , phi ) result( y )
|
|
||||||
integer(kind=idp), intent( in ) :: l , m
|
|
||||||
real(kind=dp), intent( in ) :: theta , phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( l )
|
|
||||||
case (1)
|
|
||||||
y = P_1m( m , theta ) * cos(m*phi)
|
|
||||||
case (2)
|
|
||||||
y = P_2m( m , theta ) * cos(m*phi)
|
|
||||||
case (3)
|
|
||||||
y = P_3m( m , theta ) * cos(m*phi)
|
|
||||||
case (4)
|
|
||||||
y = P_4m( m , theta ) * cos(m*phi)
|
|
||||||
case (5)
|
|
||||||
y = P_5m( m , theta ) * cos(m*phi)
|
|
||||||
case (6)
|
|
||||||
y = P_6m( m , theta ) * cos(m*phi)
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')&
|
|
||||||
&'order of spherical harmonics not implemented', l
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
end select
|
|
||||||
end function Re_Y_lm
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Im_Y_lm( l , m , theta , phi ) result( y )
|
|
||||||
integer(kind=idp), intent( in ) :: l , m
|
|
||||||
real(kind=dp), intent( in ) :: theta , phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( l )
|
|
||||||
case (1)
|
|
||||||
y = P_1m( m , theta ) * sin(m*phi)
|
|
||||||
case (2)
|
|
||||||
y = P_2m( m , theta ) * sin(m*phi)
|
|
||||||
case (3)
|
|
||||||
y = P_3m( m , theta ) * sin(m*phi)
|
|
||||||
case (4)
|
|
||||||
y = P_4m( m , theta ) * sin(m*phi)
|
|
||||||
case (5)
|
|
||||||
y = P_5m( m , theta ) * sin(m*phi)
|
|
||||||
case (6)
|
|
||||||
y = P_6m( m , theta ) * sin(m*phi)
|
|
||||||
case default
|
|
||||||
write(errmesg,'(a,i0)')&
|
|
||||||
&'order of spherical harmonics not implemented',l
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
end select
|
|
||||||
end function Im_Y_lm
|
|
||||||
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Y_1m( m , theta , phi ) result( y )
|
|
||||||
integer(kind=idp),intent( in ):: m
|
|
||||||
real(kind=dp),intent( in ):: theta , phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( m )
|
|
||||||
case (-1)
|
|
||||||
y = 0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)*cos(phi)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y = 0.5_dp*sqrt(3.0_dp/PI)*cos(theta)
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y = -0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)*cos(phi)
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(a,i0)') 'in y_1m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
end select
|
|
||||||
end function Y_1m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Y_2m(m,theta,phi) result(y)
|
|
||||||
integer(kind=idp),intent(in):: m
|
|
||||||
real(kind=dp),intent(in):: theta,phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case (m)
|
|
||||||
case (-2)
|
|
||||||
y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*cos(2.0_dp*phi)
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)*cos(theta)*cos(phi)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y=0.25_dp*sqrt(5.0_dp/PI)&
|
|
||||||
&*(3.0_dp*cos(theta)**2-1.0_dp)
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y=-0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)*cos(theta)*cos(phi)
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*cos(2.0_dp*phi)
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in y_2m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function Y_2m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Y_3m(m,theta,phi) result(y)
|
|
||||||
integer(kind=idp), intent(in) :: m
|
|
||||||
real(kind=dp), intent(in) :: theta,phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case (m)
|
|
||||||
case (-3)
|
|
||||||
y=0.125_dp*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*cos(3.0_dp*phi)
|
|
||||||
|
|
||||||
case (-2)
|
|
||||||
y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*cos(theta)*cos(2.0_dp*phi)
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=0.125_dp*sqrt(21.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp)*cos(phi)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y=0.25_dp*sqrt(7.0_dp/PI)&
|
|
||||||
&*(5.0_dp*cos(theta)**3-3.0_dp*cos(theta))
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y=-0.125_dp*sqrt(21.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp)*cos(phi)
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*cos(theta)*cos(2.0_dp*phi)
|
|
||||||
|
|
||||||
case (3)
|
|
||||||
y=-0.125_dp*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*cos(3.0_dp*phi)
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in y_3m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function Y_3m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Y_4m(m,theta,phi) result(y)
|
|
||||||
integer(kind=idp), intent(in) :: m
|
|
||||||
real(kind=dp), intent(in) :: theta,phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case (m)
|
|
||||||
case (-4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4*cos(4.0_dp*phi)
|
|
||||||
|
|
||||||
case (-3)
|
|
||||||
y=0.375_dp*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*cos(theta)*cos(3.0_dp*phi)
|
|
||||||
|
|
||||||
case (-2)
|
|
||||||
y=0.375_dp*sqrt(5.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2&
|
|
||||||
&*(7.0_dp*cos(theta)**2-1)*cos(2.0_dp*phi)
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=0.375_dp*sqrt(5.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(7.0_dp*cos(theta)**3&
|
|
||||||
&-3.0_dp*cos(theta))*cos(phi)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y=(3.0_dp/16.0_dp)/sqrt(PI)&
|
|
||||||
&*(35.0_dp*cos(theta)**4&
|
|
||||||
&-30.0_dp*cos(theta)**2+3.0_dp)
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y=-0.375_dp*sqrt(5.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(7.0_dp*cos(theta)**3&
|
|
||||||
&-3.0_dp*cos(theta))*cos(phi)
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y=0.375_dp*sqrt(5.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*(7.0_dp*cos(theta)**2-1.0_dp)&
|
|
||||||
&*cos(2*phi)
|
|
||||||
|
|
||||||
case (3)
|
|
||||||
y=-0.375_dp*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*cos(theta)*cos(3.0_dp*phi)
|
|
||||||
|
|
||||||
case (4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4*cos(4.0_dp*phi)
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(a,i0)')'in y_4m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function Y_4m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function Y_5m(m,theta,phi) result(y)
|
|
||||||
integer(kind=idp), intent(in) :: m
|
|
||||||
real(kind=dp), intent(in) :: theta,phi
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case (m)
|
|
||||||
case (-5)
|
|
||||||
y=(3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)&
|
|
||||||
&*sin(theta)**5*cos(5*phi)
|
|
||||||
|
|
||||||
case (-4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4*cos(theta)*cos(4*phi)
|
|
||||||
|
|
||||||
case (-3)
|
|
||||||
y=(1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*(9*cos(theta)**2-1.0_dp)*cos(3*phi)
|
|
||||||
|
|
||||||
case (-2)
|
|
||||||
y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*(3*cos(theta)**3-cos(theta))*cos(2*phi)
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=(1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))&
|
|
||||||
&*sin(theta)*(21.0_dp*cos(theta)**4&
|
|
||||||
&-14.0_dp*cos(theta)**2+1.0_dp)*cos(phi)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y=(1.0_dp/16.0_dp)/sqrt(11.0_dp/PI)&
|
|
||||||
&*(63.0_dp*cos(theta)**5-70.0_dp*cos(theta)**3&
|
|
||||||
&+15.0_dp*cos(theta))
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y=(-1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))&
|
|
||||||
&*sin(theta)*(21.0_dp*cos(theta)**4&
|
|
||||||
&-14.0_dp*cos(theta)**2+1.0_dp)*cos(phi)
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*(3*cos(theta)**3-cos(theta))*cos(2*phi)
|
|
||||||
|
|
||||||
case (3)
|
|
||||||
y=(-1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*(9.0_dp*cos(theta)**2-1.0_dp)&
|
|
||||||
&*cos(3.0_dp*phi)
|
|
||||||
|
|
||||||
case (4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4*cos(theta)*cos(4.0_dp*phi)
|
|
||||||
|
|
||||||
case (5)
|
|
||||||
y=(-3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)&
|
|
||||||
&*sin(theta)**5*cos(5.0_dp*phi)
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in y_5m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function Y_5m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function P_1m( m , theta ) result( y )
|
|
||||||
! >Function returns the value of the corresponding normalized Associated legendre polynom for l=1 and given m and theta
|
|
||||||
integer(kind=idp),intent( in ):: m
|
|
||||||
real(kind=dp),intent( in ):: theta
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( m )
|
|
||||||
case (-1)
|
|
||||||
y = 0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y = 0.5_dp*sqrt(3.0_dp/PI)*(cos(theta)-1.0_dp) ! -1 is subtracted to shift so that for theta=0 y=0
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y = -0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in p_1m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function P_1m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function P_2m( m , theta ) result( y )
|
|
||||||
! >Function returns the value of the corresponding normalized Associated legendre polynom for l=2 and given m and theta
|
|
||||||
integer(kind=idp),intent(in):: m
|
|
||||||
real(kind=dp),intent(in):: theta
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( m )
|
|
||||||
case (-2)
|
|
||||||
y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)*cos(theta)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y = (3.0_dp*cos(theta)**2-1.0_dp)
|
|
||||||
y = y - 2.0_dp !2.0 is subtracted to shift so that for theta=0 y=0
|
|
||||||
y = y * 0.25_dp*sqrt(5.0_dp/PI) ! normalize
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y = -0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)*cos(theta)
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y = 0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in p_2m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function P_2m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function P_3m( m , theta ) result( y )
|
|
||||||
! >Function returns the value of the corresponding normalized Associated legendre polynom for l=3 and given m and theta
|
|
||||||
integer(kind=idp), intent(in) :: m
|
|
||||||
real(kind=dp), intent(in) :: theta
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( m )
|
|
||||||
case (-3)
|
|
||||||
y=0.125_dp*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3
|
|
||||||
|
|
||||||
case (-2)
|
|
||||||
y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*cos(theta)
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=0.125_dp*sqrt(21.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(5*cos(theta)**2-1.0_dp)
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y=(5.0_dp*cos(theta)**3-3*cos(theta))
|
|
||||||
y=y-2.0_dp ! 2.0 is subtracted to shift so that for theta=0 y=0
|
|
||||||
y=y*0.25_dp*sqrt(7.0_dp/PI) ! normalize
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y=-0.125_dp*sqrt(21.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp)
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y=0.25*sqrt(105.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*cos(theta)
|
|
||||||
|
|
||||||
case (3)
|
|
||||||
y=-0.125*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in p_3m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function P_3m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function P_4m(m,theta) result(y)
|
|
||||||
! >Function returns the value of the corresponding normalized Associated legendre polynom for l=4 and given m and theta
|
|
||||||
integer(kind=idp), intent(in) :: m
|
|
||||||
real(kind=dp), intent(in) :: theta
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( m )
|
|
||||||
case (-4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4
|
|
||||||
|
|
||||||
case (-3)
|
|
||||||
y=0.375*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*cos(theta)
|
|
||||||
|
|
||||||
case (-2)
|
|
||||||
y=0.375*sqrt(5.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*(7*cos(theta)**2-1)
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=0.375*sqrt(5.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(7*cos(theta)**3-3*cos(theta))
|
|
||||||
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y=(35*cos(theta)**4-30*cos(theta)**2+3)
|
|
||||||
y = y - 8.0_dp !8.0 is subtracted to shift so that for theta=0 y=0
|
|
||||||
y = y * (3.0_dp/16.0_dp)/sqrt(PI)
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y=-0.375*sqrt(5.0_dp/(PI))&
|
|
||||||
&*sin(theta)*(7*cos(theta)**3-3*cos(theta))
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y=0.375*sqrt(5.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*(7*cos(theta)**2-1)
|
|
||||||
|
|
||||||
case (3)
|
|
||||||
y=-0.375*sqrt(35.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*cos(theta)
|
|
||||||
|
|
||||||
case (4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in p_4m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function P_4m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function P_5m(m,theta) result(y)
|
|
||||||
! >Function returns the value of the corresponding normalized Associated legendre polynom for l=5 and given m and theta
|
|
||||||
integer(kind=idp), intent(in) :: m
|
|
||||||
real(kind=dp), intent(in) :: theta
|
|
||||||
real(kind=dp) y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( m )
|
|
||||||
case (-5)
|
|
||||||
y=(3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)&
|
|
||||||
&*sin(theta)**5
|
|
||||||
|
|
||||||
case (-4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4*cos(theta)
|
|
||||||
|
|
||||||
case (-3)
|
|
||||||
y=(1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*(9*cos(theta)**2-1.0_dp)
|
|
||||||
|
|
||||||
case (-2)
|
|
||||||
y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*(3*cos(theta)**3-cos(theta))
|
|
||||||
|
|
||||||
case (-1)
|
|
||||||
y=(1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))&
|
|
||||||
&*sin(theta)*(21*cos(theta)**4-14*cos(theta)**2+1)
|
|
||||||
|
|
||||||
|
|
||||||
case (0)
|
|
||||||
y = (63*cos(theta)**5-70*cos(theta)**3+15*cos(theta))
|
|
||||||
y = y - 8.0_dp !8.0 is subtracted to shift so that for theta=0 y=0
|
|
||||||
y = y * (1.0_dp/16.0_dp)/sqrt(11.0_dp/PI)
|
|
||||||
|
|
||||||
case (1)
|
|
||||||
y=(-1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))&
|
|
||||||
&*sin(theta)*(21*cos(theta)**4-14*cos(theta)**2+1)
|
|
||||||
|
|
||||||
case (2)
|
|
||||||
y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))&
|
|
||||||
&*sin(theta)**2*(3*cos(theta)**3-cos(theta))
|
|
||||||
|
|
||||||
case (3)
|
|
||||||
y=(-1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)&
|
|
||||||
&*sin(theta)**3*(9*cos(theta)**2-1.0_dp)
|
|
||||||
|
|
||||||
case (4)
|
|
||||||
y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)&
|
|
||||||
&*sin(theta)**4*cos(theta)
|
|
||||||
|
|
||||||
case (5)
|
|
||||||
y=(-3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)&
|
|
||||||
&*sin(theta)**5
|
|
||||||
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in p_5m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function P_5m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
function P_6m(m,theta) result(y)
|
|
||||||
! >Function returns the value of the corresponding normalized Associated legendre polynom for l=6 and given m and theta
|
|
||||||
integer(kind=idp), intent(in) :: m
|
|
||||||
real(kind=dp), intent(in) :: theta
|
|
||||||
real(kind=dp):: y
|
|
||||||
character(len=70) :: errmesg
|
|
||||||
select case ( m )
|
|
||||||
case (-6)
|
|
||||||
y = (1.0_dp/64.0_dp)*sqrt(3003.0_dp/PI)&
|
|
||||||
&* sin(theta)**6
|
|
||||||
case (-5)
|
|
||||||
y = (3.0_dp/32.0_dp)*sqrt(1001.0_dp/PI)&
|
|
||||||
&* sin(theta)**5&
|
|
||||||
&* cos(theta)
|
|
||||||
case (-4)
|
|
||||||
y= (3.0_dp/32.0_dp)*sqrt(91.0_dp/(2.0_dp*PI))&
|
|
||||||
&* sin(theta)**4&
|
|
||||||
&* (11*cos(theta)**2 - 1 )
|
|
||||||
case (-3)
|
|
||||||
y= (1.0_dp/32.0_dp)*sqrt(1365.0_dp/PI)&
|
|
||||||
&* sin(theta)**3&
|
|
||||||
&* (11*cos(theta)**3 - 3*cos(theta) )
|
|
||||||
case (-2)
|
|
||||||
y= (1.0_dp/64.0_dp)*sqrt(1365.0_dp/PI)&
|
|
||||||
&* sin(theta)**2&
|
|
||||||
&* (33*cos(theta)**4 - 18*cos(theta)**2 + 1 )
|
|
||||||
case (-1)
|
|
||||||
y= (1.0_dp/16.0_dp)*sqrt(273.0_dp/(2.0_dp*PI))&
|
|
||||||
&* sin(theta)&
|
|
||||||
&* (33*cos(theta)**5 - 30*cos(theta)**3 + 5*cos(theta) )
|
|
||||||
case (0)
|
|
||||||
y = 231*cos(theta)**6 - 315*cos(theta)**4 + 105*cos(theta)**2-5
|
|
||||||
y = y - 16.0_dp !16.0 is subtracted to shift so that for theta=0 y=0
|
|
||||||
y = y * (1.0_dp/32.0_dp)*sqrt(13.0_dp/PI)
|
|
||||||
case (1)
|
|
||||||
y= -(1.0_dp/16.0_dp)*sqrt(273.0_dp/(2.0_dp*PI))&
|
|
||||||
&* sin(theta)&
|
|
||||||
&* (33*cos(theta)**5 - 30*cos(theta)**3 + 5*cos(theta) )
|
|
||||||
case (2)
|
|
||||||
y= (1.0_dp/64.0_dp)*sqrt(1365.0_dp/PI)&
|
|
||||||
&* sin(theta)**2&
|
|
||||||
&* (33*cos(theta)**4 - 18*cos(theta)**2 + 1 )
|
|
||||||
case (3)
|
|
||||||
y= -(1.0_dp/32.0_dp)*sqrt(1365.0_dp/PI)&
|
|
||||||
&* sin(theta)**3&
|
|
||||||
&* (11*cos(theta)**3 - 3*cos(theta) )
|
|
||||||
case (4)
|
|
||||||
y= (3.0_dp/32.0_dp)*sqrt(91.0_dp/(2.0_dp*PI))&
|
|
||||||
&* sin(theta)**4&
|
|
||||||
&* (11*cos(theta)**2 - 1 )
|
|
||||||
case (5)
|
|
||||||
y= -(3.0_dp/32.0_dp)*sqrt(1001.0_dp/PI)&
|
|
||||||
&* sin(theta)**5 * cos(theta)
|
|
||||||
case (6)
|
|
||||||
y = (1.0_dp/64.0_dp)*sqrt(3003.0_dp/PI)&
|
|
||||||
&* sin(theta)**6
|
|
||||||
case default
|
|
||||||
write(errmesg,'(A,i0)')'in p_6m given m not logic, ', m
|
|
||||||
error stop 'error in spherical harmonics' !error stop errmesg
|
|
||||||
|
|
||||||
end select
|
|
||||||
end function P_6m
|
|
||||||
!----------------------------------------------------------------------------------------------------
|
|
||||||
|
|
||||||
end module
|
|
||||||
|
|
@ -66,9 +66,9 @@
|
||||||
|
|
||||||
|
|
||||||
! init eigenvector matrix
|
! init eigenvector matrix
|
||||||
TYPES = int(p(pst(1,32)))
|
TYPES = int(p(pst(1,33)))
|
||||||
|
|
||||||
BLK = int(p(pst(1,32)+1)) ! BLOCK IF TYPE IS 3
|
BLK = int(p(pst(1,33)+1)) ! BLOCK IF TYPE IS 3
|
||||||
u = 0.d0
|
u = 0.d0
|
||||||
vx=0.0d0
|
vx=0.0d0
|
||||||
skip=.false.
|
skip=.false.
|
||||||
|
|
|
||||||
|
|
@ -32,14 +32,14 @@
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
if (pst(2,32) .ne. 2) then
|
if (pst(2,33) .ne. 2) then
|
||||||
|
|
||||||
write(*,*) "Error in Paramater Keys, TYPE_CAL should be 2 parameter", pst(2,32)
|
write(*,*) "Error in Paramater Keys, TYPE_CAL should be 2 parameter", pst(2,33)
|
||||||
stop
|
stop
|
||||||
end if
|
end if
|
||||||
|
|
||||||
TYPES = int(p(pst(1,32)))! TYPE OF THE CALCULATION
|
TYPES = int(p(pst(1,33)))! TYPE OF THE CALCULATION
|
||||||
BLK= int(p(pst(1,32)+1))! BLOCK IF TYPE IS 3
|
BLK= int(p(pst(1,33)+1))! BLOCK IF TYPE IS 3
|
||||||
write(*,*) "TYPE of calculation:",TYPES
|
write(*,*) "TYPE of calculation:",TYPES
|
||||||
|
|
||||||
pt=1
|
pt=1
|
||||||
|
|
@ -109,7 +109,7 @@
|
||||||
call one_dia_upper(mat_z,y(1:ntot,pt))
|
call one_dia_upper(mat_z,y(1:ntot,pt))
|
||||||
else
|
else
|
||||||
write(*,*) "Error in TYPE of calculationss",TYPES
|
write(*,*) "Error in TYPE of calculationss",TYPES
|
||||||
write(*,*) "the value:,", p(pst(1,32))
|
write(*,*) "the value:,", p(pst(1,33))
|
||||||
stop
|
stop
|
||||||
end if
|
end if
|
||||||
pt=pt+1
|
pt=pt+1
|
||||||
|
|
|
||||||
|
|
@ -10,7 +10,7 @@ contains
|
||||||
character(len=1) prefix(4)
|
character(len=1) prefix(4)
|
||||||
parameter (prefix=['N','P','A','S'])
|
parameter (prefix=['N','P','A','S'])
|
||||||
!character (len=20) key(4,25)
|
!character (len=20) key(4,25)
|
||||||
integer,parameter:: np=32
|
integer,parameter:: np=33
|
||||||
character(len=16) parname(np)
|
character(len=16) parname(np)
|
||||||
integer i,j
|
integer i,j
|
||||||
! Defining keys for potential
|
! Defining keys for potential
|
||||||
|
|
@ -63,11 +63,12 @@ contains
|
||||||
parname(27)='LZPE1E2O1'
|
parname(27)='LZPE1E2O1'
|
||||||
parname(28)='LZPE1E2O2'
|
parname(28)='LZPE1E2O2'
|
||||||
parname(29)='LZPA2E1O1'
|
parname(29)='LZPA2E1O1'
|
||||||
parname(30)='LZPA2E2O2'
|
parname(30)='LZPA2E1O2'
|
||||||
parname(31)='LZPA2E2O1'
|
parname(31)='LZPA2E2O1'
|
||||||
|
parname(32)='LZPA2E2O2'
|
||||||
|
|
||||||
|
|
||||||
parname(32)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE
|
parname(33)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE
|
||||||
|
|
||||||
do i=1,np
|
do i=1,np
|
||||||
do j=1,4
|
do j=1,4
|
||||||
|
|
|
||||||
|
|
@ -12,7 +12,7 @@
|
||||||
y=0.0d0
|
y=0.0d0
|
||||||
!y(1)=mx(4,4)+mx(5,5)
|
!y(1)=mx(4,4)+mx(5,5)
|
||||||
|
|
||||||
do i=2,3
|
do i=1,ndiab
|
||||||
y(1)=y(1)+mx(i,i)
|
y(1)=y(1)+mx(i,i)
|
||||||
y(2)=y(2)+my(i,i)
|
y(2)=y(2)+my(i,i)
|
||||||
enddo
|
enddo
|
||||||
|
|
@ -298,8 +298,8 @@
|
||||||
double precision, intent(in) :: p(:)
|
double precision, intent(in) :: p(:)
|
||||||
integer, intent(in) :: id_write
|
integer, intent(in) :: id_write
|
||||||
integer :: type_calc, blk
|
integer :: type_calc, blk
|
||||||
type_calc = int(p(pst(1,32)))
|
type_calc = int(p(pst(1,33)))
|
||||||
blk = int(p(pst(1,32)+1))
|
blk = int(p(pst(1,33)+1))
|
||||||
|
|
||||||
if (type_calc ==1) then
|
if (type_calc ==1) then
|
||||||
write(id_write,*) "Type of calculation: TRACE"
|
write(id_write,*) "Type of calculation: TRACE"
|
||||||
|
|
|
||||||
|
|
@ -64,7 +64,8 @@ module diab_mod
|
||||||
|
|
||||||
id=id+1 ! 4
|
id=id+1 ! 4
|
||||||
! order 2
|
! order 2
|
||||||
e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 5 p
|
|
||||||
|
e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 3 p
|
||||||
+p(pst(1,id)+2)*(xs*xb-ys*yb)
|
+p(pst(1,id)+2)*(xs*xb-ys*yb)
|
||||||
id =id+1 ! 5
|
id =id+1 ! 5
|
||||||
e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
|
e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
|
||||||
|
|
@ -388,13 +389,13 @@ module diab_mod
|
||||||
|
|
||||||
! w and z of E''
|
! w and z of E''
|
||||||
! order 1
|
! order 1
|
||||||
id = id
|
id = id ! 22
|
||||||
e(2,2) = e(2,2) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb
|
e(2,2) = e(2,2) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb
|
||||||
e(3,3) = e(3,3) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
|
e(3,3) = e(3,3) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
|
||||||
e(2,3) = e(2,3) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb
|
e(2,3) = e(2,3) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb
|
||||||
|
|
||||||
! order 2
|
! order 2
|
||||||
id = id +1
|
id = id +1 ! 23
|
||||||
do i =1,3
|
do i =1,3
|
||||||
e(2,2) = e(2,2) + p(pst(1,id)+(i-1))*v2(i+3)
|
e(2,2) = e(2,2) + p(pst(1,id)+(i-1))*v2(i+3)
|
||||||
e(3,3) = e(3,3) - p(pst(1,id)+(i-1))*v2(i+3)
|
e(3,3) = e(3,3) - p(pst(1,id)+(i-1))*v2(i+3)
|
||||||
|
|
@ -404,13 +405,13 @@ module diab_mod
|
||||||
! W and Z of E'
|
! W and Z of E'
|
||||||
! order 1
|
! order 1
|
||||||
|
|
||||||
id = id +1
|
id = id +1 ! 24
|
||||||
e(4,4) = e(4,4) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb
|
e(4,4) = e(4,4) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb
|
||||||
e(5,5) = e(5,5) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
|
e(5,5) = e(5,5) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
|
||||||
e(4,5) = e(4,5) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb
|
e(4,5) = e(4,5) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb
|
||||||
|
|
||||||
! order 2
|
! order 2
|
||||||
id = id +1
|
id = id +1 ! 25
|
||||||
do i =1,3
|
do i =1,3
|
||||||
e(4,4) = e(4,4) + p(pst(1,id)+(i-1))*v2(i+3)
|
e(4,4) = e(4,4) + p(pst(1,id)+(i-1))*v2(i+3)
|
||||||
e(5,5) = e(5,5) - p(pst(1,id)+(i-1))*v2(i+3)
|
e(5,5) = e(5,5) - p(pst(1,id)+(i-1))*v2(i+3)
|
||||||
|
|
@ -420,7 +421,7 @@ module diab_mod
|
||||||
! Pseudo of E' and E''
|
! Pseudo of E' and E''
|
||||||
! it must have odd power of b
|
! it must have odd power of b
|
||||||
|
|
||||||
id = id +1
|
id = id +1 !26
|
||||||
! order 0
|
! order 0
|
||||||
e(2,4) = e(2,4)
|
e(2,4) = e(2,4)
|
||||||
e(3,5) = e(3,5)
|
e(3,5) = e(3,5)
|
||||||
|
|
@ -428,13 +429,13 @@ module diab_mod
|
||||||
e(3,4) = e(3,4) - b*(p(pst(1,id)))
|
e(3,4) = e(3,4) - b*(p(pst(1,id)))
|
||||||
|
|
||||||
! order 1
|
! order 1
|
||||||
id = id +1
|
id = id +1 !27
|
||||||
e(2,4) = e(2,4) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
|
e(2,4) = e(2,4) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
|
||||||
e(3,5) = e(3,5) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
|
e(3,5) = e(3,5) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
|
||||||
e(2,5) = e(2,5) - b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
|
e(2,5) = e(2,5) - b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
|
||||||
e(3,4) = e(3,4) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
|
e(3,4) = e(3,4) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
|
||||||
! order 2
|
! order 2
|
||||||
id = id +1
|
id = id +1 !28
|
||||||
do i=1,3
|
do i=1,3
|
||||||
e(2,4) = e(2,4) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
|
e(2,4) = e(2,4) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
|
||||||
e(3,5) = e(3,5) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
|
e(3,5) = e(3,5) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
|
||||||
|
|
@ -445,11 +446,12 @@ module diab_mod
|
||||||
|
|
||||||
! the coupling between A2' and E''
|
! the coupling between A2' and E''
|
||||||
! order 1
|
! order 1
|
||||||
id = id +1
|
id = id +1 !29
|
||||||
e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)*xb))
|
e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)*xb))
|
||||||
e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)*yb))
|
e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)*yb))
|
||||||
|
|
||||||
id = id +1
|
|
||||||
|
id = id +1 !30
|
||||||
! order 2
|
! order 2
|
||||||
do i=1,3
|
do i=1,3
|
||||||
e(1,2) = e(1,2) + b*(p(pst(1,id)+(i-1)))*v2(i)
|
e(1,2) = e(1,2) + b*(p(pst(1,id)+(i-1)))*v2(i)
|
||||||
|
|
@ -459,11 +461,11 @@ module diab_mod
|
||||||
! the coupling of A2' and E'
|
! the coupling of A2' and E'
|
||||||
|
|
||||||
! order 1
|
! order 1
|
||||||
id = id +1
|
id = id +1 !31
|
||||||
e(1,2) = e(1,2) + (p(pst(1,id))*xs + p(pst(1,id)*xb))
|
e(1,2) = e(1,2) + (p(pst(1,id))*xs + p(pst(1,id)*xb))
|
||||||
e(1,3) = e(1,3) - (p(pst(1,id))*ys + p(pst(1,id)*yb))
|
e(1,3) = e(1,3) - (p(pst(1,id))*ys + p(pst(1,id)*yb))
|
||||||
|
|
||||||
id = id +1
|
id = id +1 ! 32
|
||||||
! order 2
|
! order 2
|
||||||
do i=1,3
|
do i=1,3
|
||||||
e(1,2) = e(1,2) + (p(pst(1,id)+(i-1)))*v2(i)
|
e(1,2) = e(1,2) + (p(pst(1,id)+(i-1)))*v2(i)
|
||||||
|
|
|
||||||
|
|
@ -1,38 +0,0 @@
|
||||||
!*** Relevant parameters for the analytic model
|
|
||||||
!*** offsets:
|
|
||||||
!*** offsets(1): morse equilibrium (N-H, Angström)
|
|
||||||
!*** offsets(2): reference angle (H-N-H)
|
|
||||||
!*** offsets(3): --
|
|
||||||
!*** pat_index: vector giving the position of the
|
|
||||||
!*** various coordinates (see below)
|
|
||||||
!*** ppars: polynomial parameters for tmcs
|
|
||||||
!*** vcfs: coefficients for V expressions.
|
|
||||||
!*** wzcfs: coefficients for W & Z expressions.
|
|
||||||
!*** ifc: inverse factorials.
|
|
||||||
|
|
||||||
integer matdim
|
|
||||||
parameter (matdim=5) ! matrix is (matdim)x(matdim)
|
|
||||||
|
|
||||||
real*8 offsets(2)
|
|
||||||
integer pat_index(maxnin)
|
|
||||||
|
|
||||||
! NH3 params
|
|
||||||
parameter (offsets=[1.0228710942d0,120.d0])
|
|
||||||
|
|
||||||
!##########################################################################
|
|
||||||
! coordinate order; the first #I number of coords are given to the
|
|
||||||
! ANN, where #I is the number of input neurons. The position i in
|
|
||||||
! pat_index corresponds to a coordinate, the value of pat_index(i)
|
|
||||||
! signifies its position.
|
|
||||||
!
|
|
||||||
! The vector is ordered as follows:
|
|
||||||
! a,xs,ys,xb,yb,b,rs**2,rb**2,b**2,
|
|
||||||
! es*eb, es**3, eb**3,es**2*eb, es*eb**2
|
|
||||||
! ri**2 := xi**2+yi**2 = ei**2; ei := (xi,yi), i = s,b
|
|
||||||
!
|
|
||||||
! parts not supposed to be read by ANN are marked by ';' for your
|
|
||||||
! convenience.
|
|
||||||
!##########################################################################
|
|
||||||
! a,rs**2,rb**2,es*eb,es**3,eb**3,es**2*eb,es*eb**2,b**2 #I=9 (6D)
|
|
||||||
parameter (pat_index=[1,2,3,4,5,6,7,8,9,10,11,12,13,14])
|
|
||||||
!**************************************************************************
|
|
||||||
|
|
@ -1,260 +0,0 @@
|
||||||
subroutine cart2int(cart,qint)
|
|
||||||
implicit none
|
|
||||||
! This version merges both coordinate transformation routines into
|
|
||||||
! one. JTmod's sscales(2:3) are ignored.
|
|
||||||
! This is the first version to be compatible with one of my proper 6D fits
|
|
||||||
! Time-stamp: <2024-10-22 13:52:59 dwilliams>
|
|
||||||
|
|
||||||
! Input (cartesian, in Angström)
|
|
||||||
! cart(:,1): N
|
|
||||||
! cart(:,1+i): Hi
|
|
||||||
! Output
|
|
||||||
! qint(i): order defined in JTmod.
|
|
||||||
! Internal Variables
|
|
||||||
! no(1:3): NO distances 1-3
|
|
||||||
! pat_in: temporary coordinates
|
|
||||||
! axis: main axis of NO3
|
|
||||||
include 'nnparams.incl'
|
|
||||||
include 'JTmod.incl'
|
|
||||||
|
|
||||||
|
|
||||||
real*8 cart(3,4),qint(maxnin)
|
|
||||||
|
|
||||||
real*8 no(3), r1, r2, r3
|
|
||||||
real*8 v1(3), v2(3), v3(3)
|
|
||||||
real*8 n1(3), n2(3), n3(3), tr(3)
|
|
||||||
real*8 ortho(3)
|
|
||||||
real*8 pat_in(maxnin)
|
|
||||||
logical ignore_umbrella,dbg_umbrella
|
|
||||||
logical dbg_distances
|
|
||||||
|
|
||||||
!.. Debugging parameters
|
|
||||||
!.. set umbrella to 0
|
|
||||||
parameter (ignore_umbrella=.false.)
|
|
||||||
! parameter (ignore_umbrella=.true.)
|
|
||||||
|
|
||||||
!.. break if umbrella is not 0
|
|
||||||
parameter (dbg_umbrella=.false.)
|
|
||||||
! parameter (dbg_umbrella=.true.)
|
|
||||||
|
|
||||||
!.. break for tiny distances
|
|
||||||
parameter (dbg_distances=.false.)
|
|
||||||
! parameter (dbg_distances=.true.)
|
|
||||||
|
|
||||||
integer k
|
|
||||||
|
|
||||||
!.. get N-O vectors and distances:
|
|
||||||
do k=1,3
|
|
||||||
v1(k)=cart(k,2)-cart(k,1)
|
|
||||||
v2(k)=cart(k,3)-cart(k,1)
|
|
||||||
v3(k)=cart(k,4)-cart(k,1)
|
|
||||||
enddo
|
|
||||||
no(1)=norm(v1,3)
|
|
||||||
no(2)=norm(v2,3)
|
|
||||||
no(3)=norm(v3,3)
|
|
||||||
|
|
||||||
!.. temporarily store displacements
|
|
||||||
do k=1,3
|
|
||||||
pat_in(k)=no(k)-offsets(1)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do k=1,3
|
|
||||||
v1(k)=v1(k)/no(1)
|
|
||||||
v2(k)=v2(k)/no(2)
|
|
||||||
v3(k)=v3(k)/no(3)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!.. compute three normal vectors for the ONO planes:
|
|
||||||
call xprod(n1,v1,v2)
|
|
||||||
call xprod(n2,v2,v3)
|
|
||||||
call xprod(n3,v3,v1)
|
|
||||||
|
|
||||||
do k=1,3
|
|
||||||
tr(k)=(n1(k)+n2(k)+n3(k))/3.d0
|
|
||||||
enddo
|
|
||||||
r1=norm(tr,3)
|
|
||||||
do k=1,3
|
|
||||||
tr(k)=tr(k)/r1
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! rotate trisector
|
|
||||||
call rot_trisec(tr,v1,v2,v3)
|
|
||||||
|
|
||||||
!.. determine trisector angle:
|
|
||||||
if (ignore_umbrella) then
|
|
||||||
pat_in(7)=0.0d0
|
|
||||||
else
|
|
||||||
pat_in(7)=pi/2.0d0 - acos(scalar(v1,tr,3))
|
|
||||||
pat_in(7)=sign(pat_in(7),cart(1,2))
|
|
||||||
endif
|
|
||||||
|
|
||||||
!.. molecule now lies in yz plane, compute projected ONO angles:
|
|
||||||
v1(1)=0.d0
|
|
||||||
v2(1)=0.d0
|
|
||||||
v3(1)=0.d0
|
|
||||||
r1=norm(v1,3)
|
|
||||||
r2=norm(v2,3)
|
|
||||||
r3=norm(v3,3)
|
|
||||||
do k=2,3
|
|
||||||
v1(k)=v1(k)/r1
|
|
||||||
v2(k)=v2(k)/r2
|
|
||||||
v3(k)=v3(k)/r3
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! make orthogonal vector to v3
|
|
||||||
ortho(1)=0.0d0
|
|
||||||
ortho(2)=v3(3)
|
|
||||||
ortho(3)=-v3(2)
|
|
||||||
|
|
||||||
!.. projected ONO angles in radians
|
|
||||||
pat_in(4)=get_ang(v2,v3,ortho)
|
|
||||||
pat_in(5)=get_ang(v1,v3,ortho)
|
|
||||||
|
|
||||||
pat_in(6)=dabs(pat_in(5)-pat_in(4))
|
|
||||||
|
|
||||||
!.. account for rotational order of atoms
|
|
||||||
if (pat_in(4).le.pat_in(5)) then
|
|
||||||
pat_in(5)=2*pi-pat_in(4)-pat_in(6)
|
|
||||||
else
|
|
||||||
pat_in(4)=2*pi-pat_in(5)-pat_in(6)
|
|
||||||
endif
|
|
||||||
|
|
||||||
pat_in(4)=rad2deg*pat_in(4)-offsets(2)
|
|
||||||
pat_in(5)=rad2deg*pat_in(5)-offsets(2)
|
|
||||||
pat_in(6)=rad2deg*pat_in(6)-offsets(2)
|
|
||||||
pat_in(7)=rad2deg*pat_in(7)
|
|
||||||
|
|
||||||
call genANN_ctrans(pat_in)
|
|
||||||
|
|
||||||
qint(:)=pat_in(:)
|
|
||||||
|
|
||||||
contains
|
|
||||||
|
|
||||||
!-------------------------------------------------------------------
|
|
||||||
! compute vector product n1 of vectors v1 x v2
|
|
||||||
subroutine xprod(n1,v1,v2)
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
real*8 n1(3), v1(3), v2(3)
|
|
||||||
|
|
||||||
n1(1) = v1(2)*v2(3) - v1(3)*v2(2)
|
|
||||||
n1(2) = v1(3)*v2(1) - v1(1)*v2(3)
|
|
||||||
n1(3) = v1(1)*v2(2) - v1(2)*v2(1)
|
|
||||||
|
|
||||||
end subroutine
|
|
||||||
|
|
||||||
|
|
||||||
!-------------------------------------------------------------------
|
|
||||||
! compute scalar product of vectors v1 and v2:
|
|
||||||
real*8 function scalar(v1,v2,n)
|
|
||||||
implicit none
|
|
||||||
integer i, n
|
|
||||||
real*8 v1(*), v2(*)
|
|
||||||
|
|
||||||
scalar=0.d0
|
|
||||||
do i=1,n
|
|
||||||
scalar=scalar+v1(i)*v2(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end function
|
|
||||||
|
|
||||||
!-------------------------------------------------------------------
|
|
||||||
! compute norm of vector:
|
|
||||||
real*8 function norm(x,n)
|
|
||||||
implicit none
|
|
||||||
integer i, n
|
|
||||||
real*8 x(*)
|
|
||||||
|
|
||||||
norm=0.d0
|
|
||||||
do i=1,n
|
|
||||||
norm=norm+x(i)**2
|
|
||||||
enddo
|
|
||||||
norm=sqrt(norm)
|
|
||||||
|
|
||||||
end function
|
|
||||||
|
|
||||||
!-------------------------------------------------------------------
|
|
||||||
|
|
||||||
subroutine rot_trisec(tr,v1,v2,v3)
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
real*8 tr(3),v1(3),v2(3),v3(3)
|
|
||||||
|
|
||||||
real*8 vrot(3)
|
|
||||||
real*8 rot_ax(3)
|
|
||||||
real*8 cos_phi,sin_phi
|
|
||||||
|
|
||||||
! evaluate cos(-phi) and sin(-phi), where phi is the angle between
|
|
||||||
! tr and (1,0,0)
|
|
||||||
cos_phi=tr(1)
|
|
||||||
sin_phi=dsqrt(tr(2)**2+tr(3)**2)
|
|
||||||
|
|
||||||
if (sin_phi.lt.1.0d-12) then
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
|
|
||||||
! determine rotational axis
|
|
||||||
rot_ax(1) = 0.0d0
|
|
||||||
rot_ax(2) = tr(3)
|
|
||||||
rot_ax(3) = -tr(2)
|
|
||||||
! normalize
|
|
||||||
rot_ax=rot_ax/sin_phi
|
|
||||||
|
|
||||||
! now the rotation can be done using Rodrigues' rotation formula
|
|
||||||
! v'=v*cos(p) + (k x v)sin(p) + k (k*v) (1-cos(p))
|
|
||||||
! for v=tr k*v vanishes by construction:
|
|
||||||
|
|
||||||
! check that the rotation does what it should
|
|
||||||
call rodrigues(vrot,tr,rot_ax,cos_phi,sin_phi)
|
|
||||||
if (dsqrt(vrot(2)**2+vrot(3)**2).gt.1.0d-12) then
|
|
||||||
write(6,*) "ERROR: BROKEN TRISECTOR"
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
tr=vrot
|
|
||||||
|
|
||||||
call rodrigues(vrot,v1,rot_ax,cos_phi,sin_phi)
|
|
||||||
v1=vrot
|
|
||||||
call rodrigues(vrot,v2,rot_ax,cos_phi,sin_phi)
|
|
||||||
v2=vrot
|
|
||||||
call rodrigues(vrot,v3,rot_ax,cos_phi,sin_phi)
|
|
||||||
v3=vrot
|
|
||||||
|
|
||||||
|
|
||||||
end subroutine
|
|
||||||
|
|
||||||
!-------------------------------------------------------------------
|
|
||||||
|
|
||||||
subroutine rodrigues(vrot,v,axis,cos_phi,sin_phi)
|
|
||||||
implicit none
|
|
||||||
real*8 vrot(3),v(3),axis(3)
|
|
||||||
real*8 cos_phi,sin_phi
|
|
||||||
|
|
||||||
real*8 ortho(3)
|
|
||||||
|
|
||||||
call xprod(ortho,axis,v)
|
|
||||||
vrot = v*cos_phi + ortho*sin_phi
|
|
||||||
> + axis*scalar(axis,v,3)*(1-cos_phi)
|
|
||||||
|
|
||||||
end subroutine
|
|
||||||
|
|
||||||
!-------------------------------------------------------------------
|
|
||||||
|
|
||||||
real*8 function get_ang(v,xaxis,yaxis)
|
|
||||||
implicit none
|
|
||||||
! get normalized [0:2pi) angle from vectors in the yz plane
|
|
||||||
real*8 v(3),xaxis(3),yaxis(3)
|
|
||||||
|
|
||||||
real*8 phi
|
|
||||||
|
|
||||||
real*8 pi
|
|
||||||
parameter (pi=3.141592653589793d0)
|
|
||||||
|
|
||||||
phi=atan2(scalar(yaxis,v,3),scalar(xaxis,v,3))
|
|
||||||
if (phi.lt.0.0d0) then
|
|
||||||
phi=2*pi+phi
|
|
||||||
endif
|
|
||||||
get_ang=phi
|
|
||||||
|
|
||||||
end function
|
|
||||||
|
|
||||||
end subroutine cart2int
|
|
||||||
|
|
@ -1,88 +0,0 @@
|
||||||
!-------------------------------------------------------------------
|
|
||||||
! Time-stamp: "2024-10-09 13:33:50 dwilliams"
|
|
||||||
|
|
||||||
subroutine genANN_ctrans(pat_in)
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
include 'nnparams.incl'
|
|
||||||
include 'JTmod.incl'
|
|
||||||
|
|
||||||
real*8 pat_in(maxnin)
|
|
||||||
|
|
||||||
real*8 raw_in(maxnin),off_in(maxnin),ptrans_in(7)
|
|
||||||
real*8 r0
|
|
||||||
real*8 a,b,xs,ys,xb,yb
|
|
||||||
|
|
||||||
integer k
|
|
||||||
|
|
||||||
off_in(1:7)=pat_in(1:7)
|
|
||||||
r0=offsets(1)
|
|
||||||
|
|
||||||
! transform primitives
|
|
||||||
! recover raw distances from offset coords
|
|
||||||
do k=1,3
|
|
||||||
raw_in(k)=off_in(k)+offsets(1)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do k=1,3
|
|
||||||
ptrans_in(k)=off_in(k)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! rescale ONO angles
|
|
||||||
ptrans_in(4)=deg2rad*off_in(4)
|
|
||||||
ptrans_in(5)=deg2rad*off_in(5)
|
|
||||||
ptrans_in(6)=deg2rad*off_in(6)
|
|
||||||
! rescale umbrella
|
|
||||||
ptrans_in(7)=off_in(7)*deg2rad
|
|
||||||
|
|
||||||
! compute symmetry coordinates
|
|
||||||
|
|
||||||
! A (breathing)
|
|
||||||
a=(ptrans_in(1)+ptrans_in(2)+ptrans_in(3))/dsqrt(3.0d0)
|
|
||||||
! ES
|
|
||||||
call prim2emode(ptrans_in(1:3),xs,ys)
|
|
||||||
! EB
|
|
||||||
call prim2emode(ptrans_in(4:6),xb,yb)
|
|
||||||
! B (umbrella)
|
|
||||||
b=ptrans_in(7)
|
|
||||||
|
|
||||||
! overwrite input with output
|
|
||||||
|
|
||||||
pat_in(pat_index(1))=a ! 1
|
|
||||||
pat_in(pat_index(2))=xs
|
|
||||||
pat_in(pat_index(3))=ys
|
|
||||||
pat_in(pat_index(4))=xb
|
|
||||||
pat_in(pat_index(5))=yb
|
|
||||||
pat_in(pat_index(6))=b
|
|
||||||
! totally symmetric monomials
|
|
||||||
pat_in(pat_index(7))=xs**2 + ys**2 ! 2
|
|
||||||
pat_in(pat_index(8))=xb**2 + yb**2 ! 3
|
|
||||||
pat_in(pat_index(9))=b**2 ! 9
|
|
||||||
pat_in(pat_index(10))=xs*xb+ys*yb ! 4
|
|
||||||
! S^3, B^3
|
|
||||||
pat_in(pat_index(11))=xs*(xs**2-3*ys**2) ! 5
|
|
||||||
pat_in(pat_index(12))=xb*(xb**2-3*yb**2) ! 6
|
|
||||||
! S^2 B, S B^2
|
|
||||||
pat_in(pat_index(13))=xb*(xs**2-ys**2) - 2*yb*xs*ys ! 7
|
|
||||||
pat_in(pat_index(14))=xs*(xb**2-yb**2) - 2*ys*xb*yb ! 8
|
|
||||||
|
|
||||||
do k=11,14
|
|
||||||
pat_in(pat_index(k))=tanh(0.1d0*pat_in(pat_index(k)))*10.0d0
|
|
||||||
enddo
|
|
||||||
|
|
||||||
contains
|
|
||||||
|
|
||||||
subroutine prim2emode(prim,ex,ey)
|
|
||||||
implicit none
|
|
||||||
! Takes a 2D-vector prim and returns the degenerate modes x and y
|
|
||||||
! following our standard conventions.
|
|
||||||
|
|
||||||
real*8 prim(3),ex,ey
|
|
||||||
|
|
||||||
ex=(2.0d0*prim(1)-prim(2)-prim(3))/dsqrt(6.0d0)
|
|
||||||
ey=(prim(2)-prim(3))/dsqrt(2.0d0)
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
end subroutine
|
|
||||||
|
|
@ -1,43 +0,0 @@
|
||||||
!**** Declarations
|
|
||||||
|
|
||||||
real*8 pi
|
|
||||||
real*8 hart2eV, eV2hart
|
|
||||||
real*8 hart2icm, icm2hart
|
|
||||||
real*8 eV2icm, icm2eV
|
|
||||||
real*8 deg2rad, rad2deg
|
|
||||||
integer maxnin,maxnout
|
|
||||||
|
|
||||||
!**********************************************************
|
|
||||||
!**** Parameters
|
|
||||||
!*** maxnin: max. number of neurons in input layer
|
|
||||||
!*** maxnout: max. number of neurons in output layer
|
|
||||||
|
|
||||||
parameter (maxnin=14,maxnout=15)
|
|
||||||
|
|
||||||
!**********************************************************
|
|
||||||
!**** Numerical Parameters
|
|
||||||
!*** infty: largest possible double precision real value.
|
|
||||||
!*** iinfty: largest possible integer value.
|
|
||||||
|
|
||||||
! 3.14159265358979323846264338327950...
|
|
||||||
parameter (pi=3.1415926536D0)
|
|
||||||
|
|
||||||
!**********************************************************
|
|
||||||
!**** Unit Conversion Parameters
|
|
||||||
!*** X2Y: convert from X to Y.
|
|
||||||
!***
|
|
||||||
!*** hart: hartree
|
|
||||||
!*** eV: electron volt
|
|
||||||
!*** icm: inverse centimeters (h*c/cm)
|
|
||||||
!****
|
|
||||||
!*** deg: degree
|
|
||||||
!*** rad: radians
|
|
||||||
|
|
||||||
parameter (hart2icm=219474.69d0)
|
|
||||||
parameter (hart2eV=27.211385d0)
|
|
||||||
parameter (eV2icm=hart2icm/hart2eV)
|
|
||||||
parameter (icm2hart=1.0d0/hart2icm)
|
|
||||||
parameter (eV2hart=1.0d0/hart2eV)
|
|
||||||
parameter (icm2eV=1.0d0/eV2icm)
|
|
||||||
parameter (deg2rad=pi/180.0d0)
|
|
||||||
parameter (rad2deg=1.0d0/deg2rad)
|
|
||||||
|
|
@ -400,10 +400,10 @@
|
||||||
|
|
||||||
call print_plotheader(id_plot,coord,npt,set)
|
call print_plotheader(id_plot,coord,npt,set)
|
||||||
do i=1,npt
|
do i=1,npt
|
||||||
write(id_plot,"(ES16.8,*(3(ES16.8),:))")
|
! write(id_plot,"(ES16.8,*(3(ES16.8),:))")
|
||||||
> x(coord,i), ymod(:,i),y(:,i),(wt(:,i))
|
! > x(coord,i), ymod(:,i),y(:,i),(wt(:,i))
|
||||||
! write(id_plot,"(2ES16.8,*(3(ES16.8),:))")
|
write(id_plot,"(2ES16.8,*(3(ES16.8),:))")
|
||||||
! > x(coord,i), x(coord+1,i),y(:,i)
|
> x(coord,i), x(coord+1,i),y(:,i)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
end subroutine
|
end subroutine
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue