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