source/wham/src-M-homology/read_ref_str.F

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