source/wham/src/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.21 .or. itype(i+1).eq.21) then
  56 #else
  57         if (itype(i).eq.21) then
  58 #endif
  59           itel(i)=0
  60 #ifdef PROCOR
  61         else if (itype(i+1).ne.20) then
  62 #else
  63         else if (itype(i).ne.20) then
  64 #endif
  65           itel(i)=1
  66         else
  67           itel(i)=2
  68         endif
  69       enddo
  70       call read_bridge
  71
  72       if (with_dihed_constr) then
  73
  74       read (inp,*) ndih_constr
  75       if (ndih_constr.gt.0) then
  76         read (inp,*) ftors
  77         write (iout,*) 'FTORS',ftors
  78         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
  79         write (iout,*)
  80      &   'There are',ndih_constr,' constraints on phi angles.'
  81         do i=1,ndih_constr
  82           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
  83         enddo
  84         do i=1,ndih_constr
  85           phi0(i)=deg2rad*phi0(i)
  86           drange(i)=deg2rad*drange(i)
  87         enddo
  88       endif
  89
  90       endif
  91
  92       nnt=1
  93       nct=nres
  94       if (itype(1).eq.21) nnt=2
  95       if (itype(nres).eq.21) nct=nct-1
  96       write(iout,*) 'NNT=',NNT,' NCT=',NCT
  97 c Read distance restraints
  98       if (constr_dist.gt.0) then
  99         if (refstr) call read_ref_structure(*11)
 100         call read_dist_constr
 101         call hpb_partition
 102       endif
 103
 104       call setup_var
 105       call init_int_table
 106       if (ns.gt.0) then
 107         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 108      &  ' disulfide-bridging cysteines.'
 109         write (iout,'(20i4)') (iss(i),i=1,ns)
 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       write (iout,'(a)')
 122       return
 123    11 stop "Error reading reference structure"
 124       end
 125 c-----------------------------------------------------------------------------
 126       logical function seq_comp(itypea,itypeb,length)
 127       implicit none
 128       integer length,itypea(length),itypeb(length)
 129       integer i
 130       do i=1,length
 131         if (itypea(i).ne.itypeb(i)) then
 132           seq_comp=.false.
 133           return
 134         endif
 135       enddo
 136       seq_comp=.true.
 137       return
 138       end
 139 c-----------------------------------------------------------------------------
 140       subroutine read_bridge
 141 C Read information about disulfide bridges.
 142       implicit real*8 (a-h,o-z)
 143       include 'DIMENSIONS'
 144       include 'DIMENSIONS.ZSCOPT'
 145       include 'COMMON.IOUNITS'
 146       include 'COMMON.GEO'
 147       include 'COMMON.VAR'
 148       include 'COMMON.INTERACT'
 149       include 'COMMON.NAMES'
 150       include 'COMMON.CHAIN'
 151       include 'COMMON.FFIELD'
 152       include 'COMMON.SBRIDGE'
 153 C Read bridging residues.
 154       read (inp,*) ns,(iss(i),i=1,ns)
 155       print *,'ns=',ns
 156       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 157 C Check whether the specified bridging residues are cystines.
 158       do i=1,ns
 159         if (itype(iss(i)).ne.1) then
 160           write (iout,'(2a,i3,a)')
 161      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 162      &   ' can form a disulfide bridge?!!!'
 163           write (*,'(2a,i3,a)')
 164      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 165      &   ' can form a disulfide bridge?!!!'
 166          stop
 167         endif
 168       enddo
 169 C Read preformed bridges.
 170       if (ns.gt.0) then
 171       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 172       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 173       if (nss.gt.0) then
 174         nhpb=nss
 175 C Check if the residues involved in bridges are in the specified list of
 176 C bridging residues.
 177         do i=1,nss
 178           do j=1,i-1
 179             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 180      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 181               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 182      &      ' contains residues present in other pairs.'
 183               write (*,'(a,i3,a)') 'Disulfide pair',i,
 184      &      ' contains residues present in other pairs.'
 185               stop
 186             endif
 187           enddo
 188           do j=1,ns
 189             if (ihpb(i).eq.iss(j)) goto 10
 190           enddo
 191           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 192    10     continue
 193           do j=1,ns
 194             if (jhpb(i).eq.iss(j)) goto 20
 195           enddo
 196           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 197    20     continue
 198           dhpb(i)=dbr
 199           forcon(i)=fbr
 200         enddo
 201         do i=1,nss
 202           ihpb(i)=ihpb(i)+nres
 203           jhpb(i)=jhpb(i)+nres
 204         enddo
 205       endif
 206       endif
 207       return
 208       end
 209 c------------------------------------------------------------------------------
 210       subroutine read_angles(kanal,iscor,energ,iprot,*)
 211       implicit real*8 (a-h,o-z)
 212       include 'DIMENSIONS'
 213       include 'DIMENSIONS.ZSCOPT'
 214       include 'COMMON.INTERACT'
 215       include 'COMMON.SBRIDGE'
 216       include 'COMMON.GEO'
 217       include 'COMMON.VAR'
 218       include 'COMMON.CHAIN'
 219       include 'COMMON.IOUNITS'
 220       character*80 lineh
 221       read(kanal,'(a80)',end=10,err=10) lineh
 222       read(lineh(:5),*,err=8) ic
 223       read(lineh(6:),*,err=8) energ
 224       goto 9
 225     8 ic=1
 226       print *,'error, assuming e=1d10',lineh
 227       energ=1d10
 228       nss=0
 229     9 continue
 230       read(lineh(18:),*,end=10,err=10) nss
 231       IF (NSS.LT.9) THEN
 232         read (lineh(20:),*,end=10,err=10)
 233      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 234       ELSE
 235         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 236         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 237      &    I=9,NSS),iscor
 238       ENDIF
 239 c      print *,"energy",energ," iscor",iscor
 240       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 241       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 242       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 243       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 244       do i=1,nres
 245         theta(i)=deg2rad*theta(i)
 246         phi(i)=deg2rad*phi(i)
 247         alph(i)=deg2rad*alph(i)
 248         omeg(i)=deg2rad*omeg(i)
 249       enddo
 250       return
 251    10 return1
 252       end
 253 c-------------------------------------------------------------------------------
 254       subroutine read_dist_constr
 255       implicit real*8 (a-h,o-z)
 256       include 'DIMENSIONS'
 257       include 'COMMON.CONTROL'
 258       include 'COMMON.CHAIN'
 259       include 'COMMON.IOUNITS'
 260       include 'COMMON.SBRIDGE'
 261       integer ifrag_(2,100),ipair_(2,100)
 262       double precision wfrag_(100),wpair_(100)
 263       character*500 controlcard
 264 c      write (iout,*) "Calling read_dist_constr"
 265 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 266 c      call flush(iout)
 267       call card_concat(controlcard,.true.)
 268       call readi(controlcard,"NFRAG",nfrag_,0)
 269       call readi(controlcard,"NPAIR",npair_,0)
 270       call readi(controlcard,"NDIST",ndist_,0)
 271       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 272       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 273       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 274       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 275       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 276       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 277       write (iout,*) "IFRAG"
 278       do i=1,nfrag_
 279         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 280       enddo
 281       write (iout,*) "IPAIR"
 282       do i=1,npair_
 283         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 284       enddo
 285       call flush(iout)
 286       if (.not.refstr .and. nfrag_.gt.0) then
 287         write (iout,*)
 288      &  "ERROR: no reference structure to compute distance restraints"
 289         write (iout,*)
 290      &  "Restraints must be specified explicitly (NDIST=number)"
 291         stop
 292       endif
 293       if (nfrag_.lt.2 .and. npair_.gt.0) then
 294         write (iout,*) "ERROR: Less than 2 fragments specified",
 295      &   " but distance restraints between pairs requested"
 296         stop
 297       endif
 298       call flush(iout)
 299       do i=1,nfrag_
 300         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 301         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 302      &    ifrag_(2,i)=nstart_sup+nsup-1
 303 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 304         call flush(iout)
 305         if (wfrag_(i).gt.0.0d0) then
 306         do j=ifrag_(1,i),ifrag_(2,i)-1
 307           do k=j+1,ifrag_(2,i)
 308             write (iout,*) "j",j," k",k
 309             ddjk=dist(j,k)
 310             if (constr_dist.eq.1) then
 311             nhpb=nhpb+1
 312             ihpb(nhpb)=j
 313             jhpb(nhpb)=k
 314               dhpb(nhpb)=ddjk
 315             forcon(nhpb)=wfrag_(i)
 316             else if (constr_dist.eq.2) then
 317               if (ddjk.le.dist_cut) then
 318                 nhpb=nhpb+1
 319                 ihpb(nhpb)=j
 320                 jhpb(nhpb)=k
 321                 dhpb(nhpb)=ddjk
 322                 forcon(nhpb)=wfrag_(i)
 323               endif
 324             else
 325               nhpb=nhpb+1
 326               ihpb(nhpb)=j
 327               jhpb(nhpb)=k
 328               dhpb(nhpb)=ddjk
 329               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 330             endif
 331             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 332      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 333           enddo
 334         enddo
 335         endif
 336       enddo
 337       do i=1,npair_
 338         if (wpair_(i).gt.0.0d0) then
 339         ii = ipair_(1,i)
 340         jj = ipair_(2,i)
 341         if (ii.gt.jj) then
 342           itemp=ii
 343           ii=jj
 344           jj=itemp
 345         endif
 346         do j=ifrag_(1,ii),ifrag_(2,ii)
 347           do k=ifrag_(1,jj),ifrag_(2,jj)
 348             nhpb=nhpb+1
 349             ihpb(nhpb)=j
 350             jhpb(nhpb)=k
 351             forcon(nhpb)=wpair_(i)
 352             dhpb(nhpb)=dist(j,k)
 353             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 354      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 355           enddo
 356         enddo
 357         endif
 358       enddo
 359       do i=1,ndist_
 360         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 361      &     ibecarb(i),forcon(nhpb+1)
 362         if (forcon(nhpb+1).gt.0.0d0) then
 363           nhpb=nhpb+1
 364           if (ibecarb(i).gt.0) then
 365             ihpb(i)=ihpb(i)+nres
 366             jhpb(i)=jhpb(i)+nres
 367           endif
 368           if (dhpb(nhpb).eq.0.0d0)
 369      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 370         endif
 371       enddo
 372       do i=1,nhpb
 373           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 374      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 375       enddo
 376       call flush(iout)
 377       return
 378       end