From: Cezary Czaplewski Date: Tue, 7 Aug 2018 15:00:12 +0000 (+0200) Subject: NMR log-exp from Adam X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=61fae8832983d4a71eb8efc042d3e0ba5087740c NMR log-exp from Adam --- diff --git a/source/unres/src_MD-M/COMMON.SBRIDGE b/source/unres/src_MD-M/COMMON.SBRIDGE index 469add2..0df27ec 100644 --- a/source/unres/src_MD-M/COMMON.SBRIDGE +++ b/source/unres/src_MD-M/COMMON.SBRIDGE @@ -2,16 +2,24 @@ integer ns,nss,nfree,iss common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss, & ns,nss,nfree,iss(maxss) - double precision dhpb,dhpb1,forcon,fordepth,xlscore,wboltzd - integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb,irestr_type + double precision dhpb,dhpb1,forcon,fordepth,xlscore,wboltzd, + & dhpb_peak,dhpb1_peak,forcon_peak,fordepth_peak,scal_peak + integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb,ibecarb_peak,npeak, + & ipeak, + & irestr_type,irestr_type_peak,ihpb_peak,jhpb_peak,nhpb_peak common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim), & fordepth(maxdim),xlscore(maxdim),wboltzd, & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),irestr_type(maxdim), & nhpb + common /NMRpeaks/ dhpb_peak(maxdim),dhpb1_peak(maxdim), + & forcon_peak(maxdim),fordepth_peak(maxdim),scal_peak, + & ihpb_peak(maxdim),jhpb_peak(maxdim),ibecarb_peak(maxdim), + & irestr_type_peak(maxdim),ipeak(2,maxdim),npeak,nhpb_peak double precision weidis common /restraints/ weidis - integer link_start,link_end - common /links_split/ link_start,link_end + integer link_start,link_end,link_start_peak,link_end_peak + common /links_split/ link_start,link_end,link_start_peak, + & link_end_peak double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss logical dyn_ss,dyn_ss_mask common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres), 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 badacd2..332f1ee 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -5192,7 +5192,7 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' - dimension ggg(3) + dimension ggg(3),ggg_peak(3,20) ehpb=0.0D0 do i=1,3 ggg(i)=0.0d0 @@ -5200,8 +5200,67 @@ C 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 +c & " constr_dist",constr_dist," link_start_peak",link_start_peak, +c & " link_end_peak",link_end_peak + if (link_end.eq.0.and.link_end_peak.eq.0) return + do i=link_start_peak,link_end_peak + ehpb_peak=0.0d0 +c print *,"i",i," link_end_peak",link_end_peak," ipeak", +c & ipeak(1,i),ipeak(2,i) + do ip=ipeak(1,i),ipeak(2,i) + ii=ihpb_peak(ip) + jj=jhpb_peak(ip) + dd=dist(ii,jj) + iip=ip-ipeak(1,i)+1 +C iii and jjj point to the residues for which the distance is assigned. + if (ii.gt.nres) then + iii=ii-nres + jjj=jj-nres + else + iii=ii + jjj=jj + endif + aux=rlornmr1(dd,dhpb_peak(ip),dhpb1_peak(ip),forcon_peak(ip)) + aux=dexp(-scal_peak*aux) + ehpb_peak=ehpb_peak+aux + fac=rlornmr1prim(dd,dhpb_peak(ip),dhpb1_peak(ip), + & forcon_peak(ip))*aux/dd + do j=1,3 + ggg_peak(j,iip)=fac*(c(j,jj)-c(j,ii)) + enddo + if (energy_dec) write (iout,'(a6,3i5,6f10.3,i5)') + & "edisL",i,ii,jj,dd,dhpb_peak(ip),dhpb1_peak(ip), + & forcon_peak(ip),fordepth_peak(ip),ehpb_peak + enddo +c write (iout,*) "ehpb_peak",ehpb_peak," scal_peak",scal_peak + ehpb=ehpb-fordepth_peak(ipeak(1,i))*dlog(ehpb_peak)/scal_peak + do ip=ipeak(1,i),ipeak(2,i) + iip=ip-ipeak(1,i)+1 + do j=1,3 + ggg(j)=ggg_peak(j,iip)/ehpb_peak + enddo + ii=ihpb_peak(ip) + jj=jhpb_peak(ip) +C iii and jjj point to the residues for which the distance is assigned. + if (ii.gt.nres) then + iii=ii-nres + jjj=jj-nres + else + iii=ii + jjj=jj + endif + 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,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo + enddo + enddo do i=link_start,link_end C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a C CA-CA distance used in regularization of structure. diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index d789d82..11a8267 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -1526,3 +1526,25 @@ cd & " lim_dih",lim_dih #endif return end +c------------------------------------------------------------------------------ + subroutine NMRpeak_partition + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SBRIDGE' + include 'COMMON.IOUNITS' + include 'COMMON.SETUP' +#ifdef MPI + call int_bounds(npeak,link_start_peak,link_end_peak) + write (iout,*) 'Processor',fg_rank,' CG group',kolor, + & ' absolute rank',MyRank, + & ' npeak',npeak,' link_start_peak=',link_start_peak, + & ' link_end_peak',link_end_peak +#else + link_start_peak=1 + link_end_peak=npeak +#endif + return + end diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index add54af..a94bb53 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -53,6 +53,24 @@ C Print restraint information & irestr_type(i) enddo endif + if (npeak.gt.0) then + write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)') + & "The following",npeak, + & " NMR peak restraints have been imposed:", + & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V", + & " score"," type"," ipeak" + do i=1,npeak + do j=ipeak(1,i),ipeak(2,i) + write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j), + & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j), + & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j) + enddo + enddo + write (iout,*) "The ipeak array" + do i=1,npeak + write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i) + enddo + endif #ifdef MPI endif #endif @@ -998,11 +1016,14 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) enddo call contact(.true.,ncont_ref,icont_ref,co) endif -c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup + endif +c print *, "A TU" + write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr -c write (iout,*) "After read_dist_constr nhpb",nhpb + write (iout,*) "After read_dist_constr nhpb",nhpb if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp + call NMRpeak_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co if (pdbref) then @@ -1018,7 +1039,6 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i) enddo endif - endif write (iout,*) "calling read_saxs_consrtr",nsaxs if (nsaxs.gt.0) call read_saxs_constr @@ -2548,9 +2568,9 @@ c------------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.SBRIDGE' - integer ifrag_(2,100),ipair_(2,100) - double precision wfrag_(100),wpair_(100) - character*500 controlcard + integer ifrag_(2,100),ipair_(2,1000) + double precision wfrag_(100),wpair_(1000) + character*5000 controlcard logical normalize,next integer restr_type double precision xlink(4,0:4) / @@ -2566,6 +2586,10 @@ c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) next=.true. + npeak=0 + ipeak=0 + nhpb_peak=0 + DO WHILE (next) call card_concat(controlcard) @@ -2579,21 +2603,24 @@ c call flush(iout) call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) if (restr_type.eq.10) & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0) + if (restr_type.eq.12) + & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0) call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0) call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) normalize = index(controlcard,"NORMALIZE").gt.0 write (iout,*) "WBOLTZD",wboltzd -c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ -c write (iout,*) "IFRAG" -c do i=1,nfrag_ -c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) -c enddo -c write (iout,*) "IPAIR" -c do i=1,npair_ -c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) -c enddo + write (iout,*) "SCAL_PEAK",scal_peak + write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ + write (iout,*) "IFRAG" + do i=1,nfrag_ + write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + enddo + write (iout,*) "IPAIR" + do i=1,npair_ + write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) + enddo if (nfrag_.gt.0) write (iout,*) & "Distance restraints as generated from reference structure" do i=1,nfrag_ @@ -2653,13 +2680,14 @@ c write (iout,*) "j",j," k",k endif do j=ifrag_(1,ii),ifrag_(2,ii) do k=ifrag_(1,jj),ifrag_(2,jj) + ddjk=dist(j,k) if (restr_type.eq.1) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk - forcon(nhpb)=wfrag_(i) + forcon(nhpb)=wpair_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 @@ -2667,7 +2695,7 @@ c write (iout,*) "j",j," k",k ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk - forcon(nhpb)=wfrag_(i) + forcon(nhpb)=wpair_(i) endif else if (restr_type.eq.3) then nhpb=nhpb+1 @@ -2675,14 +2703,14 @@ c write (iout,*) "j",j," k",k ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk - forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) + forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", + & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", + write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo @@ -2692,9 +2720,42 @@ c write (iout,*) "j",j," k",k c print *,ndist_ write (iout,*) "Distance restraints as read from input" do i=1,ndist_ - if (restr_type.eq.11) then + if (restr_type.eq.12) then + read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), + & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), + & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), + & fordepth_peak(nhpb_peak+1),npeak +c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), +c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), +c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), +c & fordepth_peak(nhpb_peak+1),npeak + if (forcon_peak(nhpb_peak+1).le.0.0d0.or. + & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle + nhpb_peak=nhpb_peak+1 + irestr_type_peak(nhpb_peak)=12 + if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i + ipeak(2,npeak)=i +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), + & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), + & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), + & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) +#else + write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), + & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), + & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), + & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) +#endif + if (ibecarb_peak(nhpb_peak).gt.0) then + ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres + jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres + endif + else if (restr_type.eq.11) then read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), - & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle nhpb=nhpb+1 @@ -2713,7 +2774,7 @@ c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif - else if (constr_dist.eq.10) then + else if (restr_type.eq.10) then c Cross-lonk Markov-like potential call card_concat(controlcard) call readi(controlcard,"ILINK",ihpb(nhpb+1),0)