AFM+tube+kcc
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index 215fb6c..21c0197 100644 (file)
@@ -55,6 +55,8 @@ C FG slaves as WEIGHTS array.
           weights_(17)=wbond
           weights_(18)=scal14
           weights_(21)=wsccor
+          weights_(22)=wtube
+
 C FG Master broadcasts the WEIGHTS_ array
           call MPI_Bcast(weights_(1),n_ene,
      &        MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
@@ -81,6 +83,7 @@ C FG slaves receive the WEIGHTS array
           wbond=weights(17)
           scal14=weights(18)
           wsccor=weights(21)
+          wtube=weights(22)
         endif
         time_Bcast=time_Bcast+MPI_Wtime()-time00
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
@@ -142,8 +145,10 @@ C the shielding factor is set this factor is describing how each
 C peptide group is shielded by side-chains
 C the matrix - shield_fac(i) the i index describe the ith between i and i+1
 C      write (iout,*) "shield_mode",shield_mode
-      if (shield_mode.gt.0) then
+      if (shield_mode.eq.1) then
        call set_shield_fac
+      else if  (shield_mode.eq.2) then
+       call set_shield_fac2
       endif
 c      print *,"Processor",myrank," left VEC_AND_DERIV"
       if (ipot.lt.6) then
@@ -302,6 +307,15 @@ C      print *,"za lipidami"
       else if (selfguide.gt.0) then
         call AFMvel(Eafmforce)
       endif
+      if (TUBElog.eq.1) then
+C      print *,"just before call"
+        call calctube(Etube)
+       elseif (TUBElog.eq.2) then
+        call calctube2(Etube)
+       else
+       Etube=0.0d0
+       endif
+
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
 #endif
@@ -346,6 +360,7 @@ C
       energia(22)=eliptran
       energia(23)=Eafmforce
       energia(24)=ethetacnstr
+      energia(25)=Etube
 c    Here are the energies showed per procesor if the are more processors 
 c    per molecule then we sum it up in sum_energy subroutine 
 c      print *," Processor",myrank," calls SUM_ENERGY"
@@ -440,6 +455,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       eliptran=energia(22)
       Eafmforce=energia(23)
       ethetacnstr=energia(24)
+      Etube=energia(25)
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
@@ -447,7 +463,7 @@ cMS$ATTRIBUTES C ::  proc_proc
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
      & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
-     & +ethetacnstr
+     & +ethetacnstr+wtube*Etube
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
      & +wang*ebe+wtor*etors+wscloc*escloc
@@ -456,7 +472,7 @@ cMS$ATTRIBUTES C ::  proc_proc
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
      & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
      & +Eafmforce
-     & +ethetacnstr
+     & +ethetacnstr+wtube*Etube
 #endif
       energia(0)=etot
 c detecting NaNQ
@@ -566,6 +582,8 @@ c      enddo
      &                 +wturn3*gshieldc_t3(j,i)
      &                 +wturn4*gshieldc_t4(j,i)
      &                 +wel_loc*gshieldc_ll(j,i)
+     &                +wtube*gg_tube(j,i)
+
 
 
         enddo
@@ -589,6 +607,8 @@ c      enddo
      &                 +wcorr*gshieldc_ec(j,i)
      &                 +wturn4*gshieldc_t4(j,i)
      &                 +wel_loc*gshieldc_ll(j,i)
+     &                +wtube*gg_tube(j,i)
+
 
 
         enddo
@@ -745,11 +765,7 @@ C          print *,gradafm(1,13),"AFM"
      &                 +wturn4*gshieldc_loc_t4(j,i)
      &                 +wel_loc*gshieldc_ll(j,i)
      &                 +wel_loc*gshieldc_loc_ll(j,i)
-
-
-
-
-
+     &                +wtube*gg_tube(j,i)
 
 #else
           gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
@@ -782,9 +798,7 @@ C          print *,gradafm(1,13),"AFM"
      &                 +wturn4*gshieldc_loc_t4(j,i)
      &                 +wel_loc*gshieldc_ll(j,i)
      &                 +wel_loc*gshieldc_loc_ll(j,i)
-
-
-
+     &                +wtube*gg_tube(j,i)
 
 
 #endif
@@ -799,6 +813,7 @@ C          print *,gradafm(1,13),"AFM"
      &                 +wturn3*gshieldx_t3(j,i)
      &                 +wturn4*gshieldx_t4(j,i)
      &                 +wel_loc*gshieldx_ll(j,i)
+     &                 +wtube*gg_tube_sc(j,i)
 
 
 
@@ -1004,6 +1019,7 @@ c-------------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.SBRIDGE'
+      include 'COMMON.CONTROL'
       double precision kfac /2.4d0/
       double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/
 c      facT=temp0/t_bath
@@ -1039,6 +1055,11 @@ c      facT=2*temp0/(t_bath+temp0)
 #endif
        stop 555
       endif
+      if (shield_mode.gt.0) then
+       wscp=weights(2)*fact
+       wsc=weights(1)*fact
+       wvdwpp=weights(16)*fact
+      endif
       welec=weights(3)*fact
       wcorr=weights(4)*fact3
       wcorr5=weights(5)*fact4
@@ -1093,14 +1114,16 @@ C------------------------------------------------------------------------
       eliptran=energia(22)
       Eafmforce=energia(23) 
       ethetacnstr=energia(24)
+      etube=energia(25)
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,
      &  escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
      &  ecorr,wcorr,
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
-     &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
+     &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,
      &  ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
+     &  etube,wtube,
      &  etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
@@ -1128,6 +1151,7 @@ C------------------------------------------------------------------------
      & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
      & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
      & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
+     & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/
      & 'ETOT=  ',1pE16.6,' (total)')
 
 #else
@@ -1138,6 +1162,7 @@ C------------------------------------------------------------------------
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
      &  ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
+     &  etube,wtube,
      &  etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
@@ -1164,6 +1189,7 @@ C------------------------------------------------------------------------
      & 'UCONST=',1pE16.6,' (Constraint energy)'/ 
      & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
      & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
+     & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/
      & 'ETOT=  ',1pE16.6,' (total)')
 #endif
       return
