+! Lipid transfer energy function
+ subroutine Eliptransfer(eliptran)
+!C this is done by Adasko
+!C print *,"wchodze"
+!C structure of box:
+!C water
+!C--bordliptop-- buffore starts
+!C--bufliptop--- here true lipid starts
+!C lipid
+!C--buflipbot--- lipid ends buffore starts
+!C--bordlipbot--buffore ends
+ real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip
+ integer :: i
+ eliptran=0.0
+! print *, "I am in eliptran"
+ do i=ilip_start,ilip_end
+!C do i=1,1
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1).or.(i.eq.nres))&
+ cycle
+
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),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
+ if ((positi.gt.bordlipbot) &
+ .and.(positi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (positi.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((positi-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+!C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+
+!C print *,"doing sccale for lower part"
+!C print *,i,sslip,fracinbuf,ssgradlip
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+!C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+!C print *, "doing sscalefor top part"
+!C print *,i,sslip,fracinbuf,ssgradlip
+ else
+ eliptran=eliptran+pepliptran
+!C print *,"I am in true lipid"
+ endif
+!C else
+!C eliptran=elpitran+0.0 ! I am in water
+ endif
+ if (energy_dec) write(iout,*) i,"eliptran=",eliptran,positi,sslip
+ enddo
+! here starts the side chain transfer
+ do i=ilip_start,ilip_end
+ if (itype(i,1).eq.ntyp1) cycle
+ 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)
+!C print *,positi,bordlipbot,buflipbot
+ if ((positi.gt.bordlipbot) &
+ .and.(positi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (positi.lt.buflipbot) then
+ fracinbuf=1.0d0- &
+ ((positi-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
+ gliptranx(3,i)=gliptranx(3,i) &
+ +ssgradlip*liptranene(itype(i,1))
+ gliptranc(3,i-1)= gliptranc(3,i-1) &
+ +ssgradlip*liptranene(itype(i,1))
+!C print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0- &
+ ((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
+ gliptranx(3,i)=gliptranx(3,i) &
+ +ssgradlip*liptranene(itype(i,1))
+ gliptranc(3,i-1)= gliptranc(3,i-1) &
+ +ssgradlip*liptranene(itype(i,1))
+!C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i,1))
+!C print *,"I am in true lipid"
+ endif
+ endif ! if in lipid or buffor
+!C else
+!C eliptran=elpitran+0.0 ! I am in water
+ if (energy_dec) write(iout,*) i,"eliptran=",eliptran
+ enddo
+ return
+ end subroutine Eliptransfer
+!----------------------------------NANO FUNCTIONS
+!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)
+ real(kind=8),dimension(3) :: vectube
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
+ sc_aa_tube,sc_bb_tube
+ integer :: i,j,iti
+ 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 for UNRES
+ do i=itube_start,itube_end
+!C lets ommit dummy atoms for now
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
+!C now calculate distance from center of tube and direction vectors
+ xmin=boxxsize
+ ymin=boxysize
+! Find minimum distance in periodic box
+ 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,1)
+!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
+ 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 subroutine calctube
+!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 7) allocate matrices
+
+
+!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)
+ real(kind=8),dimension(3) :: vectube
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,&
+ sstube,ssgradtube,sc_aa_tube,sc_bb_tube
+ integer:: i,j,iti
+ 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,1).eq.ntyp1).or.(itype(i+1,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)
+!C 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
+!C print *,ssgradtube, sstube,tubetranene(itype(i,1))
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!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,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!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,1)
+!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)
+!C 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
+!C print *,ssgradtube, sstube,tubetranene(itype(i,1))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!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,1))
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!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,1))
+!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 subroutine calctube2
+!=====================================================================================================================================
+ subroutine calcnano(Etube)
+ real(kind=8),dimension(3) :: vectube
+
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,&
+ sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact
+ integer:: i,j,iti,r
+
+ Etube=0.0d0
+! print *,itube_start,itube_end,"poczatek"
+ 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,1).eq.ntyp1).or.(itype(i+1,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)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=dabs(vectube(1)-tubecenter(1))
+ yminact=dabs(vectube(2)-tubecenter(2))
+ zminact=dabs(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.0d0+dcavtubpep*rdiff6*rdiff6)
+ enecavtube(i)= &
+ (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep) &
+ /denominator
+ enecavtube(i)=0.0
+ faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff)) &
+ *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff) &
+ +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0) &
+ /denominator**2.0d0
+!C faccav=0.0
+!C fac=fac+faccav
+!C 667 continue
+ endif
+ if (energy_dec) write(iout,*),i,rdiff,enetube(i),enecavtube(i)
+ 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.0d0
+!C Lets not jump over memory as we use many times iti
+ iti=itype(i,1)
+!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)=dmod((c(1,i+nres)),boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=dmod((c(2,i+nres)),boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ vectube(3)=dmod((c(3,i+nres)),boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=dabs(vectube(1)-tubecenter(1))
+ yminact=dabs(vectube(2)-tubecenter(2))
+ zminact=dabs(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
+ 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.0d0
+ faccav=0.0d0
+ else
+ denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6)
+ enecavtube(i+nres)= &
+ (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti)) &
+ /denominator
+!C enecavtube(i)=0.0
+ faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff)) &
+ *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(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
+ if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres),enecavtube(i+nres)
+ enddo
+
+
+
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) &
+ +enecavtube(i+nres)
+ enddo
+! do i=1,20
+! print *,"begin", i,"a"
+! do r=1,10000
+! rdiff=r/100.0d0
+! rdiff6=rdiff**6.0d0
+! sc_aa_tube=sc_aa_tube_par(i)
+! sc_bb_tube=sc_bb_tube_par(i)
+! enetube(i)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+! denominator=(1.0d0+dcavtub(i)*rdiff6*rdiff6)
+! enecavtube(i)= &
+! (bcavtub(i)*rdiff+acavtub(i)*dsqrt(rdiff)+ccavtub(i)) &
+! /denominator
+
+! print '(5(f10.3,1x))',rdiff,enetube(i),enecavtube(i),enecavtube(i)+enetube(i)
+! enddo
+! print *,"end",i,"a"
+! enddo
+!C print *,"ETUBE", etube
+ return
+ end subroutine calcnano
+
+!===============================================
+!--------------------------------------------------------------------------------
+!C first for shielding is setting of function of side-chains
+
+ subroutine set_shield_fac2
+ real(kind=8) :: div77_81=0.974996043d0, &
+ div4_81=0.2222222222d0
+ real (kind=8) :: dist_pep_side,dist_side_calf,dist_pept_group, &
+ scale_fac_dist,fac_help_scale,VofOverlap,VolumeTotal,costhet,&
+ short,long,sinthet,costhet_fac,sh_frac_dist,rkprim,cosphi, &
+ sinphi,cosphi_fac,pep_side0pept_group,cosalfa,fac_alfa_sin
+!C the vector between center of side_chain and peptide group
+ real(kind=8),dimension(3) :: pep_side_long,side_calf, &
+ pept_group,costhet_grad,cosphi_grad_long, &
+ cosphi_grad_loc,pep_side_norm,side_calf_norm, &
+ sh_frac_dist_grad,pep_side
+ integer i,j,k
+!C write(2,*) "ivec",ivec_start,ivec_end
+ do i=1,nres
+ fac_shield(i)=0.0d0
+ do j=1,3
+ grad_shield(j,i)=0.0d0
+ enddo
+ enddo
+ do i=ivec_start,ivec_end
+!C do i=1,nres-1
+!C if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
+ ishield_list(i)=0
+ if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
+!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,1).eq.ntyp1).or.(itype(k,1).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=sqrt(dist_pep_side)
+ dist_pept_group=sqrt(dist_pept_group)
+ dist_side_calf=sqrt(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 print *,ishield_list(i),i
+!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
+ 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,1))
+ long=long_r_sidechain(itype(k,1))
+ 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,itype(i,1),fac_shield(i)
+ enddo
+ return
+ end subroutine set_shield_fac2
+!----------------------------------------------------------------------------
+! SOUBROUTINE FOR AFM
+ subroutine AFMvel(Eafmforce)
+ use MD_data, only:totTafm
+ real(kind=8),dimension(3) :: diffafm
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i
+!C Only for check grad COMMENT if not used for checkgrad
+!C totT=3.0d0
+!C--------------------------------------------------------
+!C print *,"wchodze"
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! totTafm=3.0
+ Eafmforce=0.5d0*forceAFMconst &
+ *(distafminit+totTafm*velAFMconst-afmdist)**2
+!C Eafmforce=-forceAFMconst*(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ gradafm(i,afmbeg-1)=forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ enddo
+! print *,'AFM',Eafmforce,totTafm*velAFMconst,afmdist
+ return
+ end subroutine AFMvel
+!---------------------------------------------------------
+ subroutine AFMforce(Eafmforce)
+
+ real(kind=8),dimension(3) :: diffafm
+! real(kind=8) ::afmdist
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! print *,afmdist,distafminit
+ Eafmforce=-forceAFMconst*(afmdist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/afmdist
+ gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/afmdist
+ enddo
+!C print *,'AFM',Eafmforce
+ return
+ end subroutine AFMforce
+