X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=ea2fe9e313792eb6e44c511643a8ab390014394e;hb=c987a84c8742686406ea6b8807c9132098dc567e;hp=205f78a063dafb7ffcffd877c463e97e5151ad67;hpb=d24fd92f8762d3b48a94c699ef90e0ab37d3c204;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 205f78a..ea2fe9e 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -6895,7 +6895,7 @@ c c - do i=1,19 + do i=1,max_template distancek(i)=9999999.9 enddo @@ -6914,9 +6914,13 @@ 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 c write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii) - 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 @@ -6955,7 +6959,15 @@ 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)/ @@ -7110,8 +7122,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) ! waga_angle rmvd from Gaussian argument +#endif c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i) gdih(k)=dexp(kat3) kat2=kat2+gdih(k) @@ -7138,7 +7153,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) ! waga_angle rmvd +#endif c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle sum_sgdih=sum_sgdih+sgdih enddo @@ -7220,7 +7239,7 @@ c utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta? gtheta(k)=dexp(utheta_i) ! + min_utheta_i? - gutheta_i=gutheta_i+dexp(utheta_i) ! Sum of Gaussians (pk) + gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk) c Gradient for single Gaussian restraint in subr Econstr_back c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1) c @@ -7292,7 +7311,9 @@ c c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? c uscdiffk(k)=usc_diff(i) guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff - guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk) +c write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i), +c & " guscdiff2",guscdiff2(k) + guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk) c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j), c & xxref(j),yyref(j),zzref(j) enddo @@ -11050,6 +11071,42 @@ C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist end c---------------------------------------------------------------------------- + double precision function sscale2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) +c write (2,*) "r",r," r_cut",r_cut," r0",r0," rlamb",rlamb +c write (2,*) "rr",rr + if(rr.lt.r_cut-rlamb) then + sscale2=1.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0) + else + sscale2=0d0 + endif + return + end +C----------------------------------------------------------------------- + double precision function sscalgrad2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) + if(rr.lt.r_cut-rlamb) then + sscalgrad2=0.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + if (r.ge.r0) then + sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb + else + sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb + endif + else + sscalgrad2=0.0d0 + endif + return + end +c---------------------------------------------------------------------------- subroutine e_saxs(Esaxs_constr) implicit none include 'DIMENSIONS' @@ -11084,6 +11141,7 @@ c & sigma2CACA,sigma2CASC,sigma2SCCA,sigma2SCSC,expCACA,expCASC, & expSCCA,expSCSC,CASCgrad,SCCAgrad,SCSCgrad,aux,auxC,auxC1, & auxX,auxX1,CACAgrad,Cnorm + double precision sss2,ssgrad2,rrr,sscalgrad2,sscale2 double precision dist external dist c SAXS restraint penalty function @@ -11177,22 +11235,46 @@ c SC SC enddo ! k #else dijCACA=dist(i,j) - sigma2CACA=0.25d0/(restok(itype(j))**2+restok(itype(i))**2) + sigma2CACA=scal_rad**2*0.25d0/ + & (restok(itype(j))**2+restok(itype(i))**2) + + IF (saxs_cutoff.eq.0) THEN do k=1,nsaxs dk = distsaxs(k) expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2) Pcalc(k) = Pcalc(k)+expCACA + CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA + do l=1,3 + aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + enddo ! l + enddo ! k + ELSE + rrr = saxs_cutoff*2.0d0/dsqrt(sigma2CACA) + do k=1,nsaxs + dk = distsaxs(k) +c write (2,*) "ijk",i,j,k + sss2 = sscale2(dijCACA,rrr,dk,0.3d0) + if (sss2.eq.0.0d0) cycle + ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0) + if (energy_dec) write(iout,'(a4,3i5,5f10.4)') + & 'saxs',i,j,k,dijCACA,rrr,dk,sss2,ssgrad2 + expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)*sss2 + Pcalc(k) = Pcalc(k)+expCACA #ifdef DEBUG write(iout,*) "i j k Pcalc",i,j,Pcalc(k) #endif - CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA + CACAgrad = -sigma2CACA*(dijCACA-dk)*expCACA+ + & ssgrad2*expCACA/sss2 do l=1,3 c CA CA aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA - PgradC(k,l,i) = PgradC(k,l,i)-aux - PgradC(k,l,j) = PgradC(k,l,j)+aux + PgradC(k,l,i) = PgradC(k,l,i)+aux + PgradC(k,l,j) = PgradC(k,l,j)-aux enddo ! l enddo ! k + ENDIF #endif enddo ! j enddo ! iint @@ -11239,9 +11321,10 @@ c CA CA do k=1,nsaxs Cnorm = Cnorm + Pcalc(k) enddo - Esaxs_constr = dlog(Cnorm) + Esaxs_constr = dlog(Cnorm)-wsaxs0 do k=1,nsaxs - Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) + if (Pcalc(k).gt.0.0d0) + & Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) #ifdef DEBUG write (iout,*) "k",k," Esaxs_constr",Esaxs_constr #endif @@ -11256,7 +11339,8 @@ c CA CA auxX=0.0d0 auxX1=0.d0 do k=1,nsaxs - auxC = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k) + if (Pcalc(k).gt.0) + & auxC = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k) auxC1 = auxC1+PgradC(k,l,i) #ifdef ALLSAXS auxX = auxX +Psaxs(k)*PgradX(k,l,i)/Pcalc(k)