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=d44e78fefa5263e7ba8b02dfcee4806d5bf6c335;hb=5ba2a6014fbd859202ed28885b66520ac09bb77d;hp=c0670a7cce6bea6218e894826890bd6c1f60a925;hpb=688de817003e3e5f38e18eee73f13c25419c7f00;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 c0670a7..d44e78f 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -28,6 +28,12 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.TIME1' include 'COMMON.SPLITELE' include 'COMMON.SHIELD' + double precision fac_shieldbuf(maxres), + & grad_shield_locbuf(3,maxcontsshi,-1:maxres), + & grad_shield_sidebuf(3,maxcontsshi,-1:maxres), + & grad_shieldbuf(3,-1:maxres) + integer ishield_listbuf(maxres), + &shield_listbuf(maxcontsshi,maxres) #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -162,34 +168,50 @@ C#define DEBUG & grad_shield_side(1,1,i),grad_shield_loc(1,1,i) enddo #endif - call MPI_Allgatherv(fac_shield(ivec_start),ivec_count(fg_rank1), - & MPI_DOUBLE_PRECISION,fac_shield(1),ivec_count(0),ivec_displ(0), + call MPI_Allgatherv(fac_shield(ivec_start), + & ivec_count(fg_rank1), + & MPI_DOUBLE_PRECISION,fac_shieldbuf(1),ivec_count(0), + & ivec_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM,IERR) call MPI_Allgatherv(shield_list(1,ivec_start), & ivec_count(fg_rank1), - & MPI_I50,shield_list(1,1),ivec_count(0), + & MPI_I50,shield_listbuf(1,1),ivec_count(0), & ivec_displ(0), & MPI_I50,FG_COMM,IERR) call MPI_Allgatherv(ishield_list(ivec_start), & ivec_count(fg_rank1), - & MPI_INTEGER,ishield_list(1),ivec_count(0), + & MPI_INTEGER,ishield_listbuf(1),ivec_count(0), & ivec_displ(0), & MPI_INTEGER,FG_COMM,IERR) call MPI_Allgatherv(grad_shield(1,ivec_start), & ivec_count(fg_rank1), - & MPI_UYZ,grad_shield(1,1),ivec_count(0), + & MPI_UYZ,grad_shieldbuf(1,1),ivec_count(0), & ivec_displ(0), & MPI_UYZ,FG_COMM,IERR) call MPI_Allgatherv(grad_shield_side(1,1,ivec_start), & ivec_count(fg_rank1), - & MPI_SHI,grad_shield_side(1,1,1),ivec_count(0), + & MPI_SHI,grad_shield_sidebuf(1,1,1),ivec_count(0), & ivec_displ(0), & MPI_SHI,FG_COMM,IERR) call MPI_Allgatherv(grad_shield_loc(1,1,ivec_start), & ivec_count(fg_rank1), - & MPI_SHI,grad_shield_loc(1,1,1),ivec_count(0), + & MPI_SHI,grad_shield_locbuf(1,1,1),ivec_count(0), & ivec_displ(0), & MPI_SHI,FG_COMM,IERR) + do i=1,nres + fac_shield(i)=fac_shieldbuf(i) + ishield_list(i)=ishield_listbuf(i) + do j=1,3 + grad_shield(j,i)=grad_shieldbuf(j,i) + enddo !j + do j=1,ishield_list(i) + shield_list(j,i)=shield_listbuf(j,i) + do k=1,3 + grad_shield_loc(k,j,i)=grad_shield_locbuf(k,j,i) + grad_shield_side(k,j,i)=grad_shield_sidebuf(k,j,i) + enddo !k + enddo !j + enddo !i #ifdef DEBUG write(iout,*) "after reduce fac_shield reduce" do i=1,nres @@ -374,6 +396,8 @@ C print *,"just before call" call calctube(Etube) elseif (TUBElog.eq.2) then call calctube2(Etube) + elseif (TUBElog.eq.3) then + call calcnano(Etube) else Etube=0.0d0 endif @@ -649,7 +673,26 @@ c enddo 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 @@ -704,7 +747,7 @@ c call flush(iout) #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 @@ -880,7 +923,42 @@ C print *,gradafm(1,13),"AFM" 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 @@ -932,7 +1010,7 @@ c#undef DEBUG 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) @@ -3614,6 +3692,7 @@ 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) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) if (ymedi.lt.0) ymedi=ymedi+boxysize @@ -12073,20 +12152,38 @@ C simple Kihara potential 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) @@ -12106,12 +12203,12 @@ C calculte rdiffrence between r and r0 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 @@ -12123,7 +12220,10 @@ C now direction of gg_tube vector 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 @@ -12131,13 +12231,29 @@ 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) @@ -12149,6 +12265,7 @@ C now calculte the distance 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 @@ -12156,10 +12273,10 @@ 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 @@ -12167,8 +12284,8 @@ C now direction of gg_tube vector 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 @@ -12208,13 +12325,14 @@ C simple Kihara potential 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 @@ -12242,12 +12360,12 @@ C calculte rdiffrence between r and r0 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 @@ -12259,7 +12377,8 @@ C now direction of gg_tube vector 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 @@ -12337,11 +12456,11 @@ 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) & *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 @@ -12354,8 +12473,224 @@ C now direction of gg_tube vector &+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)=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)+ccavtub(iti)) + & /denominator + faccav=(acavtub(iti)*(1.0+bcavtub(iti)/2.0/sqrt(rdiff)) + & *denominator-acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff) + & +ccavtub(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