X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=1d241159ebdc692cc85d242be6868114c8c4c1e7;hb=d2088a7676356c745c2cff79ca5a7c2c6025e19c;hp=4a763c20c4506c21ab6916269d8b2848e0eab645;hpb=e837223cef419e16dfd7fc71baee35cd807096e9;p=unres4.git

diff --git a/source/unres/energy.F90 b/source/unres/energy.F90
index 4a763c2..1d24115 100644
--- a/source/unres/energy.F90
+++ b/source/unres/energy.F90
@@ -74,7 +74,7 @@
 ! amino-acid residue.
 !      common /precomp1/
       real(kind=8),dimension(:,:),allocatable :: mu,muder,Ub2,Ub2der,&
-       Ctobr,Ctobrder,Dtobr2,Dtobr2der      !(2,maxres)
+       Ctobr,Ctobrder,Dtobr2,Dtobr2der,gUb2      !(2,maxres)
       real(kind=8),dimension(:,:,:),allocatable :: EUg,EUgder,CUg,&
        CUgder,DUg,Dugder,DtUg2,DtUg2der      !(2,2,maxres)
 ! This common block contains vectors and matrices dependent on two
@@ -87,6 +87,7 @@
       real(kind=8),dimension(:,:,:,:),allocatable :: Ug2DtEUgder,&
        DtUg2EUgder      !(2,2,2,maxres)
 !      common /rotat_old/
+      real(kind=8),dimension(4) :: gmuij,gmuij1,gmuij2,gmuji1,gmuji2
       real(kind=8),dimension(:),allocatable :: costab,sintab,&
        costab2,sintab2      !(maxres)
 ! This common block contains dipole-interaction matrices and their 
@@ -293,7 +294,7 @@
 !          print *,"Processor",myrank," BROADCAST iorder"
 ! FG master sets up the WEIGHTS_ array which will be broadcast to the 
 ! FG slaves as WEIGHTS array.
-         ! weights_(1)=wsc
+          weights_(1)=wsc
           weights_(2)=wscp
           weights_(3)=welec
           weights_(4)=wcorr
@@ -609,9 +610,9 @@
              .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
              .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
 #endif
-            write(iout,*),"just befor eelec call"
+!            print *,"just befor eelec call"
             call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-!         write (iout,*) "ELEC calc"
+!            print *, "ELEC calc"
          else
             ees=0.0d0
             evdw1=0.0d0
@@ -662,12 +663,26 @@
 ! Calculate the virtual-bond-angle energy.
 !       write(iout,*) "in etotal afer edis",ipot
 
-      if (wang.gt.0.0d0) then
-        call ebend(ebe,ethetacnstr)
+!      if (wang.gt.0.0d0) then
+!        call ebend(ebe,ethetacnstr)
+!      else
+!        ebe=0
+!        ethetacnstr=0
+!      endif
+      if (wang.gt.0d0) then
+       if (tor_mode.eq.0) then
+         call ebend(ebe)
+       else
+!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+         call ebend_kcc(ebe)
+       endif
       else
-        ebe=0
-        ethetacnstr=0
+        ebe=0.0d0
       endif
+      ethetacnstr=0.0d0
+      if (with_theta_constr) call etheta_constr(ethetacnstr)
+
 !       write(iout,*) "in etotal afer ebe",ipot
 
 !      print *,"Processor",myrank," computed UB"
@@ -681,12 +696,27 @@
 ! Calculate the virtual-bond torsional energy.
 !
 !d    print *,'nterm=',nterm
-      if (wtor.gt.0) then
-       call etor(etors,edihcnstr)
+!      if (wtor.gt.0) then
+!       call etor(etors,edihcnstr)
+!      else
+!       etors=0
+!       edihcnstr=0
+!      endif
+      if (wtor.gt.0.0d0) then
+         if (tor_mode.eq.0) then
+           call etor(etors)
+         else
+!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+           call etor_kcc(etors)
+         endif
       else
-       etors=0
-       edihcnstr=0
+        etors=0.0d0
       endif
+      edihcnstr=0.0d0
+      if (ndih_constr.gt.0) call etor_constr(edihcnstr)
+!c      print *,"Processor",myrank," computed Utor"
+
 !      print *,"Processor",myrank," computed Utor"
        
 !
@@ -801,6 +831,7 @@
        eespp=0.0d0
       endif
 !      write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
+!      print *,"before ecatcat"
       if (nfgtasks.gt.1) then
       if (fg_rank.eq.0) then
       call ecatcat(ecationcation)
@@ -2689,18 +2720,168 @@
 !      include 'COMMON.FFIELD'
       real(kind=8) :: auxvec(2),auxmat(2,2)
       integer :: i,iti1,iti,k,l
-      real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2
+      real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2,cost1,sint1,&
+       sint1sq,sint1cub,sint1cost1,b1k,b2k,aux
 !       print *,"in set matrices"
 !
 ! Compute the virtual-bond-torsional-angle dependent quantities needed
 ! to calculate the el-loc multibody terms of various order.
 !
 !AL el      mu=0.0d0
+   
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
 #else
       do i=3,nres+1
 #endif
