From 464dca1bbb7c488dfa727068a8df0fb0fa48e25c Mon Sep 17 00:00:00 2001 From: jean paul nshuti Date: Wed, 8 Oct 2025 15:52:19 +0200 Subject: [PATCH] remove PES --- src/model/new_PES/.ctrans.f90.un~ | Bin 24607 -> 0 bytes src/model/new_PES/.xy_stretch_lib.f90.un~ | Bin 13220 -> 0 bytes src/model/new_PES/.xyumb_stretch_lib.f90.swp | Bin 12288 -> 0 bytes src/model/new_PES/.xyumb_stretch_lib.f90.un~ | Bin 16728 -> 0 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 - 15 files changed, 5543 deletions(-) delete mode 100644 src/model/new_PES/.ctrans.f90.un~ delete mode 100644 src/model/new_PES/.xy_stretch_lib.f90.un~ delete mode 100644 src/model/new_PES/.xyumb_stretch_lib.f90.swp delete mode 100644 src/model/new_PES/.xyumb_stretch_lib.f90.un~ delete mode 100644 src/model/new_PES/d3h_umb_stretch_lib.f90 delete mode 100644 src/model/new_PES/keys.f90 delete mode 100644 src/model/new_PES/pes_ctrans.f90 delete mode 100644 src/model/new_PES/pes_model.f90 delete mode 100644 src/model/new_PES/plot_surface.f90 delete mode 100644 src/model/new_PES/select_monom_mod.f90 delete mode 100644 src/model/new_PES/sphericalharmonics_mod.f90 delete mode 100644 src/model/new_PES/surface_mod.f90 delete mode 100644 src/model/new_PES/xy_stretch_lib.f90 delete mode 100644 src/model/new_PES/xy_stretch_lib.f90~ delete mode 100644 src/model/new_PES/xyumb_stretch_lib.f90 diff --git a/src/model/new_PES/.ctrans.f90.un~ b/src/model/new_PES/.ctrans.f90.un~ deleted file mode 100644 index 91f036ba43e0c5635c63714455278e96835ad46a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 24607 zcmeHPUymEfaX$x%9XS3I+wp(LZeA?yaF@GUlB?BfZ8_2Dbk?0uCuOB`U|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 diff --git a/src/model/new_PES/.xy_stretch_lib.f90.un~ b/src/model/new_PES/.xy_stretch_lib.f90.un~ deleted file mode 100644 index 9327c3e936e39e527e6c027ff3d28e3865cf11ae..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 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 diff --git a/src/model/new_PES/.xyumb_stretch_lib.f90.un~ b/src/model/new_PES/.xyumb_stretch_lib.f90.un~ deleted file mode 100644 index 168ce0965bd90d159e8efef6ec2717b6a10b47a2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 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) diff --git a/src/model/new_PES/d3h_umb_stretch_lib.f90 b/src/model/new_PES/d3h_umb_stretch_lib.f90 deleted file mode 100644 index 6d1cf0a..0000000 --- a/src/model/new_PES/d3h_umb_stretch_lib.f90 +++ /dev/null @@ -1,82 +0,0 @@ -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 deleted file mode 100644 index b66ee29..0000000 --- a/src/model/new_PES/keys.f90 +++ /dev/null @@ -1,530 +0,0 @@ - !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 deleted file mode 100644 index d4b7ba5..0000000 --- a/src/model/new_PES/pes_ctrans.f90 +++ /dev/null @@ -1,671 +0,0 @@ -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % SUBROUTINE CTRANS(...) -! % -! % M. Vossel 21.03.2023 -! % -! % Routine to transform symmetryinput coordinates to symmetrized -! % coordinates. Distances Are discribet by Morse coordinates or -! % TMC depending on Set Parameters in the Genetic Input. -! % -! % input variables -! % q: -! % q(1): H1x -! % q(2): y -! % q(3): z -! % q(4): H2x -! % q(5): y -! % q(6): z -! % q(7): H3x -! % q(8): y -! % q(9): z -! -! -! -! % Internal variables: -! % t: primitive coordinates (double[qn]) -! % t(1): -! % t(2): -! % t(3): -! % t(4): -! % t(5): -! % t(6): -! % t(7): -! % t(8): -! % t(9): -! % t: dummy (double[qn]) -! % p: parameter vector -! % npar: length of parameter vector -! % -! % Output variables -! % s: symmetrized coordinates (double[qn]) -! % s(1): CH-symetric streatch -! % s(2): CH-asymetric streatch-ex -! % s(3): CH-asymetric streatch-ey -! % s(4): CH-bend-ex -! % s(5): CH-bend-ey -! % s(6): CH-umbrella -! % s(7): CH-umbrella**2 -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -module ctrans_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 deleted file mode 100644 index 4acc849..0000000 --- a/src/model/new_PES/pes_model.f90 +++ /dev/null @@ -1,982 +0,0 @@ -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 deleted file mode 100644 index 47af619..0000000 --- a/src/model/new_PES/sphericalharmonics_mod.f90 +++ /dev/null @@ -1,575 +0,0 @@ -! Module contains the spherical harmonics up to l=5 m=-l,..,0,..,l listed on https://en.wikipedia.org/wiki/Table_of_spherical_harmonics from 19.07.2022 -! the functions are implementde by calling switch case function for given m or l value and return the corresdpondig value for given theta and phi -! the functions are split for diffrent l values and are named by P_lm. -! example for l=1 and m=-1 the realpart of the spherical harmonic for given theta and phi -! is returned by calling Re_Y_lm(1,-1,theta,phi) which itself calls the corresponding function P_1m(m,theta) and multilpies it by cos(m*phi) to account for the real part of exp(m*phi*i) -! Attention the legendre polynoms are shifted to account for the missing zero order term in spherical harmonic expansions -module sphericalharmonics_mod - use accuracy_constants, only: dp, idp - implicit none - real(kind=dp), parameter :: PI = 4.0_dp * atan( 1.0_dp ) -contains - !---------------------------------------------------------------------------------------------------- - function Y_lm( l , m , theta , phi ) result( y ) - integer(kind=idp), intent( in ) :: l , m - real(kind=dp), intent( in ) :: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( l ) - case (1) - y = Y_1m( m , theta , phi ) - case (2) - y = Y_2m( m , theta , phi ) - case (3) - y = Y_3m( m , theta , phi ) - case (4) - y = Y_4m( m , theta , phi ) - case (5) - y = Y_5m( m , theta , phi ) - case default - write(errmesg,'(A,i0)')& - &'order of spherical harmonics not implemented', l - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Y_lm - !---------------------------------------------------------------------------------------------------- - function Re_Y_lm( l , m , theta , phi ) result( y ) - integer(kind=idp), intent( in ) :: l , m - real(kind=dp), intent( in ) :: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( l ) - case (1) - y = P_1m( m , theta ) * cos(m*phi) - case (2) - y = P_2m( m , theta ) * cos(m*phi) - case (3) - y = P_3m( m , theta ) * cos(m*phi) - case (4) - y = P_4m( m , theta ) * cos(m*phi) - case (5) - y = P_5m( m , theta ) * cos(m*phi) - case (6) - y = P_6m( m , theta ) * cos(m*phi) - case default - write(errmesg,'(A,i0)')& - &'order of spherical harmonics not implemented', l - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Re_Y_lm - !---------------------------------------------------------------------------------------------------- - function Im_Y_lm( l , m , theta , phi ) result( y ) - integer(kind=idp), intent( in ) :: l , m - real(kind=dp), intent( in ) :: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( l ) - case (1) - y = P_1m( m , theta ) * sin(m*phi) - case (2) - y = P_2m( m , theta ) * sin(m*phi) - case (3) - y = P_3m( m , theta ) * sin(m*phi) - case (4) - y = P_4m( m , theta ) * sin(m*phi) - case (5) - y = P_5m( m , theta ) * sin(m*phi) - case (6) - y = P_6m( m , theta ) * sin(m*phi) - case default - write(errmesg,'(a,i0)')& - &'order of spherical harmonics not implemented',l - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Im_Y_lm - - !---------------------------------------------------------------------------------------------------- - !---------------------------------------------------------------------------------------------------- - function Y_1m( m , theta , phi ) result( y ) - integer(kind=idp),intent( in ):: m - real(kind=dp),intent( in ):: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-1) - y = 0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)*cos(phi) - - case (0) - y = 0.5_dp*sqrt(3.0_dp/PI)*cos(theta) - - case (1) - y = -0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)*cos(phi) - - case default - write(errmesg,'(a,i0)') 'in y_1m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Y_1m - !---------------------------------------------------------------------------------------------------- - function Y_2m(m,theta,phi) result(y) - integer(kind=idp),intent(in):: m - real(kind=dp),intent(in):: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-2) - y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(2.0_dp*phi) - - case (-1) - y=0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta)*cos(phi) - - case (0) - y=0.25_dp*sqrt(5.0_dp/PI)& - &*(3.0_dp*cos(theta)**2-1.0_dp) - - case (1) - y=-0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta)*cos(phi) - - case (2) - y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(2.0_dp*phi) - - case default - write(errmesg,'(A,i0)')'in y_2m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_2m - !---------------------------------------------------------------------------------------------------- - function Y_3m(m,theta,phi) result(y) - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-3) - y=0.125_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(3.0_dp*phi) - - case (-2) - y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta)*cos(2.0_dp*phi) - - case (-1) - y=0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp)*cos(phi) - - case (0) - y=0.25_dp*sqrt(7.0_dp/PI)& - &*(5.0_dp*cos(theta)**3-3.0_dp*cos(theta)) - - case (1) - y=-0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp)*cos(phi) - - case (2) - y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta)*cos(2.0_dp*phi) - - case (3) - y=-0.125_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(3.0_dp*phi) - - case default - write(errmesg,'(A,i0)')'in y_3m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_3m - !---------------------------------------------------------------------------------------------------- - function Y_4m(m,theta,phi) result(y) - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(4.0_dp*phi) - - case (-3) - y=0.375_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta)*cos(3.0_dp*phi) - - case (-2) - y=0.375_dp*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2& - &*(7.0_dp*cos(theta)**2-1)*cos(2.0_dp*phi) - - case (-1) - y=0.375_dp*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7.0_dp*cos(theta)**3& - &-3.0_dp*cos(theta))*cos(phi) - - case (0) - y=(3.0_dp/16.0_dp)/sqrt(PI)& - &*(35.0_dp*cos(theta)**4& - &-30.0_dp*cos(theta)**2+3.0_dp) - - case (1) - y=-0.375_dp*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7.0_dp*cos(theta)**3& - &-3.0_dp*cos(theta))*cos(phi) - - case (2) - y=0.375_dp*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(7.0_dp*cos(theta)**2-1.0_dp)& - &*cos(2*phi) - - case (3) - y=-0.375_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta)*cos(3.0_dp*phi) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(4.0_dp*phi) - - case default - write(errmesg,'(a,i0)')'in y_4m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_4m - !---------------------------------------------------------------------------------------------------- - function Y_5m(m,theta,phi) result(y) - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-5) - y=(3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5*cos(5*phi) - - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta)*cos(4*phi) - - case (-3) - y=(1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9*cos(theta)**2-1.0_dp)*cos(3*phi) - - case (-2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta))*cos(2*phi) - - case (-1) - y=(1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21.0_dp*cos(theta)**4& - &-14.0_dp*cos(theta)**2+1.0_dp)*cos(phi) - - case (0) - y=(1.0_dp/16.0_dp)/sqrt(11.0_dp/PI)& - &*(63.0_dp*cos(theta)**5-70.0_dp*cos(theta)**3& - &+15.0_dp*cos(theta)) - - case (1) - y=(-1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21.0_dp*cos(theta)**4& - &-14.0_dp*cos(theta)**2+1.0_dp)*cos(phi) - - case (2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta))*cos(2*phi) - - case (3) - y=(-1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9.0_dp*cos(theta)**2-1.0_dp)& - &*cos(3.0_dp*phi) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta)*cos(4.0_dp*phi) - - case (5) - y=(-3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5*cos(5.0_dp*phi) - - case default - write(errmesg,'(A,i0)')'in y_5m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_5m - !---------------------------------------------------------------------------------------------------- - function P_1m( m , theta ) result( y ) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=1 and given m and theta - integer(kind=idp),intent( in ):: m - real(kind=dp),intent( in ):: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-1) - y = 0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta) - - case (0) - y = 0.5_dp*sqrt(3.0_dp/PI)*(cos(theta)-1.0_dp) ! -1 is subtracted to shift so that for theta=0 y=0 - - case (1) - y = -0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta) - - case default - write(errmesg,'(A,i0)')'in p_1m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_1m - !---------------------------------------------------------------------------------------------------- - function P_2m( m , theta ) result( y ) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=2 and given m and theta - integer(kind=idp),intent(in):: m - real(kind=dp),intent(in):: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-2) - y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2 - - case (-1) - y=0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta) - - case (0) - y = (3.0_dp*cos(theta)**2-1.0_dp) - y = y - 2.0_dp !2.0 is subtracted to shift so that for theta=0 y=0 - y = y * 0.25_dp*sqrt(5.0_dp/PI) ! normalize - - case (1) - y = -0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta) - - case (2) - y = 0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2 - - case default - write(errmesg,'(A,i0)')'in p_2m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_2m - !---------------------------------------------------------------------------------------------------- - function P_3m( m , theta ) result( y ) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=3 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-3) - y=0.125_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3 - - case (-2) - y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta) - - case (-1) - y=0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5*cos(theta)**2-1.0_dp) - - case (0) - y=(5.0_dp*cos(theta)**3-3*cos(theta)) - y=y-2.0_dp ! 2.0 is subtracted to shift so that for theta=0 y=0 - y=y*0.25_dp*sqrt(7.0_dp/PI) ! normalize - - case (1) - y=-0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp) - - case (2) - y=0.25*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta) - - case (3) - y=-0.125*sqrt(35.0_dp/PI)& - &*sin(theta)**3 - - case default - write(errmesg,'(A,i0)')'in p_3m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_3m - !---------------------------------------------------------------------------------------------------- - function P_4m(m,theta) result(y) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=4 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4 - - case (-3) - y=0.375*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta) - - case (-2) - y=0.375*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(7*cos(theta)**2-1) - - case (-1) - y=0.375*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7*cos(theta)**3-3*cos(theta)) - - - case (0) - y=(35*cos(theta)**4-30*cos(theta)**2+3) - y = y - 8.0_dp !8.0 is subtracted to shift so that for theta=0 y=0 - y = y * (3.0_dp/16.0_dp)/sqrt(PI) - - case (1) - y=-0.375*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7*cos(theta)**3-3*cos(theta)) - - case (2) - y=0.375*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(7*cos(theta)**2-1) - - case (3) - y=-0.375*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4 - - case default - write(errmesg,'(A,i0)')'in p_4m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_4m - !---------------------------------------------------------------------------------------------------- - function P_5m(m,theta) result(y) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=5 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-5) - y=(3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5 - - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta) - - case (-3) - y=(1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9*cos(theta)**2-1.0_dp) - - case (-2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta)) - - case (-1) - y=(1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21*cos(theta)**4-14*cos(theta)**2+1) - - - case (0) - y = (63*cos(theta)**5-70*cos(theta)**3+15*cos(theta)) - y = y - 8.0_dp !8.0 is subtracted to shift so that for theta=0 y=0 - y = y * (1.0_dp/16.0_dp)/sqrt(11.0_dp/PI) - - case (1) - y=(-1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21*cos(theta)**4-14*cos(theta)**2+1) - - case (2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta)) - - case (3) - y=(-1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9*cos(theta)**2-1.0_dp) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta) - - case (5) - y=(-3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5 - - case default - write(errmesg,'(A,i0)')'in p_5m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_5m - !---------------------------------------------------------------------------------------------------- - function P_6m(m,theta) result(y) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=6 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp):: y - character(len=70) :: errmesg - select case ( m ) - case (-6) - y = (1.0_dp/64.0_dp)*sqrt(3003.0_dp/PI)& - &* sin(theta)**6 - case (-5) - y = (3.0_dp/32.0_dp)*sqrt(1001.0_dp/PI)& - &* sin(theta)**5& - &* cos(theta) - case (-4) - y= (3.0_dp/32.0_dp)*sqrt(91.0_dp/(2.0_dp*PI))& - &* sin(theta)**4& - &* (11*cos(theta)**2 - 1 ) - case (-3) - y= (1.0_dp/32.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**3& - &* (11*cos(theta)**3 - 3*cos(theta) ) - case (-2) - y= (1.0_dp/64.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**2& - &* (33*cos(theta)**4 - 18*cos(theta)**2 + 1 ) - case (-1) - y= (1.0_dp/16.0_dp)*sqrt(273.0_dp/(2.0_dp*PI))& - &* sin(theta)& - &* (33*cos(theta)**5 - 30*cos(theta)**3 + 5*cos(theta) ) - case (0) - y = 231*cos(theta)**6 - 315*cos(theta)**4 + 105*cos(theta)**2-5 - y = y - 16.0_dp !16.0 is subtracted to shift so that for theta=0 y=0 - y = y * (1.0_dp/32.0_dp)*sqrt(13.0_dp/PI) - case (1) - y= -(1.0_dp/16.0_dp)*sqrt(273.0_dp/(2.0_dp*PI))& - &* sin(theta)& - &* (33*cos(theta)**5 - 30*cos(theta)**3 + 5*cos(theta) ) - case (2) - y= (1.0_dp/64.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**2& - &* (33*cos(theta)**4 - 18*cos(theta)**2 + 1 ) - case (3) - y= -(1.0_dp/32.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**3& - &* (11*cos(theta)**3 - 3*cos(theta) ) - case (4) - y= (3.0_dp/32.0_dp)*sqrt(91.0_dp/(2.0_dp*PI))& - &* sin(theta)**4& - &* (11*cos(theta)**2 - 1 ) - case (5) - y= -(3.0_dp/32.0_dp)*sqrt(1001.0_dp/PI)& - &* sin(theta)**5 * cos(theta) - case (6) - y = (1.0_dp/64.0_dp)*sqrt(3003.0_dp/PI)& - &* sin(theta)**6 - case default - write(errmesg,'(A,i0)')'in p_6m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_6m - !---------------------------------------------------------------------------------------------------- - -end module diff --git a/src/model/new_PES/surface_mod.f90 b/src/model/new_PES/surface_mod.f90 deleted file mode 100644 index 885c02c..0000000 --- a/src/model/new_PES/surface_mod.f90 +++ /dev/null @@ -1,71 +0,0 @@ -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 deleted file mode 100644 index 9b30191..0000000 --- a/src/model/new_PES/xy_stretch_lib.f90 +++ /dev/null @@ -1,82 +0,0 @@ -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~ deleted file mode 100644 index db8cc78..0000000 --- a/src/model/new_PES/xy_stretch_lib.f90~ +++ /dev/null @@ -1,82 +0,0 @@ -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 deleted file mode 100644 index cc8073e..0000000 --- a/src/model/new_PES/xyumb_stretch_lib.f90 +++ /dev/null @@ -1,83 +0,0 @@ - -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