X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fenergy_p_new.F;h=0e50c7f00151b05edb9e450bdc686f007104c99c;hb=147d86ad190bc360fa4b0c9a4c7d6858070c506f;hp=a3e2b12dea69d757601571c3bc7f8b810da13427;hpb=ca5951f307725bdad7b4e6a390c4cdfde1b4e17a;p=unres.git diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F index a3e2b12..0e50c7f 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -116,13 +116,29 @@ c write(iout,*) "TEST_ENE1 constr_homology=",constr_homology 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 @@ -2948,16 +2964,23 @@ C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' + include 'DIMENSIONS.FREE' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' + 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 @@ -2976,79 +2999,81 @@ 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. +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 (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 +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 (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 -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) @@ -3056,7 +3081,6 @@ C Cartesian gradient in the SC vectors (ghpbx). enddo endif enddo - ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -3200,7 +3224,12 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d 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)) then + nexl=nexl+1 + cycle + endif distance(k)=odl(k,ii)-dij c write (iout,*) "distance(",k,") =",distance(k) c @@ -3220,7 +3249,17 @@ c endif enddo - min_odl=minval(distancek) +c min_odl=minval(distancek) + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + & min_odl=distancek(kk) + enddo c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij @@ -3228,11 +3267,20 @@ c write (iout,* )"min_odl",min_odl 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 & (2*(sigma_odl(i,j,k))**2)) + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c c For Gaussian-type Urestr @@ -3283,6 +3331,7 @@ c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) c & *waga_dist)+min_odl c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist c + if(.not.l_homo(k,ii)) cycle if (waga_dist.ge.0.0d0) then c For Gaussian-type Urestr c @@ -3373,7 +3422,7 @@ c write (iout,*) idihconstr_start_homo,idihconstr_end_homo do i=idihconstr_start_homo,idihconstr_end_homo kat2=0.0d0 c betai=beta(i,i+1,i+2,i+3) - betai = phi(i+3) + betai = phi(i) c write (iout,*) "betai =",betai do k=1,constr_homology dih_diff(k)=pinorm(dih(k,i)-betai) @@ -3382,8 +3431,11 @@ c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= 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) @@ -3411,7 +3463,11 @@ c ---------------------------------------------------------------------- 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