尝试计算两个粒子之间的 RDF - 收到错误消息

问：

“我正在尝试计算A和B粒子（AA、AB和BB）之间的rdf。我收到一条错误消息，说 “文件 ~\anaconda3\Lib\site-packages\matplotlib\axes_base.py：504 在_plot_args raise ValueError（f“x 和 y 必须具有相同的第一维，但”

ValueError： x 和 y 必须具有相同的第一维，但具有形状 （1，） 和 （100，）”

谁能告诉我我错过了什么

这是我的代码：

import numpy as np
from homework_helper_functions import *
import matplotlib.pyplot as plt


 
def calc_my_rdf(pos, atype, len_box, r_max=5, dr=0.05, rdf_type='AA'):
   """
"Parameters
----------
"pos : Nx3 array
   "An array of the (unwrapped) xyz positions of N simulation beads in a polymer.
atype: Nx1 vector
   A vector of the nanoparticle(np) types for all N nps in the simulation. The order matches the order of the pos array.
   Type 1 = A and Type 2 = B
len_box: float
   The length of one side of the simuulation box (this is loaded in with the read_dump_np function)
rmax: float
   A number giving the maximum r value we calculate our rdf up to (keeping this at 5 will give the same results as LAMMPS from HW2a)
dr: float
   The width of our bins (keeping this at 0.05 will give the same results as LAMMPS from HW2a)
rdf_type: string
   The type of rdf you want to calculate ('AA', 'AB', 'BB', 'all')" 




Returns
-------
r_vals : M-length vector of the r-values we use to compute the rdf
   (M = 100 if you use the rmax = 5, dr = 0.05)
   (each r value is the middle of its bin in bin_radius)
rdf : M-length vector of the rdf at each r-value
"""
   # Initialize some values/arrays you might need
   nbins = int(r_max / dr)
   counts = [0] * nbins
   rdf = [0] * nbins




   # make array of bin radii from rmax and drs
   radii = np.linspace(0, r_max + dr, nbins + 1)
   r_vals = (nbins) / len(counts)
   volumes = (4 / 3) * 3.14 * radii ** 3
   shell_vols = np.zeros(len(volumes) - 1)
   for i in range(len(volumes) - 1):
       shell_vols[i] = volumes[i + 1] - volumes[i]
       
   # masking
   # initiate counts of each bin "done"




   # Select particles from pos that match the rdf_type you want to calculate
   mask = atype == 1
   pos_A = pos[mask]
   pos_B = pos[~mask]
   if rdf_type == "AA":
       pos0 = pos_A
       pos1 = pos_A
   elif rdf_type == "BB":
       pos0 = pos_B
       pos1 = pos_B
   elif rdf_type == "AB":
       pos0 = pos_A
       pos1 = pos_B
   elif rdf_type == "all":
       pos0 = pos
       pos1 = pos
   # you can use an if statement (based on 'rdf_type') to do this, or else you might have to modify the code slighlty. done
   num_nps_type1 = len(pos0)
   num_nps_type2 = len(pos1)
   for i in range(num_nps_type1):  # len - length of array = len(pos0) len(pos1),,, p1 = poso[i] p2 = pos1[j]
       for j in range(num_nps_type2):
           if j < i:
               continue
           p1 = pos0[i]
           p2 = pos1[j]
           r = pbc_distance(p1, p2, len_box)
           # determine which bin this r goes in, and update counts
           if r > r_max:
               continue


           bin_no = int(r / dr)
           counts[bin_no] += 1 if j == i else 2


   # Normalize the bin counts and store rdf
   # total volume = box volume p(r) counts/volume, p(b) = total counts/total volume
   p_r = np.array(counts) / np.array(shell_vols)
   p_b = num_nps_type2 * num_nps_type1 / len_box ** 3
   rdf = p_r / p_b
   return r_vals, rdf



# change this to where your files are stored
directory = r'C:\Users\Kelsey\.spyder-py3\HW2'




files = load_all_trajectories(directory)

all_rdfs = []




# You can change which files to loop through (i.e. for testing your code on a single file)
# by changing "files[:1]:" to "files:" in the following line, so that the script only
# loops through every file in "files" rather than only the first.
for file in files:
   print(file)
   # load the data from the LAMMPS trajectory
   box_lens, timesteps, xpos, ypos, zpos, atype, natoms = read_dump_np(file)




   # Specify the timestep(s) we are interested in-- This is a list of the indices of the relevant timesteps in the "timesteps" array
   ts_for_analysis = [-1]  # Change to range(len(timesteps)) if only interested in all timesteps
   rdfs_file = []
   for i, ts in enumerate(ts_for_analysis):
       print(ts)
       # get the data for the current timestep
       len_box = box_lens[i]
       x = xpos[ts, :][:, np.newaxis]
       y = ypos[ts, :][:, np.newaxis]
       z = zpos[ts, :][:, np.newaxis]
       pos = np.hstack((x, y, z))




       rdf_ts = []
       r_vals, rdf_aa = calc_my_rdf(pos, atype, len_box, rdf_type='AA')  # calculate the rdfs using the function we defined above
       r_vals, rdf_ab = calc_my_rdf(pos, atype, len_box, rdf_type='AB')
       r_vals, rdf_bb = calc_my_rdf(pos, atype, len_box, rdf_type='BB')
       r_vals, rdf_all = calc_my_rdf(pos, atype, len_box, rdf_type='all')
       rdf_all_types = np.hstack((np.array(rdf_aa)[:, np.newaxis], np.array(rdf_ab)[:, np.newaxis], np.array(rdf_bb)[:, np.newaxis],np.array(rdf_all)[:, np.newaxis]))
       rdf_ts.append(rdf_all_ty, pes)
       rdfs_file.append(rdf_ts)
       all_rdfs.append(rdfs_file)


plt.plot(r_vals, rdf_all)
plt.xlabel('Distance (r)')
plt.ylabel('RDF')
plt.title(f'RDF for {rdf_type} interactions')
plt.show()

            # Save the RDF for the current rdf_type
rdf_filename = f'{file}_{rdf_type}_rdf.npy'
np.save(rdf_filename, rdf)

rdfs_file.append(rdf)
all_rdfs.append(rdfs_file)

all_rdfs = np.array(all_rdfs)
np.save(r"C:\Users\Kelsey\.spyder-py3\HW2\all_rdfs.npy", all_rdfs)

我的代码在上面列出

`