Dear experts,
I am adding activity to my Gaussian polymer, monomers connected with springs with infinite extensibility, along the bond vector of ith monomer.
I am contemplating good ways to do so, and I think of doing it using atom-type variable. I have to calculate on per-atom basis the dx and dy from n and n+1th bead for nth bead. And then, using trigonometry, I will add force in the direction of the bond vector. But I am not able to calculate the bond vector in first place. Any suggestions will be very helpful for me.
I don’t think that there is much of a chance for that. This sounds more like you need to implement a custom bond style.
I think I am having a very bad habit of writing strange and computationally inefficient ways to employ required simulation demand. I made a workaround to propel the polymer along its contour, and I am pasting it below. It was more to address the question, whether I can do that without writing custom styles.
#tangential force
label addF
variable nbead loop 1 199 #(200 bead polymer)
variable nextb equal ${nbead}+1
variable dx${nbead} atom (x[${nextb}]-x[${nbead}])*(id==${nbead})
variable dy${nbead} atom (y[${nextb}]-y[${nbead}])*(id==${nbead})
variable Fx${nbead} atom ${fp}*v_dx${nbead}/((!(sqrt(v_dx${nbead}*v_dx${nbead}+v_dy${nbead}*v_dy${nbead}))*1.0)+sqrt(v_dx${nbead}*v_dx${nbead}+v_dy${nbead}*v_dy${nbead}))
variable Fy${nbead} atom ${fp}*v_dy${nbead}/((!(sqrt(v_dx${nbead}*v_dx${nbead}+v_dy${nbead}*v_dy${nbead}))*1.0)+sqrt(v_dx${nbead}*v_dx${nbead}+v_dy${nbead}*v_dy${nbead}))
fix addF${nbead} all addforce v_Fx${nbead} v_Fy${nbead} 0.0
next nbead
jump SELF addF
Try this:
compute a property/atom xu yu id # relevant info
#endfinder
compute molchunk all chunk/atom molecule ids once
compute molmaxid all reduce/chunk molchunk max c_a[3]
compute s_maxid all chunk/spread/atom molchunk c_molmaxid
variable is_end atom c_s_maxid==id
#parity chunking id: 1 2 3 4 5 6 ...
variable parity atom id%2 # 1 0 1 0 1 0 ...
variable oddpair atom id+v_parity # 2 2 4 4 6 6 ...
variable evnpair atom id-v_parity # 0 2 2 4 4 6 ...
compute oddchunk all chunk/atom v_oddpair ids once
compute evnchunk all chunk/atom v_evnpair ids once
# process unwrapped displacements
compute oddsum all reduce/chunk oddchunk sum c_a[1] c_a[2]
compute evnsum all reduce/chunk evnchunk sum c_a[1] c_a[2]
compute s_oddsum all chunk/spread/atom oddsum[*]
compute s_evnsum all chunk/spread/atom evnsum[*]
variable dx atom (v_parity*s_oddsum[1]+(1-v_parity)*s_evnsum[1])-xu
variable dy atom (v_parity*s_oddsum[2]+(1-v_parity)*s_evnsum[2])-yu
variable dr atom sqrt(v_dx*v_dx+v_dy*v_dy+1e-6*v_is_end) # avoid divBy0
1 Like
It is a very bad practice to post the same question in different places and, even worse, in an unrelated thread.
Besides that, your script ends with:
print "All done"
So, everything is going to be fine.