Correction of D-AA to parmread

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 91a7c94..a20ec7b 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -6329,6 +6329,9 @@ c
           else
             do j=1,nbi
               diff=vbld(i+nres)-vbldsc0(j,iti) 
+            if (energy_dec)  write (iout,*)
+     &      "estr sc",i,iti,vbld(i+nres),vbldsc0(j,iti),diff,
+     &      AKSC(j,iti),AKSC(j,iti)*diff*diff
               ud(j)=aksc(j,iti)*diff
               u(j)=abond0(j,iti)+0.5d0*ud(j)*diff
             enddo
@@ -6662,6 +6665,7 @@ C
       etheta=0.0D0
       do i=ithet_start,ithet_end
 c        print *,i,itype(i-1),itype(i),itype(i-2)
+C        if (itype(i-1).eq.ntyp1) cycle
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
 C        print *,i,theta(i)
@@ -12337,10 +12341,31 @@ 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,boxysize)
-          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+C      vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+C          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+C      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+C          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+      xmin=boxxsize
+      ymin=boxysize
+        do j=-1,1
+         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
       vectube(1)=vectube(1)-tubecenter(1)
       vectube(2)=vectube(2)-tubecenter(2)
 
@@ -12359,14 +12384,61 @@ C calculte rdiffrence between r and r0
       rdiff=tub_r-tubeR0
 C and its 6 power
       rdiff6=rdiff**6.0d0
+C THIS FRAGMENT MAKES TUBE FINITE
+        positi=mod((c(3,i)+c(3,i+1))/2.0d0,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)=enetube(i)+sstube*tubetranenepep
+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 sccale for lower part"
+        elseif (positi.gt.buftubetop) then
+         fracinbuf=1.0d0-
+     &((bordtubetop-positi)/tubebufthick)
+         sstube=sscalelip(fracinbuf)
+         ssgradtube=sscagradlip(fracinbuf)/tubebufthick
+         enetube(i)=enetube(i)+sstube*tubetranenepep
+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)=enetube(i)+sstube*tubetranenepep
+C         print *,"I am in true lipid"
+        endif
+        else
+C          sstube=0.0d0
+C          ssgradtube=0.0d0
+        cycle
+        endif ! if in lipid or buffor
+
 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
+       enetube(i)=enetube(i)+sstube*
+     &(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
+     &       6.0d0*pep_bb_tube)/rdiff6/rdiff*sstube
 C       write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
 C     &rdiff,fac
 
@@ -12375,6 +12447,11 @@ C now direction of gg_tube vector
         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
+         gg_tube(3,i)=gg_tube(3,i)
+     &+ssgradtube*enetube(i)/sstube/2.0d0
+         gg_tube(3,i-1)= gg_tube(3,i-1)
+     &+ssgradtube*enetube(i)/sstube/2.0d0
+
         enddo
 C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
 C        print *,gg_tube(1,0),"TU"

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 78d0afb..dc040cb 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -1038,22 +1038,26 @@ C
         write (iout,*) 'Type',i
         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
         endif
-c        B1(1,i)  = b(3)
-c        B1(2,i)  = b(5)
-c        B1(1,-i) = b(3)
-c        B1(2,-i) = -b(5)
+        b(3,-i)= b(3,i)
+        b(5,-i)=-b(5,i)
+        b(4,-i)=-b(4,i)
+        b(2,-i)= b(2,i)     
+        B1(1,i)  = b(3,i)
+        B1(2,i)  = b(5,i)
+        B1(1,-i) = b(3,i)
+        B1(2,-i) = -b(5,i)
 c        b1(1,i)=0.0d0
 c        b1(2,i)=0.0d0
-c        B1tilde(1,i) = b(3)
-c        B1tilde(2,i) =-b(5)
-c        B1tilde(1,-i) =-b(3)
-c        B1tilde(2,-i) =b(5)
+        B1tilde(1,i) = b(3,i)
+        B1tilde(2,i) =-b(5,i)
+        B1tilde(1,-i) =-b(3,i)
+        B1tilde(2,-i) =b(5,i)
 c        b1tilde(1,i)=0.0d0
 c        b1tilde(2,i)=0.0d0
-c        B2(1,i)  = b(2)
-c        B2(2,i)  = b(4)
-c        B2(1,-i)  =b(2)
-c        B2(2,-i)  =-b(4)
+        B2(1,i)  = b(2,i)
+        B2(2,i)  = b(4,i)
+        B2(1,-i)  =b(2,i)
+        B2(2,-i)  =-b(4,i)
 cc        B1tilde(1,i) = b(3,i)
 cc        B1tilde(2,i) =-b(5,i)
 C        B1tilde(1,-i) =-b(3,i)

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index b5f5705..6c943cc 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -9188,3 +9188,692 @@ C      write(2,*) "TU",rpp(1,1),short,long,buff_shield
       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=itube_start,itube_end
