c*******************************************************
c  out_energy(n,ntime)
c
      subroutine  out_energy(n,ntime)
      
      implicit double precision(a-h, o-z)

      parameter(maxa=500,maxn=maxa*3,ncmax=5)
      double precision e_elec,e_bend,e_tor
      double precision xx(ncmax*maxn),fe12(3),r_re(3),cost(maxa)
        common/matrix/x(maxa,3),bangle(maxa),rz(maxa)
       common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma,
     +amass,am(maxn),coef,dtime2
	common/tor_param/dr,cg,phi0,dlk,phii,subseg
	common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
        common/euler/alpha(maxa),beta(maxa),gama(maxa),phi(maxa)
	
       common/debyehu/econst,debye,f_art,ddna,rmax,rmax2,rmin,nc,nc3
 
      
c e_bend
       write(*,*) 'energy calculation!'
       e_bend=0
	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))
	
        e_bend=e_bend+bendc*0.5*bangle(i)*bangle(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))
          e_bend=e_bend+bendc*0.5*bangle(n)*bangle(n)
          write(*,*) 'e_bend=',e_bend

c e_tor
       e_tor=0
       do i=1,n
          e_tor=e_tor+cg*0.5*phi(i)*phi(i)
       enddo
c 	do i=1,n-1
c	e_tor=e_tor+0.5*cg*(phi(i)-phi(i+1))*(phi(i)-phi(i+1))
c	enddo
c	e_tor=e_tor+0.5*cg*(phi(1)-phi(n))*(phi(1)-phi(n))
       write(*,*) 'e_tor=',e_tor
	write(*,*) 'phi(1,2,n):',phi(1),phi(2),phi(n)
c e_elec
c update xx(ncmax*maxn)  xx(i,k,j)  i=1:n  k=1:nc  j=1:3      
	nc33=nc3+3
      e_elec=0
      ii=0
      do i=1,n
         i2=i+1
         if (i.eq.n) then
            i2=1
         endif
         do k=1,nc
            do j=1,3
               ii=ii+1
               xx(ii)=x(i,j)+(x(i2,j)-x(i,j))*(k-1)/nc
            enddo
         enddo
      enddo 
c      write(*,*) 'xx updated!'
      ncount=0
      do i1=1,n-2
        i2max=min(n,n-2+i1)
        i2=i1+2
          do 2000 while (i2.le.i2max)
             ncount=ncount+1
             do k1=1,nc
                do k2=1,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)**2
                    enddo
                      r_rel=dsqrt(r_rel)
                     yy=debye*r_rel
                     etmp=econst*exp(-yy)/r_rel
                   e_elec=e_elec+etmp
                   if (etmp.gt.1e-3) then
                      write(111,*) i1,i2,k1,k2,etmp
                   endif
                   
                 enddo
              enddo
              i2=i2+1
             
 2000      enddo
c          write(*,*) i1,e_elec
       enddo
       write(*,*) 'ncount=',ncount

       e_elec2=0
      do i1=1,n-2
        i2max=min(n,n-2+i1)
        i2=i1+2
          do 4000 while (i2.le.i2max)
             do k1=1,nc
                do k2=1,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)**2
                    enddo
                      r_rel=dsqrt(r_rel)
                   if (r_rel.lt.rmax) then
                      yy=debye*r_rel
                   e_elec2=e_elec2+econst*exp(-yy)/r_rel
                  endif
                 enddo
              enddo
              i2=i2+1
             
 4000      enddo
c          write(*,*) i1,e_elec
       enddo
       
       

       write(787,*) e_bend,e_elec,e_elec2
	call flush(787)
	

       return
       end