+        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          if (itype(i-2,1).eq.0) then 
+          iti = nloctyp
+          else
+          iti = itype2loc(itype(i-2,1))
+          endif
+        else
+          iti=nloctyp
+        endif
+!c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+          iti1 = itype2loc(itype(i-1,1))
+        else
+          iti1=nloctyp
+        endif
+!        print *,i,itype(i-2,1),iti
+#ifdef NEWCORR
+        cost1=dcos(theta(i-1))
+        sint1=dsin(theta(i-1))
+        sint1sq=sint1*sint1
+        sint1cub=sint1sq*sint1
+        sint1cost1=2*sint1*cost1
+!        print *,"cost1",cost1,theta(i-1)
+!c        write (iout,*) "bnew1",i,iti
+!c        write (iout,*) (bnew1(k,1,iti),k=1,3)
+!c        write (iout,*) (bnew1(k,2,iti),k=1,3)
+!c        write (iout,*) "bnew2",i,iti
+!c        write (iout,*) (bnew2(k,1,iti),k=1,3)
+!c        write (iout,*) (bnew2(k,2,iti),k=1,3)
+        k=1
+!        print *,bnew1(1,k,iti),"bnew1"
+        do k=1,2
+          b1k=bnew1(1,k,iti)+(bnew1(2,k,iti)+bnew1(3,k,iti)*cost1)*cost1
+!          print *,b1k
+!          write(*,*) shape(b1) 
+!          if(.not.allocated(b1)) print *, "WTF?"
+          b1(k,i-2)=sint1*b1k
+!
+!             print *,b1(k,i-2)
+
+          gtb1(k,i-2)=cost1*b1k-sint1sq*&
+                   (bnew1(2,k,iti)+2*bnew1(3,k,iti)*cost1)
+!             print *,gtb1(k,i-2)
+
+          b2k=bnew2(1,k,iti)+(bnew2(2,k,iti)+bnew2(3,k,iti)*cost1)*cost1
+          b2(k,i-2)=sint1*b2k
+!             print *,b2(k,i-2)
+
+          gtb2(k,i-2)=cost1*b2k-sint1sq*&
+                   (bnew2(2,k,iti)+2*bnew2(3,k,iti)*cost1)
+!             print *,gtb2(k,i-2)
+
+        enddo
+!        print *,b1k,b2k
+        do k=1,2
+          aux=ccnew(1,k,iti)+(ccnew(2,k,iti)+ccnew(3,k,iti)*cost1)*cost1
+          cc(1,k,i-2)=sint1sq*aux
+          gtcc(1,k,i-2)=sint1cost1*aux-sint1cub*&
+                   (ccnew(2,k,iti)+2*ccnew(3,k,iti)*cost1)
+          aux=ddnew(1,k,iti)+(ddnew(2,k,iti)+ddnew(3,k,iti)*cost1)*cost1
+          dd(1,k,i-2)=sint1sq*aux
+          gtdd(1,k,i-2)=sint1cost1*aux-sint1cub*&
+                   (ddnew(2,k,iti)+2*ddnew(3,k,iti)*cost1)
+        enddo
+!        print *,"after cc"
+        cc(2,1,i-2)=cc(1,2,i-2)
+        cc(2,2,i-2)=-cc(1,1,i-2)
+        gtcc(2,1,i-2)=gtcc(1,2,i-2)
+        gtcc(2,2,i-2)=-gtcc(1,1,i-2)
+        dd(2,1,i-2)=dd(1,2,i-2)
+        dd(2,2,i-2)=-dd(1,1,i-2)
+        gtdd(2,1,i-2)=gtdd(1,2,i-2)
+        gtdd(2,2,i-2)=-gtdd(1,1,i-2)
+!        print *,"after dd"
+
+        do k=1,2
+          do l=1,2
+            aux=eenew(1,l,k,iti)+eenew(2,l,k,iti)*cost1
+            EE(l,k,i-2)=sint1sq*aux
+            gtEE(l,k,i-2)=sint1cost1*aux-sint1cub*eenew(2,l,k,iti)
+          enddo
+        enddo
+        EE(1,1,i-2)=EE(1,1,i-2)+e0new(1,iti)*cost1
+        EE(1,2,i-2)=EE(1,2,i-2)+e0new(2,iti)+e0new(3,iti)*cost1
+        EE(2,1,i-2)=EE(2,1,i-2)+e0new(2,iti)*cost1+e0new(3,iti)
+        EE(2,2,i-2)=EE(2,2,i-2)-e0new(1,iti)
+        gtEE(1,1,i-2)=gtEE(1,1,i-2)-e0new(1,iti)*sint1
+        gtEE(1,2,i-2)=gtEE(1,2,i-2)-e0new(3,iti)*sint1
+        gtEE(2,1,i-2)=gtEE(2,1,i-2)-e0new(2,iti)*sint1
+!        print *,"after ee"
+
+!c        b1tilde(1,i-2)=b1(1,i-2)
+!c        b1tilde(2,i-2)=-b1(2,i-2)
+!c        b2tilde(1,i-2)=b2(1,i-2)
+!c        b2tilde(2,i-2)=-b2(2,i-2)
+#ifdef DEBUG
+        write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
+        write(iout,*)  'b1=',(b1(k,i-2),k=1,2)
+        write(iout,*)  'b2=',(b2(k,i-2),k=1,2)
+        write (iout,*) 'theta=', theta(i-1)
+#endif
+#else
+        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          iti = itype2loc(itype(i-2,1))
+        else
+          iti=nloctyp
+        endif
+!c        write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti
+!c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+          iti1 = itype2loc(itype(i-1,1))
+        else
+          iti1=nloctyp
+        endif
+!        print *,i,iti
+        b1(1,i-2)=b(3,iti)
+        b1(2,i-2)=b(5,iti)
+        b2(1,i-2)=b(2,iti)
+        b2(2,i-2)=b(4,iti)
+        do k=1,2
+          do l=1,2
+           CC(k,l,i-2)=ccold(k,l,iti)
+           DD(k,l,i-2)=ddold(k,l,iti)
+           EE(k,l,i-2)=eeold(k,l,iti)
+          enddo
+        enddo
+#endif
+        b1tilde(1,i-2)= b1(1,i-2)
+        b1tilde(2,i-2)=-b1(2,i-2)
+        b2tilde(1,i-2)= b2(1,i-2)
+        b2tilde(2,i-2)=-b2(2,i-2)
+!c
+        Ctilde(1,1,i-2)= CC(1,1,i-2)
+        Ctilde(1,2,i-2)= CC(1,2,i-2)
+        Ctilde(2,1,i-2)=-CC(2,1,i-2)
+        Ctilde(2,2,i-2)=-CC(2,2,i-2)
+!c
+        Dtilde(1,1,i-2)= DD(1,1,i-2)
+        Dtilde(1,2,i-2)= DD(1,2,i-2)
+        Dtilde(2,1,i-2)=-DD(2,1,i-2)
+        Dtilde(2,2,i-2)=-DD(2,2,i-2)
+      enddo
+#ifdef PARMAT
+      do i=ivec_start+2,ivec_end+2
+#else
+      do i=3,nres+1
+#endif
+
 !      print *,i,"i"
         if (i .lt. nres+1) then
           sin1=dsin(phi(i))
@@ -2773,37 +2954,43 @@
            if (itype(i-2,1).eq.0) then
           iti=ntortyp+1
            else
-          iti = itortyp(itype(i-2,1))
+          iti = itype2loc(itype(i-2,1))
            endif
         else
-          iti=ntortyp+1
+          iti=nloctyp
         endif
 !        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         if (i.gt. nnt+1 .and. i.lt.nct+1) then
            if (itype(i-1,1).eq.0) then
