X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC-NEWCORR%2Fread_ref_str.F;fp=source%2Fwham%2Fsrc-NEWSC-NEWCORR%2Fread_ref_str.F;h=4b5618120c3e17a5982ebb082152418628ce8a40;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git

diff --git a/source/wham/src-NEWSC-NEWCORR/read_ref_str.F b/source/wham/src-NEWSC-NEWCORR/read_ref_str.F
new file mode 100644
index 0000000..4b56181
--- /dev/null
+++ b/source/wham/src-NEWSC-NEWCORR/read_ref_str.F
@@ -0,0 +1,165 @@
+      subroutine read_ref_structure(*)
+C
+C Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
+C angles.
+C
+      implicit none
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.LOCAL'
+      include 'COMMON.NAMES'
+      include 'COMMON.CHAIN'
+      include 'COMMON.FFIELD'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.HEADER'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CONTACTS1'
+      include 'COMMON.PEPTCONT'
+      include 'COMMON.TIME1'
+      include 'COMMON.COMPAR'
+      character*4 sequence(maxres)
+      integer rescode
+      double precision x(maxvar)
+      integer itype_pdb(maxres)
+      logical seq_comp
+      integer i,j,k,nres_pdb,iaux
+      double precision ddsc,dist
+      integer ilen
+      external ilen
+C
+      nres0=nres
+      write (iout,*) "pdbref",pdbref
+      if (pdbref) then
+        read(inp,'(a)') pdbfile
+        write (iout,'(2a,1h.)') 'PDB data will be read from file ',
+     &    pdbfile(:ilen(pdbfile))
+        open(ipdbin,file=pdbfile,status='old',err=33)
+        goto 34 
+  33    write (iout,'(a)') 'Error opening PDB file.'
+        return1
+  34    continue
+        do i=1,nres
+          itype_pdb(i)=itype(i)
+        enddo
+        call readpdb(.true.)
+        do i=1,nres
+          iaux=itype_pdb(i)
+          itype_pdb(i)=itype(i)
+          itype(i)=iaux
+        enddo
+        close (ipdbin)
+        nres_pdb=nres
+        nres=nres0
+        nstart_seq=nnt
+        if (nsup.le.(nct-nnt+1)) then
+          do i=0,nct-nnt+1-nsup
+            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),
+     &        nsup)) then
+              do j=nnt+nsup-1,nnt,-1
+                do k=1,3
+                  cref(k,nres+j+i)=cref(k,nres_pdb+j)
+                enddo
+              enddo
+              do j=nnt+nsup-1,nnt,-1
+                do k=1,3
+                  cref(k,j+i)=cref(k,j)
+                enddo
+                phi_ref(j+i)=phi_ref(j)
+                theta_ref(j+i)=theta_ref(j)
+                alph_ref(j+i)=alph_ref(j)
+                omeg_ref(j+i)=omeg_ref(j)
+              enddo
+#ifdef DEBUG
+              do j=nnt,nct
+                write (iout,'(i5,3f10.5,5x,3f10.5)') 
+     &            j,(cref(k,j),k=1,3),(cref(k,j+nres),k=1,3)
+              enddo
+#endif
+              nstart_seq=nnt+i
+              nstart_sup=nnt+i
+              goto 111
+            endif
+          enddo
+          write (iout,'(a)') 
+     &            'Error - sequences to be superposed do not match.'
+          return1
+        else
+          do i=0,nsup-(nct-nnt+1)
+            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),
+     &        nct-nnt+1)) 
+     &      then
+              nstart_sup=nstart_sup+i
+              nsup=nct-nnt+1
+              goto 111
+            endif
+          enddo 
+          write (iout,'(a)') 
+     &            'Error - sequences to be superposed do not match.'
+        endif
+  111   continue
+        write (iout,'(a,i5)') 
+     &   'Experimental structure begins at residue',nstart_seq
+      else
+        call read_angles(inp,*38)
+        goto 39
+   38   write (iout,'(a)') 'Error reading reference structure.'
+        return1
+   39   call chainbuild     
+        nstart_sup=nnt
+        nstart_seq=nnt
+        nsup=nct-nnt+1
+        do i=1,2*nres
+          do j=1,3
+            cref(j,i)=c(j,i)
+          enddo
+        enddo
+      endif
+      nend_sup=nstart_sup+nsup-1
+      do i=1,2*nres
+        do j=1,3
+          c(j,i)=cref(j,i)
+        enddo
+      enddo
+      do i=1,nres
+        do j=1,3
+          dc(j,nres+i)=cref(j,nres+i)-cref(j,i)
+        enddo
+        if (itype(i).ne.10) then
+          ddsc = dist(i,nres+i)
+          do j=1,3
+            dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
+          enddo
+        else
+          do j=1,3
+            dc_norm(j,nres+i)=0.0d0
+          enddo
+        endif
+c        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
+c         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
+c         dc_norm(3,nres+i)**2
+        do j=1,3
+          dc(j,i)=c(j,i+1)-c(j,i)
+        enddo
+        ddsc = dist(i,i+1)
+        do j=1,3
+          dc_norm(j,i)=dc(j,i)/ddsc
+        enddo
+      enddo
+c      print *,"Calling contact"
+      call contact(.true.,ncont_ref,icont_ref(1,1),
+     &  nstart_sup,nend_sup)
+c      print *,"Calling elecont"
+      call elecont(.true.,ncont_pept_ref,
+     &   icont_pept_ref(1,1),
+     &   nstart_sup,nend_sup)
+       write (iout,'(a,i3,a,i3,a,i3,a)')
+     &    'Number of residues to be superposed:',nsup,
+     &    ' (from residue',nstart_sup,' to residue',
+     &    nend_sup,').'
+      return
+      end