etot=wsc*(evdw+fact(6)*evdw_t)+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
etot=wsc*(evdw+fact(6)*evdw_t)+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(19)=esccor
energia(20)=edihcnstr
energia(21)=evdw_t
+c if (dyn_ss) call dyn_set_nss
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
include 'COMMON.ENEPS'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
logical lprn
common /srutu/icall
integer icant
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)
+C /06/28/2013 Adasko: In case of dyn_ss - dynamic disulfide bond
+C formation no electrostatic interactions should be calculated. If it
+C would be allowed NaN would appear
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+C /06/28/2013 Adasko: dyn_ss_mask is logical statement wheather this Cys
+C residue can or cannot form disulfide bond. There is still bug allowing
+C Cys...Cys...Cys bond formation
+ call dyn_ssbond_ene(i,j,evdwij)
+C /06/28/2013 Adasko: dyn_ssbond_ene is dynamic SS bond foration energy
+C function in ssMD.F
+ 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)
dscj_inv=vbld_inv(j+nres)
c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj,
c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
c & aux*e2/eps(itypi,itypj)
+c write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
C
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
- include 'DIMENSIONS.ZSCOPT'
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 (.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 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
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
+c & +akct*deltad*deltat12
& +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 & " deltat12",deltat12," eij",eij
+ write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
+ & " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
+ & " deltat12",deltat12," eij",eij,"ebr",ebr
ed=2*akcm*deltad+akct*deltat12
pom1=akct*deltad
pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
double precision u(3),ud(3)
+ logical :: lprn=.false.
estr=0.0d0
do i=nnt+1,nct
diff = vbld(i)-vbldp0
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
+ if (lprn)
+ & write (iout,*) 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
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
usum=usum+uprod1
usumsqder=usumsqder+ud(j)*uprod2
enddo
-c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
-c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
+ if (lprn)
+ & write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
+ & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
estr=estr+uprod/usum
do j=1,3
gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
etheta=0.0D0
c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
do i=ithet_start,ithet_end
+ if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+ &(itype(i).eq.ntyp1)) cycle
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3) then
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
enddo
else
phii=0.0d0
- ityp1=nthetyp+1
+ ityp1=ithetyp(itype(i-2))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres) then
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp(itype(i))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
enddo
enddo
10 continue
- if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
- & i,theta(i)*rad2deg,phii*rad2deg,
+c lprn1=.true.
+ if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)')
+ & 'ebe',i,theta(i)*rad2deg,phii*rad2deg,
& phii1*rad2deg,ethetai
+c lprn1=.false.
etheta=etheta+ethetai
+
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
gloc(nphi+i-2,icg)=wang*dethetai
escloc = escloc + sumene
c write (2,*) "escloc",escloc
if (.not. calc_grad) goto 1
-#ifdef DEBUG
+
+#ifdef DEBUG2
C
C This section to check the numerical derivatives of the energy of ith side
C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
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,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)
+
+ 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)
+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.21).or.
+ & (itype(i-1).eq.21)))
+ & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+ & .or.(itype(i-2).eq.21)))
+ & .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))
+ & 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
+#ifdef DEBUG
+ esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
+#endif
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
+c do i=1,nres
+c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc(i,icg)
+c enddo
return
end
c------------------------------------------------------------------------------
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