c force_electro_linear.f
c force_electro_linear(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_linear(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-1
         i2=i+1
         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))*k/nc
            enddo
         enddo
      enddo

      do 2010 i1=1,n-interval-1
        i2=i1+interval
        i2max=n-1
          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
                             ii2=i2*3-3
                          ii1=i1*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
                             gc(ii1)=gc(ii1)-art(j)
                             gc(ii2)=gc(ii2)+art(j)
                          enddo
                          goto 2020
                         else
                    yy=debye*r_rel
                    fe_rel=econst*exp(-yy)*(yy+1)/(r_rel*r_rel*r_rel)
                    
                             ii2=i2*3-3
                          ii1=i1*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
                   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
                             ii2=i2*3
                          ii1=i1*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+j)=gc(ii2+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
                          goto 2020
                        else
                     yy=debye*r_rel
                     fe_rel=econst*exp(-yy)*(yy+1)/(r_rel*r_rel*r_rel)

                    
                             ii2=i2*3
                          ii1=i1*3
                          do j=1,3
                             fe12(j)=fe_rel*r_re(j)
                             ii1=ii1+1
                             ii2=ii2+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+j)=gc(ii2+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
                      endif
                    endif
                 enddo
                enddo
 2020          i2=i2+1
             endif
 2000     enddo
 2010   enddo                          

        return
        end