C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+cmc
+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"
#ifdef TIMING
time01=MPI_Wtime()
eello_turn4=0.0d0
endif
else
-c write (iout,*) "Soft-spheer ELEC potential"
- call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
- & eello_turn4)
+ 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
Uconst=0.0d0
Uconst_back=0.0d0
endif
+C 01/27/2015 added by adasko
+C the energy component below is energy transfer into lipid environment
+C based on partition function
+ if (wliptran.gt.0) then
+ call Eliptransfer
+ endif
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
energia(17)=estr
energia(20)=Uconst+Uconst_back
energia(21)=esccor
+ energia(22)=eliptrans
c Here are the energies showed per procesor if the are more processors
c per molecule then we sum it up in sum_energy subroutine
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
+ if (dyn_ss) call dyn_set_nss
c print *," Processor",myrank," left SUM_ENERGY"
#ifdef TIMING
time_sumene=time_sumene+MPI_Wtime()-time00
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ energia(22)=eliptrans
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
#endif
energia(0)=etot
c detecting NaNQ
#endif
#ifdef MPI
include 'mpif.h'
+#endif
double precision gradbufc(3,maxres),gradbufx(3,maxres),
& glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
-#endif
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
& wturn6*gcorr6_turn(j,i)+
& wsccor*gsccorc(j,i)
& +wscloc*gscloc(j,i)
+ & +wliptran*gliptranc(j,i)
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& wturn6*gcorr6_turn(j,i)+
& wsccor*gsccorc(j,i)
& +wscloc*gscloc(j,i)
+ & +wliptran*gliptranc(j,i)
#endif
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
enddo
#ifdef DEBUG
include 'COMMON.CALC'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
logical lprn
integer xshift,yshift,zshift
evdw=0.0D0
ind=0
C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0
C we have the original box)
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
do i=iatsc_s,iatsc_e
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
if (yi.lt.0) yi=yi+boxysize
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
+C xi=xi+xshift*boxxsize
+C yi=yi+yshift*boxysize
+C zi=zi+zshift*boxzsize
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+C xj=xj-xi
+C yj=yj-yi
+C zj=zj-zi
c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
c write (iout,*) "j",j," dc_norm",
c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
gg(3)=zj*fac
C Calculate angular part of the gradient.
call sc_grad
- endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
- enddo ! zshift
- enddo ! yshift
- enddo ! xshift
+C enddo ! zshift
+C enddo ! yshift
+C enddo ! xshift
c write (iout,*) "Number of loop steps in EGB:",ind
cccc energy_dec=.false.
return
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
dimension ggg(3)
-cd write(iout,*) 'In EELEC_soft_sphere'
+C write(iout,*) 'In EELEC_soft_sphere'
ees=0.0D0
evdw1=0.0D0
eel_loc=0.0d0
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)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=mod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=mod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
dxj=dc(1,j)
dyj=dc(2,j)
dzj=dc(3,j)
- xj=c(1,j)+0.5D0*dxj-xmedi
- yj=c(2,j)+0.5D0*dyj-ymedi
- zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ isubchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(sqrt(rij))
+ sssgrad=sscagrad(sqrt(rij))
if (rij.lt.r0ijsq) then
evdw1ij=0.25d0*(rij-r0ijsq)**2
fac=rij-r0ijsq
evdw1ij=0.0d0
fac=0.0d0
endif
- evdw1=evdw1+evdw1ij
+ evdw1=evdw1+evdw1ij*sss
C
C Calculate contributions to the Cartesian gradient.
C
- ggg(1)=fac*xj
- ggg(2)=fac*yj
- ggg(3)=fac*zj
+ ggg(1)=fac*xj*sssgrad
+ ggg(2)=fac*yj*sssgrad
+ ggg(3)=fac*zj*sssgrad
do k=1,3
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
-C Return atom into box, boxxsize is size of box in x dimension
-c 184 continue
-c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
-c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
-c if ((xmedi.gt.((0.5d0)*boxxsize)).or.
-c & (xmedi.lt.((-0.5d0)*boxxsize))) then
-c go to 184
-c endif
-c 185 continue
-c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
-c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-cC Condition for being inside the proper box
-c if ((ymedi.gt.((0.5d0)*boxysize)).or.
-c & (ymedi.lt.((-0.5d0)*boxysize))) then
-c go to 185
-c endif
-c 186 continue
-c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
-c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-cC Condition for being inside the proper box
-c if ((zmedi.gt.((0.5d0)*boxzsize)).or.
-c & (zmedi.lt.((-0.5d0)*boxzsize))) then
-c go to 186
-c endif
xmedi=mod(xmedi,boxxsize)
if (xmedi.lt.0) xmedi=xmedi+boxxsize
ymedi=mod(ymedi,boxysize)
num_cont_hb(i)=num_conti
enddo ! i
C Loop over all neighbouring boxes
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
if (ymedi.lt.0) ymedi=ymedi+boxysize
zmedi=mod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
+C xmedi=xmedi+xshift*boxxsize
+C ymedi=ymedi+yshift*boxysize
+C zmedi=zmedi+zshift*boxzsize
-C Return atom into box, boxxsize is size of box in x dimension
+C Return tom into box, boxxsize is size of box in x dimension
c 164 continue
c if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
c if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
enddo ! j
num_cont_hb(i)=num_conti
enddo ! i
- enddo ! zshift
- enddo ! yshift
- enddo ! xshift
+C enddo ! zshift
+C enddo ! yshift
+C enddo ! xshift
c write (iout,*) "Number of loop steps in EELEC:",ind
cd do i=1,nres
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"
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ isubchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
c 174 continue
c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
c go to 176
c endif
C endif !endPBC condintion
- xj=xj-xmedi
- yj=yj-ymedi
- zj=zj-zmedi
+C xj=xj-xmedi
+C yj=yj-ymedi
+C zj=zj-zmedi
rij=xj*xj+yj*yj+zj*zj
sss=sscale(sqrt(rij))
r0_scp=4.5d0
cd print '(a)','Enter ESCP'
cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
do i=iatscp_s,iatscp_e
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
if (yi.lt.0) yi=yi+boxysize
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
-
+C xi=xi+xshift*boxxsize
+C yi=yi+yshift*boxysize
+C zi=zi+zshift*boxzsize
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
c c endif
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+C xj=xj-xi
+C yj=yj-yi
+C zj=zj-zi
rij=xj*xj+yj*yj+zj*zj
r0ij=r0_scp
enddo ! iint
enddo ! i
- enddo !zshift
- enddo !yshift
- enddo !xshift
+C enddo !zshift
+C enddo !yshift
+C enddo !xshift
return
end
C-----------------------------------------------------------------------------
dimension ggg(3)
evdw2=0.0D0
evdw2_14=0.0d0
+c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
cd print '(a)','Enter ESCP'
cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
do i=iatscp_s,iatscp_e
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
if (yi.lt.0) yi=yi+boxysize
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
-
+c xi=xi+xshift*boxxsize
+c yi=yi+yshift*boxysize
+c zi=zi+zshift*boxzsize
+c print *,xi,yi,zi,'polozenie i'
C Return atom into box, boxxsize is size of box in x dimension
c 134 continue
c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
c go to 134
c endif
c 135 continue
+c print *,xi,boxxsize,"pierwszy"
+
c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
C Condition for being inside the proper box
c & (zj.lt.((-0.5d0)*boxzsize))) then
c go to 176
c endif
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+c print *,xj,yj,zj,'polozenie j'
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+c print *,rrij
sss=sscale(1.0d0/(dsqrt(rrij)))
+c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz'
+c if (sss.eq.0) print *,'czasem jest OK'
+ if (sss.le.0.0d0) cycle
sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
- if (sss.gt.0.0d0) then
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- endif !endif for sscale cutoff
+c endif !endif for sscale cutoff
enddo ! j
enddo ! iint
enddo ! i
- enddo !zshift
- enddo !yshift
- enddo !xshift
+c enddo !zshift
+c enddo !yshift
+c enddo !xshift
do i=1,nct
do j=1,3
gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
iii=ii
jjj=jj
endif
-cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
+c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
& iabs(itype(jjj)).eq.1) then
+cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
+ if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
+ if (ii.gt.nres
+ & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+>>>>>>> f5379d3246c4bd95e946c4d35d4a1c13e329c4cb
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+ endif
cd write (iout,*) "eij",eij
else
C Calculate the distance between the two points and its difference from the
C target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ dd=dist(ii,jj)
+ rdis=dd-dhpb(i)
C Get the force constant corresponding to this distance.
- waga=forcon(i)
+ waga=forcon(i)
C Calculate the contribution to energy.
- ehpb=ehpb+waga*rdis*rdis
+ ehpb=ehpb+waga*rdis*rdis
C
C Evaluate gradient.
C
- fac=waga*rdis/dd
+ fac=waga*rdis/dd
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
- do j=1,3
- ggg(j)=fac*(c(j,jj)-c(j,ii))
- enddo
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
C If this is a SC-SC distance, we need to calculate the contributions to the
C Cartesian gradient in the SC vectors (ghpbx).
- if (iii.lt.ii) then
+ if (iii.lt.ii) then
do j=1,3
ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
- endif
+ endif
cgrad do j=iii,jjj-1
cgrad do k=1,3
cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k)
cgrad enddo
cgrad enddo
- do k=1,3
- ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
- ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
- enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
endif
enddo
ehpb=0.5D0*ehpb
C NO vbldp0 is the equlibrium lenght of spring for peptide group
diff = vbld(i)-vbldp0
endif
- if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+ if (energy_dec) write (iout,'(a7,i5,4f7.3)')
& "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
do j=1,3
nbi=nbondterm(iti)
if (nbi.eq.1) then
diff=vbld(i+nres)-vbldsc0(1,iti)
- if (energy_dec) write (iout,*)
+ if (energy_dec) write (iout,*)
& "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
& AKSC(1,iti),AKSC(1,iti)*diff*diff
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
C o o C
C \ /l\ /j\ / C
C \ / \ / \ / C
-C o| o | | o |o C
+C o| o | | o |o C
C \ j|/k\| \ |/k\|l C
C \ / \ \ / \ C
C o o C
-C i i C
-C C
+C i i C
+C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
C AL 7/4/01 s1 would occur in the sixth-order moment,
double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C C
+C C
C Parallel Antiparallel C
C C
-C o o C
+C o o C
C /l\ / \ /j\ C
C / \ / \ / \ C
C /| o |o o| o |\ C
& auxvec1(2),auxmat1(2,2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C C
+C C
C Parallel Antiparallel C
C C
C o o C
C / \ / \ / \ C
C /| o |o o| o |\ C
C \ j|/k\| \ |/k\|l C
-C \ / \ \ / \ C
+C \ / \ \ / \ C
C o \ o \ C
C i i C
-C C
+C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C 4/7/01 AL Component s1 was removed, because it pertains to the respective
return
end
+CCC----------------------------------------------
+ subroutine Eliptransfer(eliptran)
+ 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'
+C structure of box:
+C water
+C--bordliptop-- buffore starts
+C--bufliptop--- here true lipid starts
+C lipid
+C--buflipbot--- lipid ends buffore starts
+C--bordlipbot--buffore ends
+ eliptran=0.0
+ do i=1,nres
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+ if ((mod(c(3,i),boxzsize).gt.bordlipbot)
+ &.and.(mod(c(3,i),boxzsize).lt.bordliptop)) then
+C the energy transfer exist
+ if (mod(c(3,i),boxzsize).lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((mod(c(3,i),boxzsize)-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ ssslip=sscale(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+C print *,"doing sccale for lower part"
+ elseif (mod(c(3,i),boxzsize).gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-mod(c(3,i),boxzsize))/lipbufthick)
+ ssslip=sscale(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+ print *, "doing sscalefor top part"
+ else
+ eliptran=eliptran+1.0d0
+ print *,"I am in true lipid"
+ endif
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ endif
+ enddo
+C now multiply all by the peptide group transfer factor
+ eliptran=eliptran*pepliptran
+C now the same for side chains
+ do i=1,nres
+c for each residue check if it is in lipid or lipid water border area
+ if ((mod(c(3,i+nres),boxzsize).gt.bordlipbot)
+ & .and.(mod(c(3,i+nres),boxzsize).lt.bordliptop)) then
+C the energy transfer exist
+ if (mod(c(3,i+nres),boxzsize).lt.buflipbot) then
+ fracinbuf=1.0d0-
+ & ((mod(c(3,i+nres),boxzsize)-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ ssslip=sscale(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i))
+ print *,"doing sccale for lower part"
+ elseif (mod(c(3,i+nres),boxzsize).gt.bufliptop) then
+ fracinbuf=1.0d0-
+ &((bordliptop-mod(c(3,i+nres),boxzsize))/lipbufthick)
+ ssslip=sscale(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i))
+ print *, "doing sscalefor top part"
+ else
+ eliptran=eliptran+liptranene(itype(i))
+ print *,"I am in true lipid"
+ endif
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ enddo
projects / unres.git / blobdiff