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