[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