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