X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=6c60e98dfc180682913be231305b7686444a1ffe;hb=aa10e88041f16954d81953a1ecfbfc956e089f72;hp=6727cd8655a84a47b1151f8d8eeedb118d168b56;hpb=b657faa3790a5b86aeeaa6f88aa849df89121b23;p=unres4.git

diff --git a/source/unres/energy.F90 b/source/unres/energy.F90
index 6727cd8..6c60e98 100644
--- a/source/unres/energy.F90
+++ b/source/unres/energy.F90
@@ -1,4 +1,4 @@
-      module energy
+	      module energy
 !-----------------------------------------------------------------------------
       use io_units
       use names
@@ -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
@@ -377,6 +378,13 @@
           wscbase=weights(46)
           wscpho=weights(47)
           wpeppho=weights(48)
+!      welpsb=weights(28)*fact(1)
+!
+!      wcorr_nucl= weights(37)*fact(1)
+!     wcorr3_nucl=weights(38)*fact(2)
+!     wtor_nucl=  weights(35)*fact(1)
+!     wtor_d_nucl=weights(36)*fact(2)
+
         endif
         time_Bcast=time_Bcast+MPI_Wtime()-time00
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
@@ -609,9 +617,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 +670,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 +703,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"
        
 !
@@ -763,6 +800,8 @@
         call AFMforce(Eafmforce)
       else if (selfguide.gt.0) then
         call AFMvel(Eafmforce)
+      else
+        Eafmforce=0.0d0
       endif
       endif
       if (tubemode.eq.1) then
@@ -801,6 +840,7 @@
        eespp=0.0d0
       endif
 !      write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
+!      print *,"before ecatcat",wcatcat
       if (nfgtasks.gt.1) then
       if (fg_rank.eq.0) then
       call ecatcat(ecationcation)
@@ -821,7 +861,7 @@
       epeppho=0.0
       endif
 !      call ecatcat(ecationcation)
-!      print *,"after ebend", ebe_nucl
+!      print *,"after ebend", wtor_nucl 
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
 #endif
@@ -885,8 +925,8 @@
 !    Here are the energies showed per procesor if the are more processors 
 !    per molecule then we sum it up in sum_energy subroutine 
 !      print *," Processor",myrank," calls SUM_ENERGY"
-      energia(41)=ecation_prot
-      energia(42)=ecationcation
+      energia(42)=ecation_prot
+      energia(41)=ecationcation
       energia(46)=escbase
       energia(47)=epepbase
       energia(48)=escpho
@@ -1011,8 +1051,8 @@
       etors_d_nucl=energia(36)
       ecorr_nucl=energia(37)
       ecorr3_nucl=energia(38)
-      ecation_prot=energia(41)
-      ecationcation=energia(42)
+      ecation_prot=energia(42)
+      ecationcation=energia(41)
       escbase=energia(46)
       epepbase=energia(47)
       escpho=energia(48)
@@ -1156,6 +1196,11 @@
       wtor=weights(13)*fact(1)
       wtor_d=weights(14)*fact(2)
       wsccor=weights(21)*fact(1)
+      welpsb=weights(28)*fact(1)
+      wcorr_nucl= weights(37)*fact(1)
+      wcorr3_nucl=weights(38)*fact(2)
+      wtor_nucl=  weights(35)*fact(1)
+      wtor_d_nucl=weights(36)*fact(2)
 
       return
       end subroutine rescale_weights
@@ -1224,8 +1269,8 @@
       etors_d_nucl=energia(36)
       ecorr_nucl=energia(37)
       ecorr3_nucl=energia(38)
-      ecation_prot=energia(41)
-      ecationcation=energia(42)
+      ecation_prot=energia(42)
+      ecationcation=energia(41)
       escbase=energia(46)
       epepbase=energia(47)
       escpho=energia(48)
@@ -1300,11 +1345,11 @@
         ecorr,wcorr,&
         ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
         eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,&
-        ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,     &
+        ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforce,     &
         etube,wtube, &
         estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,&
-        evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb&
-        evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl&
+        evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,&
+        evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,&
         etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
         ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat,  &
         escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
@@ -1982,8 +2027,8 @@
 !          write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
