source/wham/src-M/molread_zs.F

   1       subroutine molread(*)
   2 C
   3 C Read molecular data.
   4 C
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7       include 'DIMENSIONS.ZSCOPT'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.GEO'
  10       include 'COMMON.VAR'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.LOCAL'
  13       include 'COMMON.NAMES'
  14       include 'COMMON.CHAIN'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.SBRIDGE'
  17       include 'COMMON.TORCNSTR'
  18       include 'COMMON.CONTROL'
  19       character*4 sequence(maxres)
  20       integer rescode
  21       double precision x(maxvar)
  22       character*320 controlcard,ucase
  23       dimension itype_pdb(maxres)
  24       logical seq_comp
  25       call card_concat(controlcard,.true.)
  26       call reada(controlcard,'SCAL14',scal14,0.4d0)
  27       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
  28       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
  29       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
  30       r0_corr=cutoff_corr-delt_corr
  31       call readi(controlcard,"NRES",nres,0)
  32       iscode=index(controlcard,"ONE_LETTER")
  33       if (nres.le.0) then
  34         write (iout,*) "Error: no residues in molecule"
  35         return1
  36       endif
  37       if (nres.gt.maxres) then
  38         write (iout,*) "Error: too many residues",nres,maxres
  39       endif
  40       write(iout,*) 'nres=',nres
  41 C Read sequence of the protein
  42       if (iscode.gt.0) then
  43         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  44       else
  45         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  46       endif
  47 C Convert sequence to numeric code
  48       do i=1,nres
  49         itype(i)=rescode(i,sequence(i),iscode)
  50       enddo
  51       write (iout,*) "Numeric code:"
  52       write (iout,'(20i4)') (itype(i),i=1,nres)
  53       do i=1,nres-1
  54 #ifdef PROCOR
  55         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
  56 #else
  57         if (itype(i).eq.ntyp1) then
  58 #endif
  59           itel(i)=0
  60 #ifdef PROCOR
  61         else if (iabs(itype(i+1)).ne.20) then
  62 #else
  63         else if (iabs(itype(i)).ne.20) then
  64 #endif
  65           itel(i)=1
  66         else
  67           itel(i)=2
  68         endif
  69       enddo
  70        write (iout,*) "ITEL"
  71        do i=1,nres-1
  72          write (iout,*) i,itype(i),itel(i)
  73        enddo
  74       call read_bridge
  75
  76       if (with_dihed_constr) then
  77
  78       read (inp,*) ndih_constr
  79       if (ndih_constr.gt.0) then
  80 C        read (inp,*) ftors
  81 C        write (iout,*) 'FTORS',ftors
  82         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
  83      &   i=1,ndih_constr)
  84         write (iout,*)
  85      &   'There are',ndih_constr,' constraints on phi angles.'
  86         do i=1,ndih_constr
  87           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
  88      &    ftors(i)
  89         enddo
  90         do i=1,ndih_constr
  91           phi0(i)=deg2rad*phi0(i)
  92           drange(i)=deg2rad*drange(i)
  93         enddo
  94       endif
  95
  96       endif
  97
  98       nnt=1
  99       nct=nres
 100       if (itype(1).eq.ntyp1) nnt=2
 101       if (itype(nres).eq.ntyp1) nct=nct-1
 102       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 103       call setup_var
 104       call init_int_table
 105       if (ns.gt.0) then
 106         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 107      &  ' disulfide-bridging cysteines.'
 108         write (iout,'(20i4)') (iss(i),i=1,ns)
 109        if (dyn_ss) then
 110           write(iout,*)"Running with dynamic disulfide-bond formation"
 111        else
 112         write (iout,'(/a/)') 'Pre-formed links are:'
 113         do i=1,nss
 114           i1=ihpb(i)-nres
 115           i2=jhpb(i)-nres
 116           it1=itype(i1)
 117           it2=itype(i2)
 118          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 119      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 120      &    dhpb(i),ebr,forcon(i)
 121         enddo
 122       endif
 123       endif
 124       write (iout,'(a)')
 125       if (ns.gt.0.and.dyn_ss) then
 126           do i=nss+1,nhpb
 127             ihpb(i-nss)=ihpb(i)
 128             jhpb(i-nss)=jhpb(i)
 129             forcon(i-nss)=forcon(i)
 130             dhpb(i-nss)=dhpb(i)
 131           enddo
 132           nhpb=nhpb-nss
 133           nss=0
 134           call hpb_partition
 135           do i=1,ns
 136             dyn_ss_mask(iss(i))=.true.
 137           enddo
 138       endif
 139       return
 140       end
 141 c-----------------------------------------------------------------------------
 142       logical function seq_comp(itypea,itypeb,length)
 143       implicit none
 144       integer length,itypea(length),itypeb(length)
 145       integer i
 146       do i=1,length
 147         if (itypea(i).ne.itypeb(i)) then
 148           seq_comp=.false.
 149           return
 150         endif
 151       enddo
 152       seq_comp=.true.
 153       return
 154       end
 155 c-----------------------------------------------------------------------------
 156       subroutine read_bridge
 157 C Read information about disulfide bridges.
 158       implicit real*8 (a-h,o-z)
 159       include 'DIMENSIONS'
 160       include 'DIMENSIONS.ZSCOPT'
 161       include 'COMMON.IOUNITS'
 162       include 'COMMON.GEO'
 163       include 'COMMON.VAR'
 164       include 'COMMON.INTERACT'
 165       include 'COMMON.NAMES'
 166       include 'COMMON.CHAIN'
 167       include 'COMMON.FFIELD'
 168       include 'COMMON.SBRIDGE'
 169 C Read bridging residues.
 170       read (inp,*) ns,(iss(i),i=1,ns)
 171       print *,'ns=',ns
 172       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 173 C Check whether the specified bridging residues are cystines.
 174       do i=1,ns
 175         if (itype(iss(i)).ne.1) then
 176           write (iout,'(2a,i3,a)')
 177      &   'Do you REALLY think that the residue ',
 178      &    restyp(itype(iss(i))),i,
 179      &   ' can form a disulfide bridge?!!!'
 180           write (*,'(2a,i3,a)')
 181      &   'Do you REALLY think that the residue ',
 182      &    restyp(itype(iss(i))),i,
 183      &   ' can form a disulfide bridge?!!!'
 184          stop
 185         endif
 186       enddo
 187 C Read preformed bridges.
 188       if (ns.gt.0) then
 189       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 190       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 191       if (nss.gt.0) then
 192         nhpb=nss
 193 C Check if the residues involved in bridges are in the specified list of
 194 C bridging residues.
 195         do i=1,nss
 196           do j=1,i-1
 197             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 198      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 199               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 200      &      ' contains residues present in other pairs.'
 201               write (*,'(a,i3,a)') 'Disulfide pair',i,
 202      &      ' contains residues present in other pairs.'
 203               stop
 204             endif
 205           enddo
 206           do j=1,ns
 207             if (ihpb(i).eq.iss(j)) goto 10
 208           enddo
 209           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 210    10     continue
 211           do j=1,ns
 212             if (jhpb(i).eq.iss(j)) goto 20
 213           enddo
 214           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 215    20     continue
 216           dhpb(i)=dbr
 217           forcon(i)=fbr
 218         enddo
 219         do i=1,nss
 220           ihpb(i)=ihpb(i)+nres
 221           jhpb(i)=jhpb(i)+nres
 222         enddo
 223       endif
 224       endif
 225       return
 226       end
 227 c------------------------------------------------------------------------------
 228       subroutine read_angles(kanal,iscor,energ,iprot,*)
 229       implicit real*8 (a-h,o-z)
 230       include 'DIMENSIONS'
 231       include 'DIMENSIONS.ZSCOPT'
 232       include 'COMMON.INTERACT'
 233       include 'COMMON.SBRIDGE'
 234       include 'COMMON.GEO'
 235       include 'COMMON.VAR'
 236       include 'COMMON.CHAIN'
 237       include 'COMMON.IOUNITS'
 238       character*80 lineh
 239       read(kanal,'(a80)',end=10,err=10) lineh
 240       read(lineh(:5),*,err=8) ic
 241       read(lineh(6:),*,err=8) energ
 242       goto 9
 243     8 ic=1
 244       print *,'error, assuming e=1d10',lineh
 245       energ=1d10
 246       nss=0
 247     9 continue
 248       read(lineh(18:),*,end=10,err=10) nss
 249       IF (NSS.LT.9) THEN
 250         read (lineh(20:),*,end=10,err=10)
 251      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 252       ELSE
 253         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 254         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 255      &    I=9,NSS),iscor
 256       ENDIF
 257 c      print *,"energy",energ," iscor",iscor
 258       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 259       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 260       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 261       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 262       do i=1,nres
 263         theta(i)=deg2rad*theta(i)
 264         phi(i)=deg2rad*phi(i)
 265         alph(i)=deg2rad*alph(i)
 266         omeg(i)=deg2rad*omeg(i)
 267       enddo
 268       return
 269    10 return1
 270       end