call calctube(Etube)
elseif (TUBElog.eq.2) then
call calctube2(Etube)
+ elseif (TUBElog.eq.3) then
+ call calcnano(Etube)
else
Etube=0.0d0
endif
enddo
enddo
- j=1
- i=0
- print *,"KUPA2",gradbufc(j,i),wsc*gvdwc(j,i),
- & wscp*gvdwc_scp(j,i),gvdwc_scpp(j,i),
- & welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i),
- & wel_loc*gel_loc_long(j,i),
- & wcorr*gradcorr_long(j,i),
- & wcorr5*gradcorr5_long(j,i),
- & wcorr6*gradcorr6_long(j,i),
- & wturn6*gcorr6_turn_long(j,i),
- & wstrain*ghpbc(j,i)
- & ,wliptran*gliptranc(j,i)
- & ,gradafm(j,i)
- & ,welec*gshieldc(j,i)
- & ,wcorr*gshieldc_ec(j,i)
- & ,wturn3*gshieldc_t3(j,i)
- & ,wturn4*gshieldc_t4(j,i)
- & ,wel_loc*gshieldc_ll(j,i)
- & ,wtube*gg_tube(j,i)
+C j=1
+C i=0
+C print *,"KUPA2",gradbufc(j,i),wsc*gvdwc(j,i),
+C & wscp*gvdwc_scp(j,i),gvdwc_scpp(j,i),
+C & welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i),
+C & wel_loc*gel_loc_long(j,i),
+C & wcorr*gradcorr_long(j,i),
+C & wcorr5*gradcorr5_long(j,i),
+C & wcorr6*gradcorr6_long(j,i),
+C & wturn6*gcorr6_turn_long(j,i),
+C & wstrain*ghpbc(j,i)
+C & ,wliptran*gliptranc(j,i)
+C & ,gradafm(j,i)
+C & ,welec*gshieldc(j,i)
+C & ,wcorr*gshieldc_ec(j,i)
+C & ,wturn3*gshieldc_t3(j,i)
+C & ,wturn4*gshieldc_t4(j,i)
+C & ,wel_loc*gshieldc_ll(j,i)
+C & ,wtube*gg_tube(j,i)
#else
do i=0,nct
do j=1,3
C 5) add COMMONs
C 6) add to zerograd
+
+C#-------------------------------------------------------------------------------
+C This subroutine is to mimic the histone like structure but as well can be
+C utilizet to nanostructures (infinit) small modification has to be used to
+C make it finite (z gradient at the ends has to be changes as well as the x,y
+C gradient has to be modified at the ends
+C The energy function is Kihara potential
+C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6)
+C 4eps is depth of well sigma is r_minimum r is distance from center of tube
+C and r0 is the excluded size of nanotube (can be set to 0 if we want just a
+C simple Kihara potential
+ subroutine calcnano(Etube)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+ double precision tub_r,vectube(3),enetube(maxres*2)
+ Etube=0.0d0
+ do i=itube_start,itube_end
+ enetube(i)=0.0d0
+ enetube(i+nres)=0.0d0
+ enddo
+C first we calculate the distance from tube center
+C first sugare-phosphate group for NARES this would be peptide group
+C for UNRES
+ do i=itube_start,itube_end
+C lets ommit dummy atoms for now
+ if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+C now calculate distance from center of tube and direction vectors
+ xmin=boxxsize
+ ymin=boxysize
+ zmin=boxzsize
+
+ do j=-1,1
+ vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ zminact=abs(vectube(3)-tubecenter(3))
+
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ if (zmin.gt.zminact) then
+ zmin=zminact
+ ztemp=vectube(3)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(3)=ztemp
+
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+ vectube(3)=vectube(3)-tubecenter(3)
+
+C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+C as the tube is infinity we do not calculate the Z-vector use of Z
+C as chosen axis
+C vectube(3)=0.0d0
+C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+ vectube(3)=vectube(3)/tub_r
+C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+C and its 6 power
+ rdiff6=rdiff**6.0d0
+C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6
+C write(iout,*) "TU13",i,rdiff6,enetube(i)
+C print *,rdiff,rdiff6,pep_aa_tube
+C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+C now we calculate gradient
+ fac=(-12.0d0*pep_aa_tube/rdiff6-
+ & 6.0d0*pep_bb_tube)/rdiff6/rdiff
+C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+C &rdiff,fac
+
+C now direction of gg_tube vector
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+ gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+ enddo
+ enddo
+
+ do i=itube_start,itube_end
+ enecavtube(i)=0.0
+C Lets not jump over memory as we use many times iti
+ iti=itype(i)
+C lets ommit dummy atoms for now
+ if ((iti.eq.ntyp1)
+C in UNRES uncomment the line below as GLY has no side-chain...
+C .or.(iti.eq.10)
+ & ) cycle
+ xmin=boxxsize
+ ymin=boxysize
+ zmin=boxzsize
+ do j=-1,1
+ vectube(1)=mod((c(1,i+nres)),boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i+nres)),boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ vectube(3)=mod((c(3,i+nres)),boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ zminact=abs(vectube(3)-tubecenter(3))
+
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ if (zmin.gt.zminact) then
+ zmin=zminact
+ ztemp=vectube(3)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(3)=ztemp
+
+C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2),
+C & tubecenter(2)
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+ vectube(3)=vectube(3)-tubecenter(3)
+C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+ vectube(3)=vectube(3)/tub_r
+
+C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+C and its 6 power
+ rdiff6=rdiff**6.0d0
+C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ sc_aa_tube=sc_aa_tube_par(iti)
+ sc_bb_tube=sc_bb_tube_par(iti)
+ enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+C 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
+C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
+ if (acavtub(iti).eq.0.0d0) go to 667
+ denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
+ enecavtube(i)=
+ & acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff)+cavtub(iti))
+ & /denominator
+ faccav=(acavtub(iti)*(1.0+bcavtub(iti)/2.0/sqrt(rdiff))
+ & *denominator-acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff)
+ & +cavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti))
+ & /denominator**2.0d0
+ fac=fac+faccav
+ 667 continue
+ do j=1,3
+ gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ enddo
+C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
+C do i=itube_start,itube_end
+C enecav(i)=0.0
+C iti=itype(i)
+C if (acavtub(iti).eq.0.0) cycle
+
+
+
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)
+ 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
+
projects / unres.git / commitdiff