-            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
-            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_cut=sscale_ele(1.0d0/(rij))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij))
 !            print *,sss_ele_cut,sss_ele_grad,&
 !            1.0d0/(rij),r_cut_ele,rlamb_ele
             if (sss_ele_cut.le.0.0) cycle
@@ -2045,7 +2090,7 @@
             fac=rij*fac
 !            print *,'before fac',fac,rij,evdwij
             fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
-            /sigma(itypi,itypj)*rij
+            *rij
 !            print *,'grad part scale',fac,   &
 !             evdwij*sss_ele_grad/sss_ele_cut &
 !            /sigma(itypi,itypj)*rij
@@ -2689,18 +2734,174 @@
 !      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
+!         write(iout,*) "i,",molnum(i)
+!         print *, "i,",molnum(i),i,itype(i-2,1)
+        if (molnum(i).eq.1) then
+          iti = itype2loc(itype(i-2,1))
+        else
+          iti=nloctyp
+        endif
+        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 +2974,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,44 +3025,44 @@
             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
 !        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
           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 +3683,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 +4125,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 +4360,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 +4773,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 +4835,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 +4858,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 +5007,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 +5026,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 +5074,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 +5184,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 +5207,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 +6251,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 +6457,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 +7345,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 +7427,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 +7592,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 +7689,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 +8954,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 +11228,7 @@
 #ifdef TIMING
       time01=MPI_Wtime()
 #endif
+!#define DEBUG
 #ifdef DEBUG
       write (iout,*) "sum_gradient gvdwc, gvdwx"
       do i=1,nres
@@ -10718,18 +11262,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
@@ -10797,7 +11341,7 @@
                      +wcorr3_nucl*gradcorr3_nucl(j,i) +&
                      wcatprot* gradpepcat(j,i)+ &
                      wcatcat*gradcatcat(j,i)+   &
-                     wscbase*gvdwc_scbase(j,i)  &
+                     wscbase*gvdwc_scbase(j,i)+ &
                      wpepbase*gvdwc_pepbase(j,i)+&
                      wscpho*gvdwc_scpho(j,i)+&
                      wpeppho*gvdwc_peppho(j,i)
@@ -11028,7 +11572,7 @@
                      +gradafm(j,i) &
                      +wliptran*gliptranc(j,i) &
                      +welec*gshieldc(j,i) &
-                     +welec*gshieldc_loc(j,) &
+                     +welec*gshieldc_loc(j,i) &
                      +wcorr*gshieldc_ec(j,i) &
                      +wcorr*gshieldc_loc_ec(j,i) &
                      +wturn3*gshieldc_t3(j,i) &
@@ -12240,7 +12784,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
@@ -13496,8 +14040,8 @@
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
-            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_cut=sscale_ele(1.0d0/(rij))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij))
             sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
             if (sss_ele_cut.le.0.0) cycle
             if (sss.lt.1.0d0) then
@@ -13553,7 +14097,7 @@
               sigder=fac*sigder
               fac=rij*fac
               fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
-            /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij  &
+              *rij-sss_grad/(1.0-sss)*rij  &
             /sigmaii(itypi,itypj))
 !              fac=0.0d0
 ! Calculate the radial part of the gradient
@@ -13787,8 +14331,8 @@
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
             sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
-            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
-            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_cut=sscale_ele(1.0d0/(rij))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij))
             if (sss_ele_cut.le.0.0) cycle
 
             if (sss.gt.0.0d0) then
@@ -13844,7 +14388,7 @@
               sigder=fac*sigder
               fac=rij*fac
               fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
-            /sigma(itypi,itypj)*rij+sss_grad/sss*rij  &
+            *rij+sss_grad/sss*rij  &
             /sigmaii(itypi,itypj))
 
 !              fac=0.0d0
@@ -16123,20 +16667,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 +20548,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 +22689,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 +22704,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 +22723,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 +22813,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 +22853,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 +22898,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
@@ -22289,12 +22911,29 @@
             vcatprm(k)=catprm(k,inum)
             enddo
             dASGL=catprm(7,inum)
