text: from __future__ import division
  from mmtbx.monomer_library import pdb_interpretation
  from mmtbx import monomer_library
  import mmtbx.monomer_library.server
  from iotbx import pdb
  from cctbx.geometry_restraints.lbfgs import lbfgs as geometry_restraints_lbfgs
  import scitbx.lbfgs
  import libtbx.load_env
  import math
  import sys, os
  
  
  def add_ring_h(site_0,site_1,site_2,d0,alpha,beta):
    #                 H (xh,yh,zh)
    #                 |
    #                 | d0
    #                 |
    #      /_alpha    C (x0,y0,z0)   /_beta
    #                / \
    #            d1 /   \ d2
    #              /     \
    #  (x1,y1,z1) C       C (x2,y2,z2)
    #             |       |
    #             |       |
    #             |       |
    #             C       C
    #              \     /
    #               \   /
    #                \ /
    #                 C
    #
    x0 = site_0[0]
    y0 = site_0[1]
    z0 = site_0[2]
    x1 = site_1[0]
    y1 = site_1[1]
    z1 = site_1[2]
    x2 = site_2[0]
    y2 = site_2[1]
    z2 = site_2[2]
    a1 = x2-x0
    a2 = y2-y0
    a3 = z2-z0
    b1 = x1-x0
    b2 = y1-y0
    b3 = z1-z0
    c1 = b2*a3-a2*b3
    c2 = -(b1*a3-a1*b3)
    c3 = b1*a2-a1*b2
    R = -x1*c1-y1*c2-z1*c3
    d1 = math.sqrt(b1*b1+b2*b2+b3*b3)
    d2 = math.sqrt(a1*a1+a2*a2+a3*a3)
    A = d1*d0*math.cos(alpha*math.pi/180.)
    B = d2*d0*math.cos(beta*math.pi/180.)
    P = -(B+a1*x0+a2*y0+a3*z0)
    Q = -(A+b1*x0+b2*y0+b3*z0)
    t1 = b2-b1*a2/a1
    t2 = b3-b1*a3/a1
    t3 = Q-b1*P/a1
    f1 = c2-c1*a2/a1
    f2 = c3-c1*a3/a1
    f3 = R-c1*P/a1
    zh = (f1*t3/t1-f3)/(f2-f1*t2/t1)
    yh = -1./t1*(zh*t2+t3)
    xh = -1./a1*(a2*yh+a3*zh+P)
    return (xh,yh,zh)
  
  
  def add_ca_ha1_and_ha2(site_0,site_1,site_2,d0,alpha,beta):
    #            HA1      HA2
    #              \      /
    #            d0 \    / d0
    #                \  /
    #      /_alpha    CA (site_0)   /_beta
    #               /    \
    #           d1 /      \ d2
    #             /        \
    #   (site_1) C          N (site_2)
    #
    x0 = site_0[0]
    y0 = site_0[1]
    z0 = site_0[2]
    x1 = site_1[0]
    y1 = site_1[1]
    z1 = site_1[2]
    x2 = site_2[0]
    y2 = site_2[1]
    z2 = site_2[2]
    a1 = x2-x0
    a2 = y2-y0
    a3 = z2-z0
    b1 = x1-x0
    b2 = y1-y0
    b3 = z1-z0
    d1 = math.sqrt(b1*b1+b2*b2+b3*b3)
    d2 = math.sqrt(a1*a1+a2*a2+a3*a3)
    A = d1*d0*math.cos(alpha*math.pi/180.)
    B = d2*d0*math.cos(beta*math.pi/180.)
    P = -(B+a1*x0+a2*y0+a3*z0)
    Q = -(A+b1*x0+b2*y0+b3*z0)
    B1 = (a3/a1-b3/b1)/(b2/b1-a2/a1)
    B2 = (P/a1-Q/b1)/(b2/b1-a2/a1)
    B3 = a3/a1+a2/a1*B1
    B4 = P/a1+a2/a1*B2
    D1 = 1.+B1*B1+B3*B3
    D2 = 2.*B3*(B4+x0)+2.*B1*(B2-y0)-2.*z0
    D3 = (B4+x0)*(B4+x0)+(B2-y0)*(B2-y0)+z0*z0-d0*d0
    discr = math.sqrt(D2*D2-4.*D1*D3)
    coeff = (-D2 - discr)/(2.*D1)
    xh1 = -B4-B3*coeff
    yh1 = B2+B1*coeff
    zh1 = coeff
    coeff = (-D2 + discr)/(2.*D1)
    xh2 = -B4-B3*coeff
    yh2 = B2+B1*coeff
    zh2 = coeff
    return ((xh1,yh1,zh1),(xh2,yh2,zh2))
  
  def add_ca_ha1_or_ha2(site_0,site_1,site_2,site_3,d0,alpha,beta,gamma):
    #     HA1 (or HA2)    CB (or HA1 or HA2)
    #             \ /_gamma/
    #           d0 \      / d3
    #               \    /
    #      /_alpha    CA (site_0)   /_beta
    #               /    \
    #           d1 /      \ d2
    #             /        \
    #   (site_1) C          N (site_2)
    #
    x0 = site_0[0]
    y0 = site_0[1]
    z0 = site_0[2]
    x1 = site_1[0]
    y1 = site_1[1]
    z1 = site_1[2]
    x2 = site_2[0]
    y2 = site_2[1]
    z2 = site_2[2]
    x3 = site_3[0]
    y3 = site_3[1]
    z3 = site_3[2]
    a1 = x2-x0
    a2 = y2-y0
    a3 = z2-z0
    b1 = x1-x0
    b2 = y1-y0
    b3 = z1-z0
    c1 = x3-x0
    c2 = y3-y0
    c3 = z3-z0
    d1 = math.sqrt(b1*b1+b2*b2+b3*b3)
    d2 = math.sqrt(a1*a1+a2*a2+a3*a3)
    d3 = math.sqrt(c1*c1+c2*c2+c3*c3)
    A = d1*d0*math.cos(alpha*math.pi/180.)
    B = d2*d0*math.cos(beta*math.pi/180.)
    C = d3*d0*math.cos(gamma*math.pi/180.)
    P = -(B+a1*x0+a2*y0+a3*z0)
    Q = -(A+b1*x0+b2*y0+b3*z0)
    R = -(C+c1*x0+c2*y0+c3*z0)
    t1 = b2-b1*a2/a1
    t2 = b3-b1*a3/a1
    t3 = Q-b1*P/a1
    f1 = c2-c1*a2/a1
    f2 = c3-c1*a3/a1
    f3 = R-c1*P/a1
    zh = (f1*t3/t1-f3)/(f2-f1*t2/t1)
    yh = -1./t1*(zh*t2+t3)
    xh = -1./a1*(a2*yh+a3*zh+P)
    return (xh,yh,zh)
  
  def add_hhh(site_0,site_1,d0,alpha,beta):
    x0 = site_0[0]
    y0 = site_0[1]
    z0 = site_0[2]
    x1 = site_1[0]
    y1 = site_1[1]
    z1 = site_1[2]
    a1 = x1-x0
    a2 = y1-y0
    a3 = z1-z0
    d1 = math.sqrt(a1*a1+a2*a2+a3*a3)
    A = d1*d0*math.cos(alpha*math.pi/180.)
    P = -(A+a1*x0+a2*y0+a3*z0)
    b1 = z0+P/a3
    b2 = a1/a3
    b3 = a2/a3
    c1 = 1.+b2*b2
    c2 = 1.+b3*b3
    c3 = 2.*b1*b2-2.*x0
    c4 = 2.*b1*b3-2.*y0
    c5 = 2.*b2*b3
    R = -(d0*d0-b1*b1-x0*x0-y0*y0)
    e1 = c5*c5-4.*c1*c2
    e2 = 2.*c3*c5-4.*c1*c4
    e3 = c3*c3-4.*c1*R
    y1 = (-e2+math.sqrt(e2*e2-4.*e1*e3)) / (2.*e1)
    y2 = (-e2-math.sqrt(e2*e2-4.*e1*e3)) / (2.*e1)
    x1 = -(c3+c5*y1)/(2.*c1)
    x2 = -(c3+c5*y2)/(2.*c1)
    z1 = -(P+a1*x1+a2*y1) / a3
    z2 = -(P+a1*x2+a2*y2) / a3
    return (x1,y1,z1), add_ca_ha1_and_ha2(site_0,site_1,(x1,y1,z1),d0,alpha,beta)
    #print x2,y2,z2
    #add_ca_ha1_and_ha2(site_0,site_1,(x2,y2,z2),d0,alpha,alpha)
  
  def add_arg_like_h(site_0,site_1,site_2,d,alpha,flag):
    F = site_0[0]
    G = site_0[1]
    H = site_0[2]
    a1,a2,a3 = site_1[0]-site_0[0],site_1[1]-site_0[1],site_1[2]-site_0[2]
    b1,b2,b3 = site_2[0]-site_1[0],site_2[1]-site_1[1],site_2[2]-site_1[2]
    c1,c2,c3 = site_2[0]-site_0[0],site_2[1]-site_0[1],site_2[2]-site_0[2]
    d1 = math.sqrt(a1*a1+a2*a2+a3*a3)
    A = a2/a1
    B = a3/a1
    C = (-d*d1*math.cos(alpha*math.pi/180.)-F*a1-G*a2-H*a3)/a1
    AA = A*A
    BB = B*B
    CC = C*C
    FF = F*F
    GG = G*G
    HH = H*H
    M = d*d
    f1 = a2*c3-a3*c2
    f2 = -(a1*c3-c1*a3)
    f3 = a1*c2-c1*a2 ;
    f4 = -F*f1 - G*f2 - H*f3
    P = f2/f1
    Q = f3/f1
    R = f4/f1
    PP = P*P
    QQ = Q*Q
    RR = R*R
    xh1 = (-2.*(A-P)*(C*(PP+QQ)+BB*(P*(G-F*P)+R)+AA*(Q*(H-F*Q)+R)-A*(C*P+B*H*P+  \
         B*G*Q-2.*B*F*P*Q+H*P*Q-G*QQ+P*R)+B*(H*PP-Q*(C+G*P+R)))+(-B*P+A*Q)*      \
         math.sqrt(-4.*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ))*((FF+ \
         GG+HH-M)*(A-P)**2+2.*C*(A-P)*(G-F*P)+CC*(1.+PP)+2.*(A*F-G)*(A-P)*R-2.*C*\
         (1.+A*P)*R+(1.+AA)*R**2)+4.*(-2.*A*H*P+H*PP+C*Q-G*P*Q+A*(G+(C+F)*P)*Q-  \
         Q*R+B*((A-P)*(-G+F*P)-C*(1.+PP)+R+A*P*R)+AA*(H-Q*(F+R)))**2))/(2.*(A-   \
         P)*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ)))
    yh1 = (-2.*(A-P)*(-(C*P)+H*P*Q-G*QQ-B*(-2.*G*Q+P*(H+(C+F)*Q))+P*R+A*((B-Q)*    \
         (H-F*Q)+C*(1.+QQ)-(1.+B*Q)*R)+BB*(-G+P*(F+R)))+(B-Q)*math.sqrt(-4.*(PP+   \
         BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ))*((FF+GG+HH-M)*(A-P)**2+2.*\
         C*(A-P)*(G-F*P)+CC*(1.+PP)+2.*(A*F-G)*(A-P)*R-2.*C*(1.+A*P)*R+(1.+AA)*RR)+\
         4.*(-2.*A*H*P+H*PP+C*Q-G*P*Q+A*(G+(C+F)*P)*Q-Q*R+B*((A-P)*(-G+F*P)-C*     \
         (1.+PP)+R+A*P*R)+AA*(H-Q*(F+R)))**2))/(2.*(A-P)*(PP+BB*(1.+PP)-2.*B*Q+QQ- \
         2.*A*(P+B*P*Q)+AA*(1.+QQ)))
    zh1 = -(-2.*H*(A-P)**2-2.*(C-(A*F-G)*(A-P)+A*C*P)*Q+2.*(1.+AA)*Q*R+2.*B*((A-P)*\
         (G-F*P)+C*(1.+PP)-(1.+A*P)*R)+math.sqrt(-4.*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*\
         (P+B*P*Q)+AA*(1.+QQ))*((FF+GG+HH-M)*(A-P)**2+2.*C*(A-P)*(G-F*P)+CC*(1.+   \
         PP)+2.*(A*F-G)*(A-P)*R-2.*C*(1.+A*P)*R+(1.+AA)*RR)+4.*(-2.*A*H*P+H*PP+C*Q-\
         G*P*Q+A*(G+(C+F)*P)*Q-Q*R+B*((A-P)*(-G+F*P)-C*(1.+PP)+R+A*P*R)+AA*(H-Q*(F+\
         R)))**2))/(2.*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ)))
    if(flag > 1):
      xh2 = (-2.*(A-P)*(C*(PP+QQ)+BB*(P*(G-F*P)+R)+AA*(Q*(H-F*Q)+R)-A*(C*P+B*H*P+  \
           B*G*Q-2.*B*F*P*Q+H*P*Q-G*QQ+P*R)+B*(H*PP-Q*(C+G*P+R)))-(-B*P+A*Q)*      \
           math.sqrt(-4.*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ))*((FF+ \
           GG+HH-M)*(A-P)**2+2.*C*(A-P)*(G-F*P)+CC*(1.+PP)+2.*(A*F-G)*(A-P)*R-2.*C*\
           (1.+A*P)*R+(1.+AA)*R**2)+4.*(-2.*A*H*P+H*PP+C*Q-G*P*Q+A*(G+(C+F)*P)*Q-  \
           Q*R+B*((A-P)*(-G+F*P)-C*(1.+PP)+R+A*P*R)+AA*(H-Q*(F+R)))**2))/(2.*(A-   \
           P)*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ)))
      yh2 = (-2.*(A-P)*(-(C*P)+H*P*Q-G*QQ-B*(-2.*G*Q+P*(H+(C+F)*Q))+P*R+A*((B-Q)*    \
           (H-F*Q)+C*(1.+QQ)-(1.+B*Q)*R)+BB*(-G+P*(F+R)))-(B-Q)*math.sqrt(-4.*(PP+   \
           BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ))*((FF+GG+HH-M)*(A-P)**2+2.*\
           C*(A-P)*(G-F*P)+CC*(1.+PP)+2.*(A*F-G)*(A-P)*R-2.*C*(1.+A*P)*R+(1.+AA)*RR)+\
           4.*(-2.*A*H*P+H*PP+C*Q-G*P*Q+A*(G+(C+F)*P)*Q-Q*R+B*((A-P)*(-G+F*P)-C*     \
           (1.+PP)+R+A*P*R)+AA*(H-Q*(F+R)))**2))/(2.*(A-P)*(PP+BB*(1.+PP)-2.*B*Q+QQ- \
           2.*A*(P+B*P*Q)+AA*(1.+QQ)))
      zh2 = -(-2.*H*(A-P)**2-2.*(C-(A*F-G)*(A-P)+A*C*P)*Q+2.*(1.+AA)*Q*R+2.*B*((A-P)*\
           (G-F*P)+C*(1.+PP)-(1.+A*P)*R)-math.sqrt(-4.*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*\
           (P+B*P*Q)+AA*(1.+QQ))*((FF+GG+HH-M)*(A-P)**2+2.*C*(A-P)*(G-F*P)+CC*(1.+   \
           PP)+2.*(A*F-G)*(A-P)*R-2.*C*(1.+A*P)*R+(1.+AA)*RR)+4.*(-2.*A*H*P+H*PP+C*Q-\
           G*P*Q+A*(G+(C+F)*P)*Q-Q*R+B*((A-P)*(-G+F*P)-C*(1.+PP)+R+A*P*R)+AA*(H-Q*(F+\
           R)))**2))/(2.*(PP+BB*(1.+PP)-2.*B*Q+QQ-2.*A*(P+B*P*Q)+AA*(1.+QQ)))
    if(flag == 1):
      return (xh1,yh1,zh1)
    if(flag == 2):
      return (xh1,yh1,zh1),(xh2,yh2,zh2)
  
  def run(file_name):
    print "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
    mon_lib_srv = monomer_library.server.server()
    ener_lib = monomer_library.server.ener_lib()
    processed = monomer_library.pdb_interpretation.process(
      mon_lib_srv           = mon_lib_srv,
      ener_lib              = ener_lib,
      file_name             = file_name,
      keep_monomer_mappings = True,
      log                   = sys.stdout).all_chain_proxies
    total_missing = 0.0
    still_missing = 0.0
    still_missing_h = []
    atom_number = 0
    file_name_ =  os.path.basename(file_name)+"_h"
    file = open(file_name_,"w")
  
    for monomer_mapping in processed.all_monomer_mappings:
      atom_number = write_atoms(monomer_mapping,
                                atom_number,
                                file)
      bond_list = monomer_mapping.monomer.bond_list
      angle_list = monomer_mapping.monomer.angle_list
      missing_h = monomer_mapping.missing_hydrogen_atoms.keys()
      residue_name = monomer_mapping.monomer.chem_comp.three_letter_code
      print
      print "Residue name: ",
      print monomer_mapping.residue_name, monomer_mapping.monomer.chem_comp.three_letter_code
      print "Missing hydrogen atoms: ", missing_h
      print
      total_missing += len(missing_h)
      unknown_h = []
      n_missed = len(missing_h)-1
      while n_missed > 0:
        n_missed -= 1
        h = missing_h[0]
        for bl in bond_list:
          if(h == bl.atom_id_1 or h == bl.atom_id_2):
            for atom in (bl.atom_id_1,bl.atom_id_2):
              if(h != atom): target_atom_name = atom
            for atom_name,atom in monomer_mapping.expected_atoms.items():
              if(atom_name == target_atom_name):
                target_site = atom.xyz
            bond_dist = bl.value_dist
            format="missing %4s: bond: %4s %4s bond distance = %5.3f"
            print format % (h, bl.atom_id_1, bl.atom_id_2, bl.value_dist)
        angles = []
        for al in angle_list:
          if(h == al.atom_id_1 or h == al.atom_id_2 or h == al.atom_id_3):
            format="        %4s: angle: %4s %4s %4s = %5.3f"
            print format % (h, al.atom_id_1,al.atom_id_2,al.atom_id_3,al.value_angle)
            angles.append([al.value_angle,al.atom_id_1,al.atom_id_2,al.atom_id_3])
        bonded_to_target_site = []
        #for bl in bond_list:
        #  if(target_atom_name in (bl.atom_id_1,bl.atom_id_2)):
        #    for atom in (bl.atom_id_1,bl.atom_id_2):
        #      if(atom != target_atom_name):
        #        bonded_to_target_site.append(atom)
        #print "target_site, bonded_to_target_site ", target_atom_name, bonded_to_target_site
  ### ring hydrogens: TYR,PHE,HIS,TRP
        if(residue_name in ["PHE","TYR","HIS","TRP","ARG"]):
          ring_h = {
           "HE1":{"PHE":["CE1","CD1","CZ"],
                  "TYR":["CE1","CD1","CZ"],
                  "HIS":["CE1","ND1","NE2"],
                  "TRP":["NE1","CD1","CE2"]},
           "HE2":{"PHE":["CE2","CD2","CZ"],
                  "TYR":["CE2","CD2","CZ"],
                  "HIS":["NE2","CE1","CD2"]},
           "HE3":{"TRP":["CE3","CZ3","CD2"]},
           "HD1":{"PHE":["CD1","CE1","CG"],
                  "TYR":["CD1","CE1","CG"],
                  "HIS":["ND1","CG","CE1"],
                  "TRP":["CD1","CG","NE1"]},
           "HD2":{"PHE":["CD2","CE2","CG"],
                  "TYR":["CD2","CE2","CG"],
                  "HIS":["CD2","NE2","CG"]},
           "HH2":{"TRP":["CH2","CZ2","CZ3"]},
           "HZ" :{"PHE":["CZ","CE1","CE2"]},
           "HZ2":{"TRP":["CZ2","CE2","CH2"]},
           "HZ3":{"TRP":["CZ3","CH2","CE3"]},
           "HE" :{"ARG":["NE","CZ","CD"]}
                   }
          if(h in ring_h.keys()):
            if(residue_name in ring_h[h].keys()):
              targets = ring_h[h][residue_name]
              print "Building:", h, " ..."
              site_0 = None
              site_1 = None
              site_2 = None
              for atom_name,atom in monomer_mapping.expected_atoms.items():
                if(atom_name == targets[0]):
                  site_0 = atom.xyz
                if(atom_name == targets[1]):
                  site_1 = atom.xyz
                if(atom_name == targets[2]):
                  site_2 = atom.xyz
              alpha = angles[0][0]
              beta  = angles[1][0]
              assert site_0 is not None and site_1 is not None and site_2 is not None
              xyz = add_ring_h(site_0,site_1,site_2,bond_dist,alpha,beta)
              atom_number = write_atoms(monomer_mapping,
                                        atom_number,
                                        file,
                                        new_atom_name        = h,
                                        new_atom_coordinates = xyz)
              missing_h.remove(h)
  ### CA hydrogen(s): HA or/and HA1 or/and HA2
        if(h == "HA" or h == "HA1" or h == "HA2"):
          build_one = False
          build_two = False
          if("HA" in missing_h):
            build_one = True
          if("HA1" in missing_h and "HA2" in missing_h):
            build_two = True
          if(("HA1" in missing_h or "HA2" in missing_h) and build_two == False):
            build_one = True
          assert build_one == True and build_two == False or \
                 build_one == False and build_two == True
          if(build_two == True):
            print "Building: HA1 and HA2"
            site_0 = None
            site_1 = None
            site_2 = None
            for atom_name,atom in monomer_mapping.expected_atoms.items():
              if(atom_name == "CA"):
                site_0 = atom.xyz
              if(atom_name == "C"):
                site_1 = atom.xyz
              if(atom_name == "N"):
                site_2 = atom.xyz
            if(site_0 is not None and site_1 is not None and site_2 is not None):
              for angle in angles:
                if("C" in angle and "CA" in angle): alpha = angle[0]
                if("N" in angle and "CA" in angle): beta  = angle[0]
              xyz = add_ca_ha1_and_ha2(site_0,site_1,site_2,bond_dist,alpha,beta)
              atom_number = write_atoms(monomer_mapping,
                                        atom_number,
                                        file,
                                        new_atom_name        = "HA1",
                                        new_atom_coordinates = xyz[0])
              atom_number = write_atoms(monomer_mapping,
                                        atom_number,
                                        file,
                                        new_atom_name        = "HA2",
                                        new_atom_coordinates = xyz[1])
              missing_h.remove("HA1")
              missing_h.remove("HA2")
          if(build_one == True):
            print "Building: ", h
            site_0 = None
            site_1 = None
            site_2 = None
            site_3 = None
            for atom_name,atom in monomer_mapping.expected_atoms.items():
              if(atom_name == "CA"):
                site_0 = atom.xyz
              if(atom_name == "C"):
                site_1 = atom.xyz
              if(atom_name == "N"):
                site_2 = atom.xyz
              if(atom_name == "CB" or atom_name == "HA1" or atom_name == "HA2"):
                site_3 = atom.xyz
                cb_or_ha1_or_ha2 = atom_name
            if(site_0 is not None and site_1 is not None and site_2 is not None):
              for angle in angles:
                if("C" in angle and "CA" in angle): alpha  = angle[0]
                if("N" in angle and "CA" in angle): beta   = angle[0]
                if("CA" in angle and cb_or_ha1_or_ha2 in angle): gamma = angle[0]
              xyz = add_ca_ha1_or_ha2(site_0,site_1,site_2,site_3,bond_dist,alpha,beta,gamma)
              atom_number = write_atoms(monomer_mapping,
                                        atom_number,
                                        file,
                                        new_atom_name        = h,
                                        new_atom_coordinates = xyz)
              missing_h.remove(h)
  ### CB hydrogen(s): HB1 or/and HB2 and similar
        cb_like_h = {
         "HB1":[["CB","CA","CG","HB2"],["CB","CA","OG","HB2"],["CB","CA","SG","HB2"]],
         "HB2":[["CB","CA","CG","HB1"],["CB","CA","OG","HB1"],["CB","CA","SG","HB1"]],
         "HG1":[["CG","CB","CD","HG2"],],
         "HG2":[["CG","CB","CD","HG1"],],
         "HD1":[["CD","CG","NE","HD2"],["CD","CG","CE","HD2"]],
         "HD2":[["CD","CG","NE","HD1"],["CD","CG","CE","HD1"]],
         "HE1":[["CE","CD","NZ","HE2"],],
         "HE2":[["CE","CD","NZ","HE1"],],
         "HG" :[["CG","CB","CD1","CD2"],],
         "HB" :[["CB","CG2","CA","CG1"],["CB","CG2","CA","OG1"]]
                    }
        if(h in cb_like_h.keys()):
          if(h in missing_h):
            targets_ = cb_like_h[h]
            for targets in targets_:
              build_one = False
              build_two = False
              if(h in missing_h and targets[3] in missing_h):
                build_two = True
              if((h in missing_h or targets[3] in missing_h) and build_two == False):
                build_one = True
              if(build_one == True and build_two == False or \
                                           build_one == False and build_two == True):
                if(build_two == True):
                  print "Building: ", h," and ",targets[3]
                  site_0 = None
                  site_1 = None
                  site_2 = None
                  for atom_name,atom in monomer_mapping.expected_atoms.items():
                    if(atom_name == targets[0]):
                      site_0 = atom.xyz
                    if(atom_name == targets[1]):
                      site_1 = atom.xyz
                    if(atom_name == targets[2]):
                      site_2 = atom.xyz
                  if(site_0 is not None and site_1 is not None and site_2 is not None):
                    for angle in angles:
                      if(targets[1] in angle and targets[0] in angle): alpha = angle[0]
                      if(targets[0] in angle and targets[2] in angle): beta  = angle[0]
                    xyz = add_ca_ha1_and_ha2(site_0,site_1,site_2,bond_dist,alpha,beta)
                    atom_number = write_atoms(monomer_mapping,
                                        atom_number,
                                        file,
                                        new_atom_name        = h,
                                        new_atom_coordinates = xyz[0])
                    atom_number = write_atoms(monomer_mapping,
                                        atom_number,
                                        file,
                                        new_atom_name        = targets[3],
                                        new_atom_coordinates = xyz[1])
                    missing_h.remove(h)
                    missing_h.remove(targets[3])
                if(build_one == True):
                  print "Building: ", h
                  site_0 = None
                  site_1 = None
                  site_2 = None
                  site_3 = None
                  alpha  = None
                  beta   = None
                  gamma  = None
                  for atom_name,atom in monomer_mapping.expected_atoms.items():
                    if(atom_name == targets[0]):
                      site_0 = atom.xyz
                    if(atom_name == targets[1]):
                      site_1 = atom.xyz
                    if(atom_name == targets[2]):
                      site_2 = atom.xyz
                    if(atom_name == h or atom_name == targets[3]):
                      site_3 = atom.xyz
                      cb_or_ha1_or_ha2 = atom_name
                  if(site_0 is not None and site_1 is not None and site_2 is not None and site_3 is not None):
                    for angle in angles:
                      if(targets[1] in angle and targets[0] in angle): alpha = angle[0]
                      if(targets[0] in angle and targets[2] in angle): beta  = angle[0]
                      if(targets[0] in angle and cb_or_ha1_or_ha2 in angle): gamma = angle[0]
                    if(alpha is not None and beta is not None and gamma is not None):
                      xyz = add_ca_ha1_or_ha2(site_0,site_1,site_2,site_3,bond_dist,alpha,beta,gamma)
                      atom_number = write_atoms(monomer_mapping,
                                        atom_number,
                                        file,
                                        new_atom_name        = h,
                                        new_atom_coordinates = xyz)
                      missing_h.remove(h)
  ### HT1,HT2,HT3 like hydrogens; rotational optimization is necessary
        hhh_like_h = {
         "HB1" :[["CB","CA","HB2","HB3"],],
         "HB2" :[["CB","CA","HB1","HB3"],],
         "HB3" :[["CB","CA","HB1","HB2"],],
         "HD11":[["CD1","CG","HD12","HD13"],["CD1","CG1","HD12","HD13"]],
         "HD12":[["CD1","CG","HD11","HD13"],["CD1","CG1","HD11","HD13"]],
         "HD13":[["CD1","CG","HD11","HD12"],["CD1","CG1","HD11","HD12"]],
         "HD21":[["CD2","CG","HD22","HD23"],],
         "HD22":[["CD2","CG","HD21","HD23"],],
         "HD23":[["CD2","CG","HD21","HD22"],],
         "HG21":[["CG2","CB","HG22","HG23"],],
         "HG22":[["CG2","CB","HG21","HG23"],],
         "HG23":[["CG2","CB","HG21","HG22"],],
         "HG11":[["CG1","CB","HG12","HG13"],],
         "HG12":[["CG1","CB","HG11","HG13"],],
         "HG13":[["CG1","CB","HG11","HG12"],],
         "HZ1" :[["NZ","CE","HZ2","HZ3"],],
         "HZ2" :[["NZ","CE","HZ1","HZ3"],],
         "HZ3" :[["NZ","CE","HZ1","HZ2"],]
                     }
        hhh_residues = ["ALA","LEU","ILE","VAL","THR","LYS"]
        if(h in hhh_like_h.keys() and residue_name in hhh_residues and h in missing_h):
          targets = hhh_like_h[h]
          for target in targets:
            site_0 = None
            site_1 = None
            for atom_name,atom in monomer_mapping.expected_atoms.items():
              if(atom_name == target[0]):
                site_0 = atom.xyz
              if(atom_name == target[1]):
                site_1 = atom.xyz
            alpha = angles[0][0]
            beta  = angles[2][0]
            if(site_0 is not None and site_1 is not None):
              print "Building:", h, " ..."
              xyz = add_hhh(site_0,site_1,bond_dist,alpha,beta)
              for atom_name, coordinates in zip((h,target[2],target[3]),(xyz[0],xyz[1][0],xyz[1][1])):
                atom_number = write_atoms(monomer_mapping,
                                          atom_number,
                                          file,
                                          new_atom_name        = atom_name,
                                          new_atom_coordinates = coordinates)
              if(h in missing_h): missing_h.remove(h)
              if(target[2] in missing_h): missing_h.remove(target[2])
              if(target[3] in missing_h): missing_h.remove(target[3])
  ### add ARG-like NH's:
        arg_like_h = {
         "HH11":[["NH1","CZ","NH2","HH12",2],],
         "HH12":[["NH1","CZ","NH2","HH11",2],],
         "HH21":[["NH2","CZ","NH1","HH22",2],],
         "HH22":[["NH2","CZ","NH1","HH21",2],],
         "HG"  :[["OG","CB","CA","HG",1],["SG","CB","CA","HG",1]] ,
         "HH"  :[["OH","CZ","CE1","HH",1],]
                     }
        residue_names = ["ARG","SER","TYR","CYS"]
        if(h in arg_like_h.keys() and residue_name in residue_names):
          targets = arg_like_h[h]
          for target in targets:
            site_0 = None
            site_1 = None
            site_2 = None
            for atom_name,atom in monomer_mapping.expected_atoms.items():
              if(atom_name == target[0]):
                site_0 = atom.xyz
              if(atom_name == target[1]):
                site_1 = atom.xyz
              if(atom_name == target[2]):
                site_2 = atom.xyz
            alpha = angles[0][0]
            if(site_0 is not None and site_1 is not None and site_2 is not None):
              print "Building:", h, " ..."
              xyz = add_arg_like_h(site_0,site_1,site_2,bond_dist,alpha,flag = target[4])
              if(target[4] == 1):
                atom_number = write_atoms(monomer_mapping,
                                          atom_number,
                                          file,
                                          new_atom_name        = h,
                                          new_atom_coordinates = xyz)
              if(target[4] == 2):
                atom_number = write_atoms(monomer_mapping,
                                          atom_number,
                                          file,
                                          new_atom_name        = h,
                                          new_atom_coordinates = xyz[0])
                atom_number = write_atoms(monomer_mapping,
                                          atom_number,
                                          file,
                                          new_atom_name        = target[3],
                                          new_atom_coordinates = xyz[1])
              if(h in missing_h): missing_h.remove(h)
              if(target[3] in missing_h): missing_h.remove(target[3])
  
  ### UNKNOWN H:
        if(h in missing_h):
          print "Unknown hydrogen type: ", h
          missing_h.remove(h)
          unknown_h.append(h)
        n_missed = len(missing_h)
  
      print
      print "Still missing: ", residue_name,unknown_h
      still_missing_h.append((residue_name,unknown_h))
      still_missing += len(unknown_h)
  
    print "Build ", int(total_missing - still_missing),
    print " from ", int(total_missing), " % = ",
    print (total_missing - still_missing)*100./total_missing
    print
    for item in still_missing_h:
      print item
  
    file.close()
    file = open(file_name_,"r")
    processed_pdb_file = pdb_interpretation.process(
      mon_lib_srv = mon_lib_srv,
      ener_lib = ener_lib,
      file_name = file_name_,
      strict_conflict_handling = False,
      crystal_symmetry = processed.pdb_inp.crystal_symmetry(),
      force_symmetry = True,
      log = sys.stdout)
    file.close()
    assert processed_pdb_file.xray_structure() is not None
    assert processed_pdb_file.geometry_restraints_manager() is not None
    xray_structure = processed_pdb_file.xray_structure()
    geometry_restraints_manager = processed_pdb_file.geometry_restraints_manager()
    regularize_model(xray_structure = xray_structure,
                     geometry_restraints_manager = geometry_restraints_manager,
                     max_iterations = 50000)
    processed_pdb_file.all_chain_proxies.pdb_atoms.set_xyz(
      new_xyz=xray_structure.sites_cart())
    processed_pdb_file.all_chain_proxies.pdb_hierarchy.write_pdb_file(
      file_name="out.pdb")
  
  def write_atoms(monomer_mapping,
                  atom_number,
                  file_object,
                  new_atom_name = None,
                  new_atom_coordinates = None,
                  crystal_symmetry = None):
    assert new_atom_name is None and new_atom_coordinates is None or \
           new_atom_name is not None and new_atom_coordinates is not None
    if(crystal_symmetry is not None):
      print >> file_object, pdb.format_cryst1_record(
                                               crystal_symmetry=crystal_symmetry)
      print >> file_object, pdb.format_scale_records(
                                          unit_cell=crystal_symmetry.unit_cell())
    if(new_atom_name is None and new_atom_coordinates is None):
      for atom_name,atom in monomer_mapping.expected_atoms.items():
        assert atom_name.strip() == atom.name.strip()
        atom_number += 1
        orig = atom.serial
        atom.serial = "%5d" % atom_number
        print >> file_object, atom.format_atom_record()
        atom.serial = orig
    else:
      atom = monomer_mapping.expected_atoms.values()[0]
      atom_number += 1
      orig = atom.serial, atom.name, atom.xyz, atom.occ, atom.b
      atom.serial = "%5d" % atom_number
      atom.name = new_atom_name
      atom.xyz = new_atom_coordinates
      atom.occ = 1
      atom.b = 0
      print >> file_object, atom.format_atom_record()
      atom.serial, atom.name, atom.xyz, atom.occ, atom.b = orig
    return atom_number
  
  def regularize_model(xray_structure,
                       geometry_restraints_manager,
                       max_iterations):
    sites_cart = xray_structure.sites_cart()
    minimized = geometry_restraints_lbfgs(
      sites_cart                  = sites_cart,
      geometry_restraints_manager = geometry_restraints_manager,
      lbfgs_termination_params    = scitbx.lbfgs.termination_parameters(
        max_iterations = max_iterations))
    xray_structure.set_sites_cart(sites_cart)
    print
    print "Energies at start of minimization:"
    minimized.first_target_result.show()
    print
    print "Energies after minimization:"
    minimized.final_target_result.show()
  
  
  if (__name__ == "__main__"):
    files = ["arg.pdb","lys.pdb","ala.pdb","gly.pdb","his.pdb","ile.pdb","leu.pdb",
             "phe.pdb","tyr.pdb","trp.pdb","thr.pdb","val.pdb","cys.pdb","ser.pdb"]
    if(len(sys.argv) > 1):
      run(sys.argv[1])
    else:
      pdb_dir = libtbx.env.find_in_repositories(
        relative_path="phenix_regression/hydrogens",
        test=os.path.isdir)
      if (pdb_dir is None):
        print "Skipping build_hydrogens.run(): input files not available"
      else:
        for file in files:
          run(os.path.join(pdb_dir, file))
  