end
c----------------------------------------------------------------------------
+ double precision function sscale2(r,r_cut,rlamb)
+ double precision r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscale2=1.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(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,rlamb)
+ double precision r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscalgrad2=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscalgrad2=0.0d0
+ endif
+ return
+ end
+c----------------------------------------------------------------------------
subroutine e_saxs(Esaxs_constr)
implicit none
include 'DIMENSIONS'
& 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
dijCACA=dist(i,j)
sigma2CACA=scal_rad**2*0.25d0/
& (restok(itype(j))**2+restok(itype(i))**2)
+ rrr = 2.0d0/dsqrt(sigma2CACA)
+ sss2 = sscale2(dijCACA,rrr,0.3d0)
+ if (sss2.ne.0.0d0) then
+ ssgrad2 = sscalgrad2(dijCACA,rrr,0.3d0)
do k=1,nsaxs
dk = distsaxs(k)
- expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+ 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,j) = PgradC(k,l,j)+aux
enddo ! l
enddo ! k
+ endif ! sss
#endif
enddo ! j
enddo ! iint