#ifdef MPI
include "mpif.h"
double precision weights_(n_ene)
+ integer IERR
+ integer status(MPI_STATUS_SIZE)
#endif
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
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
call set_shield_fac
else if (shield_mode.eq.2) then
call set_shield_fac2
+ if (nfgtasks.gt.1) then
+C#define DEBUG
+#ifdef DEBUG
+ write(iout,*) "befor reduce fac_shield reduce"
+ do i=1,nres
+ write(2,*) "fac",itype(i),fac_shield(i),grad_shield(1,i)
+ write(2,*) "list", shield_list(1,i),ishield_list(i),
+ & 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_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_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_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_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_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_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
+ write(2,*) "fac",itype(i),fac_shield(i),grad_shield(1,i)
+ write(2,*) "list", shield_list(1,i),ishield_list(i),
+ & grad_shield_side(1,1,i),grad_shield_loc(1,1,i)
+ enddo
+#endif
+C#undef DEBUG
+ endif
+#ifdef DEBUG
+ do i=1,nres
+ write(iout,*) fac_shield(i),ishield_list(i),i,grad_shield(1,i)
+ do j=1,ishield_list(i)
+ write(iout,*) "grad", grad_shield_side(1,j,i),
+ & grad_shield_loc(1,j,i)
+ enddo
+ enddo
+#endif
endif
c print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
C print *,"przed lipidami"
if (wliptran.gt.0) then
call Eliptransfer(eliptran)
+ else
+ eliptran=0.0d0
endif
C print *,"za lipidami"
if (AFMlog.gt.0) then
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
endif
- if (TUBElog.gt.0) then
+ if (TUBElog.eq.1) then
C print *,"just before call"
call calctube(Etube)
+ elseif (TUBElog.eq.2) then
+ call calctube2(Etube)
else
Etube=0.0d0
endif
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
include 'COMMON.SPLITELE'
+ include 'COMMON.SHIELD'
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
if (ymedi.lt.0) ymedi=ymedi+boxysize
zmedi=mod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ zmedi2=mod(zmedi,boxzsize)
+ if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize
+ if ((zmedi2.gt.bordlipbot)
+ &.and.(zmedi2.lt.bordliptop)) then
+C the energy transfer exist
+ if (zmedi2.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zmedi2-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zmedi2.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0d0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0d0
+ endif
num_conti=0
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
if (ymedi.lt.0) ymedi=ymedi+boxysize
zmedi=mod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
-
+ zmedi2=mod(zmedi,boxzsize)
+ if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize
+ if ((zmedi2.gt.bordlipbot)
+ &.and.(zmedi2.lt.bordliptop)) then
+C the energy transfer exist
+ if (zmedi2.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zmedi2-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zmedi2.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
num_conti=num_cont_hb(i)
c write(iout,*) "JESTEM W PETLI"
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (ymedi.lt.0) ymedi=ymedi+boxysize
zmedi=mod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ if ((zmedi.gt.bordlipbot)
+ &.and.(zmedi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zmedi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zmedi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zmedi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+C print *,sslipi,"TU?!"
C xmedi=xmedi+xshift*boxxsize
C ymedi=ymedi+yshift*boxysize
C zmedi=zmedi+zshift*boxzsize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
xj_safe=xj
yj_safe=yj
el2=el2*fac_shield(i)**2*fac_shield(j)**2
eesij=(el1+el2)
ees=ees+eesij
+C FOR NOW SHIELD IS NOT USED WITH LIPSCALE
+C & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
else
fac_shield(i)=1.0
fac_shield(j)=1.0
eesij=(el1+el2)
ees=ees+eesij
+ &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+C print *,"TUCC",(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
endif
evdw1=evdw1+evdwij*sss
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+C print *,sslipi,sslipj,lipscale**2,
+C & (sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
C
#ifdef SPLITELE
facvdw=-6*rrmij*(ev1+evdwij)*sss
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
facel=-3*rrmij*(el1+eesij)
+ &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
fac1=fac
erij(1)=xj*rmij
erij(2)=yj*rmij
C & +grad_shield(k,j)*eesij/fac_shield(j)
enddo
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
+
+ gelc_long(3,i)=gelc_long(3,i)+
+ & ssgradlipi*eesij/2.0d0*lipscale**2
*
* Loop over residues i+1 thru j-1.
cgrad enddo
if (sss.gt.0.0) then
ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
else
ggg(1)=0.0
ggg(2)=0.0
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
enddo
+C Lipidic part for scaling weight
+ gvdwpp(3,j)=gvdwpp(3,j)+
+ & sss*ssgradlipj*evdwij/2.0d0*lipscale**2
+ gvdwpp(3,i)=gvdwpp(3,i)+
+ & sss*ssgradlipi*evdwij/2.0d0*lipscale**2
+
*
* Loop over residues i+1 thru j-1.
*
#else
C MARYSIA
facvdw=(ev1+evdwij)*sss
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
facel=(el1+eesij)
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
cgrad enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
do k=1,3
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
enddo
+ gvdwpp(3,j)=gvdwpp(3,j)+
+ & sss*ssgradlipj*evdwij/2.0d0*lipscale**2
+ gvdwpp(3,i)=gvdwpp(3,i)+
+ & sss*ssgradlipi*evdwij/2.0d0*lipscale**2
+
#endif
*
* Angular part
do k=1,3
ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
& fac_shield(i)**2*fac_shield(j)**2
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
& +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
& *fac_shield(i)**2*fac_shield(j)**2
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
gelc(k,j)=gelc(k,j)
& +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
& *fac_shield(i)**2*fac_shield(j)**2
+ & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
endif
eel_loc_ij=eel_loc_ij
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Now derivative over eel_loc
if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
& (shield_mode.gt.0)) then
& +a32*gmuij1(3)
& +a33*gmuij1(4))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
c write(iout,*) "derivative over thatai"
c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
c & a33*gmuij1(4)
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& geel_loc_ij*wel_loc
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
c Derivative over j residue
geel_loc_ji=a22*gmuji1(1)
gloc(nphi+j,icg)=gloc(nphi+j,icg)+
& geel_loc_ji*wel_loc
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
geel_loc_ji=
& +a22*gmuji2(1)
gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
& geel_loc_ji*wel_loc
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
#endif
cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
& (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
& +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gel_loc_loc(j-1)=gel_loc_loc(j-1)+
& (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
& +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
do l=1,3
ggg(l)=(agg(l,1)*muij(1)+
& agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
cgrad ghalf=0.5d0*ggg(l)
cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
enddo
+ gel_loc_long(3,j)=gel_loc_long(3,j)+
+ & ssgradlipj*eel_loc_ij/2.0d0*lipscale/
+ & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+ gel_loc_long(3,i)=gel_loc_long(3,i)+
+ & ssgradlipi*eel_loc_ij/2.0d0*lipscale/
+ & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
cgrad do k=i+1,j2
cgrad do l=1,3
cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
& aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
& aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
& aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
& aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
enddo
ENDIF
& dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
& num_conti,j1,j2
j=i+2
-c write (iout,*) "eturn3",i,j,j1,j2
+C xj=(c(1,j)+c(1,j+1))/2.0d0
+C yj=(c(2,j)+c(2,j+1))/2.0d0
+ zj=(c(3,j)+c(3,j+1))/2.0d0
+C xj=mod(xj,boxxsize)
+C if (xj.lt.0) xj=xj+boxxsize
+C yj=mod(yj,boxysize)
+C if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+C sslipj=0.0
+C ssgradlipj=0.0d0
+
+C write (iout,*) "eturn3",i,j,j1,j2
a_temp(1,1)=a22
a_temp(1,2)=a23
a_temp(2,1)=a32
call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
if (shield_mode.eq.0) then
- fac_shield(i)=1.0
- fac_shield(j)=1.0
+ fac_shield(i)=1.0d0
+ fac_shield(j)=1.0d0
C else
C fac_shield(i)=0.4
C fac_shield(j)=0.6
endif
- eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+C if (j.eq.78)
+C & write(iout,*) i,j,fac_shield(i),fac_shield(j)
+ eello_turn3=eello_turn3+
+C & 1.0*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+ &0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
- eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+ eello_t3=
+ &0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
-C#ifdef NEWCORR
+#ifdef NEWCORR
C Derivatives in theta
gloc(nphi+i,icg)=gloc(nphi+i,icg)
& +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
& +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
& *fac_shield(i)*fac_shield(j)
-C#endif
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+#endif
C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
C Derivatives in shield mode
call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Derivatives in gamma(i+1)
call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
call transpose2(auxmat2(1,1),auxmat3(1,1))
gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Cartesian derivatives
do l=1,3
c ghalf1=0.5d0*agg(l,1)
gcorr3_turn(l,i)=gcorr3_turn(l,i)
& +0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
a_temp(1,1)=aggi1(l,1)!+agg(l,1)
a_temp(1,2)=aggi1(l,2)!+agg(l,2)
gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
a_temp(1,1)=aggj(l,1)!+ghalf1
a_temp(1,2)=aggj(l,2)!+ghalf2
a_temp(2,1)=aggj(l,3)!+ghalf3
gcorr3_turn(l,j)=gcorr3_turn(l,j)
& +0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
a_temp(2,1)=aggj1(l,3)
gcorr3_turn(l,j1)=gcorr3_turn(l,j1)
& +0.5d0*(pizda(1,1)+pizda(2,2))
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
enddo
+ gshieldc_t3(3,i)=gshieldc_t3(3,i)+
+ & ssgradlipi*eello_t3/4.0d0*lipscale
+ gshieldc_t3(3,j)=gshieldc_t3(3,j)+
+ & ssgradlipj*eello_t3/4.0d0*lipscale
+ gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+
+ & ssgradlipi*eello_t3/4.0d0*lipscale
+ gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+
+ & ssgradlipj*eello_t3/4.0d0*lipscale
+
+C print *,ssgradlipi,ssgradlipj,eello_t3,sslipi,sslipj
return
end
C-------------------------------------------------------------------------------
cd call checkint_turn4(i,a_temp,eello_turn4_num)
c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
c write(iout,*)"WCHODZE W PROGRAM"
+ zj=(c(3,j)+c(3,j+1))/2.0d0
+C xj=mod(xj,boxxsize)
+C if (xj.lt.0) xj=xj+boxxsize
+C yj=mod(yj,boxysize)
+C if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+C if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+
a_temp(1,1)=a22
a_temp(1,2)=a23
a_temp(2,1)=a32
endif
eello_turn4=eello_turn4-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
eello_t4=-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
gloc(nphi+i,icg)=gloc(nphi+i,icg)
& -(gs13+gsE13+gsEE1)*wturn4
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
& -(gs23+gs21+gsEE2)*wturn4
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
& -(gs32+gsE31+gsEE3)*wturn4
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
c & gs2
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Derivatives in gamma(i+1)
call transpose2(EUgder(1,1,i+2),e2tder(1,1))
call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Derivatives in gamma(i+2)
call transpose2(EUgder(1,1,i+3),e3tder(1,1))
call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Cartesian derivatives
C Derivatives of this turn contributions in DC(i+2)
if (j.lt.nres-1) then
ggg(l)=-(s1+s2+s3)
gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
enddo
endif
C Remaining derivatives of this turn contribution
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
a_temp(1,1)=aggi1(l,1)
a_temp(1,2)=aggi1(l,2)
a_temp(2,1)=aggi1(l,3)
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
a_temp(1,1)=aggj(l,1)
a_temp(1,2)=aggj(l,2)
a_temp(2,1)=aggj(l,3)
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
a_temp(2,1)=aggj1(l,3)
c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3
gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
& *fac_shield(i)*fac_shield(j)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
enddo
+ gshieldc_t4(3,i)=gshieldc_t4(3,i)+
+ & ssgradlipi*eello_t4/4.0d0*lipscale
+ gshieldc_t4(3,j)=gshieldc_t4(3,j)+
+ & ssgradlipj*eello_t4/4.0d0*lipscale
+ gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+
+ & ssgradlipi*eello_t4/4.0d0*lipscale
+ gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+
+ & ssgradlipj*eello_t4/4.0d0*lipscale
return
end
C-----------------------------------------------------------------------------
include 'COMMON.IOUNITS'
include 'COMMON.SHIELD'
include 'COMMON.INTERACT'
+ include 'COMMON.LOCAL'
+
C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
double precision div77_81/0.974996043d0/,
&div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
-
+
C the vector between center of side_chain and peptide group
double precision pep_side(3),long,side_calf(3),
&pept_group(3),costhet_grad(3),cosphi_grad_long(3),
&cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+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
C the line belowe needs to be changed for FGPROC>1
- do i=1,nres-1
- if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ do i=ivec_start,ivec_end
+C do i=1,nres-1
+C if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
ishield_list(i)=0
+ if ((itype(i).eq.ntyp1).and.itype(i+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
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
VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
enddo
fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
-C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+C write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i)
enddo
return
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)=(c(1,i)+c(1,i+1))/2.0d0-tubecenter(1)
- vectube(2)=(c(2,i)+c(2,i+1))/2.0d0-tubecenter(2)
+ 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 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
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)-tubecenter(1)
- vectube(2)=c(2,i+nres)-tubecenter(2)
+ 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
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
C 5) add COMMONs
C 6) add to zerograd
+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)
+ 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=1,2*nres
+ enetube(i)=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
+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
+ 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)
+ do i=1,nres
+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...
+ & .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)
+ 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
+ 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 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+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 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+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+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=1,2*nres
+ Etube=Etube+enetube(i)
+ 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 / blobdiff