source/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 '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,kkk
  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         do kkk=1,nperm
  57         nres_pdb=nres
  58         nres=nres0
  59         nstart_seq=nnt
  60         if (nsup.le.(nct-nnt+1)) then
  61           do i=0,nct-nnt+1-nsup
  62             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),
  63      &        nsup)) then
  64               do j=nnt+nsup-1,nnt,-1
  65                 do k=1,3
  66                   cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
  67                 enddo
  68               enddo
  69               do j=nnt+nsup-1,nnt,-1
  70                 do k=1,3
  71                   cref(k,j+i,kkk)=cref(k,j,kkk)
  72                 enddo
  73                 phi_ref(j+i)=phi_ref(j)
  74                 theta_ref(j+i)=theta_ref(j)
  75                 alph_ref(j+i)=alph_ref(j)
  76                 omeg_ref(j+i)=omeg_ref(j)
  77               enddo
  78 #ifdef DEBUG
  79               do j=nnt,nct
  80                 write (iout,'(i5,3f10.5,5x,3f10.5)')
  81      &            j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
  82               enddo
  83 #endif
  84               nstart_seq=nnt+i
  85               nstart_sup=nnt+i
  86               goto 111
  87             endif
  88           enddo
  89           write (iout,'(a)')
  90      &            'Error - sequences to be superposed do not match.'
  91           return1
  92         else
  93           do i=0,nsup-(nct-nnt+1)
  94             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),
  95      &        nct-nnt+1))
  96      &      then
  97               nstart_sup=nstart_sup+i
  98               nsup=nct-nnt+1
  99               goto 111
 100             endif
 101           enddo
 102           write (iout,'(a)')
 103      &            'Error - sequences to be superposed do not match.'
 104         endif
 105         enddo
 106   111   continue
 107         write (iout,'(a,i5)')
 108      &   'Experimental structure begins at residue',nstart_seq
 109       else
 110         call read_angles(inp,*38)
 111         goto 39
 112    38   write (iout,'(a)') 'Error reading reference structure.'
 113         return1
 114    39   call chainbuild
 115         kkk=1
 116         nstart_sup=nnt
 117         nstart_seq=nnt
 118         nsup=nct-nnt+1
 119         do i=1,2*nres
 120           do j=1,3
 121             cref(j,i,kkk)=c(j,i)
 122           enddo
 123         enddo
 124       endif
 125       nend_sup=nstart_sup+nsup-1
 126       do i=1,2*nres
 127         do j=1,3
 128           c(j,i)=cref(j,i,kkk)
 129         enddo
 130       enddo
 131       do i=1,nres
 132         do j=1,3
 133           dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
 134         enddo
 135         if (itype(i).ne.10) then
 136           ddsc = dist(i,nres+i)
 137           do j=1,3
 138             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
 139           enddo
 140         else
 141           do j=1,3
 142             dc_norm(j,nres+i)=0.0d0
 143           enddo
 144         endif
 145 c        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
 146 c         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
 147 c         dc_norm(3,nres+i)**2
 148         do j=1,3
 149           dc(j,i)=c(j,i+1)-c(j,i)
 150         enddo
 151         ddsc = dist(i,i+1)
 152         do j=1,3
 153           dc_norm(j,i)=dc(j,i)/ddsc
 154         enddo
 155       enddo
 156 c      print *,"Calling contact"
 157       call contact(.true.,ncont_ref,icont_ref(1,1),
 158      &  nstart_sup,nend_sup)
 159 c      print *,"Calling elecont"
 160       call elecont(.true.,ncont_pept_ref,
 161      &   icont_pept_ref(1,1),
 162      &   nstart_sup,nend_sup)
 163        write (iout,'(a,i3,a,i3,a,i3,a)')
 164      &    'Number of residues to be superposed:',nsup,
 165      &    ' (from residue',nstart_sup,' to residue',
 166      &    nend_sup,').'
 167       return
 168       end