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