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)) 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
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 & -(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)
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
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
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine Eliptransfer(eliptran)
implicit real*8 (a-h,o-z)
& 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
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)
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)
+ 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
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
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)
projects / unres.git / blobdiff