+ 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+
+ & 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
+ enddo ! l
+ enddo ! k
+ ENDIF
+#endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+#ifdef MPI
+ if (nfgtasks.gt.1) then
+ call MPI_Reduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION,
+ & MPI_SUM,king,FG_COMM,IERR)
+ if (fg_rank.eq.king) then
+ do k=1,nsaxs
+ Pcalc(k) = Pcalc_(k)
+ enddo
+ endif
+ call MPI_Reduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ if (fg_rank.eq.king) then
+ do i=1,nres
+ do l=1,3
+ do k=1,nsaxs
+ PgradC(k,l,i) = PgradC_(k,l,i)
+ enddo
+ enddo
+ enddo
+ endif
+#ifdef ALLSAXS
+ call MPI_Reduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ if (fg_rank.eq.king) then
+ do i=1,nres
+ do l=1,3
+ do k=1,nsaxs
+ PgradX(k,l,i) = PgradX_(k,l,i)
+ enddo
+ enddo
+ enddo
+ endif
+#endif
+ endif
+#endif
+#ifdef MPI
+ if (fg_rank.eq.king) then
+#endif
+ Cnorm = 0.0d0
+ do k=1,nsaxs
+ Cnorm = Cnorm + Pcalc(k)
+ enddo
+ Esaxs_constr = dlog(Cnorm)-wsaxs0
+ do k=1,nsaxs
+ 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
+ enddo
+#ifdef DEBUG
+ write (iout,*) "Cnorm",Cnorm," Esaxs_constr",Esaxs_constr
+#endif
+ do i=nnt,nct
+ do l=1,3
+ auxC=0.0d0
+ auxC1=0.0d0
+ auxX=0.0d0
+ auxX1=0.d0
+ do k=1,nsaxs
+ 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)
+ auxX1 = auxX1+PgradX(k,l,i)
+#endif
+ enddo
+ gsaxsC(l,i) = auxC - auxC1/Cnorm
+#ifdef ALLSAXS
+ gsaxsX(l,i) = auxX - auxX1/Cnorm
+#endif
+c write (iout,*) "l i",l,i," gradC",wsaxs*(auxC - auxC1/Cnorm),
+c * " gradX",wsaxs*(auxX - auxX1/Cnorm)
+ enddo
+ enddo
+#ifdef MPI
+ endif
+#endif
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine e_saxsC(Esaxs_constr)
+ implicit none
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+#ifdef MPI
+ include "mpif.h"
+ include "COMMON.SETUP"
+ integer IERR
+#endif
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.DERIV'
+ include 'COMMON.CONTROL'
+ include 'COMMON.NAMES'
+ include 'COMMON.FFIELD'
+ include 'COMMON.LANGEVIN'
+c
+ double precision Esaxs_constr
+ integer i,iint,j,k,l
+ double precision PgradC(3,maxres),PgradX(3,maxres),Pcalc,logPtot
+#ifdef MPI
+ double precision gsaxsc_(3,maxres),gsaxsx_(3,maxres),logPtot_
+#endif
+ double precision dk,dijCASPH,dijSCSPH,
+ & sigma2CA,sigma2SC,expCASPH,expSCSPH,
+ & CASPHgrad,SCSPHgrad,aux,auxC,auxC1,
+ & auxX,auxX1,Cnorm
+c SAXS restraint penalty function
+#ifdef DEBUG
+ write(iout,*) "------- SAXS penalty function start -------"
+ write (iout,*) "nsaxs",nsaxs," isaxs_start",isaxs_start,
+ & " isaxs_end",isaxs_end
+ write (iout,*) "nnt",nnt," ntc",nct
+ do i=nnt,nct
+ write(iout,'(a6,i5,3f10.5,5x,2f10.5)')
+ & "CA",i,(C(j,i),j=1,3),pstok,restok(itype(i))
+ enddo
+ do i=nnt,nct
+ write(iout,'(a6,i5,3f10.5)')"CSaxs",i,(Csaxs(j,i),j=1,3)
+ enddo
+#endif
+ Esaxs_constr = 0.0d0
+ logPtot=0.0d0
+ do j=isaxs_start,isaxs_end
+ Pcalc=0.0d0
+ do i=1,nres
+ do l=1,3
+ PgradC(l,i)=0.0d0
+ PgradX(l,i)=0.0d0
+ enddo
+ enddo
+ do i=nnt,nct
+ dijCASPH=0.0d0
+ dijSCSPH=0.0d0
+ do l=1,3
+ dijCASPH=dijCASPH+(C(l,i)-Csaxs(l,j))**2
+ enddo
+ if (itype(i).ne.10) then
+ do l=1,3
+ dijSCSPH=dijSCSPH+(C(l,i+nres)-Csaxs(l,j))**2
+ enddo
+ endif
+ sigma2CA=2.0d0/pstok**2
+ sigma2SC=4.0d0/restok(itype(i))**2
+ expCASPH = dexp(-0.5d0*sigma2CA*dijCASPH)
+ expSCSPH = dexp(-0.5d0*sigma2SC*dijSCSPH)
+ Pcalc = Pcalc+expCASPH+expSCSPH
+#ifdef DEBUG
+ write(*,*) "processor i j Pcalc",
+ & MyRank,i,j,dijCASPH,dijSCSPH, Pcalc
+#endif
+ CASPHgrad = sigma2CA*expCASPH
+ SCSPHgrad = sigma2SC*expSCSPH
+ do l=1,3
+ aux = (C(l,i+nres)-Csaxs(l,j))*SCSPHgrad
+ PgradX(l,i) = PgradX(l,i) + aux
+ PgradC(l,i) = PgradC(l,i)+(C(l,i)-Csaxs(l,j))*CASPHgrad+aux
+ enddo ! l
+ enddo ! i
+ do i=nnt,nct
+ do l=1,3
+ gsaxsc(l,i)=gsaxsc(l,i)+PgradC(l,i)/Pcalc
+ gsaxsx(l,i)=gsaxsx(l,i)+PgradX(l,i)/Pcalc
+ enddo
+ enddo
+ logPtot = logPtot - dlog(Pcalc)
+c print *,"me",me,MyRank," j",j," logPcalc",-dlog(Pcalc),
+c & " logPtot",logPtot
+ enddo ! j
+#ifdef MPI
+ if (nfgtasks.gt.1) then
+c write (iout,*) "logPtot before reduction",logPtot
+ call MPI_Reduce(logPtot,logPtot_,1,MPI_DOUBLE_PRECISION,
+ & MPI_SUM,king,FG_COMM,IERR)
+ logPtot = logPtot_
+c write (iout,*) "logPtot after reduction",logPtot
+ call MPI_Reduce(gsaxsC(1,1),gsaxsC_(1,1),3*nres,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ if (fg_rank.eq.king) then
+ do i=1,nres
+ do l=1,3
+ gsaxsC(l,i) = gsaxsC_(l,i)
+ enddo
+ enddo
+ endif
+ call MPI_Reduce(gsaxsX(1,1),gsaxsX_(1,1),3*nres,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+ if (fg_rank.eq.king) then
+ do i=1,nres
+ do l=1,3
+ gsaxsX(l,i) = gsaxsX_(l,i)
+ enddo
+ enddo
+ endif
+ endif
+#endif
+ Esaxs_constr = logPtot
+ return
+ end