def get_center_of_mass(trajectory, atom_selection='all'):
return md.compute_center_of_mass(trajectory.atom_slice(atom_selection))
def angle_between_vectors(v1, v2):
v1_u = v1 / np.linalg.norm(v1)
v2_u = v2 / np.linalg.norm(v2)
return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
def base_opening_angle(trajectory, angI, angII, angIII, get_center_of_mass, angle_between_vectors):
print(trajectory.n_frames)
comI = get_center_of_mass(trajectory, angI)
comII = get_center_of_mass(trajectory, angII)
comIII = get_center_of_mass(trajectory, angIII)
theta = list()
for frame in range(trajectory.n_frames):
vec1 = np.subtract(comI[frame], comII[frame])
vec2 = np.subtract(comIII[frame], comII[frame])
theta.append(angle_between_vectors(vec1,vec2))
return np.asarray(theta)
theta = paths.MDTrajFunctionCV(
name='theta',
f=base_opening_angle,
topology=engine.topology,
angI=angI_atoms,
angII=angII_atoms,
angIII=angIII_atoms,
get_center_of_mass=get_center_of_mass,
angle_between_vectors=angle_between_vectors)
for t in tqdm(trajectories):
print(t, len(t), t.to_mdtraj().n_frames)
print(theta(t))
# comI = get_center_of_mass(t.to_mdtraj(), angI_atoms)
# comII = get_center_of_mass(t.to_mdtraj(), angII_atoms)
# comIII = get_center_of_mass(t.to_mdtraj(), angIII_atoms)
# theta = list()
# for frame in range(t.to_mdtraj().n_frames):
# vec1 = np.subtract(comI[frame], comII[frame])
# vec2 = np.subtract(comIII[frame], comII[frame])
# theta.append(angle_between_vectors(vec1,vec2))
# if not len(t) == np.asanyarray(theta).shape[0]:
# print('Uh oh')
Trajectory[20] 20 20
20
[0.47252255 0.48211991 0.48338283 0.45771182 0.4409368 0.43005326
0.40664393 0.41614949 0.40332052 0.40805739 0.40741695 0.42707964
0.43450178 0.42754318 0.40565154 0.31015819 0.35187514 0.34593174
0.3893361 0.38687415]
Trajectory[19] 19 19
19
[0.47252255 0.48211991 0.48338283 0.45771182 0.4409368 0.43005326
0.40664393 0.41614949 0.40332052 0.4160452 0.52950192 0.4253559
0.4181968 0.41631153 0.43509364 0.40095428 0.43764869 0.32323956
0.29988872]
Trajectory[18] 18 18
18
[0.47252255 0.48211991 0.48338283 0.45771182 0.4409368 0.43005326
0.40664393 0.41614949 0.40332052 0.4160452 0.52950192 0.48810919
0.40606379 0.43931374 0.40811438 0.4223151 0.40136018 0.43337016]
Trajectory[15] 15 15
1
[0.19723178 0.26085034 0.14685368] [1.25004246 3.75457833 0.17888462]
python version 3.9.16openpathsampling version 1.5.2mdtraj version 1.9.7I have a collection of accepted TPS trajectories and want add a CV for analysis followig:
https://mattermodeling.stackexchange.com/questions/6583/how-is-the-atom-selection-going-wrong-when-creating-the-collective-variable.
For some of them, the number of frames is incorrectly given as 1 inside MDTrajFunctionCV! (Which, of course, results in a
TypeError: iteration over a 0-d array).The number of frames in the OPS trajectory,
len(trajectory), andtrajectory.to_mdtraj().n_framesreturn the same correct length. My functions on their own work correctly; the issue seems to be somewhere inMDTrajFunctionCV.