@@ -1772,6 +1798,7 @@ C lipbufthick is thickenes of lipid buffore
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+C      write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj)
 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)')
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
 C      if (ssgradlipj.gt.0.0d0) print *,"??WTF??"
@@ -2048,6 +2075,7 @@ C lipbufthick is thickenes of lipid buffore
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5') 
 C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
+C      write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -2775,28 +2803,28 @@ c      write(iout,*) 'nphi=',nphi,nres
 #endif
 #ifdef NEWCORR
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
-          iti = itortyp(itype(i-2))
+          iti = itype2loc(itype(i-2))
         else
-          iti=ntortyp+1
+          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 = itortyp(itype(i-1))
+          iti1 = itype2loc(itype(i-1))
         else
-          iti1=ntortyp+1
+          iti1=nloctyp
         endif
 c        write(iout,*),i
-        b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0)
+        b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0d0)
      &           +bnew1(2,1,iti)*dsin(theta(i-1))
-     &           +bnew1(3,1,iti)*dcos(theta(i-1)/2.0)
+     &           +bnew1(3,1,iti)*dcos(theta(i-1)/2.0d0)
         gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
      &             +bnew1(2,1,iti)*dcos(theta(i-1))
      &             -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
 c     &           +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
 c     &*(cos(theta(i)/2.0)
-        b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0)
+        b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0d0)
      &           +bnew2(2,1,iti)*dsin(theta(i-1))
-     &           +bnew2(3,1,iti)*dcos(theta(i-1)/2.0)
+     &           +bnew2(3,1,iti)*dcos(theta(i-1)/2.0d0)
 c     &           +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
 c     &*(cos(theta(i)/2.0)
         gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
@@ -2835,15 +2863,15 @@ c       write (iout,*) 'theta=', theta(i-1)
        enddo
 #else
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
-          iti = itortyp(itype(i-2))
+          iti = itype2loc(itype(i-2))
         else
-          iti=ntortyp+1
+          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 = itortyp(itype(i-1))
+          iti1 = itype2loc(itype(i-1))
         else
-          iti1=ntortyp+1
+          iti1=nloctyp
         endif
         b1(1,i-2)=b(3,iti)
         b1(2,i-2)=b(5,iti)
@@ -2932,15 +2960,15 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         endif
 c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
-          iti = itortyp(itype(i-2))
+          iti = itype2loc(itype(i-2))
         else
-          iti=ntortyp
+          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 = itortyp(itype(i-1))
+          iti1 = itype2loc(itype(i-1))
         else
-          iti1=ntortyp
+          iti1=nloctyp
         endif
 cd        write (iout,*) '*******i',i,' iti1',iti
 cd        write (iout,*) 'b1',b1(:,iti)
@@ -2953,8 +2981,8 @@ c        if (i .gt. iatel_s+2) then
           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
-c          write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti),
-c     &    EE(1,2,iti),EE(2,2,iti)
+c          write(iout,*) "co jest kurwa", iti, EE(1,1,i),EE(2,1,i),
+c     &    EE(1,2,iti),EE(2,2,i)
           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))
 c          write(iout,*) "Macierz EUG",
@@ -2989,18 +3017,24 @@ c     &    eug(2,2,i-2)
 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
           if (itype(i-1).le.ntyp) then
-            iti1 = itortyp(itype(i-1))
+            iti1 = itype2loc(itype(i-1))
           else
-            iti1=ntortyp
+            iti1=nloctyp
           endif
         else
-          iti1=ntortyp
+          iti1=nloctyp
         endif
         do k=1,2
           mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
         enddo
-C        write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
-c        write (iout,*) 'mu ',mu(:,i-2),i-2
+#ifdef MUOUT
+        write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1),
+     &     rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2),
+     &       -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2),
+     &       dsqrt(b2(1,i-1)**2+b2(2,i-1)**2)
+     &      +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2),
+     &      ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itype2loc(iti))
+#endif
 cd        write (iout,*) 'mu1',mu1(:,i-2)
 cd        write (iout,*) 'mu2',mu2(:,i-2)
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
@@ -3287,11 +3321,11 @@ c      endif
 #endif
 #endif
 cd      do i=1,nres
-cd        iti = itortyp(itype(i))
+cd        iti = itype2loc(itype(i))
 cd        write (iout,*) i
 cd        do j=1,2
 cd        write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)') 
-cd     &  (EE(j,k,iti),k=1,2),(Ug(j,k,i),k=1,2),(EUg(j,k,i),k=1,2)
+cd     &  (EE(j,k,i),k=1,2),(Ug(j,k,i),k=1,2),(EUg(j,k,i),k=1,2)
 cd        enddo
 cd      enddo
       return
@@ -3409,21 +3443,23 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
 C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
       do i=iturn3_start,iturn3_end
-        if (i.le.1) cycle
+c        if (i.le.1) cycle
 C        write(iout,*) "tu jest i",i
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((i+4).gt.nres)
-     & .or.((i-1).le.0)
+C Adam: Unnecessary: handled by iturn3_end and iturn3_start
+c     & .or.((i+4).gt.nres)
+c     & .or.((i-1).le.0)
 C end of changes by Ana
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i+3).eq.ntyp1) cycle
-        if(i.gt.1)then
-          if(itype(i-1).eq.ntyp1)cycle
-        end if
-        if(i.LT.nres-3)then
-          if (itype(i+4).eq.ntyp1) cycle
-        end if
+C Adam: Instructions below will switch off existing interactions
+c        if(i.gt.1)then
+c          if(itype(i-1).eq.ntyp1)cycle
+c        end if
+c        if(i.LT.nres-3)then
+c          if (itype(i+4).eq.ntyp1) cycle
+c        end if
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -3445,17 +3481,17 @@ C end of changes by Ana
         num_cont_hb(i)=num_conti
       enddo
       do i=iturn4_start,iturn4_end
