for i_o in range(1, 4):
                id_o = 'OP' + str(i_o)
                id_h = 'H' + id_o
                if id_o in atoms and id_h in atoms:
                    atom_o = mol.GetAtom(atoms[id_o]['position'])