#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*fact(1)*ees+wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +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
#else
etot=wsc*evdw+wscp*evdw2+welec*fact(1)*(ees+evdw1)
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +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
energia(18)=estr
energia(19)=esccor
energia(20)=edihcnstr
+cc if (dyn_ss) call dyn_set_nss
c detecting NaNQ
i=0
#ifdef WINPGI
& +wturn4*fact(3)*gel_loc_turn4(i)
& +wturn3*fact(2)*gel_loc_turn3(i)
& +wturn6*fact(5)*gel_loc_turn6(i)
- & +wel_loc*fact(2)*gel_loc_loc(i)+
+ & +wel_loc*fact(2)*gel_loc_loc(i)
& +wsccor*fact(1)*gsccor_loc(i)
enddo
endif
cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
c endif
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
dscj_inv=vbld_inv(j+nres)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
logical lprn
common /srutu/icall
integer icant
c if (icall.gt.0) lprn=.true.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,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
+c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
dscj_inv=vbld_inv(j+nres)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! SSBOND
enddo ! j
enddo ! iint
enddo ! i
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
common /srutu/ icall
logical lprn
integer icant
c if (icall.gt.0) lprn=.true.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+C in case of diagnostics write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
dscj_inv=vbld_inv(j+nres)
sig0ij=sigma(itypi,itypj)
r0ij=r0(itypi,itypj)
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
C Uncomment following three lines for SC-p interactions
c xj=c(1,nres+j)-xi
c yj=c(2,nres+j)-yi
C
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
- include 'sizesclu.dat'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
dimension ggg(3)
ehpb=0.0D0
-cd print *,'edis: nhpb=',nhpb,' fbr=',fbr
-cd print *,'link_start=',link_start,' link_end=',link_end
+cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
+cd write(iout,*)'link_start=',link_start,' link_end=',link_end
if (link_end.eq.0) return
do i=link_start,link_end
C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
iii=ii
jjj=jj
endif
+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. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+ if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
+ if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+ & iabs( itype(jjj)).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+cd write (iout,*) "eij",eij
+ endif
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ dd=dist(ii,jj)
+ rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+ waga=forcon(i)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
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)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "alph nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ 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 write (iout,*) "alpha reg",dd,waga*rdis*rdis
C
C Evaluate gradient.
C
- fac=waga*rdis/dd
+ fac=waga*rdis/dd
+ endif
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
- do j=iii,jjj-1
+ endif
do k=1,3
- ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
enddo
- enddo
endif
enddo
ehpb=0.5D0*ehpb
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
- itypi=itype(i)
+ itypi=iabs(itype(i))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=dsc_inv(itypi)
- itypj=itype(j)
+ itypj=iabs(itype(j))
dscj_inv=dsc_inv(itypj)
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
deltat12=om2-om1+2.0d0
cosphi=om12-om1*om2
eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
- & +akct*deltad*deltat12
+ & +akct*deltad*deltat12+ebr
& +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
c
do i=nnt,nct
- iti=itype(i)
+ iti=iabs(itype(i))
if (iti.ne.10) then
nbi=nbondterm(iti)
if (nbi.eq.1) then
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)
+ ichir1=isign(1,itype(i-2))
+ ichir2=isign(1,itype(i))
+ if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1))
+ if (itype(i).eq.10) ichir2=isign(1,itype(i-1))
+ if (itype(i-1).eq.10) then
+ itype1=isign(10,itype(i-2))
+ ichir11=isign(1,itype(i-2))
+ ichir12=isign(1,itype(i-2))
+ itype2=isign(10,itype(i))
+ ichir21=isign(1,itype(i))
+ ichir22=isign(1,itype(i))
+ endif
c if (i.gt.ithet_start .and.
c & (itel(i-1).eq.0 .or. itel(i-2).eq.0)) goto 1215
c if (i.gt.3 .and. (i.le.4 .or. itel(i-3).ne.0)) then
C In following comments this theta will be referred to as t_c.
thet_pred_mean=0.0d0
do k=1,2
- athetk=athet(k,it)
- bthetk=bthet(k,it)
+ athetk=athet(k,it,ichir1,ichir2)
+ bthetk=bthet(k,it,ichir1,ichir2)
+ if (it.eq.10) then
+ athetk=athet(k,itype1,ichir11,ichir12)
+ bthetk=bthet(k,itype2,ichir21,ichir22)
+ endif
thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
enddo
c write (iout,*) "thet_pred_mean",thet_pred_mean
thet_pred_mean=thet_pred_mean*ss+a0thet(it)
c write (iout,*) "thet_pred_mean",thet_pred_mean
C Derivatives of the "mean" values in gamma1 and gamma2.
- dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
- dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
+ dthetg1=(-athet(1,it,ichir1,ichir2)*y(2)
+ &+athet(2,it,ichir1,ichir2)*y(1))*ss
+ dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2)
+ & +bthet(2,it,ichir1,ichir2)*z(1))*ss
+ if (it.eq.10) then
+ dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2)
+ &+athet(2,itype1,ichir11,ichir12)*y(1))*ss
+ dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2)
+ & +bthet(2,itype2,ichir21,ichir22)*z(1))*ss
+ endif
if (theta(i).gt.pi-delta) then
call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
& E_tc0)
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
+ 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(itype(i-1))
+ ityp2=ithetyp((itype(i-1)))
do k=1,nntheterm
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
#else
phii=phi(i)
#endif
- ityp1=ithetyp(itype(i-2))
+ ityp1=ithetyp((itype(i-2)))
do k=1,nsingle
cosph1(k)=dcos(k*phii)
sinph1(k)=dsin(k*phii)
#else
phii1=phi(i+1)
#endif
- ityp3=ithetyp(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)
do i=loc_start,loc_end
it=itype(i)
if (it.eq.10) goto 1
- nlobit=nlob(it)
+ nlobit=nlob(iabs(it))
c print *,'i=',i,' it=',it,' nlobit=',nlobit
c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
theti=theta(i+1)-pipol
do iii=-1,1
do j=1,nlobit
- expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
+ expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
cd print *,'j=',j,' expfac=',expfac
escloc_i=escloc_i+expfac
do k=1,3
dersc12=0.0d0
do j=1,nlobit
- expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
+ expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin)
escloc_i=escloc_i+expfac
do k=1,2
dersc(k)=dersc(k)+Ax(k,j)*expfac
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=itype(i)
+ it=iabs(itype(i))
if (it.eq.10) goto 1
c
C Compute the axes of tghe local cartesian coordinates system; store in
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)
xx = xx + x_prime(j)*dc_norm(j,i+nres)
yy = yy + y_prime(j)*dc_norm(j,i+nres)
zz = zz + z_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
enddo
xxtab(i)=xx
yytab(i)=yy
zztab(i)=zz
+
C
C Compute the energy of the ith side cbain
C
c write (2,*) "xx",xx," yy",yy," zz",zz
- it=itype(i)
+ it=iabs(itype(i))
do j = 1,65
x(j) = sc_parmin(j,it)
enddo
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))
do i=1,ndih_constr
itori=idih_constr(i)
phii=phi(itori)
- difi=phii-phi0(i)
+ difi=pinorm(phii-phi0(i))
if (difi.gt.drange(i)) then
difi=difi-drange(i)
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
endif
-! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+c write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
+c & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
enddo
-! write (iout,*) 'edihcnstr',edihcnstr
+ write (iout,*) 'edihcnstr',edihcnstr
return
end
c------------------------------------------------------------------------------
etors=0.0D0
do i=iphi_start,iphi_end
if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
+ if (iabs(itype(i)).eq.20) then
+ iblock=2
+ else
+ iblock=1
+ endif
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
phii=phi(i)
gloci=0.0D0
C Regular cosine and sine terms
- do j=1,nterm(itori,itori1)
- v1ij=v1(j,itori,itori1)
- v2ij=v2(j,itori,itori1)
+ do j=1,nterm(itori,itori1,iblock)
+ v1ij=v1(j,itori,itori1,iblock)
+ v2ij=v2(j,itori,itori1,iblock)
cosphi=dcos(j*phii)
sinphi=dsin(j*phii)
etors=etors+v1ij*cosphi+v2ij*sinphi
C
cosphi=dcos(0.5d0*phii)
sinphi=dsin(0.5d0*phii)
- do j=1,nlor(itori,itori1)
+ do j=1,nlor(itori,itori1,iblock)
vl1ij=vlor1(j,itori,itori1)
vl2ij=vlor2(j,itori,itori1)
vl3ij=vlor3(j,itori,itori1)
gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
enddo
C Subtract the constant term
- etors=etors-v0(itori,itori1)
+ etors=etors-v0(itori,itori1,iblock)
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,
- & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
+ & (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
1215 continue
enddo
! 6/20/98 - dihedral angle constraints
edihcnstr=0.0d0
+c write (iout,*) "Dihedral angle restraint energy"
do i=1,ndih_constr
- print *,"i",i
itori=idih_constr(i)
phii=phi(itori)
- difi=phii-phi0(i)
+ difi=pinorm(phii-phi0(i))
+c write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
+c & rad2deg*difi,rad2deg*phi0(i),rad2deg*drange(i)
if (difi.gt.drange(i)) then
difi=difi-drange(i)
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+c write (iout,*) 0.25d0*ftors*difi**4
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+c write (iout,*) 0.25d0*ftors*difi**4
endif
-! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
enddo
-! write (iout,*) 'edihcnstr',edihcnstr
+c write (iout,*) 'edihcnstr',edihcnstr
return
end
c----------------------------------------------------------------------------
phii1=phi(i+1)
gloci1=0.0D0
gloci2=0.0D0
+ iblock=1
+ if (iabs(itype(i+1)).eq.20) iblock=2
C Regular cosine and sine terms
- do j=1,ntermd_1(itori,itori1,itori2)
- v1cij=v1c(1,j,itori,itori1,itori2)
- v1sij=v1s(1,j,itori,itori1,itori2)
- v2cij=v1c(2,j,itori,itori1,itori2)
- v2sij=v1s(2,j,itori,itori1,itori2)
+ do j=1,ntermd_1(itori,itori1,itori2,iblock)
+ v1cij=v1c(1,j,itori,itori1,itori2,iblock)
+ v1sij=v1s(1,j,itori,itori1,itori2,iblock)
+ v2cij=v1c(2,j,itori,itori1,itori2,iblock)
+ v2sij=v1s(2,j,itori,itori1,itori2,iblock)
cosphi1=dcos(j*phii)
sinphi1=dsin(j*phii)
cosphi2=dcos(j*phii1)
gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
- do k=2,ntermd_2(itori,itori1,itori2)
+ do k=2,ntermd_2(itori,itori1,itori2,iblock)
do l=1,k-1
- v1cdij = v2c(k,l,itori,itori1,itori2)
- v2cdij = v2c(l,k,itori,itori1,itori2)
- v1sdij = v2s(k,l,itori,itori1,itori2)
- v2sdij = v2s(l,k,itori,itori1,itori2)
+ v1cdij = v2c(k,l,itori,itori1,itori2,iblock)
+ v2cdij = v2c(l,k,itori,itori1,itori2,iblock)
+ v1sdij = v2s(k,l,itori,itori1,itori2,iblock)
+ v2sdij = v2s(l,k,itori,itori1,itori2,iblock)
cosphi1p2=dcos(l*phii+(k-l)*phii1)
cosphi1m2=dcos(l*phii-(k-l)*phii1)
sinphi1p2=dsin(l*phii+(k-l)*phii1)
c lprn=.true.
c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
esccor=0.0D0
- do i=iphi_start,iphi_end
+ do i=itau_start,itau_end
esccor_ii=0.0D0
- itori=itype(i-2)
- itori1=itype(i-1)
+ isccori=isccortyp(itype(i-2))
+ isccori1=isccortyp(itype(i-1))
phii=phi(i)
+cccc Added 9 May 2012
+cc Tauangle is torsional engle depending on the value of first digit
+c(see comment below)
+cc Omicron is flat angle depending on the value of first digit
+c(see comment below)
+
+
+ do intertyp=1,3 !intertyp
+cc Added 09 May 2012 (Adasko)
+cc Intertyp means interaction type of backbone mainchain correlation:
+c 1 = SC...Ca...Ca...Ca
+c 2 = Ca...Ca...Ca...SC
+c 3 = SC...Ca...Ca...SCi
gloci=0.0D0
- do j=1,nterm_sccor
- v1ij=v1sccor(j,itori,itori1)
- v2ij=v2sccor(j,itori,itori1)
- cosphi=dcos(j*phii)
- sinphi=dsin(j*phii)
+ if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
+ & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or.
+ & (itype(i-1).eq.ntyp1)))
+ & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+ & .or.(itype(i-2).eq.ntyp1)))
+ & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
+ & (itype(i-1).eq.ntyp1)))) cycle
+ if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
+ if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1))
+ & cycle
+ do j=1,nterm_sccor(isccori,isccori1)
+ v1ij=v1sccor(j,intertyp,isccori,isccori1)
+ v2ij=v2sccor(j,intertyp,isccori,isccori1)
+ cosphi=dcos(j*tauangle(intertyp,i))
+ sinphi=dsin(j*tauangle(intertyp,i))
esccor=esccor+v1ij*cosphi+v2ij*sinphi
- gloci=gloci+fact*j*(v2ij*cosphi-v1ij*sinphi)
+ gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
+ gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
+c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
+c &gloc_sc(intertyp,i-3,icg)
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,intertyp,itori,itori1),j=1,6)
+ & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
+ gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+ enddo !intertyp
enddo
+
return
end
c------------------------------------------------------------------------------
C Set lprn=.true. for debugging
lprn=.false.
eturn6=0.0d0
+ ecorr6=0.0d0
#ifdef MPL
n_corr=0
n_corr1=0
cd write(2,*)'ijkl',i,j,i+1,j1
if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3
& .or. wturn6.eq.0.0d0))then
-cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
- ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
-cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
-cd & 'ecorr6=',ecorr6
+c write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
+c ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
+c write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
+c & 'ecorr6=',ecorr6, wcorr6
cd write (iout,'(4e15.5)') sred_geom,
cd & dabs(eello4(i,j,i+1,j1,jj,kk)),
cd & dabs(eello5(i,j,i+1,j1,jj,kk)),
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
projects / unres.git / blobdiff