-        if (i.le.1) cycle
+        if (i.lt.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((i+5).gt.nres)
-     & .or.((i-1).le.0)
+c     & .or.((i+5).gt.nres)
+c     & .or.((i-1).le.0)
 C end of changes suggested by Ana
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
-     &    .or. itype(i+5).eq.ntyp1
-     &    .or. itype(i).eq.ntyp1
-     &    .or. itype(i-1).eq.ntyp1
+c     &    .or. itype(i+5).eq.ntyp1
+c     &    .or. itype(i).eq.ntyp1
+c     &    .or. itype(i-1).eq.ntyp1
      &                             ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
 CTU KURWA
       do i=iatel_s,iatel_e
 C        do i=75,75
-        if (i.le.1) cycle
+c        if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((i+2).gt.nres)
-     & .or.((i-1).le.0)
+c     & .or.((i+2).gt.nres)
+c     & .or.((i-1).le.0)
 C end of changes by Ana
-     &  .or. itype(i+2).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
+c     &  .or. itype(i+2).eq.ntyp1
+c     &  .or. itype(i-1).eq.ntyp1
      &                ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -3575,14 +3611,14 @@ C I TU KURWA
         do j=ielstart(i),ielend(i)
 C          do j=16,17
 C          write (iout,*) i,j
-         if (j.le.1) cycle
+C         if (j.le.1) cycle
           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((j+2).gt.nres)
-     & .or.((j-1).le.0)
+c     & .or.((j+2).gt.nres)
+c     & .or.((j-1).le.0)
 C end of changes by Ana
-     & .or.itype(j+2).eq.ntyp1
-     & .or.itype(j-1).eq.ntyp1
+c     & .or.itype(j+2).eq.ntyp1
+c     & .or.itype(j-1).eq.ntyp1
      &) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
@@ -3643,6 +3679,7 @@ C 13-go grudnia roku pamietnego...
       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
      &                   0.0d0,1.0d0,0.0d0,
      &                   0.0d0,0.0d0,1.0d0/
+       integer xshift,yshift,zshift
 c          time00=MPI_Wtime()
 cd      write (iout,*) "eelecij",i,j
 c          ind=ind+1
@@ -3782,6 +3819,7 @@ cd     &      xmedi,ymedi,zmedi,xj,yj,zj
               write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') 
      &'evdw1',i,j,evdwij
      &,iteli,itelj,aaa,evdw1
+              write (iout,*) sss
               write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
      &fac_shield(i),fac_shield(j)
           endif
@@ -4662,6 +4700,7 @@ C        fac_shield(j)=0.6
      &  *fac_shield(i)*fac_shield(j)
         eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
+C#ifdef NEWCORR
 C Derivatives in theta
         gloc(nphi+i,icg)=gloc(nphi+i,icg)
      &  +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
@@ -4669,7 +4708,7 @@ C Derivatives in theta
         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
 
 C        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 C Derivatives in shield mode
@@ -4824,9 +4863,9 @@ c        write(iout,*)"WCHODZE W PROGRAM"
         a_temp(1,2)=a23
         a_temp(2,1)=a32
         a_temp(2,2)=a33
-        iti1=itortyp(itype(i+1))
-        iti2=itortyp(itype(i+2))
-        iti3=itortyp(itype(i+3))
+        iti1=itype2loc(itype(i+1))
+        iti2=itype2loc(itype(i+2))
+        iti3=itype2loc(itype(i+3))
 c        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))
@@ -5832,6 +5871,7 @@ C       Checking if it involves dummy (NH3+ or COO-) group
          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
 C YES   vbldpDUM is the equlibrium length of spring for Dummy atom
         diff = vbld(i)-vbldpDUM
+        if (energy_dec) write(iout,*) "dum_bond",i,diff 
          else
 C NO    vbldp0 is the equlibrium lenght of spring for peptide group
         diff = vbld(i)-vbldp0
@@ -5845,6 +5885,7 @@ C NO    vbldp0 is the equlibrium lenght of spring for peptide group
 c        write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
 c        endif
       enddo
+      
       estr=0.5d0*AKP*estr+estr1
 c
 c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
@@ -7465,11 +7506,12 @@ C The rigorous attempt to derive energy function
       include 'COMMON.TORCNSTR'
       include 'COMMON.CONTROL'
       logical lprn
-      double precision thybt1(maxtermkcc),thybt2(maxtermkcc)
+c      double precision thybt1(maxtermkcc),thybt2(maxtermkcc)
 C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
 C      print *,"wchodze kcc"
+      if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode
       if (tor_mode.ne.2) then
       etors=0.0D0
       endif
@@ -7489,60 +7531,86 @@ c     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
         sumnonchebyshev=0.0d0
         sumchebyshev=0.0d0
 C to avoid multiple devision by 2
-        theti22=0.5d0*theta(i)
+c        theti22=0.5d0*theta(i)
 C theta 12 is the theta_1 /2
 C theta 22 is theta_2 /2
-        theti12=0.5d0*theta(i-1)
+c        theti12=0.5d0*theta(i-1)
 C and appropriate sinus function
-        sinthet2=dsin(theta(i))
         sinthet1=dsin(theta(i-1))
+        sinthet2=dsin(theta(i))
         costhet1=dcos(theta(i-1))
         costhet2=dcos(theta(i))
+c Cosines of halves thetas
+        costheti12=0.5d0*(1.0d0+costhet1)
+        costheti22=0.5d0*(1.0d0+costhet2)
 C to speed up lets store its mutliplication
-         sint1t2=sinthet2*sinthet1        
+        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
+        etori=0.0d0
         do j=1,nterm_kcc(itori,itori1)
 
+          nval=nterm_kcc_Tb(itori,itori1)
           v1ij=v1_kcc(j,itori,itori1)
           v2ij=v2_kcc(j,itori,itori1)
+c          write (iout,*) "i",i," j",j," v1",v1ij," v2",v2ij
 C v1ij is c_n and d_n in euation above
           cosphi=dcos(j*phii)
           sinphi=dsin(j*phii)