-          iti1=ntortyp+1
+          iti1=nloctyp
            else
-          iti1 = itortyp(itype(i-1,1))
+          iti1 = itype2loc(itype(i-1,1))
            endif
         else
-          iti1=ntortyp+1
+          iti1=nloctyp
         endif
 !          print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1)
 !d        write (iout,*) '*******i',i,' iti1',iti
-!d        write (iout,*) 'b1',b1(:,iti)
-!d        write (iout,*) 'b2',b2(:,iti)
+!        write (iout,*) 'b1',b1(:,iti)
+!        write (iout,*) 'b2',b2(:,i-2)
 !d        write (iout,*) 'Ug',Ug(:,:,i-2)
 !        if (i .gt. iatel_s+2) then
         if (i .gt. nnt+2) then
-          call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2))
-          call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
+          call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2))
+#ifdef NEWCORR
+          call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2))
+!c          write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
+#endif
+
+          call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2))
+          call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2))
           if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) &
           then
-          call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
-          call matmat2(DD(1,1,iti),Ug(1,1,i-2),DUg(1,1,i-2))
-          call matmat2(Dtilde(1,1,iti),Ug2(1,1,i-2),DtUg2(1,1,i-2))
-          call matvec2(Ctilde(1,1,iti1),obrot(1,i-2),Ctobr(1,i-2))
-          call matvec2(Dtilde(1,1,iti),obrot2(1,i-2),Dtobr2(1,i-2))
+          call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2))
+          call matmat2(DD(1,1,i-2),Ug(1,1,i-2),DUg(1,1,i-2))
+          call matmat2(Dtilde(1,1,i-2),Ug2(1,1,i-2),DtUg2(1,1,i-2))
+          call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2))
+          call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2))
           endif
         else
           do k=1,2
@@ -2818,8 +3005,8 @@
             enddo
           enddo
         endif
-        call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2))
-        call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
+        call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2))
+        call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2))
         do k=1,2
           muder(k,i-2)=Ub2der(k,i-2)
         enddo
@@ -2828,34 +3015,34 @@
           if (itype(i-1,1).eq.0) then
            iti1=ntortyp+1
           elseif (itype(i-1,1).le.ntyp) then
-            iti1 = itortyp(itype(i-1,1))
+            iti1 = itype2loc(itype(i-1,1))
           else
-            iti1=ntortyp+1
+            iti1=nloctyp
           endif
         else
-          iti1=ntortyp+1
+          iti1=nloctyp
         endif
         do k=1,2
-          mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
+          mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
         enddo
-!        if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2)
-!        if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,iti1)
-!        if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2)
+        if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2)
+        if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,i-1)
+        if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2)
 !d        write (iout,*) 'mu1',mu1(:,i-2)
 !d        write (iout,*) 'mu2',mu2(:,i-2)
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
         then  
-        call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2))
-        call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2))
-        call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
-        call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
-        call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
+        call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
+        call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2))
+        call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
+        call matvec2(Ctilde(1,1,i-1),obrot_der(1,i-2),Ctobrder(1,i-2))
+        call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2))
 ! Vectors and matrices dependent on a single virtual-bond dihedral.
-        call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
+        call matvec2(DD(1,1,i-2),b1tilde(1,i-1),auxvec(1))
         call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) 
         call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) 
-        call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2))
-        call matvec2(CC(1,1,iti1),Ub2der(1,i-2),CUgb2der(1,i-2))
+        call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2))
+        call matvec2(CC(1,1,i-1),Ub2der(1,i-2),CUgb2der(1,i-2))
         call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2))
         call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2))
         call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2))
@@ -3476,6 +3663,7 @@
       real(kind=8),dimension(2,2) :: acipa !el,a_temp
 !el      real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
       real(kind=8),dimension(4) :: muij
+      real(kind=8) :: geel_loc_ij,geel_loc_ji
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
                     dist_temp, dist_init,rlocshield,fracinbuf
       integer xshift,yshift,zshift,ilist,iresshield
@@ -3917,6 +4105,15 @@
             do l=1,2
               kkk=kkk+1
               muij(kkk)=mu(k,i)*mu(l,j)
+#ifdef NEWCORR
+             gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
+!c             write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
+             gmuij2(kkk)=gUb2(k,i)*mu(l,j)
+             gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
+!c             write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
+             gmuji2(kkk)=mu(k,i)*gUb2(l,j)
+#endif
+
             enddo
           enddo  
 !d         write (iout,*) 'EELEC: i',i,' j',j
@@ -4143,6 +4340,61 @@
            enddo
            endif
 
+#ifdef NEWCORR
+         geel_loc_ij=(a22*gmuij1(1)&
+          +a23*gmuij1(2)&
+          +a32*gmuij1(3)&
+          +a33*gmuij1(4))&
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
+!c         write(iout,*) "derivative over thatai"
+!c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
+!c     &   a33*gmuij1(4) 
+         gloc(nphi+i,icg)=gloc(nphi+i,icg)+&
+           geel_loc_ij*wel_loc
+!c         write(iout,*) "derivative over thatai-1" 
+!c         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
+!c     &   a33*gmuij2(4)
+         geel_loc_ij=&
+          a22*gmuij2(1)&
+          +a23*gmuij2(2)&
+          +a32*gmuij2(3)&
+          +a33*gmuij2(4)
+         gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
+           geel_loc_ij*wel_loc&
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
+
+!c  Derivative over j residue
+         geel_loc_ji=a22*gmuji1(1)&
+          +a23*gmuji1(2)&
+          +a32*gmuji1(3)&
+          +a33*gmuji1(4)
+!c         write(iout,*) "derivative over thataj" 
+!c         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+!c     &   a33*gmuji1(4)
+
+        gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
+           geel_loc_ji*wel_loc&
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
+
+         geel_loc_ji=&
+          +a22*gmuji2(1)&
+          +a23*gmuji2(2)&
+          +a32*gmuji2(3)&
+          +a33*gmuji2(4)
+!c         write(iout,*) "derivative over thataj-1"
+!c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+!c     &   a33*gmuji2(4)
+         gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
+           geel_loc_ji*wel_loc&
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+#endif
 
 !          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 !           eel_loc_ij=0.0
@@ -4501,7 +4753,9 @@
 !      include 'COMMON.CONTROL'
       real(kind=8),dimension(3) :: ggg
       real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,&
