Getting the charge of a single atom, per loop in MD Analysis

216 views Asked by At

I have been trying to use the partial charge of one particular ion to go through a calculation within mdanalysis. I have tried(This is just a snippet from the code that I know is throwing the error):

Cl = u.select_atoms('resname CLA and prop z <= 79.14') 
Lz = 79.14                          #Determined from system set-up
Q_sum = 0
COM = 38.42979431152344             #Determined from VMD
file_object1 = open(fors, 'a')
print(dcd, file = file_object1)
for ts in u.trajectory[200:]:
    frame = u.trajectory.frame
    time = u.trajectory.time
    for coord in Cl.positions:
        q= Cl.total_charge(Cl.position[coord][2])
        coords = coord - (Lz/COM)
        q_prof = q * (coords + (Lz / 2)) / Lz
        Q_sum = Q_sum + q_prof
        print(q)

But I keep getting an error associated with this.

How would I go about selecting this particular atom as it goes through the loop to get the charge of it in MD Analysis? Before I was setting q to equal a constant and the code ran fine so I know it is only this line that is throwing the error:

q = Cl.total_charge(Cl.position[coord][2])

Thanks for the help!

1

There are 1 answers

0
Confused_scientist On

I figured it out with:

def Q_code(dcd, topo):
Lz = u.dimensions[2]
Q_sum = 0
count = 0
CLAs = u.select_atoms('segid IONS or segid PROA or segid PROB or segid MEMB')
ini_frames = -200
n_frames = len(u.trajectory[ini_frames:])
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
with open('Q_av.txt', 'a') as f:
    print('The Q_av for {}   is   {}'.format(s, Q_av), file = f)
return Q_av