-          sint1t2n=sint1t2**j
-          sumnonchebyshev=
-     &                    sint1t2n*(v1ij*cosphi+v2ij*sinphi)
-          actval=sint1t2n*(v1ij*cosphi+v2ij*sinphi)
+          sint1t2n1=sint1t2n
+          sint1t2n=sint1t2n*sint1t2
+          sumth1tyb1=tschebyshev(1,nval,v11_chyb(1,j,itori,itori1),
+     &        costheti12)
+          gradth1tyb1=-0.5d0*sinthet1*gradtschebyshev(0,nval-1,
+     &        v11_chyb(1,j,itori,itori1),costheti12)
+c          write (iout,*) "v11",(v11_chyb(k,j,itori,itori1),k=1,nval),
+c     &      " sumth1tyb1",sumth1tyb1," gradth1tyb1",gradth1tyb1
+          sumth2tyb1=tschebyshev(1,nval,v21_chyb(1,j,itori,itori1),
+     &        costheti22)
+          gradth2tyb1=-0.5d0*sinthet2*gradtschebyshev(0,nval-1,
+     &        v21_chyb(1,j,itori,itori1),costheti22)
+c          write (iout,*) "v21",(v21_chyb(k,j,itori,itori1),k=1,nval),
+c     &      " sumth2tyb1",sumth2tyb1," gradth2tyb1",gradth2tyb1
+          sumth1tyb2=tschebyshev(1,nval,v12_chyb(1,j,itori,itori1),
+     &        costheti12)
+          gradth1tyb2=-0.5d0*sinthet1*gradtschebyshev(0,nval-1,
+     &        v12_chyb(1,j,itori,itori1),costheti12)
+c          write (iout,*) "v12",(v12_chyb(k,j,itori,itori1),k=1,nval),
+c     &      " sumth1tyb2",sumth1tyb2," gradth1tyb2",gradth1tyb2
+          sumth2tyb2=tschebyshev(1,nval,v22_chyb(1,j,itori,itori1),
+     &        costheti22)
+          gradth2tyb2=-0.5d0*sinthet2*gradtschebyshev(0,nval-1,
+     &        v22_chyb(1,j,itori,itori1),costheti22)
+c          write (iout,*) "v22",(v22_chyb(k,j,itori,itori1),k=1,nval),
+c     &      " sumth2tyb2",sumth2tyb2," gradth2tyb2",gradth2tyb2
 C          etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi)
 C          if (energy_dec) etors_ii=etors_ii+
 C     &                v1ij*cosphi+v2ij*sinphi
 C glocig is the gradient local i site in gamma
-          glocig=j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n
+          actval1=v1ij*cosphi*(1.0d0+sumth1tyb1+sumth2tyb1)
+          actval2=v2ij*sinphi*(1.0d0+sumth1tyb2+sumth2tyb2)
+          etori=etori+sint1t2n*(actval1+actval2)
+          glocig=glocig+
+     &        j*sint1t2n*(v2ij*cosphi*(1.0d0+sumth1tyb2+sumth2tyb2)
+     &        -v1ij*sinphi*(1.0d0+sumth1tyb1+sumth2tyb1))
 C now gradient over theta_1
-          glocit1=actval/sinthet1*j*costhet1
-          glocit2=actval/sinthet2*j*costhet2
+          glocit1=glocit1+
+     &       j*sint1t2n1*costhet1*sinthet2*(actval1+actval2)+
+     &       sint1t2n*(v1ij*cosphi*gradth1tyb1+v2ij*sinphi*gradth1tyb2)
+          glocit2=glocit2+
+     &       j*sint1t2n1*sinthet1*costhet2*(actval1+actval2)+
+     &       sint1t2n*(v1ij*cosphi*gradth2tyb1+v2ij*sinphi*gradth2tyb2)
 
 C now the Czebyshev polinominal sum
-        do k=1,nterm_kcc_Tb(itori,itori1)
-         thybt1(k)=v1_chyb(k,j,itori,itori1)
-         thybt2(k)=v2_chyb(k,j,itori,itori1)
+c        do k=1,nterm_kcc_Tb(itori,itori1)
+c         thybt1(k)=v1_chyb(k,j,itori,itori1)
+c         thybt2(k)=v2_chyb(k,j,itori,itori1)
 C         thybt1(k)=0.0
 C         thybt2(k)=0.0
-        enddo 
-        sumth1thyb=tschebyshev
-     &         (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theti12)**2)
-        gradthybt1=gradtschebyshev
-     &         (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1),
-     &        dcos(theti12)**2)
-     & *dcos(theti12)*(-dsin(theti12))
-        sumth2thyb=tschebyshev
-     &         (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theti22)**2)
-        gradthybt2=gradtschebyshev
-     &         (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
-     &         dcos(theti22)**2)
-     & *dcos(theti22)*(-dsin(theti22))
+c        enddo 
 C        print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2,
 C     &         gradtschebyshev
 C     &         (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
@@ -7550,22 +7618,19 @@ C     &         dcos(theti22)**2),
 C     &         dsin(theti22)
 
 C now overal sumation
-         etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev
 C         print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb
+        enddo ! j
+        etors=etors+etori
 C derivative over gamma
-         gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
-     &   *(1.0d0+sumth1thyb+sumth2thyb)
+        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*(1.0d0+sumth1thyb+sumth2thyb)+
-     &   sumnonchebyshev*gradthybt1)
+        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*(1.0d0+sumth1thyb+sumth2thyb)+
-     &   sumnonchebyshev*gradthybt2)
-       enddo
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2
+        if (lprn) 
+     &    write (iout,*) i-2,i-1,itype(i-2),itype(i-1),itori,itori1,
+     &       theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori
       enddo
-     
 C        gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
 ! 6/20/98 - dihedral angle constraints
       if (tor_mode.ne.2) then
@@ -7614,7 +7679,8 @@ C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
 C      print *,"wchodze kcc"
-      if (tormode.ne.2) etheta=0.0D0
+      if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode
+      if (tor_mode.ne.2) etheta=0.0D0
       do i=ithet_start,ithet_end
 c        print *,i,itype(i-1),itype(i),itype(i-2)
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
@@ -7627,6 +7693,8 @@ c        print *,i,itype(i-1),itype(i),itype(i-2)
          enddo
          sumth1thyb=tschebyshev
      &         (1,nbend_kcc_Tb(iti),thybt1(1),costhet)
+        if (lprn) write (iout,*) i-1,itype(i-1),iti,theta(i)*rad2deg,
+     &    sumth1thyb
         ihelp=nbend_kcc_Tb(iti)-1
         gradthybt1=gradtschebyshev
      &         (0,ihelp,thybt1(1),costhet)