-        e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2
+        e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,gpizda1,&
+       gpizda2,auxgmat1,auxgmatt1,auxgmat2,auxgmatt2
+
       real(kind=8),dimension(2) :: auxvec,auxvec1
 !el      real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
       real(kind=8),dimension(2,2) :: auxmat3 !el, a_temp
@@ -4561,8 +4815,15 @@
 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
 !d        call checkint_turn3(i,a_temp,eello_turn3_num)
         call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
+        call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1))
+        call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1))
         call transpose2(auxmat(1,1),auxmat1(1,1))
+        call transpose2(auxgmat1(1,1),auxgmatt1(1,1))
+        call transpose2(auxgmat2(1,1),auxgmatt2(1,1))
         call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
+        call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
+        call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
+
         if (shield_mode.eq.0) then
         fac_shield(i)=1.0d0
         fac_shield(j)=1.0d0
@@ -4577,6 +4838,18 @@
 
         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
                'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
+!C#ifdef NEWCORR
+!C Derivatives in theta
+        gloc(nphi+i,icg)=gloc(nphi+i,icg) &
+       +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3&
+        *fac_shield(i)*fac_shield(j)
+        gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)&
+       +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3&
+        *fac_shield(i)*fac_shield(j)
+!C#endif
+
+
+
           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. &
        (shield_mode.gt.0)) then
 !C          print *,i,j     
@@ -4714,8 +4987,14 @@
 !      include 'COMMON.CONTROL'
       real(kind=8),dimension(3) :: ggg
       real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,&
-        e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2
-      real(kind=8),dimension(2) :: auxvec,auxvec1
+        e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,& 
+        gte1t,gte2t,gte3t,&
+        gte1a,gtae3,gtae3e2, ae3gte2,&
+        gtEpizda1,gtEpizda2,gtEpizda3
+
+      real(kind=8),dimension(2) :: auxvec,auxvec1,auxgEvec1,auxgEvec2,&
+       auxgEvec3,auxgvec
+
 !el      real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
       real(kind=8),dimension(2,2) :: auxmat3 !el a_temp
 !el      real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
@@ -4727,7 +5006,7 @@
 !el local variables
       integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield
       real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,&
-         rlocshield
+         rlocshield,gs23,gs32,gsE13,gs13,gs21,gsE31,gsEE1,gsEE2,gsEE3
       
       j=i+3
 !      if (j.ne.20) return
@@ -4775,22 +5054,63 @@
         a_temp(1,2)=a23
         a_temp(2,1)=a32
         a_temp(2,2)=a33
-        iti1=itortyp(itype(i+1,1))
-        iti2=itortyp(itype(i+2,1))
-        iti3=itortyp(itype(i+3,1))
+        iti1=i+1
+        iti2=i+2
+        iti3=i+3
 !        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
         call transpose2(EUg(1,1,i+1),e1t(1,1))
         call transpose2(Eug(1,1,i+2),e2t(1,1))
         call transpose2(Eug(1,1,i+3),e3t(1,1))
+!C Ematrix derivative in theta
+        call transpose2(gtEUg(1,1,i+1),gte1t(1,1))
+        call transpose2(gtEug(1,1,i+2),gte2t(1,1))
+        call transpose2(gtEug(1,1,i+3),gte3t(1,1))
+
         call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
         call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
+        call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1))
+        call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1))
+!c       auxalary matrix of E i+1
+        call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1))
         s1=scalar2(b1(1,iti2),auxvec(1))
+!c derivative of theta i+2 with constant i+3
+        gs23=scalar2(gtb1(1,i+2),auxvec(1))
+!c derivative of theta i+2 with constant i+2
+        gs32=scalar2(b1(1,i+2),auxgvec(1))
+!c derivative of E matix in theta of i+1
+        gsE13=scalar2(b1(1,i+2),auxgEvec1(1))
+
         call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
+        call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1))
         call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
-        s2=scalar2(b1(1,iti1),auxvec(1))
+!c auxilary matrix auxgvec of Ub2 with constant E matirx
+        call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
+!c auxilary matrix auxgEvec1 of E matix with Ub2 constant
+        call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1))
+        s2=scalar2(b1(1,i+1),auxvec(1))
+!c derivative of theta i+1 with constant i+3
+        gs13=scalar2(gtb1(1,i+1),auxvec(1))
+!c derivative of theta i+2 with constant i+1
+        gs21=scalar2(b1(1,i+1),auxgvec(1))
+!c derivative of theta i+3 with constant i+1
+        gsE31=scalar2(b1(1,i+1),auxgEvec3(1))
+
         call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
+        call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1))
+!c ae3gte2 is derivative over i+2
+        call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1))
+
         call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
+        call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1))
+!c i+2
+        call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1))
+!c i+3
+        call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1))
+
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
+        gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
+        gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
+        gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
         if (shield_mode.eq.0) then
         fac_shield(i)=1.0
         fac_shield(j)=1.0
@@ -4844,7 +5164,21 @@
 !           print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1)
            enddo
            endif
-
+#ifdef NEWCORR
+        gloc(nphi+i,icg)=gloc(nphi+i,icg)&
+                       -(gs13+gsE13+gsEE1)*wturn4&
+       *fac_shield(i)*fac_shield(j)
+        gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)&
+                         -(gs23+gs21+gsEE2)*wturn4&
+       *fac_shield(i)*fac_shield(j)
+
+        gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)&
+                         -(gs32+gsE31+gsEE3)*wturn4&
+       *fac_shield(i)*fac_shield(j)
+
+!c         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
+!c     &   gs2
+#endif
         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
            'eturn4',i,j,-(s1+s2+s3)
 !d        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
@@ -4853,7 +5187,7 @@
         call transpose2(EUgder(1,1,i+1),e1tder(1,1))
         call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
         call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
-        s1=scalar2(b1(1,iti2),auxvec(1))
+        s1=scalar2(b1(1,i+1),auxvec(1))
         call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) &
@@ -5897,7 +6231,7 @@
       end subroutine theteng
 #else
 !-----------------------------------------------------------------------------
-      subroutine ebend(etheta,ethetacnstr)
+      subroutine ebend(etheta)
 !
 ! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
 ! angles gamma and its derivatives in consecutive thetas and gammas.
@@ -6103,34 +6437,6 @@
       enddo
 !-----------thete constrains
 !      if (tor_mode.ne.2) then
