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.FREE'
  10       include 'DIMENSIONS.COMPAR'
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.GEO'
  13       include 'COMMON.VAR'
  14       include 'COMMON.INTERACT'
  15       include 'COMMON.LOCAL'
  16       include 'COMMON.NAMES'
  17       include 'COMMON.CHAIN'
  18       include 'COMMON.FFIELD'
  19       include 'COMMON.SBRIDGE'
  20       include 'COMMON.HEADER'
  21       include 'COMMON.CONTROL'
  22       include 'COMMON.CONTACTS1'
  23       include 'COMMON.PEPTCONT'
  24       include 'COMMON.TIME1'
  25       include 'COMMON.COMPAR'
  26       character*4 sequence(maxres)
  27       integer rescode
  28       double precision x(maxvar)
  29       integer itype_pdb(maxres)
  30       logical seq_comp
  31       integer i,j,k,nres_pdb,iaux
  32       double precision ddsc,dist
  33       integer ilen
  34       external ilen
  35 C
  36       nres0=nres
  37       write (iout,*) "pdbref",pdbref
  38       if (pdbref) then
  39         read(inp,'(a)') pdbfile
  40         write (iout,'(2a,1h.)') 'PDB data will be read from file ',
  41      &    pdbfile(:ilen(pdbfile))
  42         open(ipdbin,file=pdbfile,status='old',err=33)
  43         goto 34
  44   33    write (iout,'(a)') 'Error opening PDB file.'
  45         return1
  46   34    continue
  47         do i=1,nres
  48           itype_pdb(i)=itype(i)
  49         enddo
  50         call readpdb(.true.)
  51         do i=1,nres
  52           iaux=itype_pdb(i)
  53           itype_pdb(i)=itype(i)
  54           itype(i)=iaux
  55         enddo
  56         close (ipdbin)
  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)=cref(k,nres_pdb+j)
  67                 enddo
  68               enddo
  69               do j=nnt+nsup-1,nnt,-1
  70                 do k=1,3
  71                   cref(k,j+i)=cref(k,j)
  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),k=1,3),(cref(k,j+nres),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   111   continue
 106         write (iout,'(a,i5)')
 107      &   'Experimental structure begins at residue',nstart_seq
 108       else
 109         call read_angles(inp,*38)
 110         goto 39
 111    38   write (iout,'(a)') 'Error reading reference structure.'
 112         return1
 113    39   call chainbuild
 114         nstart_sup=nnt
 115         nstart_seq=nnt
 116         nsup=nct-nnt+1
 117         do i=1,2*nres
 118           do j=1,3
 119             cref(j,i)=c(j,i)
 120           enddo
 121         enddo
 122       endif
 123       nend_sup=nstart_sup+nsup-1
 124       do i=1,2*nres
 125         do j=1,3
 126           c(j,i)=cref(j,i)
 127         enddo
 128       enddo
 129       do i=1,nres
 130         do j=1,3
 131           dc(j,nres+i)=cref(j,nres+i)-cref(j,i)
 132         enddo
 133         if (itype(i).ne.10) then
 134           ddsc = dist(i,nres+i)
 135           do j=1,3
 136             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
 137           enddo
 138         else
 139           do j=1,3
 140             dc_norm(j,nres+i)=0.0d0
 141           enddo
 142         endif
 143 c        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
 144 c         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
 145 c         dc_norm(3,nres+i)**2
 146         do j=1,3
 147           dc(j,i)=c(j,i+1)-c(j,i)
 148         enddo
 149         ddsc = dist(i,i+1)
 150         do j=1,3
 151           dc_norm(j,i)=dc(j,i)/ddsc
 152         enddo
 153       enddo
 154 c      print *,"Calling contact"
 155       call contact(.true.,ncont_ref,icont_ref(1,1),
 156      &  nstart_sup,nend_sup)
 157 c      print *,"Calling elecont"
 158       call elecont(.true.,ncont_pept_ref,
 159      &   icont_pept_ref(1,1),
 160      &   nstart_sup,nend_sup)
 161        write (iout,'(a,i3,a,i3,a,i3,a)')
 162      &    'Number of residues to be superposed:',nsup,
 163      &    ' (from residue',nstart_sup,' to residue',
 164      &    nend_sup,').'
 165       return
 166       end