X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=26f7a010f9c68724fb19d602be95b74d61548e0c;hb=855057911aaa7d3d083c300f4a9b1e7323f5d7da;hp=2c40d287b060a928c69452908b30c6b190f9f43a;hpb=3e42b5f53e0ec03cd03009f3b857931ce7ee25c2;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 2c40d28..26f7a01 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 @@ -3201,7 +3221,6 @@ C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' - include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' @@ -3212,9 +3231,13 @@ C include 'COMMON.IOUNITS' dimension ggg(3) ehpb=0.0D0 -cd print *,'edis: nhpb=',nhpb,' fbr=',fbr -cd print *,'link_start=',link_start,' link_end=',link_end -C write(iout,*) link_end, "link_end" + do i=1,3 + ggg(i)=0.0d0 + enddo +C write (iout,*) ,"link_end",link_end,constr_dist +cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr +c write(iout,*)'link_start=',link_start,' link_end=',link_end, +c & " constr_dist",constr_dist if (link_end.eq.0) return do i=link_start,link_end C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a @@ -3229,120 +3252,92 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif +c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, +c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. -C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. C & iabs(itype(jjj)).eq.1) then -C write(iout,*) constr_dist,"const" - if (.not.dyn_ss .and. i.le.nss) then - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. - & iabs(itype(jjj)).eq.1) then - call ssbond_ene(iii,jjj,eij) - ehpb=ehpb+2*eij - endif !ii.gt.neres - else if (ii.gt.nres .and. jj.gt.nres) then -c Restraints from contact prediction +cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds + if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond - additional check + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then + call ssbond_ene(iii,jjj,eij) + ehpb=ehpb+2*eij + endif +cd write (iout,*) "eij",eij +cd & ' waga=',waga,' fac=',fac +! else if (ii.gt.nres .and. jj.gt.nres) then + else +C Calculate the distance between the two points and its difference from the +C target distance. dd=dist(ii,jj) - if (constr_dist.eq.11) then -C ehpb=ehpb+fordepth(i)**4.0d0 -C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) - ehpb=ehpb+fordepth(i)**4.0d0 - & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) - fac=fordepth(i)**4.0d0 - & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd -C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, -C & ehpb,fordepth(i),dd -C write(iout,*) ehpb,"atu?" -C ehpb,"tu?" -C fac=fordepth(i)**4.0d0 -C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd - else - if (dhpb1(i).gt.0.0d0) then - ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + if (irestr_type(i).eq.11) then + ehpb=ehpb+fordepth(i)!**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)!**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +c if (energy_dec) write (iout,'(a6,2i5,6f10.3,i5)') +c & "edisL",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i), +c & ehpb,irestr_type(i) + else if (irestr_type(i).eq.10) then +c AL 6//19/2018 cross-link restraints + xdis = 0.5d0*(dd/forcon(i))**2 + expdis = dexp(-xdis) +c aux=(dhpb(i)+dhpb1(i)*xdis)*expdis+fordepth(i) + aux=(dhpb(i)+dhpb1(i)*xdis*xdis)*expdis+fordepth(i) +c write (iout,*)"HERE: xdis",xdis," expdis",expdis," aux",aux, +c & " wboltzd",wboltzd + ehpb=ehpb-wboltzd*xlscore(i)*dlog(aux) +c fac=-wboltzd*(dhpb1(i)*(1.0d0-xdis)-dhpb(i)) + fac=-wboltzd*xlscore(i)*(dhpb1(i)*(2.0d0-xdis)*xdis-dhpb(i)) + & *expdis/(aux*forcon(i)**2) +c if (energy_dec) write(iout,'(a6,2i5,6f10.3,i5)') +c & "edisX",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i), +c & -wboltzd*xlscore(i)*dlog(aux),irestr_type(i) + else if (irestr_type(i).eq.2) then +c Quartic restraints + ehpb=ehpb+forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) +c if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') +c & "edisQ",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i), +c & forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)),irestr_type(i) fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd -c write (iout,*) "beta nmr", -c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) else - dd=dist(ii,jj) +c Quadratic restraints rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis -c write (iout,*) "beta reg",dd,waga*rdis*rdis + ehpb=ehpb+0.5d0*waga*rdis*rdis +c if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') +c & "edisS",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i), +c & 0.5d0*waga*rdis*rdis,irestr_type(i) C C Evaluate gradient. C fac=waga*rdis/dd - endif !end dhpb1(i).gt.0 - endif !end const_dist=11 + endif +c Calculate Cartesian gradient do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo - do j=1,3 - ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) - ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) - enddo - do k=1,3 - ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) - ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) - enddo - else !ii.gt.nres -C write(iout,*) "before" - dd=dist(ii,jj) -C write(iout,*) "after",dd - if (constr_dist.eq.11) then - ehpb=ehpb+fordepth(i)**4.0d0 - & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) - fac=fordepth(i)**4.0d0 - & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd -C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) -C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd -C print *,ehpb,"tu?" -C write(iout,*) ehpb,"btu?", -C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i) -C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, -C & ehpb,fordepth(i),dd - else - if (dhpb1(i).gt.0.0d0) then - ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) - fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd -c write (iout,*) "alph nmr", -c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) - else - rdis=dd-dhpb(i) -C Get the force constant corresponding to this distance. - waga=forcon(i) -C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis -c write (iout,*) "alpha reg",dd,waga*rdis*rdis -C -C Evaluate gradient. -C - fac=waga*rdis/dd - endif - endif - - do j=1,3 - ggg(j)=fac*(c(j,jj)-c(j,ii)) - enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). - if (iii.lt.ii) then - do j=1,3 - ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) - ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) - enddo - endif - do j=iii,jjj-1 + if (iii.lt.ii) then + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + endif do k=1,3 - ghpbc(k,j)=ghpbc(k,j)+ggg(k) + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) enddo - enddo endif enddo - if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -3485,8 +3480,12 @@ c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d j = jres_homo(ii) dij=dist(i,j) c write (iout,*) "dij(",i,j,") =",dij + nexl=0 do k=1,constr_homology - if(.not.l_homo(k,ii)) cycle + if(.not.l_homo(k,ii)) then + nexl=nexl+1 + cycle + endif distance(k)=odl(k,ii)-dij c write (iout,*) "distance(",k,") =",distance(k) c @@ -3524,7 +3523,15 @@ c write (iout,* )"min_odl",min_odl write (iout,*) "distancek",(distancek(k),k=1,constr_homology) write (iout,* )"min_odl",min_odl #endif +#ifdef OLDRESTR odleg2=0.0d0 +#else + if (waga_dist.ge.0.0d0) then + odleg2=nexl + else + odleg2=0.0d0 + endif +#endif do k=1,constr_homology c Nie wiem po co to liczycie jeszcze raz! c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ @@ -3680,8 +3687,11 @@ c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= c & -(6.28318-dih_diff(i,k)) c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)= c & 6.28318+dih_diff(i,k) - +#ifdef OLD_DIHED kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#else + kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) +#endif c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i) gdih(k)=dexp(kat3) kat2=kat2+gdih(k) @@ -3709,7 +3719,11 @@ c ---------------------------------------------------------------------- sum_gdih=kat2 sum_sgdih=0.0d0 do k=1,constr_homology +#ifdef OLD_DIHED sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd +#else + sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) +#endif c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle sum_sgdih=sum_sgdih+sgdih enddo @@ -4446,6 +4460,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -8647,4 +8662,418 @@ C sscagrad=0.0d0 C endif return end +c---------------------------------------------------------------------------- + double precision function sscale2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) +c write (2,*) "r",r," r_cut",r_cut," r0",r0," rlamb",rlamb +c write (2,*) "rr",rr + if(rr.lt.r_cut-rlamb) then + sscale2=1.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0) + else + sscale2=0d0 + endif + return + end +C----------------------------------------------------------------------- + double precision function sscalgrad2(r,r_cut,r0,rlamb) + implicit none + double precision r,gamm,r_cut,r0,rlamb,rr + rr = dabs(r-r0) + if(rr.lt.r_cut-rlamb) then + sscalgrad2=0.0d0 + else if(rr.le.r_cut.and.rr.ge.r_cut-rlamb) then + gamm=(rr-(r_cut-rlamb))/rlamb + if (r.ge.r0) then + sscalgrad2=gamm*(6*gamm-6.0d0)/rlamb + else + sscalgrad2=-gamm*(6*gamm-6.0d0)/rlamb + endif + else + sscalgrad2=0.0d0 + endif + return + end +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 sss2,ssgrad2,rrr,sscalgrad2,sscale2 + double precision dist + external dist +c SAXS restraint penalty function +#ifdef DEBUG + write(iout,*) "------- SAXS penalty function start -------" + write (iout,*) "nsaxs",nsaxs + write (iout,*) "Esaxs: iatsc_s",iatsc_s," iatsc_e",iatsc_e + write (iout,*) "Psaxs" + do i=1,nsaxs + write (iout,'(i5,e15.5)') i, Psaxs(i) + enddo +#endif + Esaxs_constr = 0.0d0 + do k=1,nsaxs + Pcalc(k)=0.0d0 + do j=1,nres + do l=1,3 + PgradC(k,l,j)=0.0d0 + PgradX(k,l,j)=0.0d0 + enddo + enddo + enddo + do i=iatsc_s,iatsc_e + if (itype(i).eq.ntyp1) cycle + do iint=1,nint_gr(i) + do j=istart(i,iint),iend(i,iint) + if (itype(j).eq.ntyp1) cycle +#ifdef ALLSAXS + dijCACA=dist(i,j) + dijCASC=dist(i,j+nres) + dijSCCA=dist(i+nres,j) + dijSCSC=dist(i+nres,j+nres) + sigma2CACA=2.0d0/(pstok**2) + sigma2CASC=4.0d0/(pstok**2+restok(itype(j))**2) + sigma2SCCA=4.0d0/(pstok**2+restok(itype(i))**2) + sigma2SCSC=4.0d0/(restok(itype(j))**2+restok(itype(i))**2) + do k=1,nsaxs + dk = distsaxs(k) + expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2) + if (itype(j).ne.10) then + expCASC = dexp(-0.5d0*sigma2CASC*(dijCASC-dk)**2) + else + endif + expCASC = 0.0d0 + if (itype(i).ne.10) then + expSCCA = dexp(-0.5d0*sigma2SCCA*(dijSCCA-dk)**2) + else + expSCCA = 0.0d0 + endif + if (itype(i).ne.10 .and. itype(j).ne.10) then + expSCSC = dexp(-0.5d0*sigma2SCSC*(dijSCSC-dk)**2) + else + expSCSC = 0.0d0 + endif + Pcalc(k) = Pcalc(k)+expCACA+expCASC+expSCCA+expSCSC +#ifdef DEBUG + write(iout,*) "i j k Pcalc",i,j,Pcalc(k) +#endif + CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA + CASCgrad = sigma2CASC*(dijCASC-dk)*expCASC + SCCAgrad = sigma2SCCA*(dijSCCA-dk)*expSCCA + SCSCgrad = sigma2SCSC*(dijSCSC-dk)*expSCSC + do l=1,3 +c CA CA + aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux +c CA SC + if (itype(j).ne.10) then + aux = CASCgrad*(C(l,j+nres)-C(l,i))/dijCASC + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + PgradX(k,l,j) = PgradX(k,l,j)+aux + endif +c SC CA + if (itype(i).ne.10) then + aux = SCCAgrad*(C(l,j)-C(l,i+nres))/dijSCCA + PgradX(k,l,i) = PgradX(k,l,i)-aux + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + endif +c SC SC + if (itype(i).ne.10 .and. itype(j).ne.10) then + aux = SCSCgrad*(C(l,j+nres)-C(l,i+nres))/dijSCSC + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + PgradX(k,l,i) = PgradX(k,l,i)-aux + PgradX(k,l,j) = PgradX(k,l,j)+aux + endif + enddo ! l + enddo ! k +#else + dijCACA=dist(i,j) + sigma2CACA=scal_rad**2*0.25d0/ + & (restok(itype(j))**2+restok(itype(i))**2) + IF (saxs_cutoff.eq.0) THEN + do k=1,nsaxs + dk = distsaxs(k) + expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2) + Pcalc(k) = Pcalc(k)+expCACA + CACAgrad = sigma2CACA*(dijCACA-dk)*expCACA + do l=1,3 + aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA + PgradC(k,l,i) = PgradC(k,l,i)-aux + PgradC(k,l,j) = PgradC(k,l,j)+aux + enddo ! l + enddo ! k + ELSE + rrr = saxs_cutoff*2.0d0/dsqrt(sigma2CACA) + do k=1,nsaxs + dk = distsaxs(k) +c write (2,*) "ijk",i,j,k + sss2 = sscale2(dijCACA,rrr,dk,0.3d0) + if (sss2.eq.0.0d0) cycle + ssgrad2 = sscalgrad2(dijCACA,rrr,dk,0.3d0) + expCACA = dexp(-0.5d0*sigma2CACA*(dijCACA-dk)**2)*sss2 + Pcalc(k) = Pcalc(k)+expCACA +#ifdef DEBUG + write(iout,*) "i j k Pcalc",i,j,Pcalc(k) +#endif + CACAgrad = -sigma2CACA*(dijCACA-dk)*expCACA+ + & ssgrad2*expCACA/sss2 + do l=1,3 +c CA CA + aux = CACAgrad*(C(l,j)-C(l,i))/dijCACA + PgradC(k,l,i) = PgradC(k,l,i)+aux + PgradC(k,l,j) = PgradC(k,l,j)-aux + enddo ! l + enddo ! k + ENDIF +#endif + enddo ! j + enddo ! iint + enddo ! i +#ifdef MPI + if (nfgtasks.gt.1) then + call MPI_Reduce(Pcalc(1),Pcalc_(1),nsaxs,MPI_DOUBLE_PRECISION, + & MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do k=1,nsaxs + Pcalc(k) = Pcalc_(k) + enddo + endif + call MPI_Reduce(PgradC(k,1,1),PgradC_(k,1,1),3*maxsaxs*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do i=1,nres + do l=1,3 + do k=1,nsaxs + PgradC(k,l,i) = PgradC_(k,l,i) + enddo + enddo + enddo + endif +#ifdef ALLSAXS + call MPI_Reduce(PgradX(k,1,1),PgradX_(k,1,1),3*maxsaxs*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + if (fg_rank.eq.king) then + do i=1,nres + do l=1,3 + do k=1,nsaxs + PgradX(k,l,i) = PgradX_(k,l,i) + enddo + enddo + enddo + endif +#endif + endif +#endif +#ifdef MPI + if (fg_rank.eq.king) then +#endif + Cnorm = 0.0d0 + do k=1,nsaxs + Cnorm = Cnorm + Pcalc(k) + enddo + Esaxs_constr = dlog(Cnorm)-wsaxs0 + do k=1,nsaxs + if (Pcalc(k).gt.0.0d0) + & Esaxs_constr = Esaxs_constr - Psaxs(k)*dlog(Pcalc(k)) +#ifdef DEBUG + write (iout,*) "k",k," Esaxs_constr",Esaxs_constr +#endif + enddo +#ifdef DEBUG + write (iout,*) "Cnorm",Cnorm," Esaxs_constr",Esaxs_constr +#endif + do i=nnt,nct + do l=1,3 + auxC=0.0d0 + auxC1=0.0d0 + auxX=0.0d0 + auxX1=0.d0 + do k=1,nsaxs + if (Pcalc(k).gt.0) + & auxC = auxC +Psaxs(k)*PgradC(k,l,i)/Pcalc(k) + auxC1 = auxC1+PgradC(k,l,i) +#ifdef ALLSAXS + auxX = auxX +Psaxs(k)*PgradX(k,l,i)/Pcalc(k) + auxX1 = auxX1+PgradX(k,l,i) +#endif + enddo + gsaxsC(l,i) = auxC - auxC1/Cnorm +#ifdef ALLSAXS + gsaxsX(l,i) = auxX - auxX1/Cnorm +#endif +c write (iout,*) "l i",l,i," gradC",wsaxs*(auxC - auxC1/Cnorm), +c * " gradX",wsaxs*(auxX - auxX1/Cnorm) + enddo + enddo +#ifdef MPI + endif +#endif + return + end +c---------------------------------------------------------------------------- + subroutine e_saxsC(Esaxs_constr) + implicit none + include 'DIMENSIONS' + 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