source/maxlik/src-Fmatch/molread_zs.F

   1       subroutine molread_zs(molnum)
   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.ALLPROT'
  18       include 'COMMON.TORCNSTR'
  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",molnum
  35         stop
  36       endif
  37       write(iout,*) 'molnum:',molnum,' nres=',nres
  38 C Read sequence of the protein
  39       if (iscode.gt.0) then
  40         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  41       else
  42         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  43       endif
  44 C Convert sequence to numeric code
  45       do i=1,nres
  46         itype(i)=rescode(i,sequence(i),iscode,.true.)
  47       enddo
  48       write (iout,*) "Numeric code:"
  49       write (iout,'(20i4)') (itype(i),i=1,nres)
  50       do i=1,nres-1
  51 #ifdef PROCOR
  52         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
  53 #else
  54         if (itype(i).eq.ntyp1) then
  55 #endif
  56           itel(i)=0
  57 #ifdef PROCOR
  58         else if (itype(i+1).ne.20) then
  59 #else
  60         else if (itype(i).ne.20) then
  61 #endif
  62           itel(i)=1
  63         else
  64           itel(i)=2
  65         endif
  66       enddo
  67       call read_bridge
  68       nnt=1
  69       nct=nres
  70       if (itype(1).eq.ntyp1) nnt=2
  71       if (itype(nres).eq.ntyp1) nct=nct-1
  72       write(iout,*) 'NNT=',NNT,' NCT=',NCT
  73       call setup_var
  74       call init_int_table
  75       call seq2chains(nres,itype,nchain,chain_length,chain_border)
  76       write(iout,*) "Protein",molnum," nres",nres," nchain",nchain
  77       do i=1,nchain
  78         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
  79      &    chain_border(2,i)
  80       enddo
  81       call chain_symmetry(nchain,nres,itype,chain_border,
  82      &    chain_length,npermchain,tabpermchain)
  83       call store_molinfo(molnum)
  84       if (ns_zs(molnum).gt.0) then
  85         write (iout,'(/a,i3,a)') 'The chain contains',ns_zs(molnum),
  86      &  ' disulfide-bridging cysteines.'
  87         write (iout,'(20i4)') (iss_zs(i,molnum),i=1,ns_zs(molnum))
  88         write (iout,'(/a/)') 'Pre-formed links are:'
  89         do i=1,nss_zs(1,molnum)
  90           i1=ihpb_zs(i,1,molnum)-nres_zs(molnum)
  91           i2=jhpb_zs(i,1,molnum)-nres_zs(molnum)
  92           it1=itype_zs(i1,molnum)
  93           it2=itype_zs(i2,molnum)
  94          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
  95      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
  96      &    dhpb_zs(i,molnum),ebr,forcon_zs(i,molnum)
  97         enddo
  98       endif
  99       write (iout,*) "Protein:",molnum," leaving MOLREAD_ZS"
 100       call flush(iout)
 101       return
 102       end
 103 c-----------------------------------------------------------------------------
 104       logical function seq_comp(itypea,itypeb,length)
 105       implicit none
 106       integer length,itypea(length),itypeb(length)
 107       integer i
 108       do i=1,length
 109         if (itypea(i).ne.itypeb(i)) then
 110           seq_comp=.false.
 111           return
 112         endif
 113       enddo
 114       seq_comp=.true.
 115       return
 116       end
 117 c-----------------------------------------------------------------------------
 118       subroutine read_bridge
 119 C Read information about disulfide bridges.
 120       implicit real*8 (a-h,o-z)
 121       include 'DIMENSIONS'
 122       include 'DIMENSIONS.ZSCOPT'
 123       include 'COMMON.IOUNITS'
 124       include 'COMMON.GEO'
 125       include 'COMMON.VAR'
 126       include 'COMMON.INTERACT'
 127       include 'COMMON.NAMES'
 128       include 'COMMON.CHAIN'
 129       include 'COMMON.FFIELD'
 130       include 'COMMON.SBRIDGE'
 131       include 'COMMON.ALLPROT'
 132 C Read bridging residues.
 133       read (inp,*) ns,(iss(i),i=1,ns)
 134 c      print *,'ns=',ns
 135       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 136 C Check whether the specified bridging residues are cystines.
 137       do i=1,ns
 138         if (itype(iss(i)).ne.1) then
 139           write (iout,'(2a,i3,a)')
 140      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 141      &   ' can form a disulfide bridge?!!!'
 142           write (*,'(2a,i3,a)')
 143      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 144      &   ' can form a disulfide bridge?!!!'
 145          stop
 146         endif
 147       enddo
 148 C Read preformed bridges.
 149       if (ns.gt.0) then
 150       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 151       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 152       if (nss.gt.0) then
 153         nhpb=nss
 154 C Check if the residues involved in bridges are in the specified list of
 155 C bridging residues.
 156         do i=1,nss
 157           do j=1,i-1
 158             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 159      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 160               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 161      &      ' contains residues present in other pairs.'
 162               write (*,'(a,i3,a)') 'Disulfide pair',i,
 163      &      ' contains residues present in other pairs.'
 164               stop
 165             endif
 166           enddo
 167           do j=1,ns
 168             if (ihpb(i).eq.iss(j)) goto 10
 169           enddo
 170           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 171    10     continue
 172           do j=1,ns
 173             if (jhpb(i).eq.iss(j)) goto 20
 174           enddo
 175           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 176    20     continue
 177           dhpb(i)=dbr
 178           forcon(i)=fbr
 179         enddo
 180         do i=1,nss
 181           ihpb(i)=ihpb(i)+nres
 182           jhpb(i)=jhpb(i)+nres
 183         enddo
 184       endif
 185       endif
 186       return
 187       end
 188 c------------------------------------------------------------------------------
 189       subroutine store_molinfo(molnum)
 190       implicit real*8 (a-h,o-z)
 191       include "DIMENSIONS"
 192       include 'DIMENSIONS.ZSCOPT'
 193       include "COMMON.ALLPROT"
 194       include 'COMMON.CHAIN'
 195       include 'COMMON.SBRIDGE'
 196       include "COMMON.VAR"
 197       include "COMMON.INTERACT"
 198       include "COMMON.LOCAL"
 199       include "COMMON.IOUNITS"
 200       nres_zs(molnum)=nres
 201       nnt_zs(molnum)=nnt
 202       nct_zs(molnum)=nct
 203 c      write (iout,*) "store_molinfo nres",nres_zs(molnum),
 204 c     & " nnt",nnt_zs(molnum)," nct",nct_zs(molnum)
 205       iatsc_s_zs(molnum)=iatsc_s
 206       iatsc_e_zs(molnum)=iatsc_e
 207       iatel_s_zs(molnum)=iatel_s
 208       iatel_e_zs(molnum)=iatel_e
 209       iturn3_start_zs(molnum)=iturn3_start
 210       iturn3_end_zs(molnum)=iturn3_end
 211       iturn4_start_zs(molnum)=iturn4_start
 212       iturn4_end_zs(molnum)=iturn4_end
 213       iatscp_s_zs(molnum)=iatscp_s
 214       iatscp_e_zs(molnum)=iatscp_e
 215       loc_start_zs(molnum)=loc_start
 216       loc_end_zs(molnum)=loc_end
 217       ithet_start_zs(molnum)=ithet_start
 218       ithet_end_zs(molnum)=ithet_end
 219       iphi_start_zs(molnum)=iphi_start
 220       iphi_end_zs(molnum)=iphi_end
 221       itau_start_zs(molnum)=itau_start
 222       itau_end_zs(molnum)=itau_end
 223       nvar_zs(molnum)=nvar
 224       nphi_zs(molnum)=nphi
 225       ntheta_zs(molnum)=ntheta
 226       nside_zs(molnum)=nside
 227       do i=1,nres
 228         ialph_zs(i,1,molnum)=ialph(i,1)
 229         ialph_zs(i,2,molnum)=ialph(i,2)
 230       enddo
 231       do i=1,nres
 232         nint_gr_zs(i,molnum)=nint_gr(i)
 233         nscp_gr_zs(i,molnum)=nscp_gr(i)
 234         do j=1,nint_gr(i)
 235           istart_zs(i,j,molnum)=istart(i,j)
 236           iend_zs(i,j,molnum)=iend(i,j)
 237         enddo
 238         itype_zs(i,molnum)=itype(i)
 239         itel_zs(i,molnum)=itel(i)
 240         ielstart_zs(i,molnum)=ielstart(i)
 241         ielend_zs(i,molnum)=ielend(i)
 242 c        write (iout,*) "i",i," nscp_gr",nscp_gr(i),
 243 c     &   " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 244         do j=1,nscp_gr(i)
 245           iscpstart_zs(i,j,molnum)=iscpstart(i,j)
 246           iscpend_zs(i,j,molnum)=iscpend(i,j)
 247         enddo
 248       enddo
 249       ns_zs(molnum)=ns
 250 c      nss_zs(molnum)=nss
 251 c      nfree_zs(molnum)=nfree
 252 c      do i=1,ns
 253 c        iss_zs(i,molnum)=iss(i)
 254 c      enddo
 255 c      do i=1,nss
 256 c        ihpb_zs(i,molnum)=ihpb(i)
 257 c        jhpb_zs(i,molnum)=jhpb(i)
 258 c        dhpb_zs(i,molnum)=dhpb(i)
 259 c        forcon_zs(i,molnum)=forcon(i)
 260 c      enddo
 261       link_start_zs(molnum)=link_start
 262       link_end_zs(molnum)=link_end
 263       nchain_zs(molnum)=nchain
 264       npermchain_zs(molnum)=npermchain
 265       do i=1,nchain
 266         chain_border_zs(1,i,molnum)= chain_border(1,i)
 267         chain_border_zs(2,i,molnum)= chain_border(2,i)
 268         chain_length_zs(i,molnum)=chain_length(i)
 269       enddo
 270       do i=1,npermchain
 271         do j=1,nchain
 272           tabpermchain_zs(j,i,molnum)=tabpermchain(j,i)
 273         enddo
 274       enddo
 275       return
 276       end
 277 c------------------------------------------------------------------------------
 278       subroutine restore_molinfo(molnum)
 279       implicit real*8 (a-h,o-z)
 280       include "DIMENSIONS"
 281       include 'DIMENSIONS.ZSCOPT'
 282       include "COMMON.ALLPROT"
 283       include "COMMON.SBRIDGE"
 284       include "COMMON.CHAIN"
 285       include "COMMON.VAR"
 286       include "COMMON.INTERACT"
 287       include "COMMON.TORCNSTR"
 288       include "COMMON.LOCAL"
 289       ndih_constr=0
 290       nres=nres_zs(molnum)
 291       nnt=nnt_zs(molnum)
 292       nct=nct_zs(molnum)
 293       iatsc_s=iatsc_s_zs(molnum)
 294       iatsc_e=iatsc_e_zs(molnum)
 295       iatel_s=iatel_s_zs(molnum)
 296       iatel_e=iatel_e_zs(molnum)
 297       iturn3_start=iturn3_start_zs(molnum)
 298       iturn3_end=iturn3_end_zs(molnum)
 299       iturn4_start=iturn4_start_zs(molnum)
 300       iturn4_end=iturn4_end_zs(molnum)
 301       iatscp_s=iatscp_s_zs(molnum)
 302       iatscp_e=iatscp_e_zs(molnum)
 303       ithet_start=ithet_start_zs(molnum)
 304       ithet_end=ithet_end_zs(molnum)
 305       iphi_start=iphi_start_zs(molnum)
 306       iphi_end=iphi_end_zs(molnum)
 307       itau_start=itau_start_zs(molnum)
 308       itau_end=itau_end_zs(molnum)
 309       loc_start=loc_start_zs(molnum)
 310       loc_end=loc_end_zs(molnum)
 311       nvar=nvar_zs(molnum)
 312       nphi=nphi_zs(molnum)
 313       ntheta=ntheta_zs(molnum)
 314       nside=nside_zs(molnum)
 315       do i=1,nres
 316         ialph(i,1)=ialph_zs(i,1,molnum)
 317         ialph(i,2)=ialph_zs(i,2,molnum)
 318       enddo
 319       do i=1,nres
 320         nint_gr(i)=nint_gr_zs(i,molnum)
 321         nscp_gr(i)=nscp_gr_zs(i,molnum)
 322         do j=1,nint_gr(i)
 323           istart(i,j)=istart_zs(i,j,molnum)
 324           iend(i,j)=iend_zs(i,j,molnum)
 325         enddo
 326         itype(i)=itype_zs(i,molnum)
 327         itel(i)=itel_zs(i,molnum)
 328         ielstart(i)=ielstart_zs(i,molnum)
 329         ielend(i)=ielend_zs(i,molnum)
 330         do j=1,nscp_gr(i)
 331           iscpstart(i,j)=iscpstart_zs(i,j,molnum)
 332           iscpend(i,j)=iscpend_zs(i,j,molnum)
 333         enddo
 334       enddo
 335       ns=ns_zs(molnum)
 336 c      nss=nss_zs(molnum)
 337 c      nfree=nfree_zs(molnum)
 338 c      do i=1,ns
 339 c        iss(i)=iss_zs(i,molnum)
 340 c      enddo
 341 c      do i=1,nss
 342 c        ihpb(i)=ihpb_zs(i,molnum)
 343 c        jhpb(i)=jhpb_zs(i,molnum)
 344 c        dhpb(i)=dhpb_zs(i,molnum)
 345 c        forcon(i)=forcon_zs(i,molnum)
 346 c      enddo
 347       link_start=link_start_zs(molnum)
 348       link_end=link_end_zs(molnum)
 349       nchain=nchain_zs(molnum)
 350       npermchain=npermchain_zs(molnum)
 351       do i=1,nchain
 352         chain_border(1,i)=chain_border_zs(1,i,molnum)
 353         chain_border(2,i)=chain_border_zs(2,i,molnum)
 354         chain_length(i)=chain_length_zs(i,molnum)
 355       enddo
 356       do i=1,npermchain
 357         do j=1,nchain
 358           tabpermchain(j,i)=tabpermchain_zs(j,i,molnum)
 359         enddo
 360       enddo
 361       return
 362       end
 363 c------------------------------------------------------------------------------
 364       subroutine read_angles(kanal,iscor,energ,iprot,*)
 365       implicit real*8 (a-h,o-z)
 366       include 'DIMENSIONS'
 367       include 'DIMENSIONS.ZSCOPT'
 368       include 'COMMON.INTERACT'
 369       include 'COMMON.SBRIDGE'
 370       include 'COMMON.GEO'
 371       include 'COMMON.VAR'
 372       include 'COMMON.CHAIN'
 373       include 'COMMON.IOUNITS'
 374       character*80 lineh
 375       read(kanal,'(a80)',end=10,err=10) lineh
 376       read(lineh(:5),*,err=8) ic
 377       read(lineh(6:),*,err=8) energ
 378       goto 9
 379     8 ic=1
 380       print *,'error, assuming e=1d10',lineh
 381       energ=1d10
 382       nss=0
 383     9 continue
 384       read(lineh(18:),*,end=10,err=10) nss
 385       IF (NSS.LT.9) THEN
 386         read (lineh(20:),*,end=10,err=10)
 387      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 388       ELSE
 389         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 390         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 391      &    I=9,NSS),iscor
 392       ENDIF
 393 c      print *,"energy",energ," iscor",iscor
 394       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 395       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 396       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 397       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 398       do i=1,nres
 399         theta(i)=deg2rad*theta(i)
 400         phi(i)=deg2rad*phi(i)
 401         alph(i)=deg2rad*alph(i)
 402         omeg(i)=deg2rad*omeg(i)
 403       enddo
 404       return
 405    10 return1
 406       end