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