X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Freadrtns.F;h=06777ec96681579c450704d1d1303b25b341cea2;hb=c9dddd7cdc2cbf55f6fa8df9404eea7ecfe59605;hp=1b4c91180d5f44fc83604cff958d30e65bc69680;hpb=3dbe5ceeea4b858bc25fa1469237e697c0bf293f;p=unres.git diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 1b4c911..06777ec 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -19,7 +19,7 @@ C #ifdef MPL include 'COMMON.INFO' #endif - integer i + integer i,i1,i2,it1,it2 read (INP,'(a80)') titel call card_concat(controlcard) @@ -35,6 +35,11 @@ C pdbref=(index(controlcard,'PDBREF').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) + with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 + call readi(controlcard,'CONSTR_DIST',constr_dist,0) + write (iout,*) "with_dihed_constr ",with_dihed_constr, + & " CONSTR_DIST",constr_dist + call flush(iout) min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) @@ -87,6 +92,7 @@ C include 'COMMON.CONTROL' include 'COMMON.CONTACTS' include 'COMMON.TIME1' + include 'COMMON.TORCNSTR' #ifdef MPL include 'COMMON.INFO' #endif @@ -96,7 +102,7 @@ C double precision x(maxvar) integer itype_pdb(maxres) logical seq_comp - integer i,j + integer i,j,kkk,i1,i2,it1,it2 C C Body C @@ -115,6 +121,7 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'WTURN4',wturn4,1.0D0) call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) + call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) @@ -125,6 +132,51 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) if (index(weightcard,'SOFT').gt.0) ipot=6 + call reada(weightcard,"D0CM",d0cm,3.78d0) + call reada(weightcard,"AKCM",akcm,15.1d0) + call reada(weightcard,"AKTH",akth,11.0d0) + call reada(weightcard,"AKCT",akct,12.0d0) + call reada(weightcard,"V1SS",v1ss,-1.08d0) + call reada(weightcard,"V2SS",v2ss,7.61d0) + call reada(weightcard,"V3SS",v3ss,13.7d0) + call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,"ATRISS",atriss,0.301D0) + call reada(weightcard,"BTRISS",btriss,0.021D0) + call reada(weightcard,"CTRISS",ctriss,1.001D0) + call reada(weightcard,"DTRISS",dtriss,1.001D0) + write (iout,*) "ATRISS=", atriss + write (iout,*) "BTRISS=", btriss + write (iout,*) "CTRISS=", ctriss + write (iout,*) "DTRISS=", dtriss + dyn_ss=(index(weightcard,'DYN_SS').gt.0) + do i=1,maxres + dyn_ss_mask(i)=.false. + enddo + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(weightcard,"HT",Ht,0.0D0) + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) + Ht=Ht/wsc-0.25*eps(1,1) + akcm=akcm*wstrain/wsc + akth=akth*wstrain/wsc + akct=akct*wstrain/wsc + v1ss=v1ss*wstrain/wsc + v2ss=v2ss*wstrain/wsc + v3ss=v3ss*wstrain/wsc + else + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + endif + write (iout,'(/a)') "Disulfide bridge parameters:" + write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr + write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm + write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct + write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, + & ' v3ss:',v3ss + C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 @@ -146,9 +198,10 @@ C 12/1/95 Added weight for the multi-body term WCORR weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 + weights(19)=wsccor write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3, - & wturn4,wturn6 + & wturn4,wturn6,wsccor 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -166,7 +219,9 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)') + if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' @@ -229,6 +284,26 @@ C Convert sequence to numeric code print *,'Call Read_Bridge.' call read_bridge +C this fragment reads diheadral constrains + if (with_dihed_constr) then + + read (inp,*) ndih_constr + if (ndih_constr.gt.0) then + read (inp,*) ftors + write (iout,*) 'FTORS',ftors +C ftors is the force constant for torsional quartic constrains + read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) + write (iout,*) + & 'There are',ndih_constr,' constraints on phi angles.' + do i=1,ndih_constr + write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) + enddo + do i=1,ndih_constr + phi0(i)=deg2rad*phi0(i) + drange(i)=deg2rad*drange(i) + enddo + endif ! endif ndif_constr.gt.0 + endif ! with_dihed_constr nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT @@ -299,14 +374,54 @@ c endif nstart_sup=nnt nstart_seq=nnt nsup=nct-nnt+1 + kkk=1 do i=1,2*nres do j=1,3 - cref(j,i)=c(j,i) + cref(j,i,kkk)=c(j,i) enddo enddo endif call contact(.true.,ncont_ref,icont_ref) endif + if (ns.gt.0) then +C write (iout,'(/a,i3,a)') +C & 'The chain contains',ns,' disulfide-bridging cysteines.' + write (iout,'(20i4)') (iss(i),i=1,ns) + if (dyn_ss) then + write(iout,*)"Running with dynamic disulfide-bond formation" + else + write (iout,'(/a/)') 'Pre-formed links are:' + do i=1,nss + i1=ihpb(i)-nres + i2=jhpb(i)-nres + it1=itype(i1) + it2=itype(i2) + write (iout,'(2a,i3,3a,i3,a,3f10.3)') + & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i), + & ebr,forcon(i) + enddo + write (iout,'(a)') + endif + endif + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + ihpb(i-nss)=ihpb(i) + jhpb(i-nss)=jhpb(i) + forcon(i-nss)=forcon(i) + dhpb(i-nss)=dhpb(i) + enddo + nhpb=nhpb-nss + nss=0 + call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. + enddo + endif +c Read distance restraints + if (constr_dist.gt.0) then + call read_dist_constr + call hpb_partition + endif return end c----------------------------------------------------------------------------- @@ -351,10 +466,12 @@ C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then write (iout,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' #ifdef MPL call mp_stopall(error_msg) @@ -395,8 +512,8 @@ C bridging residues. enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 20 continue - dhpb(i)=dbr - forcon(i)=fbr +C dhpb(i)=dbr +C forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres @@ -477,6 +594,25 @@ c---------------------------------------------------------------------------- read (rekord(iread:),*) wartosc return end +C---------------------------------------------------------------------- + subroutine multreadi(rekord,lancuch,tablica,dim,default) + implicit none + integer dim,i + integer tablica(dim),default + character*(*) rekord,lancuch + character*80 aux + integer ilen,iread + external ilen + do i=1,dim + tablica(i)=default + enddo + iread=index(rekord,lancuch(:ilen(lancuch))//"=") + if (iread.eq.0) return + iread=iread+ilen(lancuch)+1 + read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) + 10 return + end + c---------------------------------------------------------------------------- subroutine card_concat(card) include 'DIMENSIONS' @@ -578,3 +714,131 @@ C #endif return end + subroutine read_dist_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' + integer ifrag_(2,100),ipair_(2,100) + double precision wfrag_(100),wpair_(100) + character*500 controlcard + logical lprn /.true./ + write (iout,*) "Calling read_dist_constr" +C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +C call flush(iout) + write(iout,*) "TU sie wywalam?" + call card_concat(controlcard) + write (iout,*) controlcard + call flush(iout) + call readi(controlcard,"NFRAG",nfrag_,0) + call readi(controlcard,"NPAIR",npair_,0) + call readi(controlcard,"NDIST",ndist_,0) + call reada(controlcard,'DIST_CUT',dist_cut,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) + 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 + call flush(iout) + do i=1,nfrag_ + if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup + if (ifrag_(2,i).gt.nstart_sup+nsup-1) + & ifrag_(2,i)=nstart_sup+nsup-1 +c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + call flush(iout) + if (wfrag_(i).gt.0.0d0) then + do j=ifrag_(1,i),ifrag_(2,i)-1 + do k=j+1,ifrag_(2,i) + write (iout,*) "j",j," k",k + ddjk=dist(j,k) + if (constr_dist.eq.1) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + else if (constr_dist.eq.2) then + if (ddjk.le.dist_cut) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + endif + else + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) + endif + if (lprn) + & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + do i=1,npair_ + if (wpair_(i).gt.0.0d0) then + ii = ipair_(1,i) + jj = ipair_(2,i) + if (ii.gt.jj) then + itemp=ii + ii=jj + jj=itemp + endif + do j=ifrag_(1,ii),ifrag_(2,ii) + do k=ifrag_(1,jj),ifrag_(2,jj) + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + forcon(nhpb)=wpair_(i) + dhpb(nhpb)=dist(j,k) + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + do i=1,ndist_ + if (constr_dist.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) +C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", +C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + else + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + endif + if (forcon(nhpb+1).gt.0.0d0) then + nhpb=nhpb+1 + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + endif +C endif + enddo + call hpb_partition + call flush(iout) + return + end