@@ -7635,7 +7703,7 @@ C        print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0)
         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
      &   gradthybt1*sinthet*(-0.5d0)
       enddo
-      if (tormode.ne.2) then
+      if (tor_mode.ne.2) then
       ethetacnstr=0.0d0
 C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
       do i=ithetaconstr_start,ithetaconstr_end
@@ -8677,7 +8745,7 @@ c     & ' eij',eij,' eesij',ees0pij,ees0mij,' and ',k,l
 c     & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' energy=',ekont*ees,
 c     & 'gradcorr_long'
 C Calculate the multi-body contribution to energy.
-c      ecorr=ecorr+ekont*ees
+C      ecorr=ecorr+ekont*ees
 C Calculate multi-body contributions to the gradient.
       coeffpees0pij=coeffp*ees0pij
       coeffmees0mij=coeffm*ees0mij
@@ -8831,9 +8899,9 @@ C---------------------------------------------------------------------------
      &  auxmat(2,2)
       iti1 = itortyp(itype(i+1))
       if (j.lt.nres-1) then
-        itj1 = itortyp(itype(j+1))
+        itj1 = itype2loc(itype(j+1))
       else
-        itj1=ntortyp
+        itj1=nloctyp
       endif
       do iii=1,2
         dipi(iii,1)=Ub2(iii,i)
@@ -8921,16 +8989,16 @@ cd      write (iout,*) "a_chujkl",((a_chuj(iii,jjj,kk,k),iii=1,2),jjj=1,2)
       if (l.eq.j+1) then
 C parallel orientation of the two CA-CA-CA frames.
         if (i.gt.1) then
-          iti=itortyp(itype(i))
+          iti=itype2loc(itype(i))
         else
-          iti=ntortyp
+          iti=nloctyp
         endif
-        itk1=itortyp(itype(k+1))
-        itj=itortyp(itype(j))
+        itk1=itype2loc(itype(k+1))
+        itj=itype2loc(itype(j))
         if (l.lt.nres-1) then
-          itl1=itortyp(itype(l+1))
+          itl1=itype2loc(itype(l+1))
         else
-          itl1=ntortyp
+          itl1=nloctyp
         endif
 C A1 kernel(j+1) A2T
 cd        do iii=1,2
@@ -9074,17 +9142,17 @@ C End vectors
       else
 C Antiparallel orientation of the two CA-CA-CA frames.
         if (i.gt.1) then
-          iti=itortyp(itype(i))
+          iti=itype2loc(itype(i))
         else
-          iti=ntortyp
+          iti=nloctyp
         endif
-        itk1=itortyp(itype(k+1))
-        itl=itortyp(itype(l))
-        itj=itortyp(itype(j))
+        itk1=itype2loc(itype(k+1))
+        itl=itype2loc(itype(l))
+        itj=itype2loc(itype(j))
         if (j.lt.nres-1) then
-          itj1=itortyp(itype(j+1))
+          itj1=itype2loc(itype(j+1))
         else 
-          itj1=ntortyp
+          itj1=nloctyp
         endif
 C A2 kernel(j-1)T A1T
         call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
@@ -9422,9 +9490,9 @@ cd      endif
 cd      write (iout,*)
 cd     &   'EELLO5: Contacts have occurred for peptide groups',i,j,
 cd     &   ' and',k,l
-      itk=itortyp(itype(k))
-      itl=itortyp(itype(l))
-      itj=itortyp(itype(j))
+      itk=itype2loc(itype(k))
+      itl=itype2loc(itype(l))
+      itj=itype2loc(itype(j))
       eello5_1=0.0d0
       eello5_2=0.0d0
       eello5_3=0.0d0
@@ -9493,7 +9561,7 @@ C Cartesian gradient
 c      goto 1112
 c1111  continue
 C Contribution from graph II 
-      call transpose2(EE(1,1,itk),auxmat(1,1))
+      call transpose2(EE(1,1,k),auxmat(1,1))
       call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
       vv(1)=pizda(1,1)+pizda(2,2)
       vv(2)=pizda(2,1)-pizda(1,2)
@@ -9574,7 +9642,7 @@ C Cartesian gradient
 cd        goto 1112
 C Contribution from graph IV
 cd1110    continue
-        call transpose2(EE(1,1,itl),auxmat(1,1))
+        call transpose2(EE(1,1,l),auxmat(1,1))
         call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
@@ -9647,7 +9715,7 @@ C Cartesian gradient
 cd        goto 1112
 C Contribution from graph IV
 1110    continue
-        call transpose2(EE(1,1,itj),auxmat(1,1))
+        call transpose2(EE(1,1,j),auxmat(1,1))
         call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
@@ -9944,7 +10012,7 @@ C       o     o       o     o                                                  C
 C       i             i                                                        C
 C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-      itk=itortyp(itype(k))
+      itk=itype2loc(itype(k))
       s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
       s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
       s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
@@ -10236,16 +10304,16 @@ C 4/7/01 AL Component s1 was removed, because it pertains to the respective
 C           energy moment and not to the cluster cumulant.
       iti=itortyp(itype(i))
       if (j.lt.nres-1) then
-        itj1=itortyp(itype(j+1))
+        itj1=itype2loc(itype(j+1))
       else
-        itj1=ntortyp
+        itj1=nloctyp
       endif
-      itk=itortyp(itype(k))
-      itk1=itortyp(itype(k+1))
+      itk=itype2loc(itype(k))
+      itk1=itype2loc(itype(k+1))
       if (l.lt.nres-1) then
-        itl1=itortyp(itype(l+1))
+        itl1=itype2loc(itype(l+1))
       else
-        itl1=ntortyp
+        itl1=nloctyp
       endif
 #ifdef MOMENT
       s1=dip(4,jj,i)*dip(4,kk,k)
@@ -10254,7 +10322,7 @@ C           energy moment and not to the cluster cumulant.
       s2=0.5d0*scalar2(b1(1,k),auxvec(1))
       call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1))
       s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
