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