rename the nncoords and invariants then update the makefile
This commit is contained in:
jean paul nshuti
parent
b23f225ec8
commit
55ee0ac472
2
Makefile
2
Makefile
|
|
@ -100,7 +100,7 @@ genann_objects = geNNetic.o iNNterface.o puNNch.o \
|
|||
|
||||
# Objects depending on model
|
||||
#mod_objects = nnmodel.o nncoords.o
|
||||
pln_objects = invariants_no3.o nncoords_no3.o genetic_param_nh3_planar.o model_nh3.o nnadia_nh3.o
|
||||
pln_objects = invariants_nh3.o nncoords_nh3.o genetic_param_nh3_planar.o model_nh3.o nnadia_nh3.o
|
||||
|
||||
|
||||
# Objects belonging to internal library
|
||||
|
|
|
|||
|
|
@ -0,0 +1,82 @@
|
|||
Module invariants_mod
|
||||
implicit none
|
||||
contains
|
||||
!----------------------------------------------------
|
||||
subroutine invariants(a,xs,ys,xb,yb,b,inv)
|
||||
implicit none
|
||||
!include "nnparams.incl"
|
||||
double precision, intent(in) :: a, xs, ys, xb, yb, b
|
||||
double precision, intent(out) :: inv(3)
|
||||
double precision:: invar(24)
|
||||
complex(8) :: q1, q2
|
||||
LOGICAL,PARAMETER:: debg =.TRUE.
|
||||
integer :: i
|
||||
! express the coordinate in complex
|
||||
|
||||
q1 = dcmplx(xs, ys)
|
||||
q2 = dcmplx(xb, yb)
|
||||
|
||||
! compute the invariants
|
||||
invar(24) = a
|
||||
invar(23) =b**2
|
||||
|
||||
! INVARIANTS OF KIND II
|
||||
!------------------------
|
||||
|
||||
invar(1) = dreal( q1 * conjg(q1) ) ! r11
|
||||
invar(2) = dreal( q1 * conjg(q2) ) ! r12
|
||||
invar(3) = dreal( q2 * conjg(q2) ) ! r22
|
||||
invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2
|
||||
|
||||
|
||||
!INVATIANTS OF KIND III
|
||||
!------------------------
|
||||
|
||||
invar(5) = dreal( q1 * q1 * q1 ) ! r111
|
||||
invar(6) = dreal( q1 * q1 * q2 ) ! r112
|
||||
invar(7) = dreal( q1 * q2 * q2 ) ! r122
|
||||
invar(8) = dreal( q2 * q2 * q2 ) ! r222
|
||||
invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2
|
||||
invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2
|
||||
invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2
|
||||
invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222
|
||||
|
||||
! INVARIANTS OF KIND IV
|
||||
!-------------------------
|
||||
|
||||
invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 ))
|
||||
invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 ))
|
||||
invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 ))
|
||||
invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 ))
|
||||
|
||||
! INVARIANTS OF KIND V
|
||||
!----------------------
|
||||
|
||||
invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 ))
|
||||
invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 ))
|
||||
invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 ))
|
||||
invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 ))
|
||||
invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 ))
|
||||
invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 ))
|
||||
|
||||
! the only non zero invariant for bend pure cuts
|
||||
|
||||
inv(1) = invar(1)
|
||||
inv(2) = invar(5)
|
||||
inv(3) = invar(9)
|
||||
|
||||
if (debg) then
|
||||
write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4)
|
||||
write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12)
|
||||
write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16)
|
||||
write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22)
|
||||
|
||||
write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES"
|
||||
write(14,*)"---------------------------------------"
|
||||
write(14,*)"xs =",dreal(q1), "ys=",dimag(q1)
|
||||
endif
|
||||
! modify the invariants to only consider few of them
|
||||
!
|
||||
!invar(13:22)=0.0d0
|
||||
end subroutine invariants
|
||||
end module invariants_mod
|
||||
|
|
@ -0,0 +1,716 @@
|
|||
|
||||
!--------------------------------------------------------------------------------------
|
||||
|
||||
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
! % 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
|
||||
implicit none
|
||||
!include 'only_model.incl'
|
||||
include 'nnparams.incl'
|
||||
! precalculate pi, 2*pi and angle to radian conversion
|
||||
double precision, parameter :: pii = 4.00d0*datan(1.00d0)
|
||||
double precision, parameter :: pi2 = 2.00d0*pii
|
||||
double precision, parameter :: ang2rad = pii/180.00d0
|
||||
! precalculate roots
|
||||
double precision, parameter:: sq2 = 1.00d0/dsqrt(2.00d0)
|
||||
double precision, parameter:: sq3 = 1.00d0/dsqrt(3.00d0)
|
||||
double precision, parameter:: sq6 = 1.00d0/dsqrt(6.00d0)
|
||||
! change distances for equilibrium
|
||||
double precision, parameter :: dchequi = 1.0228710942d0 ! req NH3+
|
||||
|
||||
contains
|
||||
subroutine ctrans(q)
|
||||
!use dim_parameter, only: qn
|
||||
use invariants_mod, only: invariants
|
||||
integer k !running indices
|
||||
double precision, intent(inout) :: q(maxnin) !given coordinates
|
||||
double precision :: s(maxnin) !output coordinates symmetry adapted and scaled
|
||||
double precision :: t(maxnin) !output coordinates symmetry adapted but not scaled
|
||||
! ANN Variables
|
||||
!double precision, optional, intent(out) :: invariants(:)
|
||||
! kartesian coordianates copy from MeF+ so substitute c by n and removed f
|
||||
double precision ch1(3), ch2(3), ch3(3), c_atom(3)
|
||||
double precision nh1(3), nh2(3), nh3(3)
|
||||
double precision zaxis(3), xaxis(3), yaxis(3)
|
||||
double precision ph1(3), ph2(3), ph3(3)
|
||||
! primitive coordinates
|
||||
double precision dch1, dch2, dch3 !nh-distances
|
||||
double precision umb !Umbrella Angle from xy-plane
|
||||
|
||||
! Symmetry coordinates
|
||||
double precision aR !a1-modes H-Dist.,
|
||||
double precision exR, exAng !ex components H-Dist., H-Ang.
|
||||
double precision eyR, eyAng !ey components H-Dist., H-Ang.
|
||||
double precision inv(3)
|
||||
! debugging
|
||||
logical, parameter :: dbg = .false.
|
||||
|
||||
! initialize coordinate vectors
|
||||
s = 0.0d0
|
||||
t = 0.0d0
|
||||
|
||||
! write kartesian coords for readability
|
||||
c_atom(1:3) = 0.0d0 ! N-atom at origin
|
||||
do k = 1, 3
|
||||
ch1(k) = q(k )
|
||||
ch2(k) = q(k + 3)
|
||||
ch3(k) = q(k + 6)
|
||||
end do
|
||||
q=0.d0
|
||||
|
||||
! 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
|
||||
! call invariants and get them
|
||||
! 24 invariants
|
||||
|
||||
call invariants(0.0d0,exR,eyR,exAng,eyAng,umb,inv)
|
||||
q(1:3)=inv(1:3)
|
||||
q(4) = aR
|
||||
q(5) = exR
|
||||
q(6) = eyR
|
||||
q(7) = exAng
|
||||
q(8) = -1.0d0*eyAng
|
||||
q(9) = umb
|
||||
! pairwise distances as second coordinate set
|
||||
!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
|
||||
! RETURN Q AS INTERNAL COORD
|
||||
end subroutine ctrans
|
||||
! subroutine ctrans(q)
|
||||
! use invariants_mod, only: invariants
|
||||
! implicit none
|
||||
! double precision, intent(inout):: q(maxnin)
|
||||
! double precision:: invar(24)
|
||||
! double precision:: a,b,esx,esy,ebx,eby
|
||||
! a=q(1)
|
||||
! esx=q(2)
|
||||
! esy=q(3)
|
||||
! ebx=q(4)
|
||||
! eby=q(5)
|
||||
! b=q(6)
|
||||
! call invariants(a,esx,esy,ebx,eby,b,invar)
|
||||
!
|
||||
! q(1:24)=invar(1:24)
|
||||
! q(25)=esx
|
||||
! q(26)=esy
|
||||
! q(27)=ebx
|
||||
! q(28)=eby
|
||||
! q(29)=b
|
||||
! end subroutine ctrans
|
||||
subroutine pair_distance(q, r)
|
||||
double precision, intent(in) :: q(9)
|
||||
double precision, intent(out) :: r(6)
|
||||
double precision :: atom(3, 4)
|
||||
integer :: n, k, count
|
||||
|
||||
!atom order: H1 H2 H3 N
|
||||
atom(:, 1:3) = reshape(q, [3, 3])
|
||||
atom(:, 4) = (/0.00d0, 0.00d0, 0.00d0/)
|
||||
|
||||
! 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)
|
||||
double precision, intent(in),dimension(3) :: x
|
||||
double precision, intent(in),dimension(11) :: p
|
||||
integer, intent(in),dimension(2) :: pst
|
||||
integer :: k
|
||||
double precision, dimension(3) :: s
|
||||
double precision, 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
|
||||
! double precision, intent(in) :: s(qn)
|
||||
! double precision, intent(out) :: inv_out(:)
|
||||
! ! double precision, parameter :: ck = 1.00d0, dk = 1.00d0/ck ! scaling for higher order invariants
|
||||
! double precision 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]
|
||||
! double precision 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
|
||||
! double precision function ylm(theta, phi, l, m)
|
||||
! implicit none
|
||||
! double precision theta, phi
|
||||
! integer l, m
|
||||
! ylm = plm2(dcos(theta), l, m)*cos(m*phi) - plm2(1.00d0, l, m)
|
||||
! end function ylm
|
||||
!!----------------------------------------------------------
|
||||
! double precision function plm2(x, l, n)
|
||||
! implicit none
|
||||
! double precision x
|
||||
! integer l, m, n
|
||||
!
|
||||
! double precision pmm, p_mp1m, pllm
|
||||
! integer 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*ll - 1)/dsqrt(dble(ll**2 - m**2))*p_mp1m& ! compute P(m+2,m) up to P(l,m) recursively
|
||||
! &- dsqrt(dble((ll - 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
|
||||
!----------------------------------------------------------------------------------------------------
|
||||
double precision function fac(i)
|
||||
integer i
|
||||
select case (i)
|
||||
case (0)
|
||||
fac = 1.00d0
|
||||
case (1)
|
||||
fac = 1.00d0
|
||||
case (2)
|
||||
fac = 2.00d0
|
||||
case (3)
|
||||
fac = 6.00d0
|
||||
case (4)
|
||||
fac = 24.00d0
|
||||
case (5)
|
||||
fac = 120.00d0
|
||||
case (6)
|
||||
fac = 720.00d0
|
||||
case (7)
|
||||
fac = 5040.00d0
|
||||
case (8)
|
||||
fac = 40320.00d0
|
||||
case (9)
|
||||
fac = 362880.00d0
|
||||
case (10)
|
||||
fac = 3628800.00d0
|
||||
case (11)
|
||||
fac = 39916800.00d0
|
||||
case (12)
|
||||
fac = 479001600.00d0
|
||||
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)
|
||||
double precision, intent(in) :: x
|
||||
double precision, intent(in) :: p(11)
|
||||
integer, intent(in) :: pst(2)
|
||||
double precision :: t
|
||||
if (pst(2) == 11) then
|
||||
t = 1.00d0 - 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 , intent(in) :: ii !1 for CCL and 2 for CCH
|
||||
double precision, intent(in) :: x !coordinate
|
||||
double precision, intent(in) :: p(11) !parameter-vector
|
||||
|
||||
integer i !running index
|
||||
|
||||
double precision r !equilibrium distance
|
||||
double precision gaus !gaus part of f
|
||||
double precision poly !polynom part of f
|
||||
double precision skew !tanh part of f
|
||||
|
||||
double precision f !prefactor of exponent and returned value
|
||||
|
||||
integer 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.50d0*p(3)*(dtanh(dabs(p(4))*(r - p(6))) + 1.00d0)
|
||||
|
||||
! set up gaussian function:
|
||||
gaus = dexp(-dabs(p(5))*(r - p(6))**2)
|
||||
|
||||
! set up power series:
|
||||
poly = 0.00d0
|
||||
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)
|
||||
double precision, intent(in) :: a(3), b(3)
|
||||
double precision :: 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)
|
||||
double precision, intent(in) :: v(:)
|
||||
double precision :: r(size(v))
|
||||
r = v/norm(v)
|
||||
end function normalized
|
||||
|
||||
pure function norm(v) result(n)
|
||||
double precision, intent(in) :: v(:)
|
||||
double precision 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)
|
||||
double precision, intent(in) :: x(:), n(:), r0(:)
|
||||
double precision :: 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)
|
||||
double precision, intent(in) :: x(:), n(:), r0(:)
|
||||
double precision 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
|
||||
double precision, intent(in) :: q(:)
|
||||
integer :: i
|
||||
if (all(abs(q) <= epsilon(0.00d0))) 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.00d0))) 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)
|
||||
double precision, intent(in) :: a(3), z(3)
|
||||
double precision :: r(3, 3)
|
||||
double precision :: alpha
|
||||
double precision :: 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.00d0, 0.00d0)
|
||||
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)
|
||||
double precision, intent(in) :: alpha, beta, gamma
|
||||
double precision :: rotor(3, 3)
|
||||
rotor = 0.00d0
|
||||
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)
|
||||
double precision, intent(in) :: a(3), b(3), c(3)
|
||||
double precision :: n(3)
|
||||
double precision :: 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)
|
||||
double precision, intent(in) :: q1, q2, q3
|
||||
character, intent(in) :: sym
|
||||
double precision :: s
|
||||
select case (sym)
|
||||
case ('a')
|
||||
s = (q1 + q2 + q3)*sq3
|
||||
case ('x')
|
||||
s = sq6*(2.00d0*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)
|
||||
double precision, intent(in) :: ph1(3), ph2(3), ph3(3)
|
||||
double precision, intent(in) :: x_axis(3), y_axis(3)
|
||||
double precision, intent(out) :: ex, ey
|
||||
double precision :: x1, y1, alpha1
|
||||
double precision :: x2, y2, alpha2
|
||||
double precision :: 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.00d0*ang2rad
|
||||
alpha2 = alpha2 !- 120.00d0*ang2rad
|
||||
alpha3 = alpha3 !- 120.00d0*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)
|
||||
double precision, intent(in) :: nh1(3), nh2(3), nh3(3)
|
||||
double precision, intent(in) :: n(3)
|
||||
double precision :: umb
|
||||
double precision :: 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.00d0 - 90.00d0*ang2rad
|
||||
end function construct_umbrella
|
||||
|
||||
pure subroutine construct_sphericals&
|
||||
&(theta, phi, cf, xaxis, yaxis, zaxis)
|
||||
double precision, intent(in) :: cf(3), xaxis(3), yaxis(3), zaxis(3)
|
||||
double precision, intent(out) :: theta, phi
|
||||
double precision :: 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)
|
||||
! double precision, intent(in) :: internal(6)
|
||||
! double precision, intent(out) :: kart(9)
|
||||
! double precision :: h1x, h1y, h1z
|
||||
! double precision :: h2x, h2y, h2z
|
||||
! double precision :: h3x, h3y, h3z
|
||||
! double precision :: dch0, dch1, dch2, dch3
|
||||
! double precision :: a1, a2, a3, wci
|
||||
!
|
||||
! kart = 0.00d0
|
||||
! 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
|
||||
!**** Define coordinate transformation applied to the input before fit.
|
||||
!***
|
||||
!***
|
||||
!****Conventions:
|
||||
!***
|
||||
!*** ctrans: subroutine transforming a single point in coordinate space
|
||||
|
||||
subroutine trans_in(pat_in,ntot)
|
||||
use ctrans_mod, only: ctrans
|
||||
implicit none
|
||||
|
||||
include 'nnparams.incl'
|
||||
!include 'nndbg.incl'
|
||||
!include 'nncommon.incl'
|
||||
|
||||
double precision pat_in(maxnin,maxpats)
|
||||
integer ntot
|
||||
|
||||
integer j
|
||||
|
||||
do j=1,ntot
|
||||
call ctrans(pat_in(:,j))
|
||||
! FIRST ELEMENT OF PAT-IN ARE USED BY NEURON NETWORK
|
||||
write(62,'(6f16.8)') pat_in(4:9,j)
|
||||
enddo
|
||||
end
|
||||
Loading…
Reference in New Issue