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
C
C Calculate the SC local energy.
C
+C print *,"TU DOCHODZE?"
call esc(escloc)
c print *,"Processor",myrank," computed USC"
C
else
esccor=0.0d0
endif
+C print *,"PRZED MULIt"
c print *,"Processor",myrank," computed Usccorr"
C
C 12/1/95 Multi-body terms
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
+C print *,"przed lipidami"
+ if (wliptran.gt.0) then
+ call Eliptransfer(eliptran)
+ endif
+C print *,"za lipidami"
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
energia(17)=estr
energia(20)=Uconst+Uconst_back
energia(21)=esccor
+ energia(22)=eliptran
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)
+ eliptran=energia(22)
#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'
- double precision gradbufc(3,maxres),gradbufx(3,maxres),
- & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
#endif
+ double precision gradbufc(3,0:maxres),gradbufx(3,0:maxres),
+ & glocbuf(4*maxres),gradbufc_sum(3,0:maxres),gloc_scbuf(3,0:maxres)
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
+ & +wliptran*gliptranc(j,i)
+
enddo
enddo
#else
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
+ & +wliptran*gliptranc(j,i)
enddo
enddo
#endif
& 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
do i=1,4*nres
glocbuf(i)=gloc(i,icg)
enddo
-#define DEBUG
+c#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc_sc before reduce"
do i=1,nres
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
do i=1,nres
do j=1,3
gloc_scbuf(j,i)=gloc_sc(j,i,icg)
call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
time_reduce=time_reduce+MPI_Wtime()-time00
-#define DEBUG
+c#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc_sc after reduce"
do i=1,nres
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
#ifdef DEBUG
write (iout,*) "gloc after reduce"
do i=1,4*nres
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ eliptran=energia(22)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
& edihcnstr,ebr*nss,
- & Uconst,etot
+ & Uconst,eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ebr*nss,Uconst,etot
+ & ebr*nss,Uconst,eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
+C have you changed here?
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
evdwij=e1+e2
rij=1.0D0/r_inv_ij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
+C have you changed here?
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
evdwij=e_augm+e1+e2
call sc_angular
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
+C have you changed here?
fac=(rrij*sigsq)**expon2
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
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
yi=c(2,nres+i)
zi=c(3,nres+i)
C Return atom into box, boxxsize is size of box in x dimension
- 134 continue
- if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
- if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+c 134 continue
+c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
+c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
C Condition for being inside the proper box
- if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 134
- endif
- 135 continue
- if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
- if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
+c & (xi.lt.((xshift-0.5d0)*boxxsize))) then
+c go to 134
+c endif
+c 135 continue
+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
- if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
- & (yi.lt.((yshift-0.5d0)*boxysize))) then
- go to 135
- endif
- 136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
+c if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
+c & (yi.lt.((yshift-0.5d0)*boxysize))) then
+c go to 135
+c endif
+c 136 continue
+c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
C Condition for being inside the proper box
- if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 136
- endif
+c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
+c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
+c go to 136
+c endif
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ 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
yj=c(2,nres+j)
zj=c(3,nres+j)
C Return atom J into box the original box
- 137 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+c 137 continue
+c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
C Condition for being inside the proper box
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 137
- endif
- 138 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+c if ((xj.gt.((0.5d0)*boxxsize)).or.
+c & (xj.lt.((-0.5d0)*boxxsize))) then
+c go to 137
+c endif
+c 138 continue
+c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
C Condition for being inside the proper box
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 138
- endif
- 139 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+c if ((yj.gt.((0.5d0)*boxysize)).or.
+c & (yj.lt.((-0.5d0)*boxysize))) then
+c go to 138
+c endif
+c 139 continue
+c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
C Condition for being inside the proper box
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 139
- endif
-
+c if ((zj.gt.((0.5d0)*boxzsize)).or.
+c & (zj.lt.((-0.5d0)*boxzsize))) then
+c go to 139
+c endif
+ 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-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)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
+C here to start with
+C if (c(i,3).gt.
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
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)
C
C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
do i=iturn3_start,iturn3_end
+ if (i.le.1) cycle
+C write(iout,*) "tu jest i",i
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+2).eq.ntyp1
& .or. itype(i+3).eq.ntyp1
-c & .or. itype(i-1).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
& .or. itype(i+4).eq.ntyp1
& ) cycle
dxi=dc(1,i)
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
- 184 continue
- if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
- if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
- if ((xmedi.gt.((0.5d0)*boxxsize)).or.
- & (xmedi.lt.((-0.5d0)*boxxsize))) then
- go to 184
- endif
- 185 continue
- if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
- if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-C Condition for being inside the proper box
- if ((ymedi.gt.((0.5d0)*boxysize)).or.
- & (ymedi.lt.((-0.5d0)*boxysize))) then
- go to 185
- endif
- 186 continue
- if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
- if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-C Condition for being inside the proper box
- if ((zmedi.gt.((0.5d0)*boxzsize)).or.
- & (zmedi.lt.((-0.5d0)*boxzsize))) then
- go to 186
- endif
+ 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
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
+ if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+3).eq.ntyp1
& .or. itype(i+4).eq.ntyp1
& .or. itype(i+5).eq.ntyp1
+ & .or. itype(i).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
& ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
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
- 194 continue
- if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
- if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
+c 194 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
- if ((xmedi.gt.((0.5d0)*boxxsize)).or.
- & (xmedi.lt.((-0.5d0)*boxxsize))) then
- go to 194
- endif
- 195 continue
- if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
- if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+c if ((xmedi.gt.((0.5d0)*boxxsize)).or.
+c & (xmedi.lt.((-0.5d0)*boxxsize))) then
+c go to 194
+c endif
+c 195 continue
+c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
+c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
C Condition for being inside the proper box
- if ((ymedi.gt.((0.5d0)*boxysize)).or.
- & (ymedi.lt.((-0.5d0)*boxysize))) then
- go to 195
- endif
- 196 continue
- if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
- if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
+c if ((ymedi.gt.((0.5d0)*boxysize)).or.
+c & (ymedi.lt.((-0.5d0)*boxysize))) then
+c go to 195
+c endif
+c 196 continue
+c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
C Condition for being inside the proper box
- if ((zmedi.gt.((0.5d0)*boxzsize)).or.
- & (zmedi.lt.((-0.5d0)*boxzsize))) then
- go to 196
- endif
+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
+ zmedi=mod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=num_cont_hb(i)
call eelecij(i,i+3,ees,evdw1,eel_loc)
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
do i=iatel_s,iatel_e
+ if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+2).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
& ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
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
- 164 continue
- if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
- if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
- if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 164
- endif
- 165 continue
- if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
- if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+ 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
+C xmedi=xmedi+xshift*boxxsize
+C ymedi=ymedi+yshift*boxysize
+C zmedi=zmedi+zshift*boxzsize
+
+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
C Condition for being inside the proper box
- if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
- & (ymedi.lt.((yshift-0.5d0)*boxysize))) then
- go to 165
- endif
- 166 continue
- if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
- if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
+c if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
+c & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
+c go to 164
+c endif
+c 165 continue
+c if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
+c if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
C Condition for being inside the proper box
- if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 166
- endif
+c if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
+c & (ymedi.lt.((yshift-0.5d0)*boxysize))) then
+c go to 165
+c endif
+c 166 continue
+c if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+c if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
+cC Condition for being inside the proper box
+c if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
+c & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
+c go to 166
+c endif
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
-c write (iout,*) i,j,itype(i),itype(j)
+C write (iout,*) i,j
+ if (j.le.1) cycle
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
& .or.itype(j+2).eq.ntyp1
+ & .or.itype(j-1).eq.ntyp1
&) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
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
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
+ 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
- 174 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+c 174 continue
+c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
C Condition for being inside the proper box
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 174
- endif
- 175 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+c if ((xj.gt.((0.5d0)*boxxsize)).or.
+c & (xj.lt.((-0.5d0)*boxxsize))) then
+c go to 174
+c endif
+c 175 continue
+c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
C Condition for being inside the proper box
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 175
- endif
- 176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+c if ((yj.gt.((0.5d0)*boxysize)).or.
+c & (yj.lt.((-0.5d0)*boxysize))) then
+c go to 175
+c endif
+c 176 continue
+c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
C Condition for being inside the proper box
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 176
- endif
+c if ((zj.gt.((0.5d0)*boxzsize)).or.
+c & (zj.lt.((-0.5d0)*boxzsize))) then
+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)
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
C Return atom into box, boxxsize is size of box in x dimension
- 134 continue
- if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
- if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
-C Condition for being inside the proper box
- if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 134
- endif
- 135 continue
- if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
- if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+c 134 continue
+c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
+c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
C Condition for being inside the proper box
- if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
- & (yi.lt.((yshift-0.5d0)*boxysize))) then
- go to 135
- endif
- 136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
+c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
+c & (xi.lt.((xshift-0.5d0)*boxxsize))) then
+c go to 134
+c endif
+c 135 continue
+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
- if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 136
- endif
+c if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
+c & (yi.lt.((yshift-0.5d0)*boxysize))) then
+c go to 135
+c c endif
+c 136 continue
+c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
+cC Condition for being inside the proper box
+c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
+c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
+c go to 136
+c endif
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ 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)
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- 174 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 174
- endif
- 175 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+c 174 continue
+c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
C Condition for being inside the proper box
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 175
- endif
- 176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+c if ((xj.gt.((0.5d0)*boxxsize)).or.
+c & (xj.lt.((-0.5d0)*boxxsize))) then
+c go to 174
+c endif
+c 175 continue
+c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+cC Condition for being inside the proper box
+c if ((yj.gt.((0.5d0)*boxysize)).or.
+c & (yj.lt.((-0.5d0)*boxysize))) then
+c go to 175
+c endif
+c 176 continue
+c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
C Condition for being inside the proper box
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 176
- endif
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+c if ((zj.gt.((0.5d0)*boxzsize)).or.
+c & (zj.lt.((-0.5d0)*boxzsize))) then
+c go to 176
+ 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-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
+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)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ 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
- 134 continue
- if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
- if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+c 134 continue
+c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
+c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
C Condition for being inside the proper box
- if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 134
- endif
- 135 continue
- if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
- if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
+c & (xi.lt.((xshift-0.5d0)*boxxsize))) then
+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
- if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
- & (yi.lt.((yshift-0.5d0)*boxysize))) then
- go to 135
- endif
- 136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
+c if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
+c & (yi.lt.((yshift-0.5d0)*boxysize))) then
+c go to 135
+c endif
+c 136 continue
+c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
C Condition for being inside the proper box
- if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 136
- endif
+c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
+c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
+c go to 136
+c endif
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- 174 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 174
- endif
- 175 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+ 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
+c 174 continue
+c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
C Condition for being inside the proper box
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 175
- endif
- 176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+c if ((xj.gt.((0.5d0)*boxxsize)).or.
+c & (xj.lt.((-0.5d0)*boxxsize))) then
+c go to 174
+c endif
+c 175 continue
+c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+cC Condition for being inside the proper box
+c if ((yj.gt.((0.5d0)*boxysize)).or.
+c & (yj.lt.((-0.5d0)*boxysize))) then
+c go to 175
+c endif
+c 176 continue
+c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
C Condition for being inside the proper box
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 176
- endif
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+c if ((zj.gt.((0.5d0)*boxzsize)).or.
+c & (zj.lt.((-0.5d0)*boxzsize))) then
+c go to 176
+c endif
+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
+C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C & 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
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)
+ 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'
+C print *,"wchodze"
+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=ilip_start,ilip_end
+ if (itype(i).eq.ntyp1) cycle
+
+ positi=(mod((c(3,i)+c(3,i+1)),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+C print *,i
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+ if ((positi.gt.bordlipbot)
+ &.and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+C print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+ print *, "doing sscalefor top part"
+ else
+ eliptran=eliptran+pepliptran
+ print *,"I am in true lipid"
+ endif
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ endif
+ enddo
+C print *, "nic nie bylo w lipidzie?"
+C now multiply all by the peptide group transfer factor
+C eliptran=eliptran*pepliptran
+C now the same for side chains
+ do i=ilip_start,ilip_end
+ if (itype(i).eq.ntyp1) cycle
+ 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)
+ if ((positi.gt.bordlipbot)
+ & .and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))/2.0d0
+ gliptranc(3,i-1)=
+ &+ssgradlip*liptranene(itype(i))
+ print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-
+ &((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))/2.0d0
+ gliptranc(3,i-1)=
+ &+ssgradlip*liptranene(itype(i))
+ print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i))
+ print *,"I am in true lipid"
+ endif
+ endif ! if in lipid or buffor
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ enddo
+ return
+ end
projects / unres.git / blobdiff