# -*- mode:python;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
import os
import sys
import inspect
newpath = os.path.dirname(inspect.getfile(
inspect.currentframe())) # script directory
if not newpath in sys.path:
sys.path.append(newpath)
import string
import re
import gzip
import itertools
from sym_util import *
[docs]def homogenizechains(sel='all'):
chains = cmd.get_chains()
# cmd.remove("hydro")
# cmd.remove("resn HOH")
# cmd.remove("(HET and not resn MSE+CSW)")
for c1, c2 in ((x, y) for x in chains for y in chains if x < y):
sel1 = "(%s) and chain %s and name ca" % (sel, c1)
sel2 = "(%s) and chain %s and name ca" % (sel, c2)
#
a = cmd.get_model("%s" % (sel1))
b = cmd.get_model("%s" % (sel2))
sa = "".join([getname1(x.resn) for x in a.atom])
sb = "".join([getname1(x.resn) for x in b.atom])
if sa == sb:
continue
ra = [myint(x.resi) for x in a.atom]
rb = [myint(x.resi) for x in b.atom]
# if max(ra) - min(ra) + 1 != len(ra): print "missing residue numbers",max(ra),min(ra),len(ra)
# if max(rb) - min(rb) + 1 != len(rb): print "missing residue
# numbers",rb
mla, mua, mlb, mub = lcs(sa, sb)
bla, bua, blb, bub = lcs(sa[:mla], sb[:mlb])
ala, aua, alb, aub = lcs(sa[mua + 1:], sb[mub + 1:])
ra = ra[mla:(mua + 1)]
rb = rb[mlb:(mub + 1)]
if len(ra[bla:(bua + 1)]) > 10:
ra = ra[bla:(bua + 1)] + ra[mla:(mua + 1)] + ra[ala:(aua + 1)]
rb = rb[blb:(bub + 1)] + rb[mlb:(mub + 1)] + rb[alb:(aub + 1)]
if len(ra[ala:(aua + 1)]) > 10:
ra += ra[ala:(aua + 1)]
rb += rb[alb:(aub + 1)]
for i in getres("%s" % (sel1), False):
if not i in ra:
print "remove resi", i, "from '" + sel1 + "'"
cmd.remove("%s and resi %i" % (sel1, i))
for i in getres("%s" % (sel2), False):
if not i in rb:
print "remove resi", i, "from '" + sel2 + "'"
cmd.remove("%s and resi %i" % (sel2, i))
return False
# def pickandfixchains(N=None,sel="all"):
# # find chains
# # homogenize all pairs until fixed
# if not N: N = cmd.getnum_chains(sel)
# cc = []
# for c in getchain(sel):
# cc.append((-cmd.select("%s and chain %s and name CA"%(sel,c)),c))
# cc.sort()
# chains = [x[1] for x in cc[:N]]
# done = False
# count = 0
# while not done:
# if count > 10: break
# count += 1
# done = True;
# random.shuffle(chains)
# for i in range(1,len(chains)):
# done = done and homogenizechains(sel,chains[0],chains[i])
# print chains
# if N is 2: alignc2(sel,"name ca",chains=chains)
# if N is 3: alignc3(sel,"name ca",chains=chains)
# if N is 4: alignc4(sel,"name ca",chains=chains)
# if N is 5: alignc5(sel,"name ca",chains=chains)
# chains.sort()
# return chains[0]
[docs]def processhomomers():
o = open("log", 'w')
for n in (2, 3, 4, 5):
files = glob.glob("c%ipdb/*.pdb.gz" % n)
random.shuffle(files)
for f in files:
o.write(f + "\n")
o.flush()
cmd.delete("all")
try:
cmd.load(f)
c = pickandfixchains(n)
cmd.save("c%ia/" % n + f[3:-3], "chain %s" % c)
except:
print "fail on", f
o.close()
[docs]def iscontig(sel):
m = cmd.get_model(sel + " and name N+CA+C").atom
for i in range(1, len(m)):
if (xyz.Vec(m[i - 1].coord) - xyz.Vec(m[i].coord)).length() > 1.8:
return False
return True
[docs]def procCdat(N=3, lfile=None, biod="/data/biounit", outd=None):
if lfile is None:
lfile = os.path.expanduser(
"~/Dropbox/project/sym_comp/meta/C%i.dat" % N)
if outd is None:
outd = os.path.expanduser(".//C%i" % N)
print outd
Nnobio = 0
Nok = 0
Nbig = 0
Nnsym = 0
Nnomxatm = 0
Nhomogen = 0
for fn in open(lfile).readlines():
try:
fn = fn.split()[0].strip()
pdb = fn[3:7]
bnum = int(fn[-1:])
if os.path.exists(outd + "/" + pdb + "_" + str(bnum) + "_sub1.pdb") or os.path.exists(outd + "/" + pdb + "_" + str(bnum) + "_sub1.pdb.gz"):
Nok += 1
continue
fname = biod + "/" + fn
if not os.path.exists(fname):
fname += ".gz"
if not os.path.exists(fname):
Nnobio += 1
print "no file", fname
continue
cmd.delete("all")
cmd.load(fname, 'm')
cmd.remove("resn HOH")
cmd.remove('not alt a+""')
#hf = cmd.select("(HET and not resn MSE+CSW)",state=1) / cmd.select("ALL",state=1)
# if hf > 0.1: continue
#cmd.remove("(HET and not resn MSE+CSW)")
if cmd.select('all', state=N) != 0:
for i in range(1, N + 1):
cmd.create("sub%i" % i, "m", i, 1)
else:
cc = chaincount("m")
if len(cc) < N:
Nnsym += 1
print "ERROR:", fname, " symmetry error"
continue
for i in range(1, N + 1):
cmd.create("sub%i" % i, "m and chain %s" %
(cc[-i][1]), 1, 1)
for i in range(1, N + 1):
if iscontig("sub%i" % i):
cmd.create("mxatm", "sub%i" % i)
break
if cmd.select("mxatm") < 50:
Nnomxatm += 1
print "ERROR: mxatm < 50"
continue
if cmd.select("name CA and mxatm") > 500:
Nbig += 1
print "ERROR:", fname, " more than 500 res"
continue
chains = ["sub%i" % i for i in range(1, N + 1)]
done = False
count = 0
while not done and count < 50:
done = True
random.shuffle(chains)
for i in range(len(chains)):
for j in range(i + 1, len(chains)):
done = done and homogenizechains(chains[i], chains[j])
count += 1
if count >= 50:
Nhomogen += 1
print "ERROR: error homogenizing"
continue
if cmd.select("sub1") < 50:
Nnomxatm += 1
print "ERROR: less than 50 atoms"
continue
cm = com("sub*")
for i in range(1, N + 1):
trans("sub%i" % i, -cm)
a = [cmd.get_model("sub%i and name CA" %
i).atom for i in range(1, N + 1)]
axis = xyz.Vec(0, 0, 0)
for i in range(len(a[0])):
axis1 = xyz.Vec(0, 0, 0)
for j in range(N):
axis1 += xyz.Vec(a[j][i].coord)
if axis1.length() > 0.0001 and axis.dot(axis1) < 0:
axis1 *= -1
axis += axis1
axis.normalize()
for i in range(1, N + 1):
alignaxis("sub%i" % i, xyz.Vec(
0, 0, 1), axis, xyz.Vec(0, 0, 0))
# cmd.create("final1","mxatm")
# cmd.create("final2","mxatm")
# cmd.create("final3","mxatm")
# cmd.align("final1","sub1")
# cmd.align("final2","sub2")
# cmd.align("final3","sub3")
# return
if not os.path.exists(outd):
os.mkdir(outd)
cmd.align("mxatm", "sub1")
cmd.save(outd + "/" + pdb + "_" + str(bnum) + "_sub1.pdb", "mxatm")
Nok += 1
print "SUCCESS on", N, fname
except (pymol.CmdException, Exception) as e:
print "EXCEPTION!!!!!!", e
print fname
pass
print "DONE"
print Nok, Nbig, Nnsym, Nnobio, Nnomxatm, Nhomogen
[docs]def procD2dat(lfile=None, biod="/data/biounit", outd=None):
N = 4
if lfile is None:
lfile = os.path.expanduser('~/Dropbox/project/sym_comp/meta/D2.dat')
if outd is None:
outd = os.path.expanduser("./D2")
print outd
Nnobio = 0
Nok = 0
Ncontact = 0
Nbig = 0
Nnsym = 0
Nnomxatm = 0
Nhomogen = 0
for fn in open(lfile).readlines():
try:
fn = fn.strip()
fn = fn.split()[0]
pdb = fn[-9:-5]
print fn, pdb
bnum = int(fn[-1:])
if os.path.exists(outd + "/" + pdb + "_" + str(bnum) + "_sub1.pdb") or os.path.exists(outd + "/" + pdb + "_" + str(bnum) + "_sub1.pdb.gz"):
Nok += 1
continue
fname = biod + "/" + fn
if not os.path.exists(fname):
fname += ".gz"
if not os.path.exists(fname):
Nnobio += 1
print "can't find", fname
continue
# print pdb,bnum,fname
cmd.delete("all")
cmd.load(fname, 'm')
cmd.remove("resn HOH")
cmd.remove('not alt a+""')
#hf = cmd.select("(HET and not resn MSE+CSW)",state=1) / cmd.select("ALL",state=1)
# if hf > 0.1: continue
#cmd.remove("(HET and not resn MSE+CSW)")
if cmd.select('all', state=4) != 0:
for i in range(1, N + 1):
cmd.create("sub%i" % i, "m", i, 1)
elif cmd.select('all', state=2) != 0:
cc = chaincount('m')
if len(cc) < 2:
Nnsym += 1
continue
cmd.create("sub1", "m and chain %s" % (cc[0][1]), 1, 1)
cmd.create("sub2", "m and chain %s" % (cc[1][1]), 1, 1)
cmd.create("sub3", "m and chain %s" % (cc[0][1]), 2, 1)
cmd.create("sub4", "m and chain %s" % (cc[1][1]), 2, 1)
else:
cc = chaincount("m")
if len(cc) < N:
sym = cmd.get_symmetry("m")
if sym[6] == "I 2 2 2":
trans('m', xyz.Vec(0, -sym[1], 0))
print pdb
# elif sym[6] == "P 21 21 2" and len(cc)==2:
# trans('m',xyz.Vec(-sym[0]/2.0,-sym[1]/2.0,0))
# cmd.create("sub1","m and chain %s and not (HET and not resn MSE+CSW)"%(cc[0][1]),1,1)
# cmd.create("sub2","m and chain %s and not (HET and not resn MSE+CSW)"%(cc[1][1]),1,1)
# cmd.create("sub3","m and chain %s and not (HET and not resn MSE+CSW)"%(cc[0][1]),1,1)
# cmd.create("sub4","m and chain %s and not (HET and not resn MSE+CSW)"%(cc[1][1]),1,1)
# rot("sub3",xyz.Vec(0,0,1),180,xyz.Vec(0,0,0))
# rot("sub4",xyz.Vec(0,0,1),180,xyz.Vec(0,0,0))
elif sym[6] in ('C 1 2 1', 'P 21 21 21', 'P 62 2 2', 'P 64 2 2', 'P 65 2 2', 'P 63 2 2', 'P 61 2 2', 'C 2 2 21'):
Nnsym += 1
# if pid != "1y2k_2": return
continue
else:
Nnsym += 1
continue # return
cmd.save(outd + "/" + pdb + "_" +
str(bnum) + "_sub1.pdb", "m")
Nok += 1
continue
else:
for i in range(1, N + 1):
cmd.create("sub%i" % i, "m and chain %s" %
(cc[-i][1]), 1, 1)
for i in range(1, N + 1):
if iscontig("sub%i" % i):
cmd.create("mxatm", "sub%i" % i)
break
if cmd.select("name CA and mxatm") > 250:
Nbig += 1
continue
if cmd.select("mxatm") < 50:
Nnomxatm += 1
continue
chains = ["sub%i" % i for i in range(1, N + 1)]
done = False
count = 0
while not done and count < 50:
done = True
random.shuffle(chains)
for i in range(len(chains)):
for j in range(i + 1, len(chains)):
done = done and homogenizechains(chains[i], chains[j])
count += 1
if count >= 50:
Nhomogen += 1
continue
if cmd.select("sub1") < 50:
Nnomxatm += 1
continue
cm = com("sub*")
for i in range(1, N + 1):
trans("sub%i" % i, -cm)
a = [cmd.get_model("sub%i and name CA" %
i).atom for i in range(1, N + 1)]
a1 = xyz.Vec(0, 0, 0)
for i in range(len(a[0])):
axis1 = xyz.Vec(a[0][i].coord) + xyz.Vec(a[1][i].coord)
if axis1.length() > 0.0001 and a1.dot(axis1) < 0:
axis1 *= -1
a1 += axis1
a1.normalize()
for i in range(1, N + 1):
alignaxis("sub%i" % i, xyz.Vec(1, 0, 0), a1, xyz.Vec(0, 0, 0))
a = [cmd.get_model("sub%i and name CA" %
i).atom for i in range(1, N + 1)]
a1 = xyz.Vec(0, 0, 0)
for i in range(len(a[0])):
axis1 = xyz.Vec(a[0][i].coord) + xyz.Vec(a[2][i].coord)
if axis1.length() > 0.0001 and a1.dot(axis1) < 0:
axis1 *= -1
a1 += axis1
a1.normalize()
for i in range(1, N + 1):
alignaxis("sub%i" % i, xyz.Vec(0, 1, 0), a1, xyz.Vec(0, 0, 0))
cmd.align("mxatm", "sub1")
cmd.create("final2", "mxatm")
cmd.create("final3", "mxatm")
cmd.create("final4", "mxatm")
rot('final2', xyz.Vec(1, 0, 0), 180, xyz.Vec(0, 0, 0))
rot('final3', xyz.Vec(0, 1, 0), 180, xyz.Vec(0, 0, 0))
rot('final4', xyz.Vec(0, 0, 1), 180, xyz.Vec(0, 0, 0))
n1 = cmd.select('mxatm within 4 of final2')
n2 = cmd.select('mxatm within 4 of final3')
n3 = cmd.select('mxatm within 4 of final4')
if n1 < 10 and n2 < 10 and n3 < 10:
Ncontact += 1
continue
if not os.path.exists(outd):
os.mkdir(outd)
cmd.save(outd + "/" + pdb + "_" + str(bnum) + "_sub1.pdb", "mxatm")
Nok += 1
# return
except:
print "exception!"
print Nok, Nbig, Nnsym, Ncontact, Nnobio, Nnomxatm, Nhomogen
[docs]def prepare_c2_nmr(pattern, outdir=None):
if None is outdir:
outdir = os.path.dirname(pattern)
if not os.path.exists(outdir):
os.mkdir(outdir)
for i in glob.glob(pattern):
try:
name = os.path.basename(i).split(".")[0]
cmd.delete('all')
cmd.load(i)
cmd.remove("hydro")
cmd.remove("resn hoh")
cmd.remove("not chain A+B")
Nres = cmd.select(name + " and chain A and name CA")
cmd.split_states(name)
cmd.delete(name)
for obj in cmd.get_object_list("all"):
s1, s2 = ("{0} and chain {1}".format(obj, c) for c in "AB")
n1, n2 = (cmd.select(x) for x in (s1, s2))
# homogenizechains(s1,s2)
if cmd.select(s1) < 0.8 * min(n1, n2):
print "ERROR", obj, n1, n2, cmd.select(s1)
raise Exception()
if not 0 == alignc2(obj, tgtaxis=xyz.Vec(0.816496579408716, 0, 0.57735027133783)):
print "ERROR"
raise Exception()
if not os.path.exists(outdir + "/" + obj[:-5]):
os.mkdir(outdir + "/" + obj[:-5])
cmd.save(outdir + "/" + obj[:-5] + "/" + obj + "A.pdb", s1)
cmd.save(outdir + "/" + obj[:-5] + "/" + obj + "B.pdb", s2)
print "success", i, Nres
except:
print "caught exception"
[docs]def make_cryst1_P432(fn):
print "makecryst1_P432", fn
cmd.load(fn, 'work_prot')
sele = 'work_prot'
fn = os.path.basename(fn)
if fn.endswith(".gz"):
fn = fn[:-3]
if fn.endswith(".pdb"):
fn = fn[:-4]
fn += "_cryst1.pdb"
A2 = Vec(-1.000, 0.000, 1.000)
I2 = Vec(0.750, 0.875, 0.000)
A3 = Vec(1.000, 1.000, -1.000)
I3 = Vec(0.000, 0.500, 1.000)
Dsel = sele + ' and chain A+D'
Tsel = sele + ' and chain A+B+C'
oldc2 = com(Dsel)
oldc3 = com(Tsel)
olda2 = c2axis(Dsel, chains=('A', 'D'))
olda3 = c3axis(Tsel)
print olda2.angle_degrees(olda3)
print A2.angle_degrees(A3)
assert abs(olda2.angle_degrees(olda3) - A2.angle_degrees(A3)) < 0.01
Xaln = xyz.alignvectors(olda2, olda3, A2, A3)
Xaln = Xaln # center trimer, then align axes, then move so a2 intersects x==0
xform(sele, Xaln)
trans(sele, -com(Tsel))
newc2 = com(Dsel)
newa2 = c2axis(Dsel, chains=('A', 'D'))
newa3 = c3axis(Tsel)
assert newa2.angle_degrees(A2) < 0.0001
assert newa3.angle_degrees(A3) < 0.0001
# cmd.hide('ev')
# x1 = xyz.rotation_around_degrees(A2,180,newc2)
# x2 = xyz.rotation_around_degrees(A3,120, V0 )
# x3 = xyz.rotation_around_degrees(A3,240, V0 )
# axes = set()
# cmd.create('work_minimal',sele+" and chain A and name N+CA+C+SG")
# for i,X in enumerate( xyz.expand_xforms((x1,x2,x3),9) ):
# cmd.create("sub%i"%i,"work_minimal")
# xform("sub%i"%i,X)
# axes.add( X.rotation_axis()[0].key() )
# print axes
# cmd.show('rib')
# cmd.show('sph','name SG')
# return
d = slide_to_make_lines_intersect(A3, A2, com(Dsel), I2, V0)
trans(sele, d * A3)
d = slide_to_make_lines_intersect(I2, A3, com(Tsel), I3, V0)
trans(sele, d * I2)
newc2 = com(Dsel)
print "c2", newc2
newc3 = com(Tsel)
assert A3.lineangle(newa3) < 0.0001
assert A2.lineangle(newa2) < 0.0001
assert line_line_distance(A2, newc2, I2, Vec(0, 0, 0)) < 0.0001
assert line_line_distance(A3, newc3, I3, Vec(0, 0, 0)) < 0.0001
cellsize = abs(newc2.x / I2.x)
print "\nAXIS ERROR is:", newa2.lineangle(A2) + newa3.lineangle(A3)
print "CELL SIZE", cellsize, '\n'
cmd.save(".tmp.pdb", sele + " and chain A")
with open(fn, 'w') as out:
out.write("CRYST1 %7.3f %7.3f %7.3f 90.00 90.00 90.00 I 21 3\n" % (
(cellsize,) * 3))
os.system("cat .tmp.pdb >> %s" % fn)
os.system("rm .tmp.pdb")
[docs]def make_cryst1_I213(fn):
print "makecryst1_i213", fn
cmd.load(fn, 'work_prot')
sele = 'work_prot'
fn = os.path.basename(fn)
if fn.endswith(".gz"):
fn = fn[:-3]
if fn.endswith(".pdb"):
fn = fn[:-4]
fn += "_cryst1.pdb"
Dsel = sele + ' and chain A+D'
Tsel = sele + ' and chain A+B+C'
oldc2 = com(Dsel)
oldc3 = com(Tsel)
olda2 = c2axis(Dsel, chains=('A', 'D'))
olda3 = c3axis(Tsel)
Xaln = xyz.alignvectors(olda2, olda3, Vec(0, 0, 1), Vec(1, 1, 1))
Xaln = Xaln # center trimer, then align axes, then move so a2 intersects x==0
xform(sele, Xaln)
trans(sele, -com(Tsel))
newc2 = com(Dsel)
trans(sele, xyz.Vec(-newc2.x))
print com(Dsel)
assert abs(com(Dsel).x) < xyz.SQRTEPS
newa2 = c2axis(Dsel, chains=('A', 'D'))
newa3 = c3axis(Tsel)
newc3 = com(Tsel)
assert xyz.Vec(1, 1, 1).lineangle(newc3) < 0.0001
assert xyz.Vec(1, 1, 1).lineangle(newa3) < 0.0001
assert xyz.Vec(0, 0, 1).lineangle(newa2) < 0.0001
cellsize = abs((newc2.y - newc2.x)) * 4.0
print "\nAXIS ERROR is:", newa2.lineangle(xyz.Vec(0, 0, 1)) + newa3.lineangle(xyz.Vec(1, 1, 1))
print "CELL SIZE", cellsize, '\n'
cmd.save(".tmp.pdb", sele + " and chain A")
with open(fn, 'w') as out:
out.write("CRYST1 %7.3f %7.3f %7.3f 90.00 90.00 90.00 I 21 3\n" % (
(cellsize,) * 3))
os.system("cat .tmp.pdb >> %s" % fn)
os.system("rm .tmp.pdb")
# x1 = xyz.rotation_around_degrees(newa2,180,newc2)
# x2 = xyz.rotation_around_degrees(newa3,120,newc3)
# x3 = xyz.rotation_around_degrees(newa3,240,newc3)
# for i,X in enumerate( xyz.expand_xforms((x1,x2,x3),6) ):
# cmd.create("sub%i"%i,sele+" and chain A and not sub*")
# xform("sub%i"%i,X)
# cx,cy,cz = get_cell_bounds_orthogonal_only((x1,x2),12,com(sele+' and chain A'))
# d = abs((c2.y-c2.x)*4.0)
# cmd.hide('ev','not chain A')
# for i,X in enumerate( xyzfind_identities((x1,x2,x3),9) ):
# cmd.create("sub%ii"%i,sele+" and chain A and not sub*")
# xform("sub%ii"%i,X)
# cmd.show('spheres',"sub%ii"%i)
# (0.374996,0.500000,0.250000)
# def process_xtal_dir(d):
# for fn in glob.glob(d+"/*_I213_*.pdb"):
# makecryst1_i213(fn)
[docs]def make_cryst1_23(fn, a2in, i2, a3in, i3):
print "makecryst1_i213", fn
cmd.load(fn, 'work_prot')
sele = 'work_prot'
fn = os.path.basename(fn)
if fn.endswith(".gz"):
fn = fn[:-3]
if fn.endswith(".pdb"):
fn = fn[:-4]
fn += "_cryst1.pdb"
Dsel = sele + ' and name CA and chain A+D'
Tsel = sele + ' and name CA and chain A+B+C'
olda2 = c2axis(Dsel, chains=('A', 'D'))
olda3 = c3axis(Tsel)
axis_angle_error = abs(a2in.lineangle_degrees(
a3in) - olda2.lineangle_degrees(olda3))
if axis_angle_error > 0.01:
print "WARNING: axis angle error is", axis_angle_error
print "target axis ang", a2in.lineangle_degrees(a3in)
print "actual axis ang", olda2.lineangle_degrees(olda3)
else:
print "INFO: axis angle error is", axis_angle_error
for a2i in (0, 1, 2):
a2 = rotation_matrix_degrees() * a2in
for a3 in (-a3in, a3in):
oldc2 = com(Dsel)
oldc3 = com(Tsel)
olda2 = c2axis(Dsel, chains=('A', 'D'))
olda3 = c3axis(Tsel)
if abs(olda2.angle(olda3) - a2.angle(a3)) > 0.001:
olda2 = -olda2
print olda2.angle(olda3)
print a2.angle(a3)
assert abs(olda2.angle(olda3) - a2.angle(a3)) < 0.001
Xaln = xyz.alignvectors(olda2, olda3, a2, a3)
xform(sele, Xaln)
trans(sele, -com(Tsel))
newa2 = c2axis(Dsel, chains=('A', 'D'))
newa3 = c3axis(Tsel)
assert a2.lineangle(newa2) < 0.001
assert a3.lineangle(newa3) < 0.001
# newc3 = com(Tsel)
# assert newc3.length() < 0.00001
# shift in unit cell
# replace me
newc2 = com(Dsel)
# if newc2.x > 0:
# print "skip c2.x > 0"
# continue
print "found c2.x < 0"
trans(sele, -newc2.x)
cellsize = abs((newc2.y - newc2.x)) * 4.0
# fix me!
print com(Dsel)
print a2
print i2 * cellsize
assert point_line_distance(com(Dsel), a2, i2 * cellsize) < 0.0001
assert point_line_distance(com(Tsel), a3, i3 * cellsize) < 0.0001
print "\nAXIS ERROR is:", newa2.lineangle(a2) + newa3.lineangle(a3)
print "CELL SIZE", cellsize, '\n'
# save
cmd.save(".tmp.pdb", sele + " and chain A")
with open(fn, 'w') as out:
out.write("CRYST1 %7.3f %7.3f %7.3f 90.00 90.00 90.00 I 21 3\n" % (
(cellsize,) * 3))
os.system("cat .tmp.pdb >> %s" % fn)
os.system("rm .tmp.pdb")
return
# x1 = xyz.rotation_around_degrees(newa2,180,newc2)
# x2 = xyz.rotation_around_degrees(newa3,120,newc3)
# x3 = xyz.rotation_around_degrees(newa3,240,newc3)
# for i,X in enumerate( xyz.expand_xforms((x1,x2,x3),6) ):
# cmd.create("sub%i"%i,sele+" and chain A and not sub*")
# xform("sub%i"%i,X)
# cx,cy,cz = get_cell_bounds_orthogonal_only((x1,x2),12,com(sele+' and chain A'))
# d = abs((c2.y-c2.x)*4.0)
# cmd.hide('ev','not chain A')
# for i,X in enumerate( xyzfind_identities((x1,x2,x3),9) ):
# cmd.create("sub%ii"%i,sele+" and chain A and not sub*")
# xform("sub%ii"%i,X)
# cmd.show('spheres',"sub%ii"%i)
# (0.374996,0.500000,0.250000)
# def process_xtal_dir(d):
# for fn in glob.glob(d+"/*_I213_*.pdb"):
# makecryst1_i213(fn)
# def makecryst1_i213(fn): return make_cryst1_23(fn,Uz,Uy/4,Vec(1,1,1),V0)
[docs]def nulltest():
"""
>>> print "foo"
foo
"""
[docs]def load_tests(loader, tests, ignore):
tests.addTests(doctest.DocTestSuite())
return tests
if __name__ == '__main__':
import doctest
r = doctest.testmod()
print r