X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=4c48bc0ebc7bcf62b84632ad8f398b20e235b4ad;hb=04e787c1d9ff8dde198693024cac7de49d1435b1;hp=c511fb5fd7c7d43f84b6194b41aed29a64ff0a66;hpb=8990ec14fa8e20cd2830e9e482631510fd058dcd;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index c511fb5..4c48bc0 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 @@ -112,6 +130,7 @@ C Set up the time limit (caution! The time must be input in minutes!) write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE", & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) + read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0 call readi(controlcard,'SYM',symetr,1) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 @@ -997,10 +1016,12 @@ 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 if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co @@ -1017,11 +1038,9 @@ 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 - - + if (constr_homology.gt.0) then call read_constr_homology if (indpdb.gt.0 .or. pdbref) then @@ -1068,6 +1087,7 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb if (nhpb.gt.0) call hpb_partition + if (peak.gt.0) call NMRpeak_partition c write (iout,*) "After read_dist_constr nhpb",nhpb c call flush(iout) if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 @@ -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) @@ -2868,6 +2929,7 @@ c Alternative: reading from input call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) + dist1cut=(index(controlcard,'DIST1CUT').gt.0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0) @@ -2922,6 +2984,10 @@ c enddo enddo + if (read_homol_frag) then + call read_klapaucjusz + else + do k=1,constr_homology read(inp,'(a)') pdbfile @@ -3001,6 +3067,19 @@ c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut + if (dist1cut .and. k.gt.1) then + ii=ii+1 + if (l_homo(1,ii)) then + ii_in_use(ii)=1 + l_homo(k,ii)=.true. + ires_homo(ii)=i + jres_homo(ii)=j + odl(k,ii)=distal + sigma_odl(k,ii)=sigma_odl(1,ii) + else + l_homo(k,ii)=.false. + endif + else if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then @@ -3040,6 +3119,7 @@ c & constr_homology ii=ii+1 l_homo(k,ii)=.false. endif + endif enddo enddo lim_odl=ii @@ -3192,6 +3272,9 @@ c enddo enddo endif + + endif ! .not. klapaucjusz + if (constr_homology.gt.0) call homology_partition if (constr_homology.gt.0) call init_int_table c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end @@ -3367,3 +3450,257 @@ c r1 = prng_next(me) #endif return end + +c ----------------------------------------------------------------- + subroutine read_klapaucjusz + + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.GEO' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + character*256 fragfile + integer ninclust(maxclust),inclust(max_template,maxclust), + & nresclust(maxclust),iresclust(maxres,maxclust) + + character*2 kic2 + character*24 model_ki_dist, model_ki_angle + character*500 controlcard + integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp + logical lprn /.true./ + integer ilen + external ilen + logical liiflag +c +c + double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision, dimension (max_template,maxres) :: rescore + double precision, dimension (max_template,maxres) :: rescore2 + character*24 pdbfile,tpl_k_rescore + +c +c For new homol impl +c + include 'COMMON.VAR' +c + call getenv("FRAGFILE",fragfile) + open(ientin,file=fragfile,status="old",err=10) + read(ientin,*) constr_homology,nclust + l_homo = .false. + sigma_theta=0.0 + sigma_d=0.0 + sigma_dih=0.0 +c Read pdb files + do k=1,constr_homology + read(ientin,'(a)') pdbfile + write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', + & pdbfile(:ilen(pdbfile)) + open(ipdbin,file=pdbfile,status='old',err=33) + goto 34 + 33 write (iout,'(a,5x,a)') 'Error opening PDB file', + & pdbfile(:ilen(pdbfile)) + stop + 34 continue + unres_pdb=.false. + call readpdb_template(k) + do i=1,nres + rescore(k,i)=0.2d0 + rescore2(k,i)=1.0d0 + enddo + enddo +c Read clusters + do i=1,nclust + read(ientin,*) ninclust(i),nresclust(i) + read(ientin,*) (inclust(k,i),k=1,ninclust(i)) + read(ientin,*) (iresclust(k,i),k=1,nresclust(i)) + enddo +c +c Loop over clusters +c + do l=1,nclust + do ll = 1,ninclust(l) + + k = inclust(ll,l) + do i=1,nres + idomain(k,i)=0 + enddo + do i=1,nresclust(l) + if (nnt.gt.1) then + idomain(k,iresclust(i,l)+1) = 1 + else + idomain(k,iresclust(i,l)) = 1 + endif + enddo +c +c Distance restraints +c +c ... --> odl(k,ii) +C Copy the coordinates from reference coordinates (?) + do i=1,2*nres + do j=1,3 + c(j,i)=chomo(j,i,k) +c write (iout,*) "c(",j,i,") =",c(j,i) + enddo + enddo + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + if (waga_dist.ne.0.0d0) then + ii=0 + do i = nnt,nct-2 + do j=i+2,nct + + x12=c(1,i)-c(1,j) + y12=c(2,i)-c(2,j) + z12=c(3,i)-c(3,j) + distal=dsqrt(x12*x12+y12*y12+z12*z12) +c write (iout,*) k,i,j,distal,dist2_cut + + if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 + & .and. distal.le.dist2_cut ) then + + ii=ii+1 + ii_in_use(ii)=1 + l_homo(k,ii)=.true. + +c write (iout,*) "k",k +c write (iout,*) "i",i," j",j," constr_homology", +c & constr_homology + ires_homo(ii)=i + jres_homo(ii)=j + odl(k,ii)=distal + if (read2sigma) then + sigma_odl(k,ii)=0 + do ik=i,j + sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) + enddo + sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) + if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = + & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) + else + if (odl(k,ii).le.dist_cut) then + sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) + else +#ifdef OLDSIGMA + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) +#else + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) +#endif + endif + endif + sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) + else + ii=ii+1 +c l_homo(k,ii)=.false. + endif + enddo + enddo + lim_odl=ii + endif +c +c Theta, dihedral and SC retraints +c + if (waga_angle.gt.0.0d0) then + do i = nnt+3,nct + if (idomain(k,i).eq.0) then +c sigma_dih(k,i)=0.0 + cycle + endif + dih(k,i)=phiref(i) + sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2)+rescore(k,i-3))/4.0 +c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), +c & " sigma_dihed",sigma_dih(k,i) + if (sigma_dih(k,i).ne.0) + & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) + enddo + lim_dih=nct-nnt-2 + endif + + if (waga_theta.gt.0.0d0) then + do i = nnt+2,nct + if (idomain(k,i).eq.0) then +c sigma_theta(k,i)=0.0 + cycle + endif + thetatpl(k,i)=thetaref(i) + sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2))/3.0 + if (sigma_theta(k,i).ne.0) + & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + enddo + endif + + if (waga_d.gt.0.0d0) then + do i = nnt,nct + if (itype(i).eq.10) cycle + if (idomain(k,i).eq.0 ) then +c sigma_d(k,i)=0.0 + cycle + endif + xxtpl(k,i)=xxref(i) + yytpl(k,i)=yyref(i) + zztpl(k,i)=zzref(i) + sigma_d(k,i)=rescore(k,i) + if (sigma_d(k,i).ne.0) + & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) + if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 + enddo + endif + enddo ! l + enddo ! ll +c +c remove distance restraints not used in any model from the list +c shift data in all arrays +c + if (waga_dist.ne.0.0d0) then + ii=0 + liiflag=.true. + do i=nnt,nct-2 + do j=i+2,nct + ii=ii+1 + if (ii_in_use(ii).eq.0.and.liiflag) then + liiflag=.false. + iistart=ii + endif + if (ii_in_use(ii).ne.0.and..not.liiflag.or. + & .not.liiflag.and.ii.eq.lim_odl) then + if (ii.eq.lim_odl) then + iishift=ii-iistart+1 + else + iishift=ii-iistart + endif + liiflag=.true. + do ki=iistart,lim_odl-iishift + ires_homo(ki)=ires_homo(ki+iishift) + jres_homo(ki)=jres_homo(ki+iishift) + ii_in_use(ki)=ii_in_use(ki+iishift) + do k=1,constr_homology + odl(k,ki)=odl(k,ki+iishift) + sigma_odl(k,ki)=sigma_odl(k,ki+iishift) + l_homo(k,ki)=l_homo(k,ki+iishift) + enddo + enddo + ii=ii-iishift + lim_odl=lim_odl-iishift + endif + enddo + enddo + endif + + return + 10 stop "Error infragment file" + end +c----------------------------------------------------------------------