-             do k=1,3
-                vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
-                valpha(k)=c(k,i)
-                vcat(k)=c(k,j)
-              enddo
-                      do k=1,3
+!             do k=1,3
+!                vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+                vcm(1)=(cm1(1)/cm1mag)*dASGL+xi
+                vcm(2)=(cm1(2)/cm1mag)*dASGL+yi
+                vcm(3)=(cm1(3)/cm1mag)*dASGL+zi
+
+!                valpha(k)=c(k,i)
+!                vcat(k)=c(k,j)
+                if (subchap.eq.1) then
+                 vcat(1)=xj_temp
+                 vcat(2)=yj_temp
+                 vcat(3)=zj_temp
+                 else
+                vcat(1)=xj_safe
+                vcat(2)=yj_safe
+                vcat(3)=zj_safe
+                 endif
+                valpha(1)=xi-c(1,i+nres)+c(1,i)
+                valpha(2)=yi-c(2,i+nres)+c(2,i)
+                valpha(3)=zi-c(3,i+nres)+c(3,i)
+
+!              enddo
+        do k=1,3
           dx(k) = vcat(k)-vcm(k)
         enddo
         do k=1,3
@@ -22307,7 +22946,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)
@@ -22316,7 +22963,7 @@
         wquad1=wquad1/wh2o
         wquad2 = vcatprm(4)
         wquad2=wquad2/wh2o
-        wquad2p = 1-wquad2
+        wquad2p = 1.0d0-wquad2
         wvan1 = vcatprm(5)
         wvan2 =vcatprm(6)
         opt = dx(1)**2+dx(2)**2
@@ -22327,11 +22974,11 @@
         rsixp = rfourp*rsecp
         reight=rsixp*rsecp
         Ir = 1.0d0/rs
-        Irsecp = 1/rsecp
+        Irsecp = 1.0d0/rsecp
         Irthrp = Irsecp/rs
         Irfourp = Irthrp/rs
-        Irsixp = 1/rsixp
-        Ireight=1/reight
+        Irsixp = 1.0d0/rsixp
+        Ireight=1.0d0/reight
         Irtw=Irsixp*Irsixp
         Irthir=Irtw/rs
         Irfourt=Irthir/rs
@@ -22442,11 +23089,27 @@
             vcatprm(k)=catprm(k,inum)
             enddo
             dASGL=catprm(7,inum)
-             do k=1,3
-                vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
-                valpha(k)=c(k,i)
-                vcat(k)=c(k,j)
-              enddo
+!             do k=1,3
+!                vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+!                valpha(k)=c(k,i)
+!                vcat(k)=c(k,j)
+!              enddo
+                vcm(1)=(cm1(1)/cm1mag)*dASGL+xi
+                vcm(2)=(cm1(2)/cm1mag)*dASGL+yi
+                vcm(3)=(cm1(3)/cm1mag)*dASGL+zi
+                if (subchap.eq.1) then
+                 vcat(1)=xj_temp
+                 vcat(2)=yj_temp
+                 vcat(3)=zj_temp
+                 else
+                vcat(1)=xj_safe
+                vcat(2)=yj_safe
+                vcat(3)=zj_safe
+                endif
+                valpha(1)=xi-c(1,i+nres)+c(1,i)
+                valpha(2)=yi-c(2,i+nres)+c(2,i)
+                valpha(3)=zi-c(3,i+nres)+c(3,i)
+
 
         do k=1,3
           dx(k) = vcat(k)-vcm(k)
@@ -22460,7 +23123,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)
@@ -22555,6 +23221,11 @@
             dscmag = 0.0d0
             do k=1,3
               dscvec(k) = c(k,i+nres)-c(k,i)
+! TU SPRAWDZ???
+!              dscvec(1) = xj
+!              dscvec(2) = yj
+!              dscvec(3) = zj
+
               dscmag = dscmag+dscvec(k)*dscvec(k)
             enddo
             dscmag3 = dscmag
@@ -22576,7 +23247,10 @@
            else
             rcal = 0.0d0
             do k=1,3
-              r(k) = c(k,j)-c(k,i+nres)
+!              r(k) = c(k,j)-c(k,i+nres)
+              r(1) = xj
+              r(2) = yj
+              r(3) = zj
               rcal = rcal+r(k)*r(k)
             enddo
             ract=sqrt(rcal)
@@ -24720,6 +25394,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
@@ -24733,24 +25417,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
@@ -24759,6 +25478,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
@@ -24767,6 +25495,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
@@ -24797,7 +25534,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)
@@ -24842,10 +25579,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!-------------------------------------------------------------------
@@ -25870,4 +26612,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