c write(iout,*) "TEST_ENE1 ehomology_constr=",ehomology_constr
C BARTEK for dfa test!
- if (wdfa_dist.gt.0) call edfad(edfadis)
+ if (wdfa_dist.gt.0) then
+ call edfad(edfadis)
+ else
+ edfadis=0.0d0
+ endif
c write(iout,*)'edfad is finished!', wdfa_dist,edfadis
- if (wdfa_tor.gt.0) call edfat(edfator)
+ if (wdfa_tor.gt.0) then
+ call edfat(edfator)
+ else
+ edfator=0.0d0
+ endif
c write(iout,*)'edfat is finished!', wdfa_tor,edfator
- if (wdfa_nei.gt.0) call edfan(edfanei)
+ if (wdfa_nei.gt.0) then
+ call edfan(edfanei)
+ else
+ edfanei=0.0d0
+ endif
c write(iout,*)'edfan is finished!', wdfa_nei,edfanei
- if (wdfa_beta.gt.0) call edfab(edfabet)
+ if (wdfa_beta.gt.0) then
+ call edfab(edfabet)
+ else
+ edfabet=0.0d0
+ endif
c write(iout,*)'edfab is finished!', wdfa_beta,edfabet
c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
+ include 'COMMON.IOUNITS'
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
- 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 !end dhpb1(i).gt.0
- endif !end const_dist=11
- 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)
-C write(iout,*) "after",dd
- 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
-C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i))
-C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd
-C print *,ehpb,"tu?"
-C write(iout,*) ehpb,"btu?",
-C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i)
-C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
-C & ehpb,fordepth(i),dd
- else
- if (dhpb1(i).gt.0.0d0) then
- ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ 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
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+c if (energy_dec) write (iout,'(a6,2i5,6f10.3,i5)')
+c & "edisL",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+c & 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)
+c if (energy_dec) write(iout,'(a6,2i5,6f10.3,i5)')
+c & "edisX",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+c & -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))
+c if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)')
+c & "edisQ",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+c & 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
+c if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)')
+c & "edisS",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+c & 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
do k=1,3
ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
enddo
endif
enddo
- if (constr_dist.ne.11) ehpb=0.5D0*ehpb
return
end
C--------------------------------------------------------------------------
j = jres_homo(ii)
dij=dist(i,j)
c write (iout,*) "dij(",i,j,") =",dij
+ nexl=0
do k=1,constr_homology
- if(.not.l_homo(k,ii)) cycle
+ if(.not.l_homo(k,ii)) then
+ nexl=nexl+1
+ cycle
+ endif
distance(k)=odl(k,ii)-dij
c write (iout,*) "distance(",k,") =",distance(k)
c
write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
write (iout,* )"min_odl",min_odl
#endif
+#ifdef OLDRESTR
odleg2=0.0d0
+#else
+ if (waga_dist.ge.0.0d0) then
+ odleg2=nexl
+ else
+ odleg2=0.0d0
+ endif
+#endif
do k=1,constr_homology
c Nie wiem po co to liczycie jeszcze raz!
c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
c & -(6.28318-dih_diff(i,k))
c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
c & 6.28318+dih_diff(i,k)
-
+#ifdef OLD_DIHED
kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#else
+ kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i)
+#endif
c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
gdih(k)=dexp(kat3)
kat2=kat2+gdih(k)
sum_gdih=kat2
sum_sgdih=0.0d0
do k=1,constr_homology
+#ifdef OLD_DIHED
sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+#else
+ sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i)
+#endif
c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
sum_sgdih=sum_sgdih+sgdih
enddo