source/maxlik/src-Fmatch_safe/readrtns_compar.F

   1       subroutine read_ref_structure(iprot,*)
   2 C
   3 C Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
   4 C angles.
   5 C
   6       implicit none
   7       include 'DIMENSIONS'
   8       include 'DIMENSIONS.ZSCOPT'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.GEO'
  11       include 'COMMON.VAR'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.NAMES'
  15       include 'COMMON.CHAIN'
  16       include 'COMMON.FFIELD'
  17       include 'COMMON.SBRIDGE'
  18       include 'COMMON.HEADER'
  19       include 'COMMON.CONTROL'
  20       include 'COMMON.TIME1'
  21       include 'COMMON.COMPAR'
  22       include 'COMMON.VMCPAR'
  23       include 'COMMON.EREF'
  24       character*4 sequence(maxres)
  25       integer rescode
  26       double precision x(maxvar)
  27       integer itype_pdb(maxres)
  28       logical seq_comp
  29       integer i,j,k,ib,iref,iprot,nres_pdb
  30       double precision ddsc,dist
  31       double precision efree_temp
  32       integer ilen
  33       external ilen
  34 C
  35       write (2,*) "PDBREF ",pdbref
  36       nres0=nres
  37       DO IB=1,NBETA(1,iprot)
  38       if (pdbref) then
  39         iref=0
  40         write (iout,'(2a,1h.)') 'PDB data will be read from file ',
  41      &    pdbfile(ib,iprot)(:ilen(pdbfile(ib,iprot)))
  42         open(ipdbin,file=pdbfile(ib,iprot),status='old',err=33)
  43         goto 34
  44   33    write (iout,'(a)') 'Error opening PDB file.'
  45         return1
  46   34    continue
  47         nres_pdb=nres
  48         DO IREF=1,MAXREF
  49         write (2,*) "Reading reference structure",iref
  50         call readpdb(print_refstr,iprot,efree_temp,*44)
  51         write (2,*) "After readpdb"
  52         call flush(iout)
  53         if (nres_pdb.eq.0) then
  54           exit
  55         endif
  56 C Copy the coordinates to reference coordinates
  57         efreeref(iref,ib,iprot)=efree_temp
  58         do i=1,2*nres
  59           do j=1,3
  60             cref(j,i,iref,ib,iprot)=c(j,i)
  61           enddo
  62         enddo
  63         do i=1,nres
  64           phi_ref(i,iref,ib,iprot)=phi(i)
  65           theta_ref(i,iref,ib,iprot)=theta(i)
  66           xxref(i,iref,ib,iprot)=xxtab(i)
  67           yyref(i,iref,ib,iprot)=yytab(i)
  68           zzref(i,iref,ib,iprot)=zztab(i)
  69         enddo
  70         do i=1,nres
  71           itype_pdb(i)=itype(i)
  72         enddo
  73         nres_pdb=nres
  74         nres=nres0
  75         nstart_seq=nnt
  76         if (nsup(iprot).le.(nct-nnt+1)) then
  77           do i=0,nct-nnt+1-nsup(iprot)
  78             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup(iprot)),
  79      &        nsup(iprot))) then
  80               do j=nnt+nsup(iprot)-1,nnt,-1
  81                 do k=1,3
  82                   cref(k,j+i,iref,ib,iprot)=cref(k,j,iref,ib,iprot)
  83                   cref(k,nres+j+i,iref,ib,iprot)=
  84      &                             cref(k,nres_pdb+j,iref,ib,iprot)
  85                 enddo
  86                 phi_ref(j+i,iref,ib,iprot)=phi_ref(j,iref,ib,iprot)
  87                 theta_ref(j+i,iref,ib,iprot)=theta_ref(j,iref,ib,iprot)
  88                 xxref(j+i,iref,ib,iprot)=xxref(j,iref,ib,iprot)
  89                 yyref(j+i,iref,ib,iprot)=yyref(j,iref,ib,iprot)
  90                 zzref(j+i,iref,ib,iprot)=zzref(j,iref,ib,iprot)
  91               enddo
  92 #ifdef DEBUG
  93               write (iout,*) "Free energy",efreeref(iref,ib,iprot)
  94               do j=nnt,nct
  95                 write (iout,'(i5,3f10.5,5x,3f10.5)')
  96      &            j,(cref(k,j,iref,ib,iprot),k=1,3),
  97      &              (cref(k,j+nres,iref,ib,iprot),k=1,3)
  98               enddo
  99 #endif
 100               nstart_seq=nnt+i
 101               nstart_sup(iprot)=nnt+i
 102               goto 111
 103             endif
 104           enddo
 105           write (iout,'(a)')
 106      &            'Error - sequences to be superposed do not match.'
 107           return1
 108         else
 109           do i=0,nsup(iprot)-(nct-nnt+1)
 110             if (seq_comp(itype(nnt),itype_pdb(nstart_sup(iprot)+i),
 111      &        nct-nnt+1))
 112      &      then
 113               nstart_sup(iprot)=nstart_sup(iprot)+i
 114               nsup(iprot)=nct-nnt+1
 115               goto 111
 116             endif
 117           enddo
 118           write (iout,'(a)')
 119      &            'Error - sequences to be superposed do not match.'
 120         endif
 121   111   continue
 122         write (iout,'(a,i5)')
 123      &   'Experimental structure begins at residue',nstart_seq
 124         ENDDO ! IREF
 125    44   continue
 126         nref(ib,iprot)=iref-1
 127         write (iout,*) "NREF",ib,iprot,nref(ib,iprot)
 128         call flush(iout)
 129       else
 130         DO IREF=1,MAXREF
 131         call read_angles(inp,*38)
 132         goto 39
 133    38   nref(ib,iprot)=iref-1
 134         exit
 135    39   call chainbuild
 136         nstart_sup(iprot)=nnt
 137         nstart_seq=nnt
 138         nsup(iprot)=nct-nnt+1
 139         do i=1,2*nres
 140           do j=1,3
 141             cref(j,i,iref,ib,iprot)=c(j,i)
 142           enddo
 143         enddo
 144         phi_ref(i,iref,ib,iprot)=phi(i)
 145         theta_ref(i,iref,ib,iprot)=theta(i)
 146         xxref(i,iref,ib,iprot)=xxtab(i)
 147         yyref(i,iref,ib,iprot)=yytab(i)
 148         zzref(i,iref,ib,iprot)=zztab(i)
 149         ENDDO ! IREF
 150         nref(ib,iprot)=iref-1
 151       endif
 152       nend_sup(iprot)=nstart_sup(iprot)+nsup(iprot)-1
 153 c      do i=1,2*nres
 154 c        do j=1,3
 155 c         c(j,i)=cref(j,i,iprot)
 156 c        enddo
 157 c      enddo
 158       do i=1,nres
 159 c        do j=1,3
 160 c          dc(j,nres+i)=cref(j,nres+i,iprot)-cref(j,i,iprot)
 161 c        enddo
 162 c        if (itype(i).ne.10) then
 163 c          ddsc = dist(i,nres+i)
 164 c          do j=1,3
 165 c            dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
 166 c          enddo
 167 c        else
 168 c          do j=1,3
 169 c            dc_norm(j,nres+i)=0.0d0
 170 c          enddo
 171 c        endif
 172 c        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
 173 c         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
 174 c         dc_norm(3,nres+i)**2
 175         do j=1,3
 176           dc(j,i)=c(j,i+1)-c(j,i)
 177         enddo
 178         ddsc = dist(i,i+1)
 179         do j=1,3
 180           dc_norm(j,i)=dc(j,i)/ddsc
 181         enddo
 182       enddo
 183 c      print *,"Calling contact"
 184 c      call contact(print_contact,ncont_ref(iprot),icont_ref(1,1,iprot),
 185 c     &  nstart_sup(iprot),nend_sup(iprot))
 186 c      print *,"Calling elecont"
 187 c      call elecont(print_contact,ncont_pept_ref(iprot),
 188 c     &   icont_pept_ref(1,1,iprot),
 189 c     &   nstart_sup(iprot),nend_sup(iprot))
 190 c       write (iout,'(a,i3,a,i3,a,i3,a)')
 191 c     &    'Number of residues to be superposed:',nsup(iprot),
 192 c     &    ' (from residue',nstart_sup(iprot),' to residue',
 193 c     &    nend_sup(iprot),').'
 194       ENDDO ! IB
 195       close (ipdbin)
 196       return
 197       end