From: Cezary Czaplewski Date: Sat, 27 Jan 2018 21:41:47 +0000 (+0100) Subject: wham SAXS cutoff X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=a1e2f7505712e67c4a485fc9b56e8cdfdc25c84e wham SAXS cutoff --- diff --git a/source/wham/src-M/COMMON.CONTROL b/source/wham/src-M/COMMON.CONTROL index 9c2ac68..5a2e954 100644 --- a/source/wham/src-M/COMMON.CONTROL +++ b/source/wham/src-M/COMMON.CONTROL @@ -5,7 +5,8 @@ real*8 waga_homology real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut, & dist2_cut - real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs),scal_rad + real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs), + & scal_rad,wsaxs0,saxs_cutoff logical refstr,pdbref,punch_dist,print_rms,caonly,verbose, & merge_helices,bxfile,cxfile,histfile,entfile,zscfile, & rmsrgymap,with_dihed_constr,check_conf,histout,out1file, @@ -21,4 +22,6 @@ common /homol/ waga_homology(maxR), & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut, & iset,ihset,l_homo(max_template,maxdim) - common /saxsretr/ Psaxs,distsaxs,csaxs,scal_rad,nsaxs,saxs_mode + common /saxsretr/ Psaxs,distsaxs,csaxs,Wsaxs0,scal_rad, + & saxs_cutoff,nsaxs,saxs_mode + diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index fdc98e4..73ac0f7 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -8669,6 +8669,42 @@ C endif return 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' @@ -8704,6 +8740,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 @@ -8727,8 +8764,10 @@ c SAXS restraint penalty function enddo enddo do i=iatsc_s,iatsc_e + if (itype(i).eq.ntyp1) cycle do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + if (itype(j).eq.ntyp1) cycle #ifdef ALLSAXS dijCACA=dist(i,j) dijCASC=dist(i,j+nres) @@ -8796,22 +8835,43 @@ c SC SC #else dijCACA=dist(i,j) sigma2CACA=scal_rad**2*0.25d0/ - & (restok(itype(j))**2+restok(itype(i))**2) + & (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) + 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 @@ -8858,9 +8918,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 @@ -8875,7 +8936,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) diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F index 55c3fec..df062e5 100644 --- a/source/wham/src-M/molread_zs.F +++ b/source/wham/src-M/molread_zs.F @@ -361,6 +361,11 @@ c SAXS distance distribution do i=1,nsaxs write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i) enddo + Wsaxs0=0.0d0 + do i=1,nsaxs + Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i)) + enddo + write (iout,*) "Wsaxs0",Wsaxs0 else c SAXS "spheres". do i=1,nsaxs diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F index e2eca89..78f2521 100644 --- a/source/wham/src-M/readrtns.F +++ b/source/wham/src-M/readrtns.F @@ -129,8 +129,9 @@ C endif call readi(controlcard,'NSAXS',nsaxs,0) call readi(controlcard,'SAXS_MODE',saxs_mode,0) call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0) + call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0) write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE", - & SAXS_MODE," SCAL_RAD",scal_rad + & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff C /06/28/2013 Adasko: dyn_ss is keyword allowing to break and create bond C disulfide bond. Note that in conterary to dynamics this in C CONTROLCARD. The bond are read in molread_zs.F