From cee21a9e2a76369eca7c837bb89c35e0f36871c1 Mon Sep 17 00:00:00 2001 From: lgsmith Date: Mon, 18 Jan 2021 07:10:00 -0600 Subject: [PATCH] OCF tool (#39) * added noWeights subclass to Weights to provide correct behavior for uniform weighting of frames. * fixed operator[] decs for noWeights. * fixed issue with transition from noWeights to uniform weights vector * I was working on these on the master. Whoops. This is to replace mods to the weights class with a derived class that does the things we need. * Added trajectory accepting weights c oding, for weights applications that need current_frame defined in the no weights use-case. * Brought Weights code up from the past since we are now pursuing the full subclass strategy. * Fixed some issues with delimiters and definitions. Introduced _frameWeight ILO the pervasive 1.0s. * UniformWeight builds. Now to add it to OptionsFramework. * Trying to achieve polymorphism the wrong way. Need to do the virtual functions thing, then use a pointer in weights_options. Perhaps the function-pointer approach would have been better? * vtable problem was coming from missing lines in sconscript. Code now builds. Next, see if we can add consts back in. * Changed several functions back to returing const, as they had before I made them virtual inside weights. * Weights vector swap added. Now we have a setter, a getter, and a replacer. * added error msg reporting differences in lengths for replacement weights vector. * drafting bootstrap code. * bootsy collins ref too obscure. Wrote header-only template class for making vector of observations. * Added ignore for pycharm configs. * edited rdf code to take advantage of weights new default behavior. * updated rdf to compute weights using weights default. * reordered things in uniform weight class to get rid of compiler warnings. * changed weights opts back to holding a reference to weights, rather than pointer, for consistency with master. * Use unique pointer so weights becomes an interface within weightsOptions. * added arrows to rdf and xyrdf. They now build. * May have erroneously switched an error msg line but nothing else. * needed to include correct template type in upcast in postConditions. * Adding zeroth pass changes. Switching to loos.hpp functions next. * Reorienting things as classes. * Tool remains untested, and classes need to be put in hierarchy. * regoofing classes for computing OCF in various cases. * Added uvdot to Coord.hpp, unit vector dot product. Working on class hierarchy for two OCF classes. * Abandoning class hierarchy idea. * Repetitive boilerplate version. Uncertain if it compiles. * Fixed some small typos, and booboos. updated AG to use 'uvdot' * Added testing config. Issue with weights default. Need to go finish that branch. * Adding files to record weights from file, so weights can be base. * copied and name-changed Weights over to WeightsFromFile. * Updated UniformWeight to no longer use add_traj, since that's no longer a base-class function. * Updated Weights to have only protected and public members. Updated UniformWeight and WeightsFromFile to inherit from protected content of Weights. * Updating to test building modified Weights and Weights-options * Trying to figure out why pointer to base is not working. * This now builds. Changed name of weights to pWeights in opts framework. Made add_traj virtual method of base. Amended all weights classes to respect this. * Trying fix object slicing via make unique as input to unique_ptr constructor. * Was experiencing object slicing due to make_unique's template parameter. * I appear to have been using the variadic argument of make_unique incorrectly. Modded code looks much better. * Added a comment reflecting my new lived experience. * Reordered contents of weights to get rid of -wreorder * Got rid of unneeded virtual calls * Removed extraneous files in prep for push * Apparently the public methods of Weights all need to be virtual. * fixed failure to accumulate OCFs. Fixed interaction with weights. * Moved ocf.cpp to Tools, changed sconscripts accordingly. * Fixed eroneous weights issue * trying to improve the code and also to accurately perform the correlation length sum. * added compute function that includes repetitive elements of if-elsehiemers. * Finished making compute function to reduce repetitive code. * Rewrote IO, and mostly ready to test now. Adding one flag option to do every possible offset... * Builds now. Output is modified. * Got rid of integer sign warnings. * Have something that appears to build, compute everything I want. Now onto the fullhelp. * Added fullhelp msg. * Undid small accidental changes and inadvertent includes. * Used function pointers to move t-loop outside of if-clauses. Tested on polypeptide system found no obvious issues. * Added the requested comment to lines 255-256 Co-authored-by: Louis Smith --- Packages/User/noes | Bin 0 -> 59488 bytes Tools/SConscript | 2 +- Tools/ocf.cpp | 324 ++++++++++++++++++++++++++++++++++++++++++++ src/AtomicGroup.cpp | 16 +++ src/AtomicGroup.hpp | 3 + src/Coord.hpp | 18 +++ 6 files changed, 362 insertions(+), 1 deletion(-) create mode 100755 Packages/User/noes create mode 100644 Tools/ocf.cpp diff --git a/Packages/User/noes b/Packages/User/noes new file mode 100755 index 0000000000000000000000000000000000000000..a0026d82e21733844bac695fcee762d9f27f7399 GIT binary patch literal 59488 zcmeFa34D`P)<6DqgUZqb7p5vo!020CYE#6G^esf+%PRZQ2IfrjrHAxK-OC znHa5(3-37dzKo8F;!&Zivg+-}l^mpX5oJxG?j7|DWII z=cD01_uO;OJ@?#m&OLWYdac*4iQlAv*0*P|P^Ejkf6s+qPFGz}|Lg^U%J5f4S z8i9N&gYoiWo}iRJM+9j#O&9PKq*SjknFC_e2?9?|RSqdl6J_xd?_C_Krjle}0IEm4 zHr&s$NDb^oO}hnSeaG-}wH~oY^)_zc@QoV`>_ttN2>#Vn?H~1xKGmW;HZ5mVBuPyP zm!8TGytVY+9c0;I>1{#~YO2<|81;zHAN%M9;ylq_BDt&;9a2-#ymW4(b7?`|xsA2w zG&-Aot>?5B7MxR%m+f(7=dykieZrkRcRp*0##@Mjl^>5h!Kgf|Vd7y~f zf24QPBc(6ioq0OSX5d53(~ubS#}-Qy^vUBJ(zi(KCa37j{~+C3Ja1%)UdlA4PcoiR zoSf7q&Fi~S>H_^S#>OPXDyiCd|%^-^$GwGiAU$I1oA(apxsYK$EV{>z|Y+Y?29F~0HSU!8#ds|oOH63A^0f`NGS(-Y921O8Do z>=~JW|CtHs+@3%VXC|QYW&-@I1pWG-1n^pLZWzWr=9vV3;QmqZ<#TfaIsax%eEi?Q zekV&AQsw2m@#TVtrjZ=p)5ZbM3-~_;{yLgO@l4@~bPD((Iz2S8;(3oJQaS7{9y}ib z4#B%GlTMNh-5jn*|LGYEeB#p} z=zuJH?n0nNbT*21ReAnGjBid%yQ$!l@N?E1SbbHV=L$LoK_`ooms$n?eL~J)gFQ|m z|8;^upMYnhgOen^)O{lds5bs21VZpO0iPuBFB!r6)mF;^PYD6F0-yTZCiD=^u?K3- z_?HC#DxbFif%?@ZIqF z{j9&xhX)J{PnD-E==TV@p_}Y^0{TPk#@gd|LLYiWzq$qf$wEG{`Z+@I(?vuz-JGwALwx~_;{%$Ms8jqw|*fv)n2}W9Ei^~G2?zd;b|s+?Cv!HXnPR+(W9;ZN@N=v*T1r1xAc|y!Xt(P+4){hWs3iv^<4 zHwo*pMwiQj6_{kZta`!R{H4fwP5CWuSH0WbWOKE6ovvn&Y3d5Q+i73g=@8-= zW~r{Q)jHgcdZ)+ha93B9HM*J|)r2A0Z1qj9W}C-rcYAF%jv62FvblNGdR~dw)l^X9 za=B|g<`rOMZk4xSg`>vna$8i6ay$)ox1-k9;sq>H%WO7R*)3IhHnX|h#Ouzf^5)lg zSG71SHRi%vmu;o9*3oRPq5i4V3Tho*yR*?W&C_D5arv6PAna^vX|%Nn+9p$_!|im{ zI%`T@t!RUIVvU*#oXsm-HFlDU?65-urTI0jttOM%Sy@SK%$*ZqZkpTSX|^{xR20@V zRy^$zhqJ!H>nXI?*4n&o`!e&qxE-ysnR2R8%w(~OGlEePtX#wgXBGy^vC`#UZptgQ zdz>{@jz-{vK@Hm`60x&pw%g@vakiSBB#kQ3Uk zI#COyZFE3mH7!MWP*sRCpp}{68ZwJvJJQJrnKF}k-W=%gG@Fg-hADR`jSDQo?QE{M z)Kqy*rW&X-jXdN8J%tX}2sP8_Pd@a=S<5*6oJv4UZWQ&3W0 z4->U}9inI|6w_QuLutwrph|nqa(lf)(}HG%a`SyHwW!V%Y=NWERpazJ9Ucpe3aSk? z)^fv09Nkz(#v!_XzQ&RvY~r%c*-UC4qvWxOu+c^bRba*$P4vgQhjK36u>1U}3aSyQ z$-|77tkuJmC`NKesa0Y1(F_^28iwar_!_-VjjJKa9hDarnpWj4w9Q)z4R<&!N5gcL zH?PI+bX(?B<=gPT88(w+>id6=2iI^g|WtGJnQAE%w#Oj8@|o`~Pz zJ}-8zs{_Sc$4GDKedE zhM5H%`OfBAM{6QY%2h{yRhg@Kg~RPZAT-NPu_Uh%(-neEIz25eb@t$+ zqO}Y$ENY)jaRw&rcr`OP0#!~+Gr}4vXk}Ecm2b|lP3c8 z@Zey<#Fj`l&nhf@n&{}tlwB|_JBOB1yd-Arq(KQsBR{mjLIkOaM}}h;xp?Rg=Mf(K zTvHRWb#_eTHk__GY_P3*Zv&Gk_fo^?OwEDvB2aTR!wk&(>**vexRgHBhMLP&0C?>Zrb*?a(C_Znr zdpx%A9LpO-H$qsix@DUUO5 z3)&?rf|LrrMrDx+jXOBB+gf~Xhi!$^?e*Ck(NJw;rCPESi|%F{)}S>F0@vP1Q;FAM zbF|htm}Z~`X8i)2&4p>#frY$~Mx3<=#DhqQuTTWGrP@}VgCL>6Ryn6?s*T!5u;pm3 zZH%{g`CbBoL)29{W?g9l$!ZH+#&C|k4Q|&;NS`WrykrL84?`7!v+hy>#Zys3r7EnJ zu;(Zj|J~)=ypf zk;2wspu7(Aye1M36=vkLQvFgepYJXtOSoEm!zM13(eukAcS68pq^m`J|Fu()KSsJyv+Vr6&fio=@zIfvfxG zeKBzLJYY`@Ts^1hkAbWE&PoitQtU_fJQyuc^?ao-2CkmRC^2w#U;ao8Ts{Ag#D2ay zUg2}57`S?#mlXq7&mnSR;Oe_Qje0TSUp!IxNvy{4xd{QyoEvdKk9iO!Nn0M zD^SnH2;R(5*}v*JBf;$&xIyR@!K*d!N)6nsfp=@*Rt;P|cP2U)Y2Z1J7+4O+%508P ziwZSxnz7ZVSOcedS$$*;97nw2r&0r_w$*2$22MJmK1(!k;z@lPG;r0nsjNiYv64fIF9tf&l(LpBMgrqQPtf4^XyCulz;|lk>RB%p_i5lKY4G=G;Odzu758i4>Rn%gD;oGI1Q5@l27amv zM0!L6pQwRL!X9aC4H|g520lpx*K6RDHSkOg{4@>Rpn;#RfoEypXK3I#8u))`;Ds9a z6b-ys13yy(mo@OSH1J9d{A>+;p$2}A2EIfC&(^>jH1KmZ@D>f+q=C0;;JF%jn+Be* zfv?fP3pDU94SbpgzE%TQ?^00lIt~0h4gPu!e7XkSt%0Affy?a&(&a$f6Kw`b?(Fd< z5B16I+tRm6L-}7{V~~bU{}%u1CKlt3@(olo)IWsh^iL_DMEO3J-%ojRVWA$D-$!|J zU7>E4f0Od$vO?=vemCXGRfW1({sqdDiwd=|{1(cS`i5FqeiP-%C54u-{GTXKQ(UN$ z%99HUb+i0slqc5{TF3J9C{HdY)Wz}_Ql3Jl zP#epaQl4B-sD<@^4a}LV(aZmfuZzavh;AmVbfr7AWPc9== z$?^|Uo?Jz!nC0)LJh_Na4$J?R^5hyq2A02-^5haidY1n+<+CX-vHbOvCl?SJ{FaSB zb{WIz&9Hh0SY6nP^7$c+GKpm*6{mM#a&3Xvyfr5$C3R94B$z!ea(fUs~! z*+MF4Q491akWtXOqpVf3`ZugW0YD8q+PfG^Z~IzG@yNaH>(~o>U9Y~PcbUqczah4_ z!zWI{_yLnK2$UV>ERfsKjo0Uhd)aRWD`vE6&IIX+t5mv^8vl&iA4qOE#4(2yw`HA-E_CGy%T^M;?B3`*Ll2F(XDWr?S^eDNI&2uun z?fgn~Lzt9mV!^Zvyao*3)tPeeY8qpR0%_8x${FyN7@+UwQ(`VeF*MQrF;KFUQB)Qv z2b-#GwD$-zbfLS?Db(bAG->ho%22h17?>x2Y2aP6^v>23cNYIhOK^5($&T5X5@0z5 zL#7Y~v(r&hSkgP2c&EZjkuG7We|8qi8-UzEHClmjLmQRgLFI*(;0@jM+Izzq24k;X zk=M16*LC8xq?c8p9IBB*(_``9(8H6Rl=R-vhlCzQ&LXAx&>E;U=2-n7DOpUukWR01 z3VVH*UXR7AmAqT8G7_(WJR{s4p*xZ5=+Sk2{tFW1)9XM|_aw7j(oL)XA!9#4a&L!` z^pzfq|6A)>_Ymi6;6IV*rgRtwky`!z^5!F{6wyqUH^DX1qmzsJbT`}r>eSkcUosV3 z2W%NG&Kh9e!IHG?Gu<`F`v1sUsna#~Sc`g9y53URzuD^FQr9JSKJSIu=gNMQ9B`-0 zorh)rFVnqggYC1@ljQa-dbxcwM*r`=*DX*3>y%!({i8uSSY8N!(P2D7tp?g@Kq#Az z{vb;d2FnB?TP0$}AQoV1Oe(Vk@1(~JsodD0VgzqAV(Snje2C@d>*ZhxT6snGzX+X` z{bP(okr0$R#0PTl86#P=9DJM$)c>aJzsguC`@Kek;yprQd4|&qs`JPrB&=YX%|&MM zjFAjY4tS|fhmmz6$WyByzH&k z{-pAveY(|V1d`=YMgLH_zZP@Wj4!BL@F33^yHHEH`Y`2rrlfsaQn~+w^1!*u4hD{C z&EF?a*-7eDQnXEXJ*nBiEg1hE@6|LZ)1ZB-(+fHD7VD7*U>pXv5~FdOu?KllT1x2w zTu*Sa9)c4AdgR_>155ptP^{aGeGJyYQUNy*KGkFS_MR+~bgLZr9Y^kb%R73oeR+D4 zuH#_@69n%8rVO}Vz~=+5n^;bT(YSHN&1}24lX!dDQALw zhK+h)qH;U3X0Jq9^Hjj zbRDOVJc476m5W}G{Z?a^ZuJvb3djg@vIgNAXnM&~HWIs#D8!FetSMe&rni{{6*I<> zG5fScH*SQ*zaZ{F10FTd>c9aWcYpym5F;&e(Kil^@PgHku0OBqcoTMJ4a`f2k|xW& zB}_lzoU-7gvY@9Wd#wJ8`+=ayy``kw@QxQN$PXglPkB^2LMh;bfFA)|&C~lJD=T`# zcZ(db8Z+kvkN3+(2>YscRq8tKB`F8zKv8$dQ~oOZKarmmLo`MP_ZXkQ%Fn9Y%acBl ziF<1>?*%q$)0m%`L#(HuCBga@u)dX8FWn>iFYK2CbHI8TV;y+~d9YsAPi(2I0zL>h z*ezqMBTw&w^-sY%e0;|@R0Hysrwm5qP5z@$-mIjQ`x}kL_{d;-b427;|8|T23whE` zIQO{Vu*zL&^U{CkJy zk+F_R@;yWOVA`=X9bzdExG+;L+N|sR2y4`#dFcbc=XQ)yw6uSdrR%r`hVWDtnhRy% zjYMQ={|@j;5bq`ie++Ql#15dqIbxD7EDxYd_{d!CShThE6EX)ff)p4BCPTbhF*RiB z#uE_c*oNxN=Z3CFsH>YudH9nQ;=PjduIv0Ai-J?+;MjRs(=z_y){0)xbsC8O%Jk5$ zfKivN@)KKsX1M;bG4;QIRZ#s?V(M#Ur|VotGGOxCjWxDbeveT7dFdgNm!^Gk1tcZr z3POV*&gC1-Khp|maK?U2ew%3hu5>XGdEtnd)7hvx9UyH75@s<{>_Mdgu7&2Y79rAw zXd{z^)wZ#p5?$v_*uSU~{&+D>?*BZ1f2NuK8OkZql)pd@9E{0Oj zf^H&auEj6a>9StK=OZ?)=>f0PRiC~bQmE5qe~k(Gtxt7V^na>bc=|}dKGmJQ?Ni;E zdp^}&x<_~Y!yp9`-=O+_eBNawGC7G;vN(w^Kw>vYpvuBiw#{KAwt_?-Nc8Ehzl@Xk z6xBZf8VEVnRPmSak_3CK7EtjB=V>i<$%4iZ3; z$+?LtCA!5Glx@}HpO#FB>3vKvi?7zxzdq!%LKum2{=rCIp`eH7bh7uew*&VeM(1US zu^&?C*R@ZDdIN(EW+X$ISem|-`K zNz!$`2s}9eU~p_JcDqm}2ZUvtMsdqr4-`v)g6NN^9*Li7UMIk)c>!voG+c8Kc0e|L z1yFdqfwE%gYOnzzL5hEQ@lYiU!ZU^p#6tx`l?Z#d$7#Ttq4RxgHV8Xv#MBMv)?z6- zpzCnZyzd{7H$`l8?8VsBvrc(WrjX?$*}n&&9og%aO*6?ZbpO*Mw-3&g*UD42$Vt1E zt@!|6_e~pyO?NYW3lHl*9VqbXN-cJ%$d(Y4z@V^rfypxO**?IcGpr9gjQylp2rP;z zGbl20C`dj)_~VM8*O(4bjYz^!_F(*v!)NRxGPFn5oeAFry@rA4)7lB7yJ8EJs@hnn zR3Bh!gjWEtI|C&+p$#ay8`gt52xddkZzi@iia(R`#2`3Lmx~ZS>pG@`ai;y)a!R)b zmdBWY94I3bVBwLjb1vo-EuSNX^ZC~7C_Wz=7|!P&py*cDFj`q#;5~0xW{{NKz{I-Y z`B~UKgQy@l%&0q*lX*wcMY5iUInp0vkPb6&hjAS$TK%;0B!>0M$01VM+l=c;_Siw= z9$r$hm=Lxr<9X$W35@1Y`H+%KmEqgrZLTsFTd{tU1V@VC zsGo6sNZA9T=pY=IQUMfd>5y_9AqVppt0YBB4-o@_^I^k>l`e%#kNg@*k7gIDx`s1( zS`3rFJ^G$K)9?9F@xZ;;n&i z^nJnRVocV`i=rS$`D8k6{E}d43%FLEns<5@dadjHQyQpTw~h9W7hPVmxa9JZD@trz z-q2|PSg`X;?_%bKhfpTwg(YK9 zu9+8X=m2?lKdatB)ulM|!(}KYFy@C-K|wn|j29H90Y&BSOGWvRDvY47vV)&WtnSlZ z8;#eZExL|>Uz|LdVtTIja- zU=P`gSWNXDyFVHlQiZFD7vb}l&~&7Y*RI757LV6HgQ$N*ymsIFBs3APJ@5gK*T#&b z?tD+Yb|UmOG#V}Nc#Y=s*Dy777+c0A#A}EvvJly2A-dAv#zGa47~6_!K>vt)`YT=lsjhO{R6mRqvB|>M@gXv3eVi4$3FFaaG>MG`m=$mk*mQVz@n_{;n7Qs0vf_oQ&J!K*z$phf6#YFz>t-yT?!u~8A2)9o{Fa#NV z8k)dg*dT=#Aj9LfZJ{80Z$qyo9=Fvcj@w=XuCPz~uC~@bA47ez&p$D83GK5T#oRvg zuK|UyeSYd|q}sJ8ju*AL;gU6x+I~39&;snU4ebuOU~ouU*#=W2vK#m-fQt|x5A6y) zPTMoPLc1v09hySG-J!40O(e-^5y@n17a}=}?HX<+_e(!w#NDq$yb1Tag4_eGl@zNe z$Whz{B#phwNf=J-C=KkuKQg%?xPO{yRQJD8S-)>+Xy7qc9A}&VX=M$SKdwCX54fHj z$dozRt%?cxftLmlifqTJP8L*FnTvHK-g|(gY{7aM?`!D&N%EY5D;wmZtFZT!rWAl8 z@-pT1$e+RUaK|tiu;bu8(Sj+PPV3b3KQwJ=-)0!t2P)WWptGiZfT6$7 z(vGTfdp~6kDW@>PPxpZdVzYAyt0hQPF~hfpuD}~jv6(<-9itkHkWm&;HP)@Q#QsVg zcZber*w_F<+{?(}>@d`dp#q|k#@pFQ1>k0f@;;{Zf#;NidxwVFUl~&FWz202*$F>D zhc8<#Ma1>lD0o(8^bRR|fgM`ONNybng_g0mUZoY_Fn>=YqddvU(b_E?qoJ%oYveyb z_mDytg!!~>BfhuN`~Vi0QXO$NMRJfqU6}#g&;`6j89Q*A-~e0wMSu7HPWd$iGH^T7 z4n*UmJv^@93N) zs`pqKg@dApsb=tmx)CT3=HGCs0W}UOy#$SKEosWHSm7aM9vBbZ#zi-PkbGtg=P-npIZFWh1JJNHk9BxjM9f9IWXga%aISB@Fq+0zrBJ=x21>ND6Q)< z9klp&T_OAT$?fkCR#sbsX>S%vKp$Jfk{JhpHFQEV>VD4GrNJ40KxRH-l>a1xMMq>m z>=io%XHpiY31crG>W{=1mv5z8H<3K}vlC)4|1w5q3J%adU^`9zLF`V>l#_PIyS}3% zsk9_Pa10St*O#&M+Gb2aha}&*a{G)&7vpU{@q-2SjQz+dr+|G$Mi+LZB>^vqs<%Of zq|Y!>O)Tf=mxtab8kA-{SwgO7{N~BGuQzp|YTsV$yU3yPb!9 z$IupU=WjsMIt5bvS||HomHo#XSC%7#WxrYSxl43U9&cO&V}s7Jfb;V79jUYn1}X@Y zIZfXc7#BT7$Xg@wGT~?e5$qww^9~K>@i+scfE9;h{x?}{+p9R=5*h`{1AkQaF*|pA z*ABc1aZjEGi#Y)}Nev=zqxmFJt?rjC&h&+mE0>Ft1eO z%0gt4e3zpyW@!yh)Y#mD-O4ibB{{Y)Gtiez9__ULZ-5iS3360z4bWdG`(d} zk4|1<{TJ)6Vul$^y^ExqjtFySNpP_-2?r4f_m6+x;(r_al9dQbFD*v=c`2gramAUs zadgjOr#z*PHoq*oM|WUR&f=esgY}AZ3*Bv@eTDFf)aqA!r5G)Px*FBUh%*n3 zJT%m`7M-D~3iq8bIk2Fk=?$nAH2DdnO4lCD)0sdDjmI1JA9CPYJ)65j$pnKxP77UR zkVuPV|45Xk%fT`%-RM@6ffHL4y7TX0{Q&+`iSP$Oks5CRI?`upE?gA)Epm}Dt|CkO zCa$jv&W2upDEFqa1va=*CK3bWQ=VFZ@1W6720_=!R>sq2p}z|FvYw8Bg{jlXMCm=* z3mhQa7A#=VR^^%~UP#}V;{3BJF+(G*fjWa6xB~0gf5=I1c#mNl6BNA5?cc!nOfasT zj<+OV6L+HUd=xu?jq`>1315mqy_`7V4pf~l!v1qA2`{eyobo!AaL+t&ci4Zv<-LRX z&waQitGxCa#(NzQ2islgN#0{VQ|+$`1`L6#2S(nJi%4i-z{qb$zdV(jKWW0CS{_t28n$&R7FyKvIYRnWWMv9K0qhrFFB%=Q2N0ALeFvD_q&0&Fu!cU5OV4`|~z3x?Rg|w6=RE;Pn`i=gwIEXYNC(4gx!H@XgQv%djF#lUnC4nK3wR}ai zSzPy5imAvR+qWWG-3uFfo1vd@I|(|N|0GO@csNW@&dyylx-gm)Hp;vz%|Wd8e*w?x ze{TS5d$E5_{*Lxv<{F4qcs>lUEr-THbLe>|w6h}M$g~75*ISC-_FiEP=1=_^^Vg^@ z;*3DEu?LOw|AXF-*&&S&X9U90;baj95Wd$1fYifU^r^3vYGVC~%07)hq9&;6Uy#Cx z&9l%)5V6ox6{hMV>78G5emNhJ4cEn6nHcw~>Zqf&?lhEQ$vK-RQ;?|rhvNOj{)7Dg zhy902?o|%GIyAHhuK8I^<&r9l-|p6`*!T4R2vdX~=s(7OGx{F!*9Q3O9PY1M(7YU2 z#KvPJJL~qQFB0>8bh|!<%luSoQ`r9~^g!4@jpfij5q}2LD%A$*4m+D9G>v|P?m?H( z=s}#@(Se`$WY{ptgRLVMT^?Gyi`0Q!CTu1~&Hly!!Lhe}NxfIkxO{0C{eg``y9XZE zS||JgjXxXsZ>y{UR|}NC!BX^&*8vT`7$gSts()m9{?1RQ_U?@%Vb_0-q79(Giv2#^ba_xh(E!7YHRS?93FSCq`3R368uQO zN=L0>ldGAW+&6OY+G4oFh2??M-?`9$?Ut1YWE+hJM07i8On1!gi;i?A^kOb&=Qg;?<{vwJcP}4d)KZ8;r=)dq9wO zPJ`|2e2GP5H&I)Gb~*yH2JT_|9_SVPb1k;9#HJK&nqZ?OO|b)+ZA$ga0mdp<)(eS- z{}67?D0I9+-PJ4iinny-F7|e%fo}($!Wx`zE&37!^}MFCkO=#?4?Kn0Z{WGGyx6Qo zB2&--tN(+6M}`-@EQ)ZchWi6r=h4y?iyWJQ^jbKsi*S)g>&QluPr>RX#*c-#A6tOq z1=Wugvzu+)11;i9bBZIP<{jT3g#V*@<5cOQ(#3!{nnpf2Hm;C$-A z@c;zOtR-7y5n#c}NOiz~OA- zQZCw~>sSpHv<7Q%%Z<1-f?Mkp+}>4Etih|cW(^$6{ai)Sw+J}A3!yk!?3g;@))5Wy zx#cZ>j|#3^JD2;g5}`Z^lL7`WgvqJ-{uB`);sw?pxHUFML=->t*h{oO6f-{rV#1Kq zDJ<|{Ou~Mi)FOOnhXFf}d%24k4wHxMA0&~W<4h`r>=LmTmlHeIu)}yxs_ zB^r{c>-ZCGeFefedR@otXh^!Q<85Y;Et%{Dbpf|pcU0@xs6)_c#33~B&0hajgaw*W%B?tY{++n#9g z9~#&s2yntf4D~k~pRI+)QIga^M6u2@xZyS^{f?z6}E~uT1>l^GFUtABU$ROx_ zz@|Fd?goF0$-U6(Fd3y(01j!F9$ygbl+9 zR|vv}Xu<`8aJH7Pfe6!gvWVm3wS?grVmbbNCl^IlG~ss$jOwo;4T~a+2zLs?PxBIp zqJ3MT#*&3zlQ)CA^5H;=g?(S(8(%hTn8#V1!PjATbfmh^WbxxZ+c?_BV@iR|1@YYv z9uWkvWbMVbKWJj6ITWsv0vdw~Fbu4s{Q;Pej62jQE<>^SI+gkjs!jwmRmHNuOb*|5 ztA@R!P9eK@9l%aqZ<*Nj5?{$;brAoQ6^83jbOg>*=&BG{hTR62Lq_kyR!W9mqiw4d za1rcgE4!2oEe9{)nemDj&UEKnd|7NZgm1q-$+*Ep%61x25BxwHYCX0OhhafI=)GAV zt~Y|!TWmCdtyCj;!rBmutX``94^hk}Xcc%XsFIyzg<~IG=dFa;zAZ-`*R8SsUDuhy zx(50@_{Yql>olPn`T37vhpL~Co)>hTccLtMZqRi~qTZF-dfGWk*LfLhj`lt#4}d^; z?$>qh0{_%#AuuBIRrDOE>->;4OZE5dkFNh~QGXjS64&1zu77t_ecxhg`+lk$-k!)g*iNE` zz88I|>-ak+I5m*42J5r5M&atcg?8z^V{qw@ zU2pvg^{Fmq0jsf}>Vorfb|}m0vOb8q*zHApuIuoDb5Xz3iFZq2KBlrmbSE8D=@w^| zfMUVls*n@K+F0C7KabQLH>-1~`?TV;`me;%G-h}542Uh@T{wT(3)FRxVzUVhC)j!v zVTOe707H1!d+9PHE!tPW)Nvevp*YSuJCrO22AmM%v`r@Ln(m-c$JzZ=)^Q5B+0C?I{ki3B$>2DutzKlf| z3!|~bsggGzoQdGk;(tj_+K#EGp#tBs@GG>=#N=Cz%U`>WAUIa{d#^GsG!4l`K%s~- zcqUnl|BHcBxLxD=s)3YnB=@TRcXUBLVFb-d)?lh+2`vpLk zc-3nV(<@e}AGTM&D;av-N}JssJutn85mEx#5^MybuKX)Y;3C8gD-9N1`5t_I(Ey=C z57nuz0#jLqLH7Hw+IpTmx|4LI6Bdtph>%w5fy-tHqRTmS&0>%%Au@@=Aep`R5!=CC(t-*?We%O+zlEm)Vy-ExqQ~^pxia?&ZQX4h>l8x_T+BA%WQg4G8HL*|=HElyb;C36-s0o`w z;r%*If3RN5g2Y6;5Uf6#4kK^n+J$pUXw`P1Pq6od`ZOm>pI+lS$n8C$KFuM03hSVm zciT+UID)NLLp2EPBl3@qzw~U3DSo<}-5J^|_Fixc7i7p)e@uQjd?qYHJa&G?riZ&w z%z(kMOBe-p&sf*FQc$2vP>I^7&``!c1>(P$_GyQY+Q0Ps+kcO#AIv>CwjVccMA*u{ zag`=0^l0@jIEF^zWr;^+yFxF1gO$usU+6QOM_>+`30YD0T1=>N0Q$rYlBN~?(8b&; zp8_UXrHpGM0Jnby(+seFg?it^%AOhzX9*; zJeDTkDc^;+uv(+e;2LlSM_;zP=|PB^Zsk~>QR7M4*3v!tA?_EGkH)q_xe-@)w$t2A zv-YlkvYAp2PUs!D0PzFtoy;yr*Fn*g+zW`Zk9Lf}MlCZKxh! z8CC>F!0@5zUYb&@^o4EgVprq02YUIk&+7k5vl29r2eD#@S}Z7GI5gXzRA!5Jy<(i=sb zT%%qg<(@}JhFor)JAsOLN#h6_NI>vsc#s-apqi?Ul!L@1}T*O{~h{aPy zM}jk!yc5mhRL-IZ_Bu#N;OB?vO$m2dA48885i~E21jdi&ZyPk3i?HYIX#V~n=-th9 zUF9#@Kkj1Okbk+h7rL^D6MUk?aS(>-I(vYGz^WKqBgL3+@6uSMDJ{i0=bc_`k&L zVrl<4i!ZNqi}dMErD^evNekxbMXBx8W`FvrCS^r+_OlHAlWI3I!q%-EbX{B zjecUoRM%NUD2Rix4n=2hvQ&#h(Ds8_R{DO2awcXr%X6xm4jqCDkpIQ7;Hz|qo3T8v z@b_E&pAa)<38QSnZ+yUE5DH}fCjGzwFyRDof$Hxx;lc$LgACfoSLz52Ljw=Wd}BCF z6X6|xVSq3ysS*$U{2j(NR942}2!y^_N*S#i+bf)Zixzmj0Fqw~M?Mq>M z9@kwkQN015gwzE7-tV-1QAKW9uCj=BA zpfeCA&QqXwP()FN+hpW^ak(s$*knIav4HVzrUe397_d7V0c`n+*4Pa+RtHBJG?=_z z=Rz9Ss52g2XE}r6`|@M={hbtJuN)}F9I+D-30fh~G8q`!_)WDbUFW652kj4WyMtw9 z`fS)LAYF!RAfXB@WrbydDd^#Rs4=AA$BGCY11{7Us~zlDq8#t@q0*2eyKibHm98K< zP@~XVcs8bISdKFKa^T`Yyjb>F=qFp~{!qyvH6RB{j>!It@gWPoxL<7m<#>nD21*o^ ziW~4GyhAC{0Qf+Q^9Fg?{e}5orSu5L2&3IdU$_ z?G9bQlgYg7020OZC_9~?SS-z0IA3`fhaS>)7JBY(s-w3)!OFvVCJvwn+=Tv1gS|<cOte9 z1|C-TyHT)?7x+F>H3|1W&|lb4_N*b146KRRANIR;hkRg9Jx`(dnxBVa zie^Xen)A>nn5PJzb>ckK9u2{H=-6loHZ$pG*EA=ja3T&stbqpFWJC~(4BbX!=b^M2 z895JaM@~HtEv55NbWJ@E#Tnp5IThHCl8avPW-7mVghe1q04a`Hy~kmC9C)dV&OKW> z$cM<3klBo4}?@7w#C?XlpXDCD76RiGk27=@l04(M}-$3LTGkm@|QJIWV+95xQ zNYE)Bon@9VNtRRa@R1_%65ofGz}4W^oWmR|4nj{-K6scaA5z{yN|jH*3?KRuo{@k% zInXx{bttJGf=*P{v&vs8>qO>d zbeaJhd>a}l=f`&m2*s86x1{}~ZItKq{1A{{E5bFVQsuCgy zhWQ^zg4kcIHa?WNB0RkJckdTxbz~(4ob6j{l|0Vg|MOp|7n3OR=;Sjj7GUgrAlouHUUzy-H6 z$9a(ZJ5RnqG)eP^#s7EOkHoJx+%bWDweXz(LRwWq`xJ+*Q$Aof?Kb~2Mft@a(SHcT zV0ufnAFKasEal5-$#eYM=yVeO!o&XPPH*js%vYanB?;r^9{Pz_*acJ-5C5xM{OIU5@6D;x$N{@R%r zCHTvT2K)g!`cusK`@05vbFHDq-b{3s0?C8>zh0LCe`u=KU^mbop5}j9I$O}6>oSPH ziMh(a|Ddmd{w}to$ho&q8$7<27MI%#-n>NWG=sz#lJKXx zq)Qr{H4O$kV;zly36#(EH7#|x4X!#vgZT5S9>WYnQ~~|{R|Y&rI!`iW86?9LgUV<0 zpWPh(hd~9NnhgJ|&1wlxOspStt zL;ZN~f?87yu@0$%?*-4TkZx{nraEFg0l`x?kK)Md;@o7jg^U*Km zuS4kY5K=t?rahELnvQ-9;@0)|JeTVV0CP?Xb(icL$ZAjN4rH6hH zkDjDo&6AQ^^+_j>8JXUdG*VCS3Ha>UGBk7+A&4hFKL6V=gumv$F(#0Wl%b!Mk$ItR z^vd)$Y351O&zx$UOeCm&IX-Q7qwfR|4{abV#%B-4nyS-7uois&1{e)DJp}8(rw=fQ zkUa#WjggmP>OKV6PQV}v_7GkVKHC7p5B;!*V1J{xtwTd&7=*W>;QblkK`q!rc;Di) zF{TYY^z`us zd_NA@62R#D2=vfq-%h|512%;Ka)!P=d2t4fM-D#kqfElM_aPxZZcHx8$XuOLo?&QD zot=@BTr(=e07!}2F43cP!TF$vAL9`8oQ$3US+Y%P|9UDzt@dO~h9S9mREEAhBfVra zYjYuRUII?MHqq_mdTpB*j+P2R>PFPb0-qa*F8BEclt>9a;J8nrbe-8{M+|4Ml%dIqwn2HtykPfddn#>j(s zTc+i0HmCnN>WzUvp|8xPl2J$M%_2FLXH-g+D4PI2I)EnuKXn_&6XOmcX~#V|zGPtDh5sRc{%V|jHdT{X?B4E+6mLVh}=RJK6K~d0B9`7dukqNFdcD1MzcpR z$G@ecF?J|B8=rn^55Dj^!W{!#jPdGZIU}<@Wp;)Ebuuj(`jU+F3J|A0>tPsOz@JTh z?uQjaHloj^(S1%iAt|PdmSLSP8BH|FU*8Bn+zMV&yEx68Ngr-XfxLsMr5RZ_rkOKx zR*xvjC~O~jSw>e<$|4XiEMb;}cEd7&hHx2>4{z-m8hRdi-X~~Ldx|AvN%9}jDfAMZ z!g!aX_Pe*I>ge^4WM%s=yYHn{f_j4ez0#&n6F6RBYIJu(US6YBDF~F$w=S2-Dnq8pTy>DpfuSip02ltV&ipV5F zb9P2wQp(ry`DC_I4PM^|t+PYiR#JZlnvnBN@bgjg7idCnVE|DxSJ~(mA=gqWU!8Vi zYB1#{kRY9Efs^?o=zd5-j_?`LDJ`EZDZR-Fbc*a?4aP-Aj8shh{uag%p%=g3#VNO? z#T_2Y=*4m2rOX!%qMh#l#1G$0L_2pywG-hxVhdFn-J?>@A0L+>lP}G$tq5W!;KB}D zCqcGXFz?H4LK=#Z$@XJ{xXBMS`F zr{JXW>!e8(*(Z(mCzay<^&3f0cx1=qrY)0Zq+KmJfQWp&&uRB3rM;6X1(VaZq)JaD zr~OZ=v^hC#W2*FWa@w6K(w|c~DQ(r!LrE!bjF3JW!HbTW1HSGPc>4R2U<(7ODZ4VH zy%`)~^n=N1ca4<(Jd$U}|32v})Oav0?Xj`a?db?SpC7Af1RWf@>!IYNKP1cGdga*t z)1#_!06SCCzAlpdscHR1Nx`(Vtwqv-w6uE)rH4nPwM|cY zfLcS7-<&6XFf#4%dD3<1X|J9q-Iu-(4A)MS5K4vOZQFk`U|A3$EAIAzBF+1_9SEuPL#hoU)nqQ_mp}5j7xtzU3%g_EXs80 zmaJ>noG;xoWybpRrO&3^nzZwL>4|gB0McjKXW{*gb8ny;PvmY%dg21<_5wXsFG#_^ zo2Stk^A|;F_nap^dtMs&xqW)tj_J~-8EN-iAZ?$K_SFpO^BHNt?7kokln!5zhF{XX zZf4pOGo|jC1V>MPetwp~&l31q0zXUOX9@f)fuAMtvjl#Yz|Ru+e?S7oRlJk>d7`r- zuXcUsg9i4x?h+2^Q&aqjMeQ>v@C*w$v`yf3iL^)H=ZNxRkyeUSOY7(NpC#~ry9C;9 zjP5~qOkOSRi-pr4ox!tF*b~oT4?ZTZ!u1RJ+d7e|@Yu9p*xz7`-NxcA`HA+}eVL#j z>ZyAPv1u-T&z&9>PQUN2&2JNjS8DmDana>FWAYz}{3#qwdYOJU7ta`;Na^SD7txHA z7vT95MR~*h2KJ^sCnUg042u0YQLditVVjaYUwk*@q8B>1rU$!2?Ab{_Ife&YkL+o| zHatD*Io&6MztJ4b&gkg%1d&%0?0T^$^AQ6}MKjV`K@Zzp?7!>*@HbX_M9qG6%UtOSpP5%41^mkEZ=<^N5q*w%L{g-L6P4ynU~|+iL6{5 z7eW5tmTcT$h)b*BTn?Tk(gh-|6X^<(t`_MXBE3(fPl|MhNM9G}hax>J(h)hlo$(?) zL!<>Foh8x*BCQkY3X!fB=^Y}yPoz(ZbcaY^7wLy0JuK1@bio48c#)nV(gKmr66pew z)`@h5NLP#W4w2p`(kDf_L!_^Z^h1#z7U_sw(Y{E}5NUx(dFj9J-}G)n6tUT5W#<{P z<}dX%dwqs!*?HMH=a_sfXS&8zn4OcCZJNU2QPtvRQs9jWRRhytq)5q9bzFFg@qhH}i4 zf?Fr5Ki3wez`Lq`JQ_Y;Y7_oF8vY9@o_|S!w^sdeG(KjN82zO36wD%qxNvp0iBO^W zE>-$uC<5VcG;btB&Pzl*-X{>QOaLDc@Ekg1C=@yyPTAErO23wW%)ksnElqbGA2e7Ig+!sraw|7Czr z!c8AFUOZNGxDD{h7)Ny-rN$X_-;LnwY#fp9N6_<*fLAW!c`(PGF9DwizWf{yel~+=NI3ycU#N*h z_Yc^(3p-Wg=MKOpbJVDR*D*Nw6zh5swzK;I@!Cyga*HCtGA9F0{p}V5OZPYE=>uFZ z4fm_>F!)Ks#w7&!I9@gC-;s!3nI6_~Fx`!zrwZ_ReBPP>{zwA&_5|?#3E;;f@Fsrx zmU2Gm{trEK0guOjo4{}D;COVWfu6?zC;OGw8(2Pz{R3$_Wze|AseggPcfRjAc{U)J4 z(zMa>`ELN6VB3mYw2wPSNG4R2>j0k zToUVMH8MCE3sy#7_&wncq&&dmj!Oc#UC=QIeY;n*dtUpK!gv9|Z%C*X~(> z6QAA99Fgus(L;azoa{H&e|`)&&1>qu$Rv&~9j8k`A8?|t?nfAeo?Mgwe@Oy(E8sNl zJ!0IuMZ32rz~2NowHvEfwODw^lf#h&bZRiM#>4*$;KYv_f2#KTvVgZdU|<}o_V`Z$ zPv5}d-GcuU;53Mix-SGY_M8Ja@!t~T_f`pbx7hzt_476XkM&D^fYZFt^J`AmA?Wu5 zPVMfD;d2Z+8?V2401`iyS8xJ!CzPHB0nb^%^K}1>o)eFYPbWVC{5OD;UP(b-zlbMN zUjqF741c_|{(b|)lm-5qtiQw8NB@CwN71mS9Q`2q^ow{~^#h9&=*bGenH*x`nVS>f zcME)VzwQ{p!9KvLzv{jrvh4Yq!N*I62MsK*>eX@BV50F-_vh7ppGv@)p2zry+X0U^ zj*kmE>Un@l=Y<6LA%R~h_)+!w7nr)J-M;ypZ-g=I$pt)KySw363~}Otn+2U75of+9 z=)5K1vHGUNMp-=imn48+1-L=b)FyWTPVr&vzQjF%Q@_-GR@J^<5b#(&$3uWjALvUF z^a#7i0X!a`Zh_xc%hA~#T)bu=R3bXX3piYDcNyUE=vleVsLC@8FgH%EPfWRMI#Q}83m7aSsfsTWo zcME+}@ecw{{Hy!acXD(o2OGNtSI@Hw1bjB&lQ6#eZUcj)3wVvdSNGZb1UvwEJUL6v zE{7+(rlm!)T{gGMJI!XRX>B!`Ou5=GKk24 z?TxltNX%oi`&y+M7k%~6>!{71R!}r8K9UXJMs(WjZnu4v4c}38uafHA_9lm|*4Nau z3RN^Y8?h3FYO~FnS5jfNndg?&8jk2+p0@St!C6LzjSU% zg{2HNX3w2(Gs}XmTs{wysID6R?L%nY_kV+NTCD@$19h!3lPtxz6AS3?so?vev#5tF zUGC+kf|B}rh|}(Mh@vu+wt`9XEz$B^R#I*Gfriaelik@IX;e$Y*=lyOp*0oRY~BX9 zYo)ExRqw2^IoxiS8wyigJlkq1EnB$IHZ?nsX#7AY-Ht|w-Qyt6V1RbpQm96=t;X)H zX(%oRnUYeA&6GWrC{=lj9IZ7DCIRZpu%_}FT`rF)r#k$d)^d_XlcTAoWtBj-E>^6_=rdWhWS;^H?=U14`S>Ce`AY)$qSA@bVB$`wso zrn8ICsp^U{Rl6+J6}E`-qkvSun#fT1IQIq4<`ph?_5`e0F(3&psY zay)Qupwr@YTR01Oq%vl+xoVoN7XNw+R0gwY=T_y|z=x~FW6HBFKpSMK ztRYL~tflZ8K$aOMoMLqj41~t7&^B)=9I3-$ISP+!uF8k6E3#Se=YYIU?D)#Z1Q5yb{A0rN$8R(i&Th z$H$xvP)xZsO)Vs-O8SN}`ukG^o5ut@jr#QktwVKV-YMq`I7dNS>eok{cJh5HX40%C~%jK^1FdqTA z6_7I|rGh5fJosYw|1ZjwZ?n)m@qOA=W#tNHtF?RWHsPdgOFbU-{jzF{)^)@TE$L-g z#ed2u=Ce_R&J{rCn5PL3*)WGOH3uW&BvHatgcE@)@YJ}S%%gA*LJZd|x79Q(x7Fd# zVn+0VInpXuqsQ!w<(`jJ#F5q5J>KDdETZhDsWlCbn&p-J_W@QhL(W@@$i_?>7ZJse z@DJ#EUZm^Q3+9?~SP*4j>Re$W_gd3v$0&GLwK$+kbuP2kH&ScD%`!YlUs46E2yBPBD`*qFszLG((KEEPm&q3X^^ff;_Hs zuKc=2pQjlNW5y*?V4Rg}GDp$2FN>nt({u^c%Zyx>9gk=uqkAy0gr_T8{{2YhR zK~p0Xo_EFM@mi`%hfhZU{lIht4C=wrxq4F(bP2%#4Mr{e35CjTKEzDoMb?Whsxnh# zrurwYhuRp)lqY`f1HTnw!L*{Ou0<=)VFeW9L=6BRrl~DHuX%Xv2ML*M)f8~B8OvgA z$}#_#$b(IgD$76GOXRoM-I(@7pd@B^Z68b~jmyNr&qavgplK1Kgf&RCk`4dG z6z=cXToTj&c*4wy=nYh(#kDfV>CJCmiQl5J&2xBsO%&Udm0?MMC1k`O^gp$+dv#WkFC^L)Z>8Z#$APTcWm;oDV$b79y1w;&7|EXKv?;c z6`yHetdsh&rj5#sZMP$y`+{*CoyQQCr44Lz5#)NDKM>pUz|reyb*hmaHXcMM$OWh| zFCJI1Ds@;JCyF!GwVZZF%T*`uRb{T`6_~-W!kA@ep%pD;)qN*z711RzjGFBsRuH>4 z^b_Pyw)?M4_uPnSLEPeVJ8UbQZm-XdMNrfx5H_zE=TXj&TRxA#Y3nLxB8vFB;oE0< z%n|S)4zeQ|h`BD(Lf*;-yO%{j=pdHGF2wi_YzB$AG|m<=yaxgXzL6}jEo>XajKyLl zZKA2<^8?=88$g%pC@Vtr^Rgy677Vk zo50cgV76M1%hrJ3tmB(vapKjf4R+502e#|5o_6r01e%V$4O$X6LZMN6SYJfu20Pi`_1pNTn) zF+O#zdvuLODU4~j*BQ2mS0mFnx)-&b|0762_u=|27hy-D7CY(mdyg0a?zx;e9&yoL z9tOnm*xRlq|m z5nbidOwH#h+BcE1J*%3$_N7R@Zk{%%Zy05VyG6>zFMc|*aei_R4j)AhX9L+wea^<( zbMSkj*_4AbZIH5SS2cqkPrYtlBEpHNj16Ttv9VJH@z&Dlm9kmHosHM*dKaELJp?tyvGV_E<7{0b3Zf{!RS+zOEej8{S=3|} zgUMtt8U=&Juo%Q(v6wb&H;4$L&7@+p8Egi@Vio-dyl3b9hS`e)>zR*x=HA(T55qg> zO*P8{>fA!~B)NjW?GD)0$kO3?IMXZ9Y59ldgNp)sH*2Z&!xd~q`WJm&Ph zT0T*IA^7G(wldZ7?|&AnEg%0#MFb??1$`PgN_>1KAoWrNV^I3pHb1&`@N)-dLGm(t zk@)xnLG)I{=ks5g&+ms$1=Kql+t#&%=-15^O(;1pi9J??Ee`7)C z=Q2Z2d1C%rD!5}KFvz<0;p6+czb*m0K7YMG%;Q|+7`OTH2la1C;Hws!M%~~{thfB? zBf2dl?sK+)$9OKtnxlUFMt%4xe+SVx@BGVY9{y6AXS){DB{E@bI{(n3s6ouG#)tYj jm+miWpAh}$>A7?s(mz1*5-UySfBfT*5-~`GH0=EWr#Na9 literal 0 HcmV?d00001 diff --git a/Tools/SConscript b/Tools/SConscript index 0556e9f77..e3162de94 100644 --- a/Tools/SConscript +++ b/Tools/SConscript @@ -35,7 +35,7 @@ apps = apps + ' traj2pdb merge-traj center-molecule contact-time perturb-structu apps = apps + ' big-svd kurskew periodic_box area_per_lipid residue-contact-map' apps = apps + ' cross-dist fcontacts serialize-selection transition_contacts fixdcd smooth-traj membrane_map packing_score' apps = apps + ' mops dibmops xtcinfo model-meta-stats verap lipid_survival multi-rmsds rms-overlap' -apps = apps + ' esp_mesh dihedrals rna_suites' +apps = apps + ' esp_mesh dihedrals rna_suites ocf' list = [] diff --git a/Tools/ocf.cpp b/Tools/ocf.cpp new file mode 100644 index 000000000..3f77c899a --- /dev/null +++ b/Tools/ocf.cpp @@ -0,0 +1,324 @@ +#include +#include + +using namespace std; +using namespace loos; + +namespace opts = loos::OptionsFramework; +namespace po = loos::OptionsFramework::po; + +const string fullHelpMessage = + // clang-format off +"SYNOPSIS \n" +" \n" +"This tool is designed to compute the Orientational Correlation Function, in the\n" +" style of its use in polymer chemistry contexts such as in Ullner, M. & \n" +"Woodward, C. E. Orientational Correlation Function and Persistence Lengths of \n" +"Flexible Polyelectrolytes. Macromolecules 35, 1437–1445 (2002) and more \n" +"specifically as it was used in Plumridge, A., Andresen, K. & Pollack, L. \n" +"Visualizing Disordered Single-Stranded RNA: Connecting Sequence, Structure, and\n" +" Electrostatics. J. Am. Chem. Soc. 142, 109–119 (2020). The user may specify \n" +"abstracted 'bond vectors' between links in the polymer chain using single atom \n" +"selectors or group selectors. \n" +" \n" +"DESCRIPTION \n" +" \n" +"This tool uses the definition of the orientational correlation function from \n" +"Ullner & Woodward to estimate how correlated links in a polymer chain are. This\n" +" is done by looking at the normalized projection of the i bond-vector in the \n" +"chain onto the i+o bond-vector in the chain, for all offsets o between 1 and a \n" +"max offset specified by the user (default is -1, or all possible). Each bond \n" +"vector is defined as being a link between a point on a certain residue and a \n" +"point on a neighboring residue. These could literally be a bond vector, if the \n" +"points are atoms bonded to one another, or it could be a 'coarse grained' \n" +"linkage between two monomers in the chain. For example, in the CA default for \n" +"proteins, each residue is being treated as a link, with the link position at \n" +"the alpha carbon; in such a coarsening of the polypeptide chain the vector \n" +"between CAs becomes the chain bond vector. \n" +" \n" +"Thus, in the CA example, it would be the projection of the vector between CA_i \n" +"and CA_i+1, v_i, onto the vector between CA_{i+o} and CA_{i+o+1}. These \n" +"projections are averaged across all possible i for each o requested, then \n" +"reported as a list of correlations as a function of offset. It is also possible\n" +" for anticorrelations to be exhibited by this function--for example, a pretty \n" +"solid antiparallel beta sheet would likely produce vectors that are pointed in \n" +"opposite directions but are nearly coplanar. \n" +" \n" +"The notion here is that stiffer chains have a persistence of orientation, which\n" +" is quantified by the projection of these vectors. Thus, a 'length' is also \n" +"defined; it is the average length of the bond-vectors, multiplied by the \n" +"average correlation between bond vectors summed over all pairs. \n" +" \n" +"The tool writes the results of the requested analysis to stdout as JSON. The \n" +"JSON has the following tags: \"mean ocfs\", \"variance of means\", \"mean \n" +"variances\", \"mean projections summed\", and \"mean bondlength\". The first three \n" +"are all arrays with lengths equal to the number of offsets analyzed. The \"mean \n" +"ocfs\" are the normalized projection vectors averaged over all pairs of bond \n" +"vectors with a given offset, and then across each frame analyzed. The \"variance\n" +" of means\" is the variance in each such mean across the whole trajectory. The \n" +"\"mean variances\" are the variances at each offset, averaged across all analyzed\n" +" frames. Finally, the \"mean projections summed\" and \"mean bondlength\" when \n" +"multiplied together should correspond to l_OCF as given in Plumridge et al. \n" +" \n" +"The idea behind reporting the mean variances and the variances of the mean is \n" +"that the one reports on the variability of each chain projection across the \n" +"trajectory, whereas the other reports on how variable the projections that are \n" +"associated with each offset are within a given frame, on average. If the \n" +"variability within an offset is high on average, but the variability of that \n" +"offset is low across the trajectory, it implies a strong conformational \n" +"preference, or static disorder/glassy behavior. This may mean the chain is not \n" +"really sampling different conformations, even if it is exhibiting correlation \n" +"die-off as a function of offset. \n" +" \n" +"EXAMPLES \n" +" \n" +"ocf model traj > ocf_traj.json \n" +" \n" +"This will look for either alpha carbons or phosphorus atoms within the entire \n" +"model, then use their ordering in the model to draw vectors between each such \n" +"atom and the next one in the chain. It will compute the ocf on each frame in \n" +"the trajectory. It will do this for all possible offsets. \n" +" \n" +"ocf --bond-atom-selection 'name ~= \"C*\'\"' --center-of-mass --residue-centroids\n" +" \ \n" +"model traj > ocf_sugar_carbons_com.json \n" +" \n" +"This will use the centers of mass of atoms matching the regex 'C*\'' (any \n" +"carbon with a single quote at the end, which hopefully amounts to sugar \n" +"carbons) as the points between which to draw bond-vectors for each link in the \n" +"chain. It will then proceed to compute the OCF as normal for these. To use the \n" +"centroids instead of the COM, elide the --center-of-mass flag. \n" +" \n" +"ocf --selection 'resid < 31' --bond-atom-selection 'name ~= \"C*\'\"' --center-\n" +"of-mass --residue-centroids \ \n" +"model traj > ocf_sugar_carbons_com.json \n" +" \n" +"Like the above, but only look for the bond atom selection within the first 30 \n" +"residues in the model. \n" +" \n" +"POTENTIAL COMPLICATIONS \n" +" \n" +"Be careful with selection strings; results that are only subtly wrong could \n" +"emerge from a string that grabs atoms or groups you're not expecting. While \n" +"this is always a good thing to be careful about when analyzing trajectories, \n" +"the peril here (because the selections are being split internally across either\n" +" residues or contiguous sections within the bond atom selection) seems great \n" +"indeed. \n" +" \n" +"The '--group-centroids' flag shouldn't be used unless you're after treating a \n" +"collection of atoms that is trans-residue according to how your model defines \n" +"residues. If you do need this functionality, make sure your model has \n" +"connectivity, or find a way to add it. Using the '--infer-connectivity' flag to\n" +" do this is applying a simple distance cutoff to decide where the chemical \n" +"bonds are in your system from the first frame, which should be treated with \n" +"caution.\n" +; +// clang-format on + +class ToolOptions : public opts::OptionsPackage { +public: + ToolOptions() {} + // clang-format off + void addGeneric(po::options_description& o) { + o.add_options() + ("bond-atom-selection,B", po::value(&bond_atom_selection)-> + default_value("name == 'CA' || name == 'P'"), + "Selection of atoms to compute the OCF across") + ("max-offset,M", po::value(&max_offset)->default_value(-1), + "Consider all |i - j| up to this value. -1 considers all possible offsets.") + ("group-centroids", po::bool_switch(&group_centroids)->default_value(false), + "If thrown, split bond-atom-selection by molecule and compute BVs between centroids.") + ("residue-centroids", po::bool_switch(&residue_centroids)->default_value(false), + "Split bond-atom-selection by residue, then track centroids for bond-vectors.") + ("center-of-mass,c", po::bool_switch(&com)->default_value(false), + "Instead of using centroids, use centers of mass for groups/residues.") + ("infer-connectivity", po::value(&bondlength)->default_value(-1), + "Infer connectivity using provided distance for models lacking this. ALERT: " + "uses hard distance cutoff on first frame of traj to infer connectivity. " + "Only does this for values greater than zero.") + ; + } + // clang-format on + string print() const { + ostringstream oss; + oss << boost::format("bond_atom_selection=%s,max_offset=%d,group_centroids=" + "%b,bondlength=%d,residue_centroids%b,com=%b") % + bond_atom_selection % max_offset % group_centroids % bondlength % + residue_centroids % com; + return (oss.str()); + } + + bool postConditions(po::variables_map &map) { + if (group_centroids && residue_centroids) { + cerr << "ERROR: --group-centroids and --residue-centroids flags are " + "mutually exclusive.\n"; + return (false); + } else if (com && !(group_centroids || residue_centroids)) { + cerr << "ERROR: --center-of-mass must be used with --group-centroids or" + "--residue-centroids.\n"; + return (false); + } else + return (true); + } + string bond_atom_selection; + int max_offset; + bool group_centroids; + bool residue_centroids; + bool com; + float bondlength; +}; +const string indent = " "; + +inline void ag_bond_vectors(AtomicGroup &chain, vector &bond_vectors) { + for (uint i = 0; i < chain.size() - 1; i++) + bond_vectors[i] = chain[i]->coords() - chain[i + 1]->coords(); +} + +inline void centroid_bond_vectors(vector &chain, + vector &bond_vectors) { + for (uint i = 0; i < chain.size() - 1; i++) + bond_vectors[i] = chain[i].centroid() - chain[i + 1].centroid(); +} + +inline void com_bond_vectors(vector &chain, + vector &bond_vectors) { + for (uint i = 0; i < chain.size() - 1; i++) + bond_vectors[i] = chain[i].centerOfMass() - chain[i + 1].centerOfMass(); +} +// this is the work to be done inside the traj loop, that is, per-frame. +inline void compute_ocf_bondlength(uint max_offset, + vector &bond_vectors, + greal &accum_ocfs, vector &mean_ocfs, + vector &var_ocfs, + vector &accum_sq_mean, greal weight, + greal &bl_accum) { + for (uint offset_idx = 0; offset_idx < max_offset; offset_idx++) { + uint offset = offset_idx + 1; + greal accumulated_bvproj = 0; + greal accumulated_square = 0; + greal bvproj = 0; + for (uint i = 0; i < bond_vectors.size() - offset; i++) { + bvproj = bond_vectors[i].uvdot(bond_vectors[i + offset]); + accumulated_bvproj += bvproj; + accumulated_square += bvproj * bvproj; + } + accum_ocfs += accumulated_bvproj * weight; + greal mean_ocf_atoffset = + accumulated_bvproj / (bond_vectors.size() - offset) * weight; + mean_ocfs[offset_idx] += mean_ocf_atoffset; + var_ocfs[offset_idx] += + accumulated_square * weight / (bond_vectors.size() - offset) - + mean_ocf_atoffset * mean_ocf_atoffset; + accum_sq_mean[offset_idx] += mean_ocf_atoffset * mean_ocf_atoffset; + } + for (auto bond : bond_vectors) + bl_accum += bond.length() * weight; +} + +int main(int argc, char *argv[]) { + + // parse the command line options + string hdr = invocationHeader(argc, argv); + opts::BasicOptions *bopts = new opts::BasicOptions(fullHelpMessage); + opts::BasicSelection *sopts = new opts::BasicSelection("all"); + opts::MultiTrajOptions *mtopts = new opts::MultiTrajOptions; + opts::WeightsOptions *wopts = new opts::WeightsOptions; + ToolOptions *topts = new ToolOptions; + + opts::AggregateOptions options; + options.add(bopts).add(sopts).add(mtopts).add(wopts).add(topts); + if (!options.parse(argc, argv)) + exit(-1); + + cout << "# " << hdr << "\n"; + // establish system, and subsystems + AtomicGroup model = mtopts->model; + if (model.hasBonds()) { + } else if (topts->bondlength > 0) + model.findBonds(topts->bondlength); + else + throw(LOOSError( + "Model does not appear to have chemical connectivity, and " + "infer-connectivity has not been set to a positive value.\n")); + if (topts->max_offset == 0) + throw( + LOOSError("You asked for an offset of zero, which is not possible.\n")); + AtomicGroup scope = selectAtoms(model, sopts->selection); + pTraj traj = mtopts->trajectory; + // move unique ptr to Weights into main function ownership for ease of use. + auto weights = move(wopts->pWeights); + // Attach trajectory to weights + weights->addTraj(traj); + // initialize max offset at top level, define either with user input + // or as a function of chain, below. + vector chain; + // Define pointer to the function that will update the bond sites -- will be + // either com_bond_vectors or centroid_bond_vectors, depending on user input + void (*bv_getter)(vector &, vector &); + // determine what points to use for the centers of each link in the chain + // based on user input + if (topts->com) { + bv_getter = &com_bond_vectors; + if (topts->group_centroids) { + chain = scope.splitByMolecule(topts->bond_atom_selection); + } else if (topts->residue_centroids) { + chain = selectAtoms(scope, topts->bond_atom_selection).splitByResidue(); + } + } else { + bv_getter = ¢roid_bond_vectors; + if (topts->group_centroids) { + chain = scope.splitByMolecule(topts->bond_atom_selection); + } else if (topts->residue_centroids) { + chain = selectAtoms(scope, topts->bond_atom_selection).splitByResidue(); + } else { + chain = selectAtoms(scope, topts->bond_atom_selection).splitByMolecule(); + } + } + // Now figure out how many bond vectors and offsets we're tracking + vector bond_vectors(chain.size() - 1, 0); + vector mean_ocfs(chain.size() - 2, 0); + vector var_ocfs(chain.size() - 2, 0); + vector accum_sq_mean(chain.size() - 2, 0); + greal accum_ocf = 0; + greal bondlength = 0; + uint max_offset; + if (topts->max_offset > 0) + max_offset = topts->max_offset; + else if (topts->max_offset < 0) + max_offset = bond_vectors.size() - 1; + // loop over trajectory + for (auto frame_index : mtopts->frameList()) { + traj->readFrame(frame_index); + traj->updateGroupCoords(scope); + // get frame weights; defaults to zero + const double weight = weights->get(); + weights->accumulate(); + bv_getter(chain, bond_vectors); + compute_ocf_bondlength(max_offset, bond_vectors, accum_ocf, mean_ocfs, + var_ocfs, accum_sq_mean, weight, bondlength); + } + bondlength /= bond_vectors.size(); + // create the JSON report, written to stdout. + cout << "{\n" + indent + "\"mean ocfs\": [\n"; + for (auto i : mean_ocfs) + cout << indent + indent << i / weights->totalWeight() << ",\n"; + cout << indent + "],\n"; + cout << indent + "\"variance of means\": [\n"; + greal mean_ocf; + for (uint i = 0; i < mean_ocfs.size(); i++) { + mean_ocf = mean_ocfs[i] / weights->totalWeight(); + cout << indent + indent + << accum_sq_mean[i] / weights->totalWeight() - mean_ocf * mean_ocf + << ",\n"; + } + cout << indent + "],\n"; + cout << indent + "\"mean variances\": [\n"; + for (auto i : var_ocfs) + cout << indent + indent << i / weights->totalWeight() << ",\n"; + cout << indent + "],\n"; + cout << indent + "\"mean projections summed\": " + << accum_ocf / weights->totalWeight() << ",\n"; + cout << indent + "\"mean bondlength\": " + << bondlength / weights->totalWeight() << "\n"; + cout << "}\n"; +} diff --git a/src/AtomicGroup.cpp b/src/AtomicGroup.cpp index f80827dbe..12f19edee 100644 --- a/src/AtomicGroup.cpp +++ b/src/AtomicGroup.cpp @@ -843,6 +843,22 @@ namespace loos { } + // this function takes a bond pair offset for the whole AG. + // Returns the unit vector projection across all such bond pairs. + const greal AtomicGroup::ocf(uint offset){ + greal part_ocf = 0; + GCoord bv1, bv2; + for (auto i = 0; i < atoms.size() - offset - 1; i++){ + // compute bond vector between atom and its next neighbor. + bv1 = atoms[i]->coords() - atoms[i+1]->coords(); + // compute bond vector between offset atom and its next neighbor. + bv2 = atoms[i + offset]->coords() - atoms[i + offset + 1]->coords(); + // accumulate dot product of unit vectors + part_ocf += bv1.uvdot(bv2); + } + return(part_ocf / (atoms.size() - offset - 1)); + } + void AtomicGroup::reimage() { if (!(isPeriodic())) throw(LOOSError("trying to reimage a non-periodic group")); diff --git a/src/AtomicGroup.hpp b/src/AtomicGroup.hpp index 132d52835..77d36b035 100644 --- a/src/AtomicGroup.hpp +++ b/src/AtomicGroup.hpp @@ -434,6 +434,9 @@ namespace loos { box.box(GCoord(x,y,z)); } + //! compute OCF for all atom-pairs in AG of distance offset from one another + const greal ocf(uint offset); + //! Provide access to the underlying shared periodic box... loos::SharedPeriodicBox sharedPeriodicBox() const { return(box); } diff --git a/src/Coord.hpp b/src/Coord.hpp index 6890e175b..a1437c8b9 100644 --- a/src/Coord.hpp +++ b/src/Coord.hpp @@ -349,6 +349,24 @@ namespace loos { return(s); } + //! Unit vector dot product (cosine of projective angle) + T uvdot(const Coord& rhs) const { + // define quantities to accumulate with first step of loop. + T s = v[0] * rhs.v[0]; + T d = v[0] * v[0]; + T drhs = rhs.v[0] * rhs.v[0]; + + for (uint i = 1; i < MAXCOORD; ++i){ + // dot product + s += v[i] * rhs.v[i]; + // length2 of this + d += v[i] * v[i]; + // length2 of rhs + drhs += rhs.v[i] * rhs.v[i]; + } + return(s / sqrt(d * drhs)); + } + T operator*(const Coordrhs) const {