From 15082646ca52fc41e1e897cfd62c7fc182dcc327 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@users.sourceforge.net>
Date: Fri, 11 Feb 2011 09:11:15 +0100
Subject: [PATCH] Acceleration of multi-det

---
 bin/to_ezfio.exe        | Bin 96866 -> 97535 bytes
 src/det.irp.f           | 173 ++++++++++++++++++++++------------
 src/eplf_function.irp.f | 203 +++++++++++++++++++++++++++++-----------
 src/mo.irp.f            |  10 --
 src/overlap.irp.f       |   8 +-
 src/test_1d.irp.f       |   2 +-
 6 files changed, 269 insertions(+), 127 deletions(-)

diff --git a/bin/to_ezfio.exe b/bin/to_ezfio.exe
index 2d1d7450166b127868818a11446f1e4c868c1bc0..222b347afd41dd56de145bbd6da66e8afe21ee9b 100755
GIT binary patch
delta 8181
zcmZvBWl$YJv*o=IT!Tw+4esv2-Q696I|R76yItI!i+h0J?(Ptr;7;((SNpbhU+tUe
zAJbFSJvCKbb9&C%y@UC9f~gY}`s`IrDU?U6Y%E0y0Ni#10A$HMrnrq3LLQ``L0j%G
z0!*u(rQxUWs}J)^kQ|JsxC@A3ik(GBlfQlMT8&C&!_t|^-J21L%c6a-o^B<E7yNcO
zbUCbGrM<C7<T+=cM_8^WUIB|edgNgdwO++Kq3<!r2^H$ZPY17!t1nrabC(1+z(cwb
zvQd1kAS&MJo23U6bLV)l)u;*RyW<%&bjvsQ8m6$OJdu!04Xbd6yiUi5X&|xEK+$9D
z;!5<;$}SDXpyuO(WM1mw<*^`t5J=uYf1@-%R<U;5HaG1H%ckK6(M%LjeNB*9_+Z<j
z0Y;CP_RO7QD0O2~VI1)DvCQpGxeN*GjrTG$*G`iCYW16cYPT!}4*LnH!D;zy89w^>
zb%M1Pi)D-6uy{|=wfK&S)%aMq43~f9WplbZIZ0$Pj<{SOC&U@x9E+xr?#6DJC!`G3
z_p+LG_n91%tSmZf6&eqF@5XCOov(eCl&=Kc8nHM>a$SM>n|>IlpdC#RfeY0HjnYEV
zk4)NvschHmlgXZAwog8wD~6$H?#8)jK~lHwnnTrSm*>al-_kjkcw!mMJtnDRBDql;
zZ3bc~bKT2uweeB|u~|SuLdmLy;TOZl&kW3+?C5ER3*EJEGjv^TM-Bn+L&Zdwor@eb
z6OQx73MYaOsd&+5c9vm>6LJFfcDj|hT;Y?dr2uS<4W7&%Mx8=Xw;iTv3Jq3a^sU>o
zeCI0C5q2iA9`3aaQL65NQ<Kp}+T)^BK7VsX>k!4avDiDhD8J6APeKQ87yeJ%^b_(&
zE_i+L8Uy*QdO{Uos!fY$1FPVz?z^|J70<WxFjym_1vE1;!_Jz;qboRU%bJg`iMW)j
z5Zw3NI@o1HF{2+K9X;QTp7d%G&C0v32)7gb^F}0#rym%hQ?ttCehSd&R9nfkSaF-G
z{TH|4m#;R@icK#1)KM-|edDw&o>V#`8U4c6SIB-(4V_3)tM3IUqkj~oi6J3pVP+0A
z`*u=dkd(8q{UOECr<xRPnN^}_?8%qSS?%kl*1cIaXH`6q`#3RXxlSOIEou+<bu@|7
zMgC=TQSDVeG_a=$DbDMYbM<r`BBVaZJgZOi3uJddPvNc0U88fIlCZ<-kxLe{Ny}%r
zP%R7FFDq&#<8J>{CB^`?QOVJYIc{s|B>Wclo}!72<@{o~>XZPKM>-cn=9a-@aV;IK
zNW@k{ttAO+`XPH3&zL-GJlc_JyxGd6HO}d#oTMGr^3@|OtyvUW>$mZ=ae?IQnx;7&
za)3Yt^&q~>lL&8LM$Q{D1OBuIJzRj*U*qyEtn}P>r~m-Pmz6))>6Oo3op(WtUK59$
zuQk0NB+7EopRfP`06gHn#e!dszut8p4FE9W1^^nnO-mI%uLbi9{VqN+=OY3D*g^mR
zivPS1;2C!yT2Xc#GKl<UWW!;Z<!ZC4b${{0t7|AgDg9!2fs_;Ub5<3_7wSt$32FF@
zN=C${9^0s|08KZBpLiLLqKYm>rZ5i_o#Zb31PS+jLp8o>;<VEyv-dq5MHfqEZDED{
z?;IHbPl5ekGRKe%b>s7;L7grC=ZpzP3f#rmFG|<brLy@X_8^&{7IhN&>pfv+Gk5oD
zPM7oO1mQM*HNDhzcr>X7+)S)WeFNPa_8<=_s;$)sp44)@)KT>HB@xrRyr~FnWrnZW
zY_ggJX*C1FdgaKB#JhQ;or?YxjKL*^jyU=%FOC5`gvs8)nfbkjv+*e$x>|2|GMqb?
z1;jR3T2vBn)*w}0{;n_CkF|q|p$K<>29$cOqJs&Z4p*f<4(P?TY}9l=no!Z+Lc=$I
zTHhqdQiC9uv*pDR*5dFH>Ua4Vk&BDweIa3_yu6$|L(BaL3E{E{ho6u;Sb@+Xpa&_K
z0{`sUlYRUxK>0P;Hro`z%|>D-3HD7%17IC{{aN|340L4CT-nt7%NoYCk!bkU8_Ztq
z!>fl2Kov(x+cNpy5z8*9o$PEAJ~~cZA{rcv9N&$j+JJQAJP;l}mbxR2!@T0q-##a;
zKFqbL9H>Co+Lru^`thVsb{eUZohg5gh+r>sGmt+rbycT5tsdVeSQV}DkfEvAy*+p7
z>GL*l07AFOq+uKow7F~${&A-iF6jQiWi$JS#P+izh9l<p5h6eDal&EU$ERDy7ii{c
zV#9X;0#mMI(&*RJ(XYw>Wj5i0=Z$4=v6OBPaFll$NC)ipy!QCGb2)$M%{@CJ%>gm<
zqV42-xLn-jda@oUM4$G?3e+834FZL_L88B{)PGya|CjAf3w<<dZH<hr0-*(FIBv%(
zGti<e4L|quxxP%4Fzi!u@@frKdjuV%#`g0U{G85nEjk?QN<Db?{$|P>-CnSPSjk)8
z5gU2H!BcrG)0*_Q=|H}Adtb<Ct_RNzb|BkYSI9VO@hy{XLw9B}quy@@VL684q&v8-
z<sqG8pRA{w&e}U7@r6Dh-8-(Gh*`GRAX;o<FO5>xr;*z)890defTmI-5Xxv3aSZWZ
z%cfOz%;f30@Io+B_CI6N6Yee%Z&J`ce1p=nVUJK19<_}}`HJepq3JSO#I-$Fwf&vM
zK~0bf+?}B2IW&f&G#nWLQm-~AkisuLIW`hVz?E9-1j@)qr(|0A$nH7jlkGa0)r^#t
zr;L*iS6uKk&SGV#rTJJ@sL_APGgDmb5|-dqgfG~zBxn}chPX~1$<6cqjPqO1_GL<x
zdQ9Rg%Ocfht;#fQQ@U=+WU%GnLnT0kB6KFzQ9+IJa(A$nHOx2wZP$}-hYz7uRT7Fu
z<yd#F)Xxll{oT8~;lE&S4rNrEQpj^UpgO*U*{mPC+x)yBmfq3;SpmJpdp$)V3uXW6
z3jA=MUeIg7F6q(kmQ$FXS#2#sp>uzurMocPS;$$>6Gx2(R$4$<?J|{)3ai7mg4Ipn
zJM!_|k&2TvU$MeL?`|+Vl<m^gqKBdT+CnZ>lb<)UuFg;@iFcffPlTH#nq(UcI$}^W
zyhKluznAA&HhMxh>1<r*rpWl{yo<UhtSTC~VdvZF@@y-?XGv?l4!;bCEJyEKFGmPf
zDb~IRY|=Q?h<CLacz2~{J-0x<|6TQvO^AVYt^#e4a;Z;)_+acy3NInWu?tOu<Xt;M
z3@<VQPRIEb_XOkjgDRD<J1AHWAI;qEEf)&w<uWACQGqd1CF*xib4E;7_VqN1gpgFk
zAGS^inPvrfneF~JWjfi~OAF&*yAAIvV<u$D20S^uPk{s-<W$ppvh;YDCz5GB*oyYe
z%Bt$$6Re;>!&ixUVogK-Z7l{xRcr8(xFQ;-yR*a$a|E}W)^{Ut+o3A8V_6xtZtThS
zg9Qv&oE46ZmL{p2b#INn*L{Y~DpWr-g6W1!)@5COb3eh=J@_ByS@O1g)0(q|-G+!j
zzi8|$h;8F+-sZ@p?JvMKKlQ$0dWmvKw_x0BGsh%>cGqz(7*zcvwrYc%sQJ?3%y~Cc
zTt_v<X5d&IAg8_sO$`B6M5v-}ItrFMm{@}oN9)hSYQKhuTAIrOpu#1VRSU}2qY8lj
zpA@tUHgvv+6_iZ#GF;WoGyUoJ<%%dZqUml%O(gnx6YU6DWU&g9)+RG>E6Q0u-;!-R
z{XIwtdcAsh_9w@rV;5oKvE(`ta#=$WV~_6@x_}vVzd2B#HG}8Q%f5ToI$HKjRUOt|
zuyC#JWRI{ET&3>t#bS{3ErP>$9c7A3Us5m|OrT@rI&swwp-o7wuygKD{?4Z-qFvfF
zM_ap%Mpowmx%aL=Jq!C}%KwXO82f=vuzBniL_nt%%G;y4c%aQ62fI5o*{m=iJ*ZmM
zt~bs$dmwI|OP7w(GCkjdInf?HS|2O3%D<i%y0gJ_5aCm+iPX8OR9AU(xN+2IQ+(Ux
z)b_Uw#&f+rVW@4nv@nRaMvpf|pBMuX%Yq{OP$Qu8<Tq^+0_8QJq0}hmPtxEYdpfc#
zP`fB+G$#>!55%B+pqMJb@kG)m_2$hHb(%8RZYR^6WAB8(XFaCq*DhD)z2^W6PiH)}
z+*Meg!e#VlPHFDj^zhHoP4&*usQp}-p`M2+O^DIl)>tqvX<RA;bJXV@S>lO{{H6@x
z8(v*2Xm?A?>K^jYF`DB3U6vLfl1LaWkYI$~2F&qqjUcB6+J2`g2&W!TigM1}hDa3k
zci9RI90Vw3&cYqp?)t|tru*Vs-xIX$`!*y|@(=okKx04w#~McI=wKhTE8<)RvR^|*
zF~ij^VY`X&@vu(!sO4Ei0I6_IBlBjBdUdtT+xVg8F`E_IHXWa@gwUW$C(WV~s5H0f
zW3<a5;1*TAz51`%cw3>1_TsQsnQcUz#9hG41b3^Djoys)xQnp0BEE&LH{x_T3Q7!a
zOBM>gl(+H~<rl?1mm704vqZZEG9z^3dJMJ9i7n-(Tr$0Ow;MsK%T<ZX!U$*)-rmvH
z<2kil>JMc5UwzU%NL<JAQ_amNAWK+{aRumSw@u<&t$1eGuCD0bm)p~#cr8A~KapKG
zYBpX$DxvYEUNSt<2!Cz=8tTOtPL@hIC)d&Y%3PXV{`GO1{mAt_kr&elVU;OvbR>+t
zadCS8W9=lhMPiU+lsKRDetrKW0irj=b<RwUayynK-Vx>hayhAJkdk4r209K>m$6(P
zuUPNMQH?~Rj2_-Xma523*Ch^bbQtgbk{a#p5Tutpifnl@T5C}3JwktRh)j?@f0KP~
zKo{KAB8r8-z;*OH^^QqPFGco(*73)?eIc356;$^|kjHwKwp9s8lu40uy>kxyd6V6(
z!GE<_AMb!#*P(a<gXhy)KWK(rYTX;h6Uv#smXt?(TctRyts_c&qO+6tV_;OAr!b!>
z#wGX}Fpc%rs$Hn(dur!qfDsk8d4nEd%h_lD)%Kh`#4M_4(?4bV?IP3l0Fd?-w?o>X
zu_yGgGveta6z{}_Z<!ZA84g*Uvy!d#i>zYv7w5SF0qDRLUkvjFOOQ5-!2zZ-`<BcZ
zVHtNQLjpz20{()saV#p8VT>ADz_1)zd$B&3Rm9Kn&k(+e-UYd?=&Wa=$`n$rB@6V5
zpA2-hAA1~esa&%64kOPb1zt6eY&oY7!=9&(B_2AA7CzwL3ojG5Bz8{6U4XiBKLW%u
zrd18YXRHX-_J_-2FHqH>&7Nchwla@(fQUDutXlz{wnmJ=As!q3ybQ3^Dkf#S-pMK!
zi$zV4a36>E>j9t0n*(FQG;aUd2HvxQ9;%_UpYMu8TtpH={nXaT{orV5udB|~xIr1l
zNnO37^86d{4Fz3RWGd$?*1f}Y8l@c>c87;q<z{>0vN4i<IY?>uiKkPb5MtW??3;C&
zl|gK}a*)JapzQWNS6hmet(Gvo3iScBs_$;tv26|E=rVWtJt4$HDZxKJlF?(*tyy&L
zMu&sne47>tx)f%vFoMd`Ww@rXr!%tMSFm4Md+~14`%tY#83kR&9`9u_D?GF+Zk6AH
zJLB7U=O#FgK-;Mi9!((W(SweCT?IdkR#|kZyy;9NhMaeVcD(F_@$Yk|kFZvr&oS~k
z8Jl*c%lyKp_a8oY*QHVn>*CEQ!33kBU<HyD2e-Er;jDKA0K2Z%Y=+`jPeVN0XpyB&
zVK+|#8sl|!tx|ZO3A4a)^Rp&)y>f$NPhy16(uEv2Ao*;?&=1w*Hz7j9@*n74jE;W#
z$>mnv#lxDx1MVql28kkWBR$85Ds#3j^7zq`j8L+=2!D$*UM7=!FtW}*c}g*~flJz!
zm(3}|H?dPOZS4Ae^-OYQbT(<s;$8?(<`<61(^?UZ$5Jb?qcmdcqCa`c-{Kq$nX<aZ
zQH>8GgL>S(^?_tWLi2M=rKhLczCPQ8#;PYL&&NMKC1AU{H|}a`e$WtZ5jSZ-tyxDW
zMWnrn+uZau59X7`l;{X3zDm5RpW&O$;ymRzaEya=s;|f=B%|}3)%+E!jY(;;jdE3x
z2)dV7hgLS0TxkQf7go$+ctlrXss&_>&4xh_ptAfh((KsTlq-icEa{_#8+J{M-)xa!
z|D-OKA<mHu5}XZg+V2cq{r-1~Ya@x2TXKxKuc)6*4Rj~OtB9e+6Jm!KatZ&cd3r%j
z<EWe9$Kg091%kZLSk-P+G2vrn*(w$^=9*HS3oj#CCj#Gp5kp~M3#VQa{+%4e=`ZbM
z1jQVR1`KDF<ri^bf^^R+q{J%A{=m5vA}HUrZt0#}!PgE0>4}MnJ4Jsd0$Gj#kXL0G
z5g4(88zz#CCkg5RiwP#D?W^9b(1)~JHq@|0A-|_C+BMD#3ZE22ejgKB?2FS-rA{Xy
z{bLDv7ht#8xb)|5qL>?cEG556!b;akLGdW2cJ7iQ^#x}#^gj?@(XF&IYFt<ZW}e0H
zIQc_NFBJ0AlAqs;Dp!mD1~MWR5*sr%BH-NBe`DTnWvZu*Pd%?^+tf|5&B80S4PvL)
zdSYrE=3KH~+wE(pU2hTyU*u(+vawf`bJf;~47E2mU>e`8kzu93Fbv_QZr*6+1i3K7
z^p8!^yY%{V`=Zyf;X$c-HJ15%c+5!txnzkq8H($cT+8of8AFzc4a!$W3tEAZ?|uuT
zMNuXDcv#+I-EM1-Sglcvqk?PZFT|%z7v$FVRi|X5H8m+$3Na%6l$e^Sj!sB_+Y_oM
zU<t@OcJNu<I9d_DZ<-{x6xBMo2lYa3e!`xOz-fqoedFIa9ibi@!M4)JgM=VYj28bv
z1vRF;6J4B>$ZF*fHsRN$mB#J6%A8O$qIBK-e5-j<tG#)Sd{oOwsqD3Fe){o6oKsTl
zU>SU97`FA}Ejh|go}KxY<2##^Vq8r;og{D#>7|4hne^ZkTrcheec*Qn3mPQkx(r<@
zQvx}cFwT+D;o}!@u2nx&b=*|xkF()4bw*uQcK8LGq*4`mRIjcWtLK3unW!%K>w>x+
zT2;jiFS*xU{|5VHSAppEIY%|eXKEvhw!bkKaUXE?Z)?riuT}d58`8^pX2*Nj7-1|=
zAgblRJReJ-0UcBv=^jPuKxzzUK>@F2e#Lrsl{Tr;usqSs$s$IH#%j!{ACMmQo5k@H
zO-cxdVoR|JmCH}Y&-AH{)fL*KzS>HXZzeroBn{L0&J!UbTQ1c)YogE3G<WJ=86+Lc
zLAGCdaNH*!d+=J+b6P%Ch2q&AP)G;EQ|rtw)G`w85}>oIR1^8)>f<i#o>A@ze~MY?
ze5L;U7NuRSAUtL0>2KGfce{L!>Ve~Q-EltDFv_eZp37>J^05$9(#WSzRO<m}WHcTg
zb6~V2S;!BB4O}4`Hmf|n_S%7&^_qBv_Fr~-$8DUNJP81`dj|_%ZefA2p#XsGPXNF_
zb_U>LVeW2bVdP+L;pWQZ<Xxh!;JD6$+;ytu#1yp_dsMv2n2y%3*l~u26-c`7#;}ki
zAfq_{_N?NPO|_Xs<iRE_BqZcN;z^QL)Hp7CLUkgD?D&g4leGfYOmIQ+D_eZp57&0*
zb}*vO_gE0r>#8jDur&vDH^#55&G?wj-euHN@CdaVbl1X&HhuN2(HZKPh(3JBy0EZ)
zg&UIKiQ9vPa$4OM3&48Db$iuC1EQ&f{L)P(Y(Toc1rLG}sBOhB;f+b}l|3iJ&so`l
zO&UaKUSq<rwuleKf0GZbVfRZ882}6f-~Hh;gMfN*7`V)0vlEymp*+7x?9j5$6uXw;
zhsMy;ZWje{RKWeU-q~XOidmhcYkj%3E<t6O^nxknPmW|aO#I8{Psy=pe5DvtX{al#
zEmIy=&W~VyMFbI!<^m?!HmE2SF&i@t?WJ<nwX1{&F^$-OepuA_aDjX_9?W~e9+{|1
z^nH-Z@2VAwrlOiI4!&ck$NpyymhcG{QKp+edInnNCz+o&2pv@L3mJC2<>HpK6L=TY
z(Gq;V(F>HS=;-s#3`UHVxWnch4Jrn*H;5rapHQc+O2qcfa+a-zj_(pz_aDbLf7?;>
zGDTgC?nMqnGrqV+Q@USgpGGrsWZ*FLKRXln;3}#Z6K-+Y^B6<IZ@nms32s9zOvSJy
zm@>^Axob;F@bS^Woo0Ba;aL_VJKEV1b|TKt7XEaWGPuIrTB8W0gQ6LC2^JKoU0Vx)
z1pvk{0RW_bCFJI4WZ`3J<M^NQ(cX4k=R)@fe<%Y9K=SKnlq#eZ3lB!!{O<|0Qi>`J
zBBZw$9edXudx>)$c@F$kj!0V@4TqskG=As`Sl2(XRg+a0IFUD{bJXe7E4?QR;bC>;
z?gIr9I3pd=3N2eJ%{i<8RGU{;H9aS5F;HiPhD-8)wZ))~bbTd@+DhC;gS;+JNM|ut
zc=Z=Dg67+Wq|;4P>$8kVAw(70>{8+C7e*L2S0Wb|6k;R7nYIsgG1{D=PrN->l8<Be
zo76ze!I)pE`4rA&$GMxyUb97MSrYOd=>NW>L=;#WXCwEoOp%2|U|CSgkBCNY9eZU|
zIbYHrrOx%PXmT;+^b)sylnl!IRH>O?6|=lCK&zwzTvY6%vh|LCGwcyG8dGQ)kB&)4
z)%8R_HQJ>?95o5K<w@5@*e0j=nFNL{(SarWHwrG3E)*<*qLQk2D>Io?BOxjlh%4Kx
zJvpfJckQMnX@x3+ao#d<-7dP#aGBf&C~|Pq_^%BO-(g|q!CqK9=1K-StXqU!mg#%_
zAbjJ}S?lTi_a*_V07kR~Af<m-Zz8%2Kf(U+^%Vzk?wb=0WkCoxWp)o%R2F%W{dnei
zwxvfl9?=2`fg8#);&;b~KToyYR!cSUq%0D^eFl7>A~;|X=C?0Gn`^+Y%<asxG^TF;
ze+>OV3hFH%zaQ#ZyiLq_Oti}mx$`tT=<s8Y{Ye~qSoc;?y8nHpKixKuh5>JdvzcVs
zZJT9Mjv+BfhPia|S5<juLdz~;@;8Z<9UB1wG?+dQenXd~yW~bQ`b^^(3AMk4Us*h~
zSOK&R7zy$OKqzln&NDR;_J}AMYo$HghscmjG4n%O=&8;yQcNLM-};8Hw<%sQ*bd3$
z-(*aIV<AO~&G@_+0&9uLjqc!FH6RX1<a@kigUA_r?tx(!wSAb_H?wHuwSz<ye8i-v
z4Rg|Pl06TTf#m=T80kuyW}fid1K!-Bon7F5x&hDH#_gAk{vofi?@^0Z20bc?N($lm
zyhBJxBPKu`)@D*Q&C>V$tyN!$=RMX{*CAm;52kMh{(P<3JS&?_zZXUW{sAn=pm&>a
z3u-^5>iQzdZavjl>eN(jl$EONhl_sp9C}(`uG~=8Rb|Go7onO#FEuK+lgMvc`ET2U
zeO^AI3xjj8V+PkrPqaWIZqZn@gxspyWx6#+07<e>E3L0Fy`S$8(=+X_nJ#|o?}?{s
z<5-5o<nw0PI8zXsN`4u7PrFdiF$DX3L-I9@{R^UB{iz7NCt@>~N0y03bG(l^Co;Ln
zgOp?cuSkrJKWo@J36NgO`sSze04)|_;gy1^X(fM@!>hMD&Z);Q>1{Yn8nt~bWL;cb
zgBlH*dLN7$bZmFXO)x1d490CN^Wt9ApX#t0erIYjp-5)3{t8DDsOA#@xr24`y#v?N
z5`-sK+~DzE=p}`?<u7>E20D=!Oc;3bcy~yzCe4TN3HBmSyTnyMkuK&K$L_H+`v@%|
zD0}S;c_NA<oJ;+w+E;hBkf~bhjIJEx=pn$kZ2R7F=SSxNbfZHF^gOze4#Hr|?-3(7
zy3jWa4_)4lfjQM%T&_@1NuvzxIO14TmZdWO1Mwh5===lT0|)q1t3ygUY>=21S#N3|
z|5Ya`Rt7%LJ-qp%CXzmVPSJip+>7eNSp0#fDNGC4o9CfQ7;huRc56nGzQ2`b#=;jp
z;2~6QbCymS4r8W2n(|IuDI6c0j8|{^jVpbT$plG5DE+51Ys@U@HR5pnexmbmf~&Xs
z-1-Ek>04VgX;YvJf9%npV`i5EPv5W4Bva9+Q})bwU64&<ReR;gC&7utyX86OnQ@Y9
zUEksF&O3{ib=@EBkWW8NA@m~L**D7$FiZ#DPNY>>mvG-`S}^yJ2%cgYUI+tLtB0=@
zzRH0AOpetCCuhNd8lrp~tP7CgnT6}DUdv*IWE=v$e%_c<GBAZ4GO<RSIUpi&OC9yk
zo^oGu?*AosMzJ8wk99HGiwbI;<Sx?o13lc+t-t<0Q`wf4RcR&ipR=aFyO;1_u&caI
z<BU4Rbe`@kw1+9Uyb|7)`$=VOGe~SU*tJ!Kf-X`2vE_0p<lqRQyynS%T#0!v7EzJk
z!X#d;W#&AHZT|U+8oaHe`f*_O*-8?V=CDyKcMmsRfV;I(K#IoOp~h;-=pO;#zz!56
zY;AD?J^@kwo9qnkbORbgVVflPnv#Nr-2Xuqr*ra-3ne(t9f(Z&|0z)a(Chz<^T6?%
ze^EI95u|vbBq#nw5rJtvfSgbkOaIwVCRsrO2YCF;_;i|_bN&T9@}Kpm8}OY6@Ea8J
zTXNky5m?g`NC&h72YLdffqCFDPvAG;2N=~0C;|O~3k2um0|~*-UO+nNH6kDwL;@rO
zcmMN%FW?g|AQvzkOy&&~1SWuuyn(!uN31}=|AAFyIUp1k^#9h@{--AH|K5RG0r=lq
zI&gtEkP*sS2|VKs#3TKm(EI<!MyP<{eSq8)tlIzY$uR$r9snTvHy7Xl04G=AKM(PL
E0cioVFaQ7m

delta 7509
zcmZvhbxa(<x9@kcmM-qD1s2!hi<aW<wm7A@J8a8BvEs70yIXO$;_hCoIK?RxxV^vJ
zo0pgOM*cWw=A7gtlYD08`|<~M`x><dB8IgWtqFl<JCpLh1pr#R5-KRk>aD~)X<#x%
ze#fu%Z)R07yu$ruy)-!$8hNA)x*AO(i(qYikSEQ<L^WE=QDYAwyEeLL@Od9l3l5#*
z=gRNDO_Si1o)igJJCK^CB8lb|m9<wOGL{zWPDB=}^L|mRsoNi3S@DcuJG;1K)3-;D
zvg9@u(fV~{E+_wrZf;jl#QkwM09Fsf`t4mqnQ2ooxIsL)z(v~A&|1v@P{9Y2tsZ4|
zdFC#?qWUP8n-ir{y=j}3Ou@AdeW9F+a8X-<P@(N@T1}uB0>CF8e1jhI>&l~o&3DD_
zIBJ{d?sr1xX<57<@DG=z1L?aI=!tlb>YNu(7GFmm@Qrd;6L4*?LJM|OV49iN9NeY{
z2F1i8OONZ5l?m|@<0?2MMnu6b02gv>?G$%at86i~?bgR-&;>3X{y~vl;xh6Q+76cY
zhL*7Gw4Pugj?J5b78!h{b3UW5?;sD1=J>9XqrW8P^8M1zb5uUNWz>$}@AIJe{+1hv
z6sVtxgwVKmRTaIDbbYvcfNjZTHj+rSoKBi0jYwogtTn%tN}TChd{rGI+aHw<B!3hu
zpVNMXW*@PiSMuT{L+84xpQczknt#=|Ja6YyR#(mQRgG~Dm?|GZZj;o)E$pqo?vE*o
zI@lYOW$}lOYZL(p@!Gl5LfQ56x_I#=6PXC|B4-^Ql-hGZ^&_8>VRpvBq#jzfNqy6S
zEruJq)L+i#iZ;kppTWeV<>1?k54iV<cWsY*uxhTJj8U`K^jf&^4aYwP&}4VSp;soR
z@5{_0f%)bQf`H6FKwbPBy3a+cPIe!)R!Ylr@q1{jWj*xXub0fo4?V!o%HhHV*`hH}
zuXJXSe5Ukc!{k{AY_D2`Kv<SO=-L4!k?6git0w5VjM0`%-Fa1fSUx`Dk_sLX0&(Vd
zuNhg!Szx{7^zdHqjb3tqtsg`A-xiC|`qbdnkw>>=#@17!*J{*}$CHaN{A$xMXRz{O
z)8RY!lC>C^-H4V|ERi>x%k&z9Tcf@y5dv)Xa|U*|z&o>G(VgLNEEK@(8T0m+?(brK
zEH#5tv8ShKE;SMOPVa8+7PK!S;_Y8&8Fzt|r=W+FARY&u>q4+Od^*8Nh?}2vOspfb
zL(tc-oO;4!E@tX$cSD;=b(P)U@1Ny(KlEBXCk)WD3wsRxPBGe=I%B34ZCsIoGgQk3
z?X9*Nu`f=kOUj#*<`ylH8Z4kiiGP7Z;5S4vc7}STT!2fOGy>V4J><Q&)}FV#HI|~f
zl@x%ipYUAYpC-OxAOI4EmTB=1i?XAgiy4s44tM(E?>4j-Sbd9mw=(W~cR2_O)@v6Y
zN&ui94gk<6R8tZoxV(Yb`872dFmrkPRojg!s6{`6(!GT(ZkqI}O(lvLFgi#hQSzWh
zuK$)f)lw`c-XvVCKm$J?T#1|6O${WZcSF7nZL?k8cwIYUeLsO{e}Kn9Nyc@jwHSX-
zyai-T@&An#=*^g~_dMQaa}s}2+E31fAAP{Gx&NriFcs6tC}CI@V|{-C%SP0;w$GEb
z@G}lkuThn<M>+(DN)Lm^N3^Oq*wOJM+w$?w?cZl8XkkQJhB?^La6hPie9s0>rcfYd
zV8v9_-Ogs$`bvZi9I}Kmb;Tw3XiN)$s!E;&Br>9g+a;B!OsfxM6{yO~1HV-yJ2$Jv
z@O>+b&!c}^nvn^`ns{7<ZI5-MJh*QcZ}APthx^|h<-c4K_NF&5vr8B}xxs##HF2bg
z7Q_!zAN%pKfu-Ky!f6R1Wzyyz9{l3D+!87(Q7RWV&Pcntcm+bhK1Zc)z|ZU_b;%Oy
zj}5{;5$^!De+$XLdUVPLxD5r8FPek^0}FqYxu?T$^-30p26hA3mnHjAPT0ntQaN;D
z0rfMoI=s+pX~l6zDJI#{XvOTZPWIE3DSh8gbx7h#uMfq=&s<k4!T14C<EnENAHiI$
zvVfLp!U((LUS}UnOn1oFoT!*lvrS6MFkX}!Tqn)Pw!|}^?>u>Ao)l1&G<)83{g??A
zcM*D@x-SoMQC)$xolH_a&6=rZb|@RTKZF|er}wY+lex;8Y+Mf3E(ARLk3M%RmWLw5
zC4wS<rCYixSh<c^{?Ax3Q++)71fxgKd;+2e>`VPoGE^}f3~^_AezUypNJk3TaaXGf
zn+O@lJW7r?mVVXK><Ug=HPX$Qbi%nOCA=razlZ&sr1;JE^^-3rNdk4Qp4|X!zVY@t
z+kJRms&0SH#G05JGe}5Yd4kZToRZeZYUj5~$=TJfwUQQBC##mgjaa=2b4PGfkgdvz
zqlLK;W;BcMOr<5PDLw4gzrzN!Ow>DGXc%_qskz_6zR8To@jhcT1>J{zz%UoJXf@*&
z_Djky<bl6P&wX_@NG114>M6akHCDf>TuR?7#Dhn^qEIf~13@0Uk?qD{KFAc25JHUh
za#(LPLvrPdiI$7P;M3nf^0m3?!`<tMr`T?t6^y(&)w&f6=?FMgRH4>o$#~QC(cMbu
zs+Q{T9_Q{_YHdURl^GqKHf!pj`qvg7*aNX@|3DnSlLX!@U9nYa*33P*O_<<X3hN~^
z%U+!vPnlt-J)g``{9dR2+02_Czfv?ku9VfI@uk2`sUO%??~jbG{>(H)S6OWGGi598
zPInYc95TMv9V;s$In!*fWQ-(>iiUrMcYmOe(@BKWK{LtI;+L$@o!N#2f-YwBU`Hz7
z%vnh|JuCQJ-gsZ~HBoy|Bf1w^3&9h@wR3fM((JnV-j!nCaf&M0SZ_^U=%Jl-aO`3y
zr!Cl8!*%j0Dl^xD6aFl(`|Xn{T)|iz6tOBR=<(~BKCaLFvn{C|M*9A2i%og%9gQ|g
zF9-O;I#Ou?q4+?qxUr<5JV(4dEY^osD0#Bo*6sZrhoB$4OxByK<U!Q|w=&*-kt1K4
z5gPljOJ-YqxJzjlrLd9K;1(9GuwipPq+e#Q(&Jt=Qk-_g&{p7Clg_DEtMB@BDsuUC
zF6@@f^e#M|Rh{#U80xfF)X5G41vX~lz0H?gJE9f{!<?^Iku)ql5Mljyc;7kHDmXTC
zkAB6%rVzHV?HF_oa7iqa^hOSOf}JX7R{G8I+s=e`n$zs4GZO}M-X@Rg8=i#EqS+Pp
zW|0{M?du&BTUTZO{Y1o+EOUCvK<-5)ngW8yu?^y(N#+cmfrqb}XWk69wZ<_xhiDqI
zClbt)?WT5f`|th$RMt|%I5-v$>UgBKF_q?9y&GOdHVpcUkU-penk47SFCFv9)|D*{
zUBDn2<)hUKeLX7y`xmxN{~BEftzqbV&WOulXEepV@-wiaBl}OWe)Q>7QbSsg7ig+G
zkkBV#s!8U~-!7Ix=)0stpz{j~6bCcD@;zh$lwWY)7$Op>WfBLwE2qf_k1nLwk0#zT
zs!PJTmi$J3ENlX<WZSY0mBumFNzlJXk>S4j)A)Rm=;T^^TFIQ#Z8rT~>)^1_+~za2
zB)$&VYaz|d=FklbBzn@I_`sEonRmE2>{Kaj8Cvyso80c*SSl@^jvVp2%`cC(CcHCY
zZ{FA9kdDsmq2Li1Ml2{gHL5i(-4i}FW>FL4C+V8$?}I-=(z{O0)bTw3=f=2jj+@A1
z*@*QVBTrwoahQV|H4=C8KyuKiE|R{;V^N2@Kp`^ffr5zkjinWfp;7H_+5{0!GIJ&0
z$T2bB>k3M{0O*OZDD9V}3DqH&kPZdaN}jG)Bid$IHQF99p@r!DFU@YQ4*1*}wR5(P
zh&m;AA>6d$H$}k?4Gj(hu8#DsC-fTp)t;7Sf_JP!c^l>PNYxAY6=m&0cYo~34n^51
zIG#WCtb1}4aG!lIEuo1x)=@!4ga#wYEDz=S5Yki&3r1KX$%}rjI3%`cwRG`KR~WK(
zR2w%}eucyZvm<8V{)E1p-Wmrj^arN*@eb1%HZzirH;a6Q-ffpuH)B~pF|cD`aRv|L
zwbMDzf(+kd&ddw$f9j@lus4$ItBJDyS@@2jvUeMkOs1!^P??~%LH_ZKQQh(;o@lhv
zJ9=SRA~RK=wMH)Y68*475VI2c+Pp_wwU-ESA>^F|qsd?!;mY8(8z;hI&?@E=h|L>L
z(<lp;2?A(E*)-2Mopj_XLJmc<ANsr?XHv04FA?@MBog4L#^6t`wZFSFg51GU+4d~V
zpmv`Rj!W1NJ~K{rq#|@QLovrEicPHAc6J&sbuEOgW6OO_5m+X1sxAcQ@;3>kOs6o{
zYVs{uJF>PTYv&YCRSAX8&Bc5vS)K(h#7yjJyi#?Ah;OBC4LeQqL#j?zFR!gM@=nO9
zwe+Am6@e&vXhIFe)6$aOTXK&fVidiK{zS%l0fNepLpmIp{;I=%p3`mMZGf7J^Cf%h
zbz<?bxO93o7v4k=rf`Crg9AKOa^7{$B)bZ>x^r`MF_6>b#qA9X47k2f_`wX}(G6t7
zh`}lQJx$<rXw3~Y!Q8YvR&Q8t06XG&hH@MRoKTuzb3AM&sAkHSL~t>NyWiW>tN*)C
zh*M&o4L#<UU$MIS&Wd-6bsabj9Mh@PYs6t057Ubte%<L<E{)`=S0p!%Wf?y474m@j
zHEFWAe6hS(O8*H?!q*(#ab9o`O38QSo>X?W<4NsNg8ca&$%pMTu%00acQeSU?@+X=
zS>z9+gRksUwBv8=*+yy%Q|>uTb}93{MTc%745u9Z@Q*95Ifn*Ht<!ZMc=G(dsw{A(
zCxY_V-j7Bl>LZ0p_>maw+K5bUtI)t4_Xts59g9TF;d$xCg{3wz{mK?heJx8c7ozT{
zM!l7q_d=xT9_f@BRi!93U1Lavq+a=<nO&yRIM`XqrQtaAa&*zyt%{<XXU|+9)sCq7
zcNQ_XNe7NJMERF6lf~udFJ&LJ_micEOOxqU(o#Os;H0YteOtjO&n@?-a~xEJeF@mP
zFNj5T4UQ9rlf9{!&h3wMIRgInKBPq?j_zZqdkF<;hYF?mmV3HrcS<MW^$)gU&}k!N
z_HJMqgS9gsN;fd|U+4)YZ=gkP>(BR&PnB~V<%LS;cp1LA-p`TfSRKRwqVKBGKj#?Q
zMijG=)L9mvGYOo&3p%WUzazbYy~4x+t4fA*XrY!z#QmU<wzd*x?!NvOcPWi`RuP9}
zUG@>ahQv^~_M_Ou5%uzqc}ga^bl1s94{KEn_UbsUIuGr_GI>j<U<;Qw$EZ)sv6d7i
z{{p&gt9XuhuO#|V3vWo?#_4i=FnfXd9aQPNBIhNUeJ<vWyeM^OlR%<2j7>M&=_-JC
zRvM$t=#qgT!f_n`82n-11KDl9=b1CTe8TRIJSwROp@fyO{3^7-<|C;PzTn9(^M*{7
z0lIw84W5hhr(X;=il&;849`VMdx}x*rI?0(;xV4bCP&|C|D1Q=2KZkU4nHpplY2bJ
z-{A_#rls!xZNr<6zQty4f%zkh%1c$Zb<npDuJql$6wZp%62tLbeUyElHnvVwvg4Rp
z>(j@teSdxIcs7X<+qe$1N-t?NQ=@PBU3ZAQ8da5j`EgbpX!_xpMR=BZZvrIkfIHOs
z7<#6%gJPe4G(sR5<HP2h`YqGb8uT6XG($^^B%e4K6h2$}J6oa71P1N{d0-<|?^VXw
z;oY1RpDfwNHe=7!vEqeck&>*7Di|n7Gh0lT$VkYu0%3e^T|Rhv<_VcwXJM#!bp&@0
zY2S8h+VWiKq8w73SjTvo%O5SrTwu}xX=KLzm;6BE7LD&xwPm=Taa+-!HuBc=*sLdN
zEn8B6&5lm2GzEPhVFoKJ?g3Bfs^@;lN(#ot)W&yJLpgVd^o@CsbIZPCzeWWKatC$1
z>AGyNGPJpogdn(<FY-rtS6Frg#DZ;wu4^{R`j2WlAvAc|=0e1F<UExSV>kF}l-E{r
z=T3N)Bz@HvivcH+EWHB2^g6vZi3#2Q<bLv*U=Am^o;wW>wlN-$beLffIH-VR#CS2I
zT1PzCj3#u|7xai`74`2aP(2|}E)?JIS2~ceNh?*M7@2A_Z}Itc-9^fZ8#2reOLc2G
zh=r|hnofo1H}*NLwF1-@DOX*dF{xD{t&vCN8p8cmXo|X7S|V?|m?sdqt7XZ2DX}SJ
zTJ-gA?<)d#*qrsYY9!=FPVR>9>9YzSFyeX1K#`K_=Z&E0KwlqR;I$|ZV4vc15_K7%
zL;|d}u?H9F8bYc0wv2j4F6Lc7wc!dK?BSXz=IkPtX&Vf9AERL=%ZX9g8_@APoOmJ>
z)Db4nvDFn~S`e0-gMP-IqOqwVJ?jwr1C6^;m_BC-rgmG9)PWKh##^DJD?~zN6F~bd
zjZ%vHj!|t?KHx=t-d^KxFCnTrc1`j!3g}@em-x~Up_lk$<W?`$kzL`tGKS{6p46mr
zf5h?*nH`^li`nv`gF#l`Q3lbk2L9qgp^A$BUj8}`gq6!;K?<Q4XD@kz!PzQV3nIwO
z>K9=P%v5|05~fRv0c37p))W^JvB|#<js3x8@gp1AVJiFH@@w7)Nwtz!k|#(-n}O1a
zXr-~RN4xVBr!}K~dtIg<@$ePOX$-uvEOIt67CcgdTCO>!UQV2#Aq8cYu&5-TtzLF2
z-p{J;0v+~eZYSBkE7rL6IpvoRraC^x$ez1!7)xH_Gv%3wCBD^e1TP1MjuhIy<&1l1
zm%hRa*4-_3VcQ_)3WlDdrtD?QTaSd~uTY=$z{S6%x;{_+kEKcq^QPZd{z9<w$mnw~
zj}_fGrO!K_<z@JbQV~9B???riWEW)#^8u$JPjlpGv6oLj^yoUiT#9+Zw=vW2q)ujw
zVd6eBM29#@-@Z+eo$5T!8~@R+vLXc;KnJ|$*E@S^M&i}m7VPc0$|gtOBk=Ah^?FR2
z44G-G@6u0&{z^h7=EVsw(=5xa4i@{T&vOye<+kZR&sh)nPjxyZ&pir7CS=<hIK%y*
zA<R<&0oQq<fTj|GBtM!ZW+G2XiRV>e*bKYxe1pPU<ctA?d``7oO(}6GOlOPq&yke}
zE1nJ+li^Gap3r)7SklPyf>wp4SOru5Gmi-tJ*VQ&TR9yXr<Mh_Q<{J$d|#IKv6yu`
z$#+<$k*)iFrWl50DuzMxQW%YwxCN3Z2^BU+i^-TSm|>|U3s{|z<mN7(@fRlBS~9lK
zz=W_s+46-V>W|eT;{2<Xe~vw0iGZ5#*F7zkcM;5VlPP>7u;cC-eafn6-8Oi2^`7+L
zrK%zQbS3~UaT2@3J%8if{yWG<)cP6NWheLosy`Y(`rJgM4x~Y1+K>7zdXpFV5{L!>
zjNt<SAOHa1?qp)+YyHWI!`Y`$Z{2B?9|wjo)&!El-mjj#Ri`<eyEW+&36x8uF9s9Q
z#QNN|VQ)~zFm*g0-sbE&^sE{#bjF?)+~q~lZMaL>YD5%$DKsoiPzX_Y=cGr_Lw8)V
zk?^XH@!}N9ccPP2k&D@pq(o_3`Mmosa@hK3?zB$CR}`gf0$5m{*X|vJ_2nJWC&Re3
zCevUinEd#P5KE5hU3_flD_1XST&d@2t@==_{v+nGN#NI}){j-^{ODm|>-$U;s&MRx
zvn5D4F~qh-zuyvsWxbta5lk0J*^;%%px&^jaTw$J%J5u9hrFOKnQAEd)7sivA70;(
z!gMEBsT=!cyuvs^@u_Hd1wK`Je5z&P()Yh(FcztKHU@=JAr?8F`T4CLwseN$^!Y({
z{NMa_aaO<BiN%Vnab*Zs<DKH|lF5TyFnD6~29x>OdgMrSic%(;7F>SxHTYNQmloP+
z$5>>=?^I_jtY__?=6<P!?_&h&ue|#bbf_pAKeaTQTS0n`YFIyZ@e_q>7TjQ==46nC
z2kTDu8v>vOcJw6S)RD1Y?s)r<aZU!d5XP(Z8sd?0g={NXHXQ32j(@UV4F4(ll-r%V
zXJB?ZE(UqIxf>g~p+d=%`<%f*-+hQz(+bGo1bSgj>0TKQIv!rO-9pV6m4wHRPp*Y(
z#EZB&cqc|fuKFR9DK9aTyZG7^cH?&xFw9>Tuk6YKvT97X=J*2Fr;JDWIBVa(VT7bM
zwn;1Rp?=5TOrc=u<xt((^8w>Du3yI2HwsKdXWGEu5h?$SERxtX9suMFDM;{qdi{Er
zdEQwCINr_5m<x$2UuMFyeG``oBJCe@d;Ou8=_{uWzm?eAtS>(bhGOIOKJ%QFXu>RX
zI9^ko24O5<EwRqtlQ6xd1zx)o`byb+7B!K{^0UNQpo;22s*}D4-<`ZW<awLM<w+vh
z2C1W!mehz#$JqUCvBrV%hpeUoD1Y;o*g?=(Y;aa^P4u&oB{7HnnNReXl?C4(=D<-s
z%U0P?@ffX5^m8wo>(XZ#C5NHJq*9o8i5EZhl`#bM0(yn1zZKOn_U&M?_JAa{AXglH
z`I{@eZ0j7;0qCkQfA21@L#Ys(udJ@_sx*9<ML}H57hLGPAAdpH&a%%99n{iw4(mgv
zKyltJLf;%1^!CLqiLAt=^%_vz&qwiJM|meE@t_ybp)>7>e4n!0ru-Ax@DR4c>p3lZ
z5puRBvl?cUEnMX69i|lHV$@upRU#&5&}PT_{$u9B{=$wx;!ll(NX=L8Q<&#gh8(^c
zzuP#Do2v-}ESq$&KQbzu6xr!oO4gyp^)Ku<L847J#D-`mSDCud`2EELqOKsXFa&xF
zkCIjw3Kx5eniUp4m3QSU$S}f0rI&(Ym)LUJOh5HPc97)W=sPRhA<(49Nj7`BRgRH2
zT1VNo+PfPVcClscS5U8?wBL(0DV+&w_G7##AAnYP1t<hNQ^y3Q&Dvc$lEFuzU}7bt
zjw-@5(Vn;(M^VHc{uFu((Ty_`k!T0`=#_|pwG-1)-gHl$&e_Aea}gNl1&)5xEbYve
zi9Tq6MQ!)T7X!1d-~^osj@pK~zah$IhxR-{R+v@;%Dg4Y*Hi=aA*;6}w|oecCZ~cF
zw4nF8v^{b^3GbU}2<UnH=jgf08VyV<+#}+#sZSJ_St0zL_#!nOV|}$IPd%m@w%&?{
zQBmbV>FVdT5xS&pe+FQd-Jx^y-4cK^5aL1N_%o$ih1fy7n1r=<@Xs=T{<3cDw$e98
zfjAz==6d<}xC7-kW9|22{5=$2d`TE<J?k=CE3qhcgR*A_vTbT2os3V|djlSOjUwNm
zn~q%L3#I*ksqEiA<p_KF{-HOBJhx%zXjs9x#m^m8=w=;kNAiaCb-)}{oS7%#tOkXG
z+!!6pWP`f=m*o1lykCU6eq5X+7K~n7@7?tbV<tN7Cdh!SBm}W~m@f%8n1&s2&JDA^
zg$<xiD-k0PXI0lYyJ?9T4L%E+NJq;?+JRT`DiNS#U4nL4=i$K7<HBOO{S|y6zq2JV
z?-1)gF27uTk=Vw^)GR%Ub%CWywI6bo=J~F<bKu3F>647wOgv;T-mmi`IZeFz3+n#H
zxLd{Bi>1uFrEkv#IC&}UOrmw07GxuatR8=?(7M_So(-qyCNP*(#zJdEe>m~@IPQNJ
zm1Xd8MB^{#jCe+FisS`Kk~a}~0RI3n|1UidLE!~7MFN^7)KC&79GkyEWO)HukZ4{0
z+eU180m*3oA3fn8(Ehiberc-u&(44Hf)EnR*#AboLWFq(-yz8_{3or*5(u&74P*hX
zBS1btS)_vtgq9Cb8tM1ngt2EzM5E8YB{jrPAD|pE90iClMFSEeq<#N6cMy)gKpEss
zVj!ZO6iAMk{%>a|;;%1I8Muj1@B{J#?GRu5fDoWNqRS5`gg(j*1pJ4z2r+*kJJK6v
zgo8hjgy#Rp<o{Cwiz?#ZPmBQlwD$i5|HHJJBI-YW0{|5NI^q=ofa3!EXTA9^evD@Y

diff --git a/src/det.irp.f b/src/det.irp.f
index 967cc7f..0276d4f 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -32,80 +32,70 @@ BEGIN_PROVIDER  [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
 
 END_PROVIDER
 
-
-BEGIN_PROVIDER [ integer*1, det_exc, (det_num, det_num, 2) ]
+BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ]
  implicit none
- BEGIN_DOC  
-! Degree of excitation between two determinants. Indices are alpha, beta
-! The sign is the phase factor
+ BEGIN_DOC
+! Occupation numbers of molecular orbitals
  END_DOC
 
- integer :: p
- do p=1,2
+ call get_mo_basis_mo_occ(mo_occ)
 
-  integer :: k, l
-  do l=1,det_num
-   det_exc(l,l,p) = 0
-   do k=l+1,det_num
-    det_exc(k,l,p) = 0
+END_PROVIDER
 
-    ! Excitation degree
-    integer :: i, j
-    do i=1,elec_num_2(p)-mo_closed_num
 
-     logical :: found
-     found = .False.
-     do j=1,elec_num_2(p)-mo_closed_num
-      if (det(j,l,p) == det(i,k,p)) then
-        found = .True.
-        exit
-      endif
-     enddo
-     if (.not.found) then
-       det_exc(k,l,p) += 1
-     endif
 
-    enddo
+integer function det_exc(k,l,p)
+ implicit none
+! Degree of excitation between two determinants. Indices are alpha, beta
+! The sign is the phase factor
 
-   enddo
+ integer :: k,l,p
+
+ integer :: i, j
+ det_exc = 0
+ do i=1,elec_num_2(p)-mo_closed_num
+
+  logical :: found
+  found = .False.
+  do j=1,elec_num_2(p)-mo_closed_num
+   if (det(j,l,p) == det(i,k,p)) then
+     found = .True.
+!    exit
+   endif
   enddo
+  if (.not.found) then
+    det_exc += 1
+  endif
 
  enddo
 
  ! Phase
 
- do l=1,det_num
-  do k=l+1,det_num
-   integer :: nperm
-   nperm = 0
-   do p=1,2
-    integer :: buffer(0:mo_num-mo_closed_num)
-    do i=1,elec_num_2(p)-mo_closed_num
-     buffer(i) = det(i,k,p)
-    enddo
-    do i=1,elec_num_2(p)-mo_closed_num
-     if (buffer(i) /= det(i,l,p)) then
-      integer :: m
-      m=elec_num_2(p)-mo_closed_num
-      do j=i+1,elec_num_2(p)-mo_closed_num
-       if (buffer(i) == det(j,l,p)) then  ! found
-        m=j
-        exit
-       endif 
-      enddo
-      buffer(0) = buffer(i)
-      buffer(i) = det(m,l,p)
-      buffer(m) = buffer(0)
-      nperm += 1
-     endif
-    enddo
-    det_exc(k,l,p) *= (1-2*mod( nperm, 2 ))
-    det_exc(l,k,p) = det_exc(k,l,p)
-   enddo
-  enddo
+ integer :: nperm
+ nperm = 0
+ integer :: buffer(0:mo_num-mo_closed_num)
+ do i=1,elec_num_2(p)-mo_closed_num
+  buffer(i) = det(i,k,p)
  enddo
+ do i=1,elec_num_2(p)-mo_closed_num
+  if (buffer(i) /= det(i,l,p)) then
+   integer :: m
+   m=elec_num_2(p)-mo_closed_num
+   do j=i+1,elec_num_2(p)-mo_closed_num
+    if (buffer(i) == det(j,l,p)) then  ! found
+     m=j
+     exit
+    endif 
+   enddo
+   buffer(0) = buffer(i)
+   buffer(i) = det(m,l,p)
+   buffer(m) = buffer(0)
+   nperm += 1
+  endif
+ enddo
+ det_exc *= (1-2*mod( nperm, 2 ))
 
-END_PROVIDER
+end
 
 subroutine get_single_excitation(k,l,m,n,p)
  implicit none
@@ -198,3 +188,70 @@ subroutine get_double_excitation(k,l,m,n,r,s,p)
 
 end
 
+BEGIN_PROVIDER [ real, ci_mo, (mo_num,mo_num,3) ]
+ implicit none
+ BEGIN_DOC  
+! Spin Density matrix in the AO basis
+ END_DOC
+ integer :: i,j,k,l,m,ispin, ik,il
+ do ispin=1,3
+
+  do j=1,mo_num
+   do i=1,mo_num
+    ci_mo(i,j,ispin) = 0.
+   enddo
+  enddo
+
+  do l=1,det_num
+   do m=1,det_num
+    real :: factor
+    factor = 2.*det_coef(l)*det_coef(m) 
+    do il=1,mo_closed_num
+     do ik=1,mo_closed_num
+      ci_mo(ik,il,ispin) += factor
+     enddo
+    enddo
+   enddo
+  enddo
+
+ enddo
+
+
+ do l=1,det_num
+  do m=1,det_num
+   factor = det_coef(l)*det_coef(m) 
+   do ispin=1,2
+    do j=mo_closed_num+1,elec_num_2(ispin)
+     ik = det(j-mo_closed_num,l,ispin)
+     do il=1,mo_closed_num
+      ci_mo(ik,il,ispin) += factor
+      ci_mo(il,ik,ispin) += factor
+     enddo
+     do i=mo_closed_num+1,elec_num_2(ispin)
+      il = det(i-mo_closed_num,m,ispin)
+      ci_mo(ik,il,ispin) += factor
+      ci_mo(il,ik,ispin) += factor
+     enddo
+    enddo
+   enddo
+
+   ispin=3
+   do j=mo_closed_num+1,elec_num_2(1)
+    ik = det(j-mo_closed_num,l,1)
+    do il=1,mo_closed_num
+     ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) 
+     ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) 
+    enddo
+    do i=mo_closed_num+1,elec_num_2(2)
+     il = det(i-mo_closed_num,m,2)
+     ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) 
+     ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) 
+    enddo
+   enddo
+
+  enddo
+ enddo
+
+END_PROVIDER
+
+
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index c3df9d1..28c9bac 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -68,6 +68,7 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
  enddo
 END_PROVIDER
 
+
  BEGIN_PROVIDER [ double precision, eplf_up_up ]
 &BEGIN_PROVIDER [ double precision, eplf_up_dn ]
  implicit none
@@ -84,40 +85,151 @@ END_PROVIDER
 
  do i=1,mo_closed_num
   do j=1,mo_closed_num
-    eplf_up_up += mo_value_p(i)* (  &
-      mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
-      mo_value_p(j)*mo_eplf_integral_matrix(j,i) )
-    eplf_up_dn += mo_value_p(i)*mo_value_p(i)* &
-        mo_eplf_integral_matrix(j,j)
+    double precision :: temp
+    temp = mo_value_prod_p(i,i)*mo_eplf_integral_matrix(j,j)
+    eplf_up_up += temp - mo_value_prod_p(j,i)*mo_eplf_integral_matrix(j,i) 
+    eplf_up_dn += temp
   enddo
  enddo
  eplf_up_up *= 2.d0
  eplf_up_dn *= 2.d0 
 
+ integer :: k,l,m
+
+ do m=1,eplf_factor_num_max
+  i=eplf_factor_indice(1,m)
+  j=eplf_factor_indice(2,m)
+  k=eplf_factor_indice(3,m)
+  l=eplf_factor_indice(4,m)
+  temp = mo_value_prod_p(i,j)*mo_eplf_integral_matrix(k,l)
+  eplf_up_up += eplf_factor_value(1,m)*temp
+  eplf_up_dn += eplf_factor_value(2,m)*temp
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, eplf_factor_num_max ]
+ implicit none
+ BEGIN_DOC  
+! Number of factors containing the Slater rules 
+ END_DOC
+
+ eplf_factor_num_max = 0
+
+ integer :: k,l
+ integer :: exc(3), nact, nact2, p, p2
+ integer :: det_exc
+ do k=1,det_num
+  do l=k,det_num
+   exc(1) = det_exc(k,l,1)
+   exc(2) = det_exc(k,l,2)
+   exc(4) = exc(1)*exc(2)
+   exc(1) = abs(exc(1))
+   exc(2) = abs(exc(2))
+   exc(3) = exc(1)+exc(2)
+
+   do p=1,2
+     p2 = 1+mod(p,2)
+     nact  = elec_num_2(p) -mo_closed_num
+     nact2 = elec_num_2(p2)-mo_closed_num
+     if ( exc(3) == 0 ) then
+      eplf_factor_num_max += 2*nact*mo_num
+     else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
+      eplf_factor_num_max += 2*mo_num
+     else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
+      eplf_factor_num_max += 2
+     else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
+      eplf_factor_num_max += 1
+     endif
+   enddo
+
+  enddo
+ enddo
+ 
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, eplf_factor_indice, (4,eplf_factor_num_max) ]
+&BEGIN_PROVIDER [ real, eplf_factor_value, (2,eplf_factor_num_max) ]
+ implicit none
+ BEGIN_DOC  
+! Compact representation of eplf factors
+ END_DOC
+
+ integer :: i,j,k,l,m
+
+ m=1
+ do i=1,mo_num
+  do j=1,mo_num
+   do k=1,mo_num
+    do l=1,mo_num
+     if ( (eplf_factor(1,l,k,j,i) /= 0.).or. &
+          (eplf_factor(2,l,k,j,i) /= 0.) ) then
+       eplf_factor_indice(1,m) = l
+       eplf_factor_indice(2,m) = k
+       eplf_factor_indice(3,m) = j
+       eplf_factor_indice(4,m) = i
+       eplf_factor_value(1,m) = eplf_factor(1,l,k,j,i)
+       eplf_factor_value(2,m) = eplf_factor(2,l,k,j,i)
+       m += 1
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+ FREE eplf_factor
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ real, eplf_factor, (2,mo_num,mo_num,mo_num,mo_num) ]
+ implicit none
+ BEGIN_DOC  
+! Factors containing the Slater rules 
+ END_DOC
+
+ integer :: i, j
+
  integer :: k,l,m,n,p,p2
  integer :: ik,il,jk,jl
- double precision :: phase,dtemp(2)
+ real :: phase
  integer :: exc(4), nact, nact2
+ real :: det_kl
+ integer :: det_exc
+
+ do m=1,2
+  do i=1,mo_num
+   do j=1,mo_num
+    do k=1,mo_num
+     do l=1,mo_num
+      eplf_factor(m,l,k,j,i) = 0.
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
 
  PROVIDE det
- PROVIDE elec_num_2
- PROVIDE mo_value_prod_p
 
  do k=1,det_num
   do l=k,det_num
 
-   exc(1) = abs(det_exc(k,l,1))
-   exc(2) = abs(det_exc(k,l,2))
+   exc(1) = det_exc(k,l,1)
+   exc(2) = det_exc(k,l,2)
+   exc(4) = exc(1)*exc(2)
+   exc(1) = abs(exc(1))
+   exc(2) = abs(exc(2))
    exc(3) = exc(1)+exc(2)
-   exc(4) = det_exc(k,l,1)*det_exc(k,l,2)
    if (exc(4) /= 0) then
     exc(4) = exc(4)/abs(exc(4))
    else
     exc(4) = 1
    endif
+   phase = dble(exc(4))
+
+   det_kl = phase*det_coef(k)*det_coef(l)
+   if (k /= l) then
+    det_kl += det_kl
+   endif
 
-   dtemp(1) = 0.d0
-   dtemp(2) = 0.d0
    do p=1,2
      p2 = 1+mod(p,2)
      nact  = elec_num_2(p) -mo_closed_num
@@ -129,30 +241,27 @@ END_PROVIDER
 
        do i=1,mo_closed_num
         ! Closed-open shell interactions
-        dtemp(1) += (  &
-           mo_value_prod_p(jl,jk)*mo_eplf_integral_matrix(i,i) - &
-           mo_value_prod_p(i,jk)*mo_eplf_integral_matrix(i,jl) )
-        dtemp(2) += mo_value_prod_p(jl,jk)*mo_eplf_integral_matrix(i,i) 
+        eplf_factor(1,jk,jl,i,i) += det_kl
+        eplf_factor(2,jk,jl,i,i) += det_kl
+        eplf_factor(1,i,jl,jk,i) -= det_kl
 
         !- Open-closed shell interactions
-        dtemp(1) += (  &
-           mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(i,jl)*mo_eplf_integral_matrix(i,jk) )  
-        dtemp(2) += mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) 
+        eplf_factor(1,i,i,jk,jl) += det_kl
+        eplf_factor(2,i,i,jk,jl) += det_kl
+        eplf_factor(1,jk,i,i,jl) -= det_kl
        enddo
 
        !- Open-open shell interactions
        do m=1,nact
          ik = det(m,k,p)
          il = det(m,l,p)
-         dtemp(1) += (  &
-           mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
+         eplf_factor(1,ik,il,jk,jl) += det_kl
+         eplf_factor(1,jk,il,ik,jl) -= det_kl
        enddo
        do m=1,nact2
          ik = det(m,k,p2)
          il = det(m,l,p2)
-         dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk) 
+         eplf_factor(2,ik,il,jk,jl) += det_kl
        enddo
 
       enddo
@@ -164,60 +273,46 @@ END_PROVIDER
 
       do i=1,mo_closed_num
        !- Open-closed shell interactions
-       dtemp(1) += (  &
-           mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(i,i) - &
-           mo_value_prod_p(i,ik)*mo_eplf_integral_matrix(i,il) )
-       dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(i,i)
+       eplf_factor(1,ik,il,i,i) += det_kl
+       eplf_factor(2,ik,il,i,i) += det_kl
+       eplf_factor(1,i,il,ik,i) -= det_kl
 
        !- Closed-open shell interactions
-       dtemp(1) += (  &
-           mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(i,jl)*mo_eplf_integral_matrix(i,jk) )  
-       dtemp(2) += mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk)
+       eplf_factor(1,i,i,jk,jl) += det_kl
+       eplf_factor(2,i,i,jk,jl) += det_kl
+       eplf_factor(1,jk,i,i,jl) -= det_kl
       enddo
 
       !- Open-open shell interactions
       do m=1,nact
        jk = det(m,k,p)
        jl = det(m,l,p)
-       dtemp(1) += (  &
-           mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
-           mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
+       eplf_factor(1,ik,il,jk,jl) += det_kl
+       eplf_factor(1,jk,il,ik,jl) -= det_kl
       enddo
       do m=1,nact2
        jk = det(m,k,p2)
        jl = det(m,l,p2)
-       dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk) 
+       eplf_factor(2,ik,il,jk,jl) += det_kl
       enddo
 
      else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
     
      ! Consider only the double excitations of same-spin electrons
       call get_double_excitation(k,l,ik,il,jk,jl,p)
-    
-      dtemp(1) += (  &
-             mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
-             mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
+      eplf_factor(1,ik,il,jk,jl) += det_kl
+      eplf_factor(1,jk,il,ik,jl) -= det_kl
     
      else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
     
-     ! Consider only the double excitations of opposite-spin electrons
-     call get_single_excitation(k,l,ik,il,p)
-     call get_single_excitation(k,l,jk,jl,p2)
-    
-      dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk)
+      ! Consider only the double excitations of opposite-spin electrons
+      call get_single_excitation(k,l,ik,il,p)
+      call get_single_excitation(k,l,jk,jl,p2)
+      eplf_factor(2,ik,il,jk,jl) += det_kl
     
      endif
    enddo
 
-   phase = dble(exc(4))
-   eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1) 
-   eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2)
-   if (k /= l) then
-    eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1)
-    eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2)
-   endif
-
   enddo
  enddo
  
diff --git a/src/mo.irp.f b/src/mo.irp.f
index a83c022..4b6fac9 100644
--- a/src/mo.irp.f
+++ b/src/mo.irp.f
@@ -43,16 +43,6 @@ BEGIN_PROVIDER [ integer, mo_num ]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ]
- implicit none
- BEGIN_DOC
-! Occupation numbers of molecular orbitals
- END_DOC
-
- call get_mo_basis_mo_occ(mo_occ)
-
-END_PROVIDER
-
 BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ]
  implicit none
  BEGIN_DOC
diff --git a/src/overlap.irp.f b/src/overlap.irp.f
index 3027a31..c51f37b 100644
--- a/src/overlap.irp.f
+++ b/src/overlap.irp.f
@@ -155,10 +155,10 @@ double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
  xp = xp*inv_p
 
  c = a*b*inv_p*(xa-xb)**2
-!if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
-! primitive_overlap_oneD = 0.d0
-! return
-!endif
+ if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
+  primitive_overlap_oneD = 0.d0
+  return
+ endif
 
  c =  dexp(-c)
 
diff --git a/src/test_1d.irp.f b/src/test_1d.irp.f
index 79e3112..ef845f6 100644
--- a/src/test_1d.irp.f
+++ b/src/test_1d.irp.f
@@ -10,7 +10,7 @@ subroutine run
  point(1) = 0.
  point(2) = 0.
  integer :: i
- do i=-40,40
+ do i=-40,60
   point(3) = real(i)/20.
   TOUCH point
   print *,  point(3), eplf_value_p, eplf_up_up, eplf_up_dn