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