& constr_homology,homol_nset,nsaxs,saxs_mode
real*8 waga_homology
real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut,
- & dist2_cut, scal_rad
+ & dist2_cut, scal_rad, saxs_cutoff
real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs)
logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
& sideadd,lsecondary,read_cart,unres_pdb,
& constr_homology,homol_nset,read2sigma,start_from_model
common /homol/ waga_homology(maxprocs/20),
& waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut
- common /saxsretr/ Psaxs,distsaxs,csaxs,scal_rad,nsaxs,saxs_mode
+ common /saxsretr/ Psaxs,distsaxs,csaxs,scal_rad,saxs_cutoff,
+ & nsaxs,saxs_mode
C... minim = .true. means DO minimization.
C... energy_dec = .true. means print energy decomposition matrix
sscalgrad2=0.0d0
else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then
gamm=(rr-(r_cut-rlamb))/rlamb
- sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb
- if (r.lt.r0) sscalgrad2=-sscalgrad2
+ 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
dijCACA=dist(i,j)
sigma2CACA=scal_rad**2*0.25d0/
& (restok(itype(j))**2+restok(itype(i))**2)
- rrr = 2.0d0/dsqrt(sigma2CACA)
+ 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.ne.0.0d0) then
ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0)
- if (energy_dec) write(iout,'(a4,3i5,4f10.4)')
- & 'saxs',i,j,k,dijCACA,rrr,dk,sss2
+ 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
endif ! sss
enddo ! k
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
call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
call readi(controlcard,'SYM',symetr,1)
call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours