fi
elif [ "$1" == "1l2y_micro" ]; then
- refe="37.527"
+ refe="74.2623"
stat=`awk '{if ( $1 != "#" ) {n++;a=a+$5;a2=a2+$5^2}}END{print a/n,sqrt((a2-a^2/n)/n)}' $file_stat`
array=(${stat// / })
echo 'average total energy ' ${array[0]}
integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,symetr,
& constr_dist,shield_mode,tor_mode,
& constr_homology,homol_nset,
- & iset,ihset
+ & iset,ihset,nsaxs,saxs_mode
real*8 waga_homology
real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut,
& dist2_cut
+ real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs)
logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
& lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
& with_dihed_constr,with_theta_constr,out1file,
& waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut,
& iset,ihset,l_homo(max_template,maxdim),
& print_homology_restraints,print_homology_models
+ common /saxsretr/ Psaxs,distsaxs,csaxs,nsaxs,saxs_mode
& gshieldc, gshieldc_loc, gshieldx_ec, gshieldc_ec,
& gshieldc_loc_ec, gshieldx_t3,gshieldc_t3,gshieldc_loc_t3,
& gshieldx_t4, gshieldc_t4,gshieldc_loc_t4,gshieldx_ll,
- & gshieldc_ll, gshieldc_loc_ll
+ & gshieldc_ll, gshieldc_loc_ll,gsaxsC,gsaxsX,
+ & gliptranc,gliptranx
integer nfl,icg
& gshieldx_ll(3,-1:maxres), gshieldc_ll(3,-1:maxres),
& gshieldc_loc_ll(3,-1:maxres),
& gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres),
+ & gsaxsC(3,-1:maxres),gsaxsX(3,-1:maxres),
& gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres),
& gradcorr5(3,maxres),gradcorr6(3,maxres),
& gel_loc(3,maxres),gcorr3_turn(3,maxres),gcorr4_turn(3,maxres),
double precision wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
& wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,wturn6,
& wvdwpp,wbond,weights,scal14,scalscp,cutoff_corr,delt_corr,
- & r0_corr,wliptran
+ & r0_corr,wliptran,wsaxs
integer ipot,n_ene_comp,rescale_mode
common /ffield/ wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
& wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,wturn6,
- & wvdwpp,wbond,wliptran,weights(max_ene),scalscp,
+ & wvdwpp,wbond,wliptran,wsaxs,weights(max_ene),scalscp,
& scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp,
& rescale_mode
common /potentials/ potname(5)
c integer ithetaconstr_start_homo,ithetaconstr_end_homo
c
integer nresn,nyosh,nnos
- common /back_constr/ uconst_back,
+ common /back_constr/ uconst_back,uscdiff,
& dutheta,dugamma,duscdiff,duscdiffx
common /homrestr/ odl,dih,sigma_dih,sigma_odl,
& lim_odl,lim_dih,ires_homo,jres_homo,link_start_homo,
--- /dev/null
+ double precision scal_fric,rwat,etawat,gamp,
+ & gamsc(ntyp1),stdfp,stdfsc(ntyp),stdforcp(MAXRES),
+ & stdforcsc(MAXRES),pstok,restok(ntyp+1),cPoise,Rb
+ common /langevin/ pstok,restok,gamp,gamsc,
+ & stdfp,stdfsc,stdforcp,stdforcsc,rwat,etawat,cPoise,Rb
+ double precision IP,ISC(ntyp+1),mp,
+ & msc(ntyp+1)
+ common /inertia/ IP,ISC,MP,MSC
& sigc0,dsc,dsc_inv,bsc,censc,gaussc,dsc0,vbl,vblinv,vblinv2,
& vbl_cis,vbl0,vbld_inv
integer nlob,loc_start,loc_end,ithet_start,ithet_end,
- & iphi_start,iphi_end,itau_start,itau_end
+ & iphi_start,iphi_end,itau_start,itau_end,isaxs_start,isaxs_end
C Parameters of the virtual-bond-angle probability distribution
common /thetas/ a0thet(-ntyp:ntyp),athet(2,-ntyp:ntyp,-1:1,-1:1)
& ,bthet(2,-ntyp:ntyp,-1:1,-1:1),
C Virtual-bond lenghts
common /peptbond/ vbl,vblinv,vblinv2,vbl_cis,vbl0
common /indices/ loc_start,loc_end,ithet_start,ithet_end,
- & iphi_start,iphi_end,itau_start,itau_end
+ & iphi_start,iphi_end,itau_start,itau_end,isaxs_start,isaxs_end
C Inverses of the actual virtual bond lengths
common /invlen/ vbld_inv(maxres2)
C Maximum number of templates in homology-modeling restraints
integer max_template
parameter(max_template=25)
+C Maximum number of bins in SAXS restraints
+ integer MaxSAXS
+ parameter (MaxSAXS=1000)
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)
& +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),
include 'COMMON.INTERACT'
include 'COMMON.FFIELD'
include 'COMMON.IOUNITS'
+ integer xshift,yshift,zshift
dimension ggg(3)
evdw2=0.0D0
evdw2_14=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))
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 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
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+#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=0.25d0/(restok(itype(j))**2+restok(itype(i))**2)
+ do k=1,nsaxs
+ dk = distsaxs(k)
+ expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+ Pcalc(k) = Pcalc(k)+expCACA
+#ifdef DEBUG
+ write(iout,*) "i j k Pcalc",i,j,Pcalc(k)
+#endif
+ CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+ 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
+ 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)
+ do k=1,nsaxs
+ 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
+ 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
& "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
& "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
& "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
- & "EAFM","ETHETC",""/
+ & "EAFM","ETHETC","ESAXS"/
data wname /
& "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
& "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
& "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
- & "WLIPTRAN","WAFM","WTHETC",""/
- data nprint_ene /22/
+ & "WLIPTRAN","WAFM","WTHETC","WSAXS"/
+ data nprint_ene /23/
data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
- & 16,15,17,20,21,24,22,23,0/
+ & 16,15,17,20,21,24,22,23,25/
end
c---------------------------------------------------------------------------
subroutine init_int_table
#ifdef MPL
include 'COMMON.INFO'
#endif
+ include 'COMMON.CONTROL'
include 'COMMON.CHAIN'
include 'COMMON.INTERACT'
include 'COMMON.LOCAL'
call int_bounds(nres-3,itau_start,itau_end)
itau_start=itau_start+3
itau_end=itau_end+3
+ call int_bounds(nsaxs,isaxs_start,isaxs_end)
if (lprint) then
write (iout,*) 'Processor:',MyID,
& ' loc_start',loc_start,' loc_end',loc_end,
iphi_end=nct
itau_start=4
itau_end=nres
+ isaxs_start=1
+ isaxs_end=nsaxs
+ write (iout,*) "OSAXS_START",isaxs_start," ISAXS_END",isaxs_end
#endif
return
end
include 'COMMON.SBRIDGE'
include 'COMMON.SCCOR'
include 'COMMON.SCROT'
+ include 'COMMON.LANGEVIN'
character*1 t1,t2,t3
character*1 onelett(-2:2) /"p","a","G","A","P"/
logical lprint
dimension blower(3,3,maxlob)
- double precision ip,mp
C
C Body
C
vblinv=1.0D0/vbl
vblinv2=vblinv*vblinv
#ifdef CRYST_BOND
- read (ibond,*) vbldp0,vbldpdum,akp
+ read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
do i=1,ntyp
nbondterm(i)=1
- read (ibond,*) vbldsc0(1,i),aksc(1,i)
+ read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
dsc(i) = vbldsc0(1,i)
if (i.eq.10) then
dsc_inv(i)=0.0D0
endif
enddo
#else
- read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
+ read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
do i=1,ntyp
read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
- & j=1,nbondterm(i))
+ & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
dsc(i) = vbldsc0(1,i)
if (i.eq.10) then
dsc_inv(i)=0.0D0
double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc,
& ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
& eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
- & evdw_t
+ & evdw_t,esaxs
integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
double precision qfree,sumprob,eini,efree,rmsdev
character*80 bxname
estr=enetb(18,i)
esccor=enetb(19,i)
edihcnstr=enetb(20,i)
+ esaxs=enetb(25,i)
c#ifdef SPLITELE
c etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
c &ft(1)*welec*ees+wvdwpp*evdw1
c enddo
write (iout,*) "nlist",nlist
#endif
+ write (iout,'(80(1h-)/i4,a,f6.3,a,f5.1/80(1h-))')
+ & nlist," conformations found within",prob_limit,
+ & " probablity cut_off at temperature ",1.0d0/(1.987d-3*beta_h(ib))
return
end
!--------------------------------------------------
include 'COMMON.CHAIN'
include 'COMMON.HEADER'
include 'COMMON.FFIELD'
- include 'COMMON.FREE'
include 'COMMON.INTERACT'
include "COMMON.SPLITELE"
include 'COMMON.SHIELD'
+ include 'COMMON.FREE'
+ include 'COMMON.LANGEVIN'
character*320 controlcard,ucase
#ifdef MPL
include 'COMMON.INFO'
print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
print_homology_models=
& index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
+ call readi(controlcard,'NSAXS',nsaxs,0)
+ call readi(controlcard,'SAXS_MODE',saxs_mode,0)
+ write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+ & SAXS_MODE
if (min_var) iopt=1
return
end
call reada(weightcard,'WTORD',wtor_d,1.0D0)
call reada(weightcard,'WANG',wang,1.0D0)
call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+ call reada(weightcard,'WSAXS',wsaxs,0.0D0)
call reada(weightcard,'SCAL14',scal14,0.4D0)
call reada(weightcard,'SCALSCP',scalscp,1.0d0)
call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
if (nstart.lt.nnt) nstart=nnt
if (nend.gt.nct .or. nend.eq.0) nend=nct
write (iout,*) "nstart",nstart," nend",nend
+ write (iout,*) "calling read_saxs_consrtr",nsaxs
+ if (nsaxs.gt.0) call read_saxs_constr
nres0=nres
if (constr_homology.gt.0) then
call read_constr_homology(print_homology_restraints)
return
end
+c-------------------------------------------------------------------------------
+ subroutine read_saxs_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ double precision cm(3)
+c read(inp,*) nsaxs
+ write (iout,*) "Calling read_saxs nsaxs",nsaxs
+ call flush(iout)
+ if (saxs_mode.eq.0) then
+c SAXS distance distribution
+ do i=1,nsaxs
+ read(inp,*) distsaxs(i),Psaxs(i)
+ enddo
+ Cnorm = 0.0d0
+ do i=1,nsaxs
+ Cnorm = Cnorm + Psaxs(i)
+ enddo
+ write (iout,*) "Cnorm",Cnorm
+ do i=1,nsaxs
+ Psaxs(i)=Psaxs(i)/Cnorm
+ enddo
+ write (iout,*) "Normalized distance distribution from SAXS"
+ do i=1,nsaxs
+ write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
+ enddo
+ else
+c SAXS "spheres".
+ do i=1,nsaxs
+ read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
+ enddo
+ do j=1,3
+ cm(j)=0.0d0
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ cm(j)=cm(j)+Csaxs(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/nsaxs
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ Csaxs(j,i)=Csaxs(j,i)-cm(j)
+ enddo
+ enddo
+ write (iout,*) "SAXS sphere coordinates"
+ do i=1,nsaxs
+ write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
+ enddo
+ endif
+ return
+ end
c====-------------------------------------------------------------------
subroutine read_constr_homology(lprn)
& buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
& buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+ double precision forceAFMconst, distafminit
common /afm/ forceAFMconst, distafminit,afmend,afmbeg,
& velAFMconst,
& totTafm
integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
- & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
- & constr_homology,homol_nset
+ & iprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
+ & 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
+ real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs)
logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
& sideadd,lsecondary,read_cart,unres_pdb,
& vdisulf,searchsc,lmuca,dccart,extconf,out1file,
& 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,nsaxs,saxs_mode
C... minim = .true. means DO minimization.
C... energy_dec = .true. means print energy decomposition matrix
- double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gelc_long
+ double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gelc_long,
& gvdwpp,gel_loc,gel_loc_long,gvdwc_scpp,gliptranc,gliptranx,
& gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gloc_x,dtheta,dphi,dalpha,
& domega,gscloc,gsclocx,gradcorr,gradcorr_long,gradcorr5_long,
- & gradcorr6_long,gcorr6_turn_long,gvdwx
+ & gradcorr6_long,gcorr6_turn_long,gvdwx,gradafm,gradxorr,gradcorr5,
+ &gradcorr6,gcorr3_turn,gcorr4_turn,gcorr6,gcorr6_turn,gradb,gradbx,
+ & gel_loc_loc,gel_loc_turn3,gel_loc_turn4,gel_loc_turn6,gcorr_loc,
+ & g_corr5_loc,g_corr6_loc,gsccorc,gsccorx,gsccor_loc,gsaxsC,gsaxsX
integer nfl,icg
common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres),
& gradx(3,-1:maxres,2),gradc(3,-1:maxres,2),gvdwx(3,-1:maxres),
& gliptranx(3,-1:maxres),
& gradafm(3,-1:maxres),
& gradx_scp(3,-1:maxres),gvdwc_scp(3,-1:maxres),
- & ghpbx(3,-1:maxres),
- & ghpbc(3,-1:maxres),gloc(maxvar,2),gradcorr(3,-1:maxres),
+ & ghpbx(3,-1:maxres),ghpbc(3,-1:maxres),
+ & gsaxsC(3,-1:maxres),gsaxsX(3,-1:maxres),
+ & gloc(maxvar,2),gradcorr(3,-1:maxres),
& gradcorr_long(3,-1:maxres),gradcorr5_long(3,-1:maxres),
& gradcorr6_long(3,-1:maxres),gcorr6_turn_long(3,-1:maxres),
& gradxorr(3,-1:maxres),gradcorr5(3,-1:maxres),
C calculations and defines weights of various energy terms.
C 12/1/95 wcorr added
C-----------------------------------------------------------------------
- integer n_ene_comp,rescale_mode
+ integer n_ene_comp,rescale_mode,ipot
+ double precision wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,
+ & wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
+ & wturn6,wvdwpp,weights,wliptran,wsaxs,temp0,
+ & scal14,cutoff_corr,delt_corr,r0_corr
common /ffield/ wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,
& wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
- & wturn6,wvdwpp,weights(n_ene),wliptran,temp0,
+ & wturn6,wvdwpp,weights(n_ene),wliptran,wsaxs,temp0,
& scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp,
& rescale_mode
common /potentials/ potname(5)
& ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
& iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
& iint_end,iphi1_start,iphi1_end,itau_start,itau_end,ilip_start,
- & ilip_end,
+ & ilip_end,isaxs_start,isaxs_end,
& ibond_displ(0:max_fg_procs-1),ibond_count(0:max_fg_procs-1),
& ithet_displ(0:max_fg_procs-1),ithet_count(0:max_fg_procs-1),
& iphi_displ(0:max_fg_procs-1),iphi_count(0:max_fg_procs-1),
& iint_end,iphi1_start,iphi1_end,iint_count,iint_displ,ivec_displ,
& ivec_count,iset_displ,itau_start,itau_end,
& iset_count,ibond_displ,ibond_count,ithet_displ,ithet_count,
- & iphi_displ,iphi_count,iphi1_displ,iphi1_count,ilip_start,ilip_end
+ & iphi_displ,iphi_count,iphi1_displ,iphi1_count,
+ & ilip_start,ilip_end,isaxs_start,isaxs_end
C Inverses of the actual virtual bond lengths
common /invlen/ vbld_inv(maxres2)
common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,
& ns,nss,nfree,iss(maxss)
double precision dhpb,dhpb1,forcon
- integer ihpb,jhpb,nhpb,idssb,jdssb
+ integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb
common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
& ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb
double precision weidis
parameter (max_pool=10)
C Number of energy components
integer n_ene,n_ene2
- parameter (n_ene=24,n_ene2=2*n_ene)
+ parameter (n_ene=25,n_ene2=2*n_ene)
C Number of threads in deformation
integer max_thread,max_thread2
parameter (max_thread=4,max_thread2=2*max_thread)
C Maximum number of templates in homology-modeling restraints
integer max_template
parameter(max_template=25)
+C Maximum number of bins in SAXS restraints
+ integer MaxSAXS
+ parameter (MaxSAXS=1000)
common /gucio/ cm
integer itime
logical ovrtim
+ integer nharp,iharp(4,maxres/3)
c
#ifdef MPI
if (ilen(tmpdir).gt.0)
t_enegrad=0.0d0
t_sdsetup=0.0d0
write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
+ write (iout,'(/a)')
+ & "Cartesian coordinates of the initial structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+ & (c(j,ires+nres),j=1,3)
+ enddo
+ write (iout,'(/a)')
+ & "Initial dC vectors of the chain"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(dc(j,ires),j=1,3),
+ & (dc(j,ires+nres),j=1,3)
+ enddo
+ write (iout,'(/a)')
+ & "Initial dC_norm vectors of the chain"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(dc_norm(j,ires),j=1,3),
+ & (dc_norm(j,ires+nres),j=1,3)
+ enddo
#ifdef MPI
tt0=MPI_Wtime()
#else
include 'COMMON.NAMES'
include 'COMMON.REMD'
real*8 energia_long(0:n_ene),
- & energia_short(0:n_ene),vcm(3),incr(3)
+ & energia_short(0:n_ene),energia(0:n_ene),vcm(3),incr(3)
double precision cm(3),L(3),xv,sigv,lowb,highb
double precision varia(maxvar)
character*256 qstr
dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
enddo
enddo
+ call etotal(energia(0))
endif
c call chainbuild_cart
+ write (iout,*) "Initial energies"
+ call enerprint(energia(0))
write (iout,*) "PREMINIM ",preminim
- if(iranconf.eq.0 .and. preminim) then
+ if(iranconf.ne.0 .or. preminim) then
if (overlapsc) then
write (iout,*) 'Calling OVERLAP_SC'
call overlap_sc(fail)
C Loop over i,i+2 and i,i+3 pairs of the peptide groups
C
do i=iturn3_start,iturn3_end
- if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
- & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
- & .or. itype(i-1).eq.ntyp1
- & .or. itype(i+4).eq.ntyp1
- & ) cycle
+c AL 7/8/16 CHUJ DUPA I KAMIENI KUPA. PRZECIEZ TO BYLO KURWA MAC W INNYCH
+C WERSJACH DAWNO DO CHUJA JEBANEGO POPRAWIONE!!!
+c Wylaczamy oddzialywanie 1-3 tylko wtedy gdy ktoras grupa peptydowa
+c jest dummy a to oznacza, ze reszta i lub i+1 lub i+2 lub i+3 jest dummy
+c reszta i0-1 do tego nie nalezy!
+c if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+c & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
+c & .or. itype(i-1).eq.ntyp1
+c & .or. itype(i+4).eq.ntyp1
+c & ) cycle
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+ & .or. itype(i+2).eq.ntyp1
+ & .or. itype(i+3).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
+c JAK WYZEJ!!!
+c if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+c & .or. itype(i+3).eq.ntyp1
+c & .or. itype(i+4).eq.ntyp1
+c & .or. itype(i+5).eq.ntyp1
+c & .or. itype(i-1).eq.ntyp1
+c & ) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+3).eq.ntyp1
& .or. itype(i+4).eq.ntyp1
- & .or. itype(i+5).eq.ntyp1
- & .or. itype(i-1).eq.ntyp1
- & ) cycle
+ & ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
- & .or. itype(i+2).eq.ntyp1
- & .or. itype(i-1).eq.ntyp1
- &) cycle
+C PRZECIEZ TU ODDZIALUUJA GRUPY PEPTYDOWE MIEDZY RESZTAMI I I+1 oraz j j+1
+c po co sprawdzac typy reszt i-1 oraz i-2?
+c if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+c & .or. itype(i+2).eq.ntyp1
+c & .or. itype(i-1).eq.ntyp1
+c &) cycle
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
- & .or.itype(j+2).eq.ntyp1
- & .or.itype(j-1).eq.ntyp1
- &) cycle
+c if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+c & .or.itype(j+2).eq.ntyp1
+c & .or.itype(j-1).eq.ntyp1
+cc &) cycle
+ if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
call eelecij_scale(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
weights_(17)=wbond
weights_(18)=scal14
weights_(21)=wsccor
+ weights_(25)=wsaxs
C FG Master broadcasts the WEIGHTS_ array
call MPI_Bcast(weights_(1),n_ene,
& MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
wbond=weights(17)
scal14=weights(18)
wsccor=weights(21)
+ wsaxs=weights(25)
endif
time_Bcast=time_Bcast+MPI_Wtime()-time00
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
cd write (iout,*) "multibody_hb ecorr",ecorr
endif
+c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode
+ 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 print *,"Processor",myrank," computed Ucorr"
C
C If performing constraint dynamics, call the constraint energy
energia(22)=eliptran
energia(23)=Eafmforce
energia(24)=ehomology_constr
+ energia(25)=Esaxs_constr
c Here are the energies showed per procesor if the are more processors
c per molecule then we sum it up in sum_energy subroutine
c print *," Processor",myrank," calls SUM_ENERGY"
eliptran=energia(22)
Eafmforce=energia(23)
ehomology_constr=energia(24)
+ esaxs_constr=energia(25)
+c write (iout,*) "sum_energy esaxs_constr",esaxs_constr,
+c & " wsaxs",wsaxs
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
- & +wliptran*eliptran+Eafmforce
+ & +wsaxs*esaxs_constr+wliptran*eliptran+Eafmforce
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
- & +wliptran*eliptran
+ & +wsaxs*esaxs_constr+wliptran*eliptran
& +Eafmforce
#endif
energia(0)=etot
#endif
#ifdef DEBUG
write (iout,*) "sum_gradient gvdwc, gvdwx"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)')
+ do i=0,nres
+ write (iout,'(i3,3e15.5,5x,3e15.5)')
& i,(gvdwx(j,i),j=1,3),(gvdwc(j,i),j=1,3)
enddo
call flush(iout)
#endif
+#ifdef DEBUG
+ write (iout,*) "sum_gradient gsaxsc, gsaxsx"
+ do i=0,nres
+ write (iout,'(i3,3e15.5,5x,3e15.5)')
+ & i,(gsaxsc(j,i),j=1,3),(gsaxsx(j,i),j=1,3)
+ enddo
+ call flush(iout)
+#endif
#ifdef MPI
C FG slaves call the following matching MPI_Bcast in ERGASTULUM
if (nfgtasks.gt.1 .and. fg_rank.eq.0)
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wsaxs*gsaxsc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wsaxs*gsaxsc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,-1,-1
+ do i=nres-2,0,-1
+c do i=nres-2,nnt,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
#ifdef DEBUG
write (iout,*) "gradbufc after summing"
do i=1,nres
- write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+ write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3)
enddo
call flush(iout)
#endif
#endif
#ifdef DEBUG
write (iout,*) "gradbufc"
- do i=1,nres
- write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+ do i=0,nres
+ write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3)
enddo
call flush(iout)
#endif
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,-1,-1
+ do i=nres-2,0,-1
+c do i=nres-2,nnt,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
c enddo
#ifdef DEBUG
write (iout,*) "gradbufc after summing"
- do i=1,nres
- write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+ do i=0,nres
+ write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3)
enddo
call flush(iout)
#endif
#endif
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
- & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
- & wsccor*gsccorx(j,i)
+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)
+ & +wsaxs*gsaxsx(j,i)
+ & +wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
& +wliptran*gliptranx(j,i)
enddo
#ifdef DEBUG
write (iout,*) "gradc gradx gloc"
do i=1,nres
- write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)')
+ write (iout,'(i5,3e15.5,5x,3e15.5,5x,e15.5)')
& i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg)
enddo
#endif
Uconst=energia(20)
esccor=energia(21)
ehomology_constr=energia(24)
+ esaxs_constr=energia(25)
eliptran=energia(22)
Eafmforce=energia(23)
#ifdef SPLITELE
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
- & edihcnstr,ehomology_constr, ebr*nss,
+ & edihcnstr,ehomology_constr,esaxs_constr*wsaxs, ebr*nss,
& Uconst,eliptran,wliptran,Eafmforce,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)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ehomology_constr,ebr*nss,Uconst,
+ & ehomology_constr,esaxs_constr*wsaxs,ebr*nss,Uconst,
& eliptran,wliptran,Eafmforc,
& etot
10 format (/'Virtual-chain energies:'//
& '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)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
return
end
+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.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.CONTROL'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ include 'COMMON.FFIELD'
+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 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=0.25d0/(restok(itype(j))**2+restok(itype(i))**2)
+ do k=1,nsaxs
+ dk = distsaxs(k)
+ expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+ Pcalc(k) = Pcalc(k)+expCACA
+#ifdef DEBUG
+ write(iout,*) "i j k Pcalc",i,j,Pcalc(k)
+#endif
+ CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+ 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
+ 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)
+ do k=1,nsaxs
+ 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
+ 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.GEO'
+ include 'COMMON.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.CONTROL'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ include 'COMMON.FFIELD'
+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
+
+ do i=nnt,nct
+ print *,MyRank,"C",i,(C(j,i),j=1,3)
+ enddo
+ do i=nnt,nct
+ print *,MyRank,"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
+ if (itype(i).eq.ntyp1) cycle
+ 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
double precision weights_(n_ene)
#endif
include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
include 'COMMON.IOUNITS'
double precision energia(0:n_ene)
include 'COMMON.FFIELD'
weights_(17)=wbond
weights_(18)=scal14
weights_(21)=wsccor
+ weights_(25)=wsaxs
C FG Master broadcasts the WEIGHTS_ array
call MPI_Bcast(weights_(1),n_ene,
& MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
wbond=weights(17)
scal14=weights(18)
wsccor=weights(21)
+ wsaxs=weights(25)
endif
call MPI_Bcast(dc(1,1),6*nres,MPI_DOUBLE_PRECISION,
& king,FG_COMM,IERR)
double precision weights_(n_ene)
#endif
include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
include 'COMMON.IOUNITS'
double precision energia(0:n_ene)
include 'COMMON.FFIELD'
weights_(17)=wbond
weights_(18)=scal14
weights_(21)=wsccor
+ weights_(25)=wsaxs
C FG Master broadcasts the WEIGHTS_ array
call MPI_Bcast(weights_(1),n_ene,
& MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
wbond=weights(17)
scal14=weights(18)
wsccor=weights(21)
+ wsaxs=weights(25)
endif
c write (iout,*),"Processor",myrank," BROADCAST weights"
call MPI_Bcast(c(1,1),maxres6,MPI_DOUBLE_PRECISION,
else
esccor=0.0d0
endif
+
+c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode
+ 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
C Put energy components into an array
C
energia(17)=estr
energia(19)=edihcnstr
energia(21)=esccor
+ energia(25)=Esaxs_constr
c write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
call flush(iout)
call sum_energy(energia,.true.)
do ires=1,ioverlap_last
i=ioverlap(ires)
iti=iabs(itype(i))
- if (iti.ne.10) then
+ if (iti.ne.10 .and. iti.lt.ntyp1 .and.iti.gt.0) then
nsi=0
fail=.true.
do while (fail.and.nsi.le.maxsi)
do i=iatsc_s,iatsc_e
itypi=iabs(itype(i))
itypi1=iabs(itype(i+1))
+ if (itypi.gt.ntyp .or. itypi.le.0) cycle
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
itypj=iabs(itype(j))
+ if (itypj.gt.ntyp .or. itypj.le.0) cycle
dscj_inv=dsc_inv(itypj)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
include 'COMMON.IOUNITS'
include 'COMMON.HEADER'
include 'COMMON.SBRIDGE'
- character*32 tytul,fd
+ character*50 tytul,fd
character*3 zahl
character*6 res_num,pom,ucase
#ifdef AIX
line2=' '
endif
if (print_compon) then
+#ifdef DEBUG
+ write (iout,*) "itime",itime," temperature",t_bath,
+ & " potential energy",potE,potEcomp(0)
+ call enerprint(potEcomp)
+#endif
if(itime.eq.0) then
write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
& ",100a12)"
time00=MPI_Wtime()
#endif
icg=1
+#ifdef DEBUG
+ write (iout,*) "Before sum_gradient"
+ do i=1,nres-1
+ write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
+ write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
+ enddo
+#endif
call sum_gradient
#ifdef TIMING
#endif
gel_loc_long(j,i)=0.0d0
ghpbc(j,i)=0.0D0
ghpbx(j,i)=0.0D0
+ gsaxsc(j,i)=0.0D0
+ gsaxsx(j,i)=0.0D0
gcorr3_turn(j,i)=0.0d0
gcorr4_turn(j,i)=0.0d0
gradcorr(j,i)=0.0d0
& "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
& "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
& "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB ","EVDWPP ",
- & "ESTR ","EVDW2_14 ","UCONST ", " ","ESCCOR",
- & "Eliptran","Eafmforce","Ehomology"/
+ & "ESTR ","EVDW2_14 ","UCONST ", "EDIHC ","ESCCOR",
+ & "Eliptran","Eafmforce","Ehomology","ESAXS"/
data wname /
& "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
& "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
- & "WSTRAIN","WVDWPP","WBOND","SCAL14"," "," ","WSCCOR",
- & "Wliptran"," ","EHOMO"/
- data nprint_ene /21/
+ & "WSTRAIN","WVDWPP","WBOND","SCAL14","WDIHCSN"," ","WSCCOR",
+ & "Wliptran","WAFM ","EHOMO","WSAXS"/
+ data nprint_ene /25/
data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
- & 21,24,22,23,0/
+ & 21,24,22,23,20,25/
end
c---------------------------------------------------------------------------
subroutine init_int_table
call int_bounds(nres-2,ithet_start,ithet_end)
ithet_start=ithet_start+2
ithet_end=ithet_end+2
+ call int_bounds(nsaxs,isaxs_start,isaxs_end)
+c isaxs_start=isaxs_start+nnt-1
+c isaxs_end=isaxs_end+nnt-1
+ write (iout,*) me," isaxs_start",isaxs_start,
+ & " isaxs_end",isaxs_end
call int_bounds(nct-nnt-2,iturn3_start,iturn3_end)
iturn3_start=iturn3_start+nnt
iphi_start=iturn3_start+2
iaux=iphi1_end-iphi1_start+1
call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
& MPI_INTEGER,FG_COMM,IERROR)
- do i=0,maxprocs-1
+ write (iout,*) "ielstart before zeroing out",max_fg_procs
+ call flush (iout)
+ do i=0,max_fg_procs-1
do j=1,maxres
ielstart_all(j,i)=0
ielend_all(j,i)=0
enddo
enddo
+ write (iout,*) "ielstart zeroed out"
+ call flush (iout)
call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
& iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
iphi1_end=nres
idihconstr_start=1
idihconstr_end=ndih_constr
+ isaxs_start=1
+ isaxs_end=nsaxs
iphid_start=iphi_start
iphid_end=iphi_end-1
itau_start=4
C
if (lprint) then
write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
+ do iblock=1,2
do i=1,nthetyp+1
do j=1,nthetyp+1
do k=1,nthetyp+1
enddo
enddo
enddo
+ enddo
call flush(iout)
endif
c write (2,*) "Start reading THETA_PDB",ithep_pdb
close (itorp)
if (lprint) then
write (iout,'(/a/)') 'Torsional constants:'
+ do iblock=1,2
do i=1,ntortyp
do j=1,ntortyp
write (iout,*) 'ityp',i,' jtyp',j
enddo
enddo
enddo
+ enddo
endif
C
v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
&(1+vlor3sccor(k,i,j)**2)
enddo
- v0sccor(i,j,iblock)=v0ijsccor
+ v0sccor(l,i,j)=v0ijsccor
enddo
enddo
enddo
#endif
if (lprint) then
write (iout,'(/a/)') 'Torsional constants:'
+ do l=1,maxinter
do i=1,nsccortyp
do j=1,nsccortyp
write (iout,*) 'ityp',i,' jtyp',j
enddo
enddo
enddo
+ enddo
endif
C
C Set up the time limit (caution! The time must be input in minutes!)
read_cart=index(controlcard,'READ_CART').gt.0
call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ write (iout,*) "constr_dist",constr_dist
+ call readi(controlcard,'NSAXS',nsaxs,0)
+ call readi(controlcard,'SAXS_MODE',saxs_mode,0)
+ write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+ & SAXS_MODE
call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
call readi(controlcard,'SYM',symetr,1)
call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
call reada(weightcard,'TEMP0',temp0,300.0d0)
call reada(weightcard,'WLT',wliptran,0.0D0)
+ call reada(weightcard,'WSAXS',wsaxs,1.0D0)
if (index(weightcard,'SOFT').gt.0) ipot=6
C 12/1/95 Added weight for the multi-body term WCORR
call reada(weightcard,'WCORRH',wcorr,1.0D0)
weights(17)=wbond
weights(18)=scal14
weights(21)=wsccor
+ weights(25)=wsaxs
if(me.eq.king.or..not.out1file)
& write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
nsi=nsi+1
enddo
+c AL 7/10/16
+c Calculalte only the coordinates of the current sidechain; no need to rebuild
+c whole chain
+ call locate_side_chain(i)
if(fail) write(iout,*)'Adding sidechain failed for res ',
& i,' after ',nsi,' trials'
endif
enddo
C 10/03/12 Adam: Recalculate coordinates with new side chain positions
- call chainbuild
+c call chainbuild
endif
endif
if (indpdb.eq.0) then
phibound(2,i)=pi
enddo
read (inp,*) ndih_constr
+ write (iout,*) "ndish_constr",ndih_constr
if (ndih_constr.gt.0) then
read (inp,*) ftors
read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
enddo
endif
endif
+ write (iout,*) "calling read_saxs_consrtr",nsaxs
+ if (nsaxs.gt.0) call read_saxs_constr
if (constr_homology.gt.0) then
end
c-------------------------------------------------------------------------------
+ subroutine read_saxs_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ double precision cm(3)
+c read(inp,*) nsaxs
+ write (iout,*) "Calling read_saxs nsaxs",nsaxs
+ call flush(iout)
+ if (saxs_mode.eq.0) then
+c SAXS distance distribution
+ do i=1,nsaxs
+ read(inp,*) distsaxs(i),Psaxs(i)
+ enddo
+ Cnorm = 0.0d0
+ do i=1,nsaxs
+ Cnorm = Cnorm + Psaxs(i)
+ enddo
+ write (iout,*) "Cnorm",Cnorm
+ do i=1,nsaxs
+ Psaxs(i)=Psaxs(i)/Cnorm
+ enddo
+ write (iout,*) "Normalized distance distribution from SAXS"
+ do i=1,nsaxs
+ write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
+ enddo
+ else
+c SAXS "spheres".
+ do i=1,nsaxs
+ read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
+ enddo
+ do j=1,3
+ cm(j)=0.0d0
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ cm(j)=cm(j)+Csaxs(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/nsaxs
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ Csaxs(j,i)=Csaxs(j,i)-cm(j)
+ enddo
+ enddo
+ write (iout,*) "SAXS sphere coordinates"
+ do i=1,nsaxs
+ write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
+ enddo
+ endif
+ return
+ end
+c-------------------------------------------------------------------------------
subroutine read_dist_constr
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
itypi=iabs(itype(i))
+ if (itypi.gt.ntyp .or. itypi.le.0) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
ind=ind+1
itypj=iabs(itype(j))
+ if (itypj.gt.ntyp .or. itypj.le.0) cycle
dscj_inv=dsc_inv(itypj)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
& allihpb(maxdim),alljhpb(maxdim),
& newnss,newihpb(maxdim),newjhpb(maxdim)
logical found
- integer i_newnss(max_fg_procs),displ(max_fg_procs)
+ integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
allnss=0
integer iscode,indpdb,outpdb,outmol2,icomparfunc,pdbint,
& ensembles,constr_dist,symetr,
& constr_homology,homol_nset,
- & iset,ihset
+ & iset,ihset,nsaxs,saxs_mode
real*8 waga_homology
real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut,
& dist2_cut
+ real*8 Psaxs(maxsaxs),distsaxs(maxsaxs),CSAXS(3,maxsaxs)
logical refstr,pdbref,punch_dist,print_rms,caonly,verbose,
& merge_helices,bxfile,cxfile,histfile,entfile,zscfile,
& rmsrgymap,with_dihed_constr,check_conf,histout,out1file,
common /homol/ waga_homology(maxR),
& waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut,
& iset,ihset,l_homo(max_template,maxdim)
-
+ common /saxsretr/ Psaxs,distsaxs,csaxs,nsaxs,saxs_mode
c integer ithetaconstr_start_homo,ithetaconstr_end_homo
c
integer nresn,nyosh,nnos
- common /back_constr/ uconst_back,
+ common /back_constr/ uconst_back,uscdiff,
& dutheta,dugamma,duscdiff,duscdiffx
common /homrestr/ odl,dih,sigma_dih,sigma_odl,
& lim_odl,lim_dih,ires_homo,jres_homo,link_start_homo,
--- /dev/null
+ double precision scal_fric,rwat,etawat,gamp,
+ & gamsc(ntyp1),stdfp,stdfsc(ntyp),stdforcp(MAXRES),
+ & stdforcsc(MAXRES),pstok,restok(ntyp+1),cPoise,Rb
+ common /langevin/ pstok,restok,gamp,gamsc,
+ & stdfp,stdfsc,stdforcp,stdforcsc,rwat,etawat,cPoise,Rb
+ double precision IP,ISC(ntyp+1),mp,
+ & msc(ntyp+1)
+ common /inertia/ IP,ISC,MP,MSC
C Maximum number of templates in homology-modeling restraints
integer max_template
parameter(max_template=25)
+C Maximum number of bins in SAXS restraints
+ integer MaxSAXS
+ parameter (MaxSAXS=1000)
c Maximum number of structures in the database, energy components, proteins,
c and structural classes
c#ifdef JUBL
- parameter (maxstr=200000,max_ene=25,maxprot=7,maxclass=5000)
+ parameter (maxstr=200000,max_ene=26,maxprot=7,maxclass=5000)
parameter (maxclass1=10)
c Maximum number of structures to be dealt with by one processor
parameter (maxstr_proc=10000)
& " the value read in: ",energia(0),eini," point",
& iii+1,indstart(me1)+iii," T",
& 1.0d0/(1.987D-3*beta_h(ib,ipar))
- call pdbout(indstart(me1)+iii,
- & 1.0d0/(1.987D-3*beta_h(ib,ipar)),
- &energia(0),eini,0.0d0,0.0d0)
+#ifdef DEBUG
+c call pdbout(indstart(me1)+iii,
+c & 1.0d0/(1.987D-3*beta_h(ib,ipar)),
+c &energia(0),eini,0.0d0,0.0d0)
+ write (iout,*) "wsaxs",wsaxs
call enerprint(energia(0),fT)
+#endif
errmsg_count=errmsg_count+1
if (errmsg_count.gt.maxerrmsg_count)
& write (iout,*) "Too many warning messages"
endif
endif
potE(iii+1,iparm)=energia(0)
- do k=1,22
+ do k=1,max_ene
enetb(k,iii+1,iparm)=energia(k)
enddo
#ifdef DEBUG
c call enerprint(energia(0),fT)
#endif
#ifdef DEBUG
- write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
- write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
+ write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
+ & ((c(l,k+nres),l=1,3),k=nnt,nct)
write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
endif
+
+ write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr
+ if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
+ call e_saxs(Esaxs_constr)
+ 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_ENE1 constr_homology=",constr_homology
if (constr_homology.ge.1) then
call e_modeller(ehomology_constr)
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +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
+ & +wbond*estr+wsccor*fact(1)*esccor+wsaxs*esaxs_constr
#else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
& +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
+ & +wbond*estr+wsccor*fact(1)*esccor+wsaxs*esaxs_constr
#endif
energia(0)=etot
energia(1)=evdw
energia(20)=edihcnstr
energia(21)=evdw_t
energia(22)=ehomology_constr
+ energia(26)=esaxs_constr
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
#ifdef MPL
c endif
#endif
+#define DEBUG
#ifdef DEBUG
call enerprint(energia,fact)
#endif
+#undef DEBUG
if (calc_grad) then
C
C Sum up the components of the Cartesian gradient.
edihcnstr=energia(20)
estr=energia(18)
ehomology_constr=energia(22)
+ esaxs_constr=energia(26)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
& wvdwpp,
& 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,
+ & esaxs_constr*wsaxs,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
C endif
return
end
-
+c----------------------------------------------------------------------------
+ subroutine e_saxs(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(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 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
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+#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=0.25d0/(restok(itype(j))**2+restok(itype(i))**2)
+ do k=1,nsaxs
+ dk = distsaxs(k)
+ expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)
+ Pcalc(k) = Pcalc(k)+expCACA
+#ifdef DEBUG
+ write(iout,*) "i j k Pcalc",i,j,Pcalc(k)
+#endif
+ CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA
+ 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
+ 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)
+ do k=1,nsaxs
+ 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
+ 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
& gshieldc, gshieldc_loc, gshieldx_ec, gshieldc_ec,
& gshieldc_loc_ec, gshieldx_t3,gshieldc_t3,gshieldc_loc_t3,
& gshieldx_t4, gshieldc_t4,gshieldc_loc_t4,gshieldx_ll,
- & gshieldc_ll, gshieldc_loc_ll
+ & gshieldc_ll, gshieldc_loc_ll,gsaxsC,gsaxsX
integer nfl,icg
logical calc_grad
& gshieldx_ll(3,-1:maxres), gshieldc_ll(3,-1:maxres),
& gshieldc_loc_ll(3,-1:maxres),
& gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres),
+ & gsaxsC(3,-1:maxres),gsaxsX(3,-1:maxres),
& gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres),
& gradcorr5(3,maxres),gradcorr6(3,maxres),
& gel_loc(3,maxres),gcorr3_turn(3,maxres),gcorr4_turn(3,maxres),
double precision wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
& wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
& wturn6,wvdwpp,wbond,weights,scal14,cutoff_corr,delt_corr,
- & r0_corr,wliptran
+ & r0_corr,wliptran,wsaxs
integer ipot,n_ene_comp
common /ffield/ wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
& wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
- & wturn6,wvdwpp,wbond,wliptran,weights(max_ene),
+ & wturn6,wvdwpp,wbond,wliptran,wsaxs,weights(max_ene),
& scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp
common /potentials/ potname(5)
character*3 potname
& iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
& ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
& iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
- & iint_end,iphi1_start,iphi1_end,itau_start,itau_end
+ & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,
+ & isaxs_start,isaxs_end
common /peptbond/ vbl,vblinv,vblinv2,vbl_cis,vbl0
common /indices/ loc_start,loc_end,ithet_start,ithet_end,
& iphi_start,iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
& ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
& iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
- & iint_end,iphi1_start,iphi1_end,itau_start,itau_end
+ & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,
+ & isaxs_start,isaxs_end
C Inverses of the actual virtual bond lengths
common /invlen/ vbld_inv(maxres2)
& "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
& "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
& "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
- & "EAFM","ETHETC","EMPTY"/
+ & "EAFM","ETHETC","ESHIELD","ESAXS"/
data wname /
& "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
& "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
& "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
- & "WLIPTRAN","WAFM","WTHETC","WSHIELD"/
+ & "WLIPTRAN","WAFM","WTHETC","WSHIELD","WSAXS"/
data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
& 1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
- & 0.0d0,0.0,0.0d0,0.0d0,0.0d0,0.0d0/
+ & 0.0d0,0.0,0.0d0,0.0d0,0.0d0,0.0d0,0.0d0/
data nprint_ene /22/
data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
- & 16,15,17,20,21,24,22,23,1/
+ & 16,15,17,20,21,24,22,23,26/
end
c---------------------------------------------------------------------------
subroutine init_int_table
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
#ifdef MPI
include 'mpif.h'
#endif
#ifdef MP
include 'COMMON.INFO'
#endif
+ include 'COMMON.CONTROL'
include 'COMMON.CHAIN'
include 'COMMON.INTERACT'
include 'COMMON.LOCAL'
call int_bounds(nres-3,itau_start,itau_end)
itau_start=itau_start+3
itau_end=itau_end+3
+ call int_bounds(nsaxs,isaxs_start,isaxs_end)
if (lprint) then
write (iout,*) 'Processor:',MyID,
& ' loc_start',loc_start,' loc_end',loc_end,
iphi_end=nct
itau_start=4
itau_end=nres
+ isaxs_start=1
+ isaxs_end=nsaxs
#endif
return
end
double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
& eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
- & escloc,ehomology_constr,
+ & escloc,ehomology_constr,esaxs,
& ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
& eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
ehomology_constr=enetb(22,i,iparm)
+ esaxs=enetb(26,i,iparm)
#ifdef SPLITELE
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr+ehomology_constr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
#else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr+ehomology_constr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
#endif
#ifdef MPI
Fdimless_(i)=
if (itype(1).eq.ntyp1) nnt=2
if (itype(nres).eq.ntyp1) nct=nct-1
write(iout,*) 'NNT=',NNT,' NCT=',NCT
+ write (iout,*) "calling read_saxs_consrtr",nsaxs
+ if (nsaxs.gt.0) call read_saxs_constr
+
if (constr_homology.gt.0) then
c write (iout,*) "About to call read_constr_homology"
c call flush(iout)
return
10 return1
end
+c-------------------------------------------------------------------------------
+ subroutine read_saxs_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ double precision cm(3)
+c read(inp,*) nsaxs
+ write (iout,*) "Calling read_saxs nsaxs",nsaxs
+ call flush(iout)
+ if (saxs_mode.eq.0) then
+c SAXS distance distribution
+ do i=1,nsaxs
+ read(inp,*) distsaxs(i),Psaxs(i)
+ enddo
+ Cnorm = 0.0d0
+ do i=1,nsaxs
+ Cnorm = Cnorm + Psaxs(i)
+ enddo
+ write (iout,*) "Cnorm",Cnorm
+ do i=1,nsaxs
+ Psaxs(i)=Psaxs(i)/Cnorm
+ enddo
+ write (iout,*) "Normalized distance distribution from SAXS"
+ do i=1,nsaxs
+ write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
+ enddo
+ else
+c SAXS "spheres".
+ do i=1,nsaxs
+ read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
+ enddo
+ do j=1,3
+ cm(j)=0.0d0
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ cm(j)=cm(j)+Csaxs(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/nsaxs
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ Csaxs(j,i)=Csaxs(j,i)-cm(j)
+ enddo
+ enddo
+ write (iout,*) "SAXS sphere coordinates"
+ do i=1,nsaxs
+ write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
+ enddo
+ endif
+ return
+ end
c====-------------------------------------------------------------------
subroutine read_constr_homology
include 'COMMON.SCROT'
include 'COMMON.FREE'
include 'COMMON.CONTROL'
+ include 'COMMON.LANGEVIN'
character*1 t1,t2,t3
character*1 onelett(4) /"G","A","P","D"/
character*1 toronelet(-2:2) /"p","a","G","A","P"/
external ilen
character*16 key
integer iparm
- double precision ip,mp
character*6 res1
C
C Body
wvdwpp=ww(16)
wbond=ww(18)
wsccor=ww(19)
+ wsaxs=ww(26)
endif
c and Stokes' radii of the peptide group and side chains
c
#ifdef CRYST_BOND
- read (ibond,*) vbldp0,vbldpdum,akp
+ read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
do i=1,ntyp
nbondterm(i)=1
- read (ibond,*) vbldsc0(1,i),aksc(1,i)
+ read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
dsc(i) = vbldsc0(1,i)
if (i.eq.10) then
dsc_inv(i)=0.0D0
endif
enddo
#else
- read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
+ read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
do i=1,ntyp
read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
- & j=1,nbondterm(i))
+ & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
dsc(i) = vbldsc0(1,i)
if (i.eq.10) then
dsc_inv(i)=0.0D0
refstr = index(controlcard,'REFSTR').gt.0
pdbref = index(controlcard,'PDBREF').gt.0
dyn_ss=(index(controlcard,'DYN_SS').gt.0)
+ call readi(controlcard,'NSAXS',nsaxs,0)
+ call readi(controlcard,'SAXS_MODE',saxs_mode,0)
+ write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+ & SAXS_MODE
C /06/28/2013 Adasko: dyn_ss is keyword allowing to break and create bond
C disulfide bond. Note that in conterary to dynamics this in
C CONTROLCARD. The bond are read in molread_zs.F
ww_all(17,iparm)=wbond
ww_all(19,iparm)=wsccor
ww_all(22,iparm)=wliptran
+ ww_all(25,iparm)=wsaxs
c Store bond parameters
vbldp0_all(iparm)=vbldp0
akp_all(iparm)=akp
wbond=ww_all(17,iparm)
wsccor=ww_all(19,iparm)
wliptran=ww_all(22,iparm)
+ wsaxs=ww_all(25,iparm)
c Restore bond parameters
vbldp0=vbldp0_all(iparm)
akp=akp_all(iparm)
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
& escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
& eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
- & ehomology_constr
+ & ehomology_constr,esaxs
integer ind_point(maxpoint),upindE,indE
do iparm=1,nParmSet
#ifdef DEBUG
write (iout,'(2i5,21f8.2)') i,iparm,
- & (enetb(k,i,iparm),k=1,22)
+ & (enetb(k,i,iparm),k=1,26)
#endif
call restore_parm(iparm)
#ifdef DEBUG
write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
- & wtor_d,wsccor,wbond
+ & wtor_d,wsccor,wbond,wsaxs
#endif
do ib=1,nT_h(iparm)
if (rescale_mode.eq.1) then
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
ehomology_constr=enetb(22,i,iparm)
+ esaxs=enetb(26,i,iparm)
#ifdef DEBUG
write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
& evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
- & etors,etors_d,eello_turn3,eello_turn4,esccor
+ & etors,etors_d,eello_turn3,eello_turn4,esccor,esaxs
#endif
#ifdef SPLITELE
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wsaxs*esaxs
#else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wsaxs*esaxs
#endif
#ifdef DEBUG
write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
do iparm=1,nParmSet
#ifdef DEBUG
write (iout,'(2i5,21f8.2)') i,iparm,
- & (enetb(k,i,iparm),k=1,22)
+ & (enetb(k,i,iparm),k=1,26)
#endif
call restore_parm(iparm)
#ifdef DEBUG
write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
- & wtor_d,wsccor,wbond
+ & wtor_d,wsccor,wbond,wsaxs
#endif
do ib=1,nT_h(iparm)
if (rescale_mode.eq.1) then
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
ehomology_constr=enetb(22,i,iparm)
+ esaxs=enetb(26,i,iparm)
#ifdef SPLITELE
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr+ehomology_constr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
#else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr+ehomology_constr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
#endif
etot=etot-entfac(i)/beta_h(ib,iparm)
esccor=enetb(19,t,iparm)
edihcnstr=enetb(20,t,iparm)
ehomology_constr=enetb(22,t,iparm)
+ esaxs=enetb(26,t,iparm)
if (homol_nset.gt.1)
& ehomology_constr=waga_homology(ihset)*ehomology_constr
do k=0,nGridT
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr+ehomology_constr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
& +ftprim(1)*wtor*etors+
& ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr+ehomology_constr
+ & +wbond*estr+ehomology_constr+wsaxs*esaxs
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
& +ftprim(1)*wtor*etors+
& ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
projects / unres.git / commitdiff