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       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       call init_int_table
 264       if (ns.gt.0) then
 265         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 266      &  ' disulfide-bridging cysteines.'
 267         write (iout,'(20i4)') (iss(i),i=1,ns)
 268        if (dyn_ss) then
 269           write(iout,*)"Running with dynamic disulfide-bond formation"
 270        else
 271         write (iout,'(/a/)') 'Pre-formed links are:'
 272         do i=1,nss
 273           i1=ihpb(i)-nres
 274           i2=jhpb(i)-nres
 275           it1=itype(i1)
 276           it2=itype(i2)
 277          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 278      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 279      &    dhpb(i),ebr,forcon(i)
 280         enddo
 281       endif
 282       endif
 283       write (iout,'(a)')
 284       if (ns.gt.0.and.dyn_ss) then
 285           do i=nss+1,nhpb
 286             ihpb(i-nss)=ihpb(i)
 287             jhpb(i-nss)=jhpb(i)
 288             forcon(i-nss)=forcon(i)
 289             dhpb(i-nss)=dhpb(i)
 290           enddo
 291           nhpb=nhpb-nss
 292           nss=0
 293           call hpb_partition
 294           do i=1,ns
 295             dyn_ss_mask(iss(i))=.true.
 296           enddo
 297       endif
 298       write (iout,*) "calling read_saxs_consrtr",nsaxs
 299       if (nsaxs.gt.0) call read_saxs_constr
 300
 301       return
 302       end
 303 c-----------------------------------------------------------------------------
 304       logical function seq_comp(itypea,itypeb,length)
 305       implicit none
 306       integer length,itypea(length),itypeb(length)
 307       integer i
 308       do i=1,length
 309         if (itypea(i).ne.itypeb(i)) then
 310           seq_comp=.false.
 311           return
 312         endif
 313       enddo
 314       seq_comp=.true.
 315       return
 316       end
 317 c-----------------------------------------------------------------------------
 318       subroutine read_bridge
 319 C Read information about disulfide bridges.
 320       implicit real*8 (a-h,o-z)
 321       include 'DIMENSIONS'
 322       include 'DIMENSIONS.ZSCOPT'
 323       include 'COMMON.IOUNITS'
 324       include 'COMMON.GEO'
 325       include 'COMMON.VAR'
 326       include 'COMMON.INTERACT'
 327       include 'COMMON.NAMES'
 328       include 'COMMON.CHAIN'
 329       include 'COMMON.FFIELD'
 330       include 'COMMON.SBRIDGE'
 331 C Read bridging residues.
 332       read (inp,*) ns,(iss(i),i=1,ns)
 333       print *,'ns=',ns
 334       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 335 C Check whether the specified bridging residues are cystines.
 336       do i=1,ns
 337         if (itype(iss(i)).ne.1) then
 338           write (iout,'(2a,i3,a)')
 339      &   'Do you REALLY think that the residue ',
 340      &    restyp(itype(iss(i))),i,
 341      &   ' can form a disulfide bridge?!!!'
 342           write (*,'(2a,i3,a)')
 343      &   'Do you REALLY think that the residue ',
 344      &    restyp(itype(iss(i))),i,
 345      &   ' can form a disulfide bridge?!!!'
 346          stop
 347         endif
 348       enddo
 349 C Read preformed bridges.
 350       if (ns.gt.0) then
 351       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 352       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 353       if (nss.gt.0) then
 354         nhpb=nss
 355 C Check if the residues involved in bridges are in the specified list of
 356 C bridging residues.
 357         do i=1,nss
 358           do j=1,i-1
 359             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 360      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 361               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 362      &      ' contains residues present in other pairs.'
 363               write (*,'(a,i3,a)') 'Disulfide pair',i,
 364      &      ' contains residues present in other pairs.'
 365               stop
 366             endif
 367           enddo
 368           do j=1,ns
 369             if (ihpb(i).eq.iss(j)) goto 10
 370           enddo
 371           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 372    10     continue
 373           do j=1,ns
 374             if (jhpb(i).eq.iss(j)) goto 20
 375           enddo
 376           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 377    20     continue
 378           dhpb(i)=dbr
 379           forcon(i)=fbr
 380         enddo
 381         do i=1,nss
 382           ihpb(i)=ihpb(i)+nres
 383           jhpb(i)=jhpb(i)+nres
 384         enddo
 385       endif
 386       endif
 387       return
 388       end
 389 c------------------------------------------------------------------------------
 390       subroutine read_angles(kanal,iscor,energ,iprot,*)
 391       implicit real*8 (a-h,o-z)
 392       include 'DIMENSIONS'
 393       include 'DIMENSIONS.ZSCOPT'
 394       include 'COMMON.INTERACT'
 395       include 'COMMON.SBRIDGE'
 396       include 'COMMON.GEO'
 397       include 'COMMON.VAR'
 398       include 'COMMON.CHAIN'
 399       include 'COMMON.IOUNITS'
 400       character*80 lineh
 401       read(kanal,'(a80)',end=10,err=10) lineh
 402       read(lineh(:5),*,err=8) ic
 403       read(lineh(6:),*,err=8) energ
 404       goto 9
 405     8 ic=1
 406       print *,'error, assuming e=1d10',lineh
 407       energ=1d10
 408       nss=0
 409     9 continue
 410       read(lineh(18:),*,end=10,err=10) nss
 411       IF (NSS.LT.9) THEN
 412         read (lineh(20:),*,end=10,err=10)
 413      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 414       ELSE
 415         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 416         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 417      &    I=9,NSS),iscor
 418       ENDIF
 419 c      print *,"energy",energ," iscor",iscor
 420       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 421       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 422       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 423       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 424       do i=1,nres
 425         theta(i)=deg2rad*theta(i)
 426         phi(i)=deg2rad*phi(i)
 427         alph(i)=deg2rad*alph(i)
 428         omeg(i)=deg2rad*omeg(i)
 429       enddo
 430       return
 431    10 return1
 432       end
 433 c-------------------------------------------------------------------------------
 434       subroutine read_saxs_constr
 435       implicit real*8 (a-h,o-z)
 436       include 'DIMENSIONS'
 437       include 'DIMENSIONS.ZSCOPT'
 438       include 'DIMENSIONS.FREE'
 439 #ifdef MPI
 440       include 'mpif.h'
 441 #endif
 442       include 'COMMON.CONTROL'
 443       include 'COMMON.CHAIN'
 444       include 'COMMON.IOUNITS'
 445       include 'COMMON.SBRIDGE'
 446       include 'COMMON.SAXS'
 447       double precision cm(3)
 448 c      read(inp,*) nsaxs
 449       write (iout,*) "Calling read_saxs nsaxs",nsaxs
 450       call flush(iout)
 451       if (saxs_mode.eq.0) then
 452 c SAXS distance distribution
 453       do i=1,nsaxs
 454         read(inp,*) distsaxs(i),Psaxs(i)
 455       enddo
 456       Cnorm = 0.0d0
 457       do i=1,nsaxs
 458         Cnorm = Cnorm + Psaxs(i)
 459       enddo
 460       write (iout,*) "Cnorm",Cnorm
 461       do i=1,nsaxs
 462         Psaxs(i)=Psaxs(i)/Cnorm
 463       enddo
 464       write (iout,*) "Normalized distance distribution from SAXS"
 465       do i=1,nsaxs
 466         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
 467       enddo
 468       Wsaxs0=0.0d0
 469       do i=1,nsaxs
 470         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
 471       enddo
 472       write (iout,*) "Wsaxs0",Wsaxs0
 473       else
 474 c SAXS "spheres".
 475       do i=1,nsaxs
 476         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
 477       enddo
 478       do j=1,3
 479         cm(j)=0.0d0
 480       enddo
 481       do i=1,nsaxs
 482         do j=1,3
 483           cm(j)=cm(j)+Csaxs(j,i)
 484         enddo
 485       enddo
 486       do j=1,3
 487         cm(j)=cm(j)/nsaxs
 488       enddo
 489       do i=1,nsaxs
 490         do j=1,3
 491           Csaxs(j,i)=Csaxs(j,i)-cm(j)
 492         enddo
 493       enddo
 494       write (iout,*) "SAXS sphere coordinates"
 495       do i=1,nsaxs
 496         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
 497       enddo
 498       endif
 499       return
 500       end