include 'COMMON.CONTROL'
dimension ggg(3)
ehpb=0.0D0
+ do i=1,3
+ ggg(i)=0.0d0
+ enddo
+C write (iout,*) ,"link_end",link_end,constr_dist
cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
-cd write(iout,*)'link_start=',link_start,' link_end=',link_end
+c write(iout,*)'link_start=',link_start,' link_end=',link_end,
+c & " constr_dist",constr_dist
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
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.
+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
+ 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
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 (constr_dist.eq.11) then
- ehpb=ehpb+fordepth(i)**4.0d0
- & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
- fac=fordepth(i)**4.0d0
- & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
- if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)')
- & "edisl",ii,jj,
- & fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
- & fordepth(i),dd
- else
- 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
- 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
+cd & ' waga=',waga,' fac=',fac
+! else if (ii.gt.nres .and. jj.gt.nres) then
+ else
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
- if (constr_dist.eq.11) then
- ehpb=ehpb+fordepth(i)**4.0d0
+ if (irestr_type(i).eq.11) then
+ ehpb=ehpb+fordepth(i)!**4.0d0
& *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
- fac=fordepth(i)**4.0d0
+ fac=fordepth(i)!**4.0d0
& *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
- if (energy_dec) write (iout,'(a6,2i5,f15.6,2f8.3)')
- 7 "edisl",ii,jj,
- & fordepth(i)**4.0d0*rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)),
- & fordepth(i),dd
-c if (energy_dec)
-c & write (iout,*) fac
- else
- if (dhpb1(i).gt.0.0d0) then
- ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ if (energy_dec) write (iout,'(a6,2i5,6f10.3,i5)')
+ & "edisL",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+ & ehpb,irestr_type(i)
+ else if (irestr_type(i).eq.10) then
+c AL 6//19/2018 cross-link restraints
+ xdis = 0.5d0*(dd/forcon(i))**2
+ expdis = dexp(-xdis)
+c aux=(dhpb(i)+dhpb1(i)*xdis)*expdis+fordepth(i)
+ aux=(dhpb(i)+dhpb1(i)*xdis*xdis)*expdis+fordepth(i)
+c write (iout,*)"HERE: xdis",xdis," expdis",expdis," aux",aux,
+c & " wboltzd",wboltzd
+ ehpb=ehpb-wboltzd*xlscore(i)*dlog(aux)
+c fac=-wboltzd*(dhpb1(i)*(1.0d0-xdis)-dhpb(i))
+ fac=-wboltzd*xlscore(i)*(dhpb1(i)*(2.0d0-xdis)*xdis-dhpb(i))
+ & *expdis/(aux*forcon(i)**2)
+ if (energy_dec) write(iout,'(a6,2i5,6f10.3,i5)')
+ & "edisX",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+ & -wboltzd*xlscore(i)*dlog(aux),irestr_type(i)
+ else if (irestr_type(i).eq.2) then
+c Quartic restraints
+ ehpb=ehpb+forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)')
+ & "edisQ",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+ & forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)),irestr_type(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
+c Quadratic restraints
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,*) "alpha reg",dd,waga*rdis*rdis
+ ehpb=ehpb+0.5d0*waga*rdis*rdis
+ if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)')
+ & "edisS",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+ & 0.5d0*waga*rdis*rdis,irestr_type(i)
C
C Evaluate gradient.
C
fac=waga*rdis/dd
endif
- 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
+c Calculate Cartesian gradient
+ 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
- do j=1,3
- ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
- ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
- enddo
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
endif
cgrad do j=iii,jjj-1
cgrad do k=1,3
enddo
endif
enddo
- if (constr_dist.ne.11) ehpb=0.5D0*ehpb
-c do i=1,nres
-c write (iout,*) "ghpbc",i,(ghpbc(j,i),j=1,3)
-c enddo
return
end
C--------------------------------------------------------------------------
projects / unres.git / blobdiff