waterdynamics. AngularDistribution

[AQ27na]

I am using the waterdynamics.AngularDistribution function. And I want to study the water molecules in certain area (8.55 <= y <= 11.25), here is the python content. But there is an error:

 line 623, in _getCosTheta
    OHVector0 = H1t0 - Ot0
                ~~~~~^~~~~
ValueError: operands could not be broadcast together with shapes (93,3) (94,3)

Below is my python content. how can I fix this problem? And I have tried if I only select water without the y range it can work. So I guess, if it is that water molecules are not complete in my chosen and how can I do to make sure the water molecules are complete? Thanks a lot.

frame_begin = 5000
frame_end = 5001

#orientation analysis
AD_analysis = AD(
    u,
    'prop y >= 8.55 and prop y <= 11.25 and resname water',
    bins=40,
    axis='y'
)
AD_analysis.run(start=frame_begin, stop=frame_end)
#The first two columns are P(cos(theta)) vs cos(theta) for OH vector,
#the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector,
#and thirds two columns are P(cos(theta)) vs cos(theta) for dipole vector
with open ('c:/Users/AnqiDuan/Desktop/orientation.txt', 'w') as file:
    for bin in range(bins):
        print(
            "{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(
                AD_analysisOH=AD_analysis.graph[0][bin],
                AD_analysisHH=AD_analysis.graph[1][bin],
                AD_analysisDip=AD_analysis.graph[2][bin]
            ), file=file
        )

[orbeckst]

Did you manually check the selection 'prop y >= 8.55 and prop y <= 11.25 and resname water' and does it result in broken water molecules?

If this is the case, try using the same keyword to select complete residues

'same resid as (prop y >= 8.55 and prop y <= 11.25 and resname water)'

Double-check manually that these selections give you what you want, i.e. if your water model contains 3 atoms, check that you have 3 atoms for each resid

w = u.select_atoms('same resid as (prop y >= 8.55 and prop y <= 11.25 and resname water)')
w_groups = w.atoms.groupby("resids")     # print w_groups to understand this dict!
print([resid for resid,ag in w_groups.items() if ag.n_atoms != 3])

The last line will print a list of resids for which you did NOT have all three atoms in the w group: this list should be empty [] if everything worked. Try it out with your orginal selection and the one with same.