c force_linear.f
c add the streching and bending forces to gc(i)

c force_linear(n,gc,ntime)
c gc(maxn) using to store the force_linear vector
c refurbish force_linear.f only i=1..n is used now.
c
c parameters used: ss1,bendc
c note: when use r(i,i+1) don't use r(i+1,i) instead (may not be updated)
	subroutine force_linear(n,gc,ntime)
	
	implicit double precision(a-h, o-z)
        parameter(maxa=500,maxn=maxa*3)
 	double precision gc(maxn),A(maxa,3),B(maxa,3),cost(maxa)
        common/matrix/x(maxa,3),bangle(maxa),rz(maxa)
        common/matrix2/ r(maxa,maxa)
	common/ddx/d_x(maxa,3)
       common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma,
     +amass,am(maxn),coef,dtime2
        common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)

c calculate the stretching and bending forces
c 1. calculate bangle(i)  1..n-1
	do i=1,n-1
	 cost(i)=0;  ii=i+1
	 do j=1,3
	 cost(i)=cost(i)+u(i,j)*u(ii,j)
	 enddo
	bangle(i)=acos(cost(i))
	enddo

c 2. calculate A    A(0)=A(n-1)=A(n)=0 ends for linear DNA
	do i=1,n-2
         ii=i+1
	 if (bangle(i).eq.0) then
	  tmp=1/r(i,ii)
	 else 
	  tmp=bangle(i)/r(i,ii)/sin(bangle(i))
	 endif
	 do j=1,3
	 A(i,j)=tmp*(u(ii,j)-u(i,j)*cost(i))
	 enddo
	enddo
        
	 do j=1,3
	 A(n,j)=0; A(n-1,j)=0 
	 enddo
	
c 3.calculate B    B(0)=B(1)=B(n)=0
	do i=2,n-1
         ii=i-1  
	 if (bangle(ii).eq.0) then
	  tmp=1/r(i,i+1)
	 else
	  tmp=bangle(ii)/r(i,i+1)/sin(bangle(ii))
	 endif
	 do j=1,3
	 B(i,j)=tmp*(u(ii,j)-u(i,j)*cost(ii))
	 enddo
	enddo

	 do j=1,3
	 B(n,j)=0; B(1,j)=0
	 enddo

c calculate force  fs+fb
c	write(*,*) 'ss1,bendc,u,A,f',ss1,bendc,u(1,1),A(1,1),f(1,1)
	ii=3
	do i=2,n-1
	 i3=i-1; i2=i+1
            if (i.eq.n) then
             i2=1
            endif
	  do j=1,3
	    ii=ii+1 
c	    ii=3*i+j-3
	    gc(ii)=gc(ii)+ss1*((r(i,i2)-r0)*u(i,j)-(r(i3,i)-r0)*u(i3,j))
            gc(ii)=gc(ii)-bendc*(A(i,j)-A(i3,j)+B(i,j)-B(i3,j))
c	 gc(i,j)=gc(i,j)+ss1*((r(i,i+1)-r0)*u(i,j)-(l(i-1)-r0)*u(i-1,j))
c        gc(i,j)=gc(i,j)-bendc*(A(i,j)-A(i-1,j)+B(i,j)-B(i-1,j))
	   enddo
	enddo

c for bead i=n
	do j=1,3 
	   ii=ii+1
	    gc(ii)=gc(ii)-ss1*(r(n-1,n)-r0)*u(n-1,j)
            gc(ii)=gc(ii)+bendc*B(n-1,j)
	enddo

c for bead i=1
	do j=1,3 
	    gc(j)=gc(j)+ss1*(r(1,2)-r0)*u(1,j)
            gc(j)=gc(j)-bendc*A(1,j)
	enddo

	return
	end