X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Freadrtns.F;h=6bb388ab9286a16ba7437c3c50415991afad8290;hb=058b4e72d61c149f66570d580bda9b47464d9df6;hp=1680dab476c9d11d3a1cdb2d90ab0b4f2cea15b2;hpb=7006db11d703450fb82ce20cfc53f74614334362;p=unres.git diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 1680dab..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' @@ -91,7 +92,8 @@ C long axis of side chain min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) - call readi(controlcard,'NCUT',ncut,1) + call readi(controlcard,'NCUT',ncut,0) + call readi(controlcard,'NCLUST',nclust,5) call readi(controlcard,'SYM',symetr,1) write (iout,*) 'sym', symetr call readi(controlcard,'NSTART',nstart,0) @@ -102,7 +104,8 @@ C long axis of side chain lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) - call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) + if (ncut.gt.0) + & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) call readi(controlcard,'IOPT',iopt,2) lside = index(controlcard,"SIDE").gt.0 efree = index(controlcard,"EFREE").gt.0 @@ -125,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 @@ -166,8 +173,7 @@ C C Read weights of the subsequent energy terms. call card_concat(weightcard) write(iout,*) weightcard -C call reada(weightcard,'WSC',wsc,1.0d0) - write(iout,*) wsc + call reada(weightcard,'WSC',wsc,1.0d0) call reada(weightcard,'WLONG',wsc,wsc) call reada(weightcard,'WSCP',wscp,1.0d0) call reada(weightcard,'WELEC',welec,1.0D0) @@ -186,6 +192,7 @@ C call reada(weightcard,'WSC',wsc,1.0d0) 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) @@ -414,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) @@ -955,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) @@ -968,6 +1037,7 @@ c====------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'COMMON.NAMES' include 'COMMON.HOMRESTR' c c For new homol impl @@ -988,13 +1058,15 @@ c & sigma_odl_temp(maxres,maxres,max_template) integer ilen external ilen logical lprn - logical unres_pdb + logical unres_pdb,liiflag c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision rescore_tmp,x12,y12,z12,rescore2_tmp, + & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 + double precision, dimension (max_template,maxres) :: rescore3 character*24 tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -1048,8 +1120,6 @@ cd call flush(iout) c c New c - lim_theta=0 - lim_xx=0 c c Reading HM global scores (prob not required) c @@ -1102,7 +1172,7 @@ c tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. - call readpdb + call readpdb_template(k) do i=1,2*nres do j=1,3 crefjlee(j,i)=c(j,i) @@ -1145,14 +1215,18 @@ c c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore open (ientin,file=tpl_k_rescore,status='old') if (nnt.gt.1) rescore(k,1)=0.0d0 - do irec=nnt,maxdim ! loop for reading res sim + do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, - & idomain_tmp + & rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + rescore3(k,i_tmp)=rescore3_tmp + write(iout,'(a7,i5,3f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -1252,7 +1326,8 @@ c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2)*rescore(k,i-3) ! right expression ? c Instead of res sim other local measure of b/b str reliability possible - sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*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)) c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) enddo lim_dih=nct-nnt-2 @@ -1284,14 +1359,14 @@ c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 - sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + if (sigma_theta(k,i).ne.0) + & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif - lim_theta=nct-nnt-1 if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -1317,15 +1392,14 @@ c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i) - sigma_d(k,i)=rescore(k,i) ! right expression ? - sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) + sigma_d(k,i)=rescore3(k,i) ! right expression ? + if (sigma_d(k,i).ne.0) + & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) c read (ientin,*) sigma_d(k,i) ! 1st variant - if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right? enddo - lim_xx=nct-nnt+1 endif enddo c @@ -1334,22 +1408,34 @@ 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) then - do ki=ii,lim_odl-1 - ires_homo(ki)=ires_homo(ki+1) - jres_homo(ki)=jres_homo(ki+1) - ii_in_use(ki)=ii_in_use(ki+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+1) - sigma_odl(k,ki)=sigma_odl(k,ki+1) - l_homo(k,ki)=l_homo(k,ki+1) + 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-1 - lim_odl=lim_odl-1 + ii=ii-iishift + lim_odl=lim_odl-iishift endif enddo enddo @@ -1374,17 +1460,19 @@ cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" - do i=nnt+3,lim_dih - write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i), + do i=nnt+3,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + & (rad2deg*dih(ki,i), & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" - do i=nnt+2,lim_theta - write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i), + do i=nnt+2,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + & (rad2deg*thetatpl(ki,i), & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" - do i=nnt,lim_xx + do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)