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