X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC-NEWCORR%2Freadpdb.f;fp=source%2Fwham%2Fsrc-NEWSC-NEWCORR%2Freadpdb.f;h=0b824766544120991e491ed45bf8e69de10f51f5;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/wham/src-NEWSC-NEWCORR/readpdb.f b/source/wham/src-NEWSC-NEWCORR/readpdb.f new file mode 100644 index 0000000..0b82476 --- /dev/null +++ b/source/wham/src-NEWSC-NEWCORR/readpdb.f @@ -0,0 +1,219 @@ + subroutine readpdb +C Read the PDB file and convert the peptide geometry into virtual-chain +C geometry. + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.CONTROL' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + character*3 seq,atom,res + character*80 card + double precision sccor(3,20) + integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old + double precision dcj + integer rescode + ibeg=1 + ishift1=0 + do i=1,10000 + read (ipdbin,'(a80)',end=10) card + if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10 +C Fish out the ATOM cards. + if (index(card(1:4),'ATOM').gt.0) then + read (card(14:16),'(a3)') atom + if (atom.eq.'CA' .or. atom.eq.'CH3') then +C Calculate the CM of the preceding residue. + if (ibeg.eq.0) call sccenter(ires,iii,sccor) +C Start new residue. + ires_old=ires+ishift-ishift1 + read (card(23:26),*) ires +c print *,"ires_old",ires_old," ires",ires + if (card(27:27).eq."A" .or. card(27:27).eq."B") then +c ishift1=ishift1+1 + endif + read (card(18:20),'(a3)') res + if (ibeg.eq.1) then + ishift=ires-1 + if (res.ne.'GLY' .and. res.ne. 'ACE') then + ishift=ishift-1 + itype(1)=21 + endif + ibeg=0 + else + ishift=ishift+ires-ires_old-1 + endif + ires=ires-ishift+ishift1 + if (res.eq.'ACE') then + ity=10 + else + itype(ires)=rescode(ires,res,0) + endif + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + write (iout,'(2i3,2x,a,3f8.3)') + & ires,itype(ires),res,(c(j,ires),j=1,3) + iii=1 + do j=1,3 + sccor(j,iii)=c(j,ires) + enddo +c write (*,*) card(23:27),ires,itype(ires) + else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and. + & atom.ne.'N ' .and. atom.ne.'C ') then + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif + endif + enddo + 10 write (iout,'(a,i5)') ' Nres: ',ires +C Calculate the CM of the last side chain. + call sccenter(ires,iii,sccor) + nres=ires + nsup=nres + nstart_sup=1 + if (itype(nres).ne.10) then + nres=nres+1 + itype(nres)=21 + do j=1,3 + dcj=c(j,nres-2)-c(j,nres-3) + c(j,nres)=c(j,nres-1)+dcj + c(j,2*nres)=c(j,nres) + enddo + endif + do i=2,nres-1 + do j=1,3 + c(j,i+nres)=dc(j,i) + enddo + enddo + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + if (itype(1).eq.21) then + nsup=nsup-1 + nstart_sup=2 + do j=1,3 + dcj=c(j,4)-c(j,3) + c(j,1)=c(j,2)-dcj + c(j,nres+1)=c(j,1) + enddo + endif +C Copy the coordinates to reference coordinates + do i=1,2*nres + do j=1,3 + cref(j,i)=c(j,i) + enddo + enddo +C Calculate internal coordinates. + do ires=1,nres + write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') + & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + enddo + call flush(iout) + call int_from_cart(.true.,.true.) + do i=1,nres + phi_ref(i)=phi(i) + theta_ref(i)=theta(i) + alph_ref(i)=alph(i) + omeg_ref(i)=omeg(i) + enddo + ishift_pdb=ishift + return + end +c--------------------------------------------------------------------------- + subroutine int_from_cart(lside,lprn) + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + character*3 seq,atom,res + character*80 card + double precision sccor(3,20) + integer rescode + double precision dist,alpha,beta,di + integer i,j,iti + logical lside,lprn + if (lprn) then + write (iout,'(/a)') + & 'Internal coordinates calculated from crystal structure.' + if (lside) then + write (iout,'(8a)') ' Res ',' dvb',' Theta', + & ' Phi',' Dsc_id',' Dsc',' Alpha', + & ' Omega' + else + write (iout,'(4a)') ' Res ',' dvb',' Theta', + & ' Phi' + endif + endif + do i=2,nres + iti=itype(i) + write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1) + if (itype(i-1).ne.21 .and. itype(i).ne.21 .and. + & (dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0)) then + write (iout,'(a,i4)') 'Bad Cartesians for residue',i + stop + endif + theta(i+1)=alpha(i-1,i,i+1) + if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) + enddo + if (itype(1).eq.21) then + do j=1,3 + c(j,1)=c(j,2)+(c(j,3)-c(j,4)) + enddo + endif + if (itype(nres).eq.21) then + do j=1,3 + c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3)) + enddo + endif + if (lside) then + do i=2,nres-1 + do j=1,3 + c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)) + enddo + iti=itype(i) + di=dist(i,nres+i) + if (iti.ne.10) then + alph(i)=alpha(nres+i,i,maxres2) + omeg(i)=beta(nres+i,i,maxres2,i+1) + endif + if (lprn) + & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1), + & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di, + & rad2deg*alph(i),rad2deg*omeg(i) + enddo + else if (lprn) then + do i=2,nres + iti=itype(i) + write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1), + & rad2deg*theta(i),rad2deg*phi(i) + enddo + endif + return + end +c--------------------------------------------------------------------------- + subroutine sccenter(ires,nscat,sccor) + implicit none + include 'DIMENSIONS' + include 'COMMON.CHAIN' + integer ires,nscat,i,j + double precision sccor(3,20),sccmj + do j=1,3 + sccmj=0.0D0 + do i=1,nscat + sccmj=sccmj+sccor(j,i) + enddo + dc(j,ires)=sccmj/nscat + enddo + return + end