-      ethetacnstr=0.0d0
-      print *,ithetaconstr_start,ithetaconstr_end,"TU"
-      do i=ithetaconstr_start,ithetaconstr_end
-        itheta=itheta_constr(i)
-        thetiii=theta(itheta)
-        difi=pinorm(thetiii-theta_constr0(i))
-        if (difi.gt.theta_drange(i)) then
-          difi=difi-theta_drange(i)
-          ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
-          gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
-         +for_thet_constr(i)*difi**3
-        else if (difi.lt.-drange(i)) then
-          difi=difi+drange(i)
-          ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
-          gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
-         +for_thet_constr(i)*difi**3
-        else
-          difi=0.0
-        endif
-       if (energy_dec) then
-        write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", &
-         i,itheta,rad2deg*thetiii, &
-         rad2deg*theta_constr0(i),  rad2deg*theta_drange(i), &
-         rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, &
-         gloc(itheta+nphi-2,icg)
-        endif
-      enddo
-!      endif
 
       return
       end subroutine ebend
@@ -7019,7 +7325,7 @@
       end subroutine etor_d
 #else
 !-----------------------------------------------------------------------------
-      subroutine etor(etors,edihcnstr)
+      subroutine etor(etors)
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.VAR'
@@ -7101,8 +7407,156 @@
 !       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo
 ! 6/20/98 - dihedral angle constraints
-      edihcnstr=0.0d0
-!      do i=1,ndih_constr
+      return
+      end subroutine etor
+!C The rigorous attempt to derive energy function
+!-------------------------------------------------------------------------------------------
+      subroutine etor_kcc(etors)
+      double precision c1(0:maxval_kcc),c2(0:maxval_kcc)
+      real(kind=8) :: etors,glocig,glocit1,glocit2,sinthet1,&
+       sinthet2,costhet1,costhet2,sint1t2,sint1t2n,phii,sinphi,cosphi,&
+       sint1t2n1,sumvalc,gradvalct1,gradvalct2,sumvals,gradvalst1,&
+       gradvalst2,etori
+      logical lprn
+      integer :: i,j,itori,itori1,nval,k,l
+
+      if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode
+      etors=0.0D0
+      do i=iphi_start,iphi_end
+!C ANY TWO ARE DUMMY ATOMS in row CYCLE
+!c        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+!c     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1))  .or.
+!c     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+        if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 &
+           .or. itype(i,1).eq.ntyp1 .or. itype(i-3,1).eq.ntyp1) cycle
+        itori=itortyp(itype(i-2,1))
+        itori1=itortyp(itype(i-1,1))
+        phii=phi(i)
+        glocig=0.0D0
+        glocit1=0.0d0
+        glocit2=0.0d0
+!C to avoid multiple devision by 2
+!c        theti22=0.5d0*theta(i)
+!C theta 12 is the theta_1 /2
+!C theta 22 is theta_2 /2
+!c        theti12=0.5d0*theta(i-1)
+!C and appropriate sinus function
+        sinthet1=dsin(theta(i-1))
+        sinthet2=dsin(theta(i))
+        costhet1=dcos(theta(i-1))
+        costhet2=dcos(theta(i))
+!C to speed up lets store its mutliplication
+        sint1t2=sinthet2*sinthet1
+        sint1t2n=1.0d0
+!C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma)
+!C +d_n*sin(n*gamma)) *
+!C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2))) 
+!C we have two sum 1) Non-Chebyshev which is with n and gamma
+        nval=nterm_kcc_Tb(itori,itori1)
+        c1(0)=0.0d0
+        c2(0)=0.0d0
+        c1(1)=1.0d0
+        c2(1)=1.0d0
+        do j=2,nval
+          c1(j)=c1(j-1)*costhet1
+          c2(j)=c2(j-1)*costhet2
+        enddo
+        etori=0.0d0
+
+       do j=1,nterm_kcc(itori,itori1)
+          cosphi=dcos(j*phii)
+          sinphi=dsin(j*phii)
+          sint1t2n1=sint1t2n
+          sint1t2n=sint1t2n*sint1t2
+          sumvalc=0.0d0
+          gradvalct1=0.0d0
+          gradvalct2=0.0d0
+          do k=1,nval
+            do l=1,nval
+              sumvalc=sumvalc+v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l)
+              gradvalct1=gradvalct1+ &
+                (k-1)*v1_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l)
+              gradvalct2=gradvalct2+ &
+                (l-1)*v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1)
+            enddo
+          enddo
+          gradvalct1=-gradvalct1*sinthet1
+          gradvalct2=-gradvalct2*sinthet2
+          sumvals=0.0d0
+          gradvalst1=0.0d0
+          gradvalst2=0.0d0
+          do k=1,nval
+            do l=1,nval
+              sumvals=sumvals+v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l)
+              gradvalst1=gradvalst1+ &
+                (k-1)*v2_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l)
+              gradvalst2=gradvalst2+ &
+                (l-1)*v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1)
+            enddo
+          enddo
+          gradvalst1=-gradvalst1*sinthet1
+          gradvalst2=-gradvalst2*sinthet2
+          if (lprn) write (iout,*)j,"sumvalc",sumvalc," sumvals",sumvals
+          etori=etori+sint1t2n*(sumvalc*cosphi+sumvals*sinphi)
+!C glocig is the gradient local i site in gamma
+          glocig=glocig+j*sint1t2n*(sumvals*cosphi-sumvalc*sinphi)
+!C now gradient over theta_1
+         glocit1=glocit1+sint1t2n*(gradvalct1*cosphi+gradvalst1*sinphi)&
+        +j*sint1t2n1*costhet1*sinthet2*(sumvalc*cosphi+sumvals*sinphi)
+         glocit2=glocit2+sint1t2n*(gradvalct2*cosphi+gradvalst2*sinphi)&
+        +j*sint1t2n1*sinthet1*costhet2*(sumvalc*cosphi+sumvals*sinphi)
+        enddo ! j
+        etors=etors+etori
+        gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
+!C derivative over theta1
+        gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*glocit1
+!C now derivative over theta2
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2
+        if (lprn) then
+         write (iout,*) i-2,i-1,itype(i-2,1),itype(i-1,1),itori,itori1,&
+            theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori
+          write (iout,*) "c1",(c1(k),k=0,nval), &
+         " c2",(c2(k),k=0,nval)
+        endif
+      enddo
+      return
+       end  subroutine etor_kcc
+!------------------------------------------------------------------------------
+
+        subroutine etor_constr(edihcnstr)
+      real(kind=8) :: etors,edihcnstr
+      logical :: lprn
+!el local variables
+      integer :: i,j,iblock,itori,itori1
+      real(kind=8) :: phii,gloci,v1ij,v2ij,cosphi,sinphi,&
+                   vl1ij,vl2ij,vl3ij,pom1,difi,etors_ii,pom,&
+                   gaudih_i,gauder_i,s,cos_i,dexpcos_i
+
+      if (raw_psipred) then
+        do i=idihconstr_start,idihconstr_end
+          itori=idih_constr(i)
+          phii=phi(itori)
+          gaudih_i=vpsipred(1,i)
+          gauder_i=0.0d0
+          do j=1,2
+            s = sdihed(j,i)
+            cos_i=(1.0d0-dcos(phii-phibound(j,i)))/s**2
+            dexpcos_i=dexp(-cos_i*cos_i)
+            gaudih_i=gaudih_i+vpsipred(j+1,i)*dexpcos_i
+          gauder_i=gauder_i-2*vpsipred(j+1,i)*dsin(phii-phibound(j,i)) &
+                 *cos_i*dexpcos_i/s**2
+          enddo
+          edihcnstr=edihcnstr-wdihc*dlog(gaudih_i)
+          gloc(itori-3,icg)=gloc(itori-3,icg)-wdihc*gauder_i/gaudih_i
+          if (energy_dec) &
+          write (iout,'(2i5,f8.3,f8.5,2(f8.5,2f8.3),f10.5)') &
+          i,itori,phii*rad2deg,vpsipred(1,i),vpsipred(2,i),&
+          phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,vpsipred(3,i),&
+          phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,&
+          -wdihc*dlog(gaudih_i)
+        enddo
+      else
+
       do i=idihconstr_start,idihconstr_end
         itori=idih_constr(i)
         phii=phi(itori)
