do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
c dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
chi1=chi(itypi,itypj)
evdw=evdw+evdwij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
& 'evdw',i,j,evdwij,' ss'
+C triple bond artifac removal
+ do k=j+1,iend(i,iint)
+C search over all next residues
+ if (dyn_ss_mask(k)) then
+C check if they are cysteins
+C write(iout,*) 'k=',k
+ call triple_ssbond_ene(i,j,k,evdwij)
+C call the energy function that removes the artifical triple disulfide
+C bond the soubroutine is located in ssMD.F
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,'tss'
+ endif!dyn_ss_mask(k)
+ enddo! k
ELSE
+C cycle
ind=ind+1
itypj=itype(j)
c dscj_inv=dsc_inv(itypj)
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
if (i.gt.3) then
#ifdef OSF
phii=phi(i)
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
dthett=thet_pred_mean*ssd
thet_pred_mean=thet_pred_mean*ss+a0thet(it)
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)
enddo
endif
if (i.lt.nres) then
+
+ if (iabs(itype(i+1)).eq.20) iblock=2
+ if (iabs(itype(i+1)).ne.20) iblock=1
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
sinph2(k)=0.0d0
enddo
endif
- 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)
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)
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
Cc diagnostics - remove later
xx1 = dcos(alph(2))
yy1 = dsin(alph(2))*dcos(omeg(2))
- zz1 = -dsin(alph(2))*dsin(omeg(2))
+ zz1 = -dsign(1.0, dfloat(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
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
- etors_ii=0.0D0
+ if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1) cycle
+ etors_ii=0.0D0
+ 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 (energy_dec) write (iout,'(a6,i5,0pf7.3)')
- & 'etor',i,etors_ii-v0(itori,itori1)
+ & 'etor',i,etors_ii-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,iblock),j=1,6),
+ & (v2(j,itori,itori1,iblock),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
enddo
lprn=.false.
c lprn=.true.
etors_d=0.0D0
+c write(iout,*) "a tu??"
do i=iphid_start,iphid_end
+ if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
itori2=itortyp(itype(i))
phii1=phi(i+1)
gloci1=0.0D0
gloci2=0.0D0
- 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)
+ 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,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)
enddo
gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
-c write (iout,*) "gloci", gloc(i-3,icg)
enddo
return
end
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
-c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
+c write (iout,*) "EBACK_SC_COR",itau_start,itau_end
esccor=0.0D0
do i=itau_start,itau_end
esccor_ii=0.0D0
+ if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
phii=phi(i)
c 3 = SC...Ca...Ca...SCi
gloci=0.0D0
if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
- & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
- & (itype(i-1).eq.21)))
+ & (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.21)))
+ & .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)
+ & .or.(itype(i).eq.ntyp1)))
& .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
- & (itype(i-1).eq.21)))) cycle
- if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
- if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
+ & (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+ & (itype(i-3).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)
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,intertyp,itori,itori1),j=1,6)
- & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
+ & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,
+ & (v1sccor(j,intertyp,isccori,isccori1),j=1,6)
+ & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6)
gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
enddo !intertyp
enddo
projects / unres.git / blobdiff