From 733a5d41720a117e7654d25c16b4a01a1cdf0973 Mon Sep 17 00:00:00 2001 From: jean paul nshuti Date: Wed, 8 Oct 2025 15:50:50 +0200 Subject: [PATCH] first commit of the model of nh3+ --- src/model/.giosaveVT1ebr | 0 src/model/JTmod.incl | 38 + src/model/adia.f90 | 98 + src/model/ctrans.f90 | 367 +++ src/model/data_transform.f90 | 40 + src/model/keys.f90 | 88 + src/model/matrix_form.f90 | 301 +++ src/model/model.f90 | 375 +++ src/model/new_PES/.ctrans.f90.un~ | Bin 0 -> 24607 bytes src/model/new_PES/.xy_stretch_lib.f90.un~ | Bin 0 -> 13220 bytes src/model/new_PES/.xyumb_stretch_lib.f90.swp | Bin 0 -> 12288 bytes src/model/new_PES/.xyumb_stretch_lib.f90.un~ | Bin 0 -> 16728 bytes src/model/new_PES/d3h_umb_stretch_lib.f90 | 82 + src/model/new_PES/keys.f90 | 530 ++++ src/model/new_PES/pes_ctrans.f90 | 671 +++++ src/model/new_PES/pes_model.f90 | 982 ++++++++ src/model/new_PES/plot_surface.f90 | 92 + src/model/new_PES/select_monom_mod.f90 | 2293 ++++++++++++++++++ src/model/new_PES/sphericalharmonics_mod.f90 | 575 +++++ src/model/new_PES/surface_mod.f90 | 71 + src/model/new_PES/xy_stretch_lib.f90 | 82 + src/model/new_PES/xy_stretch_lib.f90~ | 82 + src/model/new_PES/xyumb_stretch_lib.f90 | 83 + src/model/nnparams.incl | 43 + src/model/old_keys.f90 | 104 + src/model/pyramidal_model.f90 | 357 +++ src/model/temp | 110 + src/model/temp.keys | 540 +++++ src/model/weight.f | 50 + src/model/write.f | 757 ++++++ 30 files changed, 8811 insertions(+) create mode 100644 src/model/.giosaveVT1ebr create mode 100644 src/model/JTmod.incl create mode 100644 src/model/adia.f90 create mode 100644 src/model/ctrans.f90 create mode 100644 src/model/data_transform.f90 create mode 100644 src/model/keys.f90 create mode 100644 src/model/matrix_form.f90 create mode 100644 src/model/model.f90 create mode 100644 src/model/new_PES/.ctrans.f90.un~ create mode 100644 src/model/new_PES/.xy_stretch_lib.f90.un~ create mode 100644 src/model/new_PES/.xyumb_stretch_lib.f90.swp create mode 100644 src/model/new_PES/.xyumb_stretch_lib.f90.un~ create mode 100644 src/model/new_PES/d3h_umb_stretch_lib.f90 create mode 100644 src/model/new_PES/keys.f90 create mode 100644 src/model/new_PES/pes_ctrans.f90 create mode 100644 src/model/new_PES/pes_model.f90 create mode 100644 src/model/new_PES/plot_surface.f90 create mode 100644 src/model/new_PES/select_monom_mod.f90 create mode 100644 src/model/new_PES/sphericalharmonics_mod.f90 create mode 100644 src/model/new_PES/surface_mod.f90 create mode 100644 src/model/new_PES/xy_stretch_lib.f90 create mode 100644 src/model/new_PES/xy_stretch_lib.f90~ create mode 100644 src/model/new_PES/xyumb_stretch_lib.f90 create mode 100644 src/model/nnparams.incl create mode 100644 src/model/old_keys.f90 create mode 100644 src/model/pyramidal_model.f90 create mode 100644 src/model/temp create mode 100644 src/model/temp.keys create mode 100644 src/model/weight.f create mode 100644 src/model/write.f diff --git a/src/model/.giosaveVT1ebr b/src/model/.giosaveVT1ebr new file mode 100644 index 0000000..e69de29 diff --git a/src/model/JTmod.incl b/src/model/JTmod.incl new file mode 100644 index 0000000..77e5f69 --- /dev/null +++ b/src/model/JTmod.incl @@ -0,0 +1,38 @@ +!*** 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]) +!************************************************************************** diff --git a/src/model/adia.f90 b/src/model/adia.f90 new file mode 100644 index 0000000..552cc09 --- /dev/null +++ b/src/model/adia.f90 @@ -0,0 +1,98 @@ + module adia_mod + implicit none + contains + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % SUBROUTINE ADIA(N,P,NPAR,ymod,v,u,SKIP) +! % +! % determines the adiabatic energies by diagonalizing diabatic matrix. +! % The Eingenvalues are sorted according to the best fitting ordering +! % of the CI vectors. +! % +! % ATTENTION: The interface has changed. To sort by the ci's, +! % the datavalue of the current points are given +! % +! % input variables: +! % n: number of point (int) +! % p: parameter evector(double[npar]) +! % npar: number of parameters (int) +! % skip: .false. if everything should be done +! % +! % output variables: +! % ymod: firtst nstat energies and than nci*ndiab ci's (double[ntot]) +! % v: eigenvalues (double[ndiab]) +! % u: eigenvectors (double[ndiab,ndiab]) +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + subroutine adia(n,p,npar,ymod,vx,u,skip) + use dim_parameter,only: ndiab,nstat,ntot,nci,pst + use data_module,only: q_m,x1_m,x2_m,y_m + use diabmodel, only:diab + use data_matrix + !use dipole, only: diab + implicit none + + integer i,j !running indices + + integer iref !getting correction or refference + + double precision e(ndiab,ndiab) !full diabatic matrix + double precision mx(ndiab,ndiab) + double precision my(ndiab,ndiab) + double precision vxs,vys,vxb,vyb + + integer n !current point + integer npar !number of parameters + double precision p(npar) !parameters + + double precision u(ndiab,ndiab),ut(ndiab,ndiab) !ci-vectors + double precision ymod(ntot) !fitted data + double precision vx(ndiab),vy(nstat) !eigen values + double precision,allocatable,dimension(:,:):: mat + logical skip,dbg + parameter (dbg=.false.) +! lapack variables + integer,parameter :: lwork = 1000 + double precision work(lwork) + integer info + integer TYPES, BLK ! TYPE OF THE CALCULATION + ! variabke for dgemm + + double precision,dimension(ndiab,ndiab):: ex,ey + double precision:: alpha + integer:: lda,ldb,beta,ldc + double precision,dimension(ndiab,ndiab):: temp1,temp2 + call diab(ex,ey,n,x1_m(:,n),x2_m(:,n),p) +! init eigenvector matrix + + u = 0.d0 + skip=.false. + ymod=0.0d0 + call Full_diab_upper(ex,ey,ymod) + !write(*,'(16f10.4)') ex +! ymod(1)=ex(1,1) +! ymod(2)=ex(1,2) +! ymod(3)=ex(1,3) +! ymod(4)=ex(1,4) +! ymod(5)=ex(2,2) +! ymod(6)=ex(2,3) +! ymod(7)=ex(2,4) +! ymod(8)=ex(3,3) +! ymod(9)=ex(3,4) +! ymod(10)=ex(4,4) + + end subroutine + subroutine matrix_mult(C,A,B,N) + implicit none + integer:: n,i,j,k + double precision,dimension(n,n):: A,B,C + do i = 1, n ! Rows of C + do j = 1, n ! Columns of C + C(i,j) = 0.0 ! Initialize element + do k = 1, n ! Dot product + C(i,j) = C(i,j) + A(i,k) * B(k,j) + end do + end do + end do + end subroutine + + end module adia_mod diff --git a/src/model/ctrans.f90 b/src/model/ctrans.f90 new file mode 100644 index 0000000..1aec059 --- /dev/null +++ b/src/model/ctrans.f90 @@ -0,0 +1,367 @@ +module ctrans_mod + use dim_parameter, only: qn + contains +!! subroutine ctrans + + subroutine ctrans(q,x1,x2) + implicit none + include 'nnparams.incl' + include 'JTmod.incl' + double precision,intent(in):: q(qn) + double precision,intent(out):: x1(qn),x2(qn) + double precision:: cart(3,4),qint(maxnin) + cart(:,1)=0.0d0 + cart(1:3,2:4) = reshape([ q(1:9) ], shape(cart(1:3,2:4))) + call cart2int(cart,qint) + x1(1:qn)=qint(1:qn) + !x1(6)=-x1(6) + x2(1:qn)=0.0d0 !qint(1:qn) + !x1=q + end subroutine ctrans + + + 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 + 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 + + + + + end subroutine + 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 module ctrans_mod + diff --git a/src/model/data_transform.f90 b/src/model/data_transform.f90 new file mode 100644 index 0000000..20043b3 --- /dev/null +++ b/src/model/data_transform.f90 @@ -0,0 +1,40 @@ +! |HYr}z4FSh{q&XU^AY^ye?I=Z{66{nfB*TfYsuk1YklRDRey{B zA7$x(;wj0n;QycZCxKFuKnuTr0l!7U7><9n;rP{;sieScIG$-s5mTie~O&UU9;1%1PV z9*m}wG)vrMjgJON0jQ7rscCFGqUq*IE$aa>Jy?zeAEzj)@)ZXio>Ox2M z6p_3xbmXQHC95B4R6f4y`}Jl9VQ2BcmBOREwte+PL^00OVLD0AtHT-fM^i}q$yxRu z74&l+N4cI*#>&}zjM6^BfRI+iL@4G>golH*M=-A<>D*7BYnBnCxKbTH9HOs*PW7j;% zVn3Xsm6=xK^NL5o&=jKkckXRsdSmWmW)+h>i6^~N&~8=$IS~LDFF>m7i6nAZ9; z!%=@az-$mxdWc_WAOdewSiHE`o91!v;;=W$XdEX+3)XFLvG4lhmP`9%!2UL1)8TlK z_R@*VMp=TgZ)qRKdD4pqz3BjU*tv1qa-H>Y>f)@A|0ie^mN}36X`ErivvXL=bd&+T z3z8VeLDU~NTY-A*?}LcD<#x8V5BuZwehii(cnkdCtZiO~)&X8e@y02DudN0KZym4& zgmHSUy|*n`Krv0O^eP{XCbZEjSO&+%8Jkuvul0+wd=e=iWcf0*ZXK;#_VhBe-8$NC ziFQ?orFV)Riib0Q2M6ksvuQd=kMeXn1Yhmd(&_h339GGiI(s*_J3ANvAo2N44CEwE zGs=YHQB03=bJwb2(OJtuizZN?OiTfkijq*)MCL zmtYqqHcINzhmT8K+>ghTgqq)X1*U~WB2QjGdh)(VAHlq4qsg4S`LWyYyyLzf1QfgV zR!8FrW?nqtpom<2*XCms?Q@hKE2a+C(>OxGKCfBF`UXo2zfFbokJkUpu1h=$n{R!Q5SQ+ciOSPUHy9)Pewz01ij9f$O8YM-~RxRQQ`UUg87Sh})a6 z_6T?HZ8G%TYJo~^p15Ro(833rb&X^OqQfyPJfwLaxRwT_3;>81m`+*!RI-xnB z&%#rHKxy5rY!a4{W=m;4%TuUtP2~$m(}@A=$49W96Z-ehEA5_WYDnxeTG1cDWWgb{ zTo=$^&_aAt<)1cDbA~^t`ASgJDQnz@6P*rgyJ77vaX+EB2uWOjgx`F=p#2a}o=}L}Fz%{aB^CQ>muEs#FW3Vpg4#IR}6qgf&SHMROyx z*m9Vz{19OBc~!YFZp-lDXrzlPKVblAvg5q6rajB^8qYIogxh(IXBRbP+k2I^8S{vw z#D?cwxoPM!NzJNT#3}LSTZ)UBbt{NvBjIk7)1F>a=fA7|?uNhD;x$Uuj0fi5BDv6r zwA|XPf^hL-gbP7fl(qHO%|tPzuNTgAb&E;?JPUyHRZ(}gT5|ytfb@h>cw=r%jl5vyD9@1y#ef(Y2o3$ z*D#+WF~(+AVG=om;F$OUaw46U6X8+XJ4%R-1XHi11e-Pm8u0<@O%IBwWA&Y~9y%WabPb?82e=@> zXMGDh1Q4-=rPrmD&IhI5-bcGE1)C`QoX5XmWxBTdlrP?6@iqXb&}8~c(4d@|NIxZ3j9>cFz`!2q|;hhHK2|U{s|_Y^j7G|f2dd`#IJM0ALh z6Qob>wcWiAGGXo=0Iccg3gP{hYwwHFL(Zr1IH58;!Q}2rYRVjg2NmUyu4xq0vS`Nt zy+#K4a6ONwK~2T-td}SkuHA9%u50hOPSd&ckxPaWSeycpRYVrB zqGuZxUZM$$ZZgCvOW=%4AOgIMaHx0|r=+^S zdZ#T0C5Yb+5WL|(7u@)!!OQ1$;1^e#7S#Zq77{8C+BQpc*#P~Y7wN(Ph6hEq6qR*7 zQAJi(DC>0=UMHkQN#IV}*rlB%0p+O_i#Y#JYGDjx*_VRS@SJNH2u^_;l#19j9JA!= zoT5vR`Au_LAs>8TSY;NL8XS3nN#D(jf&@AQb_1|do)`mm%l4Nhr)l1|@m|`TM3Dg< z8y}Ht(42&o0a(zSfPP(-1oA2S=Mv^7rpTL1gu3k&^*oYz6vS+1mVKNng>}M*6@p5c zgW&LpyPohcl$*^93Phq=$bVg>jY=PMuvQ{l)&_mBu8O!UejsRjo1z7kRj?#&K2YFg z0Z^kI$Vcdh@LnfE`gM4=VVk)9csNGCnSa#u;}P~8X#TJsu~1y+GgpdcLckZZp-5wa zjv09JswO@?$s&n2;-dlyViq^yLF`)xG~-(WNTdjYTvRRl6nHve>P?vRn98RpUQwr3B z>`P2+5S*aWERosU+M=fF=Lv*o=YyQV0A8eK_lksa+%21k?vLRtnKwm_Bu5O8_NeS1 zv<=Ip{9u5C7N_Yk_c4$~*%rSHmpO@Vaiz?JZRv!|@aJMMj0A>BinQgB$sKY?%(c|_ z4!=x!1$5@9!1ysM2t;HZf<+psLdb=~rqKbH4XG<*7ZKl*Rb!* z450rjn++wf(MDS!%N4smVV{JO@7F^Q!N0z1qU~YDCEtUd`DW)@W8Jf4kW#^77Wvs* z(2)%W{+u>hjXtNDhyWcQz~Z!Zao zz=Rxdw}~6|;d2#R3lh9Y3mASj8n(k3IbVA#>eWlQ8d9nBpj**9#oRxFCD`bQEf{XJ z-H(!INjktr*c5BCeww%}PI{*WrVIAVatTF2;2`7?-*s-7j@0Uw{>13We7FTMpjlOm z1*$sUH2;wM3m?)HQWzqBrL4iLtxpu4u0WDTwo|XAUcX+d89@oP>nYkn0?T22PJ<;G zz-rn`OwuwF^+=5=GRlRTqtjSJ-A!1k1jitCTO@+{8s+MwPjG?<8K}g;dEbhDFxZ5} z5x@P`V?mLHg4yWwILc_b5#z(y(W>uQ4LV_+HiPEYwV(lFuwNv5O}=c15I6nf9WWVAIp)73g9OhJ94tvuZca_t`^afS*oyEeMbMIM zk}yz|gDNyi0}E7CKxMh}a-gzutb^*&p*mGr)Yn1v+)o`;mio3ZEP3tdiK;rRWB@GW zDq-85a?8BZGs%|Mfp4q;zCnqZn<^<^-`-hHdXLZ+7naaF+uaq!CA~l2AlKpO?soR( zN4ade>F{B30s5A~Ea@vk$xszThH>rA#slVYOT@!&8vqNE~)zXu9k zcLRq80sRzPqRi#RJWvX!DJuzpYy+g6k_%~I1grQLWZ5MdY7v2-(q$IX4smIAF6w};@)`oojZ@;eDIEYLoNqERqo$@^Ly^~yAKcE zymK3O2fphbeE;r&dwlP%d;H*`yZ_+e(cJ*+!GXK~gGYDYegi~5W<>GT+}0I$kW8j| zR#4z_aEXGLOfw#>UpUndypHFqgvaC=j|LaWq2oXioZbnNOLFEdDxM zgZKzH7%6bJhNoGtqFPmblNHrGGUKG5rWB3=ji*LnsOfh2$@YE;g@EugJ=v^Q`T2#0sE zTX`JiO#8`X&qiHj8OBLo<@mE6t2}AA4}hz_bheMpK|Acmrn!4q_(SLvVtE9c2hitv zdOd#b0iSfXWf!yrPSbq|n}+Ysoq>q8W54dfQ)gl;J|vyxR3A@~{V%7;F?IlHol=P* zY-JV35Or&krr}8Nd&4vTc9nSuq@%Fj;?!x)JSi-geD)JqdzWL0)gwrw6!B(x!XN1K zRAm?m`KHRh;SpOoKu=U#R(OV8?!gr?QR#%G45`5|ndN zN*Xph%fuI<>x4N~6&i;vPb;8+e&x2*&50lDZwd*MO5GJywzLUMCFf#TU}?|e!IYDY zPnp>agqM*)Sb0wqR5c+DipIMl`+Eg&+39dP^p(0j0T2rfl%8Ziuh&Z@v&m9E>pJ!5 zRxqRV(ewpUH|$R3#J5d)`lAwBP7Go3Akbz&LD-Dz7>V&kP)^dH=|{N6b}{`x0TO!VSf9qls&2%8-o* zEq!m493SIO12am1zV4~-ZDbxk1Yslh<(Ak(9*2VE&(nm6eC{c78}~!TKYl_HZ}8jZ z7`p*aX85z)eD8;ypFDXSUkpb1kH^ssH1d;3Ynt!9Cot{idp}Wf{tC(CiXZb<2{3ka zE98Vup$b2@8gh0^b&`wy?5>ucXNgQh7UqjoIZJYsu?wY4O=qOPgXHIPF7B(WX)7M$Vi%(>kPy z$FANPSvfFg8bXQ|T2jcMSS}ytr*2XfBeP402+ju5;bED18plxO&jtd?8|YYCjJdzEubJ+&-izg`1`875cos;-i>Wt3$AfI^H{35Hs zhX0i#SCVvpIVNqbTQbS?md!C9YtytT&zP}OoV;X_iJ^)t$s#$!^VWy37P6%gsC`S# zZm=YJXlJV*%Gv5SdA9l?Z+sx?%vw2ebS4A2{I=)OH82t1i^T;zeLguIo$z)~8TZ6? z6#^%{bMD)j9l(Ve#NKvZ8j+2h>89iGXQdmy4Tor6YmX9|`w|Dl=0F)U6TD4eF+5y{ zr&QiW&bRlS^e=4Z%@0km;OEwY3K_SHXVwwTXDwHkER!C1yx!Dp#0=U3Jo(JGbw1S2 zNXbz=oewka@Xw0!2N?8+ap+E%2aq8#lRIWvfkF}7Ay3y+8My9HPJz$qTv5?g{R1qY z4_9s5UMJr0M%Ho#!&_PBi!(?ix?RcZ&ozd zGulrcY^v%W4OFS8fZ#*>9Gr5*RhZsWTq*XeR-7y(gBHL`t*qfC`*~q^u5hO`JmzCm z13A>^<(W>!WhVxB_74QMw1>>y;iiS`OrJH~w{e+ZT%-dW($!FmaiV4O?2G}sBU!k+ zl97w_W-Asoys%NbB4($OaBA?|x7z*9y=A#K+kL=>A|pMV6U33hE*LCVNwSq&@))K2 zDZ(x|s7slpcP*bv-uMi(@0C^f67;6*%<_IbPx6rr4l}yU;jljp`1GMJ)Hkmr)TD;7 za;^rz@OG!`^)bu$zR&v<>YzVBIs|0Jp1dnsl2g}{%+2PsV)*2 zm&Hz!sz`ka^diq8^uhvbJ{QOx9!}+|`K{d0<(e6HnMw;R%PsSpUfb)AcfM!_rrKfK zZ}u3RN4x@q_G{K$W|*w0(eRAepfMvwz?wtDdK#0YG-R@@pVtB_tQGF(wV)N&_F66` zWqXxCEi__78)O5YKqMedc4);a_nA@8hrXb+XTgP$GtAgCXz>{a8O5bcizRZX%XE&5 z25h;D*LdLFEZ7jL1A{oA{F>YB7$wVmqfa zXo3&qFN!lpDOEJRHDi~!4XrraUALE>tb@WPDd0l}wTi7bXy%a%8^sG$)}sK!@`Ueu zi#+FRZX`Xe??~X=0g$jfiq7SDGig9jwJn`twO4ACt@*n> z3D%tn(|k(35+Kc>-WT{Uu7&BzWiU|x3;Y)t-XmECgXwXB{{ln(+HMFG&P?;bDPqXV zt0b6I&^(FHSWwGt?|;U4EL^;ZytVU`u7O=F1)C3or`B1ZhQ2IPTMD*Ft-C-CyZl*pe=>(tZ=L-$=hfia zoY+eE?i{`wz*kmqBoN-<-oT0QW0N28BB0mf-3XBQF5PIP=2TwAKM}=&Axwy!qQiW;LQ|%ikjKRt9Ib` zN}rej7^=nx-If+ zfdne3B-rIJ&he^&CN7@w9eS{8H+5wi%y-JY8u2cSy{QR9c~O(4aZlmzBNbep(TZHH z+zRR)!^c?kqkrTbE9~8E0ORpvS4Jgv)~|Ogaf~2qL_>V}CL}bpF}t-u#D_K_4hCc? za9+s)>6>fFY5JM5EhVX&vT+2}nMK4u&0D6dmRdc>+sb$%^Cd z_({5}#=!>DBk0RF%=fMwG-*SQ`z{%=(>BKJw2e_aZ8I~p45UMl4naDrK)M9!5~NFz zB@0%(*ao}ve8)}rW+uRzckt}X;+9sN797_5SuX_;A6u)jI~a;GJG8HANXi^Jd+{Rn z-Osw~rJ!JYkPF6IYbBR8fDW2jVGlQMYwXD>b&o*i77%J2gvZ)32EZu0Du>uPB9@qIxYuyst#_JD zk&)!^;C+X<29vn&mAe<*1~JuYxD|I20AjiTr9^uiQMrqV=Dn@u*;OWA65Dgx@mwnU5-x-3h~pO!2u2M84vWQwW1MQm z;;Pl$gwS$L9G$K~%Uy7opIWip^;Uxv5|O*y7~LAQ+||ao)r#dpwwf%)J~dMgg%<-2 VBtD;KHQpdleG%~WKl#Pq{1+WM+ARP8 literal 0 HcmV?d00001 diff --git a/src/model/new_PES/.xy_stretch_lib.f90.un~ b/src/model/new_PES/.xy_stretch_lib.f90.un~ new file mode 100644 index 0000000000000000000000000000000000000000..9327c3e936e39e527e6c027ff3d28e3865cf11ae GIT binary patch literal 13220 zcmeI(J8u&~5C`y!2@oKJcYyFZ0SRX?k0j)QB_bgxQ}L9bh%WI3o#dAso1~?oLwp1Z zx_kgC8X5{jLB}URbflt-1oOY_*y9wpmxS$RrOB<&$=!N>+4bFe{WPw;e^Px=`|0Jr z{`j!|@#33_+4Xn7SKoe`sq}wo4fh38rGanbgZI~-d!F}@V%qopMiiEVQaxYrOK~Nt zF2%KKP^3p;xm;Tao8f#pD)@1=8C9D>t<}u?#iC!QzuiL-Q1noor9i=!!+ZI8ZLt^6 z=DfM@c$P|TW0U;EQ@i=`Un%dVGc9bvl)EViF&=XP?^OyNdpRRcw?ukZj=g!4 z`kjrv{I<(Y2e6`YHk@y5sM{KQVb~{{j>{@F6$x=$CDWO-j=G((7livnv%!0ZBDogU z>t2x=CsiW7myq^s+?q_A^y6f>?%no{tGNOO=YAZ_6i|;*;KoRAYmfw|g?Jy!9y>~B z8>hgXCwPaJTJsH>Y>2B7O+CcTxusTPF;MH9U=CB9!4B=8$b!p4fX^i0 zWjfgi1p{6XfYGb4JXZ;ujd&$kE+qG%tO&J!(N2h}qDY9_s!j2HilX1{iaQle6J`Yr`DQfY}~qVRu9{8B_U~JP5E) zhgoF7Wg*LXIn45091pWQy{iQ0QSkymt-K;3ZmTxlsgybRu1~u+eb{YoIiYz>yaGt8 zu1JX6s!ew?t8Vv727?RCgt)ERY;R|^ zl?#vzEhq;?^TDjNA|Y<8Hrny5v|?$Qk%i<*(RwJWtVoF4s!et*t88b@pn(Ns@-d5f z0YFw*x0>SU1iO%S4>Xz;w!IeAxWe(YI1!B1RV2i1)#iFDt1B<@HLPG95={rOvWkSb zt=d>`W@R-Cw~ZmjVFZYjmlZK`R@db?&SHdf;sR3P$7e4GrxKpnvYon1s0 zToyttrM;z;?WZLf3v&y$oVWGMMx;K^6&OB-^_gT|K^*$ zlc(L>eBmNrT5K^K4>NZ9*RQuO;qq90j=iR{?c69TIr1p&p4iWhG-Gd+hK*Eq6C;D| zq$?ZGCSl%{+GLIWB;9GOJ%9F0V|{gPZS_*{^11UboqKMzaXyJ;BTK_Z7Z-7(Kghf7 zATy~n;dT(K_F`wLd34XN)G6Q;VC*X_9CyGlUI0c*nP64NYQ@|lPZ-E?K z04KqX`x*NMd<#AUSAYZ&coCcc*Y0ENQ}8wj!8&*xJO<{$&3hU94O|D`fveyEya)PV z9kjq7_rM1H46cF$@D6wtJOdc`09w8aGC;a+0CzYAoB~dPJ3|4fBaY2lc`W(VZY@TI zS>8?)>@CxhE3HhB_t7QU3J8lm-*0+ukhAlGl7+{YYVpTc#52Eid>OPZA0v+E54P*4ZgPCTj4hDuw%0YbaHX zhC1>;{Bj6avQxlFsoUj-{>yD~upQ9eeSx z6>kw^Ud+UvtW#e&Y9ovOpxUfAWi_T+*$uO`dAN-@gmBJ7rz_K<6|IwTF@CzzCNc0x zM@qE$8Ll%U44p$|OGAair|hLIscx)!d9$lwvt;#}^uUhYACH62$3sMcECiztRJRwa zP#LZhZRd2JNiM=L$Cw5|n2>AO@#ynJ$Ae`a^{5irs3aIpz+o}w zO%*(f+EhfVy2F-{L&L_$AKJvKtTDQgFI87oTci^Sh!X_ICxrEPHJg#r-U7aiO=i5@ z6LHWLCRKgv%_cZS!Q1zF%jbOn|0sdW^v1chV5gq74NWR^=2dpKpcW=c8YwM|ZADCHH4Z&FBD1YeU>feq7_*itUA`Z-$eQLxG3xQa7f;58A5#`$_FFTwuU1=X`qp(J!cu70C&|5w!c+;L zd&@K!qz4`r*Hjxv5x#67_q_UM-RJcUe3TOfw!>?z)wIJG>qv1zx+GQW6I~+cNuGvs zWxp!^H!}nt)kub%W`Y%yk(W@sSfRZgrK1&v#hwy2VYaN86}E&5$U)D#S$Qy3`Ia|@ zHHJ8?X0ur~S+gz2q9?)~z9&+IGL{(?j+U(3Ii1Je5ff{Z4I3skVnLn%r}KZ5E7jz; z%i&=@*y%M)t;!rLLQY=KNS9el6jh;|xjhM0ri({282-=-W#*Z~{9=J``vTYjp=xdNamdCb2s`fCQ>GDt^j(WBDnCfiicQN)4OD=JC literal 0 HcmV?d00001 diff --git a/src/model/new_PES/.xyumb_stretch_lib.f90.un~ b/src/model/new_PES/.xyumb_stretch_lib.f90.un~ new file mode 100644 index 0000000000000000000000000000000000000000..168ce0965bd90d159e8efef6ec2717b6a10b47a2 GIT binary patch literal 16728 zcmeI(NpI6Y6bJB7%DyjUDVu@PrVt9HtdS5xY&XPyt=z^GONkRD4wNI`f*S&E9Jq2r zs>B6Odp^F|RJzW_he=w8S63%b^*nxY|6%R?$&u+3j|+vuWvWAt<1{_DTC6k%OHReFc|pss z2gNcy@`KO|!lEAxI_0u6Ib3oY^o~7IokQC`swmjnT9}x6zo{F~&Ml0-#IsboIvevB zPobLfoeI5Bi&zFGrrsVg9EHO420BWddK~f49YfU7`WXVaP#0mBA)4_n?A7xBGC)M- zBoX=qIw%s~Q>$0nRnM8cJ=zMJUO0YZwCayVc5v#ZYKwj2GVD-0n--}p*`4Bcm)L<| zx6unKPMROQ=CL=y@(k7a*r$HTf@C4T3*7HC9qcSs+xPnD3zAC34Fk949CrqOcEX{C z+Z<}SH@$(VD}v6jMFE2qpS1vwl@eKyEJS(HPClVLQMHsiI+WWPX1g{Ucr-}$!_ii|dEFiNijBdj zCjx9Ciw~G_&(xmCG9%I?z&9MwasbedIYDB^mr(@Ziu;@rMOhKs5YY!tWI3+jh-RXI zJwP#gqJD+Pi>cZIKt;3J5Z@BoM~;SmWz%nFXb}@wmau3*hSwFbPx}_KAX$iMI>}^v z5}1zCK;(CeE^f72A9q7{tm?%pF!WVtmS7qYFF2R7e*o0V%Mv2n%;G(gm$$RNmehr& znY>5;1*BD1mXv*ZmhR!ax@Pe-i3`l->}Y^iT$T{oW)|+Dytr5xPujwA1&bffYReKL z+sv{(nAi3fw4hwc!Uwa`vV_PsiS|p{fNCT!?d(nq$yF?RD61?>h-@=U_CQ|QS!ZKF zu4cglSz%d1WSd#A!+Bwo-?DIA!(xZAy0V1GHnUv!=XFiYX~DRbg$`n6WeJgOX0h(e z%Q|zV0*XB>atNy`ONeYUOEvB8pBW2LGnSu#xQ+!5U`1sKk!=#{*R;X#-n?BEcXin% zMHsGUal;p;{n)1~#>j$XA=7VZrlq{5j1+{p;ma2A;=;gA>`+oC8yn{G9}h)wJMw@OaC)Vbr*UeSi$Miwt*Rg)z|wwa~bpF`7$Z%+VTFAEs3 zipdfp+swk;l>;-jgxw|asN7P4HxPqTJC&z@$}hPc1D!9BXd3dGGU zZU8GPONeYU3pH(sn?#L=xG>zp;)bu9vV_PsvrN;5xJgqxYXRX_7B+Y>V=-M= zB26M)&d zIV@0s0?`bOv@5dAh%^atMOs(fl7?tIBFwh4@EmH7P&KOvYREDp(j>#x+z`JwI7ro= zPP1=&cj(nSsM-R|GFAM>NarF*f@mSYwH)vm9S!M6!IICv-uDSkvwDXQqylP8Gj*zo zY={@)?d80?>4+t&m#KpH@1~kS?4xQ6fECqdLwuWs`d3p;-bwKvz*kMhE$snf)482r OvLrj|A5JyjJo*7aa$6w) literal 0 HcmV?d00001 diff --git a/src/model/new_PES/d3h_umb_stretch_lib.f90 b/src/model/new_PES/d3h_umb_stretch_lib.f90 new file mode 100644 index 0000000..6d1cf0a --- /dev/null +++ b/src/model/new_PES/d3h_umb_stretch_lib.f90 @@ -0,0 +1,82 @@ +module d3h_umb_stretch_lib + use accuracy_constants, only: dp,idp + implicit none + private + public eval_surface, init_surface,eval_matrix + real(dp), dimension(:), allocatable :: p + contains +subroutine eval_surface(e, w, u, x1) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + use dim_parameter, only: ndiab + implicit none + real(dp), dimension(:, :), intent(out) :: w, u + real(dp), dimension(:), intent(out) :: e + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + real(dp), allocatable, dimension(:, :) :: Mat + + !coordinate transformation if needed + call ctrans(x1, s, t) + + block + ! lapack variables + integer(kind=idp), parameter :: lwork = 1000 + real(kind=dp) work(lwork) + integer(kind=idp) info + !evaluate model + call diab(w, 1, x1, s, t, p, size(p, 1)) + allocate (Mat, source=w) + call dsyev('V', 'U', ndiab, Mat, ndiab, e, work, lwork, info) + u(:, :) = Mat(:, :) + deallocate (Mat) + end block + +end subroutine eval_surface + +subroutine eval_matrix(w,x1) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + implicit none + real(dp), dimension(:, :), intent(out) :: w + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + + !coordinate transformation if needed + call ctrans(x1, s, t) + call diab(w, 1, x1, s, t, p, size(p, 1)) +end subroutine eval_matrix + +subroutine init_surface() + use dim_parameter, only: ndiab, nstat, ntot, nci ,qn + use parameterkeys, only: parameterkey_read + use fileread_mod, only: get_datfile, internalize_datfile + use io_parameters, only: llen + use accuracy_constants, only: dp + implicit none + character(len=llen), allocatable, dimension(:) :: infile + + qn = 9 + ndiab = 4 + nstat = 4 + nci = 4 + ntot = ndiab + nstat + nci + + block + character(len=:),allocatable :: datnam + integer :: linenum + datnam = 'planar+pyramidal.par.save' + ! datnam = 'umbstr.par.save' + call internalize_datfile(datnam, infile, linenum, llen) + end block + + !read parameters from file + block + real(dp), dimension(:), allocatable :: p_spread + integer,dimension(:),allocatable :: p_act + integer :: npar + real(dp), parameter :: facspread = 1.0_dp, gspread = 1.0_dp + call parameterkey_read(infile, size(infile, 1), p, p_act, p_spread, npar, gspread, facspread) + end block +end subroutine init_surface +end module d3h_umb_stretch_lib diff --git a/src/model/new_PES/keys.f90 b/src/model/new_PES/keys.f90 new file mode 100644 index 0000000..b66ee29 --- /dev/null +++ b/src/model/new_PES/keys.f90 @@ -0,0 +1,530 @@ + !Start of Element: EV 3 + !Start of Element: A2V 18 + !Start of Element: A1V 33 + !Start of Element: EWZ 48 + !Start of Element: A1EWZ 60 + !Start of Element: A2EQWZ 72 + !Start of Element: A2A1Q 78 + module keys_mod + implicit none +contains + subroutine init_keys + use io_parameters, only: key + +key(1,1) = 'NEXITEN:' +key(2,1) = 'PEXITEN:' +key(3,1) = 'AEXITEN:' +key(4,1) = 'SEXITEN:' + + +key(1,2) = 'NTMC_CH:' +key(2,2) = 'PTMC_CH:' +key(3,2) = 'ATMC_CH:' +key(4,2) = 'STMC_CH:' + + +key(1,3) = 'NEVA1:' +key(2,3) = 'PEVA1:' +key(3,3) = 'AEVA1:' +key(4,3) = 'SEVA1:' + + +key(1,4) = 'NEVU:' +key(2,4) = 'PEVU:' +key(3,4) = 'AEVU:' +key(4,4) = 'SEVU:' + + +key(1,5) = 'NEVE1:' +key(2,5) = 'PEVE1:' +key(3,5) = 'AEVE1:' +key(4,5) = 'SEVE1:' + + +key(1,6) = 'NEVE2:' +key(2,6) = 'PEVE2:' +key(3,6) = 'AEVE2:' +key(4,6) = 'SEVE2:' + + +key(1,7) = 'NEVA1U:' +key(2,7) = 'PEVA1U:' +key(3,7) = 'AEVA1U:' +key(4,7) = 'SEVA1U:' + + +key(1,8) = 'NEVA1E1:' +key(2,8) = 'PEVA1E1:' +key(3,8) = 'AEVA1E1:' +key(4,8) = 'SEVA1E1:' + + +key(1,9) = 'NEVA1E2:' +key(2,9) = 'PEVA1E2:' +key(3,9) = 'AEVA1E2:' +key(4,9) = 'SEVA1E2:' + + +key(1,10) = 'NEVUE1:' +key(2,10) = 'PEVUE1:' +key(3,10) = 'AEVUE1:' +key(4,10) = 'SEVUE1:' + + +key(1,11) = 'NEVUE2:' +key(2,11) = 'PEVUE2:' +key(3,11) = 'AEVUE2:' +key(4,11) = 'SEVUE2:' + + +key(1,12) = 'NEVE1E2:' +key(2,12) = 'PEVE1E2:' +key(3,12) = 'AEVE1E2:' +key(4,12) = 'SEVE1E2:' + + +key(1,13) = 'NEVA1UE1:' +key(2,13) = 'PEVA1UE1:' +key(3,13) = 'AEVA1UE1:' +key(4,13) = 'SEVA1UE1:' + + +key(1,14) = 'NEVA1UE2:' +key(2,14) = 'PEVA1UE2:' +key(3,14) = 'AEVA1UE2:' +key(4,14) = 'SEVA1UE2:' + + +key(1,15) = 'NEVA1E1E2:' +key(2,15) = 'PEVA1E1E2:' +key(3,15) = 'AEVA1E1E2:' +key(4,15) = 'SEVA1E1E2:' + + +key(1,16) = 'NEVUE1E2:' +key(2,16) = 'PEVUE1E2:' +key(3,16) = 'AEVUE1E2:' +key(4,16) = 'SEVUE1E2:' + + +key(1,17) = 'NEVA1UE1E2:' +key(2,17) = 'PEVA1UE1E2:' +key(3,17) = 'AEVA1UE1E2:' +key(4,17) = 'SEVA1UE1E2:' + + +key(1,18) = 'NA2VA1:' +key(2,18) = 'PA2VA1:' +key(3,18) = 'AA2VA1:' +key(4,18) = 'SA2VA1:' + + +key(1,19) = 'NA2VU:' +key(2,19) = 'PA2VU:' +key(3,19) = 'AA2VU:' +key(4,19) = 'SA2VU:' + + +key(1,20) = 'NA2VE1:' +key(2,20) = 'PA2VE1:' +key(3,20) = 'AA2VE1:' +key(4,20) = 'SA2VE1:' + + +key(1,21) = 'NA2VE2:' +key(2,21) = 'PA2VE2:' +key(3,21) = 'AA2VE2:' +key(4,21) = 'SA2VE2:' + + +key(1,22) = 'NA2VA1U:' +key(2,22) = 'PA2VA1U:' +key(3,22) = 'AA2VA1U:' +key(4,22) = 'SA2VA1U:' + + +key(1,23) = 'NA2VA1E1:' +key(2,23) = 'PA2VA1E1:' +key(3,23) = 'AA2VA1E1:' +key(4,23) = 'SA2VA1E1:' + + +key(1,24) = 'NA2VA1E2:' +key(2,24) = 'PA2VA1E2:' +key(3,24) = 'AA2VA1E2:' +key(4,24) = 'SA2VA1E2:' + + +key(1,25) = 'NA2VUE1:' +key(2,25) = 'PA2VUE1:' +key(3,25) = 'AA2VUE1:' +key(4,25) = 'SA2VUE1:' + + +key(1,26) = 'NA2VUE2:' +key(2,26) = 'PA2VUE2:' +key(3,26) = 'AA2VUE2:' +key(4,26) = 'SA2VUE2:' + + +key(1,27) = 'NA2VE1E2:' +key(2,27) = 'PA2VE1E2:' +key(3,27) = 'AA2VE1E2:' +key(4,27) = 'SA2VE1E2:' + + +key(1,28) = 'NA2VA1UE1:' +key(2,28) = 'PA2VA1UE1:' +key(3,28) = 'AA2VA1UE1:' +key(4,28) = 'SA2VA1UE1:' + + +key(1,29) = 'NA2VA1UE2:' +key(2,29) = 'PA2VA1UE2:' +key(3,29) = 'AA2VA1UE2:' +key(4,29) = 'SA2VA1UE2:' + + +key(1,30) = 'NA2VA1E1E2:' +key(2,30) = 'PA2VA1E1E2:' +key(3,30) = 'AA2VA1E1E2:' +key(4,30) = 'SA2VA1E1E2:' + + +key(1,31) = 'NA2VUE1E2:' +key(2,31) = 'PA2VUE1E2:' +key(3,31) = 'AA2VUE1E2:' +key(4,31) = 'SA2VUE1E2:' + + +key(1,32) = 'NA2VA1UE1E2:' +key(2,32) = 'PA2VA1UE1E2:' +key(3,32) = 'AA2VA1UE1E2:' +key(4,32) = 'SA2VA1UE1E2:' + + +key(1,33) = 'NA1VA1:' +key(2,33) = 'PA1VA1:' +key(3,33) = 'AA1VA1:' +key(4,33) = 'SA1VA1:' + + +key(1,34) = 'NA1VU:' +key(2,34) = 'PA1VU:' +key(3,34) = 'AA1VU:' +key(4,34) = 'SA1VU:' + + +key(1,35) = 'NA1VE1:' +key(2,35) = 'PA1VE1:' +key(3,35) = 'AA1VE1:' +key(4,35) = 'SA1VE1:' + + +key(1,36) = 'NA1VE2:' +key(2,36) = 'PA1VE2:' +key(3,36) = 'AA1VE2:' +key(4,36) = 'SA1VE2:' + + +key(1,37) = 'NA1VA1U:' +key(2,37) = 'PA1VA1U:' +key(3,37) = 'AA1VA1U:' +key(4,37) = 'SA1VA1U:' + + +key(1,38) = 'NA1VA1E1:' +key(2,38) = 'PA1VA1E1:' +key(3,38) = 'AA1VA1E1:' +key(4,38) = 'SA1VA1E1:' + + +key(1,39) = 'NA1VA1E2:' +key(2,39) = 'PA1VA1E2:' +key(3,39) = 'AA1VA1E2:' +key(4,39) = 'SA1VA1E2:' + + +key(1,40) = 'NA1VUE1:' +key(2,40) = 'PA1VUE1:' +key(3,40) = 'AA1VUE1:' +key(4,40) = 'SA1VUE1:' + + +key(1,41) = 'NA1VUE2:' +key(2,41) = 'PA1VUE2:' +key(3,41) = 'AA1VUE2:' +key(4,41) = 'SA1VUE2:' + + +key(1,42) = 'NA1VE1E2:' +key(2,42) = 'PA1VE1E2:' +key(3,42) = 'AA1VE1E2:' +key(4,42) = 'SA1VE1E2:' + + +key(1,43) = 'NA1VA1UE1:' +key(2,43) = 'PA1VA1UE1:' +key(3,43) = 'AA1VA1UE1:' +key(4,43) = 'SA1VA1UE1:' + + +key(1,44) = 'NA1VA1UE2:' +key(2,44) = 'PA1VA1UE2:' +key(3,44) = 'AA1VA1UE2:' +key(4,44) = 'SA1VA1UE2:' + + +key(1,45) = 'NA1VA1E1E2:' +key(2,45) = 'PA1VA1E1E2:' +key(3,45) = 'AA1VA1E1E2:' +key(4,45) = 'SA1VA1E1E2:' + + +key(1,46) = 'NA1VUE1E2:' +key(2,46) = 'PA1VUE1E2:' +key(3,46) = 'AA1VUE1E2:' +key(4,46) = 'SA1VUE1E2:' + + +key(1,47) = 'NA1VA1UE1E2:' +key(2,47) = 'PA1VA1UE1E2:' +key(3,47) = 'AA1VA1UE1E2:' +key(4,47) = 'SA1VA1UE1E2:' + + +key(1,48) = 'NEWZE1:' +key(2,48) = 'PEWZE1:' +key(3,48) = 'AEWZE1:' +key(4,48) = 'SEWZE1:' + + +key(1,49) = 'NEWZE2:' +key(2,49) = 'PEWZE2:' +key(3,49) = 'AEWZE2:' +key(4,49) = 'SEWZE2:' + + +key(1,50) = 'NEWZE1A1:' +key(2,50) = 'PEWZE1A1:' +key(3,50) = 'AEWZE1A1:' +key(4,50) = 'SEWZE1A1:' + + +key(1,51) = 'NEWZE2A1:' +key(2,51) = 'PEWZE2A1:' +key(3,51) = 'AEWZE2A1:' +key(4,51) = 'SEWZE2A1:' + + +key(1,52) = 'NEWZE1U:' +key(2,52) = 'PEWZE1U:' +key(3,52) = 'AEWZE1U:' +key(4,52) = 'SEWZE1U:' + + +key(1,53) = 'NEWZE2U:' +key(2,53) = 'PEWZE2U:' +key(3,53) = 'AEWZE2U:' +key(4,53) = 'SEWZE2U:' + + +key(1,54) = 'NEWZE1A1U:' +key(2,54) = 'PEWZE1A1U:' +key(3,54) = 'AEWZE1A1U:' +key(4,54) = 'SEWZE1A1U:' + + +key(1,55) = 'NEWZE2A1U:' +key(2,55) = 'PEWZE2A1U:' +key(3,55) = 'AEWZE2A1U:' +key(4,55) = 'SEWZE2A1U:' + + +key(1,56) = 'NEWZE1E2:' +key(2,56) = 'PEWZE1E2:' +key(3,56) = 'AEWZE1E2:' +key(4,56) = 'SEWZE1E2:' + + +key(1,57) = 'NEWZE1E2A1:' +key(2,57) = 'PEWZE1E2A1:' +key(3,57) = 'AEWZE1E2A1:' +key(4,57) = 'SEWZE1E2A1:' + + +key(1,58) = 'NEWZE1E2U:' +key(2,58) = 'PEWZE1E2U:' +key(3,58) = 'AEWZE1E2U:' +key(4,58) = 'SEWZE1E2U:' + + +key(1,59) = 'NEWZE1E2A1U:' +key(2,59) = 'PEWZE1E2A1U:' +key(3,59) = 'AEWZE1E2A1U:' +key(4,59) = 'SEWZE1E2A1U:' + + +key(1,60) = 'NA1EWZE1:' +key(2,60) = 'PA1EWZE1:' +key(3,60) = 'AA1EWZE1:' +key(4,60) = 'SA1EWZE1:' + + +key(1,61) = 'NA1EWZE2:' +key(2,61) = 'PA1EWZE2:' +key(3,61) = 'AA1EWZE2:' +key(4,61) = 'SA1EWZE2:' + + +key(1,62) = 'NA1EWZE1A1:' +key(2,62) = 'PA1EWZE1A1:' +key(3,62) = 'AA1EWZE1A1:' +key(4,62) = 'SA1EWZE1A1:' + + +key(1,63) = 'NA1EWZE2A1:' +key(2,63) = 'PA1EWZE2A1:' +key(3,63) = 'AA1EWZE2A1:' +key(4,63) = 'SA1EWZE2A1:' + + +key(1,64) = 'NA1EWZE1U:' +key(2,64) = 'PA1EWZE1U:' +key(3,64) = 'AA1EWZE1U:' +key(4,64) = 'SA1EWZE1U:' + + +key(1,65) = 'NA1EWZE2U:' +key(2,65) = 'PA1EWZE2U:' +key(3,65) = 'AA1EWZE2U:' +key(4,65) = 'SA1EWZE2U:' + + +key(1,66) = 'NA1EWZE1A1U:' +key(2,66) = 'PA1EWZE1A1U:' +key(3,66) = 'AA1EWZE1A1U:' +key(4,66) = 'SA1EWZE1A1U:' + + +key(1,67) = 'NA1EWZE2A1U:' +key(2,67) = 'PA1EWZE2A1U:' +key(3,67) = 'AA1EWZE2A1U:' +key(4,67) = 'SA1EWZE2A1U:' + + +key(1,68) = 'NA1EWZE1E2:' +key(2,68) = 'PA1EWZE1E2:' +key(3,68) = 'AA1EWZE1E2:' +key(4,68) = 'SA1EWZE1E2:' + + +key(1,69) = 'NA1EWZE1E2A1:' +key(2,69) = 'PA1EWZE1E2A1:' +key(3,69) = 'AA1EWZE1E2A1:' +key(4,69) = 'SA1EWZE1E2A1:' + + +key(1,70) = 'NA1EWZE1E2U:' +key(2,70) = 'PA1EWZE1E2U:' +key(3,70) = 'AA1EWZE1E2U:' +key(4,70) = 'SA1EWZE1E2U:' + + +key(1,71) = 'NA1EWZE1E2A1U:' +key(2,71) = 'PA1EWZE1E2A1U:' +key(3,71) = 'AA1EWZE1E2A1U:' +key(4,71) = 'SA1EWZE1E2A1U:' + + +key(1,72) = 'NA2EQWZE1U:' +key(2,72) = 'PA2EQWZE1U:' +key(3,72) = 'AA2EQWZE1U:' +key(4,72) = 'SA2EQWZE1U:' + + +key(1,73) = 'NA2EQWZE2U:' +key(2,73) = 'PA2EQWZE2U:' +key(3,73) = 'AA2EQWZE2U:' +key(4,73) = 'SA2EQWZE2U:' + + +key(1,74) = 'NA2EQWZE1UA1:' +key(2,74) = 'PA2EQWZE1UA1:' +key(3,74) = 'AA2EQWZE1UA1:' +key(4,74) = 'SA2EQWZE1UA1:' + + +key(1,75) = 'NA2EQWZE2UA1:' +key(2,75) = 'PA2EQWZE2UA1:' +key(3,75) = 'AA2EQWZE2UA1:' +key(4,75) = 'SA2EQWZE2UA1:' + + +key(1,76) = 'NA2EQWZE1E2U:' +key(2,76) = 'PA2EQWZE1E2U:' +key(3,76) = 'AA2EQWZE1E2U:' +key(4,76) = 'SA2EQWZE1E2U:' + + +key(1,77) = 'NA2EQWZE1E2UA1:' +key(2,77) = 'PA2EQWZE1E2UA1:' +key(3,77) = 'AA2EQWZE1E2UA1:' +key(4,77) = 'SA2EQWZE1E2UA1:' + + +key(1,78) = 'NA2A1QU:' +key(2,78) = 'PA2A1QU:' +key(3,78) = 'AA2A1QU:' +key(4,78) = 'SA2A1QU:' + + +key(1,79) = 'NA2A1QUA1:' +key(2,79) = 'PA2A1QUA1:' +key(3,79) = 'AA2A1QUA1:' +key(4,79) = 'SA2A1QUA1:' + + +key(1,80) = 'NA2A1QUE1:' +key(2,80) = 'PA2A1QUE1:' +key(3,80) = 'AA2A1QUE1:' +key(4,80) = 'SA2A1QUE1:' + + +key(1,81) = 'NA2A1QUE2:' +key(2,81) = 'PA2A1QUE2:' +key(3,81) = 'AA2A1QUE2:' +key(4,81) = 'SA2A1QUE2:' + + +key(1,82) = 'NA2A1QUA1E1:' +key(2,82) = 'PA2A1QUA1E1:' +key(3,82) = 'AA2A1QUA1E1:' +key(4,82) = 'SA2A1QUA1E1:' + + +key(1,83) = 'NA2A1QUA1E2:' +key(2,83) = 'PA2A1QUA1E2:' +key(3,83) = 'AA2A1QUA1E2:' +key(4,83) = 'SA2A1QUA1E2:' + + +key(1,84) = 'NA2A1QUE1E2:' +key(2,84) = 'PA2A1QUE1E2:' +key(3,84) = 'AA2A1QUE1E2:' +key(4,84) = 'SA2A1QUE1E2:' + + +key(1,85) = 'NA2A1QUA1E1E2:' +key(2,85) = 'PA2A1QUA1E1E2:' +key(3,85) = 'AA2A1QUA1E1E2:' +key(4,85) = 'SA2A1QUA1E1E2:' + +key(1,86) = 'NCORECORE:' +key(2,86) = 'PCORECORE:' +key(3,86) = 'ACORECORE:' +key(4,86) = 'SCORECORE:' + + +end subroutine init_keys +end module diff --git a/src/model/new_PES/pes_ctrans.f90 b/src/model/new_PES/pes_ctrans.f90 new file mode 100644 index 0000000..d4b7ba5 --- /dev/null +++ b/src/model/new_PES/pes_ctrans.f90 @@ -0,0 +1,671 @@ +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % 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_pes_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 + +contains + subroutine ctrans_pes(q, s, t, invariants) + use dim_parameter, only: qn + integer(idp) k !running indices + real(dp), intent(in) :: q(qn) !given coordinates + real(dp), intent(out) :: s(qn) !output coordinates symmetry adapted and scaled + real(dp), intent(out) :: t(qn) !output coordinates symmetry adapted but not scaled +! ANN Variables + real(dp), optional, intent(out) :: invariants(:) + ! 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 = 0.0_dp + do k = 1, 3 + ch1(k) = q(k) + ch2(k) = q(k + 3) + ch3(k) = q(k + 6) + 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(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)) + call Hplane_pairdistances(ph1,ph2,ph3,t(7:9)) + + 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 + end subroutine ctrans_pes + + 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/) + + ! distance 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 + + subroutine Hplane_pairdistances(ph1,ph2,ph3, r) + real(dp), intent(in),dimension(3) :: ph1,ph2,ph3 + real(dp), intent(out) :: r(3) + real(dp) :: x(3) + x = ph1-ph2 + r(1) = norm(x) + x = ph2-ph3 + r(2) = norm(x) + x = ph3-ph1 + r(3) = norm(x) + end subroutine Hplane_pairdistances + + 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_pes_mod diff --git a/src/model/new_PES/pes_model.f90 b/src/model/new_PES/pes_model.f90 new file mode 100644 index 0000000..4acc849 --- /dev/null +++ b/src/model/new_PES/pes_model.f90 @@ -0,0 +1,982 @@ +module diab_pes + use dim_parameter, only: qn, ndiab, pst + use accuracy_constants, only: dp, idp + implicit none + logical :: debug = .false. + ! real(dp),parameter :: nuclear_energy_shift = 11.76027390_dp + real(dp), parameter :: nuclear_energy_shift = 0.881380_dp + real(dp), parameter :: ang2bohr = 1.0/0.52917721067_dp + real(dp), parameter :: a1_asymptote = 0.205024_dp*3._dp +!-------------------------------------------------------------------- +contains +!-------------------------------------------------------------------- + subroutine pote(e, n, q, s, t, p, npar) + integer(idp), intent(in) :: npar !number of parameters + integer(idp), intent(in) :: n + real(dp), intent(out) :: e(ndiab, ndiab) + real(dp), intent(in) :: q(qn), s(qn), t(qn) !< transformed coordinates + real(dp), intent(in), contiguous :: p(:) + call model_matrix(e, s, t, p) + end Subroutine + + !---------------------------------------------------------------------------------------------------- + ! calculate core repulsion potential depending on the shortest distance r and allow it only at values smaller than xmax + function core_core_rmin(r, p) result(v) + real(dp), intent(in) :: r(:) + real(dp), intent(in) :: p(:) + real(dp) :: v + real(dp) :: shift, width, rmin +if(size(p,1) /= 3) error stop 'Expected 3 parameters in core_core_min()' + rmin = minval(r) + v = p(1)/abs(rmin) + + shift = p(2) + width = p(3) + v = v*(1._dp - smootherstep(rmin, shift, width)) + end function core_core_rmin + + ! calculating the nuclear repulsion energy for a given set of pairwise distances and given nuclear charges + function nuclear_repulsion(r, charge) result(v) + real(dp), intent(in) :: r(:) + real(dp), intent(in) :: charge(:) + real(dp) :: v + v = sum(charge(:)/abs(r(:))) + end function nuclear_repulsion + + pure function smootherstep(x, shift, width) result(s) + real(dp), intent(in) :: x, shift, width + real(dp) :: s, q + q = (x - shift)/width + if (x <= 0._dp) then + s = 0._dp + else if (x >= 1._dp) then + s = 1._dp + else + s = 6._dp*q**5 - 15*q**4 + 10*q**3 + end if + end function smootherstep + + function pairwise_charge(charge) result(v) + real(dp), intent(in) :: charge(:) + real(dp), dimension(:), allocatable :: v + integer :: n, k, count + allocate (v(size(charge)*(size(charge) - 1)/2), source=0._dp) + count = 0 + do n = 1, size(charge, 1) + do k = n + 1, size(charge, 1) + count = count + 1 + v(count) = charge(n)*charge(k) + end do + end do + end function pairwise_charge + + function nuclear_repulsion_model(r, p) result(v) + real(dp), dimension(:), intent(in) :: r + real(dp), dimension(:), intent(in) :: p + real(dp) :: v + integer :: key, start, ende + key = 86 + if (pst(2, key) == 0) then + v = 0._dp + return + end if + start = pst(1, key) + ende = start + pst(2, key) - 1 + v = nuclear_repulsion(r(1:6)*ang2bohr, pairwise_charge(p(start+1:ende))) - p(start) + ! v = core_core_rmin(r([1, 2, 4]), p(start:ende)) + end function nuclear_repulsion_model + +!---------------------------------------------------------------------------------------------------- +! 3 + + case (2) + monom = v_aa_monom(a, b, 3)*w_e_monom(x, y, 1) ! order 3,1 => 4 + case (3) + monom = v_aa_monom(a, b, 2)*w_e_monom(x, y, 1) ! order 3,1 => 4 + case (4) + monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 2) ! order 2,2 => 4 + + case (5) + monom = v_aa_monom(a, b, 6)*w_e_monom(x, y, 1) ! order 4,1 => 5 + case (6) + monom = v_aa_monom(a, b, 5)*w_e_monom(x, y, 1) ! order 4,1 => 5 + case (7) + monom = v_aa_monom(a, b, 4)*w_e_monom(x, y, 1) ! order 4,1 => 5 + + case (8) + monom = v_aa_monom(a, b, 3)*w_e_monom(x, y, 2) ! order 3,2 => 5 + case (9) + monom = v_aa_monom(a, b, 2)*w_e_monom(x, y, 2) ! order 3,2 => 5 + case (10) + monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 2) ! order 2,3 => 5 + + case (11) + monom = v_aa_monom(a, b, 10)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (12) + monom = v_aa_monom(a, b, 9)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (13) + monom = v_aa_monom(a, b, 8)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (14) + monom = v_aa_monom(a, b, 7)*w_e_monom(x, y, 1) ! order 5,1 => 6 + + case (15) + monom = v_aa_monom(a, b, 6)*w_e_monom(x, y, 2) ! order 4,2 => 6 + case (16) + monom = v_aa_monom(a, b, 5)*w_e_monom(x, y, 2) ! order 4,2 => 6 + case (17) + monom = v_aa_monom(a, b, 4)*w_e_monom(x, y, 2) ! order 4,2 => 6 + + case (18) + monom = v_aa_monom(a, b, 3)*w_e_monom(x, y, 3) ! order 3,3 => 6 + case (19) + monom = v_aa_monom(a, b, 2)*w_e_monom(x, y, 3) ! order 3,3 => 6 + case (20) + monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 4) ! order 2,4 => 6 + case (21) + monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 5) ! order 2,4 => 6 + + case default + error stop 'called case of w_aae_monom not implemented' + monom = 0.0_dp + end select + end function w_aae_monom + + function z_aae_monom(a, b, x, y, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: a, b, x, y + integer(kind=idp), intent(in) :: o + select case (o) + case (1) + monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 1) ! order 2,1 => 3 + + case (2) + monom = v_aa_monom(a, b, 3)*z_e_monom(x, y, 1) ! order 3,1 => 4 + case (3) + monom = v_aa_monom(a, b, 2)*z_e_monom(x, y, 1) ! order 3,1 => 4 + case (4) + monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 2) ! order 2,2 => 4 + + case (5) + monom = v_aa_monom(a, b, 6)*z_e_monom(x, y, 1) ! order 4,1 => 5 + case (6) + monom = v_aa_monom(a, b, 5)*z_e_monom(x, y, 1) ! order 4,1 => 5 + case (7) + monom = v_aa_monom(a, b, 4)*z_e_monom(x, y, 1) ! order 4,1 => 5 + + case (8) + monom = v_aa_monom(a, b, 3)*z_e_monom(x, y, 2) ! order 3,2 => 5 + case (9) + monom = v_aa_monom(a, b, 2)*z_e_monom(x, y, 2) ! order 3,2 => 5 + case (10) + monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 2) ! order 2,3 => 5 + + case (11) + monom = v_aa_monom(a, b, 10)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (12) + monom = v_aa_monom(a, b, 9)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (13) + monom = v_aa_monom(a, b, 8)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (14) + monom = v_aa_monom(a, b, 7)*z_e_monom(x, y, 1) ! order 5,1 => 6 + + case (15) + monom = v_aa_monom(a, b, 6)*z_e_monom(x, y, 2) ! order 4,2 => 6 + case (16) + monom = v_aa_monom(a, b, 5)*z_e_monom(x, y, 2) ! order 4,2 => 6 + case (17) + monom = v_aa_monom(a, b, 4)*z_e_monom(x, y, 2) ! order 4,2 => 6 + + case (18) + monom = v_aa_monom(a, b, 3)*z_e_monom(x, y, 3) ! order 3,3 => 6 + case (19) + monom = v_aa_monom(a, b, 2)*z_e_monom(x, y, 3) ! order 3,3 => 6 + case (20) + monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 4) ! order 2,4 => 6 + case (21) + monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 5) ! order 2,4 => 6 + + case default + error stop 'called case of z_aae_monom not implemented' + monom = 0.0_dp + end select + end function z_aae_monom + + function w_ee_monom(x1, y1, x2, y2, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: x1, y1, x2, y2 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) ! order 2 + monom = x1*x2 - y1*y2 + + case (2) ! order 3 + monom = x1**2*x2 + x2*y1**2 + case (3) + monom = x1*y2**2 + x1*x2**2 + case (4) + monom = 2.0_dp*x1*y1*y2 + x1**2*x2 - x2*y1**2 + case (5) + monom = x1*x2**2 + 2.0_dp*x2*y1*y2 - x1*y2**2 + + case (6) ! order 4 + monom = x1**3*x2 - 3.0_dp*x1**2*y1*y2& + &- 3.0_dp*x1*x2*y1**2 + y1**3*y2 + case (7) + monom = x1**2*x2**2 - x1**2*y2**2 - x2**2*y1**2 + y1**2*y2**2& + &- 4.0_dp*x1*x2*y1*y2 + case (8) + monom = x2**3*x1 - 3.0_dp*x1*x2*y2**2& + &- 3.0_dp*x2**2*y1*y2 + y2**3*y1 + case (9) + monom = x1**3*x2 + 3.0_dp*x1**2*y1*y2& + &- 3.0_dp*x1*x2*y1**2 - y1**3*y2 + case (10) + monom = x2**3*x1 - 3.0_dp*x1*x2*y2**2& + &+ 3.0_dp*x2**2*y1*y2 - y2**3*y1 + case (11) + monom = x1**3*x2 - x1**2*y1*y2 + x1*x2*y1**2 - y1**3*y2 + case (12) + monom = x1**2*x2**2 + x1**2*y2**2 - x2**2*y1**2& + &- y1**2*y2**2 + case (13) + monom = x1**2*x2**2 - x1**2*y2**2 + x2**2*y1**2& + &- y1**2*y2**2 + case (14) + monom = x1*x2**3 + x1*x2*y2**2 - x2**2*y1*y2 - y2**3*y1 + + case (15) ! order 5 + monom = (x1*y2**4 + 4.0_dp*x2*y1*y2**3 - 6.0_dp*x1*x2**2*y2**2& + &- 4.0_dp*x2**3*y1*y2 + x1*x2**4) + case (16) + monom = (2.0_dp*x1*y1*y2**3 + 3.0_dp*x2*y1**2*y2**2& + &- 3.0_dp*x1**2*x2*y2**2 - 6.0_dp*x1*x2**2*y1*y2& + &- x2**3*y1**2 + x1**2*x2**3) + case (17) + monom = (3.0_dp*x1*y1**2*y2**2 - x1**3*y2**2& + &+ 2.0_dp*x2*y1**3*y2& + &- 6.0_dp*x1**2*x2*y1*y2 - 3.0_dp*x1*x2**2*y1**2& + &+ x1**3*x2**2) + case (18) + monom = (4.0_dp*x1*y1**3*y2 - 4.0_dp*x1**3*y1*y2 + x2*y1**4& + &- 6.0_dp*x1**2*x2*y1**2 + x1**4*x2) + case (19) + monom = -(y2**2 + x2**2)*(x1*y2**2 - 2.0_dp*x2*y1*y2 - x1*x2**2) + case (20) + monom = x1*(y2**2 + x2**2)**2 + case (21) + monom = -(2.0_dp*x1*y1*y2**3 - 3.0_dp*x2*y1**2*y2**2& + &+ 3.0_dp*x1**2& + &*x2*y2**2 - 6.0_dp*x1*x2**2*y1*y2& + &+ x2**3*y1**2 - x1**2*x2**3) + case (22) + monom = x2*(y1**2 + x1**2)*(y2**2 + x2**2) + case (23) + monom = -(y1**2 + x1**2)*(x1*y2**2 - 2.0_dp*x2*y1*y2 - x1*x2**2) + case (24) + monom = x1*(y1**2 + x1**2)*(y2**2 + x2**2) + case (25) + monom = x2*(y1**2 + x1**2)**2 + case (26) + monom = (y1**2 + x1**2)*(2.0_dp*x1*y1*y2 - x2*y1**2 + x1**2*x2) + + case (27) ! order 6 + monom = (y1*y2**5 + 5.0_dp*x1*x2*y2**4 - 10.0_dp*x2**2*y1*y2**3& + &- 10.0_dp*x1*x2**3*y2**2 + 5.0_dp*x2**4*y1*y2 + x1*x2**5) + case (28) + monom = (y2**2 + x2**2)*(y1*y2**3 - 3.0_dp*x1*x2*y2**2& + &- 3.0_dp*x2**2*y1*y2 + x1*x2**3) + case (29) + monom = (y1**2 + x1**2)*(y2**2 - 2.0_dp*x2*y2 - x2**2)& + &*(y2**2 + 2.0_dp*x2*y2 - x2**2) + case (30) + monom = (y1*y2 - x1*y2 - x2*y1 - x1*x2)*(y1*y2 + x1*y2& + &+ x2*y1 - x1*x2)*(y2**2 + x2**2) + case (31) + monom = (y1**2 + x1**2)*(y1*y2**3 - 3.0_dp*x1*x2*y2**2& + &- 3.0_dp*x2**2*y1*y2 + x1*x2**3) + case (32) + monom = (y1**3*y2 - 3.0_dp*x1**2*y1*y2& + &- 3.0_dp*x1*x2*y1**2 + x1**3*x2)*(y2**2 + x2**2) + case (33) + monom = (y1**2 + x1**2)*(y1*y2 - x1*y2 - x2*y1 - x1*x2)& + &*(y1*y2 + x1*y2 + x2*y1 - x1*x2) + case (34) + monom = (y1**2 - 2.0_dp*x1*y1 - x1**2)& + &*(y1**2 + 2.0_dp*x1*y1 - x1**2)& + &*(y2**2 + x2**2) + case (35) + monom = (y1**2 + x1**2)*(y1**3*y2 - 3.0_dp*x1**2*y1*y2& + &- 3.0_dp*x1*x2*y1**2 + x1**3*x2) + case (36) + monom = (y1**5*y2 - 10.0_dp*x1**2*y1**3*y2 + 5.0_dp*x1**4*y1*y2& + &+ 5.0_dp*x1*x2*y1**4 - 10.0_dp*x1**3*x2*y1**2 + x1**5*x2) + case (37) + monom = (y1*y2 - x1*x2)*(y2**2 + x2**2)**2 + case (38) + monom = (y2**2 + x2**2)*(y1*y2**3 + 3.0_dp*x1*x2*y2**2& + &- 3.0_dp*x2**2*y1*y2 - x1*x2**3) + case (39) + monom = (y1**2 + x1**2)*(y2 - x2)*(y2 + x2)*(y2**2 + x2**2) + case (40) + monom = (y1 - x1)*(y1 + x1)*(y2**2 + x2**2)**2 + case (41) + monom = (y1*y2**2 - x1*y2**2 + 2.0_dp*x2*y1*y2 + 2.0_dp*x1*x2*y2& + &- x2**2*y1 + x1*x2**2)*(y1*y2**2 + x1*y2**2 - 2.0_dp*x2*y1*y2& + &+ 2.0_dp*x1*x2*y2 - x2**2*y1 - x1*x2**2) + case (42) + monom = (y1**3*y2 - 3.0_dp*x1**2*y1*y2& + &+ 3.0_dp*x1*x2*y1**2 - x1**3*x2)*(y2**2 + x2**2) + case (43) + monom = (y1**2 + x1**2)*(y1*y2 - x1*x2)*(y2**2 + x2**2) + case (44) + monom = (y1**2 + x1**2)*(y1*y2**3 + 3.0_dp*x1*x2*y2**2& + &- 3.0_dp*x2**2*y1*y2 - x1*x2**3) + case (45) + monom = (y1**2*y2 - 2.0_dp*x1*y1*y2 - x1**2*y2 + x2*y1**2& + &+ 2.0_dp*x1*x2*y1 - x1**2*x2)& + &*(y1**2*y2 + 2.0_dp*x1*y1*y2 - x1**2*y2& + &- x2*y1**2 + 2.0_dp*x1*x2*y1 + x1**2*x2) + case (46) + monom = (y1 - x1)*(y1 + x1)*(y1**2 + x1**2)*(y2**2 + x2**2) + case (47) + monom = (y1**2 + x1**2)**2*(y2 - x2)*(y2 + x2) + case (48) + monom = (y1**2 + x1**2)*(y1**3*y2 - 3.0_dp*x1**2*y1*y2& + &+ 3.0_dp*x1*x2*y1**2 - x1**3*x2) + case (49) + monom = (y1**2 + x1**2)**2*(y1*y2 - x1*x2) + case default + error stop 'called case of w_ee_monom not implemented' + monom = 0.0_dp + end select + end function w_ee_monom + + function z_ee_monom(x1, y1, x2, y2, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: x1, y1, x2, y2 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) ! order 2 + monom = -x1*y2 - x2*y1 + + case (2) ! order 3 + monom = y1**2*y2 + x1**2*y2 + case (3) + monom = y1*y2**2 + x2**2*y1 + case (4) + monom = y1**2*y2 + 2.0_dp*x1*x2*y1 - x1**2*y2 + case (5) + monom = y2**2*y1 + 2.0_dp*x1*x2*y2 - x2**2*y1 + + case (6) ! order 4 + monom = x1**3*y2 + 3.0_dp*x1**2*x2*y1& + &- 3.0_dp*x1*y1**2*y2 - x2*y1**3 + case (7) + monom = 2.0_dp*x1**2*x2*y2 + 2.0_dp*x1*x2**2*y1& + &- 2.0_dp*x1*y1*y2**2 - 2.0_dp*x2*y1**2*y2 + case (8) + monom = x2**3*y1 + 3.0_dp*x1*x2**2*y2& + &- 3.0_dp*x2*y1*y2**2 - x1*y2**3 + case (9) + monom = x1**3*y2 - 3.0_dp*x1**2*x2*y1 - 3.0_dp*x1*y1**2*y2& + &+ x2*y1**3 + case (10) + monom = x2**3*y1 - 3.0_dp*x1*x2**2*y2 - 3.0_dp*x2*y1*y2**2& + &+ x1*y2**3 + case (11) + monom = -x1**3*y2 - x1**2*x2*y1 - x1*y1**2*y2 - x2*y1**3 + case (12) + monom = -2.0_dp*x1*x2**2*y1 - 2.0_dp*x1*y1*y2**2 + case (13) + monom = -2.0_dp*x1**2*x2*y2 - 2.0_dp*x2*y1**2*y2 + case (14) + monom = -y2**3*x1 - x1*x2**2*y2 - x2*y1*y2**2 - y1*x2**3 + + case (15) ! order 5 + monom = -(y1*y2**4 - 4.0_dp*x1*x2*y2**3 - 6.0_dp*x2**2*y1*y2**2& + &+ 4.0_dp*x1*x2**3*y2 + x2**4*y1) + case (16) + monom = -(y1**2*y2**3 - x1**2*y2**3 - 6.0_dp*x1*x2*y1*y2**2& + &- 3.0_dp*x2**2*y1**2*y2 + 3.0_dp*x1**2*x2**2*y2& + &+ 2.0_dp*x1*x2**3*y1) + case (17) + monom = -(y1**3*y2**2 - 3.0_dp*x1**2*y1*y2**2& + &- 6.0_dp*x1*x2*y1**2*y2 + 2.0_dp*x1**3*x2*y2 - x2**2*y1**3& + &+ 3.0_dp*x1**2*x2**2*y1) + case (18) + monom = -(y1**4*y2 - 6.0_dp*x1**2*y1**2*y2 + x1**4*y2& + &- 4.0_dp*x1*x2*y1**3 + 4.0_dp*x1**3*x2*y1) + case (19) + monom = (y2**2 + x2**2)*(y1*y2**2 + 2.0_dp*x1*x2*y2 - x2**2*y1) + case (20) + monom = y1*(y2**2 + x2**2)**2 + case (21) + monom = (y1**2*y2**3 - x1**2*y2**3 + 6.0_dp*x1*x2*y1*y2**2& + &- 3.0_dp*x2**2*y1**2*y2 + 3.0_dp*x1**2*x2**2*y2& + &- 2.0_dp*x1*x2**3*y1) + case (22) + monom = (y1**2 + x1**2)*y2*(y2**2 + x2**2) + case (23) + monom = (y1**2 + x1**2)*(y1*y2**2 + 2.0_dp*x1*x2*y2 - x2**2*y1) + case (24) + monom = y1*(y1**2 + x1**2)*(y2**2 + x2**2) + case (25) + monom = (y1**2 + x1**2)**2*y2 + case (26) + monom = (y1**2 + x1**2)*(y1**2*y2 - x1**2*y2 + 2.0_dp*x1*x2*y1) + + case (27) ! order 6 + monom = (x1*y2**5 - 5.0_dp*x2*y1*y2**4 - 10.0_dp*x1*x2**2*y2**3& + &+ 10.0_dp*x2**3*y1*y2**2 + 5.0_dp*x1*x2**4*y2 - x2**5*y1) + case (28) + monom = -(y2**2 + x2**2)*(x1*y2**3 + 3.0_dp*x2*y1*y2**2& + &- 3.0_dp*x1*x2**2*y2 - x2**3*y1) + case (29) + monom = -4.0_dp*x2*(y1**2 + x1**2)*y2*(y2 - x2)*(y2 + x2) + case (30) + monom = -2.0_dp*(x1*y2 + x2*y1)*(y1*y2 - x1*x2)*(y2**2 + x2**2) + case (31) + monom = -(y1**2 + x1**2)*(x1*y2**3 + 3.0_dp*x2*y1*y2**2& + &- 3.0_dp*x1*x2**2*y2 - x2**3*y1) + case (32) + monom = -(3.0_dp*x1*y1**2*y2 - x1**3*y2 + x2*y1**3& + &- 3.0_dp*x1**2*x2*y1)*(y2**2 + x2**2) + case (33) + monom = -2.0_dp*(y1**2 + x1**2)*(x1*y2 + x2*y1)*(y1*y2 - x1*x2) + case (34) + monom = -4.0_dp*x1*y1*(y1 - x1)*(y1 + x1)*(y2**2 + x2**2) + case (35) + monom = -(y1**2 + x1**2)*(3.0_dp*x1*y1**2*y2 - x1**3*y2 + x2*y1**3& + &- 3.0_dp*x1**2*x2*y1) + case (36) + monom = -(5.0_dp*x1*y1**4*y2& + &- 10.0_dp*x1**3*y1**2*y2 + x1**5*y2& + &- x2*y1**5 + 10.0_dp*x1**2*x2*y1**3 - 5.0_dp*x1**4*x2*y1) + case (37) + monom = (x1*y2 + x2*y1)*(y2**2 + x2**2)**2 + case (38) + monom = -(y2**2 + x2**2)*(x1*y2**3 - 3.0_dp*x2*y1*y2**2& + &- 3.0_dp*x1*x2**2*y2 + x2**3*y1) + case (39) + monom = 2.0_dp*x1*y1*(y2**2 + x2**2)**2 + case (40) + monom = 2.0_dp*x2*(y1**2 + x1**2)*y2*(y2**2 + x2**2) + case (41) + monom = -2.0_dp*(x1*y2**2 - 2.0_dp*x2*y1*y2 - x1*x2**2)& + &*(y1*y2**2 + 2.0_dp*x1*x2*y2 - x2**2*y1) + case (42) + monom = (3.0_dp*x1*y1**2*y2 - x1**3*y2 - x2*y1**3& + &+ 3.0_dp*x1**2*x2*y1)*(y2**2 + x2**2) + case (43) + monom = (y1**2 + x1**2)*(x1*y2 + x2*y1)*(y2**2 + x2**2) + case (44) + monom = -(y1**2 + x1**2)*(x1*y2**3 - 3.0_dp*x2*y1*y2**2& + &- 3.0_dp*x1*x2**2*y2 + x2**3*y1) + case (45) + monom = 2.0_dp*(2.0_dp*x1*y1*y2 - x2*y1**2 + x1**2*x2)& + &*(y1**2*y2 - x1**2*y2 + 2.0_dp*x1*x2*y1) + case (46) + monom = 2.0_dp*x1*y1*(y1**2 + x1**2)*(y2**2 + x2**2) + case (47) + monom = 2.0_dp*x2*(y1**2 + x1**2)**2*y2 + case (48) + monom = (y1**2 + x1**2)*(3.0_dp*x1*y1**2*y2 - x1**3*y2 - x2*y1**3& + &+ 3.0_dp*x1**2*x2*y1) + case (49) + monom = (y1**2 + x1**2)**2*(x1*y2 + x2*y1) + case default + error stop 'called case of z_ee_monom not implemented' + monom = 0.0_dp + end select + end function z_ee_monom + + function w_aee_monom(a, x1, y1, x2, y2, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: a, x1, y1, x2, y2 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) + monom = a*w_ee_monom(x1, y1, x2, y2, 1) ! order 3 + + case (2) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 1) ! order 4 + case (3) + monom = a*w_ee_monom(x1, y1, x2, y2, 2) ! order 4 + case (4) + monom = a*w_ee_monom(x1, y1, x2, y2, 3) ! order 4 + case (5) + monom = a*w_ee_monom(x1, y1, x2, y2, 4) ! order 4 + case (6) + monom = a*w_ee_monom(x1, y1, x2, y2, 5) ! order 4 + + case (7) + monom = a**3*w_ee_monom(x1, y1, x2, y2, 1) ! order 5 + case (8) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 2) ! order 5 + case (9) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 3) ! order 5 + case (10) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 4) ! order 5 + case (11) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 5) ! order 5 + case (12) + monom = a*w_ee_monom(x1, y1, x2, y2, 6) ! order 5 + case (13) + monom = a*w_ee_monom(x1, y1, x2, y2, 7) ! order 5 + case (14) + monom = a*w_ee_monom(x1, y1, x2, y2, 8) ! order 5 + case (15) + monom = a*w_ee_monom(x1, y1, x2, y2, 9) ! order 5 + case (16) + monom = a*w_ee_monom(x1, y1, x2, y2, 10) ! order 5 + case (17) + monom = a*w_ee_monom(x1, y1, x2, y2, 11) ! order 5 + case (18) + monom = a*w_ee_monom(x1, y1, x2, y2, 12) ! order 5 + case (19) + monom = a*w_ee_monom(x1, y1, x2, y2, 13) ! order 5 + case (20) + monom = a*w_ee_monom(x1, y1, x2, y2, 14) ! order 5 + + case (21) + monom = a**4*w_ee_monom(x1, y1, x2, y2, 1) ! order 6 + case (22) + monom = a**3*w_ee_monom(x1, y1, x2, y2, 2) ! order 6 + case (23) + monom = a**3*w_ee_monom(x1, y1, x2, y2, 3) ! order 6 + case (24) + monom = a**3*w_ee_monom(x1, y1, x2, y2, 4) ! order 6 + case (25) + monom = a**3*w_ee_monom(x1, y1, x2, y2, 5) ! order 6 + case (26) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 6) ! order 6 + case (27) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 7) ! order 6 + case (28) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 8) ! order 6 + case (29) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 9) ! order 6 + case (30) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 10) ! order 6 + case (31) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 11) ! order 6 + case (32) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 12) ! order 6 + case (33) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 13) ! order 6 + case (34) + monom = a**2*w_ee_monom(x1, y1, x2, y2, 14) ! order 6 + case (35) + monom = a*w_ee_monom(x1, y1, x2, y2, 15) ! order 6 + case (36) + monom = a*w_ee_monom(x1, y1, x2, y2, 16) ! order 6 + case (37) + monom = a*w_ee_monom(x1, y1, x2, y2, 17) ! order 6 + case (38) + monom = a*w_ee_monom(x1, y1, x2, y2, 18) ! order 6 + case (39) + monom = a*w_ee_monom(x1, y1, x2, y2, 19) ! order 6 + case (40) + monom = a*w_ee_monom(x1, y1, x2, y2, 20) ! order 6 + case (41) + monom = a*w_ee_monom(x1, y1, x2, y2, 21) ! order 6 + case (42) + monom = a*w_ee_monom(x1, y1, x2, y2, 22) ! order 6 + case (43) + monom = a*w_ee_monom(x1, y1, x2, y2, 23) ! order 6 + case (44) + monom = a*w_ee_monom(x1, y1, x2, y2, 24) ! order 6 + case (45) + monom = a*w_ee_monom(x1, y1, x2, y2, 25) ! order 6 + case (46) + monom = a*w_ee_monom(x1, y1, x2, y2, 26) ! order 6 + + case default + error stop 'called case of w_aee_monom not implemented' + monom = 0.0_dp + end select + end function w_aee_monom + + function z_aee_monom(a, x1, y1, x2, y2, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: a, x1, y1, x2, y2 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) + monom = a*z_ee_monom(x1, y1, x2, y2, 1) ! order 3 + + case (2) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 1) ! order 4 + case (3) + monom = a*z_ee_monom(x1, y1, x2, y2, 2) ! order 4 + case (4) + monom = a*z_ee_monom(x1, y1, x2, y2, 3) ! order 4 + case (5) + monom = a*z_ee_monom(x1, y1, x2, y2, 4) ! order 4 + case (6) + monom = a*z_ee_monom(x1, y1, x2, y2, 5) ! order 4 + + case (7) + monom = a**3*z_ee_monom(x1, y1, x2, y2, 1) ! order 5 + case (8) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 2) ! order 5 + case (9) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 3) ! order 5 + case (10) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 4) ! order 5 + case (11) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 5) ! order 5 + case (12) + monom = a*z_ee_monom(x1, y1, x2, y2, 6) ! order 5 + case (13) + monom = a*z_ee_monom(x1, y1, x2, y2, 7) ! order 5 + case (14) + monom = a*z_ee_monom(x1, y1, x2, y2, 8) ! order 5 + case (15) + monom = a*z_ee_monom(x1, y1, x2, y2, 9) ! order 5 + case (16) + monom = a*z_ee_monom(x1, y1, x2, y2, 10) ! order 5 + case (17) + monom = a*z_ee_monom(x1, y1, x2, y2, 11) ! order 5 + case (18) + monom = a*z_ee_monom(x1, y1, x2, y2, 12) ! order 5 + case (19) + monom = a*z_ee_monom(x1, y1, x2, y2, 13) ! order 5 + case (20) + monom = a*z_ee_monom(x1, y1, x2, y2, 14) ! order 5 + + case (21) + monom = a**4*z_ee_monom(x1, y1, x2, y2, 1) ! order 6 + case (22) + monom = a**3*z_ee_monom(x1, y1, x2, y2, 2) ! order 6 + case (23) + monom = a**3*z_ee_monom(x1, y1, x2, y2, 3) ! order 6 + case (24) + monom = a**3*z_ee_monom(x1, y1, x2, y2, 4) ! order 6 + case (25) + monom = a**3*z_ee_monom(x1, y1, x2, y2, 5) ! order 6 + case (26) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 6) ! order 6 + case (27) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 7) ! order 6 + case (28) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 8) ! order 6 + case (29) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 9) ! order 6 + case (30) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 10) ! order 6 + case (31) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 11) ! order 6 + case (32) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 12) ! order 6 + case (33) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 13) ! order 6 + case (34) + monom = a**2*z_ee_monom(x1, y1, x2, y2, 14) ! order 6 + case (35) + monom = a*z_ee_monom(x1, y1, x2, y2, 15) ! order 6 + case (36) + monom = a*z_ee_monom(x1, y1, x2, y2, 16) ! order 6 + case (37) + monom = a*z_ee_monom(x1, y1, x2, y2, 17) ! order 6 + case (38) + monom = a*z_ee_monom(x1, y1, x2, y2, 18) ! order 6 + case (39) + monom = a*z_ee_monom(x1, y1, x2, y2, 19) ! order 6 + case (40) + monom = a*z_ee_monom(x1, y1, x2, y2, 20) ! order 6 + case (41) + monom = a*z_ee_monom(x1, y1, x2, y2, 21) ! order 6 + case (42) + monom = a*z_ee_monom(x1, y1, x2, y2, 22) ! order 6 + case (43) + monom = a*z_ee_monom(x1, y1, x2, y2, 23) ! order 6 + case (44) + monom = a*z_ee_monom(x1, y1, x2, y2, 24) ! order 6 + case (45) + monom = a*z_ee_monom(x1, y1, x2, y2, 25) ! order 6 + case (46) + monom = a*z_ee_monom(x1, y1, x2, y2, 26) ! order 6 + + case default + error stop 'called case of z_aee_monom not implemented' + monom = 0.0_dp + end select + end function z_aee_monom + + function w_aaae_monom(a, b, c, x, y, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: a, b, c, x, y + integer(kind=idp), intent(in) :: o + select case (o) + case (1) + monom = v_aaa_monom(a, b, c, 1)*w_e_monom(x, y, 1) ! order 3,1 => 4 + + case (2) + monom = v_aaa_monom(a, b, c, 4)*w_e_monom(x, y, 1) ! order 4,1 => 5 + case (3) + monom = v_aaa_monom(a, b, c, 3)*w_e_monom(x, y, 1) ! order 4,1 => 5 + case (4) + monom = v_aaa_monom(a, b, c, 2)*w_e_monom(x, y, 1) ! order 4,1 => 5 + case (5) + monom = v_aaa_monom(a, b, c, 1)*w_e_monom(x, y, 2) ! order 3,2 => 5 + + case (6) + monom = v_aaa_monom(a, b, c, 10)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (7) + monom = v_aaa_monom(a, b, c, 9)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (8) + monom = v_aaa_monom(a, b, c, 8)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (9) + monom = v_aaa_monom(a, b, c, 7)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (10) + monom = v_aaa_monom(a, b, c, 6)*w_e_monom(x, y, 1) ! order 5,1 => 6 + case (11) + monom = v_aaa_monom(a, b, c, 5)*w_e_monom(x, y, 1) ! order 5,1 => 6 + + case (12) + monom = v_aaa_monom(a, b, c, 4)*w_e_monom(x, y, 2) ! order 4,2 => 6 + case (13) + monom = v_aaa_monom(a, b, c, 3)*w_e_monom(x, y, 2) ! order 4,2 => 6 + case (14) + monom = v_aaa_monom(a, b, c, 2)*w_e_monom(x, y, 2) ! order 4,2 => 6 + case (15) + monom = v_aaa_monom(a, b, c, 1)*w_e_monom(x, y, 3) ! order 3,3 => 6 + + case default + error stop 'called case of w_aaae_monom not implemented' + monom = 0.0_dp + end select + end function w_aaae_monom + + function z_aaae_monom(a, b, c, x, y, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: a, b, c, x, y + integer(kind=idp), intent(in) :: o + select case (o) + case (1) + monom = v_aaa_monom(a, b, c, 1)*z_e_monom(x, y, 1) ! order 3,1 => 4 + + case (2) + monom = v_aaa_monom(a, b, c, 4)*z_e_monom(x, y, 1) ! order 4,1 => 5 + case (3) + monom = v_aaa_monom(a, b, c, 3)*z_e_monom(x, y, 1) ! order 4,1 => 5 + case (4) + monom = v_aaa_monom(a, b, c, 2)*z_e_monom(x, y, 1) ! order 4,1 => 5 + case (5) + monom = v_aaa_monom(a, b, c, 1)*z_e_monom(x, y, 2) ! order 3,2 => 5 + + case (6) + monom = v_aaa_monom(a, b, c, 10)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (7) + monom = v_aaa_monom(a, b, c, 9)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (8) + monom = v_aaa_monom(a, b, c, 8)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (9) + monom = v_aaa_monom(a, b, c, 7)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (10) + monom = v_aaa_monom(a, b, c, 6)*z_e_monom(x, y, 1) ! order 5,1 => 6 + case (11) + monom = v_aaa_monom(a, b, c, 5)*z_e_monom(x, y, 1) ! order 5,1 => 6 + + case (12) + monom = v_aaa_monom(a, b, c, 4)*z_e_monom(x, y, 2) ! order 4,2 => 6 + case (13) + monom = v_aaa_monom(a, b, c, 3)*z_e_monom(x, y, 2) ! order 4,2 => 6 + case (14) + monom = v_aaa_monom(a, b, c, 2)*z_e_monom(x, y, 2) ! order 4,2 => 6 + case (15) + monom = v_aaa_monom(a, b, c, 1)*z_e_monom(x, y, 3) ! order 3,3 => 6 + + case default + error stop 'called case of z_aaae_monom not implemented' + monom = 0.0_dp + end select + end function z_aaae_monom + + function w_eee_monom(x1, y1, x2, y2, x3, y3, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: x1, y1, x2, y2, x3, y3 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) !order 3 + monom = x1*x2*x3 + y1*y2*x3 + y1*x2*y3 - x1*y2*y3 + case (2) + monom = x1*x2*x3 + y1*y2*x3 - y1*x2*y3 + x1*y2*y3 + case (3) + monom = x1*x2*x3 - y1*y2*x3 + y1*x2*y3 + x1*y2*y3 + + case (4) !order 4 + monom = (x1**2)*x2*x3 - (y1**2)*x2*x3 - 2.d0*x1*y1*y2*x3& + &- 2.0_dp*x1*y1*x2*y3 - (x1**2)*y2*y3 + (y1**2)*y2*y3 + case (5) + monom = (x1**2)*x2*x3 + (y1**2)*x2*x3 - (x1**2)*y2*y3& + &- (y1**2)*y2*y3 + case (6) + monom = (x1**2)*x2*x3 - (y1**2)*x2*x3 + 2.0_dp*x1*y1*y2*x3& + &- 2.0_dp*x1*y1*x2*y3 + (x1**2)*y2*y3 - (y1**2)*y2*y3 + case (7) + monom = (x1**2)*x2*x3 - (y1**2)*x2*x3 - 2.0_dp*x1*y1*y2*x3& + &+ 2.0_dp*x1*y1*x2*y3 + (x1**2)*y2*y3 - (y1**2)*y2*y3 + case (8) + + monom = (x2**2)*x1*x3 - (y2**2)*x1*x3 - 2.0_dp*x2*y2*y1*x3& + &- 2.0_dp*x2*y2*x1*y3 - (x2**2)*y1*y3 + (y2**2)*y1*y3 + case (9) + monom = (x2**2)*x1*x3 + (y2**2)*x1*x3 - (x2**2)*y1*y3& + &- (y2**2)*y1*y3 + case (10) + monom = (x2**2)*x1*x3 - (y2**2)*x1*x3 + 2.0_dp*x2*y2*y1*x3& + &- 2.0_dp*x2*y2*x1*y3 + (x2**2)*y1*y3 - (y2**2)*y1*y3 + case (11) + monom = (x2**2)*x1*x3 - (y2**2)*x1*x3 - 2.0_dp*x2*y2*y1*x3& + &+ 2.0_dp*x2*y2*x1*y3 + (x2**2)*y1*y3 - (y2**2)*y1*y3 + case (12) + monom = (x3**2)*x2*x1 - (y3**2)*x2*x1 - 2.0_dp*x3*y3*y2*x1& + &- 2.0_dp*x3*y3*x2*y1 - (x3**2)*y2*y1 + (y3**2)*y2*y1 + case (13) + monom = (x3**2)*x2*x1 + (y3**2)*x2*x1 - (x3**2)*y2*y1& + &- (y3**2)*y2*y1 + case (14) + monom = (x3**2)*x2*x1 - (y3**2)*x2*x1 + 2.0_dp*x3*y3*y2*x1& + &- 2.0_dp*x3*y3*x2*y1 + (x3**2)*y2*y1 - (y3**2)*y2*y1 + case (15) + monom = (x3**2)*x2*x1 - (y3**2)*x2*x1 - 2.0_dp*x3*y3*y2*x1& + &+ 2.0_dp*x3*y3*x2*y1 + (x3**2)*y2*y1 - (y3**2)*y2*y1 + case default + error stop 'called case of w_eee_monom not implemented' + monom = 0.0_dp + end select + end function w_eee_monom + + function z_eee_monom(x1, y1, x2, y2, x3, y3, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: x1, y1, x2, y2, x3, y3 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) ! order 3 + monom = -y1*x2*x3 + x1*y2*x3 + x1*x2*y3 + y1*y2*y3 + case (2) + monom = y1*x2*x3 - x1*y2*x3 + x1*x2*y3 + y1*y2*y3 + case (3) + monom = y1*x2*x3 + x1*y2*x3 - x1*x2*y3 + y1*y2*y3 + + case (4) ! order 4 + monom = (2.0*x1*y1*x2*x3 + (x1**2)*y2*x3 + (x1**2)*x2*y3& + &- 2.0*x1*y1*y2*y3 - (y1**2)*y2*x3 - (y1**2)*x2*y3) + case (5) + monom = (-(x1**2)*y2*x3 - (x1**2)*x2*y3 - (y1**2)*y2*x3& + &- (y1**2)*x2*y3) + case (6) + monom = -(2.0*x1*y1*x2*x3 - (x1**2)*y2*x3 + (x1**2)*x2*y3& + &+ 2.0*x1*y1*y2*y3 + (y1**2)*y2*x3 - (y1**2)*x2*y3) + case (7) + monom = -(2.0*x1*y1*x2*x3 + (x1**2)*y2*x3 - (x1**2)*x2*y3& + &+ 2.0*x1*y1*y2*y3 - (y1**2)*y2*x3 + (y1**2)*x2*y3) + + case (8) + monom = (2.d0*x2*y2*x1*x3 + (x2**2)*y1*x3 + (x2**2)*x1*y3& + &- 2.d0*x2*y2*y1*y3 - (y2**2)*y1*x3 - (y2**2)*x1*y3) + case (9) + monom = (-(x2**2)*y1*x3 - (x2**2)*x1*y3 - (y2**2)*y1*x3& + &- (y2**2)*x1*y3) + case (10) + monom = -(2.d0*x2*y2*x1*x3 - (x2**2)*y1*x3 + (x2**2)*x1*y3& + &+ 2.d0*x2*y2*y1*y3 + (y2**2)*y1*x3 - (y2**2)*x1*y3) + case (11) + monom = -(2.d0*x2*y2*x1*x3 + (x2**2)*y1*x3 - (x2**2)*x1*y3& + &+ 2.d0*x2*y2*y1*y3 - (y2**2)*y1*x3 + (y2**2)*x1*y3) + + case (12) + monom = (2.d0*x3*y3*x2*x1 + (x3**2)*y2*x1 + (x3**2)*x2*y1& + &- 2.d0*x3*y3*y2*y1 - (y3**2)*y2*x1 - (y3**2)*x2*y1) + case (13) + monom = (-(x3**2)*y2*x1 - (x3**2)*x2*y1 - (y3**2)*y2*x1& + &- (y3**2)*x2*y1) + case (14) + monom = -(2.d0*x3*y3*x2*x1 - (x3**2)*y2*x1 + (x3**2)*x2*y1& + &+ 2.d0*x3*y3*y2*y1 + (y3**2)*y2*x1 - (y3**2)*x2*y1) + case (15) + monom = -(2.d0*x3*y3*x2*x1 + (x3**2)*y2*x1 - (x3**2)*x2*y1& + &+ 2.d0*x3*y3*y2*y1 - (y3**2)*y2*x1 + (y3**2)*x2*y1) + case default + error stop 'called case of z_eee_monom not implemented' + monom = 0.0_dp + end select + end function z_eee_monom + + function w_aaee_monom(a, b, x1, y1, x2, y2, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: a, b, x1, y1, x2, y2 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 1) ! order 2,2 => 4 + + case (2) + monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 + case (3) + monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 + + case (4) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 2) ! order 2,3 => 5 + case (5) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 3) ! order 2,3 => 5 + case (6) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 4) ! order 2,3 => 5 + case (7) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 5) ! order 2,3 => 5 + + case (8) + monom = v_aa_monom(a, b, 6)*w_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 + case (9) + monom = v_aa_monom(a, b, 5)*w_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 + case (10) + monom = v_aa_monom(a, b, 4)*w_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 + + case (11) + monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 + case (12) + monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 + case (13) + monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 + case (14) + monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 + case (15) + monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 + case (16) + monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 + case (17) + monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 + case (18) + monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 + + case (19) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 6) ! order 2,4 => 6 + case (20) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 7) ! order 2,4 => 6 + case (21) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 8) ! order 2,4 => 6 + case (22) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 9) ! order 2,4 => 6 + case (23) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 10) ! order 2,4 => 6 + case (24) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 11) ! order 2,4 => 6 + case (25) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 12) ! order 2,4 => 6 + case (26) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 13) ! order 2,4 => 6 + case (27) + monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 14) ! order 2,4 => 6 + + case default + error stop 'called case of w_aaee_monom not implemented' + monom = 0.0_dp + end select + end function w_aaee_monom + + function z_aaee_monom(a, b, x1, y1, x2, y2, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: a, b, x1, y1, x2, y2 + integer(kind=idp), intent(in) :: o + select case (o) + case (1) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 1) ! order 2,2 => 4 + + case (2) + monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 + case (3) + monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 + + case (4) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 2) ! order 2,3 => 5 + case (5) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 3) ! order 2,3 => 5 + case (6) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 4) ! order 2,3 => 5 + case (7) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 5) ! order 2,3 => 5 + + case (8) + monom = v_aa_monom(a, b, 6)*z_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 + case (9) + monom = v_aa_monom(a, b, 5)*z_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 + case (10) + monom = v_aa_monom(a, b, 4)*z_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 + + case (11) + monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 + case (12) + monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 + case (13) + monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 + case (14) + monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 + case (15) + monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 + case (16) + monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 + case (17) + monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 + case (18) + monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 + + case (19) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 6) ! order 2,4 => 6 + case (20) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 7) ! order 2,4 => 6 + case (21) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 8) ! order 2,4 => 6 + case (22) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 9) ! order 2,4 => 6 + case (23) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 10) ! order 2,4 => 6 + case (24) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 11) ! order 2,4 => 6 + case (25) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 12) ! order 2,4 => 6 + case (26) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 13) ! order 2,4 => 6 + case (27) + monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 14) ! order 2,4 => 6 + + case default + error stop 'called case of z_aaee_monom not implemented' + monom = 0.0_dp + end select + end function z_aaee_monom + + function q_u_monom(u, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: u + integer(kind=idp), intent(in) :: o + ! only odd orders of u are allowed + monom = u**(2*(o - 1) + 1) + end function q_u_monom + + function q_ua_monom(u, a, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: u, a + integer(kind=idp), intent(in) :: o + ! only odd orders of u are allowed + select case (o) + case (1) !order 2 + monom = u*a + case (2) !order 3 + monom = u*a**2 + case (3) !order 4 + monom = u**3*a + case (4) + monom = u*a**3 + case (5) ! order 5 + monom = u**3*a**2 + case (6) + monom = u*a**4 + case (7) ! order 6 + monom = u**5*a + case (8) + monom = u**3*a**3 + case (9) + monom = u*a**5 + case (10) !order 7 + monom = u**5*a**2 + case (11) + monom = u**3*a**4 + case (12) + monom = u*a**6 + case (13) ! order 8 + monom = u**7*a + case (14) + monom = u**5*a**3 + case (15) + monom = u**3*a**5 + case (16) + monom = u*a**7 + case default + error stop 'called case of q_ua_monom not implemented' + monom = 0.0_dp + end select + end function q_ua_monom + + function q_ue_monom(u, x, y, o) result(monom) + real(kind=dp) :: monom + real(kind=dp), intent(in) :: u, x, y + integer(kind=idp), intent(in) :: o + select case (o) + case (1) !order 3 + monom = u*v_e_monom(x, y, 1) + case (2) !order 4 + monom = u*v_e_monom(x, y, 2) + case (3) !order 5 + monom = u**3*v_e_monom(x, y, 1) + case (4) + monom = u*v_e_monom(x, y, 3) + case (5) !order 6 + monom = u**3*v_e_monom(x, y, 2) + case (6) + monom = u*v_e_monom(x, y, 4) + case (7) !order 7 + monom = u**5*v_e_monom(x, y, 1) + case (8) + monom = u**3*v_e_monom(x, y, 3) + case (9) + monom = u*v_e_monom(x, y, 5) + case (10) + monom = u*v_e_monom(x, y, 6) + case default + error stop 'called case of q_ue not implemented' + monom = 0.0_dp + end select + end function q_ue_monom + + function q_uae_monom(u, a, x, y, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, a, x, y + integer(idp), intent(in) :: o + select case (o) + case (1) !order 4 + monom = u*v_ae_monom(a, x, y, 1) + case (2) !order 5 + monom = u*v_ae_monom(a, x, y, 2) + case (3) + monom = u*v_ae_monom(a, x, y, 3) + case (4) !order 6 + monom = u**3*v_ae_monom(a, x, y, 1) + case (5) + monom = u*v_ae_monom(a, x, y, 4) + case (6) + monom = u*v_ae_monom(a, x, y, 5) + case (7) + monom = u*v_ae_monom(a, x, y, 6) + case (8) !order 7 + monom = u**3*v_ae_monom(a, x, y, 2) + case (9) + monom = u**3*v_ae_monom(a, x, y, 3) + case (10) + monom = u*v_ae_monom(a, x, y, 7) + case (11) + monom = u*v_ae_monom(a, x, y, 8) + case (12) + monom = u*v_ae_monom(a, x, y, 9) + case (13) + monom = u*v_ae_monom(a, x, y, 10) + case default + error stop 'called case of q_uae not implemented' + monom = 0.0_dp + end select + end function q_uae_monom + + function q_uee_monom(u, x1, y1, x2, y2, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, x1, y1, x2, y2 + integer(idp), intent(in) :: o + select case (o) + case (1) !order 3 + monom = u*v_ee_monom(x1, y1, x2, y2, 1) + case (2) !order 4 + monom = u*v_ee_monom(x1, y1, x2, y2, 2) + case (3) + monom = u*v_ee_monom(x1, y1, x2, y2, 3) + case (4) !order 5 + monom = u**3*v_ee_monom(x1, y1, x2, y2, 1) + case (5) + monom = u*v_ee_monom(x1, y1, x2, y2, 4) + case (6) + monom = u*v_ee_monom(x1, y1, x2, y2, 5) + case (7) + monom = u*v_ee_monom(x1, y1, x2, y2, 6) + case (8) + monom = u*v_ee_monom(x1, y1, x2, y2, 7) + case (9) ! order 6 + monom = u**3*v_ee_monom(x1, y1, x2, y2, 2) + case (10) + monom = u**3*v_ee_monom(x1, y1, x2, y2, 3) + case (11) + monom = u*v_ee_monom(x1, y1, x2, y2, 8) + case (12) + monom = u*v_ee_monom(x1, y1, x2, y2, 9) + case (13) + monom = u*v_ee_monom(x1, y1, x2, y2, 10) + case (14) + monom = u*v_ee_monom(x1, y1, x2, y2, 11) + case (15) + monom = u*v_ee_monom(x1, y1, x2, y2, 12) + case (16) + monom = u*v_ee_monom(x1, y1, x2, y2, 13) + case (17) + monom = u*v_ee_monom(x1, y1, x2, y2, 14) + case (18) + monom = u*v_ee_monom(x1, y1, x2, y2, 15) + case default + error stop 'called case of q_uee not implemented' + monom = 0.0_dp + end select + end function q_uee_monom + + function q_uaee_monom(u, a, x1, y1, x2, y2, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, a, x1, y1, x2, y2 + integer(idp), intent(in) :: o + select case (o) + case (1) !order 4 + monom = u*v_aee_monom(a, x1, y1, x2, y2, 1) + case (2) !order 5 + monom = u*v_aee_monom(a, x1, y1, x2, y2, 2) + case (3) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 3) + case (4) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 4) + case (5) !order 6 + monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 1) + case (6) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 5) + case (7) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 6) + case (8) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 7) + case (9) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 8) + case (10) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 9) + case (11) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 10) + case (12) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 11) + case (13) !order 7 + monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 2) + case (14) + monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 3) + case (15) + monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 4) + case (16) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 12) + case (17) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 13) + case (18) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 14) + case (19) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 15) + case (20) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 16) + case (21) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 17) + case (22) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 18) + case (23) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 19) + case (24) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 20) + case (25) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 21) + case (26) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 22) + case (27) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 23) + case (28) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 24) + case (29) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 25) + case (30) + monom = u*v_aee_monom(a, x1, y1, x2, y2, 26) + case default + error stop 'called case of q_uaee not implemented' + monom = 0.0_dp + end select + end function q_uaee_monom + + function qw_ue_monom(u, x, y, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, x, y + integer(idp), intent(in) :: o + select case (o) + case (1) !order 2 + monom = u*w_e_monom(x, y, 1) + case (2) !order 3 + monom = u*w_e_monom(x, y, 2) + case (3) !order 4 + monom = u**3*w_e_monom(x, y, 1) + case (4) + monom = u*w_e_monom(x, y, 3) + case (5) !order 5 + monom = u**3*w_e_monom(x, y, 2) + case (6) + monom = u*w_e_monom(x, y, 4) + case (7) + monom = u*w_e_monom(x, y, 5) + case (8) !order 6 + monom = u**5*w_e_monom(x, y, 1) + case (9) + monom = u**3*w_e_monom(x, y, 3) + case (10) + monom = u*w_e_monom(x, y, 6) + case (11) + monom = u*w_e_monom(x, y, 7) + case (12) !order 7 + monom = u**5*w_e_monom(x, y, 2) + case (13) + monom = u**3*w_e_monom(x, y, 4) + case (14) + monom = u**3*w_e_monom(x, y, 5) + case (15) + monom = u*w_e_monom(x, y, 8) + case (16) + monom = u*w_e_monom(x, y, 9) + case default + error stop 'called case of qw_ue_monom not implemented' + monom = 0.0_dp + end select + end function qw_ue_monom + + function qz_ue_monom(u, x, y, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, x, y + integer(idp), intent(in) :: o + select case (o) + case (1) !order 2 + monom = u*z_e_monom(x, y, 1) + case (2) !order 3 + monom = u*z_e_monom(x, y, 2) + case (3) !order 4 + monom = u**3*z_e_monom(x, y, 1) + case (4) + monom = u*z_e_monom(x, y, 3) + case (5) !order 5 + monom = u**3*z_e_monom(x, y, 2) + case (6) + monom = u*z_e_monom(x, y, 4) + case (7) + monom = u*z_e_monom(x, y, 5) + case (8) !order 6 + monom = u**5*z_e_monom(x, y, 1) + case (9) + monom = u**3*z_e_monom(x, y, 3) + case (10) + monom = u*z_e_monom(x, y, 6) + case (11) + monom = u*z_e_monom(x, y, 7) + case (12) !order 7 + monom = u**5*z_e_monom(x, y, 2) + case (13) + monom = u**3*z_e_monom(x, y, 4) + case (14) + monom = u**3*z_e_monom(x, y, 5) + case (15) + monom = u*z_e_monom(x, y, 8) + case (16) + monom = u*z_e_monom(x, y, 9) + case default + error stop 'called case of qz_ue_monom not implemented' + monom = 0.0_dp + end select + end function qz_ue_monom + + function qw_uee_monom(u, x1, y1, x2, y2, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, x1, y1, x2, y2 + integer(idp), intent(in) :: o + select case (o) + case (1) !order 3 + monom = u*w_ee_monom(x1, y1, x2, y2, 1) + case (2) !order 4 + monom = u*w_ee_monom(x1, y1, x2, y2, 2) + case (3) + monom = u*w_ee_monom(x1, y1, x2, y2, 3) + case (4) + monom = u*w_ee_monom(x1, y1, x2, y2, 4) + case (5) + monom = u*w_ee_monom(x1, y1, x2, y2, 5) + case (6) !order 5 + monom = u**3*w_ee_monom(x1, y1, x2, y2, 1) + case (7) + monom = u*w_ee_monom(x1, y1, x2, y2, 6) + case (8) + monom = u*w_ee_monom(x1, y1, x2, y2, 7) + case (9) + monom = u*w_ee_monom(x1, y1, x2, y2, 8) + case (10) + monom = u*w_ee_monom(x1, y1, x2, y2, 9) + case (11) + monom = u*w_ee_monom(x1, y1, x2, y2, 10) + case (12) + monom = u*w_ee_monom(x1, y1, x2, y2, 11) + case (13) + monom = u*w_ee_monom(x1, y1, x2, y2, 12) + case (14) + monom = u*w_ee_monom(x1, y1, x2, y2, 13) + case (15) + monom = u*w_ee_monom(x1, y1, x2, y2, 14) + case default + error stop 'called case of qw_uee_monom not implemented' + monom = 0.0_dp + end select + end function qw_uee_monom + + function qz_uee_monom(u, x1, y1, x2, y2, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, x1, y1, x2, y2 + integer(idp), intent(in) :: o + select case (o) + case (1) !order 3 + monom = u*z_ee_monom(x1, y1, x2, y2, 1) + case (2) !order 4 + monom = u*z_ee_monom(x1, y1, x2, y2, 2) + case (3) + monom = u*z_ee_monom(x1, y1, x2, y2, 3) + case (4) + monom = u*z_ee_monom(x1, y1, x2, y2, 4) + case (5) + monom = u*z_ee_monom(x1, y1, x2, y2, 5) + case (6) !order 5 + monom = u**3*z_ee_monom(x1, y1, x2, y2, 1) + case (7) + monom = u*z_ee_monom(x1, y1, x2, y2, 6) + case (8) + monom = u*z_ee_monom(x1, y1, x2, y2, 7) + case (9) + monom = u*z_ee_monom(x1, y1, x2, y2, 8) + case (10) + monom = u*z_ee_monom(x1, y1, x2, y2, 9) + case (11) + monom = u*z_ee_monom(x1, y1, x2, y2, 10) + case (12) + monom = u*z_ee_monom(x1, y1, x2, y2, 11) + case (13) + monom = u*z_ee_monom(x1, y1, x2, y2, 12) + case (14) + monom = u*z_ee_monom(x1, y1, x2, y2, 13) + case (15) + monom = u*z_ee_monom(x1, y1, x2, y2, 14) + case default + error stop 'called case of qz_uee_monom not implemented' + monom = 0.0_dp + end select + end function qz_uee_monom + + function qw_uaee_monom(u, a, x1, y1, x2, y2, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, a, x1, y1, x2, y2 + integer(idp), intent(in) :: o + select case (o) + case (1) !order 4 + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 1) + case (2) !order 5 + monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 1) + case (3) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 2) + case (4) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 3) + case (5) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 4) + case (6) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 5) + case (7) !order 6 + monom = q_ua_monom(u, a, 3)*w_ee_monom(x1, y1, x2, y2, 1) + case (8) + monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 2) + case (9) + monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 3) + case (10) + monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 4) + case (11) + monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 5) + case (12) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 6) + case (13) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 7) + case (14) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 8) + case (15) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 9) + case (16) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 10) + case (17) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 11) + case (18) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 12) + case (19) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 13) + case (20) + monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 14) + case default + error stop 'called case of qw_uaee_monom not implemented' + monom = 0.0_dp + end select + end function qw_uaee_monom + + function qz_uaee_monom(u, a, x1, y1, x2, y2, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, a, x1, y1, x2, y2 + integer(idp), intent(in) :: o + select case (o) + case (1) !order 4 + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 1) + case (2) !order 5 + monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 1) + case (3) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 2) + case (4) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 3) + case (5) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 4) + case (6) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 5) + case (7) !order 6 + monom = q_ua_monom(u, a, 3)*z_ee_monom(x1, y1, x2, y2, 1) + case (8) + monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 2) + case (9) + monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 3) + case (10) + monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 4) + case (11) + monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 5) + case (12) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 6) + case (13) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 7) + case (14) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 8) + case (15) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 9) + case (16) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 10) + case (17) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 11) + case (18) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 12) + case (19) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 13) + case (20) + monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 14) + case default + error stop 'called case of qz_uaee_monom not implemented' + monom = 0.0_dp + end select + end function qz_uaee_monom + + function qw_uae_monom(u, a, x, y, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, a, x, y + integer(idp), intent(in) :: o + select case (o) + case (1) ! order 3 + monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 1) + case (2) ! order 4 + monom = q_ua_monom(u, a, 2)*w_e_monom(x, y, 1) + case (3) + monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 2) + case (4) !order 5 + monom = q_ua_monom(u, a, 3)*w_e_monom(x, y, 1) + case (5) + monom = q_ua_monom(u, a, 4)*w_e_monom(x, y, 1) + case (6) + monom = q_ua_monom(u, a, 2)*w_e_monom(x, y, 2) + case (7) + monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 3) + case (8) !order 6 + monom = q_ua_monom(u, a, 5)*w_e_monom(x, y, 1) + case (9) + monom = q_ua_monom(u, a, 6)*w_e_monom(x, y, 1) + case (10) + monom = q_ua_monom(u, a, 3)*w_e_monom(x, y, 2) + case (11) + monom = q_ua_monom(u, a, 4)*w_e_monom(x, y, 2) + case (12) + monom = q_ua_monom(u, a, 2)*w_e_monom(x, y, 3) + case (13) + monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 4) + case (14) + monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 5) + case default + error stop 'called case of qw_uae not implemented' + monom = 0.0_dp + end select + end function qw_uae_monom + + function qz_uae_monom(u, a, x, y, o) result(monom) + real(dp) :: monom + real(dp), intent(in) :: u, a, x, y + integer(idp), intent(in) :: o + select case (o) + case (1) ! order 3 + monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 1) + case (2) ! order 4 + monom = q_ua_monom(u, a, 2)*z_e_monom(x, y, 1) + case (3) + monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 2) + case (4) !order 5 + monom = q_ua_monom(u, a, 3)*z_e_monom(x, y, 1) + case (5) + monom = q_ua_monom(u, a, 4)*z_e_monom(x, y, 1) + case (6) + monom = q_ua_monom(u, a, 2)*z_e_monom(x, y, 2) + case (7) + monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 3) + case (8) !order 6 + monom = q_ua_monom(u, a, 5)*z_e_monom(x, y, 1) + case (9) + monom = q_ua_monom(u, a, 6)*z_e_monom(x, y, 1) + case (10) + monom = q_ua_monom(u, a, 3)*z_e_monom(x, y, 2) + case (11) + monom = q_ua_monom(u, a, 4)*z_e_monom(x, y, 2) + case (12) + monom = q_ua_monom(u, a, 2)*z_e_monom(x, y, 3) + case (13) + monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 4) + case (14) + monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 5) + case default + error stop 'called case of qz_uae not implemented' + monom = 0.0_dp + end select + end function qz_uae_monom + +end module select_monom_mod diff --git a/src/model/new_PES/sphericalharmonics_mod.f90 b/src/model/new_PES/sphericalharmonics_mod.f90 new file mode 100644 index 0000000..47af619 --- /dev/null +++ b/src/model/new_PES/sphericalharmonics_mod.f90 @@ -0,0 +1,575 @@ +! 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 diff --git a/src/model/new_PES/surface_mod.f90 b/src/model/new_PES/surface_mod.f90 new file mode 100644 index 0000000..885c02c --- /dev/null +++ b/src/model/new_PES/surface_mod.f90 @@ -0,0 +1,71 @@ +module surface_mod + use accuracy_constants, only: dp + implicit none + private + public eval_surface + contains + subroutine eval_surface(e, w, u, x1, p) + use ctrans_pes_mod, only: ctrans_pes + use diab_pes, only: pote + use accuracy_constants, only: dp, idp + use dim_parameter, only: ndiab + implicit none + real(dp), dimension(:, :), intent(out) :: w, u + real(dp), dimension(:), intent(out) :: e + real(dp), dimension(:), intent(in) :: x1, p + real(dp), dimension(size(x1, 1)) :: s, t + real(dp), allocatable, dimension(:, :) :: Mat + + !coordinate transformation if needed + call ctrans_pes(x1, s, t) + write(11,'(*(f18.8))') s + block + ! lapack variables + integer(kind=idp), parameter :: lwork = 1000 + real(kind=dp) work(lwork) + integer(kind=idp) info + !evaluate model + call pote(w, 1, x1, s, t, p, size(p, 1)) + allocate (Mat, source=w) + call dsyev('V', 'U', ndiab, Mat, ndiab, e, work, lwork, info) + u(:, :) = Mat(:, :) + deallocate (Mat) + end block + + end subroutine eval_surface + +! !subroutine init_surface(p) +! use dim_parameter, only: ndiab, nstat, ntot, nci ,qn +! use parameterkeys, only: parameterkey_read +! use fileread_mod, only: get_datfile, internalize_datfile +! use io_parameters, only: llen +! use accuracy_constants, only: dp +! implicit none +! real(dp), dimension(:), allocatable, intent(out) :: p +! character(len=llen), allocatable, dimension(:) :: infile +! +! qn = 9 +! ndiab = 4 +! nstat = 4 +! nci = 4 +! ntot = ndiab + nstat + nci +! +! block +! character(len=:),allocatable :: datnam +! integer :: linenum +! !get parameter file +! call get_datfile(datnam) +! !internalize datfile +! call internalize_datfile(datnam, infile, linenum, llen) +! end block +! +! !read parameters from file +! block +! real(dp), dimension(:), allocatable :: p_spread +! integer,dimension(:),allocatable :: p_act +! integer :: npar +! real(dp), parameter :: facspread = 1.0_dp, gspread = 1.0_dp +! call parameterkey_read(infile, size(infile, 1), p, p_act, p_spread, npar, gspread, facspread) +! end block +! end subroutine init_surface + end module surface_mod diff --git a/src/model/new_PES/xy_stretch_lib.f90 b/src/model/new_PES/xy_stretch_lib.f90 new file mode 100644 index 0000000..9b30191 --- /dev/null +++ b/src/model/new_PES/xy_stretch_lib.f90 @@ -0,0 +1,82 @@ +module xy_stretch_lib + use accuracy_constants, only: dp,idp + implicit none + private + public eval_surface, init_surface,eval_matrix + real(dp), dimension(:), allocatable :: p + contains +subroutine eval_surface(e, w, u, x1) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + use dim_parameter, only: ndiab + implicit none + real(dp), dimension(:, :), intent(out) :: w, u + real(dp), dimension(:), intent(out) :: e + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + real(dp), allocatable, dimension(:, :) :: Mat + + !coordinate transformation if needed + call ctrans(x1, s, t) + + block + ! lapack variables + integer(kind=idp), parameter :: lwork = 1000 + real(kind=dp) work(lwork) + integer(kind=idp) info + !evaluate model + call diab(w, 1, x1, s, t, p, size(p, 1)) + allocate (Mat, source=w) + call dsyev('V', 'U', ndiab, Mat, ndiab, e, work, lwork, info) + u(:, :) = Mat(:, :) + deallocate (Mat) + end block + +end subroutine eval_surface + +subroutine eval_matrix(w,x1) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + implicit none + real(dp), dimension(:, :), intent(out) :: w + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + + !coordinate transformation if needed + call ctrans(x1, s, t) + call diab(w, 1, x1, s, t, p, size(p, 1)) +end subroutine eval_matrix + +subroutine init_surface() + use dim_parameter, only: ndiab, nstat, ntot, nci ,qn + use parameterkeys, only: parameterkey_read + use fileread_mod, only: get_datfile, internalize_datfile + use io_parameters, only: llen + use accuracy_constants, only: dp + implicit none + character(len=llen), allocatable, dimension(:) :: infile + + qn = 9 + ndiab = 4 + nstat = 4 + nci = 4 + ntot = ndiab + nstat + nci + + block + character(len=:),allocatable :: datnam + integer :: linenum + datnam = 'xy_stretch.par.save' + ! datnam = 'umbstr.par.save' + call internalize_datfile(datnam, infile, linenum, llen) + end block + + !read parameters from file + block + real(dp), dimension(:), allocatable :: p_spread + integer,dimension(:),allocatable :: p_act + integer :: npar + real(dp), parameter :: facspread = 1.0_dp, gspread = 1.0_dp + call parameterkey_read(infile, size(infile, 1), p, p_act, p_spread, npar, gspread, facspread) + end block +end subroutine init_surface +end module xy_stretch_lib diff --git a/src/model/new_PES/xy_stretch_lib.f90~ b/src/model/new_PES/xy_stretch_lib.f90~ new file mode 100644 index 0000000..db8cc78 --- /dev/null +++ b/src/model/new_PES/xy_stretch_lib.f90~ @@ -0,0 +1,82 @@ +module xy_stretch_lib + use accuracy_constants, only: dp,idp + implicit none + private + public eval_surface, init_surface,eval_matrix + real(dp), dimension(:), allocatable, intent(out) :: p + contains +subroutine eval_surface(e, w, u, x1, p) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + use dim_parameter, only: ndiab + implicit none + real(dp), dimension(:, :), intent(out) :: w, u + real(dp), dimension(:), intent(out) :: e + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + real(dp), allocatable, dimension(:, :) :: Mat + + !coordinate transformation if needed + call ctrans(x1, s, t) + + block + ! lapack variables + integer(kind=idp), parameter :: lwork = 1000 + real(kind=dp) work(lwork) + integer(kind=idp) info + !evaluate model + call diab(w, 1, x1, s, t, p, size(p, 1)) + allocate (Mat, source=w) + call dsyev('V', 'U', ndiab, Mat, ndiab, e, work, lwork, info) + u(:, :) = Mat(:, :) + deallocate (Mat) + end block + +end subroutine eval_surface + +subroutine eval_matrix(w,x1, p) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + implicit none + real(dp), dimension(:, :), intent(out) :: w + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + + !coordinate transformation if needed + call ctrans(x1, s, t) + call diab(w, 1, x1, s, t, p, size(p, 1)) +end subroutine eval_matrix + +subroutine init_surface(p) + use dim_parameter, only: ndiab, nstat, ntot, nci ,qn + use parameterkeys, only: parameterkey_read + use fileread_mod, only: get_datfile, internalize_datfile + use io_parameters, only: llen + use accuracy_constants, only: dp + implicit none + character(len=llen), allocatable, dimension(:) :: infile + + qn = 9 + ndiab = 4 + nstat = 4 + nci = 4 + ntot = ndiab + nstat + nci + + block + character(len=:),allocatable :: datnam + integer :: linenum + datnam = 'xy_stretch.par.save' + ! datnam = 'umbstr.par.save' + call internalize_datfile(datnam, infile, linenum, llen) + end block + + !read parameters from file + block + real(dp), dimension(:), allocatable :: p_spread + integer,dimension(:),allocatable :: p_act + integer :: npar + real(dp), parameter :: facspread = 1.0_dp, gspread = 1.0_dp + call parameterkey_read(infile, size(infile, 1), p, p_act, p_spread, npar, gspread, facspread) + end block +end subroutine init_surface +end module xy_stretch_lib diff --git a/src/model/new_PES/xyumb_stretch_lib.f90 b/src/model/new_PES/xyumb_stretch_lib.f90 new file mode 100644 index 0000000..cc8073e --- /dev/null +++ b/src/model/new_PES/xyumb_stretch_lib.f90 @@ -0,0 +1,83 @@ + +module xyumb_stretch_lib + use accuracy_constants, only: dp,idp + implicit none + private + public eval_surface, init_surface,eval_matrix + real(dp), dimension(:), allocatable :: p + contains +subroutine eval_surface(e, w, u, x1) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + use dim_parameter, only: ndiab + implicit none + real(dp), dimension(:, :), intent(out) :: w, u + real(dp), dimension(:), intent(out) :: e + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + real(dp), allocatable, dimension(:, :) :: Mat + + !coordinate transformation if needed + call ctrans(x1, s, t) + + block + ! lapack variables + integer(kind=idp), parameter :: lwork = 1000 + real(kind=dp) work(lwork) + integer(kind=idp) info + !evaluate model + call diab(w, 1, x1, s, t, p, size(p, 1)) + allocate (Mat, source=w) + call dsyev('V', 'U', ndiab, Mat, ndiab, e, work, lwork, info) + u(:, :) = Mat(:, :) + deallocate (Mat) + end block + +end subroutine eval_surface + +subroutine eval_matrix(w,x1) + use ctrans_mod, only: ctrans + use diabmodel, only: diab + implicit none + real(dp), dimension(:, :), intent(out) :: w + real(dp), dimension(:), intent(in) :: x1 + real(dp), dimension(size(x1, 1)) :: s, t + + !coordinate transformation if needed + call ctrans(x1, s, t) + call diab(w, 1, x1, s, t, p, size(p, 1)) +end subroutine eval_matrix + +subroutine init_surface() + use dim_parameter, only: ndiab, nstat, ntot, nci ,qn + use parameterkeys, only: parameterkey_read + use fileread_mod, only: get_datfile, internalize_datfile + use io_parameters, only: llen + use accuracy_constants, only: dp + implicit none + character(len=llen), allocatable, dimension(:) :: infile + + qn = 9 + ndiab = 4 + nstat = 4 + nci = 4 + ntot = ndiab + nstat + nci + + block + character(len=:),allocatable :: datnam + integer :: linenum + !datnam = 'xy_stretch.par.save' + datnam = 'umbstr.par.save' + call internalize_datfile(datnam, infile, linenum, llen) + end block + + !read parameters from file + block + real(dp), dimension(:), allocatable :: p_spread + integer,dimension(:),allocatable :: p_act + integer :: npar + real(dp), parameter :: facspread = 1.0_dp, gspread = 1.0_dp + call parameterkey_read(infile, size(infile, 1), p, p_act, p_spread, npar, gspread, facspread) + end block +end subroutine init_surface +end module xyumb_stretch_lib diff --git a/src/model/nnparams.incl b/src/model/nnparams.incl new file mode 100644 index 0000000..635c33c --- /dev/null +++ b/src/model/nnparams.incl @@ -0,0 +1,43 @@ +!**** 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) diff --git a/src/model/old_keys.f90 b/src/model/old_keys.f90 new file mode 100644 index 0000000..3d5e2a3 --- /dev/null +++ b/src/model/old_keys.f90 @@ -0,0 +1,104 @@ +module keys_mod + implicit none +contains + subroutine init_keys + use io_parameters, only: key + character(len=1) prefix(4) + parameter (prefix=['N','P','A','S']) + !character (len=20) key(4,56) + + character(len=16) parname(56) + integer i,j + + ! the electronic states are ordered as: A2" E' and A1' + ! the name convention here is : A2 E1 A1 + + ! Naming + !-------------------- + ! V: V-TERM OR diagonal term + ! J: Jahn teller coupling term in E + ! P: pseudo jahn teller between As and E + + ! S: it involves the symmetric term of x**2+y**2 + ! N: It does not involve symmetric term + + + ! diagonal term for 4 states + !parname( 1)='VA2N0' ! order 0 + !parname( 2)='VA2S0' ! order 0 with S + !parname( 3)='VE1N0' ! order 0 witH N + !parname( 4)='VE1S0' ! order 0 with S + !parname( 5)='VA1N0' ! order 0 with N + !parname( 6)='VA1S0' ! order 0 with S + parname( 7)='VA2N1' ! order 1 + parname( 8)='VA2S1' ! order 1 with S + parname( 9)='VE1N1' ! order 1 witH N + parname(10)='VE1S1' ! order 1 with S + parname(11)='VA1N1' ! order 1 with N + parname(12)='VA1S1' ! order 1 with S + parname(13)='VA2N2' ! order 2 + parname(14)='VA2S2' ! order 2 with S ! only 2 term + parname(15)='VE1N2' ! order 2 witH N + parname(16)='VE1S2' ! order 2 with S + parname(17)='VA1N2' ! order 2 with N + parname(18)='VA1S2' ! order 2 with S + !parname(19)='VA2N3' ! order 3 + !parname(20)='VA2S3' ! order 3 with S + !parname(21)='VE1N3' ! order 3 witH N + !parname(22)='VE1S3' ! order 3 with S + !parname(23)='VA1N3' ! order 3 with N + !parname(24)='VA1S3' ! order 3 with S + + ! Jahn teller within E + + parname(25)='JE1N0' ! order 0 with N + parname(26)='JE1S0' ! order 0 with S + parname(27)='JE1N1' ! order 1 with N + parname(28)='JE1S1' ! order 1 with S + parname(29)='JE1N2' ! order 2 + parname(30)='JE1S2' ! order 2 + parname(31)='JE1N3' ! order 3 ! this has 8 terms + !parname(32)='JE1S3' ! order 3 ! i do not have this term + + ! PSeudo Jahn teller couplings + + ! coupling of A2 with other + !parname(33)='PA2E1N0' ! order 0 ! is not there + ! parname(34)='PA2E1S0' ! order 0 ! is not there + ! parname(35)='PA2A1N0' ! ORDER 0 + ! parname(36)='PA2A1S0' ! order 0 + !parname(37)='PA2E1N1' ! order 1 + !parname(38)='PA2E1S1' ! order 1 + parname(39)='PA2A1N1' ! order 1 + parname(40)='PA2A1S1' ! order 1 + !parname(41)='PA2E1N2' ! order 2 + !parname(42)='PA2E1S2' ! order 2 + parname(43)='PA2A1N2' + parname(44)='PA2A1S2' + parname(45)='PA2E1N3' ! order 3 + !parname(46)='PA2E1S3' ! order 3 + parname(47)='PA2A1N3' + !parname(48)='PA2A1S3' + + + + ! coupling of A1 with other + ! A2 with A1 is already included above + + parname(49)='PE1A1N0' ! order 1 + parname(50)='PE1A1S0' + !parname(51)='PE1A1N1' + !parname(52)='PE1A1S1' + !parname(53)='PE1A1N2' + !parname(54)='PE1A1S2' + parname(55)='PE1A1N3' + !parname(56)='PE1A1S3' + + do i=1,56 + do j=1,4 + key(i, j)=prefix(j)//trim(parname(i))//':' + enddo + enddo + end subroutine + +end module keys_mod diff --git a/src/model/pyramidal_model.f90 b/src/model/pyramidal_model.f90 new file mode 100644 index 0000000..206a487 --- /dev/null +++ b/src/model/pyramidal_model.f90 @@ -0,0 +1,357 @@ +module diabmodel + use dim_parameter,only:qn,ndiab,pst + use accuracy_constants, only:dp,idp + implicit none + logical :: debug=.false. + contains + + subroutine diab(ex,ey,n,x1,x2,p) + use ctrans_mod, only:ctrans + integer,intent(in),optional :: n ! number of parameter & nmbr of points \ + integer id + integer key,i,j + double precision, intent(in)::x1(qn),x2(qn) + double precision, contiguous,intent(in):: p(:)! array containing parameters + double precision, intent(out)::ex(ndiab,ndiab),ey(ndiab,ndiab) + key =87 + call diab_x(ex,x1,x2,key,p) + !ey=0.0d0 + call diab_y(ey,x1,x2,key,p) + end subroutine + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine diab_x(e,q,t,key,p) + + + real(dp),intent(in)::q(qn),t(qn) + real(dp),intent(out)::e(:,:) + integer(idp),intent(in)::key + real(dp),intent(in),contiguous::p(:) + integer(idp) id,i,j + real(dp) tmp_v,xs,xb,ys,yb,a,b,ss,sb,v3_vec(8) + xs=q(2) + ys=q(3) + xb=q(4) + yb=q(5) + a=q(1) + b=q(6) + + ss=xs**2+ys**2 ! totaly symmetric term + sb=xb**2+yb**2 + + v3_vec( 1) = xs*(xs**2-3*ys**2) + v3_vec( 2) = xb*(xb**2-3*yb**2) + v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys + v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb + v3_vec( 5) = ys*(3*xs**2-ys**2) + v3_vec( 6) = yb*(3*xb**2-yb**2) + v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys + v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb + + + + e=0.0d0 + + + + id=key !1 + ! V-term + ! order 1 + e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + id=id+1 !2 + e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + e(3,3)=e(3,3)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + id=id+1 !3 + e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + ! order 2 + id=id+1 !4 + e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& + +p(pst(1,id)+2)*(xs*xb-ys*yb) + 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) +& + p(pst(1,id)+2)*(xs*xb-ys*yb) + e(3,3)=e(3,3)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) +& + p(pst(1,id)+2)*(xs*xb-ys*yb) + id=id+1 !6 + e(4,4)=e(4,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) + & + p(pst(1,id)+2)*(xs*xb-ys*yb) + ! order 3 + id=id+1 !7 + e(1,1)=e(1,1)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb + b**2* & + (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb) + id=id+1 !8 + e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb+ b**2* & + (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb) + e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb+ b**2* & + (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb) + id=id+1 !9 + e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb + b**2* & + (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb) + + ! JAHN TELLER COUPLING W AND Z + ! order 0 + id=id+1 !10 + e(2,2)=e(2,2)+p(pst(1,id)) + e(3,3)=e(3,3)-p(pst(1,id)) + ! order 1 + id=id+1 !11 + e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + e(3,3)=e(3,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb + e(2,3)=e(2,3)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb + ! order 2 + id=id+1 !12 + e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb + & + b**2*(p(pst(1,id)+5) +b**4*(p(pst(1,id)+6)) + b**6*(p(pst(1,id)+7)) + & + b**8*(p(pst(1,id)+8))) + e(3,3)=e(3,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb) -& + b**8*(p(pst(1,id)+5)) + e(2,3)=e(2,3)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+ & + p(pst(1,id)+2)*(xs*yb+xb*ys) + ! order 3 + + id=id+1 !13 + do i=1,4 + j=i-1 + e(2,2)=e(2,2)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) + e(3,3)=e(3,3)-(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) + e(2,3)=e(2,3)+(-p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i+4) + enddo + + e(2,2)=e(2,2)+p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb + e(3,3)=e(3,3)-(p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb) + e(2,3)=e(2,3)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb + ! PSEUDO JAHN TELLER + + ! A2 ground state coupled with E + ! ################################################### + ! ################################################### + + ! order 0 + id=id+1 !14 + e(1,2)=e(1,2)+b*p(pst(1,id)) + + ! order 1 + id=id+1 !15 + e(1,2)=e(1,2)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) + e(1,3)=e(1,3)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) + ! order 2 + id=id+1 !16 + e(1,2)=e(1,2)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& + +p(pst(1,id)+2)*(xs*xb-ys*yb)) + e(1,3)=e(1,3)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& + +p(pst(1,id)+2)*(xs*yb+xb*ys)) + + !! THE COUPLING OF A2 WITH A1 + !#################################################### + !#################################################### + ! order 1 + id=id+1 !17 + e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) + id=id+1 !18 + e(1,4)=e(1,4)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& + +p(pst(1,id)+2)*(xs*xb-ys*yb)) + + + !!! THE COUPLING OF A1 WITH E + !!#################################################### + !#################################################### + ! order 0 + id=id+1 !19 + e(2,4)=e(2,4)+p(pst(1,id)) + + ! order 1 + id=id+1 !20 + e(2,4)=e(2,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + e(3,4)=e(3,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb + + ! order 2 + id=id+1 !21 + e(2,4)=e(2,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb) + e(3,4)=e(3,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & + -p(pst(1,id)+2)*(xs*yb+xb*ys) + + !! End of the model + + e(2,1)=e(1,2) + e(3,1)=e(1,3) + e(3,2)=e(2,3) + e(4,1)=e(1,4) + e(4,2)=e(2,4) + e(4,3)=e(3,4) + end subroutine diab_x + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! THE Y COMPONENT OF DIPOLE + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine diab_y(e,q,t,key,p) + !integer(idp), intent(in)::npar + real(dp),intent(in)::q(qn),t(qn) + real(dp),intent(out)::e(:,:) + integer(idp),intent(in):: key + real(dp),intent(in),contiguous::p(:) + integer(idp) id,i,j + real(dp) tmp_v,ys,xb,a,b,xs,yb,ss,sb,v3_vec(8) + xs=q(2) + ys=q(3) + xb=q(4) + yb=q(5) + a=q(1) + b=q(6) + + ss=xs**2+ys**2 ! totaly symmetric term + sb=xb**2+yb**2 + + v3_vec( 1) = xs*(xs**2-3*ys**2) + v3_vec( 2) = xb*(xb**2-3*yb**2) + v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys + v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb + v3_vec( 5) = ys*(3*xs**2-ys**2) + v3_vec( 6) = yb*(3*xb**2-yb**2) + v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys + v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb + + + + + + e=0.0d0 + ! V-term + id=key !1 + e(1,1)=e(1,1)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb + + id=id+1 !2 + 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 + id=id+1 !3 + e(4,4)=e(4,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb + ! order 2 + id=id+1 !4 + e(1,1)=e(1,1)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & + -p(pst(1,id)+2)*(xs*yb+xb*ys) + id=id+1 !5 + e(2,2)=e(2,2)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & + -p(pst(1,id)+2)*(xs*yb+xb*ys) + + + e(3,3)=e(3,3)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & + -p(pst(1,id)+2)*(xs*yb+xb*ys) + id=id+1 !6 + e(4,4)=e(4,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & + -p(pst(1,id)+2)*(xs*yb+xb*ys) + ! order 3 + id=id+1 !7 + e(1,1)=e(1,1)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* & + (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb) + id=id+1 !8 + e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb+b**2* & + (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb) + e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* & + (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb) + id=id+1 !9 + e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* & + (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb) + + ! V- term + totally symmetric coord a + + ! JAHN TELLER COUPLING TERM + ! order 0 + id=id+1 !10 + e(2,3)=e(2,3)+p(pst(1,id)) + ! order 1 + + id=id+1 !11 + 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(2,3)=e(2,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb + !id=id+1 !12 + ! order 2 + id=id+1 !12 + e(2,2)=e(2,2)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+p(pst(1,id)+2)*(xs*yb+xb*ys) + e(3,3)=e(3,3)-p(pst(1,id))*2*xs*ys-p(pst(1,id)+1)*2*xb*yb-p(pst(1,id)+2)*(xs*yb+xb*ys) + e(2,3)=e(2,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)) & + -p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb +& + b**8*(p(pst(1,id)+5)) + ! order 3 + id=id+1 !13 + do i=1,4 + j=i-1 + e(2,2)=e(2,2)+(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) + e(3,3)=e(3,3)-(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) + e(2,3)=e(2,3)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) + enddo + e(2,2)=e(2,2)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb + e(3,3)=e(3,3)+p(pst(1,id)+8)*ys*ss+p(pst(1,id)+9)*yb*sb + e(2,3)=e(2,3)-p(pst(1,id)+8)*xs*ss-p(pst(1,id)+1)*xb*sb + ! PSEUDO JAHN TELLER + ! ORDER 0 + ! THE COUPLING OF A2 GROUND STATE WITH E + ! ################################################### + ! ################################################### + ! order 0 + id=id+1 !14 + e(1,3)=e(1,3)-b*(p(pst(1,id))) + ! order 1 + id=id+1 !15 + e(1,2)=e(1,2)-b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) + e(1,3)=e(1,3)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) + + ! order 2 + id=id+1 !16 + e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& + +p(pst(1,id)+2)*(xs*yb+xb*ys)) + e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& + +p(pst(1,id)+2)*(xs*xb-ys*yb)) + ! THE COUPLING OF A2 WITH A1 + !#################################################### + !#################################################### + ! order 1 + id=id+1 !17 + e(1,4)=e(1,4)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) + ! order 2 + id=id+1 !18 + e(1,4)=e(1,4)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& + +p(pst(1,id)+2)*(xs*yb+xb*ys)) + + + ! THE COUPLING OF A1 WITH E + !#################################################### + !#################################################### + ! order 0 + id=id+1 !19 + e(3,4)=e(3,4)-p(pst(1,id)) + ! order 1 + id=id+1 !20 + e(2,4)=e(2,4)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb + e(3,4)=e(3,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + ! order 2 + id=id+1 !21 + e(2,4)=e(2,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb) & + +p(pst(1,id)+2)*(xs*yb+xb*ys) + e(3,4)=e(3,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb) + ! end of the model + e(2,1)=e(1,2) + e(3,1)=e(1,3) + e(3,2)=e(2,3) + e(4,1)=e(1,4) + e(4,2)=e(2,4) + e(4,3)=e(3,4) +end subroutine diab_y + subroutine copy_2_lower_triangle(mat) + real(dp), intent(inout) :: mat(:, :) + integer :: m, n +! write lower triangle of matrix symmetrical + do n = 2, size(mat, 1) + do m = 1, n - 1 + mat(n, m) = mat(m, n) + end do + end do + end subroutine copy_2_lower_triangle + +end module diabmodel diff --git a/src/model/temp b/src/model/temp new file mode 100644 index 0000000..b76671c --- /dev/null +++ b/src/model/temp @@ -0,0 +1,110 @@ + NVA2N1: + PVA2N1: + AVA2N1: + SVA2N1: + + NVE1N1: + PVE1N1: + AVE1N1: + SVE1N1: + + NVA1N1: + PVA1N1: + AVA1N1: + SVA1N1: + + NVA2N2: + PVA2N2: + AVA2N2: + SVA2N2: + + NVE1N2: + PVE1N2: + AVE1N2: + SVE1N2: + + NVA1N2: + PVA1N2: + AVA1N2: + SVA1N2: + + NVA2N3: + PVA2N3: + AVA2N3: + SVA2N3: + + NVE1N3: + PVE1N3: + AVE1N3: + SVE1N3: + + NVA1N3: + PVA1N3: + AVA1N3: + SVA1N3: + + NJE1N0: + PJE1N0: + AJE1N0: + SJE1N0: + + NJE1N1: + PJE1N1: + AJE1N1: + SJE1N1: + + NJE1N2: + PJE1N2: + AJE1N2: + SJE1N2: + + NJE1N3: + PJE1N3: + AJE1N3: + SJE1N3: + + NPA2E1N0: + PPA2E1N0: + APA2E1N0: + SPA2E1N0: + + NPA2E1N1: + PPA2E1N1: + APA2E1N1: + SPA2E1N1: + + NPA2E1N2: + PPA2E1N2: + APA2E1N2: + SPA2E1N2: + + NPA2A1N1: + PPA2A1N1: + APA2A1N1: + SPA2A1N1: + + NPA2A1N2: + PPA2A1N2: + APA2A1N2: + SPA2A1N2: + + NPE1A1N0: + PPE1A1N0: + APE1A1N0: + SPE1A1N0: + + NPE1A1N1: + PPE1A1N1: + APE1A1N1: + SPE1A1N1: + + NPE1A1N2: + PPE1A1N2: + APE1A1N2: + SPE1A1N2: + + NTYPE_CAL: + PTYPE_CAL: + ATYPE_CAL: + STYPE_CAL: + diff --git a/src/model/temp.keys b/src/model/temp.keys new file mode 100644 index 0000000..4a02c6a --- /dev/null +++ b/src/model/temp.keys @@ -0,0 +1,540 @@ + NEXITEN: + PEXITEN: + AEXITEN: + SEXITEN: + + NTMC_CH: + PTMC_CH: + ATMC_CH: + STMC_CH: + + NEVA1: + PEVA1: + AEVA1: + SEVA1: + + NEVU: + PEVU: + AEVU: + SEVU: + + NEVE1: + PEVE1: + AEVE1: + SEVE1: + + NEVE2: + PEVE2: + AEVE2: + SEVE2: + + NEVA1U: + PEVA1U: + AEVA1U: + SEVA1U: + + NEVA1E1: + PEVA1E1: + AEVA1E1: + SEVA1E1: + + NEVA1E2: + PEVA1E2: + AEVA1E2: + SEVA1E2: + + NEVUE1: + PEVUE1: + AEVUE1: + SEVUE1: + + NEVUE2: + PEVUE2: + AEVUE2: + SEVUE2: + + NEVE1E2: + PEVE1E2: + AEVE1E2: + SEVE1E2: + + NEVA1UE1: + PEVA1UE1: + AEVA1UE1: + SEVA1UE1: + + NEVA1UE2: + PEVA1UE2: + AEVA1UE2: + SEVA1UE2: + + NEVA1E1E2: + PEVA1E1E2: + AEVA1E1E2: + SEVA1E1E2: + + NEVUE1E2: + PEVUE1E2: + AEVUE1E2: + SEVUE1E2: + + NEVA1UE1E2: + PEVA1UE1E2: + AEVA1UE1E2: + SEVA1UE1E2: + + NA2VA1: + PA2VA1: + AA2VA1: + SA2VA1: + + NA2VU: + PA2VU: + AA2VU: + SA2VU: + + NA2VE1: + PA2VE1: + AA2VE1: + SA2VE1: + + NA2VE2: + PA2VE2: + AA2VE2: + SA2VE2: + + NA2VA1U: + PA2VA1U: + AA2VA1U: + SA2VA1U: + + NA2VA1E1: + PA2VA1E1: + AA2VA1E1: + SA2VA1E1: + + NA2VA1E2: + PA2VA1E2: + AA2VA1E2: + SA2VA1E2: + + NA2VUE1: + PA2VUE1: + AA2VUE1: + SA2VUE1: + + NA2VUE2: + PA2VUE2: + AA2VUE2: + SA2VUE2: + + NA2VE1E2: + PA2VE1E2: + AA2VE1E2: + SA2VE1E2: + + NA2VA1UE1: + PA2VA1UE1: + AA2VA1UE1: + SA2VA1UE1: + + NA2VA1UE2: + PA2VA1UE2: + AA2VA1UE2: + SA2VA1UE2: + + NA2VA1E1E2: + PA2VA1E1E2: + AA2VA1E1E2: + SA2VA1E1E2: + + NA2VUE1E2: + PA2VUE1E2: + AA2VUE1E2: + SA2VUE1E2: + + NA2VA1UE1E2: + PA2VA1UE1E2: + AA2VA1UE1E2: + SA2VA1UE1E2: + + NA1VA1: + PA1VA1: + AA1VA1: + SA1VA1: + + NA1VU: + PA1VU: + AA1VU: + SA1VU: + + NA1VE1: + PA1VE1: + AA1VE1: + SA1VE1: + + NA1VE2: + PA1VE2: + AA1VE2: + SA1VE2: + + NA1VA1U: + PA1VA1U: + AA1VA1U: + SA1VA1U: + + NA1VA1E1: + PA1VA1E1: + AA1VA1E1: + SA1VA1E1: + + NA1VA1E2: + PA1VA1E2: + AA1VA1E2: + SA1VA1E2: + + NA1VUE1: + PA1VUE1: + AA1VUE1: + SA1VUE1: + + NA1VUE2: + PA1VUE2: + AA1VUE2: + SA1VUE2: + + NA1VE1E2: + PA1VE1E2: + AA1VE1E2: + SA1VE1E2: + + NA1VA1UE1: + PA1VA1UE1: + AA1VA1UE1: + SA1VA1UE1: + + NA1VA1UE2: + PA1VA1UE2: + AA1VA1UE2: + SA1VA1UE2: + + NA1VA1E1E2: + PA1VA1E1E2: + AA1VA1E1E2: + SA1VA1E1E2: + + NA1VUE1E2: + PA1VUE1E2: + AA1VUE1E2: + SA1VUE1E2: + + NA1VA1UE1E2: + PA1VA1UE1E2: + AA1VA1UE1E2: + SA1VA1UE1E2: + + NEWZE1: + PEWZE1: + AEWZE1: + SEWZE1: + + NEWZE2: + PEWZE2: + AEWZE2: + SEWZE2: + + NEWZE1A1: + PEWZE1A1: + AEWZE1A1: + SEWZE1A1: + + NEWZE2A1: + PEWZE2A1: + AEWZE2A1: + SEWZE2A1: + + NEWZE1U: + PEWZE1U: + AEWZE1U: + SEWZE1U: + + NEWZE2U: + PEWZE2U: + AEWZE2U: + SEWZE2U: + + NEWZE1A1U: + PEWZE1A1U: + AEWZE1A1U: + SEWZE1A1U: + + NEWZE2A1U: + PEWZE2A1U: + AEWZE2A1U: + SEWZE2A1U: + + NEWZE1E2: + PEWZE1E2: + AEWZE1E2: + SEWZE1E2: + + NEWZE1E2A1: + PEWZE1E2A1: + AEWZE1E2A1: + SEWZE1E2A1: + + NEWZE1E2U: + PEWZE1E2U: + AEWZE1E2U: + SEWZE1E2U: + + NEWZE1E2A1U: + PEWZE1E2A1U: + AEWZE1E2A1U: + SEWZE1E2A1U: + + NA1EWZE1: + PA1EWZE1: + AA1EWZE1: + SA1EWZE1: + + NA1EWZE2: + PA1EWZE2: + AA1EWZE2: + SA1EWZE2: + + NA1EWZE1A1: + PA1EWZE1A1: + AA1EWZE1A1: + SA1EWZE1A1: + + NA1EWZE2A1: + PA1EWZE2A1: + AA1EWZE2A1: + SA1EWZE2A1: + + NA1EWZE1U: + PA1EWZE1U: + AA1EWZE1U: + SA1EWZE1U: + + NA1EWZE2U: + PA1EWZE2U: + AA1EWZE2U: + SA1EWZE2U: + + NA1EWZE1A1U: + PA1EWZE1A1U: + AA1EWZE1A1U: + SA1EWZE1A1U: + + NA1EWZE2A1U: + PA1EWZE2A1U: + AA1EWZE2A1U: + SA1EWZE2A1U: + + NA1EWZE1E2: + PA1EWZE1E2: + AA1EWZE1E2: + SA1EWZE1E2: + + NA1EWZE1E2A1: + PA1EWZE1E2A1: + AA1EWZE1E2A1: + SA1EWZE1E2A1: + + NA1EWZE1E2U: + PA1EWZE1E2U: + AA1EWZE1E2U: + SA1EWZE1E2U: + + NA1EWZE1E2A1U: + PA1EWZE1E2A1U: + AA1EWZE1E2A1U: + SA1EWZE1E2A1U: + + NA2EQWZE1U: + PA2EQWZE1U: + AA2EQWZE1U: + SA2EQWZE1U: + + NA2EQWZE2U: + PA2EQWZE2U: + AA2EQWZE2U: + SA2EQWZE2U: + + NA2EQWZE1UA1: + PA2EQWZE1UA1: + AA2EQWZE1UA1: + SA2EQWZE1UA1: + + NA2EQWZE2UA1: + PA2EQWZE2UA1: + AA2EQWZE2UA1: + SA2EQWZE2UA1: + + NA2EQWZE1E2U: + PA2EQWZE1E2U: + AA2EQWZE1E2U: + SA2EQWZE1E2U: + + NA2EQWZE1E2UA1: + PA2EQWZE1E2UA1: + AA2EQWZE1E2UA1: + SA2EQWZE1E2UA1: + + NA2A1QU: + PA2A1QU: + AA2A1QU: + SA2A1QU: + + NA2A1QUA1: + PA2A1QUA1: + AA2A1QUA1: + SA2A1QUA1: + + NA2A1QUE1: + PA2A1QUE1: + AA2A1QUE1: + SA2A1QUE1: + + NA2A1QUE2: + PA2A1QUE2: + AA2A1QUE2: + SA2A1QUE2: + + NA2A1QUA1E1: + PA2A1QUA1E1: + AA2A1QUA1E1: + SA2A1QUA1E1: + + NA2A1QUA1E2: + PA2A1QUA1E2: + AA2A1QUA1E2: + SA2A1QUA1E2: + + NA2A1QUE1E2: + PA2A1QUE1E2: + AA2A1QUE1E2: + SA2A1QUE1E2: + + NA2A1QUA1E1E2: + PA2A1QUA1E1E2: + AA2A1QUA1E1E2: + SA2A1QUA1E1E2: + + NCORECORE: + PCORECORE: + ACORECORE: + SCORECORE: + + NVA2N1: + PVA2N1: + AVA2N1: + SVA2N1: + + NVE1N1: + PVE1N1: + AVE1N1: + SVE1N1: + + NVA1N1: + PVA1N1: + AVA1N1: + SVA1N1: + + NVA2N2: + PVA2N2: + AVA2N2: + SVA2N2: + + NVE1N2: + PVE1N2: + AVE1N2: + SVE1N2: + + NVA1N2: + PVA1N2: + AVA1N2: + SVA1N2: + + NVA2N3: + PVA2N3: + AVA2N3: + SVA2N3: + + NVE1N3: + PVE1N3: + AVE1N3: + SVE1N3: + + NVA1N3: + PVA1N3: + AVA1N3: + SVA1N3: + + NJE1N0: + PJE1N0: + AJE1N0: + SJE1N0: + + NJE1N1: + PJE1N1: + AJE1N1: + SJE1N1: + + NJE1N2: + PJE1N2: + AJE1N2: + SJE1N2: + + NJE1N3: + PJE1N3: + AJE1N3: + SJE1N3: + + NPA2E1N0: + PPA2E1N0: + APA2E1N0: + SPA2E1N0: + + NPA2A1N1: + PPA2A1N1: + APA2A1N1: + SPA2A1N1: + + NPA2E1N2: + PPA2E1N2: + APA2E1N2: + SPA2E1N2: + + NPA2A1N2: + PPA2A1N2: + APA2A1N2: + SPA2A1N2: + + NPA2E1N3: + PPA2E1N3: + APA2E1N3: + SPA2E1N3: + + NPA2A1N3: + PPA2A1N3: + APA2A1N3: + SPA2A1N3: + + NPE1A1N0: + PPE1A1N0: + APE1A1N0: + SPE1A1N0: + + NPE1A1N2: + PPE1A1N2: + APE1A1N2: + SPE1A1N2: + + NPE1A1N3: + PPE1A1N3: + APE1A1N3: + SPE1A1N3: + diff --git a/src/model/weight.f b/src/model/weight.f new file mode 100644 index 0000000..99164c5 --- /dev/null +++ b/src/model/weight.f @@ -0,0 +1,50 @@ +! hybrid,wt_en2ci,wt_en,wt_ci + implicit none +! data arrays and their dimensions + double precision wt(ntot,numdatpt),y(ntot,numdatpt) +! loop index + integer i,j,k,n + + do i=1,numdatpt + wt(1,i)=1.d0 + enddo + + call norm_weight(wt,ntot,numdatpt) + + end + +!---------------------------------------------------------------------------------------------------- +! (q,x1,x2,y,wt,par,p_act,p_spread,nset,npar, + > flag,lauf) + use adia_mod, only: adia + use dim_parameter,only: qn,ntot,numdatpt,ndiab + use ctrans_mod,only: ctrans + implicit none +! IN: variables + integer lauf + integer flag !< 0= initial output 1=fit not converged 2= Fit Converged, 3= max iteration reached + integer npar,nset + double precision par(npar,nset),p_spread(npar) + integer p_act(npar) + double precision q(qn,numdatpt),x1(qn,numdatpt),x2(qn,numdatpt) + double precision y(ntot,numdatpt),wt(ntot,numdatpt) + +! INTERNAL: Variables + integer,parameter :: id_out = 20 , std_out = 6 + integer pt + integer i, id_print + double precision, allocatable :: ymod(:,:) + double precision, allocatable :: ew(:,:) + double precision, allocatable :: ev(:,:,:) + + logical skip + + allocate(ymod(ntot,numdatpt)) + allocate(ew(ndiab,numdatpt)) + allocate(ev(ndiab,ndiab,numdatpt)) + + skip=.false. + +! get Model Outputs for all geometries for current best parameter set par(:,1) + do pt=1,numdatpt + call adia(pt,par(1:npar,1),npar,ymod(1:ntot,pt), + > ew(1:ndiab,pt),ev(1:ndiab,1:ndiab,pt),skip) + call ctrans(q(:,pt),x1(:,pt),x2(:,pt)) + enddo + +! Initial write print everything you want to see before the fit and return + if(flag.eq.0) then + call print_parameterstate(std_out,par(:,1),p_act,npar) + call print_ErrorSummary(std_out,y,ymod,wt) +! print Data into the plotfiles + return + endif +! open output files for individual makro iterations + call open_outfile(id_out,lauf) +! print Data into the plotfiles + call print_plotfiles(x1,y,wt,ymod) + +! print Genetic output into files + do i=1, 2 + if (i.eq.1) then + id_print= std_out + else + id_print= id_out + endif + write(id_print,'("Writing Iteration: ",i4)') lauf + write(id_print,block_line) +! write data information only in outfile + if(i.eq.2) then + call print_data(id_print,x1,y,ymod,wt) + call print_Set_Errors(id_print,y,ymod,wt) + endif + call print_parameterblock + > (id_print,par(:,1),p_act,p_spread,npar) + call print_ErrorSummary(id_print,y,ymod,wt) + + enddo + + call print_fortranfile(par(:,1),npar) + + ! write the type of calc at the end of the output + + + close (id_out) + deallocate(ymod,ev,ew) + end subroutine +!---------------------------------------------------------------------------------------------------- +! + subroutine print_Set_Errors(id_out,y, ymod, wt) + use io_parameters,only: llen + use dim_parameter,only: ndata,nstat,ntot,numdatpt,sets + integer , intent(in) :: id_out + double precision, intent(in) :: y(ntot,numdatpt), + > ymod(ntot,numdatpt), wt(ntot,numdatpt) + integer :: set, setpoint, pt + double precision :: Set_rms(sets,ntot), Set_num(sets,ntot) + double precision :: Total_rms, Total_Energy_rms,Energy_rms(nstat) + character(len=llen) fmt + write(id_out,'(A)') 'Errors in icm for individual Sets' // + > '(specified by sets: and npoints:)' + write(id_out,'(A5,3A16)')'Set','Total', + > 'Total_Energy', 'Energy[nstat]' + write(id_out,sep_line) + write(fmt,'("(I5,2f16.1,",I2,"f16.1)")') nstat + Set_rms = 0.d0 + pt = 0 + do set=1, sets + do setpoint=1, ndata(set) + pt = pt + 1 + where(wt(:,pt) > 0.d0) + Set_rms(set,:) = Set_rms(set,:)+(ymod(:,pt)-y(:,pt))**2 + Set_num(set,:) = Set_num(set,:) + 1 + end where + enddo + Total_rms + > = dsqrt(sum(Set_rms(set,:)) + > / (sum(Set_num(set,:)))) + Total_Energy_rms + > = dsqrt(sum(Set_rms(set,1:nstat)) + > / (sum(Set_num(set,1:nstat)))) + Energy_rms(1:nstat) + > = dsqrt(Set_rms(set,1:nstat) + > / (Set_num(set,1:nstat))) + write(id_out,fmt) set, Total_rms*h2icm, Total_Energy_rms*h2icm, + > Energy_rms(1:nstat)*h2icm + enddo + write(id_out,block_line) + write(id_out,*) '' + end subroutine print_Set_Errors + +!---------------------------------------------------------------------------------------------------- +! wt(:,sstart:ssend),ymod(:,sstart:ssend), + > ndata(set),id_plot,set) + else + call print_plotcoord(plot_coord(set), + > x(:,sstart:ssend),y(:,sstart:ssend), + > wt(:,sstart:ssend),ymod(:,sstart:ssend), + > ndata(set),id_plot,set) + endif + close(id_plot) + enddo + + end subroutine + +!---------------------------------------------------------------------------------------------------- + subroutine print_plotwalk(x,y,wt,ymod,npt,id_plot,set) + use dim_parameter,only: qn,ntot + use io_parameters,only: llen + implicit none +! IN: variables + integer id_plot,npt,set + double precision x(qn,npt),y(ntot,npt),ymod(ntot,npt),wt(ntot,npt) +! INTERNAL: variables + double precision xdiff(qn),walktime + double precision walknorm +! loop control + integer i,j + + character(len=llen) fmt + j=ntot-1 + + call print_plotheader(id_plot,0,npt,set) + + call getwalknorm(x,walknorm,npt) + walktime = 0.d0 + do i=1,npt + if(i.gt.1) then + xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1) + walktime = walktime + dsqrt(sum(xdiff(1:qn)**2))/walknorm + endif + write(id_plot,"(ES16.8,*(3(ES16.8),:))") + > walktime ,ymod(:,i),y(:,i),(wt(:,i)) + enddo + + end subroutine + +!---------------------------------------------------------------------------------------------------- + subroutine print_plotcoord(coord,x,y,wt,ymod,npt,id_plot,set) + use dim_parameter,only: qn,ntot + use io_parameters,only: llen + implicit none +! IN: variables + integer, intent(in) :: id_plot,npt,set,coord + double precision, intent(in) :: x(qn,npt),y(ntot,npt) + double precision, intent(in) :: ymod(ntot,npt),wt(ntot,npt) +! loop control + integer i + + call print_plotheader(id_plot,coord,npt,set) + do i=1,npt + write(id_plot,"(ES16.8,*(3(ES16.8),:))") + > x(coord,i), ymod(:,i),y(:,i),(wt(:,i)) +! write(id_plot,"(2ES16.8,*(3(ES16.8),:))") +! > x(coord,i), x(coord+1,i),y(:,i) + enddo + + end subroutine + +!---------------------------------------------------------------------------------------------------- + subroutine print_plotheader(id_plot,coord,npt,set) + use dim_parameter,only: qn,ntot + use io_parameters,only: llen + implicit none + integer, intent(in) :: id_plot,npt,set,coord + + character(len=llen) fmt + + write(id_plot,'("#SET: ",i5)') set + write(id_plot,'("#OUTPUT VALUES",i4)') ntot + write(id_plot,'("#DATA POINTS: ",i4)') npt + if (coord.le.0) then + write(id_plot,'("#t(x) = WALK")') + else + write(id_plot,'("#t(x) = x(",I0,")")') coord + endif + write(id_plot,'("#UNIT: hartree")') + write(id_plot,'()') + write(id_plot,'("#",A15)',advance='no') "t(x)" + write(fmt,'("(3(7X,A9,",I3,"(16x)))")') ntot-1 + write(id_plot,fmt) 'ymod(p,x)','y(x) ','wt(x) ' + + + end subroutine + +!---------------------------------------------------------------------------------------------------- +! 'ERROR: No rule for Outputfile naming for MAXIT >= 1000' + stop + endif + + open (id_out,file=outname) + + end subroutine + +!---------------------------------------------------------------------------------------------------- +! =2 +! <@param id_out specifies the file in which the Parameters are Printed +! <@param p vector containing one set of parameter values +! <@param p_act vector containing the active state 0 (inactive) or 1 (active) for each parameter +! <@param p_spread vector containing the spreads for each parameter +! <@param npar lenght of the parmeter vectors (p,p_act,p_spread) +! <@TODO extract subroutine for printing the multiline values, would make this more readable + subroutine print_parameterblock(id_out,p,p_act,p_spread,npar) + use dim_parameter,only: pst, facspread + use io_parameters,only: key, parkeynum,parkeylen,llen + implicit none +! IN: Variables + integer id_out,npar,p_act(npar) + double precision p(npar),p_spread(npar) + +! INTERNAL: variables +! loop index + integer i,k,l,t,n !< internal variables for loops and positions in parameter vectors + +! number of values per line, values must be atleast 2 set this to personal preference + integer, parameter :: np=5,nsp=5 + + character(len=llen) fmt + + +! Write header for Parameter block + 1 format('!',200('=')) + write(id_out,1) + write(id_out,'(A2,5x,A11,i3)') '! ','PARAMETER: ',npar + write(id_out,1) + +! loop over all Parameter Keys + do i = 1, parkeynum +! save start and end of parameter block for specific key + k = pst(1,i) + l = pst(1,i)+pst(2,i)-1 +! print only used keys with atleast one parameter + if(pst(2,i).gt.0) then + write(fmt,'("(a",I3,"'' ''i3)")') parkeylen + write(id_out,fmt) adjustl(key(1,i)), pst(2,i) + +! write the actual parameters -> subroutine print_parameterlines()? + if(l-k.le.(np-1)) then + write(fmt,'("(a",I3,"'' ''",I3,"g24.15)")') parkeylen,np + write(id_out,fmt) key(2,i),(p(n), n=k,l) + + else +! start of multi line parameter print, number of values per line specified by np + write(fmt,'("(a",I3,"'' ''",I3,"g24.15'' &'')")') + $ parkeylen,np + write(id_out,fmt) key(2,i),(p(n), n=k,k+(np-1)) + + t=k+np +! write continuation lines till left parameters fit on last line + do while(t.le.l) + if(l-t.le.(np-1)) then + write(fmt,'("(",I3,"x'' ''",I3,"g24.15)")') + $ parkeylen,np + write(id_out,fmt) (p(n), n=t, l) + + else + write(fmt,'("(",I3,"x'' ''",I3,"g24.15'' &'')")') + $ parkeylen,np + write(id_out,fmt) (p(n), n=t, t+(np-1)) + + endif + t=t+np + enddo + + endif !-> end subroutine print_parameterlines + +! write parameter active state in one line + write(fmt,'("(a",I3,"'' ''","50i3)")') parkeylen + write(id_out,fmt) key(3,i),(p_act(n),n=k,l) + +! write the spreads for each parameter + if(l-k.le.(np-1)) then + write(fmt,'("(a",I3,"'' ''",I3,"g24.8)")') parkeylen,nsp + write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,l) + + else +! start of multiline spread values + write(fmt,'("(a",I3,"'' ''",I3,"g24.8'' &'')")') + $ parkeylen,nsp + write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,k + > +(np-1)) + + t=k+nsp +! write continuation lines till left spreads fit on last line + do while(t.le.l) + if(l-t.le.(np-1)) then + write(fmt,'("(",I3,"x'' ''",I3,"g24.8)")') + $ parkeylen,nsp + write(id_out,fmt) (p_spread(n)/facspread, n=t, l) + else + write(fmt,'("(",I3,"x'' ''",I3,"g24.8'' &'')")') + $ parkeylen,nsp + write(id_out,fmt) (p_spread(n)/facspread, n=t, t + > +(np-1)) + + endif + t=t+np + enddo + + endif +! print empty line between diffrent parameter blocks for better readability + write(id_out,'(" ")') + endif + + enddo + + end subroutine + +!---------------------------------------------------------------------------------------------------- +! adjustl('y(x)-ymod(x)'),adjustl('weight'), + > adjustl('x(1:qn_read) ') + write(id_out,sep_line) +! loop over all datapoints for each set and count the actual datapointnumber with point + point=0 + do i=1,sets + write(id_out,18) 'Set: ', i + do j=1,ndata(i) + write(id_out,18) 'Point: ', j + point=point+1 +! print all data for one datapoint + call print_datapoint(id_out,x(:,point),y(:,point), + > ymod(:,point),wt(:,point)) + write(id_out,sep_line) + + enddo + enddo +! write end of data statement and two seperating lines + write(id_out,block_line) + write(id_out,*) '' + end subroutine +!---------------------------------------------------------------------------------------------------- +! wt(k), x(1:qn_read) + enddo + enddo + + end subroutine + + + + end module write_mod