@@ -7118,13 +7572,13 @@
         else
           difi=0.0
         endif
-!d        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
-!d     &    rad2deg*phi0(i),  rad2deg*drange(i),
-!d     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
       enddo
-!d       write (iout,*) 'edihcnstr',edihcnstr
+
+      endif
+
       return
-      end subroutine etor
+
+      end subroutine etor_constr
 !-----------------------------------------------------------------------------
       subroutine etor_d(etors_d)
 ! 6/23/01 Compute double torsional energy
@@ -7215,6 +7669,75 @@
       return
       end subroutine etor_d
 #endif
+
+      subroutine ebend_kcc(etheta)
+      logical lprn
+      double precision thybt1(maxang_kcc),etheta
+      integer :: i,iti,j,ihelp
+      real (kind=8) :: sinthet,costhet,sumth1thyb,gradthybt1
+!C Set lprn=.true. for debugging
+      lprn=energy_dec
+!c     lprn=.true.
+!C      print *,"wchodze kcc"
+      if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode
+      etheta=0.0D0
+      do i=ithet_start,ithet_end
+!c        print *,i,itype(i-1),itype(i),itype(i-2)
+        if ((itype(i-1,1).eq.ntyp1).or.itype(i-2,1).eq.ntyp1 &
+       .or.itype(i,1).eq.ntyp1) cycle
+        iti=iabs(itortyp(itype(i-1,1)))
+        sinthet=dsin(theta(i))
+        costhet=dcos(theta(i))
+        do j=1,nbend_kcc_Tb(iti)
+          thybt1(j)=v1bend_chyb(j,iti)
+        enddo
+        sumth1thyb=v1bend_chyb(0,iti)+ &
+         tschebyshev(1,nbend_kcc_Tb(iti),thybt1(1),costhet)
+        if (lprn) write (iout,*) i-1,itype(i-1,1),iti,theta(i)*rad2deg,&
+         sumth1thyb
+        ihelp=nbend_kcc_Tb(iti)-1
+        gradthybt1=gradtschebyshev(0,ihelp,thybt1(1),costhet)
+        etheta=etheta+sumth1thyb
+!C        print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0)
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)-wang*gradthybt1*sinthet
+      enddo
+      return
+      end subroutine ebend_kcc
+!c------------
+!c-------------------------------------------------------------------------------------
+      subroutine etheta_constr(ethetacnstr)
+      real (kind=8) :: ethetacnstr,thetiii,difi
+      integer :: i,itheta
+      ethetacnstr=0.0d0
+!C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+      do i=ithetaconstr_start,ithetaconstr_end
+        itheta=itheta_constr(i)
+        thetiii=theta(itheta)
+        difi=pinorm(thetiii-theta_constr0(i))
+        if (difi.gt.theta_drange(i)) then
+          difi=difi-theta_drange(i)
+          ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+          gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
+         +for_thet_constr(i)*difi**3
+        else if (difi.lt.-drange(i)) then
+          difi=difi+drange(i)
+          ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+          gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
+          +for_thet_constr(i)*difi**3
+        else
+          difi=0.0
+        endif
+       if (energy_dec) then
+        write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",&
+         i,itheta,rad2deg*thetiii,&
+         rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),&
+         rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,&
+         gloc(itheta+nphi-2,icg)
+        endif
+      enddo
+      return
+      end subroutine etheta_constr
+
 !-----------------------------------------------------------------------------
       subroutine eback_sc_corr(esccor)
 ! 7/21/2007 Correlations between the backbone-local and side-chain-local
@@ -8411,9 +8934,9 @@
 
       iti1 = itortyp(itype(i+1,1))
       if (j.lt.nres-1) then
-        itj1 = itortyp(itype(j+1,1))
+        itj1 = itype2loc(itype(j+1,1))
       else
-        itj1=ntortyp+1
+        itj1=nloctyp
       endif
       do iii=1,2
         dipi(iii,1)=Ub2(iii,i)
@@ -10685,6 +11208,7 @@
 #ifdef TIMING
       time01=MPI_Wtime()
 #endif
+!#define DEBUG
 #ifdef DEBUG
       write (iout,*) "sum_gradient gvdwc, gvdwx"
       do i=1,nres
@@ -10718,18 +11242,18 @@
 !        write (iout,'(i5,3f10.5,2x,f10.5)') 
 !     &  i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i)
 !      enddo
