## FIRST PART

 # fetch, remove hydrogens and heteroatoms
 # and align NMR and X-RAY structures
 fetch 1TVG, xray
 fetch 1XPW, nmr
 remove resn hoh
 remove hetatm
 align nmr, xray

 # create beautiful pic
 color salmon, xray
 color teal, nmr
 split_states nmr, 
 set cartoon_transparency, 0.6, nmr*
 bg_color white
 set ray_trace_color, black            # outline color
 set ray_trace_mode,  1                # turn outline on
 ray 600, 600                          # for better quality
 png all.png


## SECOND PART

 # select backbone
 select bb, (name c+o+n+ca and xray)

 # select the 1st hydrogen bond
 select don, (elem n and xray and bb)
 select acc, (elem o,n and xray and bb)
 dist hbonds, don, acc, 3.5, mode=2
 hide labels
 select hb1, (i. 64 or i. 103) and xray

 # select the 2nd & the 3rd hydrogen bonds
 select sc_don, (elem n and xray and not bb)
 select sc_acc, (elem o,n and xray and not bb)
 dist hbonds_sc, sc_don, sc_acc, 3.5, mode=2
 hide labels
 select hb2, (i. 66 or i. 101) and xray
 select hb3, (i. 95 and i. 97) and xray

 # settings for beautiful pics
 h_add xray
 set transparency, 0.7
 set cartoon_transparency, 0.6
 set stick_radius, 0.18
 set sphere_scale, 0.2
 set label_size, 20
 set label_bg_color, white
 set ray_trace_color, black
 set ray_trace_mode=1
 set ray_trace_gain, 1

 # create pics (ex. for hbond #3)
 hide everything
 show sticks, hb3
 zoom hb3
 show cartoon, xray
 show spheres, hb3
 dist hbonds3, i. 95 and n. NE2 and xray, i. 97 and n. OE1 and xray, mode=2
 select th, (xray) and not ((xray) within 20 of hb3)    # for faster render
 hide everything, th					# hide everything outside pic

 # calculate the number, min, max, median of hbonds in NMR models
 # (ex. for hb #1)
 from numpy import median

 do = "i. 64 and n. N and "
 ac = "i. 103 and n. O and "
 radius = 3.5
 distances = []
 a = ["nmr_000" + str(i) for i in range(1,10)]
 a = a + ["nmr_00" + str(i) for i in range(10,21)]
 for model in a:
 	sele1 = do + model 
 	sele2 = ac + model
 	distances.append(cmd.distance("name", sele1, sele2, mode=2, cutoff=radius))

 print(distances)
 print(len([i for i in distances if i > 0]))
 print(min([i for i in distances if i > 0]))
 print(max([i for i in distances if i > 0]))
 print(median([i for i in distances if i > 0]))