From: Cezary Czaplewski Date: Sun, 17 Dec 2017 21:52:36 +0000 (+0100) Subject: saxs and adam's corrections to multichain X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=d24fd92f8762d3b48a94c699ef90e0ab37d3c204;p=unres.git saxs and adam's corrections to multichain --- diff --git a/ctest/prota_unres_energy_check_mult.sh b/ctest/prota_unres_energy_check_mult.sh index 565675a..8831cfe 100755 --- a/ctest/prota_unres_energy_check_mult.sh +++ b/ctest/prota_unres_energy_check_mult.sh @@ -57,7 +57,7 @@ elif [ "$1" == "1l2y_MIN_REGULAR_INT" ]; then 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]} diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL index ec85f0b..bc4444e 100644 --- a/source/cluster/wham/src-M/COMMON.CONTROL +++ b/source/cluster/wham/src-M/COMMON.CONTROL @@ -2,10 +2,11 @@ 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, @@ -23,3 +24,4 @@ & 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 diff --git a/source/cluster/wham/src-M/COMMON.DERIV b/source/cluster/wham/src-M/COMMON.DERIV index 4e116e1..c1b9727 100644 --- a/source/cluster/wham/src-M/COMMON.DERIV +++ b/source/cluster/wham/src-M/COMMON.DERIV @@ -7,7 +7,8 @@ & 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 @@ -29,6 +30,7 @@ & 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), diff --git a/source/cluster/wham/src-M/COMMON.FFIELD b/source/cluster/wham/src-M/COMMON.FFIELD index fa85436..8b0f907 100644 --- a/source/cluster/wham/src-M/COMMON.FFIELD +++ b/source/cluster/wham/src-M/COMMON.FFIELD @@ -6,11 +6,11 @@ C----------------------------------------------------------------------- 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) diff --git a/source/cluster/wham/src-M/COMMON.HOMRESTR b/source/cluster/wham/src-M/COMMON.HOMRESTR index 95ea932..5c23caf 100644 --- a/source/cluster/wham/src-M/COMMON.HOMRESTR +++ b/source/cluster/wham/src-M/COMMON.HOMRESTR @@ -27,7 +27,7 @@ c 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, diff --git a/source/cluster/wham/src-M/COMMON.LANGEVIN b/source/cluster/wham/src-M/COMMON.LANGEVIN new file mode 100644 index 0000000..982bde9 --- /dev/null +++ b/source/cluster/wham/src-M/COMMON.LANGEVIN @@ -0,0 +1,8 @@ + 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 diff --git a/source/cluster/wham/src-M/COMMON.LOCAL b/source/cluster/wham/src-M/COMMON.LOCAL index d92ce8e..8563e53 100644 --- a/source/cluster/wham/src-M/COMMON.LOCAL +++ b/source/cluster/wham/src-M/COMMON.LOCAL @@ -2,7 +2,7 @@ & 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), @@ -39,6 +39,6 @@ C Parameters of the side-chain probability distribution 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) diff --git a/source/cluster/wham/src-M/DIMENSIONS b/source/cluster/wham/src-M/DIMENSIONS index 1fef045..2965259 100644 --- a/source/cluster/wham/src-M/DIMENSIONS +++ b/source/cluster/wham/src-M/DIMENSIONS @@ -73,3 +73,6 @@ C Maximum number of terms in SC bond-stretching potential 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) diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index 968a9f1..340af75 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -121,7 +121,15 @@ c print *,ecorr,ecorr5,ecorr6,eturn6 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) @@ -143,7 +151,7 @@ c write (iout,*) "ft(6)",fact(6),wliptran,eliptran & +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 @@ -153,7 +161,7 @@ c write (iout,*) "ft(6)",fact(6),wliptran,eliptran & +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 @@ -165,7 +173,7 @@ c write (iout,*) "ft(6)",fact(6),wliptran,eliptran & +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) @@ -175,7 +183,7 @@ c write (iout,*) "ft(6)",fact(6),wliptran,eliptran & +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 @@ -212,6 +220,7 @@ c write (iout,*) "ft(6)",fact(6),wliptran,eliptran 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 @@ -398,6 +407,7 @@ C------------------------------------------------------------------------ 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, @@ -407,8 +417,8 @@ c ethetacnstr=energia(24) & 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)'/ @@ -431,6 +441,7 @@ c ethetacnstr=energia(24) & '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 @@ -440,7 +451,7 @@ c ethetacnstr=energia(24) & 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)'/ @@ -463,6 +474,7 @@ c ethetacnstr=energia(24) & '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 @@ -1182,6 +1194,7 @@ C 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 @@ -2127,6 +2140,7 @@ C 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), @@ -3542,6 +3556,7 @@ C include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.IOUNITS' + integer xshift,yshift,zshift dimension ggg(3) evdw2=0.0D0 evdw2_14=0.0d0 @@ -9333,7 +9348,7 @@ C gradient po costhet & )*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))*( @@ -9501,7 +9516,7 @@ C gradient po costhet &*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)) @@ -9666,4 +9681,351 @@ C eliptran=elpitran+0.0 ! I am in water 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 diff --git a/source/cluster/wham/src-M/initialize_p.F b/source/cluster/wham/src-M/initialize_p.F index 40d035c..4a6aa4c 100644 --- a/source/cluster/wham/src-M/initialize_p.F +++ b/source/cluster/wham/src-M/initialize_p.F @@ -260,15 +260,15 @@ c------------------------------------------------------------------------- & "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 @@ -281,6 +281,7 @@ c--------------------------------------------------------------------------- #ifdef MPL include 'COMMON.INFO' #endif + include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.LOCAL' @@ -516,6 +517,7 @@ C Partition local interactions 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, @@ -541,6 +543,9 @@ C Partition local interactions 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 diff --git a/source/cluster/wham/src-M/parmread.F b/source/cluster/wham/src-M/parmread.F index 8de6af5..c9617ec 100644 --- a/source/cluster/wham/src-M/parmread.F +++ b/source/cluster/wham/src-M/parmread.F @@ -17,11 +17,11 @@ C 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 @@ -34,10 +34,10 @@ C Assign virtual-bond length 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 @@ -46,10 +46,10 @@ C Assign virtual-bond length 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 diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F index 1ce5ba9..22db88f 100644 --- a/source/cluster/wham/src-M/probabl.F +++ b/source/cluster/wham/src-M/probabl.F @@ -21,7 +21,7 @@ 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 @@ -191,6 +191,7 @@ c write (iout,*) evdw 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 @@ -318,6 +319,9 @@ c & " indend",indend(iproc) 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 !-------------------------------------------------- diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 9e80118..6bb388a 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -13,10 +13,11 @@ C 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' @@ -127,6 +128,10 @@ C long axis of side chain 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 @@ -187,6 +192,7 @@ C Read weights of the subsequent energy terms. 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) @@ -415,6 +421,8 @@ C enddo 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) @@ -956,6 +964,66 @@ C endif 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) diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index dd53a8e..753b832 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -21,6 +21,7 @@ & 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 diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index 35fb03d..be1739c 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -1,9 +1,10 @@ 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, @@ -18,5 +19,6 @@ & 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 diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV index 8082b0a..407acfd 100644 --- a/source/unres/src_MD-M/COMMON.DERIV +++ b/source/unres/src_MD-M/COMMON.DERIV @@ -1,8 +1,11 @@ - 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), @@ -12,8 +15,9 @@ & 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), diff --git a/source/unres/src_MD-M/COMMON.FFIELD b/source/unres/src_MD-M/COMMON.FFIELD index a63fe78..44aa937 100644 --- a/source/unres/src_MD-M/COMMON.FFIELD +++ b/source/unres/src_MD-M/COMMON.FFIELD @@ -3,10 +3,14 @@ C The following COMMON block selects the type of the force field used in 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) diff --git a/source/unres/src_MD-M/COMMON.LOCAL b/source/unres/src_MD-M/COMMON.LOCAL index 5d1ced7..631f593 100644 --- a/source/unres/src_MD-M/COMMON.LOCAL +++ b/source/unres/src_MD-M/COMMON.LOCAL @@ -41,7 +41,7 @@ C Virtual-bond lenghts & 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), @@ -57,6 +57,7 @@ C Virtual-bond lenghts & 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) diff --git a/source/unres/src_MD-M/COMMON.SBRIDGE b/source/unres/src_MD-M/COMMON.SBRIDGE index 91dd2cd..f866aa7 100644 --- a/source/unres/src_MD-M/COMMON.SBRIDGE +++ b/source/unres/src_MD-M/COMMON.SBRIDGE @@ -3,7 +3,7 @@ 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 diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index 635b007..29c29b5 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -95,7 +95,7 @@ C Max. number of conformations in the pool 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) @@ -142,3 +142,6 @@ C to master; depends on nstex / ntwx ratio 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) diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 070476d..84c5a18 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -36,6 +36,7 @@ c------------------------------------------------ common /gucio/ cm integer itime logical ovrtim + integer nharp,iharp(4,maxres/3) c #ifdef MPI if (ilen(tmpdir).gt.0) @@ -51,6 +52,33 @@ c 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 @@ -1443,7 +1471,7 @@ c Set up the initial conditions of a MD simulation 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 @@ -1632,10 +1660,13 @@ c Removing the velocity of the center of mass 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) diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index c0b4c84..4e7691d 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -1452,11 +1452,19 @@ C 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) @@ -1478,12 +1486,17 @@ C 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) @@ -1509,10 +1522,13 @@ c 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) @@ -1531,10 +1547,11 @@ c 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 diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 2e80cf1..205f78a 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -55,6 +55,7 @@ C FG slaves as WEIGHTS array. 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) @@ -81,6 +82,7 @@ C FG slaves receive the WEIGHTS array 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 @@ -266,6 +268,16 @@ cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 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 @@ -334,6 +346,7 @@ C 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" @@ -428,6 +441,9 @@ cMS$ATTRIBUTES C :: proc_proc 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 @@ -435,7 +451,7 @@ cMS$ATTRIBUTES C :: proc_proc & +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 @@ -443,7 +459,7 @@ cMS$ATTRIBUTES C :: proc_proc & +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 @@ -502,12 +518,20 @@ cMS$ATTRIBUTES C :: proc_proc #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) @@ -547,7 +571,8 @@ c enddo & 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) @@ -565,7 +590,8 @@ c enddo & 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) @@ -625,7 +651,8 @@ c enddo 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 @@ -633,7 +660,7 @@ 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) + write (iout,'(i3,3e15.5)') i,(gradbufc(j,i),j=1,3) enddo call flush(iout) #endif @@ -641,8 +668,8 @@ c enddo #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 @@ -655,7 +682,8 @@ c enddo 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 @@ -672,8 +700,8 @@ c 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 @@ -731,8 +759,9 @@ c enddo #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 @@ -929,7 +958,7 @@ c #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 @@ -1032,6 +1061,7 @@ C------------------------------------------------------------------------ Uconst=energia(20) esccor=energia(21) ehomology_constr=energia(24) + esaxs_constr=energia(25) eliptran=energia(22) Eafmforce=energia(23) #ifdef SPLITELE @@ -1041,7 +1071,7 @@ C------------------------------------------------------------------------ & 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)'/ @@ -1065,6 +1095,7 @@ C------------------------------------------------------------------------ & '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)'/ @@ -1078,7 +1109,7 @@ C------------------------------------------------------------------------ & 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:'// @@ -1102,6 +1133,7 @@ C------------------------------------------------------------------------ & '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)'/ @@ -11017,3 +11049,354 @@ C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist 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 diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F index 24ab8dd..0f67d19 100644 --- a/source/unres/src_MD-M/energy_split-sep.F +++ b/source/unres/src_MD-M/energy_split-sep.F @@ -15,6 +15,7 @@ cMS$ATTRIBUTES C :: proc_proc double precision weights_(n_ene) #endif include 'COMMON.SETUP' + include 'COMMON.CONTROL' include 'COMMON.IOUNITS' double precision energia(0:n_ene) include 'COMMON.FFIELD' @@ -65,6 +66,7 @@ C FG slaves as WEIGHTS array. 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) @@ -91,6 +93,7 @@ C FG slaves receive the WEIGHTS array 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) @@ -252,6 +255,7 @@ cMS$ATTRIBUTES C :: proc_proc double precision weights_(n_ene) #endif include 'COMMON.SETUP' + include 'COMMON.CONTROL' include 'COMMON.IOUNITS' double precision energia(0:n_ene) include 'COMMON.FFIELD' @@ -303,6 +307,7 @@ C FG slaves as WEIGHTS array. 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) @@ -329,6 +334,7 @@ C FG slaves receive the WEIGHTS array 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, @@ -436,6 +442,17 @@ C 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 @@ -463,6 +480,7 @@ 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.) diff --git a/source/unres/src_MD-M/gen_rand_conf.F b/source/unres/src_MD-M/gen_rand_conf.F index 6caa718..4b6bbc0 100644 --- a/source/unres/src_MD-M/gen_rand_conf.F +++ b/source/unres/src_MD-M/gen_rand_conf.F @@ -781,7 +781,7 @@ c overlapping residues left, or false otherwise (success) 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) @@ -842,6 +842,7 @@ c print *,'>>overlap_sc nnt=',nnt,' nct=',nct 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) @@ -854,6 +855,7 @@ c 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) diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F index 3ad3912..ec3ae46 100644 --- a/source/unres/src_MD-M/geomout.F +++ b/source/unres/src_MD-M/geomout.F @@ -156,7 +156,7 @@ C format. 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 @@ -514,6 +514,11 @@ C print *,'A CHUJ',potEcomp(23) 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)" diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F index 7295d2e..d0709d2 100644 --- a/source/unres/src_MD-M/gradient_p.F +++ b/source/unres/src_MD-M/gradient_p.F @@ -271,6 +271,13 @@ c 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 @@ -365,6 +372,8 @@ C 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 diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index 446cad5..d789d82 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -333,16 +333,16 @@ c------------------------------------------------------------------------- & "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 @@ -628,6 +628,11 @@ C Partition local interactions 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 @@ -741,12 +746,16 @@ c nlen=nres-nnt+1 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, @@ -1169,6 +1178,8 @@ c write (iout,*) "MPI_ROTAT2",MPI_ROTAT2 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 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index a8f75a7..bbfcabb 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -358,6 +358,7 @@ C Control printout of the coefficients of virtual-bond-angle potentials 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 @@ -392,6 +393,7 @@ C enddo enddo enddo + enddo call flush(iout) endif c write (2,*) "Start reading THETA_PDB",ithep_pdb @@ -658,6 +660,7 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) 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 @@ -673,6 +676,7 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo enddo enddo + enddo endif C @@ -879,7 +883,7 @@ cc maxinter is maximum interaction sites 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 @@ -888,6 +892,7 @@ cc maxinter is maximum interaction sites #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 @@ -902,6 +907,7 @@ cc maxinter is maximum interaction sites enddo enddo enddo + enddo endif C diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 69c0816..f18e6eb 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -98,6 +98,11 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file 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 @@ -619,6 +624,7 @@ C Read weights of the subsequent energy terms. 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) @@ -642,6 +648,7 @@ C 12/1/95 Added weight for the multi-body term WCORR 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, @@ -803,12 +810,16 @@ c print *,'Finished reading pdb data' 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 @@ -869,6 +880,7 @@ C 8/13/98 Set limits to generating the dihedral angles 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) @@ -998,6 +1010,8 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb enddo endif endif + write (iout,*) "calling read_saxs_consrtr",nsaxs + if (nsaxs.gt.0) call read_saxs_constr if (constr_homology.gt.0) then @@ -2449,6 +2463,67 @@ CCCC NOW PROPERTIES FOR AFM 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' diff --git a/source/unres/src_MD-M/sc_move.F b/source/unres/src_MD-M/sc_move.F index c552ee0..af6e99d 100644 --- a/source/unres/src_MD-M/sc_move.F +++ b/source/unres/src_MD-M/sc_move.F @@ -720,6 +720,7 @@ c if (icall.eq.0) lprn=.true. 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) @@ -765,6 +766,7 @@ C 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) diff --git a/source/unres/src_MD-M/ssMD.F b/source/unres/src_MD-M/ssMD.F index 3c89ef5..7c3b6b6 100644 --- a/source/unres/src_MD-M/ssMD.F +++ b/source/unres/src_MD-M/ssMD.F @@ -650,7 +650,7 @@ c Local variables & 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 diff --git a/source/wham/src-M/COMMON.CONTROL b/source/wham/src-M/COMMON.CONTROL index 4f4fb65..8c940ab 100644 --- a/source/wham/src-M/COMMON.CONTROL +++ b/source/wham/src-M/COMMON.CONTROL @@ -1,10 +1,11 @@ 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, @@ -20,4 +21,4 @@ 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 diff --git a/source/wham/src-M/COMMON.HOMRESTR b/source/wham/src-M/COMMON.HOMRESTR index 95ea932..5c23caf 100644 --- a/source/wham/src-M/COMMON.HOMRESTR +++ b/source/wham/src-M/COMMON.HOMRESTR @@ -27,7 +27,7 @@ c 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, diff --git a/source/wham/src-M/COMMON.LANGEVIN b/source/wham/src-M/COMMON.LANGEVIN new file mode 100644 index 0000000..982bde9 --- /dev/null +++ b/source/wham/src-M/COMMON.LANGEVIN @@ -0,0 +1,8 @@ + 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 diff --git a/source/wham/src-M/DIMENSIONS b/source/wham/src-M/DIMENSIONS index d2d54bc..2fe913c 100644 --- a/source/wham/src-M/DIMENSIONS +++ b/source/wham/src-M/DIMENSIONS @@ -146,3 +146,6 @@ C Maximum number of terms in SC bond-stretching potential 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) diff --git a/source/wham/src-M/DIMENSIONS.ZSCOPT b/source/wham/src-M/DIMENSIONS.ZSCOPT index 6519938..a99d1ed 100644 --- a/source/wham/src-M/DIMENSIONS.ZSCOPT +++ b/source/wham/src-M/DIMENSIONS.ZSCOPT @@ -3,7 +3,7 @@ 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) diff --git a/source/wham/src-M/enecalc1.F b/source/wham/src-M/enecalc1.F index 4fb7c9d..5618f74 100644 --- a/source/wham/src-M/enecalc1.F +++ b/source/wham/src-M/enecalc1.F @@ -208,10 +208,13 @@ c call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) & " 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" @@ -227,7 +230,7 @@ c call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) 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 @@ -236,8 +239,8 @@ c call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) 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) diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 6b893c4..de97b0b 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -109,6 +109,18 @@ c print *,ecorr,ecorr5,ecorr6,eturn6 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) @@ -126,7 +138,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +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) @@ -135,7 +147,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +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 @@ -170,6 +182,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(20)=edihcnstr energia(21)=evdw_t energia(22)=ehomology_constr + energia(26)=esaxs_constr c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -189,9 +202,11 @@ c detecting NaNQ #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. @@ -293,6 +308,7 @@ C------------------------------------------------------------------------ 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, @@ -301,7 +317,9 @@ C------------------------------------------------------------------------ & 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)'/ @@ -324,6 +342,7 @@ C------------------------------------------------------------------------ & '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 @@ -333,7 +352,7 @@ C------------------------------------------------------------------------ & 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)'/ @@ -356,6 +375,7 @@ C------------------------------------------------------------------------ & '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 @@ -8648,4 +8668,355 @@ C sscagrad=0.0d0 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 diff --git a/source/wham/src-M/include_unres/COMMON.DERIV b/source/wham/src-M/include_unres/COMMON.DERIV index 7f8ddfb..e461169 100644 --- a/source/wham/src-M/include_unres/COMMON.DERIV +++ b/source/wham/src-M/include_unres/COMMON.DERIV @@ -8,7 +8,7 @@ & 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 @@ -29,6 +29,7 @@ & 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), diff --git a/source/wham/src-M/include_unres/COMMON.FFIELD b/source/wham/src-M/include_unres/COMMON.FFIELD index 6c432a9..08be325 100644 --- a/source/wham/src-M/include_unres/COMMON.FFIELD +++ b/source/wham/src-M/include_unres/COMMON.FFIELD @@ -6,11 +6,11 @@ C----------------------------------------------------------------------- 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 diff --git a/source/wham/src-M/include_unres/COMMON.LOCAL b/source/wham/src-M/include_unres/COMMON.LOCAL index 3e68e82..88a984b 100644 --- a/source/wham/src-M/include_unres/COMMON.LOCAL +++ b/source/wham/src-M/include_unres/COMMON.LOCAL @@ -42,12 +42,14 @@ C Virtual-bond lenghts & 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) diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F index 06e556e..76d561f 100644 --- a/source/wham/src-M/initialize_p.F +++ b/source/wham/src-M/initialize_p.F @@ -267,30 +267,32 @@ c------------------------------------------------------------------------- & "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' @@ -526,6 +528,7 @@ C Partition local interactions 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, @@ -551,6 +554,8 @@ C Partition local interactions iphi_end=nct itau_start=4 itau_end=nres + isaxs_start=1 + isaxs_end=nsaxs #endif return end diff --git a/source/wham/src-M/make_ensemble1.F b/source/wham/src-M/make_ensemble1.F index 7e2f88c..38838ec 100644 --- a/source/wham/src-M/make_ensemble1.F +++ b/source/wham/src-M/make_ensemble1.F @@ -23,7 +23,7 @@ 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 @@ -163,6 +163,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -172,7 +173,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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) @@ -182,7 +183,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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)= diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F index 3480671..55c3fec 100644 --- a/source/wham/src-M/molread_zs.F +++ b/source/wham/src-M/molread_zs.F @@ -137,6 +137,9 @@ C enddo 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) @@ -324,6 +327,68 @@ c print *,"energy",energ," iscor",iscor 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 diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F index 963212c..4315d20 100644 --- a/source/wham/src-M/parmread.F +++ b/source/wham/src-M/parmread.F @@ -22,6 +22,7 @@ C 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"/ @@ -35,7 +36,6 @@ C external ilen character*16 key integer iparm - double precision ip,mp character*6 res1 C C Body @@ -140,6 +140,7 @@ c wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) + wsaxs=ww(26) endif @@ -212,10 +213,10 @@ c Read the virtual-bond parameters, masses, and moments of inertia 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 @@ -224,10 +225,10 @@ c 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 diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F index 271fca8..51216d9 100644 --- a/source/wham/src-M/readrtns.F +++ b/source/wham/src-M/readrtns.F @@ -126,6 +126,10 @@ C endif 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 diff --git a/source/wham/src-M/store_parm.F b/source/wham/src-M/store_parm.F index f44b6f4..64b5889 100644 --- a/source/wham/src-M/store_parm.F +++ b/source/wham/src-M/store_parm.F @@ -41,6 +41,7 @@ c Store weights 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 @@ -316,6 +317,7 @@ c Restore weights 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) diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index 718bca4..536dae6 100644 --- a/source/wham/src-M/wham_calc1.F +++ b/source/wham/src-M/wham_calc1.F @@ -88,7 +88,7 @@ c parameter (MaxHdim=200000) 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 @@ -238,13 +238,13 @@ c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet) 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 @@ -324,10 +324,11 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -339,7 +340,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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) @@ -349,7 +350,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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), @@ -603,13 +604,13 @@ c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet) 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 @@ -689,6 +690,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -698,7 +700,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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) @@ -708,7 +710,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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) @@ -932,6 +934,7 @@ c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18) 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 @@ -1053,7 +1056,7 @@ c write (iout,*) "ftbis",ftbis & +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+ @@ -1076,7 +1079,7 @@ c write (iout,*) "ftbis",ftbis & +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+