From 12fe614ddf5854f1b18109a0d77fbca9cb8cc7bc Mon Sep 17 00:00:00 2001 From: Oliver Henrich Date: Tue, 5 Sep 2017 10:48:38 +0100 Subject: [PATCH] Added sequence-dependent stacking and base-pairing strength --- doc/src/PDF/USER-CGDNA-overview.pdf | Bin 3429502 -> 3428837 bytes doc/src/pair_oxdna.txt | 36 +++++--- doc/src/pair_oxdna2.txt | 29 ++++--- .../examples/oxDNA/duplex1/input.duplex1 | 8 +- .../oxDNA/duplex1/log.24Mar17.duplex1.g++.1 | 8 +- .../oxDNA/duplex1/log.24Mar17.duplex1.g++.4 | 8 +- .../examples/oxDNA/duplex2/input.duplex2 | 8 +- .../oxDNA/duplex2/log.24Mar17.duplex2.g++.1 | 8 +- .../oxDNA/duplex2/log.24Mar17.duplex2.g++.4 | 8 +- .../examples/oxDNA2/duplex1/input.duplex1 | 8 +- .../oxDNA2/duplex1/log.24Mar17.duplex1.g++.1 | 8 +- .../oxDNA2/duplex1/log.24Mar17.duplex1.g++.4 | 8 +- .../examples/oxDNA2/duplex2/input.duplex2 | 8 +- .../oxDNA2/duplex2/log.24Mar17.duplex2.g++.1 | 8 +- .../oxDNA2/duplex2/log.24Mar17.duplex2.g++.4 | 8 +- src/USER-CGDNA/pair_oxdna_hbond.cpp | 80 ++++++++++++------ src/USER-CGDNA/pair_oxdna_hbond.h | 2 + src/USER-CGDNA/pair_oxdna_stk.cpp | 72 +++++++++++----- src/USER-CGDNA/pair_oxdna_stk.h | 2 + 19 files changed, 196 insertions(+), 121 deletions(-) diff --git a/doc/src/PDF/USER-CGDNA-overview.pdf b/doc/src/PDF/USER-CGDNA-overview.pdf index cada7eaebf80121698d5ef4812c04f0864b77d98..e329877bc2cc5953574ea8fbd65789158a050879 100644 GIT binary patch delta 16790 zcmagkQ*dVCwl3fp9ox2Tr(@f;ZT&&VwyjRbwr$(C?VWYjK2_)LRL$pNR(<#HsCSI7 zzq7!Hg|oo^4ZyCEpd1|RsrhTnR0x;7z;Hl7Kz(^%W%lqKdCj$vz&9Y^{_3DMn4zih zz&}8Kfc^mc0p6S%57C|g#>~P_n0h6h=ho)4k|>A>76u6%7opa!+VJ*D?X=yk zuiyDbl}LjgteI51_~@2gLKHc#N?fQQRoh0d-VNz#4}UIDy{-4@;rn$Qjvv!QZ;zO5CS(t_ zX(_&Z-PrF_feoqp?b_1UP)YyxW(cS>xj}Z+bZD-Q=*+cyTUywRO4MV)+IVhrB%W9r zeFBTFa(l;NX<)4=ue$DgBu7_t=#(=I$Om`bs?x;1V6{RZ-(*A@<8+Pz$q;$9r3l|< z-?=LxT~I;XsU!rkcQJ&|_#09(8+_H)Q7Ncgz&Z*4C@ZRP-cgZ(w_EwDumo^exFk>{ zg<#=~>A7ZPRgHg6wiLtx{^end-n4yLnc5n);ReikWV?ORR2aV~?bW_^`tMP>)9xOG zjBBs`Ulnwl<3i)v;2FMwQ74*%PB@(N7- zGlvjrCdP=8`xR!A)ZVIQ z{MIOJ00A1RVMPI2P(Pu~VH@(r`M_!n2h-wi@ArDRjL?%nM`aoTLQ7|Q<4WP1-{T0u zN`n@T2(nYX$blbfDbk@)gul$;Tdobs{Fc-X-@l71VxgPW9=?W z1`%tvsLoX~5M_+r(@~G31uGr;&{GmOlzp8Y3i* z^bG+I^*W5A=tK7gzWdToMAzJxzdX0v-emDIieSU_g|bZ953>T$R4>YV5TyHd(;SYd z(qs-iy>GAReDmFV>1~MqvRmRx@t?C?KFxG?8S0IqssjO%KQ#)vSE0iU4p@5&`x>y7 zz^*Jt&1SL%A`YL2(DBKA)XZw-OF*0^0<2Jb*nRsL{ZlDxTE~N{kZBrNf14 zw6y>kNM8bI1wZ>|#bhfe(SZM9MQCMqQCkYH@pX7k>LITR=s81verR!?ytzgpZ;rS> z1H~Y_8)O$w22n0tu%28t0d4+qeAk{EuPVCfE391sABDYB;C^CqAV&`ja1IM*nrJS` zZXBAIn{b1WHmo3(7@GneyPSbw#uo>hsWLyL@^b^A_*`-R)PokO!njWUcg_TrS!Dy5nH> zkaGgy@H|6ddQryfB)Cz9O1G!!vfP+GH+aoRIC>I2`m2*Jk`P=F*scBZG0kJc(Pp}f z>|(Gh-LLxQ#Jj4rrrirxW|Y~tSt`+*dtaUFJ3u#87y2GuT-FdJ{bfBYT0x8~0mtl{ zGG7UD;0@0s@w)KroG5EVbn7p5P$)S%B;g6*r5b>8WL6{p4SjlJT~W;W=h5_y_W;d1 z8l!TM-=fShT=cCWs}$%vya z)mi91J=}bO_zY)X(_gvwAI2@>4f2K3E-%-fxttSsDU6(y7)^wS8&F#8&HSjwO>8q@ z7349|W=QvMZF5xmL`#BjiDp7Lmgxn@mr4Y`iE0Ll;7m;vJ)A`uHjM z%cpD>+c^8aHSxt0U-y9((qVSK%hIb(_r0y;Vq(rj1-O9r!+-y{#zGJ{;|Fo4%xiT0&je&b8)$}CEK z0bxjXg%nys)Fh#tVbp`c_H8Xq=gM#dg?ne#qi&~ToNG**;8zN)!uB*BtxqWjBff_P z4LKVL;QUyJZXS#7z902OTvBrKBHDW*b0J}Q*Poz|nNI0}{v?2DPYg)iWrajZ?Wco8 z1LI_6`QQ1bQA0m=pB)~sEg#>^ZLznO@?G4+FJtc>fxamk0MQwG&PUlqjx<%PG`^L7 zQOZE2^hCPfDm&|jKZGi(R3PysLZ)?rF}wH9fduGY?|MDrAIEUSj$oBXUk?>U3KdPK zi=h2gFGJh(xF;ZlO+N05E2?-S_vF+nS4Y{h98C<6YGOm`{_Dqy)T(Qj(358Nx_ueD zy>Yxs9MdP;5bhXkQ7dPFT7P>9e0gB0ka(#B)-9?y$Joql%iwBlNtziFCKn)>0G@l3 zYvZld;LD-%IU||EJt~K492{Y~O+yN1RbAs9-nQzeQH*%%m%(|t1q z%NGx*pb2#2yiy6JFOaW=s(1O#pAN~W*f-?m++e^NEv@w1%wDd9>76BJPSGWcKSNCQ zccI>i=BH+7!G5R4P=IB4Ba8f05BW^xm7(x{{?0_HX%_fDrsS)~mWVSNhxRaKj99Ex zzqa?A-DCY8;dqEQs3*B-Zt5_&9K71irki)b!_~f$NiQBeIIh<`ChUA+vF92)VpHfw zyzb>;8_6vh%+i`#P~*QNQ*C72$5Os-3D4F&oKzT?gBRHq$}gR) zwN;w+^4V?FB$=QXTkTp*h~qAR7SRL#@y=o*4q2@`+iY;NDAGaWQ#(I|Z)eA=Lhx9k z$>tQ17a6)6CTxPsWz8=e3oDUM(9DumDBMM8pb9QXdt$*B%|`|WYG`;+7MQpaj%6$A zc|}nNemO#W+FSU?Tbq^R_@kgV=#UHW>ZpRtvuwFf*t$LnG*q3TG2tVu*@6KiQ_5uD z7m*dCTr2TncwH5;wRGuNhj3t$+;>?pE4j@>N}d`|BAf83RYP6YZRqVXIT1wZd)9d; zL?pNNYL6WE+Ovp5V6XOUuvxq7i$RTjjn!_`2cvd6$2$Fo+&Zi|_#gu!!Nkt|-lt zr}R!gD3bWSi>HfR<`uWtK6_!kLqJzuq?XKM$?Kkn3)ahn7B4LgYVtB2T%jpp|Nb?4 zak_#nc%_hV(6d2MkM}Wrh1x09i#)C;ZsnpAh!f%;B84W|W^~4mV4RAuDG^IpKJvvK zNClg-4p_g}HrYc)iQ5Itp+o03IS?%JKF*MklB|0B_WANYu^2LMG_aJW0O6G~+2rt$ z46UC}Mv&4b^hE;q9n3U$xi@EwB;@VKjMMS)LD9SY*{S*3#e+kf*i)@TW_5`sf!Wur zN@Rr|@95xFI)f(Ny510eS3@w2&9t;@<&$ft$0ZD*)MD2wW9vx!_-bX> zhU@s|pgd0lVRH?8fImYB<1pR2oa=+GA*%|B-%XU*95Y*Vtm_y{9+Dqy&F?^Gw0P}F z#@YhKe7;IO9rd+(`&~~hzyW~#Z{&UPM4OrS%BRbVkI97NQ7!*1A+apoENu=*ayg1X z+3)h^V1AoNm)`*niDmnlCkxHVmqjs?i8p1ki2tz23}9eT{Z@P-HcWxyFnt2wlHIw9 z{{_^eUX0F+=m}<6+ zQ0dw0p$nINn5|lZ6;oL%?qYDxvYTE*EJTY*%YOCD52h#k^SGH5T z>07$jg!DF*?hLS5$)m}GDHChEq%dMm>o8Ji@;+fR({G-R`R3)9u&3(PA5PO}8H*E!oF^Bz1itj^lI4@DFCFLX8-FynMw)9$;I)1 z<5X+ZXR#B&Z{V+$-IVr46A80a1|{2?7`!NV}lhkC~=i&>w`% zb;X!KL3%0We7xwzhHtJVBl1DUA>~eGBBi-rqlR;~Uc;Kf+g|o??UBwma*S0szPG7wvvTY>H}yNfo&Jvz zx6CP>s49>qG;0Q={o-@unW$?p*BmJy-AoExkuetE!ae#|m1a}^B{mx#Q zGnsT_%0_}*R{7U#D|DhmmETlrhNvOx*lsJdGE?^0N1TymK`tgwg?afB5BVu1-{MVw|nh+Ck~0gHs)bpa+8CrUvFD@73tiOh*N*_ zG9OWK_HjfLb7cy4H6R&DJN)LA`n4`BBHLU24bN~HYR5C+quP9%20H_b0kPQYCnN+c z2=Yb(RuW~yZx>V-6d&DI1}9m!yr=1l=KA8i9WFdJ5y>bQ)Xwu9C8=!~x{9)+HWR%V<*s_B3~grmzMhKP z#@^d#ho}T#5QPf;5Had70B{bSr)fU__>K0&^|Hj>jks`;7_L@ldy&~nxWCwp;OhDv zx_UmHhz_{9;pq~@h{~KW%5hd{6Kd;`vAM6bp^}Kipt)@A^+oIHM$Bq6)aC^Aex*BJ z1&vgY)9CA+f0nsuMEczz+a#^Ku4EgTE84FuNT=H_b~{y2>m`Z>00(;1Ax&}n-ka!$ zU_2_=^_Vg5NW&riKrN{I#^pBFB}syT2gr^KzzzAqa~Z+J`0v0iWxk9ZK#p7uS4|#b!5?E(gp2}LUO~V%wnL|;_+O}@ivAD=k%VC5_CLG*|8mJck<2A= z-}X1exbbc?gWF0yxEjMj7E^PGOO6%6f}fi+z4VV_f*QR$7L`q=$jg_=n+b_sh5wv- zZC^g=PK(n}Vo{+x^Pw*YR+1oXu!$~X9RU8stYd@2nm^Qb4L8j&O0i%?CPO^Ye5Y8t z2Qy3glM!jq0zx0R3%`^GWh4{;cy*-F_|EDh{;0wy2wdpARVuSDEA>~@;ICpi6r>tE zMzf6rmhDsU@d6Q9|M@@?nw%B6kxiH+>K)7dy50^`zx{5Me)lH!n2z*u;)`k}Y>-BQ z@W>VL*Te`<4`lYNlREB|&8nddi@v8!mX8)wiV*PbS;or{BLj6!$vzL#!7WByqi=!X*Xrxgxc z`RHj}&yaB~h)7n6N+p1j#hHtyZzdfbCXO_?PcCAo8`P&tzkJ`FJ-#0<*4M*>*o>f5 z&M*UR_n#YVJrrxVTo&$ePfBS~@}sA+m7LEUxgM6>t~XRf54k~?WkLV`c_A|-e3gxk z4CT>a+A_zCpw!h9>EogTP?~Y0th{^ZARO+{1nO?0Q^w%exGCpoLj_uHoOc{4=M&?4 z{d=dHhHJumc6+~e_I`lGD|wDHd^y1X+@>w*5!elE63kJqwc3m&z?W+i2Wjvf7 zKCF)W4&#>+IeqjfhBYAgoKY{C1VerqDdbsBpl?8HfXPO&x|{UpwV}PL6%W&XDa~pDg#lRrKgd3q1G4iro;W6GE8}RCFTP} zYIwz=)t`#M?8PeF7QK6?X^sB5)tF97p#`&jty=iC!*)c{lByyngQJZ8Jnh`hZThxw z*``E#rhu&3+Xwlqvx*&e>E&s5`f;8k<^M+8)0@i(W(Fr8Nb}kE4Z7JT{&WW4EpA7@=kkSFpHApjafC-^8-k zVIF+pkB6t^Z;*-9XM(Pc!KX2~Vw!n~RiuaQK|BPSsag8>Qj5J*m0`5ZFk#swKEo&h zg?(;jH(hZd!-5guIFQa$=-^U38k1bKOVaLTb3kr={wvSFya&jroys@9xQY?5Qu*vl zZMrGHlueUsFgd{qfFLZTjKrBBy*FR(zT*Sid?)Z{1`=7r7z0k}VeJ8rt$-8vXe;S;0@CUyX4jfNp`<6r( z#>bC5HrnRIO3CZtR9#5zPbvURsFn?YW7S7}Y!#5At2VO(TXE~|D}9kU`;??~5(&g;Pk z9psI9oHbP_H78;bo5>ljy=NW1FRPDBc9Uqrdg2i*ipUd3$q>dz9TWl1a-$j3d_mwp ziI5A6pB` zo9ZtD)~JP&f>>UJ9tG8mAc_+pI}mWnxv9cRfC{fi&f zFUO87dt-=*^B6P$VV(_;`#pZf$|l>iWE+kclM~s2r3( zJpu9yN_4d&uoQC6=mXATLMLcRe;*vJmJZG5M(GJRBQ~4ULd67Vf>33FK<4{X02yD% zKKc;ti1tovaM%CSQ;9s}g%A2FtrCcEjc3qXY;E6S(^$kNcGzrCj_TP?+aWX`0PI#X zH-%sa!8JA2!8+lI0cHJzFjJdD+?1G$-axRgX{c>(&0h!71kT)H29WTa?105>rw?W4 zApF8M?YyU0cYp=_vHu&Pcu+3e^yb8cMh^hpSwMYVqj*)}hYBJS4j6I?y@2^NFwdP+ zN7phAyFze|v=qR-XZpH`1WI@(7R0o%Mu4qN5k7%(5AYWO%N|7IG&vX2JP&12-&Rde~XtEp?4?0>J`{JU-^9zOd1uXAk(*EVCnB#I6MZmCZj^e?O=lgMI z1>}*MXOzVBggllHuf`zR{44adOYPUC{+l6oB+2iXZcz+G7+%b?>{(Du z1ZHU4WJK6MXBjbfM?0c5H1%kEd?GYo42~o1Q19ZX=uASC62ws8%%;$sZmHpLsX#NY z1gz1Ium3Hw@B|^zL0MS;ckmw6Q2bAi<(I1;FAizJn?Rd(^f%UI<*3Cl)BZ{WvwLlOvzy7PU(+US zH}lgl(Ud85Rg4CmhK&8)>E7AJfuN~9sWxmrP8|*D(c3%MJvGtx67ZDzDR8}A#G=WQ zGenH$cMzCZBfgv)=IM8XeNBL`!~PkG;QHv6N~*b#ON&(1rQ=P!&77jbzI6pAUGO>T z%4%!9Zs&=PCvX)}a%Jybj~*de!W3D<*SI7|%TA0kP)Us2A*Qq@2T=A_I?G&Ov2V zDIT+u2~#?@l9p7F60-DUkttDsVvRe2GQ5i{22}NFfZ(ub3HVo#`^b>nn>G5ux4X7K zibi1sfJN5!uA!?u!jtWkMZP`T)*dbRTq2#ddsHw;5NbOYgEN?Zw*^VPJ5Q+}Noly> z-FH|j9)lq;eQ9mk>s_guK4-RK0jVEoA-h^hX&+voJDkC8|7W@qHS*9Wl?a}v3hQ#(mtvFO;CsI)|kjeltCm*K2i$4 ze7xPr3L?Ll23~&fba7B!oG_ezCGy0oKzshC>HmYylrS&OW)^h#;sprJy0PMbu>Wn7 zpfXI2l)2~gd_}LOBWaWt=V!TsPJZ?6n3AbJxe3_e06-CFxCuSj5wDjPx9(*OYG#1c zm@v&#S$+@URnv7TCyVXfU((NG=eQ?7qpTDxwpmVl9&&&0?RZ0UxmY^6S*Ax|;bLq& zq);_M7&6QK!A#SO#pPztfWu5BRvRpnfXrL(owrly_4gc_kUOe{ke724G^xagFFG)D zp1Ni;17O-jFU-oj@--X7#xMv)AdGxW9?%jpNgUoP zsh@`Z0JALv@;Yln&NDPuQCekMaY2}%t@H|P!<~Kqv@Eacb9!3*2&hM_n=+W{e%X0tBfA;5B*X}6%ntwak383jE~~l2-`G5rZNLQ9 ze4kb)Bv6uRp`vlwvrPE}R%oBig&X<~{9golwpkwwMFkkZA>-aX47S%zAXDisGvcyz`g6dO1hcM>dvIeDzI&>W`#Db3C_LC* z9akZpr10gf$`l;XuCa_;(R7A@AYf<*fxh-0e)!Njwmn!70fO0L^s?RB<0=YU5LcXL zbtVMVu^t{daCnLUQ~_MG71Fz>cFAq|-?c>wKE*FD5gRy)4TX^OA_L8Ab*magqklpyvffMboHy{v|;6*+eu|DbYuv;J;t2h(G)FlW{x5 zBSVBGoiwV9&^EA|+QiiI1L2=-@f ze))Byg3IMuf7aCV4rw0ZrfT>M-$D|0hX9gTx_! zK>dLJ@#_c757-}YKj42LGzW0D94C7zQ4*eejSMR~u9LiNltG z_9#4cc<5Kl+sbajG7|pjLfpZB*Due`yw1-`xuJ&A4I@~_8w{bHL!hUS2f@BbyE32? z8F$mdW(Q=XMvDyQTKv6Y3L{@bB*&BbB>Iw)2NazHkV??7;jS{MhFzk_&A5BL0C$KO0BRGQqYY;8I)u^}jHANXwEt>57$NG3Fd&&kX2+QOiSSNf!CZyGG=6giEzXke zfCYq0<+(bRYa1Dxp;!P7R{V}Mr=>v?FHc)AIAyd8{zXeAN}FFL3)?R|ClreLK&l({ z3{5^bqZ6p#H!B?Q=OJVP8cfju#G1O$1}Hp<36K(y%$x>>^Y_o6eZ~+{JT6#Lq8NWx zt%P3BKzK}OA@V>+@S<5MMMyJ{cR?s-<_G`q8Ann^&-wV-kv)e%LXx2f3*+a<*N0mxQHj?k09bWcn~X(j*gaBPfP*ObEn;mw}j%0C>`aSU2GBYU&5=pRuG*#5eFrxC&1N4Rv*eH!g^c+MrCh;6Z^Z*i~ zBLDL&X33i`zDRlv7Y)4>+BI|S#_5a>T7Yovq~9CW3O9YmKisLc+RNyHI+CuA+lr3$ zZ_^&R@q7aN5_P;!{J$%HdK{`>reX{lTCl5nT{8um#TgsrA(N zY{g}62MY-IF3larxm;w&1Apc;Z2$~NL!+ud{Aid49fkO)qR?)+RL*l_VV_O4y8CV9 zR~c^brm|L5f#P?KQGq8$PsNtG@^^dFgllx=`PEvI1F9~lc;en3{zdAfUQ0F8)JXXS zsHCnV@yM#kp>Y^o)U160Ha%)^8Hdpb(_rB+5aPKT_JcQ8E0d|~9|TKLBL%tNfw#L))7-$Ei$ z!SN`KFeqT*N#rh^NRDx#V|>PGvGun#W1%G9Ln*UjAt7LKybaLVkGw_BAmU6Q{HwGI zo*cAu*YQD5nhbtuduqV;D*#3a1*H6;+Yn6la6OuSSN;rD{)suuYeWwmNh1x?9GHYS zwheVC;`hJ@=GjufSHKBEyuH!Xu9e!v(0;k>p9H9E(V+S#Qz=`NKeyH%zk&vVHSrbk z32~ad@G!YwkIUwX&sRMOET>uqTYAH{pP2ZB`=%V~2=7JB-yvz+9sur>R&v^-ty5J; zI;GuY$M|Yw@T~0_0k#x>6LkyUG@H;Y_?gyl$YM%4 zGN{}w5`9Tz$d(|Ee*lJMEugiCuYqq2rCZ|iC>=;`3ajv0tIvxIA3xH% zv^iv1DgjXiw3#5LQG~je{Z_aonw;+}*t+%tYr#_;P8j^eHttc$$PixO`C#IpsNIIv zbj_3dmkZlYDw_lHO$4^}@&w%SikhVz%5H2@Uj#x|6Fdm({)vMH^P zOPVIvjBlS)f!h9h6gLm$CE0A%*9Jq$kxBNHYM)3yV0pgvwNlW@UvC9I`)7C<0uEoG zCqR?>#H|xJdFymG2>^Mhkx;0`!C$6}VpX_r+E1C)P*lulQf?{itCK6+Z%iy$VUYDl zAr{hEH30jPR{-gx!eGcgv-I$VJ|)_iq*i_C{j>i1daUi+ZWUJtP)?|!Ku$Jutw($2 zQkSbhQPzvsF|JRHnALv?46Ms zgYG7FaKvAFCbu_5)pm_SII@3s>j9A6is|byv=`!sR{$r2UFLtj?S$z(Cx$!bfNu;X z!A-Tn&Ap)d#O^Yo>tL-v6UQeW_OVf{CN6OPv*qd*%CrO_mxL(-g<8nj>(`oAqWy8VZo4<$ z@22nCOF&uU3n!JnZ2LzVd-9Cy#GWO+=y<+?MT`a~03LptO7>e$Fn6@(Go7v=woo1VDT?Gs@ z+ZY{ zDqsY)mAKPdH*^Jkv(1G`s5!F=5YkaJ#uCt^Mims~Bu@0W>jHtiyJ4pS>uM;^Bx&`; zTJMyvS%8$tZyTgOo|CcWN-7>T!PtW5`pf#)I^?5dO!1dX5C+4;GH-mqLBx2ELs{So zHcFZXMSudtp@KowUI3_!l;F7xVP#wZBOn|O*WD{X#^_b5;y%nl&yh4=c@Hf~Dz8Uh zIg__IiHR{NTpu*@?ZzBBGN+FKYHc$+PkuRg`={<{D94nbeAHDJbSacVo((>8-@K_H zbXO^2CUOQIoT00NLMhX(<;Ddow$iUnF2M5x|3K*mZ_1rtvk~ z_;pxdASBE3Zmmdj%u-gIOGGd}&j9Ro2%4J;!dLX=tp^*ytq<~O3hZz&3~25et$(1} zYc}KS(JtHe$lYAMI=#I*ZI65(O?BG9*!2k88P9yrE(BM@De+e_oej=w)WyTnP3cDW z6zg7tZig^<`Zpd4+%7lF_9%z7`J~j_ib-7>si*!&tUVODY(T%Xo(c8E^}6HR!)4XG z<)YgQeUsk5ao=~~RAeKxkd}EB%u&*^qx2hd=Rp9}Wwi-xf^WE@XA6FE=oYyb+kc8H zaFEgDi-3cP=HtJEuAVVssGc#>59A*xKTv<5{XqYL@dNV*)(`9-I6rWI;Qhe=LGXj{ z2hk7WA0$6We~|qk|3UGC@(0xq>K`;eXn)ZCp#Q<}gYgH`59Ve)V;0F(P)@f0?XNDi zG@aKuoBWomfA%WL%t9ACD$6RK;VACIwd-!~+>)(!v&gpA1(jsfDAU#5+rX0H3+U92 zo2A4XkIS-!wT6{bum?1sd>$`msB;+UgRaURrW(Zlmh1jZYr6V_+T6W;xQ`WZ^# zDbf(-+-{78p@Z>54}mDMN?mLJ&I1|{uLp{ z{e|N#SdnXb7-dx?<0|$+$Vx!T^JK>VgdmqB7Z625ge8OY015ZA*vf3gg8utKX1k|z4=_w)`4swwEpBmhIJ8I(vcN^5}A8%Ts)vSKTMl`%?FSwXQP z$NgoL4x%n&gkdc(%p)AI2nOgTrWxQmV|$ShS(av>3wlj>3Mqk-&w+CZN7vLZgm^-- z7KH?o<%O}!1@=T4gcE5RR}KV~&c)S1Ls}rK*>v`-5lRjC7t`mXoiduxz?Q|7@cT0< z3GMdU3yr~!{*K1=iz7BG&_*U=g)o^7YMrxKCgMCbQc|a-wFIXsV*uo<`6FnFIZ&q* zu`l>L1>!M+B36vxOV|>!uu4=0)F}L_)S=(mS5xx zjF-w^3>w<|q;_81L;@@sGM~f|E(SFja1K8|PRWQE(_!m~wlp_hm%F(?4nSXUhw;7w ze9cT!a`Vr>o>Cw(R~T+jZu&J-Qygl6=uZ!9mjE6`Xk5UvBLOnJGXzBlEG019F)>5i zHE#hmF(JLvm`hN$t#m`%rTsS1z~-LZC+-O12tiRP-4LSyB|s2)))OSwYcc0dQNZ=C z1#4CPyDs1X+y);-njfX3<0HhvY+IqzER!l7RrxtW88F)Q2aA8iAeHn*kN6JOK0(!g#y@fuPW^n>nx- z+)+^p;o`t}D2?nRa|=RX5W%iFrH|qMO&&5sZ~pt258&45!ioDz&)fhV_3a?(HuFMz zO3mhs*>}W%PY@M?1@zs z3bfyjl;2=#5{XHO@i*(8kBcwX-$ax5{9M+bK7bQKsuG(vIjv|^vKJdB;ja$z!{+hgT zbpv)_gv)rXUjp_mQxd}D=Xi}7K1gD|14y5IHs0}}B(5ggS5Z$H z$vSM9$A=TQKF1#$no9aP<~_wg5v>8-tC(n1{1-~JPY`>JH zf=i*TWX%`0K#L+5I&h9k3~wgASa*t&i6yX>EIZZWC0TGcDsA;}g|S`BU} z$^|4*DJ(LGei}67g?@ZzG2SXB{&a=?vlaTT3y~Skt2@>y&+VsLLciJGweG&&w7-tH zL;Z;|e!Jzun=NpgJEH!^#7I?L_E2i1r=P|mb$$@*Gcyv-HMbUT~oSk;pOD&>pc%`n#00Kkl6~Ca#32x3FOC!J(2(Js8ebP zhp{80pGjA2Z<+Z+yDyB)MYiTIRSqCQ{ib8?@(No~mFDDTeowcNi*T)`^Bb{Vro^aR zrDCs3jCs@Yw+e_V76s{?VwZ=EQUxng7E$kn;s(Euo>aXYL_2tS&7HQA(wR8yolQ>S1rD};J1fP&Ac8&OpWwKMOgTuO(iy;Pc`)rqTfj> zkSd5C$|hNe9@h&J8cHk_a+01F4+QnR?lN@$fqOAb73r^jB~Uki@V-quP6RmdON921 z(qai5gqvtH{URclF&1YYLPtVWt;)%q*I^5J*g4A8Mih|CCkK_}( z%k^DHY}OQWlnx>pt@)2}-ZgFhgZh+e3~0~`t)F4ftFf_M53mZA+5|L%6juy1Xf5_M zLa$n>?(3K&xaE8WDN<2`Kp8#e{l#GhOIY)@3G;RUViW9#jk<03I(TOYl>YBvhK>Juq zV$Hw5y@v=|w9w)oIsu_$5_4eGBxSlu!}%=@sP7gs)_7iB>T2^$p{%0X!_s(A!H%|{ z|6JgwqV*OZ8jy=vOsL7059oobr~0 zs^hs4-6qThCrQgzH+1jO^SO?$wk4g-mCvh5R-aW%n+oO1W6zb}&~eEaXI>yq!?~UO zD=cbq!(^l-GX_Ly8le5Twiat4TpkpllhQK2g`uDxudfT9jZ;hrlPab_M-v(dtql&C z^#)hRHucRGl<=KxBg1CMs2eyS9KX}N)}NqSk2G#ti`Sv52H${??$~scY2;Pre!F}) zyE|IEeSJ+2sx5*2qepgKC>6vwnOxCx>hUHL7VMVGV| zS_o~kVM0H`2=13Z;HbOZKeF_1r2T?3E(?HPY4QYXzVll zSP0hl-hSKP8Mp7ZK4SzgBmkcbUNe%df+0U%5=V_thUbfW7a-@5z$=63Pqyf)?(!mxtSF^<}ORsnX8# zQowFIozsnp>w_yz2~SAv-}LNo^@%CAn5oYvX?_^Vna`qPQn@;zoI4K6n0azrKL!?_ zcCOBsmAuOxueg}IZ=pONJL<)VK0H#K5gWu;!Gl}gYz_%lZ?5se#QtZD(qb~*=+Nuo z(RNU7<2B)&#|uHCewr@Tm5_6rw&WtmMF5Rl(})6VM+DmWlS`j`zvM5?h4&2E;^Nzi z%t`ifstU6~gN&5AS+R_d&2i=LujTW40eQ8T_I^zyKymd8jcauG6vvCTNOpIh-8;XD z)dkxv`BN7}&^8$Q_ZUFin5DG}0q)p!_`1aMSvn=uU#;o4#QNDp1f@peqyfTk6M)Dq z;U?@(DhzPneZ0?mz6l3n@Tk-i(pHvAP6K&h0yateWK_x9daVYjUo=%-ij@X1-iB018iHxn6gpR8 z=tt{VoRo4`mWd`=ECF}dmN0*qi1*~&{w@sm75SF>vAAb+Xzj;|*YZz4m$A#uw~yXl zUSY0jcERnkE8(Nsu1bxg^rYU?x$*m`-17D8%FOk7=%}&t@XOm?He5a4pFLf^{#q0CM0*5kEv zF{vJK{)`y!?~LBa^|`2xRnc!`t1)z? zN&HpT@>oX|ICHN~#8W;0?W=?N*(e7A5z4$+P~%lH^PqQ(CC6#^@Uobn7_U5Vh^V}% zL9kWjb=>kPxy-Btu_Q(4Psy3C^7uOA^=aYKn5S}Wj8+B8wTAfH8h7T6rcavrqWD+z zIh3}(xSa;_Fzr%&wClfgaBdi>FO07pM)*{wpS4pet_-U&FlVa23@ZsV7b7bZCo>l- zYiflIEB62MJ1xViLCHx-sbZ=@C}ZhnO32K}#6d-9?O^C);cRM1=w|BVY-w*t$WG5p z&(bU}%i01e|65Fm>$fPIh%gHyyD$gWZ%$TbRz`LqMj;kqPEjE?4t_%3|L+ssUkoy) zcIGY?gdC|g^Yk0q~7NK!s1kxPK*6U3{iZuNk_n96wAE?WI<~(;^mo_)e%9;&ANwB zW87Jw_k&FOgw0g^@yq~Mkf=wJo2Lay@yyHvlQLwn38<`r8_Sc*8gIT&?ssddyPxbq z;)JIl1eRO@_0CT=Z8re38ccpz@na*7RL}CDiLGn}brS|vi5$0V4B&Q4M{KQ^Q(iye zfsH11=b$$QV7mQI9QpL3xI0#@xqo$L#|}2L2*4Z&*!PmW^+SbBTu#XXEny;E>VYMs;A>qMZmMc((zicqLmE z)hx-<U^fKTy=Aqa(dL$~7jb%{8 z9~&tHpxP1~50#RnXMa(qp0snwWeE5<&DCujkv_eZR{jlC3g2qCkcHu{7WPCND_779 zb`zrep8(7NGyf1)t(NFYfvuw!#aoy1>u*c6hj#8HW(5=SK<+9dIh&Zdw~ zC7o(7p5!RWQIex1N5_HqDDhF^qr^uCB>u(Tj2@EsTCaaqgOVt|ZPkz@;xVgIlDN=s z)vzQkhpifs#Ep$Zj!NRmH><`ZDZj1ay(t;eIZ38o@DDLP1Yz(Y12Q2CvLOd@ArJB) z0tJ#Bx)sD?tkC@j%xSD93T19&b98cLVQmU!Ze(v_Y6^37VRCeMa%E-;Gc+fD~HeyX~v*F|6T zdTUi57R`fx{+S02XaaSQ0%zr9Pn(McrA4{y2ZaRz0RhYZ%&Q>`TnB5dj{?000|%-7 zzdJ~fO*~cCu##BNZ(!fRzd?M1Y%Pg}>P!S@=j2YCj|HUx9O4c!+N>2)O(xjbSKR2{sq5>Or;#P2O(e6i zEnGIitoHAU*t|A!ax(0mUW$}3!LqP-?=&j088h#YRmY&C)NJisW)|7bSbS~W@Z=WS z+`?f}kCX`hE>U?NIBUOSBUOwBO7D{bn zOZ}}+sn$Z>NdFvnu0S9mVwV}4AMUS>Qy^WtJ%TQZj~m9eTZyP*IHtp- zOl21sE1>-`L~Mi zA!g%Tc(PbZ&3$C{1Bq+`e9#5pg)O>PvTAAp%~?E!MPv|Vvr2eKl45Z!{JX8Ec3^jz z(ta<;$RXB~YDPNM{38~g6fp0Nl&hl&P+Pj5E^H2u+Z>tUhsV>2aXFPzY~7~)Z(UGQ z2XrsusU|B3Cud<$L6X5l3^>T2pIe5Xl>W z?A9TcvOKQ8{IEOk3f$Sd_QICb$8UE?jOKA}MlG);8hw)DbkUGmJm)?0z0Xg4@I;Uw z-f5nh#$D?#b~zY<^D%=3oz6jx`D3dLe_%)W54rGODJ#KleOCEZe-8Eumvu1_m>xHG z-+Rh~KBWqT$NjMGPyV_43YWpwK~5F~2!L=rQT`(=zGc@X9d3j8;W64()V^-waRM9D z{gS$lcgR$kLw;*lT(~*SiC0u>0mgp9Bd*F%w3E}><-AyWJ>w_Eu2~=O>+7g2&;JvK z0=Q5huM$w20@rg{%6z~*BPuQ7j1g9!{AJ*?&n8b^NLk{Ox&1)B#S=e2tElx4uurj% zDP^#U`toOEn50-7r+z6%5iLt^g2jt-vSx~v!jfPRWF$xIZ$7+gzOVqKq4rirV>!YI zWM`}U240TdBB{$cl@ed=R+oZZz=YlFqfKpj1nVHhMxKqcRkT>)ucr0!e=68qx8uUU zH+6vtn~WKq^{cm9bfH507a%)BQ`!MK~Og1OXF~!9G)u@C*>eAXrP6vb*g1r2BUIR-tE?I^6C-%>V&~g|sXSST-#1sze1`S60 zHMK^lHcBcL0)@rPks(Nj0f2Dp6~bZ{f6%^MINgAb(Z)wSTs00b^)2iaE`#aiycJ<9`NHG_~5K-1@_9=(cyGPq|{9w%|6S#k^mc&vI^oT^KXyU5mFgNNgga-@+53;(QDJgbANY!eVEL!?LM?T8yK4^IMsqh z)k@O|f2KqNwOiB+ye3g`+>3 ztJe3c7lYSI|E)$z zQVOGyr~XVIQx<2pFgQVcSk7MbwKPxgLgjqkukF|r4PrDI)j$Amf|%<_G?CQezk^o5 zmQO;cKzBI;Tuyb~AL0p^rGzf$QL;;;B5M?z)9u%z(M>1HPie|BzeJVLLE=!TBo60y zQv~Y@#zMc^wtN9g04&K8H^bTb$M^^~fyQRU=_{<0_4qw`i3hq=pQ?N>e>ZeKaHpos zXfA7$>f+TY*3juY*JwVGtKO-*MzaSO9%Qkgbji8_x-u^1-R5Y9wlPt|`|vc(?JSyN zrS~NlqVJl5MAuT?qkK>L@dKEd+XJ%Ik6TIgaZree`AJDcl$axo+bq%NVU+zK-GfGx3xwGdv zS$#cI@T+$kPbK;DjNKV~-+ippZ$t9fvxl*453?;Vo;8ny1j2WD_)0n6d<7pXkg2Ej z2W{~!)~|xa(UG_s?&_Oq%ub+X{Uun(F^ef6=}MBppf_sN^u&`Ci#X(dfosgfe1jH- zqoXj5{as9HzOMce7cz=R2D@erS1v!xR>bTxv&;VTGg;I7&=d$>77ddWW#t?DM?Y!yXYJk`rFy zB8f7;<)|NR>|zj7eytTg`(MZpMuUvPi;2Lft2&)@7MsgJ)QPVo@C;@{~$G~Yr@}j$Uz$uAUFU9Tk zEm+)r8)^#z)807so`&MvA*l8ni}tuD3~Vz}xPiQ}1jV8Ej5r=G#nx zdP)U95EG!FaAG<7ca)p>n9Rk|nKs)%Ep=KTJR4f!&7cNmE~kneXqj454a1zlLjTaR zs%znEuMh(M^(X14Pxbo0zvo8%r6tUpH2MkKDOJC8I)o2?&E(xLxd!lt3zVI9HvI86 z5ht6gP%Ws*yUnZOA0RT{umR(ajL8^H*jLQ)ijoJA8E1glFKxzU`>Cu-bW0^W>AOoN zjVq!(jNb>Tq5|u-ICB_atKMwEOYznsMLr@OU|S10O#145z0|CtCreruq~Cv+T*I^B zHmP%VOu*;pmHHXu0UH#%CF!XyW{wkFH=9CvwY77S4R`oVLaeMW^?-d&apAH#F-wtL$PLzWBa~-E#+#$QyNJ$319)g##U#(psHzoQ#h}Rf=5loPL66YhbD=~(tVf9;NtIKl5@UY0ByID{J7FUo;1#2lY zdsY;eO{0G$TT)EzwU$-$RGa}Sd>xkuc|v0I$)?3G9e;Rxw;eBnF37-YQ`J zaA{WzQPz|}7vdtbhZ4q}><>SYJQw-eT&qTyrj5k}TH(Fbg--pT{&t`fqsEpP7MDSS zj*Yg|Ek6|`gP?bU=dbr_y70mxG$p?v;8_3nP z$L&&pV`bOlb}x&~N6g@}k7K#0J`k)B2g6I8@@#~Ri?{gcxeSZwtOkb@8x50L{xxqNs4P(|X1kjnV%U#OESqgxbvhR}I zvv(qemee6M+L$8(*lOc6T3Er)51NS5cY-eQK4*p=Dv~aWHo+bViWxgFGWj(!<3X(L zqUcr+!c=nz{z)Zw;w}H;Oa$c1kX5jNk2&bWcJ@7=U*%OWn-tOvxa#R^M!cR!8t5uJ zO`FacMZ3*cJ1&o*lG^yG*V}+Ou&9h)7_6YwQfG71heI|PaPF8z*ECf*l3Tw~H0+pa zqSPNSgDA;`_o{kj2?%lsB5qZ^FnR3aGX?XI&l9vM8>7RhS!2VC7~adtwi|jH(I&n@ zA(s7PH?US!Vn6Og(VIeJ55!BW%E+v}HssA@PN@#I2aI7OccL3K^wt6EP{Z~7hI58F zJgsHdSK`|@&~PZ9<-|X(bwdER-n@?t5%fV62+Ip*Lu_urDa6I(l5H$7!=K?V1Ej~N zlQMW6HaNei(b6K$OUU$N(uW|4q4p}Rifg3p6Io>2Dmh-^t}6C!63{>*YhlRMoFs76 zTw`G(r%?70_i%!tusHx6$6_qop!ZNpG=(hq!(ci5%SDO53q4Jf0FRyf*b!+cKgS3B7&+rCwqLCAHHmdL9wbdME4`_vyUg+b=Fd- z%FCg$Mc`aRfpf}l1d9ENwod-eF1BMCm<;{9Blbyj_jgS#5kZiw%6K7Z=`m>z{R;Pw zFHbAA#~U#Ui_GhqLK>FoBJBN2o?xO81IQt~9OunyqYvnE5JB7jn%^cXG$uGVH`jkt zg>GFP*9{&t|L58{at*tA9vGsGOlazi;x5k7EAi(0_=7++lL*YnrTC?_iQ2!89DVc; zteC45$RT}_;3xOzr9jhG?`@k~qgF;hj%AgxwN>}f8F0Ft= ziPEglH!aOZB|EL)jwy|vw8XP)m)5kF_sqfnS1Ix3_Kaf%+n~x`^Kw01MUKuHfLSXM z)=5D<;yHp2J1Nsy%cY?U0ODJrFDQH?Jxe?y4*H7F^bj^(xmKY>k}U z%G}F3y0RC6;vrfa2AZW>4t{jK`WnwXg3aEn)G#$TRcd|xaCVft;pbO?9;=(gU!Csk zt(^7*F9KL*Dnyr*DxUm^El{fc%h*EI%%5F5b%c61ka1#jYIVwB1BnMwm%DrWRxM*Y zl-BcM`OR(Y#bKcj5)>xt41)Uao=oUB5UK%@2+UNFgoD;9Z=>uv#q)q%xAtZKUw;?m zQ(I13j}uJ}!q~hMT*qMmBl^*^?pfIEd+{;{Qt{J_aFY&9h16}xl8kJn_ULqg? zLx1GArDZcN0&l240D+4u|J9j`h*c;StTo4?I5DB?L)aao_9$JN#c!8}tb?dTxed@* zt0-YfsSTXO=a<;MSad044<G?eJg~-~6g`!ak7?NKs}Q@|TK3+s#Q=l{ zZoELqmp!l)j7HBf!~MvLDgir{^Gx|?9ZDA|LMB4+0=EJH#67~}{?ix#ohH9zHV=25 zCk*G8;6SoJ`OZ`f+E#?u(OBtG`_cU53R;xYc=l1~3inT| zCFHC;RC6bNux5~X{xEsWP?*8-cQ zG=7XAS|g_2ugY`o>($j+D8G<3Wh&W1yS>3%4l_WTO3y~?!=z09FN+WuGlVNoUQyTA z;UN}X{Ot^;`-M-O5XVem`XOps@NsTTCL$LW+b$(mOp`8X+OSP1+E&_BKKAsEzcs-w zHwWAGwl3;YVcHw*W7E9-97pzs_3i4Y2x7}0`UgmTBmAQuyQD;0lE>aeWH7l0^EvzH zo1p+n0sAuyFf=rQmm7bzuvoN5b6eHL!E$KiXF`EeS{L7bl}HY3-!qyKi+1qOO^i%) z^kzX}9?-aaZ1PFMe=O;aG8{Ikr#*(a+UtGPG`}h0#nvcFV-w(0CEMw}9V9GRF*Ul_ zvHEntPQ-g)tYx&{3FZ<7jLMKjQ;{qrK79ctkhtNK)?A-{dU~WcaZ92j2sidOe$ALW zcCewacd9it+vjAx0K-j@ zn@WK~G71lPS((07QFxJkE41c@eH8alI;{bReq&>iK z@M3)WicM-6j#b+0sOmj2l#@I*v|MG);OQTgIgqKs|?ge{BWX>b+^ z2@r(nPT7^`!Ayr7;?7Xq?n@F5d9;v`nT*2>tm>My&Baun=3fV?HFGw9B!?8!5?B6-5?8vmf8MSI>cUM};W{dOGMJUFakz?O_s{-K}KwL?^fLFkye$O=84B0c{<22tKvz zeWWG|^lLSf%|Rh5A`nJ>4{HLR40#Una1lrKzv^5sD|npI(JG6|WgCDWPXk`&A{GQl zf+HYhDhg<7-Mq3}aNhW<2Gl!pgTrm=Y&p-9vzE zjR_F)0dFh%=r!cmSLr$%cK;b^=okFLs+c!v_)@iFtexoo7m}IHTB^TkGh)A}0??Vf z4&adZ)uhz^k`J*?ZEy&7LtH2=x38oIi5HOL~jjLmEP88i^Zhue;J*eA~oAr zdpuFENt6T2X?98n_%r(x4VQ{{JnMP-Qj0Aumzw7DxBE5LXEQC=27Uxw&mQz{Zua** zHiPO!t)i(3MZK`T%PLyvX1>^UpJ3+r8_U8^w+(HnYn8(4U>wgD11IK(TG!t7n{<~~ zn0i$PblgA3%RWRzkxJfZrEPdTnt7@X9E7~-t=NhhQ|ac8Pecw$GFX1W!QL*eZ6Qnn zcC{Cv+k7%^53a9I^y%yw02+Y)FM;+w@ryVSQ#<6AlI5-ucp75JnPj& zaPoB0o?CxFFT8tWc2)d=a74gF4?$monGSXD-7SPK|Gk(_!P?0###2cDm*=Prvzac) zBP8$y&$J{5g9NQ$Ubj=KT`9s0gj7 z&WBsP>Z?U#y0~8A*l$C25@|}*3f$HViD##IX$Ao7zrA)cC8eO5HKwo5;rn%+zoc-f z)W+kpN0wlACU1v7k$6e1xuHIk9QT%3awxO*e4npkZ3-}1<)wQ$?BUU$zj!LA?~ksi zUl0eNQw>nMv&IZ3DY5|;G7S;vZD6=yyCwXFdI)Mc;&)CP2C1beTN>L+%2ArlN_% z@7GJFroaVReSpD)N`$F^rt8`h3*9h5)a@5UR0&GRFprXMILNIug3RcbcjCVBb;0n> zsK$iL1k99+{UuC#g1OTtVf=LKje@ti!^x+ElwJ8(k{b`%x@hRj=M4n&0mY)C#%caIxq8@| zkj5wH<+)FTmy5pQEMiLHawye>HNBei*WHV4A=#_irSqmg=^UjQpJ2D zm*5jZ8X5cigjONxd^=E6XIB<<$2J0-jUyicMh^v2;;*!;fQ9<>;8(j=O>JNY%*>7P zn-srh*8zac^hP>6>tz>@V=^5!>3w^u4kjc0bWNUpP&@r`sVUl>kjF)iPAEX3%ys+o zuAI+li>lZ21X6PQk>$>WdPfKHgtt>^J-v9_U=;=CX0;fg36=HtC<)>BhY?xfmoW)+ zbJP*S^Xd2XoK~V`OsdzxU&b5eOY7(KmRgmzQvvWalOvi|pqByvQ2sd6{r-e#zF|fp zmo%E26T^i!^#@*S{511P#vX zb^v(V;tb(0dqN(amlbr|72d20_T4GCcL%8;p?}s!=2=<4(tRYHR{dF+FBP)#a{w{a zek<%zPb?-fC(7mB*^r4|jUJ9D1BV0!v?OA}qX=?;y zKLmF?y0hQ5aSQ{ z$X^N*$tC9p_kGtY8ET6NX3(2Xh!m_-Wr*Vb<*A*NMO-`AOpHH+e{oBb(_icmsfA z2vv>3Krdpt*j)Eu7UW-bM+5S}rBoU%W~5?fau&1mnOx$EJ}F)l7zT5^FVs2Xc94<* zExSI}(q3=>1MOP3_-$I!zOd}s;f~q0}jT3k&^OBN9oafW(2>i z&VAT(j(2#fuz^F<(*_3IZr zXQ{u;|MtwVf8x;JV7|eAgZl>m4dEN&H>7XKt^dSPo=`wp|F^HE;b@|u1G@35nfh)D zLkW@Aw%T*VsN)s`oY47Q2(YeJ|EhV2%9Hd>k>rRjU3$2?ce_8T5kl-wu?%Neth3B! zn*_1fruXFyX|bCKpvOjFC`=S5%fnH{MEyDD=}31J7G0oYL! z1Q|XT9BEMkP@5)-VAnV#99R*mAZN(pd07=`ORzU-7&f+tfQUI~3Ks9hg!!?(pFt$# zrIFSG|3Ekr)u4TM6mk+XZ>aft9y=AQV?sjURiR-ifw?#XlYF6|Aryk;7a9ygQGS{-2H1g|OKRXcN}sj_#S&>u-yBaA zk|<3Cy0M!KB8I|9(wY@9U^T3af234;>`HF7vLei+OSdB`h@`spn1S51NACj(C$bp_bR zr>U2hEm-L$Dj$q_N#3*GM%M4Vme8gQ73`Ej72oOS$;O_$3=zj?6NJDJ~(HZXwC@)6k+2?L?ow> z5xLwKjdx_4$w9@jysCwFHVEDcMW*BgV^H?X`QGy(FU7VMomKBDF7v~$S<*%AF_a}j z{bmRHO#r%nzHqLRYN4-?Mw2f&PIx4hx)h7Gze?0hz=cZnm+iPs=V{7psiK8ovIh6r{Nf0jW_@leQCi`=21jTEA_3u?iU>n>x34m|Oi`6OgnI&64@iNn z2^66c3+*JiRl7Z+C_wDuvPYP_ zEdmy_xwmk8nEl#Wv(?E#UYTN4mDRCdG4DRsK9m$-WAp(6mVb0SXeurH;;Q~YV3Vto z7vMm7%-5Oen)8{O-D#x7Gk0MN)smC1=n->f1=P&{H*~-}`0-l8D7!#D~( z1&M+lA{L}@ROp=%;U03ajul+}6ku~G6>vj$M+8FL5)uU)2_`kbb+%_%q*-8I_Y;`^ z<88OWND2A0(|L;{=`>Cz@fbZnnqHb4(fLPYqd@LQhQ5KL;9BTRmkYKaiM?lZ3L2CT zWC4T}I7Y9f9b@bC;N{X^7xk^;xE5lEMkV6j_b>$MDx^{2Y32)I6qCivt~zwJW2Kq1S<5IgDjzw|5x2BR>n53pRQmU=npOw0#~ zerC)36Q8F9a{5 z-ydEsQ$u@~9$bjkZ={Dj@~w*fW%}}k*eZ7gQM$}7y$Z`%}&@9 zuFzY^?*vz$qHyUg8qTdRI|nCDwsXGCbE$U2~S@$7y&cj*$PX6&*aA zuIVS;T67w%`sD&PzSf?WbOn5~lWR*2{N=M`PxvfvInoLbC5!I2vtzeu=`OvER0RwK zoa_mr9txeVw{nns{RF&Lp9e-a}t%1~%mrq814Kh9m_f(ZB}8aH#8PD_l%!zE(*SfnhE7HtyW+3=U}MmxN`kD?kxM;5 zmewI^O)loqk6}fhP{5h~g{r{|SMOD&AH3NU`;avVmCU>8iC#ivO&Rv4Bi~ONxZD&n z_o%z`0ixauGtkQ${l3rQ@;A0bPhoL&?VzOkp&HKSLzYp6m zUPv5+q7~~;KlrKbsds$4$Oa23Q^DKY(oNlgTk3@x@ zy`N;Grtmv0t2U&@uD{z0Ym-^HO!ni0U1k26$iMd{M{G=8DCpkLxEF5EhzTBmwZzfb zi`DERn56qkdRHLuL6)y2eB;eJH+a^IeVCt(CHJR?0TgNAgo;BV_s7WliU8uRCLPbV zm1l#WZ=2uy)rI4G3Eraw4*}UrARQe53tF_l^D=!#BonOy8Klv3z6w#@4!Q%q{~2a|6GQ9jhK9Hnkiy3Iu_d1kX?ITH?VY2-|Nuo_y(MQ;K zowaf)y~vx)b)`4aEt>|n*szO?ejQs;b3lW;kXU&zZnDY2CR~0#f zT4N%HuOc}`4<}=7DHcu^>;e$oX{kc{JgiL)vv4&r`t1rl`Hg`R54^@mjkugRFu=T9frB zxU4+v@w{|;Vu;8^hAi5dsFg#u7K`Kc!KTwLFM7KYE}QEHGCqM{)1 z_257zD4KZ8i2{=dbr5(F)v%DzgdMOBu)lpO;84NmFLW}0Nat=ySo^i$AmpMggvs&J z<ybopg?amt(m{#P z8=zr>1E%IkKA`Yq3KfpijQHMPoikbe`EZX3w+%Lrc5h!ULbjFjE^dn4a_7nS)1;2A zq>gJl;)*uL>(G1n#*eyz#+gc%{o!UF%LTKYarXlFIWgB%(U$EKSJ3yqPh*AzPt>iz zNu=KN$3~#B?<~MD`?1kkaE#ap4Q%n01%EeVI0x`iIjAHUFzdXzOT1=vqp*Wb<>e}t zqi~_+<@U1MB_+0|4U*>e%7wovG#8#7@aI_A2l^rG^dq}GD7L6*oIBbdq*r)r#|ptN zq-7K1<9bo$DFoEYn~&fDY4mw8`uBCvfl#>!aW>L#TyH=D&rNaQ^{$sddE?ub@jIk_ zE{dEWdRN=7(?49HkznIn#M|1m#@QPI)QnDNLD?PG5YL3b(~O^K0sIiHxYul-mM%X` z-b1=i5bqfHWO|WayTh8E5A$v`UMa!(W`k4;SzSesipcL)Z?;qrprH?AvNgV_aX|C% zB!3d9{G0)3N<>hoEhkVfds3mfuuFOdzYQ})`J8eP%a5q2jvhhShXsoMj{bVQMeq6v z@i#~tqklFUdvJGJQ(7jTDisx>1{60sBA$mcGg{EsCGw|n+9icow_tWuWCN_#<_DIcJK_~Ol(z&Y}Q@ehvBC%(i0598!}<7W##@+5L-tiE95WAhnJ zh-4bm#g8Lz43haB{w{YiOo+IE3O{i{G!>qP5qtFjfo{T@gW6mr{9UIah<|XR0$0Lx z^9hIm6MZSs2)_hU%*N`Jc|i?HE7EMui8c+;b3QFCAX~A6@?xpwXZREEm-t)uht8!O z+@}))UoQKw2d)>vAHKsi3y|^ zVa(Lh|E#&H6f(1K(^LPg!EZPteuG2~P`ifnfNVW=;&tt^UGy5>WQ!K}TPm#2N=G9s zXnhxj-VLiaWjkBl-!4|-gj!=0B$$O>j`C5s=fd}&z7t{gQ=L^F|D^B)d@uqQr%8>& zkK(sV26>xjbO?q~MusgOZP}`IvD1zb5_x^dcUc-b-B?4kzj%5ognGJbe(0|PI(0Gm zxN5NZG~0a2$^A6W;yXC(2*p15mbCx)NEEG=xAdXtPKAU*;(2!8|)mszn^pjEc_a&3s7I^BkubPf+=Ug(ZG$nmy%GZpkHU_ z9GDutJ{(GAj*IPUjy?~ij)xfmC{`HJLPqVhh1Dfb5=xq$*B$Cu#S%Jbwo8s$xxl~& z(u(-QD^teTH%aH`HN>nSDqw=z04Ft(o7PUY?jk09z8L0JUdntw$Ez|lFDT zboit}_^9O%XpkSEULAHdtPR929gKu@8bNU>3Dw4W19#WO9@y>wldfFe=rYMl($P2|v;@;q~skvWG`7r6rCd3#6kf zBdwQO~+6nuChIoA^R5Nqpa~prW2yhXUmnu3H$SgC9 z$KF!v>DerlzGPKHvOr)2w2{#U>_itF^j(Dk3I05e+Hm-JWh49CTi-Ft{+i#Lj5CA& zEP4sn5UEC3+UOtKndG=tI?Fr8mzp{}KJ`jKDhl)E@C9;!0zI}0L2V?=RTA{*cRhz4 z*A}uCcVJi@ge?eo7ai3@J}GR#Z3s3)?y7}NfD29F8dxbEc960Fod`$Hil-);bz5}@ zhhjf|<;$YAH9xy~$>DYeHc$SBS-9JUxli0pq%BS!Vx~4K1%91nS0EKmN}rJ*7ukg{ z=_~aEmJ@@t8x7S-tF;%O}oGnKcPX`RmSRG22(ZeOd!Zmol1 z_moH$xWcy&=e2IZ$ULwpwpxjnQ!@vJD9+q`3dh6Kd*N!4RRb5>Q{qvV*f4AMUZHSJ zp9|rmV1itL$YZ;DEth%TmgdkIcm+4grv&D2ku{<|6ZkFjPnu-Zt0D!QxZ1BU@Az zhYnMYtlKJnDXy0K!ZX&;wfdL2mv&>)qGdf%pSG5A6H2ZNSTEnKtH$?w`EYjkx_tZZ zUq(3hRUa`$B4Ip;c$IK_X-Ad*5u4$1T6?sd?W8^+i0}Md9a>a_UPH4+vqD`o42)<2 zmlA{7%ZUxI#5Z#qDoHR8dGX)i z-RQrvx=-9}&VnTk_UT)P5_XdD3TTId6#f@P2HAwcLsz2?9AxPHH*We8`B)8_wbl^* zLR7*4Os5jIoMV*j_eq**;umWc!>NDPtW(g7yJobSc!F6Leu9KdN1%l=z3XA0f8%QW zznb69K8)Qu)AN)@u90dQ-i9qDVHV6Wc%Z;MMVHJO-E)w2H~5H?x5(@$ZS~_zXF|pN zhvNR(eL|)D<^O*2UnqM;r#+^>b^v$DW#l)2AlYWh-WSGek+(3KOE~->K1H=Mie(puaX2KZXoHx$|y$^Q??+`+aY! zS!i?b-%V021iofNm|wI$V}$O$?!GJ=Wz;Y%iBCJ8WwbV-*|jlC6HnP zRk0qeCrv5g64MuMU6vj#^=?-!dAB;UL9x`;6wZl=_sP8?&e>cxuD>#Xz{k_!Q78xK zEfc<(FX+I#_Hy|ZjH^B;l17n(9|@!PZu2B6*M|&oib#KrvQtytIB~o4UOb_jMw&te z_SZy2secR_hT62vl(fo?S@_G?$_K<7ujn6>hKE(=o;zmjmq^vj=Qf z5K3)m@g0iC|8@&JUIV%X-|o$0dik?4Z?!ridxvGp6MI)UUn(VXdI#(`h!;4Y&b-in z83p${fkgohdw-a)x7Q#e9=lIES=n66rG^D)wunh@oJ~d2XeUkoL>_G+^+>z{gx|@A z1MYjVlLXHHB7qw}D);>^E9BjGTdxmwdY|RG{__l{o+ee>+3UYf;-9$8g(-Ycy8k`+ zTT$=h(CBNEJWyOAZ#676Y*r905!a*CF7Jc>JdNny;w)*U>b!KZ!!v*n-av`d)ElD? z?&LqfG}EQ{Vh9dG%s@7^j4(KZ9hd!Ow_b-jXpz<~#X%1k??7cE7DM3N5VKii5jt|9 zElKe`vy@P!RSoW{IdFS15Bu1>5JU+yQu#Y>#LTqfaX)z0r}6O;;?ca1Zt3n5l-Lqe z+G6+Ad;svfGxu|e6wIPlddlPw-o@3r3An1Jz2!smAdDE|jIABRBcdFb zL?>_)@OJ`Caao0*TmN%!k|$);ON)L19rvY9(i@Z%Qc*>vc4vy{XxjU@1*UjZ@85}A z9cp_2JC0wiNwy&5?ndW6v~9~c?z{YT>qaQlxwxG*f8xJcdCl?zpJqs%} z3kx+Z9J7Lxxwwg&1qqEf9~%oRH_QJkN#m8}&_>}Rp;5QcCXu)Cups$9=S(Y*<-i7I zZEct3XaiT3VCUkJkl%Ud{F zxmlBNvvHFA_vM6gNlrP&8E5qN>wbmT=~y7w$^WA8ccUH!} zYEdO6)ZK?KB>?*Hp+p6jOII#0Mv??B_yNd9%w`-exB=irZ!o41FDq{#g;t#ROr^uV zv#tqd=eQ@8Y$Q-zWC-BoI5&aq@*&sc~cK#>O~EFZq3 zK+3oMcxcwT(ZHtw`DM7f=qJJN;abOL+R=cb8%e-$tx#4ssPW|D%y@YQLXG4 z+9Jw{qFJ5%5(cTK0$Hpr)42kcv4khn*v_%CZ}zBn4<9`ks(rNnSOjHLc=ey>wTueU&gT8 zX8C^tP5`n0fFx)CKjfe!A6w9MNP_fhR&n2y4C$lall(FLERR,"Incorrect args for pair coefficients in oxdna/hbond"); + if (narg != 27) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/hbond"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -638,36 +646,42 @@ void PairOxdnaHbond::coeff(int narg, char **arg) double a_hb8_one, theta_hb8_0_one, dtheta_hb8_ast_one; double b_hb8_one, dtheta_hb8_c_one; - epsilon_hb_one = force->numeric(FLERR,arg[2]); - a_hb_one = force->numeric(FLERR,arg[3]); - cut_hb_0_one = force->numeric(FLERR,arg[4]); - cut_hb_c_one = force->numeric(FLERR,arg[5]); - cut_hb_lo_one = force->numeric(FLERR,arg[6]); - cut_hb_hi_one = force->numeric(FLERR,arg[7]); + if (strcmp(arg[2], "seqav") != 0 && strcmp(arg[2], "seqdep") != 0) { + error->all(FLERR,"Incorrect setting, select seqav or seqdep in oxdna/hbond"); + } + if (strcmp(arg[2],"seqav") == 0) seqdepflag = 0; + if (strcmp(arg[2],"seqdep") == 0) seqdepflag = 1; - a_hb1_one = force->numeric(FLERR,arg[8]); - theta_hb1_0_one = force->numeric(FLERR,arg[9]); - dtheta_hb1_ast_one = force->numeric(FLERR,arg[10]); + epsilon_hb_one = force->numeric(FLERR,arg[3]); + a_hb_one = force->numeric(FLERR,arg[4]); + cut_hb_0_one = force->numeric(FLERR,arg[5]); + cut_hb_c_one = force->numeric(FLERR,arg[6]); + cut_hb_lo_one = force->numeric(FLERR,arg[7]); + cut_hb_hi_one = force->numeric(FLERR,arg[8]); - a_hb2_one = force->numeric(FLERR,arg[11]); - theta_hb2_0_one = force->numeric(FLERR,arg[12]); - dtheta_hb2_ast_one = force->numeric(FLERR,arg[13]); + a_hb1_one = force->numeric(FLERR,arg[9]); + theta_hb1_0_one = force->numeric(FLERR,arg[10]); + dtheta_hb1_ast_one = force->numeric(FLERR,arg[11]); - a_hb3_one = force->numeric(FLERR,arg[14]); - theta_hb3_0_one = force->numeric(FLERR,arg[15]); - dtheta_hb3_ast_one = force->numeric(FLERR,arg[16]); + a_hb2_one = force->numeric(FLERR,arg[12]); + theta_hb2_0_one = force->numeric(FLERR,arg[13]); + dtheta_hb2_ast_one = force->numeric(FLERR,arg[14]); - a_hb4_one = force->numeric(FLERR,arg[17]); - theta_hb4_0_one = force->numeric(FLERR,arg[18]); - dtheta_hb4_ast_one = force->numeric(FLERR,arg[19]); + a_hb3_one = force->numeric(FLERR,arg[15]); + theta_hb3_0_one = force->numeric(FLERR,arg[16]); + dtheta_hb3_ast_one = force->numeric(FLERR,arg[17]); - a_hb7_one = force->numeric(FLERR,arg[20]); - theta_hb7_0_one = force->numeric(FLERR,arg[21]); - dtheta_hb7_ast_one = force->numeric(FLERR,arg[22]); + a_hb4_one = force->numeric(FLERR,arg[18]); + theta_hb4_0_one = force->numeric(FLERR,arg[19]); + dtheta_hb4_ast_one = force->numeric(FLERR,arg[20]); - a_hb8_one = force->numeric(FLERR,arg[23]); - theta_hb8_0_one = force->numeric(FLERR,arg[24]); - dtheta_hb8_ast_one = force->numeric(FLERR,arg[25]); + a_hb7_one = force->numeric(FLERR,arg[21]); + theta_hb7_0_one = force->numeric(FLERR,arg[22]); + dtheta_hb7_ast_one = force->numeric(FLERR,arg[23]); + + a_hb8_one = force->numeric(FLERR,arg[24]); + theta_hb8_0_one = force->numeric(FLERR,arg[25]); + dtheta_hb8_ast_one = force->numeric(FLERR,arg[26]); b_hb_lo_one = 2*a_hb_one*exp(-a_hb_one*(cut_hb_lo_one-cut_hb_0_one))* 2*a_hb_one*exp(-a_hb_one*(cut_hb_lo_one-cut_hb_0_one))* @@ -718,6 +732,7 @@ void PairOxdnaHbond::coeff(int narg, char **arg) for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon_hb[i][j] = epsilon_hb_one; + if (seqdepflag) epsilon_hb[i][j] *= alpha[i-1][j-1]; a_hb[i][j] = a_hb_one; cut_hb_0[i][j] = cut_hb_0_one; cut_hb_c[i][j] = cut_hb_c_one; @@ -728,6 +743,7 @@ void PairOxdnaHbond::coeff(int narg, char **arg) b_hb_lo[i][j] = b_hb_lo_one; b_hb_hi[i][j] = b_hb_hi_one; shift_hb[i][j] = shift_hb_one; + if (seqdepflag) shift_hb[i][j] *= alpha[i-1][j-1]; a_hb1[i][j] = a_hb1_one; theta_hb1_0[i][j] = theta_hb1_0_one; @@ -814,7 +830,12 @@ double PairOxdnaHbond::init_one(int i, int j) error->all(FLERR,"Offset not supported in oxDNA"); } - epsilon_hb[j][i] = epsilon_hb[i][j]; + if (seqdepflag) { + epsilon_hb[j][i] = epsilon_hb[i][j] / alpha[i-1][j-1] * alpha[j-1][i-1]; + } + else { + epsilon_hb[j][i] = epsilon_hb[i][j]; + } a_hb[j][i] = a_hb[i][j]; cut_hb_0[j][i] = cut_hb_0[i][j]; cut_hb_c[j][i] = cut_hb_c[i][j]; @@ -824,7 +845,12 @@ double PairOxdnaHbond::init_one(int i, int j) b_hb_hi[j][i] = b_hb_hi[i][j]; cut_hb_lc[j][i] = cut_hb_lc[i][j]; cut_hb_hc[j][i] = cut_hb_hc[i][j]; - shift_hb[j][i] = shift_hb[i][j]; + if (seqdepflag) { + shift_hb[j][i] = shift_hb[i][j] / alpha[i-1][j-1] * alpha[j-1][i-1]; + } + else { + shift_hb[j][i] = shift_hb[i][j]; + } a_hb1[j][i] = a_hb1[i][j]; theta_hb1_0[j][i] = theta_hb1_0[i][j]; diff --git a/src/USER-CGDNA/pair_oxdna_hbond.h b/src/USER-CGDNA/pair_oxdna_hbond.h index 1c9f37bf50..028853a087 100644 --- a/src/USER-CGDNA/pair_oxdna_hbond.h +++ b/src/USER-CGDNA/pair_oxdna_hbond.h @@ -67,6 +67,8 @@ class PairOxdnaHbond : public Pair { double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast; double **b_hb8, **dtheta_hb8_c; + int seqdepflag; + virtual void allocate(); }; diff --git a/src/USER-CGDNA/pair_oxdna_stk.cpp b/src/USER-CGDNA/pair_oxdna_stk.cpp index 9a5779afad..9c3a11b9b7 100644 --- a/src/USER-CGDNA/pair_oxdna_stk.cpp +++ b/src/USER-CGDNA/pair_oxdna_stk.cpp @@ -38,6 +38,14 @@ using namespace LAMMPS_NS; using namespace MathConst; using namespace MFOxdna; +// sequence-specific stacking strength +// A:0 C:1 G:2 T:3, 5'- (i,j) -3' +static const double alpha[4][4] = +{{1.11960,1.00852,0.96950,0.99632}, + {1.01889,0.97804,1.02681,0.96950}, + {0.98169,1.05913,0.97804,1.00852}, + {0.94694,0.98169,1.01889,0.96383}}; + /* ---------------------------------------------------------------------- */ PairOxdnaStk::PairOxdnaStk(LAMMPS *lmp) : Pair(lmp) @@ -194,6 +202,7 @@ void PairOxdnaStk::compute(int eflag, int vflag) } + // a now in 5' direction, b in 3' direction atype = type[a]; btype = type[b]; @@ -665,7 +674,7 @@ void PairOxdnaStk::coeff(int narg, char **arg) { int count; - if (narg != 21) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk"); + if (narg != 22) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -691,28 +700,34 @@ void PairOxdnaStk::coeff(int narg, char **arg) double a_st1_one, cosphi_st1_ast_one, b_st1_one, cosphi_st1_c_one; double a_st2_one, cosphi_st2_ast_one, b_st2_one, cosphi_st2_c_one; - T = force->numeric(FLERR,arg[2]); + if (strcmp(arg[2], "seqav") != 0 && strcmp(arg[2], "seqdep") != 0) { + error->all(FLERR,"Incorrect setting, select seqav or seqdep in oxdna/stk"); + } + if (strcmp(arg[2],"seqav") == 0) seqdepflag = 0; + if (strcmp(arg[2],"seqdep") == 0) seqdepflag = 1; + + T = force->numeric(FLERR,arg[3]); epsilon_st_one = stacking_strength(T); - a_st_one = force->numeric(FLERR,arg[3]); - cut_st_0_one = force->numeric(FLERR,arg[4]); - cut_st_c_one = force->numeric(FLERR,arg[5]); - cut_st_lo_one = force->numeric(FLERR,arg[6]); - cut_st_hi_one = force->numeric(FLERR,arg[7]); + a_st_one = force->numeric(FLERR,arg[4]); + cut_st_0_one = force->numeric(FLERR,arg[5]); + cut_st_c_one = force->numeric(FLERR,arg[6]); + cut_st_lo_one = force->numeric(FLERR,arg[7]); + cut_st_hi_one = force->numeric(FLERR,arg[8]); - a_st4_one = force->numeric(FLERR,arg[8]); - theta_st4_0_one = force->numeric(FLERR,arg[9]); - dtheta_st4_ast_one = force->numeric(FLERR,arg[10]); - a_st5_one = force->numeric(FLERR,arg[11]); - theta_st5_0_one = force->numeric(FLERR,arg[12]); - dtheta_st5_ast_one = force->numeric(FLERR,arg[13]); - a_st6_one = force->numeric(FLERR,arg[14]); - theta_st6_0_one = force->numeric(FLERR,arg[15]); - dtheta_st6_ast_one = force->numeric(FLERR,arg[16]); - a_st1_one = force->numeric(FLERR,arg[17]); - cosphi_st1_ast_one = force->numeric(FLERR,arg[18]); - a_st2_one = force->numeric(FLERR,arg[19]); - cosphi_st2_ast_one = force->numeric(FLERR,arg[20]); + a_st4_one = force->numeric(FLERR,arg[9]); + theta_st4_0_one = force->numeric(FLERR,arg[10]); + dtheta_st4_ast_one = force->numeric(FLERR,arg[11]); + a_st5_one = force->numeric(FLERR,arg[12]); + theta_st5_0_one = force->numeric(FLERR,arg[13]); + dtheta_st5_ast_one = force->numeric(FLERR,arg[14]); + a_st6_one = force->numeric(FLERR,arg[15]); + theta_st6_0_one = force->numeric(FLERR,arg[16]); + dtheta_st6_ast_one = force->numeric(FLERR,arg[17]); + a_st1_one = force->numeric(FLERR,arg[18]); + cosphi_st1_ast_one = force->numeric(FLERR,arg[19]); + a_st2_one = force->numeric(FLERR,arg[20]); + cosphi_st2_ast_one = force->numeric(FLERR,arg[21]); b_st_lo_one = 2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))* 2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))* @@ -756,11 +771,11 @@ void PairOxdnaStk::coeff(int narg, char **arg) b_st2_one = a_st2_one*a_st2_one*cosphi_st2_ast_one*cosphi_st2_ast_one/(1-a_st2_one*cosphi_st2_ast_one*cosphi_st2_ast_one); cosphi_st2_c_one=1/(a_st2_one*cosphi_st2_ast_one); - for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon_st[i][j] = epsilon_st_one; + if (seqdepflag) epsilon_st[i][j] *= alpha[i-1][j-1]; a_st[i][j] = a_st_one; cut_st_0[i][j] = cut_st_0_one; cut_st_c[i][j] = cut_st_c_one; @@ -771,6 +786,7 @@ void PairOxdnaStk::coeff(int narg, char **arg) b_st_lo[i][j] = b_st_lo_one; b_st_hi[i][j] = b_st_hi_one; shift_st[i][j] = shift_st_one; + if (seqdepflag) shift_st[i][j] *= alpha[i-1][j-1]; a_st4[i][j] = a_st4_one; theta_st4_0[i][j] = theta_st4_0_one; @@ -849,7 +865,12 @@ double PairOxdnaStk::init_one(int i, int j) error->all(FLERR,"Offset not supported in oxDNA"); } - epsilon_st[j][i] = epsilon_st[i][j]; + if (seqdepflag) { + epsilon_st[j][i] = epsilon_st[i][j] / alpha[i-1][j-1] * alpha[j-1][i-1]; + } + else { + epsilon_st[j][i] = epsilon_st[i][j]; + } a_st[j][i] = a_st[i][j]; b_st_lo[j][i] = b_st_lo[i][j]; b_st_hi[j][i] = b_st_hi[i][j]; @@ -859,7 +880,12 @@ double PairOxdnaStk::init_one(int i, int j) cut_st_hi[j][i] = cut_st_hi[i][j]; cut_st_lc[j][i] = cut_st_lc[i][j]; cut_st_hc[j][i] = cut_st_hc[i][j]; - shift_st[j][i] = shift_st[i][j]; + if (seqdepflag) { + shift_st[j][i] = shift_st[i][j] / alpha[i-1][j-1] * alpha[j-1][i-1]; + } + else { + shift_st[j][i] = shift_st[i][j]; + } a_st4[j][i] = a_st4[i][j]; theta_st4_0[j][i] = theta_st4_0[i][j]; diff --git a/src/USER-CGDNA/pair_oxdna_stk.h b/src/USER-CGDNA/pair_oxdna_stk.h index 950c276228..74a9a6bf68 100644 --- a/src/USER-CGDNA/pair_oxdna_stk.h +++ b/src/USER-CGDNA/pair_oxdna_stk.h @@ -58,6 +58,8 @@ class PairOxdnaStk : public Pair { double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c; double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c; + int seqdepflag; + virtual void allocate(); };