-      write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
-      do i=1,nres
-        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
-         i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
-         (gvdwc_scpp(j,i),j=1,3)
-      enddo
-      write (iout,*) "gelc_long gvdwpp gel_loc_long"
-      do i=1,nres
-        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
-         i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
-         (gelc_loc_long(j,i),j=1,3)
-      enddo
+!      write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
+!      do i=1,nres
+!        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+!         i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
+!         (gvdwc_scpp(j,i),j=1,3)
+!      enddo
+!      write (iout,*) "gelc_long gvdwpp gel_loc_long"
+!      do i=1,nres
+!        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+!         i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
+!         (gelc_loc_long(j,i),j=1,3)
+!      enddo
       call flush(iout)
 #endif
 #ifdef SPLITELE
@@ -12240,7 +12764,7 @@
 !      call intcartderiv
 !      call checkintcartgrad
       call zerograd
-      aincr=2.0D-5
+      aincr=1.0D-7
       write(iout,*) 'Calling CHECK_ECARTINT.',aincr
       nf=0
       icall=0
@@ -16123,20 +16647,66 @@
 !
 ! Calculate the virtual-bond-angle energy.
 !
-      call ebend(ebe,ethetacnstr)
-!
 ! Calculate the SC local energy.
 !
       call vec_and_deriv
       call esc(escloc)
 !
+      if (wang.gt.0d0) then
+       if (tor_mode.eq.0) then
+         call ebend(ebe)
+       else
+!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+         call ebend_kcc(ebe)
+       endif
+      else
+        ebe=0.0d0
+      endif
+      ethetacnstr=0.0d0
+      if (with_theta_constr) call etheta_constr(ethetacnstr)
+
+!       write(iout,*) "in etotal afer ebe",ipot
+
+!      print *,"Processor",myrank," computed UB"
+!
+! Calculate the SC local energy.
+!
+      call esc(escloc)
+!elwrite(iout,*) "in etotal afer esc",ipot
+!      print *,"Processor",myrank," computed USC"
+!
+! Calculate the virtual-bond torsional energy.
+!
+!d    print *,'nterm=',nterm
+!      if (wtor.gt.0) then
+!       call etor(etors,edihcnstr)
+!      else
+!       etors=0
+!       edihcnstr=0
+!      endif
+      if (wtor.gt.0.0d0) then
+         if (tor_mode.eq.0) then
+           call etor(etors)
+         else
+!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+           call etor_kcc(etors)
+         endif
+      else
+        etors=0.0d0
+      endif
+      edihcnstr=0.0d0
+      if (ndih_constr.gt.0) call etor_constr(edihcnstr)
+
 ! Calculate the virtual-bond torsional energy.
 !
-      call etor(etors,edihcnstr)
 !
 ! 6/23/01 Calculate double-torsional energy
 !
+      if ((wtor_d.gt.0.0d0).and.(tor_mode.eq.0)) then
       call etor_d(etors_d)
+      endif
 !
 ! 21/5/07 Calculate local sicdechain correlation energy
 !
@@ -19958,6 +20528,24 @@
       allocate(Ug2DtEUgder(2,2,2,nres))
       allocate(DtUg2EUgder(2,2,2,nres))
 !(2,2,2,maxres)
+      allocate(b1(2,nres))      !(2,-maxtor:maxtor)
+      allocate(b2(2,nres))      !(2,-maxtor:maxtor)
+      allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
+      allocate(b2tilde(2,nres)) !(2,-maxtor:maxtor)
+
+      allocate(ctilde(2,2,nres))
+      allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
+      allocate(gtb1(2,nres))
+      allocate(gtb2(2,nres))
+      allocate(cc(2,2,nres))
+      allocate(dd(2,2,nres))
+      allocate(ee(2,2,nres))
+      allocate(gtcc(2,2,nres))
+      allocate(gtdd(2,2,nres))
+      allocate(gtee(2,2,nres))
+      allocate(gUb2(2,nres))
+      allocate(gteUg(2,2,nres))
+
 !      common /rotat_old/
       allocate(costab(nres))
       allocate(sintab(nres))
@@ -22081,7 +22669,8 @@
         rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt,      &
         opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
         opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
-        Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip
+        Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+        ndiv,ndivi
         real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
         gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
         dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
@@ -22095,17 +22684,6 @@
         ecation_prot=0.0d0
 ! first lets calculate interaction with peptide groups
         if (nres_molec(5).eq.0) return
-         wconst=78
-        wdip =1.092777950857032D2
-        wdip=wdip/wconst
-        wmodquad=-2.174122713004870D4
-        wmodquad=wmodquad/wconst
-        wquad1 = 3.901232068562804D1
-        wquad1=wquad1/wconst
-        wquad2 = 3
-        wquad2=wquad2/wconst
-        wvan1 = 0.1
-        wvan2 = 6
         itmp=0
         do i=1,4
         itmp=itmp+nres_molec(i)
@@ -22125,6 +22703,23 @@
           if (zi.lt.0) zi=zi+boxzsize
 
          do j=itmp+1,itmp+nres_molec(5)
+!           print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+         wconst=78*ndiv
+        wdip =1.092777950857032D2
+        wdip=wdip/wconst
+        wmodquad=-2.174122713004870D4
+        wmodquad=wmodquad/wconst
+        wquad1 = 3.901232068562804D1
+        wquad1=wquad1/wconst
+        wquad2 = 3
+        wquad2=wquad2/wconst
+        wvan1 = 0.1
+        wvan2 = 6
+!        itmp=0
+
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
@@ -22198,6 +22793,7 @@
         E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
              2*wvan2**6*Irsistp)
         ecation_prot = ecation_prot+E1+E2
+!        print *,"ecatprot",i,j,ecation_prot,rcpm
         dE1dr = -2*costhet*wdip*Irthrp-& 
          (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
         dE2dr = 3*wquad1*wquad2*Irfourp-     &
@@ -22237,6 +22833,9 @@
           enddo
            cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
          do j=itmp+1,itmp+nres_molec(5)
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
@@ -22279,7 +22878,10 @@
        endif
 !       enddo
 !       enddo
-         if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+! 15- Glu 16-Asp
+         if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.&
+         ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.&
+         (itype(i,1).eq.25))) then
             if(itype(i,1).eq.16) then
             inum=1
             else
@@ -22324,7 +22926,15 @@
 
 !  The weights of the energy function calculated from
 !The quantum mechanical GAMESS simulations of calcium with ASP/GLU
