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!