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