MDAnalysis how to index the charge?

Question

I have been trying to figure out this process that is very dependent on the charge. However, I am finding that, for waters, I am getting a value not equal to 0 for the sum of a water molecule charge. To investigate this, I wanted to compare the value of the charges being used by MDAnalysis and my topology file and see where the issue could be.

In my initial code I am able to call for the charge by using:

    for ts in u.trajectory[ini_frames:]:
    count += 1
    membrane = u.select_atoms('segid PROA or segid PROB or segid MEMB')
    COM = membrane.atoms.center_of_mass()[2]
    q_prof = CLAs.atoms.charges * (CLAs.positions[:,2] + (Lz/2 - COM))/Lz
    Q_instant = np.sum(q_prof)
    Q_sum += Q_instant
Q_av = Q_sum / n_frames

This has worked for me where

CLAs = u.select_atoms("segid HETA or segid PROA or segid PROB or segid MEMB or segid IONS")

So to look at where this discrepancy could come from I attempted:

def Q_code(dcd, topo):
Lz = u.dimensions[2]
Q_sum = 0
count = 0
CLAs = u.select_atoms('segid TIP3')
ini_frames = -20
n_frames = len(u.trajectory[ini_frames:])
for ts in u.trajectory[ini_frames:]:
    for i in CLAs:
        with open('Q_atom_temp.csv', 'a') as f:
            new_line = [s, i, i.charges, CV1_dict[s], CV2_dict[s]]
            writes = csv.writer(f)
            writes.writerow(new_line)
            f.close()
fields = ['Window', 'Atom', 'Charge', 'CV1', 'CV2'] with
 open('Q_atom_temp.csv', 'a') as f:
      writer = csv.writer(f)
      writer.writerow(fields) f.close()

But I am getting the error: So how would I go about making a file that will show me which atom is corresponding to which charge?

Are the indexes different? Thanks for any help!

Answer 1

It turns out this is indexed in the same manner. My problem was that I was applying things to each individual atom instead of its molecule. When I did that I got an answer that made sense.

So for anyone else reading this, you can iterate through your positions and charges using the same indexing.

for i in range(len(CLAs.positions[:,2])):
        #print(CLAs.atoms.charges[i], CLAs.atoms.names[i], CLAs.positions[i,2])
        total_charge += CLAs.atoms.charges[i]
        count += 1
        position += CLAs.positions[i,2]
        if count == 5:
            count = 0
            av_position = (position / 5)
            q_prof += total_charge * (av_position + (Lz/2 - COM))/Lz
            position = 0
            total_charge = 0

This code is looking at the individual charges, mapping it to a position of the same ion and then after the water molecule is complete, taking the charge of the molecule and multiplying it by its average position relative to the box.