-      call transpose2(EE(1,1,itk),auxmat(1,1))
+      call transpose2(EE(1,1,k),auxmat(1,1))
       call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
       vv(1)=pizda(1,1)+pizda(2,2)
       vv(2)=pizda(2,1)-pizda(1,2)
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
 C           energy moment and not to the cluster cumulant.
 cd      write (2,*) 'eello_graph4: wturn6',wturn6
-      iti=itortyp(itype(i))
-      itj=itortyp(itype(j))
+      iti=itype2loc(itype(i))
+      itj=itype2loc(itype(j))
       if (j.lt.nres-1) then
-        itj1=itortyp(itype(j+1))
+        itj1=itype2loc(itype(j+1))
       else
-        itj1=ntortyp
+        itj1=nloctyp
       endif
-      itk=itortyp(itype(k))
+      itk=itype2loc(itype(k))
       if (k.lt.nres-1) then
-        itk1=itortyp(itype(k+1))
+        itk1=itype2loc(itype(k+1))
       else
-        itk1=ntortyp
+        itk1=nloctyp
       endif
-      itl=itortyp(itype(l))
+      itl=itype2loc(itype(l))
       if (l.lt.nres-1) then
-        itl1=itortyp(itype(l+1))
+        itl1=itype2loc(itype(l+1))
       else
-        itl1=ntortyp
+        itl1=nloctyp
       endif
 cd      write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l
 cd      write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk,
       j=i+4
       k=i+1
       l=i+3
-      iti=itortyp(itype(i))
-      itk=itortyp(itype(k))
-      itk1=itortyp(itype(k+1))
-      itl=itortyp(itype(l))
-      itj=itortyp(itype(j))
+      iti=itype2loc(itype(i))
+      itk=itype2loc(itype(k))
+      itk1=itype2loc(itype(k+1))
+      itl=itype2loc(itype(l))
+      itj=itype2loc(itype(j))
 cd      write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj
 cd      write (2,*) 'i',i,' k',k,' j',j,' l',l
 cd      if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then
@@ -11077,7 +11145,7 @@ C       do i=1,1
         if (itype(i).eq.ntyp1) cycle
 
         positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
-        if (positi.le.0) positi=positi+boxzsize
+        if (positi.le.0.0) positi=positi+boxzsize
 C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
@@ -11369,7 +11437,8 @@ C now costhet_grad
        enddo
 
       VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
-     &                    /VSolvSphere_div*4.0d0
+     &                    /VSolvSphere_div
+     &                    *wshield
 C now the gradient...
 C grad_shield is gradient of Calfa for peptide groups
 C      write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist,
@@ -11432,7 +11501,7 @@ C--------------------------------------------------------------------------
       include "DIMENSIONS"
       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 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
@@ -11448,4 +11517,496 @@ C        print *, x(i+1),yy(i),i
       gradtschebyshev=aux
       return
       end
