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       if (ns.gt.0.and.dyn_ss) then
 208 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
 209 C the bond
 210           do i=nss+1,nhpb
 211 C /06/28/2013 Adasko: nss number of full SS bonds
 212             ihpb(i-nss)=ihpb(i)
 213             jhpb(i-nss)=jhpb(i)
 214             forcon(i-nss)=forcon(i)
 215             dhpb(i-nss)=dhpb(i)
 216           enddo
 217           nhpb=nhpb-nss
 218           nss=0
 219           call hpb_partition
 220           do i=1,ns
 221             dyn_ss_mask(iss(i))=.true.
 222 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
 223 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 224           enddo
 225       endif
 226       return
 227       end
 228 c------------------------------------------------------------------------------
 229       subroutine read_angles(kanal,iscor,energ,iprot,*)
 230       implicit real*8 (a-h,o-z)
 231       include 'DIMENSIONS'
 232       include 'DIMENSIONS.ZSCOPT'
 233       include 'COMMON.INTERACT'
 234       include 'COMMON.SBRIDGE'
 235       include 'COMMON.GEO'
 236       include 'COMMON.VAR'
 237       include 'COMMON.CHAIN'
 238       include 'COMMON.IOUNITS'
 239       character*80 lineh
 240       read(kanal,'(a80)',end=10,err=10) lineh
 241       read(lineh(:5),*,err=8) ic
 242       read(lineh(6:),*,err=8) energ
 243       goto 9
 244     8 ic=1
 245       print *,'error, assuming e=1d10',lineh
 246       energ=1d10
 247       nss=0
 248     9 continue
 249       read(lineh(18:),*,end=10,err=10) nss
 250       IF (NSS.LT.9) THEN
 251         read (lineh(20:),*,end=10,err=10)
 252      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 253       ELSE
 254         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 255         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 256      &    I=9,NSS),iscor
 257       ENDIF
 258 c      print *,"energy",energ," iscor",iscor
 259       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 260       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 261       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 262       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 263       do i=1,nres
 264         theta(i)=deg2rad*theta(i)
 265         phi(i)=deg2rad*phi(i)
 266         alph(i)=deg2rad*alph(i)
 267         omeg(i)=deg2rad*omeg(i)
 268       enddo
 269       return
 270    10 return1
 271       end
 272 c-------------------------------------------------------------------------------
 273       subroutine read_dist_constr
 274       implicit real*8 (a-h,o-z)
 275       include 'DIMENSIONS'
 276       include 'COMMON.CONTROL'
 277       include 'COMMON.CHAIN'
 278       include 'COMMON.IOUNITS'
 279       include 'COMMON.SBRIDGE'
 280       integer ifrag_(2,100),ipair_(2,100)
 281       double precision wfrag_(100),wpair_(100)
 282       character*500 controlcard
 283 c      write (iout,*) "Calling read_dist_constr"
 284 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 285 c      call flush(iout)
 286       call card_concat(controlcard,.true.)
 287       call readi(controlcard,"NFRAG",nfrag_,0)
 288       call readi(controlcard,"NPAIR",npair_,0)
 289       call readi(controlcard,"NDIST",ndist_,0)
 290       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 291       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 292       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 293       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 294       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 295       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 296       write (iout,*) "IFRAG"
 297       do i=1,nfrag_
 298         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 299       enddo
 300       write (iout,*) "IPAIR"
 301       do i=1,npair_
 302         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 303       enddo
 304       call flush(iout)
 305       if (.not.refstr .and. nfrag_.gt.0) then
 306         write (iout,*)
 307      &  "ERROR: no reference structure to compute distance restraints"
 308         write (iout,*)
 309      &  "Restraints must be specified explicitly (NDIST=number)"
 310         stop
 311       endif
 312       if (nfrag_.lt.2 .and. npair_.gt.0) then
 313         write (iout,*) "ERROR: Less than 2 fragments specified",
 314      &   " but distance restraints between pairs requested"
 315         stop
 316       endif
 317       call flush(iout)
 318       do i=1,nfrag_
 319         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 320         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 321      &    ifrag_(2,i)=nstart_sup+nsup-1
 322 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 323         call flush(iout)
 324         if (wfrag_(i).gt.0.0d0) then
 325         do j=ifrag_(1,i),ifrag_(2,i)-1
 326           do k=j+1,ifrag_(2,i)
 327             write (iout,*) "j",j," k",k
 328             ddjk=dist(j,k)
 329             if (constr_dist.eq.1) then
 330             nhpb=nhpb+1
 331             ihpb(nhpb)=j
 332             jhpb(nhpb)=k
 333               dhpb(nhpb)=ddjk
 334             forcon(nhpb)=wfrag_(i)
 335             else if (constr_dist.eq.2) then
 336               if (ddjk.le.dist_cut) then
 337                 nhpb=nhpb+1
 338                 ihpb(nhpb)=j
 339                 jhpb(nhpb)=k
 340                 dhpb(nhpb)=ddjk
 341                 forcon(nhpb)=wfrag_(i)
 342               endif
 343             else
 344               nhpb=nhpb+1
 345               ihpb(nhpb)=j
 346               jhpb(nhpb)=k
 347               dhpb(nhpb)=ddjk
 348               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 349             endif
 350             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 351      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 352           enddo
 353         enddo
 354         endif
 355       enddo
 356       do i=1,npair_
 357         if (wpair_(i).gt.0.0d0) then
 358         ii = ipair_(1,i)
 359         jj = ipair_(2,i)
 360         if (ii.gt.jj) then
 361           itemp=ii
 362           ii=jj
 363           jj=itemp
 364         endif
 365         do j=ifrag_(1,ii),ifrag_(2,ii)
 366           do k=ifrag_(1,jj),ifrag_(2,jj)
 367             nhpb=nhpb+1
 368             ihpb(nhpb)=j
 369             jhpb(nhpb)=k
 370             forcon(nhpb)=wpair_(i)
 371             dhpb(nhpb)=dist(j,k)
 372             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 373      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 374           enddo
 375         enddo
 376         endif
 377       enddo
 378       do i=1,ndist_
 379         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 380      &     ibecarb(i),forcon(nhpb+1)
 381         if (forcon(nhpb+1).gt.0.0d0) then
 382           nhpb=nhpb+1
 383           if (ibecarb(i).gt.0) then
 384             ihpb(i)=ihpb(i)+nres
 385             jhpb(i)=jhpb(i)+nres
 386           endif
 387           if (dhpb(nhpb).eq.0.0d0)
 388      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 389         endif
 390       enddo
 391       do i=1,nhpb
 392           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 393      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 394       enddo
 395       call flush(iout)
 396       return
 397       end