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
& 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
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
C if (j.eq.78)
C & write(iout,*) i,j,fac_shield(i),fac_shield(j)
- eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+ 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)
#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)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
#endif
C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
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-----------------------------------------------------------------------------
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
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
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 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
projects / unres.git / blobdiff