X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Fcluster%2Fwham%2Fsrc%2Fenergy_p_new.F;h=8006a41109e4b94193e6d85c0fdbe48987a71ca0;hb=f37cfff89ef75a184b0c7bbbac2ecea2efc1d157;hp=8ec2d70bd6be64bcaf757db147d44b82adfaaed6;hpb=88c9804a273df101470add9049bec50351c39fdb;p=unres.git diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F index 8ec2d70..8006a41 100644 --- a/source/cluster/wham/src/energy_p_new.F +++ b/source/cluster/wham/src/energy_p_new.F @@ -262,7 +262,7 @@ cd write (iout,*) i,g_corr5_loc(i) & +wsccor*fact(1)*gsccor_loc(i) enddo endif -cd call enerprint(energia(0),fact) +c call enerprint(energia(0),fact) cd call intout cd stop return @@ -2876,6 +2876,7 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr @@ -2908,6 +2909,12 @@ 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 @@ -2925,7 +2932,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd - endif + 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 @@ -2941,6 +2949,20 @@ C 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)) fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd @@ -2958,6 +2980,7 @@ 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 @@ -2978,7 +3001,7 @@ C Cartesian gradient in the SC vectors (ghpbx). enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -3106,7 +3129,7 @@ c include 'COMMON.SETUP' include 'COMMON.NAMES' - do i=1,19 + do i=1,max_template distancek(i)=9999999.9 enddo @@ -3126,7 +3149,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 @@ -3146,7 +3174,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 @@ -3154,11 +3192,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 @@ -3209,6 +3256,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 @@ -3299,7 +3347,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) @@ -3308,8 +3356,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) @@ -3337,7 +3388,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 @@ -3646,9 +3701,12 @@ c & waga_theta*Eval+waga_d*Erot)*waga_homology(iset) c write (iout,*) "ehomology_constr=",ehomology_constr endif -c write (iout,*) "odleg",odleg," kat",kat," Eval",Eval," Erot",Erot -c write (iout,*) "ehomology_constr",ehomology_constr -c ehomology_constr=odleg+kat+Uconst_back +#ifdef DEBUG + write (iout,*) "iset",iset," waga_homology",waga_homology(iset) + write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat, + & " Eval",waga_theta,Eval," Erot",waga_d,Erot + write (iout,*) "ehomology_constr",ehomology_constr +#endif return 748 format(a8,f12.3,a6,f12.3,a7,f12.3) @@ -4015,7 +4073,7 @@ C coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-3).ne.ntyp1) 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 @@ -5130,6 +5188,7 @@ c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor 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) @@ -5163,6 +5222,9 @@ c 3 = SC...Ca...Ca...SCi 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 @@ -5175,6 +5237,9 @@ c &gloc_sc(intertyp,i-3,icg) & ,(v2sccor(j,intertyp,itori,itori1),j=1,6) gsccor_loc(i-3)=gsccor_loc(i-3)+gloci enddo !intertyp +#ifdef DEBUG + write (iout,*) "i",i,(tauangle(j,i),j=1,3),esccor_ii +#endif enddo return