c force_electro_circular.f c force_electro_circular(n,gc,ntime) c c gc(maxn) using to store the force vector c add the force_electrostatic interaction in this subroutine c saving at 1~n beads ii=3*i+j-3 c parameters used: c interval(>=1), rmax2,nc,econst,debye,f_art,ddna subroutine force_electro_circular(n,gc,ntime) implicit double precision(a-h, o-z) integer n parameter(maxa=500,maxn=maxa*3,ncmax=5) double precision fe12(3),r_re(3) double precision gc(maxn),art(3) double precision xx(ncmax*maxn) common/matrix/x(maxa,3),bangle(maxa),rz(maxa) common/matrix2/ r(maxa,maxa) common/debyehu/econst,debye,f_art,ddna,rmax,rmax2,rmin,nc,nc3 common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma, +amass,am(maxn),coef,dtime2 c inteval selection depends on interval=interl/r0*0.5 nc33=nc3+3 c update xx(ncmax*maxn) xx(i,k,j) i=1:n k=1:nc j=1:3 ii=0 do i=1,n i2=i+1 if (i.eq.n) then i2=1 endif do k=0,nc-1 do j=1,3 c ii=3*nc*i+3*k+j-nc3-3 ii=ii+1 xx(ii)=x(i,j)+(x(i2,j)-x(i,j))*1.0*k/nc enddo enddo enddo do 2010 i1=1,n-interval c it's only for circular i2max=min(n,n-interval+i1) i2=i1+interval do 2000 while (i2.le.i2max) c calculate rii2 (i1,i2) and compare with rmax2 rii2=0 do j=1,3 r_re(j)=x(i2,j)-x(i1,j) rii2=r_re(j)*r_re(j)+rii2 enddo r_rel=dsqrt(rii2) if (r_rel.ge.rmax2) then itemp=(r_rel-rmax2)/r0 i2=i2+max(1,itemp) else do j=1,3 fe12(j)=0 enddo c calculate the interaction between two subcharge c if the distance is less than rmax c first the k1=1 k2=1 the r_re() has been calculated if (r_rel.le.rmax) then if (r_rel.le.ddna) then c excluded volume potential write(*,*) 'need a push' write(*,*) 'ntime:',ntime write(*,*) 'i1,i2,k1,k2' write(*,*) i1,i2,1,1 write(*,*) 'r_rel:',r_rel write(*,*) 'x(i1)',x(i1,1),x(i1,2),x(i1,3) write(*,*) 'x(i2)',x(i2,1),x(i2,2),x(i2,3) if (i2.eq.n) then iflag=0 else iflag=1 endif ii1=i1*3-3 ii2=i2*3-3 do j=1,3 c can be modified into more realistic potential c problem. dx (timestep=450ps) ~ 0.22 (10nm) art(j)=f_art*r_re(j)/r_rel ii1=ii1+1 ii2=ii2+1 c gc(ii1)=gc(ii1)-art(j) gc(ii1)=gc(ii1)-art(j) c gc(ii2*iflag+j)=gc(ii2*iflag+j)+art(j) gc(ii2)=gc(ii2)+art(j) enddo c goto 2020 endif yy=debye*r_rel fe_rel=econst*exp(-yy)*(yy+1)/(r_rel*r_rel*r_rel) ii1=i1*3-3 ii2=i2*3-3 do j=1,3 fe12(j)=fe_rel*r_re(j) ii1=ii1+1 ii2=ii2+1 gc(ii1)=gc(ii1)-fe12(j) gc(ii2)=gc(ii2)+fe12(j) enddo c endif endif c then for (k1,k2)<>(1,1) do k1=1,nc if (k1.eq.1) then k2min=2 else k2min=1 endif do k2=k2min,nc r_rel=0 do j=1,3 ii1=nc3*i1+k1*3+j-nc33 ii2=nc3*i2+k2*3+j-nc33 r_re(j)=xx(ii2)-xx(ii1) r_rel=r_rel+r_re(j)*r_re(j) enddo r_rel=dsqrt(r_rel) if (r_rel.le.rmax) then if (r_rel.le.ddna) then c excluded volume potential write(*,*) 'need a push' write(*,*) 'ntime:',ntime write(*,*) 'i,i2,k1,k2' write(*,*) i1,i2,k1,k2 write(*,*) 'r_rel:',r_rel ii1=nc3*i1+k1*3+3-nc33 ii2=nc3*i2+k2*3+3-nc33 write(*,*) 'x(i1)',xx(ii1-2),xx(ii1-1),xx(ii1) write(*,*) 'x(i2)',xx(ii2-2),xx(ii2-1),xx(ii2) if (i2.eq.n) then iflag=0 else iflag=1 endif ii1=i1*3 ii2=i2*3 do j=1,3 c can be modified into more realistic potential c problem. dx (timestep=450ps) ~ 0.22 (10nm) art(j)=f_art*r_re(j)/r_rel ii1=ii1+1 gc(ii1)=gc(ii1)-art(j)*(k1-1)*1.0/nc gc(ii1-3)=gc(ii1-3)-art(j)*(nc-k1+1)*1.0/nc gc(ii2*iflag+j)=gc(ii2*iflag+j)+ + art(j)*(k2-1)*1.0/nc gc(ii2+j-3)=gc(ii2+j-3)+art(j) + *(nc-k2+1)*1.0/nc enddo c goto 2020 endif yy=debye*r_rel fe_rel=econst*exp(-yy)*(yy+1)/(r_rel*r_rel*r_rel) if (i2.eq.n) then iflag=0 else iflag=1 endif ii1=i1*3 ii2=i2*3 do j=1,3 fe12(j)=fe_rel*r_re(j) ii1=ii1+1 gc(ii1)=gc(ii1)-fe12(j)*(k1-1)*1.0/nc gc(ii1-3)=gc(ii1-3)-fe12(j)*(nc-k1+1)*1.0/nc gc(ii2*iflag+j)=gc(ii2*iflag+j)+ + fe12(j)*(k2-1)*1.0/nc gc(ii2+j-3)=gc(ii2+j-3)+fe12(j)*(nc-k2+1)*1.0/nc enddo c endif endif enddo enddo 2020 i2=i2+1 endif 2000 enddo 2010 enddo return end