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