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,nnt_old,nct_old
  31       double precision ddsc,dist
  32       integer ilen,kkk
  33       external ilen
  34 C
  35       nres0=nres
  36       nnt_old=nnt
  37       nct_old=nct
  38       write (iout,*) "pdbref",pdbref
  39       if (pdbref) then
  40         read(inp,'(a)') pdbfile
  41         write (iout,'(2a,1h.)') 'PDB data will be read from file ',
  42      &    pdbfile(:ilen(pdbfile))
  43         open(ipdbin,file=pdbfile,status='old',err=33)
  44         goto 34
  45   33    write (iout,'(a)') 'Error opening PDB file.'
  46         return1
  47   34    continue
  48         do i=1,nres
  49           itype_pdb(i)=itype(i)
  50         enddo
  51         call readpdb(.true.)
  52         do i=1,nres
  53           iaux=itype_pdb(i)
  54           itype_pdb(i)=itype(i)
  55           itype(i)=iaux
  56         enddo
  57         close (ipdbin)
  58         do kkk=1,nperm
  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,kkk)=cref(k,nres_pdb+j,kkk)
  69                 enddo
  70               enddo
  71               do j=nnt+nsup-1,nnt,-1
  72                 do k=1,3
  73                   cref(k,j+i,kkk)=cref(k,j,kkk)
  74                 enddo
  75                 phi_ref(j+i)=phi_ref(j)
  76                 theta_ref(j+i)=theta_ref(j)
  77                 alph_ref(j+i)=alph_ref(j)
  78                 omeg_ref(j+i)=omeg_ref(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,kkk),k=1,3),(cref(k,j+nres,kkk),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         enddo
 108   111   continue
 109         write (iout,'(a,i5)')
 110      &   'Experimental structure begins at residue',nstart_seq
 111       else
 112         call read_angles(inp,*38)
 113         goto 39
 114    38   write (iout,'(a)') 'Error reading reference structure.'
 115         return1
 116    39   call chainbuild
 117         kkk=1
 118         nstart_sup=nnt
 119         nstart_seq=nnt
 120         nsup=nct-nnt+1
 121         do i=1,2*nres
 122           do j=1,3
 123             cref(j,i,kkk)=c(j,i)
 124           enddo
 125         enddo
 126       endif
 127       nend_sup=nstart_sup+nsup-1
 128       do i=1,2*nres
 129         do j=1,3
 130           c(j,i)=cref(j,i,kkk)
 131         enddo
 132       enddo
 133       do i=1,nres
 134         do j=1,3
 135           dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
 136         enddo
 137         if (itype(i).ne.10) then
 138           ddsc = dist(i,nres+i)
 139           do j=1,3
 140             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
 141           enddo
 142         else
 143           do j=1,3
 144             dc_norm(j,nres+i)=0.0d0
 145           enddo
 146         endif
 147 c        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
 148 c         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
 149 c         dc_norm(3,nres+i)**2
 150         do j=1,3
 151           dc(j,i)=c(j,i+1)-c(j,i)
 152         enddo
 153         ddsc = dist(i,i+1)
 154         do j=1,3
 155           dc_norm(j,i)=dc(j,i)/ddsc
 156         enddo
 157       enddo
 158 c      print *,"Calling contact"
 159       call contact(.true.,ncont_ref,icont_ref(1,1),
 160      &  nstart_sup,nend_sup)
 161 c      print *,"Calling elecont"
 162       call elecont(.true.,ncont_pept_ref,
 163      &   icont_pept_ref(1,1),
 164      &   nstart_sup,nend_sup)
 165        write (iout,'(a,i3,a,i3,a,i3,a)')
 166      &    'Number of residues to be superposed:',nsup,
 167      &    ' (from residue',nstart_sup,' to residue',
 168      &    nend_sup,').'
 169       nnt_old=nnt
 170       nct_old=nct
 171       return
 172       end