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