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         read (inp,*) ftors
  81         write (iout,*) 'FTORS',ftors
  82         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
  83         write (iout,*)
  84      &   'There are',ndih_constr,' constraints on phi angles.'
  85         do i=1,ndih_constr
  86           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
  87         enddo
  88         do i=1,ndih_constr
  89           phi0(i)=deg2rad*phi0(i)
  90           drange(i)=deg2rad*drange(i)
  91         enddo
  92       endif
  93
  94       endif
  95
  96       nnt=1
  97       nct=nres
  98       if (itype(1).eq.ntyp1) nnt=2
  99       if (itype(nres).eq.ntyp1) nct=nct-1
 100       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 101       call setup_var
 102       call init_int_table
 103       if (ns.gt.0) then
 104         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 105      &  ' disulfide-bridging cysteines.'
 106         write (iout,'(20i4)') (iss(i),i=1,ns)
 107        if (dyn_ss) then
 108           write(iout,*)"Running with dynamic disulfide-bond formation"
 109        else
 110         write (iout,'(/a/)') 'Pre-formed links are:'
 111         do i=1,nss
 112           i1=ihpb(i)-nres
 113           i2=jhpb(i)-nres
 114           it1=itype(i1)
 115           it2=itype(i2)
 116          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 117      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 118      &    dhpb(i),ebr,forcon(i)
 119         enddo
 120       endif
 121       endif
 122       write (iout,'(a)')
 123       if (ns.gt.0.and.dyn_ss) then
 124           do i=nss+1,nhpb
 125             ihpb(i-nss)=ihpb(i)
 126             jhpb(i-nss)=jhpb(i)
 127             forcon(i-nss)=forcon(i)
 128             dhpb(i-nss)=dhpb(i)
 129           enddo
 130           nhpb=nhpb-nss
 131           nss=0
 132           call hpb_partition
 133           do i=1,ns
 134             dyn_ss_mask(iss(i))=.true.
 135           enddo
 136       endif
 137       return
 138       end
 139 c-----------------------------------------------------------------------------
 140       logical function seq_comp(itypea,itypeb,length)
 141       implicit none
 142       integer length,itypea(length),itypeb(length)
 143       integer i
 144       do i=1,length
 145         if (itypea(i).ne.itypeb(i)) then
 146           seq_comp=.false.
 147           return
 148         endif
 149       enddo
 150       seq_comp=.true.
 151       return
 152       end
 153 c-----------------------------------------------------------------------------
 154       subroutine read_bridge
 155 C Read information about disulfide bridges.
 156       implicit real*8 (a-h,o-z)
 157       include 'DIMENSIONS'
 158       include 'DIMENSIONS.ZSCOPT'
 159       include 'COMMON.IOUNITS'
 160       include 'COMMON.GEO'
 161       include 'COMMON.VAR'
 162       include 'COMMON.INTERACT'
 163       include 'COMMON.NAMES'
 164       include 'COMMON.CHAIN'
 165       include 'COMMON.FFIELD'
 166       include 'COMMON.SBRIDGE'
 167 C Read bridging residues.
 168       read (inp,*) ns,(iss(i),i=1,ns)
 169       print *,'ns=',ns
 170       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 171 C Check whether the specified bridging residues are cystines.
 172       do i=1,ns
 173         if (itype(iss(i)).ne.1) then
 174           write (iout,'(2a,i3,a)')
 175      &   'Do you REALLY think that the residue ',
 176      &    restyp(itype(iss(i))),i,
 177      &   ' can form a disulfide bridge?!!!'
 178           write (*,'(2a,i3,a)')
 179      &   'Do you REALLY think that the residue ',
 180      &    restyp(itype(iss(i))),i,
 181      &   ' can form a disulfide bridge?!!!'
 182          stop
 183         endif
 184       enddo
 185 C Read preformed bridges.
 186       if (ns.gt.0) then
 187       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 188       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 189       if (nss.gt.0) then
 190         nhpb=nss
 191 C Check if the residues involved in bridges are in the specified list of
 192 C bridging residues.
 193         do i=1,nss
 194           do j=1,i-1
 195             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 196      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 197               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 198      &      ' contains residues present in other pairs.'
 199               write (*,'(a,i3,a)') 'Disulfide pair',i,
 200      &      ' contains residues present in other pairs.'
 201               stop
 202             endif
 203           enddo
 204           do j=1,ns
 205             if (ihpb(i).eq.iss(j)) goto 10
 206           enddo
 207           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 208    10     continue
 209           do j=1,ns
 210             if (jhpb(i).eq.iss(j)) goto 20
 211           enddo
 212           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 213    20     continue
 214           dhpb(i)=dbr
 215           forcon(i)=fbr
 216         enddo
 217         do i=1,nss
 218           ihpb(i)=ihpb(i)+nres
 219           jhpb(i)=jhpb(i)+nres
 220         enddo
 221       endif
 222       endif
 223       return
 224       end
 225 c------------------------------------------------------------------------------
 226       subroutine read_angles(kanal,iscor,energ,iprot,*)
 227       implicit real*8 (a-h,o-z)
 228       include 'DIMENSIONS'
 229       include 'DIMENSIONS.ZSCOPT'
 230       include 'COMMON.INTERACT'
 231       include 'COMMON.SBRIDGE'
 232       include 'COMMON.GEO'
 233       include 'COMMON.VAR'
 234       include 'COMMON.CHAIN'
 235       include 'COMMON.IOUNITS'
 236       character*80 lineh
 237       read(kanal,'(a80)',end=10,err=10) lineh
 238       read(lineh(:5),*,err=8) ic
 239       read(lineh(6:),*,err=8) energ
 240       goto 9
 241     8 ic=1
 242       print *,'error, assuming e=1d10',lineh
 243       energ=1d10
 244       nss=0
 245     9 continue
 246       read(lineh(18:),*,end=10,err=10) nss
 247       IF (NSS.LT.9) THEN
 248         read (lineh(20:),*,end=10,err=10)
 249      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 250       ELSE
 251         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 252         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 253      &    I=9,NSS),iscor
 254       ENDIF
 255 c      print *,"energy",energ," iscor",iscor
 256       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 257       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 258       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 259       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 260       do i=1,nres
 261         theta(i)=deg2rad*theta(i)
 262         phi(i)=deg2rad*phi(i)
 263         alph(i)=deg2rad*alph(i)
 264         omeg(i)=deg2rad*omeg(i)
 265       enddo
 266       return
 267    10 return1
 268       end