source/wham/src-M-SAXS-homology/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,tperm
  21       double precision x(maxvar)
  22       character*320 controlcard,ucase
  23       dimension itype_pdb(maxres)
  24       logical seq_comp
  25       double precision secprob(3,maxdih_constr),phihel,phibet
  26       call card_concat(controlcard,.true.)
  27       call reada(controlcard,'SCAL14',scal14,0.4d0)
  28       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
  29       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
  30       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
  31       r0_corr=cutoff_corr-delt_corr
  32 C     Bartek
  33       call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
  34       call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
  35       call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
  36       call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
  37       write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
  38      &  " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
  39       r0_corr=cutoff_corr-delt_corr
  40       call readi(controlcard,"NRES",nres,0)
  41       iscode=index(controlcard,"ONE_LETTER")
  42       if (nres.le.0) then
  43         write (iout,*) "Error: no residues in molecule"
  44         return1
  45       endif
  46       if (nres.gt.maxres) then
  47         write (iout,*) "Error: too many residues",nres,maxres
  48       endif
  49       write(iout,*) 'nres=',nres
  50 C Read sequence of the protein
  51       if (iscode.gt.0) then
  52         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  53       else
  54         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  55       endif
  56 C Convert sequence to numeric code
  57       do i=1,nres
  58         itype(i)=rescode(i,sequence(i),iscode)
  59       enddo
  60       write (iout,*) "Numeric code:"
  61       write (iout,'(20i4)') (itype(i),i=1,nres)
  62       do i=1,nres-1
  63 #ifdef PROCOR
  64         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
  65 #else
  66         if (itype(i).eq.ntyp1) then
  67 #endif
  68           itel(i)=0
  69 #ifdef PROCOR
  70         else if (iabs(itype(i+1)).ne.20) then
  71 #else
  72         else if (iabs(itype(i)).ne.20) then
  73 #endif
  74           itel(i)=1
  75         else
  76           itel(i)=2
  77         endif
  78       enddo
  79        write (iout,*) "ITEL"
  80        do i=1,nres-1
  81          write (iout,*) i,itype(i),itel(i)
  82        enddo
  83       call read_bridge
  84       nnt=1
  85       nct=nres
  86       call seq2chains(nres,itype,nchain,chain_length,chain_border,
  87      &  ireschain)
  88       write(iout,*) "nres",nres," nchain",nchain
  89       do i=1,nchain
  90         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
  91      &    chain_border(2,i)
  92       enddo
  93       call chain_symmetry(nchain,nres,itype,chain_border,
  94      &    chain_length,npermchain,tabpermchain)
  95       write(iout,*) "ireschain permutations"
  96       do i=1,nres
  97         write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
  98      &    ii=1,npermchain)
  99       enddo
 100       write(iout,*) "residue permutations"
 101       do i=1,nres
 102         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 103       enddo
 104
 105       if (itype(1).eq.ntyp1) nnt=2
 106       if (itype(nres).eq.ntyp1) nct=nct-1
 107       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 108 #ifdef DFA
 109 C     Juyong:READ init_vars
 110 C     Initialize variables!
 111 C     Juyong:READ read_info
 112 C     READ fragment information!!
 113 C     both routines should be in dfa.F file!!
 114
 115       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 116      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 117        write (iout,*) "Calling init_dfa_vars"
 118        call flush(iout)
 119        call init_dfa_vars
 120        write (iout,*) 'init_dfa_vars finished!'
 121        call flush(iout)
 122        call read_dfa_info
 123        write (iout,*) 'read_dfa_info finished!'
 124        call flush(iout)
 125       endif
 126 #endif
 127       if (with_dihed_constr) then
 128
 129       read (inp,*) ndih_constr
 130       write (iout,*) "ndih_constr",ndih_constr
 131       if (ndih_constr.gt.0) then
 132          raw_psipred=.false.
 133 C        read (inp,*) ftors
 134 C        write (iout,*) 'FTORS',ftors
 135         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 136      &   i=1,ndih_constr)
 137         write (iout,*)
 138      &   'There are',ndih_constr,' restraints on gamma angles.'
 139         do i=1,ndih_constr
 140           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 141      &    ftors(i)
 142         enddo
 143         do i=1,ndih_constr
 144           phi0(i)=deg2rad*phi0(i)
 145           drange(i)=deg2rad*drange(i)
 146         enddo
 147       else if (ndih_constr.lt.0) then
 148         raw_psipred=.true.
 149         call card_concat(controlcard,.true.)
 150         call reada(controlcard,"PHIHEL",phihel,50.0D0)
 151         call reada(controlcard,"PHIBET",phibet,180.0D0)
 152         call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
 153         call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
 154         call reada(controlcard,"WDIHC",wdihc,0.591d0)
 155         write (iout,*) "Weight of the dihedral restraint term",wdihc
 156         read(inp,'(9x,3f7.3)')
 157      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 158         write (iout,*) "The secprob array"
 159         do i=nnt,nct
 160           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 161         enddo
 162         ndih_constr=0
 163         do i=nnt+3,nct
 164           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 165      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 166             ndih_constr=ndih_constr+1
 167             idih_constr(ndih_constr)=i
 168             sumv=0.0d0
 169             do j=1,3
 170               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
 171               sumv=sumv+vpsipred(j,ndih_constr)
 172             enddo
 173             do j=1,3
 174               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
 175             enddo
 176             phibound(1,ndih_constr)=phihel*deg2rad
 177             phibound(2,ndih_constr)=phibet*deg2rad
 178             sdihed(1,ndih_constr)=sigmahel*deg2rad
 179             sdihed(2,ndih_constr)=sigmabet*deg2rad
 180           endif
 181         enddo
 182         write (iout,*)
 183      &   'There are',ndih_constr,
 184      &   ' bimodal restraints on gamma angles.'
 185         do i=1,ndih_constr
 186           write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
 187      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
 188      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
 189      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
 190      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
 191      &      (vpsipred(j,i),j=1,3)
 192         enddo
 193
 194       endif
 195
 196       endif
 197       if (with_theta_constr) then
 198 C with_theta_constr is keyword allowing for occurance of theta constrains
 199       read (inp,*) ntheta_constr
 200 C ntheta_constr is the number of theta constrains
 201       if (ntheta_constr.gt.0) then
 202 C        read (inp,*) ftors
 203         read (inp,*) (itheta_constr(i),theta_constr0(i),
 204      &  theta_drange(i),for_thet_constr(i),
 205      &  i=1,ntheta_constr)
 206 C the above code reads from 1 to ntheta_constr
 207 C itheta_constr(i) residue i for which is theta_constr
 208 C theta_constr0 the global minimum value
 209 C theta_drange is range for which there is no energy penalty
 210 C for_thet_constr is the force constant for quartic energy penalty
 211 C E=k*x**4
 212 C        if(me.eq.king.or..not.out1file)then
 213          write (iout,*)
 214      &   'There are',ntheta_constr,' constraints on phi angles.'
 215          do i=1,ntheta_constr
 216           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 217      &    theta_drange(i),
 218      &    for_thet_constr(i)
 219          enddo
 220 C        endif
 221         do i=1,ntheta_constr
 222           theta_constr0(i)=deg2rad*theta_constr0(i)
 223           theta_drange(i)=deg2rad*theta_drange(i)
 224         enddo
 225 C        if(me.eq.king.or..not.out1file)
 226 C     &   write (iout,*) 'FTORS',ftors
 227 C        do i=1,ntheta_constr
 228 C          ii = itheta_constr(i)
 229 C          thetabound(1,ii) = phi0(i)-drange(i)
 230 C          thetabound(2,ii) = phi0(i)+drange(i)
 231 C        enddo
 232       endif ! ntheta_constr.gt.0
 233       endif! with_theta_constr
 234       if (constr_homology.gt.0) then
 235 c       write (iout,*) "About to call read_constr_homology"
 236 c       call flush(iout)
 237         call read_constr_homology
 238 c       write (iout,*) "Exit read_constr_homology"
 239 c       call flush(iout)
 240         if (indpdb.gt.0 .or. pdbref) then
 241           do i=1,2*nres
 242             do j=1,3
 243               c(j,i)=crefjlee(j,i)
 244               cref(j,i)=crefjlee(j,i)
 245             enddo
 246           enddo
 247         endif
 248 #ifdef DEBUG
 249         write (iout,*) "Array C"
 250         do i=1,nres
 251           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
 252      &      (c(j,i+nres),j=1,3)
 253         enddo
 254         write (iout,*) "Array Cref"
 255         do i=1,nres
 256           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
 257      &      (cref(j,i+nres),j=1,3)
 258         enddo
 259 #endif
 260 #ifdef DEBUG
 261        call int_from_cart1(.false.)
 262        call sc_loc_geom(.false.)
 263        do i=1,nres
 264          thetaref(i)=theta(i)
 265          phiref(i)=phi(i)
 266          write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
 267        enddo
 268        do i=1,nres-1
 269          do j=1,3
 270            dc(j,i)=c(j,i+1)-c(j,i)
 271            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 272          enddo
 273        enddo
 274        do i=2,nres-1
 275          do j=1,3
 276            dc(j,i+nres)=c(j,i+nres)-c(j,i)
 277            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 278          enddo
 279        enddo
 280 #endif
 281       else
 282         homol_nset=0
 283       endif
 284
 285       call setup_var
 286       call init_int_table
 287       if (ns.gt.0) then
 288         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 289      &  ' disulfide-bridging cysteines.'
 290         write (iout,'(20i4)') (iss(i),i=1,ns)
 291        if (dyn_ss) then
 292           write(iout,*)"Running with dynamic disulfide-bond formation"
 293        else
 294         write (iout,'(/a/)') 'Pre-formed links are:'
 295         do i=1,nss
 296           i1=ihpb(i)-nres
 297           i2=jhpb(i)-nres
 298           it1=itype(i1)
 299           it2=itype(i2)
 300          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 301      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 302      &    dhpb(i),ebr,forcon(i)
 303         enddo
 304       endif
 305       endif
 306       write (iout,'(a)')
 307       if (ns.gt.0.and.dyn_ss) then
 308           do i=nss+1,nhpb
 309             ihpb(i-nss)=ihpb(i)
 310             jhpb(i-nss)=jhpb(i)
 311             forcon(i-nss)=forcon(i)
 312             dhpb(i-nss)=dhpb(i)
 313           enddo
 314           nhpb=nhpb-nss
 315           nss=0
 316           call hpb_partition
 317           do i=1,ns
 318             dyn_ss_mask(iss(i))=.true.
 319           enddo
 320       endif
 321       write (iout,*) "calling read_saxs_consrtr",nsaxs
 322       if (nsaxs.gt.0) call read_saxs_constr
 323
 324       return
 325       end
 326 c-----------------------------------------------------------------------------
 327       logical function seq_comp(itypea,itypeb,length)
 328       implicit none
 329       integer length,itypea(length),itypeb(length)
 330       integer i
 331       do i=1,length
 332         if (itypea(i).ne.itypeb(i)) then
 333           seq_comp=.false.
 334           return
 335         endif
 336       enddo
 337       seq_comp=.true.
 338       return
 339       end
 340 c-----------------------------------------------------------------------------
 341       subroutine read_bridge
 342 C Read information about disulfide bridges.
 343       implicit real*8 (a-h,o-z)
 344       include 'DIMENSIONS'
 345       include 'DIMENSIONS.ZSCOPT'
 346       include 'COMMON.IOUNITS'
 347       include 'COMMON.GEO'
 348       include 'COMMON.VAR'
 349       include 'COMMON.INTERACT'
 350       include 'COMMON.NAMES'
 351       include 'COMMON.CHAIN'
 352       include 'COMMON.FFIELD'
 353       include 'COMMON.SBRIDGE'
 354 C Read bridging residues.
 355       read (inp,*) ns,(iss(i),i=1,ns)
 356       print *,'ns=',ns
 357       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 358 C Check whether the specified bridging residues are cystines.
 359       do i=1,ns
 360         if (itype(iss(i)).ne.1) then
 361           write (iout,'(2a,i3,a)')
 362      &   'Do you REALLY think that the residue ',
 363      &    restyp(itype(iss(i))),i,
 364      &   ' can form a disulfide bridge?!!!'
 365           write (*,'(2a,i3,a)')
 366      &   'Do you REALLY think that the residue ',
 367      &    restyp(itype(iss(i))),i,
 368      &   ' can form a disulfide bridge?!!!'
 369          stop
 370         endif
 371       enddo
 372 C Read preformed bridges.
 373       if (ns.gt.0) then
 374       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 375       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 376       if (nss.gt.0) then
 377         nhpb=nss
 378 C Check if the residues involved in bridges are in the specified list of
 379 C bridging residues.
 380         do i=1,nss
 381           do j=1,i-1
 382             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 383      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 384               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 385      &      ' contains residues present in other pairs.'
 386               write (*,'(a,i3,a)') 'Disulfide pair',i,
 387      &      ' contains residues present in other pairs.'
 388               stop
 389             endif
 390           enddo
 391           do j=1,ns
 392             if (ihpb(i).eq.iss(j)) goto 10
 393           enddo
 394           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 395    10     continue
 396           do j=1,ns
 397             if (jhpb(i).eq.iss(j)) goto 20
 398           enddo
 399           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 400    20     continue
 401           dhpb(i)=dbr
 402           forcon(i)=fbr
 403         enddo
 404         do i=1,nss
 405           ihpb(i)=ihpb(i)+nres
 406           jhpb(i)=jhpb(i)+nres
 407         enddo
 408       endif
 409       endif
 410       return
 411       end
 412 c------------------------------------------------------------------------------
 413       subroutine read_angles(kanal,iscor,energ,iprot,*)
 414       implicit real*8 (a-h,o-z)
 415       include 'DIMENSIONS'
 416       include 'DIMENSIONS.ZSCOPT'
 417       include 'COMMON.INTERACT'
 418       include 'COMMON.SBRIDGE'
 419       include 'COMMON.GEO'
 420       include 'COMMON.VAR'
 421       include 'COMMON.CHAIN'
 422       include 'COMMON.IOUNITS'
 423       character*80 lineh
 424       read(kanal,'(a80)',end=10,err=10) lineh
 425       read(lineh(:5),*,err=8) ic
 426       read(lineh(6:),*,err=8) energ
 427       goto 9
 428     8 ic=1
 429       print *,'error, assuming e=1d10',lineh
 430       energ=1d10
 431       nss=0
 432     9 continue
 433       read(lineh(18:),*,end=10,err=10) nss
 434       IF (NSS.LT.9) THEN
 435         read (lineh(20:),*,end=10,err=10)
 436      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 437       ELSE
 438         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 439         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 440      &    I=9,NSS),iscor
 441       ENDIF
 442 c      print *,"energy",energ," iscor",iscor
 443       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 444       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 445       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 446       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 447       do i=1,nres
 448         theta(i)=deg2rad*theta(i)
 449         phi(i)=deg2rad*phi(i)
 450         alph(i)=deg2rad*alph(i)
 451         omeg(i)=deg2rad*omeg(i)
 452       enddo
 453       return
 454    10 return1
 455       end
 456 c-------------------------------------------------------------------------------
 457       subroutine read_saxs_constr
 458       implicit real*8 (a-h,o-z)
 459       include 'DIMENSIONS'
 460       include 'DIMENSIONS.ZSCOPT'
 461       include 'DIMENSIONS.FREE'
 462 #ifdef MPI
 463       include 'mpif.h'
 464 #endif
 465       include 'COMMON.CONTROL'
 466       include 'COMMON.CHAIN'
 467       include 'COMMON.IOUNITS'
 468       include 'COMMON.SBRIDGE'
 469       double precision cm(3)
 470 c      read(inp,*) nsaxs
 471       write (iout,*) "Calling read_saxs nsaxs",nsaxs
 472       call flush(iout)
 473       if (saxs_mode.eq.0) then
 474 c SAXS distance distribution
 475       do i=1,nsaxs
 476         read(inp,*) distsaxs(i),Psaxs(i)
 477       enddo
 478       Cnorm = 0.0d0
 479       do i=1,nsaxs
 480         Cnorm = Cnorm + Psaxs(i)
 481       enddo
 482       write (iout,*) "Cnorm",Cnorm
 483       do i=1,nsaxs
 484         Psaxs(i)=Psaxs(i)/Cnorm
 485       enddo
 486       write (iout,*) "Normalized distance distribution from SAXS"
 487       do i=1,nsaxs
 488         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
 489       enddo
 490       Wsaxs0=0.0d0
 491       do i=1,nsaxs
 492         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
 493       enddo
 494       write (iout,*) "Wsaxs0",Wsaxs0
 495       else
 496 c SAXS "spheres".
 497       do i=1,nsaxs
 498         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
 499       enddo
 500       do j=1,3
 501         cm(j)=0.0d0
 502       enddo
 503       do i=1,nsaxs
 504         do j=1,3
 505           cm(j)=cm(j)+Csaxs(j,i)
 506         enddo
 507       enddo
 508       do j=1,3
 509         cm(j)=cm(j)/nsaxs
 510       enddo
 511       do i=1,nsaxs
 512         do j=1,3
 513           Csaxs(j,i)=Csaxs(j,i)-cm(j)
 514         enddo
 515       enddo
 516       write (iout,*) "SAXS sphere coordinates"
 517       do i=1,nsaxs
 518         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
 519       enddo
 520       endif
 521       return
 522       end