+        enetube(i)=0.0d0
+        enetube(i+nres)=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=itube_start,itube_end
+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
+      xmin=boxxsize
+      ymin=boxysize
+        do j=-1,1
+         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+       
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
+      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)
+C        print *,gg_tube(1,0),"TU"
+
+
+       do i=itube_start,itube_end
+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
+      xmin=boxxsize
+      ymin=boxysize
+        do j=-1,1
+         vectube(1)=mod((c(1,i+nres)),boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i+nres)),boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
+C          write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2),
+C     &     tubecenter(2)
+      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=itube_start,itube_end
+          Etube=Etube+enetube(i)+enetube(i+nres)
+        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=itube_start,itube_end
+        enetube(i)=0.0d0
+        enetube(i+nres)=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=itube_start,itube_end
+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
+C      vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+C          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+C      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+C          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+      xmin=boxxsize
+      ymin=boxysize
+        do j=-1,1
+         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
+      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 THIS FRAGMENT MAKES TUBE FINITE
+        positi=mod((c(3,i)+c(3,i+1))/2.0d0,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)=enetube(i)+sstube*tubetranenepep
+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 sccale for lower part"
+        elseif (positi.gt.buftubetop) then
+         fracinbuf=1.0d0-
+     &((bordtubetop-positi)/tubebufthick)
+         sstube=sscalelip(fracinbuf)
+         ssgradtube=sscagradlip(fracinbuf)/tubebufthick
+         enetube(i)=enetube(i)+sstube*tubetranenepep
+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)=enetube(i)+sstube*tubetranenepep
+C         print *,"I am in true lipid"
+        endif
+        else
+C          sstube=0.0d0
+C          ssgradtube=0.0d0
+        cycle
+        endif ! if in lipid or buffor
+
+C for vectorization reasons we will sumup at the end to avoid depenence of previous
+       enetube(i)=enetube(i)+sstube*
+     &(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*sstube
+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
+         gg_tube(3,i)=gg_tube(3,i)
+     &+ssgradtube*enetube(i)/sstube/2.0d0
+         gg_tube(3,i-1)= gg_tube(3,i-1)
+     &+ssgradtube*enetube(i)/sstube/2.0d0
+
+        enddo
+C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
+C        print *,gg_tube(1,0),"TU"
+        do i=itube_start,itube_end
+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),boxysize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+
+      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))
+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 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=itube_start,itube_end
+          Etube=Etube+enetube(i)+enetube(i+nres)
+        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 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 calcnano(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),
+     & enecavtube(maxres*2)
+      Etube=0.0d0
+      do i=itube_start,itube_end
+        enetube(i)=0.0d0
+        enetube(i+nres)=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=itube_start,itube_end
+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
+      xmin=boxxsize
+      ymin=boxysize
+      zmin=boxzsize
+
+        do j=-1,1
+         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+         vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+         vectube(3)=vectube(3)+boxzsize*j
+
+
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+         zminact=abs(vectube(3)-tubecenter(3))
+
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+           if (zmin.gt.zminact) then
+             zmin=zminact
+             ztemp=vectube(3)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
+      vectube(3)=ztemp
+
+      vectube(1)=vectube(1)-tubecenter(1)
+      vectube(2)=vectube(2)-tubecenter(2)
+      vectube(3)=vectube(3)-tubecenter(3)
+
+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
+C      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
+      vectube(3)=vectube(3)/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
+         if (acavtubpep.eq.0.0d0) then
+C go to 667
+         enecavtube(i)=0.0
+         faccav=0.0
+         else
+         denominator=(1.0+dcavtubpep*rdiff6*rdiff6)
+         enecavtube(i)=
+     &   (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep)
+     &   /denominator
+         enecavtube(i)=0.0
+         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff))
+     &   *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)
+     &   +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0)
+     &   /denominator**2.0d0
+C         faccav=0.0
+C         fac=fac+faccav
+C 667     continue
+         endif
+C         print *,"TUT",i,iti,rdiff,rdiff6,acavtubpep,denominator,
+C     &   enecavtube(i),faccav
+C         print *,"licz=",
+C     & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+CX         print *,"finene=",enetube(i+nres)+enecavtube(i)
+         
+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
+
+       do i=itube_start,itube_end
+        enecavtube(i)=0.0 
+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
+      xmin=boxxsize
+      ymin=boxysize
+      zmin=boxzsize
+        do j=-1,1
+         vectube(1)=mod((c(1,i+nres)),boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i+nres)),boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+         vectube(3)=mod((c(3,i+nres)),boxzsize)
+         vectube(3)=vectube(3)+boxzsize*j
+
+
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+         zminact=abs(vectube(3)-tubecenter(3))
+
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+           if (zmin.gt.zminact) then
+             zmin=zminact
+             ztemp=vectube(3)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
+      vectube(3)=ztemp
+
+C          write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2),
+C     &     tubecenter(2)
+      vectube(1)=vectube(1)-tubecenter(1)
+      vectube(2)=vectube(2)-tubecenter(2)
+      vectube(3)=vectube(3)-tubecenter(3)
+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
+      vectube(3)=vectube(3)/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       enetube(i+nres)=0.0d0
+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       fac=0.0
+C now direction of gg_tube vector
+C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
+         if (acavtub(iti).eq.0.0d0) then
+C go to 667
+         enecavtube(i+nres)=0.0
+         faccav=0.0
+         else
+         denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
+         enecavtube(i+nres)=
+     &   (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+     &   /denominator
+C         enecavtube(i)=0.0
+         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff))
+     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)
+     &   +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0)
+     &   /denominator**2.0d0
+C         faccav=0.0
+         fac=fac+faccav
+C 667     continue
+         endif
+C         print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator,
+C     &   enecavtube(i),faccav
+C         print *,"licz=",
+C     & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+C         print *,"finene=",enetube(i+nres)+enecavtube(i)
+         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
+C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
+C        do i=itube_start,itube_end
+C        enecav(i)=0.0        
+C        iti=itype(i)
+C        if (acavtub(iti).eq.0.0) cycle
+        
+
+
+        do i=itube_start,itube_end
+          Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i)
+     & +enecavtube(i+nres)
+        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
+