Answers.Body: <p>Not sure but here <a href="https://biopython.org/docs/1.75/api/Bio.PDB.Residue.html" rel="nofollow noreferrer">https://biopython.org/docs/1.75/api/Bio.PDB.Residue.html</a>   Bio.PDB.Residue module Your can find :</p>
  <p>add(self, atom) :</p>
  <pre><code>      Add an Atom object.
  </code></pre>
  <p>Checks for adding duplicate atoms, and raises a PDBConstructionException if so.</p>
  <p>Can’t try it just now but who knows it may work</p>
  <p>Googling around I found out that if you plan to add Hydrogens to PDB files [I believe nowadays PDBx is the repository grade standard] you could do thar following PDB format rules. I am copying them from UCSF Chimera website (<a href="https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens" rel="nofollow noreferrer">https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens</a>):</p>
  <pre><code>
  Hydrogen Atoms
  
  In brief, conventions for hydrogen atoms in version 3.0 PDB format are as follows:
  
  Hydrogen atom records follow the records of all other atoms of a particular residue.
  A hydrogen atom name starts with H. The next part of the name is based on the name of the connected nonhydrogen atom. For example, in amino acid residues, H is followed by the remoteness indicator (if any) of the connected atom, followed by the branch indicator (if any) of the connected atom; if more than one hydrogen is connected to the same atom, an additional digit is appended so that each hydrogen atom will have a unique name. Hydrogen atoms in standard nucleotides and amino acids (other than the rarely seen HXT) are named according to the IUPAC recommendations (Markley et al., Pure Appl Chem 70:117 (1998)). Names of hydrogen atoms in HETATM residues are determined in a similar fashion.
  If the name of a hydrogen has four characters, it is left-justified starting in column 13; if it has fewer than four characters, it is left-justified starting in column 14. 
  
  In the following excerpt from entry 1vm3, atom H is attached to atom N. Atom HA is attached to atom CA; the remoteness indicator A is the same for these atoms. Two hydrogen atoms are connected to CB, one is connected to CG, three are connected to CD1, and three are connected to CD2.
  
  ATOM     10  N   LEU A   2       4.595   6.365   3.756  1.00  0.00           N
  ATOM     11  CA  LEU A   2       4.471   5.443   2.633  1.00  0.00           C
  ATOM     12  C   LEU A   2       5.841   5.176   2.015  1.00  0.00           C
  ATOM     13  O   LEU A   2       6.205   4.029   1.755  1.00  0.00           O
  ATOM     14  CB  LEU A   2       3.526   6.037   1.578  1.00  0.00           C
  ATOM     15  CG  LEU A   2       2.790   4.919   0.823  1.00  0.00           C
  ATOM     16  CD1 LEU A   2       3.803   3.916   0.262  1.00  0.00           C
  ATOM     17  CD2 LEU A   2       1.817   4.196   1.769  1.00  0.00           C
  ATOM     18  H   LEU A   2       4.169   7.246   3.704  1.00  0.00           H
  ATOM     19  HA  LEU A   2       4.063   4.514   2.992  1.00  0.00           H
  ATOM     20  HB2 LEU A   2       2.804   6.675   2.065  1.00  0.00           H
  ATOM     21  HB3 LEU A   2       4.099   6.623   0.873  1.00  0.00           H
  ATOM     22  HG  LEU A   2       2.234   5.353   0.004  1.00  0.00           H
  ATOM     23 HD11 LEU A   2       4.648   4.447  -0.148  1.00  0.00           H
  ATOM     24 HD12 LEU A   2       3.334   3.331  -0.516  1.00  0.00           H
  ATOM     25 HD13 LEU A   2       4.137   3.260   1.052  1.00  0.00           H
  ATOM     26 HD21 LEU A   2       0.941   3.892   1.216  1.00  0.00           H
  ATOM     27 HD22 LEU A   2       1.522   4.860   2.568  1.00  0.00           H
  ATOM     28 HD23 LEU A   2       2.296   3.323   2.188  1.00  0.00           H
  
  
  <span class="math-container">```</span>
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T15:50:26.117Z
Answers.Text: As you can see from the algorithm, adding Hydrogen is not that straightforward. :D If it were that simple, I wouldn't have posted this question!!
Answers.UserDisplayName: user366312
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-08-24T06:16:00.670Z
Answers.Text: Are you talking about the geometry of the missing H or about the protonation state of a lateral chain or about polar hyrogens ?
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-10-06T20:40:44.583Z
Answers.Text: see here for H atom names for Phenilalanine , bond lenght and angle have to be calculated from the atom coords https://www.ebi.ac.uk/pdbe-srv/pdbechem/atom/list/PHE or tabulated here https://www.ebi.ac.uk/pdbe-srv/pdbechem/bond/list/PHE
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T15:40:23.270Z
Answers.LastActivityDate: 2021-11-01T21:52:10.323Z
Answers.LastEditDate: 2021-11-01T21:52:10.323Z
Answers.UserDisplayName: pippo1980
Answers.Body: <p>The other answer is correct. But I thought I'd give a few pointers in how one can understand how to do something with a given Python module —teach a man to fish kind of thing...</p>
  <p>In a Jupyter Notebook or IPython shell (don't), you can do:</p>
  <ul>
  <li><code>help(function)</code> print out the docstring</li>
  <li><code>dir(object)</code> shows all the public and magic attributes and methods —note a class will not show the instance methods</li>
  <li><code>type(object)</code> or <code>object.__class__</code> gives you what the class is, but <code>object.__class__.__mro__</code> tells you what it inherits.</li>
  <li>in the standard library module <code>inspect</code> is the function <code>inspect.getsource(function)</code>, which for a notebook is messy due to newlines so needs to  be <code>print_code = lambda fun: print(inspect.getsource(fun))</code> and <code>print_code(fun)</code>.</li>
  </ul>
  <p>In your case you read a file...</p>
  <pre><code>from Bio import PDB
  parser = PDB.PDBParser()
  structure : PDB.Structure.Structure = parser.get_structure(&quot;tryp-cage&quot;, &quot;1L2Y.first.pdb&quot;)
  </code></pre>
  <p>So <code>parser.get_structure</code> &quot;magically&quot; makes a structure, eh? This means it can add stuff! As mentioned <code>print_code(parser.get_structure)</code> will reveal its secrets, or within it call some <code>self._private_fun()</code>, which can be inspected (<code>self</code> = <code>parser</code>) all the way down the rabbit hole. Until you find a nice constructor that declares nicely an atom instance or whatever you ever fancy searching for —that is the great thing of open source code.</p>
  <p>In the specific case of <code>PDB</code> module, everything inherits <code>Entity</code>, which has the attribute <code>child_list</code>, which is what the iterators (<code>__iter__</code> is the topmost) read without making a copy. The end of the above rabbit hole may just add a child <code>Entity</code>-inheriting class instance to this list, but it also do other things that is worth keeping an eye out.</p>
  <p>In the case of <code>residue.add(atom)</code>, creating the <code>Bio.PDB.Atom.Atom</code> instance is slightly problematic as there's no typehinting or <code>:type xxx:</code> declarations in the docstring. But doing <code>type(atom.bfactor)</code> etc. on the attributes of an existing atom will reveal how to do it. One the serial number is a PDB thing, which is sequential for the model.</p>
  <pre><code>last_serial = list(structure.get_residues())[-1].child_list[-1].serial_number + 1
  residue.add(PDB.Atom.Atom(name=' H  ', coord=np.array([1,1,1]), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number=last_serial+1,element='H'))
  </code></pre>
  <p>The coordinates will not be <code>[1,1,1]</code>, but the point that you have calculate. Say for <code>H</code>, this would be coplanar with <code>C</code>, <code>CA</code> and <code>N</code>, that is 1 and smidge Å away from <code>N</code> at an angle not quite in the middle as found in amide bonds.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T16:10:02.260Z
Answers.Text: If the other answer is correct, what do you think about the algorithm I posted? This algorithm is taken from the c++ source code I have in my hand.
Answers.UserDisplayName: user366312
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T16:11:17.377Z
Answers.Text: As you would realize, I can write Python and C++ program, but I am just confused about hydrogens. :)
Answers.UserDisplayName: user366312
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T16:45:30.790Z
Answers.Text: Biopython documentation is not that rich especially for ones like me that use google as primary/unique docs search
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T17:04:49.863Z
Answers.Text: I have expanded the answer with an example, showing that what I wrote is how one finds the correct solution without documentation. In terms of how to do the maths, you have to do that in numpy as its simple maths —I do not know if the C++ code is 100% correct as it is projecting on the line between C-CA midpoint and N, which I am not sure is the correct angle. Totally easy to verify though
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T19:09:50.633Z
Answers.Text: I would be curious about how Bio.PDB.PDBIO module (https://biopython.org/docs/1.74/api/Bio.PDB.PDBIO.html) would save such pdb with added H atoms, like : https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens or something else ?
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T20:37:52.250Z
Answers.Text: @MatteoFerla Hi, sorry to bother, would you be so kind to point me to a PDB (x-ray) into https://www.rcsb.org/ , that contains H atoms. Don't know how to search for them, or some other PDB compliant file with Hs, to play with biopython I am curious about it, I know I should learn PDBx by now but is far more complicated to me
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T21:54:23.743Z
Answers.Text: dont worry got it from http://biophysics.cs.vt.edu/ even if it marks all H as atom 0, without renumbering everything from scratch
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T22:25:23.130Z
Answers.Text: Neutron diffraction structures have protons experimentally. Say PDB:5MOP —A beauty as it has a hydronium or two. For small molecules x-ray crystallography can detect protons but the database for that is the [CCDC](https://ccdc.cam.ac.uk/)
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-02T09:53:53.190Z
Answers.Text: Do you have code for the correct naming of H atoms ? Does PDBx follow same rules ?
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-02T10:36:23.860Z
Answers.Text: All formats that support atom names use the same atom names and _should_ match the definition in the PBD. e.g. https://files.rcsb.org/ligands/view/ALA.cif Ligands are not always numbered by canonical index order, but amino acids follows a precise lettering scheme (greek) for heavy atoms —the annoying Val <-> Ile superposition issue is common to all tools. However, for protons there are two schemes: ` HB1` vs. `1HB `, the former is like the heavy atoms while the latter is more sensible as the number prefix is compatible with heavy atoms who have a number suffix, cf. ARG with `HH11` vs. `1HH1`
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-02T16:30:34.620Z
Answers.Text: —the annoying Val <-> Ile superposition issue is common to all tools ... what is that ?
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-02T16:49:05.627Z
Answers.Text: The CB atom is not chiral in valine, but the atom names make it "name chiral" (not a real word, but just like isotopic chirality): looking at valine with the HB atom forewards, you have CG2 on the side of the N (left) and CG1 on the C (right) in a predictable counterclockwise fashion. Isoleucine is chiral on the CB atom and the CG1 is on the left as it is more complex due to CD1. As a result they do not superpose well by atom names, but do if CG2 is mapped to CG1 and CG1: CG2.
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-10-12T18:17:44.817Z
Answers.Text: @MatteoFerla  AttributeError: 'Structure' object has no attribute 'get_residue' --> last_serial = list(structure.get_atoms())[-1].serial_number ?
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-01T16:05:00.680Z
Answers.LastActivityDate: 2021-11-01T16:58:43.740Z
Answers.LastEditDate: 2021-11-01T16:58:43.740Z
Answers.UserDisplayName: Matteo Ferla
Answers.Body: <p>OK my attempt, input pdb created in <strong>pymol</strong> with <code>fab AVGLWPT , AVGLWPT , hydro=0</code> and saved as <code>avgwpt.pdb</code> :</p>
  <pre><code>ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N  
  ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C  
  ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C  
  ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O  
  ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C  
  ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N  
  ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C  
  ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C  
  ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O  
  ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C  
  ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C  
  ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C  
  ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N  
  ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C  
  ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C  
  ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O  
  ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N  
  ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C  
  ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C  
  ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O  
  ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C  
  ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C  
  ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C  
  ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C  
  ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N  
  ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C  
  ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C  
  ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O  
  ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C  
  ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C  
  ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C  
  ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C  
  ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C  
  ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C  
  ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N  
  ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C  
  ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C  
  ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C  
  ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N  
  ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C  
  ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C  
  ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O  
  ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C  
  ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C  
  ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C  
  ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N  
  ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C  
  ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C  
  ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O  
  ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C  
  ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C  
  ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O  
  TER   
  END
  
  </code></pre>
  <p>pic :
  <a href="https://i.sstatic.net/NZNZY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/NZNZY.png" alt="enter image description here" /></a></p>
  <p>code :</p>
  <pre><code>#!/usr/bin/env python3
  # -*- coding: utf-8 -*-
  &quot;&quot;&quot;
  molecule create in pymol with fab AVGLWPT , AVGLWPT , hydro=0
  
  &quot;&quot;&quot;
  
  import Bio
  print('########## Biopython  VERSION ##########################################')
  print('\n\n BIOPYTHON_VERSION : ', Bio.__version__)
  print('###################################################################')
  
  from Bio.PDB import (
                      PDBParser,
                      Atom,
                      PDBIO
                      )
  
  
  structure_code = 'avgwpt'
  
  structure_name = structure_code
  
  parser_pdb= PDBParser(QUIET=False)    ### shows warnings
  
  structure_whole = parser_pdb.get_structure(structure_name, structure_name+&quot;.pdb&quot;)  
  
  
  last_serial = list(structure_whole.get_atoms())[-1].serial_number
  
  print(&quot; last_serial : &quot;,last_serial)
  
  for model in structure_whole:
      
      for chain in model:
          
          for resi in chain:
          
              print('1 : ', resi)
              
              print('2 : ',resi.id)
              
              try :
              
                  print('3 : ',chain[resi.id[1]-1])
                  
                  doable = True
                  
              except :
                  
                  print('3 : ','residue -1 missing !!!!!!!!!!!!!!!!!!!!!!1')
                  
                  doable = False
                  
              
              print('4 : ',[i for i in resi.get_atoms()])
              
              if doable:
                  
                  print('5 : ','ok')
              
              else:
                  
                  print('5 : ','NO')
                  
              
              if doable:
                  
                  
                  print(chain[resi.id[1]-1]['C'])
                  
                  print(chain[resi.id[1]-1]['C'].get_vector(), type(chain[resi.id[1]-1]['C'].get_vector()))
                  
                  print(chain[resi.id[1]]['N'].get_vector())
                  
                  c = chain[resi.id[1]-1]['C'].get_vector()
                  
                  n = chain[resi.id[1]]['N'].get_vector()
                  
                  ca = chain[resi.id[1]]['CA'].get_vector()
                  
                  nc = n-c
                  
                  nca = n-ca
                  
                  print('vector BN : ', nc)
                  
                  nc.normalize()
                  
                  print('vector BN normalize : ', nc)
                  
                  print('vector B2B : ', nca)
                  
                  nca.normalize()
                  
                  print('vector B2B normalize : ', nca)
                  
                  print('sum of nc + nca : ' , nsum := nc + nca)
                  
                  nsum.normalize()
                  
                  nsum**1.023
                  
                  nsum += n
                  
                  print(nsum)
                  
                  print(nsum.get_array(), type(nsum.get_array()))
                  
                  resi.add(Atom.Atom(name=' H  ', coord = nsum.get_array(), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number = last_serial ,element='H'))
              
                  last_serial += 1
  
  io=PDBIO() 
  
  io.set_structure(structure_whole ) 
  
  io.save(structure_name+&quot;_H_added.pdb&quot; , preserve_atom_numbering = True) 
  </code></pre>
  <p>output , saved as <code>avgwpt_H_added.pdb</code>:</p>
  <pre><code>ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N  
  ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C  
  ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C  
  ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O  
  ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C  
  ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N  
  ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C  
  ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C  
  ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O  
  ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C  
  ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C  
  ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C  
  ATOM     52  H   VAL     2      70.368  54.810  15.182  1.00  0.00           H  
  ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N  
  ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C  
  ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C  
  ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O  
  ATOM     53  H   GLY     3      74.741  54.249  15.182  1.00  0.00           H  
  ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N  
  ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C  
  ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C  
  ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O  
  ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C  
  ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C  
  ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C  
  ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C  
  ATOM     54  H   LEU     4      76.733  58.181  15.224  1.00  0.00           H  
  ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N  
  ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C  
  ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C  
  ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O  
  ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C  
  ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C  
  ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C  
  ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C  
  ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C  
  ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C  
  ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N  
  ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C  
  ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C  
  ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C  
  ATOM     55  H   TRP     5      81.105  57.620  15.224  1.00  0.00           H  
  ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N  
  ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C  
  ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C  
  ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O  
  ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C  
  ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C  
  ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C  
  ATOM     56  H   PRO     6      83.072  61.539  15.267  1.00  0.00           H  
  ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N  
  ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C  
  ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C  
  ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O  
  ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C  
  ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C  
  ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O  
  ATOM     57  H   THR     7      87.409  60.156  16.044  1.00  0.00           H  
  TER      58      THR     7                                                       
  END   
  
  </code></pre>
  <p>pic : <a href="https://i.sstatic.net/WmGL1.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/WmGL1.png" alt="enter image description here" /></a></p>
  <p>Not sure I got the code above (in question) right, can't read it , had to use comments as trace.</p>
  <p>Spent a while banging my head with : <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2724740/" rel="nofollow noreferrer">HAAD: A Quick Algorithm for Accurate Prediction of Hydrogen Atoms in Protein Structures</a></p>
  <p>and trying to figure out their alghoritm, [code (<a href="https://seq2fun.dcmb.med.umich.edu//HAAD/haad.f90" rel="nofollow noreferrer">https://seq2fun.dcmb.med.umich.edu//HAAD/haad.f90</a>) is in FORTRAN90] from the paper , but no luck, if anybody wants to help eplaining this bit:</p>
  <blockquote>
  <p>For constructing the sp2H2 and sp2H1 H-atoms, we first decide on the orientation of the conjugated plane or the aromatic ring with respect to the neighboring heavy atoms; the normal vector of the conjugated plane is determined by taking the cross product of two vectors between the heavy atoms. For the sp2H2 H-atoms (illustrated in Fig. 1b), the normal vector of the conjugated plane is the cross product of the unit vectors B→N and B2→B; then the two H-atoms are placed at positions A1 and H, which are within the conjugated plane respected to the B→N vector with the exact bond angle from CHARMM22 force field. For sp2H1 H-atoms, two conformations are possible. The first is to place the H-atom in the peptide plane as illustrated in Fig. 1b, where A1 and B1 represent the alpha carbon atoms. The position of H in this case is decided by using the same method as the one used to determine the position of the sp2H2 H-atoms while holding the trans-conformation.</p>
  </blockquote>
  <p>it refers to Fig.1 of linked freely available article</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-10-15T16:06:15.363Z
Answers.LastActivityDate: 2023-10-15T16:06:15.363Z
Answers.UserDisplayName: pippo1980
Body: <p>I know the algorithm for creating a Hydrogen atom and adding to a residue:</p>
  <pre><code>Point3d  create_hydrogen(Point3d C, Point3d N, Point3d CA, Point3d H) 
  {
    H.set(N);
    H -= C;
    H.norm();
  
    Point3d tmp2(N);
    tmp2 -= CA;
    tmp2.norm();
  
    H += tmp2;
    H.norm(N_H_BOND_LENGTH);
    H += N;
  
    return H;
  }
  
  Point3d  create_hydrogen(Residue prev_residue, Residue residue, Point3d H) 
  {
    return create_hydrogen(prev_residue.find_atom(&quot; C  &quot;), 
                                residue.find_atom(&quot; N  &quot;),
                                residue.find_atom(&quot; CA &quot;), 
                                     H);
  }
  
  Atom create_hydrogen(Residue prev_residue, Residue residue) 
  {
    Atom H = make_shared_atom(0, &quot; H  &quot;, 1);
    create_hydrogen(prev_residue, residue, H);
    residue.add(H);
    return H;
  }
  </code></pre>
  <p>I didn't find this implementation in BioPython.</p>
  <p>Is this already implemented in BioPython? If YES, in which module?</p>
  <p>If NO, what would be the implementation like in BioPython?</p>
  <p>Added image of peptide bond too try to depict the problem at hand, I believe the hydrogen in cause is the one on the peptide bond plane(pale yellow) indentifiedin purple</p>
  <p><a href="https://i.sstatic.net/IUpqe.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/IUpqe.png" alt="enter image description here" /></a></p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T06:57:47.397Z
Comments.Text: Search for Bio.PDB.Structure.Structure object in biopython code
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T09:38:14.497Z
Comments.Text: I would suggest again using already existing tools to add hydrogens. You said your PI had concerns with this, but I must say that it is not a random process like letting kids adding baubles to a Christmas tree: for sp3 heavy atoms bound by two heavy atoms and for terminal heavy atom in sp2 the hydrogens have a defined geometric position. Terminal amines, hydroxyls and methyls will require more advanced calculations, but bar for ST turns these will not contribute to SS.
Comments.UserDisplayName: Matteo Ferla
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T10:08:46.353Z
Comments.Text: @MatteoFerla, I understand. I talked to him, and he doesn't agree. This is his project and his decision is final. :) What can I do?
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T10:09:42.420Z
Comments.Text: @pippo1980, *Search for Bio.PDB.Structure.Structure object in biopython code* --- there is nothing useful in that module.
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T10:59:55.210Z
Comments.Text: Switch to openbabel ? Look how the Arpeggio guys handled the visualization of ligand protein interaction of pdb : https://github.com/harryjubb/arpeggio/blob/master/arpeggio.py
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T11:27:45.387Z
Comments.Text: @pippo1980, How about *Biotite*?
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T11:44:08.403Z
Comments.Text: Or http://thegrantlab.org/bio3d/
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T12:24:09.593Z
Comments.Text: @pippo1980, *Bio3D* is offered in `R`. That's not useful .
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T13:41:13.267Z
Comments.Text: You could refactor it in Python ;-)
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T13:45:49.510Z
Comments.Text: Meaning translate/port
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T13:51:41.953Z
Comments.Text: @pippo1980, That will take months. I don't have that much time to invest right now.
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-01T15:13:31.953Z
Comments.Text: Is it stille related to Calculate Hydrogen bonds between each two chains Biopython : like https://stackoverflow.com/questions/43030557/calculate-hydrogen-bonds-between-each-two-chains-biopython
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-03T19:40:28.047Z
Comments.Text: The code that you posted is a trivial function that adds H atom to N atom in a given distance. Do you want to reproduce just this?
Comments.UserDisplayName: marcin
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-03T20:30:56.650Z
Comments.Text: @marcin, For me it is non-trivial. I am finding to similarity in BioPython and this code.
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-03T20:31:29.477Z
Comments.Text: @marcin, *Do you want to reproduce just this?* --- For now, YES.
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-04T10:56:24.037Z
Comments.Text: First you need to understand how this code work. It's a geometrical exercise. C, N, CA and H are in the same plane. The code seems to assume that angles H-N-C and H-N-CA are equal (which is approximately true). It calculates two vectors, normalizes them, adds them and sets the length of the resulting vector. So you'd need to do the same in Python. I don't know why you want to use BioPython (and what's the point of the whole exercise), but you can implement the same in Python using, for example, Vector from Bio.PDB and then add H atom to Bio.PDB.Residue.
Comments.UserDisplayName: marcin
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-04T11:10:14.147Z
Comments.Text: @marcin, what framework do u recommend other than BioPython?
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-04T11:22:29.667Z
Comments.Text: Hard to tell without knowing more about this project. I see that many people use PyMOL for simple scripting (it's molecular viewer, but supports scripts in Python). If you'd like to add all hydrogen atoms, the most popular program is REDUCE from the Richardson lab.
Comments.UserDisplayName: marcin
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-04T11:30:14.677Z
Comments.Text: @marcin, I need a python-based Bioinformatics protein manipulation framework that is Open source and easy to learn. E.g. SeqAn has great documentation. They also published a book which is very convenient for me.
Comments.UserDisplayName: user366312
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-13T08:38:18.490Z
Comments.Text: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2724740/ HAAD: A Quick Algorithm for Accurate Prediction of Hydrogen Atoms in Protein Structures it compares and cites other alghoritm too, reduce seems to be the older one maybe still the more robust accurate
Comments.UserDisplayName: pippo1980
Tags: <python><protein-structure>
Title: How can I programmatically add a Hydrogen 'Atom' to a 'Residue' object?
ThreadText: # How can I programmatically add a Hydrogen 'Atom' to a 'Residue' object?
  
  (Asked by: username_0 on 2021-10-31)
  
  I know the algorithm for creating a Hydrogen atom and adding to a residue:
  
  ```
  Point3d  create_hydrogen(Point3d C, Point3d N, Point3d CA, Point3d H) 
  {
    H.set(N);
    H -= C;
    H.norm();
  
    Point3d tmp2(N);
    tmp2 -= CA;
    tmp2.norm();
  
    H += tmp2;
    H.norm(N_H_BOND_LENGTH);
    H += N;
  
    return H;
  }
  
  Point3d  create_hydrogen(Residue prev_residue, Residue residue, Point3d H) 
  {
    return create_hydrogen(prev_residue.find_atom(" C  "), 
                                residue.find_atom(" N  "),
                                residue.find_atom(" CA "), 
                                     H);
  }
  
  Atom create_hydrogen(Residue prev_residue, Residue residue) 
  {
    Atom H = make_shared_atom(0, " H  ", 1);
    create_hydrogen(prev_residue, residue, H);
    residue.add(H);
    return H;
  }
  
  ```
  
  I didn't find this implementation in BioPython.
  
  Is this already implemented in BioPython? If YES, in which module?
  
  If NO, what would be the implementation like in BioPython?
  
  Added image of peptide bond too try to depict the problem at hand, I believe the hydrogen in cause is the one on the peptide bond plane(pale yellow) indentifiedin purple
  
  [![enter image description here](https://i.sstatic.net/IUpqe.png)](https://i.sstatic.net/IUpqe.png)
  
  **Tags:** python, protein-structure
  
  > **Comments:**
  > **username_1**: Search for Bio.PDB.Structure.Structure object in biopython code  
  > **username_2**: I would suggest again using already existing tools to add hydrogens. You said your PI had concerns with this, but I must say that it is not a random process like letting kids adding baubles to a Christmas tree: for sp3 heavy atoms bound by two heavy atoms and for terminal heavy atom in sp2 the hydrogens have a defined geometric position. Terminal amines, hydroxyls and methyls will require more advanced calculations, but bar for ST turns these will not contribute to SS.  
  > **username_0**: @username_2, I understand. I talked to him, and he doesn't agree. This is his project and his decision is final. :) What can I do?  
  > **username_0**: @username_1, *Search for Bio.PDB.Structure.Structure object in biopython code* --- there is nothing useful in that module.  
  > **username_1**: Switch to openbabel ? Look how the Arpeggio guys handled the visualization of ligand protein interaction of pdb : https://github.com/harryjubb/arpeggio/blob/master/arpeggio.py  
  > **username_0**: @username_1, How about *Biotite*?  
  > **username_1**: Or http://thegrantlab.org/bio3d/  
  > **username_0**: @username_1, *Bio3D* is offered in `R`. That's not useful .  
  > **username_1**: You could refactor it in Python ;-)  
  > **username_1**: Meaning translate/port  
  > **username_0**: @username_1, That will take months. I don't have that much time to invest right now.  
  > **username_1**: Is it stille related to Calculate Hydrogen bonds between each two chains Biopython : like https://stackoverflow.com/questions/43030557/calculate-hydrogen-bonds-between-each-two-chains-biopython  
  > **username_3**: The code that you posted is a trivial function that adds H atom to N atom in a given distance. Do you want to reproduce just this?  
  > **username_0**: @username_3, For me it is non-trivial. I am finding to similarity in BioPython and this code.  
  > **username_0**: @username_3, *Do you want to reproduce just this?* --- For now, YES.  
  > **username_3**: First you need to understand how this code work. It's a geometrical exercise. C, N, CA and H are in the same plane. The code seems to assume that angles H-N-C and H-N-CA are equal (which is approximately true). It calculates two vectors, normalizes them, adds them and sets the length of the resulting vector. So you'd need to do the same in Python. I don't know why you want to use BioPython (and what's the point of the whole exercise), but you can implement the same in Python using, for example, Vector from Bio.PDB and then add H atom to Bio.PDB.Residue.  
  > **username_0**: @username_3, what framework do u recommend other than BioPython?  
  > **username_3**: Hard to tell without knowing more about this project. I see that many people use PyMOL for simple scripting (it's molecular viewer, but supports scripts in Python). If you'd like to add all hydrogen atoms, the most popular program is REDUCE from the Richardson lab.  
  > **username_0**: @username_3, I need a python-based Bioinformatics protein manipulation framework that is Open source and easy to learn. E.g. SeqAn has great documentation. They also published a book which is very convenient for me.  
  > **username_1**: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2724740/ HAAD: A Quick Algorithm for Accurate Prediction of Hydrogen Atoms in Protein Structures it compares and cites other alghoritm too, reduce seems to be the older one maybe still the more robust accurate  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_1 on 2021-11-01. Score: 2
  
  Not sure but here <https://biopython.org/docs/1.75/api/Bio.PDB.Residue.html> Bio.PDB.Residue module Your can find :
  
  add(self, atom) :
  
  ```
        Add an Atom object.
  
  ```
  
  Checks for adding duplicate atoms, and raises a PDBConstructionException if so.
  
  Can’t try it just now but who knows it may work
  
  Googling around I found out that if you plan to add Hydrogens to PDB files [I believe nowadays PDBx is the repository grade standard] you could do thar following PDB format rules. I am copying them from UCSF Chimera website (<https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens>):
  
  ```
  
  Hydrogen Atoms
  
  In brief, conventions for hydrogen atoms in version 3.0 PDB format are as follows:
  
  Hydrogen atom records follow the records of all other atoms of a particular residue.
  A hydrogen atom name starts with H. The next part of the name is based on the name of the connected nonhydrogen atom. For example, in amino acid residues, H is followed by the remoteness indicator (if any) of the connected atom, followed by the branch indicator (if any) of the connected atom; if more than one hydrogen is connected to the same atom, an additional digit is appended so that each hydrogen atom will have a unique name. Hydrogen atoms in standard nucleotides and amino acids (other than the rarely seen HXT) are named according to the IUPAC recommendations (Markley et al., Pure Appl Chem 70:117 (1998)). Names of hydrogen atoms in HETATM residues are determined in a similar fashion.
  If the name of a hydrogen has four characters, it is left-justified starting in column 13; if it has fewer than four characters, it is left-justified starting in column 14. 
  
  In the following excerpt from entry 1vm3, atom H is attached to atom N. Atom HA is attached to atom CA; the remoteness indicator A is the same for these atoms. Two hydrogen atoms are connected to CB, one is connected to CG, three are connected to CD1, and three are connected to CD2.
  
  ATOM     10  N   LEU A   2       4.595   6.365   3.756  1.00  0.00           N
  ATOM     11  CA  LEU A   2       4.471   5.443   2.633  1.00  0.00           C
  ATOM     12  C   LEU A   2       5.841   5.176   2.015  1.00  0.00           C
  ATOM     13  O   LEU A   2       6.205   4.029   1.755  1.00  0.00           O
  ATOM     14  CB  LEU A   2       3.526   6.037   1.578  1.00  0.00           C
  ATOM     15  CG  LEU A   2       2.790   4.919   0.823  1.00  0.00           C
  ATOM     16  CD1 LEU A   2       3.803   3.916   0.262  1.00  0.00           C
  ATOM     17  CD2 LEU A   2       1.817   4.196   1.769  1.00  0.00           C
  ATOM     18  H   LEU A   2       4.169   7.246   3.704  1.00  0.00           H
  ATOM     19  HA  LEU A   2       4.063   4.514   2.992  1.00  0.00           H
  ATOM     20  HB2 LEU A   2       2.804   6.675   2.065  1.00  0.00           H
  ATOM     21  HB3 LEU A   2       4.099   6.623   0.873  1.00  0.00           H
  ATOM     22  HG  LEU A   2       2.234   5.353   0.004  1.00  0.00           H
  ATOM     23 HD11 LEU A   2       4.648   4.447  -0.148  1.00  0.00           H
  ATOM     24 HD12 LEU A   2       3.334   3.331  -0.516  1.00  0.00           H
  ATOM     25 HD13 LEU A   2       4.137   3.260   1.052  1.00  0.00           H
  ATOM     26 HD21 LEU A   2       0.941   3.892   1.216  1.00  0.00           H
  ATOM     27 HD22 LEU A   2       1.522   4.860   2.568  1.00  0.00           H
  ATOM     28 HD23 LEU A   2       2.296   3.323   2.188  1.00  0.00           H
  
  ```
  
  ```
  
  > **Comments:**
  > **username_0**: As you can see from the algorithm, adding Hydrogen is not that straightforward. :D If it were that simple, I wouldn't have posted this question!!  
  > **username_1**: Are you talking about the geometry of the missing H or about the protonation state of a lateral chain or about polar hyrogens ?  
  > **username_1**: see here for H atom names for Phenilalanine , bond lenght and angle have to be calculated from the atom coords https://www.ebi.ac.uk/pdbe-srv/pdbechem/atom/list/PHE or tabulated here https://www.ebi.ac.uk/pdbe-srv/pdbechem/bond/list/PHE  
  
  ---
  
  ### Answer by username_2 on 2021-11-01. Score: 2
  
  The other answer is correct. But I thought I'd give a few pointers in how one can understand how to do something with a given Python module —teach a man to fish kind of thing...
  
  In a Jupyter Notebook or IPython shell (don't), you can do:
  
  * `help(function)` print out the docstring
  * `dir(object)` shows all the public and magic attributes and methods —note a class will not show the instance methods
  * `type(object)` or `object.__class__` gives you what the class is, but `object.__class__.__mro__` tells you what it inherits.
  * in the standard library module `inspect` is the function `inspect.getsource(function)`, which for a notebook is messy due to newlines so needs to be `print_code = lambda fun: print(inspect.getsource(fun))` and `print_code(fun)`.
  
  In your case you read a file...
  
  ```
  from Bio import PDB
  parser = PDB.PDBParser()
  structure : PDB.Structure.Structure = parser.get_structure("tryp-cage", "1L2Y.first.pdb")
  
  ```
  
  So `parser.get_structure` "magically" makes a structure, eh? This means it can add stuff! As mentioned `print_code(parser.get_structure)` will reveal its secrets, or within it call some `self._private_fun()`, which can be inspected (`self` = `parser`) all the way down the rabbit hole. Until you find a nice constructor that declares nicely an atom instance or whatever you ever fancy searching for —that is the great thing of open source code.
  
  In the specific case of `PDB` module, everything inherits `Entity`, which has the attribute `child_list`, which is what the iterators (`__iter__` is the topmost) read without making a copy. The end of the above rabbit hole may just add a child `Entity`-inheriting class instance to this list, but it also do other things that is worth keeping an eye out.
  
  In the case of `residue.add(atom)`, creating the `Bio.PDB.Atom.Atom` instance is slightly problematic as there's no typehinting or `:type xxx:` declarations in the docstring. But doing `type(atom.bfactor)` etc. on the attributes of an existing atom will reveal how to do it. One the serial number is a PDB thing, which is sequential for the model.
  
  ```
  last_serial = list(structure.get_residues())[-1].child_list[-1].serial_number + 1
  residue.add(PDB.Atom.Atom(name=' H  ', coord=np.array([1,1,1]), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number=last_serial+1,element='H'))
  
  ```
  
  The coordinates will not be `[1,1,1]`, but the point that you have calculate. Say for `H`, this would be coplanar with `C`, `CA` and `N`, that is 1 and smidge Å away from `N` at an angle not quite in the middle as found in amide bonds.
  
  > **Comments:**
  > **username_0**: If the other answer is correct, what do you think about the algorithm I posted? This algorithm is taken from the c++ source code I have in my hand.  
  > **username_0**: As you would realize, I can write Python and C++ program, but I am just confused about hydrogens. :)  
  > **username_1**: Biopython documentation is not that rich especially for ones like me that use google as primary/unique docs search  
  > **username_2**: I have expanded the answer with an example, showing that what I wrote is how one finds the correct solution without documentation. In terms of how to do the maths, you have to do that in numpy as its simple maths —I do not know if the C++ code is 100% correct as it is projecting on the line between C-CA midpoint and N, which I am not sure is the correct angle. Totally easy to verify though  
  > **username_1**: I would be curious about how Bio.PDB.PDBIO module (https://biopython.org/docs/1.74/api/Bio.PDB.PDBIO.html) would save such pdb with added H atoms, like : https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens or something else ?  
  > **username_1**: @username_2 Hi, sorry to bother, would you be so kind to point me to a PDB (x-ray) into https://www.rcsb.org/ , that contains H atoms. Don't know how to search for them, or some other PDB compliant file with Hs, to play with biopython I am curious about it, I know I should learn PDBx by now but is far more complicated to me  
  > **username_1**: dont worry got it from http://biophysics.cs.vt.edu/ even if it marks all H as atom 0, without renumbering everything from scratch  
  > **username_2**: Neutron diffraction structures have protons experimentally. Say PDB:5MOP —A beauty as it has a hydronium or two. For small molecules x-ray crystallography can detect protons but the database for that is the [CCDC](https://ccdc.cam.ac.uk/)  
  > **username_1**: Do you have code for the correct naming of H atoms ? Does PDBx follow same rules ?  
  > **username_2**: All formats that support atom names use the same atom names and _should_ match the definition in the PBD. e.g. https://files.rcsb.org/ligands/view/ALA.cif Ligands are not always numbered by canonical index order, but amino acids follows a precise lettering scheme (greek) for heavy atoms —the annoying Val <-> Ile superposition issue is common to all tools. However, for protons there are two schemes: ` HB1` vs. `1HB `, the former is like the heavy atoms while the latter is more sensible as the number prefix is compatible with heavy atoms who have a number suffix, cf. ARG with `HH11` vs. `1HH1`  
  > **username_1**: —the annoying Val <-> Ile superposition issue is common to all tools ... what is that ?  
  > **username_2**: The CB atom is not chiral in valine, but the atom names make it "name chiral" (not a real word, but just like isotopic chirality): looking at valine with the HB atom forewards, you have CG2 on the side of the N (left) and CG1 on the C (right) in a predictable counterclockwise fashion. Isoleucine is chiral on the CB atom and the CG1 is on the left as it is more complex due to CD1. As a result they do not superpose well by atom names, but do if CG2 is mapped to CG1 and CG1: CG2.  
  > **username_1**: @username_2  AttributeError: 'Structure' object has no attribute 'get_residue' --> last_serial = list(structure.get_atoms())[-1].serial_number ?  
  
  ---
  
  ### Answer by username_1 on 2023-10-15. Score: 1
  
  OK my attempt, input pdb created in **pymol** with `fab AVGLWPT , AVGLWPT , hydro=0` and saved as `avgwpt.pdb` :
  
  ```
  ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N  
  ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C  
  ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C  
  ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O  
  ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C  
  ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N  
  ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C  
  ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C  
  ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O  
  ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C  
  ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C  
  ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C  
  ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N  
  ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C  
  ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C  
  ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O  
  ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N  
  ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C  
  ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C  
  ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O  
  ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C  
  ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C  
  ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C  
  ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C  
  ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N  
  ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C  
  ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C  
  ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O  
  ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C  
  ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C  
  ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C  
  ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C  
  ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C  
  ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C  
  ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N  
  ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C  
  ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C  
  ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C  
  ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N  
  ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C  
  ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C  
  ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O  
  ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C  
  ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C  
  ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C  
  ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N  
  ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C  
  ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C  
  ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O  
  ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C  
  ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C  
  ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O  
  TER   
  END
  
  ```
  
  pic :
  [![enter image description here](https://i.sstatic.net/NZNZY.png)](https://i.sstatic.net/NZNZY.png)
  
  code :
  
  ```
  #!/usr/bin/env python3
  # -*- coding: utf-8 -*-
  """
  molecule create in pymol with fab AVGLWPT , AVGLWPT , hydro=0
  
  """
  
  import Bio
  print('########## Biopython  VERSION ##########################################')
  print('\n\n BIOPYTHON_VERSION : ', Bio.__version__)
  print('###################################################################')
  
  from Bio.PDB import (
                      PDBParser,
                      Atom,
                      PDBIO
                      )
  
  structure_code = 'avgwpt'
  
  structure_name = structure_code
  
  parser_pdb= PDBParser(QUIET=False)    ### shows warnings
  
  structure_whole = parser_pdb.get_structure(structure_name, structure_name+".pdb")  
  
  last_serial = list(structure_whole.get_atoms())[-1].serial_number
  
  print(" last_serial : ",last_serial)
  
  for model in structure_whole:
      
      for chain in model:
          
          for resi in chain:
          
              print('1 : ', resi)
              
              print('2 : ',resi.id)
              
              try :
              
                  print('3 : ',chain[resi.id[1]-1])
                  
                  doable = True
                  
              except :
                  
                  print('3 : ','residue -1 missing !!!!!!!!!!!!!!!!!!!!!!1')
                  
                  doable = False
                  
              
              print('4 : ',[i for i in resi.get_atoms()])
              
              if doable:
                  
                  print('5 : ','ok')
              
              else:
                  
                  print('5 : ','NO')
                  
              
              if doable:
                  
                  
                  print(chain[resi.id[1]-1]['C'])
                  
                  print(chain[resi.id[1]-1]['C'].get_vector(), type(chain[resi.id[1]-1]['C'].get_vector()))
                  
                  print(chain[resi.id[1]]['N'].get_vector())
                  
                  c = chain[resi.id[1]-1]['C'].get_vector()
                  
                  n = chain[resi.id[1]]['N'].get_vector()
                  
                  ca = chain[resi.id[1]]['CA'].get_vector()
                  
                  nc = n-c
                  
                  nca = n-ca
                  
                  print('vector BN : ', nc)
                  
                  nc.normalize()
                  
                  print('vector BN normalize : ', nc)
                  
                  print('vector B2B : ', nca)
                  
                  nca.normalize()
                  
                  print('vector B2B normalize : ', nca)
                  
                  print('sum of nc + nca : ' , nsum := nc + nca)
                  
                  nsum.normalize()
                  
                  nsum**1.023
                  
                  nsum += n
                  
                  print(nsum)
                  
                  print(nsum.get_array(), type(nsum.get_array()))
                  
                  resi.add(Atom.Atom(name=' H  ', coord = nsum.get_array(), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number = last_serial ,element='H'))
              
                  last_serial += 1
  
  io=PDBIO() 
  
  io.set_structure(structure_whole ) 
  
  io.save(structure_name+"_H_added.pdb" , preserve_atom_numbering = True) 
  
  ```
  
  output , saved as `avgwpt_H_added.pdb`:
  
  ```
  ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N  
  ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C  
  ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C  
  ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O  
  ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C  
  ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N  
  ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C  
  ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C  
  ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O  
  ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C  
  ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C  
  ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C  
  ATOM     52  H   VAL     2      70.368  54.810  15.182  1.00  0.00           H  
  ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N  
  ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C  
  ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C  
  ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O  
  ATOM     53  H   GLY     3      74.741  54.249  15.182  1.00  0.00           H  
  ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N  
  ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C  
  ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C  
  ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O  
  ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C  
  ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C  
  ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C  
  ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C  
  ATOM     54  H   LEU     4      76.733  58.181  15.224  1.00  0.00           H  
  ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N  
  ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C  
  ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C  
  ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O  
  ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C  
  ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C  
  ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C  
  ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C  
  ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C  
  ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C  
  ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N  
  ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C  
  ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C  
  ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C  
  ATOM     55  H   TRP     5      81.105  57.620  15.224  1.00  0.00           H  
  ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N  
  ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C  
  ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C  
  ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O  
  ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C  
  ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C  
  ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C  
  ATOM     56  H   PRO     6      83.072  61.539  15.267  1.00  0.00           H  
  ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N  
  ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C  
  ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C  
  ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O  
  ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C  
  ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C  
  ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O  
  ATOM     57  H   THR     7      87.409  60.156  16.044  1.00  0.00           H  
  TER      58      THR     7                                                       
  END   
  
  ```
  
  pic : [![enter image description here](https://i.sstatic.net/WmGL1.png)](https://i.sstatic.net/WmGL1.png)
  
  Not sure I got the code above (in question) right, can't read it , had to use comments as trace.
  
  Spent a while banging my head with : [HAAD: A Quick Algorithm for Accurate Prediction of Hydrogen Atoms in Protein Structures](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2724740/)
  
  and trying to figure out their alghoritm, [code (<https://seq2fun.dcmb.med.umich.edu//HAAD/haad.f90>) is in FORTRAN90] from the paper , but no luck, if anybody wants to help eplaining this bit:
  
  > For constructing the sp2H2 and sp2H1 H-atoms, we first decide on the orientation of the conjugated plane or the aromatic ring with respect to the neighboring heavy atoms; the normal vector of the conjugated plane is determined by taking the cross product of two vectors between the heavy atoms. For the sp2H2 H-atoms (illustrated in Fig. 1b), the normal vector of the conjugated plane is the cross product of the unit vectors B→N and B2→B; then the two H-atoms are placed at positions A1 and H, which are within the conjugated plane respected to the B→N vector with the exact bond angle from CHARMM22 force field. For sp2H1 H-atoms, two conformations are possible. The first is to place the H-atom in the peptide plane as illustrated in Fig. 1b, where A1 and B1 represent the alpha carbon atoms. The position of H in this case is decided by using the same method as the one used to determine the position of the sp2H2 H-atoms while holding the trans-conformation.
  
  it refers to Fig.1 of linked freely available article
  
  
  ---
  
  