+CCC----------------------------------------------
+ subroutine Eliptransfer(eliptran)
+ 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'
+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
+ eliptran=0.0
+ do i=ilip_start,ilip_end
+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.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
+ enddo
+C print *, "nic nie bylo w lipidzie?"
+C now multiply all by the peptide group transfer factor
+C eliptran=eliptran*pepliptran
+C now the same for side chains
+CV do i=1,1
+ do i=ilip_start,ilip_end
+ if (itype(i).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))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+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))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i))
+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
+ enddo
+ return
+ end
+C---------------------------------------------------------
+C AFM soubroutine for constant force
+ subroutine AFMforce(Eafmforce)
+ 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'
+ real*8 diffafm(3)
+ dist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ dist=dist+diffafm(i)**2
+ enddo
+ dist=dsqrt(dist)
+ Eafmforce=-forceAFMconst*(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist
+ gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist
+ enddo
+C print *,'AFM',Eafmforce
+ return
+ end
+C---------------------------------------------------------
+C AFM subroutine with pseudoconstant velocity
+ subroutine AFMvel(Eafmforce)
+ 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'
+ real*8 diffafm(3)
+C Only for check grad COMMENT if not used for checkgrad
+C totT=3.0d0
+C--------------------------------------------------------
+C print *,"wchodze"
+ dist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ dist=dist+diffafm(i)**2
+ enddo
+ dist=dsqrt(dist)
+ Eafmforce=0.5d0*forceAFMconst
+ & *(distafminit+totTafm*velAFMconst-dist)**2
+C Eafmforce=-forceAFMconst*(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*
+ &(distafminit+totTafm*velAFMconst-dist)
+ &*diffafm(i)/dist
+ gradafm(i,afmbeg-1)=forceAFMconst*
+ &(distafminit+totTafm*velAFMconst-dist)
+ &*diffafm(i)/dist
+ enddo
+C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
+ return
+ end
+C-----------------------------------------------------------
+C first for shielding is setting of function of side-chains
+ subroutine set_shield_fac
+ 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.0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
+ & /dist_pep_side/buff_shield*0.5
+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 print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
+ enddo
+ endif
+C if ((i.eq.3).and.(k.eq.2)) then
+C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
+C & ,"TU"
+C 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.0+short**2/dist_pep_side**2)
+C now costhet_grad
+C costhet=0.0d0
+ costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
+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.0
+ 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.0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
+
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5/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))
+
+ cosphi_grad_loc(j)=cosphi**3*0.5/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)
+ enddo
+
+ VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
+ & /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,
+C & costhet,cosphi
+C write(iout,*) "cosphi_compon",i,k,pep_side0pept_group,
+C & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k)
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)
+C gradient po costhet
+ &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
+ &-scale_fac_dist*(cosphi_grad_long(j))
+ &/(1.0-cosphi) )*div77_81
+ &*VofOverlap
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*-2.0d0
+ & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
+ & +scale_fac_dist*(cosphi_grad_long(j))
+ & *2.0d0/(1.0-cosphi))
+ & *div77_81*VofOverlap
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*cosphi_grad_loc(j)
+ & *2.0d0/(1.0-cosphi)
+ & *div77_81*VofOverlap
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*div77_81+div4_81
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+ enddo
+ return
+ end
+C--------------------------------------------------------------------------
+ double precision function tschebyshev(m,n,x,y)
+ implicit none
+ include "DIMENSIONS"
+ integer i,m,n
+ double precision x(n),y,yy(0:maxvar),aux
+c Tschebyshev polynomial. Note that the first term is omitted
+c m=0: the constant term is included
+c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=y
+ do i=2,n
+ yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i)*yy(i)
+ enddo
+ tschebyshev=aux
+ return
+ end
+C--------------------------------------------------------------------------
+ double precision function gradtschebyshev(m,n,x,y)
+ implicit none
+ 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 m=0: the constant term is included
+c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=2.0d0*y
+ do i=2,n
+ yy(i)=2*y*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i+1)*yy(i)*(i+1)
+C print *, x(i+1),yy(i),i
+ enddo
+ gradtschebyshev=aux
+ return
+ end
+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'
+ include 'COMMON.LOCAL'
+
+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 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
+C the line belowe needs to be changed for FGPROC>1
+ do i=ivec_start,ivec_end
+C do i=1,nres-1
+C if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+ if ((itype(i).eq.ntyp1).and.itype(i+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).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 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
+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,itype(i),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,boxysize)
+ if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+ 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),boxysize)
+ if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+
+ 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,boxysize)
+ if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+ 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),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=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