source/wham/src-restraints/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       if (constr_homology.gt.0) then
 105         call read_constr_homology
 106       endif
 107
 108
 109       call setup_var
 110       call init_int_table
 111       if (ns.gt.0) then
 112         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 113      &  ' disulfide-bridging cysteines.'
 114         write (iout,'(20i4)') (iss(i),i=1,ns)
 115         write (iout,'(/a/)') 'Pre-formed links are:'
 116         do i=1,nss
 117           i1=ihpb(i)-nres
 118           i2=jhpb(i)-nres
 119           it1=itype(i1)
 120           it2=itype(i2)
 121          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 122      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 123      &    dhpb(i),ebr,forcon(i)
 124         enddo
 125       endif
 126       write (iout,'(a)')
 127       return
 128    11 stop "Error reading reference structure"
 129       end
 130 c-----------------------------------------------------------------------------
 131       logical function seq_comp(itypea,itypeb,length)
 132       implicit none
 133       integer length,itypea(length),itypeb(length)
 134       integer i
 135       do i=1,length
 136         if (itypea(i).ne.itypeb(i)) then
 137           seq_comp=.false.
 138           return
 139         endif
 140       enddo
 141       seq_comp=.true.
 142       return
 143       end
 144 c-----------------------------------------------------------------------------
 145       subroutine read_bridge
 146 C Read information about disulfide bridges.
 147       implicit real*8 (a-h,o-z)
 148       include 'DIMENSIONS'
 149       include 'DIMENSIONS.ZSCOPT'
 150       include 'COMMON.IOUNITS'
 151       include 'COMMON.GEO'
 152       include 'COMMON.VAR'
 153       include 'COMMON.INTERACT'
 154       include 'COMMON.NAMES'
 155       include 'COMMON.CHAIN'
 156       include 'COMMON.FFIELD'
 157       include 'COMMON.SBRIDGE'
 158 C Read bridging residues.
 159       read (inp,*) ns,(iss(i),i=1,ns)
 160       print *,'ns=',ns
 161       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 162 C Check whether the specified bridging residues are cystines.
 163       do i=1,ns
 164         if (itype(iss(i)).ne.1) then
 165           write (iout,'(2a,i3,a)')
 166      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 167      &   ' can form a disulfide bridge?!!!'
 168           write (*,'(2a,i3,a)')
 169      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 170      &   ' can form a disulfide bridge?!!!'
 171          stop
 172         endif
 173       enddo
 174 C Read preformed bridges.
 175       if (ns.gt.0) then
 176       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 177       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 178       if (nss.gt.0) then
 179         nhpb=nss
 180 C Check if the residues involved in bridges are in the specified list of
 181 C bridging residues.
 182         do i=1,nss
 183           do j=1,i-1
 184             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 185      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 186               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 187      &      ' contains residues present in other pairs.'
 188               write (*,'(a,i3,a)') 'Disulfide pair',i,
 189      &      ' contains residues present in other pairs.'
 190               stop
 191             endif
 192           enddo
 193           do j=1,ns
 194             if (ihpb(i).eq.iss(j)) goto 10
 195           enddo
 196           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 197    10     continue
 198           do j=1,ns
 199             if (jhpb(i).eq.iss(j)) goto 20
 200           enddo
 201           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 202    20     continue
 203           dhpb(i)=dbr
 204           forcon(i)=fbr
 205         enddo
 206         do i=1,nss
 207           ihpb(i)=ihpb(i)+nres
 208           jhpb(i)=jhpb(i)+nres
 209         enddo
 210       endif
 211       endif
 212       if (ns.gt.0.and.dyn_ss) then
 213 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
 214 C the bond
 215           do i=nss+1,nhpb
 216 C /06/28/2013 Adasko: nss number of full SS bonds
 217             ihpb(i-nss)=ihpb(i)
 218             jhpb(i-nss)=jhpb(i)
 219             forcon(i-nss)=forcon(i)
 220             dhpb(i-nss)=dhpb(i)
 221           enddo
 222           nhpb=nhpb-nss
 223           nss=0
 224           call hpb_partition
 225           do i=1,ns
 226             dyn_ss_mask(iss(i))=.true.
 227 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
 228 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 229           enddo
 230       endif
 231       return
 232       end
 233 c------------------------------------------------------------------------------
 234       subroutine read_angles(kanal,iscor,energ,iprot,*)
 235       implicit real*8 (a-h,o-z)
 236       include 'DIMENSIONS'
 237       include 'DIMENSIONS.ZSCOPT'
 238       include 'COMMON.INTERACT'
 239       include 'COMMON.SBRIDGE'
 240       include 'COMMON.GEO'
 241       include 'COMMON.VAR'
 242       include 'COMMON.CHAIN'
 243       include 'COMMON.IOUNITS'
 244       character*80 lineh
 245       read(kanal,'(a80)',end=10,err=10) lineh
 246       read(lineh(:5),*,err=8) ic
 247       read(lineh(6:),*,err=8) energ
 248       goto 9
 249     8 ic=1
 250       print *,'error, assuming e=1d10',lineh
 251       energ=1d10
 252       nss=0
 253     9 continue
 254       read(lineh(18:),*,end=10,err=10) nss
 255       IF (NSS.LT.9) THEN
 256         read (lineh(20:),*,end=10,err=10)
 257      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 258       ELSE
 259         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 260         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 261      &    I=9,NSS),iscor
 262       ENDIF
 263 c      print *,"energy",energ," iscor",iscor
 264       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 265       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 266       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 267       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 268       do i=1,nres
 269         theta(i)=deg2rad*theta(i)
 270         phi(i)=deg2rad*phi(i)
 271         alph(i)=deg2rad*alph(i)
 272         omeg(i)=deg2rad*omeg(i)
 273       enddo
 274       return
 275    10 return1
 276       end
 277 c-------------------------------------------------------------------------------
 278       subroutine read_dist_constr
 279       implicit real*8 (a-h,o-z)
 280       include 'DIMENSIONS'
 281       include 'COMMON.CONTROL'
 282       include 'COMMON.CHAIN'
 283       include 'COMMON.IOUNITS'
 284       include 'COMMON.SBRIDGE'
 285       integer ifrag_(2,100),ipair_(2,100)
 286       double precision wfrag_(100),wpair_(100)
 287       character*500 controlcard
 288 c      write (iout,*) "Calling read_dist_constr"
 289 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 290 c      call flush(iout)
 291       call card_concat(controlcard,.true.)
 292       call readi(controlcard,"NFRAG",nfrag_,0)
 293       call readi(controlcard,"NPAIR",npair_,0)
 294       call readi(controlcard,"NDIST",ndist_,0)
 295       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 296       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 297       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 298       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 299       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 300       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 301       write (iout,*) "IFRAG"
 302       do i=1,nfrag_
 303         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 304       enddo
 305       write (iout,*) "IPAIR"
 306       do i=1,npair_
 307         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 308       enddo
 309       call flush(iout)
 310       if (.not.refstr .and. nfrag_.gt.0) then
 311         write (iout,*)
 312      &  "ERROR: no reference structure to compute distance restraints"
 313         write (iout,*)
 314      &  "Restraints must be specified explicitly (NDIST=number)"
 315         stop
 316       endif
 317       if (nfrag_.lt.2 .and. npair_.gt.0) then
 318         write (iout,*) "ERROR: Less than 2 fragments specified",
 319      &   " but distance restraints between pairs requested"
 320         stop
 321       endif
 322       call flush(iout)
 323       do i=1,nfrag_
 324         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 325         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 326      &    ifrag_(2,i)=nstart_sup+nsup-1
 327 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 328         call flush(iout)
 329         if (wfrag_(i).gt.0.0d0) then
 330         do j=ifrag_(1,i),ifrag_(2,i)-1
 331           do k=j+1,ifrag_(2,i)
 332             write (iout,*) "j",j," k",k
 333             ddjk=dist(j,k)
 334             if (constr_dist.eq.1) then
 335             nhpb=nhpb+1
 336             ihpb(nhpb)=j
 337             jhpb(nhpb)=k
 338               dhpb(nhpb)=ddjk
 339             forcon(nhpb)=wfrag_(i)
 340             else if (constr_dist.eq.2) then
 341               if (ddjk.le.dist_cut) then
 342                 nhpb=nhpb+1
 343                 ihpb(nhpb)=j
 344                 jhpb(nhpb)=k
 345                 dhpb(nhpb)=ddjk
 346                 forcon(nhpb)=wfrag_(i)
 347               endif
 348             else
 349               nhpb=nhpb+1
 350               ihpb(nhpb)=j
 351               jhpb(nhpb)=k
 352               dhpb(nhpb)=ddjk
 353               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 354             endif
 355             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 356      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 357           enddo
 358         enddo
 359         endif
 360       enddo
 361       do i=1,npair_
 362         if (wpair_(i).gt.0.0d0) then
 363         ii = ipair_(1,i)
 364         jj = ipair_(2,i)
 365         if (ii.gt.jj) then
 366           itemp=ii
 367           ii=jj
 368           jj=itemp
 369         endif
 370         do j=ifrag_(1,ii),ifrag_(2,ii)
 371           do k=ifrag_(1,jj),ifrag_(2,jj)
 372             nhpb=nhpb+1
 373             ihpb(nhpb)=j
 374             jhpb(nhpb)=k
 375             forcon(nhpb)=wpair_(i)
 376             dhpb(nhpb)=dist(j,k)
 377             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 378      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 379           enddo
 380         enddo
 381         endif
 382       enddo
 383       do i=1,ndist_
 384         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 385      &     ibecarb(i),forcon(nhpb+1)
 386         if (forcon(nhpb+1).gt.0.0d0) then
 387           nhpb=nhpb+1
 388           if (ibecarb(i).gt.0) then
 389             ihpb(i)=ihpb(i)+nres
 390             jhpb(i)=jhpb(i)+nres
 391           endif
 392           if (dhpb(nhpb).eq.0.0d0)
 393      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 394         endif
 395       enddo
 396       do i=1,nhpb
 397           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 398      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 399       enddo
 400       call flush(iout)
 401       return
 402       end
 403
 404
 405
 406 c====-------------------------------------------------------------------
 407
 408       subroutine read_constr_homology
 409
 410       include 'DIMENSIONS'
 411 #ifdef MPI
 412       include 'mpif.h'
 413 #endif
 414       include 'COMMON.CONTROL'
 415       include 'COMMON.CHAIN'
 416       include 'COMMON.IOUNITS'
 417       include 'COMMON.INTERACT'
 418       include 'COMMON.GEO'
 419       double precision odl_temp,sigma_odl_temp
 420       common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 421      &    sigma_odl_temp(maxres,maxres,max_template)
 422       character*2 kic2
 423       character*24 model_ki_dist, model_ki_angle
 424       character*500 controlcard
 425       integer ki, i, j, k, l
 426       logical lprn /.true./
 427
 428       call card_concat(controlcard,.true.)
 429       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 430       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 431       write (iout,*) "wga_dist",waga_dist," waga_angle",waga_angle
 432
 433       lim_odl=0
 434       lim_dih=0
 435       do i=1,nres
 436         do j=i+2,nres
 437           do ki=1,constr_homology
 438             sigma_odl_temp(i,j,ki)=0.0d0
 439             odl_temp(i,j,ki)=0.0d0
 440           enddo
 441         enddo
 442       enddo
 443       do i=1,nres-3
 444         do ki=1,constr_homology
 445           dih(ki,i)=0.0d0
 446           sigma_dih(ki,i)=0.0d0
 447         enddo
 448       enddo
 449       do ki=1,constr_homology
 450           write(kic2,'(i2)') ki
 451           if (ki.le.9) kic2="0"//kic2(2:2)
 452
 453           model_ki_dist="model"//kic2//".dist"
 454           model_ki_angle="model"//kic2//".angle"
 455         open (ientin,file=model_ki_dist,status='old')
 456         do irec=1,maxdim !petla do czytania wiezow na odleglosc
 457           read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
 458      &       sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
 459           odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
 460           sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
 461      &     sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
 462         enddo
 463  1401 continue
 464         close (ientin)
 465         open (ientin,file=model_ki_angle,status='old')
 466         do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
 467           read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
 468      &      sigma_dih(ki,i+nnt-1)
 469           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
 470           sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
 471         enddo
 472  1402 continue
 473         close (ientin)
 474       enddo
 475       ii=0
 476       write (iout,*) "nnt",nnt," nct",nct
 477       do i=nnt,nct-2
 478         do j=i+2,nct
 479           ki=1
 480 c          write (iout,*) "i",i," j",j," constr_homology",constr_homology
 481           do while (ki.le.constr_homology .and.
 482      &        sigma_odl_temp(i,j,ki).le.0.0d0)
 483 c            write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
 484             ki=ki+1
 485           enddo
 486 c          write (iout,*) "ki",ki
 487           if (ki.gt.constr_homology) cycle
 488           ii=ii+1
 489           ires_homo(ii)=i
 490           jres_homo(ii)=j
 491           do ki=1,constr_homology
 492             odl(ki,ii)=odl_temp(i,j,ki)
 493             sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
 494           enddo
 495         enddo
 496       enddo
 497       lim_odl=ii
 498       if (constr_homology.gt.0) call homology_partition
 499 c Print restraints
 500       if (.not.lprn) return
 501       write (iout,*) "Distance restraints from templates"
 502       do ii=1,lim_odl
 503         write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
 504      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
 505       enddo
 506       write (iout,*) "Dihedral angle restraints from templates"
 507       do i=nnt,lim_dih
 508         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
 509      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 510       enddo
 511 c      write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
 512 c      write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
 513
 514       return
 515       end
 516