if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
endif
-
+ if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
+ call e_saxs(Esaxs_constr)
+c write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr
+ else if (nsaxs.gt.0 .and. saxs_mode.gt.0) then
+ call e_saxsC(Esaxs_constr)
+c write (iout,*) "From EsaxsC: Esaxs_constr",Esaxs_constr
+ else
+ Esaxs_constr = 0.0d0
+ endif
c write(iout,*) "TEST_ENE",constr_homology
if (constr_homology.ge.1) then
call e_modeller(ehomology_constr)
c write(iout,*) "TEST_ENE",ehomology_constr
- write (iout,*) "ft(6)",fact(6),wliptran,eliptran
+c write (iout,*) "ft(6)",fact(6),wliptran,eliptran
#ifdef SPLITELE
if (shield_mode.gt.0) then
etot=fact(1)*wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
- & +wliptran*eliptran
+ & +wliptran*eliptran+wsaxs*esaxs_constr
else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
& +wvdwpp*evdw1
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ehomology_constr
- & +wliptran*eliptran
+ & +wliptran*eliptran+wsaxs*esaxs_constr
endif
#else
if (shield_mode.gt.0) then
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ehomology_constr
- & +wliptran*eliptran
+ & +wliptran*eliptran+wsaxs*esaxs_constr
else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
- & +wliptran*eliptran
+ & +wliptran*eliptran+wsaxs*esaxs_constr
endif
#endif
energia(20)=edihcnstr
energia(24)=ehomology_constr
energia(21)=evdw_t
+ energia(25)=Esaxs_constr
c energia(24)=ethetacnstr
energia(22)=eliptran
c detecting NaNQ
edihcnstr=energia(20)
estr=energia(18)
ehomology_constr=energia(24)
+ esaxs_constr=energia(25)
c ethetacnstr=energia(24)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
& ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
& eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
& eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
- & esccor,wsccor*fact(1),edihcnstr,ehomology_constr,ebr*nss,
- & etot
+ & esccor,wsccor*fact(1),edihcnstr,ehomology_constr,
+ & wsaxs*esaxs_constr,ebr*nss,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
+ & 'E_SAXS=',1pE16.6,' (SAXS restraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
& ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
& eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
& eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
- & edihcnstr,ehomology_constr,ebr*nss,
+ & edihcnstr,ehomology_constr,esaxs_constr*wsaxs,ebr*nss,
& etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
+ & 'E_SAXS=',1pE16.6,' (SAXS restraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
logical lprn
integer icant
external icant
+ integer xshift,yshift,zshift
evdw=0.0D0
evdw_t=0.0d0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
include 'COMMON.FFIELD'
include 'COMMON.SHIELD'
+ integer xshift,yshift,zshift
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
gcorr_loc(i)=0.0d0
enddo
do i=iatel_s,iatel_e
- if (i.le.1) cycle
+cAna if (i.le.1) cycle
if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
- & .or. ((i+2).gt.nres)
- & .or. ((i-1).le.0)
- & .or. itype(i+2).eq.ntyp1
- & .or. itype(i-1).eq.ntyp1
+cAna & .or. ((i+2).gt.nres)
+cAna & .or. ((i-1).le.0)
+cAna & .or. itype(i+2).eq.ntyp1
+cAna & .or. itype(i-1).eq.ntyp1
&) cycle
C endif
if (itel(i).eq.0) goto 1215
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
- if (j.le.1) cycle
+cAna if (j.le.1) cycle
if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
- & .or.((j+2).gt.nres)
- & .or.((j-1).le.0)
- & .or.itype(j+2).eq.ntyp1
- & .or.itype(j-1).eq.ntyp1
+cAna & .or.((j+2).gt.nres)
+cAna & .or.((j-1).le.0)
+cAna & .or.itype(j+2).eq.ntyp1
+cAna & .or.itype(j-1).eq.ntyp1
&) cycle
C endif
if (itel(j).eq.0) goto 1216
include 'COMMON.INTERACT'
include 'COMMON.FFIELD'
include 'COMMON.IOUNITS'
+ integer xshift,yshift,zshift
dimension ggg(3)
evdw2=0.0D0
evdw2_14=0.0d0
do i=idihconstr_start_homo,idihconstr_end_homo
kat2=0.0d0
c betai=beta(i,i+1,i+2,i+3)
- betai = phi(i+3)
+ betai = phi(i)
c write (iout,*) "betai =",betai
do k=1,constr_homology
dih_diff(k)=pinorm(dih(k,i)-betai)
c write (*,'(a,i2)') 'EBEND ICG=',icg
c write (iout,*) ithet_start,ithet_end
do i=ithet_start,ithet_end
- if (i.le.2) cycle
+C if (itype(i-1).eq.ntyp1) cycle
+c if (i.le.2) cycle
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
C Zero the energy function and its derivative at 0 or pi.
y(1)=0.0D0
y(2)=0.0D0
else
- if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
c icrc=0
etheta=0.0D0
c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
do i=ithet_start,ithet_end
- if (i.eq.2) cycle
+c if (i.eq.2) cycle
c print *,i,itype(i-1),itype(i),itype(i-2)
if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1)
& .or.(itype(i).eq.ntyp1)) cycle
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
enddo
else
phii=0.0d0
- ityp1=nthetyp+1
+ ityp1=ithetyp(itype(i-2))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp(itype(i))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
& )*wshield
C grad_shield_side is Cbeta sidechain gradient
grad_shield_side(j,ishield_list(i),i)=
- & (sh_frac_dist_grad(j)*-2.0d0
+ & (sh_frac_dist_grad(j)*(-2.0d0)
& *VofOverlap
& -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
&(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
&*VofOverlap
C grad_shield_side is Cbeta sidechain gradient
grad_shield_side(j,ishield_list(i),i)=
- & (sh_frac_dist_grad(j)*-2.0d0
+ & (sh_frac_dist_grad(j)*(-2.0d0)
& +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
& +scale_fac_dist*(cosphi_grad_long(j))
& *2.0d0/(1.0-cosphi))
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)
enddo
return
end
-C-------------------------------------------------------------------------------------
+c----------------------------------------------------------------------------
+ subroutine e_saxs(Esaxs_constr)
+ implicit none
+ include 'DIMENSIONS'
+#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(maxSAXS,3,maxres),
+ & PgradX(maxSAXS,3,maxres),Pcalc(maxSAXS)
+#ifdef MPI
+ double precision PgradC_(maxSAXS,3,maxres),
+ & PgradX_(maxSAXS,3,maxres),Pcalc_(maxSAXS)
+#endif
+ double precision dk,dijCACA,dijCASC,dijSCCA,dijSCSC,
+ & 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
+#ifdef DEBUG
+ write(iout,*) "------- SAXS penalty function start -------"
+ write (iout,*) "nsaxs",nsaxs
+ write (iout,*) "Esaxs: iatsc_s",iatsc_s," iatsc_e",iatsc_e
+ write (iout,*) "Psaxs"
+ do i=1,nsaxs
+ write (iout,'(i5,e15.5)') i, Psaxs(i)
+ enddo
+#endif
+ Esaxs_constr = 0.0d0
+ do k=1,nsaxs
+ Pcalc(k)=0.0d0
+ do j=1,nres
+ do l=1,3
+ PgradC(k,l,j)=0.0d0
+ PgradX(k,l,j)=0.0d0
+ enddo
+ 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)
+ dijSCCA=dist(i+nres,j)
+ dijSCSC=dist(i+nres,j+nres)
+ sigma2CACA=2.0d0/(pstok**2)
+ sigma2CASC=4.0d0/(pstok**2+restok(itype(j))**2)
+ sigma2SCCA=4.0d0/(pstok**2+restok(itype(i))**2)
+ sigma2SCSC=4.0d0/(restok(itype(j))**2+restok(itype(i))**2)
+ do k=1,nsaxs
+ dk = distsaxs(k)
+ expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+ if (itype(j).ne.10) then
+ expCASC = dexp(-0.5d0*sigma2CASC*(dijCASC-dk)**2)
+ else
+ endif
+ expCASC = 0.0d0
+ if (itype(i).ne.10) then
+ expSCCA = dexp(-0.5d0*sigma2SCCA*(dijSCCA-dk)**2)
+ else
+ expSCCA = 0.0d0
+ endif
+ if (itype(i).ne.10 .and. itype(j).ne.10) then
+ expSCSC = dexp(-0.5d0*sigma2SCSC*(dijSCSC-dk)**2)
+ else
+ expSCSC = 0.0d0
+ endif
+ Pcalc(k) = Pcalc(k)+expCACA+expCASC+expSCCA+expSCSC
+#ifdef DEBUG
+ write(iout,*) "i j k Pcalc",i,j,Pcalc(k)
+#endif
+ CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+ CASCgrad = sigma2CASC*(dijCASC-dk)*expCASC
+ SCCAgrad = sigma2SCCA*(dijSCCA-dk)*expSCCA
+ SCSCgrad = sigma2SCSC*(dijSCSC-dk)*expSCSC
+ 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
+c CA SC
+ if (itype(j).ne.10) then
+ aux = CASCgrad*(C(l,j+nres)-C(l,i))/dijCASC
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+ PgradX(k,l,j) = PgradX(k,l,j)+aux
+ endif
+c SC CA
+ if (itype(i).ne.10) then
+ aux = SCCAgrad*(C(l,j)-C(l,i+nres))/dijSCCA
+ PgradX(k,l,i) = PgradX(k,l,i)-aux
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+ endif
+c SC SC
+ if (itype(i).ne.10 .and. itype(j).ne.10) then
+ aux = SCSCgrad*(C(l,j+nres)-C(l,i+nres))/dijSCSC
+ PgradC(k,l,i) = PgradC(k,l,i)-aux
+ PgradC(k,l,j) = PgradC(k,l,j)+aux
+ PgradX(k,l,i) = PgradX(k,l,i)-aux
+ PgradX(k,l,j) = PgradX(k,l,j)+aux
+ endif
+ enddo ! l
+ enddo ! k
+#else
+ dijCACA=dist(i,j)
+ 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+
+ & 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'
+#ifdef MPI
+ include "mpif.h"
+ include "COMMON.SETUP"
+ integer IERR
+#endif
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ 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
projects / unres.git / blobdiff