goto 106
C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
105 call egbv(evdw,evdw_t)
+C write(iout,*) 'po elektostatyce'
C
C Calculate electrostatic (H-bonding) energy of the main chain.
C
- 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-C
+ 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+C write(iout,*) 'po eelec'
+
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
C
c
c Calculate the bond-stretching energy
c
+
call ebond(estr)
-c write (iout,*) "estr",estr
+C write (iout,*) "estr",estr
C
C Calculate the disulfide-bridge and other energy and the contributions
C from other distance constraints.
C Calculate the virtual-bond-angle energy.
C
call ebend(ebe)
-cd print *,'Bend energy finished.'
+C print *,'Bend energy finished.'
C
C Calculate the SC local energy.
C
call esc(escloc)
-cd print *,'SCLOC energy finished.'
+C print *,'SCLOC energy finished.'
C
C Calculate the virtual-bond torsional energy.
C
C 21/5/07 Calculate local sicdechain correlation energy
C
call eback_sc_corr(esccor)
+
+ if (wliptran.gt.0) then
+ call Eliptransfer(eliptran)
+ endif
+
C
C 12/1/95 Multi-body terms
C
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
- & +wbond*estr+wsccor*fact(1)*esccor
+ & +wbond*estr+wsccor*fact(1)*esccor+wliptran*eliptran
#else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
- & +wbond*estr+wsccor*fact(1)*esccor
+ & +wbond*estr+wsccor*fact(1)*esccor+wliptran*eliptran
#endif
energia(0)=etot
energia(1)=evdw
energia(19)=esccor
energia(20)=edihcnstr
energia(21)=evdw_t
+ energia(22)=eliptran
+
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
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*fact(2)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
#else
do i=1,nct
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
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*fact(1)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
#endif
enddo
esccor=energia(19)
edihcnstr=energia(20)
estr=energia(18)
+ eliptran=energia(22)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
& wvdwpp,
& ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
& eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
& eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
- & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
+ & esccor,wsccor*fact(1),edihcnstr,ebr*nss,eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond,
& ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
& eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
& eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
- & edihcnstr,ebr*nss,etot
+ & edihcnstr,ebr*nss,eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic 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
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
ij=icant(itypi,itypj)
eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
rij=1.0D0/r_inv_ij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
ij=icant(itypi,itypj)
eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux
& /dabs(eps(itypi,itypj))
eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
endif
if (calc_grad) then
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),15(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+C returning the ith atom to box
+ 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
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+C returning jth atom to box
+ 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.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
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
+C checking the distance
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+C finding the closest
+ 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)
c write (iout,*) i,j,xj,yj,zj
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
+ sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
+ sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+ if (sss.le.0.0) cycle
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
+
call sc_angular
sigsq=1.0D0/sigsq
sig=sig0ij*dsqrt(sigsq)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0) then
- evdw=evdw+evdwij
+ if (bb.gt.0) then
+ evdw=evdw+evdwij*sss
else
- evdw_t=evdw_t+evdwij
+ evdw_t=evdw_t+evdwij*sss
endif
ij=icant(itypi,itypj)
aux=eps1*eps2rt**2*eps3rt**2
c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
c & aux*e2/eps(itypi,itypj)
c if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
#ifdef DEBUG
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij+e_augm
else
evdw_t=evdw_t+evdwij+e_augm
do i=1,nres
num_cont_hb(i)=0
enddo
-cd print '(a)','Enter EELEC'
-cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
+C print '(a)','Enter EELEC'
+C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
do i=1,nres
gel_loc_loc(i)=0.0d0
gcorr_loc(i)=0.0d0
enddo
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (i.eq.1) then
+ if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+ & .or. itype(i+2).eq.ntyp1) cycle
+ else
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+ & .or. itype(i+2).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ &) cycle
+ endif
if (itel(i).eq.0) goto 1215
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
+ 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)
+C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (j.eq.1) then
+ if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+ & .or.itype(j+2).eq.ntyp1
+ &) cycle
+ else
+ if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+ & .or.itype(j+2).eq.ntyp1
+ & .or.itype(j-1).eq.ntyp1
+ &) cycle
+ endif
+C
+C) cycle
if (itel(j).eq.0) goto 1216
ind=ind+1
iteli=itel(i)
dx_normj=dc_norm(1,j)
dy_normj=dc_norm(2,j)
dz_normj=dc_norm(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-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
rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(sqrt(rij))
+ sssgrad=sscagrad(sqrt(rij))
rrmij=1.0D0/rij
rij=dsqrt(rij)
rmij=1.0D0/rij
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
ees=ees+eesij
- evdw1=evdw1+evdwij
-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,
-cd & xmedi,ymedi,zmedi,xj,yj,zj
+ evdw1=evdw1+evdwij*sss
+c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)')
+c &'evdw1',i,j,evdwij
+c &,iteli,itelj,aaa,evdw1
+
+C write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+c write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
+c & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
+c & 1.0D0/dsqrt(rrmij),evdwij,eesij,
+c & xmedi,ymedi,zmedi,xj,yj,zj
C
C Calculate contributions to the Cartesian gradient.
C
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)
+ facvdw=-6*rrmij*(ev1+evdwij)*sss
facel=-3*rrmij*(el1+eesij)
fac1=fac
erij(1)=xj*rmij
gelc(l,k)=gelc(l,k)+ggg(l)
enddo
enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+C ggg(1)=facvdw*xj
+C ggg(2)=facvdw*yj
+C ggg(3)=facvdw*zj
+ if (sss.gt.0.0) then
+ ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+ ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+ else
+ ggg(1)=0.0
+ ggg(2)=0.0
+ ggg(3)=0.0
+ endif
do k=1,3
ghalf=0.5D0*ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)+ghalf
enddo
enddo
#else
- facvdw=ev1+evdwij
+ facvdw=(ev1+evdwij)*sss
facel=el1+eesij
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
C Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
& +a33*muij(4)
-cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
-cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+c write (iout,'(a6,2i5,0pf7.3)')
+c & 'eelloc',i,j,eel_loc_ij
+c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
eel_loc=eel_loc+eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
if (calc_grad) then
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))
-
+C Returning the ith atom to box
+ 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
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
c yj=c(2,nres+j)-yi
c zj=c(3,nres+j)-zi
C Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+C returning the jth atom to box
+ 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
+C Finding the closest jth atom
+ 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
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+C sss is scaling function for smoothing the cutoff gradient otherwise
+C the gradient would not be continuouse
+ sss=sscale(1.0d0/(dsqrt(rrij)))
+ if (sss.le.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+e1+e2
+ evdw2_14=evdw2_14+(e1+e2)*sss
endif
evdwij=e1+e2
-c write (iout,*) i,j,evdwij
- evdw2=evdw2+evdwij
+c write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
+c & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
+c & bad(itypj,iteli)
+ evdw2=evdw2+evdwij*sss
if (calc_grad) then
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
C
- fac=-(evdwij+e1)*rrij
+ fac=-(evdwij+e1)*rrij*sss
+ fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
double precision u(3),ud(3)
estr=0.0d0
estr1=0.0d0
- write (iout,*) "distchainmax",distchainmax
+c write (iout,*) "distchainmax",distchainmax
do i=nnt+1,nct
- if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
- estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
- do j=1,3
- gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
- & *dc(j,i-1)/vbld(i)
- enddo
- if (energy_dec) write(iout,*)
- & "estr1",i,vbld(i),distchainmax,
- & gnmr1(vbld(i),-1.0d0,distchainmax)
- else
+ if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
+C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+C do j=1,3
+C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
+C & *dc(j,i-1)/vbld(i)
+C enddo
+C if (energy_dec) write(iout,*)
+C & "estr1",i,vbld(i),distchainmax,
+C & gnmr1(vbld(i),-1.0d0,distchainmax)
+C else
+ if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+ diff = vbld(i)-vbldpDUM
+ else
diff = vbld(i)-vbldp0
c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+ endif
estr=estr+diff*diff
do j=1,3
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
enddo
- endif
-
+C endif
+C write (iout,'(a7,i5,4f7.3)')
+C & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
enddo
estr=0.5d0*AKP*estr+estr1
c
nbi=nbondterm(iti)
if (nbi.eq.1) then
diff=vbld(i+nres)-vbldsc0(1,iti)
-c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
-c & AKSC(1,iti),AKSC(1,iti)*diff*diff
+C write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+C & AKSC(1,iti),AKSC(1,iti)*diff*diff
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
c write (*,'(a,i2)') 'EBEND ICG=',icg
c write (iout,*) ithet_start,ithet_end
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
+C if (itype(i-1).eq.ntyp1) cycle
+ if (i.le.2) cycle
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
C Zero the energy function and its derivative at 0 or pi.
call splinthet(theta(i),0.5d0*delta,ss,ssd)
it=itype(i-1)
ichir21=isign(1,itype(i))
ichir22=isign(1,itype(i))
endif
+ if (i.eq.3) then
+ y(1)=0.0D0
+ y(2)=0.0D0
+ else
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+ if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
icrc=0
y(1)=0.0D0
y(2)=0.0D0
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ endif
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
icrc=0
& E_theta,E_tc)
endif
etheta=etheta+ethetai
+c write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
+c & 'ebend',i,ethetai,theta(i),itype(i)
c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i),
c & rad2deg*phii,rad2deg*phii1,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
etheta=0.0D0
c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
+C if (i.eq.2) cycle
+C if (itype(i-1).eq.ntyp1) cycle
+ if (i.le.2) cycle
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
+ if (iabs(itype(i+1)).eq.20) iblock=2
+ if (iabs(itype(i+1)).ne.20) iblock=1
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
theti2=0.5d0*theta(i)
- ityp2=ithetyp(iabs(itype(i-1)))
+ ityp2=ithetyp((itype(i-1)))
do k=1,nntheterm
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+ if (i.eq.3) then
+ phii=0.0d0
+ ityp1=nthetyp+1
+ do k=1,nsingle
+ cosph1(k)=0.0d0
+ sinph1(k)=0.0d0
+ enddo
+ else
+ if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
#else
phii=phi(i)
#endif
- ityp1=ithetyp(iabs(itype(i-2)))
+ ityp1=ithetyp((itype(i-2)))
do k=1,nsingle
cosph1(k)=dcos(k*phii)
sinph1(k)=dsin(k*phii)
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ endif
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
#else
phii1=phi(i+1)
#endif
- ityp3=ithetyp(iabs(itype(i)))
+ ityp3=ithetyp((itype(i)))
do k=1,nsingle
cosph2(k)=dcos(k*phii1)
sinph2(k)=dsin(k*phii1)
c write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
c call flush(iout)
- ethetai=aa0thet(ityp1,ityp2,ityp3)
+ ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
do k=1,ndouble
do l=1,k-1
ccl=cosph1(l)*cosph2(k-l)
enddo
endif
do k=1,ntheterm
- ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
- dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
+ ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
+ dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
& *coskt(k)
if (lprn)
- & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
+ & write (iout,*) "k",k,"
+ & aathet",aathet(k,ityp1,ityp2,ityp3,iblock),
& " ethetai",ethetai
enddo
if (lprn) then
endif
do m=1,ntheterm2
do k=1,nsingle
- aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
- & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
- & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
- & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+ aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
+ & +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
+ & +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
+ & +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
ethetai=ethetai+sinkt(m)*aux
dethetai=dethetai+0.5d0*m*aux*coskt(m)
dephii=dephii+k*sinkt(m)*(
- & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
- & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+ & ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
+ & bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
dephii1=dephii1+k*sinkt(m)*(
- & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
- & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+ & eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
+ & ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
if (lprn)
& write (iout,*) "m",m," k",k," bbthet",
- & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
- & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
- & ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
- & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ & bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
+ & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
+ & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
+ & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
enddo
enddo
if (lprn)
do m=1,ntheterm3
do k=2,ndouble
do l=1,k-1
- aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
- & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
- & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
- & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+ aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
ethetai=ethetai+sinkt(m)*aux
dethetai=dethetai+0.5d0*m*coskt(m)*aux
dephii=dephii+l*sinkt(m)*(
- & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
- & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
- & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
- & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
dephii1=dephii1+(k-l)*sinkt(m)*(
- & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
- & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
- & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
- & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
if (lprn) then
write (iout,*) "m",m," k",k," l",l," ffthet",
- & ffthet(l,k,m,ityp1,ityp2,ityp3),
- & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
- & ggthet(l,k,m,ityp1,ityp2,ityp3),
- & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ & ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock),
+ & " ethetai",ethetai
write (iout,*) cosph1ph2(l,k)*sinkt(m),
& cosph1ph2(k,l)*sinkt(m),
& sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
common /sccalc/ time11,time12,time112,theti,it,nlobit
delta=0.02d0*pi
escloc=0.0D0
-c write (iout,'(a)') 'ESC'
+C write (iout,*) 'ESC'
do i=loc_start,loc_end
it=itype(i)
if (it.eq.ntyp1) cycle
if (it.eq.10) goto 1
nlobit=nlob(iabs(it))
c print *,'i=',i,' it=',it,' nlobit=',nlobit
-c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
+C write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
theti=theta(i+1)-pipol
x(1)=dtan(theti)
x(2)=alph(i)
dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
enddo
dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
-c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c & esclocbi,ss,ssd
+ write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
+ & esclocbi,ss,ssd
escloci=ss*escloci+(1.0d0-ss)*esclocbi
c escloci=esclocbi
c write (iout,*) escloci
enddo
dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c & esclocbi,ss,ssd
+c & esclocbi,ss,ssd
escloci=ss*escloci+(1.0d0-ss)*esclocbi
-c write (iout,*) escloci
+C write (iout,*) 'i=',i, escloci
else
call enesc(x,escloci,dersc,ddummy,.false.)
endif
escloc=escloc+escloci
-c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+C write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+ write (iout,'(a6,i5,0pf7.3)')
+ & 'escloc',i,escloci
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& wscloc*dersc(1)
y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
enddo
do j = 1,3
- z_prime(j) = -uz(j,i-1)
+ z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
enddo
c write (2,*) "i",i
c write (2,*) "x_prime",(x_prime(j),j=1,3)
do j = 1,3
xx = xx + x_prime(j)*dc_norm(j,i+nres)
yy = yy + y_prime(j)*dc_norm(j,i+nres)
- zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
enddo
xxtab(i)=xx
Cc diagnostics - remove later
xx1 = dcos(alph(2))
yy1 = dsin(alph(2))*dcos(omeg(2))
- zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2))
+ zz1 = -dsign(1.0d0,itype(i))*dsin(alph(2))*dsin(omeg(2))
write(2,'(3f8.1,3f9.3,1x,3f9.3)')
& alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
& xx1,yy1,zz1
c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
c write (2,*) "escloc",escloc
+c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i),
+c & zz,xx,yy
if (.not. calc_grad) goto 1
#ifdef DEBUG
C
dZZ_Ci1(k)=0.0d0
dZZ_Ci(k)=0.0d0
do j=1,3
- dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
- dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+ dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
+ & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+ dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
+ & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
enddo
dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
- & .or. itype(i).eq.ntyp1) cycle
+ if (i.le.2) cycle
+ if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+C if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+C & .or. itype(i).eq.ntyp1) cycle
if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
if (iabs(itype(i)).eq.20) then
iblock=2
c lprn=.true.
etors_d=0.0D0
do i=iphi_start,iphi_end-1
- if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
- & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (i.le.3) cycle
+C if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+C & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
+ & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
+ & (itype(i+1).eq.ntyp1)) cycle
if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0)
& goto 1215
itori=itortyp(itype(i-2))
c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
esccor=0.0D0
do i=itau_start,itau_end
+ if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
esccor_ii=0.0D0
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
-c write (iout,*) "EBACK_SC_COR",i,esccor,intertyp
- gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
+c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp,
+c & nterm_sccor(isccori,isccori1),isccori,isccori1
+c gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
if (lprn)
& write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
& restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
- & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
- gsccor_loc(i-3)=gloci
+ & (v1sccor(j,1,itori,itori1),j=1,6)
+ & ,(v2sccor(j,1,itori,itori1),j=1,6)
+c gsccor_loc(i-3)=gloci
enddo !intertyp
enddo
return
integer dimen1,dimen2,atom,indx
double precision buffer(dimen1,dimen2)
double precision zapas
- common /contacts_hb/ zapas(3,20,maxres,7),
- & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
- & num_cont_hb(maxres),jcont_hb(20,maxres)
+ common /contacts_hb/ zapas(3,ntyp,maxres,7),
+ & facont_hb(ntyp,maxres),ees0p(ntyp,maxres),ees0m(ntyp,maxres),
+ & num_cont_hb(maxres),jcont_hb(ntyp,maxres)
num_kont=num_cont_hb(atom)
do i=1,num_kont
do k=1,7
integer dimen1,dimen2,atom,indx
double precision buffer(dimen1,dimen2)
double precision zapas
- common /contacts_hb/ zapas(3,20,maxres,7),
- & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
- & num_cont_hb(maxres),jcont_hb(20,maxres)
+ common /contacts_hb/ zapas(3,ntyp,maxres,7),
+ & facont_hb(ntyp,maxres),ees0p(ntyp,maxres),
+ & ees0m(ntyp,maxres),
+ & num_cont_hb(maxres),jcont_hb(ntyp,maxres)
num_kont=buffer(1,indx+26)
num_kont_old=num_cont_hb(atom)
num_cont_hb(atom)=num_kont+num_kont_old
include 'COMMON.GEO'
logical swap
double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
- & auxvec1(2),auxvec2(1),auxmat1(2,2)
+ & auxvec1(2),auxvec2(2),auxmat1(2,2)
logical lprn
common /kutas/ lprn
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
return
end
crc-------------------------------------------------
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+ 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.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+C this is done by Adasko
+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=1,nres
+C do i=1,1
+ if (itype(i).eq.ntyp1) cycle
+
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),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.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+ 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
+C print *, "doing sscalefor top part"
+C print *,i,sslip,fracinbuf,ssgradlip
+ else
+ eliptran=eliptran+pepliptran
+C 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
+CV do i=1,1
+ do i=1,nres
+ 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)
+C print *,positi,bordlipbot,buflipbot
+ 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))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+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*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i))
+C 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
+
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
SUBROUTINE MATVEC2(A1,V1,V2)
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
scalar=sc
return
end
+C-----------------------------------------------------------------------
+ double precision function sscale(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+ if(r.lt.r_cut-rlamb) then
+ sscale=1.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale=0d0
+ endif
+ return
+ end
+C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+ double precision function sscagrad(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+ if(r.lt.r_cut-rlamb) then
+ sscagrad=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscagrad=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscagrad=0.0d0
+ endif
+ return
+ end
+C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+ double precision function sscalelip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+C if(r.lt.r_cut-rlamb) then
+C sscale=1.0d0
+C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C gamm=(r-(r_cut-rlamb))/rlamb
+ sscalelip=1.0d0+r*r*(2*r-3.0d0)
+C else
+C sscale=0d0
+C endif
+ return
+ end
+C-----------------------------------------------------------------------
+ double precision function sscagradlip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+C if(r.lt.r_cut-rlamb) then
+C sscagrad=0.0d0
+C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C gamm=(r-(r_cut-rlamb))/rlamb
+ sscagradlip=r*(6*r-6.0d0)
+C else
+C sscagrad=0.0d0
+C endif
+ return
+ end
+
+C-----------------------------------------------------------------------
projects / unres.git / blobdiff