X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=443471d7601ecc37d4dc4ece97fd323f29b76a1f;hb=654bf369a7eb3567a67f0daf38f551b640b79369;hp=e6275dd781153a082a4f18611c6814cffe1f7cd8;hpb=6305a6478c25149aad01ec9e15e0c77814851ffa;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index e6275dd..443471d 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -141,7 +141,7 @@ cmc Sep-06: egb takes care of dynamic ss bonds too 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 @@ -232,7 +232,7 @@ C#undef DEBUG 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. @@ -254,11 +254,11 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" 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. @@ -281,9 +281,9 @@ c 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 @@ -300,13 +300,13 @@ C energy function 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 @@ -325,7 +325,7 @@ C energy function 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 @@ -334,7 +334,7 @@ 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 @@ -344,7 +344,7 @@ 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 @@ -3692,13 +3692,13 @@ c if ((zmedi.gt.((0.5d0)*boxzsize)).or. 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 @@ -3757,11 +3757,11 @@ c & .or. itype(i-1).eq.ntyp1 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 @@ -3965,6 +3965,7 @@ C lipbufthick is thickenes of lipid buffore enddo enddo if (isubchap.eq.1) then +C print *,i,j xj=xj_temp-xmedi yj=yj_temp-ymedi zj=zj_temp-zmedi @@ -4173,10 +4174,14 @@ C print *,"bafter", gelc_long(1,i), gelc_long(1,j) 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. * @@ -4647,8 +4652,8 @@ c & a33*gmuji2(4) #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) @@ -5120,6 +5125,7 @@ C Derivatives in gamma(i+1) &*((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) @@ -5809,6 +5815,7 @@ C 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' @@ -6329,6 +6336,9 @@ c 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 @@ -6387,8 +6397,9 @@ C 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. @@ -6407,7 +6418,8 @@ 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 @@ -6420,6 +6432,11 @@ C Zero the energy function and its derivative at 0 or pi. 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) @@ -6661,7 +6678,9 @@ C 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) @@ -6677,7 +6696,8 @@ 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 @@ -6698,6 +6718,15 @@ C propagation of chirality for glycine type 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) @@ -8219,6 +8248,9 @@ c 3 = SC...Ca...Ca...SCi 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) @@ -11565,7 +11597,8 @@ C--bordlipbot--buffore ends 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 @@ -12337,10 +12370,31 @@ 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) @@ -12359,14 +12413,61 @@ 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)) + 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 + & 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 @@ -12375,6 +12476,11 @@ C now direction of gg_tube vector 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" @@ -12401,7 +12507,8 @@ C THIS FRAGMENT MAKES TUBE FINITE 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 @@ -12411,7 +12518,7 @@ 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)) @@ -12531,17 +12638,17 @@ C now calculate distance from center of tube and direction vectors zmin=boxzsize do j=-1,1 - vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=dmod((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)=dmod((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)=dmod((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)) + 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 @@ -12589,7 +12696,30 @@ C now we calculate gradient & 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 @@ -12598,7 +12728,7 @@ C now direction of gg_tube vector enddo do i=itube_start,itube_end - enecavtube(i)=0.0 + 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 @@ -12610,17 +12740,17 @@ C .or.(iti.eq.10) ymin=boxysize zmin=boxzsize do j=-1,1 - vectube(1)=mod((c(1,i+nres)),boxxsize) + vectube(1)=dmod((c(1,i+nres)),boxxsize) vectube(1)=vectube(1)+boxxsize*j - vectube(2)=mod((c(2,i+nres)),boxysize) + vectube(2)=dmod((c(2,i+nres)),boxysize) vectube(2)=vectube(2)+boxysize*j - vectube(3)=mod((c(3,i+nres)),boxzsize) + vectube(3)=dmod((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)) + 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 @@ -12659,23 +12789,37 @@ C for vectorization reasons we will sumup at the end to avoid depenence of previ 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) go to 667 - denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6) - enecavtube(i)= - & acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff)+cavtub(iti)) + 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 - 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)) +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 - 667 continue +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 @@ -12690,7 +12834,8 @@ C if (acavtub(iti).eq.0.0) cycle do i=itube_start,itube_end - Etube=Etube+enetube(i)+enetube(i+nres) + Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) + & +enecavtube(i+nres) enddo C print *,"ETUBE", etube return