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