+C------------------------------------------------------------------------
+C first for shielding is setting of function of side-chains
+       subroutine set_shield_fac2
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SHIELD'
+      include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+      double precision div77_81/0.974996043d0/,
+     &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+       double precision pep_side(3),long,side_calf(3),
+     &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+     &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+      do i=1,nres-1
+      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+      ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+      VolumeTotal=0.0
+      do k=1,nres
+       if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+       dist_pep_side=0.0
+       dist_side_calf=0.0
+       do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+      pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C      pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+      side_calf(j)=c(j,k+nres)-c(j,k)
+C      side_calf(j)=2.0d0
+      pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+      dist_pep_side=pep_side(j)**2+dist_pep_side
+      dist_side_calf=dist_side_calf+side_calf(j)**2
+      dist_pept_group=dist_pept_group+pept_group(j)**2
+      enddo
+       dist_pep_side=dsqrt(dist_pep_side)
+       dist_pept_group=dsqrt(dist_pept_group)
+       dist_side_calf=dsqrt(dist_side_calf)
+      do j=1,3
+        pep_side_norm(j)=pep_side(j)/dist_pep_side
+        side_calf_norm(j)=dist_side_calf
+      enddo
+C now sscale fraction
+       sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C       print *,buff_shield,"buff"
+C now sscale
+        if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient       
+        ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+        shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+        if (sh_frac_dist.gt.1.0) then
+         scale_fac_dist=1.0d0
+         do j=1,3
+         sh_frac_dist_grad(j)=0.0d0
+         enddo
+        else
+         scale_fac_dist=-sh_frac_dist*sh_frac_dist
+     &                   *(2.0d0*sh_frac_dist-3.0d0)
+         fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2)
+     &                  /dist_pep_side/buff_shield*0.5d0
+C remember for the final gradient multiply sh_frac_dist_grad(j) 
+C for side_chain by factor -2 ! 
+         do j=1,3
+         sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C         sh_frac_dist_grad(j)=0.0d0
+C         scale_fac_dist=1.0d0
+C         print *,"jestem",scale_fac_dist,fac_help_scale,
+C     &                    sh_frac_dist_grad(j)
+         enddo
+        endif
+C this is what is now we have the distance scaling now volume...
+      short=short_r_sidechain(itype(k))
+      long=long_r_sidechain(itype(k))
+      costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
+      sinthet=short/dist_pep_side*costhet
+C now costhet_grad
+C       costhet=0.6d0
+C       sinthet=0.8
+       costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4
+C       sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet
+C     &             -short/dist_pep_side**2/costhet)
+C       costhet_fac=0.0d0
+       do j=1,3
+         costhet_grad(j)=costhet_fac*pep_side(j)
+       enddo
+C remember for the final gradient multiply costhet_grad(j) 
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication  
+      pep_side0pept_group=0.0d0
+      do j=1,3
+      pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+      enddo
+      cosalfa=(pep_side0pept_group/
+     & (dist_pep_side*dist_side_calf))
+      fac_alfa_sin=1.0d0-cosalfa**2
+      fac_alfa_sin=dsqrt(fac_alfa_sin)
+      rkprim=fac_alfa_sin*(long-short)+short
+C      rkprim=short
+
+C now costhet_grad
+       cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
+C       cosphi=0.6
+       cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4
+       sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/
+     &      dist_pep_side**2)
+C       sinphi=0.8
+       do j=1,3
+         cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+     &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+     &*(long-short)/fac_alfa_sin*cosalfa/
+     &((dist_pep_side*dist_side_calf))*
+     &((side_calf(j))-cosalfa*
+     &((pep_side(j)/dist_pep_side)*dist_side_calf))
+C       cosphi_grad_long(j)=0.0d0
+        cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+     &*(long-short)/fac_alfa_sin*cosalfa
+     &/((dist_pep_side*dist_side_calf))*
+     &(pep_side(j)-
+     &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+C       cosphi_grad_loc(j)=0.0d0
+       enddo
+C      print *,sinphi,sinthet
+      VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet))
+     &                    /VSolvSphere_div
+C     &                    *wshield
+C now the gradient...
+      do j=1,3
+      grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+     &                +(sh_frac_dist_grad(j)*VofOverlap
+C  gradient po costhet
+     &       +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0*
+     &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+     &       sinphi/sinthet*costhet*costhet_grad(j)
+     &      +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+     & )*wshield
+C grad_shield_side is Cbeta sidechain gradient
+      grad_shield_side(j,ishield_list(i),i)=
+     &        (sh_frac_dist_grad(j)*-2.0d0
+     &        *VofOverlap
+     &       -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+     &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+     &       sinphi/sinthet*costhet*costhet_grad(j)
+     &      +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+     &       )*wshield        
+
+       grad_shield_loc(j,ishield_list(i),i)=
+     &       scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+     &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(
+     &       sinthet/sinphi*cosphi*cosphi_grad_loc(j)
+     &        ))
+     &        *wshield
+      enddo
+      VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+      enddo
+      fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
+C      write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+      enddo
+      return
+      end
+C-----------------------------------------------------------------------
+C-----------------------------------------------------------
+C This subroutine is to mimic the histone like structure but as well can be
+C utilizet to nanostructures (infinit) small modification has to be used to 
+C make it finite (z gradient at the ends has to be changes as well as the x,y
+C gradient has to be modified at the ends 
+C The energy function is Kihara potential 
+C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6)
+C 4eps is depth of well sigma is r_minimum r is distance from center of tube 
+C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
+C simple Kihara potential
+      subroutine calctube(Etube)
+       implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CALC'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
+      include 'COMMON.SBRIDGE'
+      double precision tub_r,vectube(3),enetube(maxres*2)
+      Etube=0.0d0
+      do i=1,2*nres
+        enetube(i)=0.0d0
+      enddo
+C first we calculate the distance from tube center
+C first sugare-phosphate group for NARES this would be peptide group 
+C for UNRES
+      do i=1,nres
+C lets ommit dummy atoms for now
+       if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+C now calculate distance from center of tube and direction vectors
+      vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+      vectube(1)=vectube(1)-tubecenter(1)
+      vectube(2)=vectube(2)-tubecenter(2)
+
+C      print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+C      print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+
+C as the tube is infinity we do not calculate the Z-vector use of Z
+C as chosen axis
+      vectube(3)=0.0d0
+C now calculte the distance
+       tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+C now normalize vector
+      vectube(1)=vectube(1)/tub_r
+      vectube(2)=vectube(2)/tub_r
+C calculte rdiffrence between r and r0
+      rdiff=tub_r-tubeR0
+C and its 6 power
+      rdiff6=rdiff**6.0d0
+C for vectorization reasons we will sumup at the end to avoid depenence of previous
+       enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6
+C       write(iout,*) "TU13",i,rdiff6,enetube(i)
+C       print *,rdiff,rdiff6,pep_aa_tube
+C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+C now we calculate gradient
+       fac=(-12.0d0*pep_aa_tube/rdiff6+
+     &       6.0d0*pep_bb_tube)/rdiff6/rdiff
+C       write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+C     &rdiff,fac
+
+C now direction of gg_tube vector
+        do j=1,3
+        gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+        gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+        enddo
+        enddo
+C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
+        do i=1,nres
+C Lets not jump over memory as we use many times iti
+         iti=itype(i)
+C lets ommit dummy atoms for now
+         if ((iti.eq.ntyp1)
+C in UNRES uncomment the line below as GLY has no side-chain...
+C      .or.(iti.eq.10)
+     &   ) cycle
+          vectube(1)=c(1,i+nres)
+          vectube(1)=mod(vectube(1),boxxsize)
+          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+          vectube(2)=c(2,i+nres)
+          vectube(2)=mod(vectube(2),boxxsize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+
+      vectube(1)=vectube(1)-tubecenter(1)
+      vectube(2)=vectube(2)-tubecenter(2)
+
+C as the tube is infinity we do not calculate the Z-vector use of Z
+C as chosen axis
+      vectube(3)=0.0d0
+C now calculte the distance
+       tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+C now normalize vector
+      vectube(1)=vectube(1)/tub_r
+      vectube(2)=vectube(2)/tub_r
+C calculte rdiffrence between r and r0
+      rdiff=tub_r-tubeR0
+C and its 6 power
+      rdiff6=rdiff**6.0d0
+C for vectorization reasons we will sumup at the end to avoid depenence of previous
+       sc_aa_tube=sc_aa_tube_par(iti)
+       sc_bb_tube=sc_bb_tube_par(iti)
+       enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6
+C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+C now we calculate gradient
+       fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+
+     &       6.0d0*sc_bb_tube/rdiff6/rdiff
+C now direction of gg_tube vector
+         do j=1,3
+          gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+          gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+         enddo
+        enddo
+        do i=1,2*nres
+          Etube=Etube+enetube(i)
+        enddo
+C        print *,"ETUBE", etube
+        return
+        end
+C TO DO 1) add to total energy
+C       2) add to gradient summation
+C       3) add reading parameters (AND of course oppening of PARAM file)
+C       4) add reading the center of tube
+C       5) add COMMONs
+C       6) add to zerograd
+
+C-----------------------------------------------------------------------
+C-----------------------------------------------------------
+C This subroutine is to mimic the histone like structure but as well can be
+C utilizet to nanostructures (infinit) small modification has to be used to 
+C make it finite (z gradient at the ends has to be changes as well as the x,y
+C gradient has to be modified at the ends 
+C The energy function is Kihara potential 
+C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6)
+C 4eps is depth of well sigma is r_minimum r is distance from center of tube 
+C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
+C simple Kihara potential
+      subroutine calctube2(Etube)
+       implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CALC'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
+      include 'COMMON.SBRIDGE'
+      double precision tub_r,vectube(3),enetube(maxres*2)
+      Etube=0.0d0
+      do i=1,2*nres
+        enetube(i)=0.0d0
+      enddo
+C first we calculate the distance from tube center
+C first sugare-phosphate group for NARES this would be peptide group 
+C for UNRES
+      do i=1,nres
+C lets ommit dummy atoms for now
+       if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+C now calculate distance from center of tube and direction vectors
+      vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+      vectube(1)=vectube(1)-tubecenter(1)
+      vectube(2)=vectube(2)-tubecenter(2)
+
+C      print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+C      print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+
+C as the tube is infinity we do not calculate the Z-vector use of Z
+C as chosen axis
+      vectube(3)=0.0d0
+C now calculte the distance
+       tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+C now normalize vector
+      vectube(1)=vectube(1)/tub_r
+      vectube(2)=vectube(2)/tub_r
+C calculte rdiffrence between r and r0
+      rdiff=tub_r-tubeR0
+C and its 6 power
+      rdiff6=rdiff**6.0d0
+C for vectorization reasons we will sumup at the end to avoid depenence of previous
+       enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6
+C       write(iout,*) "TU13",i,rdiff6,enetube(i)
+C       print *,rdiff,rdiff6,pep_aa_tube
+C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+C now we calculate gradient
+       fac=(-12.0d0*pep_aa_tube/rdiff6+
+     &       6.0d0*pep_bb_tube)/rdiff6/rdiff
+C       write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+C     &rdiff,fac
+
+C now direction of gg_tube vector
+        do j=1,3
+        gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+        gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+        enddo
+        enddo
+C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
+        do i=1,nres
+C Lets not jump over memory as we use many times iti
+         iti=itype(i)
+C lets ommit dummy atoms for now
+         if ((iti.eq.ntyp1)
+C in UNRES uncomment the line below as GLY has no side-chain...
+     &      .or.(iti.eq.10)
+     &   ) cycle
+          vectube(1)=c(1,i+nres)
+          vectube(1)=mod(vectube(1),boxxsize)
+          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+          vectube(2)=c(2,i+nres)
+          vectube(2)=mod(vectube(2),boxxsize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+
+      vectube(1)=vectube(1)-tubecenter(1)
+      vectube(2)=vectube(2)-tubecenter(2)
+C THIS FRAGMENT MAKES TUBE FINITE
+        positi=(mod(c(3,i+nres),boxzsize))
+        if (positi.le.0) positi=positi+boxzsize
+C       print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
+c for each residue check if it is in lipid or lipid water border area
+C       respos=mod(c(3,i+nres),boxzsize)
+       print *,positi,bordtubebot,buftubebot,bordtubetop
+       if ((positi.gt.bordtubebot)
+     & .and.(positi.lt.bordtubetop)) then
+C the energy transfer exist
+        if (positi.lt.buftubebot) then
+         fracinbuf=1.0d0-
+     &     ((positi-bordtubebot)/tubebufthick)
+C lipbufthick is thickenes of lipid buffore
+         sstube=sscalelip(fracinbuf)
+         ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
+         print *,ssgradtube, sstube,tubetranene(itype(i))
+         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+         gg_tube_SC(3,i)=gg_tube_SC(3,i)
+     &+ssgradtube*tubetranene(itype(i))
+         gg_tube(3,i-1)= gg_tube(3,i-1)
+     &+ssgradtube*tubetranene(itype(i))
+C         print *,"doing sccale for lower part"
+        elseif (positi.gt.buftubetop) then
+         fracinbuf=1.0d0-
+     &((bordtubetop-positi)/tubebufthick)
+         sstube=sscalelip(fracinbuf)
+         ssgradtube=sscagradlip(fracinbuf)/tubebufthick
+         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
+C     &+ssgradtube*tubetranene(itype(i))
+C         gg_tube(3,i-1)= gg_tube(3,i-1)
+C     &+ssgradtube*tubetranene(itype(i))
+C          print *, "doing sscalefor top part",sslip,fracinbuf
+        else
+         sstube=1.0d0
+         ssgradtube=0.0d0
+         enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+C         print *,"I am in true lipid"
+        endif
+        else
+C          sstube=0.0d0
+C          ssgradtube=0.0d0
+        cycle
+        endif ! if in lipid or buffor
+CEND OF FINITE FRAGMENT
+C as the tube is infinity we do not calculate the Z-vector use of Z
+C as chosen axis
+      vectube(3)=0.0d0
+C now calculte the distance
+       tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+C now normalize vector
+      vectube(1)=vectube(1)/tub_r
+      vectube(2)=vectube(2)/tub_r
+C calculte rdiffrence between r and r0
+      rdiff=tub_r-tubeR0
+C and its 6 power
+      rdiff6=rdiff**6.0d0
+C for vectorization reasons we will sumup at the end to avoid depenence of previous
+       sc_aa_tube=sc_aa_tube_par(iti)
+       sc_bb_tube=sc_bb_tube_par(iti)
+       enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6)
+     &                 *sstube+enetube(i+nres)
+C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+C now we calculate gradient
+       fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+
+     &       6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube
+C now direction of gg_tube vector
+         do j=1,3
+          gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+          gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+         enddo
+         gg_tube_SC(3,i)=gg_tube_SC(3,i)
+     &+ssgradtube*enetube(i+nres)/sstube
+         gg_tube(3,i-1)= gg_tube(3,i-1)
+     &+ssgradtube*enetube(i+nres)/sstube
+
+        enddo
+        do i=1,2*nres
+          Etube=Etube+enetube(i)
+        enddo
+C        print *,"ETUBE", etube
+        return
+        end
+C TO DO 1) add to total energy
+C       2) add to gradient summation
+C       3) add reading parameters (AND of course oppening of PARAM file)
+C       4) add reading the center of tube
+C       5) add COMMONs
+C       6) add to zerograd