X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fmolread_zs.F;fp=source%2Fwham%2Fsrc-NEWSC%2Fmolread_zs.F;h=431680df45ddbd689bd5060d605f932c11e881ed;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/wham/src-NEWSC/molread_zs.F b/source/wham/src-NEWSC/molread_zs.F new file mode 100755 index 0000000..431680d --- /dev/null +++ b/source/wham/src-NEWSC/molread_zs.F @@ -0,0 +1,378 @@ + subroutine molread(*) +C +C Read molecular data. +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.INTERACT' + include 'COMMON.LOCAL' + include 'COMMON.NAMES' + include 'COMMON.CHAIN' + include 'COMMON.FFIELD' + include 'COMMON.SBRIDGE' + include 'COMMON.TORCNSTR' + include 'COMMON.CONTROL' + character*4 sequence(maxres) + integer rescode + double precision x(maxvar) + character*320 controlcard,ucase + dimension itype_pdb(maxres) + logical seq_comp + call card_concat(controlcard,.true.) + call reada(controlcard,'SCAL14',scal14,0.4d0) + call reada(controlcard,'SCALSCP',scalscp,1.0d0) + call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0) + call reada(controlcard,'DELT_CORR',delt_corr,0.5d0) + r0_corr=cutoff_corr-delt_corr + call readi(controlcard,"NRES",nres,0) + iscode=index(controlcard,"ONE_LETTER") + if (nres.le.0) then + write (iout,*) "Error: no residues in molecule" + return1 + endif + if (nres.gt.maxres) then + write (iout,*) "Error: too many residues",nres,maxres + endif + write(iout,*) 'nres=',nres +C Read sequence of the protein + if (iscode.gt.0) then + read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) + else + read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) + endif +C Convert sequence to numeric code + do i=1,nres + itype(i)=rescode(i,sequence(i),iscode) + enddo + write (iout,*) "Numeric code:" + write (iout,'(20i4)') (itype(i),i=1,nres) + do i=1,nres-1 +#ifdef PROCOR + if (itype(i).eq.21 .or. itype(i+1).eq.21) then +#else + if (itype(i).eq.21) then +#endif + itel(i)=0 +#ifdef PROCOR + else if (itype(i+1).ne.20) then +#else + else if (itype(i).ne.20) then +#endif + itel(i)=1 + else + itel(i)=2 + endif + enddo + call read_bridge + + if (with_dihed_constr) then + + read (inp,*) ndih_constr + if (ndih_constr.gt.0) then + read (inp,*) ftors + write (iout,*) 'FTORS',ftors + 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 + + nnt=1 + nct=nres + if (itype(1).eq.21) nnt=2 + if (itype(nres).eq.21) nct=nct-1 + write(iout,*) 'NNT=',NNT,' NCT=',NCT +c Read distance restraints + if (constr_dist.gt.0) then + if (refstr) call read_ref_structure(*11) + call read_dist_constr + call hpb_partition + endif + + call setup_var + call init_int_table + if (ns.gt.0) then + write (iout,'(/a,i3,a)') 'The chain contains',ns, + & ' disulfide-bridging cysteines.' + write (iout,'(20i4)') (iss(i),i=1,ns) + 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 + endif + write (iout,'(a)') + return + 11 stop "Error reading reference structure" + end +c----------------------------------------------------------------------------- + logical function seq_comp(itypea,itypeb,length) + implicit none + integer length,itypea(length),itypeb(length) + integer i + do i=1,length + if (itypea(i).ne.itypeb(i)) then + seq_comp=.false. + return + endif + enddo + seq_comp=.true. + return + end +c----------------------------------------------------------------------------- + subroutine read_bridge +C Read information about disulfide bridges. + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.CHAIN' + include 'COMMON.FFIELD' + include 'COMMON.SBRIDGE' +C Read bridging residues. + read (inp,*) ns,(iss(i),i=1,ns) + print *,'ns=',ns + write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns) +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, + & ' can form a disulfide bridge?!!!' + write (*,'(2a,i3,a)') + & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & ' can form a disulfide bridge?!!!' + stop + endif + enddo +C Read preformed bridges. + if (ns.gt.0) then + read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) + write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) + if (nss.gt.0) then + nhpb=nss +C Check if the residues involved in bridges are in the specified list of +C bridging residues. + do i=1,nss + do j=1,i-1 + if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) + & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then + write (iout,'(a,i3,a)') 'Disulfide pair',i, + & ' contains residues present in other pairs.' + write (*,'(a,i3,a)') 'Disulfide pair',i, + & ' contains residues present in other pairs.' + stop + endif + enddo + do j=1,ns + if (ihpb(i).eq.iss(j)) goto 10 + enddo + write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' + 10 continue + do j=1,ns + if (jhpb(i).eq.iss(j)) goto 20 + enddo + write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' + 20 continue + dhpb(i)=dbr + forcon(i)=fbr + enddo + do i=1,nss + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + enddo + endif + endif + return + end +c------------------------------------------------------------------------------ + subroutine read_angles(kanal,iscor,energ,iprot,*) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.INTERACT' + include 'COMMON.SBRIDGE' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + character*80 lineh + read(kanal,'(a80)',end=10,err=10) lineh + read(lineh(:5),*,err=8) ic + read(lineh(6:),*,err=8) energ + goto 9 + 8 ic=1 + print *,'error, assuming e=1d10',lineh + energ=1d10 + nss=0 + 9 continue + read(lineh(18:),*,end=10,err=10) nss + IF (NSS.LT.9) THEN + read (lineh(20:),*,end=10,err=10) + & (IHPB(I),JHPB(I),I=1,NSS),iscor + ELSE + read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8) + read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I), + & I=9,NSS),iscor + ENDIF +c print *,"energy",energ," iscor",iscor + read (kanal,*,err=10,end=10) (theta(i),i=3,nres) + read (kanal,*,err=10,end=10) (phi(i),i=4,nres) + read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1) + read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1) + do i=1,nres + theta(i)=deg2rad*theta(i) + phi(i)=deg2rad*phi(i) + alph(i)=deg2rad*alph(i) + omeg(i)=deg2rad*omeg(i) + enddo + return + 10 return1 + end +c------------------------------------------------------------------------------- + subroutine read_dist_constr + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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 +c write (iout,*) "Calling read_dist_constr" +c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +c call flush(iout) + call card_concat(controlcard,.true.) + 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) + if (.not.refstr .and. nfrag_.gt.0) then + write (iout,*) + & "ERROR: no reference structure to compute distance restraints" + write (iout,*) + & "Restraints must be specified explicitly (NDIST=number)" + stop + endif + if (nfrag_.lt.2 .and. npair_.gt.0) then + write (iout,*) "ERROR: Less than 2 fragments specified", + & " but distance restraints between pairs requested" + stop + endif + 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 + 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_ + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) + 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)) + endif + enddo + do i=1,nhpb + write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", + & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) + enddo + call flush(iout) + return + end