[MMTK] additions to MMTK
Peter C. McCluskey
pcm@rahul.net
Wed, 27 Dec 2000 11:20:12 -0800 (PST)
I'm uncertain whether I'll ever get back to work on the molecular
specification code that I added to MMTK 2 years ago, and I suspect
that if I do I will rewrite it drastically before trying to integrate
it into the current MMTK, but I've been using some of the more basic
parts of the code in other small MMTK-based projects that I've been
writing recently, and it seems past time to get some of them folded
into the main distribution.
So I've selected a subset of the changes that I've made to the main
MMTK modules which seem likely to be of some value even if I end up
abandoning most of the molecular specification code, and I am submitting
them as patches (at the end of this message).
I'll submit some new modules and tools as soon as I find time to clean them
up.
I've also included a few new atom definitions (Atoms.tar.gz, uuencoded) and
some test scripts to check much of the code that I've added (test.tar.gz,
uuencoded).
Here's a sketch of what this set of changes adds:
Atoms:
si
al
f
k
Geometry.py:
bug fix in Plane.projectionOf
fix typo in Sphere.__init__
changed rotateDirection to allow axis to be a Vector
new classes:
Prism
Cylinder
Pyramid
Polygon
LineSegment
new methods:
in several classes:
enclosesPoint
coordList
__repr__
Circle.planeOf
Circle.distanceFrom
Line.perpendicularVector
Line.Is_Infinite
Line.asVector
_intersectLineLine
_intersectLinePolygon
_intersectLinePlane
_intersectSphereLine
_intersectCircleCircle
vectorFromPointToPlane
distanceFromPointToPlane
_intersectSphereCircle
arbitraryNormal
_intersectConeCone
findRotationAngle
angleBetweenDirections
PDB.py, ConfigIO.py
add class PDBConnectOutputFile
createNucleotideChains - add optional terminus_3 argument
Bonds.py
fix BondList.__init__ to call UserList.__init__ is list is None so that
__getslice__ works
add bond.order support
add dangling bond support
add asVector
add calcDihedralAngle
add various DihedralAngleList methods
add class DihedralConstraint
add class Dihedral, mainly just different way to construct a DihedralAngle
add findCommonAtom
ChemicalObjects
add constructBondsToAtom
add bondList, danglingBondList
AmberData
add bondAngleAvgParameters
add valence
MMTKForceField:
add defaultBondLength
add defaultBondAngle
add valence
Biopolymers
add dihedrals, danglingbonds
Universe:
add bondList, danglingBondList
Collection:
add rotateAroundPoint, danglingBondList, bondList
Database, GroupEnvironment:
add bond.order, BlueprintDanglingBond, BlueprintDihedral
begin 664 Atoms.tar.gz
M'XL("&H#23H"`T%T;VUS+G1A<@#MUD%O@C`8QG&N\BF:[.!IV!9*Y;"#R;*+
M\;0E.RY5JR,6:BAL\=NO>)%EANW@2$B>WX78XO'_OJQ6+\O9HZK56CD]6]2V
M<#.7!S?%*$V3A`2$D(32;\\SQF-"4B&25#(1L_9`<!X0&@R@<;6J"`F.FZ+W
MO=_N1ZI4A28/9.IRDV]L.0W=J5A;TQX]Y],P+)1S_@>7429E&&ZLL55[N39J
M<_#W']O7MTIM\Z9]BT:<BS"`\5A=Z5^9V_>?]O4O?/.2,TZ%/X_;_I,T1O]#
M]J],4^1E4W0'P,)T!D`:97,F.A-@7VE=7IL`D\D=>5=F=^__J;=DWVCG0DR%
M\?2_"V[??__^3_S^3YD_])L_.7\0,(G^A^Q_9QI;Y:7N]O]TR9_-HRR;=^H_
M:6/LY\_\F61(?>3]'_ZA_][]G\C._A?G_F.*_H?L_VAK7WK^_0-@>1D`<18Q
@RO^P_L75]8_0````````````````!O(%2=%@%``H````
`
end
begin 664 test.tar.gz
M'XL("%MS23H"`W1E<W0N=&%R`.Q;6Y-;QW'FJ_`KCDL/"\A8:.X75I0JF792
MJ5BD*E&<!WJ+M<(>BHB6``)@:?+?IR]SZ8,%:3EE.9<2)"Z^,V>^[IZ>GIZ>
M<W9/X_'TY9.?]S,X%:,?G@P#(#7Y+A\U#-'`)VH7`UQ';\V3P3_Y&WP>CJ?;
MPS`\V:_??K+?G[O_?_1SPOG'']_^]C?/=MOMN#ZM]A_^NCJT^M3\:^^#H?G7
M2D,,0"QHJTQ\,JA?YO]G_WS^J^'+A^/AR^\WVR_'[;MA_^'T9K>=O3[LW@[?
M?//=/P^;M_O=X31\T9M6$">O-S_\TXMZ[\7#:?]P^H?-_2@Z03S5^QQ:0'DX
MW)XVN^T26[`[@></Z_MQ=]K<C<_>W&ZVLT(Z?C@NA\,XFVVV(/S5:^B^O7T[
M#E\-5QBN>K6_^_YJ]OGP^]WMW7!Z,Z*HX3C^Y\.X78_'&=Q]M0:E&@CG^N=3
MD8O9&A4?H6=CK=:'\?8TGMEVG"]FF]?#=G<:F/)TAA&\/VRVI^'J=X?#[K"$
MN\.VT4J_X?7N87LW;+97RV&JG`3`6%?C^\UI?JT7L]GKMR<PI3MUS@-^M:8A
M+T'6X>TM=KDJYN*BO5H@;_6GP^8TSEDIMZSO=\=QCF+O@++;C]MS>5<'('^_
MV]ZA!U[>S$#^<+_9CF#I\/IN!8ZXPTL8/`\7/(#7+Y^&F^$KL.+9B^>_>_;=
MU=/99Z?=CR/Y\3"NCOM[&-#5'X^_!A74/SR]6<P^VVSOQO<X*^"S.1->:KR!
M:N$:M99F\_0&A*+*T^$#H,_(R-7M'@9Q-V=!RRIGL>"NX_OUN#\-?[B]?QAI
M1I`()J-HM!;MI([@N--F^S#"Y>%V<QQGX_L].!)],/?+0<?%<B#@*D@(0*.I
MW[$`[FN6@RW?0)D1T-A@*Q=`+M\DPXEO58&K@-4H%@_B<I6O2`$`5X$I0"?@
MY8IL!9Y`)'-F!:G:QC==`286HTVN5A.R%9F&=&9I:!=+P]NA(=T0#QI'XAO2
M#;$8A.Q/='6H*+4F6U%1%I9#J,"DBK1AHT)59FPUP`16!DA7I'65EQHAL7Z4
M8JHIQA=ML?9'S[%^6YFZJ4*3?`%%E@4)CAV!R%?D29I%.TQ%/E=4&!A:EL4@
M=*V1O6B;&$`V-\0"G42>_8/0]4:..=7T9=8'8=_4N8:\8P1B/%`"V8"H4!#1
M\%RN%B)RJK7Q7<\(Q;2./E6W^N9@#TZT#?G44*R(UT?0U0V(G*O(^X98!Z!@
M*G*A(L\>IMO4,>"RB@VU-M^0I;L.3'6I(IM(##6ZVLCZ7*HC01141<XVQ&/R
M;`**"<UL6Y%K`["-$E2SP59DB^(9P\2W7:4@XH':4'V-B,D!&;JB6%R,C:0Y
M@);@*LJZHD)!T;ZAT!"[./@F!E#4%16!0:+44.QMO*2BJ?JBKOI2TQ<;2J$A
M8&1;4>#H(T@>29I-9!1-:Z.[,5:46[^LV*\H)E<7)YS<AK*J*.6*.&ZRK6Y`
M%'E0U-AOLS7H8E=13!4EW^ZR-3#DR.F48&M,#06R(:*MJJ+84:@H\P)/J@T%
M4#851=\0#RI6&V)J5GM&.%-M`+%1LFDV^(I"KOT0:=Y"8FBC0LAJ0JK>)D3L
MC!S>KAF*5K(W@Z)<W(Q;"N^U!`L-Y<>&4D-,CUT2;:$\.,29')>31)K+`H*J
MP[+QJ*X5!+%:W')P2V;%0&H0M]$*R?+0,-O+L.QICBV>%0RZ>GOIHAHTKHY/
MXUY9G*ZQIB!,^R!NF(KW7<)<!11L&@Y%9&B.8JQ*J4`71HE.91S4'COF\J#@
M)/K4C9FFHBBFN0@=:X%5$82C*24$857J!+[(_4:UPO81$S9>M*>.:\&!55F;
M%BV&ECI68IC*=>N4D8)\-U4E@;4@TWR4X,G=%82KCY3JLTFX",+20'.29\P5
M0VEG!53?N1*27`"J?L.S-JK>>(]CS'LXXU*&8A%#W%F[2/U&M<(+H51R2UP4
MD!&V"')"D.G8$2X5;>H6,>;AXW;<!-%%-540&#N!V2+GA`(GAH:%DO:E@!56
M$"X*L*2IPR1LHV@O2XTNN`S3CIR=!/8"*X&+(+34E95/%US9E1M9X$(F!X>.
MK<1E^NFB6H$1[&S'I0;F=BUP&;XNUI%%7EBA!'9=@97D(*S3!9,@*ZQ0PCK$
MJ<Q:[(((%Q\Y*M?K]-NN(=`Q374<)2Z6YC,<ZO33:8FG-M"1+73,119C+NTU
M5J,Z*(%+AJ2+:H470NE(:@5F!9$4E'8G!9ER@9UP"%P?Z20L2H(<\QGV55`6
M0R#L!&;K8N@*&/,(L+32J0BBBVJ%Z;[`VJ'.6O("XW"BP*G$$5_X?J-$<Q*3
M0)BKM=)N!2YK+=$YL)#I'.\[YK*',5=N.M/,QHY+(5HNM.C$&$M&S?4:XQ`$
MCAWG$D=8/=4Q$TZIXV@$+D.SW;KHQ=!"ZA:ET(>)Y7,C)V&=$5A+B_!!1]F_
M0NQDPE6SZ[-&N`C"RLK4?8TNE!$WBK_T%.ML1#L+Q5K1J%0$H:E<!19<R%:0
M\2BN>KLI&R36:B1T5B]*;<28GW?0=EQ*(JS^FB#7L*'RK=1L1G>+"N8G"50'
M3'`AZ]Z?,<=1N>!G.K3W<W%7L&%!]&Q(-US[,-;%HMQFS>@H,`U'X/I`A'!L
MF`(8!9D^"X8?F^2&37E@PK@,'YUM3,>J^H@NE.C$0S`TM;GC\M")L15];/$1
M#8T%4;V3D\"YX6:%$D-&;%2QB&XXT:GXPG?K5.Q68`W5AIF*`HJC)(:@A'5.
M8.$7PF7OP].)J>4Q753-H<\:8Q:$CZ>,\QU;@5UY\$@/O5R9-=(L,)=XQM##
M1MVQ\0(79QLOR(BMP,4*+/>J4,)<C)7V$D>6!)6IS4)0Z-C2<'B8SG>+7.6B
M(.>Z(%=-1;(8`N,L,/L1BPA3*C:^*$/#,M!X@:L5L2O`2J[,FO&J8!)$P^$;
MA+G$*UAUK(H@U_U%V'(<E1M)="IFT^-$WS&?`AAS55OZE.GW-#ML!6%^;E3:
M0\>V**"AV8YMF36^,/U&M<**(2/V5K3'CET9&A;@U<%X+JM66"LLTD*!ZM;A
M2:$)\K:;:J/`PD>$N28PUG17$*X^PO-$G4[&=G$SP]<[]3W`5P.];YB\X\$W
M)\/^]G@<[ZYFX_UQG-SEERB_^JJ(P#Z/7^_\+WW_^X_C[NUX.GSXJ[_]_;/O
M?Y4*L;W_M\X/>,S5X9?WO_^#[W\__]0+X-I40V8Y_-MV<SHNA]]@_(MWP+5#
M>Q&\N__P`[W_O;_=CO5-[_.'M^-ALV;>'V#=[`Z5P%?BE?!L=AKI%1W?F:LE
M_`=KN%ZN(*O5'Z)Y:60O0Z2;V>O;]4CO/.<@]*6]*7;1E;I9TGM"Q/IF^<=V
M86X6"[P<:C=!.NLH+NR$I<]5_226^3A+G[-N9NLWX_I'?,^*[Q*Q<%O,-G"E
MZ/7N?H?Y:K/M?ESI)?^C^K:T7NLE_3]U,/Z#7,F9;TTR$:'<^>[AM!Q.A\V"
M7AZC?Y_./MM#'VA;[0^[_P`QF]WVQ>LYF8#O@U\/^V%SI)?JSW?;<;C=WI&<
MX7K8+X8O!N[)5W\_J/)&&/6NAU\/FG1C0CZ<YFM,VCSPEYN;!;^MQE?-U'$V
M>X?^*",I4;=:[XYS&[Z@"%[=C3\L(#CJO>-F.[VWF+TS7<0U>@)\!JU6M&H*
MK]GL<]X5WNG5[?:'^W'^SBR^;**6M%96Z]O[]6\W;\:[P^W]U]Q++]_9Y:3O
MK`SO]OOC_"^@@<NNX6B\RHOA[V`]X-ZS1@>T3/]L<U@#54ZM7":Z7`)Q;3[!
MHX"0BY!\0KP]ANM$)P?PIZF+V7N<7;T"]XV'(]SZ]\WIS9R%+:HW[L?M_/T"
M9]RTMOE[6!<P[@OI(,'1>[$"T@\@"CVB5^.UDTQ]D7G]<>H,@GHZ.,YO\X]E
MI^*8B7\A3/<GS2M$KRC6GX^W!]B1YZT?.V5_,I_JQDNS=K6?ZLH)$+NY3ROF
M;K.C'.2_[M^,!Q$U.(GT;X$>PT<HL_N)5WZ_V<KN&LZ<*G@75(3<`@TIV9AS
M3@$+?6C0*NEH,H0N'=3IMS4:'0Y8UGOX/T`U![$-..IDLT^9N#E$9[-**KH2
M1_?G<72\'$*8O]YKS%OOG\JD@JMN#C>N!UAW(TH2L7!=ADPK3"E<8SR;DP@H
M4S>9SS)'8O*N2[.3=&J[VQS<),%0Z[W]N)>=CUXGDV-TG-1SA.-]`L]GBWZ"
M(WU()OEHK=@G<3:<#2:J%)*#\AH<#/.B7`Y:&0PO!2)!MH-+.':##>[C-H2`
M4PM&^)"0:DT*.=C@\4$PJ+(^)N5B,'P*^^_;\#C:8`Z6`S@=;IH+-^$,@2Z%
MR:(`L6<!<M_3PGO<FG!;PN2Y.\MESW83E^/J<S30Q;(/QI*_`QJ.3Q1A-&E!
MHLQ?+"JB,QS)PR,3BZ)=O<EYU4:"$DDJF;T<2"5T%Z&_D:%/OTZ%T3_'YY>+
M1TN`Q*QNWV^.93O;O"3*S73'\8$2Y.(2W?QY>JCTV7IS^-1>I90R*5N7-$28
MPU2@<L)?)?4Y0F`O,5-DZ.%SHD->KW.\\1!0.5H-G@RK'$)VUI@0TGBM*#!!
MMH8`=,%Z"-!$>0V.BDX#$0Z78-UQ_^;C*1$E>,QP(,=B3H/KZ)2#@/4*ON$:
M<I;V.GN7Z8$5614@Z/0=C5H(.ALF;D?3@6*+'"H(V>^.,K,H5,`_(+WLCM,*
M9FKJ<&[K<&8L!)!A^>"#D@]G_(5[]N90VV`C>;@_E7[3]06])BF8N]9@/.Q.
MMZ?Q%44)EJ]\E_?UFG\YA,A=B\^'_6;]XW![^'YS.MS"86-_FK$;T9HM_K+D
M/921U'FU/NR.QTLB*1/,R3O78AB+PF"G3NY0A,\PHA]E&.ZQ%`:0XU_M[Q\F
MG7FHW^+6RQJ6`\I;#M('3'V[V?Y$[O64+(O';D_Q8)?<'#%9CEFUVE'NAFTL
MYS,B!OQ3Y53U%%;G.RNNBK*UJH_;<)$:SKCX9`>##1_&W!Y^>(>'"4V'C=J"
M!2#^ANC^?G>Z>MK.-.\VXY]>K=_>46Z\2K`EJZOEE5;+A-^P9IR_@G3Y63T4
M;Q]0P.RS`E;X8_ZRFKNLL='&7@Z:>*;!X-KM#G>_WT`1!GGA/,'T);>D976X
MO=L\''^-0_2+";6*?$D9EXH[CH\5K`(^@W6]IS>'<;S;O(7P@HT"QPNH#AO"
MCQ]N\5.O`9\57<V>_/+Y?_;W'_0+]3^/#GS"%]S'GO\9N*V'(>K@8(?3]/Q/
MV^!^>?[WM_A\_=V+;\H\#,,+]P6`[X:O\;%L^=,L!U6F@F^[TOA]K5?6Q5FG
MV6%XIB_1HC4H%:I6_*;3D:`YH)G'M`#'04LT'S710I8T"(YG]A(M^4!&1A59
MF\Z"!O>>N4LT$S+3<N*Q3;2!J!<7M4'QQ:[1B;1%[P4-VI[YQ[0('O2DC1YX
M7QNHK)R@@<@7%VE!6:)E-/;:@I%B`C2(>DY]IS3HGA#KE;>Z3(`2-#32/*89
MT!:9%E';V;QI:'MN']/T"@IJTF*]?:P-_ZCKF;NDS:.1^(R%7'-.T^C)X4)P
M*>.8AD%V;J2!MF?A$@V*?#;2A0O:P/`7YN+8$H=R,)>,=+AP+M)4H*#*^I(V
M#&7_S:6%4Z)#VTO:8`#?#EV;:3$)13HND%5RCK7+96HPN/2W%VB>YBNNC#8\
M[5EHLPI=<HD6<&Y`FU&N)`<Q`1AP,I>8-MW&*Z+!X86NE1=&XES*7")H%%P>
MQIA(FU-BO5D[S26FK0"OV"41)P*C)HOU9MTTEPB:UN021=,.VHRD^6DNZ4;R
M"H!4A%[#L6E)"]-<TFDALK8</1FI@YR`.,TEID6)=H:T1<O+55DY`6F:2SJ-
MDX$';6RDCK;3<%7(7-*-I-=RX(H059D`28O37&+:"M`QL6N<+C1AI$O37-)I
M%NVG%*0NT/(TEW0C/4XI/NI-[C$-`T_FDC[=47&>C+;2Q`1X/<TE(DHB[S@J
MEBA11M#,-)=4&AZOV9/.Y\>>Q+7XXO+8/"<ZZ_(%(]TTEU2:7RE*=!`=ESSI
MXS27V#8!4:&(M'(47!Z3NJ"E:2[I-*M02UXEW)1P>\Q!T/(TE]@V-D<%1H(<
MDCDE)9$4@IKFDD+#O3HXKA@BA[0+8@*"GN:23E.1QQ12X.6J1*6`?XDL<TFE
M@>,3:TO8!9:KRY)FI[FD:],T7PF2>*C)0=#<-)=TFL,I)2ULK#.2YJ>YI-+J
M,DV0+^FOZ%<^B`D(89I+^M@<IE`T,CJ:-YAV08O37&);*%O*QE"?4.WEL6+H
M-`QSF4NZ-J,2IR+:37$"A)%HN,PEG18#YWXNV9`FTBO.J<PEW4A=RE!C3*&)
M4,;T)'-)#V6O>(.RC2:""RLCF4LZ+>=<\F1^K`W#7.:23M..71+)HT@3^QM6
M2"\NN`1K+DZKBBJ]<YJ9YI*^WHSFA&>U*C017%A^RUQB6\)SI7I-39L<6YCF
M$E>,Q)03>75K_HY>1$F*TUQ2::YLPAGVM<"I2$YW2M-<XMH*"(9#F5="A@F1
M1N9I+NG:N+*#R'><BI)<.%E-<TFG:<K];F70)1I?JX@HR7J:2SK->*Y+@K-$
MRW)LV4QS2:7!^C(<5!KS!FJS8KJSG>:23JLYQ$=V39+E3';37%)IH=0E4#QA
MS8]&1K%U9#_-)7UL.;I2(')=F9QT29CF$N%)S?DQ9,Y@.74C-=:8,I?TZ4Z*
M-]^$E1ZZ1"1SC6TREW2:IZT#=AZL"H@6!"U,<XEKZ\V47*)2U2:-C--<XEKF
MXM4-51#NIH]H:9I+.LUE/F3:1I-&YFDNZ4;&Q,5N-/:Q2[2:YA*Q3#W7D0$3
M'M&BH.EI+G$M<]&S&_1D+BM`Y$F-1PJ92US+RL'P.<!:W[1]][M_$=18E^K7
MPV_Z\0CW;W0'Q*?%RD!#IG928ZI+=4K#EV>*5D,.K%&+Q*`QQ'FIGM,4%O/7
M<'2GXU&&?4'0Z%AK+VK#0]@U;,0X&9K>$0N:KDOUT=CHZ!=A#C./+8DIQXEY
M82_1,B8&T,9']PQ93+C$V+I4'XT-#0&:Q6^-JT),.1UK+]+H$'V=N5K04`AG
MH0U7""_51R[!%8Y'=E?C4@F:K4OU?&P1'8_?)ER@N;I4SVGTA`;"6CE5ICL)
MFJ]+]=Q(&YGF0BPT,6\VU*4ZH>'S']2"NT^J-+'"Z0P7+M$\GBXPK^1*$Q.`
M9Z#G\3%-PP9LV,A8:7(",)33)9JC9R:PU^7JDI[0-0[@>;YDI,(<^5_M75^/
MVS@.?]Y\BKQM@DL"_;7LOA7;76"!W?9PU[?!('`2SS2XF:1(,MW]^">2DD7&
M[A0]].6`^"5C6Q0EFOJ1HD3-$L(2S7#@P+SN?37&S</P!!33=M@WUWL+1&9Z
M&%(0Q8@ZZ0(98E.Q[^9Z;T&2.8H:1IUT$!^)9)7E?>N]!4D6/2'011@XT%B8
M\;.H(7KN'((*-^02?^N:N,7^%#*O)`25OC4`!B!)30Y4X#KIM80@UC?M$/!4
M(!,9308C,Q*"3&_V#4Q*P(X[\AH"B]!H;R4$,3*<:$9NGOPNH27>20@RO;?@
MP(&*O\93$$0YWD@O(:CTK4)5#LF#;40X0OM*0E"19`!7!LB2EJB*?;?*20@R
M<NJ'EM4/'%]=>0E!14ML$Q+P9:>>B:2J)`0Q5:[K9)";$6Y!0A#32:@17!M?
M)RUA`Z>J)029'KD:,,1`9C,9T\FJD1!4R-!S7:9H[[5.@G0Y!!4RVY`J:YM5
MF>$D.(X<@@R;UOHT<#(9PTGPI3D$%9%4MD[N=AB*!%3A_;!OX)MX6K!P:9XC
MR6H)07E6%6"K![G;."N.C>4Z&1H)09FLIM!8)`N67!O#)M&X=,(AR/8#A\QB
M'-6&@OQQ)L+(M(2@0N;`00+E"M3(X-@(J(V$(-N/[L:3-:W#,$2F:RLAB#42
MK$K\$$U#WEY@\QQ=.PE!120*QEDDRU-W6_.^>0E!Y0,$P!`8":@E3BI774D(
M*MQJFT324&BS\KR104)0X58#J($U19R,WR_P[U9+""IDU+<(ZFGJ7K&U'-U4
M$H)8/`@-E4KQ(!G'TTV0$,0".X'6`GQ=#[];4TL(LCT$^4`140H#7I,U$H(*
MF6O22I^B^()EL2ZCE(0@'GVJ\5>!G4:RFI%I"4%,).B]1FY^C)N1$%2X53YY
M07TC.3<K(:CTK7')P[-C9$Y"D"WA"$NHC/HW:*27$%2B3ZHB&]!X-R0#6.(0
ME"=C#:US@,^,$*3!TV-D6D*02X*,=DV1BZU31+MVBI$9"4&%FTKK5-[00J$Q
MO)%60E`)$%A#L>4&9D\03O*<FY,05,(139/7%]U@7=AH+R'(%50.M*Q%5A4V
M4W*15!*"7`\*M')D8/<IC3NV*&-TD!!4N`5-X3^3XN>665,#'^?#*#?K*!QA
M?$C30*9<.,GT8V04[85U#PI+6!9(Q0,Q/XR285QD"NX,D9O`&@DN*H>@$E>@
M&*5*,WZYOFA,(R'(L1F5IP"!'BY+&K!Y'()*J`4M#NBBK4;(M(2@$HZHTN)G
M%3)><C(C(:B08<P'\C)J-T)F)02QF$D:.%5=CS3220CB02M/`\>&$6Y>0E`9
M`4:%M'A=C9!5$H)89*VAQGGC1AH9)`25OJDTR:Q]TW.#4,NO[]_]^>Z/VTZN
MV_XOW/]%B44_/OGSF_F?TV#!]<WG/P.F:=@&=MO_]7]]_B_J4RXQR%?[$0<%
MERI^.YZVW6_[[JEP?/N\Z4Z-*V^&7<+Z]MNWV_WNO+A.9J7B_][NN\-E_[#?
M#G)9/Y[:PYF.XHTMO<YEY1FKNPX.S(TCK-W^9W,\="B`\^S</3TL('/\>%KC
M@;V+Z:'[^T)_I_R+\\MC>TK'\L+]8_O\W++[S7J7I`H/:0@+.>-VX+N4"GG^
MKM*P]QG;@AN?63OOEOK-_6)*[2]-OGNC[V$?--"=NC-0X?.4L$E=R8?WQ@(K
M?%)X[8$B"O6QFV%6`..(B9_S!>W?QGI@DS>68B+#0M""$V;U4L&[_?WDI\US
M^QE/)3ZOML?#^7)ZV1+>?3R^O1R?9_/)3[MV()(9MO+#VBV0\I?RA\U_&,JJ
M6-!)`>O+<=W&^F)-P#'6RB2>>[YK7^?F<]U^C"VF(GV5W>8K[%!M8K%=BYFC
ME#&:G_=%\8Z^!]<3KC7_X!HG2ZX>NT,7AVWW9_MY1K6<(%NJ%'C8'W:_H/CC
MU[K$`5"JNE/WJTU*&CJ9[Z(R(M4(TTD,II(HW-R/#S0^2/7O3L?+@,._,&?B
M_'9S?+E\@P70KZ*\SK]G981"44T@4Q*5^2%6/\"?F:CDX6'UI7V"H\I9'?BI
M^UMLLR.JXZ=U+!RKC7013MJ7)]3?/RCCH=0`>D$9`)B97NJZ3L=*%4*F1),3
MG*-:I(2))*GC=CM@2;HS^'2+Z4!HK[*/-2]QZJ@ERU-W>3D=$LY-)E\Y]?W5
MP]W_MV/=)^W3$PPF*$YTT`MX0O*$]_F,]/*^R'O2KA^ZO]:Y2"Y]]\;=3S:*
M/U'WD(JR4:M6KZAOL7/OUSJEG+30*7Q+MR;=&G;<2%_&R#*ZSV:&NW-*\FKU
MG"<VJP`'K:B5QA-P=+=43F8WIT36T:H,JVJ9ZUKFRI;?K"T.XU@AJO;/H$XS
MP-_IZMU'O]:K37ON5E$2BRE_BEH4QX:>TVDJH)BGZ:L5]22RIES_',Z[5ZMS
M=_EPVG6G69S`_9`&+B"J\<-:F6J;IV,+KK[H9`?:I(L9@02L_>$1N<3O1*9B
M,?VBYZFL^8ZR]I6R\UY6N0D13)&!?&'A()RK%R`/>DDR$;5_139H\GXNB6R)
M+0R?S\?S'L?\/%JF_#@+"30U/<ML_MF7+RJ:OU9)X-Z`5-EPOM-Q"+=%U`4&
M-RJ:XJ+<L5"TKY`Q3`7)=#T_'P_H9<32L3"\;F4N&]!!*E<+2!1;;6!/Z7#L
M0+%/[1EEM2&3EAYA]6UOY5CYX^53=^HI(NNHC8/WU+K4-,@1GHSZ)QDAT;3T
M-YK?&'X#!ZIHH\3Q&5E!VFQ!4ZL5?W75(S5&5MIRU6M.7]J8"XWUK'@7^$'Y
MEX?^;/"LW6;E64<8,R9AT@:4HJCCJFVI\!S'1WG2&XAE;"ZHH5J=NN?CESA?
M``?H2ZS@XS%K\JP=Q=1^?!FL(4TX*!44S-+(-*0W9HLXJ5ADVP;_D8,.L2K$
M=YB3X(R"_'PXW0OV9(;:&6N\\B;=&ZN\AX1ENJ]-M/$V-%;3/Q9A7KX;)IE?
MM3@Z[W$XB';`N2[I!(X*#EL8.9;K%M2Y7;?K=MVNVW6[;M?MNEVWZW;=KI'K
)OUNW3DH`>```
`
end
--- old/Geometry.py Tue Feb 29 05:29:58 2000
+++ Geometry.py Thu Dec 21 16:59:02 2000
@@ -11,7 +11,9 @@
and lattice objects, which define a regular arrangements of points.
"""
+import Units
from Scientific.Geometry import Vector
+from Random import randomDirection
import Numeric
# Error type
@@ -19,6 +21,8 @@
# Small number
eps = 1.e-16
+eps_radians = 1.e-6
+min_meaningfull_distance = 1.e-3*Units.Ang
#
# The base class
@@ -49,6 +53,10 @@
def hasPoint(self, point):
return self.distanceFrom(point) < eps
+ def enclosesPoint(self, point): # subclasses that enclose a volume
+ # should override this method
+ return 0
+
_intersectTable = {}
#
@@ -70,12 +78,116 @@
"""
def __init__(self, center, radius):
- self.center = centerLattice
+ self.center = center
self.radius = radius
def volume(self):
return (4.*Numeric.pi/3.) * self.radius**3
+ def __repr__(self):
+ return 'Sphere(' + `self.center` + ', ' + `self.radius` + ')'
+ __str__ = __repr__
+
+ def enclosesPoint(self, point):
+ return (point - self.center).length() < self.radius
+
+ def coordList(self, no_nulls = 0):
+ result = []
+ from Numeric import arange
+ for normal in [Vector(1,0,0), Vector(0,1,0), Vector(0,0,1),\
+ Vector(1,1,1)]:
+ v0 = normal.cross(randomDirection())
+ if result and not no_nulls:
+ result.append(None)
+ for u in arange(0, 2*Numeric.pi, Numeric.pi/8):
+ v = rotateDirection(v0, normal, u)
+ result.append(self.center + self.radius*v.normal())
+ return result
+
+
+class Prism(GeometricalObject3D):
+
+ """Prism must be constructed from a convex Polygon which is a sequence of
+ coplanar points, ordered clockwise as seen from inside prism."""
+
+ def __init__(self, polygon, height):
+ self.polygon = polygon
+ self.height = height
+ self.vector = height*polygon.normal
+
+ def volume(self):
+ raise 'not implemented'
+
+ def __repr__(self):
+ return 'Prism(' + `self.polygon` + ', ' + `self.height` + ')'
+ __str__ = __repr__
+
+ def enclosesPoint(self, point):
+ return 0#??(point - self.center).length() < self.radius
+
+ def coordList(self):
+ result = polygon[:]
+ result.append(polygon[0])
+ for pt in polygon:
+ result.append(pt + self.vector)
+ result.append(pt)
+ result.append(pt + self.vector)
+ result.append(polygon[0] + self.vector)
+ return result
+
+class Cylinder(GeometricalObject3D):
+
+ def __init__(self, center, radius, vector):
+ self.center = center # center of base
+ self.radius = radius
+ self.vector = vector # <center of top> - self.center
+ self.height = vector.length()
+
+ def volume(self):
+ return 2*Numeric.pi*self.radius*self.radius*self.height
+
+ def __repr__(self):
+ return 'Cylinder(' + `self.center` + ', ' + `self.radius`\
+ + ', ' + `self.height` + ')'
+ __str__ = __repr__
+
+ def enclosesPoint(self, point):
+ center_line = LineSegment(self.center, self.center + self.vector)
+ pt = center_line.projectionOf(point)
+ if pt is None:
+ return 0
+ return (point - pt).length() < self.radius
+
+ def coordList(self):
+ return Circle(self.center, self.vector, self.radius).coordList() + \
+ Circle(self.center + self.vector, self.vector, self.radius).coordList()
+
+#
+# All but last point of pyramid must be coplanar
+#
+
+class Pyramid(GeometricalObject3D):
+
+ def __init__(self, points):
+ self.points = points
+
+ def center(self):
+ return reduce(__add__, self.points) / len(self.points)
+
+ def boundingBox(self):
+ min = self.points[0].array
+ max = min
+ for pt in self.points[1:]:
+ r = pt.array
+ min = Numeric.minimum(min, r)
+ max = Numeric.maximum(max, r)
+ return Vector(min), Vector(max)
+
+ def volume(self):
+ base = Polygon(self.points[0:-1])
+ return base.area() * base.planeOf().distanceFrom(self.points[-1])\
+ / (len(self.points) - 1)
+
#
# Planes
#
@@ -108,11 +220,15 @@
self.normal = (v1.cross(v2)).normal()
self.distance_from_zero = self.normal*args[1]
+ def __repr__(self):
+ return 'Plane(' + str(self.normal*self.distance_from_zero) + ', ' + str(self.normal) + ')'
+ __str__ = __repr__
+
def distanceFrom(self, point):
return abs(self.normal*point-self.distance_from_zero)
def projectionOf(self, point):
- return point - self.distanceFrom(point)*self.normal
+ return point - (self.normal*point-self.distance_from_zero)*self.normal
def rotate(self, axis, angle):
point = rotatePoint(self.distance_from_zero*self.normal, axis, angle)
@@ -148,6 +264,10 @@
self.axis = axis.normal()
self.angle = angle
+ def __repr__(self):
+ return 'Cone(' + `self.center` + ', ' + `self.axis` + ',' + `self.angle` + ')'
+ __str__ = __repr__
+
def volume(self):
return None
@@ -176,9 +296,113 @@
self.normal = normal
self.radius = radius
+ def planeOf(self):
+ return Plane(self.center, self.normal)
+
+ def __repr__(self):
+ return 'Circle(' + `self.center` + ', ' + `self.normal` + ', ' + `self.radius` + ')'
+ __str__ = __repr__
+
+ def volume(self):
+ return 0.
+
+ def distanceFrom(self, point):
+ plane = self.planeOf()
+ project_on_plane = plane.projectionOf(point)
+ center_to_projection = project_on_plane - self.center
+ if center_to_projection.length() < eps:
+ return 0
+ closest_point = self.center + self.radius*center_to_projection.normal()
+ return (point - closest_point).length()
+
+ def coordList(self):
+ result = []
+ from Numeric import arange
+ v0 = self.normal.cross(randomDirection())
+ for u in arange(0, 2*Numeric.pi, Numeric.pi/64):
+ v = rotateDirection(v0, self.normal, u)
+ result.append(self.center + self.radius*v.normal())
+ return result
+
+class Polygon(GeometricalObject3D):
+
+ def __init__(self, pt_list):
+ assert(len(pt_list) >= 3)
+ self.points = pt_list
+ v1 = pt_list[2] - pt_list[1]
+ self.normal = v1.cross(pt_list[1] - pt_list[0]).normal()
+
+ def planeOf(self):
+ return Plane(self.points[0], self.normal)
+
+ def __repr__(self):
+ points_repr = `self.points[0]`
+ for pt in self.points[1:]:
+ points_repr = points_repr + ', ' + `pt`
+ return 'Polygon(' + points_repr + ')'
+ __str__ = __repr__
+
+ def __getitem__(self, index):
+ return self.points[index]
+
+ def center(self):
+ sum = Vector(0,0,0)
+ for pt in self.points:
+ sum = sum + pt
+ return sum / len(self.points)
+
+ def area(self):
+ assert(len(self.points) <= 4)
+ if len(self.points) == 4:
+ return (self.points[1] - self.points[0]).length() *\
+ (self.points[2] - self.points[1]).length()
+ v = self.points[1] - self.points[0]
+ side1 = Line(self.points[0], v)
+ ht = side1.distanceFrom(self.points[2])
+ return 0.5 * v.length() * ht
+
def volume(self):
return 0.
+ # assumes convex polygon
+
+ def pointNearest(self, point, proj_only = 0):
+ plane = self.planeOf()
+ project_on_plane = plane.projectionOf(point)
+ n = len(self.points)
+ is_inside = 1
+ nearest_dist = 1.e9
+ for i in range(n):
+ edge = LineSegment(self.points[i], self.points[(i + 1) % n])
+ v = project_on_plane - self.points[i]
+ if v.cross(edge.direction) * self.normal < 0:
+ is_inside = 0
+ near_pt = edge.pointNearest(point)
+ d = (near_pt - point).length()
+ if d < nearest_dist:
+ nearest_dist = d
+ nearest_point = near_pt
+ if is_inside:
+ return project_on_plane
+ if proj_only:
+ return None
+ return nearest_point
+
+ def projectionOf(self, point):
+ return self.pointNearest(point, 1)
+
+ def distanceFrom(self, point):
+ return (point - self.pointNearest(point)).length()
+
+ def coordList(self):
+ return self.points
+
+ def edgeVectors(self):
+ edges = [self.points[0] - self.points[-1]]
+ for i in range(1, len(self.points)):
+ edges.append(self.points[i] - self.points[i-1])
+ return edges
+
#
# Lines
#
@@ -213,9 +437,83 @@
d = d - (d*self.direction)*self.direction
return point+d
+ def perpendicularVector(self, plane):
+ return self.direction.cross(plane.normal)
+
+ def __repr__(self):
+ return 'Line(' + `self.point` + ', ' + `self.direction` + ')'
+ __str__ = __repr__
+
+ def Is_Infinite(self):
+ return 1
+
+ def asVector(self):
+ return self.point2 - self.point
+
def volume(self):
return 0.
+ def segmentNear(self, other): # for use by intersection code only
+ pt1 = self.projectionOf(other.point)
+ v = pt1 - other.point
+ if v.length() < min_meaningfull_distance:
+ angle1 = 0
+ else:
+ angle1 = v.angle(other.direction)
+ if abs(angle1 - Numeric.pi/2) < eps_radians:
+ return None
+ if angle1 > Numeric.pi/2:
+ angle1 = Numeric.pi - angle1
+ d = self.direction * (0.1 + Numeric.tan(angle1) * v.length())
+ return LineSegment(pt1 - d, pt1 + d)
+
+class LineSegment(Line):
+ def __init__(self, point1, point2):
+ Line.__init__(self, point1, point2 - point1)
+ self.point2 = point2
+
+ def BisectPoint(self):
+ return (self.point + self.point2)/2
+
+ def PerpBisector(self, plane):
+ return Line(self.BisectPoint(), self.perpendicularVector(plane))
+
+ def Is_Infinite(self):
+ return 0
+
+ def pointNearest(self, point):
+ pt = self.projectionOf(point)
+ if pt is not None:
+ return pt
+ d1 = (self.point - point).length()
+ d2 = (self.point2 - point).length()
+ if d1 < d2:
+ return self.point
+ return self.point2
+
+ def distanceFrom(self, point):
+ pt = self.projectionOf(point)
+ if pt is not None:
+ return (pt - point).length()
+ d1 = (self.point - point).length()
+ d2 = (self.point2 - point).length()
+ return min(d1, d2)
+
+ def projectionOf(self, point):
+ d = self.point-point
+ d = d - (d*self.direction)*self.direction
+ pt = point+d
+ if Within(self, pt):
+ return pt
+ return None
+
+ def asVector(self):
+ return self.point2 - self.point
+
+ def __repr__(self):
+ return 'LineSegment(' + `self.point` + ', ' + `self.point2` + ')'
+ __str__ = __repr__
+
#
# Intersection calculations
#
@@ -225,6 +523,151 @@
class1, class2 = class2, class1
_intersectTable[(class1, class2)] = (f, switch)
+#
+# a number of routines used by _intersectLineLine
+#
+
+def LinePointColinear(line, pt):
+ return line.distanceFrom(pt) < eps
+
+def Same_Dir(dir1, dir2):
+ # for 2 colinear vectors, FALSE if in opposite dir, else TRUE
+ return (dir1[0] >= 0) == (dir2[0] >= 0)\
+ and (dir1[1] >= 0) == (dir2[1] >= 0)\
+ and (dir1[2] >= 0) == (dir2[2] >= 0)
+
+def Within(line_seg, pt): # expects point colinear with line
+ v1 = pt - line_seg.point
+ v2 = pt - line_seg.point2
+ if abs(v1 * v2) < eps:
+ return 0
+ return not Same_Dir(v1, v2)
+
+def Overlap(l1, l2): # expects colinear lines
+ if l1.Is_Infinite(): return l2
+ if l2.Is_Infinite(): return l1
+ if Same_Dir(l1.direction, l2.direction):
+ l0 = l1
+ else:
+ l0 = LineSegment(l1.point2, l1.point)
+ in1 = Within(l0, l2.point)
+ if in1:
+ p1 = l2.point
+ else:
+ p1 = l0.point
+ in2 = Within(l0, l2.point2)
+ if (not in1) and (not in2):
+ if Within(l2, l0.point):
+ return l1
+ # no overlap; do they touch?
+ if (l0.point - l2.point).length() <= eps or\
+ (l0.point - l2.point2).length() <= eps:
+ return l0.point
+ if (l0.point2 - l2.point).length() <= eps or\
+ (l0.point2 - l2.point2).length() <= eps:
+ return l0.point2
+ return None
+ if in2:
+ p2 = l2.point2
+ else:
+ p2 = l0.point2
+ return LineSegment(p1, p2)
+
+def _intersectLineLine(l1, l2):
+ if LinePointColinear(l1, l2.point) \
+ and Numeric.sin(l1.direction.angle(l2.direction)) < eps_radians:
+ return Overlap(l1, l2)
+ if l2.Is_Infinite():
+ l2 = l2.segmentNear(l1)
+ if not l2:
+ return l2
+ if l1.Is_Infinite():
+ l1 = l1.segmentNear(l2)
+ if not l1:
+ return l1
+ x11 = l1.point[0]
+ y11 = l1.point[1]
+ x21 = l2.point[0]
+ y21 = l2.point[1]
+ v1 = l1.asVector()
+ x12 = float(v1[0])
+ y12 = float(v1[1])
+ v2 = l2.asVector()
+ x22 = v2[0]
+ y22 = v2[1]
+ num = x12*(y11 - y21) + y12*(x21 - x11);
+ denom = x12*y22 - y12*x22;
+ if abs(denom*1.01) > abs(num):
+ u = num/denom
+ elif abs(num) < eps: # (?) XY Perpendicular, using XZ projection num
+ z11 = l1.point[2]
+ z21 = l2.point[2]
+ z12 = float(v1[2])
+ z22 = v2[2]
+ num = x12*(z11 - z21) + z12*(x21 - x11)
+ denom = x12*z22 - z12*x22
+ if abs(denom*1.01) > abs(num):
+ u = num/denom
+ elif abs(num) < eps: # XZ Perpendicular, using YZ projection
+ num = y12*(z11 - z21) + z12*(y21 - y11)
+ denom = y12*z22 - z12*y22
+ if abs(denom*1.01) > abs(num):
+ u = num/denom
+ else:
+ u = 1.01
+ else: u = 1.01
+ else: u = 1.01
+ if u < 0 or u > 1:
+ return None
+ z12 = float(v1[2])
+ # set t to parametric distance along l1
+ if abs(y12) >= eps:
+ t = (y21 - y11 + u*y22)/y12
+ elif abs(x12) >= eps:
+ t = (x21 - x11 + u*x22)/x12
+ elif abs(z12) >= eps:
+ t = (l2.point[2] - l1.point[2] + u*v2[2])/z12
+ else: t = 0 # can l1 be a meaningfull line if this happens?
+ #print 'parametric dists %e %e %.2f %.2f %.2f' % (t,u,y12,y21,y11)
+ if not l1.Is_Infinite() and (t < 0 or t > 1):
+ return None
+ pt_from1 = l1.point + t*v1
+ pt_from2 = l2.point + u*v2
+ if (pt_from1 - pt_from2).length() < min_meaningfull_distance:
+ return pt_from1
+ else:
+ return None
+
+_addIntersectFunction(_intersectLineLine, Line, Line)
+_addIntersectFunction(_intersectLineLine, Line, LineSegment)
+_addIntersectFunction(_intersectLineLine, LineSegment, Line)
+_addIntersectFunction(_intersectLineLine, LineSegment, LineSegment)
+
+def _intersectLinePolygon(line, poly):
+ pt = line.intersectWith(poly.planeOf())
+ if pt is not None and poly.pointNearest(pt, proj_only = 1) is not None:
+ return pt
+ return None
+
+_addIntersectFunction(_intersectLinePolygon, LineSegment, Polygon)
+
+def _intersectLinePlane(line, plane):
+ v = vectorFromPointToPlane(line.point, plane)
+ if v.length() < eps:
+ return line.point
+ angle = v.angle(line.asVector())
+ if abs(angle - Numeric.pi/2) < eps:
+ return None
+ proj1 = line.point + v
+ proj2 = plane.projectionOf(line.point + line.direction)
+ x = proj1 + Numeric.tan(angle)*(proj2 - proj1)
+ if line.__class__ is LineSegment and not Within(line, x):
+ return None
+ return x
+
+_addIntersectFunction(_intersectLinePlane, LineSegment, Plane)
+
+
# Sphere with sphere
def _intersectSphereSphere(sphere1, sphere2):
@@ -253,6 +696,121 @@
_addIntersectFunction(_intersectSphereCone, Sphere, Cone)
+#
+# Sphere - Line intersection returns None or tuple of 2 points. The
+# 2 points can be identical.
+#
+
+def _intersectSphereLine(sphere, line):
+ r = sphere.radius
+ dist = line.distanceFrom(sphere.center)
+ if dist > r:
+ return None
+ projection = line.projectionOf(sphere.center)
+ drop_normal = projection - sphere.center
+ if drop_normal.length() < 1.e-9*r: # line goes through center
+ return (sphere.center + r*line.direction,\
+ sphere.center - r*line.direction)
+ scale = Numeric.sqrt(r*r - drop_normal.length()*drop_normal.length())
+ normal2 = scale*line.direction
+ return (sphere.center + (projection + normal2),
+ sphere.center + (projection - normal2))
+
+
+_addIntersectFunction(_intersectSphereLine, Sphere, Line)
+
+# Circle with Circle
+
+def vectorsAreParallel(v1,v2):
+ return v1.cross(v2).length() < eps
+
+def Circle_Intersect_2point(c1,c2,d):
+ # must have c1.radius >= c2.Radius; must have d > 0
+ dr2 = (c1.radius - c2.radius)*(c1.radius + c2.radius) # r1*r1-r2*r2
+ d1 = dr2/(2*d) + d/2
+ v = (c2.center-c1.center) # direction from center to c2
+ v1 = (d1/d)*v
+ # p is projection of intersections onto line connecting centers
+ # Point p(c1.center + v)
+ scale = Numeric.sqrt((c1.radius-d1)*(c1.radius+d1))
+ v2 = v.cross(c1.normal).normal()
+ return (c1.center + (v1 + scale*v2), c1.center + (v1 - scale*v2))
+
+def Circle_Intersect_Noncoplanar(c1, c2):
+ plane1 = c1.planeOf()
+ plane2 = c2.planeOf()
+ intersect_line = plane1.intersectWith(plane2)
+ if intersect_line is None:
+ return None
+ assert(intersect_line.__class__ is Line)
+ proj1 = intersect_line.projectionOf(c1.center)
+ dist_proj1 = (proj1 - c1.center).length()
+ scale = Numeric.sqrt(c1.radius*c1.radius - dist_proj1*dist_proj1)
+ pt1 = proj1 + scale * intersect_line.direction
+ pt2 = proj1 - scale * intersect_line.direction
+ intersect_pt1 = abs((pt1 - c2.center).length() - c2.radius) < min_meaningfull_distance
+ intersect_pt2 = abs((pt2 - c2.center).length() - c2.radius) < min_meaningfull_distance
+ if intersect_pt1 and intersect_pt2:
+ return (pt1, pt2)
+ if intersect_pt1:
+ return pt1
+ if intersect_pt2:
+ return pt2
+ return None
+
+def _intersectCircleCircle(c1, c2):
+ if not vectorsAreParallel(c1.normal,c2.normal):
+ return Circle_Intersect_Noncoplanar(c1, c2)
+ d = (c1.center - c2.center).length()
+ r2 = c1.radius+c2.radius
+ r0 = abs(c1.radius-c2.radius)
+ if d/c1.radius < eps and d/c2.radius < eps:
+ if r0 <= d:
+ return c1 # same circles
+ return None
+ if d < r0:
+ return None # 1 circle entirely inside other
+ if abs(r2/d - 1) < eps: # single point intersection
+ pt = c1.center + c1.radius/d*(c2.center - c1.center);
+ return (pt,pt)
+ if r2 < d:
+ return None
+ # if get to this point then should intersect at 2 points
+ if c1.radius > c2.radius:
+ return Circle_Intersect_2point(c1,c2,d)
+ else:
+ return Circle_Intersect_2point(c2,c1,d)
+
+_addIntersectFunction(_intersectCircleCircle, Circle, Circle)
+
+def vectorFromPointToPlane(point, plane):
+ d = (plane.normal*plane.distance_from_zero - point) * plane.normal
+ return d * plane.normal
+
+def distanceFromPointToPlane(point, plane):
+ return vectorFromPointToPlane(point,plane).length()
+
+def _intersectSphereCircle(sphere, circle):
+ vector_to_circle = vectorFromPointToPlane(sphere.center, circle.planeOf())
+ projection = sphere.center + vector_to_circle
+ d2 = vector_to_circle.length()
+ if sphere.radius < d2 + eps:
+ if abs(sphere.radius - d2) <= eps:
+ return projection
+ return None
+ radius2 = Numeric.sqrt(sphere.radius*sphere.radius - d2*d2)
+ circle2 = Circle(projection, circle.normal, radius2)
+ return _intersectCircleCircle(circle, circle2)
+
+_addIntersectFunction(_intersectSphereCircle, Sphere, Circle)
+
+def arbitraryNormal(v):
+ if abs(v[0]) < eps:
+ return v.cross(Vector(1,0,0))
+ if abs(v[1]) < eps:
+ return v.cross(Vector(0,1,0))
+ return v.cross(Vector(0,0,1))
+
# Plane with plane
def _intersectPlanePlane(plane1, plane2):
@@ -294,6 +852,34 @@
_addIntersectFunction(_intersectCirclePlane, Circle, Plane)
+# Cone with Cone - special purpose routine, returns 0 to 2 vectors
+
+def _intersectConeCone(cone1, cone2):
+ if (cone1.center - cone2.center).length() > eps:
+ raise GeomError, "Not yet implemented"
+ if cone1.angle + cone2.angle < cone1.axis.angle(cone2.axis):
+ return None
+ # intersection is a line?
+ if abs(cone1.angle + cone2.angle - cone1.axis.angle(cone2.axis)) <= eps:
+ normal = Line(cone.center,cone1.axis.cross(cone2.axis))
+ return rotateDirection(cone1.axis,normal, cone1.angle)
+ if abs(cone1.angle - Numeric.pi/2) <= eps:
+ if abs(cone2.angle - Numeric.pi/2) <= eps:
+ return _intersectPlanePlane(Plane(cone1.center, cone1.axis),
+ Plane(cone2.center, cone2.axis))
+ return _intersectConeCone(cone2, cone1)
+ cos1 = Numeric.cos(cone1.angle)
+ radius1 = cos1*Numeric.tan(cone1.angle)
+ circle = Circle(cone1.center + cos1*cone1.axis, cone1.axis, radius1)
+ plane = Plane(cone2.center + Numeric.cos(cone2.angle)*cone2.axis, cone2.axis)
+ result = _intersectCirclePlane(circle, plane)
+ if type(result) is type(()) and len(result) == 2:
+ r = (result[0] - cone1.center, result[1] - cone1.center)
+ if r[0].angle(r[1]) < eps_radians:
+ return Line(cone1.center, r[0])
+ return r
+ return result
+
#
# Rotation
#
@@ -301,12 +887,41 @@
s = Numeric.sin(angle)
c = Numeric.cos(angle)
c1 = 1-c
- axis = axis.direction
+ try:
+ axis = axis.direction
+ except AttributeError:
+ pass
return s*axis.cross(vector) + c1*(axis*vector)*axis + c*vector
def rotatePoint(point, axis, angle):
return axis.point + rotateDirection(point-axis.point, axis, angle)
+#
+# Find the angle by which one would have to rotate vector v1 about axis
+# so that when projected on to the plane to which axis is normal, it
+# would be colinear to v2. Mainly for comparing two bonds that you are
+# trying to position about the same atom.
+#
+
+def findRotationAngle(axis, v1, v2):
+ plane1 = Plane(Vector(0,0,0), axis)
+ p1 = plane1.projectionOf(v1)
+ p2 = plane1.projectionOf(v2)
+ angle = p1.angle(p2)
+ if p1.cross(p2)*axis < 0:
+ angle = -angle
+ return angle
+
+#
+# find difference between 2 directions. With default, result is -180 to +180
+# degrees. Use modulo = pi to produce a -90 to +90 range for undirected lines.
+#
+
+def angleBetweenDirections(angle1, angle2, modulo = 2*Numeric.pi):
+ r = angle1 - angle2
+ while r > modulo/2: r = r - modulo
+ while r < -modulo/2: r = r + modulo
+ return r
#
# Lattices
--- old/PDB.py Tue Feb 29 05:29:59 2000
+++ PDB.py Tue Dec 19 11:02:20 2000
@@ -170,14 +170,16 @@
chains.append(chain)
return chains
- def createNucleotideChains(self, model='all'):
+ def createNucleotideChains(self, model='all', terminus_3 = None):
"""Returns a list of NucleotideChain objects, one for each
nucleotide chain in the PDB file. The parameter |model|
has the same meaning as for the NucleotideChain constructor."""
import NucleicAcids
chains = []
for chain in self.nucleotide_chains:
- properties = {'model': model}
+ properties = {'model': model }
+ if terminus_3 is not None:
+ properties['terminus_3'] = terminus_3
if chain.segment_id != '':
properties['name'] = chain.segment_id
if chain[0].hasPhosphate():
@@ -386,6 +388,7 @@
def __init__(self, filename, subformat= None):
self.file = Scientific.IO.PDB.PDBFile(filename, 'w', subformat)
self.warning = 0
+ self.atom_no = {}
self.atom_sequence = []
def write(self, object, configuration = None, tag = None):
@@ -412,6 +415,7 @@
except AttributeError: temp = 0.
self.file.writeAtom(atom_name, p/Units.Ang,
occ, temp, atom.type.symbol)
+ self.atom_no[id(atom)] = self.file.data['serial_number']
self.atom_sequence.append(atom)
else:
self.warning = 1
@@ -444,3 +448,47 @@
Utility.warning('Some atoms are missing in the output file ' + \
'because their positions are undefined.')
self.warning = 0
+
+class PDBConnectOutputFile(PDBOutputFile):
+
+ """Like PDBOutputFile, but also writes CONECT statements."""
+
+ def __init__(self, filename, subformat = None):
+ PDBOutputFile.__init__(self, filename, subformat)
+ self.top_level = 1
+
+ def write(self, object, configuration = None, tag = None):
+ save_level = self.top_level
+ self.top_level = 0
+ PDBOutputFile.write(self, object, configuration, tag)
+ self.top_level = save_level
+ if self.top_level:
+ self.writeBonds(object, configuration)
+
+ def writeBonds(self, object, configuration):
+ if not ChemicalObjects.isChemicalObject(object):
+ for o in object:
+ self.writeBonds(o, configuration)
+ else:
+ try:
+ a1list = object.atomList()
+ except AttributeError:
+ return
+ for a in a1list:
+ a2list = a.bondedTo()
+ connect_list = []
+ for a2 in a2list:
+ a_no = self.atom_no[id(a2)]
+ connect_list.append(a_no)
+ try:
+ a_no = self.atom_no[id(a)]
+ except KeyError: pass
+ else:
+ self.writeBond(a_no, connect_list)
+
+ # writeBond probably belongs in Scientific.IO.PDB
+ def writeBond(self, number, connect_list):
+ self.file.data['serial_number'] = number
+ self.file.data['bonded'] = connect_list
+ self.file.writeLine('CONECT', self.file.data)
+
--- old/ConfigIO.py Tue Feb 29 05:29:56 2000
+++ ConfigIO.py Thu Nov 30 20:12:47 2000
@@ -330,6 +330,7 @@
_file_types[format[0]](filename, format[1])
_file_types = {'pdb': PDB.PDBOutputFile,
+ 'pdb_connect': PDB.PDBConnectOutputFile,
('vrml',): VRMLFile,
('vrml', 'wireframe'): VRMLWireframeFile,
('vrml', 'highlight'): VRMLHighlight,
--- old/Bonds.py Mon Jun 5 06:32:30 2000
+++ Bonds.py Tue Dec 26 21:55:51 2000
@@ -9,6 +9,18 @@
from UserList import UserList
import Database, Utility
+from Scientific.Geometry import Vector
+from Units import deg
+from umath import pi
+
+BondOrderNotSpecified = 0
+
+def findCommonAtom(bond1, bond2):
+ if bond1.hasAtom(bond2.a1):
+ return bond2.a1
+ if bond1.hasAtom(bond2.a2):
+ return bond2.a2
+ raise ValueError,'no common atom between ' + `bond1` + ' and ' + `bond2`
#
# Bonds
@@ -22,9 +34,14 @@
if type(blueprint) is type(()):
self.a1 = blueprint[0]
self.a2 = blueprint[1]
+ try:
+ self.order = blueprint[2]
+ except IndexError:
+ self.order = BondOrderNotSpecified
else:
self.a1 = Database.instantiate(blueprint.a1, memo)
self.a2 = Database.instantiate(blueprint.a2, memo)
+ self.order = blueprint.order
if Utility.uniqueID(self.a2) < Utility.uniqueID(self.a1):
self.a1, self.a2 = self.a2, self.a1
Utility.uniqueID.registerObject(self)
@@ -35,9 +52,24 @@
return (None,)
def __repr__(self):
- return 'Bond(' + `self.a1` + ', ' + `self.a2` + ')'
+ if self.order != BondOrderNotSpecified:
+ ostr = ', ' + `self.order`
+ else:
+ ostr = ''
+ return 'Bond(' + `self.a1` + ', ' + `self.a2` + ostr + ')'
__str__ = __repr__
+ def removeAndConvertToDangling(self, remove_atom):
+ self.direction = self.asVector(self.otherAtom(remove_atom))
+ #self.delete()
+ if remove_atom is self.a1:
+ self.a1 = self.a2
+ elif not (remove_atom is self.a2):
+ raise ValueError, "You tried to remove atom %s which is not in bond %s" % (remove_atom, self)
+ self.a2 = None
+ self.__class__ = DanglingBond
+ self.constraint = None
+
def hasAtom(self, a):
return a is self.a1 or a is self.a2
@@ -49,6 +81,20 @@
else:
raise ValueError, 'atom not in bond'
+ def setOrder(self, order):
+ self.order = order
+
+ def asVector(self,start_atom = None,configuration = None):
+ if start_atom == None:
+ start_atom = self.a1 # ambiguous direction, pick one
+ pos1 = start_atom.position(configuration)
+ if pos1 is None:
+ raise ValueError, `start_atom`+' has no position defined'
+ pos2 = self.otherAtom(start_atom).position(configuration)
+ if pos2 is None:
+ raise ValueError, `self.otherAtom(start_atom)`+' has no position defined'
+ return pos2 - pos1
+
def _graphics(self, conf, distance_fn, model, module, options):
objects = []
if model == 'ball_and_stick':
@@ -87,14 +133,74 @@
Database.registerInstanceClass(Bond.blueprintclass, Bond)
+class DanglingBond(Bond):
+ def __init__(self, blueprint, memo = None):
+ if blueprint is not None:
+ if type(blueprint) is type(()):
+ self.a1 = blueprint[0]
+ self.a2 = None
+ if len(blueprint) > 1:
+ self.constraint = blueprint[1]
+ else: self.constraint = None
+ if len(blueprint) > 2:
+ self.direction = blueprint[2]
+ else: self.direction = None
+ else:
+ self.a1 = Database.instantiate(blueprint.a1, memo)
+ self.a2 = None
+ self.constraint = Database.instantiate(blueprint.constraint, memo)
+ self.direction = Database.instantiate(blueprint.direction, memo)
+ assert(self.a1 is not None)
+
+ blueprintclass = Database.BlueprintDanglingBond
+
+ def setDirection(self, v):
+ self.direction = v
+
+ def setConstraint(self, c):
+ self.constraint = c
+
+ def fill(self, atom):
+ self.a2 = atom
+
+ def __repr__(self):
+ return 'DanglingBond(' + `self.a1` + ', ' + `self.constraint` + ')'
+ __str__ = __repr__
+
+ def __cmp__(self,other):
+ if self.a1 is other.a1:
+ if self.a2 is other.a2:
+ if (self.constraint and not other.constraint) or\
+ (not self.constraint and other.constraint):
+ return 1
+ if self.direction is other.direction and\
+ self.constraint is other.constraint:
+ return 0
+ return 1
+
+ def DanglingPosition(self):
+ try:
+ return self.a1.position() + self.direction
+ except (AttributeError, TypeError):
+ return None
+
+ def asVector(self, start_atom = None, configuration = None):
+ if isinstance(self.direction, Vector):
+ return self.direction
+ return None
+
+Database.registerInstanceClass(DanglingBond.blueprintclass, DanglingBond)
+
#
# Bond angles
#
class BondAngle:
- def __init__(self, b1, b2, ca):
+ def __init__(self, b1, b2, ca = None):
self.b1 = b1 # bond 1
self.b2 = b2 # bond 2
+ if ca is None:
+ ca = findCommonAtom(b1, b2)
self.ca = ca # common atom
if Utility.uniqueID(self.b2) < Utility.uniqueID(self.b1):
self.b1, self.b2 = self.b2, self.b1
@@ -105,25 +211,56 @@
def __repr__(self):
return 'BondAngle(' + `self.a1` +',' + `self.ca` +','+ `self.a2` +')'
- def otherBond(self, bond):
+ def radians(self):
+ return self.b1.asVector(self.ca).angle(self.b2.asVector(self.ca))
+
+ def hasBond(self, b):
+ return b is self.b1 or b is self.b2
+
+ def hasAtom(self, a):
+ return a is self.a1 or a is self.a2 or a is self.ca
+
+ def otherAtom(self, bond):
if bond is self.b1:
- return self.b2
+ return self.b2.otherAtom(self.ca)
elif bond is self.b2:
+ return self.b1.otherAtom(self.ca)
+ else:
+ raise ValueError, 'bond not in bond angle'
+
+ def otherBond(self, bond):
+ if bond == self.b1: # need == rather than is for dangling bonds
+ return self.b2
+ elif bond == self.b2:
return self.b1
else:
raise ValueError, 'bond not in bond angle'
+def calcDihedralAngle(v1, v3, cv):
+ normal1 = v1.cross(cv)
+ normal2 = cv.cross(v3)
+ try:
+ angle = normal2.angle(normal1) - pi
+ except ZeroDivisionError:
+ Utility.warning('ZeroDivisionError %s %s %s %s' \
+ % (normal1,normal2,cv,v3))
+ return None
+ if normal1.cross(normal2) * cv < 0:
+ angle = -angle
+ return angle
+
#
# Dihedral angles
#
class DihedralAngle:
- def __init__(self, ba1, ba2, cb):
+ def __init__(self, ba1, ba2, cb, angle = None):
self.ba1 = ba1 # bond angle 1
self.ba2 = ba2 # bond angle 2
# cb is the common bond, i.e. the central bond for a proper dihedral
if Utility.uniqueID(self.ba2) < Utility.uniqueID(self.ba1):
self.ba1, self.ba2 = self.ba2, self.ba1
+ self.angle = angle
self.improper = (self.ba1.ca is self.ba2.ca)
if self.improper:
self.b1 = self.ba1.otherBond(cb)
@@ -156,6 +293,191 @@
return 'Dihedral(' + `self.a1` +','+ `self.a2` +','+ \
`self.a3` +','+ `self.a4` +')'
+ def otherBond(self, bond1, bond2):
+ blist = (self.b1, self.b2, self.b3)
+ ilist = [0, 0, 0]
+ for i in range(0, 3):
+ if bond1 == blist[i] or bond2 == blist[i]: ilist[i] = 1
+ if ilist.count(1) < 2:
+ if bond1 in blist:
+ msg = 'bond ' + repr(bond2) + ' not in dihedral' + repr(blist)
+ else:
+ msg = 'bond ' + repr(bond1) + ' not in dihedral' + repr(blist)
+ raise ValueError,msg
+ if ilist.count(1) > 2: # should be impossible
+ msg = 'DihedralAngle.otherBond too many matches ' + str(ilist)
+ raise ValueError,msg
+ return blist[ilist.index(0)]
+
+ def endsIn(self, a):
+ return a is self.a4 or a is self.a1
+
+ def hasAllBonds(self):
+ return 1
+
+ def calcAngle(self):
+ try:
+ v1 = self.b1.asVector(self.a2)
+ cv = self.b2.asVector(self.a2)
+ v3 = self.b3.asVector(self.a3)
+ except:
+ return None
+ return calcDihedralAngle(v1, v3, cv)
+
+ def setAngle(self):
+ self.angle = self.calcAngle()
+ return self.angle
+
+ def constrains(self, new_bond, require_positions = 1):
+ if self.angle is None:
+ return None
+ try:
+ if new_bond.hasAtom(self.a4) and new_bond.hasAtom(self.a3):
+ if not require_positions:
+ return self
+ if Utility.isDefinedPosition(self.a1.position())\
+ and Utility.isDefinedPosition(self.a2.position())\
+ and Utility.isDefinedPosition(self.a3.position()):
+ return self
+ elif new_bond.hasAtom(self.a1) and new_bond.hasAtom(self.a2):
+ if not require_positions:
+ return self
+ if Utility.isDefinedPosition(self.a4.position())\
+ and Utility.isDefinedPosition(self.a2.position())\
+ and Utility.isDefinedPosition(self.a3.position()):
+ return self
+ except AttributeError: pass
+ return None
+
+ def hasAtom(self, a):
+ return self.a1 is a or self.a2 is a or self.a3 is a or self.a4 is a
+
+ def hasBond(self, b):
+ return self.b1 is b or self.b2 is b or self.b3 is b
+
+ def fillAtom(self, atom):
+ if self.a1 is None:
+ self.a1 = atom
+ elif self.a4 is None:
+ self.a4 = atom
+ else:
+ raise ValueError,'Neither a1 nor a4 are empty'
+
+class DihedralConstraint(DihedralAngle):
+
+ def __init__(self, b1, b3, cb, angle):
+ if b1 is not None and cb is not None:
+ ba1 = BondAngle(b1, cb, findCommonAtom(b1, cb))
+ ba2 = BondAngle(cb, b3, findCommonAtom(cb, b3))
+ DihedralAngle.__init__(self, ba1, ba2, cb)
+ else: # b3 should always be a Bond, b1 and/or cb can be None here
+ self.b1 = None
+ self.b3 = b3
+ self.ba1 = None
+ self.b2 = cb
+ if cb is not None:
+ if cb.hasAtom(b3.a1):
+ ca = b3.a1
+ else:
+ ca = b3.a2
+ self.ba2 = BondAngle(cb, b3, ca)
+ self.a2 = cb.otherAtom(ca)
+ self.a3 = ca
+ self.a4 = b3.otherAtom(ca)
+ else:
+ self.ba2 = None
+ self.a2 = None
+ self.a3 = b3.a1
+ self.a4 = b3.a2
+ self.a1 = None
+ self.normalized = None
+ self.improper = 1
+ self.angle = angle
+
+ def hasAllBonds(self):
+ return self.b1 is not None and self.b2 is not None and self.b3 is not None
+
+ def __repr__(self):
+ try:
+ constraint_str = '%f*deg' % self.angle/deg
+ except:
+ constraint_str = `self.angle`
+ return 'DihedralConstraint(' + `self.a1` +','+ `self.a2` +','+ \
+ `self.a3` +','+ `self.a4` +','+ constraint_str +')'
+ __str__ = __repr__
+
+def fetchBond(bonds_to_atom, a1, a2):
+ for b in bonds_to_atom[a1]:
+ if b.hasAtom(a2):
+ return b
+ return None
+
+#
+# Dihedral is mainly just different way to construct a DihedralAngle.
+# bonds_to_atom is a dict produced by ChemicalObjects.constructBondsToAtom,
+# needed if you want the bonds in the Dihedral to refer to pre-existing
+# bond instances.
+#
+
+class Dihedral(DihedralAngle):
+ def __init__(self, blueprint, memo = None, bonds_to_atom = None):
+ if blueprint is not None:
+ if type(blueprint) is type(()) or type(blueprint) is type([]):
+ if len(blueprint) == 4 and hasattr(blueprint[0],'a1'):
+ cb = blueprint[2]
+ ca1 = findCommonAtom(blueprint[0], cb)
+ ba1 = BondAngle(blueprint[0], cb, ca1)
+ ca2 = findCommonAtom(blueprint[1], cb)
+ ba2 = BondAngle(blueprint[1], cb, ca2)
+ angle = blueprint[3]
+ else:
+ a1 = blueprint[0]
+ a2 = blueprint[1]
+ a3 = blueprint[2]
+ a4 = blueprint[3]
+ angle = blueprint[4]
+ if bonds_to_atom is None:
+ cb = Bond((a2, a3))
+ ba1 = BondAngle(Bond((a1, a2)), cb, a2)
+ ba2 = BondAngle(cb, Bond((a3, a4)), a3)
+ else:
+ cb = fetchBond(bonds_to_atom, a2, a3)
+ ba1 = BondAngle(fetchBond(bonds_to_atom, a1, a2), cb, a2)
+ ba2 = BondAngle(cb, fetchBond(bonds_to_atom, a3, a4), a3)
+ else:
+ a2 = Database.instantiate(blueprint.a1, memo)
+ a1 = Database.instantiate(blueprint.a2, memo)
+ a3 = Database.instantiate(blueprint.a3, memo)
+ a4 = Database.instantiate(blueprint.a4, memo)
+ angle = Database.instantiate(blueprint.angle, memo)
+ if bonds_to_atom is None:
+ cb = Bond((a1, a3))
+ ba1 = BondAngle(Bond((a1, a2)), cb, a1)
+ ba2 = BondAngle(cb, Bond((a3, a4)), a3)
+ else:
+ cb = fetchBond(bonds_to_atom, a1, a3)
+ ba1 = BondAngle(fetchBond(bonds_to_atom, a1, a2), cb, a1)
+ ba2 = BondAngle(cb, fetchBond(bonds_to_atom, a3, a4), a3)
+ DihedralAngle.__init__(self, ba1, ba2, cb, angle)
+
+ def __repr__(self):
+ try:
+ angle = `self.angle/deg` + '*deg'
+ except:
+ angle = 'None'
+ if self.improper:
+ return 'ImproperDihedral(' + `self.a1` +','+ `self.a2` +','+ \
+ `self.a3` +','+ `self.a4` + ',' + angle + ')'
+ else:
+ return 'Dihedral(' + `self.a1` +','+ `self.a2` +','+ \
+ `self.a3` +','+ `self.a4` + ',' + angle + ')'
+ __str__ = __repr__
+
+
+ blueprintclass = Database.BlueprintDihedral
+
+Database.registerInstanceClass(Dihedral.blueprintclass, Dihedral)
+
#
# Bond lists
#
@@ -166,9 +488,8 @@
class BondList(UserList):
def __init__(self, list):
- if list is not None:
- UserList.__init__(self, list)
- self._clearCache()
+ UserList.__init__(self, list)
+ self._clearCache()
def __getinitargs__(self):
return (None,)
def __getstate__(self):
@@ -304,12 +625,73 @@
def __init__(self, dihedrals):
self.data = dihedrals
+ # cache entries by bond for faster searches
+ def generateMap(self):
+ bond_map = {}
+ for d in self.data:
+ for b in (d.b1, d.b2, d.b3):
+ if not bond_map.has_key(b):
+ bond_map[b] = [d]
+ else:
+ bond_map[b].append(d)
+ self.bond_map = bond_map
+
def __repr__(self):
return repr(self.data)
def __len__(self):
return len(self.data)
def __getitem__(self, i):
return self.data[i]
+ def __setitem__(self, i, val):
+ self.data[i] = val
+ def __getslice__(self, first, last):
+ return DihedralAngleList(self.data[first:last])
+
+ def searchList(self, search_bond):
+ if search_bond is None:
+ return []
+ try:
+ try:
+ return self.bond_map[search_bond]
+ except AttributeError:
+ return self.data
+ except KeyError:
+ return []
+
+ def findConstraints(self, bond, require_positions = 1):
+ result = []
+ for da in self.searchList(bond):
+ if da.constrains(bond, require_positions):
+ result.append(da)
+ return result
+
+ def findRotatesAbout(self, cmbond, bonds_to_atoms = None):
+ for da in self.data:
+ if not hasattr(da, 'b2'):
+ Utility.warning('findRotatesAbout bad dihedral ' + str(da))
+ continue
+ if da.b2 is cmbond:
+ try:
+ if bonds_to_atoms is None or \
+ ((da.b1 in bonds_to_atoms[da.b1.a1]) \
+ and (da.b2 in bonds_to_atoms[da.b2.a1]) \
+ and (da.b3 in bonds_to_atoms[da.b3.a1])):
+ return da
+ except KeyError: pass
+ return None
+
+ def __add__(self, other):
+ return DihedralAngleList(self.data+other.data)
+ __radd__ = __add__
+
+ def index(self, object):
+ return self.data.index(object)
+
+ def remove(self, object):
+ self.data.remove(object)
+
+ def append(self, object):
+ self.data.append(object)
#
# Dummy bond length database, for constraints without a force field
--- old/ChemicalObjects.py Tue Feb 29 05:29:56 2000
+++ ChemicalObjects.py Tue Dec 26 20:27:06 2000
@@ -191,6 +191,18 @@
"Returns a list containing all atoms in the object."
return self.atoms
+ def bondList(self):
+ try:
+ return self.bonds
+ except AttributeError:
+ return []
+
+ def danglingBondList(self):
+ try:
+ return self.danglingbonds
+ except AttributeError:
+ return []
+
def setPosition(self, atom, position):
if atom.__class__ is Database.AtomReference:
atom = self.atoms[atom.number]
@@ -258,6 +270,16 @@
for o in self._subunits():
n = n + len(o._distanceConstraintList())
return n
+
+ def constructBondsToAtom(self):
+ bonds_to_atom = {}
+ for a in self.atoms:
+ bonds_to_atom[a] = []
+ for b in self.bonds:
+ bonds_to_atom[b.a1].append(b)
+ bonds_to_atom[b.a2].append(b)
+ self.bonds_to_atom = bonds_to_atom
+ return bonds_to_atom
def setBondConstraints(self, universe=None):
if universe is None:
--- old/Database.py Tue Feb 29 05:29:56 2000
+++ Database.py Tue Dec 26 20:27:38 2000
@@ -152,7 +152,8 @@
import GroupEnvironment
ChemicalObjectType.__init__(self, filename, GroupEnvironment,
('atoms', 'groups', 'bonds',
- 'chain_links'))
+ 'chain_links', 'dihedrals',
+ 'danglingbonds'))
for g in self.groups:
self.atoms = self.atoms + g.atoms
self.bonds = self.bonds + g.bonds
@@ -387,17 +388,65 @@
#
class BlueprintBond:
- def __init__(self, a1, a2):
+ def __init__(self, a1, a2, order = 0):
self.a1 = a1
self.a2 = a2
+ self.order = order
is_instance_var = 1
def _blueprintCopy(self, memo):
return BlueprintBond(_blueprintCopy(self.a1, memo),
- _blueprintCopy(self.a2, memo))
+ _blueprintCopy(self.a2, memo), self.order)
object_list = 'bonds'
+
+ def __copy__(self, memo = None):
+ return self
+ __deepcopy__ = __copy__
+
+class BlueprintDanglingBond:
+
+ def __init__(self, a1, constraint, direction = None):
+ self.a1 = a1
+ self.constraint = constraint
+ self.direction = direction
+
+ is_instance_var = 1
+
+ def _blueprintCopy(self, memo):
+ return BlueprintDanglingBond(_blueprintCopy(self.a1, memo),
+ _blueprintCopy(self.constraint, memo),
+ _blueprintCopy(self.direction, memo))
+
+ object_list = 'danglingbonds'
+
+ def __copy__(self, memo = None):
+ return self
+ __deepcopy__ = __copy__
+
+#
+# The dihedral class just keeps track of the 3 bonds involved.
+#
+class BlueprintDihedral:
+
+ def __init__(self, a1, a2, a3, a4, angle):
+ self.a1 = a1
+ self.a2 = a2
+ self.a3 = a3
+ self.a4 = a4
+ self.angle = angle
+
+ is_instance_var = 1
+
+ def _blueprintCopy(self, memo):
+ return BlueprintDihedral(_blueprintCopy(self.a1, memo),
+ _blueprintCopy(self.a2, memo),
+ _blueprintCopy(self.a3, memo),
+ _blueprintCopy(self.a4, memo),
+ _blueprintCopy(self.angle, memo))
+
+ object_list = 'dihedrals'
def __copy__(self, memo = None):
return self
--- old/Collection.py Tue Feb 29 05:29:56 2000
+++ Collection.py Sun Dec 17 18:47:09 2000
@@ -50,6 +50,9 @@
except AttributeError: pass
return n
+ def bondList(self): # temp hack
+ return self.bonds
+
def degreesOfFreedom(self):
"Returns the number of mechanical degrees of freedom."
return 3*(self.numberOfAtoms()-self.numberOfFixedAtoms())
@@ -263,6 +266,12 @@
Transformation.Translation(-cm)
self.applyTransformation(t)
+ def rotateAroundPoint(self, axis, angle, point):
+ t = Transformation.Translation(point) * \
+ Transformation.Rotation(axis, angle) * \
+ Transformation.Translation(-point)
+ self.applyTransformation(t)
+
def rotateAroundOrigin(self, axis, angle):
"""Rotates the object by the given |angle| around an axis
that passes through the coordinate origin and has the given
@@ -512,6 +521,20 @@
def atomList(self):
"Returns a list containing all atoms of all objects in the collection."
lists = map(lambda o: o.atomList(), self.objectList())
+ if lists:
+ return reduce(operator.add, lists)
+ else:
+ return []
+
+ def bondList(self):
+ lists = map(lambda o: o.bondList(), self.objectList())
+ if lists:
+ return reduce(operator.add, lists)
+ else:
+ return []
+
+ def danglingBondList(self):
+ lists = map(lambda o: o.danglingBondList(), self.objectList())
if lists:
return reduce(operator.add, lists)
else:
--- old/Universe.py Tue Feb 29 05:30:00 2000
+++ Universe.py Sun Dec 17 18:44:31 2000
@@ -136,6 +136,12 @@
self._atoms = self._objects.atomList()
return self._atoms
+ def bondList(self):
+ return self._objects.bondList()
+
+ def danglingBondList(self):
+ return self._objects.danglingBondList()
+
def universe(self):
"Returns the universe itself."
return self
--- old/Biopolymers.py Tue Feb 29 05:29:56 2000
+++ Biopolymers.py Tue Dec 26 20:38:54 2000
@@ -55,9 +55,15 @@
def _setupChain(self, circular, properties, conf):
self.atoms = []
self.bonds = []
+ self.danglingbonds = []
+ self.dihedrals = []
for g in self.groups:
self.atoms = self.atoms + g.atoms
self.bonds = self.bonds + g.bonds
+ if hasattr(g, 'danglingbonds') and g.danglingbonds:
+ self.danglingbonds = self.danglingbonds + g.danglingbonds
+ if hasattr(g, 'dihedrals') and g.dihedrals:
+ self.dihedrals = self.dihedrals + g.dihedrals
for i in range(len(self.groups)-1):
link1 = self.groups[i].chain_links[1]
link2 = self.groups[i+1].chain_links[0]
@@ -67,6 +73,8 @@
link2 = self.groups[0].chain_links[0]
self.bonds.append(Bonds.Bond((link1, link2)))
self.bonds = Bonds.BondList(self.bonds)
+ self.danglingbonds = Bonds.BondList(self.danglingbonds)
+ self.dihedrals = Bonds.DihedralAngleList(self.dihedrals)
self.parent = None
self.type = None
self.configurations = {}
--- old/GroupEnvironment.py Tue Feb 29 05:29:58 2000
+++ GroupEnvironment.py Tue Dec 26 20:27:50 2000
@@ -3,8 +3,11 @@
_undocumented = 1
-from Database import BlueprintAtom, BlueprintGroup, BlueprintBond
+from Database import BlueprintAtom, BlueprintGroup, BlueprintBond, \
+ BlueprintDihedral, BlueprintDanglingBond
from Units import *
Atom = BlueprintAtom
Group = BlueprintGroup
Bond = BlueprintBond
+DanglingBond = BlueprintDanglingBond
+Dihedral = BlueprintDihedral
--- old/ForceFields/MMForceField.py Tue Feb 29 05:31:00 2000
+++ ForceFields/MMForceField.py Tue Dec 26 12:20:29 2000
@@ -32,6 +32,7 @@
ESMPForceField, ESEwaldForceField
from ForceField import CompoundForceField, ForceFieldData
from MMTK import ParticleProperties
+from MMTK.Units import deg
from Scientific.Geometry import Vector
import copy
@@ -58,6 +59,58 @@
.charge_property)
global_data.charge = charge
+ def defaultBondLength(self, bond, object):
+ try:
+ type1 = object.getAtomProperty(bond.a1, self.dataset.atom_type_property)
+ except (KeyError, AttributeError):
+ #print bond.a1, object, self.__dict__.keys()
+ raise TypeError,'atom ' + bond.a1.fullName() + ' has no type specified'
+ try:
+ type2 = object.getAtomProperty(bond.a2, self.dataset.atom_type_property)
+ except (KeyError, AttributeError):
+ raise TypeError,'atom ' + bond.a2.fullName() + ' has no type specified'
+ (l,k) = self.dataset.bondParameters(type1, type2)
+ return l
+
+ def defaultBondAngle(self, bond1, bond2, object):
+ try:
+ type1 = object.getAtomProperty(bond1.a1, self.dataset.atom_type_property)
+ except:
+ raise TypeError,'atom ' + bond1.a1.name + ' has no type specified'
+ try:
+ type2 = object.getAtomProperty(bond1.a2, self.dataset.atom_type_property)
+ except:
+ if bond1.a2 is None:
+ type2 = None
+ else:
+ raise TypeError,'atom ' + bond1.a2.name + ' has no type specified'
+ if not bond2.hasAtom(bond1.a2):
+ (type1, type2) = (type2, type1)
+ if bond1.hasAtom(bond2.a1):
+ third_atom = bond2.a2
+ else:
+ third_atom = bond2.a1
+ try:
+ type3 = object.getAtomProperty(third_atom, self.dataset.atom_type_property)
+ except:
+ if third_atom is None:
+ (a,energy) = self.dataset.bondAngleAvgParameters(type1, type2)
+ return a
+ else:
+ raise TypeError,'atom ' + third_atom.name + ' has no type specified'
+ try:
+ (a,energy) = self.dataset.bondAngleParameters(type1, type2, type3)
+ except:
+ a = 109.5*deg
+ return a
+
+ def valence(self, atom, object):
+ try:
+ type1 = object.getAtomProperty(atom, self.dataset.atom_type_property)
+ except:
+ raise TypeError,'atom ' + `atom` + ' has no type specified'
+ return self.dataset.valence(type1)
+
def _atomType(self, object, atom, global_data):
return global_data.atom_type[atom]
@@ -113,7 +166,8 @@
p = self.dataset.bondAngleParameters(t1, tc, t2)
except KeyError:
raise KeyError, 'No parameters for angle ' + `a1` + \
- '--' + `ca` + '--' + `a2`
+ '--' + `ca` + '--' + `a2` + ', types ' + `t1` + ',' +\
+ `tc` + ',' + `t2`
if p is not None and p[1] != 0.:
data.add('angles', (i1, ic, i2) + p)
@@ -161,6 +215,11 @@
def bondLengthDatabase(self, universe):
return MMBondLengthDatabase(self, universe, self.dataset)
+
+ def defaultBondLength(self,bond,object):
+ return AmberForceFieldTerm.defaultBondLength(self,bond,object)
+ def defaultBondAngle(self,bond1,bond2,object):
+ return AmberForceFieldTerm.defaultBondAngle(self,bond1,bond2,object)
def bonds(self, global_data):
return self.data.get('bonds')
--- old/ForceFields/Amber/AmberData.py Tue Feb 29 05:31:41 2000
+++ ForceFields/Amber/AmberData.py Tue Dec 26 12:21:12 2000
@@ -17,6 +17,11 @@
amber_angle_unit = Units.deg
amber_fc_angle_unit = Units.rad
+valence = { 'HC' : 1, 'HA' : 1, 'HA' : 1, 'H1' : 1, 'H2' : 1, 'H3' : 1, \
+ 'H4' : 1, 'SI' : 4, 'CT' : 4, 'CA' : 3, 'CM' : 3, 'C' : 3, \
+ 'O' : 2, 'OS' : 2, 'O2' : 1, 'OH' : 2, \
+ 'NC' : 3, 'N2' : 3, 'N*' : 3, 'P' : 4 }
+
#
# The force field parameter file
#
@@ -46,6 +51,8 @@
self._readBondParameters(file)
self.bond_angles = {}
+ self.bond_angles_average = {}
+ bond_angle_sum = {}
self._readAngleParameters(file)
self.dihedrals = {}
@@ -126,6 +133,9 @@
l[2], l[3])
def _readAngleParameters(self, file):
+ self.bond_angles_average = {}
+ bond_angle_sum = {}
+
format = FortranFormat('A2,1X,A2,1X,A2,2F10.2')
while 1:
l = FortranLine(file.readline()[:-1], format)
@@ -140,6 +150,31 @@
self.atom_types[name3],
l[3], l[4])
+ try:
+ (asum,ksum,count) = bond_angle_sum[(name1, name2)]
+ except:
+ asum = 0
+ ksum = 0
+ count = 0
+ bond_angle_sum[(name1, name2)] = (asum + l[4],
+ ksum + l[3], count + 1)
+ try:
+ (asum,ksum,count) = bond_angle_sum[name2]
+ except:
+ asum = 0
+ ksum = 0
+ count = 0
+ bond_angle_sum[name2] = (asum + l[4], ksum + l[3], count + 1)
+ for key in bond_angle_sum.keys():
+ (angle_sum, k_sum, n) = bond_angle_sum[key]
+ try:
+ (k1, k2) = key
+ except:
+ k1 = None
+ k2 = key
+ self.bond_angles_average[key] =\
+ AmberBondAngleParameters(k1, k2, None, k_sum/n, angle_sum/n)
+
def _readDihedralParameters(self, file):
format = FortranFormat('A2,1X,A2,1X,A2,1X,A2,I4,3F15.2')
append = None
@@ -239,6 +274,16 @@
return (p.angle*amber_angle_unit,
p.k*amber_energy_unit/amber_fc_angle_unit)
+ def bondAngleAvgParameters(self, name1, name2):
+ name1 = _normalizeName(name1)
+ if name2 is None:
+ p = self.bond_angles_average[name1]
+ else:
+ name2 = _normalizeName(name2)
+ p = self.bond_angles_average[(name1, name2)]
+ return (p.angle*amber_angle_unit,
+ p.k*amber_energy_unit/amber_fc_angle_unit)
+
def dihedralParameters(self, name1, name2, name3, name4):
name1 = _normalizeName(name1)
name2 = _normalizeName(name2)
@@ -311,6 +356,10 @@
return (eps, sigma, 1)
else:
raise ValueError, 'Type ' + ljp.type + ' not supported.'
+
+ def valence(self, name):
+ name = _normalizeName(name)
+ return valence[name]
#
# Atom type class
--
------------------------------------------------------------------------------
Peter McCluskey | Fed up with democracy's problems? Examine Futarchy:
http://www.rahul.net/pcm | http://hanson.gmu.edu/futarchy.pdf or .ps