cmc
c if (dyn_ss) call dyn_set_nss
-c print *,"Processor",myrank," computed USCSC"
+C print *,"Processor",myrank," computed USCSC"
#ifdef TIMING
time01=MPI_Wtime()
#endif
enddo
#endif
endif
-c print *,"Processor",myrank," left VEC_AND_DERIV"
+C print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
#ifdef SPLITELE
if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or.
eello_turn4=0.0d0
endif
else
- write (iout,*) "Soft-spheer ELEC potential"
+C write (iout,*) "Soft-spheer ELEC potential"
c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
c & eello_turn4)
endif
-c print *,"Processor",myrank," computed UELEC"
+C print *,"Processor",myrank," computed UELEC"
C
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
C
C Calculate the disulfide-bridge and other energy and the contributions
C from other distance constraints.
-cd print *,'Calling EHPB'
+C print *,'Calling EHPB'
call edis(ehpb)
-cd print *,'EHPB exitted succesfully.'
+C print *,'EHPB exitted succesfully.'
C
C Calculate the virtual-bond-angle energy.
C
ebe=0
ethetacnstr=0
endif
-c print *,"Processor",myrank," computed UB"
+C print *,"Processor",myrank," computed UB"
C
C Calculate the SC local energy.
C
C print *,"TU DOCHODZE?"
call esc(escloc)
-c print *,"Processor",myrank," computed USC"
+C print *,"Processor",myrank," computed USC"
C
C Calculate the virtual-bond torsional energy.
C
etors=0
edihcnstr=0
endif
-c print *,"Processor",myrank," computed Utor"
+C print *,"Processor",myrank," computed Utor"
C
C 6/23/01 Calculate double-torsional energy
C
else
etors_d=0
endif
-c print *,"Processor",myrank," computed Utord"
+C print *,"Processor",myrank," computed Utord"
C
C 21/5/07 Calculate local sicdechain correlation energy
C
esccor=0.0d0
endif
C print *,"PRZED MULIt"
-c print *,"Processor",myrank," computed Usccorr"
+C print *,"Processor",myrank," computed Usccorr"
C
C 12/1/95 Multi-body terms
C
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
+ enddo
+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
#ifdef TIMING
c time_allreduce=time_allreduce+MPI_Wtime()-time00
#endif
- do i=nnt,nres
+ do i=0,nres
do k=1,3
gradbufc(k,i)=0.0d0
enddo
enddo
- enddo
+ enddo
+C i=0
+C j=1
+C print *,"KUPA", gradbufc(j,i),welec*gelc(j,i),
+C & wel_loc*gel_loc(j,i),
+C & 0.5d0*wscp*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 & wbond*gradb(j,i),
+C & wcorr*gradcorr(j,i),
+C & wturn3*gcorr3_turn(j,i),
+C & wturn4*gcorr4_turn(j,i),
+C & wcorr5*gradcorr5(j,i),
+C & wcorr6*gradcorr6(j,i),
+C & wturn6*gcorr6_turn(j,i),
+C & wsccor*gsccorc(j,i)
+C & ,wscloc*gscloc(j,i)
+C & ,wliptran*gliptranc(j,i)
+C & ,gradafm(j,i)
+C & +welec*gshieldc(j,i)
+C & +welec*gshieldc_loc(j,i)
+C & +wcorr*gshieldc_ec(j,i)
+C & +wcorr*gshieldc_loc_ec(j,i)
+C & +wturn3*gshieldc_t3(j,i)
+C & +wturn3*gshieldc_loc_t3(j,i)
+C & +wturn4*gshieldc_t4(j,i)
+C & ,wturn4*gshieldc_loc_t4(j,i)
+C & ,wel_loc*gshieldc_ll(j,i)
+C & ,wel_loc*gshieldc_loc_ll(j,i)
+C & ,wtube*gg_tube(j,i)
+
+C print *,gg_tube(1,0),"TU3"
#ifdef DEBUG
write (iout,*) "gloc before adding corr"
do i=1,4*nres
call MPI_Barrier(FG_COMM,IERR)
time_barrier_g=time_barrier_g+MPI_Wtime()-time00
time00=MPI_Wtime()
- call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres,
+ call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*nres+3,
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres,
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
c & (zmedi.lt.((-0.5d0)*boxzsize))) then
c go to 196
c endif
- xmedi=mod(xmedi,boxxsize)
+ xmedi=dmod(xmedi,boxxsize)
if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
+ ymedi=dmod(ymedi,boxysize)
if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
+ zmedi=dmod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
- zmedi2=mod(zmedi,boxzsize)
+ zmedi2=dmod(zmedi,boxzsize)
if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize
if ((zmedi2.gt.bordlipbot)
&.and.(zmedi2.lt.bordliptop)) then
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
+ xmedi=dmod(xmedi,boxxsize)
if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
+ ymedi=dmod(ymedi,boxysize)
if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
+ zmedi=dmod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
if ((zmedi.gt.bordlipbot)
&.and.(zmedi.lt.bordliptop)) then
enddo
enddo
if (isubchap.eq.1) then
+C print *,i,j
xj=xj_temp-xmedi
yj=yj_temp-ymedi
zj=zj_temp-zmedi
C Lipidic part for lipscale
gelc_long(3,j)=gelc_long(3,j)+
& ssgradlipj*eesij/2.0d0*lipscale**2
-
+C if ((ssgradlipj*eesij/2.0d0*lipscale**2).ne.0.0 )
+C & write(iout,*) "WTF",j
gelc_long(3,i)=gelc_long(3,i)+
& ssgradlipi*eesij/2.0d0*lipscale**2
+C if ((ssgradlipi*eesij/2.0d0*lipscale**2).ne.0.0 )
+C & write(iout,*) "WTF",i
+
*
* Loop over residues i+1 thru j-1.
*
#endif
cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'eelloc',i,j,eel_loc_ij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2f7.3)')
+ & 'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2)
c if (eel_loc_ij.ne.0)
c & write (iout,'(a4,2i4,8f9.5)')'chuj',
c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
&*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
C Cartesian derivatives
+!DIR$ UNROLL(0)
do l=1,3
c ghalf1=0.5d0*agg(l,1)
c ghalf2=0.5d0*agg(l,2)
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
dimension ggg(3)
+ integer xshift,yshift,zshift
evdw2=0.0D0
evdw2_14=0.0d0
c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
else
do j=1,nbi
diff=vbld(i+nres)-vbldsc0(j,iti)
+ if (energy_dec) write (iout,*)
+ & "estr sc",i,iti,vbld(i+nres),vbldsc0(j,iti),diff,
+ & AKSC(j,iti),AKSC(j,iti)*diff*diff
ud(j)=aksc(j,iti)*diff
u(j)=abond0(j,iti)+0.5d0*ud(j)*diff
enddo
c time11=dexp(-2*time)
c time12=1.0d0
etheta=0.0D0
-c write (*,'(a,i2)') 'EBEND ICG=',icg
+ write (*,'(a,i2)') 'EBEND ICG=',icg
do i=ithet_start,ithet_end
+ if (i.le.2) cycle
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
C Zero the energy function and its derivative at 0 or pi.
ichir22=isign(1,itype(i))
endif
- if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+ if (i.gt.3 ) then
+ if (itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
y(1)=0.0D0
y(2)=0.0D0
endif
+ else
+ y(1)=0.0D0
+ y(2)=0.0D0
+ endif
+
if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
+ if (i.le.2) cycle
c print *,i,itype(i-1),itype(i),itype(i-2)
+C if (itype(i-1).eq.ntyp1) cycle
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
C print *,i,theta(i)
sinkt(k)=dsin(k*theti2)
enddo
C print *,ethetai
- if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+ if (i.gt.3) then
+ if (itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
sinph1(k)=0.0d0
enddo
endif
+ else
+ phii=0.0d0
+ do k=1,nsingle
+ ityp1=ithetyp((itype(i-2)))
+ cosph1(k)=0.0d0
+ sinph1(k)=0.0d0
+ enddo
+ endif
+
if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
+ if (energy_dec) write(iout,'(a9,2i4,f8.3,3i4)') "esccor",i,j,
+ & esccor,intertyp,
+ & isccori, isccori1
c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
if (lprn)
eliptran=0.0
do i=ilip_start,ilip_end
C do i=1,1
- if (itype(i).eq.ntyp1) cycle
+ if ((itype(i).eq.ntyp1).or.(itype(i+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
include 'COMMON.SBRIDGE'
double precision tub_r,vectube(3),enetube(maxres*2)
Etube=0.0d0
- do i=1,2*nres
+ 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=1,nres
+ 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
- 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
+ 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 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
+ 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+
+ 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
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 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)
C lets ommit dummy atoms for now
C in UNRES uncomment the line below as GLY has no side-chain...
C .or.(iti.eq.10)
& ) cycle
- vectube(1)=c(1,i+nres)
- 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
-
+ 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 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
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
+ 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+
+ 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(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
enddo
enddo
- do i=1,2*nres
- Etube=Etube+enetube(i)
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)
enddo
C print *,"ETUBE", etube
return
include 'COMMON.SBRIDGE'
double precision tub_r,vectube(3),enetube(maxres*2)
Etube=0.0d0
- do i=1,2*nres
+ 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=1,nres
+ 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
- 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
+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)
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))
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+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)=enetube(i)+sstube*tubetranenepep
+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)=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)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6
+ 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
+ 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
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)
- do i=1,nres
+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)
C lets ommit dummy atoms for now
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
+C print *,positi,bordtubebot,buftubebot,bordtubetop
+
if ((positi.gt.bordtubebot)
& .and.(positi.lt.bordtubetop)) then
C the energy transfer exist
C lipbufthick is thickenes of lipid buffore
sstube=sscalelip(fracinbuf)
ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
- print *,ssgradtube, sstube,tubetranene(itype(i))
+C 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 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)
+ 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+
+ 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
&+ssgradtube*enetube(i+nres)/sstube
enddo
- do i=1,2*nres
- Etube=Etube+enetube(i)
+ 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
+
+
+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),
+ & enecavtube(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)=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
+C print *,"TUT",i,iti,rdiff,rdiff6,acavtubpep,denominator,
+C & enecavtube(i),faccav
+C print *,"licz=",
+C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+CX print *,"finene=",enetube(i+nres)+enecavtube(i)
+
+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.0d0
+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)=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
+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 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
+ 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)+enecavtube(i)
+ & +enecavtube(i+nres)
enddo
C print *,"ETUBE", etube
return
projects / unres.git / blobdiff