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

c force_circular(n,gc,ntime)
c gc(maxn) using to store the force_circular vector
c refurbish force_circular.f only i=1..n is used now.

	subroutine force_circular(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)
	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
	 if (cost(i).gt.1.0) cost(i)=1.0
	bangle(i)=acos(cost(i))
	enddo
c    then i=n case
      	 cost(n)=0
	 do j=1,3
	 cost(n)=cost(n)+u(n,j)*u(1,j)
	 enddo
	 bangle(n)=acos(cost(n))
c 2. calculate A    circular DNA no ends
	do i=1,n-1
         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
         if (bangle(n).eq.0) then
	  tmp=1/r(n,1)
	 else 
	  tmp=bangle(n)/r(n,1)/sin(bangle(n))
	 endif
	 do j=1,3
	 A(n,j)=tmp*(u(1,j)-u(n,j)*cost(n))
	 enddo
	
c 3.calculate B   
	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
c B(1)       ii=n
	 if (bangle(n).eq.0) then
	  tmp=1/r(1,2)
	 else
	  tmp=bangle(n)/r(1,2)/sin(bangle(n))
	 endif
	 do j=1,3
	 B(1,j)=tmp*(u(n,j)-u(1,j)*cost(n))
	 enddo
c B(n)       
	 ii=n-1
	 if (bangle(ii).eq.0) then
	  tmp=1/r(n,1)
	 else
	  tmp=bangle(ii)/r(n,1)/sin(bangle(ii))
	 endif
	 do j=1,3
	 B(n,j)=tmp*(u(ii,j)-u(n,j)*cost(ii))
	 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=0
	do i=1,n
	 i3=i-1; i2=i+1
            if (i.eq.n) then
             i2=1
            endif
	    if (i.eq.1) then
	     i3=n
	    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))
     +  -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


	return
	end