Element indices messed up after switching from Edit to Object Mode
No reason an edit mode bmesh cannot be used with edit ops.


  I prefer NOT to use the BMesh module, since the bpy operators suit my
  needs much better.


They are not mutually exclusive.  

In answer to


  1) For every vertex in object 1a) Select this vertex (v.select = True)
  1b) Select neighbor vertices (bpy.ops.mesh.select_more()) 1c) Store
  normal angles 2) Compare normal angles within selection... etc. For
  being able to perform a bpy operator (select_more() in this case) I
  need to toggle to Edit Mode. While (as you said) for manipulating mesh
  I need to be in Object Mode. So I constantly need to toggle, right?


As a simple example, cube in edit mode each vert selected separately and default select more (one face away) called. For the cube this will result in only one vert being unselected.  Note: Select more can also be called with bpy.ops.mesh.select_more(use_face_step=False) to only select via edge connection.

import bpy
import bmesh

context = bpy.context

ob = context.edit_object
me = ob.data

bm = bmesh.from_edit_mesh(me)
print("-" * 33)
for v in bm.verts:
    bpy.ops.mesh.select_all(action='DESELECT')
    v.select_set(True)
    print("Before", v.index, [v.select for v in bm.verts])
    bpy.ops.mesh.select_more()
    print("After", v.index, [v.select for v in bm.verts])

bmesh.update_edit_mesh(me)


Produces output.

---------------------------------
Before 0 [True, False, False, False, False, False, False, False]
After 0 [True, True, True, True, True, True, False, True]
Before 1 [False, True, False, False, False, False, False, False]
After 1 [True, True, True, True, True, True, True, False]
Before 2 [False, False, True, False, False, False, False, False]
After 2 [True, True, True, True, False, True, True, True]
Before 3 [False, False, False, True, False, False, False, False]
After 3 [True, True, True, True, True, False, True, True]
Before 4 [False, False, False, False, True, False, False, False]
After 4 [True, True, False, True, True, True, True, True]
Before 5 [False, False, False, False, False, True, False, False]
After 5 [True, True, True, False, True, True, True, True]
Before 6 [False, False, False, False, False, False, True, False]
After 6 [False, True, True, True, True, True, True, True]
Before 7 [False, False, False, False, False, False, False, True]
After 7 [True, False, True, True, True, True, True, True]


Strongly recommend simply using bmesh for this kind of thing. The linked geometry of an element is readily available.

import bpy
import bmesh

context = bpy.context

ob = context.edit_object
me = ob.data

bm = bmesh.from_edit_mesh(me)
print("-" * 33)
for v in bm.verts:
    # set of all verts surounding based on linked faces
    face_neighbours = {v for f in v.link_faces 
            for v in f.verts} - {v}  
    # list of edge neighbours        
    edge_neighbours = [e.other_vert(v) for e in v.link_edges]
    print(v.index, [v.index for v in face_neighbours])
bmesh.update_edit_mesh(me)


Producing output

---------------------------------
0 [1, 2, 3, 4, 5, 7]
1 [0, 2, 3, 4, 5, 6]
2 [0, 1, 3, 5, 6, 7]
3 [0, 1, 2, 4, 6, 7]
4 [0, 1, 3, 5, 6, 7]
5 [0, 1, 2, 4, 6, 7]
6 [1, 2, 3, 4, 5, 7]
7 [0, 2, 3, 4, 5, 6]


An example to get the angle between vert and edge connected verts' normals

[v.normal.angle(n.normal) for n in edge_neighbours]


Hopefully the examples above might move you a little towards bmesh.  Consider the first example above used on a million vertex mesh, would require two million operator calls, as would the question code. 

Python performance with Blender operators

Second example requires no operators, nor does it change the selection state. Actually There is no need to update back to the mesh, or use edit mode to obtain all the info needed for this example 

Note worth mentioning for this example could be done in object mode directly from mesh data. Quick python console code example, look for faces that contain the vert index.

&gt;&gt;&gt; me = C.object.data
&gt;&gt;&gt; v = me.vertices[0]
&gt;&gt;&gt; v.index
0

&gt;&gt;&gt; set(sum([list(p.vertices) for p in me.polygons if v.index in p.vertices], [])) - {v.index}
{1, 2, 3, 4, 5, 7}


as in both examples above the only vert not in any face containing vert 0 is vert 6.