-        wh2o=78
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            ndivi=0.5
+          else
+            ndivi=1.0
+          endif
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+        wh2o=78*ndivi*ndiv
         wc = vcatprm(1)
         wc=wc/wh2o
         wdip =vcatprm(2)
@@ -22493,7 +23103,10 @@
         v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
 !  The weights of the energy function calculated from
 !The quantum mechanical GAMESS simulations of ASN/GLN with calcium
-        wh2o=78
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+        wh2o=78*ndiv
         wdip =vcatprm(2)
         wdip=wdip/wh2o
         wquad1 =vcatprm(3)
@@ -24761,6 +25374,16 @@
 !c! Compute head-head and head-tail energies for each state
 
           isel = iabs(Qi) + iabs(Qj)
+! double charge for Phophorylated! itype - 25,27,27
+!          if ((itype(i).eq.27).or.(itype(i).eq.26).or.(itype(i).eq.25)) then
+!            Qi=Qi*2
+!            Qij=Qij*2
+!           endif
+!          if ((itype(j).eq.27).or.(itype(j).eq.26).or.(itype(j).eq.25)) then
+!            Qj=Qj*2
+!            Qij=Qij*2
+!           endif
+
 !          isel=0
           IF (isel.eq.0) THEN
 !c! No charges - do nothing
@@ -24774,24 +25397,59 @@
 
           ELSE IF (isel.eq.1 .and. iabs(Qi).eq.1) THEN
 !c! Charge-nonpolar interactions
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            Qi=Qi*2
+            Qij=Qij*2
+           endif
+          if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+            Qj=Qj*2
+            Qij=Qij*2
+           endif
+
            CALL eqn(epol)
            eheadtail = epol
 !           eheadtail = 0.0d0
 
           ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN
 !c! Nonpolar-charge interactions
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            Qi=Qi*2
+            Qij=Qij*2
+           endif
+          if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+            Qj=Qj*2
+            Qij=Qij*2
+           endif
+
            CALL enq(epol)
            eheadtail = epol
 !           eheadtail = 0.0d0
 
           ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN
 !c! Charge-dipole interactions
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            Qi=Qi*2
+            Qij=Qij*2
+           endif
+          if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+            Qj=Qj*2
+            Qij=Qij*2
+           endif
+
            CALL eqd(ecl, elj, epol)
            eheadtail = ECL + elj + epol
 !           eheadtail = 0.0d0
 
           ELSE IF (isel.eq.3 .and. icharge(itypi).eq.2) THEN
 !c! Dipole-charge interactions
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            Qi=Qi*2
+            Qij=Qij*2
+           endif
+          if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+            Qj=Qj*2
+            Qij=Qij*2
+           endif
            CALL edq(ecl, elj, epol)
           eheadtail = ECL + elj + epol
 !           eheadtail = 0.0d0
@@ -24800,6 +25458,15 @@
                iabs(Qi).eq.1).and.  &
                nstate(itypi,itypj).eq.1) THEN
 !c! Same charge-charge interaction ( +/+ or -/- )
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            Qi=Qi*2
+            Qij=Qij*2
+           endif
+          if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+            Qj=Qj*2
+            Qij=Qij*2
+           endif
+
            CALL eqq(Ecl,Egb,Epol,Fisocav,Elj)
            eheadtail = ECL + Egb + Epol + Fisocav + Elj
 !           eheadtail = 0.0d0
@@ -24808,6 +25475,15 @@
                iabs(Qi).eq.1).and. &
                nstate(itypi,itypj).ne.1) THEN
 !c! Different charge-charge interaction ( +/- or -/+ )
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            Qi=Qi*2
+            Qij=Qij*2
+           endif
+          if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+            Qj=Qj*2
+            Qij=Qij*2
+           endif
+
            CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
           END IF
        END IF  ! this endif ends the "catch the gly-gly" at the beggining of Fcav
@@ -24838,7 +25514,7 @@
       use calc_data
       use comm_momo
        real (kind=8) ::  facd3, facd4, federmaus, adler,&
-         Ecl,Egb,Epol,Fisocav,Elj,Fgb
+         Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
 !       integer :: k
 !c! Epol and Gpol analytical parameters
        alphapol1 = alphapol(itypi,itypj)
@@ -24883,10 +25559,15 @@
        dGCLdOM12 = 0.0d0
        ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
        Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
-       Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb
+       debkap=debaykap(itypi,itypj)
+       Egb = -(332.0d0 * Qij *&
+        (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
 !       print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
 !c! Derivative of Egb is Ggb...
-       dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb)
+       dGGBdFGB = -(-332.0d0 * Qij * &
+       (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+       -(332.0d0 * Qij *&
+        (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
        dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
        dGGBdR = dGGBdFGB * dFGBdR
 !c!-------------------------------------------------------------------
@@ -25911,4 +26592,50 @@
        dPOLdOM2 = 0.0d0
        RETURN
       END SUBROUTINE elgrad_init
+
+      double precision function tschebyshev(m,n,x,y)
+      implicit none
+      integer i,m,n
+      double precision x(n),y,yy(0:maxvar),aux
+!c Tschebyshev polynomial. Note that the first term is omitted 
+!c m=0: the constant term is included
+!c m=1: the constant term is not included
+      yy(0)=1.0d0
+      yy(1)=y
+      do i=2,n
+        yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+      enddo
+      aux=0.0d0
+      do i=m,n
+        aux=aux+x(i)*yy(i)
+      enddo
+      tschebyshev=aux
+      return
+      end function tschebyshev
+!C--------------------------------------------------------------------------
+      double precision function gradtschebyshev(m,n,x,y)
+      implicit none
+      integer i,m,n
+      double precision x(n+1),y,yy(0:maxvar),aux
+!c Tschebyshev polynomial. Note that the first term is omitted
+!c m=0: the constant term is included
+!c m=1: the constant term is not included
+      yy(0)=1.0d0
+      yy(1)=2.0d0*y
+      do i=2,n
+        yy(i)=2*y*yy(i-1)-yy(i-2)
+      enddo
+      aux=0.0d0
+      do i=m,n
+        aux=aux+x(i+1)*yy(i)*(i+1)
+!C        print *, x(i+1),yy(i),i
+      enddo
+      gradtschebyshev=aux
+      return
+      end function gradtschebyshev
+
+
+
+
+
       end module energy