AFM+tube+kcc

diff --git a/PARAM/tube.parm b/PARAM/tube.parm
index d1783aa..182c68f 100644 (file)
--- a/PARAM/tube.parm
+++ b/PARAM/tube.parm
@@ -1,26 +1,26 @@
-   200.00 2.20
-   1.00 0.50
-   1.00 0.20
-   1.00 0.30
-   1.00 1.00
-   1.00 1.00
-   200.00 2.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
-   1.00 1.00
+   200.00 2.20 3.0
+   1.00 0.50 3.0
+   1.00 0.20 3.0
+   1.00 0.30 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   200.00 2.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
+   1.00 1.00 3.0
 

diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN
index 8443ef4..372fa85 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CHAIN
+++ b/source/unres/src_MD-M/COMMON.CHAIN
@@ -19,11 +19,13 @@
       double precision boxxsize,boxysize,boxzsize,enecut,sscut,
      & sss,sssgrad,
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick,
-     & tubecenter,tubeR0     
+     & tubecenter,tubeR0,
+     & buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick     
       common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
       common /afm/  forceAFMconst, distafminit,afmend,afmbeg,
      & velAFMconst,
      & totTafm
-      common /tube/ tubecenter(3),tubeR0
+      common /tube/ tubecenter(3),tubeR0,
+     & buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick          
 

diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT
index c43ef6a..14d92ef 100644 (file)
--- a/source/unres/src_MD-M/COMMON.INTERACT
+++ b/source/unres/src_MD-M/COMMON.INTERACT
@@ -41,6 +41,8 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters.
      & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp),
      & distchainmax,nbondterm(ntyp)
 C 01/29/15 Lipidic parameters
-      double precision   pepliptran,liptranene
+      double precision   pepliptran,liptranene,
+     &tubetranene, tubetranenepep
       common /lipid/ pepliptran,liptranene(ntyp)
+      common /tubepar/ tubetranene(ntyp), tubetranenepep
 

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 e86bb6e..21c0197 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -307,9 +307,11 @@ C      print *,"za lipidami"
       else if (selfguide.gt.0) then
         call AFMvel(Eafmforce)
       endif
-      if (TUBElog.gt.0) then
+      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
@@ -11725,8 +11727,13 @@ C for UNRES
 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)=(c(1,i)+c(1,i+1))/2.0d0-tubecenter(1)
-      vectube(2)=(c(2,i)+c(2,i+1))/2.0d0-tubecenter(2)
+      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)
 
@@ -11768,8 +11775,15 @@ C lets ommit dummy atoms for now
 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)-tubecenter(1)
-      vectube(2)=c(2,i+nres)-tubecenter(2)
+          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
@@ -11810,3 +11824,189 @@ 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
+

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index a38147f..5001b6b 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -1367,7 +1367,7 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
       epspeptube=dabs(epspeptube)
       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
-      write(iout,*) pep_aa_tube,pep_bb_tube
+      write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
       do i=1,ntyp
        read(itube,*) epssctube,sigmasctube
        sigmasctube=sigmasctube**6
@@ -1375,7 +1375,7 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
        epssctube=dabs(epssctube)
        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
-      write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i)
+      write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
       enddo
 
 #ifdef OLDSCP

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 89ddeac..e40717f 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -147,11 +147,6 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
       print *,'AFMlog',AFMlog,selfguide,"KUPA"
       call readi(controlcard,'TUBEMOD',tubelog,0)
       write (iout,*) TUBElog,"TUBEMODE"
-      if (TUBElog.gt.0) then
-       call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
-       call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
-       call reada(controlcard,"RTUBE",tubeR0,0.0d0)
-      endif
       call readi(controlcard,'IPRINT',iprint,0)
 C SHIELD keyword sets if the shielding effect of side-chains is used
 C 0 denots no shielding is used all peptide are equally despite the 
@@ -270,6 +265,16 @@ C      endif
       write(iout,*) "bufliptop=",bufliptop
       write(iout,*) "buflipbot=",buflipbot
       write (iout,*) "SHIELD MODE",shield_mode
+      if (TUBElog.gt.0) then
+       call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+       call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
+       call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+       call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
+       call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
+       call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
+       buftubebot=bordtubebot+tubebufthick
+       buftubetop=bordtubetop-tubebufthick
+      endif
       if (shield_mode.gt.0) then
       pi=3.141592d0
 C VSolvSphere the volume of solving sphere