+
+ return
+ end subroutine dyn_ssbond_ene
+!--------------------------------------------------------------------------
+ subroutine triple_ssbond_ene(resi,resj,resk,eij)
+! implicit none
+! Includes
+ use calc_data
+ use comm_sschecks
+! include 'DIMENSIONS'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ use MD_data
+! include 'COMMON.MD'
+! use MD, only: totT,t_bath
+#endif
+#endif
+ double precision h_base
+ external h_base
+
+!c Input arguments
+ integer resi,resj,resk,m,itypi,itypj,itypk
+
+!c Output arguments
+ double precision eij,eij1,eij2,eij3
+
+!c Local variables
+ logical havebond
+!c integer itypi,itypj,k,l
+ double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
+ double precision xik,yik,zik,xjk,yjk,zjk,dxk,dyk,dzk
+ double precision sig0ij,ljd,sig,fac,e1,e2
+ double precision dcosom1(3),dcosom2(3),ed
+ double precision pom1,pom2
+ double precision ljA,ljB,ljXs
+ double precision d_ljB(1:3)
+ double precision ssA,ssB,ssC,ssXs
+ double precision ssxm,ljxm,ssm,ljm
+ double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+ eij=0.0
+ if (dtriss.eq.0) return
+ i=resi
+ j=resj
+ k=resk
+!C write(iout,*) resi,resj,resk
+ itypi=itype(i,1)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ itypj=itype(j,1)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+ itypk=itype(k,1)
+ xk=c(1,nres+k)
+ yk=c(2,nres+k)
+ zk=c(3,nres+k)
+
+ dxk=dc_norm(1,nres+k)
+ dyk=dc_norm(2,nres+k)
+ dzk=dc_norm(3,nres+k)
+ dscj_inv=vbld_inv(k+nres)
+ xij=xj-xi
+ xik=xk-xi
+ xjk=xk-xj
+ yij=yj-yi
+ yik=yk-yi
+ yjk=yk-yj
+ zij=zj-zi
+ zik=zk-zi
+ zjk=zk-zj
+ rrij=(xij*xij+yij*yij+zij*zij)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+ rrik=(xik*xik+yik*yik+zik*zik)
+ rik=dsqrt(rrik)
+ rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
+ rjk=dsqrt(rrjk)
+!C there are three combination of distances for each trisulfide bonds
+!C The first case the ith atom is the center
+!C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
+!C distance y is second distance the a,b,c,d are parameters derived for
+!C this problem d parameter was set as a penalty currenlty set to 1.
+ if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
+ eij1=0.0d0
+ else
+ eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
+ endif
+!C second case jth atom is center
+ if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
+ eij2=0.0d0
+ else
+ eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
+ endif
+!C the third case kth atom is the center
+ if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
+ eij3=0.0d0
+ else
+ eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
+ endif
+!C eij2=0.0
+!C eij3=0.0
+!C eij1=0.0
+ eij=eij1+eij2+eij3
+!C write(iout,*)i,j,k,eij
+!C The energy penalty calculated now time for the gradient part
+!C derivative over rij
+ fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
+ gg(1)=xij*fac/rij
+ gg(2)=yij*fac/rij
+ gg(3)=zij*fac/rij
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,j)=gvdwx(m,j)+gg(m)
+ enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+!C now derivative over rik
+ fac=-eij1**2/dtriss* &
+ (-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xik*fac/rik
+ gg(2)=yik*fac/rik
+ gg(3)=zik*fac/rik
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+!C now derivative over rjk
+ fac=-eij2**2/dtriss* &
+ (-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)- &
+ eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xjk*fac/rjk
+ gg(2)=yjk*fac/rjk
+ gg(3)=zjk*fac/rjk
+ do m=1,3
+ gvdwx(m,j)=gvdwx(m,j)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,j)=gvdwc(l,j)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+ return
+ end subroutine triple_ssbond_ene
+
+
+
+!-----------------------------------------------------------------------------
+ real(kind=8) function h_base(x,deriv)
+! A smooth function going 0->1 in range [0,1]
+! It should NOT be called outside range [0,1], it will not work there.
+ implicit none
+
+! Input arguments
+ real(kind=8) :: x
+
+! Output arguments
+ real(kind=8) :: deriv
+
+! Local variables
+ real(kind=8) :: xsq
+
+
+! Two parabolas put together. First derivative zero at extrema
+!$$$ if (x.lt.0.5D0) then
+!$$$ h_base=2.0D0*x*x
+!$$$ deriv=4.0D0*x
+!$$$ else
+!$$$ deriv=1.0D0-x
+!$$$ h_base=1.0D0-2.0D0*deriv*deriv
+!$$$ deriv=4.0D0*deriv
+!$$$ endif
+
+! Third degree polynomial. First derivative zero at extrema
+ h_base=x*x*(3.0d0-2.0d0*x)
+ deriv=6.0d0*x*(1.0d0-x)
+
+! Fifth degree polynomial. First and second derivatives zero at extrema
+!$$$ xsq=x*x
+!$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
+!$$$ deriv=x-1.0d0
+!$$$ deriv=deriv*deriv
+!$$$ deriv=30.0d0*xsq*deriv
+
+ return
+ end function h_base
+!-----------------------------------------------------------------------------
+ subroutine dyn_set_nss
+! Adjust nss and other relevant variables based on dyn_ssbond_ij
+! implicit none
+ use MD_data, only: totT,t_bath
+! Includes
+! include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.SETUP'
+! include 'COMMON.MD'
+! Local variables
+ real(kind=8) :: emin
+ integer :: i,j,imin,ierr
+ integer :: diff,allnss,newnss
+ integer,dimension(maxdim) :: allflag,allihpb,alljhpb,& !(maxdim)(maxdim=(maxres-1)*(maxres-2)/2)
+ newihpb,newjhpb
+ logical :: found
+ integer,dimension(0:nfgtasks) :: i_newnss
+ integer,dimension(0:nfgtasks) :: displ
+ integer,dimension(maxdim) :: g_newihpb,g_newjhpb !(maxdim)(maxdim=(maxres-1)*(maxres-2)/2)
+ integer :: g_newnss
+
+ allnss=0
+ do i=1,nres-1
+ do j=i+1,nres
+ if (dyn_ssbond_ij(i,j).lt.1.0d300) then
+ allnss=allnss+1
+ allflag(allnss)=0
+ allihpb(allnss)=i
+ alljhpb(allnss)=j
+ endif
+ enddo
+ enddo
+
+!mc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ 1 emin=1.0d300
+ do i=1,allnss
+ if (allflag(i).eq.0 .and. &
+ dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
+ emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
+ imin=i
+ endif
+ enddo
+ if (emin.lt.1.0d300) then
+ allflag(imin)=1
+ do i=1,allnss
+ if (allflag(i).eq.0 .and. &
+ (allihpb(i).eq.allihpb(imin) .or. &
+ alljhpb(i).eq.allihpb(imin) .or. &
+ allihpb(i).eq.alljhpb(imin) .or. &
+ alljhpb(i).eq.alljhpb(imin))) then
+ allflag(i)=-1
+ endif
+ enddo
+ goto 1
+ endif
+
+!mc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ newnss=0
+ do i=1,allnss
+ if (allflag(i).eq.1) then
+ newnss=newnss+1
+ newihpb(newnss)=allihpb(i)
+ newjhpb(newnss)=alljhpb(i)
+ endif
+ enddo
+
+#ifdef MPI
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(newnss,g_newnss,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+ call MPI_Gather(newnss,1,MPI_INTEGER,&
+ i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_newnss(i-1)+displ(i-1)
+ enddo
+ call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,&
+ g_newihpb,i_newnss,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,&
+ g_newjhpb,i_newnss,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ if(fg_rank.eq.0) then
+! print *,'g_newnss',g_newnss
+! print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
+! print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
+ newnss=g_newnss
+ do i=1,newnss
+ newihpb(i)=g_newihpb(i)
+ newjhpb(i)=g_newjhpb(i)
+ enddo
+ endif
+ endif
+#endif
+
+ diff=newnss-nss
+
+!mc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
+! print *,newnss,nss,maxdim
+ do i=1,nss
+ found=.false.
+! print *,newnss
+ do j=1,newnss
+!! print *,j
+ if (idssb(i).eq.newihpb(j) .and. &
+ jdssb(i).eq.newjhpb(j)) found=.true.
+ enddo
+#ifndef CLUST
+#ifndef WHAM
+! write(iout,*) "found",found,i,j
+ if (.not.found.and.fg_rank.eq.0) &
+ write(iout,'(a15,f12.2,f8.1,2i5)') &
+ "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
+#endif
+#endif
+ enddo
+
+ do i=1,newnss
+ found=.false.
+ do j=1,nss
+! print *,i,j
+ if (newihpb(i).eq.idssb(j) .and. &
+ newjhpb(i).eq.jdssb(j)) found=.true.
+ enddo
+#ifndef CLUST
+#ifndef WHAM
+! write(iout,*) "found",found,i,j
+ if (.not.found.and.fg_rank.eq.0) &
+ write(iout,'(a15,f12.2,f8.1,2i5)') &
+ "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
+#endif
+#endif
+ enddo
+
+ nss=newnss
+ do i=1,nss
+ idssb(i)=newihpb(i)
+ jdssb(i)=newjhpb(i)
+ enddo
+
+ return
+ end subroutine dyn_set_nss
+! 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
+
+!-----------------------------------------------------------------------------
+#ifdef WHAM
+ subroutine read_ssHist
+! implicit none
+! Includes
+! include 'DIMENSIONS'
+! include "DIMENSIONS.FREE"
+! include 'COMMON.FREE'
+! Local variables
+ integer :: i,j
+ character(len=80) :: controlcard
+
+ do i=1,dyn_nssHist
+ call card_concat(controlcard,.true.)
+ read(controlcard,*) &
+ dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
+ enddo
+
+ return
+ end subroutine read_ssHist
+#endif
+!-----------------------------------------------------------------------------
+ integer function indmat(i,j)
+!el
+! get the position of the jth ijth fragment of the chain coordinate system
+! in the fromto array.
+ integer :: i,j
+
+ indmat=((2*(nres-2)-i)*(i-1))/2+j-1
+ return
+ end function indmat
+!-----------------------------------------------------------------------------
+ real(kind=8) function sigm(x)
+!el
+ real(kind=8) :: x
+ sigm=0.25d0*x
+ return
+ end function sigm
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ subroutine alloc_ener_arrays
+!EL Allocation of arrays used by module energy
+ use MD_data, only: mset
+!el local variables
+ integer :: i,j
+
+ if(nres.lt.100) then
+ maxconts=nres
+ elseif(nres.lt.200) then
+ maxconts=0.8*nres ! Max. number of contacts per residue
+ else
+ maxconts=0.6*nres ! (maxconts=maxres/4)
+ endif
+ maxcont=12*nres ! Max. number of SC contacts
+ maxvar=6*nres ! Max. number of variables
+!el maxdim=(nres-1)*(nres-2)/2 ! Max. number of derivatives of virtual-bond
+ maxdim=nres*(nres-2)/2 ! Max. number of derivatives of virtual-bond
+!----------------------
+! arrays in subroutine init_int_table
+!el#ifdef MPI
+!el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
+!el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
+!el#endif
+ allocate(nint_gr(nres))
+ allocate(nscp_gr(nres))
+ allocate(ielstart(nres))
+ allocate(ielend(nres))
+!(maxres)
+ allocate(istart(nres,maxint_gr))
+ allocate(iend(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(iscpstart(nres,maxint_gr))
+ allocate(iscpend(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(ielstart_vdw(nres))
+ allocate(ielend_vdw(nres))
+!(maxres)
+ allocate(nint_gr_nucl(nres))
+ allocate(nscp_gr_nucl(nres))
+ allocate(ielstart_nucl(nres))
+ allocate(ielend_nucl(nres))
+!(maxres)
+ allocate(istart_nucl(nres,maxint_gr))
+ allocate(iend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(iscpstart_nucl(nres,maxint_gr))
+ allocate(iscpend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(ielstart_vdw_nucl(nres))
+ allocate(ielend_vdw_nucl(nres))
+
+ allocate(lentyp(0:nfgtasks-1))
+!(0:maxprocs-1)
+!----------------------
+! commom.contacts
+! common /contacts/
+ if(.not.allocated(icont_ref)) allocate(icont_ref(2,maxcont))
+ allocate(icont(2,maxcont))
+!(2,maxcont)
+! common /contacts1/
+ allocate(num_cont(0:nres+4))
+!(maxres)
+ allocate(jcont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(facont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(gacont(3,maxconts,nres))
+!(3,maxconts,maxres)
+! common /contacts_hb/
+ allocate(gacontp_hb1(3,maxconts,nres))
+ allocate(gacontp_hb2(3,maxconts,nres))
+ allocate(gacontp_hb3(3,maxconts,nres))
+ allocate(gacontm_hb1(3,maxconts,nres))
+ allocate(gacontm_hb2(3,maxconts,nres))
+ allocate(gacontm_hb3(3,maxconts,nres))
+ allocate(gacont_hbr(3,maxconts,nres))
+ allocate(grij_hb_cont(3,maxconts,nres))
+!(3,maxconts,maxres)
+ allocate(facont_hb(maxconts,nres))
+
+ allocate(ees0p(maxconts,nres))
+ allocate(ees0m(maxconts,nres))
+ allocate(d_cont(maxconts,nres))
+ allocate(ees0plist(maxconts,nres))
+
+!(maxconts,maxres)
+ allocate(num_cont_hb(nres))
+!(maxres)
+ allocate(jcont_hb(maxconts,nres))
+!(maxconts,maxres)
+! common /rotat/
+ allocate(Ug(2,2,nres))
+ allocate(Ugder(2,2,nres))
+ allocate(Ug2(2,2,nres))
+ allocate(Ug2der(2,2,nres))
+!(2,2,maxres)
+ allocate(obrot(2,nres))
+ allocate(obrot2(2,nres))
+ allocate(obrot_der(2,nres))
+ allocate(obrot2_der(2,nres))
+!(2,maxres)
+! common /precomp1/
+ allocate(mu(2,nres))
+ allocate(muder(2,nres))
+ allocate(Ub2(2,nres))
+ Ub2(1,:)=0.0d0
+ Ub2(2,:)=0.0d0
+ allocate(Ub2der(2,nres))
+ allocate(Ctobr(2,nres))
+ allocate(Ctobrder(2,nres))
+ allocate(Dtobr2(2,nres))
+ allocate(Dtobr2der(2,nres))
+!(2,maxres)
+ allocate(EUg(2,2,nres))
+ allocate(EUgder(2,2,nres))
+ allocate(CUg(2,2,nres))
+ allocate(CUgder(2,2,nres))
+ allocate(DUg(2,2,nres))
+ allocate(Dugder(2,2,nres))
+ allocate(DtUg2(2,2,nres))
+ allocate(DtUg2der(2,2,nres))
+!(2,2,maxres)
+! common /precomp2/
+ allocate(Ug2Db1t(2,nres))
+ allocate(Ug2Db1tder(2,nres))
+ allocate(CUgb2(2,nres))
+ allocate(CUgb2der(2,nres))
+!(2,maxres)
+ allocate(EUgC(2,2,nres))
+ allocate(EUgCder(2,2,nres))
+ allocate(EUgD(2,2,nres))
+ allocate(EUgDder(2,2,nres))
+ allocate(DtUg2EUg(2,2,nres))
+ allocate(Ug2DtEUg(2,2,nres))
+!(2,2,maxres)
+ allocate(Ug2DtEUgder(2,2,2,nres))
+ allocate(DtUg2EUgder(2,2,2,nres))
+!(2,2,2,maxres)
+! common /rotat_old/
+ allocate(costab(nres))
+ allocate(sintab(nres))
+ allocate(costab2(nres))
+ allocate(sintab2(nres))
+!(maxres)
+! common /dipmat/
+ allocate(a_chuj(2,2,maxconts,nres))
+!(2,2,maxconts,maxres)(maxconts=maxres/4)
+ allocate(a_chuj_der(2,2,3,5,maxconts,nres))
+!(2,2,3,5,maxconts,maxres)(maxconts=maxres/4)
+! common /contdistrib/
+ allocate(ncont_sent(nres))
+ allocate(ncont_recv(nres))
+
+ allocate(iat_sent(nres))
+!(maxres)
+ allocate(iint_sent(4,nres,nres))
+ allocate(iint_sent_local(4,nres,nres))
+!(4,maxres,maxres)
+ allocate(iturn3_sent(4,0:nres+4))
+ allocate(iturn4_sent(4,0:nres+4))
+ allocate(iturn3_sent_local(4,nres))
+ allocate(iturn4_sent_local(4,nres))
+!(4,maxres)
+ allocate(itask_cont_from(0:nfgtasks-1))
+ allocate(itask_cont_to(0:nfgtasks-1))
+!(0:max_fg_procs-1)
+
+
+
+!----------------------
+! commom.deriv;
+! common /derivat/
+ allocate(dcdv(6,maxdim))
+ allocate(dxdv(6,maxdim))
+!(6,maxdim)
+ allocate(dxds(6,nres))
+!(6,maxres)
+ allocate(gradx(3,-1:nres,0:2))
+ allocate(gradc(3,-1:nres,0:2))
+!(3,maxres,2)
+ allocate(gvdwx(3,-1:nres))
+ allocate(gvdwc(3,-1:nres))
+ allocate(gelc(3,-1:nres))
+ allocate(gelc_long(3,-1:nres))
+ allocate(gvdwpp(3,-1:nres))
+ allocate(gvdwc_scpp(3,-1:nres))
+ allocate(gradx_scp(3,-1:nres))
+ allocate(gvdwc_scp(3,-1:nres))
+ allocate(ghpbx(3,-1:nres))
+ allocate(ghpbc(3,-1:nres))
+ allocate(gradcorr(3,-1:nres))
+ allocate(gradcorr_long(3,-1:nres))
+ allocate(gradcorr5_long(3,-1:nres))
+ allocate(gradcorr6_long(3,-1:nres))
+ allocate(gcorr6_turn_long(3,-1:nres))
+ allocate(gradxorr(3,-1:nres))
+ allocate(gradcorr5(3,-1:nres))
+ allocate(gradcorr6(3,-1:nres))
+ allocate(gliptran(3,-1:nres))
+ allocate(gliptranc(3,-1:nres))
+ allocate(gliptranx(3,-1:nres))
+ allocate(gshieldx(3,-1:nres))
+ allocate(gshieldc(3,-1:nres))
+ allocate(gshieldc_loc(3,-1:nres))
+ allocate(gshieldx_ec(3,-1:nres))
+ allocate(gshieldc_ec(3,-1:nres))
+ allocate(gshieldc_loc_ec(3,-1:nres))
+ allocate(gshieldx_t3(3,-1:nres))
+ allocate(gshieldc_t3(3,-1:nres))
+ allocate(gshieldc_loc_t3(3,-1:nres))
+ allocate(gshieldx_t4(3,-1:nres))
+ allocate(gshieldc_t4(3,-1:nres))
+ allocate(gshieldc_loc_t4(3,-1:nres))
+ allocate(gshieldx_ll(3,-1:nres))
+ allocate(gshieldc_ll(3,-1:nres))
+ allocate(gshieldc_loc_ll(3,-1:nres))
+ allocate(grad_shield(3,-1:nres))
+ allocate(gg_tube_sc(3,-1:nres))
+ allocate(gg_tube(3,-1:nres))
+ allocate(gradafm(3,-1:nres))
+ allocate(gradb_nucl(3,-1:nres))
+ allocate(gradbx_nucl(3,-1:nres))
+ allocate(gvdwpsb1(3,-1:nres))
+ allocate(gelpp(3,-1:nres))
+ allocate(gvdwpsb(3,-1:nres))
+ allocate(gelsbc(3,-1:nres))
+ allocate(gelsbx(3,-1:nres))
+ allocate(gvdwsbx(3,-1:nres))
+ allocate(gvdwsbc(3,-1:nres))
+ allocate(gsbloc(3,-1:nres))
+ allocate(gsblocx(3,-1:nres))
+ allocate(gradcorr_nucl(3,-1:nres))
+ allocate(gradxorr_nucl(3,-1:nres))
+ allocate(gradcorr3_nucl(3,-1:nres))
+ allocate(gradxorr3_nucl(3,-1:nres))
+
+!(3,maxres)
+ allocate(grad_shield_side(3,50,nres))
+ allocate(grad_shield_loc(3,50,nres))
+! grad for shielding surroing
+ allocate(gloc(0:maxvar,0:2))
+ allocate(gloc_x(0:maxvar,2))
+!(maxvar,2)
+ allocate(gel_loc(3,-1:nres))
+ allocate(gel_loc_long(3,-1:nres))
+ allocate(gcorr3_turn(3,-1:nres))
+ allocate(gcorr4_turn(3,-1:nres))
+ allocate(gcorr6_turn(3,-1:nres))
+ allocate(gradb(3,-1:nres))
+ allocate(gradbx(3,-1:nres))
+!(3,maxres)
+ allocate(gel_loc_loc(maxvar))
+ allocate(gel_loc_turn3(maxvar))
+ allocate(gel_loc_turn4(maxvar))
+ allocate(gel_loc_turn6(maxvar))
+ allocate(gcorr_loc(maxvar))
+ allocate(g_corr5_loc(maxvar))
+ allocate(g_corr6_loc(maxvar))
+!(maxvar)
+ allocate(gsccorc(3,-1:nres))
+ allocate(gsccorx(3,-1:nres))
+!(3,maxres)
+ allocate(gsccor_loc(-1:nres))
+!(maxres)
+ allocate(dtheta(3,2,-1:nres))
+!(3,2,maxres)
+ allocate(gscloc(3,-1:nres))
+ allocate(gsclocx(3,-1:nres))
+!(3,maxres)
+ allocate(dphi(3,3,-1:nres))
+ allocate(dalpha(3,3,-1:nres))
+ allocate(domega(3,3,-1:nres))
+!(3,3,maxres)
+! common /deriv_scloc/
+ allocate(dXX_C1tab(3,nres))
+ allocate(dYY_C1tab(3,nres))
+ allocate(dZZ_C1tab(3,nres))
+ allocate(dXX_Ctab(3,nres))
+ allocate(dYY_Ctab(3,nres))
+ allocate(dZZ_Ctab(3,nres))
+ allocate(dXX_XYZtab(3,nres))
+ allocate(dYY_XYZtab(3,nres))
+ allocate(dZZ_XYZtab(3,nres))
+!(3,maxres)
+! common /mpgrad/
+ allocate(jgrad_start(nres))
+ allocate(jgrad_end(nres))
+!(maxres)
+!----------------------
+
+! common /indices/
+ allocate(ibond_displ(0:nfgtasks-1))
+ allocate(ibond_count(0:nfgtasks-1))
+ allocate(ithet_displ(0:nfgtasks-1))
+ allocate(ithet_count(0:nfgtasks-1))
+ allocate(iphi_displ(0:nfgtasks-1))
+ allocate(iphi_count(0:nfgtasks-1))
+ allocate(iphi1_displ(0:nfgtasks-1))
+ allocate(iphi1_count(0:nfgtasks-1))
+ allocate(ivec_displ(0:nfgtasks-1))
+ allocate(ivec_count(0:nfgtasks-1))
+ allocate(iset_displ(0:nfgtasks-1))
+ allocate(iset_count(0:nfgtasks-1))
+ allocate(iint_count(0:nfgtasks-1))
+ allocate(iint_displ(0:nfgtasks-1))
+!(0:max_fg_procs-1)
+!----------------------
+! common.MD
+! common /mdgrad/
+ allocate(gcart(3,-1:nres))
+ allocate(gxcart(3,-1:nres))
+!(3,0:MAXRES)
+ allocate(gradcag(3,-1:nres))
+ allocate(gradxag(3,-1:nres))
+!(3,MAXRES)
+! common /back_constr/
+!el in energy:Econstr_back allocate((:),allocatable :: utheta,ugamma,uscdiff !(maxfrag_back)
+ allocate(dutheta(nres))
+ allocate(dugamma(nres))
+!(maxres)
+ allocate(duscdiff(3,nres))
+ allocate(duscdiffx(3,nres))
+!(3,maxres)
+!el i io:read_fragments
+! allocate((:,:,:),allocatable :: wfrag_back !(3,maxfrag_back,maxprocs/20)
+! allocate((:,:,:),allocatable :: ifrag_back !(3,maxfrag_back,maxprocs/20)
+! common /qmeas/
+! allocate(qinfrag(50,nprocs/20),wfrag(50,nprocs/20)) !(50,maxprocs/20)
+! allocate(qinpair(100,nprocs/20),wpair(100,nprocs/20)) !(100,maxprocs/20)
+ allocate(mset(0:nprocs)) !(maxprocs/20)
+ mset(:)=0
+! allocate(ifrag(2,50,nprocs/20)) !(2,50,maxprocs/20)
+! allocate(ipair(2,100,nprocs/20)) !(2,100,maxprocs/20)
+ allocate(dUdconst(3,0:nres))
+ allocate(dUdxconst(3,0:nres))
+ allocate(dqwol(3,0:nres))
+ allocate(dxqwol(3,0:nres))
+!(3,0:MAXRES)
+!----------------------
+! common.sbridge
+! common /sbridge/ in io_common: read_bridge
+!el allocate((:),allocatable :: iss !(maxss)
+! common /links/ in io_common: read_bridge
+!el real(kind=8),dimension(:),allocatable :: dhpb,forcon,dhpb1 !(maxdim) !el dhpb1 !!! nie używane
+!el integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane
+! common /dyn_ssbond/
+! and side-chain vectors in theta or phi.
+ allocate(dyn_ssbond_ij(0:nres+4,0:nres+4))
+!(maxres,maxres)
+! do i=1,nres
+! do j=i+1,nres
+ dyn_ssbond_ij(:,:)=1.0d300
+! enddo
+! enddo
+
+! if (nss.gt.0) then
+ allocate(idssb(maxdim),jdssb(maxdim))
+! allocate(newihpb(nss),newjhpb(nss))
+!(maxdim)
+! endif
+ allocate(ishield_list(nres))
+ allocate(shield_list(50,nres))
+ allocate(dyn_ss_mask(nres))
+ allocate(fac_shield(nres))
+ allocate(enetube(nres*2))
+ allocate(enecavtube(nres*2))
+
+!(maxres)
+ dyn_ss_mask(:)=.false.
+!----------------------
+! common.sccor
+! Parameters of the SCCOR term
+! common/sccor/
+!el in io_conf: parmread
+! allocate(v1sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp))
+! allocate(v2sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
+! allocate(v0sccor(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
+! allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
+! allocate(nterm_sccor(-ntyp:ntyp,-ntyp:ntyp))
+! allocate(nlor_sccor(-ntyp:ntyp,-ntyp:ntyp)) !(-ntyp:ntyp,-ntyp:ntyp)
+! allocate(vlor1sccor(maxterm_sccor,20,20))
+! allocate(vlor2sccor(maxterm_sccor,20,20))
+! allocate(vlor3sccor(maxterm_sccor,20,20)) !(maxterm_sccor,20,20)
+!----------------
+ allocate(gloc_sc(3,0:2*nres,0:10))
+!(3,0:maxres2,10)maxres2=2*maxres
+ allocate(dcostau(3,3,3,2*nres))
+ allocate(dsintau(3,3,3,2*nres))
+ allocate(dtauangle(3,3,3,2*nres))
+ allocate(dcosomicron(3,3,3,2*nres))
+ allocate(domicron(3,3,3,2*nres))
+!(3,3,3,maxres2)maxres2=2*maxres
+!----------------------
+! common.var
+! common /restr/
+ allocate(varall(maxvar))
+!(maxvar)(maxvar=6*maxres)
+ allocate(mask_theta(nres))
+ allocate(mask_phi(nres))
+ allocate(mask_side(nres))
+!(maxres)
+!----------------------
+! common.vectors
+! common /vectors/
+ allocate(uy(3,nres))
+ allocate(uz(3,nres))
+!(3,maxres)
+ allocate(uygrad(3,3,2,nres))
+ allocate(uzgrad(3,3,2,nres))
+!(3,3,2,maxres)
+
+ return
+ end subroutine alloc_ener_arrays
+!-----------------------------------------------------------------
+ subroutine ebond_nucl(estr_nucl)
+!c
+!c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds
+!c
+
+ real(kind=8),dimension(3) :: u,ud
+ real(kind=8) :: usum,uprod,uprod1,uprod2,usumsqder
+ real(kind=8) :: estr_nucl,diff
+ integer :: iti,i,j,k,nbi
+ estr_nucl=0.0d0
+!C print *,"I enter ebond"
+ if (energy_dec) &
+ write (iout,*) "ibondp_start,ibondp_end",&
+ ibondp_nucl_start,ibondp_nucl_end
+ do i=ibondp_nucl_start,ibondp_nucl_end
+ if (itype(i-1,2).eq.ntyp1_molec(2) .or. &
+ itype(i,2).eq.ntyp1_molec(2)) cycle
+! estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+! do j=1,3
+! gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
+! & *dc(j,i-1)/vbld(i)
+! enddo
+! if (energy_dec) write(iout,*)
+! & "estr1",i,vbld(i),distchainmax,
+! & gnmr1(vbld(i),-1.0d0,distchainmax)
+
+ diff = vbld(i)-vbldp0_nucl
+ if(energy_dec)write(iout,*) "estr_nucl_bb" , i,vbld(i),&
+ vbldp0_nucl,diff,AKP_nucl*diff*diff
+ estr_nucl=estr_nucl+diff*diff
+ print *,estr_nucl
+ do j=1,3
+ gradb_nucl(j,i-1)=AKP_nucl*diff*dc(j,i-1)/vbld(i)
+ enddo
+!c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
+ enddo
+ estr_nucl=0.5d0*AKP_nucl*estr_nucl
+ print *,"partial sum", estr_nucl,AKP_nucl
+
+ if (energy_dec) &
+ write (iout,*) "ibondp_start,ibondp_end",&
+ ibond_nucl_start,ibond_nucl_end
+
+ do i=ibond_nucl_start,ibond_nucl_end
+!C print *, "I am stuck",i
+ iti=itype(i,2)
+ if (iti.eq.ntyp1_molec(2)) cycle
+ nbi=nbondterm_nucl(iti)
+!C print *,iti,nbi
+ if (nbi.eq.1) then
+ diff=vbld(i+nres)-vbldsc0_nucl(1,iti)
+
+ if (energy_dec) &
+ write (iout,*) "estr_nucl_sc", i,iti,vbld(i+nres),vbldsc0_nucl(1,iti),diff, &
+ AKSC_nucl(1,iti),AKSC_nucl(1,iti)*diff*diff
+ estr_nucl=estr_nucl+0.5d0*AKSC_nucl(1,iti)*diff*diff
+ print *,estr_nucl
+ do j=1,3
+ gradbx_nucl(j,i)=AKSC_nucl(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
+ enddo
+ else
+ do j=1,nbi
+ diff=vbld(i+nres)-vbldsc0_nucl(j,iti)
+ ud(j)=aksc_nucl(j,iti)*diff
+ u(j)=abond0_nucl(j,iti)+0.5d0*ud(j)*diff
+ enddo
+ uprod=u(1)
+ do j=2,nbi
+ uprod=uprod*u(j)
+ enddo
+ usum=0.0d0
+ usumsqder=0.0d0
+ do j=1,nbi
+ uprod1=1.0d0
+ uprod2=1.0d0
+ do k=1,nbi
+ if (k.ne.j) then
+ uprod1=uprod1*u(k)
+ uprod2=uprod2*u(k)*u(k)
+ endif
+ enddo
+ usum=usum+uprod1
+ usumsqder=usumsqder+ud(j)*uprod2
+ enddo
+ estr_nucl=estr_nucl+uprod/usum
+ do j=1,3
+ gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
+ enddo
+ endif
+ enddo
+!C print *,"I am about to leave ebond"
+ return
+ end subroutine ebond_nucl
+
+!-----------------------------------------------------------------------------
+ subroutine ebend_nucl(etheta_nucl)
+ real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm
+ real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle
+ real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble
+ logical :: lprn=.true., lprn1=.false.
+!el local variables
+ integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
+ real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
+ real(kind=8) :: aux,etheta_nucl,ccl,ssl,scl,csl,ethetacnstr
+! local variables for constrains
+ real(kind=8) :: difi,thetiii
+ integer itheta
+ etheta_nucl=0.0D0
+! print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres
+ do i=ithet_nucl_start,ithet_nucl_end
+ if ((itype(i-1,2).eq.ntyp1_molec(2)).or.&
+ (itype(i-2,2).eq.ntyp1_molec(2)).or. &
+ (itype(i,2).eq.ntyp1_molec(2))) cycle
+ dethetai=0.0d0
+ dephii=0.0d0
+ dephii1=0.0d0
+ theti2=0.5d0*theta(i)
+ ityp2=ithetyp_nucl(itype(i-1,2))
+ do k=1,nntheterm_nucl
+ coskt(k)=dcos(k*theti2)
+ sinkt(k)=dsin(k*theti2)
+ enddo
+ if (i.gt.3 .and. itype(i-2,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+ phii=phi(i)
+ if (phii.ne.phii) phii=150.0
+#else
+ phii=phi(i)
+#endif
+ ityp1=ithetyp_nucl(itype(i-2,2))
+ do k=1,nsingle_nucl
+ cosph1(k)=dcos(k*phii)
+ sinph1(k)=dsin(k*phii)
+ enddo
+ else
+ phii=0.0d0
+ ityp1=nthetyp_nucl+1
+ do k=1,nsingle_nucl
+ cosph1(k)=0.0d0
+ sinph1(k)=0.0d0
+ enddo
+ endif
+
+ if (i.lt.nres .and. itype(i,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+ phii1=phi(i+1)
+ if (phii1.ne.phii1) phii1=150.0
+ phii1=pinorm(phii1)
+#else
+ phii1=phi(i+1)
+#endif
+ ityp3=ithetyp_nucl(itype(i,2))
+ do k=1,nsingle_nucl
+ cosph2(k)=dcos(k*phii1)
+ sinph2(k)=dsin(k*phii1)
+ enddo
+ else
+ phii1=0.0d0
+ ityp3=nthetyp_nucl+1
+ do k=1,nsingle_nucl
+ cosph2(k)=0.0d0
+ sinph2(k)=0.0d0
+ enddo
+ endif
+ ethetai=aa0thet_nucl(ityp1,ityp2,ityp3)
+ do k=1,ndouble_nucl
+ do l=1,k-1
+ ccl=cosph1(l)*cosph2(k-l)
+ ssl=sinph1(l)*sinph2(k-l)
+ scl=sinph1(l)*cosph2(k-l)
+ csl=cosph1(l)*sinph2(k-l)
+ cosph1ph2(l,k)=ccl-ssl
+ cosph1ph2(k,l)=ccl+ssl
+ sinph1ph2(l,k)=scl+csl
+ sinph1ph2(k,l)=scl-csl
+ enddo
+ enddo
+ if (lprn) then
+ write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,&
+ " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
+ write (iout,*) "coskt and sinkt",nntheterm_nucl
+ do k=1,nntheterm_nucl
+ write (iout,*) k,coskt(k),sinkt(k)
+ enddo
+ endif
+ do k=1,ntheterm_nucl
+ ethetai=ethetai+aathet_nucl(k,ityp1,ityp2,ityp3)*sinkt(k)
+ dethetai=dethetai+0.5d0*k*aathet_nucl(k,ityp1,ityp2,ityp3)&
+ *coskt(k)
+ if (lprn)&
+ write (iout,*) "k",k," aathet",aathet_nucl(k,ityp1,ityp2,ityp3),&
+ " ethetai",ethetai
+ enddo
+ if (lprn) then
+ write (iout,*) "cosph and sinph"
+ do k=1,nsingle_nucl
+ write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k)
+ enddo
+ write (iout,*) "cosph1ph2 and sinph2ph2"
+ do k=2,ndouble_nucl
+ do l=1,k-1
+ write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),&
+ sinph1ph2(l,k),sinph1ph2(k,l)
+ enddo
+ enddo
+ write(iout,*) "ethetai",ethetai
+ endif
+ do m=1,ntheterm2_nucl
+ do k=1,nsingle_nucl
+ aux=bbthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)&
+ +ccthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)&
+ +ddthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)&
+ +eethet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+ ethetai=ethetai+sinkt(m)*aux
+ dethetai=dethetai+0.5d0*m*aux*coskt(m)
+ dephii=dephii+k*sinkt(m)*(&
+ ccthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)-&
+ bbthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+ dephii1=dephii1+k*sinkt(m)*(&
+ eethet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)-&
+ ddthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+ if (lprn) &
+ write (iout,*) "m",m," k",k," bbthet",&
+ bbthet_nucl(k,m,ityp1,ityp2,ityp3)," ccthet",&
+ ccthet_nucl(k,m,ityp1,ityp2,ityp3)," ddthet",&
+ ddthet_nucl(k,m,ityp1,ityp2,ityp3)," eethet",&
+ eethet_nucl(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ enddo
+ enddo
+ if (lprn) &
+ write(iout,*) "ethetai",ethetai
+ do m=1,ntheterm3_nucl
+ do k=2,ndouble_nucl
+ do l=1,k-1
+ aux=ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+ ethetai=ethetai+sinkt(m)*aux
+ dethetai=dethetai+0.5d0*m*coskt(m)*aux
+ dephii=dephii+l*sinkt(m)*(&
+ -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ dephii1=dephii1+(k-l)*sinkt(m)*( &
+ -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ if (lprn) then
+ write (iout,*) "m",m," k",k," l",l," ffthet", &
+ ffthet_nucl(l,k,m,ityp1,ityp2,ityp3), &
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ggthet",&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3),&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ write (iout,*) cosph1ph2(l,k)*sinkt(m), &
+ cosph1ph2(k,l)*sinkt(m),&
+ sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
+ endif
+ enddo
+ enddo
+ enddo
+10 continue
+ if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') &
+ i,theta(i)*rad2deg,phii*rad2deg, &
+ phii1*rad2deg,ethetai
+ etheta_nucl=etheta_nucl+ethetai
+! print *,i,"partial sum",etheta_nucl
+ if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang_nucl*dephii
+ if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang_nucl*dephii1
+ gloc(nphi+i-2,icg)=wang_nucl*dethetai
+ enddo
+ return
+ end subroutine ebend_nucl
+!----------------------------------------------------
+ subroutine etor_nucl(etors_nucl)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.TORSION'
+! include 'COMMON.INTERACT'
+! include 'COMMON.DERIV'
+! include 'COMMON.CHAIN'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.FFIELD'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+ real(kind=8) :: etors_nucl,edihcnstr
+ logical :: lprn
+!el local variables
+ integer :: i,j,iblock,itori,itori1
+ real(kind=8) :: phii,gloci,v1ij,v2ij,cosphi,sinphi,&
+ vl1ij,vl2ij,vl3ij,pom1,difi,etors_ii,pom
+! Set lprn=.true. for debugging
+ lprn=.false.
+! lprn=.true.
+ etors_nucl=0.0D0
+! print *,"iphi_nucl_start/end", iphi_nucl_start,iphi_nucl_end
+ do i=iphi_nucl_start,iphi_nucl_end
+ if (itype(i-2,2).eq.ntyp1_molec(2) .or. itype(i-1,2).eq.ntyp1_molec(2) &
+ .or. itype(i-3,2).eq.ntyp1_molec(2) &
+ .or. itype(i,2).eq.ntyp1_molec(2)) cycle
+ etors_ii=0.0D0
+ itori=itortyp_nucl(itype(i-2,2))
+ itori1=itortyp_nucl(itype(i-1,2))
+ phii=phi(i)
+! print *,i,itori,itori1
+ gloci=0.0D0
+!C Regular cosine and sine terms
+ do j=1,nterm_nucl(itori,itori1)
+ v1ij=v1_nucl(j,itori,itori1)
+ v2ij=v2_nucl(j,itori,itori1)
+ cosphi=dcos(j*phii)
+ sinphi=dsin(j*phii)
+ etors_nucl=etors_nucl+v1ij*cosphi+v2ij*sinphi
+ if (energy_dec) etors_ii=etors_ii+&
+ v1ij*cosphi+v2ij*sinphi
+ gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
+ enddo
+!C Lorentz terms
+!C v1
+!C E = SUM ----------------------------------- - v1
+!C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1
+!C
+ cosphi=dcos(0.5d0*phii)
+ sinphi=dsin(0.5d0*phii)
+ do j=1,nlor_nucl(itori,itori1)
+ vl1ij=vlor1_nucl(j,itori,itori1)
+ vl2ij=vlor2_nucl(j,itori,itori1)
+ vl3ij=vlor3_nucl(j,itori,itori1)
+ pom=vl2ij*cosphi+vl3ij*sinphi
+ pom1=1.0d0/(pom*pom+1.0d0)
+ etors_nucl=etors_nucl+vl1ij*pom1
+ if (energy_dec) etors_ii=etors_ii+ &
+ vl1ij*pom1
+ pom=-pom*pom1*pom1
+ gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
+ enddo
+!C Subtract the constant term
+ etors_nucl=etors_nucl-v0_nucl(itori,itori1)
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
+ 'etor',i,etors_ii-v0_nucl(itori,itori1)
+ if (lprn) &
+ write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
+ restyp(itype(i-2,2),2),i-2,restyp(itype(i-1,2),2),i-1,itori,itori1, &
+ (v1_nucl(j,itori,itori1),j=1,6),(v2_nucl(j,itori,itori1),j=1,6)
+ gloc(i-3,icg)=gloc(i-3,icg)+wtor_nucl*gloci
+!c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
+ enddo
+ return
+ end subroutine etor_nucl
+!------------------------------------------------------------
+ subroutine epp_nucl_sub(evdw1,ees)
+!C
+!C This subroutine calculates the average interaction energy and its gradient
+!C in the virtual-bond vectors between non-adjacent peptide groups, based on
+!C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
+!C The potential depends both on the distance of peptide-group centers and on
+!C the orientation of the CA-CA virtual bonds.
+!C
+ integer :: i,j,k,iteli,itelj,num_conti,isubchap,ind
+ real(kind=8) :: dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
+ real(kind=8) :: xj,yj,zj,rij,rrmij,sss,r3ij,r6ij,evdw1,&
+ dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
+ dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,sss_grad,fac,evdw1ij
+ integer xshift,yshift,zshift
+ real(kind=8),dimension(3):: ggg,gggp,gggm,erij
+ real(kind=8) :: ees,eesij
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ ind=0
+!c
+!c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
+!c
+ print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
+ do i=iatel_s_nucl,iatel_e_nucl
+ if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+
+ do j=ielstart_nucl(i),ielend_nucl(i)
+ if (itype(j,2).eq.ntyp1_molec(2) .or. itype(j+1,2).eq.ntyp1_molec(2)) cycle
+ ind=ind+1
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+! xj=c(1,j)+0.5D0*dxj-xmedi
+! yj=c(2,j)+0.5D0*dyj-ymedi
+! zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ isubchap=0
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+!C print *,i,j
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
+
+ rij=xj*xj+yj*yj+zj*zj
+!c write (2,*)"ij",i,j," r0pp",r0pp," rij",rij," epspp",epspp
+ fac=(r0pp**2/rij)**3
+ ev1=epspp*fac*fac
+ ev2=epspp*fac
+ evdw1ij=ev1-2*ev2
+ fac=(-ev1-evdw1ij)/rij
+! write (2,*)"fac",fac," ev1",ev1," ev2",ev2," evdw1ij",evdw1ij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"evdw1ij",evdw1ij
+ evdw1=evdw1+evdw1ij
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+ do k=1,3
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ enddo
+!c phoshate-phosphate electrostatic interactions
+ rij=dsqrt(rij)
+ fac=1.0d0/rij
+ eesij=dexp(-BEES*rij)*fac
+! write (2,*)"fac",fac," eesijpp",eesij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"eesijpp",eesij
+ ees=ees+eesij
+!c fac=-eesij*fac
+ fac=-(fac+BEES)*eesij*fac
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+!c write(2,*) "ggg",i,j,ggg(1),ggg(2),ggg(3)
+!c write(2,*) "gelpp",i,(gelpp(k,i),k=1,3)
+!c write(2,*) "gelpp",j,(gelpp(k,j),k=1,3)
+ do k=1,3
+ gelpp(k,i)=gelpp(k,i)-ggg(k)
+ gelpp(k,j)=gelpp(k,j)+ggg(k)
+ enddo
+ enddo ! j
+ enddo ! i
+!c ees=332.0d0*ees
+ ees=AEES*ees
+ do i=nnt,nct
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ do k=1,3
+ gvdwpp(k,i)=6*gvdwpp(k,i)
+!c gelpp(k,i)=332.0d0*gelpp(k,i)
+ gelpp(k,i)=AEES*gelpp(k,i)
+ enddo
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ enddo
+!c write (2,*) "total EES",ees